Structure and stability stop being separable

The sparsity field is stated entirely in the symmetric positive definite case, and that is not an accident of what happened to get written first. It is the only case where the two halves of a sparse direct solve come apart cleanly, and both of the field’s central claims depend on their coming apart.

Ordering is free. The order decides the memory compares four elimination orders on the same matrix and finds factors of 1,739, 1,413, 1,354 and 1,026 entries. All four return the same answer to 3.1·10⁻¹⁶. That statement is only available because a symmetric positive definite matrix needs no pivoting for stability, so the order can be chosen purely to minimise fill and nothing numerical objects.

The symbolic phase is exact. The factor is not sparse computes the fill twice — from the elimination graph with no arithmetic at all, and by factorising and counting — and the two agree as integers on all five test structures. That is available because the order is fixed before any arithmetic happens, so there is nothing left for the numbers to decide.

Neither survives a general matrix. This essay is what happens to both.

The conflict has a size, it is a function of one entry, and at one end of the slider it is not there at all.

Across corner values of 10⁻¹², 10⁻⁹, 10⁻⁶, 10⁻³ and 1 the sparsest choice’s growth reads 1.9·10¹¹, 1.9·10⁸, 1.9·10⁵, 188 and 1 — which is 0.19 divided by the corner, exactly, over twelve decades — and its forward error reads 3·10⁻⁵, 5.3·10⁻⁹, 3·10⁻¹¹, 1.1·10⁻¹⁵ and 1.6·10⁻¹⁶. The safe choice’s growth is 1.19 at every one of the five and its error stays between 5·10⁻¹⁶ and 2·10⁻¹⁵.

So the conflict is not a property of the matrix’s sparsity pattern, which does not change along the slider. It is a property of one number, and the price of ignoring it is that number’s reciprocal.

And at a corner of 1 the safe choice is worse on both counts. The sparsest pivot gives 233 entries, a growth of 1 and an error of 1.6·10⁻¹⁶; partial pivoting gives 241 entries, a growth of 1.19 and an error of 8.9·10⁻¹⁶. Eight extra entries and five times the error, for a safeguard against a danger that is not present.

That is the honest shape of the trade this field is built on. Partial pivoting is not insurance with a premium of a few entries; it is a decision that is worth 10¹¹ at one end of a single entry’s range and worth less than nothing at the other, and nothing about the sparsity pattern says which end a matrix is at.

The matrix, and why it is built the way it is #

The demonstration needs a matrix on which the two criteria disagree, and building one honestly means building one anybody would recognise rather than one designed to produce an answer.

The starting point is the five-point grid operator — the same one the whole sparsity field counts fill on. Three things are done to it.

A drift, which scales the two forward couplings down and the two backward ones up. That is what an upwinded convection term does to the stencil, and it is what makes the matrix unsymmetric, so Cholesky is unavailable and the question this essay is about exists at all.

One ruinous diagonal entry, at node 0. A grid corner has two neighbours where an interior node has four, so it is the unique cheapest variable to eliminate; giving it a diagonal of 10⁻¹² makes it also the worst possible pivot. The conflict is therefore at the very first step, where it can be drawn.

Nothing else. The rest of the stencil is left alone, so the structure is a structure and the counts below are counts of something.

The corner’s degree is 3 — two neighbours and itself — and that is the minimum over the whole matrix, which the build asserts rather than assumes. So a fill-reducing rule of any kind picks it first, for a reason that is entirely correct on its own terms.

What eliminating it does #

The measurement is the essay, and it did not come out where the phase plan expected.

Growth. The multiplier at the first step is about 1/10⁻¹², so entries of order one become entries of order 10¹². The growth factor — the largest entry appearing in any reduced matrix, divided by the largest entry of the original — comes out at 1.9·10¹¹, and it is the same number at every grid size, because it is set by that one division.

Fill. The unpivoted elimination produces 122 entries of fill against partial pivoting’s 137. The ruinous choice really is the cheaper one, by 11%, which is what makes this a dilemma rather than a decision.

The residual. ‖PA − LU‖/‖A‖ = 3.8·10⁻¹⁷. That is not merely acceptable. It is better than every pivoted run in the sweep — partial pivoting gives 5.4·10⁻¹⁷ on the same matrix — and it is eleven decades below what the growth factor’s bound permits. The factorisation is impeccable.

The answer. Relative error 3·10⁻⁵, against 4.8·10⁻¹⁶ for the pivoted route. Four correct digits.

That combination was the surprise. The expectation, and the shape the swap that is not optional establishes, was that skipping the pivot would damage the factorisation and that the damage would show in the residual. It does not show in the residual at all. The factors reproduce the matrix and the triangular solves destroy the answer — L carries entries of 10¹², U carries entries of 10¹², and the cancellation between them in forward and back substitution takes fifteen digits with it.

So this is the site’s oldest argument, arriving in the sparsity field by a route it had not taken before. A small residual is not a small error makes the point about a solve; here it is made about a factorisation, where the residual is not merely small but the smallest number in the table.

