I'm wondering if the Shur decompositions are ordered in any way. From the Shur decomposition development thread (url is ?f=74&t=106641) it seems the algorithm being used is based on Stewarts fortran implementation, for which I've read that 2x2 blocks are ordered so that eigenvalues appear in descending order of magnitude.

I'm working through a paper, "A Schur Method for Solving Algebraic Riccati Equations" by Alan Laub, which says that the real Shur form can be arbitrarily reordered via orthogonal similarities, citing another paper "Ill-conditioned eigensystems and the computation of the Jordan canonical form" by G. Golub and J. Wilkinson. This second paper was a bit beyond my level of understanding so I couldn't quite find the appropriate section describing how to generate such a transformation. I'd like to transform the Schur decomposition such that 2x2 blocks are ordered so that eigenvalues with negative real part occur before eigen values with positive real part in order to apply the ARE solution method from Laub's paper. Does anyone know how that can be done?

Thanks for any suggestions!

After digging into the code for RealSchur.h for a couple of days, it appears that the algorithm implemented is more or less that of Stewart's HQR3 ("HQR3 and EXCHNG: Fortran Subroutines for Calculating and Ordering the Eigenvalues of a Real Upper Hessenberg Matrix" TOMS '76). I have cloned the RealSchur class as OrderedRealSchur and have had modest success in implementing Stewart's EXCHNG from his HQR3 fortran code published with his paper. I've implemented it as a method of OrderedRealSchur. The code for this method is included below.

The original fortran code is in commented blocks, followed by my "translation". I have only tested the 2x2 <-> 2x2 exchange thus far, and it appears to work. I think the trick is understanding that in Stewarts fortran code (X,Y,W) are 'shiftInfo' in Eigen, and Stewarts (P,Q,R) are 'firstHouseholderVector' in Eigen.

exchange(L,b1,b2) alters the decomposition by a similarity transformation which "swaps" a pair of 1x1 or 2x2 blocks on the diagonal of the schur decomposition, along with their eigen values. (1x1: one real eigen value, 2x2: complex eigen value pair... but you probably knew that). The method takes three parameters. L is the starting row of the blocks to swap. b1 is the size of the first block (1 or 2) and b2 is the size of the second block (1 or 2). The code for a 1x1 block exchanged with a 1x1 block looks like some row/column operations that can be written more clearly using the Eigen API. As I investigate and understand more I will try to clean those up.

Again, I have only tested the third case in the method thus far (2x2 <-> 2x2 swap).

Code: Select all

