Numerical Linear Algebra
where the algebra stops being the arithmetic
Eigenvalues, singular values, rank

An eigenvalue count that cannot be slightly wrong

Every spectral computation on this site returns floats with errors in them. Counting eigenvalues below a shift by the signs of an unpivoted elimination returns an integer, and an integer cannot be 6.9999999997 — so the answer is exactly right, or wrong by a whole eigenvalue, and where the second happens is a band of measurable width.
36 min read 31 figures Two routes to a numberMeasured, not assertedExact ground truth

Worth reading first: Symmetry is worth more than precision · A factorisation with nothing to pivot for · Rank is a decision.

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 number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 14 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.-4-2.35371-0.7074290.9388572.585144.231435.8777102468101214shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues14steps14the marks are a Jacobi decompositionand the staircase never saw one
Fig. 1 The count, drawn as a staircase, against a spectrum computed by a completely different method. Every step is at one of the marks and every mark has a step.

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.

ν(σ), the number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 6 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.-3-2.34571-1.69143-1.03714-0.3828570.2714290.9257140123456shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues6steps6the marks are a Jacobi decompositionand the staircase never saw one
Fig. 2 Six eigenvalues, where the staircase has six steps and the shape is easiest to read.

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.

The share of shifts inside a pair of eigenvalues that count it wrongly, against the pair's separationTwo eigenvalues at 1 and 1 + gap, with six others spread around them, and 100 shifts placed strictly between the pair — where the count must read 3. Down to a separation of 10⁻¹² every shift reads it correctly. At 10⁻¹³ one of 100 does not, at 10⁻¹⁴ 19 do not, and at 10⁻¹⁵ none of them reads it correctly. The boundary sits where n‖A‖u puts it — -4.405 for this matrix — because the floating-point count is the exact count of a matrix within that distance of A. This is the only place in the method where the answer can be wrong**, and it is wrong by a whole eigenvalue when it is: the failure is a miscount, not a small error.-15-13-11-9-7-5-300.250.50.751log₁₀ separation of the pairshare of shifts counted wronglyn‖A‖uthe only place it failswrong at 10⁻¹²0wrong at 10⁻¹⁴19wrong at 10⁻¹⁵100n‖A‖u-4.4wrong by a whole eigenvalueor not wrong at all
Fig. 3 The band, measured. The horizontal axis is how far apart two eigenvalues are; the vertical axis is the share of shifts between them that count them wrongly.

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 number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 8 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.-3-2.09771-1.19543-0.2931430.6091431.511432.41371012345678shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues8steps8the marks are a Jacobi decompositionand the staircase never saw one
Fig. 4 Eight eigenvalues, where the individual steps are widest and the construction is easiest to follow.

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.

ν(σ), the number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 12 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.00.51024681012shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues12steps12the marks are a Jacobi decompositionand the staircase never saw one
Fig. 5 The graded spectrum, whose eigenvalues span six decades and whose staircase is bunched at one end.

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 share of shifts inside a pair of eigenvalues that count it wrongly, against the pair's separationTwo eigenvalues at 1 and 1 + gap, with six others spread around them, and 100 shifts placed strictly between the pair — where the count must read 3. Down to a separation of 10⁻¹² every shift reads it correctly. At 10⁻¹³ one of 100 does not, at 10⁻¹⁴ 3 do not, and at 10⁻¹⁵ none of them reads it correctly. The boundary sits where n‖A‖u puts it — 4.859·10⁷ for this matrix — because the floating-point count is the exact count of a matrix within that distance of A. This is the only place in the method where the answer can be wrong**, and it is wrong by a whole eigenvalue when it is: the failure is a miscount, not a small error.-15-13-11-9-7-5-300.250.50.751log₁₀ separation of the pairshare of shifts counted wronglyn‖A‖uthe only place it failswrong at 10⁻¹²0wrong at 10⁻¹⁴3wrong at 10⁻¹⁵100n‖A‖u4.9·10⁷wrong by a whole eigenvalueor not wrong at all
Fig. 6 The same measurement on a matrix scaled by 10⁶, where the curve moves right by exactly six decades — because the band is n‖A‖u and not an absolute distance.

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.

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.

The share of shifts inside a pair of eigenvalues that count it wrongly, against the pair's separationTwo eigenvalues at 1 and 1 + gap, with six others spread around them, and 100 shifts placed strictly between the pair — where the count must read 3. Down to a separation of 10⁻¹² every shift reads it correctly. At 10⁻¹³ one of 100 does not, at 10⁻¹⁴ 0 do not, and at 10⁻¹⁵ none of them reads it correctly. The boundary sits where n‖A‖u puts it — 4807 for this matrix — because the floating-point count is the exact count of a matrix within that distance of A. This is the only place in the method where the answer can be wrong**, and it is wrong by a whole eigenvalue when it is: the failure is a miscount, not a small error.-15-13-11-9-7-5-300.250.50.751log₁₀ separation of the pairshare of shifts counted wronglyn‖A‖uthe only place it failswrong at 10⁻¹²0wrong at 10⁻¹⁴0wrong at 10⁻¹⁵100n‖A‖u4807wrong by a whole eigenvalueor not wrong at all
Fig. 7 An intermediate scaling, for reading the shift of the curve against the two extremes.

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.

The widths a proof is available at, against the precisionTwo boundaries against the number of significand bits, on a logarithmic vertical axis, with the band between them shaded. The narrow boundary runs from 0.00162 at 12 bits to 3.33·10⁻¹⁶ at 53, a small multiple of u at each. The wide boundary is 1.414214 at every precision — which is √2, the half-width at which the box first touches the line where the Jacobian is singular.9172533414910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹significand bitsbox half-widthwidest provablenarrowest provableutwo boundaries, two causesnarrowest, in units of u3widest, at 53 bits1.4√21.4one end is the arithmeticthe other is the function
Fig. 8 The other route to a proved answer on this site, from the arithmetic field.

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 six largest Ritz values, against a doubled eigenvalue at 10Distance from 10 for each of the six largest Ritz values, on a logarithmic vertical axis. The single-vector run of 16 steps has one value at 10 to 2.6·10⁻¹² and its second is 1.696 away — it has found the eigenvalue once. The block of two has two values at 10, to 1.1·10⁻⁹ and 4.7·10⁻⁹.12345610⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹Ritz value, largest firstdistance from 10one vectora block of twoa space, not a ratecopies found, one vector1copies found, block of two2Krylov dimension, one vector16the second copy is not in the spaceat any number of steps
Fig. 9 The failure the count catches, from the field’s own essay about it.

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.

What 12 members of a sequence cost, by what changes between themEvery bar is 12 solves of an 120×120 system, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n². When nothing changes the answer is reused and the cost is one factorisation. When only the right-hand side changes, one factorisation serves every member — 16.5 per cent of the independent cost, with every solve still at a relative residual of 1.54·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 4 factorisations for 12 members. When everything changes, nothing carries.nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried
Fig. 10 The shape of the argument in the field that counts operations, where a method’s value is what it declines to compute.

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.

The 14 eigenvalues of a saddle-point matrix with 10 unknowns and 4 constraints, inside their closed-form bracketsK = [[H, Aᵀ], [A, 0]] with κ(H) = 100 and κ(A) = 100. The shaded bands are the Rusten–Winther brackets, computed from four numbers — the extreme eigenvalues of H and the extreme singular values of A — before the matrix was assembled: negative eigenvalues in [-0.995, -10·10⁻⁵] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 4 below, which is Sylvester's law of inertia and not a property of this matrix: K is congruent to blkdiag(H, −AH⁻¹Aᵀ), both blocks are definite, and congruence preserves signs. The two groups are separated by a gap containing zero, at a ratio of 39.3 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it
Fig. 11 The count this routine reproduces, from the field where it is a theorem.

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.

A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 12 The integer that is a decision rather than a property, from this field’s own essay.

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.

Three perturbations of 10⁻¹⁴, three eigenvectors, one planeA circle with three radii at widely different angles, one for each perturbation, and a table of their residuals showing all three satisfy the eigenvalue equation.every radius is an eigenvectorseed 11: 2.0°seed 23: 60.5°seed 41: 18.2°none of them is wrongspread of the returned angles1spread of the returned planes1.3·10⁻¹⁵worst ‖Ax − λx‖5·10⁻¹⁵an exactly repeated eigenvalueno correct answer to fail to return
Fig. 13 What the routine declines to compute, and why that is sometimes a mercy.

At other settings

