Counting what is inside a circle

Counting what cannot be looked at is about a trace that has to be estimated because the matrix it belongs to is never formed. The estimator is Hutchinson’s: probe with random vectors, average zᵀAz, and the variance of the estimate is what the choice of probe distribution decides.

This essay is the same device with a stranger integrand and a much stranger answer. The trace being estimated is of T(z)⁻¹T′(z) for a matrix-valued function T; it is integrated around a closed contour; and what comes out is not a real number to be reported with an error bar but an integer — the number of eigenvalues inside.

Why the trace of that thing is a count #

The argument principle says that for a function f analytic inside a contour except at its zeros, (1/2πi)∮ f′/f dz counts the zeros inside, with multiplicity. Take f = det T. Its logarithmic derivative is

(det T)′/det T = tr(T(z)⁻¹T′(z)),

by Jacobi’s formula, so the count is an integral of a trace and no determinant is ever formed — which matters, because a determinant of a matrix of any size is the number that decides nothing: it overflows before it is uninformative and is wrong long before it is absent.

The integrand costs n linear solves per quadrature point, in complex arithmetic, and it is the whole computation.

The convergence is exponential, which is why anyone bothers #

On a circle, the integrand is analytic and periodic, and the trapezoidal rule on a periodic analytic function converges geometrically in the number of points rather than at the algebraic rate a Newton–Cotes rule gets on an interval. Measured, on a circle of radius two containing two eigenvalues of a delay problem:

8 points      5.7·10⁻²
16            1.5·10⁻²
32            5.3·10⁻⁴
64            3.6·10⁻⁷
128           1.4·10⁻¹³

Each doubling buys more than the last. That is the entire economic case for contour methods: a hundred and twenty-eight solves of an n × n complex system buy thirteen digits of a quantity that would otherwise require finding the eigenvalues.

The rate depends on the distance from the contour to the nearest singularity, not on how many singularities are inside. A circle that passes close to an eigenvalue has a nearly singular integrand and needs several times as many points — at radius 1.5, which runs close to an eigenvalue at 1.587, a hundred and twenty-eight points reach only 7.6·10⁻⁴. Nothing about the count is harder there; the quadrature is.

And then it is rounded, which is where the trouble is #

The count is an integer, so the natural last step is to round. The natural last step throws away the only diagnostic available.

At eight points on the radius-two circle the rounded answer is already right while the integral is 0.06 away from it — correct by luck, with two digits of margin. On a circle of radius four, which passes closer to the spectrum, eight points return a count of minus two where the answer is four.

A negative number of eigenvalues, from a routine that did not fail, warn or diverge. Nothing in the returned integer says which of those two situations it came from.

Across the four contours the wrong-count rate does not follow the error at all. At radius 1.5 the finest error is 7.58·10⁻⁴ and none of the six rules is wrong; at radius 3 the error is a thousand times smaller, at 7.21·10⁻⁷, and one is; at radius 5 the error is 1.43·10⁻⁵ and three are. The ordering by error and the ordering by correctness are different orderings, which is the whole reason the distance to the nearest integer has to be reported rather than the error. A quadrature that is accurate to seven digits and half an integer away from a whole number is worth less than one accurate to three digits and a hundredth away. The distance from the nearest integer does, it costs nothing, and it is the quantity this essay plots — for the same reason the residual the method reports is worth printing beside an iterate: the number a method returns and the number that says whether to believe it are different numbers.

The probes, and three ceilings #

Counting is not usually the goal. Beyn’s method uses the same contour and the same solves to produce the eigenvalues themselves: form the moment

A₀ = (1/2πi) ∮ T(z)⁻¹V dz

for a random n × ℓ block V, and the eigenvalues fall out of its rank and its singular vectors. The random block is doing exactly what a sketch does in the dimension does not appear — capturing a subspace by hitting it with something generic — and it inherits the same failure mode.

Measured with four eigenvalues inside: the rank is 1, 2, 3, 4, 4 at ℓ = 1 to 5.

Ceiling one: the probe count. A caller who asks for two eigenvalues gets two — both genuine, both with small residuals — and nothing reports that there were four. This is the same shape as the sketch that is spent: a random probe carries a fixed budget of information, and asking it for more than it holds returns something plausible.

Ceiling two: the number actually inside. Past ℓ = 4 the rank stops rising, which is the good case — the method has everything and extra probes cost solves for nothing.

Ceiling three, and it is the one nobody predicts: the size of the matrix. The moment is n × ℓ, so its rank cannot exceed n. On a circle of radius five containing twelve eigenvalues of a four-by-four problem, the rank stops at four however many probes are used. The first moment cannot represent more than n eigenvalues; recovering the rest needs higher moments — ∮z^k T⁻¹V dz — rather than a wider block, and a caller who responds to a saturated rank by adding probes is buying nothing at all.

