Sparsity, and what elimination costs

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.

Worth reading first: Structure and stability stop being separable · The factor is not sparse · The order decides the memory.

The essay that opened this anchor is about a fight. A sparse factorisation’s memory is decided by an ordering computed from the graph, before any number is looked at; its stability is decided by pivots computed from the values, during the factorisation; and the second can force the first to be wrong. The essay after it measured the consequence: the symbolic phase can only produce an upper bound, and the allocation has to be able to grow.

There is a family on which the fight does not happen, and the reason is one line of algebra rather than a heuristic.

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. 1 Three orderings of one sparse saddle-point matrix, with the fill each produces. The two numbers under each panel are a prediction and a count.

The family, and the line

A quasi-definite matrix — symmetric, with a positive definite (1, 1) block and a negative definite (2, 2) block — has an LDLᵀ factorisation with a diagonal D under every symmetric permutation. The constraint field’s essay on it has the proof and measures the claim against five hundred random orderings.

The consequence for this field is immediate and it is the whole page. If a factorisation exists in every order, the numeric phase has no reason ever to depart from the order the symbolic phase chose. There is no pivot search, no threshold, no delayed pivot, no dynamic reallocation — the numeric phase is a loop that fills in an array whose shape was fixed before it started.

So the ordering can be chosen for fill alone, which is what the fill-reducing heuristics were designed to do and what nothing else lets them do without a numerical veto.

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 5 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 117, minimum-degree 71, reverse-cuthill-mckee 71: 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.natural117 predicted · 117 countedminimum-degree71 predicted · 71 countedreverse Cuthill–McKee71 predicted · 71 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural117minimum-degree71reverse-cuthill-mckee71predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 2 With constraints touching five unknowns apiece, where there is more fill to reduce and the orderings separate further.

The measurement

A stiffness block of twenty-four unknowns with four constraints, each touching three consecutive variables — a bordered band, which is what a finite-element problem with multipoint constraints actually looks like. Three orderings, all computed from the sparsity pattern alone:

ordering predicted counted natural 113 113 minimum degree 63 63 reverse Cuthill–McKee 63 63

Predicted from the graph before any number is touched, and matched entry for entry by the factorisation of the regularised matrix. Not bounded — matched.

That is the claim the anchor’s earlier essay could not make. There the symbolic count is an upper bound and the numeric count comes in under it when cancellation makes a structurally nonzero entry numerically zero, or over it when a pivot moves. Here the second cannot happen, and the first is checked: the counts agree exactly, so no cancellation occurred either.

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. 3 The same comparison in the case where the two numbers are allowed to differ, from this anchor’s own earlier essay.

What “exactly” is claiming, and how it could fail

Two numbers agreeing is a weak-looking claim, and it is worth saying what would break it, because both failure modes are real elsewhere in this field.

The prediction could be too low if the numeric phase created an entry outside the allocated pattern. That is what a moved pivot does, and it is why the anchor’s earlier essay reports a bound: with dynamic pivoting an entry can appear anywhere, and a code has to be able to grow its arrays. Quasi-definiteness removes that possibility entirely, because there is no pivot search to move anything.

The prediction could be too high if a structurally nonzero entry came out numerically zero. That happens through cancellation, and it happens for real: on a matrix with symmetry or with repeated values, two contributions to a fill entry can cancel exactly. The symbolic count then over-allocates, harmlessly, and the two numbers differ.

Neither happens here. The counts agree in both directions at every ordering and every setting of the slider, which says the allocation is exact rather than merely safe — and an exact allocation is what lets a code place the factor in a fixed buffer and reuse it across a whole sequence without checking.

Columns needed by the block between a separator's two halves, against the length of the separatorThe block is 3×4, 5×6, 7×8, 11×12 across the sweep, so its size more than triples. In the separator's own numbering — the unknowns in the order they sit on the line — it needs 3, 4, 5, 6 columns at 10⁻⁸: about one more per doubling, which is a logarithm. Renumber the same Schur complement by one symmetric permutation and the same block needs 3, 5, 7, 11 — all of them at eight digits and tighter, and a multiple of the geometric count at every tolerance drawn. Nothing about the matrix changed. What changed is that a cluster of a shuffled numbering is scattered along the whole separator, so every block mixes near interactions with far ones and there is no smooth kernel left to compress. The ordering that works here is not chosen: it is the one the separator's geometry already has, which is why nested dissection leaves a compressible fill and a fill-reducing ordering chosen on the graph alone need not.0481216202402468101214unknowns on the separatorcolumns above 10⁻⁸the same matrix, renumberedin the separator's own orderthe ordering the geometry hands overseparator 73separator 236renumbered, largest11the block, largest11share of the square stored0.52the fill is totaland it is not independent
Fig. 4 The other way a symbolic count can be wrong, from the field’s own essay on it.

