The zero that is not a missing entry

Every matrix this collection has had trouble with so far has been difficult for a reason the arithmetic supplied. A Hilbert matrix arrives ill conditioned and an exact rational inverse says by how much. A Wilkinson matrix grows under elimination until the bound that is never attained is attained after all. A kernel matrix has no zeros at all and is affordable anyway. In every case the algebra was willing and the numbers were not, and the repair was arithmetic: pivot, rescale, carry more bits.

This field opens on the first matrix here whose difficulty is algebraic, and the difference shows up as an absence rather than as a large number.

Where it comes from, which is one line of calculus #

Minimise ½xᵀHx − fᵀx subject to Ax = g, with H symmetric positive definite and A of full row rank. Introduce a multiplier y for each constraint and the Lagrangian is that objective plus yᵀ(Ax − g). Setting both gradients to zero gives

Hx + Aᵀy = f Ax = g

which is one linear system in the stacked unknown, with the matrix

K = [ H Aᵀ ] [ A 0 ]

The zero block is the second derivative of the Lagrangian with respect to the multipliers, and a Lagrangian is linear in them. So the zero is a theorem. That is the whole difference between this field and the sparsity one: spmat.js stores a zero because a particular problem happened not to couple two particular unknowns, and a different problem of the same kind would not have it. Here every problem of this kind has it, in the same place, at every size.

The row is the thing to watch rather than the pivot, and it is what a single frame cannot establish.

Why it is called a saddle point #

The name is not decoration and it says what the matrix is. The Lagrangian is a function of both x and y, and at the solution it is minimised in x and maximised in y: moving along a feasible direction raises the objective, and moving the multipliers raises the penalty for infeasibility. A stationary point that is a minimum in some directions and a maximum in others is a saddle, and the second derivative at a saddle is an indefinite matrix. So the indefiniteness is not a property of a particular H and A; it is the shape of the stationary point the problem asks for.

That also says what would have to change for the matrix to be definite: the problem would have to stop having constraints. Every technique in this field is therefore a way of living with indefiniteness rather than removing it — which is why the two rebuilds in the next essay both work by eliminating one block, and why the two preconditioners two essays on both keep the zero where it is.

Where these matrices come from, which is nearly everywhere #

Four families produce this exact block structure, and they have nothing else in common.

Constrained least squares. Fit a curve required to pass through a point, or a set of coefficients required to sum to one. The (1, 1) block is AᵀA, the constraint is the requirement, and the essay on the weighted version shows what happens when the requirement is imposed by a penalty instead.

Incompressible flow. Discretise a velocity field and a pressure, require the velocity to be divergence-free, and the pressure is the multiplier: the (2, 2) block is zero because pressure does not appear in the momentum equation except through its gradient.

Network and structural equilibrium. A flow that must balance at every node, a truss whose displacements must satisfy a support condition. The constraint rows are the topology and the multipliers are the node potentials or the reaction forces.

Every step of an interior-point method for a constrained optimisation problem, which is what the fifth essay in this field is about — and where the matrix is deliberately driven towards singularity, iteration after iteration, as the algorithm works.

The list matters because it says the shape is not an artefact of writing the Lagrangian down. It is what a constraint looks like once it has been differentiated, whatever it was a constraint on.

The count is settled before anything runs #

A congruence X ↦ SᵀXS with S nonsingular does not preserve eigenvalues, and does preserve their signs. That is Sylvester’s law of inertia, and it is the most useful theorem in this field because the block elimination

[ I 0 ] [ H Aᵀ ] [ I −H⁻¹Aᵀ ] [ H 0 ] [ −AH⁻¹ I ] [ A 0 ] [ 0 I ] = [ 0 −AH⁻¹Aᵀ ]

is a congruence, written out. H is positive definite by hypothesis. AH⁻¹Aᵀ is positive definite whenever A has full row rank, because vᵀAH⁻¹Aᵀv = ‖H⁻¹ᐟ²Aᵀv‖² and Aᵀv is nonzero for v ≠ 0. So the inertia of K is

(n positive, m negative, 0 zero)

exactly, for every H, every A, every conditioning and every size. The hero figure is that statement measured: n marks above zero and m below, and the counts printed beside them.

This is stronger than it sounds. Every other spectral claim on this site is a computation whose result has a tolerance attached — an eigenvalue one vector cannot see, a rank that is a decision about a gap, a condition number that is only as meaningful as the units it was measured in. The inertia of a saddle-point matrix is a pair of integers that is known before the matrix exists.

Cholesky does not fail somewhere #

The first thing anybody does with a symmetric matrix is try a Cholesky, because it costs half what an LU costs and because it fails if and only if the matrix is not positive definite, which makes the failure a useful test rather than an inconvenience.

Run it on K and it does not fail at an unpredictable place. The leading n × n block is H, which is positive definite, so the first n pivots are positive at any conditioning whatever. The (n + 1)-th pivot is the first diagonal entry of the Schur complement the elimination has been building, and that entry is

