Numerical Linear Algebra
where the algebra stops being the arithmetic
The matrix a constraint makes

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.
39 min read 31 figures Where the algebra stops being the arithmeticExact ground truthMeasured, not asserted

Worth reading first: A factorisation with nothing to pivot for · The condition number is an amplifier · Elimination is a sequence of choices.

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.

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. 1 The spectrum of one such matrix, inside two brackets that were computed from four numbers before the matrix was assembled. The count either side of zero was known before that.

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 pivots Cholesky computes on a saddle-point matrix, and the one it stops at — row 11 of 14Every pivot Cholesky forms on K = [[H, Aᵀ], [A, 0]], in order, with the 10 taken inside H drawn as positive and the 11-th, which is not, drawn as the failure. It is not a near miss and it is not a rounding accident: the leading n × n block is positive definite, so the first 10 pivots are positive at any conditioning, and the next one is 0 − a₁ᵀH⁻¹a₁ where a₁ᵀ is the first row of A. That equals −(AH⁻¹Aᵀ)₁₁ = -328.56, and the pivot the routine actually reports is -328.56 — the same number to 13 digits. The failure has an address and a value, and both were available before the factorisation started.02468101210⁻⁵10⁻³10⁻¹10¹10³pivot, in order|pivot|negative, at row 11|S₁₁| from AH⁻¹Aᵀa failure with an addressrows before it10the pivot-329−(AH⁻¹Aᵀ)₁₁-329disagreement8.8·10⁻¹⁴it does not fail somewhereit fails at the constraint
Fig. 2 What that costs the first factorisation anybody would reach for.

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.

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) = 1. 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, -0.618] 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 0.04 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 ratio0.04the zero block is a theoremand so is the count either side of it
Fig. 3 With A orthogonal, so that σₘ = σ₁ = 1 and the negative bracket collapses to a point: the cleanest case the field has.

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.

Fill under three orderings of a sparse saddle-point matrix, predicted symbolically and then countedA stiffness block of 24 unknowns with 4 constraints each touching 3 of them. Each panel is the pattern of PKPᵀ under one ordering, with the entries the factorisation adds shown against the entries the matrix arrived with. The counts are natural 113, minimum-degree 63, reverse-cuthill-mckee 63: computed from the graph before any number was touched, and matched exactly by the factorisation of the regularised matrix, because quasi-definiteness means nothing during the numeric phase is allowed to move a pivot. That is the guarantee sparse LU cannot make — there the pivot order depends on the values, so the symbolic phase can only bound the fill and the allocation has to be able to grow.natural113 predicted · 113 countedminimum-degree63 predicted · 63 countedreverse Cuthill–McKee63 predicted · 63 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural113minimum-degree63reverse-cuthill-mckee63predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 4 One of those families as a sparsity pattern: a stiffness block with a few constraint rows across it, from the essay on what an ordering can be chosen for.

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.

ν(σ), 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.-8-5.30857-2.617140.07428572.765715.457148.14857024681012shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues12steps12the marks are a Jacobi decompositionand the staircase never saw one
Fig. 5 The same count, taken from the signs of an elimination rather than from a spectrum — the method the essay on that has to itself.

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 fifteen digits at κ(H) = 1 and to eight at κ(H) = 10¹⁰, which is the arithmetic showing through a claim that is otherwise exact.

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.

The pivots Cholesky computes on a saddle-point matrix, and the one it stops at — row 11 of 14Every pivot Cholesky forms on K = [[H, Aᵀ], [A, 0]], in order, with the 10 taken inside H drawn as positive and the 11-th, which is not, drawn as the failure. It is not a near miss and it is not a rounding accident: the leading n × n block is positive definite, so the first 10 pivots are positive at any conditioning, and the next one is 0 − a₁ᵀH⁻¹a₁ where a₁ᵀ is the first row of A. That equals −(AH⁻¹Aᵀ)₁₁ = -0.082076, and the pivot the routine actually reports is -0.082076 — the same number to 16 digits. The failure has an address and a value, and both were available before the factorisation started.02468101210⁻³10⁻²10⁻¹110¹10²pivot, in order|pivot|negative, at row 11|S₁₁| from AH⁻¹Aᵀa failure with an addressrows before it10the pivot-0.082−(AH⁻¹Aᵀ)₁₁-0.082disagreement0it does not fail somewhereit fails at the constraint
Fig. 6 The same figure with H the identity, where every entry of the matrix is of order one and the factorisation stops at exactly the same row.

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.

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. 7 The two rebuilds, priced. The next essay is about the gap between the two lower lines.

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.

