Hutch++ and the Split Spectrum
10 min read
The trace of a matrix is the sum of its diagonal entries. That sounds cheap to calculate, but sometimes we can multiply by a matrix much more cheaply than we can read its entries. We want to estimate the trace using only a few of those products.
Multiplying by random vectors gives a noisy estimate. Hutch++ makes a beautiful use of that budget: it first finds a part of the trace we can compute exactly [MMMW, SOSA ’21]. Only the remainder needs sampling.
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
Let \(\mathbf{A} \in \mathbb{R}^{d \times d}\) be a symmetric matrix that we can access by choosing a vector \(\mathbf{x} \in \mathbb{R}^d\) and receiving \(\mathbf{A}\mathbf{x}\). Often \(\mathbf{A}\) is \(\mathbf{B}^{-1}\), \(\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}\)). 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 [Hutchinson, ’89]. 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\) preserve 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 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}.\]
For error relative to the trace, 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. That last step is tight when \(\|\mathbf{A}\|_F \approx \mathrm{tr}(\mathbf{A})\), which happens when a few large eigenvalues carry most of the trace. On a flat spectrum, all \(d\) eigenvalues equal, the ratio is \(1/\sqrt{d}\) instead and Hutchinson is nearly free. Those large eigenvalues also make their eigenvectors easier to find: multiplying a random vector by \(\mathbf{A}\) amplifies its components in those directions. Hutch++ uses some of its products to learn that subspace and the rest to estimate the trace left outside it.
The two methods can be compared at the same matrix-vector query budget 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.
The Split
Hutch++ approximates the top eigenspace using matrix-vector products alone and splits its budget of \(m\) into 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.
Why the two traces add
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.
The Tail Inequality
First consider an ideal sketch that captures the top \(k\) eigenvectors exactly, where \(1 \le k < d\). The Hutchinson bound then depends on the squared eigenvalues left in the tail. 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, retaining the \(k\) largest eigenvalues.
Lemma: for any PSD \(\mathbf{A}\) and any \(k\),
\[\|\mathbf{A} - \mathbf{A}_k\|_F \le \frac{\mathrm{tr}(\mathbf{A})}{\sqrt{k}}.\]
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}}.\]
The right side does not mention the spectrum at all: one heavy eigenvalue or \(d\) equal ones, the same ceiling holds. For \(2k\) equal eigenvalues, the bound is tight up to a factor of two.
The two methods handle complementary spectra: low-rank sketches capture large leading eigenvalues, while Hutchinson is accurate on a flat remainder.
In the plot, the ideal rank-\(k\) sketch handles the teal head exactly while Hutchinson only sees the gray tail, every bar of which sits under the ceiling \(\mathrm{tr}(\mathbf{A})/k\).
The Accuracy
The tail bound lets us compare the error of the split estimator with the original random estimate. First suppose the sketch is perfect: \(\mathbf{Q}\) spans the top \(k = m/3\) eigenvectors, so \(\mathbf{M} = \mathbf{A} - \mathbf{A}_k\) exactly.
Claim: with a perfect sketch of rank \(k = m/3\),
\[\mathrm{Var}\big(\widehat{\mathrm{tr}}^{\,{+}{+}}(\mathbf{A})\big) \le \frac{18\, \mathrm{tr}(\mathbf{A})^2}{m^2}.\]
Proof
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)\).
In the actual algorithm, one power step on \(m/3\) random vectors only approximates the top eigenspace. With \(m/3\) random probes, the standard analysis of the randomized range finder [HMT, SIAM Rev. ’11] loses only a constant: 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}\]
Orthogonal projections do not increase the Frobenius norm, 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.
Deflation as a Control Variate
The same construction is a matrix control variate: it computes the surrogate trace exactly and estimates only the remainder. The scalar version, where the subtracted object is a single tunable constant, has a page of its own. Its learned rank-\(m/3\) surrogate contributes \(\mathrm{tr}(\mathbf{Q}^\top\mathbf{A}\mathbf{Q})\) exactly; Hutchinson handles the trace left over. The surrogate is learned from the matrix itself with the first third of the budget rather than supplied in advance as side knowledge. The range-finder bound makes the remainder’s Frobenius norm \(O(\mathrm{tr}(\mathbf{A})/\sqrt{m})\), so Hutchinson estimates its trace with error \(O(\mathrm{tr}(\mathbf{A})/m)\).
The Matching Lower Bound
The \(1/m\) rate is nearly optimal: [MMMW, SOSA ’21] rule out a faster rate 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 method applies wherever the target can be written as 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.