Stochastic diagonal estimation: probabilistic bounds and an improved algorithm

TL;DR

The paper develops probabilistic guarantees for stochastic diagonal estimation of a matrix A accessed via Av queries. It analyzes two estimators, Gaussian and Rademacher, deriving (ε,δ) bounds for both per-entry and full-diagonal accuracy, and demonstrates how spectrum and eigenvector structure influence performance. It shows that Rademacher estimators typically outperform Gaussian ones for small numbers of queries, and connects diagonal estimation to trace-based methods, including Hutchinson’s estimator. Building on these insights, the authors adapt Hutch++ into Diag++, achieving a transition from O(1/ε^2) to O(1/ε) convergence for diagonal estimation in PSD matrices, with robustness to spectrum. The work combines theory and experiments to provide practical guidance for accelerated, spectrum-aware diagonal estimation in large-scale numerical linear algebra settings.

Abstract

We study the problem of estimating the diagonal of an implicitly given matrix $A$. For such a matrix we have access to an oracle that allows us to evaluate the matrix vector product $Av$. For random variable $v$ drawn from an appropriate distribution, this may be used to return an estimate of the diagonal of the matrix $A$. Whilst results exist for probabilistic guarantees relating to the error of estimates of the trace of $A$, no such results have yet been derived for the diagonal. We analyse the number of queries $s$ required to guarantee that with probability at least $1-δ$ the estimates of the relative error of the diagonal entries is at most $\varepsilon$. We extend this analysis to the 2-norm of the difference between the estimate and the diagonal of $A$. We prove, discuss and experiment with bounds on the number of queries $s$ required to guarantee a probabilistic bound on the estimates of the diagonal by employing Rademacher and Gaussian random variables. Two sufficient upper bounds on the minimum number of query vectors are proved, extending the work of Avron and Toledo [JACM 58(2)8, 2011], and later work of Roosta-Khorasani and Ascher [FoCM 15, 1187-1212, 2015]. We find that, generally, there is little difference between the two, with convergence going as $O(\log(1/δ)/\varepsilon^2)$ for individual diagonal elements. However for small $s$, we find that the Rademacher estimator is superior. These results allow us to then extend the ideas of Meyer, Musco, Musco and Woodruff [SOSA, 142-155, 2021], suggesting algorithm Diag++, to speed up the convergence of diagonal estimation from $O(1/\varepsilon^2)$ to $O(1/\varepsilon)$ and make it robust to the spectrum of any positive semi-definite matrix $A$.

Figures (5)

• Figure 1: The convergence of Rademacher estimates of selected diagonal elements of the matrix "Boeing/msc10480" from the SuiteSparse matrix collection. The numerical errors are shown in solid, and the theoretical bounds in dash. For $100$ trials, we plot the median error (top left), the $67^{th}$ percentile error (top right), the $80^{th}$ percentile error (bottom left) and the $90^{th}$ percentile error (bottom right), along with the corresponding theoretical error bounds for the associated value of $\delta$. The constants $(\|\bm{A}_i\|_2^2 - A_{ii}^2)/A_{ii}^2$ associated with each index are found as $0.3868$, $0.0125$ and $4.7154$ for $i$ equal to $1$, $10$ and $1000$ respectively.
• Figure 2: Comparison of the convergence of the Rademacher and Gaussian diagonal estimators. We run 100 trials for each distribution and report the median error and the $90^{th}$ percentile error in each case, along with their respective theoretical bounds.
• Figure 3: Convergence of the Rademacher and Gaussian diagonal estimators to the actual diagonal of the matrix "Boeing/msc10480" from the SuiteSparse matrix collection. The effect of the infinite variance for the Gaussian estimator when $s\leq2$ is much more pronounced.
• Figure 4: Convergence of the Rademacher estimator on random matrices. We run $50$ trials and report the median and $90^{th}$ percentile errors for three values of $c,$$0.5$, $1$ and $1.5$, corresponding to slow, medium and fast eigenvalue decay respectively. As expected, with all matrices, convergence goes as $O(1/\varepsilon^2)$. The difference between the matrices arises from the fact that the values of $(\|A\|_F^2 - \|\bm{A}_d\|_2^2)/\|\bm{A}_d\|_2^2$ increase as $c$ increases.
• Figure 5: Relative error versus number of matrix-vector queries: we report the median and $90^{th}$ percentile after 10 trials for Synthetic Matrices with $c = 0.5,\: 1,\: 1.5$. Whilst pure stochastic estimation outperforms Diag++ for slower eigenvalue decay, Diag++ converges at roughly the same rate. Estimation through projection only however, is poor. Meanwhile Diag++ obtains high performance for faster eigenvalue decay. It is clear that this algorithm is robust to the spectrum of the matrix in question, converging for both steep and flat eigenvalue distributions.

Theorems & Definitions (27)

• lemma thmcounterlemma
• definition thmcounterdefinition
• definition thmcounterdefinition
• definition thmcounterdefinition
• definition thmcounterdefinition
• definition thmcounterdefinition
• theorem 1
• proof
• corollary thmcountercorollary
• proof