Positive (semi-)definite matrices, the five tests, Cholesky, Schur complements and the Loewner order

Before this lesson: lesson 06 (you need $S = Q\Lambda Q^\top$) and lesson 03 (determinants). Time: about 90 minutes, plus the notebook.

---

Why this matters #

A positive definite matrix is the matrix version of a positive number. Covariance matrices, kernel (Gram) matrices, the Hessian at a strict minimum, the Fisher information and finite-element stiffness matrices all have this property. Three practical questions follow. Is my matrix PD? (Is the loss a bowl here? Is my covariance valid?) How do I compute with it cheaply? (Cholesky: solve, sample, log-determinant.) How do I compare two such matrices? (The Loewner order: "the posterior uncertainty is smaller than the prior uncertainty in every direction".) Gaussian processes and Kalman filters use all three every time they condition on data, and the formula they use is a Schur complement.

1. Predict first #

2. The picture #

Start with one number $s$. The function $f(x) = s x^2$ is a parabola that opens upward when $s > 0$, is flat when $s = 0$, and opens downward when $s < 0$.

Now take a $2 \times 2$ diagonal matrix $S = \operatorname{diag}(\lambda_1, \lambda_2)$. Its energy (also called the quadratic form) is $x^\top S x = \lambda_1 x_1^2 + \lambda_2 x_2^2$.

• $\lambda_1, \lambda_2 > 0$: a bowl. The energy is positive in every direction.
• $\lambda_1 > 0 = \lambda_2$: a trough. The energy is zero along the $x_2$-axis.
• $\lambda_1 > 0 > \lambda_2$: a saddle. Up in one direction, down in another.

A general symmetric $S$ is the same picture, rotated. The spectral theorem gives $S = Q\Lambda Q^\top$ with $Q$ orthogonal (lesson 06). In the rotated coordinates $y = Q^\top x$:

So the signs of the eigenvalues decide everything. All positive: a bowl (PD). All non-negative: a bowl or a trough (PSD). Mixed signs: a saddle (indefinite). For a PD $S$, the level set $\{x : x^\top S x = 1\}$ is an ellipse (an ellipsoid in $n$ dimensions) with axes along the eigenvectors $q_i$ and semi-axis lengths $1/\sqrt{\lambda_i}$.

Interactive figure: loading…

Covariance. For a random vector $x$ with covariance $\Sigma$ and a fixed direction $w$, $\operatorname{Var}(w^\top x) = w^\top \Sigma w$. A variance cannot be negative, so covariance matrices are PSD.

3. Definitions and main results #

PD and PSD #

A symmetric matrix $S \in \mathbb{R}^{n \times n}$ is

• positive definite ($S \succ 0$) if $x^\top S x > 0$ for every $x \ne 0$;
• positive semidefinite ($S \succeq 0$) if $x^\top S x \ge 0$ for every $x$.

By the picture above, $S \succ 0$ exactly when all eigenvalues are $> 0$, and $S \succeq 0$ exactly when all eigenvalues are $\ge 0$. We always assume $S$ is symmetric. For a non-symmetric $M$, the energy only sees the symmetric part: $x^\top M x = x^\top \tfrac{1}{2}(M + M^\top)\, x$.

Sample covariance is PSD. Let $X \in \mathbb{R}^{N \times d}$ hold $N$ samples as rows, let $\bar x$ be the mean sample, and let $X_c$ be the centered data (row $j$ is $(x_j - \bar x)^\top$). The sample covariance is $S = X_c^\top X_c / (N-1)$. For any $w \in \mathbb{R}^d$:

the sample variance of the data projected on $w$. It is zero exactly when $X_c w = 0$: every centered sample is orthogonal to $w$, so the data lie in a hyperplane (for example, two features that are exact copies). The rows of $X_c$ sum to zero, so $\operatorname{rank} S \le N - 1$. With $d \ge N$ features, $S$ is always singular.

Theorem (the five tests) #

For a symmetric $S \in \mathbb{R}^{n\times n}$, these five statements are equivalent:

1. Eigenvalues: all $\lambda_i > 0$.
2. Energy: $x^\top S x > 0$ for every $x \ne 0$.
3. Gram form: $S = A^\top A$ for some matrix $A$ with independent columns.
4. Leading minors: $D_k = \det S_{[1:k,\,1:k]} > 0$ for $k = 1, \dots, n$ (the top-left $k \times k$ blocks).
5. Pivots: elimination without row exchanges gives $n$ pivots $d_1, \dots, d_n$, all $> 0$.

Proof of 1 ⇔ 2 ⇔ 3.

