Numerical Linear Algebra
where the algebra stops being the arithmetic
Sparsity, and what elimination costs

The order that was right last time

A pivot order computed once and reused across a sequence saves the symbolic phase, and the price is that a pivot which was large may now be small. Replacing it with √u·‖A‖ costs eight orders of backward error and iterative refinement recovers a factor of 8.8 of them. Divide each row by its largest entry first and the same reuse costs nothing at all.
29 min read 24 figures Where the algebra stops being the arithmeticWhose fault is itMeasured, not assertedAssertions that reject

Worth reading first: Structure and stability stop being separable · The units the matrix is measured in · The problem that arrives again.

The field on sparsity measures the trade every sparse elimination makes: the pivot that keeps the factor sparse and the pivot that keeps the elimination stable are different rows, and a threshold buys one with the other. Every number in it is about factorising one matrix.

A code that factorises one matrix is rare. A time-stepping solver, a Newton iteration or a parameter sweep factorises hundreds of matrices with the same sparsity pattern and different values, and the expensive part is not the arithmetic — it is the symbolic phase: choosing the order, allocating the structure, and on a distributed machine deciding which processor holds what. Doing that once and reusing it is worth a large constant factor, and it is what a distributed sparse solver does.

The price is that the order was chosen for a matrix that is no longer the one being factorised.

The backward error of a solve along a sequence, with the pivot order chosen fresh, kept, and kept after equilibrationA 10×10 conflict grid, 100 unknowns, with three of its rows scaled down through 9 decades as the sequence runs — a row that was ordinary becoming a row that is small in the matrix's own units. The sparsity pattern is identical at every member, so a symbolic phase computed once stays valid throughout. Choosing a fresh order each time holds the backward error at the working precision. Keeping the first member's order costs 2 replaced pivots and a backward error of 4.82·10⁻⁹. Keeping the same order after dividing each row by its largest entry costs nothing at all: 6.55·10⁻¹⁷, with no pivot replaced anywhere in the run. The ringed points on the upper curve are the members at which another pivot fell below the floor and was replaced.02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member9.5·10⁻¹⁷kept order, last member4.8·10⁻⁹equilibrated, last member6.5·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in
Fig. 1 The backward error along a twelve-member sequence with one sparsity pattern, under three policies: a fresh order each member, the first member’s order kept, and the same kept order after each row has been divided by its largest entry.

What static pivoting is

The reuse cannot include interchanges. If the code is allowed to swap rows when a pivot turns out small, then the order is not fixed, the structure it allocated is wrong, and everything the symbolic phase bought is given back.

So it does not swap. When a pivot is smaller than √u·‖A‖ it is replaced by that value, with its own sign, and the elimination carries on. √u is 1.05·10⁻⁸ in double precision.

That is a deliberate change to the matrix being factorised. What comes back is the exact factorisation of a matrix that is not A, the difference is bounded by how many pivots were replaced and by how far, and the standard answer to the obvious objection is: iterative refinement repairs it. The residual is computed against the true A, the factorisation is used only to approximate a correction, and a couple of steps bring the backward error to the working precision.

This essay is about a case where that does not happen, and about why.

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. 2 The conflict this whole field is about, from the essay that established it: the sparsest pivot and the largest pivot are different rows.
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 And what a threshold buys, from the essay that priced it. Static pivoting is that dial turned all the way to structure with a repair bolted on.

The sequence

The family is the conflict grid — a five-point Laplacian pushed out of symmetry, which is the object this field’s argument is about — with one addition: three of its rows are scaled down through nine decades as the sequence runs, each starting at its own member.

A row scaling is the smallest honest way to write down a pivot that was large and is now small. It leaves the solution of the system unchanged, leaves the sparsity pattern untouched — so a symbolic phase computed once stays valid, which is the whole premise — and moves nothing except which rows a pivoting rule would want. It is also what happens to a real matrix for real reasons: a row is an equation, and an equation written in different units is a row of a different size.

That the pattern never moves is asserted at every member rather than assumed: zero differing positions against the first member, at member 1, 5, 8 and 11.

