QR algorithm for eigenvalues

QR algorithm for eigenvalues#

import numpy as np
from numpy.linalg import qr, eig
from scipy.linalg import hessenberg

Pure QR algorithm#

# Algorithm 28.1
def qr_pure(A, maxiter=1000, tol=1e-12):
    for k in range(maxiter):
        Q, R = qr(A)
        A = R @ Q
        sub_diag = np.diag(A,-1)
        if np.abs(sub_diag).max() < tol:
            print('Converged in iterations = ', k+1)
            return A
    print('Did not converge in maxiter =', maxiter)
    return A

m = 5
B = np.random.random((m,m))
A = 0.5 * (B + B.T)
evals,evecs = eig(A)
print('Eigenvalues = ', evals)

Eigenvalues =  [ 2.64291079  0.61003094 -0.55932349 -0.26836586  0.11833976]

D = qr_pure(A.copy())
print(np.array_str(D,suppress_small=True))

Converged in iterations =  317
[[ 2.64291079 -0.         -0.          0.          0.        ]
 [ 0.          0.61003094 -0.         -0.          0.        ]
 [-0.         -0.         -0.55932349  0.         -0.        ]
 [-0.          0.          0.         -0.26836586 -0.        ]
 [ 0.          0.         -0.         -0.          0.11833976]]

Eigenvalues appear on diagonal with decreasing magnitude.

Practical QR algorithm#

We perform an initial reduction to Hessenberg form and use shifted QR with Rayleigh-quotient shifts.

# Algorithm 28.2
def qr_with_shifts(A, maxiter=1000, tol=1e-12):
    m = A.shape[0]
    H = hessenberg(A)
    for k in range(maxiter):
        mu = H[-1,-1]
        Q, R = qr(H - mu*np.eye(m))
        H = R @ Q + mu * np.eye(m)
        sub_diag = np.diag(H,-1)
        if np.abs(sub_diag).max() < tol:
            print('Converged in iterations = ', k+1)
            return H
    print('Did not converge in maxiter =', maxiter)
    return H

TODO: Add deflation step.

D = qr_with_shifts(A.copy())
print(np.array_str(D,suppress_small=True))

Converged in iterations =  109
[[ 2.64291079  0.         -0.          0.         -0.        ]
 [ 0.         -0.55932349 -0.         -0.         -0.        ]
 [ 0.         -0.          0.61003094  0.          0.        ]
 [ 0.          0.          0.         -0.26836586 -0.        ]
 [ 0.          0.          0.          0.          0.11833976]]

The convergence is usually faster than pure QR. The eigenvalues may not be ordered in any particular way on the diagonal.

Contents