1. (1 ⇒ 2) Write $S = Q\Lambda Q^\top$. For $x \ne 0$, $y = Q^\top x \ne 0$ because $Q$ is invertible. Then $x^\top S x = \sum_i \lambda_i y_i^2 > 0$: every $\lambda_i > 0$ and some $y_i \ne 0$.
2. (2 ⇒ 1) Take a real unit eigenvector $q$ with $Sq = \lambda q$ (it exists because $S$ is symmetric). Then $0 < q^\top S q = \lambda\, q^\top q = \lambda$.
3. (3 ⇒ 2) $x^\top A^\top A x = \|Ax\|^2 \ge 0$. It is $0$ only if $Ax = 0$, and independent columns force $x = 0$.
4. (1 ⇒ 3) Set $A = \Lambda^{1/2} Q^\top$. It is square and invertible, so its columns are independent, and $A^\top A = Q \Lambda^{1/2}\Lambda^{1/2} Q^\top = S$.

Why the determinant tests work, in one line. Elimination does not change the leading minors, and at the end $D_k = d_1 d_2 \cdots d_k$. So the pivots are the ratios $d_k = D_k / D_{k-1}$, and "all $D_k > 0$" is the same as "all $d_k > 0$". The proof below makes this precise and links tests 4 and 5 to the energy.

Proof of 4 ⇔ 5, and how they join 1–3

1. Elimination without row exchanges subtracts a multiple of row $j$ from a later row $i > j$. If $i \le k$, this is a row operation inside the top-left $k\times k$ block, so the block's determinant does not change. If $i > k$, the block does not change at all. So every $D_k$ stays the same during elimination.
2. At the end, the top-left $k \times k$ block of $U$ is upper triangular with diagonal $d_1, \dots, d_k$. So $D_k = d_1 d_2 \cdots d_k$, and $d_k = D_k / D_{k-1}$ (with $D_0 = 1$).
3. Hence all $D_k > 0$ ⇔ all $d_k > 0$. (If all $D_k > 0$, each pivot is a ratio of positive numbers, so elimination never meets a zero pivot.)
4. (2 ⇒ 4) For $x$ whose entries after position $k$ are zero, $x^\top S x = x_{1:k}^\top S_{[1:k,1:k]}\, x_{1:k}$. So each leading block is PD, its eigenvalues are positive, and $D_k$, their product, is positive.
5. (5 ⇒ 2) With nonzero pivots, elimination gives $S = \tilde L U$ with $\tilde L$ unit lower triangular. Symmetry and the uniqueness of the LDU factorization give $U = D\tilde L^\top$ with $D = \operatorname{diag}(d_1, \dots, d_n)$, so $S = \tilde L D \tilde L^\top$. Then $x^\top S x = \sum_j d_j \big((\tilde L^\top x)_j\big)^2 > 0$ for $x \ne 0$, because $\tilde L^\top$ is invertible. $\blacksquare$

The PSD trap: leading minors are not enough #

Tests 1–3 carry over to PSD with $\ge$: eigenvalues $\ge 0$, energy $\ge 0$, and $S = A^\top A$ for any $A$ (dependent columns give the zero eigenvalues). Test 4 does not carry over.

	PD ($\succ 0$)	PSD ($\succeq 0$)
determinant test	the $n$ leading minors are $> 0$	all $2^n - 1$ principal minors are $\ge 0$
why the shortcut fails	—	a zero leading minor hides what comes after it
counterexample	—	$T = \begin{pmatrix} 0 & 0 \\ 0 & -1 \end{pmatrix}$: leading minors $0, 0$, but $e_2^\top T e_2 = -1$

A principal minor is the determinant of the submatrix with the same set of row and column indices, for example rows and columns $\{2\}$ or $\{1, 3\}$. The leading minors are only the index sets $\{1\}, \{1, 2\}, \dots$.

Why the shortcut fails, step by step. For PD, each $D_k > 0$ and the pivot $d_k = D_k/D_{k-1}$ is defined. When some $D_{k-1} = 0$, this ratio does not exist: the leading minors stop giving information about the later pivots. In $T$, $D_1 = 0$ already, so $D_2 = 0$ says nothing about the entry $-1$. The principal minor with index set $\{2\}$ sees it directly. The $3\times 3$ version is the minors 0, 0, 0 preset in the figure: leading minors $0, 0, 0$, but the principal minor $M_{22} = -1$.

Test 5 is stated for PD: a singular PSD matrix has a zero pivot, and plain elimination then stops (pivoted Cholesky handles this case).

Interactive figure: loading…

Cholesky #

