The eigenvalue problem \(\bA\bv = \lambda\bv\) seeks the invariant directions of a linear operator. Along an eigenvector \(\bv\), the action of the matrix \(\bA\) reduces to scalar multiplication by \(\lambda\). Eigenvalues govern the long-term stability of dynamical systems, the resonant frequencies of vibrating structures, the convergence rates of iterative solvers, and the low-dimensional structure of data in principal component analysis.
In introductory algebra, eigenvalues are introduced as the roots of the characteristic polynomial \(p(\lambda) = \det(\bA - \lambda \bI) = 0\). In numerical computing, however, computing eigenvalues via characteristic polynomials is strictly avoided. By Abel’s Impossibility Theorem, no closed-form radical formula exists for \(n \geq 5\), meaning all eigenvalue algorithms must be iterative. Furthermore, polynomial root-finding is notoriously ill-conditioned, whereas matrix eigenvalue algorithms can be executed with backward stability.
Power Iteration and Shifted Inverse Iteration
When only a single dominant eigenvalue or an eigenvalue closest to a target frequency is required, vector iterations provide rapid, memory-efficient solutions.
Given an initial vector \(\bv^{(0)}\) with \(\|\bv^{(0)}\|_2 = 1\), power iteration generates successive approximations via \[
\begin{align}
\bw^{(k)} = \bA \bv^{(k-1)}, \qquad \bv^{(k)} = \frac{\bw^{(k)}}{\|\bw^{(k)}\|_2}, \qquad \lambda^{(k)} = (\bv^{(k)})^T \bA \bv^{(k)}.
\end{align}
\] If the dominant eigenvalue satisfies \(|\lambda_1| > |\lambda_2| \geq ... \geq |\lambda_n|\) and the initial vector has a nonzero component along the dominant eigenvector \(\bq_1\), then \(\bv^{(k)}\) converges to \(\bq_1\) at the linear rate \[
\begin{align}
\|\bv^{(k)} - (\pm \bq_1)\|_2 = O\left( \left| \frac{\lambda_2}{\lambda_1} \right|^k \right).
\end{align}
\]
To compute the eigenvalue closest to a target scalar \(\sigma \in \fC\), power iteration is applied to the shifted and inverted operator \((\bA - \sigma \bI)^{-1}\): \[
\begin{align}
(\bA - \sigma \bI) \bw^{(k)} = \bv^{(k-1)}, \qquad \bv^{(k)} = \frac{\bw^{(k)}}{\|\bw^{(k)}\|_2}.
\end{align}
\] The eigenvalues of \((\bA - \sigma \bI)^{-1}\) are \((\lambda_i - \sigma)^{-1}\). The eigenvalue of \(\bA\) closest to \(\sigma\) becomes the dominant eigenvalue of the shifted operator, accelerating convergence to the linear rate \(|\lambda_J - \sigma| / |\lambda_K - \sigma|\), where \(\lambda_J\) is the closest eigenvalue and \(\lambda_K\) is the second closest.
For symmetric matrices, updating the shift dynamically to the current Rayleigh quotient \(\sigma_k = (\bv^{(k)})^T \bA \bv^{(k)}\) yields Rayleigh Quotient Iteration. At each step, we solve \[
\begin{align}
(\bA - \sigma_k \bI) \bw^{(k+1)} = \bv^{(k)}, \qquad \bv^{(k+1)} = \frac{\bw^{(k+1)}}{\|\bw^{(k+1)}\|_2}, \qquad \sigma_{k+1} = (\bv^{(k+1)})^T \bA \bv^{(k+1)}.
\end{align}
\] For symmetric matrices, RQI converges cubically: the number of correct decimal digits triples at each iteration once the iterate enters the asymptotic regime.
The Practical QR Algorithm
The standard method for computing the complete spectrum of a dense matrix is the QR algorithm with shifts.
The unshifted QR iteration computes a sequence of orthogonally similar matrices:
Initialize \(\bA_0 = \bA\).
For \(k = 0, 1, 2, ...\): \begin{enumerate}
Compute the QR factorization: \(\bA_k = \bQ_k \bR_k\).
Recombine the factors in reverse order: \(\bA_{k+1} = \bR_k \bQ_k = \bQ_k^T \bA_k \bQ_k\).
\end{enumerate} Because \(\bA_{k+1}\) is orthogonally similar to \(\bA_k\), all iterates share identical eigenvalues. As \(k \to \infty\), \(\bA_k\) converges to the upper triangular Schur form \(\bT\).
To make the QR algorithm computationally efficient, standard numerical libraries implement two essential enhancements:
Preliminary Hessenberg Reduction: The matrix \(\bA\) is pre-transformed to upper Hessenberg form \(\bH_0 = \bQ_0^T \bA \bQ_0\) using Householder reflectors in \(O(n^3)\) operations once. Subsequent QR steps preserve the Hessenberg structure and require only \(O(n^2)\) flops using Givens rotations.
Wilkinson Shifts: At each step, a shift \(\mu_k\) is chosen as the eigenvalue of the trailing \(2 \times 2\) submatrix of \(\bH_k\) closest to \(h_{nn}\): \[
\begin{align}
\bH_k - \mu_k \bI = \bQ_k \bR_k, \qquad \bH_{k+1} = \bR_k \bQ_k + \mu_k \bI.
\end{align}
\] The shift accelerates subdiagonal convergence to cubic order for symmetric matrices and quadratic order for nonsymmetric matrices. When the subdiagonal entry satisfies \(|h_{n,n-1}| \leq \varepsilon_{\text{mach}}(|h_{nn}| + |h_{n-1,n-1}|)\), \(h_{nn}\) is decoupled as a converged eigenvalue and the matrix is deflated.
Compute the dominant eigenvalue and eigenvector of \(\bA = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}\) by hand using two steps of power iteration starting with \(\bv^{(0)} = (1, 0)^T\).
Apply one step of unshifted QR iteration to \(\bA = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\). Compute the QR factors explicitly and show that the off-diagonal entries decrease.
Let \(\bA = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}\). Show that unshifted QR iteration fails to converge because the factors oscillate indefinitely. Explain how introducing a shift resolves this cycle.