Gaussian Elimination

1. Overview

When the coefficient matrix $A$ is square and nonsingular, Gaussian elimination is the most fundamental direct method for finding the unique solution of $A\mathbf{x} = \mathbf{b}$ (in the singular case, the result of elimination reveals whether there is no solution or infinitely many). It applies elementary row operations to the augmented matrix $[A|\mathbf{b}]$ to transform its coefficient part into upper triangular form (forward elimination), then obtains the solution by back substitution.

Named after C. F. Gauss (1777–1855), similar techniques can be found in the ancient Chinese text The Nine Chapters on the Mathematical Art (Jiuzhang Suanshu), dating back to before the common era.

2. Forward Elimination

From the $n \times n$ coefficient matrix $A$ and the right-hand side vector $\mathbf{b}$, form the augmented matrix. For each column $k$ ($k = 1, 2, \ldots, n-1$), subtract a constant multiple of row $k$ from all rows below row $k$ to eliminate $a_{ik}$.

3. Back Substitution

Once forward elimination yields an upper triangular system $U\mathbf{x} = \mathbf{c}$ whose diagonal entries $u_{kk}$ are all nonzero, the solution is obtained starting from the last row:

4. Animated Computation Example

Separately from Figures 1 and 2 (a 5-variable system), the following interactive widget demonstrates the complete process of forward elimination and back substitution on a smaller, concrete system of three linear equations, step by step.

Example Using Augmented Matrix

C++ implementation with partial pivoting

#include <vector>
#include <cmath>
#include <stdexcept>
#include <limits>
#include <algorithm>

std::vector<double> gaussianElimination(
    std::vector<std::vector<double>> A,
    std::vector<double> b)
{
    const int n = static_cast<int>(A.size());
    // singularity threshold proportional to the largest entry (scale-invariant)
    double scale = 0.0;
    for (const auto& row : A) for (double v : row) scale = std::max(scale, std::abs(v));
    const double tol = 100 * std::numeric_limits<double>::epsilon() * scale;
    for (int k = 0; k < n - 1; ++k) {
        // partial pivot: find the row with max |A[i][k]|
        int pivot = k;
        for (int i = k + 1; i < n; ++i)
            if (std::abs(A[i][k]) > std::abs(A[pivot][k])) pivot = i;
        if (pivot != k) { std::swap(A[k], A[pivot]); std::swap(b[k], b[pivot]); }
        if (std::abs(A[k][k]) <= tol)
            throw std::runtime_error("singular matrix");
        // eliminate
        for (int i = k + 1; i < n; ++i) {
            double m = A[i][k] / A[k][k];
            for (int j = k; j < n; ++j) A[i][j] -= m * A[k][j];
            b[i] -= m * b[k];
        }
    }
    if (std::abs(A[n - 1][n - 1]) <= tol)   // check the last pivot as well
        throw std::runtime_error("singular matrix");
    // back substitution
    std::vector<double> x(n);
    for (int i = n - 1; i >= 0; --i) {
        double s = b[i];
        for (int j = i + 1; j < n; ++j) s -= A[i][j] * x[j];
        x[i] = s / A[i][i];
    }
    return x;
}

5. Pivot Selection

At step $k$ of forward elimination, if the pivot $a_{kk}^{(k)}$ is zero or close to zero, division by zero or large rounding errors may occur. To avoid this, pivot selection (pivoting) is performed.

5.1 Partial Pivoting

In column $k$, the element with the largest absolute value at or below row $k$ is selected, and the corresponding row swap is performed. Since swapping rows merely reorders the equations, the correspondence between variables is unaffected. In practice, partial pivoting is sufficient in most cases and is adopted by standard numerical libraries such as LAPACK.

5.2 Complete Pivoting (Full Pivoting)

In complete pivoting, at step $k$, the element with the largest absolute value is selected from all entries of the coefficient matrix at or below row $k$ and at or to the right of column $k$ (the right-hand-side column is not included), and both rows and columns are swapped.

A column swap corresponds to a reordering of variables. For example, swapping columns 2 and 3 interchanges the roles of $x_2$ and $x_3$. Therefore, after all elimination steps are complete, the column swap history must be traced in reverse order to restore the original variable ordering.

5.3 Pivot Selection for Rational Arithmetic

The pivoting strategies described above assume floating-point arithmetic. When performing Gaussian elimination exactly using rational numbers (fractions with arbitrary-precision integers), there is no rounding error; instead, coefficient swell becomes the concern. During elimination, the number of digits in the numerators and denominators can grow explosively, making computation extremely slow.

