The order decides the memory

The factor is not sparse establishes that elimination creates entries, that they can be counted from the graph before any arithmetic runs, and that the count grows faster than the matrix does. It leaves open the question that decides whether a problem is solvable: how much of that fill is a consequence of the matrix and how much is a consequence of the order the variables were eliminated in.

The answer is: a great deal of it is the order. And the order is free.

Why the order is free #

This is the premise everything rests on, and it takes four lines to check rather than to assume.

Renumbering the variables replaces A by PAPᵀ for a permutation P. That is a symmetric permutation: the same reordering applied to the equations and to the unknowns. It changes the matrix, it changes the factor, it changes the fill — and it does not change the problem. Solve PAPᵀy = Pb, unpermute y, and the answer is the answer.

The site checks it rather than stating it. On the 7×7 grid, factorised under a minimum-degree ordering and unpermuted, the relative error against the known solution is 3.1·10⁻¹⁶. The permuted matrix has exactly as many nonzeros as the original, which is the other half: a permutation moves entries around and does not create or destroy any.

So ordering is the rare decision that is pure gain. It costs nothing in accuracy, nothing in stability, nothing in the answer — for a symmetric positive definite matrix, where no pivoting is required for stability and the elimination order is therefore unconstrained.

The four #

Natural is the ordering the matrix arrived in — for the grid, row by row. It is the baseline and it is what a code gets by not thinking about the question.

Reverse Cuthill–McKee is a breadth-first sweep from a low-degree vertex, taking neighbours in order of increasing degree, then reversed. It optimises bandwidth rather than fill: it tries to push every nonzero close to the diagonal, on the grounds that a factor of a banded matrix is contained within the band. The reversal looks arbitrary and is not — Cuthill–McKee’s ordering and its reverse have the same bandwidth, and the reverse has a profile that is never worse and usually much better.

Minimum degree eliminates, at each step, the variable coupled to the fewest others, then adds the clique that elimination creates and repeats. It is a greedy heuristic with no optimality guarantee whatsoever, and it is what essentially every sparse direct solver has used for forty years.

Nested dissection finds a separator that splits the graph in two, numbers it last, and recurses. It is the only one of the four with an asymptotic guarantee: on a k×k grid it gives O(n log n) fill, against natural ordering’s $O(n^{1.5})$.

The numbers #

On the 12×12 grid — 144 variables, 408 entries in the matrix’s lower triangle, 10,440 in a dense factor:

The spread between best and worst is a factor of 1.7, and the spread between the best and a dense factorisation is a factor of ten. Both are worth having and the second is the one that decides feasibility.

The result that was not expected #

The table is ordered by measured fill and nested dissection is third.

It is the ordering with the proof. Its complexity bound is better than minimum degree’s — minimum degree has no bound at all — and on a grid, which is the structure nested dissection was designed for, it loses to the heuristic by 38%.

That is not a bug in the implementation, and it is not a bad grid size. It holds at every size on the slider, from 6×6 to 16×16. The crossover, where the asymptotics start to pay, is further out than this figure reaches.

The reason is that the asymptotic statement is about the leading term and the sizes anybody draws are governed by the constant. Nested dissection numbers separators last, which is optimal in the limit and wasteful at small sizes where the separator is a large fraction of the whole graph. Minimum degree makes locally good choices with no plan at all, and on a moderate grid a sequence of locally good choices is very hard to beat.

A complexity class is a statement about a limit and a matrix has a size, which is the same lesson as the bound that is never attained reaching from the elimination field: the growth factor’s bound is 5.5·10¹¹ and the measured median is 3.23. In both cases the theory is correct and the number it gives is not the number that governs the computation.

The honest position is that production solvers use both, and choose between them by trying them — which is an unsatisfying answer that happens to be the true one. Modern codes tend to use nested dissection at the top levels, where the graph is large enough for the asymptotics to matter, and switch to minimum degree on the subgraphs once they are small.

What minimum degree is actually doing, and what it costs to do it #

The heuristic that wins is worth a closer look, partly because it is the one in the libraries and partly because the plain version described above is not the one they run.

At each step it computes the degree of every remaining vertex, picks the smallest, eliminates it, and adds the resulting clique. The rationale is direct: eliminating a vertex of degree d creates at most d(d−1)/2 edges, so taking the smallest d minimises the immediate damage. It is greedy in the strict sense — it optimises the current step with no consideration of the next.

The trouble is the bookkeeping. Adding a clique changes the degrees of everything in it, so degrees must be recomputed, and the cliques accumulate until the graph being maintained is denser than the matrix. A naive implementation — which is what this site’s is, because at these sizes clarity is worth more than speed — costs more than the factorisation it is choosing the order for.

Production implementations spend essentially all their complexity on avoiding that. Quotient graphs represent a clique by a single node rather than by its edges, so the elimination graph never becomes dense. Approximate minimum degree replaces the exact degree with a bound that is cheap to update, and is the standard because the approximation costs a little fill and saves a great deal of time. Supernodes detect vertices with identical adjacency and eliminate them together, which is both faster and gives the numeric phase dense blocks to work on.

The pattern is worth noting: forty years of work on this ordering has gone almost entirely into computing it faster, and hardly any into computing a better one. That is a judgement about where the remaining value is, and the table above supports it — the spread between the four orderings is a factor of 1.7, and the spread between a good implementation and a naive one is far larger.

Bandwidth is a different objective, and it shows #