Why the symbolic phase stops being exact #

The second casualty is quieter and has a larger practical consequence.

In the symmetric case, the elimination order is decided before any arithmetic runs, so the fill can be computed from the graph. That computation is not an estimate — it is the answer, and on a matrix with an ordinary diagonal the site checks it as an integer identity: symbolic 233, numeric 233.

Under pivoting the row interchanges are decided by the values, and the values are not known until the elimination is under way. So the structure of the factor is not a function of the structure of the matrix any more. On the 5×5 grid the symbolic count is 233 and the pivoted factorisation produces 242, and the assertion that the first predicts the second is fed to the build and refused.

What is left is a bound, and it is the subject of what the symbolic phase can only bound. The short version: the union over every possible row permutation is computable in advance, it is the structure of the Cholesky factor of AᵀA, and it is loose by between 1.4 and 1.7 times on the grids drawn here.

The Markowitz idea, and why choosing is not the answer #

The obvious response to a conflict between two criteria is to pick one. Both choices are bad and the badness is asymmetric.

Pick stability and ignore fill — partial pivoting, exactly as the dense case does it — and the answer is right and the factor is 11% larger here, which sounds mild and is not. On the grids in this field the gap between a good ordering and no ordering is a factor of 1.7 and rises with size; a solver that discards the fill argument entirely gives that up, and on a large three-dimensional problem giving it up is the difference between fitting in memory and not.

Pick fill and ignore stability and the answer is the four-digit one above, or worse.

What production codes do instead is refuse to rank the criteria and constrain one with the other. Markowitz’s rule scores a candidate pivot by (r − 1)(c − 1), where r is the number of remaining entries in its row and c in its column: an estimate of how much fill eliminating it will create. Threshold pivoting then says a pivot is acceptable if it is at least τ times the largest candidate in its column, and among the acceptable ones the cheapest by Markowitz is taken.

τ = 1 is partial pivoting; τ = 0 is pure Markowitz; the libraries ship τ = 0.1. Turning that knob and watching both quantities move is a threshold between fill and growth, and this essay’s job is only to establish that the knob has to exist.

One detail about the implementation here, since it affects what the counts mean. Real codes permute both rows and columns; this site permutes rows only, and breaks ties on the row count with the column count held fixed. That changes the constants and not the argument — the conflict, the growth, the loss of symbolic exactness and the shape of the trade are all unchanged — and saying so is cheaper than implying a full Markowitz search was run.

The end of the range that is not a strategy #

τ = 0 is not the loose end of the knob. It is off the end of it, and the measurement is worth recording because the intuition points the other way.

With no stability constraint at all, the rule will accept as a pivot any nonzero entry — including an entry created by fill and then nearly cancelled, of size 10⁻¹⁶. The factorisation does not merely lose accuracy. It comes apart: ‖PA − LU‖/‖A‖ = 2.5·10¹¹, and the solve returns infinities.

So the sweep in the next two essays starts at 10⁻³ and not at zero, and every real code has a floor for the same reason. “Minimise fill subject to nothing” is not a strategy with a bad constant; it is a strategy with no lower bound on how small a number it will divide by.

What a real solver does with this, and the phase that goes missing #

The symmetric case gives a sparse solver a clean division of labour, and it is worth writing down because the general case breaks it in a specific place.

Analyse. Choose an ordering, compute the fill symbolically, allocate. No arithmetic on values. Factorise. Do the arithmetic into the space that was allocated. Solve. Two triangular systems, per right-hand side.

The value of that division is not tidiness. It is that the first phase depends only on the pattern, so a code stepping through time with a matrix whose structure never changes does the analysis once and repeats only the second and third phases. That is the standing advice at the end of the order decides the memory, and it is worth a great deal — the analysis is often the most expensive phase.

Under pivoting, the first phase cannot finish. It can choose an ordering and it can bound the allocation, and it cannot know what the factorisation will do. Codes resolve that in one of two ways and the choice is visible in their names.

Dynamic pivoting admits it and interleaves the phases: choose pivots as the numbers arrive, allocate as fill appears, and accept that the data structures grow during the factorisation. This is what threshold pivoting above describes, and it costs the reuse — a second matrix with the same pattern and different values may pivot differently, so the analysis does not transfer.

Static pivoting refuses to. Choose the ordering and the pivot sequence up front from the pattern and from the magnitudes as they stand, factorise into exactly that structure whatever the numbers turn out to be, and if a pivot comes out unusably small, replace it — with something like ε·‖A‖ — rather than swapping rows. That produces a factorisation of a matrix that is not quite the one that was asked about, which sounds indefensible and is not: the perturbation is a backward error of known size, and the answer is then repaired by iterative refinement.

