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

The regularisation that legalises every order

Perturb a saddle-point matrix's two blocks in opposite directions and it acquires a factorisation with a diagonal D under every symmetric permutation — not under a good one, under all of them. Five hundred random orderings, five hundred successes, and a growth factor that spans six orders across them.
40 min read 31 figures Where the algebra stops being the arithmeticMeasured, not assertedAssertions that reject

Worth reading first: The zero that is not a missing entry · A factorisation with nothing to pivot for · When the answer is a choice.

A sparse direct solver has two jobs that fight each other, and the sparsity field has already lost the fight once.

The ordering wants to be chosen from the sparsity pattern alone, before any number is looked at, because that is what lets the symbolic phase allocate the factor once and reuse it for a whole sequence of solves. The pivoting wants to be chosen from the numbers, because a pivot that is too small destroys the answer. A sparse LU has to compromise between them, and the threshold that governs the compromise is a knob with no good setting.

For one family of matrices the fight does not happen at all, and the family is exactly the one a constrained problem produces once it has been regularised.

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. 1 Five hundred random symmetric orderings of one matrix. Every one of them factorises. What each factorisation is worth is the horizontal axis.

Quasi-definite, and what the word covers

A symmetric matrix

[ E Aᵀ ] E symmetric positive definite, n × n [ A −F ] F symmetric positive definite, m × m

is quasi-definite, and Vanderbei’s theorem is that it has an LDLᵀ factorisation with D diagonal for every symmetric permutation. Not for a good permutation, not for one found by searching: for all of them.

The proof is two lines of the same congruence argument the field’s first essay uses. Any leading principal submatrix of a symmetric permutation of a quasi-definite matrix is itself quasi-definite; a quasi-definite matrix is nonsingular, because its inertia is (n, m, 0) and nothing is at zero; so every leading principal minor is nonzero, which is exactly the condition for an unpivoted LDLᵀ to exist.

The saddle-point matrix K is not quasi-definite: its (2, 2) block is zero, and zero is not −F for any positive definite F. Buy the theorem by perturbing both blocks —

K(δ, γ) = [ H + δI Aᵀ ] [ A −γI ]

— and the factorisation exists under every ordering, and the 2 × 2 pivots the symmetric indefinite essay had to introduce are not needed at all.

500 random symmetric orderings of a regularised saddle-point matrix, and what each factorisation reproducesA quasi-definite matrix — [[H + δI, Aᵀ], [A, −δI]] with δ = 0.01 — 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 301 of the 500. They span 1.17·10⁻¹⁶ to 6.76·10⁻¹⁵ — 2 orders — and the growth factor across them runs from 1 to 63.29 — which is 0.633/δ, 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 growth63growth × δ0.63the ordering is free to chooseand not free of consequence
Fig. 2 At a large regularisation, where the matrix is far enough from singular that no ordering is a bad one.

The count, and the control beside it

Five hundred random symmetric permutations of one matrix, factorised with no pivoting at all:

regularised at δ = 10⁻⁶ 500 of 500 factorise the same matrix at δ = 0 345 of 500

Thirty-one per cent of orderings break down outright on the matrix the problem posed. That is the control that makes the count a measurement: the theorem says something the matrix does not have for free, and the counterexamples are not rare.

The failures are not near-misses either. An unpivoted LDLᵀ fails when a pivot is exactly zero, and the zero block supplies them: permute two constraint rows to the front and the leading 2 × 2 submatrix is [[0, 0], [0, 0]]. The refusal the library publishes is the claim that a saddle-point matrix factorises under every ordering, fed the run at δ = 0 and required to fail.

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. 3 What has to be done instead when the theorem is not bought: two variables taken at once, from the elimination field.

What “every ordering” is actually claiming

Five hundred random permutations sounds like a sample and the theorem is not about a sample, so it is worth saying what the experiment is evidence for.

There are (n + m)! symmetric permutations of a 14 × 14 matrix, which is 8.7·10¹⁰, and no experiment visits a meaningful fraction of them. What the sample is good for is the control: it establishes that failures are common on the unregularised matrix — a third of the draws — so their complete absence on the regularised one is not the sample being lucky. If the failure rate were one in a thousand, five hundred draws would say nothing.

