Dynamic Trace Estimation
Year: 2021

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 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 a dynamically changing graph.

heading: 1 Introduction
text: 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 to an oracle for computing matrix-vector products Ax1, . . . , Axm between an n× n matrix A and chosen vectors x1, . . . , xm, the goal is to compute an approximation to A’s trace, tr(A) = ∑n i=1Aii. This problem arises when A’s diagonal entries cannot be accessed explicitly, usually because forming A is computationally prohibitive. As an example, consider A which is the Hessian matrix of a 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, if B is a graph adjacency matrix, tr(B3) equals six times the number of triangle in the graph [2]. ComputingA = B3 explicitly to evaluate the trace would require O(n3) time, while the matrix-vector multiplication Ax = B · (B · (Bx)) only requires O(n2) time. Similarly, in log-determinant approximation, useful in e.g. Bayesian log likelihood computation or determinantal point process (DPP) methods, we want to approximate the trace of A = log(B) [5, 17, 36]. Again, A takes O(n3) time to form explicitly, but Ax = log(B)x can be computed in roughly O(n2) time using iterative methods like the Lanczos algorithm [21]. Dynamic versions of the log-determinant estimation problem have been studied due to applications in greedy methods for DPP inference [15].
In data science and machine learning, other applications of implicit trace estimation include matrix norm and spectral sum estimation [16, 41, 31], as well as methods for eigenvalue counting [9] and spectral density estimation [45, 27]. Spectral density estimation methods typically use implicit trace estimation to estimate moments of a matrix’s eigenvalue distribution – i.e., tr(A), tr(A2), tr(A3), etc. – which can then be used to compute an approximation to that entire distribution. In deep learning, spectral density estimation is used to quickly analyze the spectra of weight matrices [33, 28] or to probe information about the Hessian matrix during optimization [13, 49]. Trace estimation has also been used for neural networks weight quantization [10, 34] and to understand training dynamics [44].
35th Conference on Neural Information Processing Systems (NeurIPS 2021).

