Sparse Gaussian elimination

Sparse Gaussian elimination solves a linear system by Gaussian elimination while exploiting its sparse matrix structure to avoid storing zeros and performing unnecessary arithmetic operations. The catch is that zeros in the input need not remain zero: eliminating a shared unknown can couple two equations that previously had no direct coupling. These new positions are fill-in, and their number can depend strongly on the order of elimination.

For a real symmetric positive-definite matrix, the usual direct method is a Cholesky decomposition. Choose an ordering, use symbolic factorization to allocate the predicted factor entries, compute their values, then solve two triangular systems. Compressed sparse storage makes the reduced storage useful. General LU decomposition may also need numerical partial pivoting, so an ordering chosen for sparsity can change during factorization.^[ref]

For a fixed ordering, let $P$ be the permutation matrix that puts the unknowns in that order. Permute both rows and columns:Here $C$ is the lower triangular matrix in the Cholesky decomposition of $B$; $z=Px$. Forward substitution and back substitution give the original solution. The simultaneous permutation preserves the symmetric positive-definite matrix property because $v^{T}Bv=(P^{T}v{)}^{T}A(P^{T}v)>0$ for every nonzero vector $v$. A row permutation alone generally loses the symmetric matrix structure.

The local mechanism is a Schur complement. If the active pivot element is $a_{kk}$, then for remaining indices $i,j$,Even when $a_{ij}=0$, two nonzero couplings through $k$ make the second term nonzero. Thus the elimination graph removes $k$ and joins its remaining neighbors into a clique. The fill-path criterion for symmetric elimination extends this local rule across several steps. The elimination tree organizes the resulting dependencies; the minimum degree algorithm and nested dissection choose orderings to control them.

This connection between graphs and linear systems was already explicit in Donald J. Rose, R. Endre Tarjan and George S. Lueker's 1976 study of fill-in and perfect elimination orderings.^[ref]

Take a star graph $G$ with center $c$ and four leaves $1,2,3,4$. Let $L_G$ be its Laplacian matrix and $I$ the identity matrix. In the order $(c,1,2,3,4)$, the original example isFor a nonzero vector $u$, its quadratic form isThus $A$ is a symmetric positive-definite matrix and either ordering below admits a Cholesky decomposition. Changing the ordering changes the work, rather than the mathematical solvability.

Center first. Eliminating $c$ with pivot element $5$ leaves the Schur complementEvery pair of leaves now has coupling $-1/5$: six new edges, hence six undirected fill-in positions. The four remaining vertices form a clique. In the Cholesky decomposition, the first column has $C_{cc}=\sqrt{5}$ and $C_{ic}=-1/\sqrt{5}$. The rest of the lower triangular matrix is dense. To check that later cancellation does not remove its entries, after $r$ leaves have also been eliminated, the active leaf block isIndeed a leaf step sends $\alpha$ to $\alpha +{\alpha }^2/(2-\alpha )$. All active off-diagonal entries remain negative. The successive pivot elements in ordinary symmetric Gaussian elimination, equivalently the squared diagonals in the Cholesky decomposition, areThe Cholesky decomposition therefore stores $5+4+3+2+1=15$ nonzero entries, including its diagonal.

Leaves first. Permute to $(1,2,3,4,c)$. Each leaf has only one remaining neighbor, so its elimination changes only the center's diagonal, reducing it by $(-1{)}^2/2=1/2$. The successive pivot elements are $2,2,2,2,3$, and the exact Cholesky decomposition isThere is no fill-in, and only $5+4=9$ nonzero entries in $C$. This is a concrete instance of leaf elimination of a tree-pattern positive-definite matrix; its general proof explains why the same idea works on every tree pattern. The minimum degree algorithm also chooses a leaf at every nontrivial stage of this particular example, but that success gives no general optimality guarantee.

The input has 13 nonzero matrix elements when both symmetric halves are counted, or nine in its lower triangle. The six added undirected edges correspond to six added positions in one triangular factor, or 12 positions if both symmetric halves of the filled pattern are displayed. Comparing all 13 input entries directly with the 15 entries of $C$ would mix counting conventions.