Three ceilings, one number, and no way to tell them apart from the result. The counting integral distinguishes them, it uses the same integrand, and it costs one more pass round the same contour.

A fourth ceiling, and it pushes the other way #

All three ceilings above cap the rank, so the decision procedure they support assumes the rank can only be too small. It can be too large, and the case is not exotic.

Take the three-by-three delay problem and the circle of radius two. Exactly one eigenvalue is inside — the counting integral returns 1.0000003 at 256 points, which is as unambiguous as an integer gets. The numerical rank of the moment with two probes is two:

The extra direction is not an eigenvalue and it is not a probe artefact. It is the quadrature’s own error, leaking a contribution from the eigenvalue at 2.1200 — which sits just outside the contour. The trapezoidal rule’s error against that singularity decays like e^(−Ng) with g = |ln(2/2.12)| = 0.0583, so the spurious singular value drops below a relative cut of 10⁻⁸ only past N = ln(10⁸)/0.0583 = 316 points. The measured crossover is between 256 and 512, which is where that number says it should be.

The tell is the shape of the table rather than any one entry: the rank falls when points are added. No ceiling can do that. A ceiling is a bound on what the moment can represent, and adding quadrature points does not change what a rank-one object can hold — so a rank that drops is a rank that was never measuring the object it was supposed to.

That is where the decision procedure misfires, and it misfires in the expensive direction. With ℓ = 2 probes and a true count of 1, the observed rank equals ℓ, which the procedure reads as rank = ℓ < N, the probe block is the binding ceiling; widen it. Widening it buys nothing at all, because the block was never the constraint. The repair is more quadrature points, and the procedure has no case that points there.

So the three ceilings are right about the limit and silent about the approach. The complete statement is that the rank is min(ℓ, N, n) once the quadrature has converged, and until then it is that plus however many nearby exterior eigenvalues the tolerance cannot exclude — which is a number set by the same g as everything else on the contour, and which is therefore knowable before the moment is formed.

It also says something about the tolerance, which the three-ceiling analysis treats as a detail. A rank cut of 10⁻⁸ on this contour is a demand for 316 quadrature points; a cut of 10⁻¹² is a demand for 474, and one of 10⁻⁴ for 158. So the tolerance and the point count are one decision rather than two, tied by g, and setting them independently is how a caller ends up paying for points that cannot help or reading a rank the points cannot support. Nothing in the returned object says which of those two happened, which is the same complaint the rounded count answers and for the same reason.

There is a cheap diagnostic in it, and it is the one this essay already argues for elsewhere. Two integers from one contour distinguishes three ceilings; three integers distinguishes four. Compute the rank at N points and at 2N, and a rank that falls is the fourth ceiling announcing itself for the cost of one more pass — the same pass the doubling was going to cost anyway. assertTheRankHasAFourthCeilingWhichIsTheQuadrature holds the whole of it: that the count inside is one, that the rank is two below 316 points and one above it, that it falls when points are added, and that the crossover is where ln(1/tol)/g puts it.

What the integrand costs, and what it is compared against #

Every quantity above rests on the trace being right, and the trace is the one thing here that could be silently wrong: a complex solve at a point where T is nearly singular, n of them per quadrature point, summed around a contour, with cancellation between contributions of opposite sign.

So it is computed twice. The first route embeds the complex system in a real one of twice the size — [[Re, −Im], [Im, Re]] against [xᵣ; xᵢ] — and calls the same pivoted elimination everything else on this site calls. The second uses the fact that, for the family drawn here, T is a symmetric matrix plus a scalar multiple of the identity, so T⁻¹T′ has eigenvalues in closed form and the trace is a sum of n divisions.

They agree to 5.6·10⁻¹⁵ at sixteen points around a contour. One is sixteen solves of an 8 × 8 real system; the other is sixteen sums of four quotients. Nothing is shared, which is what makes the agreement mean something — and the second route exists only for the restricted family, which is exactly why the family was restricted.

The cost, counted honestly #

A contour method looks expensive and the arithmetic says otherwise, provided the right thing is being counted.

Per quadrature point: one factorisation of T(z), which is complex and n × n — so, through the real embedding, one factorisation of a 2n × 2n real matrix, at eight times the cost of an n × n one. Then ℓ solves for the probe block, or n for the trace. At 32 points that is 32 factorisations, which for n in the thousands is the whole budget.

