LU Factorization with Partial Pivoting and Element Growth

Gaussian elimination is usually introduced as a sequence of row operations that turns a matrix into an upper triangular one. The description hides the only decision that matters numerically: which row to divide by at each step. Write the algorithm without that decision and it fails on matrices that are perfectly invertible, and it produces useless answers on matrices that are not.

The standard fix is partial pivoting, and its appeal is that it makes one clean promise about the entries of the computed factors. That promise is narrower than it first appears, and understanding exactly how narrow is the difference between trusting an LU solve and knowing when to check it.

Elimination written back into the matrix

The factorization computes a unit lower triangular L and an upper triangular U with PA = LU, where P is a permutation matrix recording the row interchanges. Nothing needs to be allocated for L beyond the space already occupied by A: the multipliers overwrite the entries that elimination drives to zero, and the diagonal of L is implicit.

At step k, the algorithm selects the row p with the largest absolute value in column k among rows k through n-1, swaps it into position k, and then eliminates the entries below the pivot using multipliers l_ik = a_ik / a_kk. Those multipliers are written into column k, and the trailing (n-k-1) x (n-k-1) block is updated by a rank-one subtraction.

#include <math.h>

/* In-place LU with partial pivoting on a row-major n x n matrix.
   On return A holds L below the diagonal and U on and above it.
   Returns the index of the first zero pivot, or n on success. */
int lu_factor(int n, double *A, int *piv)
{
    for (int k = 0; k < n; k++) {
        int p = k;
        for (int i = k + 1; i < n; i++)
            if (fabs(A[i * n + k]) > fabs(A[p * n + k]))
                p = i;
        piv[k] = p;
        if (A[p * n + k] == 0.0)
            return k;                       /* column is entirely zero below */
        if (p != k)
            for (int j = 0; j < n; j++) {   /* swap entire rows, not just the tail */
                double t = A[k * n + j];
                A[k * n + j] = A[p * n + j];
                A[p * n + j] = t;
            }
        double akk = A[k * n + k];
        for (int i = k + 1; i < n; i++) {
            double l = A[i * n + k] / akk;  /* |l| <= 1 follows from the search */
            A[i * n + k] = l;
            for (int j = k + 1; j < n; j++)
                A[i * n + j] -= l * A[k * n + j];
        }
    }
    return n;
}

The arithmetic is 2n^3/3 floating-point operations for the update, plus n^2/2 comparisons for the search. The comparisons are cheap per element but they are not free: they touch the same cache lines as the elimination, so the search roughly doubles the number of memory references in the unblocked form.

What the pivot search actually guarantees

Because the pivot is the largest magnitude entry in its column, every multiplier satisfies |l_ik| <= 1. That is the whole content of the guarantee, and it is worth being precise about what follows from it.

It follows that the entries of L are bounded by one in absolute value, so forming L introduces no growth of its own and the error in each multiplier is at most a couple of units of roundoff relative to the pivot. It does not follow that the entries of U stay comparable to those of A. The subtraction a_ij - l_ik * a_kj can cancel, and a small pivot at a late step divides a modest numerator by a small number, producing a large entry in U.

Consider the family of matrices with ones on the diagonal, minus one directly below each diagonal entry, and ones everywhere to the right. Elimination produces U entries that double at every step, reaching 2^(n-1) in the last column, even though every entry of the original matrix is plus or minus one. The factorization is exact in structure and catastrophic in scale: for n = 64, 2^63 exceeds the range in which the trailing entries can be added to anything else without losing all information about them.

The growth factor and what it bounds

The ratio of the largest entry appearing during elimination to the largest entry of the original matrix is the growth factor, usually written rho. The backward error of the computed factorization is proportional to it: the computed factors satisfy (A + E) = LU with a normwise bound on E of order n^2 * rho * u * ||A||, where u is the unit roundoff, about 1.1e-16 in double precision. The forward error on the solution then picks up the condition number of A on top of that.

