Imagine a matrix behind a curtain.

You are not allowed to read its diagonal. You are not allowed to form it. You are only allowed to hand it a vector and receive the product back.

This sounds like a bad way to compute a trace. The trace is the diagonal sum:

tr(A) = A_11 + A_22 + ... + A_nn

If the diagonal is hidden, what could a random vector possibly know?

The answer is one of those small pieces of randomized numerical linear algebra that feels like a magic trick until you write the expectation:

E[z z^T] = I

Then:

\[\mathbb{E}[z^\top A z] = \mathbb{E}[\mathrm{tr}(z^\top A z)] = \mathbb{E}[\mathrm{tr}(Azz^\top)] = \mathrm{tr}(A).\]

The trace can be heard through noise.

This is the Hutchinson trace estimator. Draw random probe vectors, compute the quadratic forms (z^\top A z), and average them. Hutchinson’s original 1989 paper used Rademacher probes, entries independently equal to +1 or -1, to estimate traces of influence matrices in smoothing-spline calculations.1 The same primitive now appears in large-scale statistics, scientific computing, machine learning, graph algorithms, and continuous normalizing flows.

It is not powerful because it is fancy. It is powerful because a lot of modern matrices are easier to multiply by than to inspect.

Curtain Is Normal

In a small textbook problem, the matrix is an array.

In real computation, it is often a procedure:

Hessian-vector product
Jacobian-vector product
linear solve
kernel operator
graph diffusion
implicit covariance
automatic-differentiation closure

The matrix may have (10^8) entries, but a product (Av) may be cheap because it is computed by backpropagation, a sparse operator, an iterative method, or a structured kernel. If all you need is a global scalar such as trace, effective dimension, divergence, log-determinant, or a spectral moment, materializing the whole object can be the wrong shape of work.

The trace estimator changes the access pattern:

do not ask for entries
ask for actions

That is why the formula keeps reappearing.

For a random vector with zero mean and identity covariance:

\[\mathbb{E}[z_i z_j] = \begin{cases} 1, & i=j,\\ 0, & i\ne j. \end{cases}\]

Expand the quadratic form:

\[z^\top A z = \sum_i A_{ii} z_i^2 + \sum_{i\ne j} A_{ij}z_iz_j.\]

The diagonal terms survive in expectation. The cross terms cancel.

That is the whole unbiasedness proof.

The Catch Is Variance

Unbiased is not the same as useful.

If one probe gives:

estimate = z^T A z

then the Monte Carlo average over m probes converges at the usual (1/\sqrt{m}) rate. The key question is how noisy each quadratic form is.

For a symmetric matrix and Gaussian probes, the single-probe variance is:

\[2\|A\|_F^2.\]

For Rademacher probes, it is:

\[2\left(\|A\|_F^2-\sum_i A_{ii}^2\right).\]

That second formula explains a delightful edge case: if (A) is diagonal, Rademacher Hutchinson is exact in one probe, because (z_i^2=1). Gaussian probes still fluctuate because (z_i^2) fluctuates.

So the question is not merely:

Is the estimator unbiased?

The useful question is:

Where does the variance live?

Avron and Toledo’s 2011 analysis made this practical concern explicit by studying sample counts needed for relative-error guarantees, rather than only comparing variance formulas.2 Bounds matter because production users do not ask, “Is the estimator elegant?” They ask, “How many matvecs do I need?”

A Trace Meter

The lab below builds small synthetic positive semidefinite matrices and then pretends they are implicit. The code knows the exact trace for auditing, but the estimators only use matrix-vector products.

The blue estimator is classic Hutchinson with the selected probe distribution. The green estimator is Hutch++, which spends part of the same matvec budget to learn a small dominant subspace before probing the residual.3

The default matrix is spiked: five directions carry about 76% of the trace, and the effective rank is about 5.4. With the default seed and 36 matvecs, the exact trace is 63.74. One Hutchinson run lands at 77.27, while Hutch++ lands at 64.88. Over repeated trials, the relative RMSE is about 9.5% for classic Hutchinson and 1.4% for Hutch++.

That is the modern story in miniature: first remove the loud room, then measure the quieter remainder.

I audited the lab before rewriting. The default matrix has trace 63.741, top-five trace mass 76.22%, effective rank 5.44, and Hutch++ uses all 36 allowed matvecs. The built-in checks confirm symmetry, finite positive trace, trace-diagonal consistency, and budget respect. A 432-case sweep over spectrum shape, probe distribution, matvec budget, dimension, trial count, and seed had no check failures. Hutch++ beat classic Hutchinson in 123 cases and lost or tied in 309, which is exactly the caveat the article needs: deflation helps when there is a loud subspace. In all 54 diagonal/Rademacher cases, classic Hutchinson was exact.

Hutchinson Hutch++ Exact trace Spectrum tail Matrix entries

Deterministic toy model. Matrices are explicit only for auditing and visualization; the estimators use matrix-vector products. Hutch++ uses Gaussian sketch and residual probes, while the classic estimator uses the selected probe distribution.

Try four moves.

First, keep spiked low rank and lower the budget. Hutch++ often remains stable because the random range sketch captures much of the large eigenspace. The residual trace is smaller and its Frobenius norm is much smaller.

Second, switch to flat spectrum. The advantage shrinks. There is no small room to discover, so using part of the budget to build a low-rank sketch may not pay for itself.

Third, switch to diagonal smoother with Rademacher probes. Classic Hutchinson becomes exact, because all off-diagonal variance is gone and Rademacher squares are identically one. Now switch to Gaussian probes. The same diagonal trace starts wobbling again.