The 11 eigenvalues of a saddle-point matrix with 10 unknowns and 1 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, -0.618] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 1 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 0.018 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue1 negative10 positivecounted before it was formedpositive10negative1at zero0innermost ratio0.018the zero block is a theoremand so is the count either side of it
Fig. 8 One constraint, so one eigenvalue below zero. The bracket around it is narrow because A has a single singular value and σₘ = σ₁.
The 18 eigenvalues of a saddle-point matrix with 10 unknowns and 8 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 8 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 66.4 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue8 negative10 positivecounted before it was formedpositive10negative8at zero0innermost ratio66the zero block is a theoremand so is the count either side of it
Fig. 9 Eight, so eight — and the negative bracket widens as A’s singular values spread.

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.

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. 10 Where a counted inertia can be wrong, which is the only place it can be: inside a cluster narrower than n‖A‖u.

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.

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) = 1. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 1 to 10·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 1 at the start of the sweep and 1 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: 7.341·10⁻⁷ against 2.025·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⁵10·10⁹κ(ZᵀHZ), all stops1range-space forward error7.3·10⁻⁷null-space forward error2·10⁻¹²both are the same algebraand only one squares
Fig. 11 With H the identity, so that the only conditioning in the picture belongs to the constraint.

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.

What the regularisation costs, and what 0 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 refinement0 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶0.0014refined at δ = 10⁻²0.93the perturbation is known exactlybecause the code chose it
Fig. 12 The one thing that does move the count, from the sixth essay in this field — and what it costs to move it.

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.

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. 13 Where that leads, four essays on: two ways of writing one system, with the same condition number and errors fourteen orders apart, both measured against a BigInt answer.

At other settings

