Sparsity, and what elimination costs

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.

Worth reading first: Structure and stability stop being separable · The factor is not sparse · The regularisation that legalises every order.

The freedom a symmetric factorisation does not have expected a conflict and did not find one. A symmetric indefinite factorisation permutes rows and columns together, so it cannot step around a bad pivot by choosing a different column, and the sparsest pivot and the sound one should fight harder than in the unsymmetric case. On a saddle-point matrix whose constraint rows have no diagonal entry at all, they did not fight: re-choosing the sparsest available pivot at every step held 72 entries in the factor against the natural order’s 113, with growth of 1.28 against 1.83, at every setting of the Bunch–Kaufman test. The explanation offered was deferral. The rows with no diagonal are the ones a sparsity rule eliminates last, and by then elimination has put something on their diagonal.

That essay’s rule re-reads the degree of every active row at every step, which no production code does. A sparse symmetric solver computes a fill-reducing permutation from the pattern alone, in a symbolic phase, before any number is read, and then pivots inside it — passing a planned pivot over for a later one only when the numerical test refuses it. The essay’s last section named the measurement: if the mechanism is deferral, a minimum-degree permutation of the pattern also puts the constraint rows late and should keep most of the benefit. It asked what the symbolic phase can promise.

It keeps all of the benefit, and slightly more. It does not defer anything.

Three orders on one factorisation

The factorisation is the previous essay’s symmetric LDLTLDL^{\mathsf T} with Bunch and Kaufman’s test at its classical constant, which at each step accepts a diagonal pivot, swaps in a better one, or takes a two-by-two block. Three rules decide which candidate it tries first. The natural order tries them as they come. The sparsest rule sorts the active rows by their current number of nonzeros, re-counted at every step. The static rule sorts them by a minimum-degree ranking computed once, before the factorisation, on the pattern of the matrix: repeatedly eliminate the vertex of least degree in the elimination graph and join its neighbours into a clique. It never looks at a value. In all three the test may refuse the first candidate and take a later one, and every factor is checked against the permuted matrix to rounding.

The saddle family is the previous essay’s: a tridiagonal Hessian on 24 variables and m constraint rows, each touching a few consecutive variables, with a zero constraint block. It is built at 2, 4, 6 and 8 constraints touching 2, 3, 4 or 6 variables each — sixteen matrices.

The static order loses nothing

Entries in a symmetric indefinite factor of a saddle-point matrix, three orders, constraints touching 3 variables eachA tridiagonal Hessian of 24 variables with m constraints, each touching 3 consecutive variables and with no diagonal entry of its own, factorised with Bunch–Kaufman pivoting at the classical constant: the entries in the factor against m. natural order: 78, 113, 150, 161; static minimum degree: 56, 65, 73, 80; sparsest, every step: 62, 72, 77, 81 at m = 2, 4, 6, 8. Growth: natural order 1.32, 1.83, 1.25, 1.25; static minimum degree 1.32, 1.28, 1.41, 1.29; sparsest, every step 1.32, 1.28, 1.25, 1.30.entries, touch 3m = 8, static80m = 8, sparsest every step81m = 8, natural16124684080120160200240constraints mentries in the factornatural orderstatic minimum degreesparsest, every stepclassical Bunch–Kaufman constant throughoutan order fixed in advance loses nothing
Fig. 1 Entries in the factor against the number of constraints, for the natural order, the static minimum-degree order and the step-by-step sparsest rule, with constraints touching three variables each. The dial sets how many variables each constraint touches.

With constraints touching three variables, the natural order’s factor holds 78, 113, 150 and 161 entries at 2, 4, 6 and 8 constraints; the step-by-step sparsest rule’s holds 62, 72, 77 and 81; the static order’s holds 56, 65, 73 and 80. Turn the dial: with constraints touching six variables the three are 88, 142, 177 and 223 against 71, 86, 121 and 120 against 68, 82, 118 and 129. Over all sixteen matrices the static order holds fewer entries than the natural one on every one, and no more than the step-by-step rule on fourteen. Its growth is under two on all sixteen — at most 1.63 — where the natural order reaches 4.72 on one.