Theorem. $S$ is PD if and only if $S = LL^\top$ with $L$ lower triangular with a positive diagonal. This $L$ is unique. It is the Cholesky factor of $S$.

By hand, $2\times 2$. Take $S = \begin{pmatrix} 4 & 2 \\ 2 & 3 \end{pmatrix}$ and write the unknown factor:

Match the entries one at a time:

1. Entry $(1,1)$: $\ell_{11}^2 = 4$, so $\ell_{11} = 2$.
2. Entry $(2,1)$: $\ell_{11}\ell_{21} = 2$, so $\ell_{21} = 2/2 = 1$.
3. Entry $(2,2)$: $\ell_{21}^2 + \ell_{22}^2 = 3$, so $\ell_{22}^2 = 3 - 1 = 2$ and $\ell_{22} = \sqrt 2$.

So $L = \begin{pmatrix} 2 & 0 \\ 1 & \sqrt2 \end{pmatrix}$, and the energy is a sum of squares: $x^\top S x = \|L^\top x\|^2 = (2x_1 + x_2)^2 + 2x_2^2$.

Where it breaks. Do the same for the predict matrix $\begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix}$: $\ell_{11} = 1$, $\ell_{21} = 2$, and step 3 needs $\ell_{22}^2 = 1 - 2^2 = -3$. No real number squares to $-3$, so the factorization stops. The number under the square root, $-3$, is the second pivot $1 - 2\cdot 2/1$. The energy shows the same thing: $x^\top S x = (x_1 + 2x_2)^2 - 3x_2^2$, which is $-3$ at $x = (-2, 1)$.

The general algorithm. The same matching works for any $n$.

1. For $i \ge j$, entry $(i, j)$ of $LL^\top$ is $S_{ij} = \sum_{k=1}^{j} L_{ik} L_{jk}$ (terms with $k > j$ vanish because $L_{jk} = 0$).
2. Diagonal ($i = j$): $L_{jj} = \sqrt{S_{jj} - \sum_{k<j} L_{jk}^2}$.
3. Below the diagonal ($i > j$): $L_{ij} = \big(S_{ij} - \sum_{k<j} L_{ik}L_{jk}\big) / L_{jj}$.
4. Go column by column, $j = 1, \dots, n$. Each formula uses only entries of earlier columns and the current $L_{jj}$.

Where PD is needed. From $S = \tilde L D \tilde L^\top$ (proof of the five tests), the matrix $\tilde L D^{1/2}$ is lower triangular with diagonal $\sqrt{d_j}$, so by uniqueness it is $L$. The number under the square root at step $j$ (the radicand) is exactly the $j$-th pivot $d_j$. If $S$ is not PD, some radicand is $\le 0$ and the algorithm stops. This is the standard practical PD test: np.linalg.cholesky raises LinAlgError.

Uniqueness of the Cholesky factor

Let $LL^\top = MM^\top$ with both diagonals positive. Then $M^{-1}L = M^\top L^{-\top}$. The left side is lower triangular and the right side is upper triangular, so both equal a diagonal matrix $\Delta$ with positive diagonal $L_{ii}/M_{ii}$. Then $L = M\Delta$, and $MM^\top = M\Delta^2M^\top$ gives $\Delta^2 = I$, so $\Delta = I$. So for a PD matrix and a fixed order of the variables, the Cholesky factor with a positive diagonal is unique. Other factors appear only with negative diagonal entries or reordered variables. $\blacksquare$

Cost and uses. About $n^3/3$ flops, half of LU's $2n^3/3$. After that, $Sx = b$ costs two triangular solves ($Ly = b$, then $L^\top x = y$), and $\log\det S = 2\sum_j \log L_{jj}$.

Sampling with the Cholesky factor #

To draw $x \sim \mathcal N(\mu, S)$: draw $z \sim \mathcal N(0, I)$ (independent standard normals) and set $x = \mu + Lz$. Why the covariance is $S$, in two steps.

1. $\mathbb E[Lz] = L\,\mathbb E[z] = 0$, so $\operatorname{Cov}(Lz) = \mathbb E\big[(Lz)(Lz)^\top\big]$.
2. $(Lz)(Lz)^\top = L\,zz^\top L^\top$, and the expectation is linear, so $\mathbb E\big[L\,zz^\top L^\top\big] = L\,\mathbb E[zz^\top]\,L^\top = L\,I\,L^\top = S$.

A linear map of a Gaussian vector is Gaussian, and a Gaussian is fixed by its mean and covariance. So $x = \mu + Lz$ has exactly the distribution $\mathcal N(\mu, S)$. Step 2 never used that $L$ is triangular: any $M$ with $MM^\top = S$ works, for example the symmetric square root $S^{1/2} = Q\Lambda^{1/2}Q^\top$. Cholesky is the cheapest such $M$. With $N$ samples stored as rows ($Z$ is $N \times d$): $X = \mathbf{1}\mu^\top + Z L^\top$.

