= 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.<Applications of Parallel Computers, CS267, Spring 1996/Lecture 27: Sparse Gaussian Elimination>[Applications of Parallel Computers, CS267, Spring 1996; CS267, Lecture 27, 23 April 1996, Cholesky versus LU Decomposition: distinction between numerical pivot selection and advance ordering, and phases 1–4; Sparse Matrix Storage Schemes]{ref}
For a fixed ordering, let $P$ be the <permutation matrix> that puts the unknowns in that order. Permute both rows and columns:
$$
B=PAP^T=CC^T,\qquad Cy=Pb,\qquad C^Tz=y,\qquad x=P^Tz.
$$
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^TBv=(P^Tv)^TA(P^Tv)>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$,
$$
a_{ij}^{\mathrm{new}}=a_{ij}-\frac{a_{ik}a_{kj}}{a_{kk}}.
$$
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>.<Algorithmic Aspects of Vertex Elimination on Graphs by Donald J. Rose and others>[Abstract: symmetric positive-definite systems, fill calculation for a prescribed ordering, and algorithms for special orderings; publication information, June 1976]{ref}
Back to article page
