LU Decomposition

Factoring a matrix into lower and upper triangular pieces, the standard way to solve dense linear systems.

The factorization

LU decomposition writes a square matrix A as the product of a lower triangular matrix L (with ones on its diagonal) and an upper triangular matrix U. It is Gaussian elimination recorded as a factorization: U holds the result of elimination, and L stores the multipliers used to eliminate each entry. Most implementations include row swaps, giving PA = LU with a permutation matrix P.

Why factor at all

Once A = LU is known, solving Ax = b splits into two easy triangular solves: first solve L y = b by forward substitution, then U x = y by back substitution. The factorization costs on the order of n^3, but each triangular solve costs only n^2, so solving for many right-hand sides with the same A is cheap after one factorization.

Pivoting for stability

Without care, elimination can divide by a tiny pivot and amplify roundoff. Partial pivoting swaps rows so the largest available entry becomes the pivot, keeping the multipliers bounded and the process numerically stable. This is why library routines return a permutation along with L and U.

Determinant and inverse

The determinant of A is the product of the diagonal entries of U, times the sign of the permutation, which is the standard cubic-cost way to compute determinants. The inverse, when truly needed, is found by solving AX = I column by column using the same factorization.

import numpy as np
from scipy.linalg import lu_factor, lu_solve
A = np.array([[3.0, 1.0], [1.0, 2.0]])
lu, piv = lu_factor(A)
x = lu_solve((lu, piv), np.array([9.0, 8.0]))
print(x)   # solves Ax = b

Dense sub-blocks of the large operator matrices in engineering simulations are commonly solved with pivoted LU, while sparse and structured systems use tailored variants.