JOVANA
Explore Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

The QR Algorithm

The power method hands you one eigenvalue at a time, slowly. The QR algorithm is the workhorse that quietly delivers them all — by doing nothing more exotic than factoring a matrix into Q R and multiplying the pieces back in the wrong order, over and over, until the matrix confesses its eigenvalues along the diagonal.

From one eigenvalue to all of them

In the last guide the power method gave you the single largest eigenvalue by repeatedly multiplying a vector by A and watching it line up with the dominant eigenvector, while inverse iteration let you zoom in on any eigenvalue you could already guess. Both are precious, but both are retail: they hand over eigenvalues one at a time, and getting the next one means deflating the matrix and starting again. For a dense matrix where you want all the eigenvalues — which is the everyday request in vibration analysis, quantum chemistry, stability studies — that retail approach is clumsy. We want a wholesale method.

And we must not be tempted by the textbook definition. Eigenvalues are the roots of the characteristic polynomial det(A - lambda I) = 0, so why not expand that determinant and solve the polynomial? Because, as the first guide in this rung argued, that route is a trap: the roots of a polynomial can be wildly sensitive to its coefficients, so forming the characteristic polynomial throws away accuracy you can never get back. The honest path computes eigenvalues directly from the matrix, never from its characteristic polynomial. The QR algorithm is that path.

The unreasonable trick: factor, then swap

Here is the whole idea, and it looks almost too simple to work. Start with A_0 = A. Factor it as A_0 = Q_0 R_0. Now form the next matrix by multiplying the two factors back together in the opposite order: A_1 = R_0 Q_0. Factor A_1 = Q_1 R_1, swap again to get A_2 = R_1 Q_1, and keep going. That is it. Repeatedly: factor into Q R, then recombine as R Q.

the bare QR iteration:

    A_0 = A
    repeat for k = 0, 1, 2, ... :
        Q_k, R_k = qr(A_k)        # factor:  A_k = Q_k R_k
        A_{k+1} = R_k Q_k         # swap the order back

key fact:  A_{k+1} = R_k Q_k = Q_k^T (Q_k R_k) Q_k = Q_k^T A_k Q_k

so each step is an ORTHOGONAL SIMILARITY transform of A_k --
same eigenvalues at every step, just a friendlier basis.
as k grows, A_k drifts toward upper-triangular form,
and its DIAGONAL entries converge to the eigenvalues of A.
The two-line loop and the one identity that explains why it works: A_{k+1} = Q_k^T A_k Q_k is a similarity transform, so eigenvalues are preserved, while the sequence creeps toward triangular form whose diagonal reveals them.

Why does this dull little shuffle converge to the eigenvalues? Two facts lock together. First, because Q_k is orthogonal, A_{k+1} = Q_k^T A_k Q_k is a similarity transform — A_{k+1} has exactly the same eigenvalues as A_k, and therefore as the original A. So no step ever changes the answer; the iteration only changes the basis you are looking at the matrix in. Second, and more subtly, this particular shuffle is secretly running the power method on every direction at once: the orthogonal columns of the accumulated Q's are themselves converging to the eigenvectors, sorted by the size of their eigenvalues. The upshot is that A_k slowly flattens into upper-triangular form, and a triangular matrix wears its eigenvalues openly on its diagonal.

A tiny mental picture: imagine A is a symmetric 2-by-2 matrix tilted at some angle. Each QR step rotates your coordinate axes a little toward the matrix's natural axes — its eigenvectors. Keep rotating and the off-diagonal entries shrink toward zero; when they vanish, the diagonal holds the eigenvalues. The QR algorithm is, in spirit, a patient machine for finding the right rotation, one factorization at a time.

Making it fast enough to use: Hessenberg first

The bare loop above is beautiful and unusable. Each QR factorization of a dense n-by-n matrix costs O(n^3) work, and you might need dozens of iterations, so a naive run is O(n^4) or worse — far too expensive. The fix that turns QR from a curiosity into the workhorse is a one-time preparation step: before iterating, reduce A to upper Hessenberg form — upper-triangular plus one extra nonzero sub-diagonal below the main diagonal. This is done once, exactly (up to rounding), by a sequence of orthogonal reflections, and it is itself a similarity transform, so it changes nothing about the eigenvalues.

