QR algorithm

The QR algorithm is the standard method for finding all eigenvalues of a dense matrix, and it is almost magically simple to state. Factor A_0 = A into Q_0 R_0 (orthogonal times upper-triangular), then multiply the factors back in the reverse order to get A_1 = R_0 Q_0. Repeat: factor A_1 = Q_1 R_1, set A_2 = R_1 Q_1, and so on. Astonishingly, the sequence A_0, A_1, A_2, ... converges to (block) upper-triangular form, with the eigenvalues marching onto the diagonal.

Why it works: each step is a similarity transform, A_{k+1} = R_k Q_k = Q_k^T A_k Q_k, so every A_k has the same eigenvalues as A. The iteration is secretly running a sophisticated version of the power method on all eigenvectors at once; the orthogonal factors gradually rotate the matrix into the Schur form, where eigenvalues sit on the diagonal.

Naively this would cost O(n^3) per step and converge slowly — useless. Two refinements make it the practical standard. First, reduce A once to Hessenberg form (upper-triangular plus one subdiagonal); the QR step then preserves Hessenberg structure and costs only O(n^2). Second, use shifts: instead of factoring A_k, factor A_k - mu I for a cleverly chosen mu near an eigenvalue, which dramatically accelerates convergence — to cubic near the end. Modern implementations use the implicit double-shift (Francis) step to handle complex eigenvalue pairs in real arithmetic.

For symmetric matrices the algorithm specializes beautifully: Hessenberg form becomes tridiagonal, and the symmetric QR algorithm finds all eigenvalues and eigenvectors with backward stability and remarkable speed. This is what library routines like LAPACK's syev and geev call under the hood, and it is why you should never write your own dense eigensolver.