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

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.