0 − a₁ᵀH⁻¹a₁ = −(AH⁻¹Aᵀ)₁₁ < 0

where a₁ᵀ is the first row of A. So the failure has an address — row n + 1 — and a value, and both were available before the factorisation started. The figure above prints both and the gap between them: the pivot the routine reports and the entry of AH⁻¹Aᵀ computed independently, agreeing to every bit at κ(H) = 1 and to nine digits at κ(H) = 10¹⁰, which is the arithmetic showing through a claim that is otherwise exact.

How the agreement degrades is worth a number, because the obvious guess is wrong in the useful direction:

   κ(H)      relative gap    digits agreeing
   1          0                  16+
   10²        5.88·10⁻¹⁶         15.2
   10⁴        1.41·10⁻¹⁴         13.9
   10⁶        1.34·10⁻¹²         11.9
   10⁸        3.73·10⁻¹¹         10.4
   10¹⁰       6.70·10⁻¹⁰          9.2
   10¹²       1.04·10⁻⁶           6.0

A κ(H)·u argument would say one digit lost per decade of conditioning. The fitted exponent is 0.78 — three quarters of a digit per decade — so at κ(H) = 10¹⁰ the check still carries nine digits where the linear estimate allows six.

The reason is that neither quantity is computed by an operation whose error is κ(H) times a rounding. The Cholesky pivot accumulates the Schur complement through n steps of elimination on a positive definite block, which is growth factor exactly one; the independent route solves with H once. Both inherit some of κ(H) and neither inherits all of it, and the exponent is where the two land.

The last row is outside that law by two orders, and the break is informative. Two decades past κ(H) = 10¹⁰ the gap jumps by more than a hundredfold, which is H’s own factorisation beginning to lose the Schur complement rather than the comparison losing accuracy — so the sub-linear regime has an end, it is at about κ(H) = 10¹¹, and past it the two routes are no longer two routes to one number.

What that check is for, and why it is not the usual one #

It is worth separating this comparison from the site’s ordinary two-routes habit, because it is doing a different job and the difference is the essay’s subject.

Everywhere else on this site, two routes to one number are run because neither is known to be right and agreement is evidence. Here the algebra is not in doubt: the (n+1)-th pivot is −(AH⁻¹Aᵀ)₁₁, exactly, by the congruence written out above, and no measurement could make that more or less true. What the comparison measures is not whether the claim holds but how much of it survives the arithmetic — which is the same question the whole site asks about every other identity and which is almost never asked about this one, because a claim that is exact invites being taken as exact.

That reframes what the table is evidence for. Nine digits at κ(H) = 10¹⁰ is not evidence that the inertia result holds; it is a measurement of the distance between an exact statement and a computed one, and the exponent 0.78 is the rate at which that distance opens.

The practical consequence is small and real. A routine that used the predicted pivot as a check on the factorisation — did the Cholesky stop where the algebra says it must — has a tolerance to choose, and the table says what it should be: about 10⁻¹⁶·κ(H)^0.78, which is nine digits at κ(H) = 10¹⁰ and six at 10¹². A check written with a fixed tolerance of 10⁻¹² would pass everywhere below κ(H) = 10⁷ and start firing spuriously above it, on a matrix where nothing whatever has gone wrong.

The refusal this page publishes is the reading that would otherwise survive: that a factorisation which fails is failing on the numbers, so a better pivot rule or more bits will get through. The assertion is fed the run at κ(H) = 1, where H is the identity and every entry of the matrix is of order one, and required to fail. It does. There is nothing wrong with the numbers.

Ten orders of magnitude of conditioning, and the factorisation stops in the same place every time. Row 11 of 14 at κ(H) = 1, 100, 10⁴, 10⁶, 10⁸ and 10¹⁰ alike. A numerical failure would move — a matrix ten orders harder would break something earlier, or would not break on the easy one at all. This one is a property of where the (2, 2) block starts, and no amount of conditioning puts it anywhere else.

And the pivot moves the wrong way. Its magnitude reads 0.082, 3.73, 328.6, 3.19·10⁴, 3.16·10⁶ and 3.14·10⁸ across those six stops, and the ratio per two decades of κ(H) is 45.5, 88.1, 97.1, 98.9 and 99.5 — converging on exactly a hundred. The failing pivot is asymptotically proportional to κ(H), so the harder the problem, the further from zero the offending number sits.

That rules out the reading the title is warning against. A pivot that failed because it was a small number lost in rounding would shrink towards zero as the problem worsened, and the repair would be a tolerance. This one reaches −3·10⁸, is negative at every stop, and equals −S₁₁ to every digit printed at every stop. There is nothing for a tolerance to catch, because nothing here is nearly anything.

What the repair is, and what it is not #

The essay on symmetric indefinite factorisation has the repair: allow a 2 × 2 pivot, take two variables at once, and a factorisation exists for every symmetric nonsingular matrix with a growth factor that is bounded. That essay’s headline matrix is [[0, 1], [1, 0]] — symmetric, perfectly conditioned, and with no diagonal entry to pivot on — which is the smallest saddle-point matrix there is: H is 1 × 1 and empty, A is [1].

