PaperlayerView exact revision on arXiv

Source-linked AI summary

XTrace: Making the most of every sample in stochastic trace estimation

Ethan N. Epperly, Joel A. Tropp, Robert J. Webber

arXiv:2301.07825v2math.NAcs.NA

TL;DR

Implicit trace estimation seeks the trace of a square matrix using matvec access, and existing estimators do not satisfy the exchangeability principle. This paper develops exchangeable, variance-reduced trace estimators; experiments report substantially lower errors than existing methods at fixed matvec budgets.

  • Problem

    Implicit trace estimation asks how to estimate a square matrix's trace from matvec access, while current state-of-the-art estimators do not satisfy exchangeability.

    Source: §1 · p. 2

  • Method

    The paper designs exchangeable trace estimators by symmetrizing variance-reduced estimators built from low-rank approximation and residual estimation.

    Source: §1 · p. 1, §1.2 · p. 3

  • Results

    For a fixed matvec budget, XTrace and XNysTrace reduce estimation error by several orders of magnitude compared with Hutch++.

    Source: §1.3.1 · p. 6

  • Takeaways & Limitations

    The analysis characterizes estimator performance as a function of the input matrix spectrum, and the paper also proposes exchangeable diagonal estimation with XDiag.

    Source: §1 · pp. 1–2

    The posterior error estimate can slightly underestimate the true error because individual trace estimators have a small positive correlation in practice.

    Source: §3.1 · p. 12

Abstract

from arXiv · show

The implicit trace estimation problem asks for an approximation of the trace of a square matrix, accessed via matrix-vector products (matvecs). This paper designs new randomized algorithms, XTrace and XNysTrace, for the trace estimation problem by exploiting both variance reduction and the exchangeability principle. For a fixed budget of matvecs, numerical experiments show that the new methods can achieve errors that are orders of magnitude smaller than existing algorithms, such as the Girard-Hutchinson estimator or the Hutch++ estimator. A theoretical analysis confirms the benefits by offering a precise description of the performance of these algorithms as a function of the spectrum of the input matrix. The paper also develops an exchangeable estimator, XDiag, for approximating the diagonal of a square matrix using matvecs.

1. Introduction.

The introduction frames implicit trace estimation as estimating a matrix trace from matvecs and argues that exchangeability can improve randomized estimators. It reviews Girard–Hutchinson and Hutch++ as background, noting Hutch++ reduces variance but violates exchangeability.

1.1. Stochastic trace estimators.

  • Randomized trace estimation collects matvec samples and estimates the trace from them; exchangeability requires the minimum-variance unbiased estimator to be symmetric in exchangeable test vectors.An iid family is exchangeable.Source: §1 · p. 1

1. Introduction.

  • Exchangeability motivates symmetric estimators, and the paper develops it as a design principle for randomized linear-algebra algorithms, using implicit trace estimation as its case study.The problem is to estimate the trace of a square matrix accessed through matrix–vector products.Source: §1 · p. 2

1.1.1. The Girard–Hutchinson estimator.

  • The Girard–Hutchinson estimator is unbiased under isotropic sampling, with variance decreasing at the Monte Carlo rate Θ(m^-1).The cited survey [25, §4] provides further discussion.Source: §1.1.1 · p. 2

1.1.2. The Hutch++ estimator.

  • Hutch++ combines a low-rank approximation with Girard–Hutchinson estimation of the residual; it is unbiased and has variance at most O(m^-2), but violates exchangeability.Its low-rank approximation is formed using randomized SVD.Source: §1.1.2 · p. 3

1.2. New exchangeable trace estimators.

The paper addresses Hutch++'s lack of exchangeability without combinatorially symmetrizing its sample split. It introduces XTrace and XNysTrace as exchangeable trace estimators, and extends the approach to diagonal estimation with XDiag.

1.2. New exchangeable trace estimators.

  • To avoid the combinatorial cost of symmetrizing Hutch++, the paper uses an unbalanced split and leave-one-out symmetrization, allowing each test vector to serve both approximation and residual estimation.The leave-one-out approach can be implemented at the same computational cost as Hutch++.Source: §1.2 · pp. 3–4

1.2.1. The XTrace estimator.

  • XTrace forms a low-rank approximation from all but one test vector, estimates the residual trace with the held-out vector, then averages these basic estimates.XTrace is unbiased and exchangeable.Source: §1.2.1 · p. 4

