Gaussian Elimination With Partial Pivoting

A tiny pivot, like 0.0003, in a matrix division can amplify a rounding error by a factor of 10,000. Gaussian elimination with partial pivoting fixes that by swapping rows so the largest available entry in the column becomes the pivot, and that single change turns an unstable calculation into a reliable one for most classroom-sized systems.

You already know the basic elimination steps: forward elimination to reach row-echelon form, then back-substitution. What makes a solver trustworthy in floating-point arithmetic is the pivoting strategy, not the elimination itself. Partial pivoting is the standard choice taught in numerical methods courses and used by LAPACK's dgesv routine for dense matrices up to at least 5×5 by hand. For anything larger, you call the library and let it handle the pivot decisions.

What You Will Find Here

This covers why small pivots break an otherwise correct algorithm, the exact partial pivoting procedure with pseudocode, a worked example that shows the difference with and without pivoting, what complete and scaled pivoting offer, and how production software like LAPACK actually implements it. If you need a refresher on basic Gaussian elimination steps or what Gaussian elimination is as a method, those topics have their own treatment elsewhere.

Gaussian elimination without pivoting fails numerically on perfectly invertible matrices. The cause is floating-point round-off error, which is harmless when pivot values are comparable in magnitude to other entries but catastrophic when a pivot is much smaller. A division by 0.001 when the next row's entry is 1.0 multiplies the error by 1000.

Consider a 2×2 system: 0.001x + y = 1, x + y = 2. The exact solution is x ≈ 1.000, y ≈ 0.999. Without pivoting, the first pivot is 0.001. The elimination multiplier is 1/0.001 = 1000. Subtracting 1000 times the first row from the second produces (1-1000)y ≈ 2-1000, or, -999y ≈ -998, so y ≈ 0.999, x ≈ 1.000.0. The computed y = 1.0, then x = (1-1)/0.001 = 0. The correct x is 1.0; the computed x is 0. That is a 100% relative error from one small pivot.

The problem is the choice of pivot. Partial pivoting swaps rows so the largest entry in the column, which is 1.0 from the second row, becomes the pivot. That swap eliminates the division by 0.001 entirely and the computed result matches the exact one.

This example is from Burden & Faires, Numerical Analysis (10th ed., Section 2.5 on pivoting strategies). The same principle applies to any n×n system: a pivot value that is orders of magnitude smaller than other entries in its column causes unbounded growth of rounding error. Partial pivoting bounds that growth to a factor of at most 2^{n-1} in the worst case (Golub & Van Loan, Matrix Computations, 4th ed., Ch. 3), which for a 5×5 hand-checked system is acceptable.

Partial pivoting modifies the forward elimination step. At each column k, before eliminating below the diagonal, search rows k through n for the entry with the largest absolute value. Swap that row with the current row k. Then proceed with the same elimination multiplier formula.

The algorithm in pseudocode, for an augmented matrix [A|b] of size n×(n+1):

for k = 1 to n-1
  // Partial pivoting
  max_val = |a(k,k)|
  max_row = k
  for i = k+1 to n
    if |a(i,k)| > max_val
      max_val = |a(i,k)|
      max_row = i
    end if
  end for
  if max_val == 0
    stop: system is singular or no unique solution
  end if
  // Swap rows k and max_row in A and b
  for j = k to n
    swap a(k,j) with a(max_row,j)
  end for
  swap b(k) with b(max_row)
  // Elimination
  for i = k+1 to n
    factor = a(i,k) / a(k,k)
    for j = k to n
      a(i,j) = a(i,j) - factor * a(k,j)
    end for
    b(i) = b(i) - factor * b(k)
  end for
end for

After forward elimination, the matrix is in row-echelon form (REF). Back-substitution starts from the last row and solves upward. The pivot array, which records which rows were swapped, is what LAPACK's dgesv returns as an INTEGER array of length n, where IPIV(i) tells you that row i was swapped with row IPIV(i) during the factorization PA = LU. The P matrix is the permutation matrix built from those swaps.

A zero pivot at any stage forces a row swap. If the entire column below the current row is zero, the matrix is singular and the system either has no solution or infinite solutions, which is handled on the solution types page.

Work the same system two ways to see the difference. Use 3×3 with exact fractions to show the pivot choices, then simulate 4-digit floating point on the round-off case.

System and Exact Solution

System:
x + 2y + 3z = 14
4x + 5y + 6z = 32
7x + 8y + 10z = 53
Exact solution: x=1, y=2, z=3.

Without pivoting (exact fractions): First pivot = 1. Eliminate below: subtract 4×row1 from row2, subtract 7×row1 from row3. New matrix: row2 = [0, -3, -6 | -24], row3 = [0, -6, -11 | -45]. Second pivot = -3. Subtract 2×row2 from row3: row3 = [0, 0, 1 | 3]. Back-substitute: z=3, y=(-24 + 18)/(-3)=2, x=14-4-9=1. Works exactly.

