← Alex Meterez · Notes

How do you see the spectrum of a Hessian with a hundred million rows? You can't store it, let alone diagonalize it. This primer builds up, from scratch, the tool used in A Defense of the Quadratic Model: stochastic Lanczos quadrature. We start from moment matching, show why the direct solution (Prony's method) breaks, and arrive at Lanczos as the same idea in a better basis. Then we ask how good the answer is.

1. Motivation
2. The spectrum as a measure
3. Probing with a random vector
4. Moment matching
5. Prony's trick
6. The Krylov bridge
7. Prony breaks
8. Lanczos: same subspace, better basis
9. Gram–Schmidt as a three-term recurrence
10. Reading off nodes and weights
11. How deep can we see?
12. Error bands: Gauss–Radau
13. At scale
14. Appendix: why one probe suffices

1Motivation

• Given: a large model with parameters $\theta\in\mathbb R^d$.
• Want: the spectrum of the Hessian $H_\theta\in\mathbb R^{d\times d}$.
• Problem: with $d=100\mathrm M$ in fp32, $H_\theta$ takes roughly 40 PB. A single vector $Hv$ is only 400 MB.
• Access: Hessian-vector products (HVPs), computed by autograd without ever forming $H$: $$H_\theta v=\nabla_\theta\big\langle\nabla_\theta L(\theta),\,v\big\rangle\quad(v\ \text{fixed}),$$ which costs one forward and two backward passes, about $2$–$3\times$ a gradient.

The plan:

1. Setup: the spectrum as a measure; HVPs give its moments.
2. Moment matching: the best $m$-atom summary.
3. Prony: a direct solve (which breaks).
4. Lanczos: the same idea, stable.
5. How good is it? Head vs. tail, and error bands.

Running example throughout: $\lambda_i=i^{-1.5}$ for $i=1,\dots,30$ ($d=30$), a power-law spectrum spanning a $160\times$ range.

2The spectrum as a measure

• $H\in\mathbb R^{d\times d}$ symmetric, $\;H=\sum_{i=1}^d\lambda_i\,u_iu_i^\top$.
• Empirical spectral distribution (ESD): $$\mu_H=\frac1d\sum_{i=1}^d\delta_{\lambda_i}\qquad(\delta_\lambda=\text{point mass at }\lambda).$$ Every eigenvalue counts once, with equal mass $1/d$. Equal mass is the only weighting that uses nothing but $H$ itself; any other weighting needs a preferred direction.
• Goal: estimate $\mu_H$, especially its tail.
• Moments determine the measure: anything we want is $\int f\,d\mu_H$ (e.g. the CDF, with $f=\mathbf 1[\,\cdot\le\lambda]$), and polynomials approximate $f$.

How do we compute the moments efficiently?

3Probing the spectrum with a random vector

• We want $\int\lambda^k\,d\mu_H=\frac1d\operatorname{tr}(H^k)=\frac1d\sum_{j=1}^d e_j^\top H^ke_j$. Done directly, that is $k$ HVPs per term, so $d\cdot k$ HVPs with $d\approx10^8$.
• Hutchinson: $\frac1d\operatorname{tr}(A)=\mathbb E_v\,v^\top Av$ for $v\sim\mathrm{Unif}(S^{d-1})$. So use one random probe and compute $$c_k:=v^\top H^kv.$$ Since $c_{i+j}=(H^iv)^\top(H^jv)$, a single $c_k$ needs only about $k/2$ HVPs, and $m$ HVPs give all of $c_0,\dots,c_{2m-1}$.
• Which measure has these moments? Expand $v=\sum_i(v^\top u_i)\,u_i$: $$c_k=\sum_i(v^\top u_i)^2\lambda_i^k=\int\lambda^k\,d\mu_v,\qquad \boxed{\mu_v:=\sum_i(v^\top u_i)^2\,\delta_{\lambda_i}}.$$ This is the spectrum as seen from $v$: same atoms, weights $(v^\top u_i)^2$. They sum to $\|v\|^2=1$, so $\mu_v$ is a probability measure.
• Averaging over probes recovers the ESD: $\mathbb E_v\,\mu_v=\mu_H$, because $\mathbb E\,(v^\top u_i)^2=1/d$.

