A rate that belongs to the matrix

Every improvement anybody has made to a trace estimator is the same move, made in two different places.

Hutchinson’s estimator averages zᵀAz over random probes, and the essay that derived its two variance formulas puts them side by side: one probe of ±1 entries has variance 2(‖A‖_F² − Σa_ii²), and one standard normal probe has variance 2‖A‖_F². The difference between them is an accounting. For a ±1 vector every z_i² is one, so the diagonal term Σ a_ii z_i² is not a random variable at all — it comes back exactly, once, and the sampling is over what is left.

Hutch++ does the same thing to a different part of the matrix. A third of the budget goes into a sketch, whose range is orthonormalised into Q; tr(QᵀAQ) is computed exactly at the cost of another third; and the remaining third probes (I − QQᵀ)A(I − QQᵀ). The exact head accounts for whatever share of ‖A‖_F² the leading directions hold, and the probes see only the tail.

So both improvements have one shape — account exactly for part of the squared Frobenius norm before probing, and probe only the remainder — and in both cases what the improvement is worth is the share of that norm the accounting managed to remove. That share is a number belonging to the matrix. It is not in the method, it is not in the code, and it does not appear in either method’s name.

Which makes the headline rate a strange thing to quote. Hutchinson’s exponent is genuinely a property of the estimator: nothing about the choice of random variable changes it, and the site’s own fit over five budgets returns −0.500. Hutch++'s exponent is not a number until a matrix is named, and across five matrices of one size, from one implementation, over one budget range, it reads −7.15, −4.09, −2.55, −1.03 and −0.67, the last of those on a spectrum with no decay in it at all.

The exponent is a reading of the spectrum #

The comparison is at equal products with A, which is the only unit that means anything here, because the whole premise is an operator with no entries — a subroutine that returns Av and nothing else. A comparison in flops would assume a matrix, and a comparison at equal sample counts would let Hutch++ take three times the information and report the difference as a win. The same accounting is why a hierarchical representation built from products alone is priced in applications rather than in arithmetic.

With that fixed, the only thing left free is the matrix. Hold the size at 80, hold the seeds at 32, hold the budgets at 12, 24, 48 and 96, and move the decay.

Nothing about the second run is a different algorithm. It is the same implementation, the same seeds, the same four budgets, and a matrix whose singular values fall a little faster. The exponent moved from −2.547 to −7.154.

Three matrices, three exponents — −7.154, −4.091, −2.547 — and on the same three matrices Hutchinson reads −0.280, −0.359 and −0.538. The second sequence is not a trend. It is a fit through four medians of 32 heavy-tailed draws, which is worth about a quarter of an exponent, and every value in it is a half to within that. The first sequence moves by a factor of nearly three and is not noise at all.

At 96 products with A those three matrices separate the two methods by factors of 3.3 million, 2,260 and 55. On the fastest-decaying one, matching Hutch++'s answer with plain Hutchinson at a Monte Carlo rate would take about 10¹⁵ products with A — a number with no operational meaning, offered here only as the size of the gap the spectrum opened.

What the deflation leaves can be computed in advance #

The claim that the matrix decides it is worth more than a correlation, and the mechanism is available in closed form. Hutch++ at a budget of 96 sets its sketch width to 32 and probes with the 32 products left over. If the sketch found the exact leading 32-dimensional invariant subspace — which a random sketch does not, and this is the idealisation — the probes would face a matrix whose squared Frobenius norm is the tail of the spectrum past the thirty-second value.

That tail is a property of the spectrum alone, and its square root is what multiplies the error. Computing it on the same five matrices, beside the median relative error actually measured:

The prediction is made from the singular values and nothing else, and it tracks the measurement across nearly seven orders of magnitude, always within a factor of 3.4. It is a genuinely independent route to the same number, and it is the two-routes habit applied to a convergence rate rather than to a quantity.

It also understates every time, which is the informative half. The gap is 2.65×, 3.37×, 1.84×, 1.91× and 1.36× down the column, largest exactly where the ideal deflation is most nearly perfect. The reason is that the sketch is random: its range is close to the leading subspace and is not it, and the residue of that difference is the whole of what the probes see once the true tail has gone to 10⁻¹⁵. A guarantee stated over the draw rather than for every draw is the shape a bound that holds with probability is about, and this column is where its slack shows up as a factor of three.

It is worth being explicit about which direction would have falsified the accounting. The predicted column is what the probes would face if the sketch had found the exact leading subspace, and a random sketch cannot do better than exact. A measurement sitting below its prediction would therefore mean that the deflated remainder is not what sets the error — that some other quantity is doing the work and the agreement is a coincidence of scale. Five matrices spanning seven orders, and every one of them lands above its own prediction and within a factor of 3.4 of it.

