Lanczos & Arnoldi

Eigenvalues are everywhere. The natural resonant frequencies of a bridge, the most influential pages on the web, the ground-state energy of a molecule — all of them reduce to finding the eigenvalues of some matrix $A$. The trouble is that real matrices are enormous: structural models of a car chassis can involve millions of degrees of freedom, and the matrix that describes the web graph has billions of rows.

Computing all eigenvalues of an $n \times n$ matrix costs $O(n^3)$ time and $O(n^2)$ memory. For $n = 10^6$ that is completely out of reach.

The saving insight, discovered independently by Cornelius Lanczos (1950) for symmetric matrices and extended by Walter Arnoldi (1951) to general ones, is that you almost never need all eigenvalues. You need the extremal ones — the largest, smallest, or those closest to some target. And those happen to be well-approximated by a very small Krylov subspace: the span of $\{b, Ab, A^2 b, \dots, A^{k-1} b\}$ for a starting vector $b$ and a step count $k \ll n$.

At each step you multiply by $A$ once (cheap for sparse matrices) and orthogonalize against all previous vectors. After $k$ steps you have an $n \times k$ basis $Q$ and a tiny $k \times k$ tridiagonal (Lanczos) or Hessenberg (Arnoldi) matrix $H$. The eigenvalues of $H$ — called Ritz values — converge to the true eigenvalues of $A$ at a rate that depends on how well-separated those eigenvalues are.

The algorithm belongs to the field of numerical linear algebra, and its story touches the same themes as dimensionality reduction: huge data spaces collapsed into manageable projections that still capture the essential structure.

The demo below runs the Lanczos iteration on a small $6 \times 6$ symmetric matrix whose true eigenvalues are shown on the right. Click Step to add one Krylov vector, or Run all to complete the iteration. Watch the Ritz values (orange dots) converge toward the true eigenvalues (blue lines) as the subspace grows.

<!-- {{c_layout_comment}} -->
<div id="app">
  <div id="left-panel">
    <p class="hint">{{hint_para}}</p>
    <div id="spectrum-wrap">
      <canvas id="spectrum" width="260" height="300"></canvas>
      <div id="legend">
        <span class="leg-true"></span> {{legend_true}}
        <span class="leg-ritz"></span> {{legend_ritz}}
      </div>
    </div>
  </div>
  <div id="right-panel">
    <div id="step-info">
      <span id="step-label">{{step_label}} <b id="step-num">0</b> / <b id="step-max">0</b></span>
    </div>
    <div id="ritz-list-wrap">
      <div class="list-title">{{ritz_title}}</div>
      <ul id="ritz-list"></ul>
    </div>
    <div id="true-list-wrap">
      <div class="list-title">{{true_title}}</div>
      <ul id="true-list"></ul>
    </div>
    <div class="btns">
      <button id="btn-step" type="button">{{btn_step}}</button>
      <button id="btn-run" type="button">{{btn_run}}</button>
      <button id="btn-reset" type="button" class="ghost">{{btn_reset}}</button>
    </div>
    <div id="status" class="status"></div>
  </div>
</div>