How do we get $\mu_v$ from its moments?

4Moment matching

Fix one probe $v$. Everything below is deterministic given $v$.

• Given: the moments $c_0,\dots,c_{2m-1}$ of $\mu_v$ ($m$ HVPs).
• Want: recover $\mu_v$ ($\approx\mu_H$) from them.
• Find $m$ atoms $\hat\mu_m=\sum_i\omega_i\delta_{\vartheta_i}$, with nodes $\vartheta_i$ (where) and weights $\omega_i$ (mass), such that $$\int\lambda^k\,d\hat\mu_m=\sum_{i=1}^m\omega_i\,\vartheta_i^{\,k}=c_k,\qquad k=0,\dots,2m-1.$$ The $k$-th moment of a spiky measure is just "weight times location to the $k$", summed over spikes. The $k=0$ equation says $\sum_i\omega_i=c_0=1$.
• That is $2m$ nonlinear equations (because of $\vartheta_i^{\,k}$) in $2m$ unknowns. The solution is unique among $m$-atom measures: the $m$-point Gauss quadrature rule of $\mu_v$. It integrates every polynomial of degree $\le2m-1$ exactly.

atoms

m = 1

How do we solve this nonlinear system?

5Prony's trick: nodes are roots of a polynomial

• Idea: look for the polynomial whose roots are the nodes. Its coefficients become the new unknowns: $$q(\lambda)=(\lambda-\vartheta_1)\cdots(\lambda-\vartheta_m)=a_0+a_1\lambda+\dots+a_{m-1}\lambda^{m-1}+\lambda^m.$$
• Why the $a_l$ are easy: $q(\vartheta_i)=0$ at every node. Weight by $\omega_i\vartheta_i^{\,j}$, sum over $i$, and use $\sum_i\omega_i\vartheta_i^{\,k}=c_k$: $$0=\sum_i\omega_i\,\vartheta_i^{\,j}\,q(\vartheta_i)=\sum_{l=0}^m a_l\sum_i\omega_i\,\vartheta_i^{\,j+l}=\sum_{l=0}^m a_l\,c_{j+l}\qquad(a_m=1).$$
• Taking $j=0,\dots,m-1$ gives $m$ linear equations, a Hankel system: $$M_0\,a=-(c_m,\dots,c_{2m-1})^\top,\qquad (M_0)_{jl}=c_{j+l}.$$
• Recipe: solve for $a$ → nodes $\vartheta_i$ = roots of $q$ → with the $\vartheta_i$ known, $\sum_i\omega_i\vartheta_i^{\,k}=c_k$ is linear in $\omega_i$ (a Vandermonde system) → weights $\omega_i$.

Exact on paper, but in practice $M_0$ is badly ill-conditioned. Why? Relate it back to $H$.

6The Krylov bridge: Prony's nodes are Ritz values

• Since $c_k=v^\top H^kv$, with the Krylov matrix $K=[\,v\ \ Hv\ \cdots\ H^{m-1}v\,]$: $$M_0=K^\top K,\qquad M_1=K^\top HK\qquad\big((M_1)_{jl}=c_{j+l+1}\big).$$
• Prony: the nodes are the solutions $\theta$ of the matrix pencil $M_1x=\theta M_0x$. (Each node contributes a rank-one term to both $M_0$ and $M_1$, and $\theta=\vartheta_i$ is exactly where $M_1-\theta M_0$ loses rank.)
• Substitute, and write $y=Kx\in\mathcal K_m=\operatorname{span}(v,Hv,\dots,H^{m-1}v)$: $$\underbrace{M_1x=\theta M_0x}_{\text{Prony}}\iff K^\top HKx=\theta\,K^\top Kx\iff\underbrace{K^\top(Hy-\theta y)=0}_{\text{Rayleigh–Ritz}}.$$
• Rayleigh–Ritz: ideally we'd want $Hy=\theta y$, but $\mathcal K_m$ usually contains no exact eigenvector. The best we can ask is that the residual be orthogonal to $\mathcal K_m$. The $m$ solutions $\theta$ are the Ritz values = nodes $\vartheta_i$ = atoms of $\hat\mu_m$.

