Sparsity, and what elimination costs

Two ends of the same arrow

One matrix, one row moved from the front of the elimination order to the back, and the factor goes from completely dense to no fill at all. Both factorisations are exact to rounding, and nothing numerical chose between them.

Worth reading first: The factor is not sparse · The order decides the memory.

Most of the arguments on this site are about a quantity moving by orders of magnitude while nothing in the output says so. This one is about a quantity moving from one extreme to the other — from the maximum possible to the minimum possible — on a matrix simple enough to hold in the head.

The arrowhead matrix has a dense first row, a dense first column, and a diagonal everywhere else. It is as sparse as a matrix with a full row can be: 2n − 1 entries in a matrix of n² positions.

Eliminate its first variable and the factor is completely dense. Eliminate it last and the factor has exactly the entries the matrix had. Same matrix, same arithmetic, same answer to the last bit.

The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix43 entriestip eliminated first253 entriestip eliminated last43 entries‖A − LLᵀ‖/‖A‖, tip first1.4·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 253 · sparse factor is 2n − 1 = 43one row swapped to the endnothing numerical chose between them
Fig. 1 The arrowhead matrix and its Cholesky factor computed both ways. The badge carries both residuals, because the whole argument depends on the two factorisations being equally exact. Drag the size and the dense factor grows quadratically while the sparse one grows linearly.

Two counts growing at different rates is a statement about their ratio, and the ratio is the number a caller actually pays. Across the slider:

The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix19 entriestip eliminated first55 entriestip eliminated last19 entries‖A − LLᵀ‖/‖A‖, tip first1.1·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 55 · sparse factor is 2n − 1 = 19one row swapped to the endnothing numerical chose between them
Fig. 2 Ten unknowns, the smallest the sweep draws: 55 entries against 19. The wrong order costs 2.9 times the storage — bad, and not yet alarming.
The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix31 entriestip eliminated first136 entriestip eliminated last31 entries‖A − LLᵀ‖/‖A‖, tip first8.2·10⁻¹⁷‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 136 · sparse factor is 2n − 1 = 31one row swapped to the endnothing numerical chose between them
Fig. 3 Sixteen: 136 against 31, a factor of 4.4. Six more unknowns have added 1.5 to the ratio.

And they add 1.5 to it every six unknowns, exactly. The five stops read 2.895, 4.387, 5.884, 7.382 and 8.881, and the differences are 1.492, 1.497, 1.498 and 1.499 — a straight line of slope one quarter in n. That is not a coincidence to be reported; it is the two closed forms divided:

n(n+1)/2  ÷  (2n−1)   =   n/4  +  3/8  +  (3/4)/(4n−2)

At n = 34 that evaluates to 8.5 + 0.375 + 0.0056 = 8.881, which is the measurement to four figures.

The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix55 entriestip eliminated first406 entriestip eliminated last55 entries‖A − LLᵀ‖/‖A‖, tip first1.7·10⁻¹⁶‖A − LLᵀ‖/‖A‖, tip last9.8·10⁻¹⁷dense factor is n(n+1)/2 = 406 · sparse factor is 2n − 1 = 55one row swapped to the endnothing numerical chose between them
Fig. 4 Twenty-eight: 406 entries against 55, a factor of 7.4. The dense count has passed four hundred while the sparse one has not reached sixty.
The arrowhead matrix, eliminated from each endThree sparsity plots. The first shows an arrowhead matrix with a dense first row and column. The second shows its Cholesky factor, completely dense. The third shows the factor obtained after moving the dense row to the end, which has no fill at all.the matrix67 entriestip eliminated first595 entriestip eliminated last67 entries‖A − LLᵀ‖/‖A‖, tip first5.4·10⁻¹⁷‖A − LLᵀ‖/‖A‖, tip last0dense factor is n(n+1)/2 = 595 · sparse factor is 2n − 1 = 67one row swapped to the endnothing numerical chose between them
Fig. 5 And thirty-four, the top of the slider: 595 against 67. Nearly nine times the storage for a factorisation that is exact to rounding either way.

So the penalty for the wrong ordering is not a constant to be traded off against something else. It is n/4, and it grows without bound in exactly the quantity a caller is trying to make large. A rule of thumb that said “eliminating in the wrong order costs perhaps a factor of three” would have been measured at n = 10 and would be wrong by a factor of three at n = 34.

Why each end does what it does

The rule from the factor is not sparse settles both cases immediately, which is what makes this the example everybody uses.

