Improved variants of the Hutch++ algorithm for trace estimation

1 Introduction

• Trace estimation is challenging when a matrix is available only through matrix-vector products, since explicitly computing the matrix or its diagonal can require up to O(n^3) operations.The introduction gives matrix functions such as exp(B), log(λI + B), and B^-1 as examples where matrix-vector products can be tractable.Source: §1 · p. 2
• Hutch++ reduces the matrix-vector product requirement to O(ε^-1) for relative-error estimation of positive semidefinite matrices, compared with O(ε^-2) for stochastic trace estimation alone.Its first phase forms a low-rank approximation and computes its trace; the second estimates the residual trace stochastically, with the product count prescribed in advance.Source: §1 · pp. 2–3, §1.1 · p. 3

1.1 Contributions

• The adaptive prototype guarantees a prescribed trace-error tolerance with failure probability δ while seeking near-optimal allocation between phases; A-Hutch++ is a more efficient variant.The approach does not assume positive definiteness and therefore targets absolute error; the authors report that it is never worse than Hutch++ and often outperforms it in practice.Source: §1.1 · p. 4
• For symmetric positive semidefinite matrices, Nyström++ uses one pass, retains Hutch++-like theoretical guarantees up to constants, and outperforms Single Pass Hutch++ in practice.The authors also report that it often outperforms Hutch++; Single Pass Hutch++ has nearly the same theoretical guarantees as Hutch++ but performs worse in practice.Source: §1.1 · p. 5

1.2 Notation

• The notation section defines singular values, spectral and nuclear norms, stable rank, the Moore-Penrose pseudoinverse, and standard Gaussian matrices.Source: §1.2 · p. 5

• A prototype algorithm seeks to minimize matrix-vector products while guaranteeing a bound on failure probability.It estimates the variance, equivalently the Frobenius norm, which requires additional matrix-vector products.Source: §2 · p. 5
• A-Hutch++ adaptively reuses variance-estimation matrix-vector products for trace estimation.It chooses their number in an adaptive fashion.Source: §2 · p. 5
• Although product reuse complicates the analysis, observed failure probabilities do not exceed the prescribed probability.Source: §2 · p. 5

2.1 Derivation of adaptive Hutch++

• The first phase uses 2r matrix-vector products, and the method seeks a near-optimal allocation by minimizing the total cost m(r) across both phases.The second-phase sample count M(r) is chosen to attain prescribed accuracy and success probability.Source: §2.1 · p. 6

2.1.1 Analysis of trace estimation

• A variation of the trace-estimation bound suppresses its spectral-norm term for sufficiently large sample counts; in practice, the method sets ℓ = 0 and omits the side condition, with no significant observed loss in success probability.The passage notes that the bound and lemma can be improved in Appendix A.Source: §2.1.1 · pp. 6–7

2.1.2 Finding the minimum of m(r)

• The adaptive strategy builds Q(r) column-by-column and stops when the estimated cost rises across three successive ranks; this local-minimum criterion is treated as reliable in practice.The function used to find the minimum can be updated recursively without additional matrix-vector products with A.Source: §2.1.2 · pp. 7–8

2.1.3 Estimating the Frobenius norm of the remainder

• Because computing the remainder’s Frobenius norm directly is too expensive, the method estimates it stochastically and chooses a sample count that ensures error at most ε with low failure probability.The stated relationship between sample count, α, and δ is illustrated in Figure 2.Source: §2.1.3 · pp. 8–9

2.1.4 A prototype algorithm

• |tr_adap(A) − tr(A)| ≤ ε with probability at least 1 − 2δ for the output of Algorithm 2.The algorithm can be implemented block-wise to reduce passes over A; for block size b > 1, it stops when m(r*−b) < m(r*).Source: §2.1.4 · pp. 9–10

2.2 A-Hutch++

• For fixed k, A-Hutch++ chooses alpha as large as possible because M decreases with alpha; alpha_k increases monotonically and converges to 1.Source: §2.2 · pp. 11–12
• Reusing Frobenius-norm estimation products for trace estimation makes the second-phase product count max{k, ceil(M_k)}.The resulting algorithm is Algorithm 3.Source: §2.2 · pp. 11–12
• Algorithm 3 evaluates M_k in a while loop, stops at the first k > M_k, and reuses A_rest Psi(k) to estimate tr(A_rest).Source: §2.2 · p. 12
• Dependence between Frobenius-norm and trace estimation means Lemma 2.3's proof and the same type of success guarantee do not carry over, although empirical failure probabilities remain below prescribed delta.Source: §2.2 · p. 13

2.3 Numerical experiments

• For slow singular-value decay, A-Hutch++ reached average relative error 0.001827 with 74.41 matrix-vector products, versus Hutch++ at 0.001804 with 237.7 products.For faster singular-value decay, the algorithms performed similarly; in no case did Hutch++ perform noticeably better.Source: §2.3 · p. 14, Figure 4 · p. 16

2.3.1 Synthetic matrices

• Across 100,000 synthetic-matrix trials, A-Hutch++’s empirical failure probabilities remained below the prescribed probabilities.A-Hutch++ also allocated more matrix-vector products to stochastic trace estimation for slower decay and vice versa for faster decay.Source: §2.3.1 · p. 14, Table 1 · p. 15, Table 2 · p. 18

2.3.2 Triangle counting

