Chebyshev–Hutchinson Method

• Chebyshev–Hutchinson is a randomized, matrix-free method that approximates spectral sums via Chebyshev polynomial expansion and stochastic trace estimation.
• It efficiently computes trace estimates using iterative matrix–vector products, avoiding explicit factorization of large symmetric matrices.
• Recent multilevel extensions reduce variance and computational cost, making the method scalable for high-dimensional applications like log-determinant evaluation.

The Chebyshev–Hutchinson method is a randomized matrix-free technique for approximating spectral sums of the form $\mathrm{tr}\bigl(f(A)\bigr)$, where $A\in \mathbb{R}^{d\times d}$ is large, symmetric, and $f$ is an analytic function defined on an interval containing the spectrum of $A$. By coupling Chebyshev polynomial approximation with the Hutchinson stochastic trace estimator, this approach provides unbiased, high-accuracy estimates that scale efficiently with matrix size and exploit only matrix–vector products, making it well-suited for cases where explicit formation or factorization of $A$ is infeasible. Recent developments introduce multilevel estimators that further optimize computational cost for a prescribed variance, enhancing practical efficiency and parallelism (Hallman et al., 2021, Han et al., 2016, Han et al., 2015).

1. Mathematical Foundation

The Chebyshev–Hutchinson methodology synthesizes two algorithmic paradigms:

• Chebyshev Polynomial Expansion: For any analytic function $f$ on $[a,b]$, and after linearly mapping $[a,b]$ to $[-1,1]$ via $\tau(x)=\frac{2x-(a+b)}{b-a}$, $A\in \mathbb{R}^{d\times d}$0 can be approximated by a degree-$A\in \mathbb{R}^{d\times d}$1 Chebyshev series

$A\in \mathbb{R}^{d\times d}$2

where $A\in \mathbb{R}^{d\times d}$3 denotes the $A\in \mathbb{R}^{d\times d}$4th Chebyshev polynomial. The coefficients $A\in \mathbb{R}^{d\times d}$5 admit explicit integral or discrete forms, e.g.,

$A\in \mathbb{R}^{d\times d}$6

where $A\in \mathbb{R}^{d\times d}$7 = $A\in \mathbb{R}^{d\times d}$8, $A\in \mathbb{R}^{d\times d}$9.

• Hutchinson’s Stochastic Trace Estimator: For any symmetric matrix $f$0, with $f$1 drawn i.i.d. from the Rademacher distribution ($f$2), $f$3. Averaging $f$4 such samples yields

$f$5

an unbiased estimator of $f$6, with tight high-probability error guarantees for appropriate $f$7.

Combining both, Chebyshev–Hutchinson replaces $f$8 with a Chebyshev interpolant $f$9, giving the estimator

$A$0

This estimator is unbiased for $A$1, and converges to $A$2 as $A$3 (Han et al., 2016, Han et al., 2015).

2. Algorithmic Implementation

The method requires only iterative matrix–vector multiplications. The Chebyshev–Hutchinson procedure follows:

1. Coefficient Computation: Compute Chebyshev coefficients $A$4 for $A$5 on $A$6 using the mapped Chebyshev nodes.
2. Three-Term Recurrence: For $A$7, generate $A$8, $A$9 ($A$0 is the affine-mapped $A$1 to $A$2). Iterate $A$3. Accumulate $A$4 for each random $A$5.
3. Trace Estimation: Average over $A$6 independent $A$7's to produce $A$8.

This procedure, written in pseudocode:

$A\in \mathbb{R}^{d\times d}$05 (Han et al., 2016, Hallman et al., 2021)

3. Multilevel Monte Carlo Extension

Significant variance and cost reduction is achieved by adopting a multilevel variant:

• Hierarchy of Approximations: Construct a telescoping sum with degrees $A$9. For levels $f$0, define increments

$f$1

and estimate each expectation $f$2 using $f$3 independent samples.

• Optimal Sampling: The number of samples per level $f$4 should satisfy

$f$5

where $f$6 is the variance, and $f$7 the cost per evaluation at level $f$8 (Hallman et al., 2021).

• Total Cost: For target estimator variance $f$9, multilevel cost is

$[a,b]$0

which is typically much less than single-level cost $[a,b]$1.

This multilevel structure allocates more samples to lower-cost, higher-variance increments, yielding substantial reduction in work for prescribed estimator accuracy.

4. Error Bounds and Complexity Analysis

The Chebyshev–Hutchinson method admits rigorous, non-asymptotic error and complexity guarantees:

