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

1. Introduction.

• The method first estimates the cumulative empirical spectral measure using weighted measures from random vectors, then approximates each weighted measure with polynomial quadrature.A weighted measure from a random unit vector is an unbiased estimator of the cumulative empirical spectral measure; averaging independent copies reduces variance.Source: §1 · pp. 1–2
• Different polynomial-operator choices yield qualitatively different approximations: one is piecewise constant, while another is continuous.The paper notes that either type may be preferable in different situations.Source: §1 · p. 2
• Moments of the weighted measure through degree 2k can be obtained from the Krylov subspace K_k(A, v).This subspace is spanned by v, Av, through A^k v, equivalently by p(A)v for polynomials of degree at most k.Source: §1 · p. 2
• The paper analyzes randomized matrix-free quadrature for spectrum and spectral-sum approximation, with theoretical bounds and experiments highlighting KPM and SLQ.The authors aim to compare their theoretical and practical behavior across both tasks.Source: §1 · p. 3

1.1. Notation.

• The notation section defines ∥f∥_S as the supremum of |f(x)| over x in S.It also uses bold uppercase letters for matrices, bold lowercase letters for vectors, and B∗ for the conjugate transpose.Source: §1.1 · p. 3

2.1. Past theoretical analysis.

• Earlier spectral-sum bounds controlled error for a fixed function, while prior spectrum bounds gave Wasserstein accuracy with degree O(ε^-1) and could use constant samples above the stated scale.The paper identifies bounded analytic functions as a less-studied setting and provides explicit simultaneous bounds for that class.Source: §2.1 · pp. 3–4, §2.2 · p. 4

2.2. Our contributions.

• The paper unifies spectrum-approximation analyses and derives simultaneous bounds for bounded analytic functions, alongside experiments on SLQ, KPM, and non-Chebyshev polynomial families.The authors describe the analytic-function bounds as the first of their type for trace-estimation algorithms and report that non-Chebyshev choices may offer practical benefits.Source: §2.2 · p. 4

2.3. Other algorithms.

• Hutch++ combines a low-rank approximation of f(A) with quadratic trace estimation on the remainder, achieving O(ε^-1) matrix-vector products for positive definite matrices.The cited analysis gives relative trace error 1 ± ε.Source: §2.3 · p. 4

2.4. Quadrature approximations for weighted spectral measures.

• The paper presents quadrature by interpolation, Gaussian quadrature (used by SLQ), approximation, and damped approximation (used by KPM); damping is used to ensure positivity.Gaussian quadrature integrates polynomials through degree 2s−1 exactly.Source: §2.4 · pp. 4–5
• Jackson damping gives a nonnegative kernel and a popular damped KPM variant; the paper estimates its polynomial-approximation convergence rate.Damping replaces each polynomial coefficient p_k with ρ_k p_k; the resulting approximation is weakly increasing under the stated conditions.Source: §2.4 · pp. 5–6

3. Main bounds.

• The bounds cover several instances of the protoalgorithm, prioritize simple statements, and do not analyze the ◦ = i case.The authors note that the analysis indicates how similar bounds for that case could be derived.Source: §3 · p. 6

3.1. Spectrum approximation.

• The spectrum analysis unifies KPM and SLQ bounds, slightly improves KPM constants, and trades SLQ failure-probability dependence on n for dependence on s, with slightly worse accuracy-related constants.The authors identify the unified analysis as the main advantage over the separate prior analyses.Source: §3.1 · p. 8
• For ε ≫ n^-1/2 as n → ∞, Theorem 3.11 says one test vector suffices to approximate the cumulative empirical spectral measure in Wasserstein distance to accuracy ε.Source: §3.1 · p. 10

3.2. Uniform spectral sum approximation.

• Theorem 3.15 gives high-probability spectral-sum bounds simultaneously for all functions in the specified bounded analytic class, although its sample requirement depends on s.This differs from triangle-inequality bounds that do not depend on s but apply only to a single function.Source: §3.2 · p. 11
• The uniform-bound approach does not extend directly to Jackson-damped KPM: its moments do not match, and its polynomial approximations converge only algebraically, precluding logarithmic dependence of s on ε.Source: §3.2 · p. 11

4.1. Comparison with classical quadrature.

