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

Var  =  2 · Σ_{i≠j} a_ij²  =  2(‖A‖_F² − Σ a_ii²)

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

The ratio is 1 − Σa_ii²/‖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 = Σ_i a_ii z_i²

• Σ_{i≠j} a_ij z_i z_j, and for a ±1 vector z_i² = 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_ii z_i² = Σ a_ii = tr(A) for every ±1 vector.

One probe returns the exact trace. The variance formula says the same thing — ‖A‖_F² = Σa_ii² 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_ii² = 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_ii²) 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.

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: 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.

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 σ_i²/(σ_i² + λ²), 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_i gives a_ii 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.

Write zᵀAz = Σ_i a_ii z_i² + Σ_{i≠j} a_ij z_i z_j. 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_i² = 1, so the first term is the constant tr(A) and the whole variance is the second: 2Σ_{i≠j} a_ij².
• For normal, E[z_i²] = 1 but Var[z_i²] = 2, so the first term contributes 2Σ a_ii² 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_ii². 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.