Two orderings that tie, and one that does not

Minimum degree and reverse Cuthill–McKee both give 63 on this matrix, and the natural order gives 113 — a factor of 1.8. Two things about that are worth saying.

The tie is a property of the shape rather than a coincidence. A bordered band is nearly the case both heuristics were designed for: RCM minimises bandwidth and the matrix is banded except for the border, and minimum degree eliminates the low-degree interior first, which for a band is the same order. On a two-dimensional grid the two part company decisively, and on an arrowhead they part company the other way.

The factor of 1.8 is what is being bought, and it is small because the problem is small. Fill grows superlinearly, so the ratio between a good ordering and a bad one grows with the problem: on the grids this field measures elsewhere it reaches an order of magnitude, and it is the difference between a factorisation that fits in memory and one that does not.

What the page adds is that on this family the ratio is available. A heuristic that produced a beautiful ordering the numeric phase then abandoned would be buying nothing.

Nonzeros in the Cholesky factor of the 12×12 grid Laplacian, by orderingA horizontal bar chart comparing the number of nonzero entries in the Cholesky factor under four elimination orderings, with reference marks for the matrix itself and for a dense factor.natural1739reverse Cuthill–McKee1354minimum degree1026nested dissection1413matrix: 408 entries · dense factor: 10440bandwidth 12 · 4.26× the matrixbandwidth 12 · 3.32× the matrixbandwidth 123 · 2.51× the matrixbandwidth 108 · 3.46× the matrixn = 144, five-point stencilevery ordering fills in; none avoids it
Fig. 5 The orderings ranked, in the field’s own terms.

What the border does to the graph

The constraint rows are what make this a saddle-point matrix and they are also what makes its graph interesting, so it is worth reading the panels as graphs rather than as pictures.

A constraint touching t unknowns is a vertex connected to t others. Eliminating it makes those t mutually adjacent — a clique of size t — which is t(t−1)/2 edges, most of them new. So a constraint row is an expensive vertex to eliminate early and a cheap one to eliminate late, and both fill-reducing heuristics discover that without being told: in every panel the constraint vertices are pushed towards the end of the ordering.

That is also why the fill rises with the slider. A constraint touching five unknowns creates a clique of ten edges where one touching two creates a single edge, so the cost of the border grows quadratically in how local the constraints are. A model with a few global constraints — a total that must sum to one, a mean that must vanish — is the worst case: one vertex adjacent to everything, whose elimination fills the entire remaining matrix.

The usual repair for that is to eliminate such a row last and treat it as a low-rank update, which is the same trick another field uses for a different reason.

The projected matrix, against the basis it is supposed to project ontoThree quantities against the cycle count on a logarithmic vertical axis. The disagreement between the arrowhead and the Rayleigh quotient of its own basis is at the level of rounding for two cycles and reaches 4.23·10⁻⁶ after. The basis's own orthogonality error stays at 9·10⁻¹⁶ throughout, and the run whose border entries are recomputed stays at 2.5·10⁻¹⁵.1234567810⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³cyclerelative sizea vector kept twice‖T − QᵀAQ‖‖QᵀQ − I‖border recomputedwhich model is wrong‖T − QᵀAQ‖ at cycle 21.1·10⁻¹⁵at cycle 34.2·10⁻⁶‖QᵀQ − I‖, worst9·10⁻¹⁶the basis is orthonormal throughoutand the matrix on it is not its own
Fig. 6 The extreme case, from this field’s own essay: one row adjacent to everything.

What is being given up, and it is not nothing

The theorem needs the matrix to be quasi-definite, and a saddle-point matrix is not: its (2, 2) block is zero, and zero is not negative definite. The property is bought by perturbing both blocks — δ added along the primal diagonal, −γ along the dual one — and the perturbation is a real change to the answer.

The constraint field’s essay measures what it costs: an error proportional to δ, with the constant computable, and removable by refining against the unregularised matrix as long as δ stays below σ_min(K). So the price is a few triangular solves per right-hand side, and a design constraint on δ with both ends known.

There is a second price that is easier to miss. The theorem guarantees the factorisation exists; it says nothing about the growth, and the growth across five hundred orderings spans 1 to 6.4·10⁵ with the worst residual nine orders above the best. A code that chose its ordering purely for fill and never looked at the numbers would sometimes get a factorisation that reproduces its own matrix to seven digits.

So the fight has not been abolished. It has been converted from a constraint on the ordering into a bound on the regularisation, and the second is a single scalar that can be adapted during a run.

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. 7 The spread of what the orderings are worth, which the fill counts say nothing about.

