Counting what cannot be looked at

The trace is the cheapest quantity in linear algebra. Add up the diagonal: n additions, no multiplications, exact.

It is also one of the most expensive to estimate, because the matrices whose trace is wanted are never the ones a caller has.

• tr(A⁻¹) needs n solves, and appears in every uncertainty quantification that wants the variance of an estimator.
• tr(log A) = log det A is the log-determinant, which is the normalising constant of a Gaussian and therefore in the likelihood of every Gaussian process.
• tr(f(A)) for other f — the number of eigenvalues in an interval, the Estrada index of a graph, a partition function.
• And tr(H(λ)), the effective number of parameters a regularisation leaves, which is the denominator this site has already measured in generalised cross-validation.

In every one of those, a product with the operator is available and the entries are not.

The identity #

Let z have mean zero and identity covariance — E[z] = 0, E[zzᵀ] = I. Then

E[zᵀAz]  =  E[tr(zᵀAz)]  =  E[tr(A zzᵀ)]  =  tr(A · E[zzᵀ])  =  tr(A)

for every such z, and every A. No bias, no assumption on the matrix, no asymptotics. Take s independent probes, average, done.

What it has is variance, and the variance is the whole subject.

Two distributions, and one is never worse #

For z with independent ±1 entries — a Rademacher vector — the variance of one probe is

$$\operatorname{Var} = 2 \sum_{i \ne j} a_{ij}^2 = 2 \left( \|A\|_F^2 - \sum_i a_{ii}^2 \right)$$

and for standard normals it is 2‖A‖_F². Both are exact for every n, not asymptotic.

Which makes the diagonal case a zero and the rotated case a constant fraction, and the slider says both.

At n = 8, 20, 40, 80 and 200 the ±1 variance is 0, 0, 0, 0 and 0 and the normal variance is 97.38, 265.2, 545.1, 1,105 and 2,785. Divided by n those are 12.2, 13.3, 13.6, 13.8 and 13.9 — the second formula, 2‖A‖_F², evaluated on a diagonal whose entries grow with the size, and the first formula’s Σ(i ≠ j) aᵢⱼ² having nothing in it to sum.

And the rotated case is the one that says the two formulas are the right two. Rotate the same spectrum so the matrix is no longer diagonal and the ±1 variance stops being zero: 10.05, 26.42, 57.17, 116.6 and 297.1 across the five sizes. Divided by the normal probe’s variance at the same size those are 0.103, 0.0996, 0.105, 0.106 and 0.107 — a constant tenth, which by the two formulas means 89.5% of the Frobenius mass is still on the diagonal after the rotation, at every size. The ±1 probe pays only for what is off the diagonal, and the rotation moved a tenth of the matrix there.

So a fixed number of probes gets relatively more accurate as the matrix grows, which is the opposite of most things in this collection. The normal probe’s error after sixty reads 9.9%, 5.1%, 2.9%, 0.57% and 0.86% across the five sizes: the standard error grows like √n and the trace grows like n, so their ratio falls like 1/√n. Trace estimation is a method whose accuracy improves on exactly the problems it exists for.

The ratio is 1 − Σaᵢᵢ²/‖A‖_F², a number between 0 and 1 that is a property of the matrix. So the Rademacher probe is never worse and usually better, and the reason is one line: `zᵀAz = Σᵢ aᵢᵢ zᵢ²

• Σ(i ≠ j) aᵢⱼ zᵢ zⱼ, and for a ±1 vector zᵢ² = 1`, so the diagonal contributes a constant rather than a random variable. The diagonal is estimated exactly and the off-diagonal is what is being sampled.

assertTheVarianceIsTheFormula checks both formulas against the empirical variance over 4,000 probes, and checks the ratio against the quantity computed from the matrix at 10⁻¹² — because the ratio has no sampling in it at all and is an algebraic identity between two closed forms.

The extreme case #