template<typename MatrixType>
OrderedRealSchur<MatrixType>&
  OrderedRealSchur<MatrixType>::exchange(Index L, Index b1, Index b2)
{
  Scalar P,Q,R,S,W,X,Y; //,Z;
  MatrixType& A = m_matT;
  MatrixType& V = m_matU;
  const Index N = A.rows();
  const Scalar EPS = NumTraits<Scalar>::epsilon();

  assert( L+2 < N );
  assert( b1 == 1 || b1 == 2 );
  assert( b2 == 1 || b2 == 2 );

  // exchange 1x1 and 1x1 blocks
  if( b1 == 1 && b2 == 1 ){
//      L1 = L + 1
//      Q = A(L+1,L+1) - A(L,L)
//      P = A(L,L+1)
//      R = AMAX1(ABS(P),ABS(Q))
//      IF (R.EQ.0.) RETURN
//      P = P/R
//      Q = Q/R
//      R = SQRT(P**2+Q**2)
//      P = P/R
//      Q = Q/R
    Q = A(L+1,L+1) - A(L,L);
    P = A(L,L+1);
    R = std::max( internal::abs(Q), internal::abs(P) );
    if( R == Scalar(0) ){
      m_info = InvalidInput;
      return *this;
    }
    P = P/R;
    Q = Q/R;
    R = internal::sqrt(P*P + Q*Q);
    P = P/R;
    Q = Q/R;

//      DO 10 J=L,N
//        S = P*A(L,J) + Q*A(L+1,J)
//        A(L+1,J) = P*A(L+1,J) - Q*A(L,J)
//        A(L,J) = S
//   10 CONTINUE
    for (int J = L; J < N; J++) {
      S = P*A(L,J) + Q*A(L+1,J);
      A(L+1,J) = P*A(L+1,J) - Q*A(L,J);
      A(L,J) = S;
    }

//      DO 20 I=1,L1
//        S = P*A(I,L) + Q*A(I,L+1)
//        A(I,L+1) = P*A(I,L+1) - Q*A(I,L)
//        A(I,L) = S
//   20 CONTINUE
    for(int I=0; I < L+1; I++){
      S = P*A(I,L) + Q*A(I,L+1);
      A(I,L+1) = P*A(I,L+1) - Q*A(I,L);
      A(I,L) = S;
    }

//      DO 30 I=1,N
//        S = P*V(I,L) + Q*V(I,L+1)
//        V(I,L+1) = P*V(I,L+1) - Q*V(I,L)
//        V(I,L) = S
//   30 CONTINUE
    for(int I=0; I < N; I++){
      S = P*V(I,L) + Q*V(I,L+1);
      V(I,L+1) = P*V(I,L+1) - Q*V(I,L);
      V(I,L) = S;
    }

//      A(L+1,L) = 0.
//      RETURN
    A(L+1,L) = 0;
    return *this;

// exchange 1x1 and 2x2 blocks
//   40 CONTINUE
  } else if( b1 == 1 && b2 == 2 ) {
//      X = A(L,L)
//      P = 1.
//      Q = 1.
//      R = 1.
//      CALL QRSTEP(A, V, P, Q, R, L, L+2, N, NA, NV)
    X = A(L,L);
    P = Scalar(1);
    Q = Scalar(1);
    R = Scalar(1);

    Vector3s firstHouseholderVector(P,Q,R);
    Scalar* workspace = &m_workspaceVector.coeffRef(0);
    performFrancisQRStep(L,L,L+2,true,firstHouseholderVector,workspace);

//      IT = 0
//   50 IT = IT + 1
//      IF (IT.LE.30) GO TO 60
//      FAIL = .TRUE.
//      RETURN
//   60 CONTINUE
//      P = A(L,L) - X
//      Q = A(L+1,L)
//      R = 0.
//      CALL QRSTEP(A, V, P, Q, R, L, L+2, N, NA, NV)
//      IF (ABS(A(L+2,L+1)).GT.EPS*(ABS(A(L+1,L+1))+ABS(A(L+2,L+2))))
//     * GO TO 50
//      A(L+2,L+1) = 0.
//      RETURN
    int IT = 0;
    for(; IT < m_maxIterations*N; IT++ ){
      P = A(L,L) - X;
      Q = A(L+1,L);
      R = Scalar(0);
      firstHouseholderVector = Vector3s(P,Q,R);
      performFrancisQRStep(L,L,L+2,true,firstHouseholderVector,workspace);

      if (internal::abs(A(L+2,L+1)) <=
          EPS * (internal::abs(A(L+1,L+1)) +
                 internal::abs(A(L+2,L+2)))) {
        break;
      }
    }
    if( IT == m_maxIterations*N ){
        m_info = NoConvergence;
        return *this;
    }
    A(L+2,L+1) = 0;

// exchange a 2x2 block with a 1x1 or 2x2 block
//   70 CONTINUE
  } else if( b1 == 2 ) {
//      M = L + 2
//      IF (B2.EQ.2) M = M + 1
//      X = A(L+1,L+1)
//      Y = A(L,L)
//      W = A(L+1,L)*A(L,L+1)
//      P = 1.
//      Q = 1.
//      R = 1.
//      CALL QRSTEP(A, V, P, Q, R, L, M, N, NA, NV)
    Index M = L + 2;  // index of the last row in the block
    if( b2 == 2 ) M = M + 1;

    // Note: this work is normally done in computeShift but computShift
    // uses A(M-1,M-1) through A(M,M)
    X = A(L+1,L+1);         // lower right diagonal entry
    Y = A(L,L);             // upper left diagonal entry
    W = A(L+1,L)*A(L,L+1);  // off diagonal entries

    // note: eigen values of the first 2x2 block are solutions to
    // (z -X)*(z-Y) - W = 0
    // z^2 -(X+Y)z -XY - W = 0
    // z = (X+Y) +- sqrt( (X+Y)^2 + 4*(XY+W) ) / 2
    P = Scalar(1);
    Q = Scalar(1);
    R = Scalar(1);

    Vector3s shiftInfo(X,Y,W);
    Vector3s firstHouseholderVector(P,Q,R);
    Scalar* workspace = &m_workspaceVector.coeffRef(0);
    performFrancisQRStep(L,L,M,true,firstHouseholderVector,workspace);

//      IT = 0
//   80 IT = IT + 1
//      IF (IT.LE.30) GO TO 90
//      FAIL = .TRUE.
//      RETURN
//   90 CONTINUE
//      Z = A(L,L)
//      R = X - Z
//      S = Y - Z
//      P = (R*S-W)/A(L+1,L) + A(L,L+1)
//      Q = A(L+1,L+1) - Z - R - S
//      R = A(L+2,L+1)
//      S = ABS(P) + ABS(Q) + ABS(R)
//      P = P/S
//      Q = Q/S
//      R = R/S
//      CALL QRSTEP(A, V, P, Q, R, L, M, N, NA, NV)
//      IF (ABS(A(M-1,M-2)).GT.EPS*(ABS(A(M-1,M-1))+ABS(A(M-2,M-2))))
//      GO TO 80
    int IT = 0;
    for(; IT < m_maxIterations*N; IT++ ){
      // Note: all this fortran code is done in initFrancisQRStep,
      // in the eigen version shiftInfo = [X Y W] and
      // firstHouseholderVector = [P Q R]
      Index im;
      initFrancisQRStep(L,M, shiftInfo,im,firstHouseholderVector);
      performFrancisQRStep(L,im,M,true,firstHouseholderVector,workspace);

      if (internal::abs(A(M-1,M-2)) <=
          EPS * (internal::abs(A(M-1,M-1)) +
                 internal::abs(A(M-2,M-2)))) {
        break;
      }
    }
    if( IT == m_maxIterations*N ){
      m_info = NoConvergence;
      return *this;
    }

//       A(M-1,M-2) = 0.
//       RETURN
    A(M-1,M-2) = Scalar(0);
    return *this;


  } else {
    m_info = InvalidInput;
  }

  return *this;
}