The theorem is what says all of them. The experiment says the theorem’s hypothesis is doing work, which is the part a proof cannot demonstrate and a reader is entitled to want. It also produces the thing the theorem does not have — the distribution of what those factorisations are worth — and that distribution is the next section and the reason the page is not just a citation.

120 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: 120 of 120 orderings factorise here, against 67.5 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 28 of the 120. They span 1.71·10⁻¹⁶ to 2.01·10⁻⁸ — 8 orders — and the growth factor across them runs from 1.03 to 5.112·10⁵ — which is 0.511/δ, 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.68worst growth5.1·10⁵growth × δ0.51the ordering is free to chooseand not free of consequence
Fig. 4 A smaller sample of the same experiment, where the share is unchanged and the tail is thinner because there are fewer draws to find it with.

Existence is not stability, and two numbers say so

The theorem says a factorisation exists. It says nothing about what the factorisation is worth, and across the same five hundred orderings:

growth factor 1.00 to 6.4·10⁵ residual ‖PKPᵀ − LDLᵀ‖/‖K‖ 1.6·10⁻¹⁶ to 1.0·10⁻⁷

Nine orders of residual, on a matrix every ordering of which is legal. The best ordering reproduces its own matrix to the rounding level; the worst reproduces it to seven digits, which on this site is a factorisation that has gone wrong and says so in the badge.

The refusal is the reading that the theorem invites: that a factorisation existing under every ordering is equally stable under every ordering, since the theorem does not distinguish between them. It is fed the best and worst residuals and required to reject the claim that they are within a factor of 1.5.

And the spread has a law. The worst growth factor over the orderings is 0.64/δ, an inverse relation that holds at every δ the figure is drawn at — so the regularisation that buys existence for every ordering also fixes, to within a constant, how badly the worst of them behaves. That is the whole trade in one number, and it is the reason the next section exists.

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 128 of the 500. They span 1.43·10⁻¹⁶ to 2.01 — 16 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. 5 At δ = 10⁻¹⁰, where the same five hundred orderings all still factorise and the worst residual is larger than the matrix.

What δ costs, which is exactly δ

The regularised solve is the exact answer to a different problem, so its error is proportional to δ. Measured across six decades the constant is 1,489: the relative error against the exact answer to the unregularised system is 1.49·10⁻¹¹ at δ = 10⁻¹⁴ and 1.43·10⁻³ at δ = 10⁻⁶, with the ratio error/δ constant to within a fifth of a per cent over the whole range.

That is a slope of one and not a trend, and it is the cleanest form the site’s regularisation trade-off has taken. In the regularisation field the left branch of the curve is a bias whose size is unknown, and choosing the parameter is a whole essay. Here the perturbation is known exactly — it is δ, and the code chose it — so the left branch is a straight line whose position can be computed before anything runs.

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. 6 The trade-off with no repair applied: the error is δ times a constant, over twelve decades.

And the left branch can be removed entirely

Refine against the unregularised matrix. Form the residual with K, solve the correction with the factorisation of K(δ, δ), and repeat. The perturbation was never in the residual, so it does not survive the iteration.

Measured with six steps: the refined error is at the rounding level — between 3.3·10⁻¹⁶ and 1.1·10⁻¹³ — at every δ from 10⁻¹⁴ up to 10⁻⁵, while the unrefined error at the top of that range is 1.4·10⁻². A gain of eleven orders, from three triangular solves per step against a factorisation that has already been computed.

That is what makes the regularisation practical rather than a compromise. The usual reading of a regularisation is that it trades accuracy for tractability; here it trades nothing for tractability, up to a limit, and the limit is the next section.

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. 7 With six steps of refinement, where the left branch has gone and the cliff has not moved.

Where refinement stops working, and it is not gradual

The refinement iteration’s error is multiplied at each step by roughly δ‖K⁻¹‖ = δ/σ_min(K), so it converges while that is under one and diverges above. On the figure’s matrix σ_min(K) = 6.80·10⁻⁴, and the measurement follows exactly: refined errors of 10⁻¹³ or better at every δ below 10⁻⁵, 4.9·10⁻⁷ at δ = 10⁻⁴, and 2.5·10⁻² at δ = 10⁻³.