The same holds across the test’s constant. On the four-constraint matrix at the previous essay’s five values of the constant, the static order holds 63 to 65 entries against the sparsest rule’s 70 to 74 and the natural order’s 113, with the same growth as the sparsest rule at every value: 4.16 at the loosest constants, where the test accepts nearly any pivot, and 1.28 from the classical value up.

So the question the previous essay asked — what can the symbolic phase promise — has a strong answer on this family. It can promise the whole of what re-reading the degrees buys, and it computes it once, from a pattern it already has.

It does not defer

Where each order eliminates the 8 constraint rows, among 32 stepsOn the saddle family with 8 constraints touching 3 variables each, the elimination step at which each constraint row is pivoted, for each order; the dashed line marks the start of the last 8 steps, where an order that deferred the constraint rows would put them. natural order: 3, 12, 19, 25, 26, 27, 29, 31; static minimum degree: 5, 8, 12, 16, 19, 23, 27, 31; sparsest, every step: 5, 8, 13, 16, 19, 22, 27, 30.048121620242832elimination stepnatural orderstatic minimum degreesparsest, every stepdashed: the last 8 steps beginthe constraint rows are spread through the elimination
Fig. 2 The elimination step at which each of eight constraint rows is pivoted, on the matrix with constraints touching three variables, for each order; the dashed line marks the start of the last eight steps.

Deferral makes a prediction that is easy to check: the constraint rows should be eliminated at the end. On the eight-constraint matrix the static order eliminates them at steps 5, 8, 12, 16, 19, 23, 27 and 31 of 32. One of eight is in the last eight steps. The step-by-step sparsest rule is the same shape — 5, 8, 13, 16, 19, 22, 27 and 30 — and so, in the end, is the rule whose advantage was being explained. It was never deferring. It was spreading the constraint rows through the elimination at a steady interval, one every three or four steps, which is exactly how the constraints are spread along the variables they touch.

The natural order, by contrast, puts five of its eight in the last eight steps, because the constraint rows come last in the matrix and the test takes them in turn when nothing earlier is refused. Of the three, it is the one that defers — and it is the one with twice the fill.

How many of its own variables each constraint row waits for, under the static order and the step-by-step ruleOver all 80 constraint rows of the sixteen saddle settings, the share of the variables each touches that are eliminated before it, in bins of a fifth. static minimum degree: 2, 6, 20, 35, 17, median 0.67; sparsest, every step: 1, 7, 20, 32, 20, median 0.67. A constraint row pivoted after all of its variables has a diagonal entry of its own when its turn comes.share waited forstatic minimum degree, median0.67sparsest, every step, median0.67020406080100120share of the row's own variables eliminated firstconstraint rowsunder 0.20.2–0.40.4–0.60.6–0.80.8–1static minimum degreesparsest, every step80 constraint rows, sixteen settingseach row waits for its own variables, not for the end
Fig. 3 Over all eighty constraint rows of the sixteen matrices, the share of the variables each touches that are eliminated before it, for the static order and the step-by-step rule.

What the two good orders share is local rather than global. Each constraint row is pivoted after most, but not all, of the variables it touches: the median share is two thirds under both, and 52 of the static order’s 80 rows sit between 0.6 and 1. Minimum degree does this without being told anything about constraints. A constraint row touching three variables has degree three; each of its variables has degree two in the tridiagonal Hessian plus one for each constraint touching it. The variables at the ends of a run are eliminated first, joining their neighbours and the constraint into cliques, and by the time the constraint’s own degree is the smallest, its first variables are gone and elimination has written a number where its zero diagonal was. That number is what the Bunch–Kaufman test needs, and it arrives from the variables, not from the end of the matrix.

The previous essay’s reading therefore had the right cause in the wrong place. The zero diagonal does fill in before the row is pivoted, which the previous essay measured directly. But what fills it is the row’s own variables, and waiting for them is a local ordering that interleaves the constraints with the variables rather than a global one that puts them last. An ordering that does not wait for the numbers found a family on which the ordering could be chosen for fill alone; the saddle family turns out to be another, for this reason.

The field’s other conflicts, read again

This is not the first time the field has found a conflict between fill and stability smaller than expected, and the pattern is worth naming. A threshold between fill and growth measured the unsymmetric version on a grid built so that the sparsest pivot and the sound one were different entries at every step, and there the conflict was real: a looser threshold bought fill with growth, from 172 entries and growth of 1,330 to 260 entries and growth of 1.2. The column that was never fixed found that choosing the column as well as the row changed that argument, and how few columns the search needs that a single sparsest column captured most of the full search’s benefit.