The second is worth dwelling on, because it is this site’s machinery reappearing in a place the mixed-precision field did not anticipate. The argument is exactly buying the accuracy back’s: a cheap factorisation of a nearby problem, plus refinement against accurately computed residuals, converges to the answer of the original problem provided κ times the perturbation stays under one. There the cheapness came from low precision; here it comes from a fixed pivot order. Same identity, same condition, and the same thing goes wrong when the condition fails.

What static pivoting actually repairs #

Every clause of that paragraph is checkable on the matrix this essay is about, and one of them fails. The pivot order is fixed at the natural one, ε is swept, and the answer is measured against the exact solution rather than against a better computation. κ here is 25.6 — the conflict grid is not ill conditioned, which is worth holding on to, because it means nothing below can be blamed on the condition number.

The accounting works at the bottom of the table and misses by ten orders at the top. For ε ≥ 10⁻⁹ the predicted error and the measured one agree within a factor of two, which is as well as a bound of that kind ever does. At ε = 10⁻¹⁶ it predicts 9.6·10⁻¹⁶ and the answer is wrong at the fifth digit. The residual is not the backward error of the solution, and the loss the growth factor causes passes through both columns without appearing in either — which is the essay’s own subject, one level up, applied to the argument that justifies the technique rather than to the technique.

Two things fall out of the table that the usual account does not mention.

The error is not monotone in ε. A larger deliberate perturbation gives a more accurate answer, down to ε ≈ 10⁻⁹, and gets worse again after it. The shape is not an accident: the growth factor falls exactly as 1/ε — 4.7·10¹⁰, 4.7·10⁷, 4.7·10⁴, 4.7·10², a decade of growth per three decades of ε, checked as an identity rather than read off — while the perturbation’s own contribution rises as ε. Two terms, one falling and one rising, cross near √u, and √u is where the minimum sits. So the replacement threshold is not a safety floor to be set as low as possible; it is a parameter with an optimum, and setting it low is the wrong end.

And at ε = 10⁻¹⁶ nothing is replaced at all. The pivot is 10⁻¹², which is larger, so the code takes it and the factorisation is the ruinous one from the top of this essay. Refinement recovers the answer to the rounding level in two steps anyway — as it does at every other row of the table, in between one and four steps. Which means that on this matrix what refinement repairs is the growth, not the replacement: the mechanism the argument is about is doing none of the work, and the mechanism doing the work is the one the argument treats as a finishing touch.

That does not make static pivoting wrong, and it is worth being clear about which part survives. The technique delivers what it promises — a reusable analysis phase and an answer good to the rounding level. What does not survive is the reason usually given, which prices the perturbation and is silent about the growth, and which therefore recommends the smallest ε a code can manage when the measurement recommends the largest one the condition number will bear.

That parallel is the useful way to hold the whole essay. Both halves of a sparse solve are trading a perturbation for a cost, and the currency is the same one the error field has been counting since the foundation: a backward error, multiplied by a condition number, giving a forward error. Pivoting less is a larger backward error bought with less fill. Factorising in single precision is a larger backward error bought with faster arithmetic. Static pivoting is a larger backward error bought with a reusable analysis phase. In every case the question is whether the product with κ is small, and in every case the honest report is both numbers rather than one.

What makes the sparse case harder than the mixed-precision one is that the size of the perturbation is not known in advance. u is a property of the format and is written on the box; the growth factor of a threshold-pivoted elimination is a property of the matrix and is only known once the factorisation has run. So the check has to be done afterwards, on the residual — which is why every serious sparse solver returns one, and why a code that reports only its answer is a code that has withheld the number this whole field is stated in.

What is asserted here #

The ruinous pivot is the sparsest row, checked against the minimum degree over the whole matrix rather than asserted about the corner.

And the ruinous elimination produces less fill, 122 against 137, so the fill argument’s recommendation is genuinely the cheaper one and the dilemma is real.

The unpivoted factorisation reproduces the matrix, to 3.8·10⁻¹⁷, and within a factor of ten of what pivoting achieves — which is the assertion that stops the essay from being a story about a factorisation that went wrong.

The pivoted answer is right at every position of the slider, to better than 10⁻¹¹, and the unpivoted one is wrong by a factor of at least 10⁵ at the small end of it.

The brackets exclude the middle. The growth factor here is almost exactly 1.88 divided by the corner entry, so the error ratio slides continuously from 6·10¹⁰ at the left of the slider to 1 at the right and passes through everything between. At 10⁻⁶·⁶ the unpivoted answer is wrong by 5·10⁴ times the pivoted one, which is neither a failure nor a success. Asserting a bracket over the whole range would be asserting a behaviour that does not exist, so the assertions hold at the two ends and say nothing in the middle.

The refusals: that the symbolic phase predicts the pivoted fill, which is refused; and that a factorisation reproducing its matrix to 3.8·10⁻¹⁷ gives the right answer, which is refused on the run that does exactly that.

Where the structure is all there is #

A sparse elimination trades fill against stability because it can see both. An operator supplied as a subroutine offers neither, and the trade disappears along with the ability to make it.