The 12 eigenvalues of a saddle-point matrix with 10 unknowns and 2 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 2 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 4.48 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue2 negative10 positivecounted before it was formedpositive10negative2at zero0innermost ratio4.5the zero block is a theoremand so is the count either side of it
Fig. 14 Two constraints. The gap containing zero is at its widest here, which is the easiest case for anything iterative.
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. 15 Six, for reading against the hero’s four.
The pivots Cholesky computes on a saddle-point matrix, and the one it stops at — row 11 of 14Every pivot Cholesky forms on K = [[H, Aᵀ], [A, 0]], in order, with the 10 taken inside H drawn as positive and the 11-th, which is not, drawn as the failure. It is not a near miss and it is not a rounding accident: the leading n × n block is positive definite, so the first 10 pivots are positive at any conditioning, and the next one is 0 − a₁ᵀH⁻¹a₁ where a₁ᵀ is the first row of A. That equals −(AH⁻¹Aᵀ)₁₁ = -3.1551·10⁶, and the pivot the routine actually reports is -3.1551·10⁶ — the same number to 9 digits. The failure has an address and a value, and both were available before the factorisation started.02468101210⁻⁸10⁻⁵10⁻²10¹10⁴10⁷pivot, in order|pivot|negative, at row 11|S₁₁| from AH⁻¹Aᵀa failure with an addressrows before it10the pivot-3.2·10⁶−(AH⁻¹Aᵀ)₁₁-3.2·10⁶disagreement2.5·10⁻¹⁰it does not fail somewhereit fails at the constraint
Fig. 16 At κ(H) = 10⁸, where the pivots inside H spread over eight decades and the row the failure lands on has still not moved.
The pivots Cholesky computes on a saddle-point matrix, and the one it stops at — row 11 of 14Every pivot Cholesky forms on K = [[H, Aᵀ], [A, 0]], in order, with the 10 taken inside H drawn as positive and the 11-th, which is not, drawn as the failure. It is not a near miss and it is not a rounding accident: the leading n × n block is positive definite, so the first 10 pivots are positive at any conditioning, and the next one is 0 − a₁ᵀH⁻¹a₁ where a₁ᵀ is the first row of A. That equals −(AH⁻¹Aᵀ)₁₁ = -3.189·10⁴, and the pivot the routine actually reports is -3.189·10⁴ — the same number to 11 digits. The failure has an address and a value, and both were available before the factorisation started.02468101210⁻⁷10⁻⁵10⁻³10⁻¹10¹10³10⁵pivot, in order|pivot|negative, at row 11|S₁₁| from AH⁻¹Aᵀa failure with an addressrows before it10the pivot-3.2·10⁴−(AH⁻¹Aᵀ)₁₁-3.2·10⁴disagreement3.4·10⁻¹²it does not fail somewhereit fails at the constraint
Fig. 17 And at 10⁶, for reading the two against each other.
The spectrum of P⁻¹K with S = AH⁻¹Aᵀ, exactly, at 10 unknowns and 4 constraintsP = blkdiag(H, S). With S the exact Schur complement AH⁻¹Aᵀ the preconditioned matrix has exactly three distinct eigenvalues — 1 with multiplicity n − m = 6, and (1 ± √5)/2 with multiplicity 4 each. Those are 1 − φ = -0.618034 and φ = 1.61803, the golden ratio, which arrives from λ² − λ − 1 = 0 rather than from anything anybody chose. The dashed lines are that closed form and the marks are the computed spectrum; here they are 3 distinct values and the largest distance from the closed form anywhere is 2.909·10⁻¹⁴. A minimal polynomial of degree three means a Krylov method finishes in three steps, which is what the next figure measures.00.3670080.7340171.101031.468031.835040eigenvalue of P⁻¹Kwritten down, then computeddistinct3at 16φ computed1.6off the closed form2.9·10⁻¹⁴1 − φ1φthe preconditioner's effect is a theoremand the golden ratio is in it
Fig. 18 What a preconditioner does to the two intervals, from two essays on: three points and nothing between them.
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. 19 The repair, from the elimination field: a 2 × 2 pivot taken where no 1 × 1 pivot exists.
How often Cholesky still calls a 12×12 matrix positive definite, against λ_min/λ_max in units of the format's own roundoffThree curves, one per precision, of the share of 24 seeded matrices on which the factorisation succeeds. Measured in units of each format's unit roundoff the three lie almost on top of one another, with the edge — the smallest ratio at which every seed succeeds — at 1.0, 1.0, 1.8 times u. The absolute thresholds are 6·10⁻⁸, 2.3·10⁻¹⁰ and 2·10⁻¹⁶: nine orders of magnitude apart, and the same number in the format's own units.10⁻³10⁻²10⁻¹110¹10²10³10⁴00.250.50.751λ_min / λ_max, in units of the format's own ushare of seeds that succeed24 bits32 bits53 bitsone threshold, three formats24-bit edge, absolute6·10⁻⁸32-bit edge, absolute2.3·10⁻¹⁰53-bit edge, absolute2·10⁻¹⁶in units of u, at 53 bits1.8a yes-or-no question with a precision in itand a coin flip three decades below the edge
Fig. 20 And what Cholesky’s failure means when the matrix is nearly definite rather than definitely not, which is the case this page does not cover.
Growth factor under partial pivoting to n = 40: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 21 The elimination field’s growth factor, for contrast: a difficulty that is entirely arithmetic and that a pivot rule does repair.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 22 Exact ground truth in its older form — a matrix whose entries were rationals before they were stored.
Three bases for the same null space, at 10 unknowns and 4 constraintsThe same constrained problem solved three times, differing only in which basis Z is used for the null space of A. The orthonormal basis, from a QR of Aᵀ, has κ(Z) = 1 exactly and is dense — 100 per cent of its entries are nonzero. The fundamental basis built on the first 4 columns has κ(Z) = 1.994·10⁸, and its reduced Hessian comes out at 3.801·10¹⁶, which is κ(Z)² to within a factor of 0.956 — the square attained rather than bounded. Its answer is wrong by 0.05179, against 1.07·10⁻¹⁵ for the orthonormal one. Choosing the same kind of basis by pivoting instead gives κ(Z) = 2.06, an error of 6.71·10⁻¹⁶, and the same 50 per cent density: all of the sparsity and none of the loss.κ(A) = 10⁶ throughout · κ(H) = 100 · the answer is the same answer for every basisorthonormal — κ(Z)1κ(ZᵀHZ)25.6relative error1.07·10⁻¹⁵first m basic — κ(Z)1.99·10⁸κ(ZᵀHZ)3.8·10¹⁶relative error0.0518pivoted basic — κ(Z)2.06κ(ZᵀHZ)31.9relative error6.71·10⁻¹⁶what the choice costsdensity, orthonormal1density, fundamental0.5κ(ZᵀHZ) ÷ κ(Z)², naive0.96error, pivoted choice6.7·10⁻¹⁶every one of them is a basisand one of them loses fourteen digits
Fig. 23 The next essay’s picture: three ways of removing the constraint, and the condition number each one inherits.
ν(σ), 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. 24 Twenty eigenvalues counted by an integer, which is where the inertia of this page becomes an algorithm.
How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 25 And the identity the whole site runs on, which the fifth essay in this field takes to a place where it stops being useful.
The 13 eigenvalues of a saddle-point matrix with 10 unknowns and 3 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 3 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 15.2 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue3 negative10 positivecounted before it was formedpositive10negative3at zero0innermost ratio15the zero block is a theoremand so is the count either side of it
Fig. 26 Three constraints, for reading the bracket widths against the hero’s four.
The pivots Cholesky computes on a saddle-point matrix, and the one it stops at — row 15 of 19Every pivot Cholesky forms on K = [[H, Aᵀ], [A, 0]], in order, with the 14 taken inside H drawn as positive and the 15-th, which is not, drawn as the failure. It is not a near miss and it is not a rounding accident: the leading n × n block is positive definite, so the first 14 pivots are positive at any conditioning, and the next one is 0 − a₁ᵀH⁻¹a₁ where a₁ᵀ is the first row of A. That equals −(AH⁻¹Aᵀ)₁₁ = -128.47, and the pivot the routine actually reports is -128.47 — the same number to 13 digits. The failure has an address and a value, and both were available before the factorisation started.024681012141610⁻⁴10⁻³10⁻²10⁻¹110¹10²10³10⁴pivot, in order|pivot|negative, at row 15|S₁₁| from AH⁻¹Aᵀa failure with an addressrows before it14the pivot-128−(AH⁻¹Aᵀ)₁₁-128disagreement2.2·10⁻¹⁴it does not fail somewhereit fails at the constraint
Fig. 27 A larger problem, where the failure is at row fifteen and the argument is unchanged.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 2 constraints with κ(H) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 11.43 to 1.127·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 55.28 at the start of the sweep and 55.28 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: 2.137·10⁻⁶ against 1.053·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⁵1.1·10¹¹κ(ZᵀHZ), all stops55range-space forward error2.1·10⁻⁶null-space forward error1.1·10⁻¹¹both are the same algebraand only one squares
Fig. 28 Two constraints, where the Schur complement is 2 × 2 and its conditioning still squares.
The augmented step's two condition numbers, and the one the error obeysκ₂ is a worst case over all perturbations of the same NORM, and a normwise perturbation is allowed to put its whole budget on the smallest entry of a matrix. The componentwise number — Skeel's ‖ |A⁻¹||A||x| ‖ / ‖x‖ — is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. For the augmented interior-point matrix the two are 3.044·10¹³ and 13.25 at μ = 10⁻¹², a ratio of 2.296·10¹², and the measured error is 9.434·10⁻¹⁶ — which is the second number times the unit roundoff and has nothing to do with the first. The componentwise line is still finding its level down to about μ = 10⁻⁴, where the two groups of constraints have not yet separated, and from there it does not move in the fourth digit; where it settles — 13.25 here — is set by how many constraints are active rather than by μ. Eliminating the second block to reach the condensed form destroys the distinction: there the componentwise number is 1.324·10¹⁴, tracking the normwise one, because the large entries have been summed into CᵀDC and are no longer separately identifiable.-12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²3·10¹³componentwise, same matrix13their ratio2.3·10¹²measured error9.4·10⁻¹⁶every library prints the top lineand the error obeys the third
Fig. 29 And the number that describes the error on the field’s hardest matrix, five essays on.
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. 30 And where a counted inertia stops being reliable, from the spectra field.
MINRES on a saddle-point system under three Schur-complement approximations, and under noneThe same system at 12 unknowns and 3 constraints, solved four ways. With the exact Schur complement the preconditioned matrix has three distinct eigenvalues and the residual falls to 4.31·10⁻¹⁵ in three steps, after which nothing is left to remove. Replacing S by A diag(H)⁻¹Aᵀ costs 7 steps and replacing it by a scaled AAᵀ costs 7; the unpreconditioned system takes 18. The exact preconditioner is unaffordable — forming S costs 3 solves with H and a decomposition — so its value is as the statement the cheap ones are measured against, and the measurement needs no reference solution: the distance from {1 − φ, 1, φ} is a property of the approximation alone.02468101214161810⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1MINRES steprelative residualthree stepssteps to a residual of 10⁻¹¹exact S3diag(H)7scaled AAᵀ7none18three eigenvalues, three stepsand the approximations pay for the difference
Fig. 31 An iteration on the same kind of matrix, from two essays on.

Where this goes next

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 When symmetry is not enough The matrix [[0, 1], [1, 0]] is symmetric, nonsingular and perfectly conditioned, and there is no diagonal entry to pivot on. Every factorisation restricted to symmetric interchanges and one-by-one pivots fails on it, at any depth of searching, because every entry it could search is zero. The repair is to take two variables at once. 10¹10³10⁵10⁷10⁹10¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^11partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst9the bound, 2^112048both eliminations reach the same growthand only one of them had to swap to get there A factorisation with nothing to pivot for Cholesky's growth factor is not bounded by one. It is equal to one, at every size and every condition number, and the two-line reason is why the algorithm needs no pivoting at all — not "usually gets away without it". Its only failure is the square root of a non-positive number, which is exactly the test for definiteness, and in floating point that test moves with the precision. the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ The exact answer to a nearby problem A good algorithm does not give an approximate answer to your problem. It gives the exact answer to a problem very close to yours — and once that is the definition, a wrong result has two possible authors and they can be measured apart.

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

CholeskyCondition numberConstrained minimisationExact ground truthIndefinite matrixInertiaLagrange multiplierSaddle-point systemsSchur complementSymmetric indefinite