LU Decomposition: Factoring Matrices into Lower and Upper Triangular Forms

LU decomposition is one of the most important matrix factorizations in numerical linear algebra. It expresses a square matrix A as the product of a lower triangular matrix L and an upper triangular matrix U, such that A = LU. This factorization transforms complex matrix operations into simpler triangular solves, dramatically improving computational efficiency for solving systems of equations, computing determinants, and finding matrix inverses.

The L matrix has ones on the diagonal and multipliers below it, while U is the result of Gaussian elimination — an upper triangular matrix with pivots on the diagonal. The decomposition essentially captures the entire elimination process in matrix form. Instead of performing row operations repeatedly, you factor once and then solve multiple right-hand sides cheaply.

The LU Decomposition Algorithm

The algorithm performs Gaussian elimination while recording the multipliers. For each column, you compute the multiplier for each row below the pivot: m = a[row][col] / a[pivot][pivot]. You store this multiplier in L[row][col] and subtract m times the pivot row from the current row in U. The diagonal of L is always 1, and U accumulates the elimination results.

For a 3×3 matrix, the decomposition proceeds column by column. First, you eliminate below the (1,1) pivot, recording two multipliers. Then you eliminate below the (2,2) pivot in the reduced matrix, recording one more multiplier. The result is L (with 1s on diagonal and multipliers below) and U (upper triangular).

Solving Systems with LU Decomposition

The primary advantage of LU decomposition is solving Ax = b. Once you have A = LU, you solve in two steps: first Ly = b (forward substitution), then Ux = y (back substitution). Both steps are O(n²) — far cheaper than Gaussian elimination from scratch for each new right-hand side. If you need to solve A x₁ = b₁, A x₂ = b₂, ..., you factor once and solve twice per system.

Forward substitution for Ly = b starts with y₁ = b₁ (since L[1,1] = 1), then computes y₂ = b₂ - L[2,1]y₁, and so on. Back substitution for Ux = y starts with xₙ = yₙ/U[n,n], then works upward. Both processes are straightforward and numerically stable.

When Pivoting Is Needed

In practice, the basic LU decomposition can fail if a pivot element is zero or very small. Partial pivoting addresses this by swapping rows to place the largest available element in the pivot position. The decomposition becomes PA = LU, where P is a permutation matrix recording the row swaps. This is the version used in most numerical libraries and is what our calculator implements for robustness.

Computational Advantages

• Solving multiple systems: Factor once, solve many times at O(n²) per system.
• Determinant: det(A) = det(L) × det(U) = product of U's diagonal elements.
• Inverse: Solve Ax = eᵢ for each unit vector to get columns of A⁻¹.
• Numerical stability: With partial pivoting, the algorithm is backward stable.
• Memory efficiency: L and U can overwrite A in place since L's diagonal is always 1.

Applications

LU decomposition is used in finite element analysis to solve large sparse systems, in circuit simulation (SPICE) to solve nodal equations, in computational fluid dynamics, and in any application requiring repeated solutions of linear systems with the same coefficient matrix. It is also the foundation for more advanced decompositions like Cholesky (for symmetric positive definite matrices) and QR decomposition.