Sparse Matrices and Iterative Eigensolvers

⚡ Key Takeaways

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.

💡 The Core Insight

You almost never need the matrix $A$ itself — you only need the action $v \mapsto Av$ . Every iterative eigensolver is built from this single primitive. The Krylov subspace $\mathcal{K}_m(A,v) = \mathrm{span}\{v, Av, \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.

DefinitionSparse Matrix and Storage Formats

A matrix $A \in \mathbb{R}^{n\times n}$ (or $\mathbb{C}^{n\times n}$ ) is called sparse if its number of nonzeros $\mathrm{nnz}(A)$ satisfies $\mathrm{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):

1

\quad (Av)_i = \sum_{k=\text{rowptr}[i]}^{\text{rowptr}[i+1]-1} \text{values}[k]\, v\big(\text{colind}[k]\big).

Derivation hint: Iterate over each row's stored nonzeros; the inner sum touches only stored entries, giving total work $\sum_i (\text{rowptr}[i+1]-\text{rowptr}[i]) = \mathrm{nnz}(A)$ .

DefinitionEigenvalue Problem and Rayleigh Quotient

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

2

\quad A x = \lambda x, \qquad \lambda \in \mathbb{C},\ x \in \mathbb{C}^n.

Derivation hint: $\lambda$ is a root of $\det(A-\lambda I)=0$ ; eigenvectors span the null space of $A-\lambda I$ . For Hermitian $A=A^*$ , all eigenvalues are real and the Rayleigh quotient

3

\quad \rho(x) = \frac{x^* A x}{x^* x}

satisfies $\lambda_{\min} \le \rho(x) \le \lambda_{\max}$ , with stationary points exactly at eigenvectors. Derivation hint: Diagonalize $A=Q\Lambda Q^*$ and write $x=Qy$ ; then $\rho = \sum \lambda_i |y_i|^2 / \sum |y_i|^2$ , a convex combination of eigenvalues. Setting $\nabla\rho=0$ gives $Ax=\rho(x)x$ .

DefinitionKrylov Subspace

Given $A$ and a starting vector $v_1$ with $\|v_1\|=1$ , the $m$ -th Krylov subspace is

4

\quad \mathcal{K}_m(A,v_1) = \mathrm{span}\{v_1,\, Av_1,\, A^2 v_1,\, \dots,\, A^{m-1}v_1\}.

Derivation hint: 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: $\mathcal{K}_m = \{ p(A)v_1 : \deg p \le m-1 \}$ .

TheoremArnoldi Factorization

For general $A$ , the Arnoldi process produces an orthonormal basis $V_m=[v_1,\dots,v_m]$ of $\mathcal{K}_m(A,v_1)$ and an upper Hessenberg matrix $H_m$ such that

5

\quad A V_m = V_m H_m + h_{m+1,m}\, v_{m+1} e_m^{\mathsf T},

with $H_m = V_m^* A V_m$ . The Ritz values (eigenvalues of $H_m$ ) approximate eigenvalues of $A$ . Derivation hint: Apply Gram–Schmidt to $\{v_1, Av_1, \dots\}$ ; the orthogonalization coefficients populate $H_m$ , and the new direction $v_{m+1}$ is what remains after subtracting projections of $Av_m$ . When $A$ is Hermitian, $H_m$ becomes tridiagonal $T_m$ and Arnoldi reduces to the Lanczos process:

6

\quad A V_m = V_m T_m + \beta_m v_{m+1} e_m^{\mathsf T},\qquad T_m=\mathrm{tridiag}(\beta_{j-1},\alpha_j,\beta_j).

Derivation hint: Hermitian symmetry forces $H_m=H_m^*$ ; an upper Hessenberg Hermitian matrix is tridiagonal, so each $Av_j$ orthogonalizes against only $v_{j-1}$ and $v_j$ — the three-term recurrence.

ProofLanczos Three-Term Recurrence and Ritz Convergence

Step 1 — Recurrence. Starting from $v_1$ , $\beta_0=0$ , $v_0=0$ , for $j=1,2,\dots$ :

7

\quad \alpha_j = v_j^* A v_j,\quad w_j = A v_j - \alpha_j v_j - \beta_{j-1} v_{j-1},\quad \beta_j = \|w_j\|,\quad v_{j+1}=w_j/\beta_j.

Derivation hint: This is Gram–Schmidt of $Av_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 $(\theta_i, s_i)$ of $T_m$ yield Ritz pairs $(\theta_i, y_i)$ with $y_i = V_m s_i$ .

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

8

\quad \|A y_i - \theta_i y_i\| = |\beta_m|\,|e_m^{\mathsf T} s_i|.

Derivation hint: Multiply [6] on the right by $s_i$ : $AV_m s_i = V_m T_m s_i + \beta_m v_{m+1}e_m^{\mathsf T}s_i = \theta_i y_i + \beta_m v_{m+1}(e_m^{\mathsf T}s_i)$ ; the last term has norm $|\beta_m||e_m^{\mathsf T}s_i|$ since $v_{m+1}$ is a unit vector.

Step 4 — Convergence rate (Kaniel–Paige–Saad). For the largest eigenvalue $\lambda_1 > \lambda_2 \ge \cdots \ge \lambda_n$ ,

9

\quad 0 \le \lambda_1 - \theta_1^{(m)} \le (\lambda_1-\lambda_n)\,\frac{\tan^2\angle(v_1,x_1)}{C_{m-1}(1+2\gamma)^2},\qquad \gamma=\frac{\lambda_1-\lambda_2}{\lambda_2-\lambda_n},

where $C_{m-1}$ is the degree- $(m-1)$ Chebyshev polynomial of the first kind. Derivation hint: The optimal Krylov polynomial that suppresses $\{\lambda_2,\dots,\lambda_n\}$ while preserving $\lambda_1$ is the shifted/scaled Chebyshev polynomial on $[\lambda_n,\lambda_2]$ ; its exponential growth outside that interval, $C_{m-1}(1+2\gamma)\sim \tfrac12 e^{2(m-1)\sqrt\gamma}$ , drives geometric convergence governed by the spectral gap.

∎

mathematician · 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.

🧠 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 $|\lambda|$ direction outgrows all others, but it throws away every intermediate vector. The Krylov subspace keeps those intermediate vectors. Within the $m$ -dimensional span $\{v, Av, \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 $\gamma$ 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.

Figure

Krylov convergence: Chebyshev amplification of the extremal eigenvalue

ExampleLanczos on a 4×4 tridiagonal matrix

Given. The 1D discrete Laplacian (with sign convention $A=-\Delta_h$ ),

A=\begin{pmatrix} 2&-1&0&0\\ -1&2&-1&0\\ 0&-1&2&-1\\ 0&0&-1&2\end{pmatrix},

whose exact eigenvalues are $\lambda_k = 2-2\cos\frac{k\pi}{5} = \{0.382, 1.382, 2.618, 3.618\}$ . Start vector $v_1=\tfrac12(1,1,1,1)^{\mathsf T}$ .

Steps. 1. $Av_1 = \tfrac12(1,0,0,1)^{\mathsf T}$ , so $\alpha_1 = v_1^* A v_1 = \tfrac14(1+0+0+1)=0.5$ . 2. $w_1 = Av_1 - \alpha_1 v_1 = \tfrac12(1,0,0,1) - \tfrac14(1,1,1,1) = (0.25,-0.25,-0.25,0.25)$ . Then $\beta_1=\|w_1\|=0.5$ , $v_2=(0.5,-0.5,-0.5,0.5)^{\mathsf T}$ . 3. $Av_2 = (1.5,-1,-1,1.5)^{\mathsf T}\cdot\!$ ... computing: $Av_2=(1.5,-1.0,-1.0,1.5)^{\mathsf T}$ ; $\alpha_2 = v_2^*Av_2 = 0.5(1.5)+(-0.5)(-1)+(-0.5)(-1)+0.5(1.5)=2.5$ . 4. $T_2=\begin{pmatrix}0.5 & 0.5\\ 0.5 & 2.5\end{pmatrix}$ , with Ritz values $\theta=\tfrac{3\pm\sqrt{4+1}}{2}=\{0.382,\,2.618\}$ .

Result. After just 2 Lanczos steps the Ritz values $\boxed{\theta_1=0.382,\ \theta_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.

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

Given. A 3D Poisson problem discretized on a $200\times200\times200$ grid gives $n = 8\times10^6$ unknowns; the 7-point stencil yields $\mathrm{nnz}(A)\approx 7n = 5.6\times10^7$ . We want the 10 smallest eigenvalues (lowest vibrational/diffusion modes).

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

Result. Lanczos finishes in $\boxed{\sim 7\times10^{11}\text{ flops},\ \sim 13\text{ GB}}$ versus $5\times10^{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.

⚠️ 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(nm^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.

CounterexamplePower 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 $\lambda_{\max}=1000$ , not what you want. Worse, even targeting the smallest via power iteration on $A^{-1}$ converges at rate $|\lambda_{\min}/\lambda_2|^k = (1/1.001)^k$ — agonizingly slow because the two smallest eigenvalues are clustered. The fix is shift-and-invert: apply $(A-\sigma I)^{-1}$ with $\sigma$ near the target, which maps nearby eigenvalues to huge well-separated values $1/(\lambda-\sigma)$ , 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 $\kappa(\lambda_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 $\Lambda_\epsilon(A) = \{z : \|(zI-A)^{-1}\| \ge \epsilon^{-1}\}$ are the honest diagnostic: they reveal how far eigenvalues move under perturbations of size $\epsilon$ . A small residual bounds the backward error, not the forward eigenvalue error.

connection

Iterative eigensolvers share their entire algorithmic core with iterative linear solvers. The Conjugate Gradient Method is Lanczos applied to solving $Ax=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).

🔬 Where this appears

⏳ History

- 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.

open-problem

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(\lambda)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.

Sparse Matrices and Iterative Eigensolvers

Formal Statement

Why It Holds / Derivation

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also