Lanczos Algorithm

From Arnoldi to a three-term recurrence

Arnoldi iteration builds an orthonormal Krylov basis Q = [q_1 | q_2 | ...] satisfying A·Q_k = Q_{k+1}·H̃_k, where H̃_k is upper Hessenberg. When A is symmetric, the projected matrix H_k = Q_kᵀ·A·Q_k must inherit that symmetry: H_kᵀ = (Q_kᵀ·A·Q_k)ᵀ = Q_kᵀ·Aᵀ·Q_k = Q_kᵀ·A·Q_k = H_k. So H_k is both upper Hessenberg and symmetric — which forces it to be tridiagonal. Everything outside the main diagonal and the two adjacent diagonals must be zero.

Tridiagonal structure has a dramatic consequence. The Arnoldi step

w = A·q_j - sum_{i=1}^{j} h_{ij}·q_i

becomes

w = A·q_j - α_j·q_j - β_{j-1}·q_{j-1}

because all the other h_{ij} for i < j-1 are zero. Only three terms remain. Computing q_{j+1} needs q_j (the current basis vector), q_{j-1} (the previous one), and one matrix-vector product. You can forget q_1, ..., q_{j-2} entirely — they're automatically orthogonal to q_{j+1} by the recurrence.

Two scalars carry all the new spectral information: α_j = q_jᵀ·A·q_j (the diagonal entry of T) and β_j = ‖w‖ (the sub-diagonal entry of T). Everything else in T is the symmetric reflection (β_j again, above the diagonal). The full algorithm fits in five lines.

The algorithm

q_0 = 0,  β_0 = 0,  q_1 = v / ‖v‖
for j = 1, 2, ..., m:
    w = A · q_j                                # one matrix-vector product
    α_j = q_jᵀ · w                              # diagonal entry of T
    w = w - α_j · q_j - β_{j-1} · q_{j-1}      # three-term recurrence
    β_j = ‖w‖                                   # sub-diagonal entry
    if β_j == 0: break                          # invariant subspace
    q_{j+1} = w / β_j

After m steps you have the tridiagonal matrix

T_m = [[ α_1   β_1                          ]
       [ β_1   α_2   β_2                    ]
       [       β_2   α_3   β_3              ]
       [             β_3   α_4   ...        ]
       [                   ...        β_{m-1}]
       [                         β_{m-1}  α_m ]]

Eigenvalues of T_m (Ritz values) approximate extremal eigenvalues of A. The tridiagonal eigenvalue problem itself is cheap — O(m²) via the QR algorithm specialised to tridiagonal matrices — so the total work is dominated by m sparse matrix-vector products.

Why O(N) storage matters

Arnoldi must store all k basis vectors to orthogonalise against them: O(kN) memory. For a million-dimensional problem with k = 100, that's 100 million doubles = 800 MB just for the basis — uncomfortably close to a single workstation's RAM, intolerable if you need to scale.

Lanczos keeps only three vectors: q_{j-1}, q_j, and the work vector w. Plus 2k scalars for the tridiagonal entries. The total is O(N) + O(k), independent of how many Lanczos steps you take. For a million-dimensional problem that's 24 MB for the vectors plus negligible bookkeeping — trivial on any machine. This is what makes Lanczos the standard tool for sparse symmetric eigenvalue problems with N in the millions to billions.

One caveat: if you want to recover eigenvectors (not just eigenvalues), you need either to save all the q's anyway (giving up the O(N) win) or to rerun the iteration after solving the tridiagonal eigenvalue problem (cheap because you know α and β). In practice, "two-pass" Lanczos saves α and β on the first pass, solves the tridiagonal problem to identify the desired Ritz value, then reruns the Lanczos recurrence to extract the corresponding Ritz vector.

Worked example — a 4×4 symmetric matrix

Take the symmetric matrix

A = [[4, 1, 0, 0],
     [1, 3, 1, 0],
     [0, 1, 2, 1],
     [0, 0, 1, 1]]

This is already nearly tridiagonal — a useful test case. Eigenvalues are approximately {4.86, 3.20, 1.46, 0.48}. Run Lanczos starting from q_1 = [1, 1, 1, 1]ᵀ / 2:

