Gaussian Elimination

1. Overview

Gaussian elimination is the most fundamental direct method for solving a system of linear equations $A\mathbf{x} = \mathbf{b}$. It applies elementary row operations to the augmented matrix $[A|\mathbf{b}]$ to transform it 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 the upper triangular system $U\mathbf{x} = \mathbf{c}$, the solution is obtained starting from the last row:

4. Animated Computation Example

The following interactive widget demonstrates the complete process of forward elimination and back substitution on a concrete system of three linear equations, step by step.

Example Using Augmented Matrix

C++ implementation with partial pivoting

#include <vector>
#include <cmath>
#include <stdexcept>

std::vector<double> gaussianElimination(
    std::vector<std::vector<double>> A,
    std::vector<double> b)
{
    const int n = static_cast<int>(A.size());
    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]) < 1e-12)
            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];
        }
    }
    // 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 at or below row $k$ and at or to the right of column $k$, 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, in contrast to floating-point arithmetic, it is effective to choose the nonzero element with the smallest absolute value as the pivot. More precisely, by selecting the element whose height $H(p/q) = \max(|p|, |q|)$ (where $p/q$ is in lowest terms) is minimized as the pivot, the numerators and denominators of the elimination multipliers remain small, thereby suppressing coefficient swell.

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

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$. 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 practical matrices $\rho_n = O(n)$ almost always. 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)}| < \tau \cdot \|A\|$ (with $\tau$ on the order of the machine epsilon), 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).