Why does this one preparation pay off so enormously? A near-magical structural fact: the QR step preserves Hessenberg form. If A_k is Hessenberg, so is A_{k+1}, for free. And factoring a Hessenberg matrix is cheap — only that single sub-diagonal needs to be zeroed out, which a sweep of Givens rotations (or Householder reflections) does in O(n^2) work instead of O(n^3). One O(n^3) setup, then iterations that are each only O(n^2): that is the difference between a toy and a tool. For the symmetric case the payoff is even sharper, as the next section shows.

  1. Reduce A once to upper Hessenberg form by orthogonal similarity (O(n^3), done a single time, exact up to rounding) — eigenvalues unchanged.
  2. Apply the QR step to the Hessenberg matrix; because Hessenberg structure is preserved, each step costs only O(n^2).
  3. Watch a sub-diagonal entry shrink toward zero; when it is negligible, an eigenvalue has converged in the corner — split it off (deflate) and continue on the smaller block.
  4. Repeat until every sub-diagonal entry has vanished; the diagonal now holds all the eigenvalues.

Shifts and deflation: where the speed really comes from

Hessenberg form makes each step cheap, but plain QR can still need many steps because, like the power method underneath it, its convergence rate is governed by ratios of eigenvalues — if two eigenvalues are close in size, the relevant off-diagonal entry crawls toward zero. The cure is the same one that supercharged inverse iteration: a shift. Instead of factoring A_k, factor the shifted matrix A_k - mu_k I = Q_k R_k and then set A_{k+1} = R_k Q_k + mu_k I (the shift is added back so eigenvalues stay put). A well-chosen shift mu_k that sits near an eigenvalue makes the bottom-corner sub-diagonal entry plunge — typically the convergence becomes cubic for the symmetric case, meaning the number of correct digits roughly triples each step once you are close.

How do you pick a good shift without already knowing the eigenvalue? The trailing 2-by-2 block in the bottom-right corner is itself a tiny matrix whose eigenvalues you can compute by hand; using one of them (the Wilkinson shift) gives an excellent estimate of the eigenvalue about to converge. As soon as the corner sub-diagonal entry drops below a tolerance, that bottom-right diagonal entry is locked in as a converged eigenvalue. You then deflate: peel off that row and column and run the iteration on the smaller block that remains. This is the same shift-and-deflate rhythm you saw with the power method, now industrialized — converge a corner, lock it, shrink, repeat, until the matrix is fully triangular.

The symmetric case, honesty, and where this fits

When A is symmetric, everything gets better, which is why the symmetric eigenproblem deserves its own name. The Hessenberg reduction of a symmetric matrix produces a tridiagonal matrix — nonzeros only on the diagonal and the two neighbouring diagonals — because symmetry forces the entries above the first super-diagonal to vanish too. The QR iteration then operates on a tridiagonal form at O(n) per step, and with the Wilkinson shift it converges cubically. The eigenvalues are guaranteed real, the eigenvectors orthogonal, and the whole computation is backward stable. Computing all eigenvalues of a symmetric matrix this way costs about O(n^3) overall — genuinely affordable, and the routine you reach for thousands of times without thinking.

Now the honesty that this rung insists on. Every eigenvalue the QR algorithm returns is an approximation — the iteration stops when a sub-diagonal entry is below a threshold, not when it is exactly zero, because in floating point it essentially never reaches exact zero. The good news is that QR is backward stable: the eigenvalues you get are the exact eigenvalues of a matrix very close to A. But backward stability is only half the accuracy story. The other half is conditioning: if A's eigenvalues are themselves ill-conditioned — highly sensitive to perturbations of the matrix, as non-symmetric matrices with near-defective structure can be — then even a perfect backward-stable result can be far from the true eigenvalues. Accuracy equals conditioning times stability, here as everywhere. Symmetric matrices are the happy exception: their eigenvalues are perfectly conditioned, so backward stability buys you full accuracy.

Finally, the boundary of where QR belongs. It is the right tool for dense matrices up to a few thousand on a side, where you want all the eigenvalues — because it spends O(n^3) work and O(n^2) memory touching the whole matrix. For the giant sparse matrices of the iterative-solvers rung, where you only want a handful of extreme eigenvalues and could never afford to store a dense Hessenberg form, QR is the wrong shape, and the Arnoldi and Lanczos methods take over — Krylov-subspace cousins that, like the power method, lean only on the cheap product A x. Hold that division in mind: QR for dense-and-all, Krylov for sparse-and-few. With dense eigenvalues now in hand, the next guide turns to the singular value decomposition, which — fittingly — is computed by running a close relative of this very algorithm.