The factorization
decomposition stores Gaussian elimination as two triangular matrices. With partial pivoting, a nonsingular square matrix is factored as
Nakafa
Free learning materials.For AI agents: use /llms.txt for the Nakafa content index.
The three factors have the following roles:
At elimination step , partial pivoting selects the largest-magnitude entry in column at or below the diagonal. This avoids dividing by a zero pivot and usually limits numerical error growth.
For a permutation matrix, . The equality holds only for special permutations, such as a single row swap, and is not true for every accumulated permutation.
After selecting the pivot in column , eliminate each entry below it with the multiplier
Then apply the following row update.
The multipliers are stored below the diagonal of . If two rows are swapped after earlier elimination steps, the already-filled columns of must be swapped as well. Otherwise will no longer hold.
We will track both row swaps so the final permutation agrees with the factors.
The largest entry in the first column has magnitude . Swap rows one and three, then eliminate the entry in row two with multiplier . Row three already starts with zero.
The next pivot has magnitude in row three, so swap rows two and three. Swap their already stored first-column multipliers in as well. The multiplier for the final row is then .
The resulting factors are
Multiplying the factors verifies the decomposition:
From , the equation becomes
We solve this system in two stages using triangular matrices.
First solve . Here,
For a general lower triangular matrix,
The diagonal entries of this are all one.
Then solve from the last row upward:
For the example,
Direct substitution confirms .
Dense 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 work for the two substitutions. Reusing the same factors therefore avoids repeating the factorization for every right-hand side.
The determinant also follows from the factors:
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.
This implementation solves systems only for nonsingular square matrices and treats a zero numerical pivot as a failure. Use for rectangular least-squares problems and column-pivoted or for rank-deficient problems. A factorization may still exist for a singular square matrix, but this solver cannot use it to obtain a unique solution.
Published: . Updated: .
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))