Where Lanczos & Arnoldi Matter

Introduction

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.

Try It

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.

<!-- Two-panel layout: spectrum canvas on the left, controls and lists on the right -->
<div id="app">
  <div id="left-panel">
    <p class="hint">Blue lines = true eigenvalues. Orange dots = Ritz values from the Lanczos subspace. Step through the iteration and watch them converge.</p>
    <div id="spectrum-wrap">
      <canvas id="spectrum" width="260" height="300"></canvas>
      <div id="legend">
        <span class="leg-true"></span> True eigenvalues
        <span class="leg-ritz"></span> Ritz values
      </div>
    </div>
  </div>
  <div id="right-panel">
    <div id="step-info">
      <span id="step-label">Step <b id="step-num">0</b> / <b id="step-max">0</b></span>
    </div>
    <div id="ritz-list-wrap">
      <div class="list-title">Ritz values</div>
      <ul id="ritz-list"></ul>
    </div>
    <div id="true-list-wrap">
      <div class="list-title">True eigenvalues</div>
      <ul id="true-list"></ul>
    </div>
    <div class="btns">
      <button id="btn-step" type="button">Step</button>
      <button id="btn-run" type="button">Run all</button>
      <button id="btn-reset" type="button" class="ghost">Reset</button>
    </div>
    <div id="status" class="status"></div>
  </div>
</div>

/* Global box-sizing reset and base typography */
* { 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; }
/* Spectrum canvas: vertical eigenvalue axis with true eig lines and Ritz dots */
#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; }
/* Step counter display */
#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; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Lanczos iteration demo: build a Krylov subspace K_k(A,b) = span{b, Ab, ..., A^(k-1)b},
// project A onto it to get a k×k tridiagonal T, and read off Ritz values (eigenvalues of T)
// as approximations to the true eigenvalues of A. Extreme eigenvalues converge first.

// True eigenvalues of our 6×6 symmetric demo matrix
// A = Q * diag(eigenvalues) * Q^T  where Q is a random orthogonal matrix
const TRUE_EIGS = [-7.5, -3.2, 0.8, 2.1, 5.4, 9.0];
const N = TRUE_EIGS.length; // dimension n=6

// Build A = Q * diag(eigs) * Q^T (symmetric by construction)
function buildMatrix(eigs) {
  // Construct a random orthogonal Q via Gram-Schmidt on seeded random vectors
  const n = eigs.length;
  // Fill with seeded pseudo-random values (deterministic, no Math.random)
  const raw = [];
  for (let i = 0; i < n; i++) { raw.push([]); for (let j = 0; j < n; j++) raw[i].push(seeded(i*n+j)); }
  // Orthonormalize rows to get Q
  const Q = gramSchmidt(raw, n);
  // A = Q diag(eigs) Q^T
  const A = zeros(n);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
    let s = 0; for (let k = 0; k < n; k++) s += Q[i][k] * eigs[k] * Q[j][k];
    A[i][j] = s;
  }
  return A;
}

// Deterministic pseudo-random in [-0.5, 0.5] seeded by index i
function seeded(i) {
  const x = Math.sin(i * 9301 + 49297) * 233280;
  return x - Math.floor(x) - 0.5;
}
function zeros(n) { return Array.from({length:n},()=>Array(n).fill(0)); }
function dot(a,b){ let s=0; for(let i=0;i<a.length;i++) s+=a[i]*b[i]; return s; }
function norm(v){ return Math.sqrt(dot(v,v)); }
function scale(v,c){ return v.map(x=>x*c); }
function sub(a,b){ return a.map((x,i)=>x-b[i]); }

function gramSchmidt(vecs, n) {
  const Q = [];
  for (let j = 0; j < n; j++) {
    let v = vecs[j].slice();
    for (let i = 0; i < Q.length; i++) { const p = dot(v, Q[i]); v = sub(v, scale(Q[i], p)); }
    const nr = norm(v); Q.push(scale(v, 1/nr));
  }
  return Q;
}

// Matrix-vector product w = A*v (O(n^2) for dense, O(nnz) for sparse)
function matvec(A, v) {
  const n = A.length, w = Array(n).fill(0);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) w[i] += A[i][j] * v[j];
  return w;
}