Prony = Rayleigh–Ritz on $\mathcal K_m$, in the basis $K$. Prony never mentions $H$, but that is what it is doing.

7Prony breaks

• The columns of $K$ are $v,Hv,H^2v,\dots$: power iteration. They all align with the top eigenvector $u_1$, so $K$ is nearly rank-deficient and $\kappa(M_0)=\kappa(K)^2$ explodes.
• float64 keeps about 16 digits ($\epsilon_{\rm mach}\approx10^{-16}$). Once $\kappa(M_0)\gtrsim1/\epsilon_{\rm mach}$ there are no correct digits left. In the running example that happens at $m=8$.

depth

m = 4

The subspace $\mathcal K_m$ is fine; the basis $K$ is terrible. Lanczos: same $\mathcal K_m$, orthonormal basis.

8Lanczos: same subspace, better basis

Ritz values depend only on the subspace $\mathcal K_m$, not on the basis. Fix: keep $\mathcal K_m$, and orthonormalize its basis (Gram–Schmidt).

Same subspace, same nodes. Only the basis differs. How do we compute $Q$?

9Gram–Schmidt as a three-term recurrence

• $Q=[q_0\cdots q_{m-1}]$ is an orthonormal basis of $\mathcal K_m$, and $T_m:=Q^\top HQ$ is $H$ in that basis.
• Only $q_k$ and $q_{k-1}$ need removing: by symmetry of $H$, the components along all earlier $q_j$ are already zero. So $T_m$ is tridiagonal, with $\alpha$ on the diagonal and $\beta$ off it.

In exact arithmetic that is the whole algorithm. In floating point, orthogonality is lost as Ritz values converge and "ghost" copies of $\lambda_1$ appear; at scale one stores $Q$ and fully reorthogonalizes (see At scale).

10Reading off nodes and weights

• Write the pair $(H,v)$ in the orthonormal basis $Q$ (coordinates: multiply by $Q^\top$):
  • matrix: $H\ \to\ Q^\top HQ=T_m$ ($m\times m$). A matrix maps vectors to vectors, so it gets one $Q$ for its input and one $Q^\top$ for its output.
  • probe: $v\ \to\ Q^\top v=(q_0^\top v,\ q_1^\top v,\ \dots)=e_1$, because $q_0=v$ and every other $q_k\perp v$.
• Apply the probe-measure formula to both pairs, with $T_m=P\operatorname{diag}(\vartheta)\,P^\top$: $$\begin{array}{rcl} (H,\ v) &\longrightarrow& \mu_v=\sum_i(v^\top u_i)^2\,\delta_{\lambda_i}\\[4pt] (T_m,\ e_1) &\longrightarrow& \hat\mu_m=\sum_i(e_1^\top P_{:,i})^2\,\delta_{\vartheta_i}=\boxed{\sum_i P_{1i}^2\,\delta_{\vartheta_i}} \end{array}$$
• Nodes $\vartheta_i$: eigenvalues of $T_m$. Weights $\omega_i=P_{1i}^2$: first entries of its eigenvectors, squared. (This is the Golub–Welsch algorithm.)

$m$ HVPs → $T_m$ ($m\times m$) → one small eigendecomposition → $\hat\mu_m$.

11How deep can we see? Resolution $\approx\lambda_1/m^2$

