Gaussian Elimination, in depth: row operations, partial pivoting, PA = LU, growth factors and exact elimination mod p and over GF(2)

From equations to a triangular system

Write the system Ax = b as an augmented matrix [A | b]. Three row operations leave the solution set unchanged: swapping two rows, multiplying a row by a non-zero constant, and adding a multiple of one row to another. Each is an invertible linear map applied to both sides, so any x that satisfied the old system satisfies the new one and vice versa.

Forward elimination processes columns left to right. At column k it picks a pivot row, and for every row i below it computes the multiplier m = a[i][k] / a[k][k] and subtracts m times the pivot row. After column k, every entry below the diagonal in that column is zero. After n columns the matrix is upper triangular, written U. Back substitution then solves the last equation for the last unknown, substitutes it into the row above, and climbs to the top.

Counting operations: eliminating column k touches roughly (n-k) rows of (n-k) entries, one multiply and one subtract each. Summed over k this is about 2n³/3 floating point operations. Back substitution costs about n². For n = 1,000 that is roughly 0.67 billion operations to factor and one million to substitute, which is why everything below is organised around doing the cubic part once.

Worked example with partial pivoting

Solve 2x + y - z = 8, -3x - y + 2z = -11, -2x + y + 2z = -3. In column 1 the candidates are 2, -3 and -2. Partial pivoting takes the largest magnitude, -3, so rows 1 and 2 swap. The multiplier for the new row 2 is 2 / -3 = -2/3, so we add 2/3 of row 1 to it, giving (0, 1/3, 1/3 | 2/3). For row 3 the multiplier is -2 / -3 = 2/3, giving (0, 5/3, 2/3 | 13/3).

In column 2 the candidates are now 1/3 and 5/3, so rows 2 and 3 swap and the pivot is 5/3. The multiplier for the remaining row is (1/3)/(5/3) = 1/5, which leaves (0, 0, 1/5 | -1/5). Back substitution gives z = -1, then (5/3)y + (2/3)(-1) = 13/3 so y = 3, then -3x - 3 - 2 = -11 so x = 2. Substituting (2, 3, -1) into the original three equations gives 8, -11 and -3: check every hand calculation this way.

Keep the multipliers. With the swaps recorded as a permutation P, they form the unit lower triangular L = [[1, 0, 0], [2/3, 1, 0], [-2/3, 1/5, 1]], and the final matrix is U = [[-3, -1, 2], [0, 5/3, 2/3], [0, 0, 1/5]]. Multiplying L by U reproduces the rows of A in pivot order (original rows 2, 3, 1). Notice that the multiplier -2/3 moved with its row when rows 2 and 3 swapped; forgetting that is the most common LU bookkeeping bug.

The algorithm in code

A direct implementation, with partial pivoting and a singularity test scaled to the matrix, so the threshold means the same thing for entries near 1e-6 and near 1e6:

def solve(A, b):
    n = len(A)
    M = [row[:] + [bi] for row, bi in zip(A, b)]      # augmented copy
    scale = max(abs(v) for row in A for v in row) or 1.0
    eps = 2.2e-16
    for k in range(n):
        p = max(range(k, n), key=lambda i: abs(M[i][k]))
        if abs(M[p][k]) <= n * eps * scale:
            raise ValueError(f"singular to working precision at column {k}")
        M[k], M[p] = M[p], M[k]
        for i in range(k + 1, n):
            m = M[i][k] / M[k][k]
            for j in range(k, n + 1):
                M[i][j] -= m * M[k][j]
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        s = sum(M[i][j] * x[j] for j in range(i + 1, n))
        x[i] = (M[i][n] - s) / M[i][i]
    return x

def residual(A, x, b):
    return max(abs(sum(a * xi for a, xi in zip(row, x)) - bi) for row, bi in zip(A, b))

Always check the residual ||Ax - b|| in tests. A small residual shows the algorithm did its job, not that x is accurate. For real work call LAPACK through NumPy or SciPy; this loop is far slower than a blocked factorisation.

Why pivoting matters: growth and conditioning

Take the system 1e-20 x + y = 1, x + y = 2. Its true solution is very close to x = 1, y = 1 and the matrix is well conditioned. Without pivoting, the multiplier is 1e20, the second row becomes (1 - 1e20) y = 2 - 1e20, and in double precision both sides round to -1e20, so y = 1 exactly. Back substitution then computes x = (1 - y) / 1e-20 = 0. The answer is completely wrong, not slightly wrong. Swapping rows first makes the multiplier 1e-20 and both unknowns come out as 1.

The general explanation is the growth factor: the ratio of the largest entry that appears during elimination to the largest entry of A. Rounding errors are proportional to the sizes of the numbers being subtracted, so large intermediate entries swamp the information in small ones. Partial pivoting keeps every multiplier at magnitude 1 or less, which bounds the growth factor by 2^n-1. That bound is attained: a matrix with ones on the diagonal and in the last column and -1 everywhere below the diagonal doubles its last column at every step. Running the elimination above on it gives growth 512 at n = 10 and about 5.4e8 at n = 30, which consumes most of double precision's sixteen digits.

Such matrices are rare in practice, so partial pivoting is LAPACK's default. Complete pivoting, which searches the whole remaining submatrix for the pivot, has a far smaller growth bound but costs an extra O(n³) comparisons and hurts memory locality, so it is reserved for rank-revealing uses.

