# How does the QR algorithm applied to a real matrix returns complex eigenvalues?

I'm a noob into eigenvalues algorithms, but something call my attention. QR algorithm works with real/complex matrices producing real/complex eigenvalues. However, it can not produce complex eigenvalues from a real matrix. Here a simplistic example written in Julia and derivated from here and here:

using LinearAlgebra
A = [7 3 4 11 -9 -2;
-6 4 -5 7 1 12;
-1 -9 2 2 9 1;
-8 0 -1 5 0 8;
-4 3 -5 7 2 10;
6 1 4 -11 -7 -1]
M = copy(A)

for i=1:100
global M
Q,R = LinearAlgebra.qr(M);
M=R*Q;
end

display(diag(M))
display(eigvals(A))

6-element Array{Float64,1}:
-2.8415406888480472
8.675063708533656
3.658872985794657
6.3411270142053695
0.12201942568224483
3.0444575546321087
6-element Array{Complex{Float64},1}:
2.916761509842819 + 13.248032079355992im
2.916761509842819 - 13.248032079355992im
5.000000000000005 + 6.000000000000003im
5.000000000000005 - 6.000000000000003im
1.5832384901571723 + 1.4155521348117128im
1.5832384901571723 - 1.4155521348117128im


Defining matrix A as complex, with only real components, makes no difference.

My questions are :

• what is my conceptual misunderstanding on the subject ?
• what step am I doing wrong ?
• and how to fix it ?

Thank you

• Regarding the answer of Christian Clason, you can for instance check that the lowest $2×2$ block has a trace $0.12201942568224483+3.0444575546321087=3.1664769803143535$ which is also about two times the real part of the last eigenvalue pair $2*1.5832384901571723=3.1664769803143447$. The same relation should hold for the other two pairs, but the order relation there is not as guaranteed as for the last, which is always the smallest eigenvalue in the simple QR algorithm without shifts or other bells-and-whistles. – Lutz Lehmann Oct 23 '18 at 10:57
• Real matrix may have complex eigenvalues. Real and Symmetric matrix can only have real eigenvalues. math.stackexchange.com/questions/67304/… – R zu Oct 24 '18 at 0:33

## 3 Answers

In a nutshell, the QR algorithm applied to a matrix $$A$$ is an iterative procedure that converges to the real Schur decomposition: a unitary matrix $$Q$$ and a matrix $$R$$ in block upper triangular form (see below) such that $$A = QRQ^T$$. It follows that the columns of $$Q$$ are the eigenvectors (which are the principal objects that are computed!) and that $$R$$ has the same eigenvalues as $$A$$.

The key point is the block upper triangular form, which means that $$R = \begin{pmatrix} R_{11}&&*\\&\ddots\\0&&R_{mm} \end{pmatrix},$$ where $$R_{ii}$$ are real blocks of either

• size $$1\times 1$$, in which case $$R_{ii}$$ is a (real) eigenvalue of $$A$$, or
• size $$2\times 2$$, in which case $$R_{ii}$$ has a pair of complex conjugate eigenvalues of $$A$$ (such as $$2+i$$ and $$2-i$$).

Since you can compute eigenvalues of $$2\times 2$$ matrices analytically (as roots of a quadratic polynomial), it is a cheap step to extract the complex eigenvalues from the computed (approximation of) $$R$$ in the end -- and this is what eigvals does.

So your conceptual misunderstanding is the following: Every real $$n\times n$$ matrix has $$n$$ eigenvalues, but they don't have to be real (or distinct) -- look at Richard Zhang's now deleted example, $$A=\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$$, which has the eigenvalues $$\pm i$$. Only if the matrix is symmetric are the eigenvalues guaranteed to be real (and the matrix $$R$$ to be diagonal) -- so your code only works for symmetric matrices. If the input matrix is non-symmetric, you additionally have to extract the (complex) eigenvalues by identifying the $$2\times 2$$ blocks (e.g., by checking whether a subdiagonal element is greater than a tolerance) and if so, computing the eigenvalues by a formula.

This is a bit tedious, but if you're willing to cheat a bit and use eigenvalues for the $$2\times 2$$ blocks, the following modification of your code will do it:

using LinearAlgebra
A = [7 3 4 11 -9 -2;
-6 4 -5 7 1 12;
-1 -9 2 2 9 1;
-8 0 -1 5 0 8;
-4 3 -5 7 2 10;
6 1 4 -11 -7 -1]
M = copy(A)

for i=1:100
global M
Q,R = LinearAlgebra.qr(M);
M=R*Q;
end

eig = Complex{Float64}[]
let
i=1
N=size(M,1)
while i<N
if abs(M[i+1,i])<1e-10
append!(eig,M[i,i])
i+=1
else
append!(eig,eigvals(M[i:i+1,i:i+1]))
i+=2
end
end
if length(eig)<N
append!(eig,M[N:N,N:N])
end
end

• I don't think that I understand your point. In my example, matrix M that contains the eigenvalues. How I suppose to get the complex eigenvalues from it ? ( if someone could provide a source code would be better ) – Noel Araujo Oct 23 '18 at 13:51
• Although this is the way things are generally done in practice, it's an interesting exercise to do a random complex similarity transform to the matrix (keeping the eigenvalues constant while making the matrix complex) and then do the QR algorithm with complex numbers. You'll end up with the complex eigenvalues on the diagonal of $R$ and eliminate the 2x2 blocks. – Brian Borchers Oct 23 '18 at 14:25
• One advantage of this approach is that any reasonable way of finding the eigenvalues of the 2x2 blocks ensures that the complex eigenvalues come in complex conjugate pairs as theory requires. – Brian Borchers Oct 23 '18 at 14:26
• @NoelAraujo You should look at the full matrix M , not just the diagonal diag(M) (since M is neither diagonal nor triangular). Then you'll see the 2x2 blocks. You can check my claim with, e.g., eigvals(M[1:2,1:2]). (@BrianBorchers That's a neat trick!) – Christian Clason Oct 23 '18 at 15:21
• @KutalmisB If you want complex eigenvectors (rather than a real eigenbasis), the easiest way is probably Brian's approach to work in complex arithmetic from the start. Extracting complex eigenvectors from the real Schur factorization can be done but is trickier; you can see how LAPACK does it. – Christian Clason Mar 21 '19 at 11:54

Forget about the QR algorithm, and remember what eigenvalues are - they are the roots of the characteristic polynomial for the matrix (see e.g. https://en.wikipedia.org/wiki/Characteristic_polynomial). For a real matrix of order N this is a polynomial of order N with real coefficients. But real coefficients does not mean real roots necessarily, you may have complex conjugate pairs. Hence a general real matrix may have complex eigenvalues.

I've used this equation many times to test convergence. The 11 in the first row should be -11. Running the QR algorithm without shifts, just as the matrix is, for about 60 iterations results in a matrix with the eigenvalues showing almost exactly on the diagonal , even the complex ones. The 5+-6i are in the first two lines, the 4 and 3 in the next two and finally the 1+-21 in the bottom right block. Don't think the complex values usually appear this clearly without computing the 2x2 matrix eigenvalues.