Interactive figure: loading…

The Schur complement #

The picture. Eliminate the first variable of $S = \begin{pmatrix} 4 & 2 \\ 2 & 3 \end{pmatrix}$: subtract $2/4$ times row 1 from row 2. The entry that is left is $3 - 2\cdot 2/4 = 2$. It is the second pivot, and it is $\ell_{22}^2$ from the hand computation above. Do the same with whole blocks instead of single entries, and the number $3 - 2\cdot 2/4$ becomes a matrix.

Let $M = \begin{pmatrix} A & B \\ B^\top & C \end{pmatrix}$ be symmetric, with $A$ of size $k\times k$ and invertible. The Schur complement of $A$ in $M$ is

Block elimination subtracts $B^\top A^{-1}$ times the first block row from the second one: the second block row becomes $\begin{pmatrix} 0 & C - B^\top A^{-1}B \end{pmatrix}$. So the pivots of $M$ are the pivots of $A$ followed by the pivots of $M/A$, and $\det M = \det A \cdot \det(M/A)$.

Theorem (Schur complement test). $M \succ 0 \iff A \succ 0$ and $C - B^\top A^{-1}B \succ 0$.

Proof.

1. Let $E = \begin{pmatrix} I & A^{-1}B \\ 0 & I \end{pmatrix}$. It is invertible (block triangular with identity blocks). Multiply out: $E^\top \begin{pmatrix} A & 0 \\ 0 & M/A \end{pmatrix} E = \begin{pmatrix} A & B \\ B^\top & B^\top A^{-1}B + M/A \end{pmatrix} = M$.
2. So for any $x$, with $y = Ex$: $x^\top M x = y^\top \begin{pmatrix} A & 0 \\ 0 & M/A \end{pmatrix} y$. Because $E$ is invertible, $x \ne 0 \iff y \ne 0$.
3. Hence $M \succ 0$ exactly when the block-diagonal matrix is PD, that is, when $A \succ 0$ and $M/A \succ 0$. (If $M \succ 0$, then $A \succ 0$ as a leading block, so $A^{-1}$ exists.) $\blacksquare$

With $A \succ 0$, the same argument gives $M \succeq 0 \iff M/A \succeq 0$.

Gaussian conditioning (the ML use). Split a Gaussian vector into two parts, $x = (x_a, x_b) \sim \mathcal N(\mu, \Sigma)$ with $\Sigma = \begin{pmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ba} & \Sigma_{bb} \end{pmatrix}$. After you observe $x_a$,

The conditional covariance is the Schur complement $\Sigma/\Sigma_{aa}$. It is PSD (by the test), and it is $\preceq \Sigma_{bb}$, because the subtracted term $\Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab}$ is PSD: observing $x_a$ never increases the uncertainty about $x_b$ in any direction.

Small check with Cholesky. Let $\Sigma = \begin{pmatrix} 4 & 2 \\ 2 & 3 \end{pmatrix}$ and $\mu = 0$. With $L$ from above, $x_1 = 2z_1$ and $x_2 = z_1 + \sqrt2\, z_2$. Observing $x_1$ fixes $z_1 = x_1/2$. What is left is $x_2 = x_1/2 + \sqrt2\, z_2$: mean $x_1/2 = \Sigma_{21}\Sigma_{11}^{-1}x_1$ and variance $2 = 3 - 2\cdot 2/4$, the Schur complement. A Gaussian-process posterior is this formula with $\Sigma_{aa} = K + \sigma^2 I$ (training points) and $x_b$ = function values at test points.

The Loewner order #

For symmetric $A, B$: $A \preceq B$ means $B - A \succeq 0$, that is, $x^\top A x \le x^\top B x$ for every $x$. This is the Loewner order. The PSD matrices form a convex cone: sums and non-negative multiples of PSD matrices are PSD. The cone is not a subspace: if $S \ne 0$ is PSD, then $-S$ is not.

Theorem. $\preceq$ is a partial order on symmetric matrices (reflexive, antisymmetric, transitive), and it is not a total order: some pairs cannot be compared. Example: $A = \operatorname{diag}(2, 1)$ and $B = \operatorname{diag}(1, 2)$. $B - A = \operatorname{diag}(-1, 1)$ is indefinite, so neither $A \preceq B$ nor $B \preceq A$.

Proof that $\preceq$ is a partial order