Reverse Cuthill–McKee is optimising something else, and putting it in a fill table slightly misrepresents it.

A band solver stores everything within the bandwidth, whether or not it is nonzero, and its cost is therefore set by the bandwidth alone. For such a solver RCM is exactly the right ordering, and the figure asserts what it promises: it never increases the bandwidth, and on the grid it reduces it substantially.

For a general sparse solver, which stores only what is nonzero, bandwidth is the wrong target — a matrix can have a large bandwidth and very little fill. RCM lands third-best on fill here because reducing bandwidth reduces fill as a side effect rather than as an aim.

Which is worth naming as a general trap: an ordering is only good relative to a cost model, and comparing orderings without saying which storage scheme is in use compares them at cross purposes. The figure prints the bandwidth beside each bar for that reason.

Nested dissection, and why its guarantee is the interesting one #

The ordering that comes third is the only one whose behaviour can be predicted rather than measured, and the argument for it is short enough to give.

Take the k×k grid and cut it with a vertical line of k vertices. Number everything to the left first, everything to the right second, and the separator last. Now consider what elimination does: a variable on the left and a variable on the right are not adjacent, and eliminating left-hand variables can never couple them, because every path between the two halves passes through the separator — which has not been eliminated yet.

So the two halves fill in independently, and the only coupling created between them is within the separator’s own block, which is k×k and dense. Recurse on each half and the same argument applies again. The total is a sum over levels of separator blocks, and it comes to O(n log n) entries against natural ordering’s $O(n^{1.5})$.

The construction generalises: what it needs is a small separator, a set of vertices whose removal splits the graph into pieces of bounded size. Planar graphs have separators of size O(√n), which is the theorem that makes the whole approach work for two-dimensional problems, and three-dimensional grids have separators of size $O(n^{2/3})$, which is why the exponents are worse there.

What the argument does not supply is a way to find a good separator on a graph that is not a grid. On the grid it is a line, and this site’s implementation uses that directly rather than searching. General-purpose implementations use multilevel graph partitioning to find one, and the quality of that partitioner is most of what decides whether nested dissection beats minimum degree on a particular matrix.

What no ordering does #

The refusal in the library is the belief this essay could otherwise create, and it is worth stating plainly.

No ordering removes fill. The best of the four still produces a factor 2.5 times the size of the matrix. Every ordering, at every grid size drawn, adds entries — that is asserted for all four on every frame, and an ordering that produced zero fill on this matrix would fail the build.

The arrowhead in two ends of the same arrow is the exception that makes the belief tempting: there, one ordering genuinely gives no fill at all. It is a very special structure, and reading it as the general case is the misunderstanding the refusal exists to catch.

And the underlying reason nothing better is available: finding the ordering that minimises fill is NP-hard. Every ordering in this table is a heuristic, including the one with the asymptotic guarantee — the guarantee is about a bound on the fill, not about achieving the minimum. There is no algorithm anybody expects to find that would fill the bottom row of that table with the true optimum.

What is asserted here #

Every ordering still fills in, all four, at every grid size.

Every ordering beats a dense factorisation, which is the other end of the same claim and would catch an ordering that had gone catastrophically wrong.

Some ordering beats the natural one — the claim that the decision is worth making at all.

RCM does not increase the bandwidth, which is the thing RCM is actually for.

And a symmetric permutation returns the same solution, to 3.1·10⁻¹⁶, with the permuted matrix having the same number of entries as the original.

The refusals: the claim that two orderings give the same fill must throw, and so must the claim that a good ordering eliminates fill. Both do.

The other cost the count does not show #

Fill is measured in entries, and entries are what decide whether the factor fits. They are not what decides how long the factorisation takes.

The arithmetic cost of eliminating a vertex of degree d is about d²/2 operations, so the total work is the sum of squared degrees over the elimination — a quantity that grows faster than the fill it produces. On the model problem the fill grows as $n^{1.5}$ and the work as $n^2$, and the ratio between them widens with size.

That means the ordering table understates its own case. An ordering that reduces fill by a factor of 1.7 reduces the operation count by rather more than 1.7, because the entries it avoided creating were also entries later steps would have had to process. Both numbers are reported by a real solver’s symbolic phase, and the operation count is the one to look at when time rather than memory is binding.

There is a third quantity, and on current hardware it is often the one that matters most: how much of the work happens in dense blocks. A factorisation organised into supernodes runs its arithmetic as dense matrix-matrix products, which reach a large fraction of a machine’s peak throughput, while the same operation count scattered over individual entries does not. Two orderings with identical fill and identical operation counts can differ severalfold in time for that reason alone — which is a consideration no count in this essay captures, and a reason the choice between orderings is ultimately settled by measurement.

Where this leaves a practitioner #

Three things, and the first is nearly always enough.

Use the library’s default and do not think about it. Every sparse solver applies a fill-reducing ordering automatically, usually approximate minimum degree, and the difference between it and the best available choice is tens of per cent. The difference between it and no ordering is the factor of 1.7 in the table, and rises with size.

Look at the fill count when memory is the binding constraint. The symbolic phase reports it before the numeric phase allocates, so the question “will this fit” is answerable in advance and cheaply. That is a more definite planning basis than almost anything else in this subject.

And reorder once, not every time. For a sequence of solves with the same structure and changing values — which is what a time-stepping code does — the ordering and the symbolic factorisation depend only on the pattern. They are computed once and reused, and only the numeric factorisation repeats.