Direct Methods for Linear Systems

Direct methods solve a linear system $Ax = b$ in a finite number of operations. The cornerstone is Gaussian elimination, which factors $A$ into triangular matrices. Understanding pivoting strategies and operation counts is essential for practical implementation.

Gaussian Elimination and LU Factorization

Definition5.1LU Decomposition

An LU decomposition of $A \in \mathbb{R}^{n \times n}$ is a factorization $A = LU$ where $L$ is unit lower triangular ( $\ell_{ii} = 1$ ) and $U$ is upper triangular. The system $Ax = b$ reduces to solving $Ly = b$ (forward substitution, $O(n^2)$ ) then $Ux = y$ (back substitution, $O(n^2)$ ). The factorization costs $\frac{2}{3}n^3 + O(n^2)$ flops.

Definition5.2Partial and Complete Pivoting

Partial pivoting selects the largest entry in the current column below the diagonal as pivot: $PA = LU$ where $P$ is a permutation matrix. Complete pivoting selects the largest entry in the entire remaining submatrix: $PAQ = LU$ . Partial pivoting suffices in practice and guarantees $|l_{ij}| \leq 1$ , giving a growth factor $g = \max_{ij} |u_{ij}| / \max_{ij} |a_{ij}|$ bounded by $2^{n-1}$ (though $g = O(n^{2/3})$ in practice).

Cholesky Factorization

Definition5.3Cholesky Decomposition

For a symmetric positive definite (SPD) matrix $A$ , the Cholesky factorization is $A = LL^T$ where $L$ is lower triangular with positive diagonal entries. The algorithm computes $\ell_{jj} = \sqrt{a_{jj} - \sum_{k=1}^{j-1} \ell_{jk}^2}$ and $\ell_{ij} = \frac{1}{\ell_{jj}}\left(a_{ij} - \sum_{k=1}^{j-1} \ell_{ik}\ell_{jk}\right)$ for $i > j$ . The cost is $\frac{1}{3}n^3$ flops, half that of LU. No pivoting is needed: the factorization is always numerically stable for SPD matrices.

ExampleCholesky Example

For $A = \begin{pmatrix} 4 & 2 & 2 \\ 2 & 5 & 1 \\ 2 & 1 & 6 \end{pmatrix}$ : $\ell_{11} = 2$ , $\ell_{21} = 1$ , $\ell_{31} = 1$ , $\ell_{22} = \sqrt{5 - 1} = 2$ , $\ell_{32} = (1 - 1)/2 = 0$ , $\ell_{33} = \sqrt{6 - 1 - 0} = \sqrt{5}$ . Thus $L = \begin{pmatrix} 2 & 0 & 0 \\ 1 & 2 & 0 \\ 1 & 0 & \sqrt{5} \end{pmatrix}$ .

Banded and Sparse Systems

RemarkExploiting Sparsity

For banded matrices with bandwidth $m \ll n$ , LU factorization costs $O(m^2 n)$ instead of $O(n^3)$ . Sparse direct solvers use fill-reducing orderings (minimum degree, nested dissection) to minimize the fill-in (new nonzeros created during factorization). Nested dissection achieves $O(n^{3/2})$ fill-in for 2D problems and $O(n^2)$ for 3D, compared to $O(n^2)$ and $O(n^{4/3})$ respectively for the dense case.

ExampleThomas Algorithm for Tridiagonal Systems

A tridiagonal system $a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i$ is solved in $O(n)$ operations by the Thomas algorithm (specialized LU without pivoting). Forward sweep: $c_i' = c_i / (b_i - a_i c_{i-1}')$ , $d_i' = (d_i - a_i d_{i-1}') / (b_i - a_i c_{i-1}')$ . Back substitution: $x_n = d_n'$ , $x_i = d_i' - c_i' x_{i+1}$ . This arises in cubic spline interpolation and finite difference methods.