1. Reflexive: $A - A = 0 \succeq 0$.
2. Antisymmetric: let $B - A \succeq 0$ and $A - B \succeq 0$. The symmetric matrix $D = B - A$ has all eigenvalues $\ge 0$ and $\le 0$, so all are $0$, and $D = Q\,0\,Q^\top = 0$. (Symmetry matters here: a $90°$ rotation $R$ has $x^\top R x = 0$ for all $x$ but $R \ne 0$.)
3. Transitive: $C - A = (C - B) + (B - A)$ is a sum of PSD matrices, so it is PSD.
4. Not total: the pair $\operatorname{diag}(2, 1)$, $\operatorname{diag}(1, 2)$ above. $\blacksquare$

Useful facts. If $A \preceq B$, then $C^\top A C \preceq C^\top B C$ for any matrix $C$ (apply the energy to $Cx$), and $\lambda_k(A) \le \lambda_k(B)$ for every $k$ (Courant–Fischer, lesson 06). For PD matrices, $A \preceq B$ implies $B^{-1} \preceq A^{-1}$. Geometrically, for PD $A \preceq B$ the ellipsoid $\{x : x^\top B x \le 1\}$ lies inside $\{x : x^\top A x \le 1\}$. Ridge is the simplest example: $S \preceq S + \lambda I$ for every $\lambda \ge 0$.

The covariance picture. For a covariance it is more natural to draw $\mathcal E_M = \{x : x^\top M^{-1} x \le 1\} = \{M^{1/2}z : \|z\| \le 1\}$, with half-axes $\sqrt{\lambda_i}$ along the eigenvectors. Its width in a unit direction $u$ is $\max_{\|z\| \le 1} u^\top M^{1/2} z = \|M^{1/2}u\| = \sqrt{u^\top M u}$, the standard deviation along $u$. So $A \preceq B$ says: in every direction, $\mathcal E_A$ is at most as wide as $\mathcal E_B$, that is, $\mathcal E_A \subseteq \mathcal E_B$.

Interactive figure: loading…

The generalized eigenproblem $Sx = \lambda Bx$

Linear discriminant analysis (LDA) maximizes $w^\top S_B w / w^\top S_W w$ (between-class over within-class scatter), and canonical correlation analysis (CCA) has the same form. The critical points solve $Sx = \lambda Bx$ with $S$ symmetric and $B \succ 0$. Cholesky reduces this to a symmetric problem of lesson 06:

1. Factor $B = LL^\top$ and put $y = L^\top x$, so $x = L^{-\top}y$.
2. Then $Sx = \lambda LL^\top x$ becomes $L^{-1}SL^{-\top}y = \lambda y$.
3. $C = L^{-1}SL^{-\top}$ is symmetric, so the $\lambda_i$ are real and the $y_i$ can be chosen orthonormal. The $x_i = L^{-\top}y_i$ are then $B$-orthonormal: $x_i^\top B x_j = \delta_{ij}$.

Example: $S = \begin{pmatrix} 2 & 1 \\ 1 & 3 \end{pmatrix}$, $B = \begin{pmatrix} 4 & 2 \\ 2 & 3 \end{pmatrix}$ gives $\lambda = 0.5$ and $1.25$. In SciPy: scipy.linalg.eigh(S, B).

Answers to section 1. (1) No. The eigenvalues are $3$ and $-1$, and $x = (1, -1)$ gives $x^\top S x = -2$. Cholesky stops with radicand $-3$. (2) No. $e_2^\top T e_2 = -1 < 0$. The leading-minor test works only for PD, with strict inequalities.

4. Worked example #

Cholesky, column by column.

• Column 1: $L_{11} = \sqrt 4 = 2$, $L_{21} = 2/2 = 1$, $L_{31} = 2/2 = 1$.
• Column 2: $L_{22} = \sqrt{5 - 1^2} = 2$, $L_{32} = (3 - 1\cdot 1)/2 = 1$.
• Column 3: $L_{33} = \sqrt{6 - 1^2 - 1^2} = 2$.

The five tests agree.

• Pivots = radicands: $4, 4, 4$, all $> 0$.
• Leading minors: $D_1 = 4$, $D_2 = 4\cdot 5 - 2\cdot 2 = 16$, $D_3 = \det S = (\det L)^2 = 8^2 = 64$. The ratios $4/1,\ 16/4,\ 64/16$ give the pivots back.
• Energy as a sum of squares, from $x^\top S x = \|L^\top x\|^2$: $x^\top S x = (2x_1 + x_2 + x_3)^2 + (2x_2 + x_3)^2 + (2x_3)^2$, which is $> 0$ unless $x = 0$.
• Gram form: $A = L^\top$ is triangular with a nonzero diagonal, so its columns are independent.
• Eigenvalues (numerically): $2.32,\ 2.79,\ 9.89$. Check: sum $= 15 = \operatorname{tr} S$, product $= 64 = \det S$.