Push that observation to its end. On a diagonal matrix there is no off-diagonal to sample: zᵀAz = Σ aᵢᵢ zᵢ² = Σ aᵢᵢ = tr(A) for every ±1 vector.

One probe returns the exact trace. The variance formula says the same thing — ‖A‖_F² = Σaᵢᵢ² for a diagonal matrix, so it returns exactly zero — and the two routes are asserted against each other.

The Gaussian probe on the same matrix has variance 2Σaᵢᵢ² = 545 and needs thousands of samples for two digits. The estimator that is exact on the easiest class is the one nobody’s notation distinguishes, and the difference between the two is one line in an implementation.

There is a control beside it, and it matters: the same spectrum in a general basis — QΛQᵀ with the same eigenvalues — has a Rademacher variance of 56.8. The advantage is a property of the basis, not of the spectrum, and a reader who took away “±1 probes are free” would have taken away something false.

The zero that was not zero #

A small thing worth recording, because it is this site’s own subject appearing inside its own library.

Evaluating 2(‖A‖_F² − Σaᵢᵢ²) on a diagonal matrix subtracts two quantities that are equal, and returned 1.14·10⁻¹³ where the answer is zero. The formula is written as a difference and the two terms cancel completely; what is left is the rounding of the two sums.

rademacherVariance now sums the off-diagonal entries directly — which is what the formula is counting — and returns exactly 0. It is cancellation takes the answer in a routine whose whole purpose is to compute a variance, and it is the kind of thing that would have been invisible had the assertion been written with a tolerance instead of as an equality.

The rate, and the one thing that beats it #

The relative error falls like s^(−1/2) and nothing about a better random variable changes the exponent. assertTheRateIsMinusAHalf fits it over five budgets and 40 seeds and gets −0.500.

That is the reason trace estimation has a reputation for being expensive: two more digits costs ten thousand times the work, and a Monte Carlo rate is a Monte Carlo rate whatever the variance constant is.

Hutch++ changes the exponent. Spend a third of the budget on a sketch, orthonormalise it, compute tr(QᵀAQ) exactly at another third, and Hutchinson the deflated remainder (I − QQᵀ)A(I − QQᵀ) with what is left. The deflated matrix has a far smaller Frobenius norm whenever the spectrum decays, and the variance is proportional to that norm squared.

The comparison is at equal products with A, and the distinction is not pedantry. Compared at equal sample counts, Hutch++ would be an algorithm competing against itself with three times the information, and the figure would be a picture of arithmetic.

And the flat spectrum, for the third time #

Randomisation does not create structure made this argument for the randomised SVD, and a bound that holds with probability made it for the HMT bound. This is the third method to reach it, and it is asserted here rather than referred to, because a claim that holds for a low-rank approximation is not a claim about a trace estimator.

The reason is measured too: at k = 8 out of n = 120 the leading directions hold k/n of the Frobenius norm and no more, which is what a flat spectrum means, so the deflation removes 6.7% of the variance for 67% of the budget.

One detail is worth keeping rather than smoothing away: on the decaying spectrum Hutch++ is also behind at the smallest budget. Two thirds of a budget of twelve is a sketch of four columns, which deflates almost nothing while removing eight probes. Hutch++ has a break-even and it is not at zero, and the assertion says so in both directions.

Where the break-even is #

That is asserted at one decay and one budget, and it is the quantity a caller has to know. Sweeping both — the budget at which Hutch++'s median relative error over twenty-four seeds first falls below plain Hutchinson’s, on a 60 × 60:

   decay      break-even, in products with A
   0.70            12 or below
   0.80            12 or below
   0.85            18
   0.90            48
   0.95            72
   0.99            past 144

An order of magnitude across the same slider the figures above already carry. At a decay of 0.99 — a spectrum nobody would call flat, and one whose deflation does eventually pay — Hutch++ is behind at every budget up to 2.4n products with A.

