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.
22 min read 7 figures Where the algebra stops being the arithmeticMeasured, not assumedClaims that can fail

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 and −γI on its diagonal, A and its transpose off it, with δ = 10⁻⁶ and γ = 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 in the zero block's perturbation 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.

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.

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 and −γI on its diagonal, A and its transpose off it, with δ = 10⁻⁶ and γ = 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 in the zero block's perturbation 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. 2 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 and −γI on its diagonal, A and its transpose off it, with δ = 10⁻¹⁰ and γ = 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 in the zero block's perturbation 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. 3 At δ = 10⁻¹⁰, where the same five hundred orderings all still factorise and the worst residual is larger than the matrix — 2.01, which is no correct digits at all.
500 random symmetric orderings of a regularised saddle-point matrix, and what each factorisation reproducesA quasi-definite matrix — H + δI and −γI on its diagonal, A and its transpose off it, with δ = 10⁻⁴ and γ = 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 in the zero block's perturbation 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. 4 And at δ = 10⁻⁴, six decades the other way. The same five hundred orderings, the same 100% against 69.0%, and a worst residual of 1.23·10⁻¹¹.

Six more decades, and the two columns the figure is about have still not moved by a digit while the third has moved by seven:

500 random symmetric orderings of a regularised saddle-point matrix, and what each factorisation reproducesA quasi-definite matrix — H + δI and −γI on its diagonal, A and its transpose off it, with δ = 0.01 and γ = 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 in the zero block's perturbation 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. 5 The largest δ the slider carries. Still 100% against 69.0%, and now the worst ordering’s residual is 6.76·10⁻¹⁵ — within two orders of the best one’s.

Legality does not depend on δ and accuracy depends on it enormously. Across the whole slider the regularised fraction is 100% and the bare fraction is 69.0% — the same two numbers to three figures at δ = 10⁻¹⁰, 10⁻⁸, 10⁻⁶, 10⁻⁴ and 10⁻². Whatever δ is doing for legality, it has finished doing it by the smallest value on offer.

The worst residual, over the same five stops, reads 2.01, 6.17·10⁻⁴, 1.02·10⁻⁷, 1.23·10⁻¹¹ and 6.76·10⁻¹⁵. That is eight decades of δ buying fifteen decades of residual, a slope of about −1.8 — not far off the square, and certainly not the linear relationship the legality column might suggest. The best ordering’s residual meanwhile sits at 1.2 to 1.7·10⁻¹⁶ at every stop: some orderings are fine at any δ, and what δ decides is how badly the worst one behaves.

So the two halves of the trade run in opposite directions at very different rates, and the essay has both of them measured. The section below prices the other one: the regularised solve is the exact answer to a perturbed problem, and its error is 1489·δ on that figure’s family — linear in δ, with a constant large enough to matter. Residual falling as δ⁻¹·⁸ against perturbation rising as δ means there is an interior optimum in δ and it is not at either end of this slider.

The two numbers come from different figures on different families, so this page does not multiply them together and name the optimum. What it can say is the shape: the usual instinct — pick δ as small as the arithmetic tolerates, on the grounds that a smaller perturbation is a smaller lie — is the wrong end. At δ = 10⁻¹⁰ the perturbation is negligible and the factorisation has no correct digits.

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.

The two laws multiplied

Two laws have now been measured separately — the worst growth over the orderings is 0.5/δ, and the regularisation error is 1,430·δ — and putting them together is what decides δ. The worst ordering’s residual falls as about δ⁻², so the total error against the unregularised answer is 1,430δ + c/δ²:

δ worst ordering’s residual regularisation error total
10⁻¹⁰ 2.0 1.4·10⁻⁷ 2.0
10⁻⁸ 1.0·10⁻⁴ 1.4·10⁻⁵ 1.2·10⁻⁴
10⁻⁶ 2.0·10⁻⁸ 1.4·10⁻³ 1.4·10⁻³
10⁻⁴ 7.9·10⁻¹³ 1.3·10⁻¹ 1.3·10⁻¹

A minimum at δ ≈ 10⁻⁸, and a floor of about 10⁻⁴. Four digits, on a problem a stable factorisation answers to sixteen. That is the price of reading every ordering is legal as a licence to use any of them.

The median ordering pays none of it. Its residual is 5·10⁻¹⁵ at every δ from 10⁻¹⁴ to 10⁻² — flat, twelve decades of it — so for a typical ordering the factorisation term never enters at all, the total is the regularisation error alone, and there is no optimum to find: δ should be as small as the solve allows, and 10⁻¹⁴ gives 1.5·10⁻¹¹.

So the theorem is worth having as insurance and worth nothing as a licence, and the difference between those two readings is twelve digits. A code that orders for fill and relies on the theorem only to know that the factorisation will exist pays nothing for the guarantee and should take δ as small as it can. A code that reads the theorem as permission to stop caring about the ordering acquires an optimal δ, and its best available answer is four digits.

There is a reading of the flat median that is worth separating from the arithmetic, because it says what the theorem is really buying. Five hundred random permutations of a fourteen-by-fourteen matrix are five hundred draws from a space of 8.7·10¹⁰, and the ones with a large growth factor are the ones that happen to place several constraint rows early — a structure a random draw finds occasionally and a fill-reducing ordering never produces, because it is looking at the same rows for a different reason and rejecting them. The tail of that distribution is a set of orderings nobody would choose.

So the two readings correspond to two populations rather than to two attitudes. Insurance is the statement that the factorisation exists on the ordering a code would pick anyway; licence is the statement that it exists on an ordering picked adversarially, and only the second has a price. That distinction is invisible in the theorem, which quantifies over all permutations equally, and it is the distribution the experiment produces that makes it visible — which is the reason this page is not a citation.

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 σₘᵢₙ(K) = 6.797·10⁻⁴, marked on the axis, and stops working above it: the iteration's contraction factor is δ/σₘᵢₙ 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δ = σₘᵢₙ(K)no refinement0 stepsa left branch that can be removederror ÷ δ, unrefined1489σₘᵢₙ(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.

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.

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 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. 7 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 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.

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.

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.

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.

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