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?


3 Answers 3


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:


  • $\begingroup$ A huge thanks for your reply. It seems the assigned paper is helpful. I will check it out later. Thank you :) $\endgroup$
    – jzin
    Commented May 23, 2022 at 8:39
  • $\begingroup$ The saving of complexity really comes from the Monte-Carlo estimate introduced by "probing vectors" $\epsilon$, hence highly dependent on the dimensionality. And it would be helpful to clarify that $A$ is assumed to be a $D\times D$ matrix in this post. $\endgroup$
    – Henry.L
    Commented Dec 6, 2022 at 6:01
  • $\begingroup$ Clarified that A is D x D, thanks. It is true that the advantage (or lack thereof) is highly dependent on dimensionality. If we want to calculate a log determinant, for example, for a small matrix we are much better off with matrix decomposition-based methods. Since these scale as O(D^3), however, for a large matrix stochastic Lanczos quadrature or other methods involving the Hutchison's trace estimator become highly attractive. $\endgroup$
    – HappyDog
    Commented Dec 7, 2022 at 4:32
  • $\begingroup$ Calculating the exact trace of $A^2$ actually only has a time complexity of $O(D^2)$, since it doesn’t require calculating the entire $A^2$ matrix, but just the diagonal elements, which has a $O(D^2)$ time complexity. But I guess that for calculating the trace of $A^3$, using Hutchinson's trace estimator may make sense. $\endgroup$ Commented Oct 30, 2023 at 6:08
  • 1
    $\begingroup$ This is true! I have edited the answer to reflect this. $\endgroup$
    – HappyDog
    Commented Oct 31, 2023 at 13:13

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.

  • 1
    $\begingroup$ Hi, thanks for your quick reply. Yes, maybe you are right, this trick is not useful in isolation. Actually, this trick is quite commonly used in machine learning papers (like, see section 3.1 in this paper openreview.net/pdf?id=rJxgknCcK7 ), say, it is generally used in the computation of the trace of the Jacobian matrix. But, I still can't understand why should we use this trick to do that and how can it reduce the computation complexity. I mean I can compute the Jacobian matrix regularly and sum up the diagonal items to obtain its trace. $\endgroup$
    – jzin
    Commented Sep 11, 2021 at 12:56
  • $\begingroup$ I review this question after a few months by chance. I think you're right: "this trick is not useful in isolation." I review the paper discussed above openreview.net/pdf?id=rJxgknCcK7 I found that this trick (Hutchinson's trace estimator) is applied together with the other trick (vector-Jacobian product). These two tricks work together to make sense. It dramatically reduces the complexity of the trace estimation of the Jacobian matrix. $\endgroup$
    – jzin
    Commented May 23, 2022 at 8:45

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.