The Lanczos algorithm is the symmetric specialization of the Arnoldi iteration for large, sparse real symmetric matrices \(\bA = \bA^T \in \fR^{n \times n}\). When \(\bA\) is symmetric, the upper Hessenberg projection matrix becomes symmetric tridiagonal, collapsing the \(k\)-term Gram-Schmidt orthogonalization into a three-term recurrence involving only the two immediately preceding basis vectors.
By projecting the large \(n \times n\) matrix onto an expanding \(k\)-dimensional Krylov subspace (\(k \ll n\)), the Lanczos algorithm rapidly computes extremal eigenvalues (the largest and smallest frequencies in the spectrum) using only one sparse matrix-vector product per step and minimal vector memory.
The Three-Term Lanczos Recurrence
Let \(\bA \in \fR^{n \times n}\) be real symmetric, and let \(\bv_1 \in \fR^n\) be a normalized initial vector (\(\|\bv_1\|_2 = 1\)). Set \(\bv_0 = \bzero\) and \(\beta_0 = 0\). For \(j = 1, 2, ..., k\):
Compute the matrix-vector product: \(\bw = \bA\bv_j\).
Compute the diagonal Rayleigh quotient: \(\alpha_j = \bv_j^T \bw\).
Orthogonalize against the previous two vectors: \(\bw \leftarrow \bw - \alpha_j \bv_j - \beta_{j-1} \bv_{j-1}\).
Compute the subdiagonal norm: \(\beta_j = \|\bw\|_2\).
If \(\beta_j = 0\), terminate. Otherwise, set \(\bv_{j+1} = \bw / \beta_j\).
In matrix notation, after \(k\) steps, the orthonormal Lanczos basis \(\bV_k = [\bv_1, ..., \bv_k]\) satisfies the algebraic relation \[
\begin{align}
\bA \bV_k = \bV_k \bT_k + \beta_k \bv_{k+1} \be_k^T,
\end{align}
\] where \(\bT_k \in \fR^{k \times k}\) is the symmetric tridiagonal matrix \[
\begin{align}
\bT_k = \begin{pmatrix}
\alpha_1 & \beta_1 & 0 & ... & 0 \\
\beta_1 & \alpha_2 & \beta_2 & ... & 0 \\
0 & \beta_2 & \alpha_3 & \ddots & \vdots \\
\vdots & \vdots & \ddots & \ddots & \beta_{k-1} \\
0 & 0 & ... & \beta_{k-1} & \alpha_k
\end{pmatrix}.
\end{align}
\]
In exact arithmetic, the Lanczos basis vectors are mutually orthonormal (\(\bV_k^T\bV_k = \bI_k\)). Multiplying the Lanczos relation on the left by \(\bV_k^T\) yields the Rayleigh-Ritz Galerkin projection: \[
\begin{align}
\bT_k = \bV_k^T \bA \bV_k.
\end{align}
\]
Ritz Values, Ritz Vectors, and Error Bounds
Let \(\bT_k \bs_i = \theta_i \bs_i\) be the eigendecomposition of the small symmetric tridiagonal matrix \(\bT_k\), with normalized eigenvectors \(\|\bs_i\|_2 = 1\).
The eigenvalue \(\theta_i \in \fR\) is termed the \(i\)-th Ritz value of \(\bA\) with respect to \(\mathcal{K}_k(\bA, \bv_1)\).
The vector \(\tilde{\bu}_i = \bV_k \bs_i \in \fR^n\) is termed the \(i\)-th Ritz vector.
The residual of the Ritz pair \((\theta_i, \tilde{\bu}_i)\) satisfies the explicit identity \[
\begin{align}
\|\bA \tilde{\bu}_i - \theta_i \tilde{\bu}_i\|_2 = \beta_k |s_{k,i}|,
\end{align}
\] where \(s_{k,i}\) is the final component of the \(i\)-th tridiagonal eigenvector \(\bs_i \in \fR^k\).
Multiplying the Lanczos relation \(\bA \bV_k = \bV_k \bT_k + \beta_k \bv_{k+1} \be_k^T\) on the right by \(\bs_i\) gives \[
\begin{align}
\bA (\bV_k \bs_i) = \bV_k (\bT_k \bs_i) + \beta_k \bv_{k+1} (\be_k^T \bs_i)
\implies \bA \tilde{\bu}_i = \theta_i \tilde{\bu}_i + \beta_k s_{k,i} \bv_{k+1}.
\end{align}
\] Rearranging and taking the Euclidean norm yields \[
\begin{align}
\|\bA \tilde{\bu}_i - \theta_i \tilde{\bu}_i\|_2 = \|\beta_k s_{k,i} \bv_{k+1}\|_2 = \beta_k |s_{k,i}| \|\bv_{k+1}\|_2 = \beta_k |s_{k,i}|,
\end{align}
\] because \(\|\bv_{k+1}\|_2 = 1\). This allows monitoring the convergence of every Ritz pair for \(O(k)\) operations without forming the \(n\)-dimensional Ritz vectors \(\tilde{\bu}_i\).
Loss of Orthogonality and Ghost Eigenvalues
(Paige’s theorem on loss of orthogonality) In floating-point arithmetic, the three-term recurrence suffers from a fundamental numerical instability discovered by Chris Paige. The basis vectors \(\bv_j\) do not remain orthogonal; rather, orthogonality is lost rapidly precisely as individual Ritz values converge to extreme eigenvalues of \(\bA\). As a consequence, once an eigenvalue has converged, the algorithm continues to re-discover it in subsequent iterations, generating multiple spurious copies known as ghost eigenvalues.
(Practical remedies for ghost eigenvalues) To maintain numerical reliability on large problems, scientific codes employ one of three strategies:
Full Reorthogonalization: Reorthogonalize \(\bv_{j+1}\) against all previous basis vectors \(\bv_1, ..., \bv_j\) using Modified Gram-Schmidt at every step (\(O(k^2 n)\) cost).
Selective Reorthogonalization: Orthogonalize only against converged Ritz vectors whenever loss of orthogonality is detected (\(O(k \cdot n)\) cost).
Cullum-Willoughby Identification: Run the unorthogonalized Lanczos iteration and filter out ghost eigenvalues by comparing the spectrum of \(\bT_k\) with the spectrum of the submatrix formed by deleting the first row and column of \(\bT_k\).
Execute two steps of the Lanczos algorithm by hand for the diagonal matrix \(\bA = \operatorname{diag}(4, 2, 1)\) starting with the normalized vector \(\bv_1 = \frac{1}{\sqrt{3}}(1, 1, 1)^T\). Explicitly write down \(\alpha_1, \beta_1, \alpha_2\), and the \(2 \times 2\) tridiagonal matrix \(\bT_2\).
For the \(2 \times 2\) matrix \(\bT_2\) computed above, find its Ritz values and evaluate the residual error bound using the result above.
Prove that the Lanczos algorithm applied to the shifted system \((\bA - \sigma\bI)\) produces the identical basis vectors \(\bV_k\) as the original system, with tridiagonal matrix \(\bT_k - \sigma\bI\).