More refinement steps do not move the cliff. Each step multiplies by the same factor, so below the cliff the flat region extends and above it the iteration diverges faster — which is what makes σ_min a property of the matrix rather than of the effort, and what makes the assertion a law rather than a level: the refined error is checked against plain·(δ/σ_min)^steps at every δ below the threshold, on every stop of the slider.

The practical reading is a design rule with a number in it. δ has to be large enough that the factorisation is stable — 0.64/δ is the growth — and small enough that refinement converges — δ < σ_min(K). Both ends are computable in advance, and on this matrix they leave about eight decades of room.

Iterative refinement on a factorisation of a matrix that is not the one being solvedThe last member of the sequence, solved with the first member's pivot order. Every other repair in this collection works because the outer loop recomputes the residual from the matrix, and this one does that too — and it does not work. The perturbed pivots are not a rounding, they are a change to the matrix, and the residual that would reveal it is computed with the same badly scaled rows that caused it. Six steps move the backward error from 4.82·10⁻⁹ to 5.46·10⁻¹⁰, a factor of 8.8, and stop. The equilibrated run has nothing to repair: it starts at 6.55·10⁻¹⁷.0123456710⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷refinement stepsbackward error of the solvethe kept order, in the units as giventhe same order, rows equilibratedthe repair that does not repairbefore refinement4.8·10⁻⁹after six steps5.5·10⁻¹⁰the factor gained8.8equilibrated, before6.5·10⁻¹⁷pivots replaced2the residual is recomputedand it cannot see what was done to the matrix
Fig. 8 Refinement failing for a different reason, in the field that owns it: a residual computed in the same precision as the solve.

The signs of D, which come free

An unpivoted LDLᵀ of a quasi-definite matrix has a diagonal D, and by Sylvester’s law the signs of its entries are the inertia. Across a hundred and twenty orderings of a matrix with n = 10 and m = 4, the count is ten positive and four negative in every order — the same integers the field’s first essay got from a congruence argument and from a spectrum.

That is a third route to the inertia and the only one that costs nothing: it is a by-product of a factorisation that was going to be computed anyway. The essay that turns it into an algorithm is in the spectra field, and it is where the property that the answer is an integer starts paying.

It is also a check. A code that regularises and then finds n + 1 positive pivots has a bug or a rank-deficient constraint, and it finds out for free, on every solve, without computing anything extra.

ν(σ), 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. 9 The same signs used deliberately rather than incidentally, from the spectra field.

The other thing that is bought, which is a memory layout

There is a consequence of a diagonal D that is easy to miss because it is about storage rather than about arithmetic.

A Bunch–Kaufman factorisation produces a block-diagonal D with a mixture of 1 × 1 and 2 × 2 blocks, and which is which is discovered during the factorisation. So the data structure has to accommodate either at every position, the solve has to branch, and a code that wants to reuse the factorisation’s structure across a sequence cannot: a different set of numbers gives a different arrangement of blocks.

A quasi-definite factorisation has a diagonal D and a lower-triangular L with a pattern the symbolic phase computed. There are no blocks, no branches and no discovery — the factorisation is a loop that fills in an array whose shape was fixed before it started. That is the difference between a routine that can be scheduled and one that cannot, and on a parallel machine it is the difference between a static task graph and a dynamic one.

None of that shows up in a residual, which is why it is worth stating separately from everything this page measures.

One entry 64× the rest, at 6-bit block significandsTwo bar charts: how many of 320 entries were rounded to zero, and the median entry's relative error, for the block format in two orderings and for a per-element format.entries rounded to zero, of 320block, as given310block, sorted by size22E4M3, either order0median entry's relative errorblock, as given1block, sorted0.013E4M30.022the same numbers, three waysdeleted, as given310deleted, sorted22deleted, per-element06-bit significands, blocks of 32sorting is free and changes no value
Fig. 10 A factorisation whose structure is fixed in advance, from the sparsity field, where the same property is what makes a schedule possible.

Why both blocks have to move

δ alone leaves the (2, 2) block at zero, which is not −F for any positive definite F: the theorem does not apply and the failures return. γ alone leaves E = H, which is enough when H is definite and not when it is only semidefinite — and the case where H is singular is common enough that a code cannot assume otherwise.

There is also no requirement that the two be equal, and in practice they are not. δ is a primal regularisation and is chosen against the objective’s scale; γ is a dual one and is chosen against the constraint’s. Codes pick them separately, adapt them during a run, and increase them when a factorisation reports too much growth — which is a feedback loop the 0.64/δ law makes predictable.

