Optimal Sketching for Trace Estimation
Authors: Shuli Jiang
Hai Pham
David Woodruff
Qiuyi Zhang
Venue: NeurIPS 2021 Spotlight (2021)
Keywords: Trace estimation
Query complexity
Sketching algorithms

Abstract: Matrix trace estimation is ubiquitous in machine learning applications and has traditionally relied on Hutchinson's method, which requires $O(\log(1/\delta)/\epsilon^2)$ matrix-vector product queries to achieve a $(1 \pm \epsilon)$-multiplicative approximation to $\text{trace}(A)$ with failure probability $\delta$ on positive-semidefinite input matrices $A$. Recently, the Hutch++ algorithm was proposed, which reduces the number of matrix-vector queries from $O(1/\epsilon^2)$ to the optimal $O(1/\epsilon)$, and the algorithm succeeds with constant probability. However, in the high probability setting, the non-adaptive Hutch++ algorithm suffers an extra $O(\sqrt{\log(1/\delta)})$ multiplicative factor in its query complexity. Non-adaptive methods are important, as they correspond to sketching algorithms, which are mergeable, highly parallelizable, and provide low-memory streaming algorithms as well as low-communication distributed protocols. In this work, we close the gap between non-adaptive and adaptive algorithms, showing that even non-adaptive algorithms can achieve $O(\sqrt{\log(1/\delta)}/\epsilon + \log(1/\delta))$ matrix-vector products. In addition, we prove matching lower bounds demonstrating that, up to a $\log \log(1/\delta)$ factor, no further improvement in the dependence on $\delta$ or $\epsilon$ is possible by any non-adaptive algorithm. Finally, our experiments demonstrate the superior performance of our sketch over the adaptive Hutch++ algorithm, which is less parallelizable, as well as over the non-adaptive Hutchinson's method.