Fourth, use slow power tail and increase the budget. The estimator improves, but slowly. This is the quiet lesson: Monte Carlo does not care how badly you want more digits. If the variance is there, it has to be paid down.

Hutch++ Is a Variance Budget

Hutch++ starts from the observation that not all variance deserves the same treatment.

If a positive semidefinite matrix has a few huge eigenvalues, a plain random quadratic form can be dominated by whether the probe happened to point into those directions. Averaging fixes this eventually, but eventually can be a lot of matvecs.

Hutch++ spends roughly one third of the budget on a sketch:

Y = A S
Q = orth(Y)

Then it computes the trace captured by that subspace:

tr(Q^T A Q)

and uses Hutchinson only on the residual:

(I - QQ^T) A (I - QQ^T)

That is the low-rank-plus-residual pattern. Do deterministic work where the matrix is loud, use randomness where the matrix is quieter.

Meyer, Musco, Musco, and Woodruff show that Hutch++ can obtain a ((1\pm\epsilon)) approximation for positive semidefinite matrices using (O(1/\epsilon)) matrix-vector products, improving the (O(1/\epsilon^2)) dependence of classic Hutchinson in that setting.3

The lab is not a proof of that theorem. It is the visual shape of the theorem: when spectral mass is concentrated, deflation changes the game.

Why Trace Estimation Shows Up In ML

The trace is not always the final quantity. Often it is the inside of a larger quantity.

For a positive definite matrix (K), a log-determinant can be written as:

\[\log\det(K)=\mathrm{tr}(\log K).\]

That turns a determinant problem into a trace-of-a-matrix-function problem. You still need to approximate (\log(K)v), but now matrix-vector products and Krylov methods enter the room. Ubaru, Chen, and Saad combine stochastic trace estimation with Lanczos quadrature for estimating (\mathrm{tr}(f(A))), including log-determinants, Schatten norms, trace inverses, and related quantities in machine learning and scientific computing.4

Continuous normalizing flows give a different example. The instantaneous change-of-variables formula contains the trace of a Jacobian:

\[\frac{d}{dt}\log p(z(t)) = -\mathrm{tr}\left(\frac{\partial f}{\partial z(t)}\right).\]

Computing the exact Jacobian trace can be expensive in high dimension. FFJORD uses Hutchinson’s estimator to get a scalable unbiased estimate of the log-density while allowing less restricted neural-network dynamics.5 In automatic differentiation terms, the estimator turns a trace into a vector-Jacobian product against a random vector.

That same pattern appears when people estimate Hessian traces during training, effective dimensions of kernels, uncertainty terms, graph centrality measures, or physical observables tied to traces of inverse operators.

The matrix is not always there.

The action of the matrix is.

What the Estimator Does Not Promise

The simple identity:

E[z^T A z] = tr(A)

is easy to oversell.

It does not say one sample is accurate. It does not say relative error is stable when the true trace is near zero. It does not say Hutch++ is always better than plain Hutchinson. It does not say the probe distribution is irrelevant. It does not say your matrix-vector product is exact, cheap, or deterministic.

The diagonal case in the lab is deliberately embarrassing for fancy methods. Rademacher Hutchinson wins immediately because the matrix has the one structure that makes its variance vanish. A flat spectrum can also make Hutch++ look ordinary because the sketch has no dominant room to capture.

This is a good thing. A numerical method that cannot be embarrassed by a simple counterexample is probably hiding its assumptions.

Here are the practical questions I would ask before using this estimator in a real system:

Is the target matrix symmetric or positive semidefinite?
Is the trace far from zero?
Do I need additive error or relative error?
Can I reuse probes across time, or does that bias an adaptive loop?
How noisy or approximate are the matrix-vector products?
Does the spectrum have spikes that a low-rank sketch can remove?

The estimator is a measuring device, not a spell.

A Small Philosophy of Matrix-Free Computation

There is a recurring pattern in large-scale computation:

the object is too large to inspect
but not too large to interrogate

Randomized SVD asks the matrix where its range is. Johnson-Lindenstrauss projections ask a point cloud which distances matter. Stochastic Lanczos asks an operator for spectral moments. Hutchinson asks for the diagonal sum without looking at the diagonal.

This is not randomness as decoration. It is experiment design.

A random probe is a question with a known expectation. If you know the expectation and can control the variance, you can replace a forbidden measurement with a sequence of cheap, auditable measurements.

That is the real beauty here.

The trace was never invisible. It was just hidden in quadratic forms.

  1. M. F. Hutchinson, “A Stochastic Estimator of the Trace of the Influence Matrix for Laplacian Smoothing Splines”, Communications in Statistics - Simulation and Computation 18(3), 1989. ↩

  2. Haim Avron and Sivan Toledo, “Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix”, Journal of the ACM 58(2), 2011. ↩

  3. Raphael A. Meyer, Cameron Musco, Christopher Musco, and David P. Woodruff, “Hutch++: Optimal Stochastic Trace Estimation”, SOSA 2021. ↩ ↩2

  4. Shashanka Ubaru, Jie Chen, and Yousef Saad, “Fast Estimation of tr(f(A)) via Stochastic Lanczos Quadrature”, SIAM Journal on Matrix Analysis and Applications 38(4), 2017. ↩

  5. Will Grathwohl, Ricky T. Q. Chen, Jesse Bettencourt, Ilya Sutskever, and David Duvenaud, “FFJORD: Free-Form Continuous Dynamics for Scalable Reversible Generative Models”, ICLR 2019. ↩