Randomized matrix-free quadrature for spectrum and spectral sum approximation

Tyler Chen

This is a companion piece to the paper: https://arxiv.org/abs/2204.01941

Code to replicate the figures in this document can be found on Github.

Introduction

The cumulative empirical spectral measure (CESM) $\Phi : \mathbb{R}\to [0,1]$, induced by $\mathbf{A}$, is defined by $$\Phi(x) = \Phi_{\mathbf{A}} := \sum_{i=0}^{n-1} n^{-1} \mathbb{I}[ \lambda_i \leq x].$$ Here $\mathbb{I}[ \texttt{true} ] =1$ and $\mathbb{I}[\texttt{false}] = 0$. Traces of matrix functions are related to the CESM by $$\operatorname{tr}(f(\mathbf{A})) = n \int f\,\mathrm{d}\Phi.$$ You can read more about applications of these quantities here.

In this paper, we provide a unified perspective of algorithms for randomized matrix-free quadrature. Specifically, the class of algorithms we study is characterized by the use of a Krylov subspace method to approximate independent and identically distributed samples of $\mathbf{v}^\mathsf{T}f(\mathbf{A})\mathbf{v}$, where $\mathbf{v}$ is an isotropic random vector. Since many widely used algorithms, such as stochastic Lanczos quadrature (SLQ) [1] and the kernel polynomial method (KPM) [2], fall into this class, our unified perspective allows the algorithms’ tradeoffs to be more clearly understood.

Define the weighted CESM $\Psi : \mathbb{R}\to[0,1]$, an estimator for the CESM defined by $$\Psi(x) := \mathbf{v}^\mathsf{T}\mathbb{I}(\mathbf{A}\leq x) \mathbf{v} = \sum_{i=1}^{n} |\mathbf{v}^\mathsf{T}\mathbf{u}_i|^2 \mathbb{I}(\lambda_i\leq x),$$ where $\mathbf{v}$ is a random vector satisfying $\mathbb{E}[\mathbf{v}\mathbf{v}^\mathsf{T}] = n^{-1} \mathbf{I}$; i.e., $\mathbf{v}$ is isotropic with appropriate scale.

Then, using basic properties of the trace and expectation, we have $$\begin{align*} \mathbb{E}[\Psi(x)] &= \mathbb{E}[\mathbf{v}^\mathsf{T}\mathbb{I}(\mathbf{A}\leq x) \mathbf{v}] = \mathbb{E}[\operatorname{tr}(\mathbf{v}^\mathsf{T}\mathbb{I}(\mathbf{A}\leq x) \mathbf{v})] \\& \hspace{3em} = \mathbb{E}[\operatorname{tr}(\mathbb{I}(\mathbf{A}\leq x) \mathbf{v} \mathbf{v}^\mathsf{T})] = \operatorname{tr}(\mathbb{E}[ \mathbb{I}(\mathbf{A}\leq x) \mathbf{v} \mathbf{v}^\mathsf{T}]) \\& \hspace{6em} = \operatorname{tr}( \mathbb{I}(\mathbf{A}\leq x) \mathbb{E}[ \mathbf{v} \mathbf{v}^\mathsf{T}] ) = n^{-1} \operatorname{tr}(\mathbb{I}(\mathbf{A}\leq x)) = \Phi(x). \end{align*}$$ Further, almost by definition, we see that $\int_{-\infty}^{\infty} f(x) \mathrm{d}\Psi(x) = \mathbf{v}^\mathsf{T}f(\mathbf{A})\mathbf{v}$ is an unbiased estimator for $n^{-1} \operatorname{tr}(f(\mathbf{A}))$.

The Lanczos method for matrix function approximation (Lanczos-FA) can be used to approximate $\mathbf{v}^\mathsf{T}f(\mathbf{A})\mathbf{v}$ (read more about this task here).

The prototypical algorithm

A natural approach to approximate the CESM $\Psi$ and therefore $\operatorname{tr}(f(\mathbf{A}))$ is to average several copies of $\Psi$. However, computing $\Psi$ requires knowledge of the eigenvalues of $\mathbf{A}$. Therefore, we can approximate $\Psi$ with a quadrature rule; i.e. with a simpler distribution function. The most common way to do this is by obtaining the moments of $\Psi$, and using this information to ensure that polynomials up to some degree are integrated exactly.

This is typically done with a Krylov subspace method. Indeed, for any $i,j$ with $0\leq i,j \leq k$, $$(\mathbf{A}^j\mathbf{v})^\mathsf{T}(\mathbf{A}^j\mathbf{v}) = \mathbf{v}^\mathsf{T}\mathbf{A}^{i+j} \mathbf{v} = \int_{-\infty}^{\infty} x^{i+j} \mathrm{d}\Psi(x).$$ Thus, the moments of $\Psi$ through degree $2k$ can be computed using the information in the Krylov subspace $$\mathcal{K}_{k}(\mathbf{A},\mathbf{v}) := \operatorname{span}\{\mathbf{v}, \mathbf{A}\mathbf{v}, \ldots, \mathbf{A}^{k}\mathbf{v} \} = \{ p(\mathbf{A})\mathbf{v} : \deg(p) \leq k \}.$$

This can be summarized with the following algorithm:

Algorithm. Prototypical randomized matrix free quadrature $\begin{aligned} \small{1:}&\hspace{1em} \textsf{protoalg}(\mathbf{A}, n_{\text{v}}, k, \circ) \\ \small{2:}&\hspace{1em} \textbf{for}~ \ell = 1, 2, \ldots, n_{\textup{v}}~ \textbf{do} \\ \small{3:}&\hspace{2em} \text{define (implicitly) \( \Psi_\ell \stackrel{\text{iid}}{\sim} \Psi \) by sampling \( \mathbf{v}_\ell \stackrel{\text{iid}}{\sim} \mathbf{v} \), \( \mathbb{E}[\mathbf{v}\mathbf{v}^\mathsf{T}] = n^{-1} \mathbf{I} \)} \\ \small{4:}&\hspace{2em} \text{compute the moments of \( \Psi_{\ell} \) through degree \( s \) by constructing \( \mathcal{K}_k(\mathbf{A},\mathbf{v}_\ell) \)} \\ \small{5:}&\hspace{2em} \text{approximate \( \Psi_{\ell} \) by a quadrature rule} \\ \small{6:}&\hspace{1em}\textbf{return}~\text{average of individual quadrature rules} \end{aligned}$

A unified perspective

Common algorithms like SLQ and KPM are special instances of the prototypical algorithm described above. As such, they can all be analyzed through the lens of randomized quadrature. Our analysis shows that all algorithms of this form satisfy similar theoretical bounds. However there are many practical concerns which separate the algorithms.

A number of important past works have argued that Lanczos-based approaches are unstable, and therefore should not be used without costly reorthogonalization schemes . Based on the analysis in [3] , we argue that, even without reorthogonalzation, Lanczos-based approaches tend to outperform approaches based on explicit polynomials. In fact, we go so far as to argue that explicit polynomial methods such as the kernel polynomial method should typically be implemented using Lanczos. The qualitative tradeoffs between the algorithms are further illustrated by our numerical experiments in a wite variety of numerical experiments.

References