{'fraction_non_alphanumeric': 0.0371871330217136, 'fraction_numerical': 0.011559096456284971, 'mean_word_length': 4.343257874015748, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimating the Trace of the Matrix Inverse by Interpolating from the Diagonal of an Approximate Inverse

A number of applications require the computation of the trace of a matrix that is implicitly available through a function. A common example of a function is the inverse of a large, sparse matrix, which is the focus of this paper. When the evaluation of the function is expensive, the task is computationally challenging because the standard approach is based on a Monte Carlo method which converges slowly. We present a different approach that exploits the pattern correlation, if present, between the diagonal of the inverse of the matrix and the diagonal of some approximate inverse that can be computed inexpensively. We leverage various sampling and fitting techniques to fit the diagonal of the approximation to the diagonal of the inverse. Depending on the quality of the approximate inverse, our method may serve as a standalone kernel for providing a fast trace estimate with a small number of samples. Furthermore, the method can be used as a variance reduction method for Monte Carlo in some cases. This is decided dynamically by our algorithm. An extensive set of experiments with various technique combinations on several matrices from some real applications demonstrate the potential of our method.


Introduction
Computing the trace of a matrix A that is given explicitly is a straightforward operation. However, for numerous applications we need to compute the trace of a matrix that is given implicitly by its action on a vector x, i.e., Ax. Specifically, many applications are interested in computing the trace of a function of a matrix f (A). Examples include estimating parameters in image restoration using the generalized cross-validation approach [1], exploring the inverse covariance matrix in uncertainty quantification [2,3], computing observables in lattice quantum chromodynamics (LQCD) [4], or counting triangles in large graphs [5]. The matrix A is large, and often sparse, so its action f (A)x is typically computed through iterative methods. Because it is challenging to compute f (A) explicitly, the Monte Carlo (MC) approach has become the standard method for computing the trace [6,7]. The main purpose of this paper is to develop practical numerical techniques to address the computation of the trace of the inverse of a large, sparse matrix. But our technique can also be adapted to other functions such as the trace of the logarithm (yielding the determinant) or the trace of the matrix exponential.
For small size problems, computing A −1 through a dense or sparse LDU decomposition is the most efficient and accurate approach [8]. This works well for discretizations of differential operators in low dimensions but becomes intractable in high dimensional discretizations. For larger size problems, domain decomposition and divide and conquer strategies are more tractable but still expensive [9]. In many cases, however, a low accuracy approximation is sufficient. Numerous methods have been presented to address this need for estimating the trace of the inverse of symmetric positive definite matrices through Gaussian quadratures [7,10], modified moments [11,12], and MC techniques [7,10,13,14,11,6].
MC methods for computing the trace of a matrix are based on the structure of the Hutchinson method [1], which iteratively computes an average of matrix quadratures with random vectors. Variants of MC estimators are mainly analyzed and compared based on the variance of one sample [6,15], which depends on the quality of the selected random vectors. For real matrices, choosing random vectors having each element ±1 with equal probability is known to minimize variance over all other choices of random vectors [1,6] and therefore has been widely used in many applications. For complex matrices, the same result holds for vectors with ±1, ±i elements. In [6], Avron and Toledo analyze the quality of trace estimators through three different metrics such as trace variance, (ε, δ)-approximation of the trace, and the number of random bits for different choices of random vectors. In [15], Khorasani and Ascher improve the bounds of (ε, δ)-approximation for the Hutchinson, Gaussian and unit vector estimators. However, the structure of the matrix can also be exploited to accelerate the convergence of the Hutchinson method.
There has been a number of efforts to combine the Hutchinson method with well-designed vectors based on the structure of the matrix [13,14,16,4]. In [14], they use columns of the Hadamard matrix, rather than random vectors, to systematically capture certain diagonals of the matrix. Then, the MC iteration achieves the required accuracy by continuously annihilating more diagonals with more Hadamard vectors. However, the location of the nonzeros, or of the large elements of A −1 , often does not coincide with the diagonals annihilated by the Hadamard vectors. In [16], graph coloring and probing vectors are used to identify and exploit special structures, such as bandedness or decaying properties in the elements of A −1 , to annihilate the error contribution from the largest elements. However, if the error for the chosen number of colors is large, all work has to be discarded and the probing procedure repeated until the accuracy is satisfied. In [4], the authors introduce hierarchical probing on lattices to avoid the previous problems and achieve the required accuracy in an incremental way. For all these approaches, the approximation error comes from non-zero, off-diagonal elements that have not been annihilated yet. Instead, this paper looks only at the main diagonal of A −1 .
Our motivation for focusing only on the main diagonal is that the trace of A −1 is simply a summation of a discrete, 1-D signal of either the eigenvalues or the diagonal elements of A −1 . Although we cannot compute all the diagonal elements, we may have an approximation to the whole signal from the diagonal of an approximation of A −1 (e.g., of a preconditioner). If the two diagonals have sufficiently correlated patterns, fitting methods can be used to refine the approximation both for the diagonal and the trace. Therefore, the proposed method may serve as a standalone kernel for providing a good trace estimate with a small number of samples. But it can also be viewed as a preprocessing method for stochastic variance reduction for MC in cases where the variance reduces sufficiently. This can be monitored dynamically by our method.
We present several techniques that improve the robustness of our method and implement dynamic error monitoring capabilities. Our extensive experiments show that we typically obtain trace estimates with much better accuracy than other competing methods, and in some cases the variance is sufficiently reduced to allow for further improvements through an MC.

Preliminaries
We denote by . the 2-norm of a vector or a matrix, by N the order of A, by Z an approximation of A, by D the diagonal elements of A −1 , by M the diagonal elements of Z −1 , by Tr( f (A)) the trace of the matrix f (A), and by extension, Tr(D) the sum of the elements of the vector D, by T e i ( f (A)) the MC trace estimator of f (A) using unit (orthocanonical) vectors, by T Z 2 ( f (A)) the MC trace estimator of f (A) using Rademacher vectors, by diag(.) the diagonal operator of a matrix, and by Var(.) the variance operator of a random variable or a vector.

Hutchinson trace estimator and unit vector estimator
The standard MC method to estimate the trace of the matrix inverse is due to Hutchinson [1]. It estimates the Tr(A −1 ) by averaging s quadratures with random vectors z j ∈ Z N 2 = {z(i) = ±1 with probability 0.5}, The variance of this method is given by The variance of this trace estimator is proven to be minimum over all vectors with real entries [1]. The confidence interval of a MC method reduces as O( Var(T Z 2 (A −1 ))) for the given matrix. The unit vector estimator uniformly samples s vectors from the orthocanonical basis {e 1 , . . . , e N } [6], where i j are the random indices. The variance of the unit vector estimator is given by The variance of the Hutchinson method depends on the magnitude of the off-diagonal elements. It converges in one step for diagonal matrices and rapidly if A −1 is highly diagonal dominant. On the other hand, the variance of the unit vector estimator depends only on the variance of the diagonal elements. It converges in one step if the diagonal elements are all the same and rapidly if the diagonal elements are similar. Thus, the method of choice depends on the particular matrix.

Reducing stochastic variance through matrix approximations
Given an approximation Z ≈ A, for which Z −1 and Tr(Z −1 ) are easily computable, we can decompose where E = A −1 − Z −1 . We hope that by applying the MC methods on E instead on A −1 , the variance of the underlying trace estimator, in (2) or (4), can be reduced, thereby accelerating the convergence of MC. Among many ways to obtain a Z, we focus on the following two. The first approach is when Z −1 = (LU) −1 , where the L, U matrices stem from an incomplete LU (ILU) factorization of A, one of the most commonly used preconditioners. If the ILU is sufficiently accurate, then M = diag(Z −1 ) may be a good approximation to D. To obtain the vector M without computing the entire Z −1 , we can use an algorithm described in [17]. This algorithm requires the computation of only those entries Z −1 i j for which L i j or U i j = 0. If the L,U factors are sufficiently sparse or structured, this computation can be performed efficiently (see [18] for an example in the symmetric case).
The second approach is a low rank approximation Z −1 = V Λ −1 U T , where Λ is a diagonal matrix with the k N smallest singular values of A, and U and V are the corresponding left and right singular vectors. Eigenvalues and eigenvectors can be used instead but we do not consider it in this paper. This subspace can be obtained directly by an iterative eigensolver [19,20] as a preprocessing step. The cost of this procedure is relatively small since the singular space is not needed in high accuracy. Alternatively, the space can be approximated by methods such as eigCG [21] or eigBiCG [22] while solving linear systems of equations during the MC method. This incremental approach adds only minimal overhead to the iterative linear solver but it cannot compute as many and as good quality singular vectors as the first approach. The quality of the approximation of A −1 by Z −1 depends on the separation of the computed singular space. Therefore, depending on the matrix, more singular triplets may be needed for a good low rank approximation. On the other hand, this space can also be used to deflate and thus accelerate subsequent linear systems. Once the singular space is computed, each diagonal element can be obtained with a single inner product of short vectors, M i = V (i, :) T U(i, :)/λ i , and thus it is computationally inexpensive. Finally, the incremental SVD approach requires some special algorithmic attention during our algorithm, which will be pointed out later.
A computationally inexpensive, albeit less accurate approach for computing an approximation M is based on variational bounds on the entries of A −1 [7,23]. Upper and lower bounds on the i-th diagonal entry ii are derived inexpensively since they only depend on estimates of the smallest and largest algebraic eigenvalues, λ 1 , λ N , and the entries of A. The bounds apply to both symmetric and unsymmetric matrices. For the case of a real symmetric A, we have [23], where s i j = ∑ n k=1 A ik A k j . However the bounds in (6) will not be sharp especially the upper bound [11] and the error in the approximation can be large.
In general, unless Z −1 is highly accurate, we do not expect Tr(M) to be close to Tr(A −1 ). However, the patterns of M and D often show some correlation. We demonstrate this for two example matrices, delsq50 and orsreg2205, in Figures 1 and 2 using ILU and SVD respectively. For matrix orsreg2205, both the ILU and SVD approaches return an approximate diagonal M which captures the pattern of D very well, with the M returned by ILU being slightly better than the one from SVD. For matrix delsq50, SVD clearly captures the pattern of D better than what ILU does. As in preconditioning for linear systems of equations, the appropriate approximation technique depends on the given matrix.

Comparison of different MC methods and discussion on importance sampling
Based on (2-4), we express the variance of the trace estimators T Z 2 (E) and T e i (E) as follows:  Var(T e i (E)) = N 2 s Var(diag(E)).
Figures 1 and 2 show there is potential for the variances of T Z 2 (E) or T e i (E) to be smaller than those of T Z 2 (A −1 ) and T e i (A −1 ). However, for a given matrix, we must gauge which MC method would be better, and whether the variances need further improvement. The estimator T e i (E) has the interesting property that if M = D+c, where c is a constant, then its variance in (8) is zero and we obtain the correct trace in one step. Although we cannot expect this in practice, it means that the shift observed between M and D in Figure 2 To obtain an estimator of the trace of D with importance sampling, we replace the uniform sampling of D i values with sampling with probability G i [24]. Then, instead of (3), the importance sampling estimator is: When M = cD, the variance of T IS (D) is zero and it finds the trace in one step. However, it completely fails to identify shift correlations. In general D and M have a more complex relationship that neither T IS (D) or T e i (E) can capture. This motivates our idea to explore general fitting models to approximate D.

Approximating the trace of a matrix inverse
We seek to construct a function f , such that D ≈ f (M). Then we can decompose Tr( f (M)) is trivially computed for a given f . A key difference between the approaches in (10) and (5) is that it is easier to find a fitting of the two vectors M and D if a strong pattern correlation exists between them than to fit the entire matrices A −1 and Z −1 . If Tr( f (M)) is a good approximation to Tr(A −1 ) and its accuracy can be evaluated easily, we can directly use this quantity; Otherwise, we can apply the unit vector MC estimator to compute Tr(E f it ) = Tr(D − f (M)), provided that its variance is smaller than the variances in (2), (3), (7) and (8). (3) 6: Estimate the relative trace error and, if needed, the variances for different MC methods (Section 4)

Algorithm 1 Basic algorithm for approximating
The basic description of the proposed estimator is outlined in Algorithm 1. First, our method computes an approximation M of D using one of the methods discussed in the previous section. If that method is based on a preconditioner or a low rank approximation, that preconditioner Z −1 could also be used to speed up the solution of linear systems in step 3. Second, it finds a fitting sample S f it , a set of indices that should capture the important distribution characteristics of D. Since we have no information about D, Section 3.1 discusses how to tackle this task by considering the distribution of M. Third, it computes the values of D(S f it ) by solving the corresponding linear systems using a preconditioned iterative solver. Since this is the computational bottleneck, the goal is to obtain good accuracy with far fewer fitting points than the number of vectors needed in MC. Fourth, it computes a fitting model that has sufficient predictive power to improve the diagonal approximation. This critical task is discussed in Section 3.2. Finally, since there are no a-posteriori bounds on the relative error of the trace, we use a combination of statistical approaches and heuristics to estimate it incrementally at every step. If the estimated error is not sufficiently small, our method can be followed by an MC method. For this reason, the algorithm also estimates dynamically the variances of the two different MC estimators (7) and (8) so that the one with the smallest variance is picked. The dynamic evaluation of error and variances is discussed in Section 4.

Point Identification Algorithm
We need to identify a set of indices S f it based on which we can compute a function f that fits f (M(S f it )) to D(S f it ) so that |Tr( f (M)) − Tr(D)| is minimized. It is helpful to view Tr(D) as an integral of some hypothetical one dimensional function which takes the values D(i) = D i on the discrete points i = 1, . . . , n. Our goal is to approximate n 0 D(x)dx with far fewer points than n. Because it is one dimensional, it may be surprising that Monte Carlo is the standard approach and not numerical integration. However, the effectiveness of any numerical quadrature rule relies on the smoothness of the function to be integrated, specifically on the magnitude of its higher order derivatives. In our case, our hypothetical function may have no smoothness or bounded derivatives. More practically, we deal with a set of discrete data points in D for which smoothness must be defined carefully. A typical definition of smoothness for discrete data is based on the Lipschitz continuity, i.e., |D i − D j | < c|i − j| [25]. The lower the constant c, the smoother the set of data. Hence, for a continuous function, the larger the magnitude of its first derivative, the less smooth its discretized points are. In the context of integrating an arbitrary D, the first order divided difference (i.e., the discretized first derivative) D i − D i+1 may be arbitrarily large in magnitude and therefore numerical integration may not work any better than random averaging. Even for matrices that model physically smooth phenomena (e.g., in PDEs), the matrix may be given in an ordering that does not preserve physical locality.
Consider a sorted permutation of the diagonal,D = sort(D). Obviously Tr(D) = Tr(D), butD is monotonic and, based on the definition of Lipschitz discrete smoothness, it is maximally smooth, i.e., for the smallest possible ∆ ∈ R + among all permutations of D. Monotonicity implies that, in the absence of any additional information about the data, a simple trapezoidal rule minimizes the worst case integration error [26]. In addition, if we are allowed to choose the integration points sequentially, based on the points computed so far, then a much better average case error can be obtained [27,28]. On the other hand, if bounds are known on the discrete smoothness ofD, better worst case error bounds can be established. Since D, however, is not available, we turn to its approximation M.
A close pattern correlation between M and D means that the elements of M should have a similar distribution as those of D, or thatM = sort(M) should be similar toD. Let us assume for the moment, that the index that sorts D toD sorts also M toM. In other words, we assume a complete correlation between the monotonic orderings of M and D even though their values may differ. Then, we can work on the surrogate modelM for which we can afford to identify the best quadrature points that yield the smallest error in Tr(M). These will be the ideal points for computing the fitting function f . Specifically, we need to select indices that capture the important distribution changes inM. For example, identifying minimum and maximum elements of M sets the range of approximation for f and avoids extrapolation. We also look for entries inM that deviate highly from their neighbors (where the first order divided difference ofM and hopefully ofD has a large value). This strategy has three advantages. First, between such indices the data is smooth so the integral Tr(M) should be captured well by the trapezoidal rule. Second, and more important, we can obtain a more accurate fitting function f in a piecewise manner in intervals whereD has similar behavior. Third, in this way we avoid picking points with the same valuê M i =M j which can create problems in the fitting function.
The proposed index selection method is shown in Algorithm 2. Initially the set of sampled indices, S f it , includes the indices of the extrema ofM, 1 and N. Then, for every interval (i, j), with i, j consecutive indices inŜ f it , we find the index t from i + 1 to j − 1 that minimizes the trapezoidal rule error for computing Tr(M(i : j)): The process continues until it reaches a maximum number of sampling points or until the maximum error over all intervals decreases by a fixed value, say 0.001. The value of this threshold depends on how close the patterns of M and D are correlated, a question we had postponed and we discuss next. The points selected by the previous algorithm minimize the error for computing Tr(M) with only |Ŝ f it | points. How good are these points for fittingM toD and thus obtaining a small error in Tr(D)? Consider the M and D as discrete functions from {1, . . . , n} → R and assume they are bijective in their ranges. Let G, J : {1, . . . , n} → {1, . . . , n} be the two permutation functions that sort M toM and D toD, i.e.,M(i) = M(J(i)) andD(i) = D(G(i)). Our method relies on the assumption that the distribution ofM captures the distribution ofD, or that there exists a smooth function f that can be approximated well with a low degree polynomial that fitsD Our algorithm first picks an index i forM and derives the original index k = J(i) in M. Since the permutation g is unknown, we do not know the D(m) that should be matched with M(k). Thus, the algorithm makes the simplifying assumption that g = s and computes D(k) instead. However, this corresponds to a different D( j), where G( j) = k. Therefore, using (12), we can bound the error for this mismatch by This implies that as long as the permutations J and G are locally similar, or in other words they do not shuffle the same index of M and D too far from each other, the error is small and therefore the fitting should work well. On the other hand, if M is a random permutation of D, there is no good fitting, even thoughM =D. This implies that the traces of f (M) and D may be very close but not their individual elements. It also means that its variance of E f it may not always be small (see the examples in the following section as well as in Section 4.2).
In practice, since we do not know the permutation G, it is possible that some flat area ofM is associated with important changes inD. To alleviate the effect of a possible local pattern mismatch, we instead empirically insert the midpoint of the current largest interval every k samples (we choose k = 5 in lines [11][12][13][14]. Returning to the choice of threshold in the algorithm, we see that going below a threshold helps the accuracy of Tr(M) but not necessarily of Tr(D). The reason is that the greedy point selection strategy becomes less effective in determining the local mismatching patterns. Therefore, we terminate searching for a new index based on (13) if the maximum error is less than 0.001. We found this threshold to be sufficient in our experiments. If the maximum required samples have not been generated, we continue by simply bisecting the largest intervals until maxPts is reached (in lines [16][17][18]. On exit, we compute S f it which mapsŜ f it to the original unsorted ordering. Add t intoŜ f it , numSamples = numSamples + 1 10: Push intervals (L,t) and (t, R) each with their corresponding tempErr in Q 11: if numSamples is a multiple of 5 then A typical sampling result produced by our method is shown in Figure 3. The graph on the left shows the diagonals M and D plotted in their original order. The pattern correlation between them is clear, but it is less straightforward how to pick the fitting points. Figure 3(b) shows both diagonals plotted with the order of indices that sortsM, i.e., M(J), D(J). The points picked by our algorithm based onM correspond to an almost monotonic sequence of points in D(J). The associated indices can then be used to identify a suitable set of entries in D to perform the numerical integration.

Two Fitting Models
Next we work on the monotonic sequenceM which we assume correlates well withD. Therefore, we consider a linear model and a piecewise polynomial model. The MC methods in (5) and (10)  Input : A ∈ R n×n Output : Tr(A −1 ) estimation: T f 1: Compute M using ILU or Eigendecompostion or SVD on A 2: Call Algorithm 2 to compute sample set S f it . 3 The linear model preserves the shape of M, and therefore relies exclusively on the quality of M. To take advantage of our premise that the distributionM approximates well the distribution ofD, our next fitting model is the Piecewise Cubic Hermite Spline Interpolation (PCHIP). It was proposed in [29] to construct a visually pleasing monotone piecewise cubic interpolant to monotone data. The PCHIP interpolant is only affected locally by changes in the data and, most importantly, it preserves the shape of the data and respects monotonicity. Therefore, we work onM and the indicesŜ f it = [1 = s 1 , s 2 , . . . , s k−1 , s k = N] which are given in an order such that α =M(s 1 ) ≤M(s 2 ) ≤ · · · ≤M(s k ) = β is a partition of the interval Notice that p(x) will be monotone in the subintervals where the fitting points D(J −1 (s i )) are also monotone. Therefore, as long as M is close to D in the sense of (14), integration of p(x) will be very accurate. The PCHIP model is given in Algorithm 4. The first two steps are the same as in Algorithm 3. In step 3, we apply the function unique to remove the duplicate elements ofM(Ŝ f it )) to produce a sequence of unique values as required by PCHIP. This yields a subset of the indices,Ŝ f it , which is mapped to original indices as I = J −1 (Ŝ f it ) to be used in PCHIP.   Table 1 compares the relative error for the trace, as well as the variances of D, D − M and D − p(M) when using the linear LS and PCHIP models on two test matrices, OLM5000 and KUU. In both cases, the two models can provide a trace estimate of surprising relative accuracy O(1e − 2) with only 20 fitting points. In addition, for matrix OLM5000, the standard deviation of MC on D − p(M) is reduced by a factor of 10 compared to that of MC on D (a speedup of 100 in terms of samples). However, for matrix KUU, the standard deviation of MC on D − p(M) does not improve much over MC on D. This reflects the discussion in Section 3.1, and suggests that our method can serve as a standalone kernel for giving a fast trace estimate. In addition, in some cases the method can reduce the variance significantly to accelerate a second stage MC.

Dynamic Evaluation of Variance and Relative Trace Error
Since there are no a-posteriori bounds for the accuracy of our results, we develop methods that use the information from the solution of the linear systems to incrementally estimate the trace error and the variances of the resulting approximations. This approach is also useful when M is updated with more left and right eigenvectors or singular vectors obtained from the solution of additional linear systems.

Dynamic Variance Evaluation
To decide which MC method we should use after the fitting stage or even whether it is beneficial to use the fitting process for variance reduction, we monitor incrementally the variances Var(T e i (A −1 )), Var(T e i (E f it )), Var(T Z 2 (A −1 )), and Var(T Z 2 (E)), with the aid of the cross-validation technique [30].
Our training set is the fitting sample set D(S f it ), while our test set D(S mc ) is a small random set which is independent of the fitting sample set. If we want to combine our method with MC, eventually more samples need to be computed, and thus we can pre-compute a certain number of them as the test set D(S mc ). We have used the holdout method [31], a single train-and-test experiment for some data splitting strategy since the fitting sample set is fixed.
To compute D(S f it ) or D(S mc ) , a linear system with multiple right hand sides is solved as follows: The computed column vectors x i can be used to estimate the Frobenius norm of both A −1 and E = A −1 − Z −1 [32,33]. Then Var(T Z 2 (A −1 )) and Var(T Z 2 (E)) can be estimated as follows: Var(T Z 2 (E)) ≈ 2N where E(:, i) = Ee i . Simultaneously, based on the sampled diagonal elements A −1 ii , we can also update the evaluation of Var(T e i (E f it )), Var(T e i (E)) and Var(T e i (A −1 )). Here we only show the computation of the unbiased variance estimation for Var(T e i (E f it )) by: Var(E f it (S mc )).
Note that S f it should not be used for estimating the variance of the unit vector MC estimator since these sample points are exact roots of the PCHIP function.
We implement the dynamic variance evaluation scheme in Algorithm 5. In lines 2-4 and 7-9, the approximation M can be computed using an ILU factorization at the beginning of the procedure or be updated with increasing number of singular triplets or eigenpairs. Note that if M is obtained by a partial eigendecomposition or SVD, the updated M is different in two consecutive steps, and thus Algorithm 2 will return a slightly different index set S f it which may not be incremental. This may not provide a consistent improvement of the relative trace error during the fitting progress. In line 10 of the algorithm, we force the points to be incremental between steps i and i + 1 as follows; we generate the entire set S Call Algorithm 2 to find more indices so that S f it has i fitting points 11: Call Algorithms 3 or 4 to update approximation of Tr(A −1 )

12:
Estimate variances of different MC methods based on (17), (18) and (19) 13: end for that lie the closest to the previous index set S f it . This simple scheme works quite well experimentally. Figure 5 shows how the actual and the estimated variances for three MC methods match for the test matrix RDB5000. In addition, the relative difference between different MC methods becomes clear after only a few points which facilitates not only the proper choice of MC method but also an early decision to stop if further fitting is not beneficial. In the numerical experiments section we show that these results are typical for matrices from a wide variety of applications. ) and E (denoted as Var(E)) of the matrix RDB5000 with ILU and SVD respectively.

Monitoring Relative Trace Error
As discussed in Section 3.1 and showed in Table 1, the variance of MC on the vector E f it may not be reduced if the sorting permutations ofM andD are dissimilar. Even in such cases, however, our fitting method can provide a good trace estimation with only a small number of samples. We further investigate this by comparing the elements of D and p(M) with different orderings. Figure 6 in Var(T e i (E f it )). However, Figure 6(b) reveals that the distributions of the sorted p(M) and the sorted D almost coincide; hence, the two integrals Tr(p(M)) and Tr(D) are very close.
The main obstacle for using our method as a standalone kernel for trace estimation is that there is no known way to measure or bound the relative trace error. Resorting to the confidence interval computed by the variance of a MC estimator is pessimistic as our results show (see also [6,15]). We note the similarity to the much smaller error obtained in the average case of integrating monotonic functions with adaptive quadratures versus the worst case known bounds [26,27].
Motivated by these previous research results, we show how to develop practical criteria to monitor the relative trace error in our Algorithms. Suppose at each step of the fitting process in Algorithm 5 we collect a sequence of trace estimates T i , i ∈ [1,maxPts]. Consider the trace estimations in two successive steps, (20) As long as T i converges with more fitting points, the relative difference of two successive trace estimations can serve as an approximation to the relative error. However, when the global pattern provided by M and p(M) is not fully matched to that of D, convergence of T i stagnates until enough points have been added to resolve the various local patterns. To determine whether the current relative trace error estimation can be trusted, we present our second heuristic by considering the error bound of our fitting models.
When approximating f (M) with p(M), the PCHIP Hermite cubic splines with k points on the interval [α, β], the bound on the error E(M) = f (M) − p(M) is given by [34], where h = β−α k , and f (4) ∞,[α,β] denotes the maximum value of the fourth derivative of f in the entire interval [α, β]. Since f (4) ∞,[α,β] /384 is a constant, in two successive fitting steps we have, where k i and k i+1 are the number of fitting points in two consecutive steps. We can use (22) to estimate the maximum possible improvement between two consecutive trace errors. If the i + 1 trace error estimate reduces over the i-th estimate by a factor of more than (k i /k i+1 ) 4 , we do not trust it. One caveat is that (21)   the error in other intervals. Therefore, convergence might not be fully dictated by (22). In practice, we found that the improvement ratio is between O(k i /k i+1 ) and O((k i /k i+1 ) 4 ). Therefore, if the current relative trace error estimate is determined not to be trusted, we may instead use |E i+1 (M) The choice of k depends on the quality of M. In our experiments, we use the geometric mean of the four rates, yielding k = 9/4. Recall that the corresponding ratio in MC is O( k i /k i+1 ), which is much slower than the proposed trace estimation method.
Algorithm 6 combines the two heuristics in (20) and (22) to dynamically monitor the relative trace error. It is called after step 12 of Algorithm 5. Figure 7 shows two examples of how our dynamic method provides reasonable estimates of the relative trace error.

Numerical Experiments
We run experiments on matrices that are sufficiently large to avoid problems of sampling from small spaces but can still be inverted to obtain the exact trace error. We select matrices RDB5000 and cfd1 from the University of Florida sparse matrix collection [35] and generate three test matrices from applications that appear in [7]. The Heatflow160 matrix is from the discretization of the linear heat flow problem using the simplest implicit finite difference method. The matrix Poisson150 is from 5-point central difference discretization of the 2D Poisson's equation on a square mesh. The VFH6 matrix is from the transverse vibration of a Vicsek fractal that is constructed self-similarly. We also use matrix matb5 which is a discretization of the Wilson Dirac operator on a 8 4 uniform lattice with 12 degrees of freedom at each node, using a mass near to critical. Table 2 lists these matrices along with some of their basic properties. All experiments are conducted using MATLAB 2013a. The number of fitting points increases as s = 5 : 100. The approximation M is computed by ILU with parameters type = ilutp and droptol = 1E-2, or as a low rank approximation of 2s smallest singular vectors with accuracy 1E-6 (twice the number of fitting points at each step), or by the bounds on the diagonal.

Effectiveness of the fitting models
In Figure 8, we divide the diagonal elements of the matrix RDB5000 into three contiguous sets and zoom in the details. We see that despite a good M, the linear LS model cannot scale the entire M onto D. The more flexible piecewise approach of PCHIP results in a much better fit.
In Figure 9, we look at three matrices with M generated using the SVD. The PCHIP model typically has smaller relative trace error than the LS model. We also see that as more fitting points are sampled, the relative trace error of both models decreases significantly at early stages and slowly after a certain point. This relates to the quality of M, not of the model. Typically M will approximate the global pattern of D and the two can be matched well with only a few fitting points. But if the local patterns of M and D differ, a large number of fitting points will be required. This can be seen in Figure 10. Using the permutation ofM for each case, we plotM, D, and p(M) using 100 fitting points. For the Heatflow160 matrix, p(M) approximates D well everywhere except for the small leftmost part of the plot, which allows the relative error to reach below 10 −4 before convergence slows down (Figure 9). The behavior is similar for the Poisson150. The issue is more pronounced on matrix VFH6, where M and p(M) capture the average location of D but completely miss the local pattern, which is reflected by a very slowly improving error in Figure 9. We mention that the irregularity of the relative trace errors in Figure 9 relates to the variability of successive updates of M and of the sampling indices, especially when M is of lower quality. Figure 11 demonstrates that the PCHIP model has smaller actual variance for MC on E f it than the LS model. Therefore, we only consider the PCHIP model in the rest of experiments.

Comparison between the fitting model and different MC methods
We address the question of whether the number of matrix inversions we spend on computing the fitting could have been used more efficiently in an MC method, specifically the Hutchinson method on A −1 and the Hutchinson method on E. In Table 3 we compare the relative trace error of the PCHIP model with 20 fitting points against the relative errors of the two MC methods as computed explicitly from their respective standard deviations in (2) and (7), with s = 20, divided by the actual trace of D.
When M approximates D sufficiently well, the trace from the fitted diagonal is better, and for the ILU approximations far better, than if we just use the Hutchinson method on A −1 (the first column of results). Although MC on E exploits the ILU or SVD approximation of the entire matrix (not just the diagonal that our method uses), we see that it does not always improve on MC on A −1 , and in some cases (cfd1 with ILU) it is far worse. In contrast, our diagonal fitting typically improves on MC on E. The last column shows that even with an inexpensive diagonal approximation we obtain a similar or better error than MC on A −1 . The only exception is the matrix VFH6 where, as we saw earlier, M cannot capture the pattern of D. Even then, its error is close to the errors from the MC methods and, as we show next, the best method can be identified dynamically with only a small number of samples.
If the user requires better trace accuracy than our fitting technique provides, we explore the performance of the diagonal fitting as a variance reduction for MC with unit vectors. We compute the actual values of  (2) and (7) divided by the actual trace). Var(T e i (E f it )), Var(T e i (A −1 )), Var(T e i (E)), Var(T Z 2 (A −1 )) and Var(T Z 2 (E)) for every step s = 5 : 100 and show results for three matrices in Figure 11. Note that the low rank approximation uses 2s singular vectors.

