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.
#include <cstdio>
#include <vector>
#include <cmath>
#include <random>
#include <stdexcept>
using Matrix = std::vector<std::vector<double>>;
// A = L L^T, column by column -- no linear-algebra library. Throws if A isn't PD.
Matrix cholesky(const Matrix& A) {
int n = static_cast<int>(A.size());
Matrix L(n, std::vector<double>(n, 0.0));
for (int j = 0; j < n; ++j) {
double s = 0.0;
for (int k = 0; k < j; ++k) s += L[j][k] * L[j][k];
double val = A[j][j] - s;
if (val <= 0) throw std::runtime_error("matrix is not positive definite");
L[j][j] = std::sqrt(val);
for (int i = j + 1; i < n; ++i) {
double s2 = 0.0;
for (int k = 0; k < j; ++k) s2 += L[i][k] * L[j][k];
L[i][j] = (A[i][j] - s2) / L[j][j];
}
}
return L;
}
bool isPd(const Matrix& A) {
try {
cholesky(A);
return true;
} catch (const std::runtime_error&) {
return false;
}
}
int main() {
Matrix cov = {{4.0, 2.0}, {2.0, 3.0}};
double mean[2] = {1.0, -2.0};
Matrix L = cholesky(cov);
std::mt19937 rng(0);
std::normal_distribution<double> stdNormal(0.0, 1.0);
const int nSamples = 1000;
std::vector<double> s0(nSamples), s1(nSamples);
for (int s = 0; s < nSamples; ++s) {
double z0 = stdNormal(rng), z1 = stdNormal(rng);
s0[s] = mean[0] + L[0][0] * z0 + L[0][1] * z1;
s1[s] = mean[1] + L[1][0] * z0 + L[1][1] * z1;
}
double m0 = 0.0, m1 = 0.0;
for (int s = 0; s < nSamples; ++s) { m0 += s0[s]; m1 += s1[s]; }
m0 /= nSamples; m1 /= nSamples;
double c00 = 0.0, c01 = 0.0, c11 = 0.0;
for (int s = 0; s < nSamples; ++s) {
c00 += (s0[s] - m0) * (s0[s] - m0);
c01 += (s0[s] - m0) * (s1[s] - m1);
c11 += (s1[s] - m1) * (s1[s] - m1);
}
c00 /= (nSamples - 1); c01 /= (nSamples - 1); c11 /= (nSamples - 1);
std::printf("sample_cov = [[%.3f, %.3f], [%.3f, %.3f]]\n", c00, c01, c01, c11);
std::printf("is_pd(cov) = %s\n", isPd(cov) ? "true" : "false");
return 0;
}- 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.