1. Lanczos only sees polynomials of degree $<2m$. Eigenvalues on which all those polynomials agree get lumped together.
2. Change variables to angles: $\lambda_i=\tfrac{\lambda_1}{2}(1-\cos\varphi_i)$ with $\varphi_i\in[0,\pi]$. By trig identities, every polynomial of degree $\lesssim m$ becomes a cosine sum $\sum_{k\lesssim m}b_k\cos(k\varphi)$.
3. Near $\lambda=0$ (Taylor): $\cos(k\varphi)\approx1-\tfrac{(k\varphi)^2}{2}$ and $\lambda\approx\tfrac{\lambda_1}{4}\varphi^2$, i.e. $\varphi\approx2\sqrt{\lambda/\lambda_1}$. So:
   • $\lambda\ll\lambda_1/m^2\iff\varphi\ll1/m$: all $\cos(k\varphi)\approx1$, the eigenvalue looks like $\lambda=0$ ⇒ lumped.
   • $\lambda\gg\lambda_1/m^2\iff\varphi\gg1/m$: the cosines vary ⇒ resolved.

This is a guarantee in terms of eigenvalue size, not count: everything above roughly $\lambda_1/m^2$ is resolved, whatever the spectrum. With $m=1200$ that is a dynamic range of about $1.5\times10^6$.

12Error bands: Gauss–Radau

The $m$-atom solution is unique, but measures with more atoms (including the truth) match the same $2m$ moments. So the number of eigenvalues below a point $\zeta$ is uncertain. For a query point $\zeta$:

1. Force an atom at $\zeta$: grow $T_m$ by one row and column, $$\tilde T=\begin{pmatrix}T_m & \beta_{m-1}e_m\\ \beta_{m-1}e_m^\top & \gamma\end{pmatrix},$$ with $\gamma$ chosen so that $\zeta$ is an eigenvalue. ($\beta_{m-1}$ is already computed, so this costs no extra HVPs, and $\gamma$ only affects moments of order $\ge2m+1$.)
2. Eigenvalues of $\tilde T$: $m+1$ new atoms, one of them at $\zeta$, carrying the largest lump the moments allow there.
3. Band: mass strictly left of $\zeta$ $\le F(\zeta)\le$ mass left of $\zeta$ + lump.

13At scale

Settings used in A Defense of the Quadratic Model:

• Models: $d$ from 3.3M to 1B parameters.
• HVPs: over 10M tokens (100M made no difference).
• Depth: $m=1200$ ⇒ resolves $\lambda\gtrsim\lambda_1/m^2\approx\lambda_1/1.5\mathrm M$.
• Stability: full reorthogonalization (two-pass modified Gram–Schmidt), which means storing all of $Q$: about 750 GB, sharded over 16 H100s.
• Precision: fp32 (no difference from fp64; fp16 hurts).
• Error bars: Gauss–Radau bands for finite depth; concentration for the single probe (below).

AAppendix: why one random probe suffices

The probe only affects the weights, and the weights only set the index (x-position) on the index plot, $\lambda_i$ against $i$. The question is whether those x-positions are right.

• The point at height $\lambda$, whose true index is $k=\#\{i:\lambda_i\ge\lambda\}$, is drawn at $$x=d\sum_{\lambda_i\ge\lambda}(v^\top u_i)^2\qquad\big(\text{Lanczos: }d\,(\omega_1+\dots+\omega_j)\approx x\big).$$
• Write $v=g/\|g\|$ with $g\sim\mathcal N(0,I)$, and let $h=U^\top g$. By rotation invariance $h\sim\mathcal N(0,I)$ too, so $v^\top u_i=h_i/\|h\|$ and $$x=d\,\frac{h_1^2+\dots+h_k^2}{\|h\|^2}=d\,\frac{k\pm\sqrt{2k}}{d\pm\sqrt{2d}}\approx k\Big(1\pm\sqrt{2/k}\Big).$$
• A single weight $(v^\top u_i)^2=h_i^2/\|h\|^2$ is noisy ($\chi^2_1$, about ±140%). A sum of $k$ independent $h_i^2$ concentrates: errors grow like $\sqrt k$ while the sum grows like $k$.
• Head ($k$ small): the index is off (true index 1 drawn at 4). Tail ($k\ge10^3$): essentially exact. The heights $\vartheta_j$ are unaffected.