8. Solving Sparse Linear Systems
Society for Industrial and Applied Mathematics eBooks · 2006
Solving a sparse linear system Ax = b is now a matter of putting together the methods described in Chapters 2 through 7. The typical steps are to (1) find a permutation to reduce fill-in, (2) analyze and factorize the permuted matrix A (via Cholesky, LU, or QR), (3) permute the right-hand side b, (4) solve for x using forward/backsolve or by applying Householder reflections, (5) permute the solution x, and (6) optionally perform one or two steps of iterative refinement, which can sometimes improve the solution (see Problem 8.5), even when performed in standard floating-point precision. A naive method for solving Ax = b when A is square and nonsingular is to compute the inverse A−1, or x=inv(A)*b in MATLAB. This is numerically unstable when A is ill-conditioned in the dense case and also very costly in the sparse case. The inverse normally has no zero entries at all, as shown by the following two theorems. Theorem 8.1 (Gilbert [101]). Ignoring numerical cancellation, the nonzero pattern of the solution to Ax = b, where A has a zero-free diagonal, is = ReachA. Ignoring numerical cancellation, the solution x has no zero entries if A is strong Hall. Theorem 8.2 (Gilbert [101]). The transitive closure of the directed graph of A is the graph C, where i = ReachA(i). Ignoring numerical cancellation, gives the nonzero pattern of A−1. Every edge is present in C, and A−1 has no zero entries, if A is strong Hall. 8.1 Using a Cholesky factorization When A is symmetric positive definite, the system Ax = b can be solved via Cholesky factorization. If P is the fill-reducing permutation, LLT = PAPT. The system Ax = b becomes PAPTPx = Pb. Solving Ly = Pb for y, solving LTz = y for z, and finally x = PTz results in the solution x.