Lanczos Algorithm

Lanczos algorithm is an iterative projection method in numerical linear algebra and scientific computing, devised by Cornelius Lanczos in 1950, that approximates the extremal eigenvalues and eigenvectors of a large, sparse symmetric real matrix A ∈ ℝ^n×n (or Hermitian complex matrix A ∈ ℂ^n×n). Rather than applying dense orthogonal similarity transformations directly to A—which would destroy sparsity and require O(n³) arithmetic operations and O(n²) memory—the Lanczos method constructs an orthonormal basis Q_k = [q₁, q₂, …, q_k] for the order-k Krylov subspace K_k(A, q₁) = span{q₁, Aq₁, A²q₁, …, A^k−1q₁} using only matrix-vector multiplications v ↦ Av.

Because A is symmetric, the projected Rayleigh–Ritz matrix T_k = Q_k^TAQ_k is both upper Hessenberg and symmetric, forcing T_k to be a symmetric tridiagonal matrix. As a consequence, generating each new basis vector q_j+1 requires orthogonalizing Aq_j against only the two most recent vectors q_j and q_j−1 via a three-term recurrence relation. This close structural relationship links the Lanczos process to the conjugate gradient method for symmetric positive-definite linear systems and distinguishes it from the Arnoldi iteration used inside the generalized minimal residual method for non-symmetric matrices.

Starting with a unit-norm starting vector q₁ ∈ ℝⁿ (‖q₁‖₂ = 1), setting β₁ = 0 and q₀ = 0, step j (for j = 1, 2, …, k) of the Lanczos iteration executes four vector operations:

1. Matrix-vector product: w_j = Aq_j − β_jq_j−1
2. Diagonal Rayleigh quotient: α_j = q_j^Tw_j
3. Orthogonal residual: r_j = w_j − α_jq_j
4. Off-diagonal norm and normalization: β_j+1 = ‖r_j‖₂; if β_j+1 ≠ 0, set q_j+1 = r_j / β_j+1.

In compact matrix form, after k steps the Lanczos vectors satisfy the Lanczos decomposition:

AQ_k = Q_kT_k + β_k+1q_k+1e_k^T,

where e_k is the k-th canonical unit vector of ℝ^k, and T_k ∈ ℝ^k×k is the symmetric tridiagonal matrix with main diagonal entries α₁, …, α_k and sub-/superdiagonal entries β₂, …, β_k. Each iteration requires one sparse matrix-vector product (costing O(nnz(A)) flops) plus O(n) vector arithmetic and storage for only three vectors of length n.

Ritz Values, Ritz Vectors, and Kaniel–Paige Convergence Theory

Let T_k = S Θ S^T be the spectral decomposition of the small k × k tridiagonal matrix T_k, computed in O(k²) operations via the symmetric tridiagonal QR decomposition or divide-and-conquer eigensolver. The eigenvalues θ₁, …, θ_k on the diagonal of Θ are called the Ritz values, and the lifted vectors y_i = Q_ks_i (where s_i is the i-th column of S) are the Ritz vectors. Crucially, the residual norm of each Ritz pair (θ_i, y_i) can be evaluated without ever forming y_i ∈ ℝⁿ:

‖Ay_i − θ_iy_i‖₂ = ‖(AQ_k − Q_kT_k)s_i‖₂ = β_k+1 |e_k^Ts_i| = β_k+1 |s_k,i|,

meaning that convergence of the i-th Ritz pair is monitored directly from the bottom entry s_k,i of the eigenvector of T_k. By the Kaniel–Paige–Saad convergence bound, which uses Chebyshev polynomials of degree k − 1 on the interior spectrum [λ₂, λ_n], the largest eigenvalue λ₁ is approximated with error decaying as O(e^{−4(k−1)√γ}), where γ = (λ₁ − λ₂)/(λ₂ − λ_n) is the relative spectral gap. Extremal eigenvalues at the top and bottom of the spectrum therefore emerge in k ≪ n iterations (often k ≈ 50–200 even when n > 10⁷).

