Prerequisite:

Row Operations and Echelon Forms

The key insight: the solution set of $A\mathbf{x} = \mathbf{b}$ is unchanged by elementary row operations on the augmented matrix $[A \mid \mathbf{b}]$:

  1. Swap two rows: $R_i \leftrightarrow R_j$
  2. Scale a row by a nonzero scalar: $R_i \leftarrow \alpha R_i$, $\alpha \neq 0$
  3. Add a multiple of one row to another: $R_i \leftarrow R_i + \alpha R_j$, $i \neq j$

Each operation preserves the row space of $A$ (operations 1 and 3 exactly; operation 2 up to scaling) and does not change the solution set of the system.

Definition (Row Echelon Form). A matrix is in row echelon form (REF) if:

  • All zero rows are at the bottom.
  • The leading entry (pivot) of each nonzero row is strictly to the right of the pivot in the row above.

Definition (Reduced Row Echelon Form). A matrix is in RREF if additionally:

  • Every pivot equals $1$.
  • Every pivot is the only nonzero entry in its column.

Theorem (Existence and Uniqueness of RREF). Every matrix $A$ can be reduced to a unique RREF via elementary row operations.

Existence is constructive: Gaussian elimination (below). Uniqueness is proved by showing that the pivot positions are invariants of the row space - any two RREF matrices with the same row space must be identical.

The Algorithm

Forward Elimination (Gaussian)

Start with the augmented matrix $[A \mid \mathbf{b}]$. Process columns left to right:

  1. Find the leftmost column with a nonzero entry at or below the current pivot row.
  2. Swap to bring a nonzero entry to the pivot position.
  3. Scale the pivot row to make the pivot $1$ (for RREF; skip for REF).
  4. Eliminate all entries below the pivot using row operation 3.
  5. Move to the next pivot row and repeat.

Back Substitution

After reaching REF, solve for the leading variables from the bottom up, expressing each in terms of the free variables (columns without pivots).

Step-by-Step Example

Solve the $3 \times 3$ system:

$$\begin{pmatrix} 2 & 1 & -1 \ -3 & -1 & 2 \ -2 & 1 & 2 \end{pmatrix} \mathbf{x} = \begin{pmatrix} 8 \ -11 \ -3 \end{pmatrix}$$

Augmented matrix:
[ 2   1  -1 |  8 ]
[-3  -1   2 | -11]
[-2   1   2 | -3 ]

Step 1: pivot on column 1, row 1. Pivot = 2.
R2 <- R2 + (3/2)R1:   -3 + 3 = 0, -1 + 3/2 = 1/2, 2 - 3/2 = 1/2, -11 + 12 = 1
R3 <- R3 + (1)R1:     -2 + 2 = 0,  1 + 1  = 2,  2 - 1 = 1,  -3 + 8 = 5

[ 2   1  -1 |  8 ]
[ 0  1/2 1/2|  1 ]
[ 0   2   1 |  5 ]

Step 2: pivot on column 2, row 2. Pivot = 1/2.
R3 <- R3 - 4*R2:   0, 2 - 2 = 0, 1 - 2 = -1, 5 - 4 = 1

[ 2   1  -1 |  8 ]
[ 0  1/2 1/2|  1 ]
[ 0   0  -1 |  1 ]

This is REF. Back-substitute:
  x3 = -1
  x2 = (1 - (1/2)(-1)) / (1/2) = (3/2)/(1/2) = 3
  x1 = (8 - 1*(3) - (-1)*(-1)) / 2 = (8 - 3 - 1)/2 = 2

Solution: x = (2, 3, -1).

Elementary Matrices and LU Decomposition

Each elementary row operation on $A$ corresponds to left-multiplication by an elementary matrix $E$:

  • Swap $R_i \leftrightarrow R_j$: $E$ is the identity with rows $i$ and $j$ swapped (a permutation matrix).
  • Scale $R_i$ by $\alpha$: $E$ is the identity with $E_{ii} = \alpha$.
  • $R_i \leftarrow R_i + \alpha R_j$: $E$ is the identity with $E_{ij} = \alpha$.

All elementary matrices are invertible; their inverses are elementary matrices of the same type (with the operation reversed).

LU Decomposition. Forward elimination without row swaps applies a sequence of lower-triangular elementary matrices $E_1, E_2, \ldots, E_k$ to produce an upper-triangular matrix $U$:

