How to Perform Cholesky Decomposition in Python

Key Insights

Cholesky decomposition factors a positive-definite matrix into L·L^T, offering twice the speed of LU decomposition for symmetric matrices and requiring half the storage.
NumPy’s implementation is sufficient for most use cases, but SciPy provides additional control over triangular form and better numerical stability for ill-conditioned matrices.
The decomposition is essential for efficiently generating correlated random variables, solving linear systems with symmetric coefficient matrices, and computing matrix operations in statistical modeling.

Introduction to Cholesky Decomposition

Cholesky decomposition is a specialized matrix factorization technique that decomposes a positive-definite matrix A into the product of a lower triangular matrix L and its transpose: A = L·L^T. This decomposition is named after André-Louis Cholesky, a French military officer who developed it for geodetic calculations.

You’ll encounter Cholesky decomposition frequently in numerical linear algebra, particularly when solving systems of linear equations with symmetric, positive-definite coefficient matrices. It’s approximately twice as fast as LU decomposition for these matrices and requires only half the storage since you only need to compute and store the lower triangular matrix.

The method appears in multiple domains: solving normal equations in least squares problems, generating correlated random variables for Monte Carlo simulations, computing inverse and determinant of covariance matrices in multivariate statistics, and Kalman filtering in control theory.

The critical constraint is that your matrix must be positive-definite. A matrix is positive-definite if it’s symmetric and all its eigenvalues are positive, or equivalently, if x^T·A·x > 0 for all non-zero vectors x.

Mathematical Foundation

A positive-definite matrix has special properties that make Cholesky decomposition possible. These matrices are always symmetric, invertible, and have positive eigenvalues. Common examples include covariance matrices and Gram matrices.

The decomposition produces a unique lower triangular matrix L where each element is computed from previous elements in the same row and column. The diagonal elements are always positive for a positive-definite matrix.

Compared to other factorizations, Cholesky offers significant advantages for symmetric positive-definite matrices. LU decomposition works for general square matrices but requires more computation. QR decomposition is more general and numerically stable but slower for this specific case.

Here’s how to verify a matrix is positive-definite before attempting decomposition:

import numpy as np

def is_positive_definite(A):
    """Check if a matrix is positive-definite using eigenvalues."""
    # Check symmetry
    if not np.allclose(A, A.T):
        return False, "Matrix is not symmetric"
    
    # Compute eigenvalues
    eigenvalues = np.linalg.eigvals(A)
    
    # Check if all eigenvalues are positive
    if np.all(eigenvalues > 0):
        return True, "Matrix is positive-definite"
    else:
        return False, f"Matrix has non-positive eigenvalues: {eigenvalues}"

# Example: Create a positive-definite matrix
A = np.array([[4, 2, 1],
              [2, 5, 3],
              [1, 3, 6]])

is_pd, message = is_positive_definite(A)
print(f"{message}")
print(f"Eigenvalues: {np.linalg.eigvals(A)}")

Basic Implementation with NumPy

NumPy provides a straightforward implementation through numpy.linalg.cholesky(). This function returns the lower triangular matrix L by default.

Let’s walk through a complete example:

import numpy as np

# Define a positive-definite matrix
A = np.array([[4, 2, 1],
              [2, 5, 3],
              [1, 3, 6]], dtype=float)

print("Original matrix A:")
print(A)

# Perform Cholesky decomposition
L = np.linalg.cholesky(A)

print("\nLower triangular matrix L:")
print(L)

# Verify the decomposition by reconstructing A
A_reconstructed = L @ L.T

print("\nReconstructed matrix (L·L^T):")
print(A_reconstructed)

# Check if reconstruction matches original
print("\nVerification (should be True):")
print(np.allclose(A, A_reconstructed))

# Calculate the difference
print("\nMaximum absolute difference:")
print(np.max(np.abs(A - A_reconstructed)))

The reconstruction should match the original matrix within numerical precision. The np.allclose() function handles floating-point comparison with appropriate tolerance.

Using SciPy for Advanced Cases

SciPy’s implementation offers more control and better numerical stability for challenging cases. The key difference is the ability to specify whether you want the lower or upper triangular matrix.

import numpy as np
from scipy import linalg as scipy_linalg

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

# Lower triangular (default)
L_lower = scipy_linalg.cholesky(A, lower=True)
print("SciPy lower triangular:")
print(L_lower)

# Upper triangular
U_upper = scipy_linalg.cholesky(A, lower=False)
print("\nSciPy upper triangular:")
print(U_upper)

# Verify both decompositions
print("\nVerify L·L^T:")
print(np.allclose(A, L_lower @ L_lower.T))

print("Verify U^T·U:")
print(np.allclose(A, U_upper.T @ U_upper))

# Compare with NumPy
L_numpy = np.linalg.cholesky(A)
print("\nNumPy vs SciPy lower (should be identical):")
print(np.allclose(L_numpy, L_lower))

SciPy’s implementation uses LAPACK routines that handle numerical stability better for ill-conditioned matrices. For production code dealing with real-world data, prefer SciPy.

