Matrix function trace approximation - Randomized Numerical Linear Algebra with Examples

Let $\vec{A}\in\R^{n\times n}$ be a symmetric matrix with eigendecomposition $\vec{A} = \sum_{i=1}^{n} \lambda_i \vec{u}_i \vec{u}_i^\T$ , where $\lambda_i$ are the eigenvalues and $\vec{u}_i$ are the orthonormal eigenvectors of $\vec{A}$ . Recall the matrix function

f(\vec{A}) := \sum_{i=1}^{n} f(\lambda_i) \vec{u}_i \vec{u}_i^\T.

(6.9)

A natural approach is to combine the implicit trace estimation algorithms discussed earlier in this chapter with black-box methods for approximating $\vec{x}\mapsto \vec{x}^\T f(\vec{A})\vec{x}$ or $\vec{x}\mapsto f(\vec{A})\vec{x}$ , such as the Lanczos method for matrix functions.

Stochastic Lanczos Quadrature¶

Recall that using $k-1$ matrix-vector products with $\vec{A}$ , the Lanczos method for quadratic forms produces an approximation $\Call{Lan-QF}_k(f;\vec{A},\vec{x})$ that is an approximation to $\vec{x}^\T f(\vec{A})\vec{x}$ . Stochastic Lanczos Quadrature (SLQ) simply combines this method with the Girard-Hutchinson trace estimator.

A simple application of the triangle inequality gives a bound on expected squared error of the SLQ estimator.

Proof

By the triangle inequality and since $(x+y)^2\leq 2(x^2+y^2)$ , we have

\begin{aligned} &| \tr(f(\vec{A})) - \Call{SLQ}_{k,m}(f;\vec{A}) |^2 \\&\hspace{4em}\leq 2\left| \tr(f(\vec{A})) - \frac{1}{m}\sum_{i=1}^{m} \vec{x}_i^\T f(\vec{A})\vec{x}_i \right|^2 \\&\hspace{5em}+ 2\left| \frac{1}{m}\sum_{i=1}^{m} \left( \vec{x}_i^\T f(\vec{A})\vec{x}_i - \Call{Lan-QF}_k(f;\vec{A},\vec{x}_i) \right) \right|^2. \end{aligned}

Note that

\EE\left[ \left| \tr(f(\vec{A})) - \frac{1}{m}\sum_{i=1}^{m} \vec{x}_i^\T f(\vec{A})\vec{x}_i \right|^2 \right] = \VV\left[ \widehat{\tr}_m(f(\vec{A})) \right] = \frac{2 \| f(\vec{A}) \|_\F^2}{m},

where $\widehat{\tr}_m(\cdot)$ is the Girard--Hutchinson estimator.

Next, by the triangle inequality and Theorem 1.3,

\begin{aligned} \hspace{4em}&\hspace{-4em} \left| \frac{1}{m}\sum_{i=1}^{m} \left( \vec{x}_i^\T f(\vec{A})\vec{x}_i - \Call{Lan-QF}_k(f;\vec{A},\vec{x}_i) \right) \right| \\&\leq \frac{1}{m}\sum_{i=1}^{m} \left| \vec{x}_i^\T f(\vec{A})\vec{x}_i - \Call{Lan-QF}_k(f;\vec{A},\vec{x}_i) \right| \\&\leq \frac{2}{m} \sum_{i=1}^{m} \|\vec{x}_i\|^2 \min_{\deg(p)<2k-1} \max_{x\in[\lambda_n,\lambda_1]} | f(x) - p(x) |. \end{aligned}

Hence,

\begin{aligned} \hspace{4em}&\hspace{-4em} \EE\left[ \left| \frac{1}{m}\sum_{i=1}^{m} \left( \vec{x}_i^\T f(\vec{A})\vec{x}_i - \Call{Lan-QF}_k(f;\vec{A},\vec{x}_i) \right) \right|^2 \right] \\&\leq \frac{4}{m^2}\EE\left[ \left(\sum_{i=1}^{m} \|\vec{x}_i\|^2 \right)^2 \right] \min_{\deg(p)<2k-1} \max_{x\in[\lambda_n,\lambda_1]} | f(x) - p(x) |^2 \end{aligned}

Now, note that $\sum_{i=1}^{m} \|\vec{x}_i\|^2$ is a Chi-squared random variable with $mn$ degrees of freedom. By looking on Wikipedia, we see that

\EE\left[ \left(\sum_{i=1}^{m} \|\vec{x}_i\|^2 \right)^2 \right] = 2mn + (mn)^2 \leq 3(mn)^2.

Combining all of the above gives the result.

Note that for “nice” functions $f(x)$ , the error of the best polynomial approximation decreases exponentially with $k$ Trefethen, 2019.

Variance reduction¶

As with implicit trace estimation, the accuracy of SLQ can be improved using a control variate. In particular, suppose $\widehat{\vec{F}}$ is a low-rank approximation to $f(\vec{A})$ whose trace we can compute exactly. Since

\tr(f(\vec{A})) = \tr(\widehat{\vec{F}}) + \tr(f(\vec{A}) - \widehat{\vec{F}}),

(6.10)

we can estimate the second term with SLQ; i.e. we use the estimator

\tr(\widehat{\vec{F}}) + \frac{1}{m}\sum_{i=1}^{m} \left( \Call{Lan-QF}_k(f;\vec{A},\vec{x}_i) - \vec{x}_i^\T\widehat{\vec{F}}\vec{x}_i \right).

(6.11)

The variance of this estimator scales with $\|f(\vec{A}) - \widehat{\vec{F}}\|_\F^2$ rather than $\|f(\vec{A})\|_\F^2$ , so it can be substantially more accurate when $\widehat{\vec{F}}$ is a good approximation to $f(\vec{A})$ .

Techniques for efficiently obtaining a low-rank approximation to $f(\vec{A})$ are described in Matrix function low-rank approximation. In particular, Krylov-aware methods Chen & Hallman, 2023Persson et al., 2025 were developed with this application in mind, and for operator monotone functions, a suitable $\widehat{\vec{F}}$ can be obtained using products with $\vec{A}$ only.

References¶

Chen, T. (2024). The Lanczos algorithm for matrix functions: a handbook for scientists. https://arxiv.org/abs/2410.11090
Chen, T., Trogdon, T., & Ubaru, S. (2025). Randomized Matrix-Free Quadrature: Unified and Uniform Bounds for Stochastic Lanczos Quadrature and the Kernel Polynomial Method. SIAM Journal on Scientific Computing, 47(3), A1733–A1757. 10.1137/23m1600414
Trefethen, L. N. (2019). Approximation Theory and Approximation Practice, Extended Edition. Society for Industrial. 10.1137/1.9781611975949
Chen, T., & Hallman, E. (2023). Krylov-Aware Stochastic Trace Estimation. SIAM Journal on Matrix Analysis and Applications, 44(3), 1218–1244. 10.1137/22m1494257
Persson, D., Chen, T., & Musco, C. (2025). Randomized block-Krylov subspace methods for low-rank approximation of matrix functions. https://arxiv.org/abs/2502.01888

6.3 Matrix function trace approximation

Stochastic Lanczos Quadrature¶

Variance reduction¶