Prev
Unit 6 · Eigenvalue Problems
The QR Algorithm
60m100 XP
Next
🌀Eigenvalue Problems

The QR Algorithm

Similarity transforms, QR iteration, Hessenberg reduction, shifts and deflation, eig versus eigh

London, 1959
The National Research Development Corporation is a small British government body set up to turn research into products. Its staff include Christopher Strachey, one of Britain's computing pioneers, and his assistant John Francis — 25 years old, who left Cambridge without finishing a degree. Francis is studying a method by Heinz Rutishauser for finding eigenvalues by repeatedly factoring a matrix into triangular pieces. It is elegant and unstable. Francis replaces the triangular factorisation with an orthogonal one, adds a clever way of choosing shifts, and in two papers in The Computer Journal in 1961–62 describes what he calls the QR transformation. In the Soviet Union, Vera Kublanovskaya discovers the same idea independently in 1961. Then Francis leaves numerical analysis for good. Over the following decades his algorithm becomes the way every computer computes eigenvalues; in 2000 it is named one of the ten algorithms with the greatest influence on science and engineering in the twentieth century. In 2007 Gene Golub and Frank Uhlig tracked him down, retired in Hove on the south coast of England. He had had no idea.
Two matrices AA and B=S−1ASB = S^{-1}AS are similar: they represent the same linear map in different coordinates and have the same eigenvalues. Every practical eigenvalue algorithm applies similarity transforms until the eigenvalues can simply be read off — from the diagonal of a triangular matrix (or a diagonal one, if AA is symmetric). For stability the transforms should be orthogonal, S=QS = Q, so that rounding errors are never amplified (lesson 5-2).
💡The QR iteration
Starting from A0=AA_0 = A: Ak=QkRk(factorise),Ak+1=RkQk(multiply in reverse order).A_k = Q_kR_k \quad\text{(factorise)}, \qquad A_{k+1} = R_kQ_k \quad\text{(multiply in reverse order)}. Since Rk=QkTAkR_k = Q_k^{\mathsf T}A_k, we have Ak+1=QkTAkQkA_{k+1} = Q_k^{\mathsf T}A_kQ_k: a similarity transform, so the eigenvalues never change. Remarkably, under mild conditions AkA_k converges to upper-triangular form with the eigenvalues on the diagonal in decreasing order of magnitude. The subdiagonal entry (i+1,i)(i+1, i) shrinks like ∣λi+1/λi∣k|\lambda_{i+1}/\lambda_i|^k — the QR iteration is secretly the power method applied to a whole set of vectors at once, kept orthogonal by the factorisation.
Derivation
Why QR iteration is the power method in disguise
Collect the orthogonal factors: Q^k=Q0Q1⋯Qk−1\hat Q_k = Q_0Q_1\cdots Q_{k-1} and R^k=Rk−1⋯R1R0\hat R_k = R_{k-1}\cdots R_1R_0. Two facts follow by induction.
The iterates are AA seen in rotating coordinates. Each step is Aj+1=QjTAjQjA_{j+1} = Q_j^{\mathsf T}A_jQ_j, so after kk steps Ak=Q^kTAQ^kA_k = \hat Q_k^{\mathsf T}A\hat Q_k, that is AQ^k=Q^kAkA\hat Q_k = \hat Q_kA_k.
The accumulated factors are a QR factorisation of AkA^k. It holds for k=1k = 1. If Ak=Q^kR^kA^k = \hat Q_k\hat R_k, then
Ak+1=AQ^kR^k=Q^kAkR^k=Q^k(QkRk)R^k=Q^k+1R^k+1.A^{k+1} = A\hat Q_k\hat R_k = \hat Q_kA_k\hat R_k = \hat Q_k(Q_kR_k)\hat R_k = \hat Q_{k+1}\hat R_{k+1}.
Read the first column of Ak=Q^kR^kA^k = \hat Q_k\hat R_k: since R^k\hat R_k is upper triangular, Ake1A^k\mathbf e_1 is a multiple of the first column of Q^k\hat Q_k. That column is exactly the power method started from e1\mathbf e_1, so it turns towards the dominant eigenvector, and the first column of Ak=Q^kTAQ^kA_k = \hat Q_k^{\mathsf T}A\hat Q_k approaches (λ1,0,…,0)(\lambda_1, 0, \dots, 0). The first two columns of Q^k\hat Q_k span AkA^k applied to span⁡{e1,e2}\operatorname{span}\{\mathbf e_1, \mathbf e_2\} — the power method on a plane, which settles on the plane of the two dominant eigenvectors — and so on down the matrix. The subdiagonal entry (i+1,i)(i+1, i) measures how far the ii-th subspace still is from settling, which happens at the rate ∣λi+1/λi∣k|\lambda_{i+1}/\lambda_i|^k.
📐Watching it converge
For a symmetric tridiagonal 4×44 \times 4 matrix with eigenvalues 4.7453, 3.1773, 1.8227, 0.25474.7453,\ 3.1773,\ 1.8227,\ 0.2547, the off-diagonal part shrinks from 0.960.96 after one step to 2×10−22 \times 10^{-2} after 10, 4×10−44 \times 10^{-4} after 20 and 10−710^{-7} after 40. The smallest eigenvalue appears first — it is already correct to eight digits after five steps, because its ratio 0.2547/1.8227=0.140.2547/1.8227 = 0.14 is small — while the top pair, with ratio 0.670.67, takes longest.
Plain QR iteration is too slow and too expensive — O(n3)O(n^3) per step, many steps. Three refinements make it the algorithm inside every library:
  1. Hessenberg reduction. First transform AA by Householder reflections (similarity, O(n3)O(n^3) once) to upper Hessenberg form — zeros below the first subdiagonal. QR iteration preserves this form, and each step now costs only O(n2)O(n^2). For symmetric AA the Hessenberg form is tridiagonal, and a step costs O(n)O(n).
  2. Shifts. Iterate on Ak−μkI=QkRkA_k - \mu_kI = Q_kR_k, Ak+1=RkQk+μkIA_{k+1} = R_kQ_k + \mu_kI with μk\mu_k an estimate of the bottom eigenvalue (the Wilkinson shift, from the bottom 2×22\times2 block). The last subdiagonal entry then converges quadratically, cubically for symmetric matrices — the QR algorithm has absorbed Rayleigh quotient iteration.
  3. Deflation. As soon as the last subdiagonal entry is negligible, its diagonal entry is an eigenvalue: cut it off and continue on the smaller matrix.