Why the constraint block cannot simply be ordered away

A reader who has followed the graph argument may ask why the ordering does not simply put all the constraint rows last and be done with it, since a constraint vertex is expensive to eliminate early.

It does, roughly, and that is what both heuristics discover. What it cannot do is make the constraint rows free, and the reason is the same clique argument seen from the other end. Eliminating everything else first leaves a dense m × m Schur complement among the constraint rows — the matrix AH⁻¹Aᵀ that the constraint field’s second essay forms explicitly — and that block is dense whenever any two constraints share a variable, whatever order they are eliminated in.

So the fill has a floor of m(m + 1)/2 entries that no ordering removes, and the whole exercise is about the fill above that floor. On the figure’s matrix with m = 4 the floor is ten entries out of sixty-three, and on a problem with a thousand constraints it is half a million — at which point the border, not the interior, is the memory.

That is a shape worth recognising because it changes which method wins. When m is large enough for the dense Schur block to dominate, factorising the whole matrix stops being cheaper than the null-space route, and the choice goes back to being about the two eliminations.

The Schur complement left on a 15-unknown separator, shaded by the size of each entryA 15 × 15 grid, its middle column taken as a separator, and both halves eliminated exactly. What is left on the separator is 100 per cent nonzero — the sparsity field's result, unchanged: eliminating a variable couples everything it touched, and by the end everything is coupled to everything. The shading is what that field does not measure. The entries fall away smoothly from the diagonal, because the Schur complement is a discrete Green's function — the operator mapping data on the separator to response on the separator — and away from the diagonal that is an integral operator with a smooth kernel. The outlined block is the 7 × 8 between the separator's two halves: every entry nonzero, and 5 columns describe it to eight digits. The rank structure was not put there by the elimination. It was in the differential operator before anything was discretised.the separator's 15 unknowns, in the order they sit on the lineoutlined: the block between the two halvesdense, and not independententries nonzero1the block56columns it needs5numbers stored150entries in the square225every entry is nonzeroand six columns describe them
Fig. 8 The floor, drawn: what is left among the constraint rows after everything else has gone.

Against the alternative, which is threshold pivoting

The other way to factorise a sparse saddle-point matrix accurately is to allow the numeric phase to depart from the ordering when a pivot is too small, and to accept the extra fill. This anchor’s threshold essay prices that: every setting of the threshold buys fill with growth or growth with fill, and there is no setting that is free.

The comparison is not that one is better. It is that they are different kinds of cost.

Threshold pivoting’s cost is unpredictable: how much extra fill appears depends on the values, so the allocation cannot be sized in advance, the analysis may have to be repeated, and a sequence of similar problems can have different memory requirements. Regularisation’s cost is a number chosen before the run, and refinement removes most of it.

For a sequence of solves that difference dominates everything else, and it is why interior-point codes regularise rather than threshold. One symbolic analysis at the start, dozens of numeric factorisations with identical structure, and a memory requirement known before the first one runs.

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. 9 The curve this family gets off, from the anchor’s own essay.

And the ordering survives the sequence too

The sequence field asks a related question of a different family: whether an ordering computed for one matrix is still good for the next one. Its answer there is that the pattern has to be the same, and when it is, the ordering is reusable and the pivots are not.

An interior-point sequence is the cleanest possible case of that. From one barrier parameter to the next, only the p diagonal entries of the (2, 2) block change, and no off-diagonal entry moves at all — measured, at every step. The pattern is not merely similar; it is identical.

So the symbolic phase runs once for a whole optimisation. And the numeric factorisation, which the next essay in the sequence field shows cannot be carried across even a single step, is recomputed each time into an array that never changes shape.

That is the division this page is really about: a fixed structure and moving numbers, with the structure computed once. It is what makes a sparse direct method viable inside an outer iteration, and it is available here because of one line about the signs of two blocks.

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 6 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 1.02 and 2.88·10¹¹. The reason is in the entry count: the 6 entries that moved are the 6 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 moved6off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 10 What cannot be carried across a step, from the sequence field.

What a fixed structure buys beyond memory

Predictable fill is usually argued for as an allocation question, and on a modern machine it is worth more than that.

A factorisation whose pattern is known in advance has a known elimination tree, and the tree is a task graph: which columns can be factorised independently, in what order, and with what dependencies. A code can build that schedule once, before any number exists, and then execute it — on several cores, or on a machine where the cost is words moved rather than operations counted.

A factorisation with dynamic pivoting cannot. A moved pivot changes the tree, so the schedule has to be rebuilt, the load balance changes and the communication pattern with it. Codes that do this well spend a great deal of effort on it, and the effort is the price of the numerical freedom.

