How Does Partial Pivoting Prevent Tiny Pivots in LU?

Gaussian elimination with partial pivoting prevents division by tiny numbers during LU decomposition by reordering the rows of a matrix before each elimination step. At every stage, the algorithm scans the current column from the diagonal downward, finds the entry with the largest absolute value, and swaps its row into the pivot position. This simple row exchange keeps multipliers bounded, limits floating-point rounding errors, and produces a stable factorization of a permuted matrix.

Why Tiny Pivots Cause Numerical Instability

In classical computing, numbers are stored in finite-precision floating-point formats such as IEEE 754 double precision. During standard Gaussian elimination without pivoting, the algorithm divides each subdiagonal entry \(a_{ik}\) by the diagonal pivot \(a_{kk}\) to compute the multiplier:

\[m_{ik} = \frac{a_{ik}}{a_{kk}}\]

If \(a_{kk}\) is zero, division by zero halts the algorithm. Even if \(a_{kk}\) is nonzero but very small compared to \(a_{ik}\), the multiplier \(m_{ik}\) becomes huge. When the algorithm updates the remaining submatrix using

\[a_{ij} \leftarrow a_{ij} - m_{ik} a_{kj}\]

the large product \(m_{ik} a_{kj}\) can overwhelm \(a_{ij}\). In finite-precision arithmetic, adding or subtracting numbers of vastly different magnitudes causes significant digits of \(a_{ij}\) to be lost, a phenomenon known as catastrophic cancellation or element growth.

How Partial Pivoting Works

Partial pivoting avoids small pivots by selecting the largest available entry in the current column at each step \(k\) (for \(k = 1, \dots, n-1\)). Specifically, it performs the following operations:

  1. Search the entries \(a_{ik}\) for \(i = k, k+1, \dots, n\) to find the index \(p\) that maximizes \(|a_{ik}|\).
  2. Swap row \(k\) and row \(p\) in the working matrix.
  3. Record the row swap in a permutation matrix \(P\) (or a compact pivot index vector).
  4. Compute the multipliers \(m_{ik} = a_{ik} / a_{kk}\) for \(i > k\) and eliminate the entries below the diagonal.

Because \(|a_{kk}| \ge |a_{ik}|\) for all \(i > k\) after the swap, every multiplier satisfies:

\[|m_{ik}| \le 1\]

Bounding the multipliers by \(1\) prevents them from amplifying rounding errors during the elimination updates.

The Resulting PA = LU Factorization

Swapping rows during elimination means the algorithm no longer factors the original matrix \(A\) directly. Instead, it computes an LU decomposition of a row-permuted version of \(A\):

\[PA = LU\]

Here, \(P\) is a permutation matrix whose rows are the identity matrix reordered according to the swaps, \(L\) is a unit lower triangular matrix containing the multipliers \(m_{ik}\) (all with magnitude at most \(1\)), and \(U\) is an upper triangular matrix.

When solving a linear system \(Ax = b\), the permutation is applied to the right-hand side vector \(b\) as well:

\[PAx = Pb \implies LUx = Pb\]

The system is then solved in two fast triangular stages: forward substitution for \(Ly = Pb\), followed by backward substitution for \(Ux = y\).

Computational Cost and Practical Reliability

Partial pivoting adds only \(O(n^2)\) comparisons across the entire \(O(n^3)\) elimination process, making its overhead negligible for large dense matrices. While complete pivoting (searching both rows and columns) offers slightly tighter theoretical bounds on element growth, partial pivoting is the standard choice in classical numerical libraries such as LAPACK because it is much faster and achieves excellent numerical stability in practice.