Cholesky and Positive-Definiteness
The square root of a matrix, when one exists. Push the correlation too far and the factorisation fails — which is the test.
Overview
The factorisation
For a symmetric positive-definite matrix C there is a lower triangular L with:
C = L L'L is unique if its diagonal is taken positive, and it is the nearest thing a matrix has to a square root. For a 2×2 it is short enough to write out whole:
L11 = sqrt(C11)
L21 = C21 / L11
L22 = sqrt(C22 - L21^2)That last line carries the whole story. If C22 - L21^2 is negative there is no real square root, and the algorithm stops.
Cholesky and Positive-Definiteness
Worth knowing
C = L L' with L lower triangular. It is the closest thing a matrix has to a square root.Cholesky and Positive-Definiteness
A matrix square root, what it needs to exist, and why its failure is more useful than its success.
Positive-definite
A symmetric matrix is positive definite when x' C x > 0 for every non-zero x, which is equivalent to all its eigenvalues being strictly positive.
For a covariance matrix this is not a technicality. x' C x is the variance of the combination x of your variables, and a variance cannot be negative. A matrix that fails the test is not describing any distribution that exists.
Drag the correlation control toward ±1. Past the boundary the readout reports failure, and the reason is that a correlation of 1.2 is not a thing: no pair of random variables can be more than perfectly correlated.
The factorisation failing is the test. Cholesky costs about n^3 / 3 operations, roughly half a general LU factorisation, and about a tenth of what computing the eigenvalues costs. When a routine needs to know whether a matrix is positive definite, it attempts a Cholesky and watches for the negative square root.
Generating correlated randomness
This is where the factor earns its keep. Take independent standard normals z and compute L z. The result has covariance:
Cov(Lz) = L Cov(z) L' = L I L' = L L' = CExactly the covariance you asked for. The scatter in the visualisation is generated this way, and the two orange segments are the columns of L — the axes the independent noise gets mapped onto.
That single trick underlies Monte Carlo simulation of correlated risk factors, sampling from a multivariate normal, and the reparameterisation trick in variational autoencoders.
Where else it appears
Gaussian processes. Fitting one requires solving a system with the kernel matrix and computing its log-determinant. Cholesky gives both: the solve by triangular substitution, and the log-determinant as twice the sum of the logs of L's diagonal. It is the computational core of GP regression.
Linear systems with a symmetric positive-definite matrix — twice as fast as LU and numerically stable without pivoting.
Optimisation. Newton steps need to solve with the [Hessian](jacobian_and_hessian.html), and a successful Cholesky confirms the Hessian is positive definite, which confirms the step direction is a descent direction. A failure signals a saddle, and modified Newton methods react to exactly that signal.
Whitening. L^-1 x transforms correlated data into uncorrelated data with unit variance.
When it fails on data that should be fine
An estimated covariance matrix can come out not-quite-positive-definite for reasons that are arithmetic rather than conceptual.
More variables than observations. The estimate is singular by construction: with 100 variables and 50 samples the matrix has rank at most 50.
Perfectly collinear columns. One variable is a combination of others, so some direction has genuinely zero variance.
Floating point. A matrix that is positive definite in exact arithmetic can have a tiny negative eigenvalue after rounding.
The standard repairs: add a small multiple of the identity to the diagonal — *jitter*, or *ridge*, and the same idea as [ridge regression](ridge_and_lasso_regression.html); use a shrinkage estimator such as Ledoit-Wolf; or clip the negative eigenvalues to zero and reassemble.
Cholesky by hand, and the failure
The three formulas are short enough to run on a 2×2. Take C = [[4, 2], [2, 5]]:
L₁₁ = √4 = 2. L₂₁ = 2 / L₁₁ = 1. L₂₂ = √(5 − 1²) = √4 = 2.
So L = [[2, 0], [1, 2]], and multiplying back, L L' = [[4, 2], [2, 5]] — the original. The whole story is in that last square root: C₂₂ minus what the first column already accounts for must be non-negative for a real L to exist.
Now feed it a matrix that is symmetric but not positive definite, B = [[1, 2], [2, 1]]. The last step asks for √(1 − 2²) = √(−3), which has no real value, and the factorisation stops. That is not a bug — it is the test. B fails because x = (1, −1) gives x' B x = 1 − 4 + 1 = −2 < 0, so B is not a covariance of anything: no set of variables has a negative variance. The cheap n³/3 factorisation is how a routine decides positive-definiteness, precisely because it breaks exactly when the property fails.
That is also why an estimated covariance matrix sometimes fails the factorisation when it "should" work: with more variables than observations the estimate is singular by construction, and a genuinely collinear pair gives some direction zero variance, so the last square root reaches zero or below. The standard repair — adding a small multiple of the identity, the same idea as ridge regression — nudges every eigenvalue up just enough to make the square roots real again, at the cost of admitting you changed the model.
It is worth stressing why the factorisation is the test of choice: at about n³/3 operations it is half the cost of a general LU decomposition and a fraction of computing the eigenvalues, so a routine that needs to know whether a matrix is positive definite simply attempts the Cholesky and watches for the negative square root.
Where it goes wrong
Passing a non-symmetric matrix. Most implementations read only one triangle and will silently return nonsense.
Treating failure as a bug. It is usually the data telling you something.
Adding jitter without recording it. You have changed the model. Say by how much.
Using it on a matrix that is only positive semi-definite. A zero eigenvalue gives a zero on L's diagonal, and anything that then divides by it fails.
Half a matrix, and what it costs to have one
Cholesky finds a triangular L with L @ L.T equal to your matrix -- but only if the matrix is positive definite. Here is what that condition rules out.
Check yourself
0 of …Answer without scrolling back up.
When does a Cholesky factorisation exist?
How do you generate samples with a given covariance C?
Why is attempting a Cholesky the standard positive-definiteness test?
Cholesky and Positive-Definiteness
L is unique if its diagonal is taken positive, and it is the nearest thing a matrix has to a square root. For a 2×2 it is short enough to write out whole: