Randomized estimation of spectral densities of large matrices made accurate

Abstract

Abstract. For a large Hermitian matrix A ∈ C N ×N , it is often the case that the only affordable operation is matrix-vector multiplication. In such case, randomized method is a powerful way to estimate the spectral density (or density of states) of A. However, randomized methods developed so far for estimating spectral densities only extract information from different random vectors independently, and the accuracy is therefore inherently limitedwhere Nv is the number of random vectors.In this paper we demonstrate that the "O(1/ √ Nv) barrier" can be overcome by taking advantage of the correlated information of random vectors when properly filtered by polynomials of A. Our method uses the fact that the estimation of the spectral density essentially requires the computation of the trace of a series of matrix functions that are numerically low rank. By repeatedly applying A to the same set of random vectors and taking different linear combination of the results, we can sweep through the entire spectrum of A by building such low rank decomposition at different parts of the spectrum. Under some assumptions, we demonstrate that a robust and efficient implementation of such spectrum sweeping method can compute the spectral density accurately with O(N 2 ) computational cost and O(N ) memory cost. Numerical results indicate that the new method can significantly outperform existing randomized methods in terms of accuracy. As an application, we demonstrate a way to accurately compute a trace of a smooth matrix function, by carefully balancing the smoothness of the integrand and the regularized density of states using a deconvolution procedure.

Randomized matrix-free trace and log-determinant estimators

“…Like our estimators, the spectrum-sweeping method [27,Algorithm 5] is based on a randomized low-rank approximation of A. However, it is designed to compute the trace of smooth functions of Hermitian matrices in the context of density of state estimations in quantum physics.…”

Fast Computation of Spectral Densities for Generalized Eigenvalue Problems

“…It has been proved in [19] that the expansion error in (3.6) decays as ρ −m for some constant ρ > 1. For a general matrix A whose eigenvalues are not necessarily in the interval [−1, 1], a linear transformation is first applied to A to bring its eigenvalues to the desired interval.…”

“…In order to compare with the accuracy of the DOS, the exact eigenvalues of each problem are computed with MATLAB built-in function eig. We measure the error of the approximate DOS using the relative L 1 error as proposed in [19]:…”

A posteriori error estimate for computing tr(f(A)) by using the Lanczos method

“…Computing the trace in this manner, however, requires the computation of all the eigenvalues, which is also often prohibitively expensive. Hence, various methods proposed for approximately computing tr( f ( A )) consist of the following two ingredients: Approximate the trace of f ( A ) by using the average of unbiased samples , i =1,…, N , where the u i are independent random vectors of some nature. Approximately compute the bilinear form by using some numerical technique. …”

“…Computing the trace in this manner, however, requires the computation of all the eigenvalues, which is also often prohibitively expensive. Hence, various methods proposed for approximately computing tr( f(A)) consist of the following two ingredients 1,10,[13][14][15][16][17][18][19] The various methods differ in the random mechanism of selecting the u i and the numerical technique for computing the bilinear form. Several variants of these ingredients exist (e.g., computing deterministically tr( (A)) = ∑ n i=1 e T i (A)e i rather than using random vectors u i , or even using block vectors to replace the canonical vectors e i 20 ; or using moment extrapolation for, particularly, f(t) = t with real value 21,22 ), but they are not the focus of this work.…”