Direct methods solve finite-dimensional linear systems by finite algebraic transformations. Their analysis is the model case for backward stability.
Gaussian elimination
Definition. Gaussian elimination factors a matrix into triangular factors by eliminating subdiagonal entries. In exact arithmetic, solving \(Ax=b\) reduces to forward and backward substitution.
LU factorization
Proposition. If the leading pivots are nonzero, \(A=LU\) with \(L\) unit lower triangular and \(U\) upper triangular. With row pivoting, one obtains \(PA=LU\), which greatly improves numerical reliability.
Backward stability with pivoting
Theorem. Gaussian elimination with partial pivoting computes the exact solution of \((A+\Delta A)\widehat x=b\), where \(\|\Delta A\|\) is small relative to \(\|A\|\) up to a growth factor. In most practical problems this gives stable solutions.
Pivot growth
Remark. Partial pivoting can have large worst-case growth, although such cases are rare. Complete pivoting, rook pivoting, scaling, or structure-aware methods may be used when the matrix is especially delicate.
Sparse systems
Remark. For sparse matrices, the symbolic pattern of elimination matters as much as arithmetic. Fill-in control, ordering, and exploiting symmetry are numerical-analysis questions, not merely programming details.
Gaussian elimination and LU factorization
Elimination factors a matrix as \[ PA=LU, \] where \(P\) records row pivoting, \(L\) is unit lower triangular, and \(U\) is upper triangular. The factorization costs \(2n^3/3+O(n^2)\) flops; triangular solution costs \(2n^2+O(n)\). Reusing \(LU\) makes multiple right-hand sides cheap.
Why pivoting is needed
Without pivoting, a tiny pivot causes large multipliers and magnifies roundoff. Partial pivoting chooses the largest available magnitude in the current column, ensuring \(|l_{ij}|\le1\). It is backward stable in practice, though its worst-case growth factor can be \(2^{n-1}\).
Theorem (backward error, schematic). Floating-point Gaussian elimination with partial pivoting computes factors \(\widehat L,\widehat U\) satisfying \[ PA+\Delta A=\widehat L\widehat U, \qquad |\Delta A| \le\gamma_{cn}|\widehat L||\widehat U|. \] Consequently, \[ \|\Delta A\|_\infty \lesssim nu\,\rho\,\|A\|_\infty, \qquad \rho=\frac{\|\widehat U\|_\infty}{\|A\|_\infty}, \] up to modest factors.
The proof accumulates local dot-product errors in each Schur-complement update and reinterprets them as perturbations to the original entries. Growth \(\rho\), not only multiplier size, controls the final backward error.
Forward error and iterative refinement
If the solve is backward stable, \[ \frac{\|\widehat x-x\|}{\|x\|} \lesssim\kappa(A)u \] when \(\kappa(A)u\ll1\). Iterative refinement computes a residual \(r=b-A\widehat x\), solves \(Ad=r\) with the existing factors, and updates \(\widehat x\leftarrow\widehat x+d\). Computing the residual in higher precision can recover near-working-precision accuracy even when the factors are lower precision, provided conditioning is not too severe.
Special matrices
- Symmetric positive definite \(A\) admits Cholesky \(A=LL^\top\), at half the LU cost and without pivoting.
- Symmetric indefinite matrices need symmetric pivoting, such as \(P^\top AP=LDL^\top\) with \(1\times1\) and \(2\times2\) blocks.
- Banded matrices preserve bandwidth under suitable elimination and reduce cost to \(O(nb^2)\).
- Sparse elimination creates fill; graph orderings such as nested dissection control memory and work.
Determinants and inverses
After factorization, \[ \det A=\operatorname{sgn}(P)\prod_i u_{ii}. \] For magnitude stability, sum logarithms of pivots. One should almost never form \(A^{-1}\) to solve a system: triangular solves are cheaper and have a direct backward-error interpretation.