Doing the same arithmetic for plain Hutchinson gives 4.68·10⁻², 3.47·10⁻², 2.52·10⁻², 1.93·10⁻² and 1.14·10⁻² against measurements of 5.71·10⁻², 4.59·10⁻², 2.41·10⁻², 1.69·10⁻² and 8.64·10⁻³ — ratios of 1.22, 1.32, 0.96, 0.88 and 0.76, straddling one rather than sitting above it. Hutchinson’s variance formula is exact and has no draw to be lucky in.

And the size is not the parameter #

An exponent that moves by a factor of three invites the suspicion that it is really tracking something else, and the obvious candidate is the dimension. It is not.

Holding the decay at 0.85 and running the size at 40, 80 and 160 moves Hutch++'s fitted exponent to −2.566, −2.547 and −2.606. Four times the matrix, and the third decimal place. Hutchinson on the same three reads −0.441, −0.538 and −0.611, which spans more than Hutch++ does and is entirely the fitting noise already described.

That is the honest form of the claim. The exponent is a function of the spectrum’s shape and not of its length, because the deflation depth at a given budget is a fixed number of directions and what those directions hold is decided by how fast the values fall, not by how many of them there are. At a budget of 96 the deflation removes 32 directions whether the matrix is 40 across or 160, and a geometric spectrum’s tail past its thirty-second value is the same fraction of the whole in each case, because the values beyond it are already negligible against the ones before it. A quantity that genuinely refuses to carry the dimension is the subject of the dimension does not appear, and it is worth keeping the two apart: there the independence is the theorem, and here it is a measurement that could have come out otherwise and did not.

Where the improvement is a cost #

The last stop of the sweep is the matrix with no decay at all — half its eigenvalues at +1 and the rest at −1, so the trace is small beside the Frobenius norm and there is genuinely something to estimate.

The leading 32 directions of a flat 200×200 spectrum hold 32/200 of the squared norm and no more, so the deflation removes 16% of the variance for two thirds of the budget. What is left is a method running 32 probes where the other ran 96, on 84% of the same matrix, and the measured penalty is a factor of 2.02.

The same last stop of the sweep at n = 80 reads −0.667 for Hutch++ against −0.316 for Hutchinson, and the figure above, at n = 200, reads −0.495 against −0.615. The two orderings are opposite, which is exactly what a difference smaller than the fit noise looks like and is the strongest available statement that the deflation has stopped changing the exponent. On the decaying spectra the gap was 2.0, 3.7 and 6.9.

Randomisation does not create structure made this argument for the randomised low-rank approximation and the collection has now reached it from several directions, including the sketch that is spent, where a second application of the same random map returns the zero matrix. What this sweep adds is the range. From 1.71·10⁻⁸ at a decay of 0.6 to 6.21·10⁻² on the flat spectrum above, Hutch++'s error at 96 products moves by a factor of 3.6·10⁶ — six and a half orders — while Hutchinson’s across the same matrices moves by a factor of 6.6 — and the 6.6 is itself a reading of ‖A‖_F / |tr(A)|, which falls from 0.500 to 0.163 as the decay steepens.

So the flat spectrum is not a special case; it is one end of a continuum, and the continuum is six and a half orders long. A statement of the form “Hutch++ converges at one over the budget” is true somewhere on it.

The exponent is not an exponent #

There is a sharper objection to quoting any of these numbers, and it applies to the good ones as well as the bad.

A geometric spectrum’s tail past k values falls like decay^(2k), and Hutch++'s deflation depth k is a third of the budget. So the error falls like decay^(budget/3) — exponentially in the budget, not as a power of it. Drawing that on logarithmic axes and fitting a straight line gives a slope, and the slope is a description of the window rather than of the method.

Measured on the decay-0.6 matrix, over three budget windows of the same shape:

budgets            fitted exponent
6, 12, 24, 48          −3.856
12, 24, 48, 96         −7.154
24, 48, 96, 192       −13.315

One matrix, one implementation, one set of seeds, and the exponent nearly quadruples by sliding the window one doubling at a time. At 192 products the median error is 2.63·10⁻¹⁵, which against a trace of 200 is an absolute error of 5.3·10⁻¹³ on a matrix whose Frobenius norm is 100 — about two dozen units in the last place of that norm. The deflated remainder has gone below what the arithmetic can represent, and any further fit measures rounding.

That is where this figure’s own axis stops meaning what it appears to mean, and it is the limit to carry away from the whole comparison. Hutchinson’s −1/2 extrapolates and Hutch++'s −7.15 does not. One is the exponent of a genuine power law over any budget range at all; the other is a local slope on a curve that is not a power law, and using it to predict the cost of two more digits gives an answer that is wrong in the direction of optimism. The distinction between a rate that is a constant and a rate that is an exponent is one this collection has drawn before, in a rate that does not notice the size; this is the third case, where the fitted number is neither.

