Hutch++ and the Split Spectrum
11 min read
Estimating the trace of a matrix you can only touch through matrix-vector products is one of the workhorse problems of numerical linear algebra. The standard tool is a one-line randomized estimator that Michael F. Hutchinson [Hutchinson, ’89] introduced to count the effective degrees of freedom of a smoothing spline. Its \(1/\sqrt{m}\) convergence rate then stood unimproved for three decades. Meyer, Musco, Musco & Woodruff [MMMW, SOSA ’21] improved it to \(1/m\) with a short algorithm, and proved that, up to a logarithmic factor, no matrix-vector algorithm can do better. Hutch++ is a hometown result for me: half of its author list was down the hall during my PhD at NYU.
One Random Quadratic Form
The matrix whose trace we want is usually one we never hold: \(\mathbf{A} \in \mathbb{R}^{d \times d}\) is \(\mathbf{B}^{-1}\), or \(\exp(\mathbf{B})\), or \(\log(\mathbf{B})\), some transformation of a matrix \(\mathbf{B}\) we actually store. Forming \(\mathbf{A}\) entrywise costs a full \(O(d^3)\) eigendecomposition, while multiplying a vector by it costs about as much as multiplying by \(\mathbf{B}\) (via conjugate gradient or Lanczos, for reasonably conditioned \(\mathbf{B}\)). This access model is implicit trace estimation: a symmetric matrix \(\mathbf{A}\) that we touch only by choosing a vector \(\mathbf{x} \in \mathbb{R}^d\) and receiving \(\mathbf{A}\mathbf{x}\). Probing with the \(d\) standard basis vectors reads off the diagonal and hence the trace exactly, so the game is to spend \(m \ll d\) products.
Hutchinson’s estimator uses each product on a random quadratic form. To see what one product buys, take \(d = 2\) and probe with two independent fair signs \(z_1, z_2 \in \{-1, +1\}\):
\[\begin{pmatrix} z_1 & z_2 \end{pmatrix}\begin{pmatrix} [\mathbf{A}]_{1,1} & [\mathbf{A}]_{1,2} \\ [\mathbf{A}]_{1,2} & [\mathbf{A}]_{2,2} \end{pmatrix}\begin{pmatrix} z_1 \\ z_2 \end{pmatrix} = [\mathbf{A}]_{1,1} + [\mathbf{A}]_{2,2} + 2\,[\mathbf{A}]_{1,2}\, z_1 z_2.\]
The squares \(z_1^2 = z_2^2 = 1\) hand back both diagonal entries exactly, on every draw, and the only randomness left is the off-diagonal entry times the fair sign \(z_1 z_2\). The trace arrives whole; the noise is one coin flip. The same split happens in any dimension. Draw \(\mathbf{z} \in \{-1, +1\}^d\) with independent fair signs and expand entrywise:
\[\mathbf{z}^\top \mathbf{A}\mathbf{z} = \sum_{i=1}^d [\mathbf{A}]_{i,i}\, z_i^2 + \sum_{i \neq j} [\mathbf{A}]_{i,j}\, z_i z_j.\]
Each sign satisfies \(z_i^2 = 1\), so the diagonal sum is the trace on every draw, not merely in expectation:
\[\sum_{i=1}^d [\mathbf{A}]_{i,i}\, z_i^2 = \sum_{i=1}^d [\mathbf{A}]_{i,i} = \mathrm{tr}(\mathbf{A}).\]
Each cross term has mean zero, since \(z_i\) and \(z_j\) are independent fair signs: \(\mathbb{E}[z_i z_j] = \mathbb{E}[z_i]\,\mathbb{E}[z_j] = 0\) for \(i \neq j\). So the quadratic form is unbiased for the trace: \(\mathbb{E}\big[\mathbf{z}^\top \mathbf{A}\mathbf{z}\big] = \mathrm{tr}(\mathbf{A})\).
Averaging the quadratic form over \(m\) independent probes \(\mathbf{z}_1, \dots, \mathbf{z}_m\), at one matrix-vector product each, is Hutchinson’s estimator:
\[\widehat{\mathrm{tr}}^{\,\mathrm{hutch}}(\mathbf{A}) = \frac{1}{m}\sum_{j=1}^{m} \mathbf{z}_j^\top \mathbf{A}\mathbf{z}_j.\]
All of its error is fluctuation, and the diagonal contributes none of it: the variance is entirely off-diagonal mass. An unordered pair \(\{i, j\}\) with \(i \neq j\) contributes \(2[\mathbf{A}]_{i,j}\, z_i z_j\) to the quadratic form, since \(\mathbf{A}\) is symmetric. That term has variance \(4[\mathbf{A}]_{i,j}^2\), because \(z_i z_j\) is itself a fair sign with mean zero and square one. The squared Frobenius norm is the sum of all squared entries: \(\|\mathbf{A}\|_F^2 = \sum_{i,j} [\mathbf{A}]_{i,j}^2\). Distinct pairs are uncorrelated, since any two differ in some index whose sign flips independently of everything else, so the variances add:
\[\begin{align} \mathrm{Var}\big(\mathbf{z}^\top\mathbf{A}\mathbf{z}\big) &= \sum_{\{i,j\}:\, i \neq j} 4[\mathbf{A}]_{i,j}^2 \\ &= 2 \sum_{i \neq j} [\mathbf{A}]_{i,j}^2 \\ &\le 2 \sum_{i,j} [\mathbf{A}]_{i,j}^2 = 2\, \|\mathbf{A}\|_F^2, \end{align}\]
where the second line turns the sum over unordered pairs into a sum over ordered pairs (each squared entry counted twice), and the third adds the nonnegative diagonal squares to complete the Frobenius norm. Averaging \(m\) independent probes divides this variance by \(m\):
\[\mathrm{Var}\big(\widehat{\mathrm{tr}}^{\,\mathrm{hutch}}(\mathbf{A})\big) \le \frac{2\,\|\mathbf{A}\|_F^2}{m}.\]
What a user cares about is error relative to the trace, so we put the two on one scale. Now let \(\mathbf{A}\) be positive semidefinite (PSD), the common case in applications, with eigenvalues \(\lambda_1 \ge \dots \ge \lambda_d \ge 0\). For a PSD matrix the trace is the sum of the eigenvalues, and the squared Frobenius norm is the sum of their squares:
\[\mathrm{tr}(\mathbf{A}) = \sum_{i=1}^d \lambda_i, \qquad \|\mathbf{A}\|_F^2 = \mathrm{tr}(\mathbf{A}^2) = \sum_{i=1}^d \lambda_i^2.\]
Taking the square root, the estimator’s standard deviation is at most \(\sqrt{2/m}\, \|\mathbf{A}\|_F\), so its size relative to the trace is:
\[ \frac{\sqrt{2/m}\, \|\mathbf{A}\|_F}{\mathrm{tr}(\mathbf{A})} = \sqrt{\frac{2}{m}} \cdot \frac{\sqrt{\sum_i \lambda_i^2}}{\sum_i \lambda_i} \le \sqrt{\frac{2}{m}}, \]
where the last step is because the Euclidean length of a nonnegative vector never exceeds the sum of its entries. Relative error \(\varepsilon\) therefore costs \(m = O(1/\varepsilon^2)\) products: the \(1/\sqrt{m}\) rate. But look at when that last step is tight: \(\|\mathbf{A}\|_F \approx \mathrm{tr}(\mathbf{A})\) exactly when a few large eigenvalues carry the whole trace. On a flat spectrum, all \(d\) eigenvalues equal, the ratio is \(1/\sqrt{d}\) instead and Hutchinson is nearly free. A heavy top is Hutchinson’s worst case.
The Tail Inequality
A heavy top is also exactly what a low-rank approximation captures, so suppose we could strip off the top \(k\) eigendirections and hand Hutchinson the remainder. How small is that remainder guaranteed to be? Write the eigendecomposition both ways, factored and as a sum of outer products:
\[\mathbf{A} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^\top = \sum_{i=1}^d \lambda_i \mathbf{v}_i \mathbf{v}_i^\top,\]
with \(\mathbf{V} \in \mathbb{R}^{d \times d}\) the orthonormal eigenvectors in its columns and \(\mathbf{\Lambda}\) the diagonal matrix of eigenvalues. Let \(\mathbf{A}_k = \sum_{i \le k} \lambda_i \mathbf{v}_i \mathbf{v}_i^\top\) be the best rank-\(k\) approximation, the heavy top and nothing else; the outer-product form is the one we manipulate.
Lemma: for any PSD \(\mathbf{A}\) and any \(k\),
\[\|\mathbf{A} - \mathbf{A}_k\|_F \le \frac{\mathrm{tr}(\mathbf{A})}{\sqrt{k}}.\]
The right side does not mention the spectrum at all: one heavy eigenvalue or \(d\) equal ones, the same ceiling holds.
Proof
The proof rests on two facts about the tail eigenvalues \(\lambda_{k+1}, \dots, \lambda_d\). First, \(\lambda_{k+1}\) is at most the average of the \(k\) eigenvalues above it, each of which is at least as large, and that average is at most a \(k\)-th of the whole trace:
\[\lambda_{k+1} \le \frac{1}{k}\sum_{i=1}^k \lambda_i \le \frac{1}{k}\sum_{i=1}^d \lambda_i = \frac{\mathrm{tr}(\mathbf{A})}{k}.\]
Second, in the outer-product form the tail is \(\mathbf{A} - \mathbf{A}_k = \sum_{i > k} \lambda_i \mathbf{v}_i \mathbf{v}_i^\top\), and because the \(\mathbf{v}_i\) are orthonormal its squared Frobenius norm is the sum of the squared tail eigenvalues. Bounding one factor of \(\lambda_i\) by \(\lambda_{k+1}\) and then the remaining tail sum by the full trace:
\[\begin{align} \|\mathbf{A} - \mathbf{A}_k\|_F^2 &= \sum_{i > k} \lambda_i^2 \\ &\le \lambda_{k+1} \sum_{i > k} \lambda_i \\ &\le \frac{\mathrm{tr}(\mathbf{A})}{k} \sum_{i > k} \lambda_i \\ &\le \frac{\mathrm{tr}(\mathbf{A})}{k} \cdot \mathrm{tr}(\mathbf{A}). \end{align}\]
Taking square roots gives the lemma: \(\|\mathbf{A} - \mathbf{A}_k\|_F \le \frac{\mathrm{tr}(\mathbf{A})}{\sqrt{k}}.\)
Which spectrum makes this essentially tight? \(2k\) equal eigenvalues, where the bound is off by exactly a factor of two.
Read the lemma as a statement about failure modes: heavy spectra are easy for a low-rank sketch, flat spectra are easy for Hutchinson, and every spectrum splits into one of each.
In the plot, a rank-\(k\) sketch handles the teal head exactly while Hutchinson only ever sees the gray tail, every bar of which sits under the ceiling \(\mathrm{tr}(\mathbf{A})/k\). Exploiting that complementarity takes the top \(k\) eigenvectors, which we cannot afford to compute. All we can buy is matrix-vector products.
The Split
Hutch++ buys a crude substitute for the top eigenvectors out of products alone, spending its budget of \(m\) in thirds. Draw two independent random sign matrices \(\mathbf{S}, \mathbf{G} \in \{-1, +1\}^{d \times m/3}\), the same \(\pm 1\) probes the variance bound above already covers. Spend the first third computing \(\mathbf{A}\mathbf{S}\), and let \(\mathbf{Q} \in \mathbb{R}^{d \times m/3}\) be an orthonormal basis for its column span, from a QR decomposition. This is one step of the power method from a random start: multiplying by \(\mathbf{A}\) stretches each random column along the large eigenvalues, so \(\mathbf{Q}\) comes out coarsely aligned with the top eigenvectors. Then split the trace along \(\mathbf{Q}\) and its orthogonal complement:
\[\mathrm{tr}(\mathbf{A}) = \mathrm{tr}\big(\mathbf{Q}^\top \mathbf{A} \mathbf{Q}\big) + \mathrm{tr}\big((\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\, \mathbf{A}\, (\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\big).\]
In words, the first term is the trace of \(\mathbf{A}\) inside the sketched subspace, and the second is the trace of whatever the sketch missed. The identity needs only that the trace is invariant under cyclic permutation of a product. Cycling \(\mathbf{Q}^\top\) to the back turns the first term into a projection of \(\mathbf{A}\): \(\mathrm{tr}\big(\mathbf{Q}^\top \mathbf{A}\mathbf{Q}\big) = \mathrm{tr}\big(\mathbf{A}\mathbf{Q}\mathbf{Q}^\top\big)\). The second term cycles the same way, and \(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top\) is a projection, equal to its own square (since \(\mathbf{Q}^\top\mathbf{Q} = \mathbf{I}\)), so one of the two factors collapses:
\[\begin{align} \mathrm{tr}\big((\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\mathbf{A}(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\big) &= \mathrm{tr}\big(\mathbf{A}(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)^2\big) \\ &= \mathrm{tr}\big(\mathbf{A}(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\big). \end{align}\]
Adding the two and using linearity, the projections cancel and leave the trace:
\[\begin{align} \mathrm{tr}\big(\mathbf{A}\mathbf{Q}\mathbf{Q}^\top\big) + \mathrm{tr}\big(\mathbf{A}(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\big) &= \mathrm{tr}\big(\mathbf{A}\mathbf{Q}\mathbf{Q}^\top + \mathbf{A} - \mathbf{A}\mathbf{Q}\mathbf{Q}^\top\big) \\ &= \mathrm{tr}(\mathbf{A}). \end{align}\]
Hutch++ computes the first term exactly, spending the second third of the budget on \(\mathbf{A}\mathbf{Q}\), and estimates the second term with Hutchinson, spending the last third on the probes in \(\mathbf{G}\). Equal thirds is a balance: the variance below falls like one over the sketch rank times the probe count, so starving either piece hurts, and any other constant split changes only the constant out front. (The QR and the small traces cost \(O(dm^2)\) extra arithmetic but no extra matrix-vector products, which are what we count.)
Write \(\mathbf{M}\) for the deflated matrix, whatever the sketch failed to capture: \(\mathbf{M} = (\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\,\mathbf{A}\,(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\). The whole algorithm is one exact trace plus one Hutchinson average over the columns \(\mathbf{g}_j \in \{-1, +1\}^d\) of \(\mathbf{G}\):
\[\widehat{\mathrm{tr}}^{\,{+}{+}}(\mathbf{A}) = \mathrm{tr}\big(\mathbf{Q}^\top \mathbf{A}\mathbf{Q}\big) + \frac{3}{m}\sum_{j=1}^{m/3} \mathbf{g}_j^\top \mathbf{M} \mathbf{g}_j.\]
The second term is the first section’s estimator run on \(\mathbf{M}\) in place of \(\mathbf{A}\), so its variance bound carries over verbatim: nothing in that derivation used more than symmetry. Pretend for a moment the sketch were perfect, with \(\mathbf{Q}\) spanning the top \(k = m/3\) eigenvectors, so that \(\mathbf{M} = \mathbf{A} - \mathbf{A}_k\) exactly. Then the first term is a fixed number, so all the randomness in Hutch++ lives in the second, and the lemma caps the variance of every probe:
\[\begin{align} \mathrm{Var}\big(\mathbf{g}_j^\top \mathbf{M} \mathbf{g}_j\big) &\le 2\, \|\mathbf{M}\|_F^2 \\ &= 2\, \|\mathbf{A} - \mathbf{A}_k\|_F^2 \\ &\le \frac{2\, \mathrm{tr}(\mathbf{A})^2}{k}. \end{align}\]
Averaging the \(m/3\) probes divides by \(m/3\), and \(k = m/3\) collapses the two factors:
\[\begin{align} \mathrm{Var}\big(\widehat{\mathrm{tr}}^{\,{+}{+}}(\mathbf{A})\big) &\le \frac{2\, \mathrm{tr}(\mathbf{A})^2}{k \cdot (m/3)} \\ &= \frac{2\, \mathrm{tr}(\mathbf{A})^2}{(m/3)^2} \\ &= \frac{18\, \mathrm{tr}(\mathbf{A})^2}{m^2}. \end{align}\]
Taking the square root, the standard deviation is \(O(\mathrm{tr}(\mathbf{A})/m)\): \(\sqrt{\mathrm{Var}\big(\widehat{\mathrm{tr}}^{\,{+}{+}}(\mathbf{A})\big)} \le \sqrt{18}\, \mathrm{tr}(\mathbf{A})/m\). So relative error \(\varepsilon\) costs \(m = O(1/\varepsilon)\) products, a quadratic improvement over Hutchinson’s \(O(1/\varepsilon^2)\).
Of course \(\mathbf{Q}\) is not perfect: one power step on \(m/3\) random vectors is not an eigendecomposition. The imperfection costs only constants. With \(m/3\) random probes, the standard analysis of the randomized range finder gives the following bound, with high probability, for a rank \(k\) proportional to \(m\): \(\|\mathbf{A} - \mathbf{A}\mathbf{Q}\mathbf{Q}^\top\|_F^2 \le 2\, \|\mathbf{A} - \mathbf{A}_k\|_F^2\). The deflated matrix is this same quantity, projected once more. Since \(\mathbf{A} - \mathbf{A}\mathbf{Q}\mathbf{Q}^\top = \mathbf{A}(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\), left-multiplying by the projection \(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top\) only shrinks the Frobenius norm:
\[\begin{align} \|\mathbf{M}\|_F &= \|(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\, \mathbf{A}(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\|_F \\ &\le \|\mathbf{A}(\mathbf{I} - \mathbf{Q}\mathbf{Q}^\top)\|_F \\ &= \|\mathbf{A} - \mathbf{A}\mathbf{Q}\mathbf{Q}^\top\|_F. \end{align}\]
(Do you see why a projection can never lengthen anything?) So the deflated matrix is at most twice as heavy as the ideal \(\mathbf{A} - \mathbf{A}_k\), and the variance bound loses only a constant. Running this argument with tail bounds, [MMMW, SOSA ’21] prove that for PSD \(\mathbf{A}\), \(m = O(1/\varepsilon)\) products return a \((1 \pm \varepsilon)\)-approximation of \(\mathrm{tr}(\mathbf{A})\) with any fixed success probability.
Both rates are visible in an experiment on a \(2000 \times 2000\) PSD matrix with eigenvalues \(\lambda_i = 1/i\), averaged over 100 trials.
In the plot, the fitted slopes come out to \(-0.48\) for Hutchinson and \(-1.05\) for Hutch++, each tracking its guide line. And the shape of the estimator, one exact piece plus a cheap estimate of what is left, is not special to traces.
Deflation as a Control Variate
That move has a name: a control variate is a quantity you can evaluate exactly, subtracted from the target so that only the remainder has to be estimated. The scalar version, where the subtracted object is a single tunable constant, has a page of its own. Rearranged, Hutch++ is exactly this: subtract the trace of the learned rank-\(m/3\) surrogate, computed exactly as \(\mathrm{tr}(\mathbf{Q}^\top\mathbf{A}\mathbf{Q})\), and let Hutchinson estimate only the trace that is left over. The unusual part is where it came from: learned from the matrix itself with the first third of the budget, rather than supplied in advance as side knowledge. And the lemma guarantees the deflation always helps: whatever the spectrum, the surrogate absorbs enough Frobenius norm that the leftover fluctuates like \(\mathrm{tr}(\mathbf{A})/\sqrt{k}\) rather than like \(\mathrm{tr}(\mathbf{A})\). Learning the control variate is the whole difference between \(1/\sqrt{m}\) and \(1/m\).
The Matching Lower Bound
Could a cleverer split, or something else entirely, push past \(1/m\)? No: [MMMW, SOSA ’21] rule it out by reducing from the Gap-Hamming problem of communication complexity. Two parties holding sign vectors \(\mathbf{s}, \mathbf{t} \in \{-1, +1\}^d\) can build a PSD matrix whose trace encodes their inner product: \(\mathrm{tr}(\mathbf{A}) = 2d + 2\langle \mathbf{s}, \mathbf{t} \rangle\). A sharp enough trace estimate would then decide whether \(\langle \mathbf{s}, \mathbf{t} \rangle\) is large or small (a decision known to cost \(\Omega(d)\) bits of communication), while the parties can simulate each matrix-vector query with only about \(\sqrt{d}\) bits. An algorithm with too few queries would be a communication protocol that is too good. So any algorithm issuing bounded-precision queries, however adaptively it chooses them, needs \(\Omega\big(1/(\varepsilon \log(1/\varepsilon))\big)\) of them for a \((1 \pm \varepsilon)\)-approximation. A separate information-theoretic argument gives a clean \(\Omega(1/\varepsilon)\) for non-adaptive algorithms with real-valued queries. Up to the logarithm, no matrix-vector algorithm can beat the split.
The recipe gets used wherever the number we actually want turns out to be a trace. Log-determinants are the standard case: \(\log\det(\mathbf{B}) = \mathrm{tr}(\log \mathbf{B})\) is the expensive term in a Gaussian log-likelihood, and \(\log\mathbf{B}\) is a matrix we can multiply by but never want to form. Spectral sums \(\sum_i f(\lambda_i) = \mathrm{tr}(f(\mathbf{A}))\) are the general version, counting eigenvalues in a window when \(f\) is a smoothed indicator of it. Diagonal estimation, where we want every diagonal entry rather than their sum, deflates and probes entry by entry. And when an estimator has one known worst case, it is worth spending part of the budget deleting that case exactly and estimating only what is left.