Eric Campbell 6 min read Statistics

Linear Algebra: LU Decomposition Explained

LU decomposition is a fundamental matrix factorization technique that breaks down a square matrix A into the product of two triangular matrices: a lower triangular matrix L and an upper triangular...

Key Insights

  • LU decomposition factors a matrix into lower and upper triangular matrices, enabling efficient solution of linear systems—especially when solving multiple equations with the same coefficient matrix but different right-hand sides.
  • The method is essentially Gaussian elimination with bookkeeping, achieving O(n³) complexity for decomposition but only O(n²) for each subsequent solve, making it faster than direct methods for repeated solves.
  • Partial pivoting is critical for numerical stability in real-world implementations; never use naive LU decomposition in production without considering pivot selection.

Introduction to LU Decomposition

LU decomposition is a fundamental matrix factorization technique that breaks down a square matrix A into the product of two triangular matrices: a lower triangular matrix L and an upper triangular matrix U, such that A = LU. This seemingly simple transformation is one of the workhorses of numerical linear algebra.

Why should you care? Three reasons: solving systems of linear equations efficiently, computing determinants quickly, and inverting matrices. In statistics and machine learning, you encounter these operations constantly—from fitting linear regression models to computing covariance matrices to implementing optimization algorithms.

Here’s what LU decomposition looks like for a simple 3×3 matrix:

import numpy as np

# Original matrix A
A = np.array([[2, -1, 0],
              [-1, 2, -1],
              [0, -1, 2]], dtype=float)

# After LU decomposition: A = L @ U
L = np.array([[1, 0, 0],
              [-0.5, 1, 0],
              [0, -0.667, 1]], dtype=float)

U = np.array([[2, -1, 0],
              [0, 1.5, -1],
              [0, 0, 1.333]], dtype=float)

# Verify: L @ U equals A
print("A =\n", A)
print("\nL @ U =\n", np.round(L @ U, 3))
print("\nMatch:", np.allclose(A, L @ U))

The real power emerges when you need to solve Ax = b for multiple different b vectors. Decompose once, solve many times.

The Mathematics Behind LU Decomposition

LU decomposition is Gaussian elimination in disguise. When you perform row operations to convert a matrix to row echelon form, you’re implicitly computing the U matrix. The L matrix records the multipliers you used along the way.

The process works like this: start with your matrix A. To eliminate elements below the first pivot, you multiply row 1 by some factor and subtract from rows below. Store those factors in L. Continue column by column until you have an upper triangular U.

Here’s the critical requirement: the matrix must be non-singular (invertible), and all leading principal minors must be non-zero. If you encounter a zero pivot, you need partial pivoting—swapping rows to bring a non-zero element into the pivot position.

Let’s implement this step-by-step to see the mechanics:

def lu_decomposition_verbose(A):
    """LU decomposition with detailed output."""
    n = len(A)
    L = np.eye(n)
    U = A.copy()
    
    print("Starting matrix U:")
    print(U)
    print()
    
    for k in range(n-1):
        print(f"Step {k+1}: Eliminate column {k}")
        for i in range(k+1, n):
            # Calculate multiplier
            factor = U[i, k] / U[k, k]
            L[i, k] = factor
            
            # Eliminate
            U[i, k:] = U[i, k:] - factor * U[k, k:]
            
            print(f"  Row {i}: subtract {factor:.3f} × Row {k}")
        print("U =")
        print(U)
        print()
    
    return L, U

A = np.array([[4, 3, 2],
              [2, 1, 3],
              [6, 3, 4]], dtype=float)

L, U = lu_decomposition_verbose(A)
print("Final L:")
print(L)
print("\nFinal U:")
print(U)

Partial pivoting adds a permutation matrix P, giving you PA = LU instead of A = LU. This ensures numerical stability by avoiding division by small numbers.

Implementing LU Decomposition

Doolittle’s method is the standard algorithm. It computes L with ones on the diagonal and fills in U and L simultaneously through a clever indexing scheme.

def lu_decompose(A):
    """
    Doolittle's LU decomposition without pivoting.
    Returns L and U such that A = L @ U.
    """
    n = len(A)
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    
    for i in range(n):
        # Upper triangular matrix U
        for k in range(i, n):
            sum_val = sum(L[i][j] * U[j][k] for j in range(i))
            U[i][k] = A[i][k] - sum_val
        
        # Lower triangular matrix L
        for k in range(i, n):
            if i == k:
                L[i][i] = 1  # Diagonal as 1
            else:
                sum_val = sum(L[k][j] * U[j][i] for j in range(i))
                L[k][i] = (A[k][i] - sum_val) / U[i][i]
    
    return L, U

# Test implementation
A = np.array([[2, -1, 0],
              [-1, 2, -1],
              [0, -1, 2]], dtype=float)

L, U = lu_decompose(A)
print("Custom implementation:")
print("L =\n", L)
print("\nU =\n", U)
print("\nVerification:", np.allclose(A, L @ U))

# Compare with SciPy
from scipy.linalg import lu
P, L_scipy, U_scipy = lu(A)
print("\nSciPy (with pivoting):")
print("L =\n", L_scipy)
print("\nU =\n", U_scipy)

Time complexity is O(n³) for the decomposition—same as Gaussian elimination. But here’s the win: once decomposed, solving for each new right-hand side costs only O(n²).

Compared to QR decomposition (more stable, better for least squares) and Cholesky (faster but requires positive definite matrices), LU offers a good balance of generality and efficiency for square systems.

Solving Linear Systems with LU Decomposition