Practical Applications

The most powerful application is generating correlated random variables. When you need multivariate normal samples with a specific covariance structure, Cholesky decomposition provides an efficient solution.

import numpy as np
import matplotlib.pyplot as plt

def generate_correlated_samples(mean, cov, n_samples):
    """Generate multivariate normal samples using Cholesky decomposition."""
    # Perform Cholesky decomposition
    L = np.linalg.cholesky(cov)
    
    # Generate independent standard normal samples
    n_vars = len(mean)
    z = np.random.standard_normal((n_samples, n_vars))
    
    # Transform to correlated samples: X = μ + L·Z
    samples = mean + z @ L.T
    
    return samples

# Define parameters
mean = np.array([0, 0])
cov = np.array([[1.0, 0.8],
                [0.8, 1.0]])

# Generate samples
samples = generate_correlated_samples(mean, cov, 1000)

# Verify the empirical covariance matches
empirical_cov = np.cov(samples.T)
print("Target covariance:")
print(cov)
print("\nEmpirical covariance:")
print(empirical_cov)

# Solving linear systems efficiently
# For A·x = b where A is positive-definite
A = np.array([[4, 2, 1],
              [2, 5, 3],
              [1, 3, 6]], dtype=float)
b = np.array([1, 2, 3])

# Using Cholesky decomposition
L = np.linalg.cholesky(A)

# Solve L·y = b using forward substitution
y = np.linalg.solve(L, b)

# Solve L^T·x = y using backward substitution
x = np.linalg.solve(L.T, y)

print("\nSolution to A·x = b:")
print(x)
print("Verification A·x:")
print(A @ x)

Error Handling and Edge Cases

Robust code must handle cases where decomposition fails. The most common issue is attempting to decompose a non-positive-definite matrix.

import numpy as np
from scipy import linalg

def safe_cholesky(A, regularization=1e-6):
    """
    Perform Cholesky decomposition with error handling and regularization.
    """
    try:
        # Attempt standard decomposition
        L = np.linalg.cholesky(A)
        return L, "success"
    
    except np.linalg.LinAlgError:
        # Matrix is not positive-definite
        print(f"Standard decomposition failed. Attempting regularization...")
        
        # Add small value to diagonal (Tikhonov regularization)
        A_reg = A + regularization * np.eye(A.shape[0])
        
        try:
            L = np.linalg.cholesky(A_reg)
            return L, "regularized"
        
        except np.linalg.LinAlgError:
            # Even regularization failed
            return None, "failed"

# Example: Nearly singular matrix
A_singular = np.array([[1, 0.999],
                       [0.999, 1]])

L, status = safe_cholesky(A_singular)
print(f"Decomposition status: {status}")

# Example: Non-positive-definite matrix
A_non_pd = np.array([[1, 2],
                     [2, 1]])

L, status = safe_cholesky(A_non_pd)
print(f"Non-PD matrix status: {status}")

if L is not None:
    print("Decomposition succeeded (possibly with regularization)")

Performance Considerations

Cholesky decomposition has O(n³/3) complexity, making it twice as fast as LU decomposition’s O(2n³/3) for symmetric matrices. However, the actual performance depends on matrix size and structure.

import numpy as np
import time

def benchmark_cholesky(sizes):
    """Benchmark Cholesky decomposition for different matrix sizes."""
    results = []
    
    for n in sizes:
        # Create a random positive-definite matrix
        A = np.random.randn(n, n)
        A = A @ A.T + n * np.eye(n)  # Ensure positive-definite
        
        # Time the decomposition
        start = time.time()
        L = np.linalg.cholesky(A)
        elapsed = time.time() - start
        
        results.append((n, elapsed))
        print(f"Size {n}x{n}: {elapsed:.6f} seconds")
    
    return results

# Benchmark different sizes
sizes = [100, 200, 500, 1000, 2000]
results = benchmark_cholesky(sizes)

# Compare with solving linear systems
n = 1000
A = np.random.randn(n, n)
A = A @ A.T + n * np.eye(n)
b = np.random.randn(n)

# Method 1: Direct solve
start = time.time()
x1 = np.linalg.solve(A, b)
time_direct = time.time() - start

# Method 2: Using Cholesky
start = time.time()
L = np.linalg.cholesky(A)
y = np.linalg.solve(L, b)
x2 = np.linalg.solve(L.T, y)
time_cholesky = time.time() - start

print(f"\nDirect solve: {time_direct:.6f}s")
print(f"Cholesky solve: {time_cholesky:.6f}s")
print(f"Speedup: {time_direct/time_cholesky:.2f}x")

Use Cholesky decomposition when you have symmetric positive-definite matrices, especially when solving multiple systems with the same coefficient matrix. For general matrices, stick with LU decomposition. For overdetermined systems, use QR decomposition instead.

The key to effective use is recognizing when your problem structure matches Cholesky’s requirements. Covariance matrices, normal equations, and kernel matrices in machine learning are prime candidates. Combined with proper error handling and numerical stability checks, Cholesky decomposition becomes an indispensable tool in your numerical computing toolkit.