Separate from algorithmic stability is conditioning. The condition number cond(A) measures how much the exact solution moves when A or b is perturbed. A backward-stable solver gives the exact answer to a nearby problem, so the forward error is roughly cond(A) times machine epsilon. The rule of thumb: with cond(A) near 1e10 in double precision, expect about six correct digits, however good the code. Check numpy.linalg.cond(A) (or a cheap LAPACK estimate) before trusting digits, and report it alongside results.

Factor once, solve many: PA = LU

Elimination does not depend on b except at the very end. Store the multipliers in the zeroed lower triangle, record the swaps, and you have the factorisation PA = LU at the cost of a single elimination. Each new right-hand side is then two triangular solves: Ly = Pb by forward substitution, Ux = y by back substitution, about 2n² operations instead of 2n³/3. A Newton solver that reuses its Jacobian for several steps, an implicit time-stepper, or a Kalman filter with fixed dynamics all exploit this.

import numpy as np
from scipy.linalg import lu_factor, lu_solve

A = np.array([[2., 1., -1.], [-3., -1., 2.], [-2., 1., 2.]])
lu, piv = lu_factor(A)                  # one O(n^3) factorisation (LAPACK getrf)
for b in ([8., -11., -3.], [1., 0., 0.]):
    x = lu_solve((lu, piv), np.array(b))  # O(n^2) per right-hand side (getrs)
    print(x, np.abs(A @ x - b).max())
print("cond:", np.linalg.cond(A))

Do not compute A^-1 and multiply. Forming the inverse costs about three times as much as the factorisation, is less accurate, and turns a sparse problem dense. The determinant is a free by-product: the product of U's diagonal times the sign of the permutation. Here that is (-3)(5/3)(1/5) = -1 times the sign of two swaps (+1), so det(A) = -1.

Exact elimination mod p and over GF(2)

Over a finite field there is no rounding, so any non-zero pivot will do. Modulo a prime p every non-zero element has an inverse (a^p-2 mod p, by Fermat), and reduced row-echelon form gives rank, solvability and a solution basis exactly, which is what linear codes and many contest problems need.

def rref_mod_p(M, p):
    M = [r[:] for r in M]
    rows, cols, r = len(M), len(M[0]), 0
    for c in range(cols):
        piv = next((i for i in range(r, rows) if M[i][c] % p), None)
        if piv is None:
            continue                      # free column
        M[r], M[piv] = M[piv], M[r]
        inv = pow(M[r][c], p - 2, p)
        M[r] = [v * inv % p for v in M[r]]
        for i in range(rows):
            if i != r and M[i][c] % p:
                f = M[i][c]
                M[i] = [(a - f * b) % p for a, b in zip(M[i], M[r])]
        r += 1
    return r, M                           # rank, reduced row-echelon form

Over GF(2) addition is XOR, so a row is a bitmask and eliminating a row is one instruction per 64 columns. The XOR linear basis from contest problems is this elimination done one row at a time: for each set bit from high to low, XOR out the basis vector with that leading bit, or, if none exists, store the vector there. The rank is the basis size.

How production libraries do it

Production solvers rearrange the same arithmetic for memory hierarchies. LAPACK's getrf factors a panel of columns with partial pivoting, then updates the trailing matrix with a single matrix multiply, so most flops run as cache-friendly BLAS-3 kernels. GPU libraries such as cuSOLVER offer the same getrf/getrs pair, and batched variants factor thousands of small matrices in one launch, which is how many small systems in a physics or robotics batch get solved.

Mixed precision is the modern twist: factor in a low precision that runs fast on tensor cores, then refine. Compute r = b - Ax in high precision, solve LU d = r with the cheap factors, set x += d, and repeat while the residual shrinks. It converges when cond(A) times the low-precision epsilon is well below 1; the HPL-MxP benchmark measures this pattern.

Sparse solvers reorder rows and columns first to limit fill-in, the non-zeros elimination creates. Symmetric positive definite matrices need no pivoting: use Cholesky, at half the cost.

Failure modes

• No pivoting. A tiny pivot produces huge multipliers and silent garbage, as in the 1e-20 example. Never ship elimination without at least partial pivoting.
• Absolute zero tests. Comparing pivots with == 0 in floating point misses near-singular matrices; a fixed 1e-9 breaks on badly scaled ones. Scale the tolerance by the matrix norm and size.
• Trusting a small residual. An ill-conditioned system can have a tiny residual and a solution with no correct digits. Report the condition estimate.
• Multipliers not swapped with rows. L is wrong while U looks right. Test several b vectors against a reference solver.
• Fraction blow-up. Exact rational elimination can grow numerators exponentially; use modular arithmetic or the fraction-free Bareiss variant for exact integer results.
• Forgetting the field. Modular code with a composite modulus has elements with no inverse; Fermat's inverse silently returns the wrong value.

Trade-offs

Related reading

Keep going: least squares with QR versus normal equations is covered in Linear Regression, in depth; pivoting reappears as the core step of the simplex algorithm; repeated matrix products are the subject of Matrix Exponentiation, in depth; and another O(n log n) restructuring of linear algebra is the Fast Fourier Transform.

What to do next

1. Work the 3 by 3 example above by hand, then verify it with scipy.linalg.lu_factor.
2. Run the code on the 1e-20 system with and without the row swap and watch x change from 1 to 0.
3. Build the 30 by 30 growth matrix, solve it with your code, and compare the error with numpy.linalg.solve.
4. In any code you own, replace inv(A) @ b with a factor-and-solve, reusing the factorisation across right-hand sides.
5. Add a condition estimate and a residual check to the outputs of your solver calls.
6. Implement the GF(2) basis insert and use it to compute the rank of 1,000 random 64-bit vectors.