Cholesky Decomposition

Cholesky decomposition (or Cholesky factorization) is a direct matrix factorization algorithm in computation and algorithms and numerical linear algebra that decomposes a real symmetric positive-definite matrix A ∈ R^n×n into the product of a unique lower-triangular matrix L with positive diagonal entries and its transpose L^T (or, for a complex Hermitian positive-definite matrix, its conjugate transpose L^*). Devised around 1910 by French military geodesist and artillery officer André-Louis Cholesky to solve normal equations arising from the triangulation survey of France and Crete—and published posthumously by fellow officer Commandant Benoît in 1924—the factorization requires n³/3 floating-point operations (FLOPs), half the arithmetic cost and storage of general Gaussian elimination (LU decomposition). Because positive-definiteness bounds the growth of intermediate entries, Cholesky factorization is backward stable in IEEE floating-point arithmetic without row or column pivoting, serving as the standard direct solver for normal equations in Linear Regression, covariance inversions in Kalman Filter updates, Gaussian process log-marginal likelihoods, and correlated Monte Carlo simulation.

Symmetric Positive-Definite Matrices and Uniqueness

A real n × n matrix A is symmetric positive-definite (SPD) if A^T = A and the quadratic form satisfies x^TAx > 0 for all non-zero x ∈ Rⁿ \ {0}. Equivalently, all n eigenvalues of A are real and positive (λ_i > 0), and all n leading principal minors det(A_{1:k, 1:k}) for k = 1, ..., n are positive (Sylvester's criterion).

Theorem (Cholesky Factorization): Every real symmetric positive-definite matrix A ∈ R^n×n admits a unique factorization:

A = L L^T (or in upper-triangular form A = R^T R, where R = L^T)

such that L = [L_ij] is lower triangular (L_ij = 0 for j > i) with positive diagonal entries (L_ii > 0). The determinant and log-determinant follow directly from the diagonal of L:

det(A) = ( ∏_i=1ⁿ L_ii )², and log det(A) = 2 ∑_i=1ⁿ log(L_ii)

Evaluating log det(A) via the sum of logarithms avoids floating-point overflow or underflow on high-dimensional covariance matrices.

Square-Root-Free LDL^T Variant

Factoring out the diagonal matrix D₀ = diag(L₁₁, ..., L_nn) so L = L̃D₀ (where L̃ is unit lower triangular with L̃_ii = 1) and setting D = D₀² yields the LDLT decomposition A = L̃DL̃^T, which avoids all n square-root evaluations and extends to symmetric indefinite systems.

Cholesky–Banachiewicz and Cholesky–Crout Recurrences

Equating entries A_ij = ∑_k=1^min(i,j) L_ik L_jk for i ≥ j yields explicit recurrences for the diagonal (i = j) and subdiagonal (i > j) entries:

• Diagonal Entry (j = i): L_jj = √( A_jj − ∑_k=1^j−1 L_jk² )
• Off-Diagonal Entry (i > j): L_ij = ( 1 / L_jj ) ( A_ij − ∑_k=1^j−1 L_ik L_jk )

Because A is symmetric, only the lower triangle (n(n + 1)/2 floats) is stored, and L can overwrite A in-place via row-by-row (Cholesky–Banachiewicz), column-by-column (Cholesky–Crout), or blocked Level-3 BLAS (LAPACK DPOTRF) traversal.

Operation Count (n³/3 FLOPs) and Numerical Stability

Summing the inner-product operations across columns j = 1, ..., n gives:

Total FLOPs = ∑_j=1ⁿ 2(j − 1)(n − j + 1) + O(n²) = (1/3) n³ + O(n²)

compared to (2/3)n³ FLOPs for LU decomposition and (4/3)n³ FLOPs for Householder QR. Moreover, because A_ii = ∑_k=1ⁱ L_ik², every entry of L is bounded by |L_ik| ≤ √A_ii. James H. Wilkinson proved that this bound makes Cholesky decomposition backward stable without pivoting: the computed factor L̂ satisfies L̂L̂^T = A + ΔA with ‖ΔA‖₂/‖A‖₂ = O(n ε_mach).

Forward and Backward Substitution

Once A = LL^T is factored, solving Ax = b reduces to two n²-FLOP triangular sweeps:

1. Forward Substitution (Ly = b): y_i = (b_i − ∑_k=1ⁱ⁻¹ L_ik y_k) / L_ii for i = 1, ..., n.
2. Backward Substitution (L^Tx = y): x_i = (y_i − ∑_k=i+1ⁿ L_ki x_k) / L_ii for i = n, ..., 1.

For large sparse SPD matrices, Incomplete Cholesky factorization (IC(0)) discards fill-in outside the sparsity pattern of A to precondition the Conjugate Gradient Method.

Correlated Monte Carlo Sampling and Gaussian Processes

In computational statistics, Cholesky factorization enables two fundamental operations on multivariate normal distributions N(μ, Σ):

1. Correlated Normal Sampling: Given z ∼ N(0, I_n) and Σ = LL^T, the transformation x = μ + Lz has covariance Cov(x) = L E[zz^T] L^T = Σ.
2. Gaussian Process (GP) Marginal Likelihood: Evaluating log p(y | X) = −(1/2)y^TK⁻¹y − (1/2) log det(K) − (n/2) log(2π) factors K = LL^T, solves Lα̃ = y so y^TK⁻¹y = ‖α̃‖₂², and sums ∑_i=1ⁿ log L_ii.

Why is Cholesky decomposition twice as fast as LU decomposition?

When A is symmetric (A^T = A), the upper factor is the transpose of the lower factor (U = L^T), halving the number of triangular entries to compute and cutting the arithmetic count from (2/3)n³ to (1/3)n³.

What happens if Cholesky decomposition is applied to a nearly singular covariance matrix?

Roundoff error on zero or near-zero eigenvalues can produce a negative radicand A_jj − ∑ L_jk² < 0; numerical libraries add a diagonal Tikhonov jitter A ← A + εI (ε ≈ 10⁻⁸) or apply pivoted rank-revealing Cholesky factorization.

Why should one solve L(L^Tx) = b rather than computing A⁻¹ explicitly?

Forming A⁻¹ explicitly costs an extra (2/3)n³ FLOPs and amplifies forward roundoff error by roughly κ(A) compared to solving L(L^Tx) = b via forward and backward substitution in 2n² FLOPs.

1. Benoît, C. (1924). Note sur une méthode de résolution des équations normales provenant de l'application de la méthode des moindres carrés à un système d'équations linéaires en nombre inférieur à celui des inconnues (Procédé du Commandant Cholesky). Bulletin Géodésique, 2(1), 67–77. https://doi.org/10.1007/BF03031308
2. Wilkinson, J. H. (1968). Apriori error analysis of algebraic processes. Proceedings of the International Congress of Mathematicians, 629–640. https://doi.org/10.1137/1.9781611971446
3. Higham, N. J. (1990). Analysis of the Cholesky decomposition of a semi-definite matrix. In M. G. Cox & S. Hammarling (Eds.), Reliable Numerical Computation (pp. 161–185). Oxford University Press. https://doi.org/10.1093/oso/9780198535645.003.0010
4. Davis, T. A. (2006). Direct Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics (SIAM). https://doi.org/10.1137/1.9780898718881
5. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press. https://doi.org/10.56021/9781421407944