None of that appears in a fill count and none of it appears in a residual, which is why it is worth naming separately. On the family this page is about, the schedule is as static as the allocation, and both are computed from the graph.

A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 11 The other half of what a fixed structure is worth, in the field that counts messages.

The refusal

The claim to close is the one this anchor established two essays ago and which is true everywhere except here: a symbolic phase can only bound the fill, because the numeric phase may move a pivot and every moved pivot adds entries the analysis did not allocate.

The assertion is fed the three pairs — 113 against 113, 63 against 63, 63 against 63 — and required to reject the claim that they differ. It does.

The refusal is worth reading as a boundary rather than as a correction. The earlier essay’s claim is right about sparse LU and about symmetric indefinite factorisation with threshold pivoting, which between them cover nearly everything this field studies. What this page adds is a family where it fails, and the family is characterised by an algebraic property rather than by a numerical accident — so it can be checked for rather than hoped for.

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 2 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 111, minimum-degree 59, reverse-cuthill-mckee 59: 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.natural111 predicted · 111 countedminimum-degree59 predicted · 59 countedreverse Cuthill–McKee59 predicted · 59 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural111minimum-degree59reverse-cuthill-mckee59predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 12 With the most local constraints the figure draws, where the border costs least.

Where the family is, and where it is not

The property is algebraic, so it can be checked for rather than hoped for, and it is worth knowing which problems have it.

Interior-point methods for linear and quadratic programming, which is where the regularisation is standard and where the sequence argument above applies in full.

Least squares written as an augmented system — [[I, A], [Aᵀ, 0]] — which is quasi-definite after regularising the second block and is the stable alternative to the normal equations that the least-squares field spends an essay on.

Mixed finite-element formulations for flow and for elasticity, where the (1, 1) block is a discretised operator and is definite, and where a small dual regularisation is often present already for reasons of the discretisation rather than of the arithmetic.

Where it is not: any problem whose (1, 1) block is indefinite. A saddle point of a nonconvex objective has a Hessian with negative curvature in it, and no diagonal perturbation of bounded size makes it definite — a diagonal shift δ would have to exceed the most negative eigenvalue, which is a change to the problem rather than to its conditioning. Those need the 2 × 2 pivots after all, and the fight this page describes resumes.

The curvature pᵀAp/pᵀp along conjugate gradients on a matrix with one eigenvalue at −10^-1The matrix is 40×40, positive definite apart from a single eigenvalue at -0.1. Conjugate gradients divides by pᵀAp at every step, and the quantity runs 5.43, 2.72, 1.71, 0.687 before turning negative at step 6, where it is -0.02656. That direction is a proof: one matrix–vector product from outside the iteration confirms it, and it recovers 26.6 per cent of the eigenvalue. MINRES on the same system never forms this quantity, meets nothing, and returns the solution in 37 steps.012345678-7-5-3-11357conjugate gradient steppᵀAp ⁄ pᵀpλ_min = -0.1positive: a step existsnegative: a certificate existsone matrix, two questionsstep it turns at6quotient there-0.027share of λ_min recovered0.27λ_min, by construction-0.1MINRES steps on the same system37the division that cannot be doneis the answer to a different question
Fig. 13 Negative curvature, from the field where it is useful rather than an obstacle.