Schur complement. Split after the first row and column: $A = (4)$, $B = (2, 2)$, $C = \begin{pmatrix} 5 & 3 \\ 3 & 6 \end{pmatrix}$. Then

Its leading minors are $4$ and $16$, so it is PD, and $4 > 0$: by the Schur test, $S$ is PD. Its Cholesky factor is $\begin{pmatrix} 2 & 0 \\ 1 & 2 \end{pmatrix}$, the bottom-right block of $L$.

Sampling. $x = Lz$ reads $x_1 = 2z_1$, $x_2 = z_1 + 2z_2$, $x_3 = z_1 + z_2 + 2z_3$. Two entries by hand: $\operatorname{Var}(x_2) = 1 + 4 = 5$ and $\operatorname{Cov}(x_1, x_2) = 2 \cdot 1 = 2$. ✓

A failure. Change $S_{33}$ from $6$ to $1$. Columns 1 and 2 of $L$ do not change, and the last radicand is $1 - 1 - 1 = -1$. Cholesky stops: the matrix is not PD. Indeed $D_3 = 4 \cdot 4 \cdot (-1) = -16 < 0$. (Preset S₃₃ = 1 in the five-tests figure.)

Loewner. $\lambda_{\min} I \preceq S \preceq \lambda_{\max} I$, here $2.32\, I \preceq S \preceq 9.89\, I$. But $S$ and its diagonal $D_S = \operatorname{diag}(4, 5, 6)$ are incomparable: $D_S - S \ne 0$ has trace $0$, so it has a positive and a negative eigenvalue (here $-4.70,\ 1.70,\ 3$).

5. In code #

NumPy checks the hand computations above. The notebook for this lesson has these cells and the exercises.

import numpy as np

S = np.array([[4., 2., 2.],
              [2., 5., 3.],
              [2., 3., 6.]])
L = np.linalg.cholesky(S)                  # lower triangular, positive diagonal
print(L)                                   # [[2,0,0],[1,2,0],[1,1,2]]
print(np.diag(L) ** 2)                     # the pivots: [4. 4. 4.]
print(np.array([np.linalg.det(S[:k, :k]) for k in (1, 2, 3)]).round(6))  # [ 4. 16. 64.]
print(np.linalg.eigvalsh(S).round(2))      # [2.32 2.79 9.89]

S_bad = S.copy()
S_bad[2, 2] = 1.0
try:
    np.linalg.cholesky(S_bad)
except np.linalg.LinAlgError as err:
    print("not PD:", err)                  # not PD: Matrix is not positive definite

T = np.array([[0., 0.], [0., -1.]])        # leading minors 0 and 0
print(np.linalg.eigvalsh(T))               # [-1.  0.]: not PSD

A, B, C = S[:1, :1], S[:1, 1:], S[1:, 1:]  # Schur complement of the 1x1 block
schur = C - B.T @ np.linalg.solve(A, B)
print(schur)                               # [[4. 2.] [2. 5.]]
print(np.linalg.cholesky(schur))           # [[2. 0.] [1. 2.]] = L[1:, 1:]

D_S = np.diag(np.diag(S))
print(np.linalg.eigvalsh(D_S - S).round(2))   # [-4.7  1.7  3. ]: incomparable

rng = np.random.default_rng(0)
X = rng.multivariate_normal(np.zeros(3), S, size=200_000, method="cholesky")
print(np.cov(X, rowvar=False).round(1))    # close to S

6. Where it shows up #

