Cholesky Decomposition: Efficient LDL^T Factorization for SPD Matrices

In many scientific and engineering applications, you encounter systems of linear equations where the coefficient matrix is symmetric and positive definite (SPD). Examples include finite element analysis, machine learning (e.g., Gaussian processes), and financial modeling. Solving Ax = b for such matrices can be done efficiently using Cholesky decomposition, which factors A into the product of a lower triangular matrix L and its transpose: A = LL^T. This decomposition is not only numerically stable but also requires only half the operations of LU decomposition. For large-scale problems, the LDL^T variant (where D is diagonal) avoids square roots, making it even more robust. In this tutorial, you'll learn the mathematical foundation, implement both LL^T and LDL^T in Python, and explore real-world debugging scenarios. By the end, you'll be able to apply Cholesky decomposition confidently in production code.

What is Cholesky Decomposition?

Cholesky decomposition is a matrix factorization technique applicable to symmetric positive definite (SPD) matrices. It decomposes a matrix A into the product of a lower triangular matrix L and its transpose: A = LL^T. This is analogous to taking the square root of a matrix. The decomposition is unique if L has positive diagonal entries. The LDL^T variant factors A as LDL^T, where L is unit lower triangular (diagonal of 1s) and D is diagonal. This avoids computing square roots, which can be beneficial for numerical stability. The algorithm proceeds column by column, computing the entries of L using the known entries of A and previously computed entries of L. The operation count is about n^3/3 floating point operations, half that of LU decomposition, making it efficient for large systems.

import numpy as np

def cholesky_llt(A):
    """Compute LL^T Cholesky decomposition."""
    n = A.shape[0]
    L = np.zeros_like(A)
    for j in range(n):
        # Compute diagonal entry
        L[j, j] = np.sqrt(A[j, j] - np.dot(L[j, :j], L[j, :j]))
        # Compute off-diagonal entries in column j
        for i in range(j+1, n):
            L[i, j] = (A[i, j] - np.dot(L[i, :j], L[j, :j])) / L[j, j]
    return L

# Example
A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]], dtype=float)
L = cholesky_llt(A)
print("L:\n", L)
print("LL^T:\n", L @ L.T)

LDL^T Decomposition: Avoiding Square Roots

The LDL^T decomposition is a variant that avoids square roots by factoring A as LDL^T, where L is unit lower triangular (diagonal entries are 1) and D is diagonal. This is more numerically stable when the matrix is ill-conditioned. The algorithm is similar to LL^T but computes D and L separately. The steps: for each column j, compute D[j] = A[j,j] - sum(L[j,k]^2 D[k] for k<j), then for i>j, L[i,j] = (A[i,j] - sum(L[i,k]L[j,k]*D[k] for k<j)) / D[j]. This decomposition is especially useful in applications like Kalman filters where covariance matrices may become nearly singular.

import numpy as np

def cholesky_ldlt(A):
    """Compute LDL^T decomposition."""
    n = A.shape[0]
    L = np.eye(n)
    D = np.zeros(n)
    for j in range(n):
        # Compute D[j]
        D[j] = A[j, j] - np.dot(L[j, :j]**2, D[:j])
        # Compute L[i,j] for i>j
        for i in range(j+1, n):
            L[i, j] = (A[i, j] - np.dot(L[i, :j] * L[j, :j], D[:j])) / D[j]
    return L, D

A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]], dtype=float)
L, D = cholesky_ldlt(A)
print("L:\n", L)
print("D:", D)
print("LDL^T:\n", L @ np.diag(D) @ L.T)

Solving Linear Systems with Cholesky

Once you have the Cholesky factor L (or LDL^T), solving Ax = b becomes straightforward. For LL^T: first solve Ly = b (forward substitution), then solve L^T x = y (back substitution). For LDL^T: solve Lz = b, then D^{-1} z = y, then L^T x = y. The total cost is O(n^2) after the O(n^3) factorization. This is much faster than computing the inverse. Example: solve the system from the previous sections with b = [1, 2, 3]^T.

