Stochastic Lanczos Quadrature: Scalable Spectral Estimation

• Stochastic Lanczos Quadrature is a method that approximates spectral sums by combining randomized trace estimation, Lanczos Krylov projections, and Gauss quadrature.
• It efficiently computes key quantities such as log-determinants, spectral densities, and traces while offering strong theoretical guarantees on error and convergence.
• Variants like Block-SLQ and Hutch++ enhance efficiency by reducing variance and computational complexity in large-scale linear algebra applications.

Stochastic Lanczos Quadrature (SLQ) is a randomized matrix-free algorithm for efficiently estimating spectral sums and associated quantities—such as log-determinants, spectral densities, and traces of analytic matrix functions—of large Hermitian or symmetric matrices, given only the ability to perform matrix–vector products. SLQ couples randomized trace estimation (Girard–Hutchinson or Hutch++), Krylov subspace projection (Lanczos process), and Gauss quadrature to provide strong theoretical guarantees for accuracy and computational complexity, and it forms the basis for several state-of-the-art large-scale linear algebra routines in scientific computation, machine learning, and statistical inference.

1. Mathematical Formulation and Algorithmic Framework

SLQ targets the approximation of spectral sums of the form

$$\mathrm{tr}\bigl(f(A)\bigr) = \sum_{i=1}^n f(\lambda_i)$$

where $A \in \mathbb{R}^{n \times n}$ is symmetric or Hermitian, $f$ is analytic on an ellipse containing the spectrum $\{\lambda_i\}$ of $A$, and direct access to all entries of $A$ is not assumed. Instead, the algorithm relies on matrix–vector product oracles.

The central observation is that for a random probe $v$ drawn (often as a Rademacher or standard Gaussian vector and normalized), the quadratic form $v^\top f(A) v$ is an unbiased estimator of $\frac{1}{n} \mathrm{tr}\bigl(f(A)\bigr)$, since $\mathbb{E}[v v^\top] = I/n$. Averaging over $A \in \mathbb{R}^{n \times n}$0 independent probe vectors yields the classical randomized trace estimator (Girard–Hutchinson estimator).

The innovation of SLQ lies in efficiently and deterministically approximating $A \in \mathbb{R}^{n \times n}$1 for each probe using $A \in \mathbb{R}^{n \times n}$2-step Lanczos processes:

• The Lanczos process, with starting vector $A \in \mathbb{R}^{n \times n}$3, generates an orthonormal Krylov basis $A \in \mathbb{R}^{n \times n}$4 and a tridiagonal matrix $A \in \mathbb{R}^{n \times n}$5 such that $A \in \mathbb{R}^{n \times n}$6.
• Gaussian quadrature rules constructed from the eigenpairs $A \in \mathbb{R}^{n \times n}$7 of $A \in \mathbb{R}^{n \times n}$8 give

$A \in \mathbb{R}^{n \times n}$9

where $f$0 and $f$1 is the $f$2th eigenvector of $f$3.

The SLQ estimator for the trace is then

$f$4

where $f$5, $f$6, and $f$7 arise from the $f$8th probe.

Block variants (such as BOLT) operate on blocks of orthonormal probes, leverage block Krylov spaces, and result in matrix-valued quadrature approximations with potentially improved sample complexity (Yeon et al., 18 May 2025, Zimmerling et al., 2024).

2. Theoretical Guarantees and Error Analysis

The total error of SLQ decomposes into Monte Carlo sampling error (due to random probing) and Lanczos quadrature (truncation) error:

• Sampling error: For $f$9 probes and function values $\{\lambda_i\}$0 bounded in $\{\lambda_i\}$1, concentration results (e.g., Hoeffding’s inequality or sub-Gaussian bounds) yield

$\{\lambda_i\}$2

To achieve target error $\{\lambda_i\}$3 at confidence $\{\lambda_i\}$4, one requires $\{\lambda_i\}$5 (Chen et al., 2022, Yeon et al., 18 May 2025).

• Quadrature error: For each probe, if $\{\lambda_i\}$6 is analytic in a Bernstein ellipse of radius $\{\lambda_i\}$7, the error decays exponentially:

$\{\lambda_i\}$8

with $\{\lambda_i\}$9. For $A$0 and general spectrum, careful analysis must be done to handle asymmetric node distributions (Li et al., 2023, Mbingui et al., 5 Jun 2026).

Combining both errors gives, with probability $A$1: $A$2 The precise allocation of the error budget between sampling and quadrature terms can be optimized to minimize computational work (i.e., the number of matrix–vector products) (Li et al., 2023).

For spectrum estimation (e.g., empirical spectral CDF), the SLQ approximation achieves Wasserstein and Kolmogorov–Smirnov error bounds scaling as $A$3 per probe and $A$4 over probes (Chen et al., 2021).

3. Specialized Methodologies and Practical Enhancements

3.1 Hutch++ and Low-Rank Splitting

Modern variants of SLQ, such as OSLQ (Mbingui et al., 5 Jun 2026), incorporate a two-phase procedure:

• A deterministic low-rank phase (e.g., projecting onto a dominant subspace via randomized SVD or QR decomposition applied to $A$5 operating on a random sketch $A$6).
• A stochastic residual phase (probing the orthogonal complement).

This combination, introduced in Hutch++, decomposes the trace as

$A$7

with $A$8 the projector onto the computed low-rank space. Both terms are approximated by Lanczos quadrature, using orthogonalized probes. The approach achieves lower variance and better accuracy per mat-vec cost compared to naïve Hutchinson or classical SLQ (Mbingui et al., 5 Jun 2026).

3.2 Block and Subblock Approaches