Although infinite-precision arithmetic ensures Q_k^TQ_k = I_k, Christopher C. Paige proved in his landmark 1971 dissertation that in finite-precision floating-point arithmetic, loss of orthogonality among the Lanczos vectors is not a gradual random drift; instead, it occurs precisely as the first Ritz values converge to eigenvalues of A. Once a Ritz vector y_i converges to accuracy ε_mach ‖A‖ / β_k+1, rounding errors re-inject components along y_i into subsequent vectors q_j, spawning spurious duplicate "ghost" copies of already converged eigenvalues in T_k.

Three main remedies restore numerical stability in production eigensolvers:

Quantum Many-Body Hamiltonian Ground States

In condensed-matter physics, quantum chemistry, and lattice gauge theory, the Hamiltonian matrix H of an interacting spin or fermion system on a Hilbert space of dimension n = 2^L (or 4^L) is far too large to store as a dense matrix, yet applying H to a state vector ψ requires only O(L n) operations. The Lanczos algorithm extracts the ground-state energy E₀ (lowest eigenvalue) and low-lying excitations in 50 to 300 matrix-vector products, serving as the core eigensolver inside Exact Diagonalization (ED) and the Density Matrix Renormalization Group (DMRG).

Golub–Kahan Lanczos Bidiagonalization for SVD and PCA

Applying the symmetric Lanczos process to the augmented matrix [0, A; A^T, 0] or equivalently to A^TA without explicitly forming the product yields Golub–Kahan Lanczos bidiagonalization, which constructs left and right orthonormal bases U_k and V_k reducing rectangular A ∈ ℝ^m×n to a k × k upper bidiagonal matrix B_k. This algorithm computes the top k singular values and vectors for large-scale principal component analysis, latent semantic indexing, and the LSQR iterative solver for ill-conditioned least-squares systems.

Shift-and-Invert Spectral Transformations

When interior eigenvalues near a target shift σ are needed—such as acoustic resonance modes in finite-element structural mechanics—running the Lanczos iteration on the shifted-and-inverted operator (A − σI)⁻¹ maps eigenvalues λ_i ≈ σ to huge extremal eigenvalues (λ_i − σ)⁻¹, where sparse LU decomposition or Cholesky decomposition factors A − σI once before the Lanczos loop.

How does the Lanczos algorithm relate to the conjugate gradient method?

When A is symmetric positive-definite, the conjugate gradient method for solving Ax = b computes the unique vector x_k ∈ K_k(A, r₀) that solves the projected tridiagonal system T_kz_k = ‖r₀‖₂ e₁ with x_k = Q_kz_k, updating x_k via an implicit LDL^T factorization of T_k.

Why does Lanczos fail to find the multiplicity of degenerate eigenvalues in a single run?

Because the cyclic Krylov subspace K_k(A, q₁) generated by a single starting vector q₁ can contain only one linearly independent eigenvector per distinct eigenvalue of A, single-vector Lanczos resolves only one copy of a degenerate eigenspace. Block Lanczos, which propagates a tall block Q₁ ∈ ℝ^n×b of b orthogonal vectors simultaneously, resolves eigenspaces of multiplicity up to b.

1. Lanczos, C. (1950). An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. Journal of Research of the National Bureau of Standards, 45(4), 255–282. https://doi.org/10.6028/jres.045.026
2. Paige, C. C. (1972). Computational variants of the Lanczos method for the eigenproblem. Journal of the Institute of Mathematics and Its Applications, 10(3), 373–381. https://doi.org/10.1093/imamat/10.3.373
3. Parlett, B. N., & Scott, D. S. (1979). The Lanczos algorithm with selective orthogonalization. Mathematics of Computation, 33(145), 217–238. https://doi.org/10.1090/S0025-5718-1979-0514820-3
4. Saad, Y. (2011). Numerical Methods for Large Eigenvalue Problems (Revised ed.). Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9781611970739