Sparse Matrices and Iterative Eigensolvers

Sparse matrices store only nonzero entries, enabling computation on systems with millions of unknowns. Iterative eigensolvers extract a few eigenpairs of such matrices using only matrix–vector products, without ever forming dense factorizations. Key result: Krylov subspace methods (Lanczos, Arnoldi) converge superlinearly to extremal eigenvalues, with error controlled by the spectral gap and Chebyshev polynomial growth. Applied throughout computational physics, structural mechanics, quantum chemistry, network science, and machine learning. Foundational contributions: Lanczos (1950), Arnoldi (1951), Saad, Sorensen, and the Davidson/LOBPCG lineage.

A sparse matrix is one in which the number of nonzero entries is $O (n)$ rather than $O (n^{2})$ , so that storing and operating on it scales linearly with dimension. Iterative eigensolvers compute selected eigenvalues and eigenvectors of large sparse (or merely structured) matrices by repeatedly applying the matrix to vectors and building low-dimensional approximations — typically over Krylov subspaces. They matter because direct eigensolvers (QR algorithm on a dense matrix) cost $O (n^{3})$ time and $O (n^{2})$ memory, which is utterly infeasible for the $n = 1 0^{6}$ – $1 0^{9}$ matrices arising from discretized PDEs, electronic structure, and graph Laplacians. The surprising fact: even though one needs all $n$ eigenvalues to "know" the spectrum, the few extremal eigenpairs that govern physics and stability can be obtained to high accuracy in far fewer than $n$ iterations, often with cost dominated by a handful of sparse matrix–vector products.

The Core Insight

You almost never need the matrix $A$ itself — you only need the action $v \mapsto A v$ . Every iterative eigensolver is built from this single primitive. The Krylov subspace $K_{m} (A, v) = span {v, A v, \dots, A^{m - 1} v}$ adaptively captures the directions along which $A$ amplifies most strongly, and the eigenvalues of the tiny projected matrix (the Rayleigh–Ritz values) converge rapidly to the extremal eigenvalues of $A$ . Sparsity makes the matvec cheap; Krylov projection makes few matvecs sufficient.

Formal Statement

Sparse Matrix and Storage Formats

A matrix $A \in R^{n \times n}$ (or $C^{n \times n}$ ) is called sparse if its number of nonzeros $nnz (A)$ satisfies $nnz (A) = O (n)$ , so that a matrix–vector product costs $O (n)$ work rather than $O (n^{2})$ . Common storage schemes:

COO (coordinate): triples $(i, j, a_{ij})$ .
CSR (compressed sparse row): a values array, a column-index array, and a row-pointer array rowptr of length $n + 1$ .
CSC (compressed sparse column): the column analogue.
Banded / DIA: store diagonals offset from the main diagonal.

The fundamental kernel is SpMV (sparse matrix–vector multiply):

(A v)_{i} = k = rowptr [i] \sum rowptr [i + 1] - 1 values [k] v (colind [k]) .

⟳Iterate over each row's stored nonzeros; the inner sum touches only stored entries, giving total work

\sum_{i} (rowptr [i + 1] - rowptr [i]) = nnz (A)

Eigenvalue Problem and Rayleigh Quotient

The standard eigenvalue problem seeks $(λ, x)$ , $x \neq = 0$ , with

A x = λ x, λ \in C, x \in C^{n} .

⟳

λ

is a root of

det (A - λ I) = 0

; eigenvectors span the null space of

A - λ I

For Hermitian $A = A^{*}$ , all eigenvalues are real and the Rayleigh quotient

ρ (x) = \frac{x ^{*} A x}{x ^{*} x}

satisfies

λ_{m i n} \leq ρ (x) \leq λ_{m a x}

, with stationary points exactly at eigenvectors.

⟳Diagonalize

A = Q Λ Q^{*}

and write

x = Q y

; then

ρ = \sum λ_{i} ∣ y_{i} ∣^{2} / \sum ∣ y_{i} ∣^{2}

, a convex combination of eigenvalues. Setting

\nabla ρ = 0

gives

A x = ρ (x) x

Krylov Subspace

Given $A$ and a starting vector $v_{1}$ with $∥ v_{1} ∥ = 1$ , the $m$ -th Krylov subspace is

