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.

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.

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.

The 12×12 grid Laplacian and its Cholesky factor, ordered by naturalTwo square sparsity plots side by side. The left shows the nonzeros of the matrix; the right shows the nonzeros of its Cholesky factor, with the entries created by elimination marked in a second colour.the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding
Fig. 2 The same phenomenon on a matrix with no single dominant vertex. On the grid the 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.

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.

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

Nonzeros in the Cholesky factor of the 12×12 grid Laplacian, by orderingA horizontal bar chart comparing the number of nonzero entries in the Cholesky factor under four elimination orderings, with reference marks for the matrix itself and for a dense factor.natural1739reverse Cuthill–McKee1354minimum degree1026nested dissection1413matrix: 408 entries · dense factor: 10440bandwidth 12 · 4.26× the matrixbandwidth 12 · 3.32× the matrixbandwidth 123 · 2.51× the matrixbandwidth 108 · 3.46× the matrixn = 144, five-point stencilevery ordering fills in; none avoids it
Fig. 3 And the same decision on a realistic matrix, where the answer is a factor of 1.7 rather than a factor of n. The arrowhead makes the phenomenon visible; the grid says what it is usually worth.
Fill growth under natural: the factor rises as n^1.49A log-log plot of nonzero count against matrix dimension. The matrix's own count is a straight line of slope one; the factor's is steeper; a dense factor is steeper still.10²10².³10²10³10⁴n (dimension)nonzerosdense factorfactor, n^1.49matrix, n^1.05grid Laplacians from 5×5 to 12×12fitted, not quoted
Fig. 4 The two behaviours as slopes. The arrowhead’s two orderings would be a line of slope two and a line of slope one; on the grid every ordering lies between, and none of them reaches the matrix’s own slope.
Gaussian elimination on a 4×4, one step at a timeFour copies of the same matrix: as given, and after each of the three elimination steps. The pivot in use is outlined and the entries reduced to zero are greyed.21-13-3-121-212-443-12as givenrows in the order 1 2 3 443-1201.251.252.502.51.5-30-0.5-0.52after step 1pivot 443-1202.51.5-3000.5400-0.21.4after step 2pivot 2.543-1202.51.5-3000.540003after step 3pivot 0.5‖PA − LU‖/‖A‖0largest multiplier0.75row order 4 3 2 1the pivot is chosen
Fig. 5 The arithmetic underneath the counting. Every entry this figure modifies in a position that held a zero is one unit of fill, and eliminating the arrow’s tip does it to every remaining position at once.
Elimination with and without pivoting, ε = 10⁻¹⁷The same two-by-two system solved twice. With a row swap the answer is exact; without one the upper triangular factor contains an entry of order one over epsilon and the second component of the answer is wrong.[ ε 1 ; 1 1 ] x = [ 1 ; 2 ], exact answer (1.000000, 1.000000)with partial pivoting1101U after elimination1.0000001.000000computed xbackward error 0forward error 0without10⁻¹⁷10-1·10¹⁷U after elimination0.0000001.000000computed xbackward error 0.25forward error 0.71no error is raisedgrowth 10¹⁷
Fig. 6 The decision that is not free, for contrast. A row interchange is chosen for stability and changes the answer if it is skipped; the interchange in this essay is chosen for structure and changes no answer at all.
Growth factor under partial pivoting: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 7 The other place a bound and a measurement diverge in the elimination field. There the bound is 5.5·10¹¹ and the median is 3.23; here the two orderings are both attainable and the gap between them grows without limit.