Why does Hutchinson's trace estimator reduce computation complexity?

Question

Given a matrix $A$, we want to compute its trace, in which we can use a trick name Hutchinson's trace estimator \begin{align} tr(A) = tr(A\mathbb{E}[\epsilon \epsilon^T])=\mathbb{E}[tr(A \epsilon \epsilon^T)]=\mathbb{E}[tr(\epsilon^T A \epsilon)]=\mathbb{E}[\epsilon^T A \epsilon], \end{align} where $\epsilon \sim \mathcal{N}(0,1)$ is a standard normal distribution. Then the trace of $A$ can be estimated using Monte Carlo sampling.

Some literature said, by transforming the computation into quadratic form, the computation complexity of the trace calculators can be reduced. I can't understand it. Why can such a stochastic quadratic way reduce the complexity of trace computation?

Answer 1

You are right that for calculating the trace of a matrix this does not reduce cost vs a simple calculation...but this trick is very useful when we need to compute the trace of a function of a matrix, $tr(f(A))$.

For example, consider estimating a log determinant, which can be rewritten as $tr( log(A) )$, the trace of the log of the matrix. Log determinants can be calculated efficiently using matrix decomposition methods (e.g. SVD) as long as the matrix is small, but these methods scale poorly to large matrices. Or consider (the closely related!) calculation where we need the trace of a matrix raised to a power, e.g. $tr(A^k)$. If you are calculating the trace of a function of a large matrix, Hutchison's trace estimator may dramatically reduce your computational cost.

If you are trying to calculate $tr(A^3)$, for example, where A is some square $D x D$ matrix, you could calculate $AAA$ and then take the trace of this -- for a large matrix this will be quite expensive, or you could use Hutchison's trace estimator to replace that series of matrix multiplications with a series of matrix-vector multiplications which will be substantially cheaper, as long as the number $M$ of probe vectors that you use is $M << D$.

Another use case is when forming the matrix $A$ explicitly would be very expensive but taking a matrix vector product would be much cheaper. For example, consider if we have a matrix formed from $XX^T$, where $X$ is some $N x M$ matrix and $N$ is very large, so that $A$ would be too large to fit in memory. In this case, if you need to evaluate $tr(A)$, it might be much cheaper to use matrix-vector products rather than forming $A$ explicitly.

Basically, the trick helps us reduce the cost of trace estimation when we want to estimate the trace of a function of a matrix or when we are working with matrices that we do not want to form explicitly. Estimating log determinants is a common use case -- this paper uses Hutchison's trace estimator in their approach for approximating log determinants for example:

http://proceedings.mlr.press/v37/hana15.pdf

Answer 2

I have never heard of this trick before, and it does look strange ... since the trace is a sum of the $n$ diagonal terms, while the quadratic form needs a double sum over all the $n^2$ terms. So I guess the "trick" must be used only as one idea in combination with others, for some more involved problem than just a trace.

And indeed, in your own linked reference such examples are given, so just start top read those examples with attention. But I do not find the example in that link very well written, so better to go to the reference they give Randomized algorithms for matrices and data, but that paper does not mention the trick! Search gives a lot of relevant papers, and RANDOMIZED ALGORITHMS FOR ESTIMATING THE TRACE OF AN IMPLICIT SYMMETRIC POSITIVE SEMI-DEFINITE MATRIX says, in its abstract that

these algorithms are useful in applications in which there is no explicit representation of $A$ but rather an efficient method to compute $z^T Az$ given $z$.

So this trick is not useful in isolation, but as a building block for randomized algorithms.

Answer 3

It depends a lot on what accuracy you need as well as the form of the matrix. In the application I had, there was a matrix $A=GP$ that was the product of a large (say $N\times M$) sparse matrix $G$ of genotypes and a projection matrix $P$ on to the residual space of a linear regression. This means $A$ wasn't itself sparse, but computing matrix-vector products $Ax$ or $A^TAx$ was still fast.

I needed the trace of $B=A^TA$ and $B^TB$ as part of an approximation to the distribution of a quadratic form, and they didn't need to be all that accurate. The trace of $B$ can be computed directly in $MN$ operations, but $tr(B^TB)$ requires explicitly computing all of $B$, taking $M^2N$ (assuming $M<N$, with both in the thousands). Hutchinson's trace estimator with sample size $k$ could be computed in the time needed for $2k$ matrix multiplications by $A$ (small compared to $MNk$). The error in the estimator was $O_p(k^{-1/2})$, but that was ok (with $k\sim 500$)

The rest of the approximation involved the first few (~100)leading eigenvalues of $B$ and stochastic SVD gives an algorithm whose time is dominated by doing ~100 matrix-vector multiplications by $B$, so it was desirable for the trace part of the whole computation not to be hugely slower than this.