// One Lanczos step: extend the Krylov basis by one vector
// Returns false if an invariant subspace is found (exact breakdown)
// Three-term recurrence: A*q_j = beta_{j-1}*q_{j-1} + alpha_j*q_j + beta_j*q_{j+1}
function lanczosStep(A, Qs, alphas, betas) {
  const j = Qs.length - 1;
  let w = matvec(A, Qs[j]);
  // Three-term recurrence: subtract alpha*q_j and beta_{j-1}*q_{j-1}
  if (j > 0) w = sub(w, scale(Qs[j-1], betas[betas.length-1]));
  const alpha = dot(w, Qs[j]);
  alphas.push(alpha);
  w = sub(w, scale(Qs[j], alpha));
  const beta = norm(w);
  betas.push(beta);
  if (beta < 1e-10) return false; // exact breakdown: invariant subspace found
  Qs.push(scale(w, 1/beta));
  return true;
}

// Compute eigenvalues of the k×k symmetric tridiagonal T via QR iteration
function tridiagEigs(alphas, betas_in) {
  // Simple QR iteration (sufficient for small k <= 6)
  const n = alphas.length;
  if (n === 1) return [alphas[0]];
  // Diagonal d and off-diagonal e of the tridiagonal
  const d = alphas.slice(), e = betas_in.slice(0, n-1);
  // Implicit symmetric QR step with Wilkinson shift — implicit symmetric QR with Wilkinson shift
  for (let iter = 0; iter < 30 * n; iter++) {
    // Find bottom of active submatrix by checking off-diagonal entries
    let m = n - 1;
    while (m > 0 && Math.abs(e[m-1]) < 1e-12 * (Math.abs(d[m-1]) + Math.abs(d[m]))) m--;
    if (m === 0) break;
    // Wilkinson shift: eigenvalue of bottom 2x2 closest to d[m]
    const a = d[m-1], b = e[m-1], c = d[m];
    const tr = (a + c) / 2;
    const disc = Math.sqrt((a - c) * (a - c) / 4 + b * b);
    const mu1 = tr + disc, mu2 = tr - disc;
    const mu = Math.abs(mu1 - c) < Math.abs(mu2 - c) ? mu1 : mu2;
    // Implicit QR step: chase the bulge with Givens rotations
    let x = d[0] - mu, y = e[0];
    for (let k = 0; k < m; k++) {
      const r = Math.sqrt(x * x + y * y);
      if (r < 1e-14) { if (k < m - 1) { x = e[k]; y = (k + 1 < m ? e[k + 1] : 0); } continue; }
      const cc = x / r, ss = y / r;
      if (k > 0) e[k - 1] = r;
      // Apply Givens rotation from both sides to maintain tridiagonal form
      const dk = d[k], dk1 = d[k + 1], ek = e[k];
      d[k]     = cc*cc*dk + 2*cc*ss*ek + ss*ss*dk1;
      d[k + 1] = ss*ss*dk - 2*cc*ss*ek + cc*cc*dk1;
      e[k]     = cc*ss*(dk1 - dk) + (cc*cc - ss*ss)*ek;
      if (k < m - 1) { x = e[k]; y = ss * e[k + 1]; e[k + 1] = cc * e[k + 1]; }
    }
  }
  return d.slice().sort((a, b) => a - b);
}

// Global Lanczos state: basis vectors, recurrence coefficients, Ritz history
const A = buildMatrix(TRUE_EIGS);
let Qs, alphas, betas, ritzHistory, done;

function initState() {
  // Normalized random starting vector b
  const b = Array.from({length:N},(_,i)=>seeded(i+100));
  const nb = norm(b);
  Qs = [scale(b, 1/nb)];
  alphas = []; betas = []; ritzHistory = [[]]; done = false;
}

initState();

// Canvas drawing: eigenvalue axis on the left, Ritz dots in the center
const cvs = document.getElementById('spectrum');
const ctx = cvs.getContext('2d');
const PAD = 30, W = cvs.width, H = cvs.height;
const MIN_EIG = -10, MAX_EIG = 11;

function eigToY(e) { return PAD + (H-2*PAD) * (1 - (e-MIN_EIG)/(MAX_EIG-MIN_EIG)); }