The matrix and its LU factors at τ = 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.A288 entries, 64 unknownsL + U986 entries, 8 interchanges→‖PA − LU‖/‖A‖8.8·10⁻¹⁷growth factor1.2fill created698forward error7.7·10⁻¹⁶the 8×8 grid with a drift and one ruinous cornerred marks are entries that were zero
Fig. 4 The pattern in question, from the field that measures what an elimination does to it. Nothing in this essay’s sequence changes a single position of it.
Two condition numbers of one 12×12 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 15 to 3.3·10⁶ while the componentwise condition number stays at 13.2 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.8·10¹⁶ and 1.1·10¹⁶: a matrix whose sensitivity no scaling repairs.0123456110³10⁶10⁹10¹²10¹⁵10¹⁸spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread15κ_∞ at 6 decades3.3·10⁶cond, either end13Hilbert, equilibrated1.3·10¹⁶the solution is the same at every spreadand one of these curves knows it
Fig. 5 And the units, from the essay that established that a matrix does not have one scale. The drift here is entirely a change of units.

What the reuse costs

By the last member, two pivots have fallen below the floor and been replaced.

pivots replaced ‖PA − LU‖/‖A‖ backward error after six refinements fill
a fresh order each member 0 1.6·10⁻¹⁶ 9.5·10⁻¹⁷ — 1,568
the first member’s order, kept 2 1.22·10⁻⁹ 4.8·10⁻⁹ 5.5·10⁻¹⁰ 1,460
the same order, rows equilibrated 0 9.4·10⁻¹⁷ 6.5·10⁻¹⁷ nothing to do 1,458

Two results, and the second is the one worth carrying.

Refinement does not repair the middle row. Six steps take the backward error from 4.8·10⁻⁹ to 5.5·10⁻¹⁰ — a factor of 8.8 — and then it stops falling. Every other repair in this collection’s work on sequences reaches the working precision in one or two steps. This one gets less than one decade and stalls.

And the whole thing evaporates on a row scaling. Divide each row by its largest entry, choose the order once on the equilibrated first member, and reuse it: no pivot is replaced anywhere in the run, the factorisation’s own residual stays at 9.4·10⁻¹⁷, and the backward error is at the working precision at every member with no refinement at all.

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. 6 The refinement history of the damaged run against the one that was not damaged. The flat curve is the content.

Why the floor is the problem

√u·‖A‖ is a threshold in the units of the largest row. A row that has been scaled down by nine decades has entries nine decades below that threshold, so the rule does not protect it from a small pivot — it overwrites it, with a number that has nothing to do with the equation that row represents.

That is not a rounding and it is not a perturbation of the size the analysis assumes. It is a substitution of one number for another, and the analysis behind static pivoting assumes the replacement is of order √u relative to the matrix, which it is — relative to the matrix. Relative to the row, it is a factor of 10⁹.

The refinement then fails for a reason worth stating precisely. Refinement is a fixed-point iteration whose contraction factor is how good the perturbed factorisation is as a preconditioner for the true matrix. Once a pivot has been replaced by something nine orders larger than the row it came from, that factorisation is not a good preconditioner, the contraction factor is close to one, and the iteration crawls. The failure is not that refinement is slow here. It is that its rate is nearly 1.

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. 7 The same fact one rung down, from the essay on a pivot that reads the units: a pivoting rule compares numbers, and numbers in different units are not comparable.
Iterative refinement from a 24-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 8 What refinement does when its rate is good, from the field that measured it — one or two steps to the working precision, which is what makes the stall in this essay a measurement rather than an impatience.

Where the boundary is

The damage arrives when the rows fall past the floor and not before, and the sweep says exactly where:

rows scaled down by pivots replaced backward error after refinement
0 decades 0 6.0·10⁻¹⁷
2 0 8.5·10⁻¹⁷
4 0 9.5·10⁻¹⁷
6 0 5.9·10⁻¹⁷
8 2 2.4·10⁻¹⁴
9 2 5.5·10⁻¹⁰
12 2 7.5·10⁻¹⁰

Nothing at all happens for six decades. The transition is between six and eight, which is where the row’s entries cross √u·‖A‖, and past it the refinement is not recovering. And the equilibrated column of the same sweep replaces no pivot at any setting, at any number of decades.

So the boundary is not a property of the sequence, the drift rate or the conditioning. It is the point at which a row becomes small in the units the floor is written in, and moving the units moves it out of the picture entirely.

