Gaussian Elimination

Gaussian Elimination is a foundational direct algorithm in numerical linear algebra and computation and algorithms for solving systems of linear algebraic equations Ax = b, determining the rank of a matrix, evaluating determinants, and computing the inverse of an invertible square matrix. Named in honor of German mathematician Carl Friedrich Gauss—who employed row reduction in 1809 to solve normal equations arising from celestial least squares orbit determination of the asteroid Pallas, although an equivalent tabular column-elimination procedure appeared in Chapter 8 (Fangcheng) of the Chinese mathematical classic The Nine Chapters on the Mathematical Art (Jiuzhang Suanshu, circa 150 BCE)—the algorithm systematically applies elementary row operations to reduce an augmented matrix [A | b] to upper-triangular row-echelon form.

In modern scientific computing, Gaussian elimination provides the algorithmic backbone of dense linear solvers in LAPACK and BLAS, corresponding algebraically to the triangular matrix factorization PA = LU via LU decomposition. For an n × n dense matrix A, forward elimination requires 2n³/3 + O(n²) floating-point arithmetic operations (flops), followed by n² flops for backward substitution. When paired with partial row pivoting or complete row-and-column pivoting to bound element growth factors in IEEE 754 floating-point arithmetic, Gaussian elimination achieves backward numerical stability across general non-singular systems where specialized symmetric solvers such as Cholesky decomposition or iterative Krylov methods such as the conjugate gradient method do not directly apply.

Forward Elimination and Elementary Elimination Matrices

Consider a system of n linear equations in n unknowns written in matrix form Ax = b, where A ∈ ℝ^n×n has entries a_ij (1 ≤ i, j ≤ n), x = [x₁, x₂, ..., x_n]^T ∈ ℝⁿ is the unknown column vector, and b = [b₁, b₂, ..., b_n]^T ∈ ℝⁿ is the right-hand-side vector. Gaussian elimination relies on three elementary row operations that preserve the solution set of the linear system:

1. Row Swapping (Interchange): Exchanging row i and row k (R_i ↔ R_k) to move a nonzero pivot element onto the principal diagonal.
2. Row Scaling: Multiplying row i by a nonzero scalar α ≠ 0 (R_i ← α*R_i*).
3. Row Addition (Elimination): Adding a scalar multiple of row k to row i (R_i ← R_i − m_ikR_k).

The forward elimination phase proceeds through n − 1 sequential stages indexed by k = 1, 2, ..., n − 1. Denoting the initial matrix and vector as A⁽¹⁾ = A and b⁽¹⁾ = b, at step k the diagonal entry a_kk^(k) ≠ 0 serves as the pivot. For every row i > k below the pivot row, the algorithm computes the Gauss multiplier:

Subtracting m_ik times row k from row i zeroes out all subdiagonal entries in column k (a_ik^(k+1) = 0 for i > k) and updates the remaining active submatrix entries (j = k + 1, ..., n) and right-hand-side components via the rank-1 outer-product update:

In matrix notation, step k left-multiplies A^(k) by a unit lower-triangular Gauss transformation matrix (Frobenius matrix) M_k = I − m_ke_k^T, where m_*k = [0, ..., 0, m_k+1,k, ..., m_nk]^T. After n − 1 stages, the system is transformed into upper-triangular form Ux = y, where U = A⁽ⁿ⁾ = M_n−1 ··· M₂M₁A and L = M₁⁻¹M₂⁻¹ ··· M_n−1⁻¹ is the unit lower-triangular factor storing the multipliers m_ik directly below the main diagonal (A = LU).

Back-Substitution and Gauss–Jordan Reduction

Once forward elimination produces the upper-triangular system Ux = y with diagonal pivots u_kk = a_kk^(k) ≠ 0, the unknown components x_n, x_n−1, ..., x₁ are recovered in reverse order by backward substitution:

A closely related variant introduced by Wilhelm Jordan in 1887, Gauss–Jordan elimination, eliminates entries both below and above each pivot a_kk^(k) and scales each pivot row by 1/a_kk^(k) to produce the reduced row-echelon form (RREF) [I | A⁻¹b]. However, clearing superdiagonal entries increases the arithmetic cost of solving a single linear system from 2n³/3 flops to n³ flops (a 50% overhead), so standard Gaussian elimination with back-substitution is preferred for solving Ax = b, while Gauss–Jordan reduction is reserved for explicit matrix inversion [A | I] → [I | A⁻¹] and canonical subspace basis extraction.

Floating-Point Operation Count and Determinant Evaluation

At elimination step k (1 ≤ k ≤ n − 1), operating on the trailing (n − k) × (n − k) active submatrix requires n − k divisions to compute the multipliers m_ik, followed by (n − k)² multiplications and (n − k)² subtractions to update a_ij^(k+1). Summing over k = 1 to n − 1 yields the leading cubic operation count:

Updating the right-hand-side vector b during forward elimination requires ∑_k=1..n−1 2(n − k) = n(n − 1) flops, and backward substitution requires n² flops, giving a combined cost of (2/3)n³ + O(n²) flops. As an immediate algebraic byproduct, the determinant of A equals the product of the upper-triangular diagonal pivots adjusted by the sign of the row permutation P (s = number of row swaps):

Partial, Complete, and Rook Pivoting in Finite-Precision Arithmetic

