Solving Linear Systems: Overview
Ax=b is the core of scientific computing; the choice between direct and iterative solvers depends on size, sparsity, and structure.
The central computation
Solving the linear system A x = b is the most common task in scientific computing, because discretizing differential equations, fitting models, and linearizing nonlinear problems all reduce to it. The right method depends on the size of A, whether it is sparse or dense, and whether it has structure such as symmetry or positive definiteness.
Direct methods
Direct methods compute the exact solution (up to rounding) in a finite number of steps, typically by factorizing A. Gaussian elimination and LU decomposition are the workhorses; Cholesky is a faster variant for symmetric positive-definite matrices. Direct methods are reliable but cost O(n^3) for dense matrices and can fill in sparse ones.
Iterative methods
Iterative methods generate a sequence of approximations that converge to the solution. Stationary methods (Jacobi, Gauss-Seidel) are simple; Krylov methods (conjugate gradient, GMRES) are far more powerful. Iterative methods shine for large sparse systems, where they avoid fill-in and can use the matrix only through matrix-vector products.
Choosing a solver
- Small or dense: direct factorization (LU, Cholesky).
- Large sparse and symmetric positive definite: conjugate gradient with a preconditioner.
- Large sparse and general: GMRES or BiCGSTAB with a preconditioner.
- Repeated solves with the same A: factor once, reuse the factorization.
The vast sparse systems from discretizing the field and transport equations of the breeder Hyperion are solved with preconditioned Krylov methods, where storing a dense factorization would be infeasible.