So growth matters because it is the one factor in the error bound that the algorithm itself controls. Partial pivoting keeps rho small on almost every matrix that arises from a discretized problem, typically a small constant, but there is no proof that it must: the doubling example above is a legitimate matrix, and it forces rho = 2^(n-1) under partial pivoting. The bound is attained, it is just not representative.

Two alternatives change the guarantee rather than the constant. Complete pivoting searches the entire remaining submatrix for the largest entry and swaps both rows and columns; its worst-case growth bound is far below the exponential one, and no matrix is known that comes close to it, but the search costs O(n^3) comparisons and destroys the memory access pattern of the update. Rook pivoting searches alternately along a row and a column until the entry is maximal in both, which recovers most of the stability of complete pivoting at a cost closer to partial pivoting. Both are used when the matrix is small enough that the search is affordable, or when the application cannot tolerate a bad pivot sequence.

Panels, cache lines, and the transition to matrix-matrix kernels

The loop above is written as a sequence of rank-one updates, and rank-one updates are the worst possible shape for a modern processor. Each one reads a row of length n-k and a column of length n-k and performs one multiply-add per element, so every element is touched once to do two flops: the kernel runs at whatever rate memory can deliver operands, regardless of how fast the floating-point units are. The units sit idle waiting for data.

The fix is blocking. The factorization is split into a panel of nb columns, processed with the unblocked algorithm, followed by a single update of the remaining trailing matrix expressed as a matrix-matrix product of shape (n-k-nb) x nb times nb x (n-k-nb). That product has arithmetic intensity proportional to nb, and nb in the range of a few tens of columns is enough to move the kernel from memory-bound to compute-bound on typical cache hierarchies. LAPACK’s getrf is organized exactly this way, with the panel width chosen from the cache size at build time.

The blocking changes the error analysis only in the constants, but it changes the achievable rate substantially, and it is the reason the unblocked loop and a library call on the same matrix do not run at the same speed. The panel itself remains a sequence of rank-one updates and stays memory-bound; for large n the panel work is a small fraction of the total, so the overall cost approaches the matrix-matrix rate.

Where partial pivoting is not enough

Three situations defeat the standard scheme.

The first is a singular or nearly singular matrix. Partial pivoting cannot distinguish “the column is exactly zero” from “the column is small because the matrix is ill conditioned”. The factorization completes, the triangular solves produce a vector with enormous entries, and nothing in the solve reports a problem. Detecting this requires a separate estimate of the reciprocal condition number, computed from the factors with a handful of triangular solves, or an explicit rank test on the diagonal of U.

The second is a badly scaled matrix. Pivoting by absolute magnitude implicitly assumes that all rows are measured in the same units. A row whose entries are a thousand times larger than the rest will dominate the pivot choice regardless of its numerical importance, and the multipliers for the other rows can then exceed the bound that made the method attractive in the first place. Row equilibration before factorization, scaling each row by its largest entry, restores the meaning of the comparison; the permutation it produces is different, and usually better.

The third is the sparse case. Choosing a pivot for stability may create far more fill than a purely structural ordering would, and choosing an ordering for sparsity may put a tiny entry on the diagonal. Sparse direct solvers resolve this with threshold partial pivoting, accepting a candidate pivot if it is at least some fraction of the largest entry in its column, and deferring the rest; a matrix that is structurally symmetric and positive definite avoids the conflict entirely, which is why the Cholesky factorization described in Cholesky Factorization as a Positive-Definiteness Test needs no pivoting at all. The interaction between ordering and stability is taken up in Fill-In, Elimination Ordering, and Symbolic Factorization.

For a general dense matrix, partial pivoting with a condition estimate is the default for good reason: the bound on the multipliers is cheap, the growth is small on the overwhelming majority of inputs, and the cost of the search is a small fraction of the n^3 work. The failure modes are identifiable, and each has a documented remedy that costs either comparisons or a second factorization. What partial pivoting does not provide is a certificate, and any code that treats a completed factorization as evidence of a well-posed problem is relying on a property the algorithm never promised.