On a random symmetric 8×88 \times 8 matrix, shifted QR with deflation finds all eight eigenvalues in 16 iterations — two per eigenvalue. Francis's final touch, the implicit double shift, handles complex-conjugate pairs of real nonsymmetric matrices without complex arithmetic.
Derivation
A shifted step is still a similarity transform
From Ak−μkI=QkRkA_k - \mu_kI = Q_kR_k we get Rk=QkT(Ak−μkI)R_k = Q_k^{\mathsf T}(A_k - \mu_kI). Substitute into the next iterate:
Ak+1=RkQk+μkI=QkT(Ak−μkI)Qk+μkI=QkTAkQk−μkQkTQk+μkI=QkTAkQk.A_{k+1} = R_kQ_k + \mu_kI = Q_k^{\mathsf T}(A_k - \mu_kI)Q_k + \mu_kI = Q_k^{\mathsf T}A_kQ_k - \mu_kQ_k^{\mathsf T}Q_k + \mu_kI = Q_k^{\mathsf T}A_kQ_k.
Subtracting the shift before factorising and adding it back afterwards changes nothing about the eigenvalues — only how fast they appear. Shifting by μk\mu_k makes the bottom eigenvalue of Ak−μkIA_k - \mu_kI tiny, so the ratio ∣λn−μk∣/∣λn−1−μk∣|\lambda_n - \mu_k|/|\lambda_{n-1} - \mu_k| that governs the last subdiagonal entry is tiny too.
Simplicity is prerequisite for reliability.
– Edsger W. Dijkstra, How do we tell truths that might hurt? (1975)
⚠eig or eigh?
numpy.linalg.eigh (and scipy.linalg.eigh) is for symmetric (Hermitian) matrices: it reduces to tridiagonal form, returns real eigenvalues in ascending order and orthonormal eigenvectors, and is several times faster and more accurate. numpy.linalg.eig is for general matrices: eigenvalues may be complex, eigenvectors need not be orthogonal, and they come in no particular order. A state matrix [0I−M−1K−M−1C]\begin{bmatrix} 0 & I \\ -M^{-1}K & -M^{-1}C \end{bmatrix} is not symmetric even when KK, MM and CC are — use eig, and read the natural frequency ∣λ∣|\lambda| and damping ratio −Re⁡λ/∣λ∣-\operatorname{Re}\lambda/|\lambda| off each complex pair. Calling eigh on a nonsymmetric matrix silently uses only one triangle and returns wrong answers.
Derivation
The state matrix, and reading frequency and damping from an eigenvalue
First-order form. A damped structure obeys Mx¨+Cx˙+Kx=0M\ddot{\mathbf x} + C\dot{\mathbf x} + K\mathbf x = \mathbf 0. Take positions and velocities together as the state, z=(x,x˙)\mathbf z = (\mathbf x, \dot{\mathbf x}). Its derivative is (x˙,x¨)(\dot{\mathbf x}, \ddot{\mathbf x}), and the equation of motion gives x¨=−M−1Kx−M−1Cx˙\ddot{\mathbf x} = -M^{-1}K\mathbf x - M^{-1}C\dot{\mathbf x}:
z˙=[0I−M−1K−M−1C]z.\dot{\mathbf z} = \begin{bmatrix} 0 & I \\ -M^{-1}K & -M^{-1}C \end{bmatrix}\mathbf z.
The identity block in the top right has no partner in the bottom left: the matrix is not symmetric.
One mode. For a single damped oscillator, x¨+2ζωnx˙+ωn2x=0\ddot x + 2\zeta\omega_n\dot x + \omega_n^2x = 0, the eigenvalues solve λ2+2ζωnλ+ωn2=0\lambda^2 + 2\zeta\omega_n\lambda + \omega_n^2 = 0:
λ=−ζωn±i ωn1−ζ2.\lambda = -\zeta\omega_n \pm i\,\omega_n\sqrt{1 - \zeta^2}.
Its size is ∣λ∣2=ζ2ωn2+ωn2(1−ζ2)=ωn2|\lambda|^2 = \zeta^2\omega_n^2 + \omega_n^2(1 - \zeta^2) = \omega_n^2, so ωn=∣λ∣\omega_n = |\lambda|, and the real part over the size gives the damping ratio, ζ=−Re⁡λ/∣λ∣\zeta = -\operatorname{Re}\lambda/|\lambda|. For λ=−1+28.27i\lambda = -1 + 28.27i: ωn=28.29\omega_n = 28.29 rad/s (4.504.50 Hz) and ζ=0.035\zeta = 0.035.
💡Where you will use this
Control Systems and Attitude Control: closed-loop poles of state-space models, computed by eig — QR underneath. Rocket Propulsion and Aircraft Engines: stability of linearised dynamics (pogo, combustion, rotor whirl). Structures: dense eigenproblems of reduced models after component mode synthesis. Every course: whenever a library returns eigenvalues of a dense matrix, Francis's algorithm computed them.
📚Further reading
  • Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lectures 25–29 — Hessenberg reduction and the QR algorithm without and with shifts.
  • David Watkins, 'The QR algorithm revisited', SIAM Review 50(1), 133–145 (2008) — why it works, explained through subspace iteration.
  • Gene Golub and Frank Uhlig, 'The QR algorithm: 50 years later — its genesis by John Francis and Vera Kublanovskaya and subsequent developments', IMA Journal of Numerical Analysis 29(3), 467–485 (2009) — the story of the cold open.
↓But there is a problem
We can now find every eigenvalue of any matrix. For a structure, though, the stiffness matrix is only half the story. A building does not vibrate according to KK alone — floors are heavy, columns are light, and the mass is distributed unevenly. The right equation is Kϕ=ω2MϕK\phi = \omega^2 M\phi, with a mass matrix on the right-hand side. On 10 June 2000, thousands of people walked across London's new Millennium Bridge on its opening day. Within minutes the deck was swaying sideways so much that people had to stop and hold the handrails. The bridge was closed two days later — for nearly two years.