At other settings

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 4 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 115, minimum-degree 67, reverse-cuthill-mckee 67: 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.natural115 predicted · 115 countedminimum-degree67 predicted · 67 countedreverse Cuthill–McKee67 predicted · 67 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural115minimum-degree67reverse-cuthill-mckee67predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 14 Four, midway along the slider.
The 12×12 grid Laplacian and its Cholesky factor, ordered by naturalTwo square sparsity plots side by side. The left shows the nonzeros of the matrix; the right shows the nonzeros of its Cholesky factor, with the entries created by elimination marked in a second colour.the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding
Fig. 15 The same kind of picture in the field’s own first essay.
Fill growth under natural: the factor rises as n^1.49A log-log plot of nonzero count against matrix dimension. The matrix's own count is a straight line of slope one; the factor's is steeper; a dense factor is steeper still.10²10².³10²10³10⁴n (dimension)nonzerosdense factorfactor, n^1.49matrix, n^1.05grid Laplacians from 5×5 to 12×12fitted, not quoted
Fig. 16 How fill grows with the problem, which is why the factor of 1.8 here is a small number about a small matrix.
The matrix and its LU factors at τ = 0.1Two sparsity patterns. The left is the matrix; the right is L and U together, with the entries elimination created drawn in a second colour.A156 entries, 36 unknownsL + U372 entries, 23 interchanges‖PA − LU‖/‖A‖2.9·10⁻¹⁶growth factor38fill created216forward error1.5·10⁻¹⁵the 6×6 grid with a drift and one ruinous cornerred marks are entries that were zero
Fig. 17 The unsymmetric case, where the pivots move and the pattern with them.
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. 18 The conflict this family avoids.
The Schur complement left on a 15-unknown separator, shaded by the size of each entryA 15 × 15 grid, its middle column taken as a separator, and both halves eliminated exactly. What is left on the separator is 100 per cent nonzero — the sparsity field's result, unchanged: eliminating a variable couples everything it touched, and by the end everything is coupled to everything. The shading is what that field does not measure. The entries fall away smoothly from the diagonal, because the Schur complement is a discrete Green's function — the operator mapping data on the separator to response on the separator — and away from the diagonal that is an integral operator with a smooth kernel. The outlined block is the 7 × 8 between the separator's two halves: every entry nonzero, and 5 columns describe it to eight digits. The rank structure was not put there by the elimination. It was in the differential operator before anything was discretised.the separator's 15 unknowns, in the order they sit on the lineoutlined: the block between the two halvesdense, and not independententries nonzero1the block56columns it needs5numbers stored150entries in the square225every entry is nonzeroand six columns describe them
Fig. 19 What eliminating a whole block does to a pattern, which is the alternative to factorising the matrix entire.
What one symmetric permutation does to the storage, on a matrix it does not changeThe same 256 × 256 matrix, its rows and columns renumbered by one permutation and its columns by the same one. The condition number is 24.3948 either way, to eight digits; the Frobenius norm is 6139.964 either way, to twelve. Nothing a norm can see has moved. Clustered, the partition stores 27,008 numbers. Shuffled, the admissibility test finds no admissible pair anywhere — every cluster of a shuffled numbering spans the whole interval, so every q is infinite — and the format degenerates to dense storage exactly. The rule with no test to fail does worse than that: it compresses every off-diagonal block regardless, gets ranks up to 119 out of 128, and stores 118,208 numbers — 1.80 times the matrix it was compressing. A rank-119 factorisation of a 128-column block is a more expensive way to write down the block than the block.numbers stored, 256 × 256clustered, strong27,008clustered, weak24,064the dense matrix65,536shuffled, strong65,536shuffled, weak118,208the same matrix, twiceκ, clustered24κ, shuffled24‖A‖_F, clustered6140‖A‖_F, shuffled6140shuffled weak ⁄ dense1.8the compressibility is in the numberingand the numbering is not in the matrix
Fig. 20 What the orderings cost to compute, which is not free either.
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. 21 An ordering chosen for a machine rather than for a count.
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. 22 The price of the property, and its removal.
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. 23 The matrix underneath, from the constraint field.
The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix43 entriestip eliminated first253 entriestip eliminated last43 entries‖A − LLᵀ‖/‖A‖, tip first1.4·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 253 · sparse factor is 2n − 1 = 43one row swapped to the endnothing numerical chose between them
Fig. 24 And the case where two orderings that tie here part company completely.
Fill under three orderings of a sparse saddle-point matrix, predicted symbolically and then countedA stiffness block of 30 unknowns with 6 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 182, minimum-degree 83, reverse-cuthill-mckee 83: 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.natural182 predicted · 182 countedminimum-degree83 predicted · 83 countedreverse Cuthill–McKee83 predicted · 83 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural182minimum-degree83reverse-cuthill-mckee83predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue
Fig. 25 A larger problem with more constraints, where the border’s floor is a larger share.
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. 26 The property being relied on, at a comfortable regularisation.
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 4 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.0527 and 1.49·10¹⁰. The reason is in the entry count: the 4 entries that moved are the 4 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 moved4off-diagonal0survived at σ = 0.996survived at σ = 0.106 steps0 stepsthe few entries that movedare the ones that dominate
Fig. 27 The sequence the fixed structure is for.
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. 28 And the routes that avoid the question by not factorising K at all.
ν(σ), 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.-4-2.60171-1.203430.1948571.593142.991434.38971024681012shift σν(σ)an answer that is an integershifts2000disagreements0eigenvalues12steps12the marks are a Jacobi decompositionand the staircase never saw one
Fig. 29 The signs the same factorisation produces, put to work.
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. 30 And the conditioning of the systems this ordering is computed for.

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.

Elimination treeFill-inGrowth factorLDLᵀ factorisationMinimum degreeQuasi definite matrixSaddle-point systemsSparse pivotingStatic pivotingSymbolic factorisation