import numpy as np

def solve_llt(L, b):
    # Forward substitution: Ly = b
    y = np.zeros_like(b)
    for i in range(len(b)):
        y[i] = (b[i] - np.dot(L[i, :i], y[:i])) / L[i, i]
    # Back substitution: L^T x = y
    x = np.zeros_like(b)
    for i in reversed(range(len(b))):
        x[i] = (y[i] - np.dot(L[i+1:, i], x[i+1:])) / L[i, i]
    return x

A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]], dtype=float)
b = np.array([1, 2, 3])
L = cholesky_llt(A)
x = solve_llt(L, b)
print("x:", x)
print("Check Ax:", A @ x)

Computing the Determinant and Inverse

Cholesky decomposition can also be used to compute the determinant and inverse of an SPD matrix. The determinant is the product of squares of diagonal entries of L (for LL^T) or product of D entries (for LDL^T). This is O(n) after factorization. The inverse can be computed by solving A * X = I using Cholesky, which is more efficient than direct inversion. However, in most applications, you should avoid computing the inverse explicitly.

import numpy as np

def determinant_from_cholesky(L):
    # det(A) = (prod L_ii)^2
    return np.prod(np.diag(L))**2

A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]], dtype=float)
L = cholesky_llt(A)
det = determinant_from_cholesky(L)
print("det(A):", det)
print("Check numpy:", np.linalg.det(A))

# Inverse via solving A * X = I
n = A.shape[0]
I = np.eye(n)
inv = np.zeros((n, n))
for i in range(n):
    inv[:, i] = solve_llt(L, I[:, i])
print("Inverse:\n", inv)
print("Check A * inv:\n", A @ inv)

Numerical Stability and Regularization

In floating-point arithmetic, matrices that are theoretically positive definite may become indefinite due to rounding errors. This is common in covariance matrices from sensor data. To handle this, you can add a small regularization term: A_reg = A + epsilon * I, where epsilon is a small positive number (e.g., 1e-8). This shifts eigenvalues by epsilon, ensuring positive definiteness. The choice of epsilon is a trade-off: too large distorts the solution, too small may not fix the problem. Another approach is to use the LDL^T decomposition, which is more robust. In extreme cases, you may need to fall back to a more stable but slower method like QR decomposition.

import numpy as np

# Example of a nearly singular matrix
A = np.array([[1, 1], [1, 1 + 1e-12]], dtype=float)
print("Eigenvalues:", np.linalg.eigvalsh(A))

try:
    L = np.linalg.cholesky(A)
    print("Cholesky succeeded")
except np.linalg.LinAlgError:
    print("Cholesky failed - matrix not positive definite")
    # Regularize
    epsilon = 1e-8
    A_reg = A + epsilon * np.eye(2)
    L = np.linalg.cholesky(A_reg)
    print("Regularized Cholesky succeeded")

Applications: Monte Carlo and Kalman Filters

Cholesky decomposition is essential in Monte Carlo simulations for generating correlated random variables. If you have a covariance matrix Sigma, you can generate samples from a multivariate normal distribution by computing L = chol(Sigma) and then x = L * z, where z is a vector of independent standard normals. In Kalman filters, the covariance update step involves Cholesky to maintain numerical stability. The LDL^T form is particularly useful because it avoids square roots and can be updated efficiently. These applications require robust and efficient Cholesky implementations.

import numpy as np

# Generate correlated random variables
mu = np.array([0, 0])
Sigma = np.array([[1, 0.8], [0.8, 1]])
L = np.linalg.cholesky(Sigma)

# Generate 1000 samples
z = np.random.randn(2, 1000)
x = mu[:, np.newaxis] + L @ z
print("Sample covariance:\n", np.cov(x))
print("Theoretical covariance:\n", Sigma)