Spectra from Moments

Five things you half-know are secretly the same object: the density of states, the Lanczos algorithm, Gaussian quadrature, Padé approximants, and the continued-fraction form of a Green's function. The thread that ties them is one idea — a spectral function is a probability measure, its moments are matrix–vector products, and Lanczos turns those moments into a quadrature rule for the spectrum. Once you see it, "compute the spectrum of a huge operator" stops being a dozen techniques and becomes one move.

1. A spectral function is a probability distribution

Take a large Hermitian operator $H$ with eigenvalues $E_n$ and eigenvectors $|n⟩$. The eigenvectors form a basis, so any unit vector $|v⟩$ can be expanded in them, $|v⟩={\sum }_n{c}_n|n⟩$ with coefficients $c_n=⟨n|v⟩$. The squared coefficients $|c_n{|}^2$, placed at their eigenvalues $E_n$, form a distribution over the spectrum — the spectral function:

Each spike sits at an eigenvalue $E_n$ with height $|c_n{|}^2=|⟨n|v⟩{|}^2$, the squared magnitude of that expansion coefficient. Because $|v⟩$ is normalized the heights are nonnegative and sum to one, ${\sum }_n|c_n{|}^2=1$ — so $S(E)$ is genuinely a probability distribution over the spectrum. The same object is called the local density of states (solid state), the strength function $|⟨n|\hat{O}|0⟩{|}^2$ (nuclear, with $|v⟩=\hat{O}|0⟩$), and the absorption spectrum (TDDFT). One object, many names — the difference is only which $|v⟩$ you choose.

What does it mean to put H inside a delta?

The right-hand side is an ordinary sum of spikes, but the left-hand side has an operator inside a delta function. That is shorthand from the functional calculus. For any ordinary function, the operator version is defined to keep the eigenvectors and apply the function to the eigenvalues:

This is the same rule behind every operator function you already use — time evolution $e^{-iHt}$, the thermal weight $e^{-\beta H}$, the inverse $H^{-1}$: apply the scalar function to the eigenvalues, leave the eigenvectors alone. Take the function to be the delta and you get

an operator that depends on the real dial $E$. It is zero almost everywhere; the instant $E$ lands on an eigenvalue it becomes the projector $|n⟩⟨n|$ onto that eigenstate. Sweeping the dial across all energies simply slices the identity operator up by eigenvalue:

Sandwich it between $⟨v|$ and $|v⟩$ and the projectors turn into the squared coefficients, which is exactly why the two sides of the spectral function agree:

And the bridge to your resolvent page: on the real axis this operator-delta is just the imaginary part of the Green's function (Sokhotski–Plemelj). The spectral spikes are the imaginary edge of the smooth resolvent:

2. The moments are free

A distribution has moments. The $k$-th moment of $S$ is its average of $E^k$; let the delta functions do the integral and something nice happens:

The $k$-th moment of the spectrum is just apply $H$ to $|v⟩$ a total of $k$ times, then dot with $|v⟩$ — no eigenvalues, no diagonalization, only matrix–vector products. The moments peel off the distribution's shape one feature at a time: ${\mu }_0=1$ (it is normalized), ${\mu }_1=⟨H⟩$ is the mean (its center of mass), and the first two together give the variance — the width:

3. Lanczos builds the orthogonal polynomials of S

Computing $S$ from its moments is the classical moment problem, and the tool that solves it is orthogonal polynomials. Lanczos is exactly that tool in disguise. Run Gram–Schmidt on the Krylov sequence $|v⟩,H|v⟩,H^2|v⟩,\dots$. Because $H$ is Hermitian you get a three-term recurrence, and the recurrence coefficients fill a tridiagonal Jacobi matrix:

The $j$-th Lanczos vector is $q_j=p_j(H)|v⟩$, where $p_j$ is the degree-$j$ orthonormal polynomial with respect to the measure $S$. The Lanczos $\alpha$'s and $\beta$'s are nothing but the recurrence coefficients of those polynomials. Lanczos never "knew" it was doing classical orthogonal-polynomial theory on the spectral measure — it just was.

4. It is Gaussian quadrature

Diagonalize the small $N\times N$ matrix $T_N$: eigenvalues ${\theta }_j$, eigenvectors $u_j$. Then

and this is not an approximation — it is the $N$-point Gauss quadrature rule for the measure $S$. By the defining property of Gauss quadrature, it is exact for every polynomial up to degree $2N-1$, which means it reproduces the first $2N$ moments ${\mu }_0,\dots ,{\mu }_{2N-1}$ exactly. So:

$N$ Lanczos steps — about $N$ matrix–vector products — give the unique $N$-pole spectrum that matches the true spectrum through $2N$ moments. Nothing with that many parameters does better. It also explains why Lanczos pins down the extremal eigenvalues first: the extremes dominate the high moments, and quadrature spends its accuracy exactly where the moments are largest.

5. The resolvent is a continued fraction

Your Green's function $G(z)=⟨v|(z-H{)}^{-1}|v⟩$ equals the top-left element of $(z-T{)}^{-1}$, and the $(1,1)$ element of a tridiagonal inverse is a continued fraction in the Lanczos coefficients:

Truncating at level $N$ is simultaneously the Gauss quadrature above, the $[N-1/N]$ Padé approximant of $G(z)$, and the recursion method (Haydock) for the density of states — the same truncation seen three ways. The kernel polynomial method is the cousin where you expand in Chebyshev moments instead of these power moments: same philosophy, different orthogonal family.

Why this is the spine

Every spectral task is now one move — project onto a Krylov space, read the Jacobi matrix:

• A few eigenvalues → Lanczos directly (Lanczos iteration).
• A whole density of states / strength function → the quadrature, or the continued fraction (recursion method, KPM).
• A response function $G(z)$ off the real axis → the continued fraction (the resolvent; the Liouville–Lanczos and Casida machinery).
• Resummation of a divergent series → Padé / Borel on that same continued fraction (the perturbation-methods pages).