---
title: Chebyshev–Hutchinson Method
url: https://www.emergentmind.com/topics/chebyshev-hutchinson-method
type: topic
---

# Chebyshev–Hutchinson Method

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 [2103.10516] [1606.00942] [1503.06394].

## 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}$, $f$ can be approximated by a degree-$K$ Chebyshev series
  \[
  f(x) \approx p_K(x) = \sum_{j=0}^K c_j T_j(\tau(x)),
  \]
  where $T_j$ denotes the $j$th Chebyshev polynomial. The coefficients $c_j$ admit explicit integral or discrete forms, e.g.,
  \[
  c_j = \frac{2-\delta_{j,0}}{K+1} \sum_{k=0}^K f(x_k) \cos(j \theta_k),
  \]
  where $x_k$ = $\cos(\theta_k)$, $\theta_k = \pi (k + 1/2)/(K+1)$.

- **Hutchinson’s Stochastic Trace Estimator:** For any symmetric matrix $M$, with $z$ drawn i.i.d. from the Rademacher distribution ($z \in \{\pm 1\}^d$), $\mathbb{E}[z^T M z] = \mathrm{tr}(M)$. Averaging $m$ such samples yields
  \[
  \Gamma_m = \frac{1}{m} \sum_{i=1}^{m} (z^{(i)})^T M z^{(i)},
  \]
  an unbiased estimator of $\mathrm{tr}(M)$, with tight high-probability error guarantees for appropriate $m$.

Combining both, Chebyshev–Hutchinson replaces $M=f(A)$ with a Chebyshev interpolant $p_K(A)$, giving the estimator
\[
\Gamma_m = \frac{1}{m} \sum_{i=1}^m \sum_{j=0}^K c_j (z^{(i)})^T T_j(\tau(A)) z^{(i)}.
\]
This estimator is unbiased for $\mathrm{tr}(p_K(A))$, and converges to $\mathrm{tr}(f(A))$ as $K \to \infty$ [1606.00942] [1503.06394].

## 2. Algorithmic Implementation

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

1. **Coefficient Computation:** Compute Chebyshev coefficients $c_j$ for $f$ on $[a,b]$ using the mapped Chebyshev nodes.
2. **Three-Term Recurrence:** For $j=0,\ldots, K$, generate $u_0=z$, $u_1=B z$ ($B$ is the affine-mapped $A$ to $[-1,1]$). Iterate $u_{j+1}=2 B u_j - u_{j-1}$. Accumulate $S_i = \sum_{j=0}^K c_j z^T u_j$ for each random $z$.
3. **Trace Estimation:** Average over $m$ independent $S_i$'s to produce $\Gamma_m$.

This procedure, written in pseudocode:

```python
# Input: symmetric A; function f; spectrum in [a,b]; degree K; samples m
# Output: estimate Γ ≈ tr(f(A))
Precompute Chebyshev coefficients {c_j} for f on [a,b]
for i in 1..m:
    draw z ∈ {±1}^d
    u_prev = z  # T_0(A)z
    u_curr = ((2/(b-a)) * A - ((b+a)/(b-a)) * I) @ z  # T_1(A)z
    sum_i = c_0 * (z^T u_prev) + c_1 * (z^T u_curr)
    for j in 1..K-1:
        u_next = 2 * ((2/(b-a)) * A - ((b+a)/(b-a)) * I) @ u_curr - u_prev
        sum_i += c_{j+1} * (z^T u_next)
        u_prev = u_curr
        u_curr = u_next
    record sample S_i = sum_i
Γ_m = (1/m) * sum_i S_i
```
[1606.00942] [2103.10516]

## 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 $K_1 < K_2 < \ldots < K_L = K$. For levels $\ell=1,\ldots, L$, define increments
  \[
  Q_\ell(z) = \sum_{j=K_{\ell-1}+1}^{K_\ell} c_j z^T T_j(A) z, \quad K_0=0
  \]
  and estimate each expectation $\mathbb{E}[Q_\ell(z)]$ using $N_\ell$ independent samples.
- **Optimal Sampling:** The number of samples per level $N_\ell$ should satisfy
  \[
  N_\ell = \mu \sqrt{V_\ell / C_\ell}, \quad \mu = \varepsilon^{-2}\sum_{k=1}^L \sqrt{V_k C_k}
  \]
  where $V_\ell$ is the variance, and $C_\ell$ the cost per evaluation at level $\ell$ [2103.10516].