The saddle family sits at the other end. Its hard pivots — the zero diagonal of the constraint rows — are hard only until their own variables are eliminated, and an order chosen for sparsity eliminates those variables first for its own reasons. Fill and stability ask for the same order here, which is why the static order, which knows only about fill, is never refused on the eight-constraint matrix. The conflict the unsymmetric grid had was built in; the saddle family’s absence of one is built in too, by the structure of a constraint, and it is the structure the ordering reads.

That is also why the order decides the memory could compare orderings on fill alone: on a matrix whose difficult pivots are structural, the ordering that minimises fill is also the one that makes the pivots easy, and nothing numerical needs to be consulted until the factorisation starts.

Where the plan stops fitting

The previous essay’s second question was the harder test of any explanation: a symmetric indefinite matrix with small diagonal entries scattered through it, rather than a zero block at one end, has the same shortage of good pivots and no single place to put them. Whether the sparsest order still wins there, or whether the conflict finally appears, was the measurement it proposed.

The matrix is an 8 × 8 grid Laplacian with its diagonal’s sign alternated over a sublattice, which makes it indefinite without changing its pattern, and with a random share of its diagonal entries set to 10−1210^{-12} — from none to half of them, four seeds at each share.

Entries in the factor as small pivots are scattered through a 8 × 8 indefinite grid, three ordersThe grid Laplacian with its diagonal's sign alternated over a sublattice and a share of its diagonal entries, chosen at random, set to 10⁻¹²; medians over four seeds of the entries in the factor, against the share. natural order: 519, 512, 510, 498, 495, 508; static minimum degree: 358, 399, 413, 447, 480, 530; sparsest, every step: 355, 360, 356, 374, 421, 517 at shares 0, 0.05, 0.1, 0.2, 0.3, 0.5. The static order's saving over the natural one: 31%, 22%, 19%, 10%, 3%, -4%; the step-by-step rule's: 32%, 30%, 30%, 25%, 15%, -2%.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
Fig. 4 Entries in the factor of the 8 × 8 indefinite grid against the share of diagonal entries made small, for the three orders, medians over four seeds.

With no small entries the two sparsity orders agree — 358 and 355 entries against the natural order’s 519, a saving of 31 per cent — because the pattern is all that matters and they read the same pattern. As small entries are scattered in, the static order’s saving falls steadily: 22 per cent at one in twenty, 19 at one in ten, 10 at one in five, 3 at three in ten, and at half the diagonal it holds 4 per cent more than the natural order. The step-by-step rule holds on much longer: 30 per cent through one in ten, 25 at one in five, 15 at three in ten, and it too falls to the natural order’s level only when half the diagonal is small.

The reason is visible in how often the static factorisation leaves its plan. Count the steps at which the row it pivots on is not the first remaining row in the minimum-degree ranking — a planned pivot refused by the test and a later one taken. On the eight-constraint saddle matrix there are none: the plan is followed at every step, which is why the static order cost nothing there. On the grid there are none with no small entries, and then 7, 11, 16, 21 and 25 of 64 steps as the share rises to half. Every refusal takes a row out of the order the pattern was ranked for, and the ranking, computed once, does not know. The step-by-step rule re-reads the degrees after every refusal and plans the next step from the matrix that actually remains.

Growth as small pivots are scattered through a 8 × 8 indefinite grid, three ordersThe grid Laplacian with its diagonal's sign alternated over a sublattice and a share of its diagonal entries, chosen at random, set to 10⁻¹²; medians over four seeds of the growth factor, against the share. natural order: 1.14, 1.17, 1.18, 1.40, 1.40, 1.44; static minimum degree: 1.25, 1.45, 1.45, 1.52, 1.50, 1.43; sparsest, every step: 1.27, 1.27, 1.45, 1.49, 1.44, 1.62 at shares 0, 0.05, 0.1, 0.2, 0.3, 0.5.8 × 8 gridlargest growth, any order1.600.10.20.30.40.511.251.51.752share of diagonal entries made smallgrowth factornatural orderstatic minimum degreesparsest, every stepmedians over four seedsthe growth does not separate them
Fig. 5 The growth factor of the same factorisations against the share of small diagonal entries.