The backward error of a solve along a sequence, with the pivot order chosen fresh, kept, and kept after equilibrationA 10×10 conflict grid, 100 unknowns, with three of its rows scaled down through 4 decades as the sequence runs — a row that was ordinary becoming a row that is small in the matrix's own units. The sparsity pattern is identical at every member, so a symbolic phase computed once stays valid throughout. Choosing a fresh order each time holds the backward error at the working precision. Keeping the first member's order costs 0 replaced pivots and a backward error of 9.55·10⁻¹⁷. Keeping the same order after dividing each row by its largest entry costs nothing at all: 7.21·10⁻¹⁷, with no pivot replaced anywhere in the run. The ringed points on the upper curve are the members at which another pivot fell below the floor and was replaced.02468101210⁻¹⁸10⁻¹⁶member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member1.1·10⁻¹⁶kept order, last member9.5·10⁻¹⁷equilibrated, last member7.2·10⁻¹⁷pivots replaced, kept0pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in
Fig. 9 Four decades of row scaling, where nothing happens: every policy sits at the working precision and the reuse is free.
The backward error of a solve along a sequence, with the pivot order chosen fresh, kept, and kept after equilibrationA 10×10 conflict grid, 100 unknowns, with three of its rows scaled down through 12 decades as the sequence runs — a row that was ordinary becoming a row that is small in the matrix's own units. The sparsity pattern is identical at every member, so a symbolic phase computed once stays valid throughout. Choosing a fresh order each time holds the backward error at the working precision. Keeping the first member's order costs 2 replaced pivots and a backward error of 5.27·10⁻⁹. Keeping the same order after dividing each row by its largest entry costs nothing at all: 5.33·10⁻¹⁷, with no pivot replaced anywhere in the run. The ringed points on the upper curve are the members at which another pivot fell below the floor and was replaced.02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member7.9·10⁻¹⁷kept order, last member5.3·10⁻⁹equilibrated, last member5.3·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in
Fig. 10 And twelve, where the damage has saturated — two pivots replaced, and replacing them further does not make it worse.
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 9.55·10⁻¹⁷ to 9.55·10⁻¹⁷, a factor of 1.0, and stop. The equilibrated run has nothing to repair: it starts at 7.21·10⁻¹⁷.012310⁻¹⁷10⁻¹⁵refinement stepsbackward error of the solvethe kept order, in the units as giventhe same order, rows equilibratedthe repair that does not repairbefore refinement9.5·10⁻¹⁷after six steps9.5·10⁻¹⁷the factor gained1equilibrated, before7.2·10⁻¹⁷pivots replaced0the residual is recomputedand it cannot see what was done to the matrix
Fig. 11 The refinement figure below the boundary, where there is one point and nothing to repair.

The reused order is also the sparser one

A number falls out of the same run that has nothing to do with stability and is worth having.

The equilibrated reused order leaves 1,458 entries in the factors. A freshly chosen order at the last member leaves 1,568. The order computed once is not merely usable, it is 7 per cent sparser than the one threshold pivoting picks for that member.

The reason is that threshold pivoting on a later member is choosing under a constraint the first member did not impose: it has to reject rows whose pivot candidate is below τ times the column maximum, and by the last member some of those rejections are forced by the row scalings rather than by the structure. A choice made when the matrix was well scaled is a choice made with more options available.

That is an argument for reuse this field did not have, and it points the same way as the cost argument — provided the equilibration is done.

Fill growth under reverse Cuthill–McKee: the factor rises as n^1.41A 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.41matrix, n^1.05grid Laplacians from 5×5 to 12×12fitted, not quoted
Fig. 12 Fill as the quantity a sparse elimination is judged on, from the essay that measured it. Seven per cent of it is not nothing at the sizes this matters at.
Nonzeros in the Cholesky factor of the 8×8 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.natural519reverse Cuthill–McKee428minimum degree359nested dissection434matrix: 176 entries · dense factor: 2080bandwidth 8 · 2.95× the matrixbandwidth 8 · 2.43× the matrixbandwidth 49 · 2.04× the matrixbandwidth 48 · 2.47× the matrixn = 64, five-point stencilevery ordering fills in; none avoids it
Fig. 13 And the orderings compared, from the essay that ranked them. A reused order is a fifth entry on that chart whose position depends on the units.

What equilibration is and is not

Row equilibration — divide each row by its largest entry — is usually described as a preliminary, a tidying step done before the interesting work. Everything above says it is not.

It is worth being precise about what it does and does not change.