heading: 1.1 Static Trace Estimation
text: The mostly widely used implicit trace estimation algorithm is Hutchinson’s estimator [14, 22]. Letting g1, . . . , g` ∈ Rn be random vectors with i.i.d. mean 0 and variance 1 entries (e.g., standard Gaussian or ±1 Rademachers), Hutchinson’s approximates tr(A) via the average h`(A) = 1` ∑1 i=1 g T i (Agi). This estimator requires ` matrix-vector multiplications to compute. Its variance can be shown to be O(‖A‖2F /`) and with high probability, when ` = O(1/ 2), we have the error guarantee [3, 30]:
|h`(A)− tr(A)| < ‖A‖F . (1) 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 (1) with o(1/ 2) matrix-vector products.

heading: 1.2 Dynamic Trace Estimation
text: We explore a natural and widely applicable dynamic version of the implicit trace estimation problem: given access to a matrix-vector multiplication oracle for a dynamically changing matrix A, maintain an approximation to 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 [43], edges are added to a network to optimally increase the “Estrada index” [11]. Evaluating this connectivity measure requires computing tr(exp(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 estimator to every snapshot of A as it changes over time. To achieve a guarantee like (1) for m time steps, we require O(m/ 2) 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 A1, ..., Am be n× n matrices satisfying: 1. ‖Ai‖F ≤ 1, for all i ∈ [1,m]. 2. ‖Ai+1 −Ai‖F ≤ α, for all i ∈ [1,m− 1].
Given implicit matrix-vector multiplication access to eachAi in sequence, the goal is to compute trace approximations t1, . . . , tm for tr(A1), ...., tr(Am) such that, for each i ∈ 1, . . . ,m,
P[|ti − tr(Ai)| ≥ ] ≤ δ. (2)
Above A1, . . . , Am represent different snapshots of a dynamic matrix at m time steps. We require ‖Ai‖F ≤ 1 only to simplify the form of our error bounds – no explicit rescaling is necessary for matrices with larger norm. If we assume ‖Ai‖F ≤ U for some (unknown) upper bound U , the guarantee of (2) would simply change to involve a U terms instead of . The second condition bounds how much the matrices change over time. Again for simplicity, we assume a fixed upper bound α on the difference at each time step, but the algorithms presented in this paper will be adaptive to changing gaps between Ai and Ai+1, and will perform better when these gaps are small on average. By triangle inequality, α ≤ 2, but in applications we typically have α 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 estimator repeatedly to each Ai.
Note on Matrix Functions. In many applications Ai = f(Bi) for a dynamically changing matrix B. While we may have ‖Ai+1 − Ai‖F = ‖f(Bi+1) − f(Bi)‖F ‖Bi+1 − Bi‖F for functions like the matrix exponential, this is not an immediate issue. To improve on Hutchinson’s estimator, the important requirement is simply that ‖Ai+1 − Ai‖F is small in comparison to ‖Ai+1‖F . 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 1 in the matrix-vector multiplication oracle model of computation, meaning that we consider the cost of matrix-vector products (which are the only way A1, . . . , Am can be accessed) to be significantly larger than other computational costs. We thus seek solely to minimize the number of such products used [39]. 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].

heading: 1.3 Main Result
text: Our main result is an algorithm for solving Problem 1 more efficiently than Hutchinson’s estimator: Theorem 1.1. For any , δ, α ∈ (0, 1), the DeltaShift algorithm (Algorithm 1) solves Problem 1 with
O ( m · α log(1/δ)
2 +
log(1/δ)
2 ) total matrix-vector multiplications involving A1, . . . , Am.
For large m, the first term dominates the complexity in Theorem 1.1. For comparison, a tight analysis of Hutchinson’s estimator [29] establishes that the naive approach requires O ( m · log(1/δ)/ 2 ) ,
which is worse than Theorem 1.1 by a factor of α. A natural setting is when α = O( ), in which case Algorithm 1 requires O(log(1/δ)/ ) matrix-multiplications on average over m time steps, in comparison to O(log(1/δ)/ 2) for Hutchinson’s estimator, a quadratic improvement in .
To prove Theorem 1.1, we introduce a dynamic variance reduction scheme. By linearity of trace, tr(Ai+1) = tr(Ai) + tr(∆i), where ∆i = Ai+1 −Ai. Instead of directly estimating tr(Ai+1), we combine previous estimate for tr(Ai) with an estimate for tr(∆i), computed via Hutchinson’s estimator. Each sample for Hutchinson’s estimator applied to ∆i requires just two matrix-vector multiplies: one with Ai and one with Ai+1. At the same time, when ∆i has small Frobenius norm (bounded by α), we can estimate its trace more accurately than tr(Ai+1)1. While intuitive, this approach requires care to make work. In a naive implementation, error in estimating tr(∆1), tr(∆2), . . . , compounds over time, eliminating any computational savings. To avoid this issue, we introduce a novel damping strategy that actually estimates tr(Ai+1 − (1− γ)Ai)) for a positive damping factor γ. 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 A1, . . . , Am and each Ai+1 −Ai than Problem 1.

heading: 1.4 Related Work
text: 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 [1, 12, 26, 30]. In these results, the trace of a matrix A is estimated by decomposing A = B+ ∆ where B has an easily computed trace (e.g., because it is low-rank) and ‖∆‖F ‖A‖F , so tr(∆) is more easily approximated with Hutchinson’s estimator than tr(A) directly.

heading: 2 Preliminaries
text: Notation. We let B ∈ Rm×k denote a real-valued matrix with m rows and k columns. x ∈ Rn denotes a real-valued vector with n entries. Subscripts like Bi or xj typically denote a matrix or vector in a sequence, but we use double subscripts with matrices to denote entries: Bij being the entry at the ith row and jth column. Let σ`(B) denote the `th singular value of B. ‖B‖F denotes the Frobenius norm of B, √∑ i,j B 2 ij = ∑ ` σ`(B) 2. ‖B‖∗ denotes the nuclear norm, ∑ ` σ`(B). We
let E[v] and Var[v] denote the expectation and variance of a random variable v. Hutchinson’s Estimator. Our algorithm uses Hutchinson’s trace estimator with Rademacher ±1 random variables as a subroutine. Specifically, let g1, . . . , g` ∈ Rn be independent random vectors, with each entry +1 or −1 with probability 1/2. Let A ∈ Rn×n. Hutchinson’s estimator for tr(A) is:
h`(A) = 1
` ∑̀ i=1 gTi (Agi) (3)
1While we consider the general, unstructured problem setting, we note that, if Ai has additional structure, it is not necessarily easier to estimate the trace of ∆i than that of Ai. For example, if Ai is a PSD matrix then specialized trace estimation algorithms that improve on Hutchinson’s method can be used [30]. Understanding dynamic traces estimation methods for sequences of structured matrices is a natural direction for future work.
Fact 2.1 (Hutchinson’s expectation and variance). For any positive integer ` and matrix A we have:
E[h`(A)] = tr(A), Var[h`(A)] = 2
`
( ‖A‖2F − n∑ i=1 A2ii ) ≤ 2 ` ‖A‖2F .
Fact 2.1 follows from simple calculations, found e.g. in [3]. Similar bounds hold when Hutchinson’s estimator is implemented with different random variables. For example, random Gaussians also lead to a variance bound of 2` ‖A‖2F . However, Rademachers tend to work better empirically. Given Fact 2.1, 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 E[X] = µ and variance Var[X] = σ2, for any k ≥ 1, P (|X − µ| ≥ kσ) ≤ 1/k2. Claim 2.3. For any , δ ∈ (0, 1), if ` = 2 2δ then Pr [|h`(A)− tr(A)| ≥ ‖A‖F ] ≤ δ.
The δ dependence in Claim 2.3 can be improved from 1δ to log(1/δ) via the Hanson-Wright inequality, which shows that h`(A) is a sub-exponential random variable [30, 35]. We also require HansonWright to obtain our bound involving log(1/δ). From this tighter result, Hutchinson’s yields a total matrix-vector multiplication bound of O(m · log(1/δ)/ 2) for solving Problem 1 by simply applying the estimator in sequence to A1, . . . , Am.

heading: 3 Main Algorithmic Result
text: As discussed in Section 1.2, a natural idea for solving Problem 1 with fewer than O(m/ 2) queries is to take advantage of the small differences between Ai+1 and Ai to compute a running estimate of the trace. In particular, instead of estimating tr(A1), tr(A2), . . . , tr(Am) individually using Hutchinson’s estimator, we denote ∆i = Ai −Ai−1 and use linearity of the trace to write:
tr(Aj) = tr(A1) + j∑ i=2 tr(∆i). (4)
By choosing a large `0, we can compute an accurate approximation h`0(A1) to tr(A1). Then, for j > 1 and ` `0, we can approximate tr(Aj) via the following unbiased estimator:
tr(Aj) ≈ h`0(A1) + j∑ i=2 h`(∆i) (5)
Since ‖∆i‖F ≤ α ‖Ai+1‖F , we expect to approximate tr(∆2), . . . , tr(∆m) much more accurately than tr(A2), . . . , tr(Am) directly. At the same time, the estimator in (5) only incurs a 2 factor overhead in matrix-vector multiplies in comparisons to Hutchinson’s: it requires 2 · (m − 1)` to compute h`(∆2), . . . , h`(∆m) versus (m − 1)` to compute h`(A2), . . . , h`(Am). The cost of the initial estimate h`0(A1) is necessarily higher, but can be amortized over time.

