LU Decomposition

LU decomposition, also called LU factorization, is a direct numerical method in linear algebra and scientific computing that factors a square matrix A ∈ ℝ^n×n (or ℂ^n×n) into the product of a lower triangular matrix L and an upper triangular matrix U. Introduced in modern matrix notation by Alan Turing in 1948 to analyze rounding errors in automated computing machines, LU decomposition is the matrix-algebraic formulation of Gaussian elimination. Once a matrix A is factored as PA = LU (incorporating a permutation matrix P for row pivoting), any linear system Ax = b can be solved in O(n²) arithmetic operations by two triangular sweeps: forward substitution on Ly = Pb followed by back substitution on Ux = y.

Compared with QR decomposition, which requires (4/3)n³ floating-point operations (flops) for a square matrix, LU decomposition requires only (2/3)n³ flops—half the arithmetic work—making it the default general-purpose solver for dense non-symmetric systems of linear equations, matrix inversion, determinant evaluation, and the inner linear solves of Newton's method for nonlinear systems. When A is symmetric and positive-definite, LU decomposition specializes to Cholesky decomposition (A = LL^T), which further halves the arithmetic cost to (1/3)n³ flops.

In standard Doolittle convention, L ∈ ℝ^n×n is a unit lower triangular matrix (meaning l_ii = 1 along the main diagonal and l_ij = 0 for j > i), while U ∈ ℝ^n×n is an upper triangular matrix (u_ij = 0 for i > j). At step k of Gaussian elimination (for k = 1, …, n − 1), multiples of row k are subtracted from rows i = k + 1, …, n using the multipliers l_ik = a_ik^(k) / a_kk^(k) to zero out all entries below the pivot u_kk = a_kk^(k). Recording each elimination stage as a unit lower triangular Frobenius matrix M_k yields M_n−1⋯M₂M₁A = U. Because the inverse of a Frobenius elimination matrix simply negates its subdiagonal multipliers, the product L = M₁⁻¹M₂⁻¹⋯M_n−1⁻¹ assembles all multipliers l_ik directly into their corresponding subdiagonal positions without any additional arithmetic.

A square matrix A possesses a unique unpivoted LU decomposition with unit lower triangular L if and only if all of its leading principal submatrices A_{1:k, 1:k} for k = 1, …, n − 1 are non-singular (having non-zero leading principal minors). Even when A itself is invertible—for example, the 2 × 2 permutation matrix [0, 1; 1, 0]—a zero pivot u_kk = 0 causes unpivoted elimination to halt at that step unless rows are reordered.

Doolittle, Crout, and LDU Variants

Three closely related normalizations of the triangular factors appear in numerical software:

• Doolittle algorithm: Sets the diagonal of L to unity (l_ii = 1) and leaves the pivots on the diagonal of U. This is the standard convention in LAPACK (dgetrf), MATLAB, and SciPy.
• Crout algorithm: Sets the diagonal of U to unity (u_ii = 1) and places the pivots on the diagonal of L.
• LDU decomposition: Factors U = DŨ, yielding A = LDŨ, where both L and Ũ are unit triangular and D = diag(d₁, …, d_n) is diagonal. For symmetric matrices (A = A^T), symmetry forces Ũ = L^T, giving the symmetric indefinite LDL^T factorization.

Even when every pivot u_kk is non-zero, an unpivoted pivot of tiny magnitude (|u_kk| ≪ ‖A‖) produces huge multipliers l_ik = a_ik^(k) / u_kk. In finite-precision floating-point arithmetic, multiplying row k by a large multiplier l_ik and subtracting it from row i swamps the existing entries of row i with rounded multiples of row k, destroying accuracy.

Partial, Complete, and Rook Pivoting

To bound multiplier magnitude, Gaussian elimination with partial pivoting (GEPP) searches the active column k from row k to row n at each stage and swaps row k with the row p ≥ k that maximizes |a_pk^(k)|. This produces the factorization:

PA = LU,

where P ∈ {0, 1}^n×n is a row permutation matrix, and ensures that every off-diagonal multiplier satisfies |l_ik| ≤ 1. James H. Wilkinson proved that the computed factors L̂ and Û satisfy the backward-error relation L̂Û = PA + ΔA with ‖ΔA‖_∞ / ‖A‖_∞ ≤ n³ ρ_n ε_mach, where ρ_n = max_i,j,k |a_ij^(k)| / max_i,j |a_ij| is the pivot growth factor. Although Wilkinson constructed pathological matrices where partial pivoting suffers worst-case growth ρ_n = 2ⁿ⁻¹, in practical scientific and engineering workloads ρ_n almost never exceeds O(n^1/2), making partial pivoting backward stable in practice. When stricter guarantees are required, complete pivoting (PAQ = LU, searching the entire active trailing submatrix) bounds ρ_n ≈ O(n^{(1/2) ln n}) at the cost of O(n³) comparisons, while rook pivoting searches one row and one column alternatingly to achieve near-complete stability with only O(n²) comparison overhead.

On modern cache-hierarchical CPUs and GPUs, element-by-element elimination is bottlenecked by memory bandwidth rather than floating-point throughput. Block LU decomposition partitions A into panels of width b (typically 64 to 256 columns), factors the current n × b panel using partial pivoting, and updates the trailing (n − b) × (n − b) Schur complement A₂₂ ← A₂₂ − L₂₁U₁₂ via a single Level-3 BLAS dgemm matrix-matrix multiplication.

For large sparse matrices arising from finite-element or finite-volume discretizations of partial differential equations, direct LU factorization introduces fill-in—non-zero entries appearing in L and U at positions where A was zero. Sparse direct solvers (such as SuperLU, UMFPACK, and PARDISO) reorder rows and columns prior to elimination using minimum-degree or nested-dissection graph orderings to minimize fill-in. For 3D meshes where even reordered direct LU exceeds memory budgets, incomplete LU (ILU) factorization drops fill-in entries outside a prescribed sparsity pattern (ILU(0)) or below a numerical drop tolerance (ILUT), serving as an algebraic preconditioner that accelerates Krylov subspace solvers such as the generalized minimal residual method and the biconjugate gradient stabilized method.

How does LU decomposition speed up solving multiple right-hand sides?

Factoring PA = LU costs (2/3)n³ flops once. Each subsequent right-hand side vector b^(m) requires only forward substitution (Ly = Pb^(m), costing n² flops) and back substitution (Ux = y, costing n² flops), for a total of 2n² flops per solve—orders of magnitude faster than refactoring or computing A⁻¹ explicitly.

How is the determinant of a matrix computed from its LU factorization?

Because det(L) = 1 for a unit lower triangular matrix and det(P) = (−1)^s where s is the number of row swaps in P, the determinant of A is simply det(A) = (−1)^s ∏_i=1ⁿ u_ii, computed in O(n) multiplications after factoring.

1. Turing, A. M. (1948). Rounding-off errors in matrix processes. The Quarterly Journal of Mechanics and Applied Mathematics, 1(1), 287–308. https://doi.org/10.1093/qjmam/1.1.287
2. Wilkinson, J. H. (1961). Error analysis of direct methods of matrix inversion. Journal of the ACM, 8(3), 281–330. https://doi.org/10.1145/321075.321076
3. Trefethen, L. N., & Schreiber, R. S. (1990). Average-case stability of Gaussian elimination. SIAM Journal on Matrix Analysis and Applications, 11(3), 335–360. https://doi.org/10.1137/0611023
4. Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms (2nd ed.). Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898718027