Gaussian Elimination — When Ill-Conditioned Matrices Lie

Gaussian elimination is the foundational algorithm for solving linear systems Ax = b. It runs in O(n³) and, with partial pivoting, is numerically stable for most practical problems. It is what numpy.linalg.solve, scipy.linalg.solve, and every numerical linear algebra library use under the hood (via LU decomposition).

Understanding Gaussian elimination means understanding why partial pivoting matters — without it, small pivots cause numerical catastrophe. It also means knowing how LU decomposition factorises A for efficient multiple right-hand-side solving, and why direct methods (Gauss) are preferred over iterative methods (Jacobi, Gauss-Seidel) for dense systems. Here's the thing: most engineers who reach for np.linalg.solve don't realise they're paying O(n³) every time. Factor once, solve many.

Forward Elimination with Partial Pivoting

Forward elimination reduces the augmented matrix [A|b] to upper triangular form. At each column, we first find the row with the largest absolute value in that column (partial pivoting) and swap it to the diagonal. This ensures all multipliers are ≤ 1, bounding error amplification. Then we subtract multiples of the pivot row from rows below to zero out the column. The process produces an upper triangular matrix U and a transformed right-hand side vector.

When the pivot is < 1e-12 after pivoting, the matrix is numerically singular and we raise an error. In production codes, the threshold is adjusted based on the matrix norm. Don't hardcode 1e-12 — it's a toy value. Real libraries scale the tolerance by the maximum element magnitude. If you're implementing this yourself, and you're in production, you're doing it wrong. Use LAPACK.

def gaussian_elimination(A, b):
    """Solve Ax = b using Gaussian elimination with partial pivoting."""
    import copy
    n = len(A)
    # Augmented matrix [A|b]
    M = [list(A[i]) + [b[i]] for i in range(n)]
    for col in range(n):
        # Partial pivoting: swap row with maximum pivot element
        max_row = max(range(col, n), key=lambda r: abs(M[r][col]))
        M[col], M[max_row] = M[max_row], M[col]
        pivot = M[col][col]
        if abs(pivot) < 1e-12:
            raise ValueError('Matrix is singular or nearly singular')
        # Eliminate below pivot
        for row in range(col+1, n):
            factor = M[row][col] / pivot
            for j in range(col, n+1):
                M[row][j] -= factor * M[col][j]
    # Back substitution
    x = [0.0] * n
    for i in range(n-1, -1, -1):
        x[i] = M[i][n]
        for j in range(i+1, n):
            x[i] -= M[i][j] * x[j]
        x[i] /= M[i][i]
    return x

A = [[2,1,-1],[3,3,1],[1,-1,2]]
b = [8,13,1]
solution = gaussian_elimination(A, b)
print(f'x = {[round(v,4) for v in solution]}')
# Verify with numpy
import numpy as np
print(f'numpy: {np.linalg.solve(A, b)}')

Why Partial Pivoting Matters

Without pivoting, Gaussian elimination fails when the pivot element is zero (division by zero) and becomes numerically unstable when the pivot is very small (amplifying rounding errors).

Example: A = [[0.001, 1],[1, 1]], b = [1, 2]. Without pivoting, the first pivot is 0.001 — the elimination multiplier is 1/0.001 = 1000. Floating point errors in the first row get amplified by 1000x. Partial pivoting swaps to use the row with the largest absolute value as pivot, keeping multipliers ≤ 1 and bounding error amplification. That's the difference between a correct answer and garbage. And don't think 'my matrix is well-conditioned' — it only takes one bad pivot to corrupt everything.

Back Substitution — Solving the Upper Triangular System

After forward elimination, we have an upper triangular system Ux = c. We solve from the last equation upward: x_n = c_n / U_{nn}, then for i from n-1 down to 1: x_i = (c_i - sum_{j=i+1}^{n} U_{ij} x_j) / U_{ii}. This is O(n²) — much cheaper than elimination. The numerical stability of back substitution is good because it solves a triangular system directly; however, a very small diagonal element in U can still cause large errors. Watch for tiny diagonals — they hint at singularity or extreme ill-conditioning.

def back_substitute(U, c):
    n = len(U)
    x = [0.0] * n
    for i in reversed(range(n)):
        s = c[i]
        for j in range(i+1, n):
            s -= U[i][j] * x[j]
        if abs(U[i][i]) < 1e-12:
            raise ValueError('Zero diagonal in U — singular matrix')
        x[i] = s / U[i][i]
    return x

