Sparse direct solvers: ordering, fill-in and the cost

A finite element stiffness matrix is sparse: each node couples only to the nodes it shares an element with, so a mesh of 10000 nodes produces a matrix with tens of nonzeros per row rather than tens of thousands. Storing and factorising only the nonzeros is what makes direct solution of large systems possible at all.

The complication is that factorisation creates nonzeros. Eliminating an unknown couples all of its remaining neighbours to one another, and those new couplings are fill-in — entries that are zero in the matrix and nonzero in its factor. Fill-in dominates both the memory and the arithmetic, and the amount of it depends on nothing except the order the unknowns are eliminated in.

Every number in this article is recomputed from its inputs by a script that runs on each build of this site.

The extreme case, exactlyPermalink to “The extreme case, exactly”#

Take an arrow matrix: nonzero on the diagonal, and nonzero along the whole of the first row and first column. With 1000 unknowns the lower triangle holds 1999 entries — the diagonal and the first column — which is as sparse as a connected matrix gets.

Eliminate unknown 1 first. It is coupled to all 999 others, so eliminating it couples every one of them to every other. The remaining block becomes completely dense, and the fill is

$$\frac{(n-1)(n-2)}{2} = \frac{999 \times 998}{2} = 498501 \ \text{entries}$$

which is 99.60% of a full dense triangle. Every advantage of sparsity has been given away in the first step.

Eliminate unknown 1 last. Each of the other 999 unknowns is coupled to nothing except unknown 1, so eliminating one creates no new couplings between anything. Fill-in is zero, and the factorisation costs $O(n)$.

Same matrix, same arithmetic, same answer to the last bit. 498501 entries or none, decided by the order.

The realistic casePermalink to “The realistic case”#

Real meshes are not arrows, but the mechanism is the same and the numbers are large enough to matter. A five-point stencil on a 100 by 100 grid — the structure a structured finite difference or a regular quadrilateral mesh produces — has $n = 10000$ unknowns, 19800 edges, and therefore 49600 nonzeros in the full matrix, 29800 in the lower triangle.

Order it by rows, which gives a half-bandwidth of 100, and the banded factor occupies

$$n b - \frac{b(b+1)}{2} = 10000 \times 100 - 5050 = 994950 \ \text{entries}$$

a fill ratio of 33.4 over the matrix it came from. The factorisation cost scales as $nb^2$, which is $10^8$ operations — 100 times $n^{1.5}$, the asymptotic cost that a good two-dimensional ordering achieves.

Bandwidth ordering is not a bad choice; it is a naive one, and it is the default a program gets if it does nothing. Reverse Cuthill-McKee will do better on an unstructured mesh, and nested dissection — recursively splitting the mesh with separators and numbering each separator after the pieces it separates — achieves $O(n \log n)$ fill and $O(n^{1.5})$ flops in two dimensions, which is provably the best available for this class of problem.

Why ordering is not the engineer’s problem, and isPermalink to “Why ordering is not the engineer’s problem, and is”#

Finding the ordering that minimises fill is NP-complete, so every solver uses a heuristic: minimum degree, nested dissection, or a hybrid. The heuristics are good, they are implemented in libraries, and a user should never be choosing between them by hand.

What a user should know is what the choice implies, because it explains behaviour that otherwise looks arbitrary.

Memory grows faster than the mesh. Doubling the elements in a two-dimensional model does not double the factor; it more than doubles it. A model that fits comfortably at one mesh density can exhaust memory at the next refinement, and the jump is in the factor rather than in the model.

Element numbering affects run time and not results. Two identical models with differently numbered meshes can take noticeably different times, because the ordering heuristic starts from what it is given. They should return the same answer — and if they do not, the difference is in the summation order rather than in the ordering.

The pattern is worth computing once. In a staged construction analysis or a consolidation run, the sparsity pattern does not change between phases or between time steps, only the values do. The symbolic factorisation — deciding where fill will appear and allocating for it — can be done once and reused, and in a multi-phase run that is most of the analysis. This is why KATAI 2D establishes the pattern once and repeats only the numerical factorisation against it.

Direct or iterativePermalink to “Direct or iterative”#

The alternative is not to factorise at all: iterative solvers apply the matrix repeatedly and never form a factor, so fill-in does not arise. For very large three-dimensional problems that is the only option, because the fill in three dimensions is $O(n^{4/3})$ and grows past what can be stored.

For two-dimensional geotechnical problems the case for a direct solver is different and mostly about interpretation. A direct solve returns the solution of the linear system, to round-off. An iterative solve returns a solution converged to a tolerance, and that tolerance becomes part of every result — so when a nonlinear analysis fails to converge, the question “is this the model or the linear solver” has to be asked and answered before anything else can be. Removing that question from the interpretation of a run is worth a great deal, and at these sizes it costs little.

That is a verification argument as much as a performance one. A result whose accuracy depends on a tolerance that was not reported is not reproducible in the sense that matters.