remarkRemark hypothesisHypothesis claimClaim Faster Stochastic Trace
EstimationE. Hallman ex_supplement

# Faster Stochastic Trace Estimation with a Chebyshev Product
Identity††thanks:  This research was supported in part by the National Science
Foundation through grant DMS-1745654.

Eric Hallman North Carolina State University (,
https://erhallma.math.ncsu.edu/<EMAIL_ADDRESS>
###### Abstract

Methods for stochastic trace estimation often require the repeated evaluation
of expressions of the form $z^{T}p_{n}(A)z$, where $A$ is a symmetric matrix
and $p_{n}$ is a degree $n$ polynomial written in the standard or Chebyshev
basis. We show how to evaluate these expressions using only $\lceil n/2\rceil$
matrix-vector products, thus substantially reducing the cost of existing trace
estimation algorithms that use Chebyshev interpolation or Taylor series.

###### keywords:

Chebyshev polynomials, stochastic trace estimation, spectral function,
Hutchinson’s method

15A63, 68W20, 68W25

##  1 Introduction

Given a symmetric matrix $A\in\mathbb{R}^{d\times d}$ and a function
$f:\mathbb{R}\rightarrow\mathbb{R}$, we consider the problem of estimating

(1) $\operatorname*{tr}f(A)=\sum_{i=1}^{d}f(\lambda_{i}),$

where $\lambda_{1},\ldots,\lambda_{d}$ are the eigenvalues of $A$. When $A$ is
large enough to make computing its eigenvalues impractical, one common
approach is to use Hutchinson’s method [9]. This method samples $m$
independent vectors $\\{z^{(i)}\\}_{i=1}^{m}$ from a Rademacher distribution
(entries $\pm 1$ with equal probability), and yields the estimate

(2)
$\operatorname*{tr}f(A)\approx\frac{1}{m}\sum_{i=1}^{m}\left(z^{(i)}\right)^{T}f(A)z_{i}.$

For many functions $f$ of interest (e.g., $\exp(x)$, $x^{-1}$, $x^{p/2}$, or
$\log x$), the right hand side of (2) is further simplified by approximating
$f$ by a degree $n$ polynomial $p_{n}$, most commonly through Chebyshev
interpolation or a Taylor series.

A common way to evaluate an expression of the form $z^{T}p_{n}(A)z$ is to
compute $z_{n}=p_{n}(A)z$ and return $z^{T}z_{n}$ (see [7, 6, 2] for
examples). We refer to this method as one-sided evaluation, and it will in
general require $n$ matrix-vector products (matvecs) with $A$. For polynomials
written in the standard or Chebyshev bases we show how to reduce the number of
matvecs to $\lceil n/2\rceil$. Since the cost of many existing trace
estimation algorithms is dominated by matvecs, our method will cut their
runtime approximately in half.

###  1.1 Related work

A related idea is explored in [4] where the authors exploit the symmetry of
$A$ to reduce the cost of their estimator. Our proposal is more efficient and
is more generally applicable, as their method applies only when $f(A)$ is
positive definite.

For more general background on computing matrix polynomials, see [8, Ch. 4] or
[5, Sec. 9.2]. Our method bears some resemblance to that of Paterson and
Stockmeyer [11], but aims to compute $z^{T}p_{n}(A)z$ rather than $p_{n}(A)$
itself.

A few recent papers that use Chebyshev approximations for stochastic trace
estimation are [7, 3, 6, 4], and Taylor series are used similarly in [2]. For
applications of stochastic trace estimation, see [14].

One primary competitor to Chebyshev interpolation is stochastic Lanczos
quadrature. For more information on this method and its pros and cons with
respect to using Chebyshev polynomials, see [13]. In short, the authors
suggest that Lanczos quadrature is generally superior since it converges at
twice the rate of Chebyshev interpolation. If this is true, then our methods
(which do not extend to Lanczos quadrature) should put Chebyshev interpolation
back on more or less equal footing. It may require an approximating polynomial
with twice the degree of that needed by Lanczos, but can compute it with the
same number of matvecs!

##  2 Standard basis

It was noted in [4] that expressions of the form $z^{T}A^{n}z$ can be
evaluated with $\lceil n/2\rceil$ matvecs by letting $k=\lfloor n/2\rfloor$
and computing $z_{k}=A^{k}z$, then returning $z_{k}^{T}z_{k}$ if $n$ is even
and $z_{k}^{T}Az_{k}$ if $n$ is odd. We first extend this idea to polynomials
of the form $p_{n}(x)=\sum_{j=0}^{n}\alpha_{j}x^{j}$. Algorithm 1 requires
$\lceil n/2\rceil$ matvecs and at each step $j$ needs to store only the two
most recent vectors $z_{j}$, $z_{j-1}$ in memory.

Algorithm 1 Two-sided evaluation (standard basis)

0: Symmetric $A\in\mathbb{R}^{d\times d}$, $z_{0}\in\mathbb{R}^{d}$,
polynomial coefficients $a=[\alpha_{0},\alpha_{1},\ldots,\alpha_{n}]$

0: $s=z_{0}^{T}p_{n}(A)z_{0}$

1: $s=\alpha_{0}z_{0}^{T}z_{0}$

2: for $j=1,2,\ldots,\lceil n/2\rceil$ do

3: $z_{j}=Az_{i-1}${$z_{j}=A^{j}z_{0}$}

4: $s=s+\alpha_{2j-1}z_{j-1}^{T}z_{i}$

5: if $n=2j-1$ then stop

6: $s=s+\alpha_{2j}z_{j}^{T}z_{j}$

7: end for

##  3 Chebyshev basis

We use Chebyshev polynomials of the first kind, which can be defined by the
recurrence

(3) $T_{j+1}(x)=2xT_{j}(x)-T_{j-1}(x),$

where $T_{0}(x)=1$ and $T_{1}(x)=x$. A function
$f:[-1,1]\rightarrow\mathbb{R}$ can then be approximated by a polynomial of
the form

(4) $f(x)\approx p_{n}(x)=\sum_{j=0}^{n}\alpha_{j}T_{j}(x).$

The polynomial $p_{n}$ interpolates $f$ at a set of Chebyshev nodes
$\\{x_{j}\\}_{j=0}^{n}$. Several different choices for the nodes are available
[12, 10], but as one example Trefethen [12] uses the nodes

(5) $x_{j}=\cos\frac{j\pi}{n},\quad 0\leq j\leq n$

and shows how to quickly compute the coefficients $\\{\alpha_{j}\\}_{j=0}^{n}$
by using an FFT. Our method applies as long as the polynomial in (4) is
expressed in the Chebyshev basis.

The key idea is to use the fact that Chebyshev polynomials follow the relation
[1]

(6)
$T_{j}(x)T_{k}(x)=\frac{1}{2}\left(T_{j+k}(x)+T_{|k-j|}(x)\right),\quad\forall
j,k\geq 0.$

By letting $k=j$ or $k=j+1$ in the above equation and rearranging, it follows
that for all $j\geq 0$,

(7) $T_{2j}(x)=2T_{j}(x)^{2}-T_{0}(x)=2T_{j}(x)^{2}-1$

and

(8) $T_{2j+1}(x)=2T_{j}(x)T_{j+1}(x)-T_{1}(x)=2T_{j}(x)T_{j+1}(x)-x.$

We can therefore evaluate terms of the form $z^{T}T_{2j}(A)z$ by computing
$z_{j}=T_{j}(A)z$ and returning $2z_{j}^{T}z_{j}-z^{T}z$. Similarly, we can
evaluate terms of the form $z^{T}T_{2j+1}(A)z$ by computing
$z_{j+1}=T_{j+1}(A)z$ and returning $2z_{j}^{T}z_{j+1}-z^{T}Az$.

Our method is presented in Algorithm 2. It requires $\lceil n/2\rceil$ matvecs
and at each step $j$ needs to store only the two most recent vectors $z_{j}$,
$z_{j-1}$ in memory. It should therefore take about half the time required by
one-sided evaluation.

Algorithm 2 Two-sided evaluation (Chebyshev basis)

0: Symmetric $A\in\mathbb{R}^{d\times d}$, $z_{0}\in\mathbb{R}^{d}$, Chebyshev
coefficients $a=[\alpha_{0},\alpha_{1},\ldots,\alpha_{n}]$

0: $s=z_{0}^{T}p_{n}(A)z_{0}$

1: $z_{1}=Az_{0}$

2: $\zeta_{0}=z_{0}^{T}z_{0}$

3: $\zeta_{1}=z_{0}^{T}z_{1}$

4:
$s=\alpha_{0}\zeta_{0}+\alpha_{1}\zeta_{1}+\alpha_{2}\left(2(z_{1}^{T}z_{1})-\zeta_{0}\right)$

5: for $j=2,3,\ldots,\lceil n/2\rceil$ do

6: $z_{j}=2(Az_{i-1})-z_{j-2}${$z_{j}=T_{j}(A)z_{0}$}

7: $s=s+\alpha_{2j-1}\left(2(z_{j-1}^{T}z_{j})-\zeta_{1}\right)$

8: if $n=2j-1$ then stop

9: $s=s+\alpha_{2j}\left(2(z_{j}^{T}z_{j})-\zeta_{0}\right)$

10: end for

###  3.1 Stability

It is possible that due to cancellation and roundoff error, some of the
computed values $2z_{j-1}^{T}z_{j}-\zeta_{1}$ or $2z_{j}^{T}z_{j}-\zeta_{0}$
may have poor relative accuracy. The absolute error, however, will be small
(on the order of machine precision) compared to
$2\|z_{j-1}\|_{2}\|z_{j}\|_{2}+|\zeta_{1}|$ and
$2\|z_{j}\|_{2}\|z_{j}\|_{2}+|\zeta_{0}|$, respectively. Since the first two
terms in our sum are $\alpha_{0}\zeta_{0}$ and $\alpha_{1}\zeta_{1}$, we
expect that the effect of these rounding errors will typically be minor. A
more rigorous analysis is left for future exploration.

##  4 Numerical Experiments

We tested our algorithm on synthetic data in order to verify its correctness
and to examine how closely the output of Algorithm 2 aligned with that of one-
sided evaluation. We formed a random symmetric matrix
$A\in\mathbb{R}^{5000\times 5000}$ and scaled it so that its largest and
smallest eigenvalues were $1$ and $-1$. Using the function $f(x)=\exp(10x)$,
we computed $\operatorname*{tr}f(A)$ to be approximately $7.04\times 10^{5}$.

Using the Chebyshev nodes in (5) with $n=20$, we computed a Chebyshev
interpolant $p_{n}$ to $f$. For each of 100 trials, we drew a random
Rademacher vector $z$ and used both our method and one-sided evaluation to
compute $z^{T}p_{n}(A)z$. Our method took a total of $5.42$ seconds and one-
sided evaluation required $10.70$ seconds, so as predicted two-sided
evaluation ran nearly twice as fast.

Both methods estimated $\operatorname*{tr}f(A)$ to be about $7.14\times
10^{5}$. More notably, the two outputs agreed almost perfectly with each
other: the relative error between them was $1.6\times 10^{-16}$, less than
machine precision. We then examined the relative errors between the two
methods for each evaluation $z^{T}p_{n}(A)z$, and found that the largest was
$1.6\times 10^{-14}$.

To take an even closer look, we examined the values returned by the two
algorithms for each individual term $\alpha_{j}z^{T}T_{j}(A)z$ in each sum.
That is, we compared the terms $\alpha_{j}z_{0}^{T}z_{j}$ produced by one-
sided evaluation to the corresponding terms being added in Lines 7 or 9 in
Algorithm 2. The largest relative error found between any single pair of terms
was $3.1\times 10^{-11}$, so the two methods of evaluation agreed with each
other to a level of accuracy well beyond what would be required for stochastic
trace estimation.

##  5 Conclusion

We have shown how to evaluate the expression $z^{T}p_{n}(A)z$, where $A$ is
symmetric and $p_{n}$ is a polynomial in the standard or Chebyshev basis,
using no more than $\lceil n/2\rceil$ matvecs with $A$. Our proposed method is
simple to implement and can be used for any stochastic trace estimation
technique that relies on Taylor expansions or Chebyshev interpolation. The
stability of our method remains an open question, but a small numerical
experiment suggests that its output will largely agree with that of standard
one-sided evaluation. We therefore recommend that two-sided evaluation be
incorporated into existing algorithms.

## References

  * [1] M. Abramowitz, I. A. Stegun, and R. H. Romer, Handbook of mathematical functions with formulas, graphs, and mathematical tables, 1988.
  * [2] C. Boutsidis, P. Drineas, P. Kambadur, E.-M. Kontopoulou, and A. Zouzias, A randomized algorithm for approximating the log determinant of a symmetric positive definite matrix, Linear Algebra and its Applications, 533 (2017), pp. 95–117.
  * [3] E. Di Napoli, E. Polizzi, and Y. Saad, Efficient estimation of eigenvalue counts in an interval, Numerical Linear Algebra with Applications, 23 (2016), pp. 674–692.
  * [4] E. Dudley, A. K. Saibaba, and A. Alexanderian, Monte carlo estimators for the Schatten p-norm of symmetric positive semidefinite matrices, arXiv preprint arXiv:2005.10174, (2020).
  * [5] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 4th ed., 2013.
  * [6] I. Han, D. Malioutov, H. Avron, and J. Shin, Approximating spectral sums of large-scale matrices using stochastic Chebyshev approximations, SIAM Journal on Scientific Computing, 39 (2017), pp. A1558–A1585.
  * [7] I. Han, D. Malioutov, and J. Shin, Large-scale log-determinant computation through stochastic Chebyshev expansions, in International Conference on Machine Learning, 2015, pp. 908–917.
  * [8] N. J. Higham, Functions of matrices: theory and computation, SIAM, 2008\.
  * [9] M. F. Hutchinson, A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines, Communications in Statistics-Simulation and Computation, 18 (1989), pp. 1059–1076.
  * [10] S. E. Notaris, Interpolatory quadrature formulae with Chebyshev abscissae of the third or fourth kind, Journal of computational and applied mathematics, 81 (1997), pp. 83–99.
  * [11] M. S. Paterson and L. J. Stockmeyer, On the number of nonscalar multiplications necessary to evaluate polynomials, SIAM Journal on Computing, 2 (1973), pp. 60–66.
  * [12] L. N. Trefethen, Is Gauss quadrature better than Clenshaw–Curtis?, SIAM review, 50 (2008), pp. 67–87.
  * [13] S. Ubaru, J. Chen, and Y. Saad, Fast estimation of tr(f(A)) via stochastic Lanczos quadrature, SIAM Journal on Matrix Analysis and Applications, 38 (2017), pp. 1075–1099.
  * [14] S. Ubaru and Y. Saad, Applications of trace estimation techniques, in International Conference on High Performance Computing in Science and Engineering, Springer, 2017, pp. 19–33.