To mitigate this problem, unlike in floating-point arithmetic, one can choose the pivot by looking at the size of the numerator and denominator rather than at the magnitude (absolute value). Choosing a nonzero element whose height $H(p/q) = \max(|p|, |q|)$ (where $p/q$ is in lowest terms) is small tends to keep the numerators and denominators of that step's multipliers small. This is only a local heuristic, however, and it does not guarantee that the overall coefficient swell is minimized. For example, $1/1000$ has a small absolute value but a height of $1000$, whereas $2/3$ has a larger absolute value but a height of only $3$; the smallest-height rule therefore picks $2/3$ (large in absolute value but small in height), which is the opposite of the largest-absolute-value rule used in floating point.

6. Computational Complexity

The computational complexity of Gaussian elimination for a system of $n$ linear equations is as follows.

• Forward Elimination: $\dfrac{2}{3}n^3 + O(n^2)$ floating-point operations
• Back Substitution: $n^2 + O(n)$ floating-point operations
• Overall: $O(n^3)$

7. Relationship with LU Decomposition

When elimination proceeds without row swaps (that is, when all leading principal minors of $A$ are nonzero), the forward elimination process of Gaussian elimination is equivalent to decomposing the matrix $A$ into the product of a lower triangular matrix $L$ and an upper triangular matrix $U$:

The elimination multipliers $m_{ik}$ become the subdiagonal elements of $L$, and the resulting upper triangular matrix is $U$. When partial pivoting swaps rows, a permutation matrix $P$ gives $PA = LU$; when complete pivoting also swaps columns, $PAQ = LU$ (on a row swap, the previously computed multipliers are reordered along with the rows). Once the $LU$ decomposition is computed, a solution for any different right-hand side $\mathbf{b}$ can be obtained in $O(n^2)$ operations, making it efficient when solving with the same coefficient matrix but multiple right-hand sides.

8. Pitfalls

• Growth factor and stability: The numerical stability of Gaussian elimination with partial pivoting depends on the growth factor $\rho_n = \max_{i,j,k} |a_{ij}^{(k)}| / \max_{i,j} |a_{ij}|$. The theoretical worst case is $\rho_n \le 2^{n-1}$ (attained by the Wilkinson matrix), but for most practical and random matrices it is empirically small (typically at most of order $n$), though this is not guaranteed. Complete pivoting grows more slowly, with $\rho_n \le n^{1/2} (2 \cdot 3^{1/2} \cdots n^{1/(n-1)})^{1/2}$.
• Scaling sensitivity: When rows have widely differing scales, the largest-absolute-value criterion alone can pick a poor pivot. Some implementations adopt "scaled partial pivoting", which normalizes each row by $\|A_{i,:}\|_\infty$ before applying the criterion.
• Detecting singularity: When $|a_{kk}^{(k)}| \le \tau \cdot \max_{i,j} |a_{ij}|$ (with $\tau$ a few tens to 100 times the machine epsilon, as in the implementation notes above; a rule of thumb to be chosen according to dimension, scaling, and required accuracy), the matrix is numerically singular and the solver should switch to a least-squares or SVD-based path.
• Conditioning limits accuracy: The computed solution $\hat{\mathbf{x}}$ has relative error on the order of $\kappa(A) \cdot u$ (with $u$ the unit roundoff); when $\kappa(A) \approx 1/u$, all significant digits are lost. Estimating the condition number first (e.g., LAPACK's ?gecon) is recommended.
• Iterative refinement: Compute the residual $\mathbf{r} = \mathbf{b} - A\hat{\mathbf{x}}$ in higher precision, solve $A\mathbf{d} = \mathbf{r}$ with the existing $LU$, and apply $\hat{\mathbf{x}} \leftarrow \hat{\mathbf{x}} + \mathbf{d}$ one or two times. This recovers accuracy whenever $\kappa(A) \cdot u \ll 1$.
• Dense-matrix assumption: Plain Gaussian elimination targets dense matrices. Applying it directly to a sparse matrix causes fill-in (zeros becoming nonzeros during elimination) and destroys sparsity.

9. Frequently Asked Questions

10. References

• Wikipedia, "Gaussian elimination"
• Wikipedia, "LU decomposition"
• G. H. Golub & C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins, 2013 (Chapter 3 on $LU$ factorization; §3.4 on growth factors and the Wilkinson matrix).
• L. N. Trefethen & D. Bau III, Numerical Linear Algebra, SIAM, 1997 (Lectures 20-23 on $LU$ and partial pivoting).
• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002 (§9 on backward error analysis of Gaussian elimination; §12 on iterative refinement).