Multifrontal Factorization and the Dense Front Matrix

A sparse Cholesky factorization touches a small fraction of the entries a dense one would, and yet the dense version can run several times faster per flop. The sparse code is dominated by indirect addressing, short inner loops, and irregular memory access; the dense code spends its time inside a matrix-matrix kernel that streams contiguous data at a rate close to the machine’s peak.

The multifrontal method is the resolution of that conflict. It reorganizes the sparse factorization so that nearly all of the arithmetic happens inside small dense matrices, and the sparse structure survives only as the schedule that decides which dense matrix is factorized when.

The elimination tree is the schedule

Given a fill-reducing ordering, the dependencies between columns form a tree. Column j depends on every column that contributes to it, and the parent of j is the first column after j that receives a contribution from it. The tree can be computed from the pattern alone, in time proportional to the number of nonzeros, before any arithmetic.

/* Elimination tree of a symmetric matrix given its upper triangular
   pattern in compressed sparse column form. parent[] receives the tree. */
void elimination_tree(int n, const int *col_ptr, const int *row_ind,
                      int *parent, int *ancestor)
{
    for (int k = 0; k < n; k++) { parent[k] = -1; ancestor[k] = -1; }
    for (int k = 0; k < n; k++) {
        for (int p = col_ptr[k]; p < col_ptr[k + 1]; p++) {
            int i = row_ind[p];
            while (i != -1 && i < k) {
                int next = ancestor[i];    /* walk up the path once */
                ancestor[i] = k;           /* path compression: point to k */
                if (next == -1)
                    parent[i] = k;         /* i's first parent above it */
                i = next;
            }
        }
    }
}

The tree answers the scheduling question directly: a node can be processed once all of its children are done. Leaves correspond to columns whose elimination affects nothing else, and the root corresponds to the last column eliminated, which typically has the largest front.

A front is a dense matrix with a sparse boundary

At each node of the tree, the multifrontal method assembles a small dense matrix called the front. It consists of the rows and columns that the node’s elimination touches: the node’s own column, the columns of the original matrix that have nonzeros in those rows, and any rows contributed by the children.

The assembly adds the children’s contribution blocks into the appropriate positions of the front. Each child, after being factorized, produces a Schur complement of its eliminated variables, which is a dense square block indexed by the rows that the child’s elimination coupled together. That block is exactly the update that the child would have applied to its parent in a left-looking factorization, and the multifrontal method applies it by adding it into the parent’s front.

The front is then partially factorized. Its first b1 columns, corresponding to the node’s own variables, are eliminated with a dense factorization, and the remaining part of the front is updated by the corresponding Schur complement. That Schur complement is the contribution block passed to the parent, and the eliminated columns are written out as part of the global factor.

The arithmetic for a front of size b with b1 fully summed columns is about b1^2 * b for the partial factorization and the update of the remaining block combined, since each of the b1 eliminated columns updates a panel of length b. When the whole front is eliminated at once, b1 = b and the count reduces to b^3/3, which is exactly the flop count of a dense Cholesky factorization of a b x b matrix. That identity is the point: the sparse structure has been converted into a sequence of dense problems whose sizes are determined by the ordering.

Why the front stays small

The front size is not a property of the matrix as a whole; it is a property of the separator structure at that node. For a two-dimensional grid ordered by nested dissection, the front at a node corresponding to a separator of length s has size about s, and the total work is the sum of s^3 over the recursion levels, which is O(N^{3/2}) for N unknowns.

The front is dense even though the original matrix is sparse, because the elimination of the separator couples all of its variables to each other. That density is not waste: the entries in the front are exactly the ones that the factorization needs, and the dense kernel processes them at a rate that the sparse kernel cannot approach.

The arithmetic intensity of a dense b x b factorization is proportional to b. A front of size 64 running in cache sustains a much higher fraction of peak than a front of size 8, which is why orderings that produce large separators are not penalized as much as their flop count suggests. The crossover between “more flops in a dense kernel” and “fewer flops in a sparse one” is the central tuning question of a sparse direct solver, and it is resolved by amalgamating small nodes into larger supernodes.

Supernodal and multifrontal are two schedules for the same arithmetic