Sorry for not replying so far, nice work anyway!

Regarding the upstream integration , I think that instead of adding a new class we could have a runtime option to RealSchur or a new method for post ordering, or even both.

Regarding the code, the 1 to 1 case could be simplified by observing P,Q represent a rotation in the plane, so using a JacobiRotation would allow to replace the for loops by call to applyOnThe{Left,Right}(...).

ggael wrote:Sorry for not replying so far, nice work anyway!

No worries, and thanks for getting back to me.

ggael wrote:Regarding the upstream integration , I think that instead of adding a new class we could have a runtime option to RealSchur or a new method for post ordering, or even both.

Sounds good to me. I was thinking an interface where compute() takes an additional parameter which is a callable that compares eigenvalues from two blocks and returns true if the first set of (1 or 2) eigenvalues should be ordered before the second set. I like the idea of having both options though.

ggael wrote:Regarding the code, the 1 to 1 case could be simplified by observing P,Q represent a rotation in the plane, so using a JacobiRotation would allow to replace the for loops by call to applyOnThe{Left,Right}(...).

Thanks for the hint, I will create a test and work on simplifying it.

ggael wrote:Sorry for not replying so far, nice work anyway!

Regarding the upstream integration , I think that instead of adding a new class we could have a runtime option to RealSchur or a new method for post ordering, or even both.

Regarding the code, the 1 to 1 case could be simplified by observing P,Q represent a rotation in the plane, so using a JacobiRotation would allow to replace the for loops by call to applyOnThe{Left,Right}(...).

Re-ordering the Schur decomposition is also necessary for the algorithm we use for general matrix functions, so there are more possible applications. Currently, the algorithm is only implemented in the complex case, which mean that it only needs to care about the simple 1-by-1 case. The implementation in matrix_function_permute_schur() in unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h uses JacobiRotation as Gael suggested, so you may have a look at that (I did not try to optimize the implementation though).