1.2.2. The XNysTrace estimator.

  • For positive-semidefinite matrices, XNysTrace uses a Nyström approximation in place of randomized SVD and produces an unbiased, exchangeable trace estimate.The Nyström method needs k matvecs for a rank-k approximation, while randomized SVD needs 2k.Source: §1.2.2 · pp. 4–5

1.2.3. Stochastic diagonal estimators.

  • XDiag extends the exchangeable-estimator approach to estimating the diagonal of an implicitly defined matrix.Source: §1.2.3 · p. 5

1.3. Numerical experiments.

The experiments illustrate how XTrace and XNysTrace use exchangeable estimation to improve trace estimates under a fixed matvec budget. They test spectral decay on a synthetic matrix and performance on a transverse-field Ising partition function.

1.3. Numerical experiments.

  • The experiments are presented as motivating numerical demonstrations of XTrace and XNysTrace, with further numerical work in Section 4.The first uses a synthetic matrix; the second applies trace estimation to the transverse-field Ising model's partition function.Source: §1.3 · p. 5, §1.3.1 · p. 5, §1.3.2 · p. 6

1.3.1. Exploiting spectral decay.

  • On a synthetic psd matrix with exponentially decreasing eigenvalues, XTrace and XNysTrace converge exponentially at 1.5× and 3× Hutch++'s rate, respectively, until machine precision.For a fixed matvec budget, both can reduce error by several orders of magnitude compared with Hutch++; exchangeability's variance reduction is almost as significant as using a low-rank control variate.Source: §1.3.1 · p. 6, Figure 1 · p. 6

1.3.2. Computing partition functions.

  • At m = 40 matvecs, XTrace is 240× more accurate and XNysTrace is 2400× more accurate than Hutch++ on the Ising partition function.At m = 10, all variance-reduced methods achieve errors 5 orders of magnitude smaller than Girard–Hutchinson; the eigenvalue ratio λ19/λ20 = 4600 helps explain when the methods capitalize on spectral decay.Source: §1.3.2 · p. 7, Figure 2 · p. 7

1.4. Theoretical guarantees.

The theoretical analysis gives variance guarantees for Hutch++, XTrace, and XNysTrace, including how performance relates to approximation rank and spectral decay. Under the stated normal-vector setting, the bounds support faster variance reduction than Girard–Hutchinson.

1.4. Theoretical guarantees.

  • For each method, m = O(η^−1/2) matvecs suffice to achieve the stated variance bound.The bounds involve rank-r approximations and matrix norms; the best rank-r approximation is simultaneous in Frobenius, spectral, and nuclear norms.Source: §1.4 · p. 8
  • For all three estimators, variance decreases as O(1/m^2), compared with Girard–Hutchinson's Θ(1/m) rate.Theorem 1.1 states the estimators are unbiased for m ≥ 8; the analysis uses standard normal test vectors to provide explicit constants.Source: §1.4 · pp. 7–8
  • The bounds favor XTrace and XNysTrace on matrices whose singular values decay rapidly because they permit larger approximation ranks.For exponentially decaying eigenvalues, the errors of the three methods have distinct decay behavior; the supplied passage does not provide its formula.Source: §1.4 · p. 8
  • For exponentially decaying eigenvalues, XTrace and XNysTrace converge exponentially at 1.5× and 3× Hutch++'s rate, matching Figure 1.Source: §1.4 · p. 8

1.5. Benefits.

The paper highlights accuracy, computational cost, and error estimation as benefits of XTrace and XNysTrace.

1.5. Benefits.

  • For a budget of m matvecs, XTrace and XNysTrace often yield errors orders of magnitude smaller than Hutch++.Source: §1.5 · p. 8
  • Their implementations use m matvecs plus O(m^2N) arithmetic, the same computational cost as Hutch++, and can provide reliable error estimates.Source: §1.5 · pp. 8–9

1.6. A brief history of stochastic trace estimation.

The section traces stochastic trace estimation from early randomized estimators through variance-reduced low-rank methods. XTrace and XNysTrace add exchangeability to that strategy.

