Source-linked AI summary
Hutch++: Optimal Stochastic Trace Estimation
Raphael A. Meyer, Cameron Musco, Christopher Musco, David P. Woodruff
TL;DR
The paper addresses estimating the trace of a matrix accessible only through matrix-vector products, where reducing query complexity is central. It introduces Hutch++, which combines an approximate low-rank projection with variance-reduced Hutchinson estimation, achieving O(1/ε) queries for PSD matrices instead of Hutchinson’s O(1/ε^2). The method is nearly optimal up to logarithmic factors and performs better empirically, while its relative-error theory is mainly scoped to PSD matrices.
Problem
Trace estimation seeks accurate approximations of tr(A) when A is accessible only through matrix-vector multiplication, using far fewer than d queries.
Method
Hutch++ projects onto an approximate top-eigenspace, computes that component exactly, and applies Hutchinson’s estimator to the orthogonal residual.
Results
Hutch++ achieves the same (1 ± ε) guarantee for PSD A with O(1/ε) matrix-vector queries, compared with Hutchinson’s O(1/ε^2) bound.
Takeaways & Limitations
The query complexity is nearly optimal up to logarithmic factors, and experiments show Hutch++ substantially outperforms Hutchinson on synthetic, real-world, and non-PSD applications.
Takeaways & Limitations
The relative-error guarantees mainly apply to PSD matrices; for non-PSD matrices, the paper provides generalized guarantees and empirical results.
Abstract
from arXiv · showhide
We study the problem of estimating the trace of a matrix $A$ that can only be accessed through matrix-vector multiplication. We introduce a new randomized algorithm, Hutch++, which computes a $(1 \pm ε)$ approximation to $tr(A)$ for any positive semidefinite (PSD) $A$ using just $O(1/ε)$ matrix-vector products. This improves on the ubiquitous Hutchinson's estimator, which requires $O(1/ε^2)$ matrix-vector products. Our approach is based on a simple technique for reducing the variance of Hutchinson's estimator using a low-rank approximation step, and is easy to implement and analyze. Moreover, we prove that, up to a logarithmic factor, the complexity of Hutch++ is optimal amongst all matrix-vector query algorithms, even when queries can be chosen adaptively. We show that it significantly outperforms Hutchinson's method in experiments. While our theory mainly requires $A$ to be positive semidefinite, we provide generalized guarantees for general square matrices, and show empirical gains in such applications.
1 Introduction
Trace estimation seeks accurate approximations using far fewer than d matrix-vector queries. Hutch++ reduces Hutchinson’s query complexity from O(1/ε^2) to O(1/ε) for PSD matrices by projecting onto an approximate top-eigenspace and applying Hutchinson to the residual.
- Problem: Trace estimation targets approximating tr(A) when A is accessible only through matrix-vector multiplication, using far fewer than d queries.The problem arises for implicit transformations such as powers, inverses, and matrix exponentials, where matrix-vector products can be computed more efficiently than explicitly forming A.
- Prior method: Hutchinson’s estimator achieves a (1 ± ε) guarantee for PSD A with O(1/ε^2) matrix-vector queries at constant failure probability.It averages quadratic forms from independent mean-zero, variance-one sub-Gaussian vectors.
- Hutch++: Hutch++ achieves the same (1 ± ε) guarantee for PSD A using O(1/ε) matrix-vector queries, improving the prior bound quadratically.The algorithm uses a random start matrix, computes an orthonormal basis Q for the span of AS, and estimates the residual trace with Hutchinson’s estimator.
- Variance reduction: Hutch++ reduces variance by computing the projected trace exactly and applying Hutchinson’s estimator only to the orthogonal residual.The residual variance is proportional to ∥(I − QQT)A(I − QQT)∥F^2 and is bounded by ε tr(A)^2 with m = O(1/ε).
- Optimality: The algorithm has adaptive and non-adaptive variants, while a nearly matching lower bound applies to bounded-bit-complexity query algorithms.The lower bound is obtained through a reduction from Gap-Hamming and includes a tight non-adaptive result in the real RAM model.
- Empirical results: Experiments on synthetic and real-world matrices show substantial improvement over Hutchinson and related top-eigenvalue estimators, including non-PSD applications.Applications include log determinants, the graph Estrada index, and triangle counting; the method requires only the query count m as a parameter.
- Prior work: Prior work often relied on sparsity, decomposition, or randomized SVD, whereas Hutch++ is analyzed to improve over Hutchinson even when A is far from low-rank.The analysis also shows that one power-method iteration provides a sufficiently coarse top-eigenvector approximation.
2 Preliminaries
The preliminaries define norms, traces, eigenvalue structure, and low-rank approximations used in the analysis. They focus the main guarantees on symmetric PSD matrices while also introducing non-PSD matrix notation and Hutchinson’s accuracy bound.
- Notation: For square A, tr(A) is the sum of diagonal entries, while the Frobenius norm is defined from the squared entry magnitudes.For symmetric matrices, the trace and Frobenius norm can also be expressed through eigenvalues.
- PSD matrices: The main trace-approximation results concern symmetric PSD matrices, whose eigenvalues are nonnegative and whose trace equals the ℓ1 norm of the eigenvalue vector.The Frobenius norm equals the ℓ2 norm of the eigenvalue vector.
- Low-rank approximation: For PSD A, Ak denotes the optimal rank-k approximation formed from the top k eigenvectors and eigenvalues.This approximation is used to describe spectral-tail behavior in the variance analysis.
- General matrices: For general square matrices, the singular value decomposition provides singular values and the nuclear norm used in non-PSD guarantees.The notation extends the analysis beyond symmetric PSD matrices.
- Hutchinson analysis: Hutchinson’s estimator has a high-probability accuracy bound for independent sub-Gaussian random variables and applies beyond PSD, including asymmetric matrices.Normal and ±1 random variables are both sub-Gaussian with constant parameter.
3 Complexity Analysis
Hutch++ reduces trace-estimation cost by projecting away dominant spectral components before stochastic estimation, achieving O(1/ε) queries for PSD matrices. A non-adaptive variant retains nearly the same guarantee, while the analysis extends to general square matrices through nuclear-norm bounds.
- Adaptive Hutch++: Hutchinson’s O(1/ε^2) query bound can be improved to O(1/ε) by approximately projecting away A’s largest eigenvalues before estimating the residual trace.The residual’s Frobenius norm is smaller than tr(A), reducing stochastic-estimation error.
- Adaptive Hutch++: O(1/ε) queries suffice because an approximate rank-k projection with k = O(1/ε) makes the residual accurate enough for stochastic trace estimation.The projection can be computed with O(k) matrix-vector queries using randomized linear-algebra methods.
- Adaptive Hutch++: For PSD A, Hutch++ achieves a (1 ± ε) relative-error guarantee with probability at least 1 − δ using O(log(1/δ)/ε + log(1/δ)) queries.The guarantee is stated for Algorithm 1 under the PSD assumption.
- General matrices: The same query-order guarantee extends to general square matrices, with performance expressed using the nuclear norm and improved behavior when eigenvalues decay quickly.For PSD matrices, the nuclear norm equals the trace.
- Non-adaptive Hutch++: NA-Hutch++ is a non-adaptive variant that uses O(log(1/δ)/ε) queries and achieves the same relative-error guarantee for PSD matrices.Non-adaptivity can help when block or parallel matrix-vector multiplication is faster, although experiments found slightly worse performance.
4 Lower Bounds
The paper proves that Hutch++’s O(1/ε) query complexity is nearly optimal, including for adaptive algorithms under bounded-precision queries. A separate real-RAM lower bound is tight up to constants for non-adaptive algorithms.
- Adaptive lower bound: Hutch++’s O(1/ε) query complexity is optimal up to a logarithmic factor for adaptive PSD trace-estimation algorithms with bounded-precision integer queries.The proof reduces Gap-Hamming to trace estimation.
- Adaptive lower bound: Any bounded-precision adaptive algorithm needs enough queries to achieve a (1 ± ε) trace estimate with probability greater than 2/3.The query vectors have integer entries bounded in magnitude by 2^b.
- Non-adaptive lower bound: The real-RAM lower bound is tight up to constants for non-adaptively chosen real-valued query vectors.Its proof uses a hypothesis-testing reduction involving negatively spiked covariance matrices.
- Proof strategy: The adaptive lower bound follows because a sufficiently accurate trace estimator would yield a low-communication protocol for Gap-Hamming.Gap-Hamming requires Ω(n) randomized communication, while an m-query estimator would induce communication O(m·√n(log n + b)).
5 Variance Analysis
The Gaussian Hutch++ analysis establishes unbiasedness and bounds variance by removing a low-rank component before applying Hutchinson estimation to the residual. The resulting explicit bound applies to PSD matrices and extends to general matrices through the nuclear norm.
- Variance bound: Gaussian random vectors are used because the analysis relies on an expected-error randomized low-rank approximation result.The high-probability analysis uses a different result, while the variance analysis provides explicit constants.
- Estimator construction: Gaussian-Hutch++ uses m matrix-vector products to compute a low-rank basis and estimate the residual trace with Hutchinson’s estimator.Algorithm 3 samples Gaussian S and G, forms an orthonormal basis Q from AS, and combines tr(Q^T A Q) with residual estimation.
- Gaussian-Hutch++: Gaussian-Hutch++ is unbiased: E[Gaussian-Hutch++(A)] = tr(A).The estimator decomposes the trace into a projected term and a Hutchinson estimate of the projected residual.
- Variance bound: The variance satisfies Var[Gaussian-Hutch++(A)] ≤ 16 tr^2(A)/(m−2)^2 for PSD A.The explicit constant follows from optimizing the low-rank and residual-query allocation.
- General matrices: For general square matrices, the estimator remains unbiased and its variance can be expressed using the nuclear norm rather than the trace.The PSD-specific approximation lemma is replaced by a nuclear-norm version.
6 Experimental Validation
Experiments on synthetic and real trace-estimation problems show that Hutch++ and NA-Hutch++ generally outperform baseline methods, including for some non-PSD matrices. Performance depends on spectral decay, while Hutch++ remains robust when leading eigenvalues do not dominate the trace.
- Synthetic Matrices: ε = O(m^-1) for Hutch++ and NA-Hutch++, compared with ε = O(m^-1/2) for Hutchinson’s.The experiments therefore use log-log plots to compare relative error with query count.
- Synthetic Matrices: Adaptive Hutch++ consistently performs better than non-adaptive NA-Hutch++ and avoids selecting constants c1, c2, and c3.The authors consequently use Hutch++ in their real-data experiments.
- Synthetic Matrices: Hutch++ and NA-Hutch++ typically outperform Hutchinson’s and Subspace Projection on random matrices with power law spectra.Subspace Projection performs better than Hutchinson’s mainly for very fast eigenvalue decay.
- Real Matrices: For A = exp(B), Subspace Projection and Hutch++ outperform Hutchinson’s because exponentiation produces a quickly decaying spectrum.For A = log(B + λI), Hutchinson’s performs well because the spectrum is very flat, while Hutch++ remains essentially as fast.
- Real Matrices: Hutch++ reduces variance by approximately projecting off the largest-magnitude eigenvalues, whether positive or negative.This makes it less sensitive than Subspace Projection to non-monotonic error caused by positive and negative leading eigenvalues.
- Real Matrices: Hutch++ outperforms baseline methods when A is non-PSD, including approximate triangle counting on arXiv.org and Wikipedia administrator voting networks.For these experiments, A = B^3 can be applied exactly using three matrix-vector multiplications with B.
A Proof of Lemma 2
The proof of Lemma 2 applies a sub-Gaussian quadratic-form inequality to a block-diagonal representation of Hutchinson’s estimator. This yields a high-probability trace-estimation guarantee after substituting the block matrix’s trace and norms.
- Quadratic-form reduction: The Hanson-Wright inequality applies to the vectorized query because independent normal and ±1 random variables are sub-Gaussian.The required sub-Gaussian parameter is constant for both distributions used by the estimator.
- Quadratic-form reduction: Hutchinson’s estimator with ℓ sub-Gaussian queries can be represented as a quadratic form involving a block-diagonal matrix with ℓ copies of A.Vectorizing the query matrix gives ℓ·Hℓ(A) = g^TĀg.
- Probability guarantee: For ℓ greater than a constant multiple of log(1/δ), the resulting bound holds with probability 1 − δ.The lemma states this condition for fixed constants in the sub-Gaussian concentration inequality.
- Norm substitution: The block-diagonal matrix preserves A’s spectral norm and repeats its trace and Frobenius-norm contributions ℓ times.The proof substitutes E[g^TĀg] = ℓtr(A), together with the corresponding block-matrix norms.
B Proof of Theorem 8
The proof of Theorem 8 reduces relative-error trace estimation to distinguishing two covariance models, then proves that distinguishing them requires Ω(1/ε) queries. This establishes the lower-bound contribution for PSD trace estimation.
- Reduction to testing: The lower-bound proof reduces relative-error trace estimation to a testing problem that distinguishes two PSD matrix constructions.An estimator achieving the stated relative-error guarantee would solve the testing problem with constant success probability.
- Conclusion: Combining the reduction with the testing lower bound yields the main non-adaptive trace-estimation lower bound.The proof explicitly states that the testing result and Lemma 16 together establish Theorem 8.
- Query normalization: Without loss of generality, a non-adaptive query matrix can be taken as the first m standard basis vectors.Orthonormalization and Gaussian rotational invariance justify this reduction.
- Covariance testing: The testing problem distinguishes m Gaussian samples with covariance I from samples with covariance I − zz^T, where z is a random unit vector.The two covariance cases arise with equal probability.
- Reduction to testing: The trace-estimation problem and covariance-testing problem have matching observations through a Gaussian construction, preserving the algorithm’s success probability.This is established by constructing auxiliary Gaussian matrices and using rotational invariance.
- Testing lower bound: m = Ω(1/ε) samples are necessary for the covariance-testing problem, because fewer samples cannot distinguish the two distributions with probability at least 2/3.The argument bounds divergence using the negatively spiked Wishart model.
- Divergence bound: The divergence calculation controls the testing difficulty by analyzing the overlap of independent random unit vectors and bounding the resulting beta-distribution expression.The proof uses binomial-coefficient bounds and chooses a sufficiently small constant in m = c/ε.