step j=1:
  w = A·q_1 = [2.5, 2.5, 2.0, 1.0]ᵀ
  α_1 = q_1ᵀ·w = 0.5·2.5 + 0.5·2.5 + 0.5·2.0 + 0.5·1.0 = 4.0
  w = w - 4.0·q_1 = [0.5, 0.5, 0.0, -1.0]ᵀ
  β_1 = ‖w‖ = √(0.25 + 0.25 + 0 + 1.0) = √1.5 ≈ 1.225
  q_2 = w / 1.225 ≈ [0.408, 0.408, 0.0, -0.816]ᵀ

step j=2:
  w = A·q_2 ≈ [2.041, 1.224, 0.408, -0.816]ᵀ
  α_2 = q_2ᵀ·w ≈ 0.408·2.041 + 0.408·1.224 + 0·0.408 + (-0.816)·(-0.816)
              ≈ 0.832 + 0.499 + 0 + 0.666 ≈ 1.997
  w = w - 1.997·q_2 - 1.225·q_1 ≈ [0.612, -0.204, -0.204, -0.204]ᵀ
  β_2 = ‖w‖ ≈ 0.707
  q_3 = w / 0.707 ≈ [0.866, -0.289, -0.289, -0.289]ᵀ

(continue to j=3 to fully tridiagonalise; eigenvalues of T_3 already approximate 4.86 and 0.48 well)

T_3 = [[ 4.0     1.225            ]
       [ 1.225   1.997   0.707   ]
       [         0.707   ...     ]]

Eigenvalues of T_3 (the 3×3 tridiagonal projected matrix) bracket the extremes of A's spectrum. After just three Lanczos steps the largest and smallest Ritz values are already within ~1% of the true extremes; one or two more steps gives full convergence to all four eigenvalues.

Lanczos vs Arnoldi vs CG vs MINRES

	Lanczos	Arnoldi	Conjugate Gradient	MINRES
Matrix type	Symmetric	General	Sym. positive definite	Symmetric indefinite
Solves	Eigenvalue problem	Eigenvalue problem	Linear system Ax=b	Linear system Ax=b
Per-step cost	1 mat-vec + O(N)	1 mat-vec + O(kN)	1 mat-vec + O(N)	1 mat-vec + O(N)
Storage	O(N) + O(k) scalars	O(kN)	O(N)	O(N)
Projected matrix	Tridiagonal	Hessenberg	Tridiagonal (implicit)	Tridiagonal (implicit)
Best when	Few eigenvalues, sparse sym A	Few eigenvalues, sparse general A	Linear system, sym PD A	Linear system, sym indefinite

Lanczos and CG are closely related — both are Krylov methods on symmetric matrices, and their iterates lie in the same subspace. CG implicitly does Lanczos on the same Krylov basis, just shaping the iteration around minimising the energy norm of the residual instead of extracting eigenvalues. MINRES is the symmetric-indefinite generalisation of CG; it too is mathematically equivalent to Lanczos plus a least-squares projection.

Loss of orthogonality and how to fix it

The three-term recurrence guarantees orthogonality of the q's in exact arithmetic. But in floating point, rounding errors break this property in a very specific way: once a Ritz value θ_i converges (matches an eigenvalue of A within machine precision), the corresponding Ritz vector y_i has approximately the same direction as the true eigenvector of A. The next Lanczos vector ends up with a non-trivial component along y_i — orthogonality with earlier basis vectors is lost — and the iteration starts to produce duplicates of the converged eigenvector instead of finding new ones.

Christopher Paige's 1971 PhD thesis at the University of London analysed this precisely and showed it's exactly the same phenomenon for every Lanczos run. Three fixes are now standard:

• Full reorthogonalisation. At each step, re-Gram–Schmidt against all previous basis vectors. Gives bulletproof orthogonality but kills the O(N) per-step cost — you're back to Arnoldi's O(kN). Used for small to medium problems where reliability matters more than scaling.
• Selective reorthogonalisation (Parlett–Scott, 1979). Reorthogonalise only against the converged Ritz vectors, not the whole basis. Much cheaper: typically a handful of vectors, not the full basis. Maintains good orthogonality while keeping near-O(N) cost.
• Partial reorthogonalisation (Simon, 1984). Maintain a running estimate of orthogonality loss using a recurrence on the rounding-error magnitudes. Reorthogonalise only when the loss exceeds √ε. The standard choice in modern codes — ARPACK and PRIMME use variants.