A supernode is a set of consecutive columns with identical structure below the diagonal. Factoring them together turns a sequence of rank-one updates into a sequence of matrix-matrix products, which is the same conversion the multifrontal method achieves by a different route.

The supernodal approach is left-looking: for each supernode, it gathers the updates from all previously computed supernodes that touch it, applies them, and factorizes. The data it touches is contiguous, because the supernodes are contiguous ranges of columns, and the update is a sparse-sparse matrix product that is executed as a sequence of dense calls.

The multifrontal approach is right-looking: each node pushes its contribution to its parent as soon as it is done, and the parent accumulates contributions until all children have reported. The data is a stack of fronts, and the index lists that map front positions to global row indices are explicit.

The two produce the same factor and the same flop count. They differ in memory access and in how easily they accommodate pivoting. The supernodal method has better locality for the update phase because it works with contiguous panels; the multifrontal method has a simpler dependency structure and a natural recursive formulation that maps well onto task-based runtimes.

Pivoting inside a static structure

For a symmetric positive definite matrix, no pivoting is needed, and the multifrontal structure is used exactly as described: every front is factorized with the dense Cholesky kernel, and the contribution block is symmetric.

For a general matrix, the front must be factorized with partial pivoting, and the pivot may come from a row that belongs to an ancestor’s front rather than to the current one. Moving that row into the current front changes the front’s row index set, which changes the structure that was fixed by the symbolic phase. The standard resolution is to allow the front to grow: rows that are needed for stability are appended to the front, the factorization proceeds, and the rows that could not be eliminated at this node are delayed to the parent.

Delayed pivots increase the front sizes and therefore the work and the memory. On a matrix that is strongly diagonally dominant, few pivots are delayed and the structure is close to the symbolic prediction. On an indefinite or badly scaled matrix, the delays can cascade and the factorization can become substantially more expensive than the pattern suggested. This is the fundamental tension of sparse LU: the structure wants to be fixed in advance, and stability wants to be decided during the computation.

An alternative is static pivoting, which fixes the pivot sequence in advance and replaces any pivot that turns out to be too small by a perturbation of the diagonal, of the order of the square root of the machine precision times the norm of the matrix. The factorization then solves a slightly perturbed system, and a few steps of iterative refinement recover the accuracy for the original one. The approach keeps the structure static and the kernels dense, at the cost of a solve phase that must include the refinement, and it fails when the perturbation is large enough to change the conditioning materially.

Memory, and where the parallelism lives

The peak memory of a multifrontal factorization is not the sum of all fronts. Fronts are created and destroyed in a postorder traversal, and at any moment only the fronts on the path from the current node to the root, plus their pending contribution blocks, are alive. The peak is the maximum over the traversal of that live set, which is usually dominated by the root front and the largest separators near the top of the tree.

This is a different memory profile from the supernodal method, which allocates the factor’s storage up front and has a peak that is close to the final size. The multifrontal peak can be several times the size of the factor in the worst case, which matters when the factor barely fits in memory.

The memory profile is the reason the choice between a direct factorization and an iterative method is often made on storage rather than on arithmetic. A factorization that does not fit cannot be run, whatever its flop count, and the alternative that keeps memory bounded by the matrix itself is the incomplete factorization described in Incomplete LU Preconditioning and Dropped Fill.

The parallelism comes from two places. Tree parallelism is available between subtrees that have no dependency: two children of a node can be factorized concurrently, and on a wide tree that exposes a large amount of independent work. Node parallelism comes from the dense kernels themselves: the b^3/3 flops of a front are a matrix-matrix product that a threaded BLAS can split across cores.

The two forms interact. Near the leaves the tree is wide but the fronts are small, so tree parallelism dominates and the kernels are too small to thread efficiently. Near the root the tree is narrow and the fronts are large, so node parallelism dominates and there may be only one front in flight. The critical path runs through the root, and its cost is the cost of the largest front factorizations in sequence. A solver that scales well must keep both forms busy across the whole traversal, which is why task-based implementations with dynamic scheduling are used rather than a static partition.

The fill that determines all of these sizes is set by the ordering, and the ordering is computed before the factorization begins. That dependency is why the analysis in Fill-In, Elimination Ordering, and Symbolic Factorization is the first step of every sparse direct solve, and why the quality of the ordering has a larger effect on the total cost than any decision made inside the kernels.