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

TL;DR

This paper introduces a cost-effective method for spectral density approximation using modified stochastic trace estimators with fixed randomization, enabling reuse of computations across parameter-dependent matrices without sacrificing accuracy.

Contribution

It proposes and analyzes modifications to three stochastic trace estimators that allow reuse of random vectors across different parameter values, reducing computational costs in spectral density estimation.

Findings

Constant randomization does not deteriorate estimator accuracy.

Nystr"om++ with Chebyshev approximation achieves error bounds with minimal matrix-vector products.

Numerical experiments validate the efficiency and accuracy of the proposed method.

Abstract

Stochastic trace estimation is a well-established tool for approximating the trace of a large symmetric matrix $\boldsymbol{B}$. Several applications involve a matrix that depends continuously on a parameter $t \in [a,b]$, and require trace estimates of $\boldsymbol{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 $\boldsymbol{B}(t_1), \dots, \boldsymbol{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 fixed randomization across different values of $t$, that is, every matrix $\boldsymbol{B}(t_1), \dots, \boldsymbol{B}(t_m)$ is multiplied with the 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 $\boldsymbol{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 and neural network optimization validate our findings.