It does not change the solution. Scaling a row scales the corresponding entry of the right-hand side, and the system is the same system. Nothing about the answer moves.

It does not change the conditioning in general. The collection has an essay on a condition number scaling cannot move, and it applies here: equilibration improves the conditioning of some matrices and not others, and the improvement is bounded by a factor that has nothing to do with this essay.

What it changes is which numbers are comparable. A pivoting rule compares candidates against each other and a perturbation floor compares them against a norm. Both are comparisons, both are meaningless across rows in different units, and equilibration is the statement that there is no such thing as a small row — only a small entry relative to its own row.

So the honest summary is that static pivoting is not a technique with a stability caveat. It is a technique that requires an equilibration, in the strong sense that without one it is not merely less accurate but unrepairable, and with one it is free. That is a much sharper statement than “scaling helps”, and the six orders of magnitude between the two policies at the last member is what makes it one.

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. 14 The comparison a pivoting rule makes, at a milder scaling. Nothing about the rule is wrong; what is wrong is asking it to compare numbers from different equations.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 15 And the bound that is written per entry rather than per norm, from the field that measured it — the form in which the difference between the two policies is visible at all.
Backward error with and without the structure, for Levinson and for elimination, n = 20Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 4.2·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 3.059·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.3·10⁻¹⁷Levinson, Toeplitz only10⁻¹⁰the ratio between them3.1·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 16 The quantity the whole comparison is in, once more, at the size this essay’s grid uses.

The one repair in this collection that does not work

It is worth saying plainly what this essay is doing in a run of essays about outer loops forgiving their inner objects.

Every other one has an outer loop that recomputes the residual from the matrix, and that recomputation is what makes a stale or approximate inner object survivable. A Newton step solved to two digits, a factorisation four members old, a preconditioner built for a different operator — all forgiven, because the next step measures against the real problem.

Here the outer loop does the same thing and it does not work. The residual is recomputed against the true A, in double, from the current iterate. What defeats it is that the perturbation is not an error in solving the problem, it is a change to the problem — and it is a change in a direction the residual is badly conditioned to see, because the row that was overwritten is the row whose entries are nine orders below everything else in the residual vector.

So the field’s rule has a boundary and this is it. An outer loop forgives an inner object that approximated the answer. It does not forgive an inner object that answered a different question.

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. 17 The repair working, from the field that established it — for comparison with the figure above where it does not.
A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 0.0372 from a root that is known exactly by construction. The same linear system is solved to relative residuals from 0.3 down to 10⁻¹⁴, at 109 and 1126 conjugate gradient iterations, and the resulting point is 0.005319 and 0.002497 from the root. The curve is flat below about 10⁻³: the linearisation is wrong at second order, so the step cannot land closer than the square of the distance it started at — 0.001383 — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻⁴10⁻³10⁻²10⁻¹1inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 0.00138distance before the step, 0.03721093547231126the number by each point is the iterations it costeleven decades, one landing placedistance before the step0.037its square0.0014where η = 10⁻³ lands0.0025where η = 10⁻¹⁴ lands0.0025iterations for the first354iterations for the second1126the accuracy that is thrown awaymeasured against a root that is known
Fig. 18 And the forgiveness measured at its most extreme: eleven decades of inner error, one landing place.
Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 19 The quantity all of this is in, from the essay that established the collection’s rule.

Why the symbolic phase is worth this much trouble

The essay is about an accuracy failure, so it is worth restating what is being bought, because the answer is not arithmetic and the size of it is why anybody accepts the risk.

A sparse factorisation has two halves. The numeric half is the arithmetic — the multiplications and subtractions that produce the factors — and it is the half a flop count sees. The symbolic half decides where the fill will go, allocates the structure to hold it, and on a distributed machine decides which processor owns which part of the elimination tree. It computes no entries at all.

For a matrix whose values change and whose pattern does not, the symbolic half is identical at every member and the numeric half is not. Reusing it turns a sequence of factorisations into a sequence of numeric passes over a structure allocated once, and on a distributed solver it also fixes the communication pattern — which is a second saving the flop count cannot see and which this collection’s cost field has spent a run establishing is the one that matters at scale.

That is what static pivoting is protecting, and it is why replacing a pivot with a made-up number is considered an acceptable trade rather than an obvious mistake. The measurement above does not overturn that trade. It says the trade has a precondition, that the precondition is a row scaling, and that without it the repair the trade depends on does not arrive.