Against what? For a problem with infinitely many eigenvalues there is no alternative to compare with, which is the honest answer and also the uncomfortable one: this is not a faster route to something a direct method could have done, it is the only route to a question a direct method cannot pose. The comparison worth making is with the other narrow question — Newton from a starting guess, at one factorisation per step, converging quadratically to one eigenvalue that may not be the one it was pointed at. Thirty factorisations for every eigenvalue in a region, against five for one eigenvalue near a point, and the choice is decided by how many are wanted rather than by any rate.

Two integers from one contour #

A detail worth having, since it costs nothing and is the practical use of the count.

The moment A₀ gives the eigenvalues; the count N gives how many there should be. Running both means the rank of A₀ can be compared with an independently computed integer, and the three ceilings above become distinguishable:

• rank < ℓ and rank = N — the method has everything and the probes were sufficient.
• rank = ℓ < N — the probe block is the binding ceiling; widen it.
• rank = n < N — the size is the binding ceiling; higher moments are needed and a wider block will not help.

That is a decision procedure made out of two integers, both exact, from one pass of the same factorisations. It is the same pattern an eigenvalue count that cannot be slightly wrong sets out for a symmetric matrix: an integer that is either right or wrong by a whole unit is a far better diagnostic than a real number that is approximately right, because it can be checked against another integer rather than compared against a tolerance.

What the randomisation is and is not doing #

Worth separating, because this method is randomised in a narrower way than most of the field.

The random block is not reducing a dimension. It is not a Johnson–Lindenstrauss embedding and there is no distortion bound: nothing here is claiming that a subspace is approximately preserved. What V does is guarantee, with probability one, that the columns of the moment span the whole eigenspace the contour encloses — a genericity argument rather than a concentration one.

So the quantity that behaves probabilistically is rank, not accuracy. A badly chosen V does not give a slightly worse answer; it gives an answer of the wrong dimension, which is a discrete failure and is exactly what randomisation does not create structure warns about from the other direction. The eigenvalues that do come back are as accurate as the quadrature — measured against a Lambert closed form, to the rounding level once the contour is resolved.

And the seed matters in the way an answer that changes with the seed describes: two runs with different V return the same eigenvalues to the rounding level when ℓ is large enough, and different subsets when it is not. A method whose answer changes with the seed in its cardinality is a different animal from one whose answer changes in its fifteenth digit.

Where the estimate is exact and where it is not #

One more separation, because this essay sits in a field about bounds that hold with probability.

The count is not an estimate. Once the quadrature has converged it is exact — an integer, obtained by rounding a number that is within 10⁻¹³ of it, with the distance printed. Nothing about it holds “with probability”, and no probe is involved: the count needs the trace, and the trace is computed by n solves rather than sampled.

The eigenvalues are where the randomness enters, and there the guarantee is the discrete one above. So a run of Beyn’s method contains a deterministic exact integer and a randomised set whose size is bounded by it, and the useful discipline is to compute the first and check the second against it. That is the shape a bound that holds with probability recommends in general — know which of the numbers in hand is the guaranteed one — arriving here with an unusually clean division between them.

What this does to the field’s own habit #

The randomised field on this site has a discipline: state what is guaranteed, state with what probability, and measure whether the bound is attained. Every essay in it — a bound that holds with probability onwards — reports a failure probability and a distortion, because those are the currency.

This method does not fit that currency, and the mismatch is worth recording rather than smoothing over. There is no ε, no δ, and no failure probability that shrinks with oversampling. There is a genericity argument that holds with probability one, a discrete rank that is either sufficient or not, and a deterministic integer beside it. Oversampling — the standard lever, priced in built from products alone — buys nothing here once ℓ exceeds the count, and buys everything until it does. The curve is a step rather than a slope.

That is a genuinely different shape of randomised method, and the reason it is worth its own essay in this field rather than a paragraph in the polynomial one.

Where the contour goes, which is the only real decision #

Everything else here is machinery; the placement of the contour is judgement, and it decides both of the quantities the method returns.

It sets ρ, the distance to the nearest eigenvalue, and therefore the convergence rate — measured at nine orders of difference between a circle in a gap and one across it. It sets how many eigenvalues are inside, and therefore the probe count needed. And it sets what the answer means: a contour drawn around a region of physical interest returns a physically meaningful count, and one drawn for numerical convenience returns a number about a disc.

The uncomfortable part is that all three want different things. A contour placed to avoid the spectrum converges fast and may enclose the wrong set; one placed around the interesting region may pass arbitrarily close to an eigenvalue. There is no rule, and the honest practice is the one the measurement supports: place the contour for the question, then read the distance from the nearest integer, then move it if that distance is not small.