Randomized Algorithms for Estimating the Trace of an Implicit Symmetric Positive Semi-Definite Matrix IV: Sample Bound for Hutchinson's Trace EstimatorResearch Paper
Motivation
Many computations need the trace of a matrix that is never formed explicitly and can only be applied to vectors: the trace of a matrix function such as or in statistics and lattice QCD, the Frobenius norm of an operator, or the number of triangles of a graph. The standard tool is Monte-Carlo estimation, introduced by M. F. Hutchinson (Hutchinson 1989): average quadratic forms over random sign vectors . Each sample costs one matrix–vector product, uses one random bit per entry, and needs only additions and subtractions.
Before Avron and Toledo (2011), only the variance of such estimators had been analysed. A small variance does not say how many samples guarantee a given relative error with a given probability. Avron and Toledo gave the first bounds of this kind for several estimators. This mission covers the bound for Hutchinson's estimator, their Theorem 7.1.
Setting
A Rademacher random variable takes the values and , each with probability . Let be symmetric positive semi-definite. Draw independent random vectors whose entries are independent Rademacher variables. Hutchinson's trace estimator is
A single sample is an unbiased estimator of (Lemma 2.1 of the paper, due to Hutchinson). For symmetric its variance is , twice the squared Frobenius mass of off the diagonal.
Given and , a random estimator is an -approximator of if
The rank is the number of nonzero eigenvalues of , counted with multiplicity.
Formalization targets
Goal: Theorem 7.1, sample bound for
For every symmetric positive semi-definite , every and every ,
The bound depends on only through its rank. It does not depend on the dimension , on the condition number, or on how the trace is spread over the diagonal.
Milestones
- Lemma 7.2 (Achlioptas 2001, Lemma 5). For a unit vector and , for every ,
- Per-direction bound (proof of Theorem 7.1, p. 8:11). Let and . If , then .
- From directions to the trace (proof of Theorem 7.1, p. 8:11). Write and . If for every with , then . This step is deterministic.
- Lemma 2.1 (Hutchinson). , and for symmetric , .
Significance
Theorem 7.1 gives a practitioner an explicit number of matrix–vector products after which Hutchinson's method is guaranteed to reach relative accuracy with confidence . No bound was available before, although the method had been in wide use for two decades. The bound exceeds the one the same paper proves for Gaussian test vectors () by a factor. The authors conjecture that this factor is not needed. Later work removed it: Roosta-Khorasani and Ascher (2015) proved a rank-free bound for the Rademacher case, and Cortinovis and Kressner (2022) extended sample bounds to indefinite matrices. Sample bounds of this type underlie variance-reduced estimators such as Hutch++ (Meyer et al. 2021).
The theorem is proved in the literature. To the platform's knowledge it has no machine-checked proof. Formalizing it produces a checked Rademacher concentration inequality for averages of squared linear forms (Achlioptas' lemma, which is also a core lemma of database-friendly Johnson–Lindenstrauss projections), and a checked spectral reduction from a quadratic-form estimator to its eigen-directions. Neither is currently in Mathlib.
Difficulty
The estimator is not an average of independent copies of a bounded variable with a small range: a single sample can deviate from by an amount comparable to . Chebyshev's inequality with the variance of Lemma 2.1 only gives a polynomial dependence on . Hoeffding's inequality applied to directly gives a dependence on and on the size of the entries of . Unlike the Gaussian case, the estimator cannot be written as a weighted sum of independent chi-squared variables, because rotating a Rademacher vector does not give another Rademacher vector. The central difficulty is Lemma 7.2: a tail bound for that holds uniformly over every unit direction , including directions in which is far from Gaussian (for it is a constant).
Formalization scope
Matrices are Matrix (Fin n) (Fin n) ℝ, and positive semi-definiteness is Matrix.PosSemidef. The Rademacher law is on ℝ. The sample space is Fin M → Fin n → ℝ with the product of copies of this law (hutchinsonSampleMeasure), so the law of the estimator is constructed, not assumed. . Probabilities are Measure.real, and the -approximator is Definition 4.1 verbatim. is Matrix.rank. In Lemma 7.2 the i.i.d. copies are the functions on this product space.
Corrections of printed statements, each recorded in the item's Formalization Note:
- Theorem 7.1 is printed without a range on . Its proof needs , which holds exactly when . Without a range the printed statement is false: take , , and . The condition admits , while . The goal and milestone 2 are therefore stated for .
- Lemma 2.1 is printed for an arbitrary matrix. The variance formula fails for non-symmetric : for the variance is , not . Symmetry is assumed for the variance part only.
- Proof of Theorem 7.1. The proof writes together with . These fit together only for , which is the convention of milestone 3.
For the threshold involves . Lean's Real.log 0 = 0 turns the condition into , which agrees with the paper's reading . The conclusion then holds because , so no extra hypothesis is added. A formalization that took the law of the samples as a hypothesis could make that hypothesis unsatisfiable and the theorem vacuous; the constructed product space rules this out.
A complete development needs:
- the product Rademacher measure and moment generating functions of Rademacher sums;
- a Chernoff bound for averages of i.i.d. bounded variables on a product space;
- the spectral theorem for real symmetric matrices, with rank equal to the number of nonzero eigenvalues;
- a finite union bound.
The Rademacher concentration results are reusable beyond this mission. Proofs of any milestone are welcome, as are alternative proofs of Lemma 7.2.
Selected references
- H. Avron and S. Toledo, Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix, Journal of the ACM 58(2), Article 8, 2011. https://doi.org/10.1145/1944345.1944349
- M. F. Hutchinson, A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines, Communications in Statistics – Simulation and Computation 18(3), 1059–1076, 1989. https://doi.org/10.1080/03610918908812806
- D. Achlioptas, Database-friendly random projections, Proceedings of PODS 2001, 274–281. https://doi.org/10.1145/375551.375608
- F. Roosta-Khorasani and U. Ascher, Improved bounds on sample size for implicit matrix trace estimators, Foundations of Computational Mathematics 15, 1187–1212, 2015. https://doi.org/10.1007/s10208-014-9220-1
- A. Cortinovis and D. Kressner, On randomized trace estimates for indefinite matrices with an application to determinants, Foundations of Computational Mathematics 22, 875–903, 2022. https://doi.org/10.1007/s10208-021-09525-9
- R. A. Meyer, C. Musco, C. Musco and D. P. Woodruff, Hutch++: Optimal stochastic trace estimation, SOSA 2021, 142–155. https://doi.org/10.1137/1.9781611976496.16