Randomized matrix-free quadrature: unified and uniform bounds for stochastic Lanczos quadrature and the kernel polynomial method

TL;DR

This work presents a unified theoretical framework for randomized matrix-free quadrature to approximate the spectrum and spectral sums of Hermitian matrices, focusing on Kernel Polynomial Method and Stochastic Lanczos Quadrature. It develops high-probability bounds for spectrum approximation in Wasserstein distance and extends to uniform-in-function bounds for all analytic functions on Bernstein ellipses, supported by explicit sample-degree tradeoffs. The results cover Gaussian/damped KPM settings and a broad class of analytic function tests, with detailed numerical demonstrations across sparse, smooth, and spike-structured spectra, as well as thermodynamic quantities in spin systems. The findings illuminate how choices of quadrature variant, damping, and weighting distributions affect accuracy and convergence, offering practical guidance for parameter selection and revealing the relative strengths of SLQ and KPM in diverse applications.

Abstract

We analyze randomized matrix-free quadrature algorithms for spectrum and spectral sum approximation. The algorithms studied include the kernel polynomial method and stochastic Lanczos quadrature, two widely used methods for these tasks. Our analysis of spectrum approximation unifies and simplifies several one-off analyses for these algorithms which have appeared over the past decade. In addition, we derive bounds for spectral sum approximation which guarantee that, with high probability, the algorithms are simultaneously accurate on all bounded analytic functions. Finally, we provide comprehensive and complimentary numerical examples. These examples illustrate some of the qualitative similarities and differences between the algorithms, as well as relative drawbacks and benefits to their use on different types of problems.

Figures (9)

• Figure 1: Illustration of the two main components of \ref{['alg:protoalg']}. Left: CESM $\Phi$ () and 10 independent samples of weighted CESM $\Psi$ (). Each copy is a sample from an unbiased estimator for $\Phi$ at all points $x$. Right: One sample of the weighted CESM $\Psi$ () and different approximations $\ifstrempty{}{ [\Psi]_{s}^{\circ\textup{q}} }{ [\Psi]_{s}^{\textup{{}q}} }$ to $\Psi$ based on stochastic Lanczos quadrature () and the damped kernel polynomial method (). Note that while the approximations are both induced by $\ifstrempty{}{ [\:\cdot\:]_{s}^{\circ\textup{p}} }{ [\:\cdot\:]_{s}^{\textup{{}p}} }$ for different choices of "$\circ$" they are qualitatively different: one produces a piecewise constant approximation while the other produces a continuous approximation. In different situations, one type of approximation may be preferable to the other.
• Figure 1: Errors for approximating $\int f \,\mathrm{d}\Psi = \mathbf{v}^* f(\mathbf{A}) \mathbf{v}$ when $f(x) = 1/(1+16x^2)$ for a spectrum uniformly filling $[-1,1]$ (left) and a spectrum with a gap around zero (right). Legend: Gaussian quadrature (), quadrature by interpolation (), quadrature by approximation (), Jackson's damped quadrature by approximation (), and rate $O(\rho^{-2k})$, where $\rho = (1+\sqrt{17})/4$. The behavior of the algorithms are highly dependent on the eigenvalue distribution of $\mathbf{A}$, but intuition about classical approximation theory informs our understanding. Gaussian quadrature may perform significantly better that explicit methods when the spectrum of $\mathbf{A}$ has additional structure such as gaps.
• Figure 2: Approximations to a sparse spectrum with just $12$ eigenvalues. Legend: true spectral density (), Gaussian quadrature approximation: $s=23$ (), damped quadrature by approximation: $s=500$ (). The Gaussian quadrature produces an extremely good approximation using just $12$ matrix-vector products. Even with many more matrix-vector products, damped quadrature by approximation does not have the same resolution.
• Figure 3: Wasserstien error for approximating CESM $\Phi$ supported on just 12 points. Legend: Gaussian quadrature (), Jackson's damped quadrature by approximation (), and reference for rate $O(s^{-1})$ (). The Gaussian quadrature rule converges very quickly (in just 12 steps) to the true weigthed CESM, while the damped quadrature by approximation converges at a rate of $O(s^{-1})$.
• Figure 4: Approximations to a "smooth" spectrum using quadrature by approximation with various choices of $\mu$. Legend: $\mu = \mu_{a_1,b_2}^U$ (), $\mu = \frac{1}{2} \mu_{a_1,b_1}^U + \frac{1}{2} \mu_{a_2,b_2}^U$(). A priori knowledge about properties of the spectrum allow better choices of parameters such as $\mu$.

Theorems & Definitions (31)

• Definition 2.1: Quadrature by interpolation; $\circ = \textrm{i}$
• Definition 2.2: Gaussian quadrature; $\circ = \textrm{g}$
• Definition 2.3: Quadrature by approximation; $\circ = \textrm{a}$
• Definition 2.4: damped quadrature by approximation; $\circ = \textrm{d}$
• Definition 3.1
• Definition 3.2
• Definition 3.3
• Definition 3.4
• Lemma 3.5
• Proof 1