The Nyström Method: Rebuild a Giant Kernel Matrix From a Handful of Landmark Points

Every kernel method — SVM, Gaussian processes, kernel PCA, spectral clustering — never looks at your raw features. It looks only at the n×n matrix K of pairwise similarities, K[i][j] = k(x_i, x_j). That matrix is a monster: n = 100,000 points is 10 billion entries — tens of gigabytes you can't store and an O(n³) solve you can't afford. The Nyström method buys its way out with one idea: sampling.

Sample landmarks, rebuild the whole thing

Pick m ≪ n landmark points. Compute only the thin n×m block C of kernel values from every point to the landmarks, and the tiny m×m block W among the landmarks themselves. Then reconstruct the entire matrix:

K ≈ C W⁺ Cᵀ

That's it. The whole giant is rebuilt from a strip of columns and one small pseudo-inverse. Both blocks are cheap — O(nm) and O(m²) — and you never form the full K:

function buildCW(xs, ys, land, n, m, gamma) {
  const C = new Float64Array(n * m), W = new Float64Array(m * m);
  for (let i = 0; i < n; i++)
    for (let a = 0; a < m; a++) {
      const l = land[a];
      C[i*m + a] = rbf(xs[i], ys[i], xs[l], ys[l], gamma);
    }
  for (let a = 0; a < m; a++)
    for (let b = 0; b < m; b++) {
      const la = land[a], lb = land[b];
      W[a*m + b] = rbf(xs[la], ys[la], xs[lb], ys[lb], gamma);
    }
  return { C, W };                       // never build the full n×n K
}

Why it works: kernels are nearly low-rank

A smooth RBF kernel's eigenvalues decay quickly, so the matrix is nearly low-rank — a handful of well-spread landmarks already capture almost all of it. And it's exact when the landmarks are all the points: at m = n, C = W = K and K W⁺ Kᵀ = K by the defining property of the pseudo-inverse. In between, the relative Frobenius error falls fast as m grows.

The W⁺ is a Moore–Penrose pseudo-inverse, not a plain inverse, because two nearly-identical landmarks give W a near-zero eigenvalue that would blow up. You diagonalise the tiny symmetric W, invert each eigenvalue above a tolerance, and zero the rest:

function pinvSym(W, m, jitter) {
  const A = Float64Array.from(W);
  for (let i = 0; i < m; i++) A[i*m + i] += jitter;    // ridge for stability
  const { vals, vecs } = jacobiEigen(A, m);
  let maxv = 0; for (let i = 0; i < m; i++) if (vals[i] > maxv) maxv = vals[i];
  const tol = maxv * 1e-9 + 1e-15;                      // rank cutoff
  const inv = new Float64Array(m*m);
  for (let a = 0; a < m; a++) for (let b = 0; b < m; b++) {
    let s = 0;
    for (let k = 0; k < m; k++) {
      const lam = vals[k];
      if (lam > tol) s += vecs[a*m+k] * vecs[b*m+k] / lam;   // 1/λ on kept modes
    }
    inv[a*m + b] = s;                                   // (V Λ⁺ Vᵀ)_ab
  }
  return inv;
}

Which points you pick matters

Uniform random landmarks are the classic baseline, but they clump in dense regions and waste budget. Farthest-point sampling — greedily add the point most distant from those already chosen — blankets the cloud, so for the same m it reconstructs the kernel far better.

Where it sits

Nyström is the same bet as random projections: approximate the expensive object so the algorithm scales. But it's data-dependent — it looks at your points and places landmarks where they matter — which is why it exploits a genuinely low-rank kernel so well and is exact at m = n. Random Fourier Features approximate the kernel too, but data-obliviously. In practice you rarely hand-roll it: scikit-learn ships Nystroem, which returns the feature map Φ = C W^(−1/2) so Φ Φᵀ ≈ K — feed Φ into any linear model to get kernel quality at O(n·m²) cost instead of O(n³). Gaussian-process libraries call the same idea "inducing points."

When n is huge and you need kernels, reach for Nyström. Watch the error-vs-m curve plunge toward zero, compare random against spread sampling, and read the full from-scratch build here: https://dev48v.infy.uk/ml/day55-nystrom.html