A spectrum adaptive Kernel Polynomial Method

Tyler Chen

This is a companion piece to https://doi.org/10.1063/5.0166678

Code to replicate the figures in the corresponding paper can be found on Github.

See also the spectral_density Python package.


Consider a Hermitian matrix H\mathbf{H} with eigenvalues {λi}\{\lambda_i\} and eigenvectors {ui}\{\mathbf{u}_i\}. Given a vector v\mathbf{v}, we may be interested in the Local Density of States ρ(x)=i=1dvui2δ(xλi). \rho(x) = \sum_{i=1}^{d} |\langle \mathbf{v} | \mathbf{u}_i \rangle|^2 \delta(x-\lambda_i).

We can’t efficiently compute ρ(x)\rho(x) exactly, as this is equivalent to determining all the eigenvalues. However, we can compute moments xnρ(x)dx=vHnv\int x^n \rho(x) \mathrm{d}{x} = \langle \mathbf{v}|\mathbf{H}^n|\mathbf{v}\rangle for n=0,2,,sn=0,2, \ldots, s.

The aim of the Kernel Polynomial Method (KPM) is to produce a “density” ρKPM(x)\rho_{\text{KPM}}(x) approximating ρ^(x)\hat{\rho}(x). Towards this end, let σ(x)\sigma(x) be a fixed reference density and expand ρ^(x)/σ(x)\hat{\rho}(x)/\sigma(x) as a formal polynomial series ρ^(x)σ(x)=n=0μnpn(x), \frac{\hat{\rho}(x)}{\sigma(x)} = \sum_{n=0}^{\infty} \mu_n p_n(x), where {pn}\{p_n\} are the orthonormal polynomials with respect to σ(x)\sigma(x). Using the orthonormality of the {pn}\{p_n\} with respect to σ(x)\sigma(x), we can compute μn\mu_n by μn=σ(x)ρ^(x)σ(x)pn(x)dE=ρ^(x)pn(x)dE=vpn(H)v. \mu_n = \int \sigma(x) \frac{\hat{\rho}(x)}{\sigma(x)} p_n(x) \mathrm{d}E = \int \hat{\rho}(x) p_n(x) \mathrm{d}E = \langle \mathbf{v} | p_n(\mathbf{H}) | \mathbf{v} \rangle. Thus, we see the {μn}\{\mu_n\} are the so-called (modified) moments of ρ^(x)\hat{\rho}(x) with respect to the orthogonal polynomials of σ(x)\sigma(x) and can be computed without explicit knowledge of ρ^(x)\hat{\rho}(x) using the expression above.

Computing the first ss moments naturally gives an approximation ρKPM(x)\rho_{\text{KPM}}(x) to ρ^(x)\hat{\rho}(x) defined by ρKPM(x)=σ(x)n=0sμnpn(x). \rho_{\text{KPM}}(x) = \sigma(x) \sum_{n=0}^{s} \mu_n p_n(x). When the degree of the approixmation ss\to\infty, ρKPM(x)σ(x)n=0μnpn(x)=ρ^(x). \rho_{\text{KPM}}(x) \to \sigma(x) \sum_{n=0}^{\infty} \mu_n p_n(x) = \hat{\rho}(x).

Typically, one chooses the reference density σ(x)1/1x2\sigma(x) \propto 1/\sqrt{1-x^2} to be the orthogonality measure for the Chebyshev polynomials. In this case, the modified moments vTn(H)v\langle \mathbf{v} | \mathbf{T}_n (\mathbf{H})| \mathbf{v}\rangle can be computed by computing Chebyshev polynomials in H\mathbf{H} applied to v\mathbf{v}. For such an approach to work, it is critical that H\mathbf{H} is scaled so that the spectrum is contained in [1,1][-1,1].

What’s new?

The Lanczos method for matrix function approximation (Lanczos-FA) can be used to approximate f(A)bf(\mathbf{A})\mathbf{b} (read more about this task here). For polynomials, the method is mathematically exact, and so we can use this method to compute vpn(H)v\langle \mathbf{v} | p_n(\mathbf{H}) | \mathbf{v} \rangle. In the present paper, the key observation is the expensive part of this approach (running Lanczos) is completely independent of the choice of reference density σ(x)\sigma(x). This means we can obtain lots of different KPM approxiamtions for essentially free. This avoids the need for determining how to scale H\mathbf{H} a priori and also allows us to use more “exotic” reference densities.

There is a lot of concern that the Lanczos method is unstable. Another contribution of this paper is to try and spread awareness about the way the instabilities of Lanczos actually impact methods like Lanczos-FA. For example, [1] shows that the method accurately computes (bounded) polynomial moments!


1. Knizhnerman, L.A. The Simple Lanczos Procedure: Estimates of the Error of the Gauss Quadrature Formula and Their Applications. Comput. Math. Math. Phys. 1996, 36, 1481–1492.