A Multilevel Approach to Variance Reduction in the Stochastic Estimation of the Trace of a Matrix
Pages A2536 - A2556
Abstract
The trace of a matrix function $f(A)$, most notably of the matrix inverse, can be estimated stochastically using samples $x^*f(A)x$ if the components of the random vectors $x$ obey an appropriate probability distribution. However such a Monte Carlo sampling suffers from the fact that the accuracy depends quadratically on the samples used, thus making higher precision estimation very costly. In this paper, we suggest and investigate a multilevel Monte Carlo approach which uses a multigrid hierarchy to stochastically estimate the trace. This results in a substantial reduction of the variance, so that higher precision can be obtained with much less effort. We illustrate this for the trace of the inverse using three different classes of matrices.
References
[1]
C. Alexandrou, M. Constantinou, V. Drach, K. Hadjiyiannakou, K. Jansen, G. Koutsou, A. Strelchenko, and A. Vaquero, Evaluation of disconnected quark loops for hadron structure using GPUs, Comput. Phys. Commun., 185 (2014), pp. 1370--1382, https://doi.org/10.1016/j.cpc.2014.01.009.
[2]
H. Avron and S. Toledo, Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix, J. ACM, 58 (2011), 8.
[3]
R. Babich, J. Brannick, R. C. Brower, M. A. Clark, T. A. Manteuffel, S. F. McCormick, J. C. Osborn, and C. Rebbi, Adaptive multigrid algorithm for the lattice Wilson-Dirac operator, Phys. Rev. Lett., 105 (2010), 201602.
[4]
G. Bali, S. Collins, A. Frommer, K. Kahl, I. Kanamori, B. Müller, M. Rottmann, and J. Simeth, (Approximate) low-mode averaging with a new multigrid eigensolver, in 33rd International Symposium on Lattice Field Theory (LATTICE 2015), https://pos.sissa.it/251/350/, 2016.
[5]
S. Baral, T. Whyte, W. Wilcox, and R. B. Morgan, Disconnected loop subtraction methods in lattice QCD, Comput. Phys. Commun., 241 (2019), pp. 64--79.
[6]
C. Bekas, E. Kokiopoulou, and Y. Saad, An estimator for the diagonal of a matrix, Appl. Numer. Math., 57 (2007), pp. 1214--1229.
[7]
A. H. Bentbib, M. El Ghomari, K. Jbilou, and L. Reichel, Shifted extended global Lanczos processes for trace estimation with application to network analysis, Calcolo, 58 (2021), 4, https://doi.org/10.1007/s10092-020-00395-1.
[8]
D. Braess, Towards algebraic multigrid for elliptic problems of second order, Computing, 55 (1995), pp. 379--393, https://doi.org/10.1007/BF02238488.
[9]
M. Brezina, R. Falgout, S. MacLachlan, T. Manteuffel, S. McCormick, and J. Ruge, Adaptive smoothed aggregation ($\alpha$SA) multigrid, SIAM Rev., 47 (2005), pp. 317--346, https://doi.org/10.1137/050626272.
[10]
J. Chen and Y. Saad, A posteriori error estimate for computing ${\rm tr}(f(A))$ by using the Lanczos method, Numer. Linear Algebra Appl., 25 (2018), e2170, https://doi.org/10.1002/nla.2170.
[11]
A. Cortinovis and D. Kressner, On randomized trace estimates for indefinite matrices with an application to determinants, Found. Comput. Math., 22 (2022), pp. 875--903, https://doi.org/10.1007/s10208-021-09525-9.
[12]
T. A. DeGrand and S. Schaefer, Improving meson two point functions in lattice QCD, Comput. Phys. Commun., 159 (2004), pp. 185--191, https://doi.org/10.1016/j.cpc.2004.02.006.
[13]
S. Dong and K. Liu, Stochastic estimation with ${Z}_2$ noise, Phys. Lett. B, 328 (1994), pp. 130--136.
[14]
E. Endress, C. Pena, and K. Sivalingam, Variance reduction with practical all-to-all lattice propagators, Comput. Phys. Commun., 195 (2015), pp. 35--48, https://doi.org/10.1016/j.cpc.2015.04.017.
[15]
E. Estrada, The Structure of Complex Networks: Theory and Applications, Oxford University Press, New York, 2011.
[16]
E. Estrada and D. J. Higham, Network properties revealed through matrix functions, SIAM Rev., 52 (2010), pp. 696--714, https://doi.org/10.1137/090761070.
[17]
A. Frommer, K. Kahl, S. Krieg, B. Leder, and M. Rottmann, Adaptive aggregation-based domain decomposition multigrid for the lattice Wilson--Dirac operator, SIAM J. Sci. Comput., 36 (2014), pp. A1581--A1608, https://doi.org/10.1137/130919507.
[18]
A. Frommer, C. Schimmel, and M. Schweitzer, Analysis of probing techniques for sparse approximation and trace estimation of decaying matrix functions, SIAM J. Matrix Anal. Appl., 42 (2021), pp. 1290--1318, https://doi.org/10.1137/20M1364461.
[19]
A. S. Gambhir, A. Stathopoulos, and K. Orginos, Deflation as a method of variance reduction for estimating the trace of a matrix inverse, SIAM J. Sci. Comput., 39 (2017), pp. A532--A558, https://doi.org/10.1137/16M1066361.
[20]
M. B. Giles, Multilevel Monte Carlo methods, Acta Numer., 24 (2015), pp. 259--328, https://doi.org/10.1017/S096249291500001X.
[21]
Y. Ginosar, I. Gutman, T. Mansour, and M. Schork, Estrada index and Chebyshev polynomials, Chem. Phys. Lett., 454 (2008), pp. 145--147.
[22]
L. Giusti, P. Hernandez, M. Laine, P. Weisz, and H. Wittig, Low-energy couplings of QCD from current correlators near the chiral limit, J. High Energy Phys., 2004 (2004), 013, https://doi.org/10.1088/1126-6708/2004/04/013.
[23]
G. H. Golub, M. Heath, and G. Wahba, Generalized cross-validation as a method for choosing a good ridge parameter, Technometrics, 21 (1979), pp. 215--223.
[24]
G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins Stud. Math. Sci., Johns Hopkins University Press, Baltimore, MD, 2013.
[25]
G. H. Golub and U. von Matt, Generalized cross-validation for large scale problems, J. Comput. Graph. Statist., 6 (1995), pp. 1--34.
[26]
E. Hallman and D. Troester, A multilevel approach to stochastic trace estimation, Linear Algebra Appl., 638 (2022), pp. 125--149.
[27]
I. Han, D. Malioutov, H. Avron, and J. Shin, Approximating spectral sums of large-scale matrices using stochastic Chebyshev approximations, SIAM J. Sci. Comput., 39 (2017), pp. A1558--A1585, https://doi.org/10.1137/16M1078148.
[28]
I. Han, D. Malioutov, and J. Shin, Large-scale log-determinant computation through stochastic Chebyshev expansions, in Proceedings of the International Conference on Machine Learning, 2015, pp. 908--917.
[29]
N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, 2008, https://doi.org/10.1137/1.9780898717778.
[30]
M. F. Hutchinson, A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines, Comm. Statist. Simulation Comput., 19 (1990), pp. 433--450, https://doi.org/10.1080/03610919008812866.
[31]
J. Laeuchli and A. Stathopoulos, Extending hierarchical probing for computing the trace of matrix inverses, SIAM J. Sci. Comput., 42 (2020), pp. A1459--A1485, https://doi.org/10.1137/18M1176427.
[32]
L. Lin, Y. Saad, and C. Yang, Approximating spectral densities of large matrices, SIAM Rev., 58 (2016), pp. 34--65, https://doi.org/10.1137/130934283.
[33]
Q. Liu, W. Wilcox, and R. Morgan, Polynomial Subtraction Method for Disconnected Quark Loops, preprint, https://arxiv.org/abs/1405.1763, 2014.
[34]
R. A. Meyer, C. Musco, C. Musco, and D. P. Woodruff, Hutch++: Optimal stochastic trace estimation, in Symposium on Simplicity in Algorithms (SOSA), SIAM, Philadelphia, 2021, pp. 142--155, https://doi.org/10.1137/1.9781611976496.16.
[35]
L. N. Olson and J. B. Schroder, PyAMG: Algebraic Multigrid Solvers in Python v4.0, https://github.com/pyamg/pyamg, 2018.
[36]
C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning), The MIT Press, Cambridge, MA, 2005.
[37]
E. Romero, A. Stathopoulos, and K. Orginos, Multigrid deflation for lattice QCD, J. Comput. Phys., 409 (2020), 109356, https://doi.org/10.1016/j.jcp.2020.109356.
[38]
F. Roosta Khorasani and U. Ascher, Improved bounds on sample size for implicit matrix trace estimators, Found. Comput. Math., 15 (2015), pp. 1187--1212.
[39]
H. Rue and L. Held, Gaussian Markov Random Fields: Theory And Applications, CRC Press, Boca Raton, FL, 2005.
[40]
Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM, Philadelphia, 2003, https://doi.org/10.1137/1.9780898718003.
[41]
B. Sapoval, T. Gobron, and A. Margolina, Vibrations of fractal drums, Phys. Rev. Lett., 67 (1991), pp. 2974--2977.
[42]
J. Schwinger, Gauge invariance and mass II, Phys. Rev. (2), 128 (1962), p. 2425--2429, https://doi.org/10.1103/PhysRev.128.2425.
[43]
J. Sexton and D. Weingarten, Systematic Expansion for Full QCD Based on the Valence Approximation, preprint, 1994, https://arxiv.org/abs/hep-lat/9411029, 1994.
[44]
A. Stathopoulos, J. Laeuchli, and K. Orginos, Hierarchical probing for estimating the trace of the matrix inverse on toroidal lattices, SIAM J. Sci. Comput., 35 (2013), pp. S299--S322, https://doi.org/10.1137/120881452.
[45]
J. M. Tang and Y. Saad, A probing method for computing the diagonal of a matrix inverse, Numer. Linear Algebra Appl., 19 (2012), pp. 485--501.
[46]
U. Trottenberg, C. Oosterlee, and A. Schuller, Multigrid, Academic Press, New York, 2000.
[47]
S. Ubaru, J. Chen, and Y. Saad, Fast estimation of $\mathrm{tr}(f({A}))$ via stochastic Lanczos quadrature, SIAM J. Matrix Anal. Appl., 38 (2017), pp. 1075--1099, https://doi.org/10.1137/16M1104974.
[48]
S. Ubaru and Y. Saad, Applications of trace estimation techniques, in International Conference on High Performance Computing in Science and Engineering, Springer, Cham, 2017, pp. 19--33.
[49]
W. Wilcox, Noise Methods for Flavor Singlet Quantities, preprint, https://arxiv.org/abs/hep-lat/9911013v2, 1999.
Recommendations
A Multilevel Dual Reordering Strategy for Robust Incomplete LU Factorization of Indefinite Matrices
A dual reordering strategy based on both threshold and graph reorderings is introduced to construct robust incomplete LU (ILU) factorization of indefinite matrices. The ILU matrix is constructed as a preconditioner for the original matrix to be used in ...
A sparse-sparse iteration for computing a sparse incomplete factorization of the inverse of an SPD matrix
In this paper, a method via sparse-sparse iteration for computing a sparse incomplete factorization of the inverse of a symmetric positive definite matrix is proposed. The resulting factorized sparse approximate inverse is used as a preconditioner for ...
"Wide or tall" and "sparse matrix dense matrix" multiplications
HPC '11: Proceedings of the 19th High Performance Computing SymposiaSparse matrix dense matrix (SMDM) multiplications are useful in block Krylov or block Lanczos methods. SMDM computations are AU, and VA, multiplication of a large sparse m x n matrix A by a matrix V of k rows of length m or a matrix U of k columns of ...
Comments
Please enable JavaScript to view thecomments powered by Disqus.Information & Contributors
Information
Published In
© 2022, Society for Industrial and Applied Mathematics.
Publisher
Society for Industrial and Applied Mathematics
United States
Publication History
Published: 01 January 2022
Author Tags
Qualifiers
- Research-article
Contributors
Other Metrics
Bibliometrics & Citations
Bibliometrics
Article Metrics
- 0Total Citations
- 0Total Downloads
- Downloads (Last 12 months)0
- Downloads (Last 6 weeks)0
Reflects downloads up to 13 Jan 2025