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

Question

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

Answer 1

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       

Answer 2

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.

Answer 3

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.