Cholesky Decomposition Calculator

What is Cholesky Decomposition?

The Cholesky decomposition (or Cholesky factorization) is a method of decomposing a symmetric positive definite matrix into the product of a lower triangular matrix and its transpose. Named after the French mathematician Andre-Louis Cholesky, it is one of the most efficient and numerically stable ways to solve linear systems involving symmetric positive definite (SPD) matrices.

Given a symmetric positive definite matrix A, the Cholesky decomposition finds a unique lower triangular matrix L with strictly positive diagonal entries such that:

Here L is lower triangular (all entries above the main diagonal are zero) and Lᵀ is its transpose, which is upper triangular. The matrix L is sometimes called the Cholesky factor or the square root of the matrix A, because L plays a role analogous to the square root of a positive number.

Requirements: Symmetric Positive Definite

The Cholesky decomposition exists if and only if the matrix A is symmetric positive definite. A real matrix A is symmetric positive definite when:

• Symmetric: A = Aᵀ, meaning a_ij = a_ji for all i, j. The matrix equals its own transpose.
• Positive definite: xᵀAx > 0 for every nonzero vector x. Intuitively, the matrix "curves upward" in all directions — the quadratic form it defines has a unique minimum.

If A is only positive semidefinite (xᵀAx ≥ 0, allowing equality for nonzero x), the decomposition still exists but L may have zeros on the diagonal and is no longer unique.

How to Check if a Matrix is Positive Definite

Before applying the Cholesky decomposition, you need to verify that your matrix is SPD. There are several equivalent ways to test positive definiteness:

1. Eigenvalue test: All eigenvalues of A are strictly positive. This is the definition-level check but can be expensive to compute.
2. Leading principal minors (Sylvester's criterion): Every leading principal minor (the determinant of the top-left k x k submatrix for k = 1, 2, ..., n) is strictly positive.
3. Cholesky test: Attempt the Cholesky decomposition itself. If the algorithm completes without encountering a zero or negative value under a square root, the matrix is positive definite. This is often the most practical check.
4. Diagonal dominance (sufficient condition): If A is symmetric and strictly diagonally dominant (each diagonal entry exceeds the sum of absolute values of other entries in its row), then A is positive definite.

The Cholesky Algorithm Step by Step

The entries of the lower triangular matrix L are computed column by column, from left to right. For each column j = 1, 2, ..., n:

1. Diagonal entry: Compute the diagonal element L_jj using the formula:
   L_jj = sqrt( A_jj − ∑_{k=1}^{j-1} L_jk^2 )

   The value under the square root must be positive. If it is zero or negative, the matrix is not positive definite.

2. Below-diagonal entries: For each row i > j, compute:
   L_ij = ( A_ij − ∑_{k=1}^{j-1} L_ik * L_jk ) / L_jj

This process is repeated for every column until the entire lower triangle of L is filled. The entries above the diagonal are all zero by construction.

Worked Example: 3x3 Matrix

Consider the symmetric positive definite matrix:

Column 1:

• L_11 = sqrt(4) = 2
• L_21 = 2 / 2 = 1
• L_31 = 2 / 2 = 1

Column 2:

• L_22 = sqrt(5 − 1^2) = sqrt(4) = 2
• L_32 = (3 − 1*1) / 2 = 1

Column 3:

• L_33 = sqrt(6 − 1^2 − 1^2) = sqrt(4) = 2

You can verify: LLᵀ reconstructs the original matrix A.

Computational Complexity

The Cholesky decomposition requires approximately n³/3 floating-point operations for an n x n matrix. This is roughly half the cost of LU decomposition (which needs about 2n³/3 flops) and is one of the primary reasons Cholesky is preferred whenever the matrix is known to be SPD.

In addition to being faster, the Cholesky decomposition is also more numerically stable than LU for SPD matrices because it avoids the need for pivoting. The algorithm is guaranteed not to encounter growth in the entries of L, making it inherently backward-stable.

When to Use Cholesky vs LU

Both Cholesky and LU decomposition can factor square matrices and solve linear systems. However, they serve different purposes:

• Use Cholesky when the matrix is symmetric positive definite. It is twice as fast, requires half the storage, needs no pivoting, and is more numerically stable.
• Use LU when the matrix is a general (non-symmetric or indefinite) square matrix. LU with partial pivoting handles any invertible matrix.
• Cholesky as a diagnostic: If you suspect a matrix should be SPD but are not certain, attempting a Cholesky decomposition is a fast way to test. If the algorithm fails, the matrix is not positive definite.

In practice, covariance matrices, kernel matrices (Gram matrices), normal equation matrices (AᵀA), and stiffness matrices in finite element analysis are all SPD, making Cholesky the preferred factorization in these contexts.

Connection to the Square Root of a Matrix

The Cholesky factor L is sometimes called the "matrix square root" because it satisfies A = LLᵀ, analogous to a = (sqrt(a))^2 for positive scalars. However, the Cholesky factor is not the unique symmetric positive definite square root of A (which satisfies S^2 = A with S symmetric positive definite). Instead, it is a triangular square root, which is computationally much cheaper to obtain.

The unique symmetric square root can be found via the eigendecomposition A = QDQᵀ as S = QD^{1/2}Qᵀ, but this is significantly more expensive. The Cholesky factor is preferred in practice because its triangular structure makes forward and backward substitution trivial.

Applications of Cholesky Decomposition

The Cholesky decomposition appears across many areas of applied mathematics, statistics, and engineering:

• Monte Carlo simulation: To generate correlated random samples from a multivariate normal distribution N(μ, Σ), compute L = chol(Σ) and transform independent standard normal samples z as x = μ + Lz. This is the most common method for sampling correlated random variables.
• Kalman filters: The square-root Kalman filter uses the Cholesky factor of the covariance matrix to improve numerical stability and guarantee the covariance remains positive definite across filter updates.
• Optimization and quadratic programming: In Newton's method and interior-point methods, the Hessian matrix of a convex objective is SPD near a minimizer. The Cholesky decomposition efficiently solves the Newton step H*d = −g.
• Sampling from multivariate normal distributions: Beyond Monte Carlo, Bayesian inference engines (MCMC samplers, variational inference) heavily rely on Cholesky for parameterizing and sampling from Gaussian posteriors.
• Solving linear systems: When Ax = b with A SPD, factor A = LLᵀ and solve Ly = b (forward substitution) then Lᵀx = y (backward substitution). This is the fastest direct solver for SPD systems.
• Computing determinants: det(A) = det(L)^2 = (∏ L_ii)^2. Since L is triangular, its determinant is simply the product of its diagonal entries. Taking the log avoids overflow: log det(A) = 2 ∑ log(L_ii).
• Preconditioning: Incomplete Cholesky factorizations serve as preconditioners for iterative solvers (conjugate gradient method) applied to large sparse SPD systems.
• Geodesy and surveying: Cholesky himself developed the factorization for solving least-squares problems in geodetic surveying. The method was published posthumously after his death in World War I.

Variants of the Cholesky Decomposition

• LDLᵀ decomposition: An alternative form A = LDLᵀ where L is unit lower triangular (ones on the diagonal) and D is diagonal. This avoids computing square roots and is useful when working in exact arithmetic or over fields where square roots are expensive.
• Incomplete Cholesky: A sparse approximation where fill-in entries are dropped during the factorization. Commonly used as a preconditioner for the conjugate gradient method.
• Pivoted Cholesky: A variant with row and column permutations (A = PLLᵀPᵀ) that can handle positive semi-definite matrices and produces a low-rank approximation when stopped early.