The library’s third refusal covers the reading that would make all of this unnecessary: that the regularisation is free. It is fed the sweep with no refinement and required to reject the claim that the error stays at the rounding level. The regularised matrix is a different matrix and its answer is a different answer; what refinement does is not make the perturbation harmless but remove it.

What the regularisation costs, and what 1 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 refinement1 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶2.1·10⁻⁶refined at δ = 10⁻²0.87the perturbation is known exactlybecause the code chose it
Fig. 11 With a single step, where the flat region reaches about as far as δ = σ_min/10 and no further.

What this makes possible

The regularisation is not an accuracy device and it is not a stability device. It is a scheduling device: it decouples the ordering from the numbers, so that the symbolic phase can run once for a whole sequence of solves and the numeric phase can be a pure evaluation.

That matters most in exactly the place the fifth essay in this field came from. An interior-point method solves one saddle-point system per iteration, dozens of times, with the same sparsity pattern and different numbers. Without quasi-definiteness the accurate form of that system needs dynamic pivoting and therefore a fresh symbolic analysis whenever the pivots move; with it, one ordering is computed at the start and every subsequent factorisation is a fill-in-the-blanks.

That is why the regularisation is standard in interior-point codes rather than a repair applied when something goes wrong, and it is why the next essay in the sparsity field can count fill before any number exists and have the count be exactly right.

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. 12 The consequence, counted: three orderings, a prediction made from the graph alone, and a factorisation that produces exactly that many entries.

Where the idea came from, and what it replaced

Before quasi-definiteness the accurate way to solve a sparse saddle-point system was a symmetric indefinite factorisation with threshold pivoting: allow the numeric phase to depart from the symbolic ordering when a pivot is too small, accept whatever extra fill that causes, and repeat the analysis when the departure is large. The sparsity field measured what that costs, and the answer was that the two objectives are genuinely opposed — every threshold setting buys fill with growth or growth with fill.

The regularisation does not find a better compromise on that curve. It gets off the curve, by changing the matrix so that the numeric phase has no reason to depart. What it pays instead is a perturbation of known size, which the previous two sections show can be removed. Trading an unbounded and unpredictable cost for a bounded and removable one is what makes it worth doing, and it is a shape that recurs: static pivoting does the same trade in the unsymmetric case, and pays with a perturbation that refinement sometimes cannot remove.

The difference between the two is instructive. Static pivoting perturbs whichever entries turn out to be too small, so the perturbation is discovered rather than chosen and its size is not known in advance. A regularisation perturbs a diagonal by a number the code picked before it started. The first is a repair and the second is a design, and only the second has a law like 0.64/δ attached to it.

Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹⁶, before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.2 0 interchangesno pivoting · rows scaled1 0 interchangespartial · as given1.1·10⁻¹⁶ 1 interchangepartial · rows scaled1 0 interchangesscaled partial · as given1.1·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given2.2·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given1.1·10⁻¹⁶partial, rows scaled1its relative residual5·10⁻¹⁷complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 13 The unsymmetric version of the same trade, from the sparsity field.