function draw(ritzVals) {
  ctx.clearRect(0,0,W,H);
  // Vertical axis
  ctx.strokeStyle = '#cdd9e3'; ctx.lineWidth = 1;
  ctx.beginPath(); ctx.moveTo(PAD,PAD); ctx.lineTo(PAD,H-PAD); ctx.stroke();

  // Axis tick marks and numeric labels
  ctx.fillStyle = '#888'; ctx.font = '11px ui-monospace, monospace'; ctx.textAlign = 'right';
  for (let e = -8; e <= 10; e += 2) {
    const y = eigToY(e);
    ctx.beginPath(); ctx.moveTo(PAD-4,y); ctx.lineTo(PAD,y); ctx.stroke();
    ctx.fillText(e, PAD-6, y+4);
  }
  ctx.fillStyle = '#aaa'; ctx.font = '10px system-ui'; ctx.textAlign = 'center';
  ctx.fillText('λ', PAD, PAD-8);

  // Draw horizontal dashed lines for each true eigenvalue
  ctx.strokeStyle = '#1d6fa5'; ctx.lineWidth = 2; ctx.setLineDash([4,3]);
  for (const e of TRUE_EIGS) {
    const y = eigToY(e);
    ctx.beginPath(); ctx.moveTo(PAD+4,y); ctx.lineTo(W-PAD,y); ctx.stroke();
  }
  ctx.setLineDash([]);

  // Draw orange circles for current Ritz values
  if (ritzVals && ritzVals.length > 0) {
    const cx = (PAD + W - PAD) / 2;
    ctx.fillStyle = '#f4900c';
    for (const r of ritzVals) {
      const y = eigToY(r);
      ctx.beginPath(); ctx.arc(cx, y, 5, 0, 2*Math.PI); ctx.fill();
    }
  }
}

function renderLists(ritzVals) {
  const rl = document.getElementById('ritz-list');
  const tl = document.getElementById('true-list');
  rl.innerHTML = '';
  tl.innerHTML = '';
  const EPS = 0.5; // convergence threshold for coloring
  if (ritzVals) {
    for (const r of ritzVals) {
      const li = document.createElement('li');
      const closest = TRUE_EIGS.reduce((best,e)=>Math.abs(e-r)<Math.abs(best-r)?e:best, Infinity);
      li.className = Math.abs(closest-r) < EPS ? 'converged' : 'far';
      li.textContent = r.toFixed(4);
      rl.appendChild(li);
    }
  }
  for (const e of TRUE_EIGS) {
    const li = document.createElement('li'); li.textContent = e.toFixed(4); tl.appendChild(li);
  }
}

function setStatus(t, cls) {
  const s = document.getElementById('status'); s.textContent = t; s.className = 'status ' + cls;
}

function updateStepLabel(cur) {
  document.getElementById('step-num').textContent = cur;
  document.getElementById('step-max').textContent = N;
}

function step() {
  if (done || Qs.length > N) { done = true; return; }
  const ok = lanczosStep(A, Qs, alphas, betas);
  const ritz = alphas.length >= 1 ? tridiagEigs(alphas, betas.slice(0, alphas.length-1)) : [];
  ritzHistory.push(ritz.slice());
  draw(ritz);
  renderLists(ritz);
  updateStepLabel(alphas.length);
  if (!ok || alphas.length >= N) {
    done = true;
    document.getElementById('btn-step').disabled = true;
    document.getElementById('btn-run').disabled = true;
    setStatus('✓ Lanczos complete — all eigenvalues converged.', 'ok');
  } else {
    setStatus(fmt('Added Krylov vector %length% of %N% — extreme eigenvalues converge first.', { length: alphas.length, N: N }), 'info');
  }
}

function runAll() {
  while (!done && alphas.length < N) step();
}

function reset() {
  initState();
  draw([]);
  renderLists([]);
  updateStepLabel(0);
  document.getElementById('btn-step').disabled = false;
  document.getElementById('btn-run').disabled = false;
  setStatus('', '');
}

document.getElementById('btn-step').addEventListener('click', step);
document.getElementById('btn-run').addEventListener('click', runAll);
document.getElementById('btn-reset').addEventListener('click', reset);

// Initial empty draw
draw([]);
renderLists([]);
updateStepLabel(0);
document.getElementById('step-max').textContent = N;