1.6. A brief history of stochastic trace estimation.

  • Girard introduced stochastic trace estimation, and Hutchinson proposed random sign test vectors.Girard's estimator used vectors drawn uniformly from a Euclidean sphere.Source: §1.6 · p. 9
  • Variance-reduction methods use low-rank approximations; Hutch++ [26] has a worst-case variance bound of O(1/m^2).Earlier approaches include a low-rank surrogate estimator and a hybrid low-rank plus residual estimator.Source: §1.6 · p. 9
  • XTrace and XNysTrace extend low-rank variance reduction by enforcing exchangeability and using all test vectors for both approximation and residual trace estimation.Nyström++ [28] is described as a Hutch++ version for psd matrices using Nyström approximation.Source: §1.6 · p. 9
  • Other approaches estimate traces of matrix functions through stochastic Lanczos quadrature or reuse Krylov-subspace matvecs [7].Source: §1.6 · p. 9

1.7. Reproducible research.

The paper provides optimized MATLAB implementations and experiment-reproduction code, then outlines its implementation, error-estimation, numerical-experiment, and proof sections. It also establishes notation for the paper’s matrix and norm conventions.

1.7. Reproducible research.

  • Optimized MATLAB implementations and code to reproduce the experiments are available online.The paper identifies the XTrace code repository.Source: §1.7 · p. 9

1.8. Outline.

  • The paper is organized around efficient implementations, error estimation and adaptive stopping, numerical experiments, and theoretical proofs.These topics appear in sections 2 through 5, respectively.Source: §1.8 · p. 9

1.9. Notation.

  • The notation section specifies matrix and vector conventions, matrix-deletion notation, three norms, and best rank-r approximations.The norms are spectral, Frobenius, and trace/nuclear norms.Source: §1.9 · pp. 9–10

2. Efficient implementation of exchangeable trace estimators.

The implementation section makes XTrace and XNysTrace practical through efficient linear-algebra updates, recommends test-vector normalization to improve accuracy, and introduces XDiag as an unbiased, lower-variance diagonal estimator.

2.1. Computing XTrace.

  • XTrace computes the leave-one-out orthogonal projectors using rank-one updates and has O(m^2N) cost, matching Hutch++ asymptotically.The update vectors are obtained from the full basis, with the needed linear systems dominating the update calculation.Source: §2.1 · p. 10

2.2. Computing XNysTrace.

  • XNysTrace uses rank-one updates to form its estimator with m matvecs and O(m^2N) post-processing operations.The implementation also uses numerical-stability techniques from.Source: §2.2 · pp. 10–11

2.3. Normalization of test vectors.

  • Normalizing the distinguished test vectors removes variance from their random lengths and improves accuracy over unnormalized Gaussian or uniform-sphere vectors.The procedure first uses spherically symmetric vectors to form the approximation, then normalizes vectors for the basic trace estimates.Source: §2.3 · p. 11

2.4. Diagonal estimation.

  • XDiag is an unbiased diagonal estimator with lower variance, but unlike XTrace it also requires matvecs with A* as well as A.It can be implemented in O(m^2N) operations.Source: §2.4 · p. 12

3. Error estimation and adaptive stopping.

The section reuses individual trace estimates to assess error and adapt the matvec budget to a target tolerance. Its posterior error estimate is useful in practice, though small positive correlation can make it slightly underestimate true error.

3.1. Error estimation.

  • A scaled variance of the individual trace estimates provides a posterior estimate of squared trace-estimation error.The paper analyzes this estimate in Proposition 3.1.Source: §3.1 · p. 12
  • Because individual estimates are typically slightly positively correlated, the posterior error estimate usually understates true error by a small amount.The authors nevertheless describe it as a valuable practical tool.Source: §3.1 · p. 12

3.2. Adaptive stopping.

  • A doubling strategy increases the matvec budget until the estimated error meets the requested relative tolerance, using at most twice the optimal matvec count.It maintains XTrace and XNysTrace's O(m^2N) computational cost.Source: §3.2 · pp. 12–13

4. Numerical experiments.

Across synthetic spectra and applications, the experiments compare estimation errors at fixed matvec budgets and show strong performance for the exchangeable estimators, with performance depending on spectral shape and test-vector normalization. They also assess posterior error estimates and apply the methods to quantum statistical mechanics and network centralities.

4.1. Comparison of trace estimators.

  • The synthetic comparison uses six trace estimators, random-sign test vectors, four spectral profiles, and relative error averaged over 1,000 trials.The matrix dimension is N = 1000.Source: §4.1 · p. 13
  • Across synthetic instances, XTrace has smaller errors than Hutch++, sometimes by orders of magnitude; for the exp instance, its error decays 1.5× faster.The average relative errors are shown for the spectral profiles using random-sign test vectors.Source: §4.1 · p. 14, Figure 3 · p. 14
  • For the exp instance, XNysTrace converges 2× faster than XTrace and Nyström++ and 3× faster than Hutch++, but can perform poorly on slowly decaying spectral tails.The cited error bounds depend on the trace/nuclear norm, which is sensitive to slow eigenvalue decay.Source: §4.1 · p. 14

