Cholesky Decomposition Calculator

What Is Cholesky Decomposition?

The Cholesky decomposition factors any symmetric positive definite (SPD) matrix A into the product A = LL^T, where L is a unique lower triangular matrix with strictly positive diagonal entries. Named after French military officer André-Louis Cholesky (1875–1918), who developed the technique for geodetic surveying calculations, this factorization is the cornerstone of numerical linear algebra and appears throughout scientific computing, statistics, and machine learning. Because it exploits the symmetry of A, Cholesky decomposition requires roughly half the arithmetic of a general LU factorization and is numerically the most stable direct method for SPD systems.

Prerequisites: Symmetric Positive Definite Matrices

The Cholesky decomposition exists and is unique if and only if the input matrix satisfies two conditions simultaneously:

• Symmetry: Every off-diagonal entry satisfies a_ij = a_ji, so A = A^T. The calculator enforces this by accepting a single value for each symmetric pair (e.g., A[1,2] = A[2,1]).
• Positive definiteness: All eigenvalues are strictly positive, equivalently x^TAx > 0 for every non-zero real vector x. A practical test known as Sylvester’s criterion requires all leading principal minors to be positive: for a 2×2 matrix, a₁₁ > 0 and a₁₁a₂₂ − a₁₂² > 0; for a 3×3 matrix, additionally det(A) > 0.

If either condition fails, the algorithm encounters the square root of a non-positive number during computation, and no real Cholesky factor exists. Real-world covariance matrices from multivariate data and Hessian matrices at strict local minima are naturally SPD.

The Cholesky-Banachiewicz Formulas

Given an n×n SPD matrix A, the entries of the lower triangular factor L are computed column by column using two recurrences. The diagonal entries are:

L_ii = √(a_ii − ∑_k=1ⁱ⁻¹ L_ik²)

The sub-diagonal entries (j > i) are:

L_ji = (a_ji − ∑_k=1ⁱ⁻¹ L_jkL_ik) / L_ii

All entries strictly above the main diagonal of L are zero. According to Cornell CS4210 (Golub & Van Loan, Chapter 7) and UCLA EE133A Lecture 12 on Cholesky Factorization, this algorithm requires approximately n³/3 floating-point multiply-add operations — roughly half the cost of LU factorization — making it the standard solver whenever the coefficient matrix is SPD.

Step-by-Step: 2×2 Matrix

For A = [[a₁₁, a₁₂], [a₁₂, a₂₂]], compute three non-zero entries of L in order:

• L₁₁ = √a₁₁
• L₂₁ = a₁₂ / L₁₁
• L₂₂ = √(a₂₂ − L₂₁²)

Step-by-Step: 3×3 Matrix

After completing the 2×2 block, three additional sub-diagonal entries extend the factor:

• L₃₁ = a₁₃ / L₁₁
• L₃₂ = (a₂₃ − L₃₁ ⋅ L₂₁) / L₂₂
• L₃₃ = √(a₃₃ − L₃₁² − L₃₂²)

Worked Numerical Example

Let A = [[4, 2], [2, 3]]. First verify positive definiteness: a₁₁ = 4 > 0 and det(A) = (4)(3) − (2)² = 8 > 0. Apply the formulas:

• L₁₁ = √4 = 2.0000
• L₂₁ = 2 / 2 = 1.0000
• L₂₂ = √(3 − 1²) = √2 ≈ 1.4142

Verification: [[2, 0], [1, 1.4142]] × [[2, 1], [0, 1.4142]] = [[4, 2], [2, 3]] = A. The matrix determinant equals the square of the product of diagonal entries: (2 × 1.4142)² = 8, consistent with direct computation.

Numerical Stability

The Cholesky algorithm is backward stable for SPD matrices: the computed factor L satisfies (A + E) = LL^T where the perturbation E is small relative to machine precision. No row interchanges (pivoting) are needed, which simplifies implementation and avoids the fill-in overhead that partial pivoting introduces in LU. A closely related variant, the LDL^T decomposition, avoids square roots entirely by writing A = LDL^T with D diagonal, offering a modest speed advantage on architectures where square roots are expensive.

Key Applications

• Solving linear systems: Factor A = LL^T, forward-substitute to solve Ly = b, then back-substitute to solve L^Tx = y — two O(n²) triangular solves after an O(n³/3) factorization.
• Monte Carlo simulation: Generate correlated random vectors via x = Lz, where z is a vector of independent standard normal draws; the sample covariance of x matches the target matrix A.
• Kalman filtering: Maintain numerical positive definiteness of covariance matrices across hundreds of prediction-update cycles in navigation and robotics.
• Optimization: Compute Newton search directions efficiently when the Hessian at the current iterate is confirmed positive definite.
• Gaussian process regression: Evaluate multivariate normal log-likelihoods and draw posterior samples in Bayesian machine learning pipelines, where the kernel (covariance) matrix is SPD by construction.
• Portfolio risk: Decompose asset return covariance matrices to sample correlated return scenarios and compute portfolio variance in O(n) per draw after factorization.