An eigenvalue count that cannot be slightly wrong

Everything this site has done with eigenvalues returns real numbers with errors in them. A Jacobi sweep returns a spectrum to fourteen digits. The QR algorithm returns one to thirteen. Lanczos returns one that is right about some eigenvalues and confidently wrong about a repeated one. In every case the output is a float and “is it right” is a question about a tolerance.

There is one spectral computation here whose answer is an integer.

The rule #

For a symmetric A and any shift σ, take an LDLᵀ factorisation of the shifted matrix A − σ·(the identity), with no pivoting at all — 1 × 1 pivots in the natural order. Then

ν(σ) = the number of negative entries of D = the number of eigenvalues of A strictly below σ

That is Sylvester’s law of inertia used as an algorithm. The factorisation is a congruence, A − σI = LDLᵀ with L nonsingular, and congruence preserves the signs of the eigenvalues; the eigenvalues of a diagonal matrix are its entries; so the sign pattern of D is the inertia of A − σI, which counts the eigenvalues of A on either side of σ.

The cost is one elimination — O(n³) dense, O(nnz) for a tridiagonal, and for a sparse matrix whatever a symbolic factorisation says. What is returned is not an approximation to a count. It is a count.

Why an integer is a different kind of answer #

A computed eigenvalue can be a little wrong, and how little is governed by that eigenvalue’s own condition number. A computed count cannot be a little wrong. There is no number between 6 and 7 for a rounding error to land in.

So the failure mode changes shape entirely. Where a computed spectrum degrades smoothly as the problem gets harder, a computed count is exactly right until it is wrong by one — and the question stops being “how accurate is it” and becomes “under what circumstances is it wrong at all”.

The answer, from the standard analysis, is that the floating-point count is the exact count of some matrix within a small multiple of n‖A‖u of A. So ν(σ) can only differ from the truth when σ is within that distance of an eigenvalue: a band, whose width is a property of the matrix’s norm and the format’s precision and of nothing else. Away from the band the answer is not approximately right; it is right.

What the staircase is, read carefully #

The figure repays a slow look, because it contains two computations that share nothing and one of them is drawn as a curve rather than as a set of points.

The staircase is ν evaluated at four hundred and twenty shifts across the interval. At each of them an elimination runs, the signs of its diagonal are counted, and an integer is plotted. The curve is therefore not an interpolation of anything: every point on it is a separate factorisation, and the vertical segments are where two adjacent shifts straddled an eigenvalue.

The vertical marks are the eigenvalues a Jacobi sweep returns, and the dots at (λₖ, k) are what the two computations have to agree about — the k-th eigenvalue is where the staircase reaches height k. That is the agreement, and there is a mark at every step and a step at every mark.

What the figure cannot show is the scan, which is two thousand shifts rather than four hundred and twenty and is placed at random rather than on a grid. A grid can miss a narrow feature by landing either side of it; random placement over two thousand draws does not, and the scan is what the assertion checks.

The scan, which is the claim #

A figure showing a staircase that lines up with some marks would be satisfied by a routine that was right most of the time. What is asserted here is a scan: two thousand shifts placed at random across an interval containing the whole spectrum, on each of three spectra — evenly spread, geometrically graded over six decades, and split into two clusters with a gap.

Six thousand shifts. Zero disagreements with a spectrum computed by a Jacobi sweep, which shares no code with the counting routine and no algorithm either.

That is the site’s two-routes habit in the form it takes when one of the two answers is exact. A disagreement anywhere would be a whole eigenvalue, so the comparison has no tolerance in it, and “they agree” means what it says.

Where it does fail, measured #

Two eigenvalues at 1 and 1 + gap, six others around them, and a hundred shifts placed strictly between the pair — where the count must read three:

gap shifts counted wrongly 10⁻³ 0 of 100 10⁻⁹ 0 10⁻¹² 0 10⁻¹³ 1 10⁻¹⁴ 19 10⁻¹⁵ 100

The transition is between 10⁻¹³ and 10⁻¹⁴, and n‖A‖u for this matrix is 4·10⁻¹⁵ — the boundary is where the analysis puts it, within an order.

And the failure is total when it comes. At a gap of 10⁻¹⁵ not one shift in a hundred reads three; the routine reports two or four, consistently, because the two eigenvalues it is being asked to separate are one eigenvalue as far as the format is concerned.

The refusal this page publishes is the reading that the integer makes the answer safe: a count computed in floating point is right, because a count is an integer and integers do not have rounding errors. The premise is true and the conclusion does not follow. The assertion is fed the run at a gap of 10⁻¹⁵ and required to fail.

The band scales, which is what makes it a property rather than a tolerance #