Tip first. The arrow’s tip is adjacent to every other variable. Eliminating it turns its neighbourhood — which is everything — into a clique, so after one step every pair of remaining variables is coupled. The matrix that remains is dense, and nothing after that can recover. The factor has n(n+1)/2 entries: every position in the lower triangle, the maximum a factor can have.

Tip last. Now the first n − 1 variables are eliminated first, and each of them is adjacent to exactly one other vertex — the tip. A vertex of degree one creates a clique of size one, which is no edges at all. So n − 1 eliminations create nothing, and then the tip is eliminated with nothing left to couple. The factor has 2n − 1 entries: exactly what the matrix had, and zero fill.

At n = 22 that is 253 against 43. At n = 34 it is 595 against 67. The ratio grows without limit, because one is quadratic and the other is linear.

The half that makes it an argument

The counts alone would be a fact about memory. What makes this an essay for a site whose rule is no decomposition without its residual is the other measurement, and it is on the badge.

‖A − LLᵀ‖/‖A‖ = 1.32·10⁻¹⁶ for the dense factorisation. 1.10·10⁻¹⁶ for the sparse one.

Both are at rounding. Neither is better. The expensive factorisation bought nothing — not accuracy, not stability, not a smaller residual — and the cheap one gave up nothing.

That is unusual on this site and it is worth dwelling on. Most of what is measured here is a trade: the road that squares the problem is convenience against conditioning, two Gram–Schmidts is a shorter loop against orthogonality, buying the accuracy back is speed against precision. Here there is no trade. One choice is better in the only dimension that differs, and the dimension that usually pays for such things does not move.

Nothing numerical chose between them, which is the sentence this essay exists for. The decision is structural, it is made before any value is read, and a code that does not make it deliberately has made it accidentally.

It is also why the residual badge is on this figure rather than only on the ones where a residual is in doubt. A badge that appeared only where the news was interesting would be a badge whose absence carried information, and the rule the site is named for would become a rule about presentation. Here the residual is the boring half of the measurement and printing it is what licenses the rest.

The tree, and the ordering that is genuinely perfect

The reason the tip-last ordering works is worth generalising properly, because the general statement is a clean one and it is the only place in this field where a perfect answer exists.

The arrowhead’s graph is a star: one central vertex joined to n − 1 leaves, and no other edges. A star is a tree, and every tree has a perfect elimination ordering — eliminate the leaves first, working inward. A leaf has degree one, eliminating it creates no edges, and removing it leaves a smaller tree. The induction runs all the way down.

So for a tree the fill is exactly zero under the right ordering, and the right ordering is found by a depth-first traversal in linear time. No heuristic, no NP-hardness, no measurement: a tree-structured problem is solved by elimination at the cost of storing the matrix, and that is why banded and tree-structured solvers are separate specialised routines rather than cases of the general one.

The tridiagonal matrix is the other member of the family worth naming. Its graph is a path, a path is a tree, and its natural ordering is already perfect — which is why the factor is not sparse measures the 1D Laplacian’s factor at 59 entries with no fill at all, and why incomplete Cholesky on that matrix is the complete one, a degeneracy changing the condition number on purpose has to work around.

The grid is where it stops. A grid graph contains cycles of length four everywhere, it is not a tree, and no elimination ordering avoids fill on it. Every practical two-dimensional problem sits on the wrong side of that line.

Where the arrowhead comes from

It would be easy to file this as a constructed example, and it is not. The structure occurs whenever a problem has one variable coupled to everything.

A global constraint — total mass, total charge, a Lagrange multiplier enforcing a condition over the whole domain — produces exactly this row. So does a shared parameter in a statistical model where every observation depends on one common quantity. So does a hub in a network, a bus in a power system, a single shared resource in a scheduling problem.

In every one of those the natural numbering puts the special variable first, because it is the important one and the person writing the model lists it first. That is exactly the ordering that produces the dense factor.

The natural ordering is the worst one, and it is natural for a reason that has nothing to do with elimination. A modelling convention and a numerical catastrophe, meeting because nothing connects them.

On a matrix with no single dominant vertex the same four orderings differ by a factor of 1.7 rather than by a factor of n. The arrowhead is the extreme case and the extreme case is not the representative one — it is chosen here because it makes the mechanism visible, not because a matrix in the wild costs its user a factor of nine.

The walk-back, which matters more than the example

The arrowhead is a teaching example and it teaches one thing too many. Shown a matrix where the right ordering eliminates fill entirely, the natural conclusion is that fill is an artefact of bad ordering and a good ordering removes it.