heading: 3.1 Our Approach
text: While intuitive, the problem with the approach above is that error compounds due to the sum in (5). Each h`(∆i) is roughly α/ √ ` away from tr(∆i), so after j steps we naively expect total error
O(j · α/ √ `). We can do slightly better by arguing that, due to their random nature, error actually accumulates as O( √ j · α/ √ `), but regardless, there is accumulation. One option is to “restart” the estimation process: after some number of steps q, throw out all previous trace approximations, compute an accurate estimate for tr(Aq), and for j > q construct an estimator based on tr(Aj) = tr(Aq) + ∑j−1 i=q+1 tr(∆i). While possible to analyze theoretically, this approach turns out to be difficult to implement in practice due to several competing parameters (see details in Section 5).
Instead, we introduce a more effective approach based on a damped variance reduction strategy, which is detailed in Algorithm 1, which we call DeltaShift. Instead of being based on (4), DeltaShift uses the following recursive identity involving a fixed parameter 0 ≤ γ < 1 (to be chosen later):
tr(Aj) = (1− γ)tr(Aj−1) + tr(∆̂j), where ∆̂j = Aj − (1− γ)Aj−1. (6)
Algorithm 1 DeltaShift Input: Implicit matrix-vector multiplication access to A1, ..., Am ∈ Rn×n, positive integers `0, `, damping factor γ ∈ [0, 1]. Output: t1, . . . , tm approximating tr(A1), . . . , tr(Am).
Initialize t1 ← 1`0 ∑`0 i=1 g T i A1gi, where g1, . . . , g`0 ∈ Rn are random ±1 vectors for j ← 2 to m do Draw ` random ±1 vectors g1, . . . , g` ∈ Rn z1 ← Aj−1g1, . . . , z` ← Aj−1g`, w1 ← Ajg1, . . . , w` ← Ajg` tj ← (1− γ)tj−1 + 1` ∑` i=1 g T i (wi − (1− γ)zi) end for
Given an estimate tj−1 for tr(Aj−1), DeltaShift estimates tr(Aj) by (1 − γ)tj−1 + h`(∆̂j). This approach has several useful properties: 1) if tj−1 is an unbiased estimate for tr(Aj−1), tj is an unbiased estimate for tr(Aj), 2) ‖∆̂j‖F is not much larger than ‖∆j‖F if γ is small, and 3) by shrinking tj−1 by a factor of (1 − γ) when computing tj , 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 , δ, α ∈ (0, 1), Algorithm 1 run with γ = α, `0 = O ( log(1/δ)/ 2 ) , and ` = O ( α log(1/δ)/ 2 ) solves Problem 1. In total, it requires
O ( m · α log(1/δ)
2 +
log(1/δ)
2 ) matrix-vector multiplications with A1, . . . , Am.
The full proof of Theorem 1.1 relies on the Hanson-Wright inequality, and is given in Appendix B. Here, we give a simple proof of essentially the same statement, but with a slightly weaker dependence on the failure probability δ.
Proof. Let γ = α, `0 = 2 2δ , and ` = 8α 2δ . The proof is based on an inductive analysis of the variance of tj , the algorithms estimate for tr(Aj). Specifically, we claim that that for j = 1, . . . ,m:
Var[tj ] ≤ δ 2. (7) For the base case, j = 1, (7) follows directly from Fact 2.1 because t1 is simply Hutchinson’s estimator applied to A1, and ‖A1‖F ≤ 1. For the inductive case, tj is the sum of two independent estimators, tj−1 and h`(∆̂j). 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, ‖∆̂j‖F = ‖Aj−(1−γ)Aj−1‖F ≤ ‖Aj − Aj−1‖F + γ‖Aj−1‖F ≤ 2α. Thus, by Fact 2.1, Var[h`(∆̂j)] ≤ 8`α2. Combined with the inductive assumption that Var[tj−1] ≤ δ 2, we have:
Var[tj ] = (1− γ)2Var[tj−1] + Var[h`(∆̂j)] ≤ (1− α)2δ 2 + 8α2
` ≤ (1− α)δ 2 + αδ 2 = δ 2.
This proves (7), and by Chebyshev’s inequality we thus have Pr [|tj − tr(Aj)| ≥ ] ≤ δ for all j.

heading: 3.2 Selecting γ in Practice
text: While DeltaShift is simple to implement, in practice, its performance is sensitive to the choice of γ. For the Theorem 1.1 analysis, we assume γ = α, 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 vj−1 = Var[tj−1] be the variance of our estimator for tr(Aj−1). We have that vj = (1− γ)2vj−1 + Var[h`(Aj − (1− γ)Aj−1)] ≤ (1 − γ)2vj−1 + 2` ‖Aj − (1 − γ)Aj−1‖2F . At time step j, a natural goal is to choose damping parameter γ∗ that minimizes this upper bound on the variance of tj :
γ∗ = arg min γ
[ (1− γ)2vj−1 + 2
` ‖∆̂j‖2F
] , (8)
where ∆̂j = Aj − (1 − γ)Aj−1 as before. While (8) cannot be computed directly, observing that ‖B‖2F = tr(BTB) for any matrix B, the above quantity can be estimated as ṽj = (1− γ)2ṽj−1 + 2 `h`(∆̂ T j ∆̂j), where ṽj−1 is an estimate for vj−1. The estimate h`(∆̂ T j ∆̂j) can be computed using exactly the same ` matrix-vector products with Aj and Aj−1 that are used to estimate tr(∆̂j), so there is little computational overhead. Moreover, since ∆̂Tj ∆̂j is positive semidefinite, as long as ` ≥ log(1/δ), we will obtain a relative error approximation to its trace with probability 1− δ [3]. An alternative approach to estimating vj would be to simply compute the empirical variance of the average h`(∆̂j), but this requires fixing γ. An advantage of our closed form approximation is that it can be used to analytically optimize γ. Specifically, expanding ∆̂Tj ∆̂j , we have that:
ṽj = (1− γ)2ṽj−1 + 2
`
( h`(A T j Aj)+ (1− γ)2h`(ATj Aj)− 2(1− γ)h`(ATj−1Aj) ) . (9)
Above, each estimate h` is understood to use the same set of random vectors. Taking the derivative and setting to zero, we have that the minimizer of (9), denoted γ̃∗, equals:
γ̃∗ = 1− 2h`(A
T j−1Aj)
`ṽj−1 + 2h`(ATj−1Aj−1) . (10)
This formula for γ̃∗ 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, `. For simplicity, unlike Algorithm 1, we do not use a larger number of matrix-vector multiplies when estimating A1. This leads to somewhat higher error for the first matrices in the sequence A1, . . . , Am, but error quickly falls for large j.

heading: 4 Algorithm Improvements and Lower Bound
text: In this section, we prove a lower bound showing that, in general, Theorem 1.1 is likely optimal. On the other hand, we show that, if we make a slightly stronger assumption on A1, . . . , Am and the dynamic updates ∆i = Ai+1 −Ai, an improvement on DeltaShift is possible.

heading: 4.1 Lower Bound
text: As noted, for a large number of time steps m, the matrix-vector multiplication complexity of DeltaShift is dominated by the leading term in Theorem 1.1, O(m · α log(1/δ)/ 2). We show that it is unlikely an improvement on this term can be obtained in general:
Lemma 4.1. Suppose there is an algorithm S that solves Prob. 1 with o(m · α log(1/δ)/ 2) total matrix-vector multiplies with A2, . . . , Am, and any number of matrix-vector multiplies with A1 when α = 1/(m− 1). Then there is an algorithm T that achieves (1) for a single A with o(log(1/δ)/ 2) matrix-vector multiplies.
Proof. The proof is via a direct reduction. Given a matrix A, positive integer m > 1, and parameter α = 1m−1 , construct the sequence of matrices:
A1 = 0, A2 = α ·A, . . . A1/α = (1− α)A, Am = A Since A = 0, and every A2, . . . , Am is a scaling of A, any algorithm S satisfying the assumption of Lemma 4.1 can be implemented with o ( (m− 1) · α log(1/δ)/ 2 ) = o(log(1/δ)/ 2) matrix-vector multiplications with A. Moreover, if S is run on this sequence of matrices, on the last step it outputs an approximation tm to tr(A) with Pr[|tm − tr(A)| ≥ ] ≤ δ. So algorithm T can simply simulate S on A1, . . . , Am and return its final estimate to satisfy (1).
Lemma 4.1 is a conditional lower-bound on matrix-vector query algorithms for solving Problem 1: if Hutchinson’s estimator cannot be improved for static trace estimation (and it hasn’t been for 30 years) then DeltaShift cannot be improved for dynamic trace estimation. We believe the bound could be made unconditional through a slight generalization of existing lower bounds on trace estimation in the matrix-vector multiplication model [30, 46].

heading: 4.2 Improved Algorithm
text: A recent improvement on Hutchinson’s estimator, called Hutch++, was described in [30]. For the static trace estimation problem, Hutch++ achieves a matrix-vector multiplication complexity of O(1/ ) to compute a relative error (1± ) approximation to the trace of any positive semi-definite matrix (PSD), improving on theO(1/ 2) required by Hutchinson’s. It does so via a variance reduction method (also used e.g. in [12]) which allocates some matrix-vector products to a randomized SVD algorithm which approximates the top singular vector subspace of A. This approximate subspace is projected off of A and Hutchinson’s used to estimate the trace of the remainder.
In our setting it is not realistic to assume PSD matrices – while in many applications A1, . . . , Am are all PSD, it is rarely the case the ∆1, . . . ,∆m−1 are. Nevertheless, we can take advantage of a more general bound proven in [30] for any matrix:
Fact 4.2 (Hutch++ expectation and variance). Let h++` (A) be the Hutch++ estimator of [30] applied to any matrix A with ` matrix-vector multiplications. We have:
E[h++` (A)] = tr(A), Var[h ++ ` (A)] ≤
16 `2 ‖A‖2∗.
Recall that ‖A‖∗ denotes the nuclear norm of A. Comparing to the variance of Hutchinson’s estimator from Fact 2.1, notice that the variance of Hutch++ depends on 1`2 instead of 1 ` , implying faster convergence as the number of matrix-vector products, `, increases. A trade-off is that the variance scales with ‖A‖2∗ instead of ‖A‖2F . ‖A‖2∗ is strictly larger, and possible a factor of n larger than ‖A‖2F . However, for matrices that are rank k, ‖A‖2∗ ≤ k‖A‖2F , so the norms are typically much closer for low-rank or nearly low-rank matrices. In many problems, ∆1, . . . ,∆m may have low-rank structure, in which case, an 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. ‖Ai‖∗ ≤ 1, for all i ∈ [1,m]. 2. ‖Ai+1 −Ai‖∗ ≤ α, for all i ∈ [1,m− 1]. Given matrix-vector multiplication access to each Ai in sequence, the goal is to compute trace approximations t1, . . . , tm for tr(A1), ...., tr(Am) such that, for all i, P[|ti− tr(Ai)| ≥ ] ≤ δ.
In Appendix C we prove the following result on a variant of DeltaShift that we call DeltaShift++:
Theorem 4.3. For any , δ, α ∈ (0, 1), DeltaShift++ (Algorithm 4) solves Problem 2 with
O ( m · √ α/δ + √ 1/δ ) total matrix-vector multiplications involving A1, . . . , Am.
Theorem 4.3 is stronger than Theorem 1.1 for vanilla DeltaShift in that it has a linear instead of a quadratic dependence on . In particular, its leading term scales as √ α/ 2, whereas Theorem 1.1 scaled with α/ 2. However, the result does require stronger assumptions on A1, . . . , Am and each ∆i = Ai+1 − Ai in that Problem 2 requires these matrices to have bounded nuclear norm instead of Frobenius norm. Since the nuclear norm of a matrix is strictly larger than its Frobenius norm, these requirements are stronger than those of Problem 1. As we will show in Section 5, the benefit of improved dependence often outweights the cost of these more stringent assumptions.

heading: 5 Experiments
text: We show that our proposed algorithm outperforms three alternatives on both synthetic and real-world trace estimation problems. Specifically, we evaluate the following methods:
Hutchinson. The naive method of estimating each tr(A1), . . . , tr(Am) using an independent Hutchinson’s estimator, as discussed in Section 2. NoRestart. The estimator of (5), which uses the same variance reduction strategy as DeltaShift for all j ≥ 2, but does not restart or add damping to reduce error accumulation.
Restart. The estimator discussed in Sec. 3.1, which periodically restarts the variance reduction strategy, using Hutchinson’s to obtain a fresh estimate for tr(Aj). Pseudocode is in Appendix A. DeltaShift. Our parameter free, damped variance reduction estimator detailed in Appendix A. We allocated a fixed number of matrix-vector queries, Q, to be used over all time steps 1, . . . ,m. For Hutchinson and DeltaShift, the same number of vectors Q/m was used at each step. For Restart and NoRestart, the distribution was non-uniform, and parameter selections are described in Appendix D.
Synthetic data: To simulate the dynamic setting, we generate a random matrix A ∈ Rn×n and add random perturbations for each of 100 time steps. We consider two cases: low (Fig. 1(a)) and significant (Fig. 1(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, tj , and the true trace tr(Aj). 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, so we 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 in a dynamic unweighted, undirected graphG via the fact that the number of triangles equals 16 tr(B
3), whereB is the adjacency matrix. The graph dataset we use is the Wikipedia vote network dataset with 7115 nodes [25, 24]. At each timestep we perturb the graph by adding a random k-clique, for k 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 matrix-vector products used by the estimators and the error reported as in the synthetic experiments(Appendix D). Note that for this particular application, the actual number of matrix-vector multiplications with B is 3Q, since each oracle call computes B(B(Bx)). As seen in Fig. 3, DeltaShift provides the best estimates overall.
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 a function of tr(exp(B)) for adjacency matrix B. This problem has also been explored in [7, 8, 4]. We
0
5000
10000
15000
20000
25000
30000
35000
T ot
al m
at ix
-v ec
to r
pr o
du ct
s (Q
)
Average error across 100 timesteps for 100 trials
Hutchinson NoRestart DeltaShift
use the road network data Gleich/minnesota (available at https://sparse.tamu.edu/Gleich/ minnesota). We perturb the graph over time by choosing two 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), as 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 when 100 matrix-vector products are used per time step.
The key takeaway from the experiments above is that DeltaShift and DeltaShift++ are able to obtain good dynamic trace approximations in far fewer 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 error across all time steps vs. total number of matrix-vector products(Q), for various values of Q. In order to achieve the accuracy level as DeltaShift, Hutchinson’s requires substantially more matrix-vector products.
Hessian spectral density: Finally, we evaluate the performance of DeltaShift on an application pertaining to a dynamically changing Hessian matrix, H , involved in training a neural network. As discussed in Section 1, a common goal is to approximate the spectral density of H . 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 T0, . . . , Tq , and estimate tr(T0(H)), . . . , tr(Tq(H)). Other polynomial basis sets can also be used (e.g., Legendre polynomials). Experimental details are discussed in section Appendix D, but we summarized the results here. We implement matrix vector products withH using the PyHessian library [47], and report
Table 1: Average relative error for trace of Chebyshev polynomials of Hessian.
Hutchinson NoRestart DeltaShift
T1(H) 2.5e-02 3.7e-02 1.7e-02 T2(H) 1.2e-06 1.7e-06 8.0e-07 T3(H) 4.0e-02 4.1e-02 3.1e-02 T4(H) 1.5e-06 1.7e-06 1.0e-06 T5(H) 2.1e-02 4.3e-02 1.9e-02
Table 2: Average relative error for trace of Chebyshev polynomials of Hessian.
Hutchinson NoRestart DeltaShift
T1(H) 1.9e-02 5.0e-02 1.5e-02 T2(H) 1.2e-06 2.9e-06 9.9e-07 T3(H) 7.7e-02 9.4e-02 6.1e-02 T4(H) 1.7e-06 2.8e-06 1.5e-06 T5(H) 2.1e-02 4.2e-02 1.8e-02
average error over 25 training epochs for the Hessian of a ResNet model with 269722 parameters trained it on the CIFAR-10 dataset. As it is impossible to compute the true trace of these matrices, we use Hutchinson’s estimator with a greater number of queries as placeholder for ground-truth, and compare the performance against the computed values. As can be seen in Tables 1 and 2, 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.

heading: Acknowledgements
text: 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.

Review summaries:
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. 

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.$  The 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.$  The 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.  The algorithms extend their algorithm to the case where you have a stronger nuclear norm bound assumption, and present numerical simulations demonstrating their advantages.

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.

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.