Naive Gaussian elimination without row exchanges breaks down whenever a pivot a_kk^(k) = 0 is encountered, and suffers severe numerical instability in finite-precision IEEE 754 arithmetic whenever |a_kk^(k)| is small relative to other entries in column k. For example, on the 2 × 2 system with a₁₁ = ε = 10⁻¹⁷, a₁₂ = 1, a₂₁ = 1, a₂₂ = 1, multiplier m₂₁ = 10¹⁷ causes a₂₂⁽²⁾ = 1 − 10¹⁷ to round to −10¹⁷ in double precision (machine epsilon ε_mach ≈ 1.11 × 10⁻¹⁶), yielding a 100% erroneous solution x₁ = 0 instead of x₁ ≈ 1.

James H. Wilkinson's backward error analysis (1961) established that the computed solution x̂ satisfies a perturbed system (A + ΔA)*x̂* = b with relative backward error bounded by the element growth factor ρ_n = max_i,j,k |a_ij^(k)| / max_i,j |a_ij|:

To keep the multipliers bounded (|m_ik| ≤ 1) and restrain ρ_n, practical implementations employ three pivoting strategies:

• Partial Pivoting (GEPP, PA = LU): Searches column k on and below the diagonal (i ≥ k) for the entry of maximum magnitude |a_pk^(k)| = max_k≤i≤n |a_ik^(k)| and swaps rows k and p in O(n) comparisons per step (O(n²) total search overhead). Although Wilkinson constructed pathological matrices where ρ_n = 2ⁿ⁻¹ under partial pivoting, average-case growth on physical and random matrices is ρ_n ≈ O(n^1/2 to n^2/3).
• Complete Pivoting (GECP, PAQ = LU): Searches the entire (n − k + 1) × (n − k + 1) trailing submatrix for the maximum entry and swaps both rows and columns in O(n³) comparisons, bounding ρ_n ≤ n^1/2 (2 · 3^1/2 ··· n^1/(n−1))^1/2 ≈ O(n^{(1/2) ln n}).
• Rook Pivoting: Alternates row and column scans until an entry is simultaneously maximal in its row and column, achieving complete-pivoting stability in O(n²) expected search overhead.

On modern hierarchical memory architectures, unblocked row-by-row Gaussian elimination is bottlenecked by main-memory bandwidth because each Level-2 BLAS rank-1 update (dger) performs only 2 flops per memory reference. High-performance libraries such as LAPACK (dgetrf), OpenBLAS, and Intel MKL partition the matrix into column panels of width b ≈ 64–256. After factoring a narrow n × b panel using unblocked partial pivoting, the algorithm applies Level-3 BLAS triangular solves (dtrsm) and matrix-matrix multiplications (dgemm) to update the trailing (n − b) × (n − b) submatrix at ~2b flops per memory byte—reaching over 85% of peak hardware floating-point throughput and reducible asymptotically via the Strassen algorithm.

For large sparse matrices arising from finite-element PDE discretizations or circuit simulations, naive elimination suffers from fill-in, where zero entries a_ij = 0 become nonzero during the outer-product update a_ij − m_ika_kj. Sparse direct solvers (SuperLU, UMFPACK, MUMPS, PARDISO) permute rows and columns prior to elimination using nested dissection or minimum-degree graph ordering on the elimination graph—or defer to matrix-free methods like the conjugate gradient method or singular value decomposition for rank-deficient systems—to minimize fill-in and preserve sparse complexity.

Why is Gaussian elimination with partial pivoting preferred over computing the explicit matrix inverse A⁻¹?

Computing A⁻¹ explicitly via Gauss–Jordan elimination or solving n unit-vector right-hand sides requires 2n³ flops—three times the (2/3)n³ cost of an LU factorization—and yields a larger forward residual ‖*Ax̂* − b‖₂. Once PA = LU is computed, any number of right-hand-side vectors b can be solved in just 2n² flops via forward and backward triangular substitution.

When does Gaussian elimination succeed without pivoting?

Gaussian elimination is guaranteed to proceed without zero pivots and maintain bounded element growth (ρ_n ≤ 2) without row exchanges when A is symmetric positive-definite (where Cholesky decomposition halves the work to n³/3 flops) or row/column diagonally dominant (|a_ii| > ∑_j≠i |a_ij| for all i).

How does Gaussian elimination determine the rank and nullspace of a rectangular matrix?

Applying Gaussian elimination with complete or rook pivoting to an m × n matrix A produces a row-echelon form containing r nonzero pivot rows and m − r zero rows. The number of nonzero pivots r equals rank(A), the pivot columns span the column space Range(A), and setting the n − r free variables corresponding to non-pivot columns yields a basis for the nullspace Null(A) of dimension n − r by the rank-nullity theorem.

1. Wilkinson, J. H. (1961). "Error analysis of direct methods of matrix inversion." Journal of the ACM, 8(3), 281–330. https://doi.org/10.1145/321075.321076
2. Trefethen, L. N., & Bau, D. (1997). "Growth factors of Gaussian elimination with partial pivoting." In Numerical Linear Algebra (Lectures 20–22, pp. 147–171). SIAM. https://doi.org/10.1137/1.9780898719574
3. Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms (2nd ed., Chapter 9: LU Factorization and Linear Equations). Society for Industrial and Applied Mathematics (SIAM). https://doi.org/10.1137/1.9780898718027
4. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed., Chapter 3: General Linear Systems). Johns Hopkins University Press. https://doi.org/10.56021/9781421407944
5. Grcar, J. F. (2011). "Mathematicians of Gaussian elimination." Notices of the American Mathematical Society, 58(6), 782–792. https://www.ams.org/notices/201106/rtx110600782p.pdf