It does not, and the site’s refusal is exactly that claim. On the 12×12 grid Laplacian, minimum degree — the best of the four orderings in the order decides the memory — produces 1,026 entries against the matrix’s 408. It more than doubles the storage, and it is the best available.

An assertion that a good ordering gives zero fill on the grid is fed to rejects() and must throw. It does. Without it, this essay would be free to leave a reader with a belief its own field disproves.

The correct generalisation is narrower and less satisfying: the arrowhead has a perfect elimination ordering because its graph is a tree, and a tree can always be eliminated from the leaves inward with no fill at all. Graphs that are trees, or nearly trees, have this property. Grids do not, and almost nothing from a two-dimensional discretisation does.

The class with the property has a name — chordal graphs, those with a perfect elimination ordering — and the general fact is that fill is precisely the number of edges needed to make the graph chordal. Finding the smallest such set is the NP-hard problem the whole ordering field is approximating.

What the arrowhead is really demonstrating

Not that ordering can eliminate fill, which is false in general. Two things that are true in general, and the example makes both unmissable.

The spread is unbounded. There is no constant c such that the worst ordering is at most c times the best. Here the ratio is n(n+1)/2 divided by 2n−1, which grows linearly forever. So no argument of the form “the ordering only matters by a constant factor” survives, and the question cannot be dismissed on size grounds.

And the cost is invisible in the answer. Both routes return the same solution with the same residual. A code that has quietly chosen the dense ordering produces correct results, passes every numerical test, and fails only by exhausting memory at a size that seems too small. There is no warning, because nothing is wrong — the computation is doing exactly what it was asked to do, at a cost nobody chose.

That second property is the site’s own subject in a structural rather than a numerical form. Every essay here is about a quantity that moves without the output changing shape; usually the quantity is an error and here it is a memory footprint, and the failure mode — correct answers, silently unaffordable — is the same one.

What a solver does about it

The structure is common enough that the response to it is a standard technique rather than a judgement call, and knowing the technique makes the example less of a curiosity.

Detect and postpone. A row that is much denser than the others is identified in the analysis phase and moved to the end of the ordering. Minimum degree does this automatically without being told — the tip has the highest degree in the graph, so it is the last vertex the heuristic will ever pick, and running minimum degree on the arrowhead recovers the good ordering with no special-casing.

That is a genuine point in the heuristic’s favour and worth noticing: the greedy rule that the order decides the memory measures as merely the best of four on the grid is exactly optimal here, for free, because “eliminate the least connected variable first” is precisely the leaves-inward rule a tree wants.

Or handle the dense part separately. If the dense rows are few, the matrix can be split — the sparse block factorised normally and the dense rows folded in afterwards through a small dense correction. This is what a solver does when the constraint block is a handful of rows rather than one, and the cost is a dense factorisation of the constraint block alone rather than of everything.

Both responses depend on the structure being seen. A code that hands its matrix to a solver with no analysis phase, or that assembles the constraint into the leading rows and calls a banded routine, gets the dense factor and no indication that anything else was available.

Which of the two, and when

The pair is offered above as alternatives with no criterion, and there is one: the number of dense rows. A path core of m vertices with k rows attached to all of it, counting nnzL under three orderings and under the split:

m k tips first minimum degree tips last split ordering worth split worth
12 1 91 36 36 24 2.53× 1.50×
12 3 117 65 65 29 1.80× 2.24×
12 8 182 155 155 59 1.17× 2.63×
20 1 231 60 60 40 3.85× 1.50×
20 2 252 82 82 42 3.07× 1.95×
20 3 273 105 105 45 2.60× 2.33×
20 5 315 154 154 54 2.05× 2.85×
20 8 378 235 235 75 1.61× 3.13×

Minimum degree equals tips-last at every k. The claim above is made about one dense row and it holds at eight: the greedy rule finds the optimal ordering of a k-tip arrowhead unaided, at every size and every k measured. That is a stronger statement than the paragraph makes, and worth having about a heuristic that is merely the best of four on a grid.

The two responses cross at three dense rows. Reordering is worth 3.85× at k = 1 and 1.61× at k = 8 — its value decays. Splitting is worth 1.50× at k = 1 and 3.13× at k = 8, and grows. On the twenty-vertex core they meet between k = 3 and k = 5; on the twelve, between k = 2 and k = 3.

The mechanism is short and it is the essay’s own rule applied once more. k dense rows fill each other whatever the order, because eliminating any vertex adjacent to all of them makes a clique of all of them — so the part of the fill an ordering can remove shrinks as k grows, while the part a dense correction absorbs grows as k(k+1)/2. Hence a criterion: one or two dense rows, reorder; three or more, split.