ν(σ), the number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 10 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.-3-1.84971-0.6994290.4508571.601142.751433.901710246810shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues10steps10the marks are a Jacobi decompositionand the staircase never saw one
Fig. 14 Ten eigenvalues, evenly spread.
ν(σ), the number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 20 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.-4-1.609710.7805713.170865.561147.9514310.34170369121518shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues20steps20the marks are a Jacobi decompositionand the staircase never saw one
Fig. 15 And twenty, where the steps are close enough that the marks are what makes them readable.
ν(σ), the number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 14 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.-4-2.35371-0.7074290.9388572.585144.231435.8777102468101214shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues14steps14the marks are a Jacobi decompositionand the staircase never saw one
Fig. 16 The hero’s spectrum drawn again with the scan’s own parameters.
The share of shifts inside a pair of eigenvalues that count it wrongly, against the pair's separationTwo eigenvalues at 1 and 1 + gap, with six others spread around them, and 100 shifts placed strictly between the pair — where the count must read 3. Down to a separation of 10⁻¹² every shift reads it correctly. At 10⁻¹³ one of 100 does not, at 10⁻¹⁴ 30 do not, and at 10⁻¹⁵ none of them reads it correctly. The boundary sits where n‖A‖u puts it — 4.859·10⁵ for this matrix — because the floating-point count is the exact count of a matrix within that distance of A. This is the only place in the method where the answer can be wrong**, and it is wrong by a whole eigenvalue when it is: the failure is a miscount, not a small error.-15-13-11-9-7-5-300.250.50.751log₁₀ separation of the pairshare of shifts counted wronglyn‖A‖uthe only place it failswrong at 10⁻¹²0wrong at 10⁻¹⁴30wrong at 10⁻¹⁵100n‖A‖u4.9·10⁵wrong by a whole eigenvalueor not wrong at all
Fig. 17 The band at a scaling of 10⁴.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 18 The decision this routine does not have to make, from the field’s own essay.
What a perturbation of 10⁻⁶ moves, against the gap it is applied acrossThree curves against the gap between two eigenvalues, both axes logarithmic. Two are flat; the third rises steeply as the gap closes and then levels off.10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻⁸10⁻⁵10⁻²gap between the two eigenvalueshow far it movedthe eigenvectorsthe eigenvaluestheir plane‖E‖ / gapone perturbation, three answerseigenvalue shift, spread over the sweep1plane angle, spread over the sweep1eigenvector angle, spread1.6·10⁵the dashed line is Davis–Kahan's ‖E‖/gaptwo of the three never noticed
Fig. 19 What a gap is worth to a method that computes values rather than counts.
How far a perturbation of size ε moves an eigenvalue, 8×8A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/8)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 20 Why a computed eigenvalue moves at all.
Relative error of each singular value of a bidiagonal matrix graded over 30 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 2.098·10⁻³⁰. One-sided Jacobi holds every value to 4.41·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 2.21·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 1 — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly2.1·10⁻³⁰worst, one-sided Jacobi4.4·10⁻¹⁶worst, zero-shift QR2.2·10⁻¹⁶worst, eigenvalues of BᵀB1a relative error is a ratioand the denominator is the answer
Fig. 21 And a case where an eigenvalue’s relative accuracy is available and a count would say nothing.
500 random symmetric orderings of a regularised saddle-point matrix, and what each factorisation reproducesA quasi-definite matrix — [[H + δI, Aᵀ], [A, −δI]] with δ = 10⁻⁶ — has an LDLᵀ factorisation with diagonal D for every symmetric permutation, so the ordering can be chosen for fill with no numerical veto at all. That is the count: 500 of 500 orderings factorise here, against 69.0 per cent for the same matrix with the zero block left where the problem put it. Each bar counts the orderings whose factorisation left a relative residual ‖PKPᵀ − LDLᵀ‖/‖K‖ in that decade, with the tallest holding 136 of the 500. They span 1.56·10⁻¹⁶ to 1.02·10⁻⁷ — 9 orders — and the growth factor across them runs from 1 to 6.43·10⁵ — which is 0.643/δ, an inverse law that holds at every δ this figure is drawn at. Existence is not stability, and the theorem says only the first.-18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.69worst growth6.4·10⁵growth × δ0.64the ordering is free to chooseand not free of consequence
Fig. 22 The factorisation this routine is a by-product of, on the family where it always exists.
The block diagonal of a Bunch–Kaufman factorisation of an 8×8 saddle-point matrix at τ = 10⁻⁶The matrix D from PAPᵀ = LDLᵀ, drawn as a matrix. 4 of its blocks are 2×2 — the shaded pairs — and the rest are single entries. A 2×2 block is taken exactly where no single diagonal entry is large enough to divide by safely, which on a matrix whose lower-right block is zero by construction is most of it. The badge carries the residual of all three pivot rules.D from PAPᵀ = LDLᵀ — the shaded pairs are 2×2 pivots10⁻⁶0.749······0.749·········1.2·10⁻⁶1.4······1.4·········2.1·10⁻⁶0.549······0.549·········3.6·10⁻⁶0.614······0.614·three rules, one matrix‖PAPᵀ − LDLᵀ‖, blocks5.8·10⁻¹⁷‖PAPᵀ − LDLᵀ‖, diagonal3.1·10⁻¹¹growth, blocks1.3growth, diagonal5·10⁵the zero block is what the problem saysand one rule does not need it to be nonzero
Fig. 23 And the pivot rule it deliberately does without.
The spectrum of C⁻¹T against T's own, n = 64, ρ = 0.8Two sorted spectra on a linear vertical axis. The unpreconditioned Toeplitz matrix spreads its eigenvalues across its whole range; the preconditioned operator puts 97% of them within a tenth of 1, with 0 below zero.08162432404856640123456789index, sortedeigenvalueT aloneC⁻¹Teigenvalues, not singular valuesfraction within 0.1 of 10.97eigenvalues below zero0λ_min of the preconditioner0.11the dashed line is 1, where the cluster formsthe solid line is zero
Fig. 24 Clusters, which are what the band is about resolving.
Both forms of an interior-point step: two condition numbers that climb together and two errors that do notOne Newton step of an interior-point method for a quadratic programme with 8 unknowns and 6 constraints, 3 of them active, written twice. The augmented form is [[H, Cᵀ], [C, −D⁻¹]] and the condensed form is H + CᵀDC, and they have the same solution. As the barrier parameter μ falls, the diagonal D separates — z/s runs to 1/μ on the active constraints and to μ on the inactive ones — so both condition numbers climb: 29.04 to 3.044·10¹⁵ for the augmented form and 155.8 to 2.403·10¹⁶ for the condensed one, within a factor of 7.9 of each other. The two forward errors, both measured against the exact rational solution of the system the machine is holding, do not follow: the condensed form's rises to 0.3098 and the augmented form's stays at 1.035·10⁻¹⁵ — fifteen correct digits at a condition number of 3.04·10¹⁵.-14-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μ — the barrier parametercondition number, and relative errorκ₂, condensedκ₂, augmentederror, condensederror, augmentedagainst a BigInt answerκ₂ augmented, μ = 10⁻¹⁴3·10¹⁵its relative error10⁻¹⁵κ₂ condensed2.4·10¹⁶its relative error0.31the same step, written two waysand only one of them is solvable
Fig. 25 A matrix whose inertia is worth checking on every solve, from the constraint field.
ν(σ), the number of eigenvalues below σ, counted from the signs of an unpivoted LDLᵀThe staircase is the number of negative pivots in an LDLᵀ factorisation of A − σI, taken with no pivoting at all. Congruence preserves the signs, so that count is the number of eigenvalues below σ — Sylvester's law of inertia, used as an algorithm. The vertical marks are the 8 eigenvalues a Jacobi decomposition returns, which is a completely different computation; every step of the staircase is at one of them and every one of them has a step. Over 2000 shifts placed at random across the interval the two answers disagree 0 times. That is not a statement about a tolerance: the output is a count, so it is exactly right or wrong by a whole eigenvalue, and there is nothing in between for a rounding error to land in.00.51012345678shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues8steps8the marks are a Jacobi decompositionand the staircase never saw one
Fig. 26 A graded spectrum with eight eigenvalues, where the steps bunch at one end.
The share of shifts inside a pair of eigenvalues that count it wrongly, against the pair's separationTwo eigenvalues at 1 and 1 + gap, with six others spread around them, and 200 shifts placed strictly between the pair — where the count must read 3. Down to a separation of 10⁻¹² every shift reads it correctly. At 10⁻¹³ one of 200 does not, at 10⁻¹⁴ 39 do not, and at 10⁻¹⁵ none of them reads it correctly. The boundary sits where n‖A‖u puts it — -4.405 for this matrix — because the floating-point count is the exact count of a matrix within that distance of A. This is the only place in the method where the answer can be wrong**, and it is wrong by a whole eigenvalue when it is: the failure is a miscount, not a small error.-15-13-11-9-7-5-300.250.50.751log₁₀ separation of the pairshare of shifts counted wronglyn‖A‖uthe only place it failswrong at 10⁻¹²0wrong at 10⁻¹⁴39wrong at 10⁻¹⁵200n‖A‖u-4.4wrong by a whole eigenvalueor not wrong at all
Fig. 27 Twice as many shifts a point, where the transition is measured more finely.
The 16 eigenvalues of a saddle-point matrix with 10 unknowns and 6 constraints, inside their closed-form bracketsK = [[H, Aᵀ], [A, 0]] with κ(H) = 100 and κ(A) = 100. The shaded bands are the Rusten–Winther brackets, computed from four numbers — the extreme eigenvalues of H and the extreme singular values of A — before the matrix was assembled: negative eigenvalues in [-0.995, -10·10⁻⁵] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 6 below, which is Sylvester's law of inertia and not a property of this matrix: K is congruent to blkdiag(H, −AH⁻¹Aᵀ), both blocks are definite, and congruence preserves signs. The two groups are separated by a gap containing zero, at a ratio of 19.4 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue6 negative10 positivecounted before it was formedpositive10negative6at zero0innermost ratio19the zero block is a theoremand so is the count either side of it
Fig. 28 Six negative eigenvalues to count, from the constraint field.
500 random symmetric orderings of a regularised saddle-point matrix, and what each factorisation reproducesA quasi-definite matrix — [[H + δI, Aᵀ], [A, −δI]] with δ = 10⁻⁸ — has an LDLᵀ factorisation with diagonal D for every symmetric permutation, so the ordering can be chosen for fill with no numerical veto at all. That is the count: 500 of 500 orderings factorise here, against 69.0 per cent for the same matrix with the zero block left where the problem put it. Each bar counts the orderings whose factorisation left a relative residual ‖PKPᵀ − LDLᵀ‖/‖K‖ in that decade, with the tallest holding 130 of the 500. They span 1.65·10⁻¹⁶ to 6.17·10⁻⁴ — 13 orders — and the growth factor across them runs from 1 to 6.43·10⁷ — which is 0.643/δ, an inverse law that holds at every δ this figure is drawn at. Existence is not stability, and the theorem says only the first.-18-15-12-9-6-300log₁₀ ‖PKPᵀ − LDLᵀ‖ / ‖K‖share of orderingsbestworstexistence and stabilityfactorise1unregularised0.69worst growth6.4·10⁷growth × δ0.64the ordering is free to chooseand not free of consequence
Fig. 29 And the factorisations whose signs would count them.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 3 constraints with κ(H) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 12.99 to 3.981·10¹⁰ — the square of κ(A), for the reason the normal equations square it. The null-space method solves with the reduced Hessian ZᵀHZ, whose condition number is 21.13 at the start of the sweep and 21.13 at the end: it does not contain κ(A) at all. The two forward errors, measured against a solution computed in BigInt rationals, follow their own condition numbers: 5.314·10⁻⁶ against 5.788·10⁻¹² at the far end. Both methods are correct and one of them is usable.01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares
Fig. 30 A constrained solve whose inertia this routine would check.
What the regularisation costs, and what 6 steps of refinement take backSolving [[H + δI, Aᵀ], [A, −δI]] instead of K gives the exact answer to a different problem, so its error is proportional to δ: measured at 1489·δ across six decades, which is a slope of one and not a trend. Refining against the unregularised matrix — the residual formed with K and the correction solved with the regularised factorisation — removes that term entirely, because the perturbation was never in the residual. It works while δ is below σ_min(K) = 6.797·10⁻⁴, marked on the axis, and stops working above it: the iteration's contraction factor is δ/σ_min and a fixed point needs that under one. So the trade-off curve every regularisation essay on this site has drawn — a term falling in δ against a term rising in it — has, here, a left branch that can simply be removed.-14-12-10-8-6-4-210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ δrelative error against the exact answerδ = σ_min(K)no refinement6 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶3.6·10⁻¹⁶refined at δ = 10⁻²0.62the perturbation is known exactlybecause the code chose it
Fig. 31 And the regularisation that makes the counting factorisation always exist.

Where this goes next

-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it The zero that is not a missing entry A constrained minimisation produces a matrix with a zero block, and the zero is a theorem rather than a sparsity pattern. No pivot order makes it positive definite, no precision changes that, and Cholesky does not fail somewhere on it — it fails at the first constraint row, on a number the problem already contained. 1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap. 11.522.511.52xyexactly one roota verdict, not a bound‖I − C F′(X)‖0.28width of X0.8width of K(X)0.23strictly inside is a proofand overlapping is nothing at all Proving the answer is in the box Every other method here computes a number and estimates how wrong it is. This one returns a verdict: there is exactly one solution in this box, or there is none, or — the honest third outcome — nothing can be said. Two of the three are proofs about infinitely many points from finitely many operations.

What links here

Computed from the collection, not written here: the essays that point at this one.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

BisectionCondition numberEigenvalue bracketExact ground truthInertiaLDLᵀ factorisationNumerical rankSaddle-point systemsSpectral slicingSymmetric indefinite