Numerical Linear Algebra
where the algebra stops being the arithmetic
Series

Sparse pivoting — the series

13 essays on one idea, from the one that introduces it to the one that assumes the rest.
  1. 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.

    part 1 · sparsity
  2. 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.9·10⁻¹⁶growth at τ = 0.138entries in L+U at τ = 0.1372one knob, two measurements, opposite directionsand the default is most of both

    A threshold between fill and growth

    One number decides how small a pivot an elimination will accept. At 0.001 the factor holds 172 entries and the matrix grows by 1,330; at 1 it holds 260 and grows by 1.2. The libraries ship 0.1, and the measurement says why.

    part 2 · sparsity
  3. 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

    What the symbolic phase can only bound

    Without pivoting, the fill can be computed from the graph and the count is exact — 233 predicted, 233 measured. With pivoting it is 233 predicted and 242 measured, and what survives is a bound that is right at every threshold and loose by 1.7 times at the largest grid drawn.

    part 3 · sparsity
  4. 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

    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.

    part 4 · sparsity
  5. natural113 predicted · 113 countedminimum-degree63 predicted · 63 countedreverse Cuthill–McKee63 predicted · 63 countedthe shaded entries are fill: zeros of K that the factorisation makes nonzeroallocated before the numbersnatural113minimum-degree63reverse-cuthill-mckee63predicted minus counted0the symbolic phase decides the memoryand nothing later is allowed to argue

    An ordering that does not wait for the numbers

    A sparse factorisation's memory is decided by an ordering computed from the graph, and its stability by pivots computed from the values, and the two decisions fight. On one family of matrices they do not — the ordering can be chosen for fill alone, and the fill the symbolic phase predicts is the fill the factorisation produces — exactly, not as a bound.

    part 5 · sparsity
  6. 10⁻³10⁻²10⁻¹1110¹10²10³10⁴pivot threshold τgrowth factor · entries in L + U ÷ entries in Agrowth, rows onlygrowth, rows and columnsfill ratio, rows onlyfill ratio, rows and columns8×8 grid, 64 unknowns, 288 entriesτ = 0.001: rows only — 741 entries; growth2209τ = 0.001: rows and columns — 640 entries; growth2.5τ = 0.1: rows only — 875 entries; growth38τ = 0.1: rows and columns — 640 entries; growth2.5τ = 1: rows only — 986 entries; growth1.2τ = 1: rows and columns — 659 entries; growth1.2dashed: a fixed column · solid: the column chosen toothe same threshold rule in both

    The column that was never fixed

    Every threshold-pivoting measurement so far chose the pivot row in a fixed column, and the routine's own description said that choosing the column as well would change the constants and not the argument. Measured, it changes the argument. On the 8×8 conflict grid the factor shrinks from 875 entries to 640 at the library default, and the growth factor that climbed to 2,209 as the threshold loosened stays at 2.54 at every threshold from 0.3 down to 0.001. What does most of the work is not the column but which of several equally cheap entries is taken — and on random sparse matrices, choosing the column without that makes the growth worse.

    part 6 · sparsity
  7. 80 matrices, τ = 0.1, ties to the largestno column search, median fill244one column136every column109share of the gain at one column0.804896144192240columns the search may look atmedian fill123468121624allno column search at allmedian fillmedian growthworst growth — see the captionthe fill is the quantity the width buysthe worst case is the draw

    How few columns the search needs

    A full row-and-column pivot search is quadratic in the active submatrix at every step, and no library performs one. Looking at a single sparsest column takes the median fill from 244 to 136 where the full search reaches 109 — four fifths of the benefit for a linear scan — and that share is 79, 83, 80, 89 and 92 per cent across five thresholds. The worst growth appears to favour the narrow search by a factor of six, and on the next draw it favours the wide one by two.

    part 7 · sparsity
  8. a 28 × 28 saddle-point matrixnatural order, entries113sparsest order, entries72two-by-two pivots, loosest test0and at the strictest3023466992115the constant in the pivot testentries in the factor0.10.30.640.80.95133natural ordersparsest ordergrowth, naturalgrowth, sparsestnumbers above the points: two-by-two pivots takenboth curves favour the sparser order

    The freedom a symmetric factorisation does not have

    Permuting rows and columns together leaves no column to choose, so the conflict between the sparsest pivot and the sound one should be worse rather than better. On a saddle-point matrix whose constraint rows have no diagonal entry at all, it is not there: taking the sparsest available pivot holds 70 entries against the natural order's 113 and a growth of 1.28 against 1.83 — better on both currencies at once, at every setting of the pivot test. The two-by-two blocks that make it legal cost 1.33 entries apiece.

    part 8 · sparsity
  9. 8 × 8 gridstatic saving, none small0.31static saving, three in ten0.0300.10.20.30.40.5350400450500550share of diagonal entries made smallentries in the factornatural orderstatic minimum degreesparsest, every stepmedians over four seedsthe order fixed in advance stops fitting

    An order fixed before the numbers

    On a saddle-point matrix whose constraint rows have no diagonal, taking the sparsest pivot at every step beat the natural order on fill and growth at once, and the reading was deferral: the constraint rows go last and have filled in by then. A code computes its ordering once, from the pattern. That static minimum-degree order holds fewer entries than re-reading the degrees on fourteen of sixteen settings, growth under two throughout — and it does not defer. It spreads the constraint rows through the elimination, each after two thirds of its own variables. Scatter small pivots through a grid instead and the plan is refused on up to 25 of 64 steps; its saving falls from 31 per cent to 3.

    part 9 · sparsity
  10. share of the saving recovered6 × 6, revised, worst0.978 × 8, revised, worst0.8510 × 10, revised, worst0.6600.250.50.7511.25small diagonal entriesshare of the saving recovered5%10%20%30%6 × 6, static6 × 6, revised8 × 8, static8 × 8, revised10 × 10, static10 × 10, revisedsolid: 8 × 8 · dotted: 6 × 6 · dashed: 10 × 10re-planning at refusals recovers most of it

    A plan redrawn where it was refused

    A symmetric indefinite factorisation planned once from the pattern loses its fill saving as small pivots are scattered through the diagonal — on an 8 × 8 grid it keeps 74 per cent of what re-reading the degrees at every step saves at one small pivot in twenty, and 20 per cent at three in ten. Redraw the plan for the remaining rows only when the pivoting criterion rejects the planned pivot, and it keeps 85 to 100 per cent at every share up to three in ten, reading a tenth to a quarter of the rows the step-by-step rule reads. Each redrawn plan heads off refusals that would have followed it: 3 re-plans where the fixed plan was refused 7 times, 9 where it was refused 17. Where the plan was never refused, redrawing costs nothing — and with half the diagonal small, no plan of any kind beats the natural order.

    part 10 · sparsity
  11. saving with half the diagonal small, %scattered, static, at half-4.3scattered, revised, at half-4.1clustered, static, at half13clustered, revised, at half14-0.100.10.20.3small diagonal entriessaving over the natural order5%10%20%30%50%fixed, scatteredredrawn, scatteredfixed, clusteredredrawn, clusteredsolid: clustered · dashed: scattereda region keeps the plan, not breaks it

    A small pivot with a small neighbour

    A fill-reducing plan computed from the pattern loses its saving as small pivots are scattered through a symmetric indefinite matrix, and redrawing it where the pivoting test refuses it wins the saving back. Gather the same small pivots into one region and the prediction was that the fixed plan would do worse, because a whole region's ranking is wrong at once. It does better. On 8 × 8 and 10 × 10 grids the fixed plan is refused 0 to 6 times with the small pivots clustered, against 7 to 19 scattered, and at three in ten it saves 10 and 23 per cent of the natural order's factor where scattered it saved 3 and 9. The reason is in the Bunch–Kaufman test: a small pivot whose neighbour is also small cannot be passed over for that neighbour, so the test pairs the two into a block that keeps the planned row. Both plans still save something with half the diagonal small. And the step-by-step sparsest rule, best of all when the small pivots are scattered, falls behind the redrawn plan in a cluster.

    part 11 · sparsity
  12. 8 × 8scattered, refused at 30%17along lines, refused at 30%8in a disc, refused at 30%3050100150share of the diagonal that is smallentries added to the fixed plan's factor0%5%10%20%30%50%scatteredalong linesin a discthe badge counts refusals, which the curves do not followthree layouts, one cost

    A refusal has no fixed price

    A symmetric indefinite factorisation planned once from the pattern is refused by its pivoting test where small diagonal entries sit, and the earlier essays read the refusal count as the plan's cost: scattered small pivots were refused 7 to 19 times and saved little, clustered ones 0 to 6 and saved much more. Put the small pivots along lines, each with two small neighbours, and the prediction of one refusal for every three small entries holds — between 0.25 and 0.42 on every share and grid. The rest of the reading does not. From a fifth of the diagonal up, the three layouts add within a factor of 1.6 of the same number of entries to the fixed plan's factor while their refusal counts differ by three to seven times. What moved the earlier savings was the baseline: the natural order's factor changes with the layout too, by 13 per cent at three in ten, and on the 10 × 10 grid the lines 'save' 19 per cent where the scatter saves 9 — with fixed-plan factors of 907 and 898 entries.

    part 12 · sparsity
  13. -4004080120160200share of the diagonal that is smallextra entries from interchanges5%10%20%30%50%scatteredalong rowsalong columnsin a discsolid: 10 × 10 · dashed: 8 × 8the row lines pay for the natural order's interchanges

    The baseline was the interchanges

    Small pivots along grid rows gave a fixed fill-reducing plan twice the saving over the natural order that scattered ones did, for a factor of the same size, and the reading was that the natural order did badly on rows. Turning the lines to run along columns was predicted to bring the natural order below 1,000 entries, leave the plan within 3 per cent of 907, and drop the saving to the scatter's. It misses on all three: 1,008, 977 and 3.1 per cent. What the natural order paid for on rows was its own interchanges. Offered by position, as a dense code offers it, every Bunch–Kaufman interchange moves a row out of sequence, and on rows that costs 125 entries at three in ten and 191 at half — the most of any layout at every share. Kept to the grid's index, the natural order barely notices the layout, and the row lines' saving halves to 9.0 per cent.

    part 13 · sparsity

All series