Cholesky Decomposition

What is Cholesky Decomposition?

The Cholesky decomposition is a special factorization available for symmetric positive definite (SPD) matrices — the most important class of matrices in numerical linear algebra. For any SPD matrix A, Cholesky guarantees the existence of a unique lower triangular matrix L with positive diagonal entries such that A = LLᵀ. This factorization is essentially the "square root" of a matrix: just as any positive real number x has a square root √x, any SPD matrix has a triangular "square root" L.

The elegance of Cholesky decomposition lies in its efficiency. Compared to LU decomposition (which factors A = LU for general square matrices), Cholesky requires only half the arithmetic operations — roughly n³/6 flops versus n³/3 for LU. It also needs only half the storage since the symmetry of A is fully exploited. And it is unconditionally stable without pivoting: the algorithm never encounters a zero or near-zero pivot as long as A is truly positive definite.

Cholesky decomposition is the go-to method whenever you need to solve Ax = b repeatedly with the same SPD matrix A, compute the determinant of A, sample from a multivariate Gaussian distribution, or test whether a matrix is positive definite. It appears in Kalman filters, Gaussian process regression, finite element methods, and optimization algorithms.

Symmetric Positive Definite Matrices

A matrix A is symmetric positive definite if it is square, equal to its own transpose (Aᵀ = A), and satisfies xᵀAx > 0 for every non-zero vector x. Intuitively, a positive definite matrix always "pushes" vectors away from the origin — it has no direction that it flips, shrinks to zero, or reflects through the origin.

The most natural source of SPD matrices in practice is the normal equations matrix AᵀA (or AAᵀ), which arises in least squares problems, PCA, and linear regression. If A has full column rank, then AᵀA is always SPD. Other common SPD matrices include covariance matrices in statistics (they measure variance, so xᵀΣx is the variance of the linear combination xᵀy, which must be non-negative), stiffness matrices in finite element analysis, and the Gram matrix of a set of vectors.