Once you have A = LU, solving Ax = b becomes a two-step process:

  1. Forward substitution: Solve Ly = b for y
  2. Backward substitution: Solve Ux = y for x

Both operations are O(n²), making repeated solves extremely efficient.

def forward_substitution(L, b):
    """Solve Ly = b for y."""
    n = len(b)
    y = np.zeros(n)
    
    for i in range(n):
        y[i] = b[i] - sum(L[i][j] * y[j] for j in range(i))
        y[i] /= L[i][i]
    
    return y

def backward_substitution(U, y):
    """Solve Ux = y for x."""
    n = len(y)
    x = np.zeros(n)
    
    for i in range(n-1, -1, -1):
        x[i] = y[i] - sum(U[i][j] * x[j] for j in range(i+1, n))
        x[i] /= U[i][i]
    
    return x

def solve_with_lu(A, b):
    """Solve Ax = b using LU decomposition."""
    L, U = lu_decompose(A)
    y = forward_substitution(L, b)
    x = backward_substitution(U, y)
    return x

# Example: solve a system
A = np.array([[3, 2, 1],
              [2, 3, 2],
              [1, 2, 3]], dtype=float)
b = np.array([10, 14, 14], dtype=float)

x = solve_with_lu(A, b)
print("Solution x:", x)
print("Verification Ax:", A @ x)
print("Expected b:", b)

# Benchmark: multiple solves
import time

n = 500
A_large = np.random.randn(n, n)
b_vectors = [np.random.randn(n) for _ in range(10)]

# Method 1: Direct solve each time
start = time.time()
for b in b_vectors:
    x = np.linalg.solve(A_large, b)
direct_time = time.time() - start

# Method 2: LU once, solve many
start = time.time()
from scipy.linalg import lu_factor, lu_solve
lu, piv = lu_factor(A_large)
for b in b_vectors:
    x = lu_solve((lu, piv), b)
lu_time = time.time() - start

print(f"\nDirect solve: {direct_time:.4f}s")
print(f"LU decomposition: {lu_time:.4f}s")
print(f"Speedup: {direct_time/lu_time:.2f}x")

The speedup for multiple solves is substantial—often 2-3x for even modest numbers of right-hand sides.

Practical Applications in Statistics

Linear regression is the canonical example. The normal equations (X^T X)β = X^T y require solving a linear system. With LU decomposition:

def linear_regression_lu(X, y):
    """
    Compute regression coefficients using LU decomposition.
    Solves (X^T X)β = X^T y
    """
    XTX = X.T @ X
    XTy = X.T @ y
    
    # LU decomposition approach
    from scipy.linalg import lu_factor, lu_solve
    lu, piv = lu_factor(XTX)
    beta = lu_solve((lu, piv), XTy)
    
    return beta

# Generate sample data
np.random.seed(42)
n_samples, n_features = 1000, 5
X = np.random.randn(n_samples, n_features)
true_beta = np.array([1.5, -2.0, 0.5, 3.0, -1.0])
y = X @ true_beta + np.random.randn(n_samples) * 0.5

# Compute coefficients
beta_lu = linear_regression_lu(X, y)
beta_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]

print("LU method:", beta_lu)
print("Least squares:", beta_lstsq)
print("True values:", true_beta)
print("Close match:", np.allclose(beta_lu, beta_lstsq))

Other applications include computing inverse covariance matrices for Mahalanobis distance, solving maximum likelihood estimation problems, and implementing Kalman filters.

Performance Considerations and Best Practices

Use LU decomposition when:

  • You have a square, non-singular system
  • You need to solve for multiple right-hand sides
  • The matrix isn’t particularly special (not positive definite, not sparse)

Avoid it when:

  • You have an overdetermined system (use QR instead)
  • The matrix is positive definite (use Cholesky—it’s 2x faster)
  • The matrix is sparse (use specialized sparse solvers)

Always use partial pivoting in production. The naive algorithm fails catastrophically with poor pivot choices.

# Performance comparison across matrix sizes
sizes = [100, 200, 500, 1000]
results = []

for n in sizes:
    A = np.random.randn(n, n)
    b = np.random.randn(n)
    
    # Direct solve
    start = time.time()
    x1 = np.linalg.solve(A, b)
    direct = time.time() - start
    
    # LU with SciPy
    start = time.time()
    lu, piv = lu_factor(A)
    x2 = lu_solve((lu, piv), b)
    lu_time = time.time() - start
    
    results.append((n, direct, lu_time))
    print(f"n={n:4d}: Direct={direct:.4f}s, LU={lu_time:.4f}s")

For production work, use LAPACK-backed libraries (NumPy, SciPy, Julia’s LinearAlgebra). They implement blocked algorithms and leverage BLAS for optimal cache performance.

Conclusion

LU decomposition is essential knowledge for anyone working with numerical computing. It transforms the repeated solution of linear systems from an O(n³) problem into an O(n²) one after initial decomposition. This matters enormously in iterative algorithms, optimization routines, and statistical computations.

The key takeaway: decompose once, solve many times. When you see code that repeatedly calls np.linalg.solve with the same coefficient matrix, reach for LU decomposition.

Start with library implementations—SciPy’s lu_factor and lu_solve are production-ready. Understand the algorithm so you can debug numerical issues and choose the right tool for each problem. And always remember: partial pivoting isn’t optional; it’s mandatory for numerical stability.

For deeper understanding, study the LAPACK documentation, explore blocked algorithms for large matrices, and investigate how sparse solvers exploit structure. The rabbit hole goes deep, but the fundamentals covered here will serve you well in 95% of practical applications.

Liked this? There's more.

Every week: one practical technique, explained simply, with code you can use immediately.