Stochastic trace estimation for parameter-dependent matrices applied to spectral density approximation

Abstract: Stochastic trace estimation is a well-established tool for approximating the trace of a large symmetric matrix $\mathbf{B}$. Several applications involve a matrix that depends continuously on a parameter $t \in [a,b]$, and require trace estimates of $\mathbf{B}(t)$ for many values of $t$. This is, for example, the case when approximating the spectral density of a matrix. Approximating the trace separately for each matrix $\mathbf{B}(t_1), \dots, \mathbf{B}(t_m)$ clearly incurs redundancies and a cost that scales linearly with $m$. To address this issue, we propose and analyze modifications for three stochastic trace estimators, the Girard-Hutchinson, Nystr\"om, and Nystr\"om++ estimators. Our modification uses \emph{constant} randomization across different values of $t$, that is, every matrix $\mathbf{B}(t_1), \dots, \mathbf{B}(t_m)$ is multiplied with the \emph{same} set of random vectors. When combined with Chebyshev approximation in $t$, the use of such constant random matrices allows one to reuse matrix-vector products across different values of $t$, leading to significant cost reduction. Our analysis shows that the loss of stochastic independence across different $t$ does not lead to deterioration. In particular, we show that $\mathcal{O}(\varepsilon^{-1})$ random matrix-vector products suffice to ensure an error of $\varepsilon > 0$ for Nystr\"om++, independent of low-rank properties of $\mathbf{B}(t)$. We discuss in detail how the combination of Nystr\"om++ with Chebyshev approximation applies to spectral density estimation and provide an analysis of the resulting method. This improves various aspects of an existing stochastic estimator for spectral density estimation. Several numerical experiments from electronic structure interaction, statistical thermodynamics, and neural network optimization validate our findings.