4.2. Choice of test vectors.

  • On flat and step spectra, normalized XTrace test vectors perform best among the tested distributions; differences are visible only for spectra with flat segments.Uniform-sign vectors rank next, followed by uniform-sphere vectors, while standard-normal vectors lag behind.Source: §4.2 · p. 15, Figure 4 · p. 15

4.3. Application: Quantum statistical mechanics.

  • For the quantum partition-function experiment, posterior error estimates track true errors within a factor of 3.2 across the reported trials.Figure 5 also reports XNysTrace-computed average energy per site over β, h ∈ [10^-1, 10^1] up to relative error 10^-3.Source: §4.3 · pp. 16–17, Figure 5 · p. 16

4.4. Application: Networks.

  • For yeast subgraph centrality after 200 matvecs, XDiag is five orders of magnitude more accurate than BKS and three orders more accurate than Diag++.Figure 6 compares mean errors for subgraph centralities and triangle numbers.Source: §4.4 · p. 17, Figure 6 · p. 17

5. Theoretical analysis.

The analysis proves refined variance bounds for Hutch++, XTrace, and XNysTrace, and characterizes the posterior error estimator.

5.1. XTrace variance bound.

  • XTrace is unbiased, and its variance bound depends on low-rank approximation error and uses exchangeability among its test vectors.The initial variance analysis isolates variance and covariance terms before combining their bounds.Source: §5.1 · pp. 18–19

5.2. Hutch++ variance bound.

  • Hutch++ is unbiased and receives an initial variance bound; its proof conditions on the low-rank approximation.The paper describes this bound as more elementary because it does not account for interactions between simple estimators.Source: §5.2 · pp. 19–20

5.3. XNysTrace variance bound.

  • For positive semidefinite matrices, XNysTrace is unbiased and its initial variance bound depends on randomized Nyström approximation error.The bound is developed for XNysTrace using m standard normal test vectors.Source: §5.3 · p. 20

5.4. Error bounds for low-rank approximations.

  • The analysis supplies error bounds for randomized SVD and, for positive semidefinite matrices, randomized Nyström approximation.These approximation results are used in deriving the complete variance bounds.Source: §5.4 · p. 21, §5.5 · p. 21

5.5. The complete variance bound.

  • The complete variance bounds combine estimator variance results with randomized low-rank approximation errors, using ranks ⌊m/6⌋−1 for Hutch++, ⌊m/4⌋−2 for XTrace, and ⌊m/2⌋−3 for XNysTrace.The derivation uses different approximation ranks for each estimator.Source: §5.5 · pp. 21–22

5.6. Proof of Proposition 3.1.

  • The posterior error estimator is checked against observed XTrace behavior: correlations are at most 0.06 and the average error estimate is accurate within a factor of 1.2 in one example.The example uses a matrix with exponentially decaying eigenvalues and 10^4 independent runs; the text also reports small correlations for slow singular-value decay.Source: §5.6 · p. 22

SUPPLEMENTARY MATERIAL

The supplementary material compares XTrace with adaptive Hutch++ and examines runtime trade-offs. It also describes efficient implementations of the estimators.

SM1. Comparison with Adaptive Hutch++.

  • XTrace performs similarly to adaptive Hutch++ or better by up to an order of magnitude in mean error across the tested tolerances.The comparison plots mean error and requested tolerance against the mean number of matvecs.Source: SUPPLEMENTARY MATERIAL · p. 25, Figure SM1 · p. 26

SM2. Runtime comparison.

  • At matched runtime or processing time, XTrace achieves lower error than Hutch++; XNysTrace has higher processing costs, favoring it when matvecs are expensive and XTrace when they are cheap.XTrace's processing cost is also higher than Hutch++ because it orthogonalizes a larger matrix.Source: SUPPLEMENTARY MATERIAL · pp. 25–26, Figure SM2 · p. 27

SM3. Derivation for efficient XTrace implementation.

  • The XTrace update is implemented in O(m^2N) operations, and the XNysTrace implementation is derived in a similar way.The XTrace implementation uses m/2 matvecs for each of AΩ and AQ.Source: SUPPLEMENTARY MATERIAL · p. 28

Related papers

Loading 2301.07825v2…