Growth does not separate the three. The medians sit between 1.14 and 1.62 at every share for every order, and the test’s two-by-two blocks — from none with no small entries to eleven or twelve at half — are what keep it there. So the conflict the previous essay went looking for appears, but in the currency it was not looking in: scattered small pivots do not make the sparsest order unstable, they make a sparsity order computed in advance stop being sparse.

What a symbolic phase can promise

Put in terms of what a solver can know before it reads a value:

On a matrix whose small pivots are structural — a zero block whose rows fill in from their own neighbours — a minimum-degree order computed from the pattern is followed almost without refusal — 27 refused steps over all sixteen matrices, none on the eight-constraint one — and does as well as re-reading the degrees at every step. The symbolic phase’s fill prediction is then the fill. This is the case the saddle family represents, and it is the common one for constrained optimisation.

On a matrix whose small pivots are scattered values, the pattern cannot see them, the plan is refused in proportion to how many there are, and each refusal costs fill the symbolic phase did not predict. At a share of one in ten the static factor held 16 per cent more entries than the step-by-step rule’s. A code that must promise its memory before factorising should read the refusal count as the warning: it is zero when the plan fits, it rises with the damage, and it costs nothing to keep. What survives one step of the barrier found an interior-point method’s pattern reusable for ever and its numbers not; the refusal count is how a code would know, step by step, whether the pattern’s plan is still the right one for the numbers it now has. The order that was right last time met the same phenomenon across a sequence of matrices, where an order computed for one is refused by the next; here it is refused by the same matrix’s own values.

Thirty-two rows, dense storage and exact minimum degree

Small matrices — 26 to 32 rows for the saddle family, 64 for the grid — dense storage of a sparse factorisation, and minimum degree without the approximations a production code uses. A real approximate-minimum-degree ordering differs from exact minimum degree in its tie-breaking, which on a tridiagonal Hessian with constraints is where most of the choices are; the essay’s claim is about the shape of the result, not the counts. The scattered small entries are exactly 10−1210^{-12} and the rest of the diagonal is of order one, so every small entry is refused as a pivot; a spread of magnitudes would refuse some and not others.

Why the static order holds a few entries fewer than the step-by-step rule on most of the saddle matrices — 65 against 72 on the four-constraint one, 80 against 81 on the eight — is not settled here. Both rank rows by degree, but the step-by-step rule counts the active row’s stored entries after each factorisation step, including a two-by-two block’s effect on both of its rows at once, while the minimum-degree ranking counts graph neighbours of single vertices; on a tridiagonal Hessian nearly every choice is a tie between rows of equal degree, and the two rules break those ties differently. A comparison with ties broken identically would say whether the difference is in the ties or in what is being counted.

Delayed pivoting here is unbounded: a refused candidate can be passed over for any later row. A supernodal code delays only within a supernode or to its parent, which constrains the damage differently, and whether the refusal count still predicts the fill under that constraint is not measured.

Still open: a plan that is revised, and the blocks’ real price

A plan revised after each refusal. The step-by-step rule re-reads every degree at every step; the static rule never does. A code that kept the symbolic order but recomputed it for the remaining rows only when a pivot is refused would pay for re-ordering in proportion to the refusal count — nothing on the saddle family, a quarter of the steps at half the grid’s diagonal small. The prediction with a sign is that it recovers most of the step-by-step rule’s saving at a fraction of its re-ordering work, because the refusals cluster: once one small pivot is passed over, its neighbours’ degrees are what changed.

The two-by-two pivots’ real price. The previous essay counted a block as 1.33 stored entries apiece and asked whether applying one in a solve costs the same. The scattered grid takes up to twelve blocks, far more than the saddle family’s one to three, and is the natural place to put a number on it.

Several constraint types. The saddle family’s constraints are local and alike. A matrix whose constraints couple distant variables — a global conservation law beside local ones — would give minimum degree a constraint row whose own variables are eliminated late for reasons of their own, and whether the plan is still followed without refusal there is the case that would test the local explanation hardest.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Bunch–KaufmanFill-inFill-reducing orderingGrowth factorMinimum degreeSaddle-point systemsSparse pivotingSymbolic factorisationSymmetric indefinite