A threshold between fill and growth

Every slider on this site so far improves one thing. Drag the mantissa up and the error falls. Drag the condition number down and the answer gets better. Drag the oversampling up and the randomised approximation improves. The knobs are all of the form more of this is better and costs more, and the interesting question is where the diminishing returns set in.

This one is not like that. Threshold pivoting has a single parameter τ, and moving it in either direction makes one measured quantity worse as it makes another better. There is no setting at which both are best, and the choice is therefore a genuine trade rather than a budget.

What the knob does #

At step k of the elimination, the rule looks at every remaining row that has a nonzero in column k. The largest of those candidates is m. A row is acceptable if its candidate is at least τ·m, and among the acceptable rows the one with the fewest remaining entries is taken.

τ = 1 makes only the largest acceptable, which is partial pivoting exactly as the dense case does it. τ = 0 makes everything acceptable, so the choice is made on structure alone. Anything between accepts a pivot that is somewhat smaller than the best available, in exchange for a row that couples to fewer other variables.

The bookkeeping is what makes it a real algorithm rather than a rule: the row counts change every step, since eliminating a variable adds entries to the rows that touched it, so the “fewest remaining entries” test is against a graph that is being modified as the elimination proceeds.

The two measurements, and their units #

Putting two quantities on one axis needs an argument, and here it is a short one: both are ratios.

Fill is reported as the entries in L and U together, divided by the entries in the matrix. On the 6×6 grid the matrix has 156 nonzeros; the factors have between 328 and 416 depending on τ, so the ratio runs from 2.1 to 2.7.

Growth is the largest entry appearing in any reduced matrix, divided by the largest entry in the original — the quantity the bound that is never attained is about, and the one every backward-error bound for Gaussian elimination is stated in. It runs from 1,330 down to 1.2 across the same range.

Both are dimensionless, so the axis is honest. And both are measurements rather than estimates: the growth is tracked through every step of every elimination, and the fill is counted from the factors that came out.

One detail about the growth definition, because the first version of the library got it wrong in a way that produced a plausible number. The growth factor is a statement about the reduced matrices, not about the multipliers. Including the stored L entries reported 3.8·10¹¹ where the largest entry any reduced matrix holds is 1.9·10¹¹ — a different quantity wearing the same name, since the multiplier 1/pivot is large by construction the moment a small pivot is accepted and says nothing about whether the elimination lost anything.

The numbers #

On the 6×6 grid, 36 unknowns, 156 nonzeros:

Three things in that table are worth stopping on.

The trade is monotone in trend and not step by step. Between τ = 0.01 and τ = 0.03 the growth falls from 224 to 31 while the fill rises by two entries — one extra swap improved both. And between 0.03 and 0.1 the growth goes up, from 31 to 38, while the fill also rises. The pivot sequence changes discretely, so a smooth trade-off curve is a summary rather than a law, and the site asserts the shape it actually has: fill rises with the threshold allowing for the discreteness of a pivot choice, rather than pairwise monotonicity, which would be asserting something false about a real trade.

The interchange count falls as τ rises, which is backwards from what the rule suggests — a stricter threshold should force more swapping. It does not, because with τ = 1 the largest entry is often already on the diagonal after the previous step’s elimination, and with a loose threshold the rule keeps hunting for a sparser row and moving things around. Twenty-eight interchanges at τ = 0.001 against eight at τ = 1.

Every row of the table is a usable factorisation. ‖PA − LU‖/‖A‖ is between 8·10⁻¹⁷ and 2.7·10⁻¹⁵ throughout, and every answer is right to at least 3.8·10⁻¹⁴. The trade is between memory and a margin of safety, not between memory and correctness — which is why the choice is a matter of judgement rather than a matter of the answer being wrong.

Why 0.1 #

The default is not a compromise between the two ends. It is most of both, and the measurement says so plainly.

Against the strictest setting, τ = 0.1 gives up nothing in correctness — the answer is right to 1.4·10⁻¹⁵ — and buys back 44 of the 88 entries that partial pivoting added. Against the loosest, it gives up 44 entries and buys back a factor of 35 in growth.

The reason it lands where it does is that the two curves have different shapes. Growth is steep at the left of the plot and nearly flat at the right; fill is nearly flat at the left and steep at the right. So there is a broad middle where both are close to their best, and the default sits in it. A knob whose two costs were both steep everywhere would have no such place and the choice would be much harder.

What the number is not is optimal in any provable sense. There is no theorem that says 0.1; there is a body of practice that says growth beyond a couple of orders of magnitude starts to matter and that fill beyond about a third above the minimum starts to hurt, and 0.1 satisfies both on the matrices people solve. It is a number chosen by measurement, and the measurement is the kind above.

The bound, and what it is really saying #

The reason growth is the quantity worth trading against fill — rather than, say, the smallest pivot, or the number of interchanges — is that it is the quantity the backward error bound is stated in.

Gaussian elimination with any pivoting strategy produces factors satisfying

$$\|PA - LU\| \lesssim c\,n\,\rho\,u\,\|A\|$$

where ρ is the growth factor and u the unit roundoff. So a growth factor of 38 costs a factor of 38 in the backward error bound, and a growth factor of 1,330 costs 1,330 — which at double precision means 3·10⁻¹³ rather than 2·10⁻¹⁶ and is still, for most purposes, nothing at all.

That is the real reason a loose threshold is usually safe. The bound has three decades of room in it before anything a user notices happens, and threshold pivoting spends one or two of them to save a quarter of the memory.

