Dynamic Trace Estimation
Authors: Prathamesh Dharangutte
Christopher P Musco
Venue: NeurIPS 2021 Poster (2021)
Keywords: randomized numerical linear algebra
trace estimation
dynamic algorithms

Abstract: We study a dynamic version of the implicit trace estimation problem. Given access to an oracle for computing matrix-vector multiplications with a dynamically changing matrix A, our goal is to maintain an accurate approximation to A's trace using as few multiplications as possible. We present a practical algorithm for solving this problem and prove that, in a natural setting, its complexity is quadratically better than the standard solution of repeatedly applying Hutchinson's stochastic trace estimator. We also provide an improved algorithm assuming additional common assumptions on A's dynamic updates. We support our theory with empirical results, showing significant computational improvements on three applications in machine learning and network science: tracking moments of the Hessian spectral density during neural network optimization, counting triangles and estimating natural connectivity in a dynamically changing graph.

Dynamic Trace Estimation
Prathamesh Dharangutte Dept.of Computer Science Engineering New York University
Christopher Musco
Dept.of Computer Science Engineering New York University
ptd244@nyu. edu
cmusco@nyu.edd
Abstract
We study dynamicc version of the implicit trace estimation problem. Given access to an oracle for computing matrix-vector multiplicattins with a dynamically changing matrix A, our goal is maintain trace using as few multiplications as possible. We present a practical algorithm for solving this problem and prove that, in natural setting, its complexity is quadratically better than the standard solution of repeatedly applying Hutchinson's 's stochastic trace estimator. We also provide an improved algorithm assuming slightly stronger assumptions on the dynamic matrix A. We support our theory with empirical results, showing significant computational improvements on three applications in machine learning and network science: tracking moments of the Hessian spectral density during neural network optimization, counting triangles, and estimating natural connectivity in dynamically changing graph.
accurate approximation to A's
Introduction
Implicit or "matrix-free" trace estimation is a ubiquitous computational primitive in linear algebra, which has become increasingly important in machine learning and data science. Given access oracle for computing matrix-vector products A Arm between an n x n matrix A and chosen vectors Tı, ........ Zm, the goal is compute an approximation to A's trace, tr(A) = > A,,. This problem arises when A's diagonal entries cannot be accessed explicitly, usually because forming is computationally prohibitive. As example, consider A which the Hessian matrix of loss function involving a neural network. While forming the Hessian is infeasible when the network is large, backpropagation can be used to efficiently compute Hessian-vector products 32].
In other applications, A is a matrix function of another matrix B. For example, B is a graph adjacency matrix, tr( () equals six times the number of triangle in the graph (2). Computing A B explicitly to evaluate the trace would require (nn²) time, while the matrix-vector multiplication Ax B- (B (Bx:) only requires O(n²) time. Similarly, in beg-determinant approximation, useful in e.g. Bayesian log likelihood computation or determinantal point process (DPP) methods, we want approximate the trace of A log( B3 5 17 36). Again, A takes O time form explicitly, but Ax = log B)r can be computed in roughly O(n²) time using iterative methods like the Lanczos algorithm 21. Dynamic versions of the log- deteeemmmnnn estimation problem have been studied due to applications in greedy methods for DPP inference
[[5]
In data science and machine learning, other applications of implicit trace estimation include matrix norm and spectral sum estimation 16
31 as well as methods for eigenvalue counting 9 and Spectral density estimation methods typically implicit trace estimation to estimate moments of a matrix 's eigenvalue distribution - i.e., tr(A), tr (A22), tr tr(3),,, etc. approximation to that entire distribution. In deep learning,
spectral density estimation
45 27
which can then be used to compute spectral density estimation is used to quickly analyze the spectra of weight matrices [33] (8) or to probe information about the Hessian matrix during optimization
Trace estimation has also been used for neural networks weight quantization [TO 34 and to understand training dynamics 44
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

1.1 Static Trace Estimation
The mostly widely used implicit trace estimation algorithm is Hutchinson's estimator 14 22 Letting
91..... 9e E IR" be random vectors with i.i.d mean 0 and variance I entries (e.g., standard Gaussian or 1 Rademachers), Hutchinson's is approximates tr(4) via the average (4)
Σi= of (Ag:), )
This estimator requires l matrix vectorr multiplications to compute. Its variance can be shown to be ((IAA[}//e) and with high probability, when & O(1/e"). we have the error guarantee 3 30
[A)(A) tr(A) < <||A||p.
While improvements on Hutchinson's estimator have been studied for restricted classes of matrices (positive semidefinite, sparse, nearly low-rank, etc.) [30 40 38 36 19,, the method is the best known for general matrices - no techniques achieve guarantee with o o(1/e") matrix-vector products.
1.2 Dynamic Trace Estimation
We explore natural and widely applicable dynamic version of the implicit trace estimation problem: given access to a matrix-vector multiplication oracle for dynamically changing matrix A, maintain an approximation
A's trace. This problem arises in applications involving optimization in machine learning where we need to estimate the trace of a constantly changing Hessian matrix H (or some function of H) during model training. In other applications, A is dynamic because it is repeatedly modified by some algorithmic process. E.g., in the transit planning method of 433 edges are added to network to optimally increase the "Estrada index" (II. Evaluating this connectivity measure requires computing tr (r(exp(B)) (B)), where B is the dynamically changing network adjacency matrix and exp( B) is a matrix exponential. A naive solution to the dynamic problem is to simply apply Hutchinson's 's estimator to every snapshot of A as it changes over time. To achieve a guarantee like (I) for m time steps, we require Omm/ /e) matrix- -vector multiplies. The goal of this paper is to improve on this bound when the changes to A are bounded. Formally, we abstract the problem as follows:
Problem 1 (Dynamic trace estimation). Let Al.... Am be n x n matrices satisfying:
1. | S L. for all i [m. 2. ||A+1 Alle S a. for all i i [L. 1].
Given implicit matrix-vectt vector multiplication access to each A; in sequence, the goal is to compute trace approximations t1, .tm. tr A1), ..... tr() such that, for each i 1,..., m,
Pllt tr(A (A.)) c 6
Above A....., Am represent different snapshots of a dynamic matrix at m time steps. We require ||A:|| S 1 only to simplify the form of our error bounds no explicit rescaling is necessary for matrices with larger norm. If we assume || AAllr F I U for some (unknown) upper bound U, the guarantee of (2) would simply change to involve a eU terms instead of e. The second condition bounds how much the matrices change over time. Again for simplicity, we assume a fixed upper bound a the difference at each time step, but the algorithms presented in this paper will be adaptive to changing gaps between A; and 1,+1, and will perform better when these gaps are small on average. By triangle inequality,
2, but in applications we typically have a < 1, meaning that the changes in the dynamic matrix are small relative to its Frobenius norm, If this is not the case, there is no hope to improve on the naive method of applying Hutchinson's 's estimator repeatedly to each A,.
Note on Matrix Functions. In many applications A; = [(B:) for a dynamically changing matrix B. While we may have || - A:ll F (B 1) [(B.)]
>> |B++1 B:Ilr for functions like the matrix exponential, this is not an immediate issue. To improve on Hutchinson's 's estimator, the important requirement is simply that ||Ar+1 A& is small in comparison to || As discussed in Section 5 this is typically the case for application involving matrix functions.
We will measure the complexity of any algorithm for solving Problem in the matrix-vector multiplication oracle model of computation, meaning that we consider the cost of matrix-vecte vector products (which are the only way Ai...., Am can be accessed) to be significantly larger than other computational costs. We thus seek solely to minimize the number of such products used (3 The matrix-vector oracle model has seen growing interest in recent years as it generalizes both the matrix sketching and Krylov subspace models in linear algebra, naturally captures the true computational cost of algorithms in these classes, and is amenable to proving strong lower-bounds [37 6].

---

1.3 Main Result
more efficiently than Hutchinson's 's estimator: Theorem 1.1. For any 6, E (0, 1), the DeltaShift algorithm (Algorithm) 7 solves Problem 7 with
Our main result an algorithm for solving Problem |
o (m-
log(110)
log(1/0)
?
2²
total matrix-vector multiplications involving A Am
For large m, the first term dominates the complexity in Theorem 1.1 For comparison, a tight analysis of Hutchinson's 's estimator [29] establishes that the naive approach requires O (m- (lo((/))//e)
which is worse than Theorem 1.1 by factor of A natural setting is when O(e), in which (e) mantiiiimullpppiaatinnn on average over m time steps, in comparison (Tlog(1/6)/e for Hutchinson's 's estimator, quadratic improvement in
case Algorithm| I requires Olog( (o((1//))))))))))
To prove Theorem I.I we introduce a dynamic variance reduction scheme. By linearity of trace, tr(A,,1) = tr(A,) + tr (A.). where A; = A++1 - A,. Instead of directly estimating tr(A,-1), we combine previous estimate for tr,A.) with an estimate for tr(A.). computed via Hutchinson's 's estima- tor. Each sample for Hutchinson's is estimator applied to A; requires just two matrix-vector multiplies: one with A; and one with A,+1- At the same time, when A, has small Frobenius norm (bounded by a), we can estimate its trace more accurately than tr( ((i++11) While intuitive, this approach requires naive implementation, error in estimating tr(Ai). tr(22), compounds over time, eliminating any computational savings. To avoid this issue, we introduce a novel damping strategy that actually estimates tr(A:+1 - (1 (1))) for a positive damping factor 7.
care to make work.
We compliment our main result with a nearly matching conditional lower bound: in Section 4 we argue that our DeltaShift method cannot be improved in the dynamic setting unless Hutchinson's estimator can be improved in the static setting. We also present an improvement to DeltaShift under more stringent, but commonly present, bounds on A AI....,
....,,A,, Am and each A,+1 A; than Problem
1.4 Related Work
Prior work on implicit trace estimation and applications in machine learning is discussed in the beginning of this section. While there are no other methods that improve on Hutchinson's estimator in the dynamic setting, the idea of variance reduction has found applications in other work on implicit trace estimation 12 26 30). In these results, the trace of matrix is estimated by decomposing A = B + A where B has an easily computed trace (e.g., because it is low-rank) and || < ||A||,, so tr A) more easily approximated with Hutchinson's estimator than tr(A) directly.
2 Preliminaries
Notation. We let B mg denote real-valued matrix with rows and & columns. E R" denotes a real-valued vector with n entries. Subscripts like B, or æ, typically denote matrix or vector in sequence, but we use double subscripts with matrices to denote entries: Bij being the entry at the 5th row and 7th column. Let oz () denote the hh singular value of B. || B[[ denotes the Frobenius norm of B, V E = Ee (B) [[]] denotes the nuclear norm, E, (B(). We
let E[u] and Varlu) denote the expectation and variance of random variable v.
Hutchinson's Estimator. Our algorithm uses Hutchinson's trace estimator with Rademacher #1 random variables as a subroutine. Specifically, let ....... IR" be independent random vectors, with each entry +1 or -1 with probability 1/2. Let A E R=xa, Hutchinson' 's estimator for tr(A) is:
8 he(A) Σ of (Ag:)
(3)
While we consider the general, unstructured problem setting, we note that, if A, additional structure, is not necessarily easier to estimate the trace of A, than that of A.. For example, if A, is a PSD matrix then specialized trace estimation algorithms that improve on Hutchinson's 's method can be used [30]. Understanding dynamic traces estimation methods for sequences of structured matrices natural direction for future work.
3

---

Fact 2.1 (Hutchinson's ' expectation and variance). For any positive integer i and matrix A we have:
n ||
[[h((A)] = tr(),,
Var[he(A)]] =
2 S |
i=1
[3. Similar bounds hold when Hutchinson's s estimator is implemented with different random variables. For example, random Gaussians also lead to variance bound of & ||| However, Rademachers tend to work better empirically, Given Fact
Fact 2.1 follows from simple calculations, found e.g. in
2.T Chebyshev's inequality immediately implies a concentration bound for Hutchinson's estimator.
Fact 2.2 (Chebyshev's Inequality). For a random variable X with mean [[X] = f and variance Var [X] = for any k 2 1, P ([X μ4| 2 ko) S 1/k².
Claim 2.3. For any e, E (0, 1). if 2 then Pr [[he(A) tr(A)) 2 <||A|]]] 6.
The dependence in Claim 2.3 can improved from 8 to log(1/5) via the Hanson-
Hanson-Wright Wright inequality, which shows that he(4A) is sub- exbooneepnnial random variable 80 BS]. We also require Hanson- Wright to obtain our bound involving log(1/6). From this tighter result, Hutchinson's 's yields a total matrix-vector multiplication bound of (( log(1/5)/e) for solving Problem | by simply applying the estimator in sequence to AI,..., A....,Am.. Am-
Main Algorithmic Result
As discussed in Section 1.2 natural idea for solving Problem with fewer than (m)/²²) queries
is to take advantage of the small differences between A141 and A; to compute a running estimate of the trace. In particular, instead of estimating tr tr(A)), (r(2))...... tr(An) individually using Hutchinson's 's estimator, we denote A, AA,-- A and use lincarity of the trace to write:
tr tr(4 = tr) tr(A.).
(4)
By choosing large lo, we can compute an accurate approximation he, (A1) to tr(A)). Then, for j > 1 and l < Eo. we can approximate tr(A,) via the following unbiased estimator:
tr(A,) x heg (A1) + he(A,)
(5)
Since |Aill <a < ||A,-I||r, we expect to approximate tr(A2), ..., (r(A)) much more accu-
rately than tr (r(A))....,,t ir(Am) directly. At the same time, the estimator in (5) only incurs a 2 factor overhead in metrix-vector multiplies in comparisons to Hutchinson's it requires 2 (m 1) to
compute he(A))..... he,(Am)))) versus (m 1)e to compute he(A2). (Am)). The cost of the initial estimate hgg (At) is necessarily higher, but can be amortized over time.
3.1 Our Approach
While intuitive, the problem with the approach above is that error compounds due to the sum in Each he (A,) is roughly a/ a/vl away from tr( (r(A), after steps we naively expect total error Olj O/j-a/vo Vt). We can do slightly better by arguing that, due to their random nature, error actually accumulates O(vī a/ (E), but regardless, there is accumulation. One option is "restart" the estimation process: after some number of steps q, throw out all previous trace approximations, compute accurate estimate for tr(Ag), and for j > q construct estimator based tr(A,) tr(A,) + E tr (A)) While possible to analyze thcoretically, this approach turns out to be difficult to implement practice due several competing parameters (see details in Section
5
Instead, we introduce a more effective approach based on a damped variance reduction strategy, which call DeltaShift. Instead of being based (4). DeltaShift
which is detailed in Algorithm [] uses the following recursive identity involving fixed parameter 0 S < (to be chosen later):
tr( (r,,)) = (1 - t) + tr(A)), where A = A₁ - (1 7)A,
(6)

---

Algorithm DeltaShift
Input: Implicit matrix-vector multiplication access to A1,..., Am R positive integers to, damping factor E [0, 1]. Output: , b approximating tr(Ai), tr(Am).
Initialize 1 < Eu Σ 1 of A19i, where 91...... I6o E R" are random #1 vectors
for j 2 to m do Draw & random 11 vectors 91, . E Rd 2t A 1-1 w1 t- Ag91,..., W f- A,9e
21 -- 191
t + (1 7) + E of. (w:- - (1 7)2))
end for
Given an estimate 1--1 for tr(Aj-1), DeltaShift estimates tr(A;) by (1-7)) (1- he(A,) This approach has several useful properties: 1) if t is an unbiased estimate for tr(Aj-1), tj is an unbiased estimate for tr(A,), 2) ||رر e is not much larger than ||ررر y small, and 3) by shrinking f--1 by factor of (1 - 7) when computing ig, we reduces the variance of this leading term. The last property ensures that error does not accumulate over time, leading to our main result:
Theorem 1.1 (Restated). For any 6, (0,1), Algorithm | run with & a, to O (log(1/o)fe)) and & = O (a log(1/8)/{) solves Problem 1. In total, it requires
o m log(1/0)
log(1/6)
62
matrix-vector multipicaation with A....., Am
The full proof of Theorem I.I relies on the Hanson- Hanson-Wright inequality, and is given in Appendix B Here, we give simple proof of essentially the same statement, but with a slightly weaker dependence on the failure probability
Proof. Let to 2 and é == 8a The proof is based inductive analysis of the variance of lj, the algorithms estimate for ir(A)). Specifically, we claim that that for j 1,..., m.
23
Var[t,] S 86?.
For the base case, = 1, 7 follows directly from Fact 2.1 because i is simply Hutchinson's estimator applied to A, and A1lIF S 1. For the inductive case, L is the sum of two independent estimators, 1--1 and hz (,) So, to bound its variance, we just need to bound the variance of these two terms. To address the second, note that by triangle inequality, || || ||A; - (1-1111444 (1- S ||A, Ay-1llr - 1|r 2a. Thus, by Fact 2.1 Var A, I a Combined with the inductive assumption that Var < 8c2, we have:
80²
Var[t,] (1 (Ll-7)2 Varlhe (1 (1- a)
S (1 adde² + aade de?.
This proves 7). and by Chebyshev's inequality we thus have PIt tr(A,)) 2 d S 8 for all j.
3.2 Selecting y in Practice
While DeltaShift is simple to implement, in practice, its performance sensitive to the choice of y- For the Theorem I.I analysis, we assume : a, but may not be known apriori, and may change over time. To address this issue, we describe a way to select a near optimal at each time step j (the choice may vary over time) with very little additional computational overhead. Let 1/ be the variance of our estimator for tr( (r(Aj-)). We have that vy (1 - -)) + Var[he[A, (1
V
1)4}-1)] < (1 - u3 + 11, (1 7) 11 At time step j. natural goal is to choose damping parameter 7 that minimizes this upper bound on the variance of ty:
= arg min (1 (1- 7
2 11
(8)
7
5

---

where D, A, (1 (111-) as before. While 8 cannot be computed directly, observing that |B tr(B² B) for any matrix B, the above quantity can be estimated as ëg = (1 - --) hhe (A where 8/ー is an estimate for U5-1. The estimate he ( A5) be computed using exactly the same & matrix-vector products with A₁ and Aj-1 that are used to estimate tr tr(A,): so there is little computational overhead. Moreover, since Σz A is positive semidefinite, long as 2 log (1)). we will obtain a relative error approximation its trace with probability 1- B
An alternative approach to estimating 1y would be to simply compute the empirical variance of the average he (A,), but this requires fixing 7, An advantage of our closed form approximation is that it can be used analytically optimize y. Specifically, expanding AJ Az. j. we have that:
2 8 = (1- / (hr(AfA))+ (L A,) - 2(1 - ())n))))((A)31)
(A)
(9)
Above, each estimate he is understood to use the same set of random vectors. Taking the derivative and setting zero, have that the minimizer of 9 denoted 4", equals:
2hp (AT A4,)
1
(10)
tö 2h
(47 (
A4-14
This formula for 7" motivates an essentially parameter free version of DeltaShift, which is used in our experimental evaluation (Algorithm 2 in Appendix A). The only input to the algorithm is the
number of matrix-vector multiplies used at each time step, l. For simplicity, unlike Algorithm []] we
do not use a larger number of matrix- multiplies when estimating A1. This leads to somewhat higher error for the first matrices the sequence A1,.... An. error quickly falls for large j.
Algorithm Improvements and Lower Bound
In this section, we prove lower bound showing that, in general, Theorem I.] is likely optimal. On the other hand, we show that, if we make a slightly stronger assumption on A1.-.., An and the dynamic updates A, A4+1 A,, an improvement on DeltaShift possible.
4.1 Lower Bound
As noted, for large number of time steps m, the matrix-vector multiplication complexity of DeltaShift dominated by the leading term Theorem 1.1, 0(m 0 log(1/))/e). We show that is unlikely an improvement on this term can be obtained in general:
Lemma 4.1. Suppose there is an algorithm S that solves Prob. with o(m a log(1 (1/))//))))))))))))) total matrix-vector multiplies with A2,..., Am and any number of matrix- multiplies with A when a = 1/mm - 1). Then there an algorithm T that achieves (T) for a single with o(log(1 (1/8) matrix-vector multiplies.
Proof. The proof is via a direct reduction. Given matrix A, positive integer >> 1, and parameter T- construct the sequence of matrices:
A O,
A a=
A1/o (1 (1- a)
Am A
Since A= O, and every A2...., Am is scaling of A, any algorithm S satisfying the assumption of matrix multiplications with A. Moreover, if S is run on this sequence of matrices, on the last step it outputs an approximation tm to (rA)) with Pr|itm- tr(A) 2 c) <
Lemma 4.1 can be implemented with o (mm - 1) a log(1/6)/e²) = ( =(log(1/8)/e²)
So algorithm can simply simulate
S on A Am and return its final estimate to satisfy
Lemma is conditional lower-bound on matrix- query algorithms for solving Problem
Hutchinson's is estimator cannot improved for static trace estimation (and it hasn't 't been for 30 years) then DeltaShift cannot be improved for dynamic trace estimation. We believe the bound could
be made unconditional through slight generalization of existing lower bounds on trace estimation in the matrix-vecee multiplication model [30 46).

---

4.2 Improved Algorithm
recent improvement on Hutchinson's 's estimator, called Hutch++, was described in [3). For the static trace estimation problem, Hutch++ achieves matrix-vector multiplication complexity of 0(1/e) to compute a relative enror (1 + e) approximation to the trace of any positive semi-definite matrix (PSD), improving on the O required by Hutchinson's s. It does so via a variance reduction method (also used e.g. in [2]) which allocates some metrix-vector products to a randomized SVD algorithm which approximates the top singular vector subspace of A. This approximate subspace is projected off of and Hutchinson's 's used to estimate the trace of the remainder.
In our setting it is not realistic to assume PSD matrices while in many applications AI..... Am are all PSD, it rarely the case the Ai..... Am-1 are. Nevertheless, we can take advantage of more general bound proven in [30] for any matrix:
Fact 4.2 (Hutch++ expectation and variance). Let hi (A) be the + estimator of 330 applied any matrix A with l matrix-vector multiplications. We have:
16 8 |A
Eiht (A)] = tr(A),
Varlh Varlh)((A)]
Recall that |. denotes the nuclear norm of A. Comparing to the variance of Hutchinson's 's estimator from Fact notice that the variance of Hutch++ depends on zr instead of implying faster convergence the number of matrix-vector products, f, increases. A trade-off is that the variance scales with || instead of ||
- || is strictly larger, and possible a factor of n larger than || However, for matrices that are rank k. ||A||| &||lA|??. so the norms are typically much closer for low-rank or nearly low-rank matrices. In many problems, A..... Am may have low-rank structure, in which case, alternative based on Hutch++ provides better performance.
Formally, we introduce a new variant of Problem 1 to capture this potential improvement.
Problem 2 (Dynamic trace estimation w/ Nuclear norm assumption). Let A1,..., Am satisfy:
1. || < for all i m1. 2. |A++1 Asll- S a, for all i E [1, 1].
Given matrix-vector multiplicatiin access to each A; in sequence, the goal is to compute trace approximations (1...., Im for tr(Ai), br(Am) such that, for all i, P[4 tr((:)) 2 c] 6.
In Appendix C we prove the following result on a variant of DeltaShift that we call DeltaShift++: Theorem 4.3. For any e, 6, a (0, DeltaShift++ (Algorithm 4 solves Problem 2 with
o
m Val8
V15
total metrox-vector multiplications involving A Am
Theorem is stronger than Theorem 1.I for vanilla DeltaShift in that it has linear instead of quadratic dependence on e. In particular, its leading term scales va/e², whereas Thcorem [.] scaled with az/e.. However, the result does require stronger assumptions
A An and each A; = A;+1 - A; in that Problem 2 requires these matrices to have bounded nuclear norm instead of matrix is strictly larger than its Frobenius norm, these requirements are stronger than those of Problem [ As we will show in Section 5] the benefit of improved dependence often outweights the cost of these more stringent assumptions.
of Frobenius norm. Since the nuclear
Experiments
We show that proposed algorithm outperforms three alternatives both synthetic and real-world trace estimation problems. Specifically, we evaluate the following methods:
Hutchinson. The naive method of estimating each tr(Ai), .(,A)),...,t(.... tr(Am) using an independent Hutchinson's estimator, as discussed in Section 2
NoRestart. The estimator of 55) which uses the same variance reduction strategy as DeltaShift for all j 2 2, but does not restart or add damping to reduce error accumulation.
7

---

Restart. The estimator discussed in Sec. 3.1 which periodically restarts the variance reduction strategy, using Hutchinson's 's to obtain fresh estimate for tr(A,). Pseudocode is in Appendix A DeltaShift. Our parameter free, damped variance reduction estimator detailed Appendix [A
We allocated fixed number of matrix-vecte queries, Q. to be used over all time steps 1, ,,.... For Hutchinson and DeltaShift, the same number of vectors /m was used at each step. For Restart and NoRestart, the distribution was non-uniform, and parameter selections are described in Appendix D
Average error +0.25 100 trials
Average errar 4:0.55 for 100 trials
0.05
flaciinsor
0.03
Paschinion
Recial
Rodat
0.04
Naflntas
hoBaccar
MW
MMM
0.02
0.03
0.02
0.01
0.01
0.00
40
60
80
100
40
60
100
Time step (i)
Time step (6)
(a) Synthetic data with low perturbation
Synthetic data with significant perturbation
Figure I: Comparison of DeltaShift with Hutchinson and NoRestart on synthetic data with Q - 10
Average emoor +0,250 for 100 trials
Average error 450.25er for 100 trials
Hlactiion
Nothatart
Di
Map o maniing Astiad
|=AA1-1
0.04
0.02
0.02
0.01
80
Time step (1)
Time step (+)
(a) Synthetic data with Q 2 10
(b) Synthetic data with Q 10
Figure 2: Comparison of DeltaShift with Hutchinson and NoRestart on synthetic data.
Synthetic data: To simulate the dynamic setting, we generate random matrix A E Rax= and add random perturbations for each of 100 time steps. We consider two cases: low (Fig. I(a) and significant (Fig. [(b) perturbations, the exact details of which, as well as the allocation of matrix-vector products for NoRestart and Restart, are discussed in Appendix [D] We report scaled absolute error between the estimator at time, ty. and the true trace tr(A,). As expected, Hutchinson is outperformed even by NoRestart when perturbations are small. DeltaShift performs best, and its error actually improve slightly over time. DeltaShift also performs best for the large perturbation experiment. We note that choosing the multiple parameters for the Restart method was a challenge in comparison to DeltaShift. Tuning the method becomes infeasible for larger experiments, exclude this method in our other experiments. That includes for the plots in Fig. [2] which show that DeltaShift continues to outperform Hutchinson and NoRestart for lower values of Q.
Counting triangles: Our first real-data experiment is on counting triangles a dynamic unweighted, undirected graph G via the fact that the number of triangles equals t Lir(B)), where B is the adjacency matrix. The graph dataset we use is the Wikipedia vote network dataset with 7115 nodes [25, 24]. At each timestep perturb the graph by adding random k-clique, for & chosen uniformly between 10 and 150. After 75 time steps, we start randomly deleting among the subgraphs added. We follow the same setup for number of mntrix-vector products used by the estimators and the error reported as in the synthetic matrix-vector multiplications with B is 3Q, since each oracle call computes ((((x)), As seen in Fig. 3, DeltaShift provides the best estimates overall.
D). Note that for this particular application, the actual number of
Estimation natural connectivity: To evaluate the DeltaShift++ algorithm introduced in Section 4
we address an application in [43] on estimating natural connectivity in a dynamic graph, which is function of tr ((exp(B) for adjacency matrix B. This problem has also been explored in 7 4]. We

---

Average error 25c for 50 trials
Average error =0.25a for 50 trials
Passhinn
Nifintat
0.18
Detrotee
0.03
0.06
0.02
Anmmy
0.04
0.01
40 60
100
20
40
60
Time step (1)
Time step
Graph data with Q
10
Graph data with Q 10*
Figure 3: Comparison of DeltaShift with Hutchinson and NoRestart for triangle counting experiment.
Average emror +5. for 100 trials
Average 100 timesteps for 100 trials
Toral matikuusttor products (Q)
Dendde
330100
Nothaatart
Dulu-Saii
90000
Irrinn (431-4)
mix 12 Peepp (A.))
23000
0.0006
20000
150000
0.0004
10000
0.010
Time
Average relative error
Figure 4: Error of DeltaShift and DeltaShift++ for dynamic estimation of natural connectivity.
Figure 5: Average error computational cost for synthetic data with large perturbations.
use the road network data Gleich/minnesota (available at minnesota). We perturb the graph over time by choosing nodes at random and adding an edge between them, and use the Lanczos method to approximate matrix-vector products with exp( B). We find that DeltaShift++ performs better than DeltaShift (Fig. 4 the change in exp( B) tends to be nearly low-rank, and thus have small nuclear norm (see [4] for details). Both DeltaShift and DeltaShift++ perform significantly better than naive Hutchinson's 's when 100 matrix-vector products are used per time step.
tamu. edu Gleich/
The key takeaway from the experiments above is that DeltaShift and DeltaShift++ are able obtain good dynamic trace approximations in far fewer matrix matrix-vector products compared to Hutchinson's and other methods, resulting in considerable computational savings. This is made evident in Figure 5 which plots average relative eтror across all time steps vs. total number matrix-vector products( for various values of Q. In order achieve the accuracy level as DeltaShift, Hutchinson's 's requires substantially more matrix-vector products.
Hessian spectral density: Finally, we evaluate the performance of DeltaShift on an application pertaining to dynamically changing Hessian matrix, H. involved in training a neural network. As discussed in Section[] common goal is to approximate the spectral density of !.. Most methods for doing so, like the popular Kernel Polynomial Method [45], require computing the trace of polynomials of the matrix H. We consider the sequence of Chebyshev polynomials To.... ..... Ty. and estimate tr(To (H),... .............. tr( (H)) Other polynomial basis sets can also be used (e.g., Legendre polynomials). Experimental details are discussed in section Appendix D but summarized the results here. We implement matrix vector products with H using the PyHessian library
and report
Table Average relative error for trace of Chebyshev polynomials of Hessian.
Table 2:
Average relative error for trace of Chebyshev polynomials of Hessian.
Hutchinson
NoRestart
DeltaShift
Hutchinson
NoRestart
DeltaShift
(H))
2.5e-02
3.7e-02
1.7e-02
[(H))
1.9e-02
5.0e-02
1.5e-02
(2(H)
1.2e-06
1.7e-06
8.0e-07
Te(H)
1.2e-06
2.9e-06
9.9e-07
(H(
4.0e-02
4.1e-02
3.1e-02
Ts(H)
7.7e-02
9.4e-02
6.le-02
1.Se-06
1.7e-06
1.0e-06
T((H)
1.7e-06
2.8e-06
1.5e-06
(H(
2.1e-02
4.3c-02
1.9c-02
Ts(H)
2.1e-02
4.2c-02
1.8c-02

---

average eгrоr over 25 training epochs for the Hessian of ResNet model with 269722 parameters trained it on the CIFAR-I0 10 dataset. As impossible to compute the true trace of these matrices. we use Hutchinson's s estimator with greater number of queries as placeholder for ground-truth, and compare the performance against the computed values. As can be seen in Tables [ and DeltaShift obtains uniformly better approximation to the trace values, although the improvement is small. This makes sense, as more progress on each training epoch implies a greater change in the Hessian over time, meaning is larger and thus DeltaShift's ' advantage over Hutchinson's ' is smaller.
Acknowledgements
We would like to think Cameron Musco for helpful discussions, as well as the paper referees for detailed feedback. This work was supported by NSF Award #2045590.
References
[1] Ryan P. Adams, Jeffrey Pennington, Matthew J. Johnson, Jamie Smith, Yaniv Ovadia, Brian Patton, and James Saunderson. Estimating the spectral density large implicit matrices.
arXiv:I802.00 1802.03457 2018.
[2] Haim Avron. Counting triangles in large graphs using randomized matrix trace estimation. In Proceedings of the 16th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2010.
[3] Haim Avron and Sivan Toledo. Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. Journal of the ACM, 58(2), 2011.
[4] Bernhard Beckermann, Daniel Kressner, and Marcel Schweitzer. Low-rank updates of matrix functions. SIAM Journal on Matrix Analysis and Applications, 39(1 39(1) ):539-555,, 2018.
[5] Christos Boutsidis, Petros Drineas, Prabhanjan Kambadur, and Anastasios Zouzias. A random- ized algorithm for approximating the log determinant of a symmetric positive definite matrix. Linear Algebra and its Applications, 533, 03 2015.
[6] Mark Braverman, Elad Hazan, Max Simchowitz, and Blake Woodworth. The gradient com- plexity of linear regression. In Proceedings of the 33rd Annual Conference on Computational Learning Theory (COLT), volume 125, pages 627-647 2020.
[7] Hau Chan, Leman Akoglu, and Hanghang Tong. Make it or break it: Manipulating robustness in large networks. In Proceedings of the 2014 SIAM International Conference on Data Mining, pages 325-333. SIAM, 2014.
[8] Chen Chen, Ruiyue Peng, Lei Ying, and Hanghang Tong. Network connectivity optimization: Fundamental limits and effective algorithms. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 1167-1176, 2018.
[9] Edoardo Di Napoli, Eric Polizzi, and Yousef Saad. Efficient estimation of eigenvalue counts in an interval. Numerical Linear Algebra with Applications, 2016.
[10] Zhen Dong, Zhewei Yao, Yaohui Cai, Daiyaan Arfeen, Amir Gholami, Michael W Mahoney, and Kurt Keutzer. Hawq-v2: Hessian aware trace-weighted quantization of neural networks. Advances in Neural Information Processing Systems 33 (NeurIPS), 2020.
[11] Ernesto Estrada and Naomichi Hatano. Communicability in complex networks. Phys. Rev. E. 77:036111, 1, Mar 2008
[12] Arjun Singh Gambhir, Andreas Stathopoulos, and Kostas Orginos. Deflation as method of variance reduction for estimating the trace of matrix inverse. SIAM Journal on Scientific Computing, 39(2):553--A558, 2017.
[13] Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao. An investigation into neural net opti- mization via hessian eigenvalue density. In Proceedings of the 36th International Conference on Machine Learning (ICML), volume 97, pages 2232-2241 2019.
10

---

[14] Didier Girard. Un algorithme simple et rapide pour la validation croisee géenéralisée sur des problémes de grande taille. Technical report, 1987.
[15] Insu Han, Prabhanjan Kambadur, Kyoungsoo Park, and Jinwoo Shin. Faster greedy map inference for determinantal point processes. In International Conference on Machine Learning, pages 1384-1399 1393. PMLR, 2017.
[16] Insu Han, Dmitry Maliouuov, Haim Avron, and Jinwoo Shin. Approximating the spectral sums of large-scale matrices using stochastic Chebyshev approximations. SIAM Journal on Scientific Computing, 2017.
[17] Insu Han, Dmitry Malioutov, and Jinwoo Shin. Large-scale log-determinant computation through stochastic Chebyshev expansions. In Proceedings of the 32nd International Conference on Machine Learning (ICML), pages
2015.
[18] D. L. Hanson and F. T. Wright, A bound on tail probabilities for quadratic forms in independent random variables. Ann. Math. Statist., 1079:1307-- 06 1971.
[19] Yuanyang Zhu Hanyu Li. Randomized block Krylov space methods for trace and log- determinant estimators. arXiv:2003.00022) 2020.
[20] 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.
[21] Nicholas J. Higham. Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics, 2008.
[22] Michael F. Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in Statistis-Simmllation and Computation, ((2)):333455) 1990.
[23] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
[24] Jure Leskovec, Daniel Huttenlocher, and Jon Kleinberg. Predicting positive and negative links in online social networks. In Proceedings of the 19th International Conference on World Wide Web, WWW '10, pages 641-650, New York, NY, USA, 2010. Association for Computing Machinery.
[25] Jure Leskovec, Daniel Huttenlocher, and Jon Kleinberg. Signed Networks in Social Media, pages 1361-1370. Association for Computing Machinery, New York, NY, USA, 2010.
[26] Lin Lin. Randomized estimation of spectral densities of large matrices made accurate. Nu- merische Mathematik, 136(1):183-215 2017.
[27] Lin Lin, Yousef Saad, and Chao Yang. Approximating spectral densities of large matrices. SIAM Review, 58(1):34 65, 2016.
[28] Michael Mahoney and Charles Martin. Traditional and heavy tailed self regularization in neural network models. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning (ICML), volume 97, pages 4284-4293, 2019.
[29] Per-Gunnar Martinsson and Joel A. Tropp. Randomized numerical linear algebra: Foundations and algorithms. Acta Numerica, 29:403-572, 2020.
[30] Raphael A. Meyer, Cameron Musco, Christopher Musco, and David Woodruff. Hutch++: opti- mal stochastic trace estimation. Proceedings of the 4th Symposium on Simplicity in Algorithms (SOSA), 2021.
[31] Cameron Musco, Praneeth Netrapalli, Aaron Sidford, Shashanka Ubaru, and David P. Woodruff. Spectrum approximation beyond fast matrix multiplication: Algorithms and hardness. Proceed- ings of the 9th Conference on Innovations in Theoretical Computer Science (ITCS), 2018.
11

---

[32] Barak A. Pearlmutter. Fast exact multiplication by the hessian. Neural computation, 6(1):147- 160, 1994.
[33] Jeffrey Pennington, Samuel Schoenholz, and Surya Ganguli. The emergence of spectral universality in deep networks. In Proceedings of the 21 International Conference on Artificial Intelligence and Statistics (AISTATS), pages 1924-192999 1932, 2018.
[34] Xu Qian, Victor Li, and Crews Darren. Channel-wise hessian aware trace-weighted quantization of neural networks. arXiv preprint arXiv:2008. 08284, 2020.
[35] Mark Rudelson, Roman Vershynin, et al. -tright inequality and sub-Gaussian concen- tration. Electronic Communications in Probability, 18, 2013.
[36] Arvind K. Saibaba, Alen Alexanderian, and Ilse C. F. Ipsen. Randomized matrix-frie trace and log-determinant estimators. Numerische Mathematik, 137(2)3353-33 395. 2017.
[37] Max Simchowitz, Ahmed EI Alaoui, and Benjamin Recht. Tight query complexity lower bounds for PCA via finite sample deformed Wigner law. In Proceedings of the 50th Annual ACM Symposium on Theory of Computing (STOC), pages 1249-1259, 1259, 2018.
[38] Andreas Stathopoulos, Jesse Laeuchli, and Kostas Orginos. Hierarchical probing for estimating the trace of the matrix inverse on toroidal lattices. SIAM Journal on Scientific Computing, 35(5):S299-S322 2013.
[39] Xiaoming Sun, David P. Woodruff, Guang Yang, and Jialin Zhang. Querying a matrix through matrix-vector products. In Proceedings of the 46th International Colloquium on Automata, Languages and Programming (ICALP), volume 132, pages 19,11,,, 2019.
[40] Jok M. Tang and Yousef Saad. Domain momnain-decomposttion--yyy methods for computing the diagonal of a matrix inverse. SIAM Journal on Scientific Computing, 3365):2823-2877, 2011.
[41] Shashanka Ubaru and Yousef Saad. Applications of trace estimation techniques. In High Performance Computing in Science and Engineering, pages 19-33, 2018.
[42] Martin J. Wainwright. Basic tail and concentration bounds, pages 21-57. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2019.
[43] Sheng Wang, Yuan Sun, Christopher Musco, and Zhifeng Bao. Public transport planning: When transit network connectivity meets commuting demand. In Proceedings of the 2021 ACM SIGMOD International Conference on Management of Data, 2021.
[44] Mingwei Wei and David J Schwab. How noise affects the hessian spectrum in overparameterized neural networks. arXiv preprint arXiv:1910.00195 2019.
[45] Alexander WeiBe, Gerhard Wellein, Andreas Alvermann, and Holger Fehske. The kernel polynomial method. Reviews of modern physics, 78(1):275, 2006.
[46] Karl Wimmer, Yi Wu, and Peng Zhang. Optimal query complexity for estimating the trace of a matrix. In Proceedings of the s International Colloquium on Automata, Languages and Programming (ICALP), pages 1051-1062, 2014.
[47] Zhewei Yao, Amir Gholami, Kurt Keutzer, and Michael Mahoney. Pyhessian: Neural networks through the lens of the hessian. arXiv preprint arXiv: 2019.
[48] Zhewei Yao, Amir Gholami, Kurt Keutzer, and Michael W. Mahoney. Pyhessian: Neural networks through the lens of the hessian. IEEE BigData, 2020.
[49] Zhewei Yao, Amir Gholami, Sheng Shen, Kurt Keutzer, and Michael W Mahoney. ADAHES- SIAN: An adaptive second order optimizer for machine learning. AAAI Conference on Artificial Intelligence (AAAI 202/), 2021.
12


Reviews:
[{"review_id": "qlhrOYonj_H", "reviewer": "Reviewer_jK9X", "summary": "The authors proposed a method for estimating the trace for a sequence of matrices that only changes for a small amount at each step. The method originates from the well-known stochastic trace estimator, but bring it to a new practical setting (e.g. Hessian at learning iterations, evolving network, etc). By an interesting treatment of the iterative process (i.e. damping), the resulting algorithm could cut the computational cost by a factor of the size of the tolerable error. In addition, the authors included the recent development/improvement on the original Hutchinson estimator, to get better versions of their algorithm under more assumptions.", "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": "Please address this and correct me if I'm wrong, in the inequality between line 175 and 176, $(1-\\alpha)^2\\delta\\epsilon^2+4\\alpha/\\ell\\to\\frac{3}{2}\\delta\\epsilon^2$ as $\\alpha\\to 0$. Thus the inequality doesn't seem to hold here.\n\nOriginality: I think the work is novel and this dynamic setting is probably never looked at before. The work builds on top of the well-known Hutchinson estimator, but the main contribution has great ideas that are totally independent from previous work. The paper also has good reference and sourced all the relevant work.\n\nQuality: The paper is technically sound for the most part, but please respond about my question at the top. Nonetheless the experiment results are promising and shows the algorithm works well in practice.\n\nClarity: The paper is well-written and easy to follow. The list of reference is good. One thing bothers me a little is that the authors claimed in multiple places that $tr(\\Delta_j)$ is easier to estimate and we can we can approximate $tr(\\Delta_j)$ more accurately than $tr(A_j)$. While I get authors' point in the context of this problem, the statement on its own isn't very accurate because we usually look at relative error in this kind of estimations, and $tr(\\Delta_j)$ isn't obviously easier in that sense.\n\nSignificance: I think the work is important and have many applications as suggested by authors' experiments. While I'm not sure about how much scalability problem we have in practice given the efficiency of the original estimator, the proposed work makes clear improvements. \n\nOther comment: While damping seems a great idea in general, the authors also mentioned that sometimes $\\Delta_j$ can have special properties (e.g. low-rank, diagonally concentrated after reordering, etc). In this case, the damping would destroy those structures. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There is no negative societal impact of this work.", "review_text": "", "extra_scores": {"ethical_concerns": "There is no ethical concerns from the content of this paper.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "kprcik18Ryu", "reviewer": "Reviewer_uAQs", "summary": "The authors in this paper study the problem of estimating the trace of a matrix that is changing dynamically over several iterations. Under the assumption that the diagonal entries of the matrix are not directly accessible, but multiplication with a vector $Av$ (or the quadratic form $v^{\\top}Av$), the Hutchinston estimator is the baseline algorithm for comparison. The Hutchinston estimator takes the average of $\\epsilon^{-2} \\log 1/\\delta,$ estimates $g_i^{\\top} A g_i$ where $g_i$ is a random vector with iid gaussian variables, or Rademacher ($\\pm 1$) variables, and is correct with probability $1-\\delta.$\n\nThe model is that we have a sequence of matrices $A_1, \\ldots, A_m$ each with $\\|| A_i \\|| \\le 1$ and $\\||A_{i+1} - A_{i} \\|| \\le \\alpha,$ (where the matrix norm used is the Frobenius norm) and we wish to output estimates of $tr(A_1),\\ldots,tr(A_m)$ such that each estimate is accurate to $\\pm \\epsilon.$ Using a Hutchinston estimator at each iteration will require $m \\epsilon^{-2} \\log m/\\delta$  where all the answers are accurate with probability $1-\\delta.$\n\nThe paper presents an algorithm that is able to maintain accurate estimates for such a dynamically changing $A$ using only $(m \\alpha + 1 )\\epsilon^{-2} \\log m/\\delta,$ which is essentially better by a factor of $\\alpha.$ The authors also show a lower bound claiming to show that this method cannot be improved unless the original Hutchinston method is improved.\n\nThe algorithms extend their algorithm to the case where you have a stronger nuclear norm bound assumption, and present numerical simulations demonstrating their advantages.", "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": "I think this is a nice algorithm based on a clean and simple idea that could lead to direct gains in some settings, e.g. estimating the trace of the hessian while training a neural network. The key idea is to maintain a trace estimate of a moving (exponentially) weighted average of the matrices, and to estimate the trace of the difference at each step using the Hutchinston estimator.\n\nA few comments:\n1. A nicer baseline to compare against is a Hutchinston estimator that is restarted every $1/\\alpha$ iterations.\n2. The lower bound is weaker than the claim that \"this method is the best possible\". The lower bound only shows that using the case $\\alpha \\approx 1/m,$ this problem recovers the original estimation problem for hutchinston estimator. It does not rule out running time bounds such as $m + \\sqrt{m \\alpha} \\epsilon^{-2} \\log m/\\delta.$\n3. I think your algorithm can also work in the weaker model where we only have access to the quadratic forms. I think that's the usual model for the Hutchinston estimator too.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3-4"}, "extra_text": {}}, {"review_id": "dNYr16jozAd", "reviewer": "Reviewer_hb2j", "summary": "This work considers the problem of dynamic trace estimation. Specifically, if there is a sequence of matrices $A_1, \\dots, A_n$, and the norm of $A_i - A_{i-1}$ is small relative to the norm of $A_i$, we can estimate the trace of each matrix efficiently. The paper gives an efficient algorithm, the convergence analysis and conducts relavant experiments.\n", "questions": "", "limitations": "", "rating": 6, "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": "The paper is well-organized and I don't see any major flaws in the proofs. The idea of damped variance reduction is quite good. But I worry that setting of this paper is too specific and is not very important for other types of problems.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "This method seems very specific to this problem. I don't know if it can be used in other time varying optimization problems, or would be useful for other machine learning tasks. The trace estimation problem has two special properties: The trace is linear, so we can easily divide it into subproblems; and the estimation error has a relative bound. If we have a time varying linear regression problem instead, it is hard to reasonably assume that the error can be bounded relatively rather than absolutely. Furthermore, although the problem is \"dynamic\", we are given the matrices $A_i$ and has the freedom to perform the measurement on the linear combinations of $A_i$ as we want. This is very different from most time series datasets where we can only observe data, not measure it.\n\nI not satisfied with Section 3.2, the theoretical justification of this adaptive step is lacking. This can be a main focus of the paper.\n\nOne immediate application of this dynamic trace estimation I can think of, is in the stepsize for gradient descent (quasi-Newton method). At each step, we can use the trace of the inverse Hessian (divided by dimension) as the stepsize. Then the difference of two Hessians is diminishing as the iterates get closer to the optimizer. If the difference is vanishing, how can we choose the damping parameter? And I guess the answer to this can give us a very good quasi-Newton method. It is also related to the previous problem. We must choose the damping parameter adaptively.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "9mcUOPC7C6p", "reviewer": "Reviewer_94L4", "summary": "This paper studies dynamic setting of trace estimation problem; compute the trace of a sequence of matrices efficiently. The authors consider the differences of consecutive matrices are small and utilize the linearity of trace and Hutchinson estimator, which returns trace values with only matrix-vector multiplications. In addition, in order to reduce the variance of trace estimations, they adopt a damping strategy; it controls the scale of the input matrix to estimate and reduces its variance appropriately. They propose a closed form of the optimal damping factor. Experimental results show the superiority of the proposed algorithm with various problems of trace estimation under both synthetic and real-world datasets.", "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": "This work proposes an efficient method for the dynamic trace estimation problem by utilizing well-known Hutchinson trace estimation and variance reduction technique (aka control variate). The problem is quite fundamental so that it can be applied to various of practical tasks, e.g., triangle counting, Hessian spectral density estimation and so on. The idea of Hutchinson estimator and variance reduction is not new and widely used in many problems involving computational issues, but the authors analyze with concrete theoretical results. \n\nI think that dynamic trace estimation is somewhat related to greedy algorithm of the determinant maximization problem (aka MAP inference of Determinantal Point Processes). Indeed, [1] tried to dynamically estimate log-determinant (polynomial version of Hutchinson estimator) and analyzed the variance of two distinct Hutchinson estimations when randomness is shared (see Theorem 2 in [1]). More generally, it seems that the proposed method can combine with a polynomial to extend to more various essential problems.\n\n[1] Han, Insu, et al. \"Faster greedy MAP inference for determinantal point processes.\" ICML, 2017.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The paper is quite well-written and easy to understand. But, it would be good to highlight concrete applications that I provide in the experiment section. Also, I suggest providing some details of technical terms (e.g., Hanson-Wright inequality in line 127) for non-expert readers. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper considers estimating the trace of a sequence of matrices $A_1, A_2, … A_t$ under the constraint that $||A_i - A_{i+1}||\_2$ is small. The proposed algorithm outperforms the baseline which repeatedly applies Hutchinson's stochastic trace estimator in terms of sample complexity. This paper receives unanimous support from the reviewers. Thus, I recommend acceptance. 
Decision: Accept (Poster)