Cholesky Decomposition

The Cholesky decomposition factors a symmetric positive definite (SPD) matrix $A$ as $A = LL^\top$, where $L$ is lower triangular with positive diagonal entries. It is essentially the 'square root' of an SPD matrix — analogous to writing a positive number as $x = (\sqrt{x})^2$.

A matrix is positive definite if $\mathbf{x}^\top A \mathbf{x} > 0$ for every non-zero $\mathbf{x}$. Equivalently, all eigenvalues are strictly positive. SPD matrices are ubiquitous: covariance matrices, Gram matrices $X^\top X$ (when $X$ has full column rank), and Hessians at strict minima are all SPD. Cholesky is the fastest, most numerically stable factorization available for this class.

The algorithm is specialized Gaussian elimination that exploits symmetry. For each row $i$, compute the diagonal entry $\ell_{ii} = \sqrt{a_{ii} - \sum_{k<i}\ell_{ik}^2}$ and the sub-diagonal entries via similar formulas. The process fails iff some diagonal element becomes non-positive — a clean test for positive definiteness: a matrix is SPD iff Cholesky succeeds without complex square roots.

Cholesky solves $A\mathbf{x} = \mathbf{b}$ in two triangular steps: first $L\mathbf{y} = \mathbf{b}$ (forward substitution), then $L^\top \mathbf{x} = \mathbf{y}$ (backward substitution). Each step costs $O(n^2)$; the decomposition itself is $O(n^3/3)$, half the cost of LU. That is why Cholesky is the default when solving normal equations for linear regression (since $X^\top X$ is SPD).

Beyond solving, Cholesky has distinctive ML uses. In Bayesian methods and Gaussian processes, it lets us sample from multivariate Gaussians: if $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, I)$ and $A = LL^\top$, then $L\mathbf{z} \sim \mathcal{N}(\mathbf{0}, A)$ — a single matrix-vector multiply generates correlated samples. It also computes log-determinants cheaply: $\log\det(A) = 2\sum_i \log \ell_{ii}$, needed in likelihoods and variational lower bounds.