K_{m} (A, v_{1}) = span {v_{1}, A v_{1}, A^{2} v_{1}, \dots, A^{m - 1} v_{1}} .

⟳It is the smallest

A

-invariant subspace containing

v_{1}

once

m

reaches the degree of the minimal polynomial of

v_{1}

w.r.t.

A

Polynomial characterization: $K_{m} = {p (A) v_{1} : de g p \leq m - 1}$ .

Arnoldi Factorization

For general $A$ , the Arnoldi process produces an orthonormal basis $V_{m} = [v_{1}, \dots, v_{m}]$ of $K_{m} (A, v_{1})$ and an upper Hessenberg matrix $H_{m}$ such that

A V_{m} = V_{m} H_{m} + h_{m + 1, m} v_{m + 1} e_{m}^{T},

with

H_{m} = V_{m}^{*} A V_{m}

. The Ritz values (eigenvalues of

H_{m}

) approximate eigenvalues of

A

⟳Apply Gram–Schmidt to

{v_{1}, A v_{1}, \dots}

; the orthogonalization coefficients populate

H_{m}

, and the new direction

v_{m + 1}

is what remains after subtracting projections of

A v_{m}

When $A$ is Hermitian, $H_{m}$ becomes tridiagonal $T_{m}$ and Arnoldi reduces to the Lanczos process:

A V_{m} = V_{m} T_{m} + β_{m} v_{m + 1} e_{m}^{T}, T_{m} = tridiag (β_{j - 1}, α_{j}, β_{j}) .

⟳Hermitian symmetry forces

H_{m} = H_{m}^{*}

; an upper Hessenberg Hermitian matrix is tridiagonal, so each

A v_{j}

orthogonalizes against only

v_{j - 1}

and

v_{j}

— the three-term recurrence.

Why It Holds / Derivation

Lanczos Three-Term Recurrence and Ritz Convergence

Step 1 — Recurrence. Starting from

v_{1}

β_{0} = 0

v_{0} = 0

, for

j = 1, 2, \dots

α_{j} = v_{j}^{*} A v_{j}, w_{j} = A v_{j} - α_{j} v_{j} - β_{j - 1} v_{j - 1}, β_{j} = ∥ w_{j} ∥, v_{j + 1} = w_{j} / β_{j} .

⟳This is Gram–Schmidt of

A v_{j}

against

{v_{j - 1}, v_{j}}

; symmetry guarantees orthogonality against all earlier

v_{i}

in exact arithmetic.

Step 2 — Projection. Collecting the recurrences gives [6] with

T_{m}

symmetric tridiagonal. Eigenpairs

(θ_{i}, s_{i})

T_{m}

yield Ritz pairs

(θ_{i}, y_{i})

with

y_{i} = V_{m} s_{i}

. Step 3 — Residual bound. The residual of a Ritz pair is computable for free:

∥ A y_{i} - θ_{i} y_{i} ∥ = ∣ β_{m} ∣ ∣ e_{m}^{T} s_{i} ∣.

⟳Multiply [6] on the right by

s_{i}

A V_{m} s_{i} = V_{m} T_{m} s_{i} + β_{m} v_{m + 1} e_{m}^{T} s_{i} = θ_{i} y_{i} + β_{m} v_{m + 1} (e_{m}^{T} s_{i})

; the last term has norm

∣ β_{m} ∣∣ e_{m}^{T} s_{i} ∣

since

v_{m + 1}

is a unit vector.

Step 4 — Convergence rate (Kaniel–Paige–Saad). For the largest eigenvalue

λ_{1} > λ_{2} \geq \dots \geq λ_{n}

0 \leq λ_{1} - θ_{1}^{(m)} \leq (λ_{1} - λ_{n}) \frac{tan ^{2} ∠ ( v _{1} , x _{1} )}{C _{m - 1} ( 1 + 2 γ ) ^{2}}, γ = \frac{λ _{1} - λ _{2}}{λ _{2} - λ _{n}},

where

C_{m - 1}

is the degree-

(m - 1)

Chebyshev polynomial of the first kind.

⟳The optimal Krylov polynomial that suppresses

{λ_{2}, \dots, λ_{n}}

while preserving

λ_{1}

is the shifted/scaled Chebyshev polynomial on

