LU Decomposition

What is LU Decomposition?

Suppose you need to solve the linear system Ax = b not once, but dozens or hundreds of times — each time with a different right-hand side vector b but the same coefficient matrix A. This situation arises constantly in engineering: finite-element solvers assemble the same stiffness matrix K and then solve Ku = f for many different load vectors f; circuit simulators build a conductance matrix G once and solve Gv = i for many current vectors i.

The naive approach — run Gaussian elimination from scratch each time — is wasteful. All the work of reducing A to row echelon form is identical regardless of b. LU decomposition captures that work once and for all by factoring A into a product of two triangular matrices.

The key insight is that the process of Gaussian elimination on A can be recorded as two triangular factors: a lower triangular matrix L whose diagonal entries are all 1, and an upper triangular matrix U that is the row echelon form of A. Once A = LU is computed, any system Ax = b is solved in two cheap triangular-system steps.

The L and U Factors

A lower triangular matrix L has all its nonzero entries on or below the main diagonal. The special case relevant to LU decomposition is a unit lower triangular matrix — one whose diagonal entries are all exactly 1, with arbitrary values below the diagonal and zeros above. An upper triangular matrix U has all its nonzero entries on or above the main diagonal. The LU factorization writes

The structure is elegant: L encodes what was done to reduce A, and U records the result. Together they preserve all the arithmetic of Gaussian elimination in a form that can be reused.

Connection to Gaussian Elimination

Gaussian elimination proceeds by choosing a pivot in each column and subtracting multiples of the pivot row from all rows below it. The multiplier used to eliminate the entry in row i using the pivot in row k is

m_ik = a_ik / a_kk

Ordinarily these multipliers are discarded once the elimination is done. The LU factorization saves them: the multiplier m_ik is stored in position (i, k) of L — the very slot that was zeroed out during elimination. Every elimination step below the diagonal fills in one entry of L. The remaining reduced matrix is U.

To see why this works, note that subtracting m_ik times row k from row i is equivalent to left-multiplying by an elementary matrix E_ik. Performing all such operations gives E · A = U, where E is the product of elementary matrices. The inverse of E is L: the lower triangular matrix whose off-diagonal entries are exactly the multipliers m_ik with their signs flipped back. Because the inverse of a product of lower triangular matrices is itself lower triangular, L inherits clean triangular structure.

In short: U is what Gaussian elimination produces; L remembers how it got there.

Solving Ax = b via Forward and Back Substitution

Once the factorization A = LU is in hand, solving Ax = b becomes a two-step process. Substituting A = LU into Ax = b gives LUx = b. Introducing the intermediate vector y = Ux splits this into two triangular systems:

Forward Substitution in Detail

In forward substitution, the unknowns are solved from top to bottom. Since L is unit lower triangular, the first equation is simply y₁ = b₁. The second equation is l₂₁ y₁ + y₂ = b₂, giving y₂ = b₂ − l₂₁ y₁. In general, at row i we have

yᵢ = bᵢ − Σⱼ₌₁ⁱ⁻¹ lᵢⱼ yⱼ

This costs O(n²) operations total. Back substitution for Ux = y proceeds symmetrically from the bottom up, again O(n²).

Why LU? Efficiency for Multiple Right-Hand Sides

The whole point of LU decomposition is the separation of costs. Factorizing A = LU requires approximately n³/3 multiplications and additions — a one-time O(n³) investment. Once L and U are stored, each new right-hand side b requires only a forward substitution and a back substitution, each costing O(n²).

For a large system with n = 10,000, the factorization costs roughly 3.3 × 10¹¹ flops — expensive, but done once. Each subsequent solve then costs only 2 × 10⁸ flops, roughly 1,600 times cheaper. This asymmetry is why LU decomposition, not repeated Gaussian elimination, is the algorithm of choice in every serious scientific computing library.

Partial Pivoting for Numerical Stability

The basic LU factorization assumes that no pivot is zero — that we never need to divide by zero during elimination. In practice, even nonzero pivots can be numerically dangerous if they are very small: the resulting multipliers become very large, amplifying rounding errors in floating-point arithmetic.

Partial pivoting avoids this by reordering rows before each elimination step: at column k, search rows k through n for the entry with the largest absolute value, and swap that row to the pivot position. This ensures all multipliers satisfy |m_ik| ≤ 1, which bounds error growth.

Row reordering is encoded by a permutation matrix P. Instead of factoring A directly, partial pivoting computes the factorization of a permuted version of A:

Partial pivoting is nearly always sufficient. Complete pivoting — which also permutes columns — provides even stronger stability guarantees but requires tracking column swaps and is rarely used in practice.

A Worked 3×3 Example

Consider the matrix

A = [[2, 1, 1], [4, 3, 3], [8, 7, 9]]

Step 1 — Eliminate column 1. The pivot is a₁₁ = 2. The multipliers are m₂₁ = 4/2 = 2 and m₃₁ = 8/2 = 4. Subtracting 2 × row 1 from row 2 and 4 × row 1 from row 3 gives a new matrix with zeros in column 1 below the diagonal, and stores m₂₁ = 2 and m₃₁ = 4 in L.

Step 2 — Eliminate column 2. The new pivot is 1 (the (2,2) entry after step 1). The multiplier is m₃₂ = 1/1 = 1. Subtracting 1 × row 2 from row 3 completes U, and stores m₃₂ = 1 in L.

To solve Ax = b for, say, b = [1, 3, 9]ᵀ: forward substitution gives y₁ = 1, y₂ = 3 − 2·1 = 1, y₃ = 9 − 4·1 − 1·1 = 4. Back substitution gives x₃ = 4/2 = 2, x₂ = (1 − 1·2)/1 = −1, x₁ = (1 − 1·(−1) − 1·2)/2 = 0. So x = [0, −1, 2]ᵀ.