assertTheTwoResponsesCrossAtThreeDenseRows measures the whole table, requires minimum degree to be optimal at every k, and requires the two curves to move in opposite directions and cross.

What is asserted here

Tip first gives exactly n(n+1)/2 entries, as an integer equality rather than a bound: a completely dense factor, checked as such.

Tip last gives exactly zero fill, again as an integer.

And the two differ by far more than a constant, asserted at a factor of five and rising with n.

Both factorisations are exact, both residuals below 10⁻¹³, which is the assertion that turns a memory comparison into an argument.

And the fill counts agree symbolically and numerically on both orientations, so neither number is one computation’s opinion of itself.

The refusals: the claim that fill is the same whichever end is eliminated first must throw, and so must the claim that a good ordering eliminates fill on a matrix that is not a tree. Both do.

What the figure had to be careful about

One detail of the construction, because it is the kind of thing that would quietly invalidate the comparison.

The arrowhead as described — a constant diagonal with a full first row and column — is not positive definite for any interesting n. The tip’s row has n − 1 off-diagonal entries of size 1 against a diagonal entry of 4, so diagonal dominance fails as soon as n exceeds 5, and Cholesky breaks down rather than producing a factor to count.

The fix is one line: the tip’s diagonal entry is set to n·arm + 1 rather than to the common value, so the matrix is diagonally dominant at every size on the slider and Cholesky exists. That changes the values and changes nothing about the structure, which is what the essay is about — the fill count is computed from the pattern and is identical either way.

It is worth saying because the alternative would have been to draw the figure at n = 5 and quote a ratio of 3, which is a fair number to report and a bad one to make an argument from. The claim is about a ratio that grows without limit, so the figure has to reach sizes where the growth is visible, and reaching them required the matrix to be positive definite there.

The general lesson, stated carefully

Three sentences, and the third is the one that transfers.

Elimination order does not affect the answer. It can affect the cost without limit. And it is chosen by default whenever it is not chosen deliberately — by the order variables were written down in, which is a modelling decision made for reasons that have no connection to the graph.

That last sentence is the useful one outside this subject. A default that was set for one purpose governs an unrelated quantity, nothing connects the two, and the consequence is invisible until it becomes fatal. The arrowhead is the cleanest instance because both extremes are attainable and the gap between them is unbounded — but the shape recurs wherever a structural choice and a cost are coupled through machinery neither end knows about.

The projected matrix after 3 cycles, keeping 4 vectorsA 8×8 matrix drawn as a grid. The first 4 diagonal entries are the retained Ritz values; the row and column at index 4 are their couplings to the residual vector; the remaining 4×4 block is an ordinary Lanczos tridiagonal. 22 entries are above rounding out of 64.retained Ritz values on the diagonal, their couplings in one row and column10···1.3·10⁻⁵····9.5··4.9·10⁻⁶·····9·0.000325······8.52.8·10⁻⁵···1.3·10⁻⁵4.9·10⁻⁶0.0003252.8·10⁻⁵1.830.554······0.5542.250.447······0.4471.960.427······0.4271.94the shape, countedentries above rounding22of a possible64worst entry off the arrow0a diagonal and one borderand a tridiagonal tail
Fig. 6 The same shape from a restart rather than an ordering: retained Ritz values on the diagonal, their couplings in one row and column, and a tridiagonal tail. One object coupled to everything.

The same shape, one field over

The pattern is worth one more instance, because it appears in this site’s iterative field with the same structure and a different quantity.

A rate that is known in advance notes that Gauss–Seidel’s convergence rate depends on the ordering of the unknowns — its iteration matrix is built from the lower triangle, so renumbering the variables gives a genuinely different iteration, while Jacobi’s diagonal splitting is untouched by any permutation. The exact squaring relation between the two rates holds for the natural ordering of the model problem and can fail under another.

So the same relabelling that decides whether a factor fits in memory also decides how fast a stationary iteration converges, through completely unrelated machinery, and in neither case does the answer change. Red-black ordering is adopted for parallelism and inherits both consequences without being chosen for either.

That is as far as the generalisation goes and it is far enough: a permutation is invisible to the solution and visible to nearly everything else, so an ordering picked for one reason is silently setting several other quantities. The arrowhead is the case where one of them spans everything it possibly could.

What links here

Computed from the collection, not written here: the essays that point at this one.

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

Arrowhead matrixCholeskyElimination orderFill-inFill-reducing orderingGaussian eliminationResidualSparsity