It is also why the measured residuals in the table are far below the bound. At τ = 0.001 the bound permits 3·10⁻¹³ and the measured residual is 2.7·10⁻¹⁵ — the bound is loose by a factor of a hundred, in the direction it is always loose, for the reason the bound that is never attained gives: worst-case growth requires a conspiracy among the entries, and ordinary matrices do not conspire.

Which leaves a practical rule that this site is willing to state. Do not tune τ. Measure the residual afterwards. The threshold controls a bound on a quantity that can be computed directly once the factorisation exists, and computing it costs one matrix-vector product.

What makes the slack disappear #

The rule above — do not tune it, measure the residual — is good advice on ordinary matrices and it is worth knowing where it stops being good advice, because the failure is not gradual.

The three decades of slack come from the growth factor being small in practice, and the bound that is never attained measures exactly how small: a median of 3.23 and a worst of 4.81 over thirty random 40×40 matrices, against a bound of 5.5·10¹¹. But that essay also draws Wilkinson’s matrix, which attains 2ⁿ⁻¹ exactly, at every size tested, with zero row interchanges — partial pivoting looks at it and finds nothing to swap. So the slack is a property of the matrices people usually have and not a property of the algorithm.

Two things make it thinner in the sparse case than in the dense one, and both are worth naming since neither is obvious.

A sparse elimination has more steps at which a small pivot is the only option. In a dense matrix every column has n candidates and the largest of them is usually of a reasonable size. In a sparse matrix a column may have three, and if all three are small then partial pivoting’s “largest available” is a small number and no threshold changes that. The growth then comes from the matrix rather than from the strategy, and no setting of τ helps.

And the fill argument keeps pushing towards the small ones. The Markowitz rule is actively seeking the sparsest row, which on a matrix with any structure at all correlates with the rows that have the fewest and smallest entries. So the two criteria are not merely in tension; on some structures they are anti-correlated, and the loose end of the slider is worse than a random choice would be.

The practical consequence is the one every sparse solver’s documentation states and few users read: a sparse direct solve should be followed by a residual check, and by refinement if the residual is not small. Not because the strategy is unsound, but because the quantity that bounds its error is one the strategy cannot control and the user cannot predict.

Two currencies, one identity #

It is worth putting this essay’s trade beside the other two the site has drawn, because they are the same trade in different units and the identity underneath is the one the whole site is organised on.

Forward error is bounded by the condition number times the backward error, and the condition number belongs to the problem. So every choice an algorithm makes is a choice about the backward error, and every saving is bought by making it larger.

Buying the accuracy back buys speed with precision: factorise in 24 bits, take a backward error of 10⁻⁷ instead of 10⁻¹⁶, and repair it with refinement provided κ·u stays under one. This essay buys memory with growth: accept a pivot ten times smaller than the best, take a backward error thirty-eight times larger, and repair nothing because it did not need repairing. The order decides the memory buys memory with nothing at all, because in the symmetric positive definite case the backward error does not move.

Three trades, one identity, three different currencies — and the third one is free, which is why the symmetric case is where the sparsity field was built and why this essay had to be a separate one.

The comparison also settles a question the table above raises and does not answer: whether a growth factor of 1,330 is bad. On its own the number is meaningless. Multiplied by the unit roundoff it is a backward error of 3·10⁻¹³, and multiplied by the condition number of the problem it is a forward error. For a κ of 10³ that is 3·10⁻¹⁰ and nobody cares; for a κ of 10¹³ it is 3, and the answer has no correct digits — at a threshold every library would consider aggressive but not absurd.

So the threshold cannot be chosen without knowing the conditioning, which is the sentence this site has been arriving at from a different direction in every field. The knob is not really a memory-against-safety knob. It is a memory-against-backward-error knob, and whether the backward error matters is a property of the problem rather than of the elimination.

What is asserted here #

Every threshold in the range gives a factorisation and an answer: ‖PA − LU‖/‖A‖ under 10⁻¹² and a relative error under 10⁻¹¹, at all seven values, at every grid size on the slider.

The loosest threshold produces less fill and more growth than the strictest, both directions asserted, because a trade with only one half checked is a claim about one quantity.

Partial pivoting keeps the growth factor under 20, which is the control the whole comparison needs.

The library default is under both extremes — less fill than partial pivoting, and growth within two decades of it.

And τ = 0 is asserted to fail: ‖PA − LU‖/‖A‖ above 1, which is what defines the bottom of the usable range and is why the sweep starts at 10⁻³.

The refusal: that a factorisation reproducing its matrix to 3.8·10⁻¹⁷ therefore gives the right answer. It is fed the unpivoted run from structure and stability stop being separable, which does exactly that and is wrong in the fifth digit, and it refuses.

A closing note on what a knob is for #

The site has two other sliders that are genuine parameters of an argument rather than of a presentation — the mantissa, and the condition number — and both of them have a right end. More bits is better; smaller κ is better; the slider exists so that the reader can watch a claim hold across a range rather than at a point.

This one has no right end, and that changes what dragging it is for. Watching the two curves cross is not watching a claim hold at every position; it is watching two claims, each true over part of the range, hand over to each other. The figure would be dishonest with one curve on it, in either direction: the fill curve alone recommends τ = 0.001 and the growth curve alone recommends τ = 1, and both recommendations are wrong for the same reason.

Which is the argument for the badge this site puts on every factorisation figure, arriving from an unexpected side. The rule was written so that a decomposition could not be drawn without the evidence that it is a decomposition. What it does here is stop a memory argument from being made without the stability number beside it — and the version of this figure that plots only fill, which would be the natural figure for the sparsity field to draw, is the one the rule forbids.