Dynamic evaluation of variance and relative trace error
The above results emphasize the importance of being able to assess quickly and accurately the relative differences between the variances of different methods as well as the trace error, so that we can decide whether to continue with fitting or which MC method to switch to. First we show the effectiveness of the dynamic variance evaluation algorithm for our fitting MC method on (D − p(M)) with unit vectors, and on A −1 and E with Rademacher vectors. Then, we evaluate our algorithm for estimating the relative trace error during the fitting process. Figure 12 compares the estimated variances with the actual variances of the three MC methods when increasing the number of fitting points from 5 to 100. The approximation M is computed by using ILU. We can see that the estimated values of Var(T Z 2 (A −1 )) and Var(T Z 2 (E)) converge to the actual variances after only a few sample points. The estimated value of Var(T e i (E f it )) gets close to and captures the trend of the the actual variance as the fitting samples increase. Nevertheless, the relative differences between the variances of the various MC methods are apparent almost immediately. Figure 13 shows the same experiments when the approximation M is computed by SVD. Since M is updated each step, Var(T Z 2 (E)) and Var(T e i (E f it )) change accordingly. As with ILU, Var(T Z 2 (A −1 )) and Var(T Z 2 (E) can be estimated very well in a few steps. Var(T e i (E f it )) could be underestimated but the relative variance difference between these MC methods becomes clear when the fitting points increase beyond 20. Thus we are able to determine whether the fitting process is beneficial as a variance reduction preprocessing and which is the best MC method to proceed with for the trace estimation. Figure 14 compares the estimated relative trace error with the actual one in the cases of Figure 13. We observe that the estimation is accurate as the fitting samples increase, even for cases such as VFH6 where the fitting process is not as successful. Moreover, because our algorithm is based on upper bounds on the error of a piecewise cubic polynomial, the actual relative trace error could be lower than predicted.

A large QCD problem
The trace estimator presented in this paper has the potential of improving a number of LQCD calculations, where the trace of the Dirac matrix is related to an important property of QCD called spontaneous chiral symmetry breaking [21]. In our previous work [4], we presented the method of hierarchical probing that achieves almost optimal variance reduction incrementally and inexpensively.
As shown in Figure 15(a), a low rank approximation with 200 singular vectors yields a good approximation M and an excellent fit p(M). In Figure 15(b), we see that the actual relative trace error decreases very fast to O(10 −4 ) with increasing number of fitting points and singular vectors, and can be monitored well by our dynamic trace estimation algorithm. These singular vectors can be approximated while solving the linear systems with eigCG. Interestingly, we can improve the relative trace error of hierarchical probing by two orders of magnitude. In addition, the variances of different MC methods can be estimated dynamically to allow us to continue with the best estimator if needed (Figure 15(c)).   Figure  15(b), the green square denotes the relative trace error by applying hierarchical probing technique on deflated matrix E in [4].

Conclusion and future work
A novel method has been presented to estimate the trace of the matrix inverse by exploiting the pattern correlation between the diagonal of the inverse of the matrix and some approximation. The key idea is to construct a good approximation M ≈ D through eigenvectors or some preconditioner, sample important patterns of D by using the distribution of the elements of M, and use fitting techniques to obtain a better approximation p(M) ≈ D from where we obtain a trace estimate. The proposed method can provide a fast trace estimate with 2-3 digits relative accuracy given only a few samples while may or may not improve the variance of MC. When the variance is reduced sufficiently, our method can be also used as a diagonal estimator. We also propose an effective dynamic variance evaluation algorithm to determine the MC method with the smallest variance and a dynamic relative trace error estimation algorithm without any additional costs. We demonstrated the effectiveness of these methods through a set of experiments in some real applications.