• For triangle counting, the study applied both algorithms to the cube of adjacency matrices from a Wikipedia vote network and an arXiv collaboration network.The networks had sizes 7,115 and 5,242, with tr(B^3) values 3,650,334 and 289,560, respectively; one product with A = B^3 requires three with B.Source: §2.3.2 · p. 18, Figure 5 · p. 19

2.3.3 Estrada index

• The Estrada-index experiment compared A-Hutch++ and Hutch++ on Roget’s Thesaurus semantic graph, approximating products with exp(B) using 30 Lanczos iterations.The paper defines the Estrada index as tr(exp(B)).Source: §2.3.3 · p. 19, Figure 6 · p. 19

2.3.4 Log-determinant

• For log-determinant estimation, the study compared both algorithms on symmetric positive definite matrices by applying them to A = log(B), using log det(B) = tr(log(B)).For the Thermomech TC matrix, products with log(B) were approximated using 35 Lanczos iterations.Source: §2.3.4 · p. 20, Figure 7 · p. 20

2.3.5 Trace of inverses

• For trace-of-inverse estimation, the study tested both algorithms on a 10,000 × 10,000 tridiagonal matrix and a 100 × 100-mesh Poisson discretization.Products with A = B^-1 were computed using Matlab backslash; Figure 8 presents the results.Source: §2.3.5 · pp. 20–21, Figure 8 · p. 21

3 Nystr¨om++

• For symmetric positive semidefinite A, Nyström++ uses the Nyström approximation and fits the streaming model.Its rank-r approximation requires r matrix-vector products, compared with 2r for randomized SVD; the method can allocate products to both approximation and residual trace estimation.Source: §3 · p. 22, §3.2 · p. 25

3.1 Analysis of Nystr¨om++

• Nyström++ is shown to enjoy the same theoretical guarantees as Hutch++, with analysis beginning from a Frobenius-norm error result for the Nyström approximation.Lemma 3.2 concerns a standard Gaussian matrix with 2k columns and k ≥ 5.Source: §3.1 · pp. 23–24

3.2 Adaptive Nystr¨om++

• An adaptive Nyström++ variant offers little advantage because estimating the residual norm requires extra matrix-vector products and the approach does not fit the streaming model.It would also forfeit Nyström++’s advantage of using r rather than 2r matrix-vector products for a rank-r approximation.Source: §3.2 · p. 25

3.3 Numerical results

• To stabilize the Nyström approximation, the implementation computes an eigendecomposition after adding a small shift, then removes the shift and clips eigenvalues below zero.The procedure returns the stabilized approximation in factored form without explicitly forming the Nyström approximation.Source: §3.3 · pp. 25–26
• Single Pass Hutch++ is weakest in all tested cases, while Nyström++ and Hutch++ perform similarly in many cases and Nyström++ sometimes outperforms Hutch++.The comparison includes matrices from earlier experiments and two matrices with exponential spectral decay; the results are shown in Figures 9–12.Source: §3.3 · p. 26

4 Conclusion

• A-Hutch++ often improves on Hutch++ and never requires more matrix-vector products for the same error, while Nyström++ performs significantly better than Single Pass Hutch++.The comparisons cover algebraic- and exponential-decay synthetic matrices, the Estrada index, and traces of matrix inverses.Source: §4 · p. 26, Figure 9 · p. 27, Figure 10 · p. 27
• A-Hutch++ estimates the trace of a symmetric matrix while attempting to minimize total matrix-vector products, without requiring users to set the product count in advance.It aims to meet a prescribed error tolerance with a controllable failure probability.Source: §4 · p. 26
• Nyström++ uses the Nyström approximation in a single pass, has Hutch++’s theoretical guarantees, and performs similarly to Hutch++.The paper contrasts its performance with that of Single Pass Hutch++, which Nyström++ significantly outperforms.Source: §4 · p. 26

Appendices

• The appendices prove a Hanson-Wright inequality slightly stronger than the bound in [6, Theorem 1], derive a version of Lemma 2.1, and show the bound is asymptotically optimal.Source: Appendices · p. 29

A Hanson-Wright Inequality

• Theorem A.1 bounds the tails of a Gaussian quadratic form for symmetric matrices by showing it is sub-exponential.Its proof uses sub-exponential tail bounds and the spectral decomposition of the matrix.Source: A Hanson-Wright Inequality · pp. 29–31
• Corollary A.5.1 applies Theorem A.1 to obtain a tail bound for the Gaussian stochastic trace estimator.The corollary is derived using the diagonal embedding trick [6, Theorem 1].Source: A Hanson-Wright Inequality · p. 31
• Lemma A.6 states that Theorem A.1 is slightly stronger than the corresponding tail bound in [6, Lemma 4].Source: A Hanson-Wright Inequality · p. 31
• Lemma A.7 gives a stable-rank version of the Gaussian stochastic trace-estimator bound for symmetric matrices using m matrix-vector products.It assumes m is sufficiently large.Source: A Hanson-Wright Inequality · p. 31

B Tight constants

• γ = 2 is the tight constant: Lemma B.1 gives the lower bound, while Lemma A.7 shows it is asymptotically reachable.The section states this is the smallest constant for which the guarantee holds with probability at least 1 − δ.Source: B Tight constants · p. 32
• As c → 0, the constant C in Lemma A.7 tends to 1, yielding the upper bound γ = 2.Source: B Tight constants · p. 32
• The lower-bound proof strengthens Lemma B.1 to sample means of i.i.d. zero-mean variables, establishing failure for some fixed δ* at all sufficiently large N.For symmetric A and γ < 2, Lemma B.1 states that such a δ* exists.Source: B Tight constants · pp. 32–33