That changes what the exponent comparison is worth. A better exponent guarantees an eventual crossing and says nothing at all about where it is, and the where is the whole decision for anybody with a budget. On the 0.95 figure above, the fitted exponents are −0.50 and −1.03 and the crossing is at 72 products: over most of that figure’s own budget range the method with the better exponent is the worse method.

So the practical form of the recommendation is not use Hutch++ when the spectrum decays. It is use Hutch++ when the budget is past the break-even, and the break-even moves by an order over decays that look similar. A caller who cannot see the spectrum — which is the whole premise of this essay, since a product with A is available and the entries are not — cannot read that table either, and the honest fallback is to run both at a small budget and compare, which costs one extra sketch and settles it.

That fallback is cheaper than it sounds and it is worth spelling out, because it turns an unanswerable question into a two-line one. Run Hutchinson at a budget of b and Hutch++ at the same b; the second already computes tr(QᵀAQ) exactly on its sketch, so the comparison costs nothing beyond what was going to be spent anyway. If Hutch++ is ahead, the crossing is behind and more budget will widen the gap. If it is behind, either the budget is short or the spectrum is flat, and the sketch’s own singular values say which — they are computed and discarded.

Which is the general shape this collection keeps arriving at with randomised methods. The seed is not the problem and the spectrum is, and the sketch a randomised method takes is also the cheapest available measurement of the spectrum it was betting on. Every one of these algorithms computes the evidence for whether it was the right algorithm, and none of them reports it.

What this does to a parameter-choice rule #

Choosing without knowing measured generalised cross-validation against three other rules for picking a regularisation parameter, and the step that stops mattering found its denominator to be a trace over a Krylov subspace, bounded by k + 1 whatever n is.

GCV’s whole objective is

G(λ)  =  ‖(I − H(λ))b‖²  /  (tr(I − H(λ)))²

and the denominator is exactly the kind of trace this essay is about — the quantity choosing without knowing needs and cannot afford: H(λ) = A(AᵀA + λ²I)⁻¹Aᵀ is never formed, a product with it is a solve, and n solves to get the diagonal defeats the purpose of choosing λ cheaply.

So GCV as it is actually implemented on a large problem is GCV with an estimated denominator, and the estimate is Hutchinson’s. That puts a random quantity inside the objective being minimised, which has two consequences worth naming.

The objective is noisy, so its minimiser is noisy. A rule that selects λ by minimising a function computed to two digits selects a λ known to about two digits, and the site’s measurement that GCV’s λ barely moves — a factor of 2.15 while its denominator moved by 27 — is the reason that is survivable here and would not be everywhere.

And the noise is not independent across λ. Using the same probe vectors at every λ makes the estimated objective a smooth function of λ, wrong by a consistent amount; drawing fresh probes at every λ makes it unbiased and jagged, and a minimiser of a jagged function is a minimiser of the jaggedness. Every implementation reuses the probes, which is the right choice and is a decision about correlation rather than about variance.

That is worth sitting with. The parameter-choice rule this site measured against an oracle is, at scale, a rule whose objective contains a Monte Carlo estimate held deliberately correlated so that its error cancels in the argmin. None of that appears in the formula.

Where the flat spectrum is not the pathological case #

It is worth resisting the reading that a flat spectrum is a contrived difficulty. For the traces this essay opened with, it is often the typical one, which is a bound that holds with probability read the other way.

tr(A⁻¹) for a discretised elliptic operator has a spectrum that decays like k^(−2/d) — slowly, and more slowly the higher the dimension. The influence matrix H(λ) whose trace is the effective degrees of freedom has eigenvalues σᵢ²/(σᵢ² + λ²), which are near 1 for the retained directions and near 0 for the rest: a step, not a decay, and a step has a dominant subspace exactly the size of the step, which is what makes it the one case where deflation is spectacular.

So whether to deflate is a question about the operator, and it is usually answerable in advance from what the operator is, which is the useful form of the result.

Why the trace and not the diagonal #

