This is an English translation of the Chinese original, first published on 14 July 2025 (use the 中文 button in the menu to read it). The translation was prepared with the assistance of Claude Opus 5.5 (Anthropic).
Kernel Polynomial Method (KPM)
The kernel polynomial method (KPM) uses a numerical Chebyshev polynomial expansion to reduce computational cost. Computing a spectrum the traditional way usually requires \(O(N^3)\) diagonalization, whereas KPM turns the problem into multiplying the matrix with a small set of vectors, which greatly reduces the cost. The Lanczos algorithm was developed in the same period; Lanczos is suited to finding eigenstates, while KPM is better suited to computing spectral properties and correlation functions [1].
KPM comes from expanding functions in orthogonal polynomials. Given a weight function \(w(x)\) and an interval \([a,b]\), the inner product of two integrable functions \(f,g:[a,b]\to\mathbb{R}\) is
\[\begin{equation} \langle f|g \rangle = \int_a^b w(x) f(x) g(x) dx \end{equation}\]
Gram–Schmidt orthogonalization gives orthogonal polynomials \(p_n(x)\) satisfying
\[\begin{equation} \langle p_n|p_m \rangle = \delta_{nm}/h_n \end{equation}\] where \(h_n = ⟨p_n|p_n⟩\). Any function can be expanded as \[\begin{equation} f(x) = \sum_{n=0}^\infty \alpha_n p_n(x), \quad \alpha_n = \frac{1}{h_n}\langle p_n|f \rangle \end{equation}\] In principle any set of orthogonal polynomials can be used to expand a function and hence to run KPM. Through history and experience, however, Chebyshev polynomials of the first and second kind have come to be regarded as the best choice, thanks to their excellent convergence and their similarity to the Fourier transform.
Chebyshev polynomials
Chebyshev polynomials are defined on the interval \([a,b]=[-1,1]\). The weight function of the first-kind Chebyshev polynomials \(T_n(x)\) is \(w(x)=\frac{1}{\pi\sqrt{1-x^2}}\), and that of the second kind, \(U_n(x)\), is \(w(x)=\pi\sqrt{1-x^2}\). One can check that
\[\begin{equation} \begin{split} T_n(x) &= \cos\left[n \cos^{-1}x\right] \\ U_n(x) &= \frac{\sin\left[(n+1)\cos^{-1}x\right]}{\sin\left[\cos^{-1}x\right]} \end{split} \end{equation}\]
Both kinds of Chebyshev polynomials obey recurrence relations. For the first kind,
\[\begin{equation} \begin{split} T_0(x) &= 1, \\ T_{-1}(x) &= T_1(x) = x, \\ T_{n+1}(x) &= 2xT_n(x)-T_{n-1}(x) \end{split} \end{equation}\]
and for the second kind,
\[\begin{equation} \begin{split} U_0(x) &= 1, \\ U_{-1}(x) &= 0, \\ U_{n+1}(x) &= 2xU_n(x)-U_{n-1}(x) \end{split} \end{equation}\]
There are also product relations:
\[\begin{equation} \begin{split} 2T_m(x)T_n(x) &= T_{m+n}(x)+T_{m-n}(x), \\ 2(x^2-1)U_{m-1}(x)U_{n-1}(x) &= T_{m+n}(x)-T_{m-n}(x) \end{split} \end{equation}\]
The two weight functions give two definitions of the inner product:
\[\begin{equation} \begin{split} \langle f|g \rangle_1 &= \int_{-1}^{1} \frac{f(x)g(x)}{\pi\sqrt{1-x^2}} dx, \\ \langle f|g \rangle_2 &= \int_{-1}^{1} \pi\sqrt{1-x^2} f(x)g(x) dx \end{split} \end{equation}\]
The corresponding orthogonality relations are
\[\begin{equation} \begin{split} \langle T_n|T_m \rangle_1 &= \frac{1+\delta_{n,0}}{2}\delta_{n,m} \\ \langle U_n|U_m \rangle_2 &= \frac{\pi^2}{2}\delta_{n,m} \end{split} \end{equation}\]
Define \(\phi_n(x)=\frac{T_n(x)}{\pi\sqrt{1-x^2}}\). Then, under the second inner product,
\[\begin{equation} \langle \phi_n|\phi_m \rangle_2 = \frac{1+\delta_{n,0}}{2}\delta_{nm} \end{equation}\]
One can check that any function can be expanded as
\[\begin{equation} \begin{split} f(x) &= \sum_{n=0}^{\infty} \frac{\langle f|\phi_n \rangle_2}{\langle \phi_n|\phi_n \rangle_2} \phi_n(x) \\ &= \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 + 2\sum_{n=1}^{\infty} \mu_n T_n(x)\right) \end{split} \end{equation}\]
where the coefficients \(\mu_n\) are called moments and are defined as
\[\begin{equation} \mu_n = \langle f|\phi_n \rangle_2 = \int_{-1}^{1} f(x) T_n(x) dx \end{equation}\]
Kernel polynomials and the Gibbs phenomenon
In practice the expansion has to be truncated, and truncating it directly causes a severe Gibbs phenomenon. To fix this, one usually introduces damping kernel coefficients, \(\mu_n \to g_n \mu_n\), to modify the expansion coefficients:
\[\begin{equation} \begin{split} f(x) &= \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 + 2\sum_{n=1}^{\infty} \mu_n T_n(x)\right) \\ &\to \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 g_0 + 2\sum_{n=1}^{N-1} \mu_n g_n T_n(x)\right) \end{split} \end{equation}\]
A commonly used kernel is the Jackson kernel:
\[\begin{equation} g_n^{\text{J}} = \frac{1}{N+1}\left[(N-n+1)\cos\frac{\pi n}{N+1} + \sin\frac{\pi n}{N+1}\cot\frac{\pi}{N+1}\right] \end{equation}\]
For Green's function calculations the Lorentz kernel is common:
\[\begin{equation} g_n^L = \frac{\sinh\left[\lambda\left(1-\frac{n}{N}\right)\right]}{\sinh \lambda} \end{equation}\]
where the parameter \(\lambda \in \mathbb{R}\) is usually taken between \(3\) and \(5\).
Computing the density of states
Rescaling the interval
The Chebyshev expansion is defined on \([-1,1]\), so for matrix computations the eigenvalues \(\{E_k\}\) must also lie in \([-1,1]\). This can be arranged with a linear transformation:
\[\begin{equation} \begin{split} \tilde{H} &= \frac{H-b}{a} \\ \tilde{E} &= \frac{E-b}{a} \end{split} \end{equation}\]
where \(a\) and \(b\) are fixed by the largest and smallest eigenvalues of the Hamiltonian, \(E_{\text{max}}\) and \(E_{\text{min}}\):
\[\begin{equation} a = \frac{E_{\text{max}}-E_{\text{min}}}{2-\epsilon} \\ b = \frac{E_{\text{max}}+E_{\text{min}}}{2} \end{equation}\]
The parameter \(\epsilon\) keeps the edges of the spectrum from landing exactly on, or beyond, \([-1,1]\); it is usually taken as \(0.01\).
The largest and smallest eigenvalues can be obtained with the Lanczos
algorithm (for example MATLAB's eigs):
1 | E_max = eigs(H_sparse, 1, 'largestreal') |
Discrete cosine Fourier transform
After truncating the expansion of \(f(x)\) and correcting it with a kernel, we have, at a given point \(x\), the expansion
\[\begin{equation} f(x) = \frac{1}{\pi\sqrt{1-x^2}}\left(\mu_0 g_0 + 2\sum_{n=1}^{N-1} \mu_n g_n T_n(x)\right) \end{equation}\]
The properties of Chebyshev polynomials let us cut the computing time further. Define the kernel-corrected moments as
\[\begin{equation} \tilde{\mu}_n = \mu_n g_n \end{equation}\]
The expansion can then be written as
\[\begin{equation} f(x) = \frac{1}{\pi\sqrt{1-x^2}}\left(\tilde{\mu}_0 + 2\sum_{n=1}^{N-1} \tilde{\mu}_n T_n(x)\right) \end{equation}\]
If the sampling points \(x\) are chosen as the special points \(x_k\) satisfying
\[\begin{equation} x_k = \cos\frac{\pi(k+1/2)}{\tilde{N}},\quad k=0,\ldots,(\tilde{N}-1) \end{equation}\]
then
\[\begin{equation} T_n(x_k) = \cos\left[n\cos^{-1}x_k\right] = \cos\left[\frac{n\pi(k+1/2)}{\tilde{N}}\right] \end{equation}\]
These points are called Chebyshev nodes [2]. \(\tilde{N}\) is usually taken as \(2N\). At these points, \(f(x_k)\) can be obtained by applying the type-III discrete cosine transform (DCT-III) to \(\tilde{\mu}_n\):
\[\begin{equation} y[k] = \frac{1}{2}x_0 + \sum_{n=1}^{N-1} x_n \cos\left[\frac{n\pi(k+1/2)}{N}\right] \end{equation}\]
Since the type-III cosine transform can be accelerated with the fast Fourier transform (FFT), this is faster than direct multiplication. Example MATLAB code:
1 | omega_dct = cos(pi * ((0:N_points-1)+0.5)/N_points); |
Note: the second line divides the first moment by \(\sqrt{2}\). This is because MATLAB's
dct definition differs from the type-III cosine transform
above by a factor of \(\sqrt{2}\)
between the first term and the remaining terms. Without this change, the
cosine-transform coefficients converge to zero more slowly and give a
wrong density of states.
Computing the moments
The KPM expansion can be applied to operators, for example the Hamiltonian \(\hat{H}\). When computing with KPM, two forms of moments are of interest.
The first is an expectation between two states: \[\begin{equation} \mu_n = \langle \beta | T_n(\tilde{H}) | \alpha \rangle \end{equation}\]
The second is the trace of an operator \(\hat{A}\) with the polynomial: \[\begin{equation} \mu_n = \mathrm{Tr}[A T_n(\tilde{H})] \end{equation}\]
In the first case the moments can be computed directly from the recurrence. Let \(|\alpha_n\rangle = T_n(H)|\alpha\rangle\); then
\[\begin{equation} \begin{split} |\alpha_0\rangle &= |\alpha\rangle \\ |\alpha_1\rangle &= \tilde{H}|\alpha_0\rangle \\ |\alpha_{n+1}\rangle &= 2\tilde{H}|\alpha_n\rangle - |\alpha_{n-1}\rangle \end{split} \end{equation}\]
and the moments follow as \(\mu_n = \langle \beta | \alpha_n \rangle\). The whole procedure involves only matrix multiplication, and sparse-matrix methods can make large calculations efficient. If \(|\beta\rangle = |\alpha\rangle\), it simplifies further. Using the product relation, one can check that for the even moments
\[\begin{equation} \begin{split} 2\langle \alpha_n | \alpha_n \rangle - \mu_0 &= 2\langle \alpha | T_n(\tilde{H})^\dagger T_n(\tilde{H}) | \alpha \rangle - \mu_0 \\ &= 2\langle \alpha | T_n(\tilde{H}) T_n(\tilde{H}) | \alpha \rangle - \mu_0 \\ &= \langle \alpha | \left(T_0(\tilde{H}) + T_{2n}(\tilde{H})\right) | \alpha \rangle - \mu_0 \\ &= \langle \alpha | \alpha_0 \rangle + \langle \alpha | \alpha_{2n} \rangle - \mu_0 \\ &= \mu_{2n} \end{split} \end{equation}\]
and for the odd moments
\[\begin{equation} \begin{split} 2\langle \alpha_{n+1} | \alpha_n \rangle - \mu_1 &= 2\langle \alpha | T_{n+1}(\tilde{H})^\dagger T_n(\tilde{H}) | \alpha \rangle - \mu_1 \\ &= 2\langle \alpha | T_{n+1}(\tilde{H}) T_n(\tilde{H}) | \alpha \rangle - \mu_1 \\ &= \langle \alpha | \left(T_1(\tilde{H}) + T_{2n+1}(\tilde{H})\right) | \alpha \rangle - \mu_1 \\ &= \langle \alpha | \alpha_1 \rangle + \langle \alpha | \alpha_{2n+1} \rangle - \mu_1 \\ &= \mu_{2n+1} \end{split} \end{equation}\]
So when \(|\alpha\rangle = |\beta\rangle\), the number of moment computations can be halved:
\[\begin{equation} \begin{split} \mu_{2n} &= 2\langle \alpha_n | \alpha_n \rangle - \mu_0 \\ \mu_{2n+1} &= 2\langle \alpha_{n+1} | \alpha_n \rangle - \mu_1 \end{split} \end{equation}\]
For the second case, take the density of states as an example. With eigenstates \(\tilde{H}|k\rangle = E_k|k\rangle\), the density of states is defined as
\[\begin{equation} \rho(E) = \sum_k \delta(E-E_k) = \mathrm{Tr}\left[\delta(E-\tilde{H})\right] \end{equation}\]
The \(\delta\) function can be expanded in Chebyshev polynomials, giving the moments
\[\begin{equation} \begin{split} \mu_n &= \int_{-1}^{1} \rho(E) T_n(E) dE \\ &= \int_{-1}^{1} \sum_k \delta(E-E_k) T_n(E) dE \\ &= \sum_k T_n(E_k) \end{split} \end{equation}\]
Using \(T_n(H)|k\rangle = T_n(E_k)|k\rangle\), we get
\[\begin{equation} \mu_n = \sum_k \langle k | T_n(H) | k \rangle = \mathrm{Tr}[T_n(H)] \end{equation}\]
Because Chebyshev polynomials obey a recurrence relation,
\[\begin{equation} \begin{split} |\alpha_n(k)\rangle &= T_n(H)|k\rangle \\ |\alpha_{n+1}(k)\rangle &= 2H|\alpha_n(k)\rangle - |\alpha_{n-1}(k)\rangle \end{split} \end{equation}\]
and the moments become
\[\begin{equation} \mu_n = \sum_k \langle \alpha_0(k) | \alpha_n(k) \rangle = \sum_k \langle k | T_n(H) | k \rangle \end{equation}\]
Note that the earlier formulas turned the density of states into a calculation of expansion coefficients, and the equation above shows that each expansion coefficient only requires the inner product of two vectors. Since each order of \(|\alpha_n(k)\rangle\) is really just a matrix multiplication, these steps let us skip the costly diagonalization altogether and obtain the expansion coefficients of the density of states from matrix multiplications alone.
In fact, computing the trace does not require summing over all eigenstates. Averaging over a small number of random normalized vectors \(|r\rangle\) already gives a good estimate of \(\mu_n\):
\[\begin{equation} \mu_n = \mathrm{Tr}[A T_n(H)] \approx \frac{1}{R} \sum_{r=1}^R \langle r | A T_n(H) | r \rangle \end{equation}\]
where \(|r\rangle = \frac{1}{\sqrt{D}} \sum_{i=1}^D \xi_{rk} |k\rangle\) and \(\xi_{rk}\) are complex random normalized vectors satisfying
\[\begin{equation} \begin{split} \mathbb{E}[\xi_{rk}] &= 0 \\ \mathbb{E}[\xi_{rk}\xi_{r'l}] &= 0 \\ \mathbb{E}[\xi^*_{rk}\xi_{r'l}] &= \delta_{rr'}\delta_{kl} \end{split} \end{equation}\]
If \(\xi_{rk}=e^{i\phi}\) with \(\phi\in[0,2\pi]\) drawn from a uniform or Gaussian distribution, the results are quite accurate.