Scaling the whole matrix by 10⁶ scales the eigenvalues and the band together, and the measured curve moves right by six decades exactly. Nothing in the routine has a constant in it: there is no tolerance to set, no 1e-12 anywhere, and the place where it stops working is a consequence of the matrix’s norm and the format’s precision.

The quantity behind that is n‖A‖u, and it tracks the scaling exactly.

n‖A‖u is proportional to the scale over five decades and the failure gap is not. At scalings of 10², 10³, 10⁴, 10⁵ and 10⁶ the quantity runs 4,807, 4.854·10⁴, 4.859·10⁵, 4.859·10⁶ and 4.859·10⁷ — a factor of ten per decade, to four digits from the second stop onward — while the gap at which the first miscount appears sits at 10⁻¹³ or 10⁻¹⁴ across the whole range, a single step of the sweep’s own decade grid.

That is the invariance from the other side. The absolute band moves with the norm because n‖A‖u does, and what stays fixed is the band measured against the eigenvalues it is separating — so the routine has no constant to be wrong about, and no scaling of the input can put it in a regime it was not tested in. A tolerance would have had to be re-chosen at every one of those five scalings; this has nothing to re-choose.

Stated in relative terms that is not merely a claim about one scaling — it is an invariance, and the sweep says so. The relative gap at which half of a hundred shifts miscount:

   scale of the matrix     10⁻³      10⁰      10³      10⁶      10⁹
   relative gap          3.16·10⁻¹⁵  3.16·10⁻¹⁵  3.16·10⁻¹⁵  3.16·10⁻¹⁵  3.16·10⁻¹⁵

Identical across twelve orders of magnitude in the matrix, to the resolution of the sweep. There is no absolute length anywhere in the routine and the measurement is what says so, rather than the absence of a constant in the source.

It also pins the number the analysis leaves open. n·u for this matrix is 8.9·10⁻¹⁶ and the measured boundary is 3.16·10⁻¹⁵, so the band is about 3.5·n·u. That factor is worth having: a caller deciding whether two eigenvalues are far enough apart to be counted between needs a constant and not only a shape, and 3.5 is the constant on this family.

Why a constant that stays put is unusual here #

Almost every constant this collection measures moves. The looseness of a bound moves with the problem; the ratio of a disagreement to κu moves with the vector; the fraction of a bound a perturbation attains moves with its direction. Sixteen essays into this site the default expectation for “a small multiple of n·u” is that the multiple is a distribution.

Here it is a number, and the reason is worth stating because it explains what kind of quantity is being measured. The band is not the size of an error; it is the size of the region in which two eigenvalues stop being two. Below it the count is wrong at every shift, above it at none, and the transition is a property of when a perturbation of order n‖A‖u can move an eigenvalue past its neighbour. There is no direction to be lucky in and no realisation to average over — the question is whether two numbers are further apart than the format’s spacing, and that has one answer.

Which is the same reason the count is an integer in the first place. A quantity with no continuum in it has no room for a distribution either, and the band inherits the discreteness of the thing it bounds. That is a rarer arrangement than it sounds: it is the only measurement in this collection where the constant in front of n·u came out the same on every matrix tried.

What that buys, in one line #

A caller with two eigenvalues at λ and λ + δ, on a matrix of norm ‖A‖ in a format of unit roundoff u, can decide before running anything whether the count between them is meaningful: it is, if δ ⁄ ‖A‖ > 4 n u, and it is not otherwise. No sweep, no tolerance, and no output to interpret — which is the whole reason to prefer a count to a spectrum when a count is what the question wants.

That is a different situation from the rank decision, where a threshold has to be chosen and the choice is the computation. Here the threshold is not chosen and cannot be moved; the routine simply cannot resolve two eigenvalues closer together than the format’s spacing at that magnitude, which is a statement about the format.

It is also why the routine is worth having in a higher precision when the resolution matters: doubling the significand moves the band by eight decades and the routine is otherwise unchanged. The precision slider is the knob this whole site is built around, and here it moves exactly one thing.

A bracket, rather than an estimate #

Because ν is monotone in σ, bisecting on it gives an interval guaranteed to contain the k-th smallest eigenvalue. The loop maintains ν(lo) < k ≤ ν(hi), which is an invariant about integers and therefore survives the arithmetic, and it stops when the interval is narrower than a requested width.

Measured on a fourteen-by-fourteen matrix, every one of the fourteen brackets contains the corresponding Jacobi eigenvalue and every one is narrower than 5.4·10⁻¹³ of the norm.

Nothing else on this site returns an eigenvalue with a proof attached. The interval essays return enclosures for a solve, and krawczyk.js is the only other file in the collection that proves anything. This is the second, and its proof is cheaper: an enclosure of a solve needs directed rounding throughout, and an enclosure of an eigenvalue needs a monotone integer.