• With equal matrix-vector products on a uniformly filled interval, interpolation initially converges twice as quickly as Gaussian quadrature, then matches its O(ρ^-2k) rate.This comparison uses the Runge function f(x) = 1/(1 + 16x^2).Source: §4.1 · pp. 13–14
• When the spectrum has a gap, Gaussian quadrature converges significantly faster, while quadrature by interpolation converges at its limiting rate.More broadly, Gaussian quadrature can outperform bounds based on polynomial approximation over the full spectral interval when spectra have gaps or outlying eigenvalues.Source: §4.1 · p. 14

4.2. Approximating sparse spectra.

• For the 12-point sparse spectrum, Gaussian quadrature is exact at s = 23 (k = 12), whereas damped quadrature by approximation converges at O(s^-1).The damped method has lower resolution even with s = 500 matrix-vector products; Gaussian quadrature can exploit sparsity, while fixed-distribution methods cannot.Source: §4.2 · pp. 14–15, Figure 4 · p. 15

4.3. Approximating “smooth” densities.

• For a smooth spectrum, choosing µ to reflect its known support and square-root edge behavior yields an approximation with the correct endpoint behavior and fewer interior oscillations.The tailored choice keeps the Radon–Nikodym derivative of the limiting density with respect to µ of order 1 at the spectral edges.Source: §4.3 · pp. 16–17

4.3.1. Smoothing by convolution.

• Convolution smoothing has Wasserstein error at most σ for a median-zero kernel with standard deviation σ, but too little smoothing leaves discreteness and too much blurs spectral resolution.A single smoothing scale can be problematic when different spectral regions need different amounts; the paper notes that alternative metrics and regularization merit further study.Source: §4.3.1 · pp. 17–18, Figure 6 · p. 18

4.3.2. Handling spikes.

• Explicitly accounting for a spectral spike produces a far better approximation to the continuous density, even when the reference distribution assigns the spike twice its true mass.An incorrectly estimated spike location can instead cause large oscillations.Source: §4.3.2 · pp. 19–20, Figure 7 · p. 19

4.4. Thermodyanmics of spin systems.

• Because matrix size grows exponentially with spin-site count and diagonalization is not always applicable, the experiments use Algorithm 1.1-based approaches.The passage cites Lanczos-based examples [SRS20] and Chebyshev-based examples [Sch+21].Source: §4.4 · p. 21
• For a 12-site, spin-1/2 Heisenberg ring, the experiments use k = 50 and m = 300 to approximate heat capacity across temperatures with three quadrature variants.The variants are Gaussian quadrature, quadrature by interpolation, and damped quadrature by interpolation; the same random vectors are used for all four traces.Source: §4.4 · p. 21
• Quadrature by approximation produces a nonphysical ghost dip in heat capacity, whereas Jackson damping keeps heat capacity positive at all temperatures.Positivity follows from the Cauchy–Schwarz inequality when the CESM approximation is non-decreasing; the authors note a ghost peak may be harder to identify.Source: §4.4 · p. 22
• Gaussian quadrature is simultaneously accurate for infinitely many partition-function functions, while quadrature by approximation shows larger error, especially at large β.Increasing m and s decreases error; the same qualitative trend holds for quadrature by approximation, without theoretical guarantees.Source: Figure 9 · p. 21, §4.4 · p. 22

• The paper provides a unified analysis of a general class of algorithms for spectrum and spectral sum approximation.Source: §5 · p. 22
• The bounds separately analyze approximating the CESM with a random weighted CESM and then approximating it with a polynomial quadrature rule.Source: §5 · p. 22
• Lanczos-based methods often outperform bounds based on best approximation over a fixed interval, motivating spectrum-dependent or a posteriori bounds for stopping criteria.Source: §5 · p. 22
• Illustrative numerical experiments highlight matrix-free quadrature peculiarities and the value of choosing parameters using application-specific knowledge.Source: §5 · p. 22

SM1. Introduction

• The supplement discusses selecting KPM parameters after the main computation, implementation choices, finite-precision behavior, and why Lanczos loss of orthogonality does not stop SLQ convergence.It also provides additional numerical experiments.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 26

SM2.1. Computing modified moments directly