Associate an undirected graph with a symmetric matrix's off-diagonal nonzero pattern. Eliminating vertex $i$ updates the remaining matrix by $a_{jk}\mapsto a_{jk}-a_{ji}{a}_{ik}/a_{ii}$. Only pairs of neighbors of $i$ can acquire fill-in. If the graph is a tree, remove leaves successively: each eliminated vertex has at most one remaining neighbor, so only a diagonal entry can change. Every Schur complement of a positive-definite matrix is positive definite, so no zero pivot forces a different ordering. Leaf removal therefore constructs a perfect elimination ordering and a fill-free factorization.

Symbolic factorization computes possible nonzero positions and update dependencies from a sparse matrix pattern and ordering before calculating numerical factor values. The elimination graph predicts fill and the elimination tree organizes dependencies. Numerical pivoting or exceptional cancellation can make the computed factor pattern differ from this structural prediction.

For a symmetric factorization in a fixed ordering, the elimination tree assigns the displayed parent to a column of the lower-triangular factor when that set is nonempty. A column with no later entry is a root. The tree records dependencies of descendant updates on later pivots, supports symbolic factor allocation, and exposes parallel subtrees. It describes the factorization ordering, not merely a spanning tree of the original matrix graph.

A vertex separator divides the remaining unknowns into disconnected subdomains. Nested dissection recursively eliminates each subdomain before its separator, keeping interaction localized until late in the factorization. For a regular two-dimensional mesh with $N$ unknowns, separators have size $O(\sqrt{N})$; the resulting direct factorization uses $O(N^{3/2})$ work and $O(N\log N)$ storage under the usual geometric separator assumptions.

At each step choose a remaining vertex of minimum degree in the current elimination graph. This heuristic reduces the size of the immediate Schur complement update and often reduces fill-in, but does not guarantee a globally optimal ordering. Its degrees must reflect earlier fill, rather than just the original graph.

Fill-in consists of new nonzero positions created by Gaussian elimination, absent from the original sparse matrix in the same ordering. For a Cholesky decomposition $PA{P}^T=C{C}^T$, compare $C$ with the lower triangle of $PA{P}^T$, not with the unpermuted input. Each new undirected edge of the elimination graph contributes one structural position in $C$. The sparse elimination of a four-leaf star gives six such positions when the center is eliminated first, and none when the leaves are eliminated first.

Symbolic factorization predicts possible nonzeros from the pattern alone. An initially zero $a_{ij}$ becomes nonzero in a single scalar step if both $a_{ik}$ and $a_{kj}$ are nonzero; however, at later steps several contributions can cancel. For an exact example, letEliminate indices $1,2$ in that order. The first Schur complement, in order $(2,3,4)$, isThe second pivot element updates the $(3,4)$ entry to $-1/3-(-1/2)(1)/(3/2)=0$. The remaining block is $\mathrm{diag}(7/3,7/3)$, so all four pivot elements are positive and $H$ is a symmetric positive-definite matrix. The elimination graph retains the structural edge $3-4$, but this exact choice of coefficients makes the numerical entry vanish. A small change to $H_{34}=H_{43}$ generally removes the cancellation.

Such an exact zero is different from discarding a small value computed with floating-point arithmetic. The symbolic factorization deliberately reserves space for possible entries without depending on special numerical identities. A perfect elimination ordering creates no structural fill-in at all; an ordering with cancellations need not have that property.

For a symmetric sparse matrix, its matrix graph has one vertex per unknown and an edge $i-j$ whenever the off-diagonal entry $a_{ij}$ is structurally nonzero. Simultaneous row and column permutation relabels this graph. Nonsymmetric patterns are instead represented by directed or bipartite graphs when row and column roles must be distinguished.

One step of symmetric Gaussian elimination removes the pivot vertex and connects its remaining neighbors into a clique. The rule follows from the Schur complement update $a_{ij}\leftarrow a_{ij}-a_{ik}{a}_{kj}/a_{kk}$. Repeating it predicts the generic factor sparsity; exceptional numerical cancellations can remove structurally predicted entries.

After eliminating a vertex set $E$, a generic edge between remaining vertices $i,j$ exists exactly when the original matrix graph contains a path from $i$ to $j$ whose internal vertices all belong to $E$. Inductively, the new edges at pivot $k$ concatenate two such paths through $k$; conversely split any allowed path at $k$, or retain it when it avoids $k$. This characterizes fill by graph reachability.