/* {{c_global_box}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; padding: 14px; }
#app { display: flex; gap: 1rem; padding: .5rem; flex-wrap: wrap; }
#left-panel { flex: 0 0 auto; }
#right-panel { flex: 1 1 160px; min-width: 140px; display: flex; flex-direction: column; gap: .5rem; }
.hint { font-size: .85rem; color: #444; margin: 0 0 .5rem; line-height: 1.4; max-width: 260px; }
/* {{c_canvas_comment}} */
#spectrum-wrap { position: relative; }
canvas { display: block; border: 1px solid #cdd9e3; border-radius: 8px; background: #f8fafc; }
#legend { display: flex; gap: .8rem; align-items: center; font-size: .78rem; color: #555; margin-top: .3rem; }
.leg-true { display: inline-block; width: 22px; height: 2px; background: #1d6fa5; vertical-align: middle; }
.leg-ritz { display: inline-block; width: 10px; height: 10px; border-radius: 50%; background: #f4900c; vertical-align: middle; }
/* {{c_step_info}} */
#step-info { font-size: .9rem; color: #555; }
.list-title { font-size: .78rem; font-weight: 700; color: #666; text-transform: uppercase; letter-spacing: .04em; }
ul { list-style: none; margin: 0; padding: 0; }
li { font: .82rem ui-monospace, monospace; padding: .1rem 0; }
li.converged { color: #0a7d33; }
li.far { color: #888; }
.btns { display: flex; gap: .4rem; flex-wrap: wrap; }
button { font: 600 13px system-ui, sans-serif; padding: .4rem .8rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
.status { font-size: .9rem; font-weight: 600; min-height: 1.3em; }
.status.ok { color: #0a7d33; }
.status.info { color: #1d6fa5; }

// Code not found

Notice that the extreme eigenvalues (largest and smallest) converge first — the Ritz values at the ends of the spectrum lock on quickly, while the middle ones need more steps. This is the fundamental property that makes Lanczos and Arnoldi practical: you stop as soon as the eigenvalues you care about have converged.

How fast does this work, and how hard is the underlying problem?

• Cost per step. Each Lanczos/Arnoldi step does one matrix–vector product $Av$ and one round of orthogonalization. If $A$ is sparse with $\text{nnz}$ non-zeros, the product costs $O(\text{nnz})$. Orthogonalizing against the $k$ previous vectors costs $O(kn)$. After $k$ steps the total work is $O(k(\text{nnz} + kn))$, far below the $O(n^3)$ of a dense eigensolver.

• Convergence rate. The Ritz value for the largest eigenvalue $\lambda_1$ converges at a geometric rate controlled by the ratio $(\lambda_1 - \lambda_2)/(\lambda_1 - \lambda_n)$: the more isolated $\lambda_1$ is, the faster convergence. Formally the error after $k$ steps satisfies

  $$|\lambda_1 - \theta_1^{(k)}| \leq C \cdot \left(\frac{\lambda_1 - \lambda_2}{\lambda_1 - \lambda_n}\right)^{2k}$$

  where $\theta_1^{(k)}$ is the leading Ritz value and $C$ depends on the starting vector. This is the same polynomial-approximation argument behind Chebyshev acceleration.

• Finite precision is the catch. In exact arithmetic, Lanczos produces a perfectly orthogonal basis and terminates in at most $n$ steps. In floating-point, the basis loses orthogonality ("ghost eigenvalues" appear). The fix — full re-orthogonalization (Gram–Schmidt at every step) or selective re-orthogonalization — restores reliability at the cost of $O(k^2 n)$ extra work. Modern implementations like ARPACK use implicit restarts to keep $k$ small.

• Arnoldi vs Lanczos. For a non-symmetric matrix, the three-term recurrence of Lanczos breaks down; Arnoldi keeps a full Hessenberg matrix and is more expensive ($O(k^2 n)$ orthogonalization) but always numerically stable. The eigenvalues of $A$ and $A^T$ can differ, so Arnoldi must track both left and right eigenvectors.

The eigenproblem itself is solvable in polynomial time ($O(n^3)$ via the QR algorithm), so there is no NP-hardness here. The challenge is purely about scale: squeezing useful answers out of matrices whose size makes dense methods impossible. Krylov methods are the dominant answer.

Krylov eigensolvers are invisible infrastructure in almost every field that does large-scale computation:

• Search and recommendation. Google's original PageRank is the dominant eigenvector of a billion-node web graph. The power iteration used for PageRank is the simplest Krylov method; HITS and spectral clustering use the same machinery.
• Quantum chemistry. The ground-state energy of a molecule is the smallest eigenvalue of a Hamiltonian matrix. Configuration-interaction methods produce matrices with $10^{10}$ or more rows; the Davidson method (a variant of Lanczos) extracts the few lowest eigenvalues in hours rather than centuries.
• Structural engineering. Vibration analysis of a building, airplane wing, or car chassis requires the lowest natural frequencies (eigenvalues of the stiffness/mass pencil). FEA software ships Lanczos as a built-in routine.
• Machine learning. Spectral methods for dimensionality reduction (PCA, spectral embedding, diffusion maps) need the top-$k$ eigenvectors of a covariance or graph-Laplacian matrix. At the scale of ImageNet, iterative methods are the only option.
• Signal processing and control. Stability of a dynamical system depends on the eigenvalues of its state-transition matrix. Large discretized PDEs require the rightmost (potentially destabilizing) eigenvalues — exactly what Arnoldi finds first.

The ARPACK library, built on implicitly restarted Arnoldi/Lanczos, underpins SciPy's eigs, MATLAB's eigs, Julia's Arpack.jl, and countless other scientific computing stacks. When you call eigs(A, k), you are almost certainly running the Arnoldi iteration inside.

Lanczos and Arnoldi answer a question that sounds impossible at first: how do you find the most important eigenvalues of a matrix you can barely store in memory? Their answer is to never look at the whole matrix — just multiply by it repeatedly, accumulating information in a tiny Krylov subspace, and let the Ritz values converge.

What makes these algorithms remarkable is their economy. They touch the matrix only through matrix–vector products, which for sparse $A$ costs almost nothing per step. They terminate in far fewer steps than the matrix dimension, precisely because the eigenvalues you care about are isolated enough to stand out.

From PageRank to quantum chemistry to structural resonance, the Lanczos and Arnoldi iterations quietly underlie some of the most computationally demanding tasks in science and engineering — a vivid reminder that the right mathematical structure can turn an intractable problem into a routine one. For the bigger picture of when computation becomes easy or hard, see the P vs NP article.