The order that was right last time

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.

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

What the reuse costs #

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

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.

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, and the next section is about what kind of failure it is — because the obvious account, that the perturbed factorisation is a poor preconditioner and the fixed-point iteration’s contraction factor is close to one, turns out not to survive being measured.

Refinement is sublinear, not slow #

Six steps and a factor of 8.8 is a small enough gain to describe as “it stalls”, and the natural mechanism to reach for is a fixed-point iteration whose contraction factor is close to one. That account makes a prediction — a constant ratio between successive errors — and the way to test it is to stop cutting the iteration off at six.

Run to thirty. The backward error goes 4.82·10⁻⁹, 2.32·10⁻⁹, 1.49·10⁻⁹, 1.08·10⁻⁹, 8.27·10⁻¹⁰, 6.63·10⁻¹⁰, 5.46·10⁻¹⁰ … 4.15·10⁻¹¹, and the ratios between consecutive terms go

2.08, 1.56, 1.39, 1.30, 1.25, 1.21, 1.19, 1.17, 1.16, 1.15 … 1.09.

That is not a constant near 1. It is a sequence approaching 1 from above, monotonically, once the first step’s transient is past. There is no contraction factor, because the iteration is not linear.

Nor is there a floor. The error at step thirty is still falling — the run has gained a factor of 116 rather than the six-step run’s 8.8, and nothing about the trace says where it stops. What it is nowhere near is the working precision: 4.15·10⁻¹¹ after five times the budget, against a fresh order’s 9.5·10⁻¹⁷ with no refinement at all.

Sublinear is a worse failure than slow, and it is worse in a specific way. A linear iteration with a rate of 0.99 reaches any tolerance given enough steps, and the number of steps is computable in advance from one measured ratio — so it is a cost, and a cost can be budgeted. A sublinear one never becomes geometric, so there is no ratio to measure and no budget that is the right budget. The first six steps here buy a factor of 8.8, the next six buy 2.4, and the following eighteen buy 5.5. Every extra step is worth less than the one before it, for ever.

That also disposes of the repair somebody would reach for first. Refine harder is a real option when the rate is 0.99 and it is not one here: reaching double precision from 4.8·10⁻⁹ on this trajectory is thousands of solves, each of them the cost of the thing the symbolic reuse was saving. The saving and the repair are the same arithmetic, and the repair is much larger.

Why the shape is sublinear, and what that says about the diagnosis #

The mechanism is worth naming because it is the reason equilibration is a fix rather than an improvement.

Refinement corrects the iterate towards the solution of the system it can see. The perturbed factorisation does not merely approximate A badly — it is the exact factorisation of a different matrix, one whose overwritten row asserts an equation the problem never contained. Correcting against the true residual pulls the iterate towards the true solution along the directions the factorisation can resolve, and the direction it cannot resolve is precisely the one the overwritten row spans. Progress along that direction comes only from the coupling between it and the rest, which weakens as the residual in the other directions is used up — so the gain per step decays, which is what the ratio sequence is.

Stated that way the boundary is sharper than “refinement does not repair it”. Refinement repairs an error in solving the problem; it repairs a change to the problem only through whatever coupling happens to exist, and the rate of that is not a property anybody controls. On this family the coupling is weak enough to give 1.09 per step at thirty. On another it might give 1.5. Neither is a number a code can plan around, and both are avoided entirely by the row scaling that keeps the pivot above the floor in the first place — which needs no steps, no budget and no trace to read.

Where the boundary is #

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

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.

Drawn at the stops either side, the flatness before the cliff is as striking as the cliff.

For six decades the reused order is indistinguishable from a fresh one, and sometimes better. At 1, 2, 3, 4 and 6 decades the kept order’s backward error is 8.36·10⁻¹⁷, 8.54·10⁻¹⁷, 1.01·10⁻¹⁶, 9.5·10⁻¹⁷ and 5.85·10⁻¹⁷ against fresh factorisations at 1.19·10⁻¹⁶, 1.28·10⁻¹⁶, 8.01·10⁻¹⁷ and 1.31·10⁻¹⁶ — the two are interleaved, with no pivot replaced at any of them.

The jump is at the replacement rather than at the scaling. Two further decades past the break cost another factor of six, which is real and is nothing beside the seven orders the first replaced pivot cost. What the kept ordering loses it loses all at once, on the step where it stops being the ordering the matrix would have chosen for itself.

The equilibrated column never moves at any of the seven stops: 6.05·10⁻¹⁷, 5.6·10⁻¹⁷, 5.66·10⁻¹⁷, 8.19·10⁻¹⁷, 5.74·10⁻¹⁷, 5.44·10⁻¹⁷ across one to ten decades, with no pivot replaced anywhere. That is the whole recommendation in one row, and it is worth stating as the shape rather than as a preference: the reused order has a cliff and equilibration removes the cliff rather than moving it. A policy with a threshold in it needs the threshold to be checked at run time; a policy without one does not, and on a sequence of matrices nobody is inspecting, that difference is the entire argument.

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.

One more reading of the boundary table, because it says something about how a code would find this for itself. The transition between six decades and eight is invisible from the successful side: at six decades no pivot is replaced, the backward error is 5.9·10⁻¹⁷, and every diagnostic a solver prints is the diagnostic of a healthy run. At eight, two pivots are replaced and the error is 2.4·10⁻¹⁴ — still a number most codes would accept — and at nine it is 5.5·10⁻¹⁰. So the quantity that moves first is not the accuracy but the count of replaced pivots, which crosses from zero to two while the reported error is still respectable, and which a solver knows exactly and usually does not report. A run that printed it would have had two decades of warning.

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.

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.

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.

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.