LU Decomposition with Partial Pivoting

The factorization

$\mathrm{LU}$ decomposition stores Gaussian elimination as two triangular matrices. With partial pivoting, a nonsingular square matrix $A \in \mathbb{R}^{n\times n}$ is factored as

The three factors have the following roles:

• $P$ is a permutation matrix that records row swaps
• $L$ is unit lower triangular
• $U$ is upper triangular

At elimination step $k$, partial pivoting selects the largest-magnitude entry in column $k$ at or below the diagonal. This avoids dividing by a zero pivot and usually limits numerical error growth.

For a permutation matrix, $P^{-1}=P^{\mathsf T}$. The equality $P^{-1}=P$ holds only for special permutations, such as a single row swap, and is not true for every accumulated permutation.

How elimination builds the factors

After selecting the pivot in column $k$, eliminate each entry below it with the multiplier

Then apply the following row update.

The multipliers are stored below the diagonal of $L$. If two rows are swapped after earlier elimination steps, the already-filled columns of $L$ must be swapped as well. Otherwise $PA=LU$ will no longer hold.

Factoring a Matrix with Partial Pivoting

We will track both row swaps so the final permutation agrees with the factors.

First pivot

The largest entry in the first column has magnitude $2$. Swap rows one and three, then eliminate the entry in row two with multiplier $\ell_{21}=1/2$. Row three already starts with zero.

Second pivot

The next pivot has magnitude $2$ in row three, so swap rows two and three. Swap their already stored first-column multipliers in $L$ as well. The multiplier for the final row is then $\ell_{32}=(1/2)/2=1/4$.

The resulting factors are

Multiplying the factors verifies the decomposition:

Solving a linear system

From $PA=LU$, the equation $Ax=b$ becomes

We solve this system in two stages using triangular matrices.

Forward substitution

First solve $Ly=Pb$. Here,

For a general lower triangular matrix,

The diagonal entries of this $L$ are all one.

Backward substitution

Then solve $Ux=y$ from the last row upward:

For the example,

Direct substitution confirms $Ax=b$.

Cost and reuse

Dense $\mathrm{LU}$ factorization requires approximately

This count gives the number of floating-point operations when multiplications and subtractions are counted separately. Once the factors are available, each new right-hand side requires only $O(n^2)$ work for the two substitutions. Reusing the same $\mathrm{LU}$ factors therefore avoids repeating the $O(n^3)$ factorization for every right-hand side.

The determinant also follows from the factors:

NumPy Implementation with Pivot Checks

import numpy as np


def plu_factor(matrix, tolerance=1e-12):
  """Return P, L, U such that P @ matrix == L @ U."""
  A = np.asarray(matrix)
  if np.iscomplexobj(A):
      raise ValueError("this example requires real entries")
  A = A.astype(float)
  if A.ndim != 2 or A.shape[0] != A.shape[1] or A.shape[0] == 0:
      raise ValueError("PLU factorization requires a nonempty square matrix")
  if not np.all(np.isfinite(A)):
      raise ValueError("matrix entries must be finite")
  if not np.isfinite(tolerance) or tolerance < 0:
      raise ValueError("relative tolerance must be finite and nonnegative")
  threshold = tolerance * np.linalg.norm(A, ord=np.inf)

  n = A.shape[0]
  U = A.copy()
  L = np.eye(n)
  P = np.eye(n)

  for k in range(n - 1):
      pivot = k + np.argmax(np.abs(U[k:, k]))
      if abs(U[pivot, k]) <= threshold:
          raise np.linalg.LinAlgError("matrix is singular to working precision")

      if pivot != k:
          U[[k, pivot]] = U[[pivot, k]]
          P[[k, pivot]] = P[[pivot, k]]
          L[[k, pivot], :k] = L[[pivot, k], :k]

      for i in range(k + 1, n):
          L[i, k] = U[i, k] / U[k, k]
          U[i, k:] -= L[i, k] * U[k, k:]

  if abs(U[-1, -1]) <= threshold:
      raise np.linalg.LinAlgError("matrix is singular to working precision")

  return P, L, U


def forward_substitution(L, right_hand_side):
  b = np.asarray(right_hand_side, dtype=float)
  y = np.empty_like(b)
  for i in range(len(b)):
      y[i] = (b[i] - L[i, :i] @ y[:i]) / L[i, i]
  return y


def backward_substitution(U, right_hand_side, tolerance=1e-12):
  y = np.asarray(right_hand_side, dtype=float)
  x = np.empty_like(y)
  threshold = tolerance * np.linalg.norm(U, ord=np.inf)
  for i in range(len(y) - 1, -1, -1):
      if abs(U[i, i]) <= threshold:
          raise np.linalg.LinAlgError("singular triangular system")
      x[i] = (y[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
  return x


def plu_solve(matrix, right_hand_side):
  P, L, U = plu_factor(matrix)
  b = np.asarray(right_hand_side)
  if np.iscomplexobj(b):
      raise ValueError("this example requires a real right-hand side")
  b = b.astype(float)
  if b.shape != (len(P),) or not np.all(np.isfinite(b)):
      raise ValueError("right-hand side must be a finite vector of matching size")
  y = forward_substitution(L, P @ b)
  x = backward_substitution(U, y)
  return x, P, L, U


A = np.array([[0, 2, 1], [1, 1, 0], [2, 1, 1]], dtype=float)
b = np.array([7, 3, 7], dtype=float)
x, P, L, U = plu_solve(A, b)

assert np.allclose(P @ A, L @ U)
assert np.allclose(A @ x, b)
assert np.allclose(x, np.array([1, 2, 3], dtype=float))

The checks are part of the example because factor order and row-swap bookkeeping are common implementation errors.

The pivot tolerance is relative to the matrix infinity norm. Uniformly rescaling a well-conditioned system therefore does not trigger failure just because every entry is small. The threshold is a numerical decision, not proof of exact singularity.

Scope of the method

This implementation solves systems only for nonsingular square matrices and treats a zero numerical pivot as a failure. Use $\mathrm{QR}$ for rectangular least-squares problems and column-pivoted $\mathrm{QR}$ or $\mathrm{SVD}$ for rank-deficient problems. A $\mathrm{PLU}$ factorization may still exist for a singular square matrix, but this solver cannot use it to obtain a unique solution.