Hutch++ and the Split Spectrum
11 min read
We want the trace of a large matrix, the sum of its diagonal entries, but can only access the matrix through products with vectors. For a covariance matrix, this trace is the total variance across its coordinates. Hutchinson’s estimator recovers it approximately from a small number of random matrix-vector products [Hutchinson, ’89].
Much of the estimation error can come from a few directions where the matrix is largest. Hutch++ uses some products to find a subspace that captures those directions, then computes its contribution to the trace exactly. Only the remaining contribution needs random estimation. For a positive semidefinite matrix, using \(m\) products improves the relative-error bound from order \(1/\sqrt{m}\) to order \(1/m\) [MMMW, SOSA ’21]. 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
What can a handful of matrix-vector products tell us about the trace of a matrix we never hold? 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\) 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. Hutch++ removes those few large eigenvalues before probing, and the next question is how heavy the remainder can be.
The Tail Inequality
If a sketch removes the top \(k\) eigenvalues, how heavy can what remains be? A low-rank approximation can remove the leading eigendirections that make Hutchinson expensive. The following bound controls the remainder. 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}}.\]
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, a 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 Split
How do we remove the top \(k\) eigendirections when an eigendecomposition is exactly what we cannot afford? 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. 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. 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.
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.
Deflation as a Control Variate
What has Hutch++ done, in the language of variance reduction? Hutch++ 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
Could a different split, or a different algorithm altogether, beat \(1/m\)? 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.