What is asserted

That the pattern never moves — zero differing positions at four members spread through the sequence — which is what makes the symbolic reuse legal in the first place.

That the kept order is destroyed and not repaired: two pivots replaced, a backward error seven orders above a fresh order’s, refinement gaining a factor of 8.8 and not reaching the working precision.

That equilibration makes the same reuse free: no pivot replaced at any member, the first solve already at the working precision, and six orders of magnitude between the two policies on a row scaling alone.

And that the boundary is the floor: nothing at four decades, everything at nine, with the equilibrated column flat throughout.

The refusals are the three shapes this could be misread as. The damage claimed where nothing was rescaled; the refinement failure claimed at the member the order was chosen for; and equilibration claimed to fix the fill as well as the stability, which it does not — the fill is a separate measurement and it happens to go the same way.

Growth factor under partial pivoting to n = 20: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.048121620110¹10²10³10⁴10⁵10⁶10⁷10⁸matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 20: bound 5.2·10⁵, worst 2.8
Fig. 20 The growth factor, which is what an elimination’s stability is usually judged by and which the replaced pivots send to 10⁶ here without anything else going visibly wrong.
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. 21 And the symbolic phase, whose reuse is the entire prize this essay is about.
What 12 members of a sequence cost, by what changes between themEvery bar is 12 solves of an 120×120 system, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n². When nothing changes the answer is reused and the cost is one factorisation. When only the right-hand side changes, one factorisation serves every member — 16.5 per cent of the independent cost, with every solve still at a relative residual of 1.54·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 4 factorisations for 12 members. When everything changes, nothing carries.nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried
Fig. 22 Where it sits in the ledger of what a sequence licenses.
Fill against growth as the pivot threshold moves, on the 8×8 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.1·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1875one knob, two measurements, opposite directionsand the default is most of both
Fig. 23 The trade a threshold makes, from the essay that measured it. Static pivoting is the limit of that curve, and this essay is the condition under which the limit is safe to take.
Where the same pivot rule lands on Kahan's matrix, with and without a tie tolerance, at c = 0.5Two curves of |r_nn| ÷ σ_min against n. Every trailing column norm of this matrix is exactly equal at every step, so a rule that interchanges on any strict improvement is comparing two rounding errors: it makes 0, 9, 14, 6 interchanges at n = 20, 30, 40, 50 and lands 4.5·10⁶ times nearer the truth at the largest size. It lands in the SAME place at the three largest sizes while the other rule's answer grows by orders of magnitude across them, which is what a floor set by rounding looks like rather than a rule that is finding something.18283848110²10⁴10⁶10⁸size of the matrix|r_nn| ÷ σ_minequal treated as equalstrict comparisona choice between equal numbersties found, at n = 5049interchanges, tie-tolerant0interchanges, strict6ratio between the two4.5·10⁶the improvement is real and is not a repairit is the ulps, and they are not reproducible
Fig. 24 And how a tie between two candidate pivots is broken, which is the decision being frozen when an order is reused.

Where this goes next

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. forward error, relative to a solution of exactly (1, 1)no pivoting · as given1 0 interchangesno pivoting · rows scaled1 0 interchangespartial · as given0 1 interchangepartial · rows scaled1 0 interchangesscaled partial · as given0 1 interchangescaled partial · rows scaled0 1 interchangecomplete · as given0 1 interchangecomplete · rows scaled0 1 interchangethe same problem twicepartial, as given10⁻¹⁸partial, rows scaled1its relative residual10⁻¹⁷complete, rows scaled10⁻¹⁸the two systems have the same solutionand one pivot rule cannot see it The pivot that reads the units Partial pivoting compares the entries of a column and takes the largest. Those entries carry units, so the comparison depends on them — and there is a row scaling, on the standard two-by-two that pivoting exists to fix, which makes partial pivoting perform the identical catastrophic elimination it was introduced to prevent, with no interchange at all. nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried The problem that arrives again A hundred and thirty essays have solved a system once and measured how wrong the answer was. Almost no computation is shaped like that. A solve is one step of an outer loop, its answer is an input rather than a deliverable, and four quantities this collection has treated as accuracy requirements turn out to be assets with a shelf life.

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.

Backward errorEquilibrationFill-inGrowth factorIterative refinementSparse luStatic pivotingSymbolic factorisationThreshold pivoting