[λ_{n}, λ_{2}]

; its exponential growth outside that interval,

C_{m - 1} (1 + 2 γ) \sim \frac{1}{2} e^{2 (m - 1) γ}

, drives geometric convergence governed by the spectral gap.

Cornelius Lanczos | 1893–1974

A Hungarian-American mathematician and physicist who worked as an assistant to Einstein, Lanczos introduced his "method of minimized iterations" in 1950 while at the National Bureau of Standards, seeking eigenvalues without forming the full characteristic polynomial. He was deeply suspicious of catastrophic cancellation — fittingly, since loss of orthogonality became the central practical difficulty of his own algorithm. His broader oeuvre spans the Lanczos tensor in general relativity, the Lanczos sigma factor in Fourier analysis, and the Tau method for differential equations. Walter Arnoldi published the nonsymmetric generalization in 1951 for vibration analysis of aircraft.

Building Intuition

Think of $A$ as a transformation that stretches space anisotropically. Power iteration — repeatedly applying $A$ and renormalizing — converges to the dominant eigenvector because the largest $∣ λ ∣$ direction outgrows all others, but it throws away every intermediate vector. The Krylov subspace keeps those intermediate vectors. Within the $m$ -dimensional span ${v, A v, \dots, A^{m - 1} v}$ , you can form any polynomial $p (A) v$ of degree $< m$ . The eigensolver implicitly searches for the polynomial that maximally amplifies the eigenvectors you want and damps the rest. Because polynomials can grow explosively outside a fixed interval (Chebyshev), even a modest degree separates a well-isolated extremal eigenvalue from the bulk. The spectral gap $γ$ is the lever: a large gap means a low-degree polynomial suffices, while a tiny gap (clustered eigenvalues) forces many iterations or a preconditioning/shift-invert strategy.

Krylov convergence: Chebyshev amplification of the extremal eigenvalue

Worked Examples

Lanczos on a 4×4 tridiagonal matrix