Equivalent conditions for a symmetric matrix to be positive definite include: all eigenvalues are strictly positive, all leading principal minors have positive determinant (Sylvester's criterion), and — most practically — the Cholesky algorithm completes without encountering a non-positive diagonal entry.

The Cholesky Factorization: A = LLᵀ

The Cholesky factorization writes an n × n SPD matrix A as the product of a lower triangular matrix L and its transpose Lᵀ (which is upper triangular). The diagonal entries of L are all strictly positive. The factorization is unique: there is exactly one such L for any SPD matrix A.

An alternative form, sometimes called the LDLᵀ decomposition, factors A = LDLᵀ where L is unit lower triangular (diagonal entries all equal 1) and D is diagonal with positive entries. The LDLᵀ form avoids square roots entirely — useful when floating-point square roots are expensive or when you want to work with indefinite symmetric matrices (where D may have negative entries). The standard Cholesky A = LLᵀ is related by writing D = diag(d₁, …, dₙ) and absorbing √dᵢ into each column of L.

The Cholesky Algorithm

The Cholesky algorithm proceeds column by column, computing each column of L in turn. For column j (j = 1, …, n):

1. Diagonal entry: ℓⱼⱼ = √(aⱼⱼ − Σₖ₌₁ʲ⁻¹ ℓⱼₖ²). This is the square root of the "remaining" diagonal mass after subtracting the contributions from previously computed columns. If this argument is ≤ 0, A is not positive definite.
2. Below-diagonal entries: for i = j+1, …, n, compute ℓᵢⱼ = (aᵢⱼ − Σₖ₌₁ʲ⁻¹ ℓᵢₖℓⱼₖ) / ℓⱼⱼ. These are the remaining entries in column j, scaled by the diagonal.

The algorithm has cost O(n³/6) multiplications and n square root operations. For comparison, LU decomposition costs O(n³/3) multiplications. The Cholesky algorithm is therefore about twice as fast as LU and needs only the lower triangle of A in storage. There is no need for pivoting — the diagonal entries ℓⱼⱼ are guaranteed positive as long as A is positive definite, so no zero divisors arise.

Worked Example: 3 × 3 Cholesky

Consider the SPD matrix:

Solving Linear Systems with Cholesky

Once A = LLᵀ is computed, solving Ax = b becomes a two-step triangular system solve — the same structure as LU decomposition but with only one triangle to store:

1. Forward substitution: solve Ly = b for y. Since L is lower triangular, this costs O(n²) and proceeds from the first component to the last.
2. Back substitution: solve Lᵀx = y for x. Since Lᵀ is upper triangular, this costs O(n²) and proceeds from the last component back to the first.

If you need to solve Ax = bᵢ for many right-hand sides b₁, b₂, …, bₘ with the same matrix A, you factor A = LLᵀ once (O(n³/6)) and then do m pairs of triangular solves (each O(n²)). This is especially important in optimization algorithms that repeatedly solve a system with the same Hessian, or in Bayesian inference where you compute multiple posterior draws with the same covariance matrix.

Determinant via Cholesky

The determinant of A is easy to compute once L is known: det(A) = det(L) · det(Lᵀ) = det(L)² = (ℓ₁₁ · ℓ₂₂ · ⋯ · ℓₙₙ)². Since the diagonal entries of L are all positive, this product is guaranteed to be positive — consistent with the fact that all positive definite matrices have positive determinant. In many statistical applications (e.g., evaluating a multivariate Gaussian log-likelihood), you need log det(A) = 2 Σᵢ log ℓᵢᵢ, which is numerically stable to compute from the Cholesky factor.

Testing Positive Definiteness

The Cholesky algorithm provides the most practical test for positive definiteness: attempt the factorization. If it completes without encountering a non-positive diagonal argument (i.e., aⱼⱼ − Σₖ<ⱼ ℓⱼₖ² > 0 at every step), A is positive definite. If the argument is zero or negative at any step, A is not positive definite. This is more reliable than checking eigenvalues (which requires the full eigendecomposition) or Sylvester's criterion (which requires n determinant computations).

In practice, numerical noise can make a nearly positive semidefinite matrix fail the Cholesky test even though all its theoretical eigenvalues are positive. A common remedy is to add a small regularization ε to the diagonal: compute Cholesky of A + εI. If this succeeds, A is numerically positive definite. The smallest ε for which Cholesky succeeds gives an estimate of how far A is from the boundary of positive definiteness.

Applications of Cholesky Decomposition

Solving Normal Equations

In linear least squares, the normal equations AᵀAx = Aᵀb involve the SPD matrix AᵀA (when A has full column rank). Cholesky is the standard direct method for solving these equations. However, numerical analysts often prefer QR decomposition applied directly to A, since forming AᵀA squares the condition number: κ(AᵀA) = κ(A)², amplifying round-off errors for ill-conditioned A. For well-conditioned problems, Cholesky on the normal equations is efficient and accurate.

Generating Correlated Random Variables

To generate samples from the multivariate Gaussian distribution 𝒩(μ, Σ), you compute the Cholesky factor L of the covariance matrix Σ = LLᵀ, generate a vector z of independent standard normal samples, and compute x = μ + Lz. The result has mean μ and covariance E[(Lz)(Lz)ᵀ] = LE[zzᵀ]Lᵀ = LILᵀ = Σ. This is the standard algorithm in Monte Carlo simulation, Bayesian sampling, and generative models for correlated data.

Kalman Filtering

The Kalman filter maintains a state estimate x̂ and an error covariance matrix P. At each step, P must be updated and inverted to compute the Kalman gain K. Since P is always SPD (it is a covariance matrix), Cholesky is the natural tool for this computation. Square-root Kalman filter variants propagate the Cholesky factor of P directly, avoiding explicit computation of P itself. These square-root filters are numerically more stable, ensuring that P remains SPD despite floating-point errors that could otherwise make P non-symmetric or non-positive-definite.

Gaussian Process Regression

Gaussian process (GP) regression requires solving a linear system with the kernel matrix K (which is SPD by construction) and computing log det(K). Both operations are done via the Cholesky factorization K = LLᵀ: the system K⁻¹y is solved as two triangular solves, and log det(K) = 2 Σᵢ log Lᵢᵢ. The cubic cost O(n³) of Cholesky is the main computational bottleneck of exact GP inference, motivating sparse, low-rank, and iterative approximations for large datasets.

Finite Element Methods

Structural analysis, heat transfer, and electromagnetic simulations by finite element methods (FEM) produce large sparse SPD systems Ku = f, where K is the stiffness matrix. Cholesky (or its sparse variant) is the direct solver of choice. For sparse K, the Cholesky fill-in pattern (new non-zeros introduced by the factorization) depends on the ordering of the unknowns; reordering algorithms (e.g., minimum degree, nested dissection) minimize fill-in and reduce both memory and computation.

Numerical Properties and Stability

Cholesky decomposition is numerically stable without pivoting for positive definite matrices. The growth factor (ratio of largest element in L to largest element in A) is bounded by 1, so floating-point errors in the factorization are guaranteed small. This is in contrast to LU without pivoting, which can have arbitrarily large growth factors and must use partial pivoting for reliability.

The backward error of the Cholesky algorithm satisfies (A + ΔA) = LLᵀ where ‖ΔA‖ ≤ ε_mach · ‖A‖ (up to small constants). This means the computed L is the exact Cholesky factor of a matrix A + ΔA close to A in relative terms. Combined with the condition number κ(A) = σ_max/σ_min (or equivalently λ_max/λ_min for SPD matrices), the forward error in the solution x to Ax = b satisfies ‖Δx‖/‖x‖ ≈ κ(A) · ε_mach.

Comparison with Other Decompositions

How does Cholesky compare to the other factorizations in Module 7?

• vs. LU decomposition: Cholesky is twice as fast and needs no pivoting, but requires A to be SPD. LU works for any square invertible matrix but costs 2× more and needs partial pivoting for numerical stability.
• vs. QR decomposition: QR (applied directly to A) avoids squaring the condition number and is more stable for ill-conditioned problems. QR costs O(2n³/3) — more than Cholesky — but works for non-square and rank-deficient matrices.
• vs. SVD: SVD is the most general and most informative (reveals rank, condition number, all subspaces) but costs O(n³) with a large constant — roughly 6× more than Cholesky. SVD is the choice when you need maximum numerical insight; Cholesky is the choice when A is SPD and speed matters.
• vs. Eigendecomposition: For SPD matrices, A = QΛQᵀ (eigendecomposition) and A = LLᵀ (Cholesky) are both available. Eigendecomposition provides spectral information (individual eigenvalues) but costs more than Cholesky. For solving Ax = b, Cholesky is preferred; for PCA or diagonalization, eigendecomposition is needed.