• Kernel and Gram matrices are PSD. If $K_{ij} = \varphi(x_i)^\top \varphi(x_j)$ for a feature map $\varphi$, then $K = \Phi\Phi^\top$ with rows $\varphi(x_i)^\top$, and $c^\top K c = \|\Phi^\top c\|^2 \ge 0$. A function $k(x, x')$ is a valid kernel exactly when every matrix $[k(x_i, x_j)]$ is PSD; this is the condition of Mercer's theorem. Kernel ridge regression, SVMs and Gaussian processes rely on it.
• Gaussian processes. GP regression factors $K + \sigma^2 I = LL^\top$ once, solves with two triangular solves, and gets $\log\det$ from $\operatorname{diag}(L)$. The posterior covariance at test points is a Schur complement. When Cholesky fails, add jitter such as $10^{-6} I$.
• Optimization. A PD Hessian at a stationary point means a strict local minimum. Newton solvers try a Cholesky of the (shifted) Hessian. The Gauss–Newton matrix $J^\top J$ is PSD by the Gram test.
• Data assimilation (ML for science). The error covariances $B$ and $R$ in 3D-Var or a Kalman filter must be PD. The Kalman posterior covariance $P_a = P_f - P_f H^\top (H P_f H^\top + R)^{-1} H P_f \preceq P_f$ is the Schur complement of the joint covariance of state and observation: data never increase the uncertainty in any direction. Correlated ensemble perturbations are drawn as $\mu + Lz$.
• Statistics. The Cramér–Rao bound is a Loewner inequality: an unbiased estimator has $\operatorname{Cov}(\hat\theta) \succeq I(\theta)^{-1}$, with $I(\theta)$ the Fisher information.
• Simulation. Finite-element stiffness matrices (with boundary conditions) are symmetric PD, solved by Cholesky or conjugate gradients (lesson 14).

7. Common mistakes #

• "Positive entries mean PD." $\begin{pmatrix}1&2\\2&1\end{pmatrix}$ is indefinite; $\begin{pmatrix}2&-1\\-1&2\end{pmatrix}$ is PD.
• "$\det S > 0$ means PD." $\operatorname{diag}(-1, -1)$ has determinant $1$.
• Leading minors $\ge 0$ as a PSD test. Wrong (see $T$). Use all principal minors, or eigenvalues $\ge -\text{tol}$.
• Forgetting symmetry. $M = \begin{pmatrix}1&3\\0&1\end{pmatrix}$ has eigenvalues $1, 1$, but $x = (1,-1)$ gives $x^\top M x = -1$. Test the symmetric part $\tfrac12(M + M^\top)$.
• Trusting random energy checks. They can find a negative direction but never prove PD: a singular PSD matrix passes.
• Sampling with the wrong factor. $x = L^\top z$ has covariance $L^\top L \ne S$ (in the worked example, $L^\top L$ has first row $(6, 3, 2)$). With samples as rows, use $Z L^\top$, not $Z L$.
• Checking only the Schur complement. The test needs both $A \succ 0$ and $M/A \succ 0$. For $M = \begin{pmatrix} -1 & 1 \\ 1 & 0 \end{pmatrix}$ with $A = (-1)$, the Schur complement is $0 - 1\cdot 1/(-1) = 1 > 0$, but $M$ is indefinite ($\det M = -1$), because $A < 0$.
• "$A \preceq B$ implies $A^2 \preceq B^2$." False. $A = \operatorname{diag}(2, 0)$ and $B = \begin{pmatrix}3&1\\1&1\end{pmatrix}$: $B - A = \begin{pmatrix}1&1\\1&1\end{pmatrix} \succeq 0$, but $B^2 - A^2 = \begin{pmatrix}6&4\\4&2\end{pmatrix}$ has determinant $-4 < 0$. (The square root does keep the order.)

8. Check yourself #

1. Is $\begin{pmatrix}1&2\\2&5\end{pmatrix}$ PD? Give its Cholesky factor and its energy as a sum of squares.
   Answer

   Yes: $D_1 = 1$ and $D_2 = 5 - 4 = 1$. $L = \begin{pmatrix}1&0\\2&1\end{pmatrix}$, so $x^\top S x = (x_1 + 2x_2)^2 + x_2^2$.
2. For which $a$ is $\begin{pmatrix}1&a\\a&1\end{pmatrix}$ PD? PSD?
   Answer

   The eigenvalues are $1 \pm a$. PD for $|a| < 1$, PSD for $|a| \le 1$. (It is the correlation matrix of two variables with correlation $a$.)
3. Show that $S \succ 0$ implies $S^{-1} \succ 0$.
   Answer

   $S^{-1}$ is symmetric with eigenvalues $1/\lambda_i > 0$. Or: for $x \ne 0$ put $y = S^{-1}x \ne 0$; then $x^\top S^{-1} x = y^\top S y > 0$.
4. You have $N = 5$ samples of $d = 10$ features. Can the sample covariance be PD?
   Answer

   No. $\operatorname{rank} S \le N - 1 = 4 < 10$, so $S$ has at least 6 zero eigenvalues. Regularize ($S + \lambda I$, shrinkage) before you invert it.
5. For which $c$ is $M = \begin{pmatrix}1&1&0\\1&2&1\\0&1&c\end{pmatrix}$ PD? Use the Schur complement of the top-left $2\times 2$ block.
   Answer

   $A = \begin{pmatrix}1&1\\1&2\end{pmatrix}$ is PD ($D_1 = 1$, $D_2 = 1$) and $A^{-1} = \begin{pmatrix}2&-1\\-1&1\end{pmatrix}$. With $B = (0, 1)^\top$: $B^\top A^{-1} B = 1$, so $M/A = c - 1$. $M$ is PD exactly when $c > 1$. (Check: $\det M = c - 1$.)
6. If $A \preceq B$, must $\lambda_{\max}(A) \le \lambda_{\max}(B)$? Is the converse true?
   Answer

   Yes: $\lambda_{\max}(A) = \max_{\|x\|=1} x^\top A x \le \max_{\|x\|=1} x^\top B x = \lambda_{\max}(B)$. The converse is false: $\operatorname{diag}(2,1)$ and $\operatorname{diag}(1,2)$ have the same eigenvalues but are incomparable.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Define a positive definite matrix.

   Show answer

2. Name the three PD tests that use eigenvalues, the energy, or a factor $A$.

   Show answer

3. Name the two PD tests that use determinants or elimination, and how they are linked.

   Show answer

4. Why can leading minors $\ge 0$ not prove that a matrix is PSD?

   Show answer

5. Compute the Cholesky factor of $\begin{pmatrix} 4 & 2 \\ 2 & 3 \end{pmatrix}$.

   Show answer

6. When does a Cholesky factorization $S = LL^\top$ (positive diagonal) exist, and is it unique?

   Show answer

7. Why does $x = \mu + Lz$ with $z \sim \mathcal N(0, I)$ have covariance $LL^\top$?

   Show answer

8. Why is every Gram matrix $K = \Phi\Phi^\top$ PSD?

   Show answer

9. State the Schur complement test for $M = \begin{pmatrix} A & B \\ B^\top & C \end{pmatrix}$.

   Show answer

10. What does $A \preceq B$ mean in the Loewner order?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you check that a matrix is PD in practice?" — Try a Cholesky factorization: about $n^3/3$ flops, and it fails exactly when a pivot is $\le 0$. Use eigvalsh if you need the margin $\lambda_{\min}$. Never use $\det$ alone.
• "My GP kernel matrix fails Cholesky. Why, and what do you do?" — It is numerically singular (close inputs, a very smooth kernel). Add jitter $\varepsilon I$: this moves it up in the Loewner order and lowers its condition number.
• "Derive the conditional of a Gaussian." — Block-eliminate the covariance: the conditional mean is $\mu_b + \Sigma_{ba}\Sigma_{aa}^{-1}(x_a - \mu_a)$ and the conditional covariance is the Schur complement $\Sigma_{bb} - \Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab}$.
• "Is the sum of two PD matrices PD? The product?" — The sum, yes. The product $AB$ is usually not symmetric, but its eigenvalues are real and positive: $AB$ is similar to $A^{1/2} B A^{1/2} \succ 0$.