Given. The 1D discrete Laplacian (with sign convention

A = - Δ_{h}

A = 2 - 1 00 - 1 2 - 1 0 0 - 1 2 - 1 00 - 1 2,

whose exact eigenvalues are

λ_{k} = 2 - 2 cos \frac{k π}{5} = {0.382, 1.382, 2.618, 3.618}

. Start vector

v_{1} = \frac{1}{2} (1, 1, 1, 1)^{T}

. Steps.

$A v_{1} = \frac{1}{2} (1, 0, 0, 1)^{T}$ , so $α_{1} = v_{1}^{*} A v_{1} = \frac{1}{4} (1 + 0 + 0 + 1) = 0.5$ .
$w_{1} = A v_{1} - α_{1} v_{1} = \frac{1}{2} (1, 0, 0, 1) - \frac{1}{4} (1, 1, 1, 1) = (0.25, - 0.25, - 0.25, 0.25)$ . Then $β_{1} = ∥ w_{1} ∥ = 0.5$ , $v_{2} = (0.5, - 0.5, - 0.5, 0.5)^{T}$ .
$A v_{2} = (1.5, - 1, - 1, 1.5)^{T} \cdot$ ... computing: $A v_{2} = (1.5, - 1.0, - 1.0, 1.5)^{T}$ ; $α_{2} = v_{2}^{*} A v_{2} = 0.5 (1.5) + (- 0.5) (- 1) + (- 0.5) (- 1) + 0.5 (1.5) = 2.5$ .
$T_{2} = (0.5 0.5 0.5 2.5)$ , with Ritz values $θ = \frac{3 \pm 4 + 1}{2} = {0.382, 2.618}$ .

Result. After just 2 Lanczos steps the Ritz values

θ_{1} = 0.382, θ_{2} = 2.618

already reproduce two exact eigenvalues of

A

(the symmetric starting vector aligns with two eigenmodes). Convergence to the full spectrum requires the remaining steps.

Cost comparison: dense QR vs. Lanczos for a PDE matrix

Given. A 3D Poisson problem discretized on a

200 \times 200 \times 200

grid gives

n = 8 \times 1 0^{6}

unknowns; the 7-point stencil yields

nnz (A) \approx 7 n = 5.6 \times 1 0^{7}

. We want the 10 smallest eigenvalues (lowest vibrational/diffusion modes). Steps.

Dense QR: memory $= n^{2} \times 8$ bytes $= (8 \times 1 0^{6})^{2} \times 8 \approx 5 \times 1 0^{14}$ bytes $= 512$ TB. Time $\sim 10 n^{3} \approx 5 \times 1 0^{21}$ flops. Infeasible.
Lanczos (shift-invert for smallest): one SpMV costs $\approx 2 nnz = 1.1 \times 1 0^{8}$ flops. Suppose $m = 200$ iterations to resolve 10 eigenpairs.
Total matvec work $\approx 200 \times 1.1 \times 1 0^{8} = 2.2 \times 1 0^{10}$ flops; orthogonalization $\approx 2 n m^{2} = 2 (8 \times 1 0^{6}) (200)^{2} = 6.4 \times 1 0^{11}$ flops. Memory: basis $V_{m}$ is $n \times m \times 8 = 8 \times 1 0^{6} \times 200 \times 8 = 1.3 \times 1 0^{10}$ bytes $= 13$ GB.

Result. Lanczos finishes in

\sim 7 \times 1 0^{11} flops, \sim 13 GB

versus

5 \times 1 0^{21}

flops and 512 TB for dense QR — roughly ten orders of magnitude less work and four orders less memory, the difference between a laptop-scale and an impossible computation.

Common Misconceptions

Loss of orthogonality is the silent killer

In exact arithmetic the Lanczos vectors ${v_{j}}$ are orthonormal and the three-term recurrence suffices. In floating point, once a Ritz value converges, the corresponding eigencomponent is reintroduced into later $v_{j}$ , destroying global orthogonality. The visible symptom: spurious "ghost" eigenvalues — copies of already-converged eigenvalues that appear and proliferate in $T_{m}$ . The cure is reorthogonalization: full (expensive, $O (n m^{2})$ ), selective (Parlett–Scott: reorthogonalize only when a Ritz pair converges), or partial (Simon). Never trust raw unreorthogonalized Lanczos multiplicities without the Cullum–Willoughby test.

Power iteration fails for the wrong eigenvalue distribution

Suppose you want the smallest eigenvalue of a positive-definite $A$ with spectrum ${1, 1.001, \dots, 1000}$ . Power iteration on $A$ converges to $λ_{m a x} = 1000$ , not what you want. Worse, even targeting the smallest via power iteration on $A^{- 1}$ converges at rate $∣ λ_{m i n} / λ_{2} ∣^{k} = (1/1.001)^{k}$ — agonizingly slow because the two smallest eigenvalues are clustered. The fix is shift-and-invert: apply $(A - σ I)^{- 1}$ with $σ$ near the target, which maps nearby eigenvalues to huge well-separated values $1/ (λ - σ)$ , restoring rapid convergence. This requires solving sparse linear systems each iteration — trading matvecs for factorizations.

Nonsymmetric eigenvalues can be wildly ill-conditioned

For non-normal $A$ (where $A A^{*} \neq = A^{*} A$ ), eigenvalue condition numbers $κ (λ_{i}) = 1/∣ y_{i}^{*} x_{i} ∣$ (with $y_{i}, x_{i}$ left/right eigenvectors) can be astronomically large. The Arnoldi Ritz values may then be meaningless even with a tiny residual. Pseudospectra $Λ_{ϵ} (A) = {z : ∥ (z I - A)^{- 1} ∥ \geq ϵ^{- 1}}$ are the honest diagnostic: they reveal how far eigenvalues move under perturbations of size $ϵ$ . A small residual bounds the backward error, not the forward eigenvalue error.

Connections

Iterative eigensolvers share their entire algorithmic core with iterative linear solvers. The Conjugate Gradient Method is Lanczos applied to solving $A x = b$ — same three-term recurrence, same Krylov subspace, different extraction. GMRES is Arnoldi for nonsymmetric linear systems. The Power Iteration is the degenerate one-vector ancestor of all Krylov methods, and Inverse Iteration plus Rayleigh Quotient Iteration (cubically convergent for symmetric problems) provide the refinement step. The Singular Value Decomposition of sparse $A$ is computed by Lanczos bidiagonalization on $A$ (equivalent to Lanczos on $A^{*} A$ ). Preconditioning — central to fast linear solvers — reappears as spectral transformation in shift-invert eigensolvers and in the LOBPCG and Jacobi-Davidson Method families. The convergence theory rests entirely on Chebyshev Polynomials and the Min-Max Theorem (Courant–Fischer).

Real-World Applications

Quantum chemistry / electronic structure. Density Functional Theory and Configuration Interaction require the lowest few eigenpairs of Hamiltonians with dimension up to $1 0^{10}$ ; the Davidson method (a preconditioned Krylov variant) was invented (1975) precisely for this and remains standard in codes like NWChem and Gaussian.
Structural and modal analysis. Finite-element models of bridges, aircraft, and engines compute natural frequencies via the generalized problem $K x = λ M x$ ; block Lanczos with shift-invert (as in Nastran) finds modes in a target frequency band.
Google PageRank and network centrality. The dominant eigenvector of the (sparse, $1 0^{9}$ -node) web hyperlink matrix is found by power iteration; eigenvector and Katz centrality use the same machinery on social and biological graphs.
Spectral graph partitioning and clustering. The Fiedler vector — eigenvector of the second-smallest eigenvalue of the graph Laplacian — partitions meshes and clusters data; computed by Lanczos or LOBPCG (used in METIS, scikit-learn spectral clustering).
Vibration / stability in fluid dynamics. Linear stability analysis of the Navier–Stokes operator uses Arnoldi (via ARPACK / IRAM) to find the rightmost eigenvalues governing the onset of turbulence.
Principal Component Analysis at scale. Truncated SVD of huge sparse data matrices (recommender systems, latent semantic analysis) uses Lanczos bidiagonalization (svds, irlba).
Lattice QCD. Eigenvalues of the sparse Dirac operator on a spacetime lattice probe chiral symmetry breaking and accelerate fermion inversions via deflation.

Historical Development

1929 — von Mises & Pollaczek-Geiringer formalize the power method for the dominant eigenvalue.
1950 — Cornelius Lanczos publishes the "method of minimized iterations" (NBS), introducing the symmetric three-term recurrence.
1951 — Walter Arnoldi generalizes to nonsymmetric matrices via the Hessenberg/Arnoldi factorization, motivated by aircraft flutter.
1952 — Hestenes & Stiefel introduce Conjugate Gradients, the linear-solver twin of Lanczos.
1966–1971 — Kaniel, Paige, and Saad establish sharp Chebyshev-based convergence bounds; Paige's 1971 thesis explains loss of orthogonality and rehabilitates Lanczos for finite-precision computation.
1975 — Ernest Davidson introduces the preconditioned Davidson method for quantum-chemistry Hamiltonians.
1979–1985 — Cullum & Willoughby develop ghost-eigenvalue identification; Parlett & Scott introduce selective reorthogonalization.
1992 — Sorensen develops the Implicitly Restarted Arnoldi Method (IRAM), basis of the celebrated ARPACK library (Lehoucq, Sorensen, Yang, 1998).
2001 — Knyazev publishes LOBPCG, an optimal preconditioned method for symmetric eigenproblems; Sleijpen & Van der Vorst's Jacobi–Davidson (1996) matures in parallel.

Subtleties & Limits

Several challenges remain at the research frontier. (1) Clustered and interior eigenvalues: Krylov methods excel at well-separated extremal eigenvalues but struggle when many eigenvalues cluster or lie in the spectrum's interior; contour-integral / spectral-projection methods (FEAST, Sakurai–Sugiura) compute all eigenvalues in a complex region but require many sparse linear solves and robust quadrature. (2) Non-normal and pseudospectral regimes: for highly non-normal operators, no cheap a posteriori bound certifies forward eigenvalue accuracy; computing pseudospectra at scale is itself expensive. (3) Communication-avoiding Krylov methods ( $s$ -step Lanczos/Arnoldi) reduce synchronization on exascale machines but exhibit worse numerical stability — balancing the two is open. (4) Nonlinear eigenvalue problems $T (λ) x = 0$ (arising in delay equations, photonics, fluid–structure interaction) lack a fully mature, robust large-scale iterative theory. (5) Eigensolvers for tensors and matrix functions ( $f (A) v$ , Lyapunov/Sylvester equations) are an active, only partially understood extension of the Krylov framework.