The Kernel Polynomial Method

The density of states $\rho (E)=\frac{1}{D}{\sum }_k\delta (E-E_k)$ is built from every eigenvalue of $H$. For a lattice model with millions or billions of sites, computing all $D$ of them is out of the question — and pointless, because the answer you want is a smooth curve, not a list of $10^9$ numbers.

So don't ask for the eigenvalues. Ask for the shape of their distribution directly. KPM computes $\rho (E)$ to whatever resolution you pay for, using only the one thing a sparse matrix is good at: multiplying a vector. No eigensolver, no dense storage.

Rescale $H$ so its whole spectrum lands inside $[-1,1]$ — call it $\tilde{H}=(H-b)/a$. Now $\rho$ is a function on that interval, and any such function has a Chebyshev expansion:

The expansion coefficients ${\mu }_n$ are the Chebyshev moments of the density, ${\mu }_n=\int T_n(\tilde{E})\rho (\tilde{E})d\tilde{E}$. Truncating at $N$ moments gives a resolution that sharpens like $1/N$ — more moments, finer features. The whole problem reduces to: get the moments.

Here is the key identity. Because $T_n$ is a polynomial, the moment is a trace of a matrix polynomial:

And you never form $T_n(\tilde{H})$. The Chebyshev polynomials obey a three-term recurrence, so apply it to a vector: $T_0(\tilde{H})|v⟩=|v⟩$, $T_1(\tilde{H})|v⟩=\tilde{H}|v⟩$, and $T_{n+1}(\tilde{H})|v⟩=2\tilde{H}{T}_n(\tilde{H})|v⟩-T_{n-1}(\tilde{H})|v⟩$. Each step is one sparse matrix–vector product. If that recurrence looks familiar, it should — it is the same three-term structure that makes Lanczos cheap. KPM and Lanczos are two bases for the same moment problem.

One problem remains: $\mathrm{Tr}{T}_n(\tilde{H})$ still sums a diagonal over all $D$ basis states. The escape is stochastic trace estimation. For a random vector $|r⟩$ whose components are independent with unit variance,

A handful of random probes $R$ — often just 10–30, independent of system size — gives every moment at once: run the recurrence on each $|r⟩$ and accumulate the inner products. This is the same Hutchinson trick that estimates $\log \det$ and $\mathrm{Tr}f(A)$ elsewhere, and it is why KPM scales to operators you could never even store as a dense matrix.

Cut the Chebyshev series off at $N$ terms and you get Gibbs oscillations — the reconstructed density rings around sharp features and can dip below zero, which is nonsense for a density. The fix, and the method's namesake, is to multiply each moment by a damping factor: ${\mu }_n\to g_n{\mu }_n$.

That sequence $\{g_n\}$ is the kernel. The Jackson kernel is the standard choice: it convolves the exact density with a near-Gaussian of width $\sim 1/N$, killing the ringing while keeping the result everywhere non-negative. Truncation becomes honest broadening instead of artefacts. Toggle it in the demo below and watch the wiggles vanish — that one trick is the difference between a usable spectrum and a noisy one.