- **Total Cost:** For target estimator variance $\varepsilon^2$, multilevel cost is
  \[
  W_{\rm ML} = \varepsilon^{-2} \left( \sum_{\ell=1}^L \sqrt{V_\ell C_\ell} \right)^2
  \]
  which is typically much less than single-level cost $K\,\varepsilon^{-2}$.

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 $f$ analytic in the Bernstein ellipse of parameter $\rho>1$, and $\max_{z\in E_\rho}|f(z)|\le U$,
  \[
  \|f - p_K\|_{L^\infty([-1,1])} \le \frac{4 U}{(\rho - 1)\rho^K}.
  \]
  For spectral sum estimation, this gives
  \[
  \left| \mathrm{tr}(f(A)) - \mathrm{tr}(p_K(A)) \right| \le \frac{4dU}{(\rho-1)\rho^K}.
  \]

- **Stochastic (Hutchinson) Error:** For symmetric $B\succeq 0$, to ensure with probability $\ge 1-\delta$,
  \[
  |\mathrm{tr}_m(B) - \mathrm{tr}(B)| \le \varepsilon |\mathrm{tr}(B)|,
  \]
  it suffices to take
  \[
  m \ge 6 \varepsilon^{-2} \ln\frac{2}{\delta}.
  \]
  [1606.00942] [2103.10516] [1503.06394]

- **Total Complexity:** For $m$ samples and polynomial degree $K$, cost is $O(mK \cdot \text{mat–vec}(A))$. Multilevel variants asymptotically improve this to $O(\varepsilon^{-2} (\sum_\ell \sqrt{V_\ell C_\ell})^2)$.

## 5. Representative Algorithms and Pseudocode

The Chebyshev–Hutchinson algorithm can be instantiated for a variety of matrix functions, including $\log\det(A)$, matrix inverse trace, Estrada index, nuclear norm, and triangle counting (via $f(x)=x^3$ for adjacency matrices). Pseudocode for core operations (see tables below) involves only matrix–vector operations and Chebyshev recurrence, with memory usage $O(d)$ beyond storing $A$. Specific parameter choices ($m$, $K$) are dictated by desired accuracy, spectrum, and function analyticity.

| Step                  | Operation                                       | Complexity         |
|-----------------------|-------------------------------------------------|--------------------|
| Coefficient Computation | Chebyshev expansion of $f$ (nodes/weights)     | $O(K)$             |
| Matvec Recurrence      | $T_j(A)z$ via three-term recurrence             | $O(K\cdot$ matvec$)$ |
| Stochastic Ensemble    | $m$ i.i.d. Rademacher vector samples            | $O(mK$ matvec$)$   |

[1606.00942] [2103.10516] [1503.06394]

## 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 $10^7$ dimensions [1606.00942] [2103.10516]. Key findings:

- **Variance Reduction:** For nuclear-norm estimation on “FA” matrix using $K=300$, $m=50$, the single-level standard error was $\approx 1.48$, multilevel $\approx 0.44$, implying order-of-magnitude reduction in work for the same error tolerance [2103.10516].
- **Robustness:** Chebyshev interpolants outperform Taylor expansions by factors of $5$–$10$ in accuracy for spectral sum problems [1606.00942].
- **Log-Determinant and Beyond:** For $\log\det(A)$, the approach enables tractable computation in high dimensions, bypassing cubic-cost Cholesky/SVD (see also [1503.06394]).
- **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 $\varepsilon$, $m \sim 6\varepsilon^{-2}\ln(2/\delta)$, $K$ set by analyticity through the Bernstein ellipse parameter $\rho$, with $K \approx (1/\ln\rho) \ln(4U/((\rho-1)\varepsilon L))$ (with $L = \min_{x\in[a,b]} |f(x)|$).
- **Implementation:** Rademacher estimators offer lower variance than Gaussian alternatives. Sparse matrix structure, whenever present, should be exploited.
- **Extensions:** The method applies to any $f$ analytic on a region containing $\sigma(A)$. The multilevel generalization enables variance/cost separation, allocation, and is compatible with control variates for further accuracy enhancement, as in triangle counting for graphs [2103.10516].

## References

- "A Multilevel Approach to Stochastic Trace Estimation" [2103.10516]
- "Approximating the Spectral Sums of Large-scale Matrices using Chebyshev Approximations" [1606.00942]
- "Large-scale Log-determinant Computation through Stochastic Chebyshev Expansions" [1503.06394]

Source: https://www.emergentmind.com/topics/chebyshev-hutchinson-method