<!-- Two-panel layout: spectrum canvas on the left, controls and lists on the right -->
<div id="app">
  <div id="left-panel">
    <p class="hint">Blue lines = true eigenvalues. Orange dots = Ritz values from the Lanczos subspace. Step through the iteration and watch them converge.</p>
    <div id="spectrum-wrap">
      <canvas id="spectrum" width="260" height="300"></canvas>
      <div id="legend">
        <span class="leg-true"></span> True eigenvalues
        <span class="leg-ritz"></span> Ritz values
      </div>
    </div>
  </div>
  <div id="right-panel">
    <div id="step-info">
      <span id="step-label">Step <b id="step-num">0</b> / <b id="step-max">0</b></span>
    </div>
    <div id="ritz-list-wrap">
      <div class="list-title">Ritz values</div>
      <ul id="ritz-list"></ul>
    </div>
    <div id="true-list-wrap">
      <div class="list-title">True eigenvalues</div>
      <ul id="true-list"></ul>
    </div>
    <div class="btns">
      <button id="btn-step" type="button">Step</button>
      <button id="btn-run" type="button">Run all</button>
      <button id="btn-reset" type="button" class="ghost">Reset</button>
    </div>
    <div id="status" class="status"></div>
  </div>
</div>

/* Global box-sizing reset and base typography */
* { 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; }
/* Spectrum canvas: vertical eigenvalue axis with true eig lines and Ritz dots */
#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; }
/* Step counter display */
#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; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// Lanczos iteration demo: build a Krylov subspace K_k(A,b) = span{b, Ab, ..., A^(k-1)b},
// project A onto it to get a k×k tridiagonal T, and read off Ritz values (eigenvalues of T)
// as approximations to the true eigenvalues of A. Extreme eigenvalues converge first.

// True eigenvalues of our 6×6 symmetric demo matrix
// A = Q * diag(eigenvalues) * Q^T  where Q is a random orthogonal matrix
const TRUE_EIGS = [-7.5, -3.2, 0.8, 2.1, 5.4, 9.0];
const N = TRUE_EIGS.length; // dimension n=6

// Build A = Q * diag(eigs) * Q^T (symmetric by construction)
function buildMatrix(eigs) {
  // Construct a random orthogonal Q via Gram-Schmidt on seeded random vectors
  const n = eigs.length;
  // Fill with seeded pseudo-random values (deterministic, no Math.random)
  const raw = [];
  for (let i = 0; i < n; i++) { raw.push([]); for (let j = 0; j < n; j++) raw[i].push(seeded(i*n+j)); }
  // Orthonormalize rows to get Q
  const Q = gramSchmidt(raw, n);
  // A = Q diag(eigs) Q^T
  const A = zeros(n);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
    let s = 0; for (let k = 0; k < n; k++) s += Q[i][k] * eigs[k] * Q[j][k];
    A[i][j] = s;
  }
  return A;
}

// Deterministic pseudo-random in [-0.5, 0.5] seeded by index i
function seeded(i) {
  const x = Math.sin(i * 9301 + 49297) * 233280;
  return x - Math.floor(x) - 0.5;
}
function zeros(n) { return Array.from({length:n},()=>Array(n).fill(0)); }
function dot(a,b){ let s=0; for(let i=0;i<a.length;i++) s+=a[i]*b[i]; return s; }
function norm(v){ return Math.sqrt(dot(v,v)); }
function scale(v,c){ return v.map(x=>x*c); }
function sub(a,b){ return a.map((x,i)=>x-b[i]); }

function gramSchmidt(vecs, n) {
  const Q = [];
  for (let j = 0; j < n; j++) {
    let v = vecs[j].slice();
    for (let i = 0; i < Q.length; i++) { const p = dot(v, Q[i]); v = sub(v, scale(Q[i], p)); }
    const nr = norm(v); Q.push(scale(v, 1/nr));
  }
  return Q;
}

// Matrix-vector product w = A*v (O(n^2) for dense, O(nnz) for sparse)
function matvec(A, v) {
  const n = A.length, w = Array(n).fill(0);
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) w[i] += A[i][j] * v[j];
  return w;
}