10. Practice #

The notebook for this lesson runs in Colab. It has five exercises, each with a check cell and a folded solution:

1. Sample covariance of a real dataset (load_iris): build it from the centered data, check that it is PSD, and see the rank limit $N - 1$.
2. The five tests as functions (is_pd_eigen, is_pd_minors, elimination_pivots, …), checked to agree on random PD and non-PD matrices.
3. Write cholesky column by column and make it report the step where it fails.
4. Sample correlated Gaussians with $X = \mathbf{1}\mu^\top + ZL^\top$, and check the sample covariance.
5. loewner_leq(A, B): ridge ($S \preceq S + \lambda I$), the posterior–prior pair, and an incomparable pair.

A demo section checks the Schur complement test and the Gaussian conditioning formula by sampling.

11. Watch and read #

5. Positive Definite and Semidefinite MatricesGilbert Strang · MIT 18.065 (2018), lecture 5 · 45 min

27. Positive Definite Matrices and MinimaGilbert Strang · MIT 18.06 (2010), lecture 27 · 51 min

Symmetric Matrices and Positive Definiteness | MIT 18.06SC Linear Algebra, Fall 2011MIT 18.06SC (2011) · recitation, a worked problem · 13 min

• Gilbert Strang, MIT 18.065 (2018), lecture 5: Positive definite and semidefinite matrices (textbook section I.7).
• Nick Higham, "What is a Cholesky factorization?": uniqueness, failure as a PD test, the PSD case, stability.
• Francis Bach, "Playing with positive definite matrices – I: matrix monotony and convexity": the Loewner order, and which matrix functions keep it.
• Rasmussen and Williams, Gaussian Processes for Machine Learning, free PDF: §2.2 and appendix A.2 (Gaussian conditioning, the Schur complement in GP regression).