LU decomposition, step by step

Factoring a matrix into lower and upper triangular parts, why it is faster than repeated elimination, and how Cholesky specialises it.

LU decomposition writes a square matrix as the product of a lower triangular matrix and an upper triangular one: A = L·U. The factors come straight out of Gaussian elimination — U is the row echelon form, and L records the multipliers you used to get there.

Why bother

Because the factorisation is reusable. Solving Ax = b from scratch costs about n³/3 operations. Once you have L and U, each new right-hand side costs only n²: solve Ly = b by forward substitution, then Ux = y by back substitution. With dozens of right-hand sides that is a large win.

The determinant also becomes free — it is the product of the diagonal of U, times −1 for each row swap performed.

Pivoting

Not every matrix admits a plain A = L·U: a zero in a pivot position breaks it. The fix is a permutation matrix, giving P·A = L·U. In practice pivoting is used anyway, because choosing the largest available pivot keeps rounding errors under control.

Cholesky: the symmetric case

If A is symmetric and positive definite, it factors as A = L·Lᵀ — you only need one triangular factor. Cholesky is roughly twice as fast as LU and numerically very stable, which is why it dominates in optimisation, least squares and statistics (covariance matrices are symmetric positive definite by construction).

If the algorithm ever asks for the square root of a negative number, the matrix was not positive definite — Cholesky doubles as a test for that property.