What that repair does not do is make the matrix definite. It makes the factorisation exist. Everything downstream of definiteness stays gone: no Cholesky, no conjugate gradients, no guarantee that a diagonal preconditioner is positive, no interpretation of the solve as a descent. The whole field is what has to be rebuilt once that is accepted, and the next essay starts on it.

The spectrum has a closed form, and it is a bracket #

The inertia says how many eigenvalues are on each side. Rusten and Winther’s bounds say where, from four numbers: the extreme eigenvalues μₙ ≤ μ₁ of H and the extreme singular values σₘ ≤ σ₁ of A. The negative eigenvalues lie in

[ ½(μₙ − √(μₙ² + 4σ₁²)), ½(μ₁ − √(μ₁² + 4σₘ²)) ]

and the positive ones in [μₙ, ½(μ₁ + √(μ₁² + 4σ₁²))]. Two intervals, with a gap between them that contains zero.

That is exact ground truth for a spectrum, which this collection has had twice before — for the discrete Laplacian and for a circulant — and never in this form. Those two give the eigenvalues themselves; this gives an enclosure, so what is checked is containment, and the figure additionally reports how much room the bracket leaves over. On the hero’s matrix the extreme negative eigenvalue sits at 0.96 to 1.00 of its bound.

A bracket is what a Krylov method needs to be told. The rate the condition number predicts is a statement about a spectrum in one interval; a spectrum in two, with a hole in the middle, is a different Chebyshev problem and a different rate. Two essays further into this field, that bracket becomes a preconditioner’s design.

Three routes to the same two integers #

The inertia is worth having three ways, because three routes that agree are what this site means by knowing a number.

The algebra gives it in four lines, above: a congruence to a block diagonal whose blocks are definite of known sign. The spectrum gives it by counting the eigenvalues a Jacobi rotation sweep returns, which is what the hero figure draws and what costs the most — O(n³) with a large constant, and an answer whose accuracy is a tolerance. And the pivots give it from a single unpivoted elimination: the signs of the diagonal of D in an LDLᵀ, which is Sylvester’s law used as an algorithm and is the route a sparse code would actually take.

The three agree at every conditioning the library tests, and the third is the only one that scales. It is also the only one that returns an integer directly rather than by thresholding a float, which is the subject of its own essay in the spectra field: a count cannot be slightly wrong, so a computed inertia is exactly right or wrong by a whole eigenvalue, and where the second happens is measurable.

What the condition number does here #

κ(K) is not κ(H) and not κ(A), and it is not their product either. It is set by the two brackets: the largest eigenvalue in absolute value comes from the top of the positive interval and the smallest from wherever the two intervals come nearest to zero, which is governed by σₘ — the smallest singular value of the constraint matrix.

So a saddle-point system can be badly conditioned because the objective is, or because the constraints are nearly dependent, and the two enter differently. A nearly dependent constraint is the more dangerous of the two, because it takes the matrix towards singularity rather than towards inaccuracy: at exactly dependent constraints the inertia claim above fails outright, an eigenvalue arrives at zero, and one is taken from each side. The library’s own assertion is fed that case — a third constraint row made the sum of the first two — and required to reject the count.

The two things this matrix is not #

Two readings survive the argument above and are both wrong, and naming them is cheaper than having a reader carry them into the next four essays.

It is not a matrix that a better factorisation makes definite. Definiteness is a property of the quadratic form, and every congruence preserves it. There is no similarity, no scaling, no reordering and no change of basis that turns a matrix with m negative eigenvalues into one with none: the count is an invariant of the matrix and not of how it is written down. What a factorisation can do is exist — which is what the 2 × 2 pivot buys — and what a preconditioner can do is cluster the eigenvalues without moving them across zero, which is what the two essays on preconditioning measure.

And it is not ill conditioned by construction. The hero’s matrix has a condition number of about thirty. Nothing on this page is about a hard problem; it is about a matrix that a particular routine cannot factorise for a reason that has nothing to do with how hard the problem is. Conflating the two is what makes people reach for higher precision, and the refusal this page publishes exists to close that route: at κ(H) = 1, with H the identity, the factorisation stops at the same row.

The measurement, and what it is against #

Nothing on this page estimates its own error. Every finite double is a dyadic rational — m/2ᵏ with m an integer — so a matrix of doubles is a matrix of rationals with nothing approximated and nothing chosen, and exact.js now converts one. exactSolve then solves in BigInt with no rounding anywhere, and the forward error of a float solve is a difference from the answer to the problem the machine is holding rather than from a better float answer.

That extension is what lets this field measure anything at all. Until now the site could compare against a known answer only for problems whose entries were rationals before they were stored — the Hilbert matrix, and the closed forms the model problem supplies. The constrained systems here are not of that kind, and converting the stored double instead costs nothing and gives the same guarantee. Two essays later it is what makes the field’s central claim checkable.