Positive numbers have square roots: $9 = 3 \cdot 3$. Symmetric positive-definite (SPD) matrices have a matrix analogue — the Cholesky decomposition $A = LL^\top$, where $L$ is lower triangular with a positive diagonal. It’s the fastest, most numerically stable way to work with covariance matrices, and it’s what lets you sample from a Gaussian, solve SPD systems quickly, and run the reparameterization trick in variational autoencoders.

A real-life example: generating correlated random data

You want to simulate fake-but-realistic data where height and weight are correlated — tall people tend to be heavier. You know the covariance matrix $\Sigma$. The trick: generate independent standard-normal numbers (easy), then multiply by the Cholesky factor $L$ of $\Sigma$. The output has exactly the correlation structure you asked for. Cholesky is the bridge from “independent noise” to “correlated samples.”

The decomposition

A symmetric positive-definite matrix $A$ factors uniquely as

where $L$ is lower-triangular with positive diagonal entries — the Cholesky factor. Multiplying out $LL^\top$ and matching entries gives explicit formulas; for the $3\times3$ diagonal:

and below the diagonal $l_{21} = a_{21}/l_{11}$, etc. Each entry is “back-calculated” from $A$ and the entries already found. The decomposition exists iff $A$ is SPD — a handy positive-definite test in itself.

Cholesky shapes noise into correlation

Start with a circle of independent unit vectors (uncorrelated noise). Multiply each by the Cholesky factor $L$ of a covariance matrix and the circle becomes a tilted ellipse — the samples now have the target variances and correlation. Watch the point trace the ellipse that $L$ carves out of the circle:

Why Cholesky is used everywhere

• Sampling from a Gaussian $\mathcal{N}(\boldsymbol\mu, \Sigma)$: draw $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, I)$, return $\boldsymbol\mu + L\mathbf{z}$ where $\Sigma = LL^\top$ (the demo above).
• Fast, stable solves: solving $A\mathbf{x} = \mathbf{b}$ for SPD $A$ via Cholesky is about twice as fast as generic LU and far more stable.
• Cheap determinant: $\det(A) = \det(L)^2 = \prod_i l_{ii}^2$ — no expensive general routine.
• Reparameterization trick: differentiable Gaussian sampling in VAEs relies on the $\boldsymbol\mu + L\mathbf{z}$ form so gradients flow through the sample.

NumPy

import numpy as np
 
# A symmetric positive-definite covariance matrix
A = np.array([[2.0, 1.2],
              [1.2, 1.5]])
 
L = np.linalg.cholesky(A)          # lower-triangular factor
print("L =\n", np.round(L, 4))
print("L @ L.T == A:", np.allclose(L @ L.T, A))
 
# determinant via Cholesky: product of squared diagonal
print("det via Cholesky:", round(np.prod(np.diag(L))**2, 4))
print("det direct       :", round(np.linalg.det(A), 4))
 
# generate correlated samples: mu + L @ z, with z ~ N(0, I)
rng = np.random.default_rng(0)
z = rng.standard_normal((2, 5000))
samples = L @ z                    # each column is a correlated 2-D sample
print("empirical covariance ≈ A:\n", np.round(np.cov(samples), 2))

import numpy as np
 
# A symmetric positive-definite covariance matrix
A = np.array([[2.0, 1.2],
              [1.2, 1.5]])
 
L = np.linalg.cholesky(A)          # lower-triangular factor
print("L =\n", np.round(L, 4))
print("L @ L.T == A:", np.allclose(L @ L.T, A))
 
# determinant via Cholesky: product of squared diagonal
print("det via Cholesky:", round(np.prod(np.diag(L))**2, 4))
print("det direct       :", round(np.linalg.det(A), 4))
 
# generate correlated samples: mu + L @ z, with z ~ N(0, I)
rng = np.random.default_rng(0)
z = rng.standard_normal((2, 5000))
samples = L @ z                    # each column is a correlated 2-D sample
print("empirical covariance ≈ A:\n", np.round(np.cov(samples), 2))

L =
 [[1.4142 0.    ]
 [0.8485 0.7778]]
L @ L.T == A: True
det via Cholesky: 1.56
det direct       : 1.56
empirical covariance ≈ A:
 [[1.99 1.19]
 [1.19 1.49]]

L =
 [[1.4142 0.    ]
 [0.8485 0.7778]]
L @ L.T == A: True
det via Cholesky: 1.56
det direct       : 1.56
empirical covariance ≈ A:
 [[1.99 1.19]
 [1.19 1.49]]

Why this matters for ML

• Gaussian processes and Bayesian models invert/solve SPD kernel and covariance matrices — Cholesky is the standard workhorse.
• VAEs sample latent variables via $\boldsymbol\mu + L\mathbf{z}$, keeping sampling differentiable.
• Whitening / preconditioning uses Cholesky factors to decorrelate features and speed up optimization.

🧪 Try It Yourself

Exercise 1 – Factor and verify

Exercise 2 – Cholesky as a positive-definite test

Exercise 3 – Determinant via Cholesky

Recap

• The Cholesky decomposition factors a symmetric positive-definite matrix as $A = LL^\top$ with $L$ lower-triangular, positive diagonal — a matrix square root.
• It exists iff $A$ is SPD, so a failed Cholesky is a positive-definiteness test.
• It powers Gaussian sampling ($\boldsymbol\mu + L\mathbf{z}$), fast SPD solves, cheap determinants ($\prod l_{ii}^2$), and the reparameterization trick.

Next: the general square-matrix factorization into eigenbasis coordinates — Eigendecomposition and Diagonalization.