# Example: after elimination from earlier
U = [[2, 1, -1], [0, 1.5, 2.5], [0, 0, 4]]
c = [8, 9, -4]
x = back_substitute(U, c)
print(x)  # [2.0, 3.0, -1.0]

LU Decomposition — Reusing the Factorization

Gaussian elimination actually computes an LU decomposition of A (with partial pivoting: PA = LU). The lower triangular matrix L contains the multipliers, and U is the upper triangular result. Once we have P, L, U, solving Ax = b becomes solving Ly = Pb (forward substitution) and Ux = y (back substitution). The forward elimination step (O(n³)) is done once; each new right-hand side costs only O(n²). This is how library solvers work internally.

When you call numpy.linalg.solve, it computes the LU decomposition and then solves. You can get the factors explicitly via scipy.linalg.lu. The real win: if you need to solve for 1000 different b vectors, you factor once and solve 1000 times. That's O(n³ + 1000n²) vs O(1000n³). The savings are massive. Don't be the engineer who calls solve in a loop.

import scipy.linalg
import numpy as np

A = np.array([[2,1,-1],[3,3,1],[1,-1,2]])
b1 = np.array([8,13,1])
b2 = np.array([5,10,0])

# LU decomposition
P, L, U = scipy.linalg.lu(A)
print('L:\n', L)
print('U:\n', U)

# Solve for b1 using the factors
# First solve L y = P @ b1
Pb1 = P.T @ b1  # or P @ b1 depending on convention
y = scipy.linalg.solve_triangular(L, Pb1, lower=True)
x1 = scipy.linalg.solve_triangular(U, y, lower=False)
print('Solution for b1:', x1)

# Solve for b2 without re-factoring
Pb2 = P.T @ b2
y = scipy.linalg.solve_triangular(L, Pb2, lower=True)
x2 = scipy.linalg.solve_triangular(U, y, lower=False)
print('Solution for b2:', x2)

Computational Complexity and When to Avoid Gaussian Elimination

Dense Gaussian elimination costs O(n³) flops. For n=1000, that's about 1 billion operations — fine on modern hardware. For n=10,000, a trillion operations takes minutes. Libraries like LAPACK are highly optimised (blocked algorithms, BLAS level 3) but the cubic growth remains.

For large sparse systems, iterative solvers (Conjugate Gradient, GMRES) exploit sparsity and converge in O(nnz * iterations) where nnz is number of nonzeros and iterations are often << n. The break-even point for dense vs sparse is usually around n=1000-5000, depending on sparsity. If your matrix is 99% zeros, Gaussian elimination ignores that and chews up memory storing zeros. Don't do that. Check sparsity first.

When Gaussian Elimination Fails: Ill-Conditioned Systems and Iterative Refinement

Even with partial pivoting, Gaussian elimination can produce inaccurate results for ill-conditioned systems. The condition number κ(A) = ||A|| ||A⁻¹|| measures how much errors in b or A amplify in the solution. When κ(A) ~ 10^t, you lose about t digits of precision. For double precision (16 digits), if κ=10^12, only 4 digits remain reliable.

Iterative refinement addresses this: after computing a first solution x₀ by Gaussian elimination, compute the residual r = b - Ax₀, solve Ad = r (using the same LU factors), and correct x₁ = x₀ + d. Repeat until the residual is small. This recovers full precision if the system is not too ill-conditioned. For extremely ill-conditioned systems (κ > 10^14), even iterative refinement fails, and you need regularisation (e.g., Tikhonov) or arbitrary-precision arithmetic.

In production, always compute the condition number before trusting a solve. LAPACK's DGESVX (expert driver) provides error bounds and iterative refinement automatically.

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

def iterative_refinement(A, b, max_iter=10, tol=1e-12):
    """Solve Ax = b with iterative refinement."""
    lu, piv = lu_factor(A)
    x = lu_solve((lu, piv), b)
    for i in range(max_iter):
        r = b - A @ x
        if np.linalg.norm(r) / np.linalg.norm(b) < tol:
            break
        dx = lu_solve((lu, piv), r)
        x += dx
    return x

# Example with an ill-conditioned Hilbert matrix
n = 10
A = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
cond = np.linalg.cond(A)
print(f'Condition number: {cond:.2e}')
b = np.ones(n)
x_direct = np.linalg.solve(A, b)
x_refined = iterative_refinement(A, b)
print('Direct residual:', np.linalg.norm(b - A @ x_direct) / np.linalg.norm(b))
print('Refined residual:', np.linalg.norm(b - A @ x_refined) / np.linalg.norm(b))