• Direct polynomial recurrences compute modified moments with k matrix-vector products and O(n) storage; Chebyshev identities also enable low-memory extraction through degree 2k.Chebyshev polynomials can be numerically unstable if the spectrum extends beyond [-1, 1], so the interval should contain the spectrum.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · pp. 28–29

SM2.2. Lanczos algorithm

• When µ = Ψ, the Stieltjes procedure is mathematically equivalent to Lanczos, which computes the Jacobi matrix’s diagonal and off-diagonal coefficients in k matrix-vector products.The algorithm optionally reorthogonalizes its vectors.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 29

SM2.3. Connection coefficients to compute more modified moments

• Connection coefficients transform modified moments between distributions as m = C^Tn; combined with k Lanczos iterations, they can recover moments through degree 2k when the reference Jacobi matrix is known.The coefficient computation is generally less costly than products with A.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · pp. 30–31
• Recovering quadrature nodes and weights or a Jacobi matrix from modified moments requires care because these maps may be severely ill-conditioned in finite precision.The tridiagonal matrix from Lanczos can be obtained from the modified moments.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 31

SM3.1. Quadrature by interpolation; SM3.2. Gaussian quadrature and stochastic Lanczos quadrature

• Gaussian quadrature integrates polynomials through degree 2k − 1 exactly, while degree-s interpolation quadrature integrates polynomials through degree s exactly.Using orthogonal polynomials for interpolation gives an orthogonal-column system that can be inverted through products with its transpose and a diagonal matrix.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · pp. 32–33

SM3.3.1. Evaluating spectral sums and the relation to quadrature by interpolation

• For spectral sums, the approximation's distribution function need not be computed explicitly: its polynomial projections can be known analytically or computed as a scalar problem independent of matrix size.When the projections are computed by quadrature, interpolation is equivalent to approximation with d = s + 1.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 34

SM3.4. Positivity by damping and the kernel polynomial method

• Damping preserves unit mass when ρ0 = 1 and makes approximations weakly increasing when the damping kernel is nonnegative and f is nonnegative.The standard kernel polynomial method uses Jackson-coefficient damping.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 35

SM4.1. Computational costs

• The core moment and Lanczos algorithms cost O(k(Tmv + n)) and can use O(n) storage; full reorthogonalization raises cost to O(k(Tmv + kn)).Storing the Krylov basis can also enable delayed inner products that reduce global reductions, potentially speeding computation on highly parallel machines.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · pp. 35–36

SM4.2. Stability in finite precision arithmetic

• Chebyshev-based modified moments are expected to remain close to exact when the spectrum lies in [a, b], but can become exponentially unstable when it does not.The instability arises because the associated Chebyshev polynomials of the second kind can grow exponentially outside [a, b].Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 37

SM4.2.1. The Lanczos algorithm

• For k not too large, Lanczos Gaussian quadrature computes modified moments accurately, supporting stable Chebyshev-distribution moments despite finite-precision perturbations.The guarantee concerns moments through degree 2k − 1; the associated Chebyshev polynomials are scaled to an interval slightly larger than the spectrum.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · pp. 37–38

SM4.2. Stability in finite precision arithmetic

• Lanczos-based Gaussian quadrature is expected to perform nearly as well as, and often better than, properly scaled Chebyshev approaches, without requiring prior spectral bounds.Chebyshev approaches require an interval [a, b] containing the spectrum to avoid exponential instabilities.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 38

SM6.1. Finite precision convergence of SLQ

• Without reorthogonalization, Gaussian quadrature converges more slowly but can still converge and outperform the other tested approximation methods.The comparison includes Gaussian quadrature with and without reorthogonalization, quadrature by approximation, and quadrature by interpolation.Source: Figure SM1 · p. 40
• Experiments used diagonal matrices scaled to ∥A∥2 = 1 and uniform starting vectors, including an exponential-spacing model and the California matrix [DH11].The model used n = 300, κ = 10^3, ρ = 0.85, and f(x) = 1/x; the California matrix used n = 9664 eigenvalues and f(x) = |x|.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 40
• With or without reorthogonalization, Lanczos Jacobi matrices differed by order ∥A∥2, yet modified moments for µ = µT a,b from Algorithms SM2.5 and SM2.6 differed only in the 12th and 14th digits, without noticeably impacting interpolation-quadrature convergence.Source: SUPPLEMENTARY MATERIALS: RANDOMIZED… · p. 41