The proof is conditional in the way the previous section describes — the bracket is guaranteed for a matrix within n‖A‖u of A — and that is exactly the guarantee a backward-stable computation gives about everything else, stated for once as a bracket rather than as a bound.

Eigenvalues in an interval, and nothing else #

Two counts and a subtraction give the number of eigenvalues in [a, b), which is a question that a method computing the whole spectrum can only answer by computing the whole spectrum.

That is what the routine is actually for in practice. A spectral transformation for a large sparse symmetric problem wants the eigenvalues near a shift, needs to know how many there are before it allocates, and needs to be sure it has found all of them afterwards — and a Lanczos run cannot tell it either thing. The essay on an eigenvalue that arrives twice is about exactly that failure: a method that returns a plausible spectrum with a multiplicity silently wrong.

The count is the check that catches it. Run Lanczos, get some Ritz values, then count how many eigenvalues the matrix actually has in the interval those Ritz values span. If the two numbers differ, something was missed, and the cost of finding out is one factorisation.

The cost, honestly counted #

An eigenvalue by bisection needs one factorisation per bisection step, and the figure’s brackets took about fifty steps each to reach 10⁻¹³ of the norm. Fifty dense eliminations of a fourteen-by-fourteen matrix is more arithmetic than a Jacobi sweep that returns the whole spectrum, so on a small dense matrix the method is a curiosity.

Where it pays is where the matrix is tridiagonal or sparse. A tridiagonal LDLᵀ is a recurrence of length n — three flops per entry, no fill, no storage beyond a scalar — so a count costs O(n) and fifty of them cost 50n. Against that, computing the whole spectrum of a tridiagonal matrix costs O(n²) and gives eigenvalues nobody asked for. Ask for five eigenvalues out of ten thousand and bisection is three orders cheaper.

That is why the method survives in libraries under names that mention bisection: it is the routine that answers some of the eigenvalues, in this interval, with a bracket, and it is the only one that answers all three parts of that question at once. The dense figures here are drawing a small case of something whose economics are entirely about the sparse one.

And the inertia of a saddle-point matrix, for nothing #

The constraint field’s first essay proves that In(K) = (n, m, 0) from a congruence and measures it from a spectrum. This routine counts it from a single unpivoted elimination — the third route, and the only one a sparse code would run.

Measured across three conditionings and three shapes: ten positive and four negative pivots for n = 10, m = 4, and so on, agreeing with the Jacobi spectrum in every case.

The elimination it needs is exactly the one the constraint field’s regularised matrix admits under any ordering, so on that family the inertia is a by-product of a factorisation that was going to be computed anyway. A code that solves a regularised saddle-point system and finds n + 1 positive pivots has found a rank-deficient constraint, for free, on every solve.

The library’s other refusals guard the two things a reader is most likely to over-read. One is fed the claim that bisection returns a number rather than an interval, and required to refuse: the interval is the whole point. The other is fed the claim that an unpivoted elimination of K finds no negative pivot, and required to refuse that too — it finds exactly m.

Where else on this site an answer is an integer #

It is worth collecting them, because there are exactly three and they behave alike.

A count of negative pivots, which is this page. A rank, which is not one — rank is a decision about a gap, and the integer it returns is a function of a tolerance somebody chose, so it can be made to be anything. And the number of iterations a method took, which is an integer that is exactly reproducible and reports on the arithmetic rather than on the matrix.

Only the first is an integer that the matrix determines and that the arithmetic can only get wrong by a whole unit. That is a rarer thing than it sounds, and it is the reason this routine is used as a verification tool for results computed by other means rather than as a way of computing a spectrum.

The same distinction explains why the count of a tensor’s typical rank works the way it does in the tensor field: the classification there is the sign of a polynomial in the entries, computed by four multiplications and a subtraction, and it is right or it is wrong by a whole rank. Two fields, one shape, and in both of them the fragile case is a discriminant near zero.

What it cannot do #

Three limits, and all three are why this is a specialised tool rather than a replacement.

It gives no eigenvectors. The factorisation is discarded; what survives is a sign pattern. Anything that needs a direction needs a different method.

It needs a symmetric matrix. Congruence preserves inertia for symmetric matrices, and for an unsymmetric one the eigenvalues are not real and there is nothing to count on either side of.

And an unpivoted elimination can break down — a pivot exactly zero — which the routine handles by replacing it with an infinitesimal of the right sign. That is the standard device and it is a genuine perturbation of the matrix, so it is one more contribution to the band the fourth section measures rather than a free repair.