With partial pivoting (exact fractions): At column 1, largest |entry| in column 1 is 7 from row3. Swap row1 and row3: new row1 = [7, 8, 10 | 53], row2 = [4, 5, 6 | 32], row3 = [1, 2, 3 | 14]. Pivot = 7. Eliminate: factor from row2 = 4/7, subtract from row2; factor from row3 = 1/7, subtract from row3. Row2 becomes [0, 3/7, 2/7 | 12/7], row3 becomes [0, 6/7, 11/7 | 45/7]. At column 2, largest |entry| in column 2 below row1 is 6/7 from row3. Swap row2 and row3: new row2 = [0, 6/7, 11/7 | 45/7], new row3 = [0, 3/7, 2/7 | 12/7]. Pivot = 6/7. Eliminate row3: factor = (3/7)/(6/7)=1/2. Row3 becomes [0, 0, -7/14 | -21/14] → [0,0,-1/2 | -3/2]. Back-substitute: z=3, y= (45/7-33/7)/(6/7)= (12/7)/(6/7)=2, x=(53-16-30)/7=1. Same exact solution, but different intermediate numbers. The pivot choices changed the elimination order but not the final result in exact arithmetic.

Round-Off Demonstration

Round-off demonstration (4-digit floating point): Use the same system but with pivot 0.001 in a different matrix. Let matrix be [0.001, 1, 1; 1, 1, 0] with b = [1, 2]. Without pivoting: pivot=0.001, multiplier=1000. Subtract 1000×row1 from row2: row2 becomes [0,, 999,, 1000 |, 998]. In 4-digit:, 999 is, 999.0,, 1000 is, 1000.0,, 998 is, 998.0. Back-substitute: y = (, 998)/(, 1000)=0.998, x = (1-0.998)/0.001 = 2.0. Exact y=0.999, x=1.001. Error in x: 100%. With partial pivoting: swap row2 up, pivot=1. Multiplier=0.001. Eliminate: row2 becomes [0, 0.999,, 0.001 | 0.998]. Back-substitute: y=0.998/0.999=0.999, x=(2-0.999)/1=1.001. Error ~0.1%. The swap cut the error by a factor of 1000.

This matches the analysis in Golub & Van Loan, Ch. 3: partial pivoting limits the growth factor to at most 2^{n-1}, which for n=2 is 2, while no pivoting allows unbounded growth from a single small leading entry.

Partial pivoting swaps only rows. Complete pivoting searches the entire remaining submatrix (rows k..n, columns k..n) for the largest entry, then swaps both rows and columns to bring that entry to the pivot position. Column swaps change the order of variables, so you must track the permutation. Complete pivoting reduces the theoretical growth factor to about n^{1/2} for random matrices (Golub & Van Loan), but it doubles the search cost per pivot from O(n-k) to O((n-k)^2). For classroom 5×5 systems the extra cost is negligible, but for the matrices LAPACK solves (up to thousands of rows) the extra search is not worth the marginal stability gain. Burden & Faires note that complete pivoting is rarely used in practice for that reason.

Scaled pivoting, also called partial pivoting with implicit scaling, normalises each row by its largest entry before comparing pivot candidates. It handles cases where rows have wildly different magnitudes, for example, one row in meters and another in kilometers. The search still looks only down the column, but the comparison uses scaled values. Scaled pivoting is a refinement on partial pivoting, not a separate strategy. Most numerical methods courses cover it as an extension in the same section on pivoting strategies.

LAPACK's dgesv routine solves a real general system A x = b using LU decomposition with partial pivoting. The call does three things in one shot: factor PA = LU, solve L y = P b by forward substitution, then solve U x = y by back-substitution. The pivot information comes back as an INTEGER array IPIV of length n. Each entry IPIV(i) is 1-indexed: it tells you that during the factorization, row i was swapped with row IPIV(i). To reconstruct the permutation matrix P, start with the identity and swap rows according to IPIV in order.

dgesv was first released in 1992 and has been the standard dense linear system solver in LAPACK ever since. MATLAB's backslash operator and NumPy's numpy.linalg.solve both call dgesv (or its parallel variant) under the hood. They do not expose the pivot array in the default path, but you can access it via the lower-level ''lu'' functions in both environments.

For a 5×5 system, dgesv completes the factorization in about 140 floating-point operations, most of which are elimination. The pivot decisions add about 10% overhead compared to no pivoting, but that 10% buys numerical stability for nearly all matrices that are not deliberately pathological. Wilkinson's matrix, a 12×12 triangular matrix with small off-diagonal entries, can cause partial pivoting to fail with a growth factor of 2^{n-1}, but for random matrices the growth factor averages around 10 (Golub & Van Loan).

The LAPACK documentation recommends dgesv for dense matrices up to about 10^4 rows on modern hardware. Beyond that, or for sparse matrices, you switch to dedicated solvers in scipy.sparse.linalg.splu. For a 5×5 system by hand, you can replicate the exact dgesv algorithm with the pseudocode above and check your fractions against its output.

This method suits numerical-methods students who need to understand why a textbook algorithm fails on a computer. It also suits linear algebra undergraduates who solve 3×3 and 4×4 systems by hand and want to verify that their row swaps produce the correct REF. Engineering students in a first matrix course who encounter the difference between unique, infinite, and no solution can use partial pivoting as the mechanism that reveals a singular matrix: a zero pivot that cannot be swapped means the matrix is singular, and they then check the augmented column to distinguish no solution from infinite solutions.

For solving a system of 1000 equations or more, a different approach is needed. For that scale, call a dedicated sparse solver library, scipy.sparse.linalg.splu or the equivalent in your environment, and let it handle pivoting internally. You should also skip it if you want a deep theoretical treatment of vector spaces, eigenvalues, or the LU decomposition proof, which Strang and Axler cover.

The single thing that most often goes wrong: a student swaps rows correctly but forgets to apply the same swap to the constant vector b, then back-substitutes using the wrong right-hand side. Every row operation must apply to the entire augmented row, or the solution set changes. Check the augmented column after every swap before you write the next elimination step.