Block-SLQ or BOLT (Yeon et al., 18 May 2025) replaces scalar probes with blocks of orthonormal vectors, run through block-Lanczos recurrences. Block quadrature rules act on matrix-valued functions, providing a “self-averaging” effect, improved sample complexity (error decaying as $A$9 for $A$0 the mat-vec budget), and robustness in flat-spectrum or partial-access regimes (Zimmerling et al., 2024).

3.3 Error Monitoring and Finite-Precision Robustness

A posteriori bounds can be computed directly from quadrature weights and nodes, providing certifiable accuracy post-computation (e.g., KS or Wasserstein bounds on the reconstructed empirical spectral measure) (Chen et al., 2021, Chen et al., 2022). The algorithm is numerically stable in practice, with finite-precision loss of orthogonality ameliorated via controlled or selective reorthogonalization, and stability results from the Paige–Greenbaum–Knizhnerman theory (Chen, 2024).

4. Applications and Impact

SLQ underpins scalable computations of:

• $A$1 for Gaussian process models, uncertainty quantification, and statistical physics (Mbingui et al., 5 Jun 2026, Li et al., 2023),
• matrix norms and divergences (e.g., KL, Wasserstein-2 distance for Gaussians) (Yeon et al., 18 May 2025),
• spectral sum approximations (e.g., entropy, trace exponentials, or Stieltjes transforms),
• empirical spectral density and spectrum estimation (e.g., for Hessians, Laplacians, or covariance operators) (Chen et al., 2021).

Extensive numerical experiments confirm that SLQ and its advanced variants yield substantial reductions in computational cost compared to “vanilla” Hutchinson or Chebyshev-based polynomial methods; e.g., SLQ with 30–100 mat-vecs achieves relative log-determinant errors of $A$2 to $A$3 in large-scale real-world sparse problems (Mbingui et al., 5 Jun 2026).

Block and subblock SLQ also enable efficient computation in regimes with restricted matrix access or memory (Yeon et al., 18 May 2025).

5. Recent Advances and Analytical Refinements

Recent advances establish tight probabilistic error bounds for the log-determinant and general analytic spectral sums in the presence of asymmetric quadrature nodes, relaxing earlier restrictive assumptions of symmetric node distributions (Li et al., 2023). Optimization of error allocation between the quadrature (Lanczos) and trace estimator, rather than naïvely equal partitioning, yields provably minimal total matrix–vector complexity and improved empirical performance.

Block-Gauss and block-Gauss–Radau quadrature (Zimmerling et al., 2024) provide two-sided nested bounds (in the Loewner order) for matrix-valued resolvent functions ($A$4), with monotonicity theorems and potential-theoretic extrapolation, further accelerating convergence and providing fully computable error certificates.

Randomized block enrichment, in which structured physical or random probes are added, enhances subspace coverage and reduces stagnation in applications where the spectral measure is diffuse (Zimmerling et al., 2024).

6. Complexity, Implementation, and Hyperparameter Selection

The dominant computational cost in SLQ and its variants is the product $A$5 (number of mat-vecs with $A$6). Per-probe cost is $A$7 in sparse settings and $A$8 for dense matrices, with storage $A$9 per probe or block.

Guidelines for efficient parameter selection are:

• $v$0–$v$1 typically suffices for analytic $v$2 at moderate accuracy; block methods may require fewer steps.
• $v$3 for target trace error $v$4 at confidence $v$5.
• Error budgets between sampling and quadrature should be optimized, especially for expensive MVMs (Li et al., 2023, Mbingui et al., 5 Jun 2026).
• Preconditioning via spectrum normalization (e.g., shifting and scaling $v$6 so that its spectrum clusters near 1) can significantly reduce $v$7 (Mbingui et al., 5 Jun 2026).

SLQ exhibits favorable parallelism, since probe evaluations are independent, and storage can be reused if quadrature weights/nodes are aggregated on the fly (Chen et al., 2021, Yeon et al., 18 May 2025).

7. Comparison with Related Randomized Matrix Function Methods

SLQ is contrasted with the kernel polynomial method (KPM), which uses Chebyshev expansions and Jackson damping for spectral sum approximation. SLQ (Gauss–Lanczos) achieves geometric convergence $v$8 for analytic functions, while KPM achieves $v$9 with larger constants and requires explicit spectral interval identification. SLQ is better suited for non-smooth spectral measures and provides adaptive, sample-specific quadrature weights and nodes without prior endpoint specification (Chen et al., 2022).

Block SLQ, as in BOLT and (Zimmerling et al., 2024), enjoys superior scaling with mat-vec count in certain spectral regimes, and Hutch++/OSLQ-type methods achieve variance reduction via hybrid low-rank and stochastic splitting (Mbingui et al., 5 Jun 2026).

Summary Table: Core SLQ Variants

Variant	Key Feature	Error Rate
Classic SLQ	Scalar probes, Gauss quad	$v^\top f(A) v$0
Block SLQ	Orth. block probes	$v^\top f(A) v$1 in flat spectra
Hutch++/OSLQ	Low-rank + residual split	Reduced variance; same asymptotics with lower constants
Block Gauss–Radau	Two-sided bounds	Loewner-order monotonicity, error certificates
KPM	Chebyshev expansion	$v^\top f(A) v$2 for Lipschitz $v^\top f(A) v$3, $v^\top f(A) v$4 analytic

SLQ thus forms a unified and extensible framework for high-accuracy, scalable spectral computation in constrained and large-scale linear algebraic, machine learning, and physical modeling applications (Chen et al., 2021, Chen et al., 2022, Li et al., 2023, Zimmerling et al., 2024, Yeon et al., 18 May 2025, Mbingui et al., 5 Jun 2026, Chen, 2024).