Monte Carlo Methods for Estimating the Diagonal of a Real Symmetric Matrix

Abstract

For real symmetric matrices that are accessible only through matrix vector products, we present Monte Carlo estimators for computing the diagonal elements. Our probabilistic bounds for normwise absolute and relative errors apply to Monte Carlo estimators based on random Rademacher, sparse Rademacher, normalized and unnormalized Gaussian vectors, and to vectors with bounded fourth moments. The novel use of matrix concentration inequalities in our proofs represents a systematic model for future analyses. Our bounds mostly do not depend on the matrix dimension, target different error measures than existing work, and imply that the accuracy of the estimators increases with the diagonal dominance of the matrix. An application to derivative-based global sensitivity metrics corroborates this, as do numerical experiments on synthetic test matrices. We recommend against the use in practice of sparse Rademacher vectors, which are the basis for many randomized sketching and sampling algorithms, because they tend to deliver barely a digit of accuracy even under large sampling amounts.

Approximate Euclidean lengths and distances beyond Johnson-Lindenstrauss

“…Recently, in their seminal work, Meyer, Musco, Musco, and Woodruff [25] proved a remarkable result: their Hutch++ algorithm is the first to obtain -accuracy for stochastic trace estimation while requiring only 1/ matrix-vector queries. For the related problem of estimating the diagonal elements of a matrix, which was also recently studied in depth [18,5], Baston and Nakatsukasa [5] achieved -accuracy for the total, norm-wise error of the entire diagonal using O(1/ ) matrix-vector queries, but not for each individual diagonal element, which should not be possible due to the optimality of the JL lemma [23]. It is worth noting that the squared row norms of a matrix A can be found in the diagonal of A A, therefore, our work is closely connected the diagonal estimation literature.…”