// One Lanczos step: extend the Krylov basis by one vector
// Returns false if an invariant subspace is found (exact breakdown)
// Three-term recurrence: A*q_j = beta_{j-1}*q_{j-1} + alpha_j*q_j + beta_j*q_{j+1}
function lanczosStep(A, Qs, alphas, betas) {
  const j = Qs.length - 1;
  let w = matvec(A, Qs[j]);
  // Three-term recurrence: subtract alpha*q_j and beta_{j-1}*q_{j-1}
  if (j > 0) w = sub(w, scale(Qs[j-1], betas[betas.length-1]));
  const alpha = dot(w, Qs[j]);
  alphas.push(alpha);
  w = sub(w, scale(Qs[j], alpha));
  const beta = norm(w);
  betas.push(beta);
  if (beta < 1e-10) return false; // exact breakdown: invariant subspace found
  Qs.push(scale(w, 1/beta));
  return true;
}

// Compute eigenvalues of the k×k symmetric tridiagonal T via QR iteration
function tridiagEigs(alphas, betas_in) {
  // Simple QR iteration (sufficient for small k <= 6)
  const n = alphas.length;
  if (n === 1) return [alphas[0]];
  // Diagonal d and off-diagonal e of the tridiagonal
  const d = alphas.slice(), e = betas_in.slice(0, n-1);
  // Implicit symmetric QR step with Wilkinson shift — implicit symmetric QR with Wilkinson shift
  for (let iter = 0; iter < 30 * n; iter++) {
    // Find bottom of active submatrix by checking off-diagonal entries
    let m = n - 1;
    while (m > 0 && Math.abs(e[m-1]) < 1e-12 * (Math.abs(d[m-1]) + Math.abs(d[m]))) m--;
    if (m === 0) break;
    // Wilkinson shift: eigenvalue of bottom 2x2 closest to d[m]
    const a = d[m-1], b = e[m-1], c = d[m];
    const tr = (a + c) / 2;
    const disc = Math.sqrt((a - c) * (a - c) / 4 + b * b);
    const mu1 = tr + disc, mu2 = tr - disc;
    const mu = Math.abs(mu1 - c) < Math.abs(mu2 - c) ? mu1 : mu2;
    // Implicit QR step: chase the bulge with Givens rotations
    let x = d[0] - mu, y = e[0];
    for (let k = 0; k < m; k++) {
      const r = Math.sqrt(x * x + y * y);
      if (r < 1e-14) { if (k < m - 1) { x = e[k]; y = (k + 1 < m ? e[k + 1] : 0); } continue; }
      const cc = x / r, ss = y / r;
      if (k > 0) e[k - 1] = r;
      // Apply Givens rotation from both sides to maintain tridiagonal form
      const dk = d[k], dk1 = d[k + 1], ek = e[k];
      d[k]     = cc*cc*dk + 2*cc*ss*ek + ss*ss*dk1;
      d[k + 1] = ss*ss*dk - 2*cc*ss*ek + cc*cc*dk1;
      e[k]     = cc*ss*(dk1 - dk) + (cc*cc - ss*ss)*ek;
      if (k < m - 1) { x = e[k]; y = ss * e[k + 1]; e[k + 1] = cc * e[k + 1]; }
    }
  }
  return d.slice().sort((a, b) => a - b);
}

// Global Lanczos state: basis vectors, recurrence coefficients, Ritz history
const A = buildMatrix(TRUE_EIGS);
let Qs, alphas, betas, ritzHistory, done;

function initState() {
  // Normalized random starting vector b
  const b = Array.from({length:N},(_,i)=>seeded(i+100));
  const nb = norm(b);
  Qs = [scale(b, 1/nb)];
  alphas = []; betas = []; ritzHistory = [[]]; done = false;
}

initState();

// Canvas drawing: eigenvalue axis on the left, Ritz dots in the center
const cvs = document.getElementById('spectrum');
const ctx = cvs.getContext('2d');
const PAD = 30, W = cvs.width, H = cvs.height;
const MIN_EIG = -10, MAX_EIG = 11;

function eigToY(e) { return PAD + (H-2*PAD) * (1 - (e-MIN_EIG)/(MAX_EIG-MIN_EIG)); }