Optimal Sketching for Trace Estimation
Shuli Jiang Robotics Institute
Hai Pham
Language Technologies Institute
Carnegie Mellon University shulij@andred edu
Carnegie Mellon University
htpham@cs. cmu edu
David P. Woodruff
Qiuyi (Richard) Zhang
Computer Science Department Cameeiie Mellon University dwoodruf@ca cmu. edu
Google Brain qiuyiz0google.com com
Abstract
Matrix trace estimation is ubiquitous in machine learning applications and has tra- ditionally relied on Hutchinson's 's method, which requires (llog(1/))/²² matrix- vector product queries to achieve a (1 mullipplicatiive approximation tr(AA) with failure probability on positive-semidefinic input matrices A. Recently, the Hutch++ algorithm was proposed, which reduces the number of matrix-vector queries from (1)/²) to the optimal ((1/e), and the algorithm succeeds with constant probability. However, in the high probability setting, the non- -daptive Hutch++ algorithm suffers an extra O O (logo/////// complexity. Non-adaptive methods are important, as they correspond to sketching algorithms, which are mergeable, highly parallelizable, and provide low-memory streaming algorithms well as communumniittionn distributed protocols. In this work, we close the gap between non-adaptive and adaptive algorithms, showing that even non-adaptive algorithms can achieve Ooy/ogg1 1/8)/e + log(1/8) (1/6)) matrix- vector products. In addition, we prove matching lower bounds demonstrating that, up to a log (18) factor, no further improvement in the dependence on or is possible by any non-adaptive algorithm. Finally, our experiments demonstrate the superior performance our sketch over the adaptive Hutch++ algorithm, which is less parallelizable, as well as over the non-adaptive Hutchinson's 's method.
multiplicative factor in query
Introduction
The problem of implicit matrix trace estimation arises naturally in wide range of applications [I] For example, during the training of Gaussian Process, a popular non non- parametric kernel-based method, the calculation of the marginal log-likel likelikhhhoo contains a heavy-computation term, i.e., the log determinant of the covariance matrix, log(det (K)), where K E RR" xx and n is the number of data points. The canonical way of computing log(deu (K)) is via Cholesky decomposition
K, whose time complexity is Onn²). Since log (og(det(K)) = En= log(A:), where X 's are the eigenvalues of K, one compute tr tr(log(K) (K) instead. Trace estimation combined with polynomial approximation (e.g., the Chebyshev polynomial or Stochastic Lanczos Quadrature) to log [2). or trace estimation combined with maximum entropy estimation (3) provide fast ways of estimating tr(log(K) for large-scale data. Other popular applications of implicit trace estimation include counting triangles and computing the Estrada Index in graphs 5) approximating the generalized rank of a matrix [6], and studying non-convex loss landscapes from the Hessian matrix of large neural networks (NNs) 7 g
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

To define the problem, we consider the matrix-vector product model as formalized in 19 [0], where there is a real symmetric imput matrix A E R=** that cannot be explicitly presented but one has oracle access to A via matrix-vector queries, i.e., one can obtain Aq for any desired query vector q E TR.. For example, due to tremendous amount of trainable parameters of large NNs, is often prohibitive to compute or store the entire Hessian matrix H with respect to some loss function from the parameters 77. which is often used to study the non-convex loss landscape. However, with Pearlmutter's trick [T] one can compute Hq for any chosen vector q. The goal is to efficiently estimate the trace of A, denoted by tr()), up to e error, i.e., to compute a quantity within (1=))tr((A) (A). For efficiency, such algorithms minimum number of queries q required to solve the problem is referred to as the query complexity.
randomized and succeed with probability at least 1 - 6. The
Computing matrix-vector products Aq through oracle access, however, can be costly. For exam- ple, computing Hessian-vector products Hq on large NNs takes approximately twice the time of backpropagation. When estimating the eigendensity of H, one computes tr(f(H) ()) for some density function f. and needs repeated access to the matrix- -vectrr product oracle. a result, even with Pearlmutter's s trick and distributed computation on modern GPUs, it takes 20 hours to compute the eigendensity of a single Hessian H with respect to the cross-entropy loss on the CIFAR-I0 10 dataset [1]. from a set of fixed weights for ResNet-118 D which has approximately 11 million parameters 77. Thus, it is important to understand the fundamental limits of implicit trace estimation as the query complexity in terms of the desired approximation error e and the failure probability 8.
Hutchinson's 's method 14, simple yet elegant randomized algorithm, is the ubiquitous work force for implicit trace estimation. Letting Q = (q,.... 4g E TRxx be vectors with i.i.d. Gaussian or Rademacher (i.e., +1 with equal probability) random variables, Hutchinson's method returns an estimate of tr( (A) as q Ei qf Aq: = q tr( (r( T AQ). Although Hutchinson's 's method dates back to 1990, it is surprisingly not well-understood on positive semi-definite (PSD) matrices. It was originally shown that for PSD matrices A with the qu being Gaussian random variables, order to obtain multiplicative (1 + c) approximation to tr(A) with probability at least 1 8, ((log(L/))/e) matrix- queries suffice [15].
A recent work 16 proposes a variance-reducee version of Hutchinson's 's method that shows only O(1/c) matrix-vector queries are needed to achieve (1 + c)-approximation to any PSD matrix with constant success probability, in contrast to the O(1/²²) matrix-vector queries needed for Hutchinson's s original method. The key observation is that the variance of the estimated trace in Hutchinson's method is largest when there is large gap between the top few eigenvalues and the remaining ones. Thus, by splitting the number of matrix-vector queries between approximating the top O(1/e) eigenvalues, i.e., by computing rank-O(1/e (1))) approximation to A, and performing trace estimation on the remaining part of the spectrum, one needs only O(1/e) queries in total to achieve a (1 = e) approximation to tr( (A). Furthermore, [16] shows ((1/e) queries are in fact necessary for any trace estimation algorithm, up to a logarithmic factor, for algorithms succeeding with constant success probability. While [16] mainly focuses the improvement in the query complexity with constant failure probability, we focus on the dependence on the failure probability 8.
Algorithm 1 Hutch++: Stochastic trace estima- tion with adaptive matrix-vector queries
Algorithm NA Hutthhch++: Stochastic trace esti- mation with non-adaptive matrix-vector queries
1: Input: Matrix-vector multiplication oracle for PSD matrix A E R=**, Number m of queries.
1: Input: vetri multiplication oracle for
PSD matrix A E R=x", Number of queries.
2: Output: Approximation to (r(
2: Output: Approximation to tr(A).
3: Sample S E = m and G E 1x* with i..dd
3: Fix constants e1, e2, es such that e < e2 and e1 +
02 + Cs 1.
(0, 1) entries.
4: Sample S E Ra R E R"xe2m and G E
4: Compute orthonormal basis Q n T for the span of AS via QR decomposition.
R*xram, with i.i.d. // 1) entries.
5: Z AR. W AS
5: return t tr(Q" AQ) + 3 GT
6: return
Ir((S" oym (tr GT AG) (T Z(S" 2) G))
+
AQ") A(I=)""))) G).
Achieving a low failure probability d is important in applications where failures are highly undesirable, and the low failure probability regime is well-studied in related areas such as compressed sensing [17. data stream algorithms [18] 19), distribution testing [20]. and so on. While one can always reduce the failure probability from a constant to by performing (log(1/6) independent repetitions
2

---

and taking the median, this multiplicative overhead of (1)g(1/))) can cause a huge slowdown in practice, e.g., in the examples above involving large Hessians.
Two algorithms were proposed in [6]: Hutch++ (Algorithm D.. which requires adaptively chosen matrix-vector queries and HA-HHtch+++ (Algorithm{2 2 which only requires non-adaptively chosen queries. We call the matrix- queries adaptively chosen if subsequent queries are dependent on previous queries q and observations Aq. whereas the algorithm non-adaptive if all queries can be chosen at once without any prior information about A. Note that Hutchinson's 's method uses only non-adaptive queries. [16] shows that Hutch++ can use O( O[Vlog(1/ (o//)////)///eee +
log(1/8)) adaptive matrix-vector queries to achieve (1 + E) approximation with probability least 1-ó, 6, while MA-Hutch++ can use (llo((/6)/e)) non-adaptive queries. Thus, in many parameter regimes the non-adaptive algorithm suffers an extra Vlog(1/6) multiplicative factor over the adaptive algorithm.
is important to understand the query complexity of non-adaptive algorithms for trace estimation because the advantages of non- -da-aiiit are plentiful: algorithms that require only non-adaptive queries can be easily parallelized across multiple machines while algorithms with adaptive queries are inherently sequential. Furthermore, -daatiee algorithms correspond to sketching algorithms which are the basis for many streaming algorithms with low memory [2]] or distributed protocols overhead (for an example application to low rank approximation, see [22
with
We note that there are numerous works on estimating matrix norms in a data stream [3 24 25
most of which use trace estimation subroutine.
1.1 Our Contributions
Improving the Non-adaptive Query Complexity. We give improved analysis of the query complexity of the non-adaptive trace estimation algorithm NA-Hutch++ (Algorithm (2), based on a new low-rank approximation algorithm and analysis in the high probability regime, instead applying an off-the-shelf low-rank approximation algorithm as in [16]. Instead of (llg(1/8))e) queries as shown in [16), we show that O(V/og((//)))) / + log(1 (1))) non-adaptive queries suffice to achieve a multiplicative (1 + <) approximation of the trace with probability at least 5, which matches the query complexity of the adaptive trace estimation algorithm Hutch++. Since our algorithm is non-adaptive, it can be used in subroutines in streaming and distributed settings for estimating the trace, with lower memory than was previously possible for the same failure probability.
Theorem 1.1 (Restatement of Theorem 3.I.. Let A be any PSD matrix. If NA- Hutchh+ is im- (0 Vlog(L/8) log(1/6) matrix-vector multiplication queries, then with probability 1 - 8, the oulput L of HlA-Hutch++ satisfies (1 e)tr(A) < S (1 ()rr)))(A
plemented with m 0
The improved dependence on is perhaps surprising in the non-adaptive setting, as simply repeating a constant-probabille algorithm would give an O(log(1/8)/e) dependence. Our non-adaptive algorithm as good as the best known adaptive algorithm, and much better than previous non- adaptive algorithms [16, [4]. The key difference between our analysis and the analysis in [16] is in the number of non-adaptive matrix-vector queries we need to obtain an O(1) aappoooioee rank-k approximation to A in Frobenius norm.
Specifically, to reduce the total number of matrix-vector queries, our queries are split between (1) computing Á, rank-k approximation the matrix A, and (2) performing trace estimation on A A. Let A₁ = millmankkk || A - As|| |r be the best rank-k approximation to A in Frobenius norm. For our algorithm to work, we require A - Áll S O(I)I - Aullr with probability 1 - 8. Previous results from 27 show the number of non-adaptive queries required to compute Á is O(k log(1/8)), where each query is an i.i.d. Gaussian or Rademacher vector. We prove O(k + log(1/6) non-adaptive Gaussian query vectors suffice to compute Á. Low rank approximation requires both a so-called subspace embedding and approximate matrix product guarantee (see, e.g., 128], for survey on sketching for low rank approximation), and we show both hold with the desired probability, with some case analysis, for Gaussian queries. A technical overview can be found in Section
The improvement on the number of non-adaptive queries to achieve apppaappoeee rank-k approximation has many other implications, which can be of an independent interest. For example, since low-rank approximation algorithms are extensively used in streaming algorithms suitable for low-memory settings, this new result directly improves the space complexity of the state-of-the-an
3

---

streaming algorithm for Principle Component Analysis (PCA) [29] from O(d (k (ll0)(//))) to (k+ (x+log(1/6)) for constant approximation error e. where d is the dimension of the imput.
O(d .Ok
Lower Bound. Previously, no lower bounds were known on the query complexity in terms of á in a high probability setting. In this work, we give novel matching lower bound for non-adaptive (i.e., sketching) algorithms for trace estimation, with novel techniques based on new family of hard input distributions, showing that our improved O( O /Tog(//))////// - (og(1/)) upper bound is optimal, up to log (og(1/8) factor, for any E E (0. 1). The methods previously used to prove an (2(/e) lower bound with constant success probability (up to logarithmic factors) in [16] do not apply in the high probability setting. Indeed, [16] gives two lower bound methods based on a reduction from two types of problems: (1) communication complexity problem, and (2) a distribution testing problem between clean and negatively spiked random covariance matrices. Technique (1) does not apply since there is not multi-round lower bound for the
communication problem used in [[6] that depends on One might think that since we are proving a non con-adaptive lower bound, we could use a non-adaptive lower bound for Gap-Hamming (which exists, see (T8), ), but this is wrong because even the non-adaptive lower bound in [16] uses a 2-round lower bound for HaamHHaminn, and there is no such lower bound known in terms of 6. Technique (2) also does not apply, as it involves 1/e x 1/e matrix, which can be recovered exactly with 1//e queries; further, increasing the matrix dimensions would break the lower bound as their two cases would no longer need to be distinguished. Thus, such a hard input distribution fails to show the additive (llog(1/6) term in the lower bound.
Our starting point for a hard instance is family of Wigner matrices (see Definition 2.1 shifted by an identity matrix so that they are PSD. However, due to strong concentration properties of these matrices, they can only be used to provide a lower bound of O Vlog (log(1/)/// when e < 1/ 1///log/ 1/8). Indeed, setting 8 to be constant in this case recovers the (2(1/e) lower bound shown in [16] but via completely different technique. For larger 6, we consider new distribution testing problem between clean Wigner matrices and the same distribution with a large rank-1 1 noisy PSD matrix, and then argue with probability roughly 8, all non-adaptive queries have unusually tiny correlation with this rank-1 1 matrix, thus making it indistinguishable between the two distributions. This gives the desired additive (((o)((1))) lower bound, up to a log (18) factor.
Theorem 1.2 (Restatement of aaehoreemm... I.I). Suppose A is a non-adaptive query-based algorithm that returns a (1 + e multiplicattve estimate to tr( (A) for any PSD matrix A with probability at least 1 - 8.
Then, the number of matrix-vector queries must be at least m = S ( (ooR(1/8)
log(1/5) [og(]g((1/8)
1.2 Related Work
A summary of prior work on the query complexity of trace estimation of PSD matrices is given in Table For the upper bounds, prior to the work of (30), the analysis of implicit trace estimation mainly focused on the variance of estimation with different types of query vectors. [30 gave the first upper bound on the query complexity. The work of [15] improved the bounds in [30]. On the lower bound side, although (15) gives a necessary condition on the query complexity for Gaussian query vectors, this condition does not directly translate a bound on the minimum number of query vectors. The work of [16] gives the first lower bound on the query complexity in terms of but only works for constant failure probability.
Upper Bounds
Prior Work
Query Complexity
Query Vector Type
Failure Probability
Algorithm Type
30
Ollog(1/ 1/0)/=")
Gaussian
non- aon-addaeivv
30
(rog(rank A)(d)/e))
Rademacher
non non-adaptive
1S
Ollog( 1/6)/=²)
Gaussian, Rademacher
non con-addptive adaptive non -ala-ddaptive non-adaptive
16
Ol log 1, 1/6//c // (1/8))
Gaussian, Rademacher
5
16
O log. (1/]1//8/))
Gaussian, Rademacher
6
This Work
O(ylog(1/0)) (1/6))
Gaussian
Lower Bounds
16
1111/(e ([1/(elo(1//e))
constant
adaptive non-udaptive
16
[211cc)
constant
This Work
9 / (og(1//)//
non non-adaptive
Tog 177
Table 1: Upper and lower bounds on the query complexity for trace estimation of PSD matrices.

---

2 Problem Setting
Notation. A matrix A E 1 is symmetric positive semi-definite (PSD) if it is real, symmetric and has non-negative eigenvalues. Hence, x Az 2 0 for all x E R". Let tr A) 2 1: A
denote the trace of A. Let A F (Ei= Ej=1 2 11/2 denote the Frobenius norm and ||
s1p]iv||a=1 with mean and variance o². Our analysis extensively relies on the following facts:
| denote the operator norm of A. Let M(,, o²) denote the Gaussian distribution
Definition 2.1 (Gaussian and Wigner Random Matrices), We let G ~ (() denote an X n. random Gaussian matrix with i.i.d. (O, 1) entries. We let W W(n) G + G denote an n X Wigner matrix, where G ~ (n).
Fact 2.1 (Rotational Invariance of standard Gaussian). Let R E Ruxn be an orthornormal matrix. Let g E R" be a random vector with i.i.d. (0, 1) entries. Then Rg has the same distribution as g Fact 2.2 (Upper and Lower Gaussian Tail Bounds). Letting Z ~ -N( 1) be univariate Gaussian random variable, for any t >> 0, Zr Z 2 f| = O(t exp(p - 1)
3 An Improved Analysis of NA-Hutch++
Suppose we are trying to compute a sketch so as to estimate the trace of a matrix A up to a (I + e)-factor with success probability at least 1 6. Note that we focus on the case where we make matrix-vector queries non-adaptively. For any algorithm that accomplishes this with small constant failure probability, one can simply repeat this procedure (loog(1/)) times amplify the success probability to 1 - 8. Since these queries are von-adaptive and must be presented before any observations are made, seems intuitive that the number of non-adaptive queries of NA-Hutch+* (Algorithm 2 should be ((log(1/))
as shown in [16]. In this section, we give proof sketch as to why this can be reduced to O (vv log log(1/8)/e e + log(1/6)) as stated in Theorem 3.1 All proof details are provided in the supplementary material.
Theorem 3.1. Let A be PSD matrix. If NA NA-Hutch+++ is implemented with m O( (oo((/og////// log(1 (1/6) matrix-vector multiplication queriex, then with probability 1-8, the output of NA VA-Hutch++, denoted by I, satisfies (1 - -))r(AA) (A) S 1 < (1 + e)tr(A).
NA-Hutch++ splits its matrix-vector queries between computing an 0(1) aaprrxipoeee rank-k ap- proximation Å and performing Hutchinson's estimate on the residual matrix A - Á containing the small eigenvalues. The trade-off between the rank i and the number I of queries spent on estimating the small eigenvalues is summarized in Theorem 3.2
Theorem 3.2 (Theorem 4 of (16).. Let R be PSD. E (0, N, E. Let À
and A be any matrices with tr(A) = tr(A) + tr(A) and | F & O(1)IA - AslFF where A arg mill Av ||A A&l Let H((M) denote Hutchinson's trace estimator with I queries on matrix M. For fixed constants e, C, if Z clog( ()), then with probability 1 - 5, for 7
tr( A) + (A),,, we have 12 - tr()) CV (/g(1/8) tr(
The total number of matrix-vector queries directly depends on the number of non-adaptive queries required to compute an -appapppox rank-k approximation Å. Consider S E Raxam R E Raxemm for some constants C1, 02 > as defined in Algorithm 2 and set our low rank approximation
of A to be Á A= AR(S" AR) (RR)((AA))) (AS)). The standard analysis [66] applies result from streaming low-rank approximation in [27]. which requires
O(k log lo8)1 to get A Ál < O(L) A A&IlP with probability - 8. [16] then sets & = (1/e) and O (log((//))// in Theorem 3.2 to get a (I + E) approximation to tr(A). However, the right-hand side of Theorem 3.2 suggests the optimal split between à and should be k = I. The reason [16] cannot achieve such an optimal split is due to a large number of queries to compute the ((1))appro rank k-approximation. We give an improved analysis of this result, which may be of independent interest.
To get O(1) low rank approximation error, we need the non-adaptive query matrices S, R to satisfy two properties: the subspace embedding property (see Lemma 3.3), and an approximate matrix product for orthogonal subspaces (see Lemma 4 While it is known that m O(k + log(1/8)
5

---

suffices to achieve the first property, we show that m = O( + log(1/6)) suffices to achieve the second property when S, are matrices with i.i.d. Gaussian random variables, stated in Lemma 3.4
Lemma 3.3 (Subspace Embedding (Theorem 6 of (28)). Given 8 (0, 1 and e E (0,1). Ler S E IR" be a random matrix with i...d. Gaussian random variables (N,, +). Then anyy fixed
d-dimensional subspace A E Rrndd and for r = O( + log( ())) the sollowing holds with probability 1 - 6 simultaneously for all x E R2 || |2 = (1 +
Lemma 3.4 (Approximate Matrix Product for Orthogonal Subspaces), Given 8 E (0, 2)), let U E Bonk W E RRxx be two matrices with orthonormal columns such that U 0, p 2 max(k, log(1/0)), rank(U) = k and rank(W)) p. Lel S E Brxx be a random matrix with i.i.d. Gaussian random variables N(0, t). For r O(k + log( )) the following holds with probability 1 - 6, | ST SW O(1)]]W]]..
Note that we will apply the above two lemmas with constant c. The proof intuition is as follows: consider a sketch matrix S of size with i.i.d. (O, F) random variables as in Lemma 3.4 The range of U E m corresponds to an orthonormal basis of a rank-k low rank approximation to A. and the range of W R= is the orthogonal complement. Note that both SU and SW are random matrices consisting of i..dd (o, :) random variables and thus the task is to bound the size, in Frobenius norm, of the product of two random Gaussian matrices with high probability. Intuitively, the size of the matrix product is proportional to the rank & and inversely proportional to our sketch size r. The overall failure probability 8, however, is inversely proportional k, since k grows, the matrix product involves summing over more squared Gaussian random variables, i.e., x² random variables, and thus becomes even more concentrated. We show that for k Z log(1/5) 6). sketch size of ((k) suffices since the failure probability for each x² random variable is small enough to pay a union bound over k terms. On the other hand, when k < Log(1/5),, we show that r (log(1/8) suffices for the union bound. Combining the two cases gives r = O + log(1/5)).
Having shown the above, we next show that the low rank approximation error, i.e., || 14-A] Á is upper bounded by: 1) the inflation in eigenvalues by applying sketch matrix S as in Lemma 3.3 and 2) the approximate product of the range of a low rank approximation to A and its orthogonal complement, as in Lemma 3.4 Together these show that m O(k + log(1/6)) suffices for Á to be an (1))aapproximate rank-k approximation to A with probability 1 6, as stated in Theorem 3.5 Note that in both Lemma (.3 and Lemma 3.4 the entries of the random matrix are scaled Gaussian random variables (O, 7). However, when one sets the low rank approximation as Ā = AR(ST AR) (AS)", the scale cancels and one can choose standard Gaussians in the sketching matrix for convenience as in Theorem 3.5
Theorem 3.5. Let A E IR" be an arbitrary PSD matrix. Let A arg min k.n ||A Au F be the optimal rank-k approximation to A in Frobenius norm. If S E R= and R E R xmx are random matrices with i..dd (N, 1) entries for some fixed constant c >> 0 with m O + Olgglgg(1//0 (1))). then with probability 1 - 5, the matrix A (SR) AR) (AS)" satisfies A A P S
O())IAA - AAx
This improved result enables us to choose A: = 01 llllll/lo1/////// in Theorem 3.2 and combined with Theorem 3.5 this shows that only O(y/og(L/o)/ + log(1/8)) matrix-vector queries are needed to output number in (1 + c)tr(AA) with probability 1 - 8, as we conclude in Theorem 3.1
4 Lower Bounds
In this section, we show that upper bound on the query complexity of non-adaptive trace estimation is tight, up to a factor of Olog (1)()
Theorem 4.1 (Lower Bound for Non- Non-Adaptive Queries). Let e E (0, Any algo- rithm that accesses real PSD matrix A through natrix-vector multiplication queries Aqi, Aqu,... Aqm, where qi,..., qm are real-valued, non non-adaptively chosen vectors, requires
m- -Ω /log(1/5) log(1/8) log 1og[1/55
queries to outpul an estimate such that with probability least
1 8, (1- (A st (1 e)tr(A) (A).

---

Our lower bound hinges on two separate cases: we first show an S2( (1og(1/6)/ lower bound in log(1/8) lower bound in Section 4.1 whenever O (log(1/))) Second, we show an 52( log log(1/6)
Section 4.2 that applies to any c (0. 1). Observe that for c < 1/ //log(1/)) the first lower bound holds; for 2 1/ (lv/oo/g/))) our second lower bound dominates. Therefore, combining both lower bounds implies that for every and 5, the query complexity of O(v (oog//// + log(1/6)) for non-adaptive trace estimation is tight, up to a log log(1/6) factor.
We now give proof sketch of the two lower bounds. All details are in the supplementary material. Our lower bounds crucially make use of rotational invariance of the Gaussian distribution (see 2. to argue that the first y queries are, w.l.o.g., the standard basis vectors . Eq.. Note that our queries can assumed to be orthonormal. Both lower bounds use the family of n x n Wigner matrices (see Definition 2.1) with shifted mean, i.c. W + C for some C' > 0 depending [[ olp as part of the hard input distribution. The mean shift ensures that our ultimate instance is PSD with high probability.
4.1 Case 1: Lower Bound for Small
The first lower bound is based the observation that due to rotational invariance, the not-yet-queried part of W distributed almost identically to W, up to some mean shift, conditioned on the queried known part, no matter how the queries are chosen. The sum of diagonal entries of the not-yet-queried part is Gaussian, and this still has too much deviation to determine the overall trace of the input to a (1 = <) factor when n = Vlog(1/5)/e and 1/ /lgg(()))).
Theorem 4.2 (Lower Bound for Small e). For any PSD matrix A and all = 0(1 O(1//og(1 (6)), any algorithm that succeeds with probability at least 5 in outputting an estimate such that - e)) < t < (1 + e))r(A), requires m = S2( 100-m=//0(01111110 llog(1/8)//) matrix-vector queries.
4.2
Case 2: Lower Bound for Every
The second lower bound presented in Theorem 4.3 is shown via reduction to a distribution testing problem between two distributions presented in Problem 4.4
Theorem 4.3 (Lower Bound on Non-adaptive Queries for PSD Matrices). Let 6 E (0, Any
algorithu that accesses real, PSD matrix A through matrix-vector queries Aqı Aqg........ ........ Aqm.
where qı....,qu are real-valued non-adaptively chosen vectors, requires m Q log log( log(1/8) to output an estimate such that with probability at least 1 5, (1 e)tr A) < t (1 + tr (A).
log(1/8)
In the distribution testing problem, we consider Wigner matrices W W(log(1/5) shifted by O @(vlogí (1/))((1)) The problem requires algorithm for distinguishing between sample Q from this Wigner distribution and a sample P from this distribution shifted by a random rank-1 PSD matrix. The rank-1 matrix is the outer product of random vector with itself and is chosen to provide constant factor gap between the trace of P and Q
Problem 4.4 (Hard PSD Matrix Distribution Test), Given 8 E (). ser log(1/8). Choose
E R to be an independent random vector with i...d. (0, 1) entries. Consider two distributions:
Distrribution P on matrices {C'log 212
+W -2 )) some fixed gg"
constant C' 1.
Distribution Q on matrices { (W + +2 1og
where W ~ W(n) in Definition 2.1 Let A be random matrix drawn from either P or Q with equal probability. Consider any algoritim which, for a fixed query matrix Q E Raxx, observes AQ, and guesses if A or A ~ Q with success probability at least 1 - 8.
We then show in Lemma 4.5 that any algorithm which succeeds with probability 1 5 in distinguish- Tog(1/8) ing P from @ requires 2 log (og(1 non-adaptive matrix-vector queries.
Due to rotational invariance and since queries are non-adaptive, the first g queries are the first standard unit vectors. By Fact 2.2 with probability at least [og[1/8] however, single coordinate of

---

has absolute value at most log(1/8) By independence, with probability at least ( (og(1 Y, all of the first q coordinates of g are simultancously small, and thus give the algorithm almost no information to distinguish P from Q: this probability is if q O log log(1/8]
1/g(1/5)
Lemma 4.5 (Hardness of Problem 4.4). For a non-adaptive query matrix Q E R as in Prob- lem 4.4 given 5 E (0, 2). for = log(1/)), if q ol log log (18)
1og(1/8) (8) algorithm can solve Prob-
lem 4.4 with success probability 6.
5 Experiments
Part I: Comparison of Failure Probability and Running Time We give sequential and parallel implementations of the non-adaptive trace estimation algorithm NA-Hutch++ (Algorithm 2 the adaptive algorithm Hutch++ (Algorithm[) [ and Hutchinson 's method [1]. We specifically explore the benefits of the non-adaptive algorithm in a parallel setting, where all algorithms have parallel access to matrix-vector oracle. All the code is included in the supplementary material and will be publicly released.
Metrics. We say an estimate failed if on input matrix A, the estimate returned by an algorithm falls into either case: t < (1 c)tr(A) or t >> (1+ c))r We measure the performance of each algorithm by: 1) the number of failed estimates across 100 random trials, 2) the total wall-clock time perform 100 trials with sequential execution, and 3) the total wall-clock time perform 100 trials with parallel execution.
Datasets and Applications. We consider different applications of trace estimation from synthetic to real-world datasets. In many applications, trace estimation is used to estimate not only tr (r),, but also tr(((A) for some function f IR- + R. Letting A = VEVT be the eigendecompositiin of A, we have (A) := V V(E) where f(E) denotes applying f to each of the eigenvalues. Due to the expensive computation of cigendecompositions of large matrices, the matrix- multiplication vA)v is often estimated by polynomials implicitly computed via an oracle algorithm for random v. The Lanczos algorithm is very popular choice due to its superior performance (e.g. [31 2 7
We compare the performance of our trace estimation algorithms on the following applications and datasets, and use the Lanczos algorithm as the matrix-vector oracle random vector v in some particular cases.
Fast Decay Spectrum. We first consider synthetic dataset of size 5000 with a fast decaying spectrum, following [16), which is a diagonal matrix A with i-th diagonal entry A; = 1/1². Matrices with fast decaying spectrum will cause high variance in the estimated trace of Huthinson, but low variance for Hutch++ and NA-Hutch++, The matrix-vector oracle is simply Av.
.
Graph Estrada Index. Given binary adjacency matrix A E {0. 1 of graph, the Graph Estrada Index is defined as tr(exp(A)), which measures the strength of connectivity within the graph. Following [16), we use roget's 's Thesaurus semantic grapt with 1022 nodes, which was originally studied [5), and use the Lanczos algorithm with 40 steps to approximate exp(A)v as the matrix- oracle.
Graph Triangle Counting, Given a binary adjacency matrix A E {0, 1 of a graph,
the number of triangles in the graph is 1/6 tr Ar) This is an important graph summary with numerous applications in graph-mining and social network analysis (e.g. [32 33). We use arxiv_cm. _mmm the Condense Matter collaboration network dataset from arXiv This is a common benchmark graph with 23, 133 nodes and 173, 361 triangles. The matrix-vector oracle is A Note that A in this case not necessarily a PSD matrix.
Log-likelihood Estimation for Gaussian Process. When performing maximum likelihood estimation (MLE) to optimize the hyperparameters of a kernel matrix A for Gaussian Processes, one needs compute the gradient of the log-determininant of A, which involves estimating tr A 2 Following (2), we use the precipitaaiioff which consists
OOur code is available at: wwt.iiibaaaaatppppp/tttt.ab.booa fmf www. httpa///.....aaaaaaaaaaaa/aot//d
data.

---

of the measured amount of precipitation during day collected from 5,500 weather stations in the US in 2010. We sample 1,000 data points, and construct a covariance matrix A using the RBF kernel with length scale 1. We use the Lanczos algorithm with 40 steps as in 2 to approximate -1 as the -vecoo oracle.
failed estimation across trials
synthetic
roget
arxiv
precipitation
100
100
faalee
80
86
40
40
40
01
50
130
10
50
90
230
100
300
500
700
10
50
130
querles
Wall-clock time with sequential execution
synthetic
roget
arxiv cm
precipitation
500
5000
1000
segands
220
300
3000
600
100
queries
Wall-clock time with parallel execution
synthetic
roget
arxiv cm
precipitation
2500
seconds
2000
80
1000
20
10
300
500
330
queries
Figure 1: The performance comparison of Hutch++, NA-Hutch++ and Huthinson over datasets (mean + std. across 10 random runs). The approximation error for all settings is set at 0.01. Both Hutch++ and NA-Hutch++ outperform Hutchinson in terms of failed estimates. The parallel version of the non-adaptive NA-Hutch++ significantly faster than the adaptive Hutch++, making more practical in real-world applications. Legend: Hutch++ is Hutchinson is
NA-Hutch++
and
Implementation. We use random vectors with i.idd (000 1) entries as the query vectors for all algorithms. NA-Hutch++ requires additional hyperparameters to specify how the queries are split between random matrices S. R. G (see Algorithm 2 We set suggests. For each setting, we conduct estimates across 100 trials and the mean total wall-clock time (in seconds) conducting 100 trials with one standard deviation. For all of our experiments, we fix the error parameter measure the performance of each algorithm with (10,30, 50. 130, 150) queries on synthetic, roget and precipitation and with (100, 200, 700, 800} queries on arxiv which has significantly larger size. The parallel versions are implemented using Python multiprocessing package. Due to the large size of arxiv v_ _m we use sparse apparsedot__mee dot_nk10 Python wrapper for Intel Math Kernel Library (MKL) which supports fast sparse matrix-vectu multiplications, implement the matrix-vectt vector oracle for this dataset. During the experiments, we launch a pool of 40 worker processes in our parallel exccution. All experiments are conducted machines with 40 CPU cores.
C3 and C2 random runs and report the mean number of failed
[16]
0.01 and
Results and Discussion. The results Hutch++, NA-Hutch++ and Hutchinson over the 4 datasets are presented in Figure The performance of all algorithms is consistent across different datasets with different matrix- vector oracles, and even on non-PSI) instance from arxiv_ cm. Given the same number of queries, Hutch++ and NA-Hutch++ both give significantly fewer failed estimates than Hutchinson, particularly on PSD instances. is not surprising to see that Hut thinoon fails to

---

bexe62p
Inservee
Atsorute Reletove lerer
teurt0n
Egenialees
677
0.20
Eloensoectrum
Noment
Figure 2: The eigenspectrum of the two datasets and the performance comparison of Hutch++, NA-Hutch++ and Hutchinson on maximum entropy estimation based log determinant estimation.
achieve (1 + = appproxximtno to the trace most of the time due to the high variance in its estimation, 0.01).
given a small number of queries and high accuracy requirement
For computational costs, the difference in running time of all algorithms is insignificant in our sequential execution. In our parallel execution, however, Hutch++ becomes significantly slower than the other two, NA-Hutch++ and Hutchinson, which have very little difference in their parallel running time. Hutch++ suffers from slow running time due to its adaptively chosen queries, despite the fact that Hutch++ consistently gives the least number of failed estimates.
It is not hard to that NA NA-Hutch++ gives the best trade-off between high success probability in estimating an accurate trace with only a few number of queries, and a fast parallel running time due to the use of non-adaptive queries, which makes NA HA-Hutch++ more practical on large, real-world datasets. We remark that although the Lanczos algorithm is adaptive itself, even with a sequential matrix-vectes oracle, non-adaptive trace estimation can still exploit much more parallelism than adaptive methods, as shown by our experiments.
Part II: Comparison of Performance on Log Determinant Estimation We give additional experiment to compare the performance Hutch++, NA-Hutch++ and Hutchinson on estimating tr(log( K)), for some covariance matrix K. Estimating log(det((K) is required when computing the marginal log- ligliiilihood in large-scale Gaussian Process models. Recently, B proposed a maximum entropy estimation based method for log determinant estimation, which uses Hutchinson's trace estimation as a subroutine to estimate up to the k-th moments of the eigenvalues, given fixed k. The i-th moment of the eigenvalues is E[
} tr tr'K),, where K is PSD matrix, and is the vector of eigenvalues. [ shows that their proposed approach outperforms traditional Chebyshev/Lanczan polynomials for computing log (det in terms of absolute value of the relative error, i.e., abs (estimated log determinant true log determinant)/abs(tuur log determinant).
We compare the estimated log determinant of covariance matrix with different trace estimation subroutines for estimating the moments of the eigenvalues. We use 2 PSD matrices from the UFL. Sparse Matrix Collection" bcsstk20 (size 485 x 485) and bcsstm08 (size 1074 x 1074), with varying max moments {10, ,,, 30} and 30 matrix-vector queries. We repeated each run 100 times and reported the mean estimated log determinant with each trace estimation subroutine. While an improved estimate the eigenvalue moments does not necessarily lead to an improved estimate of the log determinant, it is not hard to show that accurate moment estimation does lead to improved log determinant estimation extreme cases where the eigenspectrum of K contains a few very large eigenvalues. Such case will cause Hutchinson's method to have very large variance. while method reduces the variance by first removing the large eigenvalues. The eigenspectrums of both input matrices and the results are presented in Figure 2
Conclusion
We determine an optimal @ (ol//l/g//////e//// log(1/6)) bound on the number of queries to achieve (1 é) approximation of the trace with probability 1 for non-adaptive trace estimation algorithms, up to log 1og(1/8) factor. This involves both designing a new algorithm, as well as proving a new lower bound. We conduct experiments non-adapin non-adaptive algorithm has a higher success probability compared to Hutchinson's method for the same sketch size, and has significantly faster parallel running time compared to adaptive algorithms.
synthetic and real-world datasets and confirm that our
httpps:/sppare www.eeaaame.ommm tam..edu/
10

---

Acknowledgments and Disclosure of Funding
We would like to thank the anonymous reviewers for their feedback. We are also grateful to Raphael Meyer for many detailed comments on the lower bound proofs. D. Woodruff was supported by NSF CCF-1815840, Office of Naval Research grant N80014-18-1-266,, and Simons Investigator Award.
References
[1] Shashanka Ubaru and Yousef Saad. Applications of trace estimation techniques. In Tomáš Kozubek, Martin Čermák, Petr Tichý, Radim Blaheta, Jakub Sístek, Dalibor Lukáš, and Jiff Jaroš, editors, High Performance Computing in Science and Engineering, pages 19-33, Cham, 2018. Springer International Publishing.
[2] Kun Dong, David Eriksson, Hannes Nickisch, David Bindel, and Andrew G Wilson. Scalable log determinants for gaussian process kernel learning. In Guyon, U. Luxburg, S. Ben- gio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[3] Jack K. Fitzsimons, Diego Granziol, Kurt Cutajar, Michael A. Osborne, Maurizio Filippone, and Stephen J. Roberts. Entropic trace estimates for log determinants. In ECML/PKDD (1), volume 10534 of Lecture Notes in Computer Science, pages 323 233-338.22 Springer, 2017.
[4] Haim Avron. Counting triangles in large graphs using randomized matrix trace estimation. 08 2010.
[5] Ernesto Estrada and Naomichi Hatano. Communicability in complex networks. Phys. Rev. E, 77:036111 11 Mar 2008
[6] Yuchen Zhang, Martin Wainwright, and Michael Jordan. Distributed estimation of generalized matrix rank: Efficient algorithms and lower bounds. In International Conference on Machine Learning, pages 457-465. PMLR, 2015.
[7] Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao. An investigation into neural net opti- mization via hessian eigenvalue density. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 2232-2241. PMLR, 09-15 Jun 2019.
[8] Zhewei Yao, Amir Gholami, Kurt Keutzer, and Michael Mahoney. Pyhessian: Neural networks through the lens of the hessian, 2020.
[9] Xiaoming Sun, David P. Woodruff, Guang Yang, and Jialin Zhang. Querying matrix through matrix-vector products, ACM Trans. Algorithms, 17(4), October 2021.
[10] Cyrus Rashtchian, David P. Woodruff, and Hanlin Zhu. Vector- nector-matrix-vectov queries for solving linear algebra, statistics, and graph problems. In APPROX-RANDOM, 2020.
[11] Barak A. Pearlmutter. Fast exact multiplication by the hessian. Neural Computation, 6: 6:147-160, 1994.
[12] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
[13] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770-778, 2016
[14] Michael F. Hutchinson. stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. volume 19, page 433 1990.
[15] Farbod Roosta-Khorasani and Uri Ascher. Improved bounds on sample size for implicit matrix trace estimators. Found. Comput. Math., 15(5):1887-1212 October 2015.
[16] Raphael Meyer, Cameron Musco, Christopher Musco, and David Woodruff. Hutch++: Optimal stochastic trace estimation, 2020.

---

[17] Anna C. Gilbert, Hung Q. Ngo. Ely Porat, Atri Rudra, and Martin J. Strauss. L2/12-foreach 2013.
sparse recovery with low risk. CoRR,
[18] S. Jayram and David P. Woodruff. Optimal bounds for nohnson-lindenstraul transforms and streaming problems with sub-constant error. In Proceedings of the Twenty-Second Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2011, San Francisco, California, USA, January 23-25, 2011, pages 1-10, 2011.
[19] Akshay Kamath, Eric Price, and David P. Woodruff. A simple proof of a new set disjointness with applications to data streams, 2021.
[20] Ilias Diakonikolas, Themis Gouleakis, Daniel M. Kane, John Peebles, and Eric Price. Optimal testing of discrete distributions with high probability. CoRR, abs/2009.06544 .00400, 2020.
[21] S. Muthukrishnan. Data streams: Algorithms and applications. Found. Trends Theor. Comput. Sci., 1(2), 2005.
[22] Christos Boutsidis, David Woodruff, and Peilin Zhong. Optimal principal component analysis in distributed and streaming models. In Proceedings of the 48th Annual ACM SIGACT Symposium on Theory of Compuring, STOC 2016, Cambridge, MA, USA, June 18-21, 2016, pages 236-249, 2016.
[23] Yi Li, Huy L. Nguyen, and David P. Woodruff. On sketching matrix norms and the top singular vector. In Proceedings of the Twenty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2014, Portland, Oregon, USA, January 5-7, 2014, pages 1562-1581, 2014.
[24] Yi Li and David P. Woodruff. On approximating functions of the singular values in a stream. In Proceedings of the 48th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2016, Cambridge, MA, USA, June 18-21, 2016, pages 726-739. 2016
[25] Vladimir Braverman, Stephen R. Chestnut, Robert Krauthgamer, Yi Li, David P. Woodruff, and Lin F. Yang. Matrix norms data streams: Faster, multi-pass and row-order. In Proceedings of the 35th International Conference on Machine Learning, ICMI. 2018, Stockholmsmüssonn Stockholm, Sweden, July 10-15, 2018, pages 648-657, 2018.
[26] Vladimir Braverman, Robert Krauthgamer, Aditya Krishnan, and Roi Sinoff. Schatten norms in matrix streams: Hello sparsity, goodbye dimension. In Proceedings of the 37th International Conference on Muchine Learning, ICML, 2020, 13-18 July 2020, Virtual Event, pages 1100- 1110, 2020.
[27] Kenneth L. Clarkson and David P. Woodruff. Numerical linear algebra in the streaming model. In Proceedings of the Forty-First Annual ACM Symposium on Theory of Computing, STOC '09, page 204, New York, NY, USA, 2009. Association for Computing Machinery.
[28] David P Woodruff. Sketching as a tool for numerical linear algebra. arXiv preprint arXiv:1411.4355 2014.
[29] Christos Boutsidis, David P. Woodruff, and Peilin Zhong. Optimal principal component analysis in distributed and streaming models. In Proceedings of the Firht-Eighth Annual ACM Symposium on Theory of Computing, STOC '16, page 236-249, New York, NY, USA, 2016. Association for Computing Machinery.
[30] Haim Avron and Sivan Toledo. Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. J. ACM, 58(2), April 2011.
[31] Lin Lin, Yousef Saad, and Chao Yang. Approximating spectral densities of large matrices. SIAM Review, 58, 08 2013.
[32] Mihail N. Kolountzakis, Gary L. Miller, Richard Peng, and Charalampos E. Tsourakakis. Efficient triangle counting in large graphs via bagree-based vertex partitioning. Lecture Notes in Computer Science, page 15 24, 2010.
[33] A. Pavan, Kanat Tangwongsan, Srikanta Tirthapura, and Kun-Lung Wu. Counting and sampling triangles from a graph stream. Proc. VLDB Endow., 6(14):1870-1885 September 2013.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] We clearly state the scope to
which our theorem and lower bounds hold. See Section for an overview and Section 3
Section 4 for more details.
(c) Did you discuss any potential negative societal impacts of your work? [N/A] As our work is very fundamental and theoretical, we do not see immediate potential negative societal impacts.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section and Section for our theoretical results, including all the assumptions.
(b) Did you include complete proofs of all theoretical results? [Yes] In the main paper, we give intuitive proof sketch due to limited space. See the supplementary material for all proof details.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes] See Section for a complete description used in the experiments are public. The links the datasets are provided in the footnote in Section S All code for the experiment are included in the supplementary material. (b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section 5 for all the details on how the experiments are conducted and how the hyperparameters are chosen.
the datasets and all implementation details. All datasets
(c) Did you report error bars (e.g., with respect to the random seed after running exper- iments multiple times)? [Yes] We conduct every experiment setting for 10 random runs and report the mean and one standard deviation with respect to each metric. See Section
5 for details.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Section 5 for details on total computational time and the types of machines we use.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets... Our work does not use existing assets.
(a) If your work uses existing assets, did you cite the creators? [N/A]
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
(d) Did you discuss whether and how consent obtained from people whose data you'le 're using/curating? [N/A]
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsoureing conducted research with human subjects...
Our work does not involve crowdsourcing or human subjects.
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "jkxm3_Jflp", "reviewer": "Reviewer_5vUg", "summary": "This paper contains two main contributions. First, they provided an improved analysis on the NA-Hutch++ algorithm to bridge the gap between adaptive and non-adaptive algorithm. Secondly, they also provided a nearly matching lower-bound on the complexity of the problem, closing the gap.", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I acknowledge that I am not an expert on this particular problem.Closing the complexity gap is usually a significant contribution to any problem, which indicates that improvement should be sought elsewhere. Therefore, I believe the authors have made significant contribution to the trace estimation problem where we only have access to a linear projection oracle. \n\nThe paper is well-written. Even if I did not read the proof details, the authors have kindly provided intuitions and key idea behind the proof. \n\nMy only concern is that the problem and the application itself might still not be close enough to make direct impact on real world applications, but the theoretical advances should interest the audience. \n\nupdate: thanks for the comment", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Based on my limited knowledge, I believe that it is meaningful in some applications to break down the $1\\pm\\epsilon$-multiplicative bound into an upper bound and lower bound. I am wondering if some of the results can be refined for these error bound or whether they are nearly the same. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1"}, "extra_text": {}}, {"review_id": "UkXTnPTyiAV", "reviewer": "Reviewer_j4ZM", "summary": "This paper presents an improved analysis for Hutch++, a stochastic trace estimation method proposed recently. In particular, the paper presents a new analysis for the non-adaptive variant of the Hutch++ algorithm and show that even the non-adaptive algorithm achieves optimal dependency on the failure probability \\delta (i.e., O(\\sqrt{\\log(1/\\delta)}))for the number of samples/matvecs. The paper also presents an improved lower bound for the number of matrix-vector queries required in terms of \\delta and error tolerance \\epsilon. Some numerical results are also presented to illustrate the performance of the methods.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper presents a new analysis for the Hutch++ method that improves the non-adaptive query complexity of the algorithm. It also presents an improved lower bound for the number of matrix-vector queries required. These results are certainly interesting from theoretical analysis (TCS) perspective. The numerical results presented and the python implementation with parallel processing will be useful for practitioners. \n\nHowever, I have the following comments about the paper: \n\n1. Novelty is limited: The main concern about the paper is that the work is very incremental and novelty seems limited. The main algorithms of Hutch++ and much of the analysis have already appeared in reference [16]. The proof of the main theorem 3.1 largely depends on the results in [16] (Although, I have not gone over the proof that is in the supplement in details). The theoretical improvement is really just a \\sqrt{\\log} factor wrt. failure probability, and does not have any implications on the actual usage and implementation of the algorithm. The lower bound analysis wrt. \\delta is novel, but it is not too far from the lower bounds presented in [16] (analysis approach also appears to be similar to [16]). \n\nTherefore, it is not very clear that the publication of this paper in an AI/ML conference is warranted. To my understanding, this paper is  more suitable for Theoretical Computer Science (TCS) or Numerical Linear Algebra (NLA) venue. \n\n2. Minor comment:\nI suggest changing the title of the paper since it is almost identical to the Hutch++ paper. \n\n--------------------------\nAfter Author response:\n\nI thank the authors for their response to reviewers' comments and for clarifying things. \nOne thing I should clarify, which was perhaps confusing, is the comment \"much of the analysis has already appeared in [16]\", where I was referring to the type of analysis presented, i.e., the multiplicative error bounds and lower bound analysis results that are also developed in [16]. The proof of Theorem 3.1 in this paper does follow the low rank approximation route (similar to [16] and uses the result on splitting the trace estimation), since the algorithm uses sketching to capture the top k singular values (and subspace). This is by design and the method was perhaps designed to leverage the sketching results. To my understanding, the algorithmic advancement was mainly this, using sketching to approximate the top singular values, and then using deflation and applying the stochastic trace estimator on the deflated matrix to capture the remaining part of the spectrum. This key idea has appeared in [16]. \nHowever, the novelty in this paper seems to be the new analysis of the sketched low rank approximation (I am not sure if there are changes to the algorithm itself, I request the authors to clarify). This analysis could be of independent interest, which I agree with the authors. Maybe authors can highlight this more and center some of the discussions around this point. Authors can add bullet points in 1.1 to highlight and discuss the key differences between [16] and this paper.\n\nOverall, I agree with other reviewers that the theoretical results are strong. I am happy if the paper is accepted, since the theoretical results might be for interest for other RandNLA problems, including low rank approximation, spectral sums tr(f(A), etc.\nBut, I am still not convinced entirely that the results are of interest to broad AI readers and practitioners, since the main algorithmic advances have already appeared in [16], and the theoretical improvement of sqrt log factor wrt. to the success probability will likely have no practical consequence when the method is applied to the AI applications listed by the authors. Hence I am increasing the score to 6: Marginally above the acceptance threshold.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There does not seem to be a discussion on the limitations and potential negative societal impact of their work in the paper.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "McYtPq9Q7n3", "reviewer": "Reviewer_jykV", "summary": "The paper studies the problem of estimating the trace of an input PSD matrix with the perspective of improving the query complexity of matrix-vector products required to compute an eps-multiplicative approximation of the trace. There have been two types of randomized algorithms solving this problem: adaptive and non-adaptive.\n\nThe former uses samples that depend on previous ones and therefore, by definition, is sequential. The latter doesn't, and is therefore amenable to parallelization, clearly a useful/desirable feature to have. \n\nHowever, the prior lowest query complexity results for adaptive algorithms were lower than those for non-adaptive algorithms, by a factor of sqrt(log(1/delta)), where delta is the error probability. This paper focuses on non-adaptive algorithms for trace estimation of PSD matrices. Within the framework of this problem, it has two key contributions. \n\n1. The first is an improved analysis of a non-adaptive algorithm proposed in a previous paper. The analysis is an improvement in the sense that it is able to show that this non-adaptive algorithm in fact needs just as many queries as the best adaptive algorithm, thus shaving off a factor of sqrt(log(1/delta)) from the prior analysis. \n\n2. The second contribution is a proof that the new result obtained is in fact optimal, thereby completely resolving the question of non-adaptive query complexity of matrix-vector products for trace estimation of PSD matrices. \n\nIn the process of 1., the paper also provides the following technical contribution: a proof that to estimate the k-rank approximation of a PSD matrix with probability 1-delta, it suffices to use O(k + log(1/delta)) queries, down from the previous best query complexity of O(k log (1/delta)). The paper explains how, to achieve k-rank approximation, one needs matrices that satisfy two properties: (1) a subspace embedding property, and (2) a property about orthogonal matrix products. While property (1) was satisfiable by the improved query complexity before, it wasn't known that (2) was. This paper shows, by carefully arguing about tail bounds of chi-squared distributed variables, that (2) is also in fact true for O(k + log(1/delta)) queries. ", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "REVIEW: This is a very strong paper. Both the main technical contributions (listed in the previous field) are strong in their own right: one, giving an improved analysis of an existing algorithm, and second, showing that this is the best possible query complexity for this problem, thereby effectively closing this specific problem. Further, the k-rank approximation using O(k + log(1/delta)) queries is also an improvement of a problem so fundamental, it's taught in week 4 of an introductory machine learning class. Clearly, this paper scores highly on significance and originality. \n\nThe paper scores average on clarity: some parts of it are great! I particularly liked the succinct description of the main technical idea in lines 55-60 and of the technical contribution from lines 108 - 111. However,  ​I have some specific parts of the introduction I'd love to see re-written or explained better. \n\n1. The paper's opening paragraph (lines 20-26) gives a laundry list of applications that invoke trace estimation as a subroutine. However,  it's quite difficult to actually read five other papers to see *how* this routine actually shows up there. So, I feel the authors should \"open up\" at least one (maybe the most impactful) of these applications (in just a couple of sentences) and describe the exact use of trace estimation therein. Though I don't need convincing that trace estimation is important, I am curious to see exactly where it's applicable. \n\n2. Given that the paper measures efficiency in terms of number of queries to matrix-vector products, the authors should explain at the very outset (paragraph 2: lines 27-33) *what* is the exact connection between trace estimation and matrix-vector products.  \n\nTo  summarize, I enjoyed the technical contributions and think that the problem is well-motivated. However, I think explaining the fundamental assumptions and claimed applications of your work is critical  to making the paper accessible to people who may not have done research on this topic before. \n\nI currently rate this paper 7/10 (GOOD PAPER, ACCEPT). However, I will happily give this  paper a  CLEAR ACCEPT (SCORE 8) , CONDITIONED ON the authors addressing points 1 and 2 above. I look forward to reading the authors' response. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "None, as far as I can tell. ", "review_text": "", "extra_scores": {"ethical_concerns": "None, as far as I can tell. ", "needs_ethics_review": "No", "ethics_review_area": ["I don’t know"], "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "HGgznLu6VJ6", "reviewer": "Reviewer_H6Dn", "summary": "This paper gives improved bounds for estimating the trace of a matrix using matrix-vector products. It turns what was previous a multiplicative factor of log(1 / failure probability) into an additive one, and shows that this bound is tight up to a loglog factor. These methods are evaluated on a wide range of matrices, several of which are implicitly represented (via the Lanczos method), and a better success probability was observed.", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I believe this is a solid result: trace estimators have been well studied, and the dependence on success probability is a nature parameter for analyzing algorithmic performances. The paper essentially closes off this dependence issue by showing that a additive log(1 / failure probability) is sufficient, and also necessary up to a loglog factor. They also empirically demonstrate this improvement by showing the failure probability of this method versus previous, and raises some interesting questions about better parallel performances of such estimators.\n\nThe importance of trace estimators across numerical/ML algorithms make this paper a clear accept. On the other hand, I feel the experiments demonstrate these gains are rather limited, and occur in performance only when parallelism is brought in. However, I believe they provide natural directions to extend some of the primitives developed here, and the evidences provided in the responses indicate that significantly more can be done on this front.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "My impression is that failure probability are usually an artifact of the analysis. Outer loops that utilize these trace estimators are often ok with some randomized noise in the answer, and I feel empirical results (in either positive or negative directions) on this would be interesting.\n\nAlso, the performance gains of these methods shown in Fig 1 seems to be around 1.5x, which to me feels rather marginal. It might take more end-to-end experiments to demonstrate more pronounced gains of these methods. I'm very happy that the authors took time to provide additional empirical evaluations in their responses, and am now quite interested in running my own experiments on this topic.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper considers the fundamental problem of estimating the trace of a PSD matrix from matrix vector product queries. Prior results had a gap between the adaptive and non-adaptive versions of the problem. In this paper, the authors propose a non-adaptive algorithm with nearly the same performance as the best adaptive algorithm. The paper also makes nice contributions in establishing lower bounds on the number of queries. The reviewers all felt that this paper should be published in NeuRIPS. 
Decision: Accept (Spotlight)