At other settings

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. 14 At δ = 10⁻⁸, midway along the inverse law.
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 151 of the 500. They span 1.34·10⁻¹⁶ to 1.23·10⁻¹¹ — 5 orders — and the growth factor across them runs from 1 to 6428 — 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 growth6428growth × δ0.64the ordering is free to chooseand not free of consequence
Fig. 15 And at 10⁻⁴, where the worst residual has come back to the rounding level and the answer has moved by a tenth of a per cent.
What the regularisation costs, and what 2 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 refinement2 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶2.9·10⁻⁹refined at δ = 10⁻²0.81the perturbation is known exactlybecause the code chose it
Fig. 16 Two refinement steps, for reading the cliff’s position against the six-step figure.
What the regularisation costs, and what 10 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 refinement10 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶2.5·10⁻¹⁵refined at δ = 10⁻²0.48the perturbation is known exactlybecause the code chose it
Fig. 17 And ten, where the cliff is in exactly the same place.
Fill against growth as the pivot threshold moves, on the 6×6 gridTwo curves against the pivot threshold on a logarithmic horizontal axis. One falls steeply from left to right; the other rises gently. A vertical line marks the value libraries default to.10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L+U, ÷ entries in Agrowthfillthe library default‖PA − LU‖/‖A‖ at τ = 0.12.9·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1372one knob, two measurements, opposite directionsand the default is most of both
Fig. 18 The compromise this page avoids, from the sparsity field: fill against growth, with a knob.
The same first pivot decision, with the corner entry at 10⁻¹²Three sparsity patterns side by side: the matrix, the factor obtained by eliminating the corner first, and the factor obtained by pivoting on the largest entry instead.the matrix105 entriescorner first — sparsest227 entries, growth 1.9·10¹¹largest first — safe242 entries, growth 1.19the middle factor is the smaller one, and its answer has no correct digitsboth factorisations reproduce the matrix‖PA − LU‖/‖A‖, sparsest3.8·10⁻¹⁷‖PA − LU‖/‖A‖, pivoted5.4·10⁻¹⁷forward error, sparsest3·10⁻⁵forward error, pivoted4.8·10⁻¹⁶red marks are entries elimination createdthe fill argument and the stability argument disagree
Fig. 19 And the case that made the compromise necessary there.
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. 20 The unregularised matrix, whose inertia the signs of D reproduce.
Iterative refinement of the inverse-and-multiply solve at κ = 10^8The backward error starts at 3.3·10⁻¹¹ and falls by about κu a step — 1.8·10⁻¹⁷, 2.7·10⁻¹⁷, 2.4·10⁻¹⁷, 2.1·10⁻¹⁷, 2.1·10⁻¹⁷ — reaching the LU route's 3.1·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 1.6·10⁻⁸ to 1.2·10⁻⁹, against the LU route's 9.6·10⁻¹⁰. What refinement at working precision buys is stability, and the accuracy floor belongs to the problem.01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 3.1·10⁻¹⁷LU route: forward 9.6·10⁻¹⁰forward errorbackward errorwhat a correction buysη before refinement3.3·10⁻¹¹η after four steps2.1·10⁻¹⁷forward, unchanged1.2·10⁻⁹cost of a step, flops2048the residual is repairableand the accuracy floor is the problem's
Fig. 21 Refinement working as intended, in the field where it was introduced.
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. 22 The system this regularisation is applied to most often.
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.7312, and the pivot the routine actually reports is -3.7312 — the same number to 15 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-3.7−(AH⁻¹Aᵀ)₁₁-3.7disagreement3.6·10⁻¹⁶it does not fail somewhereit fails at the constraint
Fig. 23 What happens without it to the factorisation everyone reaches for first.
Entries in U, against the bound the symbolic phase can compute, on the 5×5 gridA row of horizontal bars. The topmost is longest and is labelled as the bound; every measured bar below it is shorter.the bound211no pivoting127τ = 0.00198τ = 0.003100τ = 0.01106τ = 0.03111τ = 0.1116τ = 0.3117τ = 1138entries in Ua bound, and its slackthe bound, from the graph alone211the worst that occurs138loose by1.5no arithmetic was done to compute the bound — only the pattern of AᵀA and its elimination graphGeorge and Ng: U fits inside chol(AᵀA), whatever the swapsprovable, cheap, and loose
Fig. 24 The symbolic phase’s prediction in the case where it can only be a bound.
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. 25 And the two routes that avoid factorising K at all.
300 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: 300 of 300 orderings factorise here, against 67.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 83 of the 300. They span 1.56·10⁻¹⁶ to 1.02·10⁻⁷ — 9 orders — and the growth factor across them runs from 1 to 5.112·10⁵ — which is 0.511/δ, 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.67worst growth5.1·10⁵growth × δ0.51the ordering is free to chooseand not free of consequence
Fig. 26 Three hundred draws rather than five hundred, where the tail is thinner and the share is the same.
Fill under three orderings of a sparse saddle-point matrix, predicted symbolically and then countedA stiffness block of 32 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 145, minimum-degree 79, reverse-cuthill-mckee 79: 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.natural145 predicted · 145 countedminimum-degree79 predicted · 79 countedreverse Cuthill–McKee79 predicted · 79 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural145minimum-degree79reverse-cuthill-mckee79predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 27 A larger sparse problem, where the orderings separate further.
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 400.1·δ 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) = 0.002082, 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 ÷ δ, unrefined400σ_min(K)0.0021refined at δ = 10⁻⁶3.5·10⁻¹⁵refined at δ = 10⁻²0.24the perturbation is known exactlybecause the code chose it
Fig. 28 With a worse-conditioned objective, where σ_min is smaller and the cliff moves left.
What a factorisation survives across one barrier step, against the step's reduction factorFrom one interior-point iteration to the next, H, C and the sparsity pattern are identical and exactly 8 entries of the matrix change — the diagonal of the (2, 2) block, and nothing off it. That is the sparsity pattern reused for free forever. The factorisation is a different question, and the answer is a threshold: carrying the decomposition of K(μ) to K(σμ) and cleaning up with iterative refinement holds at the rounding level while σ is within a per cent of one, and by σ = 0.9 it buys one step. At the schedules an interior-point method actually uses — σ between 0.5 and 0.1 — the first reuse is already at 0.915 and 2.59·10¹¹. The reason is in the entry count: the 8 entries that moved are the 8 that dominate the matrix, and they moved by a factor of 1/σ, so the relative change in K across one step at σ = 0.1 is 9 rather than the 10⁻³ the collection's reuse essays are about.-1-0.75-0.5-0.25010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴log₁₀ σ — the barrier's reduction factorresidual after one reused stepconvergedthe pattern free, the factors notentries moved8off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 29 And the sequence this regularisation is applied inside.
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. 30 Where a counted inertia fails, which the signs of D inherit.
The constraint preconditioner's spectrum, over six decades of the constraint's condition numberP = [[G, Aᵀ], [A, 0]] with G an approximation to H — here a well-conditioned matrix that is not H at all. The theorem says P⁻¹K has eigenvalue 1 with multiplicity 2m = 8 and n − m = 6 others, which are the generalised eigenvalues of the pencil (ZᵀHZ, ZᵀGZ) for any basis Z of the null space of A. A does not appear in that list, and the measurement is the flat lines: κ(A) crosses six decades along the horizontal axis and the 6 eigenvalues move by 7.14·10⁻⁶ relative, which is the arithmetic. A preconditioner for a constrained problem can decline to know anything about the constraint, because the constraint has already been inverted exactly inside it. The one quantity that does move is the drift of the eigenvalues the theorem puts at exactly one: 6.65·10⁻⁹ at κ(A) = 1 and 6.04·10⁻⁵ at 10⁶, which is κ(A) times the unit roundoff.012345610⁻³10⁻²10⁻¹110¹log₁₀ κ(A)eigenvalue of P⁻¹Kthe constraint is not in itnontrivial6at one8movement, six decades7.1·10⁻⁶drift at one6·10⁻⁵2m at onethe marks are a pencil that never saw Aand the lines are the preconditioned matrix
Fig. 31 And a preconditioner that keeps the zero block rather than perturbing it.

