Positive Definite Matrices & Cholesky
The matrix version of "positive," and the fast square-root-like factorization it enables.
On this page
Beginner: a square symmetric matrix A is positive definite if, for every nonzero vector x, the number xᵀAx comes out positive. Think of it as a generalization of "this number is positive" to matrices — it's the condition that makes the quadratic bowl from section 1.16 always curve upward, everywhere, with a single clear minimum.
Intermediate: a positive semi-definite matrix relaxes this slightly — xᵀAx ≥ 0, allowing flat directions. Every covariance matrix in existence is positive semi-definite by mathematical construction (a variance can never be negative), which is why "is this a valid covariance matrix?" and "is this matrix PSD?" are the exact same question.
Advanced: for a positive-definite matrix, a special, cheaper, more numerically stable decomposition exists: Cholesky decomposition, A = LLᵀ, where L is lower-triangular. It's the matrix equivalent of taking a square root — and just like you can't take the square root of a negative number, Cholesky only exists for positive-definite matrices, which makes attempting it a fast, reliable PD test in itself.
Equivalent tests: all eigenvalues of A are positive, or every "leading principal minor" (determinant of the top-left k×k block, for every k) is positive.
Write out A = LLᵀ entry by entry and solve column by column. For the diagonal entry Ajj, matching it against the corresponding entry of LLᵀ gives Ajj = Σₖ₌₁ʲ Ljk² (only terms up to k=j survive, since L is lower-triangular). Solving for the one new unknown, Ljj:
Every off-diagonal entry below it, Aij for i>j, similarly matches Aij = Σₖ Lik Ljk, giving:
Working column by column, left to right, every L entry only ever depends on entries computed earlier — so the whole matrix can be filled in with no guessing, no iteration, and no need to solve a system. This only works because the square root above never hits a negative number — which is precisely what positive-definiteness guarantees.
Where this is used: this is literally the algorithm np.linalg.cholesky runs — a direct, non-iterative, column-by-column computation, which is exactly why Cholesky is faster than generic eigen-decomposition or LU (section 1.13) for this specific, common case.
Test A = [[4, 2], [2, 3]] three different ways, and confirm they agree. Eigenvalues: solving det(A − λI) = 0 gives λ ≈ 5.56 and λ ≈ 1.44 — both positive, so A is positive definite. Leading minors: the top-left 1×1 block is 4 > 0, and the full 2×2 determinant is (4)(3) − (2)(2) = 8 > 0 — both positive, agreeing with the eigenvalue test. Cholesky: L = [[2, 0], [1, √2]], and multiplying LLᵀ back out reproduces A exactly — confirming the factorization exists, which it could only do because A is positive definite.
This exact pattern — Cholesky, then multiply standard normal noise by L — is how every Gaussian process library and every "sample a correlated random variable" function is implemented under the hood.
import numpy as np
cov = np.array([[4.0, 2.0], [2.0, 3.0]]) # a covariance matrix
mean = np.array([1.0, -2.0])
L = np.linalg.cholesky(cov) # fails with LinAlgError if not PD
# Turn standard-normal noise into correlated samples from N(mean, cov):
z = np.random.randn(2, 1000) # standard normal, uncorrelated
samples = mean.reshape(-1, 1) + L @ z # now has the target mean & covariance
print(np.cov(samples)) # should be close to cov
# A quick PD check anywhere in your code:
def is_pd(A):
try:
np.linalg.cholesky(A)
return True
except np.linalg.LinAlgError:
return False- Gaussian processes require their kernel (covariance) matrix to be PD by definition — Cholesky is used both to sample from the process and to compute its log-likelihood efficiently.
- Kalman filters (used in GPS, robotics, and finance) rely on Cholesky decomposition of covariance matrices at every update step.
- Physical simulation — mass and stiffness matrices in structural engineering are positive definite by physical necessity (energy can't be negative), and Cholesky is the standard solver.
- Assuming a computed covariance-like matrix is exactly PD — floating-point rounding can leave it just barely non-PD (a tiny negative eigenvalue). The standard fix is adding a small jitter term (e.g.
+ 1e-6 * I) to the diagonal before factoring. - Confusing positive semi-definite (allows zero eigenvalues, no Cholesky) with positive definite (strictly positive eigenvalues, Cholesky guaranteed) — they are not interchangeable for this purpose.
Going deeper
Cholesky is roughly twice as fast as generic LU decomposition for the same matrix, purely by exploiting symmetry — it only ever needs to compute (and store) half the matrix.
At the master level: Cholesky decomposition is the standard, numerically preferred way to compute the log-determinant needed for a multivariate Gaussian log-likelihood: log|Σ| = 2·Σᵢ log(Lᵢᵢ) — summing the log of L's diagonal — which is both faster and dramatically more numerically stable than computing the determinant directly (section 1.12), and is exactly what every serious probabilistic modeling library does.