• Chebyshev Truncation (Polynomial Approximation) Error: For $[a,b]$2 analytic in the Bernstein ellipse of parameter $[a,b]$3, and $[a,b]$4,

$[a,b]$5

For spectral sum estimation, this gives

$[a,b]$6

• Stochastic (Hutchinson) Error: For symmetric $[a,b]$7, to ensure with probability $[a,b]$8,

$[a,b]$9

it suffices to take

$[a,b]$0

(Han et al., 2016, Hallman et al., 2021, Han et al., 2015)

• Total Complexity: For $[a,b]$1 samples and polynomial degree $[a,b]$2, cost is $[a,b]$3. Multilevel variants asymptotically improve this to $[a,b]$4.

5. Representative Algorithms and Pseudocode

The Chebyshev–Hutchinson algorithm can be instantiated for a variety of matrix functions, including $[a,b]$5, matrix inverse trace, Estrada index, nuclear norm, and triangle counting (via $[a,b]$6 for adjacency matrices). Pseudocode for core operations (see tables below) involves only matrix–vector operations and Chebyshev recurrence, with memory usage $[a,b]$7 beyond storing $[a,b]$8. Specific parameter choices ($[a,b]$9, $[-1,1]$0) are dictated by desired accuracy, spectrum, and function analyticity.

Step	Operation	Complexity
Coefficient Computation	Chebyshev expansion of $[-1,1]$1 (nodes/weights)	$[-1,1]$2
Matvec Recurrence	$[-1,1]$3 via three-term recurrence	$[-1,1]$4 matvec$[-1,1]$5
Stochastic Ensemble	$[-1,1]$6 i.i.d. Rademacher vector samples	$[-1,1]$7 matvec$[-1,1]$8

(Han et al., 2016, Hallman et al., 2021, Han et al., 2015)

6. Numerical Performance and Applications

Empirical studies demonstrate that Chebyshev–Hutchinson and its multilevel extension deliver high accuracy and linear scaling on matrices with up to $[-1,1]$9 dimensions (Han et al., 2016, Hallman et al., 2021). Key findings:

• Variance Reduction: For nuclear-norm estimation on “FA” matrix using $\tau(x)=\frac{2x-(a+b)}{b-a}$0, $\tau(x)=\frac{2x-(a+b)}{b-a}$1, the single-level standard error was $\tau(x)=\frac{2x-(a+b)}{b-a}$2, multilevel $\tau(x)=\frac{2x-(a+b)}{b-a}$3, implying order-of-magnitude reduction in work for the same error tolerance (Hallman et al., 2021).
• Robustness: Chebyshev interpolants outperform Taylor expansions by factors of $\tau(x)=\frac{2x-(a+b)}{b-a}$4–$\tau(x)=\frac{2x-(a+b)}{b-a}$5 in accuracy for spectral sum problems (Han et al., 2016).
• Log-Determinant and Beyond: For $\tau(x)=\frac{2x-(a+b)}{b-a}$6, the approach enables tractable computation in high dimensions, bypassing cubic-cost Cholesky/SVD (see also (Han et al., 2015)).
• Limitations: For functions admitting accurate low-degree polynomial approximation (e.g., Estrada index), multilevel variance reduction may be marginal.

7. Practical Guidelines, Choices, and Extensions

• Parameter Tuning: For relative precision $\tau(x)=\frac{2x-(a+b)}{b-a}$7, $\tau(x)=\frac{2x-(a+b)}{b-a}$8, $\tau(x)=\frac{2x-(a+b)}{b-a}$9 set by analyticity through the Bernstein ellipse parameter $A\in \mathbb{R}^{d\times d}$00, with $A\in \mathbb{R}^{d\times d}$01 (with $A\in \mathbb{R}^{d\times d}$02).
• Implementation: Rademacher estimators offer lower variance than Gaussian alternatives. Sparse matrix structure, whenever present, should be exploited.
• Extensions: The method applies to any $A\in \mathbb{R}^{d\times d}$03 analytic on a region containing $A\in \mathbb{R}^{d\times d}$04. The multilevel generalization enables variance/cost separation, allocation, and is compatible with control variates for further accuracy enhancement, as in triangle counting for graphs (Hallman et al., 2021).

• "A Multilevel Approach to Stochastic Trace Estimation" (Hallman et al., 2021)
• "Approximating the Spectral Sums of Large-scale Matrices using Chebyshev Approximations" (Han et al., 2016)
• "Large-scale Log-determinant Computation through Stochastic Chebyshev Expansions" (Han et al., 2015)