Where this goes next

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 An ordering that does not wait for the numbers A sparse factorisation's memory is decided by an ordering computed from the graph, and its stability by pivots computed from the values, and the two decisions fight. On one family of matrices they do not — the ordering can be chosen for fill alone, and the fill the symbolic phase predicts is the fill the factorisation produces — exactly, not as a bound. 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. the matrix105 entriescorner first — sparsest227 entries, growth 1.9·10¹¹largest first — safe242 entries, growth 1.19the middle factor is the smaller one, and its answer has no correct digitsboth factorisations reproduce the matrix‖PA − LU‖/‖A‖, sparsest3.8·10⁻¹⁷‖PA − LU‖/‖A‖, pivoted5.4·10⁻¹⁷forward error, sparsest3·10⁻⁵forward error, pivoted4.8·10⁻¹⁶red marks are entries elimination createdthe fill argument and the stability argument disagree Structure and stability stop being separable The sparsest variable to eliminate on this matrix has a diagonal entry of 10⁻¹². Eliminating it produces the smaller factor, reproduces the matrix to 3.8·10⁻¹⁷ — better than pivoting does — and returns an answer wrong in the fifth digit.

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.

Bunch–KaufmanGrowth factorInertiaIterative refinementLDLᵀ factorisationQuasi definite matrixRegularisationSaddle-point systemsSymbolic factorisationSymmetric permutation