Fill-In, Elimination Ordering, and Symbolic Factorization

The nonzeros of a sparse matrix do not survive elimination. Subtracting a multiple of one row from another creates entries where the matrix had zeros, and those entries are then themselves eliminated, creating more. A matrix with five nonzeros per row can produce a factor with hundreds, and the difference between a good and a bad elimination order is the difference between a factorization that fits in memory and one that does not.

Fill is a structural phenomenon: it depends only on which positions are nonzero, not on their values. That separation is what makes it tractable. The pattern of the factor can be computed exactly, before a single floating-point operation, and the ordering can be chosen to make that pattern small.

Elimination is a graph operation

Take a symmetric matrix and its graph: one vertex per index, an edge between i and j when the entry is nonzero. Eliminating vertex k removes it from the graph and connects all of its remaining neighbors to each other, because the row operations that zero out column k subtract multiples of row k from every row that has a nonzero there, and those subtractions touch every column in which row k is nonzero.

Those new edges are fill. The rule is local and exact: eliminating k adds the clique on its neighborhood. Nothing about the numerical values enters, which is why the fill pattern can be computed from the structure alone.

The order matters enormously. Consider a matrix that is tridiagonal except for a dense last row and column, the “arrow” shape. Eliminating in the natural order, from the first index, leaves the arrow vertex connected to everything, and each elimination fills in more of the matrix, ending with a completely dense factor. Eliminating the arrow vertex first leaves the remaining matrix tridiagonal, and the factor is as sparse as the original. The same matrix, two orders, a quadratic difference in storage.

What the order costs on structured problems

For a matrix from a two-dimensional grid with N unknowns, the natural row-by-row ordering produces a bandwidth of about sqrt(N) and a factor with O(N^{3/2}) nonzeros. The work is O(N^2), and the storage is O(N^{3/2}).

A nested dissection ordering does better. The grid is recursively split by a separator, the separator is eliminated last, and the two halves are independent until then. Each level of the recursion has separators totaling O(sqrt(N)) vertices, and there are log N levels, so the factor has O(N log N) nonzeros and the work is O(N^{3/2}). Both the storage and the arithmetic improve by a factor of sqrt(N) and sqrt(N)/log N respectively.

In three dimensions the exponent worsens. A separator of a n x n x n grid with N = n^3 unknowns has O(n^2) = O(N^{2/3}) vertices, the factor has O(N^{4/3}) nonzeros, and the work is O(N^2). The exponents are the reason three-dimensional sparse direct solvers are used with care: the factor of a problem that fits in memory can fail to fit once it is factorized.

Structure	Ordering	Nonzeros in factor	Flops
Tridiagonal	natural	O(N)	O(N)
2D grid	natural	O(N^{3/2})	O(N^2)
2D grid	nested dissection	O(N log N)	O(N^{3/2})
3D grid	nested dissection	O(N^{4/3})	O(N^2)

Minimum degree: greedy and local

The obvious heuristic is to eliminate next the vertex with the fewest neighbors, on the theory that the clique created is smallest when the neighborhood is smallest. Minimum degree does exactly that, updating the degrees as the elimination proceeds.

The naive implementation is expensive because eliminating a vertex changes the adjacency of all its neighbors, and the degrees have to be updated for each of them. Maintaining explicit adjacency sets and recomputing degrees costs more than the factorization itself for large problems. The standard refinement is the quotient graph representation, which stores eliminated vertices as groups rather than updating every edge individually, and the approximate minimum degree variant, which bounds the degree estimate rather than computing it exactly and is what most libraries use in practice.

Minimum degree is local and greedy, and it has no guarantee. On the arrow matrix it eliminates the cheap-looking vertices first and produces exactly the bad order, because at the first step every vertex has degree one or two and the arrow vertex’s high degree is not yet visible. On grids it produces factors noticeably denser than nested dissection, though the gap narrows on irregular problems where no good separator is known.

Its advantage is cost and robustness. The ordering is computed in time close to the number of nonzeros, it needs no geometric information, and it works on any graph, including ones that do not come from a mesh. For matrices with irregular structure, which is most matrices outside of discretized PDEs, it is often the best available choice.

Nested dissection and the separator

Nested dissection assumes the graph has small separators: sets of vertices whose removal splits the graph into pieces of at most some fraction of the size. The algorithm finds a separator, orders it last, and recurses on the pieces.