A question worth answering directly, because the diagonal is what the trace is made of and estimating it entrywise is a different and harder problem.

e_iᵀA eᵢ gives aᵢᵢ exactly, in one product with A, so n products give the whole diagonal and the trace with it. The estimator in this essay uses far fewer than n and returns only the sum. What it is buying is precisely the s ≪ n, and what it gives up is any statement about the individual entries.

The estimator for the diagonal — the same probes, but forming z ⊙ Az and averaging entrywise rather than summing — is the natural extension, and its variance per entry is set by the off-diagonal mass in that entry’s row rather than by the whole matrix’s. So a matrix with a heavy row is a matrix whose diagonal is estimated badly in exactly one place, and there is no averaging over n to rescue it.

That asymmetry is why the trace is the quantity with a literature and the diagonal is a special request. The sum of n things each estimated to relative accuracy ε is not estimated to ε; the sum estimated directly is, because the cancellation the probes suffer is over the off-diagonal and the diagonal contributes to every probe identically.

Where the variance actually comes from #

One more decomposition, because it makes both findings in this essay obvious rather than surprising — and because randomisation does not create structure is the sentence it turns into algebra.

Write zᵀAz = Σᵢ aᵢᵢ zᵢ² + Σ(i ≠ j) aᵢⱼ zᵢ zⱼ. The first term is what the probe measures exactly in the ±1 case and noisily in the normal case; the second is the sampling.

• For ±1, zᵢ² = 1, so the first term is the constant tr(A) and the whole variance is the second: 2Σ(i ≠ j) aᵢⱼ².
• For normal, E[zᵢ²] = 1 but Var[zᵢ²] = 2, so the first term contributes 2Σ aᵢᵢ² on top.

Everything else follows. A diagonal matrix has no second term, so the ±1 variance is zero and the normal one is 2Σaᵢᵢ². A matrix with a large diagonal relative to its Frobenius norm has a large advantage; a hollow matrix has none at all, because there is no first term for either probe to get right.

That last case is worth naming: for a matrix with zero diagonal the two distributions have identical variance, and the ±1 probe’s whole advantage disappears. An adjacency matrix is such a matrix, and counting triangles in a graph is tr(A³)/6, so the one application where the trace estimator is best known is the one where the choice of probe does not matter.

What the budget is measured in #

The horizontal axis of the rate figures is products with A, and it is worth saying why that is the honest unit rather than a fussy one.

For a dense matrix a product costs n² operations and every method here costs the same per product, so products and flops agree up to a constant. For the matrices this estimator exists for they do not. A sparse operator’s product costs its number of nonzeros; a matrix-free operator’s product is a PDE solve, an FFT, or a simulation, and can cost anything at all.

So a comparison in flops would be a comparison that assumed a matrix, and the whole premise is that there is not one. A comparison in probes would be worse: it would let Hutch++ take three times as much information as Hutchinson and then report the difference as an improvement.

Counting products is also what makes the break-even visible. Hutch++'s sketch and exact head cost two thirds of its budget before a single probe is taken, so at a budget of twelve it is behind, and the figure says so rather than starting the axis past the crossing.

What is worth carrying #

An unbiased estimator is exactly right on average and wrong every time, and the whole of the engineering is in the variance rather than in the bias.

The choice of probe distribution is a decision with a closed form behind it. ±1 is never worse than normal, better by a factor that is a property of the matrix, and on a diagonal matrix it is exact.

The rate is Monte Carlo and no random variable improves the exponent. Deflation does, and only where there is a dominant subspace to deflate, and it has a break-even.

And a flat spectrum defeats it, for the third time on this site. What randomisation buys is fewer passes over the data, and on a matrix with no structure to find it buys nothing.

Another quantity computed from products alone #

A trace estimator asks the matrix only for products, which is what makes it usable on an operator nobody can look at. So does the action of a matrix function, and so does an orthogonalisation that never inverts anything.