$$E_k \cdots E_2 E_1 A = U$$

Setting $L = (E_k \cdots E_1)^{-1} = E_1^{-1} \cdots E_k^{-1}$, we get $A = LU$ where $L$ is unit lower-triangular (ones on the diagonal) and $U$ is upper-triangular. The entries of $L$ below the diagonal are exactly the multipliers $-\alpha$ used during elimination, placed in the corresponding positions.

Solving $A\mathbf{x} = \mathbf{b}$ via LU requires two triangular solves:

  1. $L\mathbf{y} = \mathbf{b}$ (forward substitution, $O(n^2)$)
  2. $U\mathbf{x} = \mathbf{y}$ (back substitution, $O(n^2)$)

After factoring once ($O(n^3)$), each new right-hand side costs only $O(n^2)$.

Pivoting for Numerical Stability

In exact arithmetic, any nonzero pivot works. In floating-point arithmetic, small pivots amplify rounding errors catastrophically.

Partial pivoting. Before eliminating column $j$, swap the current row with the row having the largest absolute value in column $j$ at or below the diagonal. This keeps multipliers $|m_{ij}| \leq 1$, bounding error growth.

Complete pivoting. Swap both rows and columns to place the global maximum remaining entry in the pivot position. More stable than partial pivoting, but more expensive (requires searching the whole submatrix). In practice, partial pivoting suffices.

With partial pivoting, the factorization becomes $PA = LU$ where $P$ is a permutation matrix recording the row swaps. LAPACK’s dgesv implements exactly this.

Complexity: Counting FLOPs

At step $k$ of forward elimination (eliminating column $k$), we perform:

  • $(n - k)$ divisions to compute multipliers
  • $(n - k)^2$ multiply-add pairs to update the submatrix

Total FLOPs for forward elimination:

$$\sum_{k=1}^{n-1} (n-k) + (n-k)^2 \approx \sum_{j=1}^{n-1} j^2 = \frac{n(n-1)(2n-1)}{6} \approx \frac{n^3}{3}$$

Back substitution adds $O(n^2)$ FLOPs, negligible by comparison. So Gaussian elimination costs $\frac{2n^3}{3}$ FLOPs (counting multiply and add separately), confirming the $O(n^3)$ complexity.

Existence, Uniqueness, and Rank via RREF

After reducing $[A \mid \mathbf{b}]$ to RREF, let $r$ be the number of pivot rows.

Existence of a solution. The system is consistent iff no pivot appears in the last column, i.e., RREF has no row of the form $[0\ 0\ \cdots\ 0 \mid c]$ with $c \neq 0$. Equivalently, $\text{rank}(A) = \text{rank}([A \mid \mathbf{b}])$.

Uniqueness of the solution. Given consistency, the solution is unique iff there are no free variables, i.e., $r = n$ ($A$ has full column rank).

Four cases for $A \in \mathbb{R}^{m \times n}$:

$r = m$, $r = n$ $r = m$, $r < n$ $r < m$, $r = n$ $r < m$, $r < n$
Unique solution for every $\mathbf{b}$ Infinitely many solutions for every $\mathbf{b}$ Unique solution or no solution No solution or infinitely many

Computing rank. The rank of $A$ equals the number of pivot columns in its RREF. Equivalently, it is the number of nonzero rows in any REF of $A$.

Examples

Gaussian elimination (with partial pivoting) is the most important algorithm in scientific computing. Its implementation in LAPACK (Linear Algebra PACKage), via dgetrf (LU factorization) and dgetrs (triangular solve), underlies essentially all dense linear system solvers, including those in NumPy (numpy.linalg.solve), SciPy, MATLAB’s \ operator, and Julia’s \.

LAPACK achieves near-peak floating-point performance by reorganizing the algorithm into block LU factorization: instead of eliminating one column at a time, it factors rectangular blocks, exposing level-3 BLAS operations (matrix-matrix multiply) that are highly cache-efficient. The block size is tuned to the cache hierarchy.

For sparse matrices (most of the $n^2$ entries are zero), direct sparse factorization methods (e.g., SuperLU, UMFPACK) permute rows and columns to minimize fill-in - nonzeros created during elimination in positions that were originally zero. For very large sparse systems, iterative methods (CG, GMRES) become preferable and avoid factorization entirely.

Read Next: