Is it possible to find complex eigenvalues with QR decomposition?

Question

I wonder if it's possible to find the complex eigenvalues with QR decomposition. I can find the real eigenvalues with QR just by doing this.

function [R] =findEig(A,k)
for i=1:k
    [Q,R] =qr(A);
    A=R*Q;
end
end

Let's say I have a matrix $A$

>> A
A =

    2.00000    3.00000    1.00000    0.50000    4.00000
    4.00000    5.00000    7.00000    0.10000    1.00000
    5.00000    3.00000    6.00000   19.20000    9.00000
    1.00000    4.00000    1.00000    4.00000    7.00000
    3.00000    1.00000    6.00000    2.00000    6.00000

Then I can run the function where I want to do 10 iterations.

>> R = findEig(A, 10)
R =

  -21.35929   -9.38466    6.85281   -4.68212   -2.85981
   0.00000    6.67519   -3.58144    3.60920    5.73361
   0.00000    0.00000    7.31262   -0.27010   -3.36928
   0.00000    0.00000    0.00000   -4.32428   -0.28359
   0.00000    0.00000    0.00000    0.00000    1.25801

>> eig(A)
ans =

   21.3589 +  0.0000i
   -1.9810 +  6.6826i
   -1.9810 -  6.6826i
    1.2580 +  0.0000i
    4.3450 +  0.0000i

As you can see. On the diagonal of matrix $R$, we can see that it has the same eigenvalues as the eigenvalues above, but $R$ does not show the complex eigenvalues. The more iteration I do, the better result I will get. But not for the complex eigenvalues.

So is there possible to find the complex eigenvalues with QR?

Answer 1

The algorithm that you've tried to use can't compute complex eigenvalues because the QR factorization of A will always be real and the next iteration will have a real A matrix. It is true that $A$ in this algorithm will converge (slowly) to a quasi-upper triangular real Schur matrix. The output $A$ will (by the nature of the similarity transformations) have the same eigenvalues as the input $A$, and you can recover the complex eigenvalues from the this quasi-triangular matrix.

Practical implementations are more sophisticated. The implicitly shifted QR algorithm with a double shift is one way to make this work. This is discussed in many textbooks on numerical linear algebra. I'd particularly recommend Fundamentals of Matrix Computations, 3rd ed., by David S. Watkins.

Your procedure is also incorrect in that it returns the last $R$ matrix rather than the last $A$ matrix. If you modify the procedure to return $A$, you'll find that the $A$ does converge to a quasi-triangular real Schur matrix with the correct complex eigenvalues.

The corrected algorithm is:

function S=myqr(A,k)
for i=1:k
  [Q,R]=qr(A);
  A=R*Q;
end
S=A;

Let's see how this works in comparison with the schur() command in MATLAB.

>> schur(A)

ans =

   21.3589   -5.6678   10.0809   -2.3429   -5.0900
         0   -1.9810    8.7273   -6.2583   -3.8582
         0   -5.1169   -1.9810   -0.3539   -1.8999
         0         0         0    1.2580   -0.3232
         0         0         0         0    4.3450
>> eig(schur(A))

ans =

  21.3589 + 0.0000i
  -1.9810 + 6.6826i
  -1.9810 - 6.6826i
   1.2580 + 0.0000i
   4.3450 + 0.0000i

The two complex eigenvalues can be computed directly from the 2x2 block in rows 2-3, columns 2-3.

>> S=schur(A); S(2:3,2:3); eig(S(2:3,2:3))

ans =

  -1.9810 + 6.6826i
  -1.9810 - 6.6826i

You can get the same result using the simple-minded QR algorithm. It's easy to confirm (using eig()) that S is converging to a quasi-triangular real Schur matrix with the correct eigenvalues.

>> S=myqr(A,100)

S =

   21.3589   -6.1095    9.8195    4.8184   -2.8601
    0.0000   -2.1410    8.7202    3.1002   -6.5951
   -0.0000   -5.1240   -1.8209    1.9923   -0.8432
    0.0000   -0.0000    0.0000    4.3450   -0.3232
    0.0000   -0.0000    0.0000    0.0000    1.2580

>> eig(S)

ans =

  21.3589 + 0.0000i
  -1.9810 + 6.6826i
  -1.9810 - 6.6826i
   4.3450 + 0.0000i
   1.2580 + 0.0000i