Where Lanczos shows up

• Quantum many-body physics. The Hubbard model, the Heisenberg model, and ab initio nuclear shell-model calculations all reduce to finding the lowest few eigenvalues of huge sparse symmetric Hamiltonians (N up to 10¹⁰ or more). Lanczos is the universal tool — sometimes block Lanczos with implicit restart, sometimes plain Lanczos with selective reorthogonalisation.
• ARPACK and PRIMME. The implementations behind MATLAB's eigs and SciPy's scipy.sparse.linalg.eigsh. ARPACK uses implicitly restarted Lanczos; PRIMME (Stathopoulos & McCombs) uses block Davidson and GD+k, both Lanczos-based.
• Graph spectra. Laplacian eigenvalues of huge networks (10⁸ nodes or more). The Fiedler vector (second-smallest eigenvalue) is computed by Lanczos with shift-and-invert; used in spectral clustering, graph partitioning, and graph drawing.
• Density functional theory. Kohn–Sham equations are eigenvalue problems on the Hamiltonian operator discretised on a grid or basis set. Codes like Quantum ESPRESSO and VASP use block Lanczos or Davidson (a Lanczos variant) to find the occupied orbitals.
• Singular value decomposition of large sparse matrices. SVD of A reduces to an eigenvalue problem on AᵀA (or the bidiagonal Golub–Kahan form). Both can be solved by Lanczos applied to the symmetric system, giving truncated SVD without forming AᵀA explicitly.
• Vibration analysis in structural mechanics. Eigenvalues of the stiffness/mass matrix pair (a generalised eigenvalue problem) give natural frequencies of buildings, bridges, aircraft. ANSYS, Abaqus, and Nastran all use block Lanczos for the modal analysis.

Common pitfalls

• Believing the orthogonality guarantee in finite precision. Plain three-term Lanczos loses orthogonality after the first Ritz value converges. Use selective or partial reorthogonalisation, or be very careful about how you interpret the output.
• Targeting interior eigenvalues without shift-and-invert. Plain Lanczos converges to extremal eigenvalues — it ignores the interior. For an interior λ near σ, apply Lanczos to (A − σI)⁻¹ and pay for a linear solve per step.
• Cluster eigenvalues confuse single-vector Lanczos. If two eigenvalues are very close, plain Lanczos may merge them into one Ritz value or miss one entirely. Block Lanczos with a block size larger than the cluster solves this — at the cost of multiple matrix-vector products per step.
• Forgetting to recompute eigenvectors. One-pass Lanczos discards the basis vectors and gives you only α and β. To recover an eigenvector you need to redo the iteration (cheap because you know the structure) or save all basis vectors during the run (sacrificing O(N) storage).
• Underestimating convergence variability. Lanczos converges at a rate set by the relative spectral gap of the target eigenvalue. Well-separated extremes converge in a few dozen steps; clustered eigenvalues can take hundreds. Don't hard-code an iteration count without a residual-based stopping test.

Cornelius Lanczos and the 1950 paper

Lanczos developed the algorithm at the U.S. National Bureau of Standards in the late 1940s, where he had emigrated after wartime work in Britain. The 1950 paper "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators" introduces both the symmetric tridiagonal recurrence and (in a related section) the non-symmetric two-sided variant now usually attributed to a separate "bi-Lanczos" tradition. The same paper proves the orthogonality and bracketing properties of Ritz values that make the method work.

The method sat almost untouched for two decades. Early computers had short mantissas (often 8–10 decimal digits) and the orthogonality loss was catastrophic in practice. Christopher Paige's 1971 PhD thesis at the University of London diagnosed exactly why Lanczos broke down in finite precision — he showed that the loss of orthogonality is not a bug but a precise mathematical consequence of Ritz value convergence — and Beresford Parlett's group at Berkeley provided the selective-reorthogonalisation fixes that made the method practical. Today Lanczos is the most-cited algorithm in computational physics; an estimated 30–40% of all CPU cycles spent on quantum chemistry codes run inside some Lanczos variant.