function draw(ritzVals) {
  ctx.clearRect(0,0,W,H);
  // Vertical axis
  ctx.strokeStyle = '#cdd9e3'; ctx.lineWidth = 1;
  ctx.beginPath(); ctx.moveTo(PAD,PAD); ctx.lineTo(PAD,H-PAD); ctx.stroke();

  // Axis tick marks and numeric labels
  ctx.fillStyle = '#888'; ctx.font = '11px ui-monospace, monospace'; ctx.textAlign = 'right';
  for (let e = -8; e <= 10; e += 2) {
    const y = eigToY(e);
    ctx.beginPath(); ctx.moveTo(PAD-4,y); ctx.lineTo(PAD,y); ctx.stroke();
    ctx.fillText(e, PAD-6, y+4);
  }
  ctx.fillStyle = '#aaa'; ctx.font = '10px system-ui'; ctx.textAlign = 'center';
  ctx.fillText('λ', PAD, PAD-8);

  // Draw horizontal dashed lines for each true eigenvalue
  ctx.strokeStyle = '#1d6fa5'; ctx.lineWidth = 2; ctx.setLineDash([4,3]);
  for (const e of TRUE_EIGS) {
    const y = eigToY(e);
    ctx.beginPath(); ctx.moveTo(PAD+4,y); ctx.lineTo(W-PAD,y); ctx.stroke();
  }
  ctx.setLineDash([]);

  // Draw orange circles for current Ritz values
  if (ritzVals && ritzVals.length > 0) {
    const cx = (PAD + W - PAD) / 2;
    ctx.fillStyle = '#f4900c';
    for (const r of ritzVals) {
      const y = eigToY(r);
      ctx.beginPath(); ctx.arc(cx, y, 5, 0, 2*Math.PI); ctx.fill();
    }
  }
}

function renderLists(ritzVals) {
  const rl = document.getElementById('ritz-list');
  const tl = document.getElementById('true-list');
  rl.innerHTML = '';
  tl.innerHTML = '';
  const EPS = 0.5; // convergence threshold for coloring
  if (ritzVals) {
    for (const r of ritzVals) {
      const li = document.createElement('li');
      const closest = TRUE_EIGS.reduce((best,e)=>Math.abs(e-r)<Math.abs(best-r)?e:best, Infinity);
      li.className = Math.abs(closest-r) < EPS ? 'converged' : 'far';
      li.textContent = r.toFixed(4);
      rl.appendChild(li);
    }
  }
  for (const e of TRUE_EIGS) {
    const li = document.createElement('li'); li.textContent = e.toFixed(4); tl.appendChild(li);
  }
}

function setStatus(t, cls) {
  const s = document.getElementById('status'); s.textContent = t; s.className = 'status ' + cls;
}

function updateStepLabel(cur) {
  document.getElementById('step-num').textContent = cur;
  document.getElementById('step-max').textContent = N;
}

function step() {
  if (done || Qs.length > N) { done = true; return; }
  const ok = lanczosStep(A, Qs, alphas, betas);
  const ritz = alphas.length >= 1 ? tridiagEigs(alphas, betas.slice(0, alphas.length-1)) : [];
  ritzHistory.push(ritz.slice());
  draw(ritz);
  renderLists(ritz);
  updateStepLabel(alphas.length);
  if (!ok || alphas.length >= N) {
    done = true;
    document.getElementById('btn-step').disabled = true;
    document.getElementById('btn-run').disabled = true;
    setStatus('✓ Lanczos complete — all eigenvalues converged.', 'ok');
  } else {
    setStatus(fmt('Added Krylov vector %length% of %N% — extreme eigenvalues converge first.', { length: alphas.length, N: N }), 'info');
  }
}

function runAll() {
  while (!done && alphas.length < N) step();
}

function reset() {
  initState();
  draw([]);
  renderLists([]);
  updateStepLabel(0);
  document.getElementById('btn-step').disabled = false;
  document.getElementById('btn-run').disabled = false;
  setStatus('', '');
}

document.getElementById('btn-step').addEventListener('click', step);
document.getElementById('btn-run').addEventListener('click', runAll);
document.getElementById('btn-reset').addEventListener('click', reset);

// Initial empty draw
draw([]);
renderLists([]);
updateStepLabel(0);
document.getElementById('step-max').textContent = N;

Result

Console

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.

The Real Complexity

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.

Where It Matters

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.

Conclusion

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.

Comments

Loading comments...

https://www.kipuhub.com/en/article/lanczos-arnoldi/Content licensed under CC BY-NC 4.0.

Cookies & privacy

Lanczos & Arnoldi

The art of listening to a matrix without reading it all