The quality of the ordering is therefore the quality of the separators. For a mesh with known geometry, a coordinate-based bisection works: sort along the longest dimension, split at the median, and take the vertices near the cut as the separator. For a graph with no coordinates, a graph partitioner such as a multilevel method must find the separator from the connectivity alone, which is a harder problem and produces weaker separators.

The recursion depth is log N, and the separator at the top level dominates the cost. That is why the top-level separator quality matters more than the deep ones, and why partitioners are usually run with a limit on how small a piece may become before the recursion switches to minimum degree. The hybrid, separators at the top and minimum degree at the bottom, is what most sparse direct solvers do.

There is a subtlety about the ordering of the separator itself. The separator’s vertices are eliminated after the two halves, but they are not independent of each other: eliminating one connects its neighbors within the separator, so the separator induces a smaller fill problem of its own, which is solved recursively. A nested dissection that ignores this produces a slightly denser factor than necessary.

Symbolic factorization predicts the pattern

Given an ordering, the nonzero pattern of the factor can be computed without any arithmetic. This is the symbolic factorization, and it has two jobs: produce the exact count of nonzeros per column so the arrays can be allocated once, and produce the elimination tree that describes the dependencies between columns.

The elimination tree has one node per column, with the parent of column j being the smallest row index greater than j that has a nonzero in column j of the factor. It is a tree even for a structurally symmetric matrix, it can be computed in time proportional to the number of nonzeros of the matrix, and it encodes which columns depend on which. A column can be computed once all of its descendants are done, which is precisely the parallelism available in the factorization.

def symbolic_eliminate(n, adj, order):
    """Return the fill count for eliminating in the given order.
    adj is a list of sets, modified in place."""
    fill = 0
    for v in order:
        nbrs = [u for u in adj[v] if u != v]
        for i, u in enumerate(nbrs):
            adj[u].discard(v)
            for w in nbrs[i + 1:]:
                if w not in adj[u]:
                    adj[u].add(w)
                    adj[w].add(u)
                    fill += 1
        adj[v].clear()
    return fill

The cost of the symbolic phase is proportional to the number of nonzeros in the factor, which is the same order as the storage, so it never dominates the numeric factorization. It does dominate when the factorization is repeated with the same pattern and different values, as in a nonlinear solver that reassembles the Jacobian but keeps its structure: the symbolic work is done once, and each subsequent numeric factorization reuses the pattern and the allocated arrays.

The elimination tree also gives the column counts, which determine the storage. A column j has nonzeros in the rows corresponding to the nodes on the path from j to the root of the subtree, and those counts can be computed by a traversal that is linear in the size of the factor. The result is an exact allocation, with no reallocation during the numeric phase.

Where the ordering meets the arithmetic

For a symmetric positive definite matrix, the ordering is purely a matter of sparsity. Any permutation preserves symmetry and definiteness, no pivoting is required, and the factorization of the permuted matrix is the permuted factorization of the original. This is why the Cholesky case is the clean one, and why Cholesky Factorization as a Positive-Definiteness Test can be applied directly to a reordered matrix.

For a general matrix, the ordering and the pivoting conflict. A fill-reducing order may place a zero or tiny entry on the diagonal, and a stability-preserving pivot may destroy the fill-reducing structure. Sparse LU solvers resolve this with threshold partial pivoting: a candidate is accepted if its magnitude is at least a fraction of the largest in its column, and otherwise the pivot is deferred to a later column. The deferral preserves stability at the cost of additional fill, and the amount depends on how badly the matrix violates diagonal dominance.

The structural prediction also assumes that no numerical cancellation occurs. In exact arithmetic, an entry that the symbolic phase predicts to be nonzero may turn out to be zero; in floating point, it will be a rounding error rather than zero. Direct solvers therefore keep the predicted structure rather than compressing it, and the factor’s storage is determined by the pattern, not by the values. Dropping the entries that turn out to be small is the idea behind incomplete factorization, which is a different trade entirely and is described in Incomplete LU Preconditioning and Dropped Fill.

The ordering decision has to be made before the factorization, and it cannot be revised cheaply afterwards. That is why the choice of ordering is the single largest lever on the cost of a sparse direct solve, larger than the kernel implementation and larger than the choice of storage format, which is described in Compressed Sparse Row Storage and the Cost of Insertion.