LU Decomposition: Factor Matrices, Solve Systems Fast

LU decomposition splits a square matrix A into a lower triangular L (zeros above the diagonal) times an upper triangular U (zeros below), so A = LU. The elimination steps of Gaussian elimination are exactly what build L and U. Once factored, solving Ax = b for any new b is cheap: solve Ly = b by forward substitution, then Ux = y by back substitution.

Real problems rarely have one right-hand side: a bridge under ten load patterns, or a circuit simulated at many time steps, would repeat the expensive elimination every time. LU does the heavy n³/3 work once, then solves each new case in only n² steps — a massive saving.

Key Concepts

• L (unit lower triangular): 1s on the diagonal, elimination multipliers below it.
• U (upper triangular): the row-echelon form produced by elimination.
• Doolittle's method: the standard construction — unit diagonal on L, pivots on U's diagonal.
• Pivoting: with row swaps it becomes PA = LU, where P records the swaps — the form used in practice for stability.

Formulas

• Factorization: A = LU (or PA = LU with pivoting), computed in about n³/3 operations — the same cost as one Gaussian elimination.
• Forward substitution: solve Ly = b from top to bottom, y_i = b_i − Σ_j<i l_ij y_j.
• Back substitution: solve Ux = y from bottom to top, x_i = (y_i − Σ_j>i u_ij x_j) / u_ii.
• Each extra solve after factorization costs only about n² operations, versus n³/3 for a fresh elimination.

Worked Example

Factor A = [[2, 1, 1], [4, 3, 3], [8, 7, 9]] by Doolittle. Multipliers: l₂₁ = 4/2 = 2, l₃₁ = 8/2 = 4. Updated block: u₂₂ = 3 − 2·1 = 1, u₂₃ = 3 − 2·1 = 1; l₃₂ = (7 − 4·1)/1 = 3; u₃₃ = 9 − 4·1 − 3·1 = 2. So L = [[1, 0, 0], [2, 1, 0], [4, 3, 1]], U = [[2, 1, 1], [0, 1, 1], [0, 0, 2]]. Solve Ax = b, b = (1, 2, 5). Forward: y₁ = 1; 2 + y₂ = 2 → y₂ = 0; 4 + y₃ = 5 → y₃ = 1. Back: 2x₃ = 1 → x₃ = 0.5; x₂ = −0.5; 2x₁ = 1 → x₁ = 0.5. Check: 2(0.5) − 0.5 + 0.5 = 1; 4(0.5) + 3(−0.5) + 3(0.5) = 2; 8(0.5) + 7(−0.5) + 9(0.5) = 5. Correct.

Common Misconceptions

• "LU decomposition always exists." Without pivoting it needs all leading pivots to be nonzero. Some matrices need row swaps (PA = LU) or cannot be factored at all if singular.
• "The factors L and U are unique." Only with a convention like Doolittle's unit diagonal on L. Without fixing one diagonal, infinitely many L–U pairs give the same A.
• "LU is faster than Gaussian elimination." For a single system the cost is identical. LU wins when the same A meets many different right-hand sides.

Real-World Applications

• Structural analysis: one stiffness matrix, many load cases — factorize once, solve for each load.
• Time-stepping simulations: implicit methods solve systems with the same matrix at every step.
• Matrix inversion: solve AX = I column by column using one factorization.
• Determinants: det(A) is the product of U's diagonal (times the sign of the row permutation).

LU decomposition factors A = LU once at the cost of one elimination, then solves each new right-hand side with two cheap triangular substitutions — the standard way to handle repeated linear systems and the backbone of professional numerical libraries.