Two questions come apart here that a single quoted exponent runs together. Where the two curves cross is a separate reading, it is a function of the same decay, and it moves by an order across decays that look similar — which is why the recommendation a caller can act on is stated in products with A rather than in exponents. How steeply the better curve falls once it is ahead is what an exponent describes, and the measurement above says that number is local. A method selection needs the first and the literature quotes the second.

The probe is the same claim, on the other axis #

The other half of the accounting behaves the same way, and it is cheaper to see because it has a closed form with no sampling in it at all.

The ±1 probe’s advantage over the normal probe is exactly √(1 − Σa_ii²/‖A‖_F²) — the square root of the share of the squared Frobenius norm that is not on the diagonal, which is the part the ±1 probe still has to sample. The formula contains no budget, no dimension and no algorithm. It is a ratio of two quantities read off the matrix.

At that size the tighter probe is worth 1.09× in error and 1.20× in variance, which on a method whose error falls like the square root of the budget means the looser probe needs 20% more probes for the same answer. The choice is still free and it is close to worthless.

The same family at a twentieth of the dimension is a different report.

Sweeping the size across the generator’s five stops, the ratio reads 0.6399, 0.7520, 0.8213, 0.8864 and 0.9150 at n = 10, 20, 40, 80 and 120 — worth 1.56× at the small end and 1.09× at the large one, or 2.44× and 1.20× as variances. The quantity moving underneath is the diagonal’s share of the squared Frobenius norm, which falls from 0.5905 to 0.1629 across the same five sizes, because this family’s singular values fall like 1/j while its dimension grows.

That is the point the two halves of this essay share. The often-quoted figure for the ±1 probe is that it halves the variance; the value measured here runs from a factor of 2.44 down to 1.20, and the extreme the derivation actually permits is a diagonal matrix, where the share is 1 and the ±1 probe is exact from a single draw. A ratio that runs from 1 to infinity has no representative value, and a matrix’s diagonal share is not something a caller who reaches for a trace estimator can see, since the reason for reaching is that the entries do not exist.

Both endpoints are reachable and neither is contrived. A diagonal share of 1 is any diagonal matrix, on which the ±1 probe returns the trace exactly from a single draw, with a variance the closed form gives as exactly zero. A share of 0 is any matrix with a zero diagonal — an adjacency matrix, for instance — where the two distributions have identical variance and the choice of probe is genuinely a detail. The advantage is bounded below by nothing and above by nothing, and what decides it is the one object the method was chosen for not having.

What follows for anything that quotes a rate #

A rate quoted without a matrix is a rate quoted about somebody else’s matrix. That is the whole of it, and the practical consequences are specific.

For a method selection, the sketch is the measurement. Hutch++ computes the singular values of its own sketch on the way to Q and discards them, and those values are the cheapest available reading of the decay that decides whether the deflation was worth its two thirds. Every randomised method here computes the evidence for whether it was the right method and reports none of it, which is the same gap that makes an answer that changes with the seed a property worth reporting rather than a nuisance to be averaged away.

For an eigenvalue count, the matrix is chosen rather than given. Counting what is inside a circle estimates the trace of a spectral projector, whose spectrum is a step: near one inside the contour and near zero outside. A step has a dominant subspace exactly the size of the step, which is the best case in the table above and is known in advance from what the operator is.

For a parameter-choice rule, the estimate is inside an argmin. The trace that choosing without knowing needs is the denominator of generalised cross-validation, and the step that stops mattering found the same quantity bounded by the Krylov dimension rather than by the size. The influence matrix’s spectrum is a step as well, so a deflating estimator is the right one there — and the reason is a fact about regularisation, not a fact about the estimator.

And for anything a sketch is asked to decide, the share is the whole answer. The half of a problem a sketch may touch draws the boundary by asking what a multiplicative guarantee can say about a quantity that is zero. The boundary here is the same one seen from the other side: a multiplicative guarantee on a probed remainder says nothing until the remainder is named, and naming it requires the spectrum that was the reason for not forming the matrix.

What is worth carrying #

One estimator here has a convergence rate and the other has a rate per matrix. Hutchinson’s −1/2 is stable across every spectrum measured and across a factor of four in the dimension; Hutch++'s runs from −0.67 to −7.15 on five matrices of one size, and the deciding parameter is the share of the squared Frobenius norm the leading directions hold.

The share can be computed from the spectrum and predicts the measurement, within a factor of 3.4 over nearly seven orders, always low, and the shortfall is the price of a random sketch rather than an error in the accounting.

A fitted exponent on a geometric decay is a description of the budget window. Sliding it one doubling took −7.154 to −13.315 on one matrix, so the number is a report and not a prediction.

And the ±1 probe’s advantage is the same quantity read on the diagonal, running from 2.44× to 1.20× in variance across five sizes of one family. Both improvements account exactly for part of ‖A‖_F² and probe the rest, both are worth the square root of what they removed, and how much that is was decided before either method was chosen.