Lesson overview | Previous part | Next part

Positive Definite Matrices: Part 12: Exercises to Appendix F: Quick Reference

12. Exercises

Exercise 1 * - Verify PD Axioms for $A^\top A + \lambda I$

Let $A \in \mathbb{R}^{m \times n}$ with $m \geq n$ and $\lambda > 0$. Define $B = A^\top A + \lambda I$.

(a) Prove $B \succ 0$ using the definition $\mathbf{x}^\top B\mathbf{x} > 0$ for $\mathbf{x} \neq \mathbf{0}$.

(b) What is the smallest eigenvalue of $B$? Express it in terms of the singular values of $A$ and $\lambda$.

(c) What happens as $\lambda \to 0$? Under what condition on $A$ does $A^\top A$ itself become PD?

(d) Show $B^{-1} \preceq \frac{1}{\lambda}I$ using the Loewner order.

(e) AI connection: Identify $B = A^\top A + \lambda I$ as the normal equations matrix for ridge regression. Why does $\lambda > 0$ guarantee $B \succ 0$ regardless of the rank of $A$?

---

Exercise 2 * - Implement Cholesky from Scratch

(a) Implement the Cholesky algorithm (Algorithm 4.2) in pure NumPy without using np.linalg.cholesky. Your function should return the lower triangular $L$ satisfying $A = LL^\top$, or raise a ValueError if $A$ is not PD.

(b) Test on $A = \begin{pmatrix}4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98\end{pmatrix}$ and verify $LL^\top = A$.

(c) Compare the result with np.linalg.cholesky - they should agree to machine precision.

(d) Test what happens when you apply your function to an indefinite matrix $A = \begin{pmatrix}1 & 3 \\ 3 & 1\end{pmatrix}$.

(e) Efficiency: What is the flop count of your implementation for $n \times n$ input? Compare to LU factorization.

---

Exercise 3 * - Sylvester's Criterion

Let $A(\alpha) = \begin{pmatrix}1 & \alpha & 0 \\ \alpha & 2 & \alpha \\ 0 & \alpha & 4\end{pmatrix}$.

(a) Compute the three leading principal minors $\Delta_1, \Delta_2, \Delta_3$ as functions of $\alpha$.

(b) Use Sylvester's criterion to find all values of $\alpha \in \mathbb{R}$ for which $A(\alpha) \succ 0$.

(c) At the boundary values of $\alpha$ (where $A$ transitions from PD to indefinite), what is the nullspace of $A$? Interpret this geometrically.

(d) Verify your answer computationally by plotting the smallest eigenvalue of $A(\alpha)$ as a function of $\alpha$.

---

Exercise 4 ** - LDL^T Factorization

(a) Implement the LDL^T algorithm (Algorithm 4.4) in NumPy. Your function returns unit lower triangular $L$ and diagonal $D$ such that $A = LDL^\top$.

(b) Test on $A = \begin{pmatrix}4 & 2 & 1 \\ 2 & 5 & 3 \\ 1 & 3 & 6\end{pmatrix}$. Verify $LDL^\top = A$.

(c) Using your factorization, solve $A\mathbf{x} = \mathbf{b}$ for $\mathbf{b} = (1, 2, 3)^\top$ by: (i) solve $L\mathbf{y} = \mathbf{b}$, (ii) solve $D\mathbf{w} = \mathbf{y}$, (iii) solve $L^\top\mathbf{x} = \mathbf{w}$.

(d) Explain why LDL^T avoids the numerical issues of Cholesky when $A$ is indefinite.

(e) AI connection: Why do second-order optimizers (like KFAC) prefer working with $D$ (the pivots) rather than $\sqrt{D}$ (the Cholesky diagonal)?

---

Exercise 5 ** - Schur Complement and Block PD Test

Let $M = \begin{pmatrix}A & B \\ B^\top & D\end{pmatrix}$ where:

$$A = \begin{pmatrix}4 & 1 \\ 1 & 3\end{pmatrix}, \quad B = \begin{pmatrix}1 \\ 2\end{pmatrix}, \quad D = \begin{pmatrix}5\end{pmatrix}.$$

(a) Compute the Schur complement $S = D - B^\top A^{-1}B$.

(b) Use Theorem 6.2 to determine whether $M \succ 0$.

(c) Compute $\det M$ using $\det M = \det A \cdot \det S$ and verify numerically.

(d) Use the block inverse formula to compute $M^{-1}$. Verify $MM^{-1} = I$.

(e) Gaussian conditioning: Interpret $M$ as a joint covariance matrix $\Sigma = \begin{pmatrix}\Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22}\end{pmatrix}$. What is the conditional covariance $\Sigma_{1|2}$? Interpret: does conditioning reduce uncertainty?

---

Exercise 6 ** - Log-Det via Cholesky and Gradient Verification

(a) Compute $\log\det A$ for $A = \begin{pmatrix}3 & 1 & 0 \\ 1 & 4 & 2 \\ 0 & 2 & 5\end{pmatrix}$ using: (i) direct eigenvalues, (ii) Cholesky diagonal formula, (iii) np.log(np.linalg.det(A)). Verify all three agree.

(b) Implement a function log_det(A) using Cholesky. Compare speed vs eigenvalue method for $n = 100, 500, 1000$ random PD matrices.

(c) Numerically verify the gradient formula $\nabla_A \log\det A = A^{-1}$ using finite differences: $\frac{f(A + hE_{ij}) - f(A)}{h} \approx (A^{-1})_{ij}$.

(d) For the Gaussian log-likelihood: given 50 samples from $\mathcal{N}(\mathbf{0}, \Sigma_{\text{true}})$ where $\Sigma_{\text{true}} = A$ from part (a), compute the log-likelihood $\ell(\Sigma) = -\frac{n}{2}\log\det\Sigma - \frac{1}{2}\text{tr}(\Sigma^{-1}S)$ where $S = \frac{1}{n}\sum_i\mathbf{x}_i\mathbf{x}_i^\top$ is the sample covariance.

---

Exercise 7 *** - Gram Matrix / Kernel Matrix and Mercer Check

(a) Given 5 points $\mathbf{x}_1, \ldots, \mathbf{x}_5 \in \mathbb{R}^3$ (design your own interesting configuration), compute the Gram matrix $G = XX^\top$.

(b) Verify $G \succeq 0$ by checking eigenvalues. Explain why $\text{rank}(G) \leq 3$.

(c) Implement and test three kernels on these points: (i) linear $k(\mathbf{x},\mathbf{z}) = \mathbf{x}^\top\mathbf{z}$, (ii) RBF $k(\mathbf{x},\mathbf{z}) = \exp(-\|\mathbf{x}-\mathbf{z}\|^2/2)$, (iii) polynomial $k(\mathbf{x},\mathbf{z}) = (\mathbf{x}^\top\mathbf{z}+1)^2$. For each, compute the $5 \times 5$ kernel matrix and verify PSD.

(d) Now try the "not-a-kernel" $k(\mathbf{x},\mathbf{z}) = \exp(-\|\mathbf{x}-\mathbf{z}\|)$ (Laplacian kernel IS a valid kernel) vs $k(\mathbf{x},\mathbf{z}) = \cos(\|\mathbf{x}-\mathbf{z}\|)$ (NOT always a valid kernel). Determine empirically whether each produces a PSD Gram matrix.

(e) AI connection: Implement a simple kernel regression: given training inputs $\mathbf{X}$ and outputs $\mathbf{y}$, compute predictions at test points $X_*$ using $\hat{\mathbf{y}}_* = K_{*n}(K_{nn} + 0.1 I)^{-1}\mathbf{y}$ (Nadaraya-Watson/GPR predictive mean).

---

Exercise 8 *** - Cholesky Reparameterization for VAE Sampling

(a) Implement a function sample_mvn(mu, L, n_samples) that samples $n$ vectors from $\mathcal{N}(\boldsymbol{\mu}, LL^\top)$ using the reparameterization trick $\mathbf{z} = \boldsymbol{\mu} + L\boldsymbol{\epsilon}$, $\boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0},I)$.

(b) Let $\boldsymbol{\mu} = (1, 2)^\top$ and $\Sigma = \begin{pmatrix}4 & 2 \\ 2 & 3\end{pmatrix}$. Compute $L = \text{chol}(\Sigma)$. Sample 1000 points, compute the sample mean and covariance, and verify they match $\boldsymbol{\mu}$ and $\Sigma$ (up to sampling noise).

(c) Simulate a VAE encoder: the encoder takes a scalar input $x$ and outputs $\boldsymbol{\mu}_\phi(x) = Wx + b$ and a lower triangular $L_\phi(x)$ (a diagonal $L$ for simplicity). Sample $\mathbf{z}$ and pass it through a simple decoder $\hat{x} = c^\top\mathbf{z}$. Compute the ELBO:

$$\text{ELBO} = \mathbb{E}_{\mathbf{z}}[\log p(x|\mathbf{z})] - D_{\text{KL}}(\mathcal{N}(\boldsymbol{\mu}, \Sigma) \| \mathcal{N}(\mathbf{0}, I)).$$

where $D_{\text{KL}}(\mathcal{N}(\boldsymbol{\mu}, \Sigma) \| \mathcal{N}(\mathbf{0},I)) = \frac{1}{2}[\text{tr}(\Sigma) + \boldsymbol{\mu}^\top\boldsymbol{\mu} - d - \log\det\Sigma]$.

(d) Verify the KL formula using both the analytical expression and Monte Carlo estimation from samples.

(e) Gradient check: Confirm that the reparameterization trick allows gradients to flow through $\boldsymbol{\mu}$ and $L$ (not $\boldsymbol{\epsilon}$) by computing $\partial\text{ELBO}/\partial\boldsymbol{\mu}$ analytically and comparing with finite differences.

---

13. Why This Matters for AI (2026 Perspective)

---

14. Conceptual Bridge

Looking Back

This section builds on a foundation laid across the previous six sections of Advanced Linear Algebra. From eigenvalues and eigenvectors (01), we borrow the spectral theorem - the fact that symmetric matrices diagonalize orthogonally - which provides the cleanest proof that positive definiteness is equivalent to positive eigenvalues. From SVD (02), we inherit the language of singular values and the connection between the spectral norm, condition number, and the stability of Cholesky. The orthogonality theory from 05 underpins the proof that $Q\Lambda Q^\top$ is symmetric (the spectral theorem) and that the PSD square root is well-defined and unique. Matrix norms from 06 give precise bounds on how a PD matrix behaves numerically: the condition number $\kappa(A) = \lambda_{\max}/\lambda_{\min}$ measures how far $A$ is from singular, and ill-conditioned PD matrices lead to numerical difficulties in Cholesky.

Most fundamentally, positive definite matrices are the matrices that have all four of the desirable properties we want from a symmetric matrix: uniquely solvable linear systems (invertible), well-defined quadratic energy (positive), computable square root (Cholesky), and monotone response to matrix operations (Loewner order). Understanding this quadrant of symmetry \times positivity is the analytical key to all probabilistic machine learning.

Looking Forward

Positive definite matrices are the entry point to the next major development: 08 (Matrix Decompositions) uses the Cholesky factorization developed here as one member of a family that includes LU (general matrices) and QR (orthogonal factorization). The Cholesky decomposition studied here in full generality is previewed (only briefly) in 08 - by the time students reach 08, the Cholesky algorithm should feel familiar.

The deeper application of positive definiteness comes later in the curriculum. In Chapter 6 (Probability Theory), covariance matrices $\Sigma \succeq 0$ appear in every multivariate distribution. In Chapter 8 (Optimization), the second-order conditions for minima require $\nabla^2\mathcal{L} \succeq 0$, and convexity of a function is equivalent to its Hessian being PSD everywhere. In Chapter 12 (Functional Analysis), Mercer's theorem and RKHS theory extend PSD kernels to infinite-dimensional inner product spaces. The PSD cone and semidefinite programming previewed in 9 appear extensively in Chapter 8 (Optimization) and in the robotics/control applications in later chapters.

Curriculum Position

POSITIVE DEFINITE MATRICES IN THE CURRICULUM
========================================================================

  Chapter 2: Linear Algebra Basics
  +-- 01 Vectors and Spaces        -> inner products, quadratic forms
  +-- 04 Determinants              -> Sylvester's criterion
  +-- 06 Vector Spaces             -> subspaces, null space
  
  Chapter 3: Advanced Linear Algebra
  +-- 01 Eigenvalues         -> spectral characterization of PD
  +-- 02 SVD                 -> singular values, condition number
  +-- 05 Orthogonality       -> spectral theorem, square root
  +-- 06 Matrix Norms        -> condition number, jitter scale
  +-- 07 Positive Definite   <== YOU ARE HERE
  |     |   Cholesky, LDL^T, Schur, log-det, Gram, SDP
  +-- 08 Matrix Decompositions  -> LU, QR, brief Cholesky recap
  
  Downstream:
  +-- Ch. 6 Probability: Gaussian covariance \Sigma \succ 0
  +-- Ch. 8 Optimization: Hessian \succ 0 at local minima; SDP
  +-- Ch. 12 Functional Analysis: Mercer, RKHS
  +-- Ch. 13 ML Math: K-FAC, VAE, GP, normalizing flows

========================================================================

The central insight to carry forward: positive definiteness is not a restriction but a characterization. A matrix that fails to be PD is not "almost OK" - it is categorically different, with no unique minimum, no valid probability interpretation, no Cholesky factorization. The condition $A \succ 0$ is the boundary between well-posedness and ill-posedness in a wide range of mathematical problems, and recognizing that boundary - and the tools to work with it (Cholesky, Schur, log-det) - is an essential skill for the working ML practitioner.

---

Appendix A: Closure Properties and Algebraic Operations

A.1 Preserving Positive Definiteness Under Operations

Understanding which operations preserve PD/PSD structure is practically important - it tells you when a composed matrix is guaranteed to be PD without recomputing.

Direct results:

The Schur product theorem. The Hadamard (element-wise) product of two PSD matrices is PSD. This is non-obvious! It is used in kernel methods: if $k_1$ and $k_2$ are valid kernels, their pointwise product $k_1(\mathbf{x},\mathbf{z}) \cdot k_2(\mathbf{x},\mathbf{z})$ is also a valid kernel.

Proof of Schur product theorem. If $A, B \succeq 0$, write $A = \sum_i \mathbf{u}_i\mathbf{u}_i^\top$ and $B = \sum_j \mathbf{v}_j\mathbf{v}_j^\top$ (Gram representations). Then $(A \odot B)_{kl} = A_{kl}B_{kl}$, which can be expressed as the Gram matrix of vectors $\mathbf{u}_i \otimes \mathbf{v}_j$ (Kronecker products of the generating vectors). Gram matrices are PSD. $\square$

Why $AB$ is not always PSD. Consider $A = \begin{pmatrix}1 & 0 \\ 0 & 2\end{pmatrix}$ and $B = \begin{pmatrix}2 & 1 \\ 1 & 1\end{pmatrix}$. Both are PD, but $AB = \begin{pmatrix}2 & 1 \\ 2 & 2\end{pmatrix}$, which is not symmetric (hence not PSD). The product of symmetric matrices is symmetric iff they commute.

For AI: Kronecker product structure is central to K-FAC: $F \approx A_l \otimes G_l$ is PSD since $A_l, G_l \succeq 0$. The Schur product theorem ensures that combining PD kernels (feature interactions) produces valid kernels.

A.2 Principal Submatrix Inheritance

Theorem A.1. If $A \succ 0$, every principal submatrix is $\succ 0$. If $A \succeq 0$, every principal submatrix is $\succeq 0$.

A principal submatrix of $A$ is obtained by deleting rows and columns with the same index set: $A_S = A[S, S]$ for $S \subseteq \{1,\ldots,n\}$.

Proof. For any $\mathbf{y} \in \mathbb{R}^{|S|}$, embed in $\mathbb{R}^n$ via $\mathbf{x}_i = y_i$ for $i \in S$ and $\mathbf{x}_i = 0$ for $i \notin S$. Then $\mathbf{y}^\top A_S \mathbf{y} = \mathbf{x}^\top A \mathbf{x} > 0$ for $\mathbf{y} \neq 0$ (since $\mathbf{x} \neq 0$). $\square$

Consequence. Every diagonal entry $A_{ii} > 0$ for $A \succ 0$ (take $S = \{i\}$). Every $2 \times 2$ principal submatrix is PD, so $A_{ii}A_{jj} > A_{ij}^2$ for all $i \neq j$ (Cauchy-Schwarz-type bound for matrix entries).

A.3 Indefinite Matrices and Saddle Points

For completeness: an indefinite symmetric matrix $A$ has both positive and negative eigenvalues. The corresponding quadratic form takes positive values in some directions and negative values in others. The level sets are hyperboloids.

Connection to saddle points in optimization. A critical point $\nabla\mathcal{L}(\boldsymbol{\theta}^*) = 0$ with indefinite Hessian $\nabla^2\mathcal{L}(\boldsymbol{\theta}^*)$ is a saddle point - not a local minimum. In neural network optimization:

• Gradient descent near a saddle point may slow down dramatically (gradient is small, Hessian is indefinite)
• The negative eigenvalue directions are "escape directions" - following them reduces the loss
• Stochastic gradient descent (SGD) with noise "escapes" saddles; deterministic gradient descent can get stuck

Sylvester's law of inertia. For a symmetric matrix $A$, the signature $(p, q, r)$ (number of positive, negative, zero eigenvalues) is invariant under congruence transformations $A \mapsto B^\top A B$ (for invertible $B$). This means no invertible change of basis can turn a saddle into a minimum - the number of descent directions is a geometric invariant of the quadratic form.

---

Appendix B: Numerical Methods and Implementation

B.1 Cholesky for Linear System Solving

A primary application of Cholesky is solving $A\mathbf{x} = \mathbf{b}$ for PD $A$.

Algorithm:

1. Factor: $A = LL^\top$ (Cholesky, $O(n^3/3)$)
2. Forward solve: $L\mathbf{y} = \mathbf{b}$ (substitute from top, $O(n^2)$)
3. Backward solve: $L^\top\mathbf{x} = \mathbf{y}$ (substitute from bottom, $O(n^2)$)

For multiple right-hand sides $B = [\mathbf{b}_1|\cdots|\mathbf{b}_k]$: factor once, solve $k$ times at cost $O(n^3/3 + kn^2)$ total.

Numerical example: Solve $\begin{pmatrix}4 & 2 \\ 2 & 3\end{pmatrix}\mathbf{x} = \begin{pmatrix}8 \\ 7\end{pmatrix}$.

Cholesky: $L = \begin{pmatrix}2 & 0 \\ 1 & \sqrt{2}\end{pmatrix}$ (since $L_{11}=2$, $L_{21}=1$, $L_{22}=\sqrt{3-1}=\sqrt{2}$).

Forward: $2y_1 = 8 \Rightarrow y_1 = 4$; $y_1 + \sqrt{2}y_2 = 7 \Rightarrow y_2 = (7-4)/\sqrt{2} = 3/\sqrt{2}$.

Backward: $\sqrt{2}x_2 = 3/\sqrt{2} \Rightarrow x_2 = 3/2$; $2x_1 + x_2 = 4 \Rightarrow x_1 = (4 - 3/2)/2 = 5/4$.

Verify: $\begin{pmatrix}4&2\\2&3\end{pmatrix}\begin{pmatrix}5/4\\3/2\end{pmatrix} = \begin{pmatrix}5+3\\5/2+9/2\end{pmatrix} = \begin{pmatrix}8\\7\end{pmatrix}$. OK

B.2 Cholesky Update and Downdate

When $A$ changes by a rank-1 update $A \leftarrow A + \mathbf{v}\mathbf{v}^\top$ (or downdate $A \leftarrow A - \mathbf{v}\mathbf{v}^\top$), it is wasteful to recompute the full Cholesky factorization. Rank-1 Cholesky update algorithms (LINPACK's dchud) update $L$ in $O(n^2)$ time.

Algorithm (rank-1 update, Gill-Golub-Murray-Saunders style):

Given $A = LL^\top$, compute $L'$ such that $L'L'^\top = A + \mathbf{v}\mathbf{v}^\top$:

x = L \ v           (forward solve: O(n^2))
for k = 1 to n:
    r = sqrt(L[k,k]^2 + x[k]^2)
    c = r / L[k,k]
    s = x[k] / L[k,k]
    L[k,k] = r
    L[k+1:n, k] = (L[k+1:n,k] + s*x[k+1:n]) / c
    x[k+1:n]    = c*x[k+1:n] - s*L[k+1:n,k]

This uses Givens rotations (-> 05: Orthogonality) to update $L$ in place.

For AI: Online learning algorithms that add one data point at a time use rank-1 Cholesky updates to maintain the Cholesky factorization of the Gram matrix. Sequential Bayesian updating (online GP regression) uses rank-1 Cholesky updates at cost $O(n^2)$ per update instead of $O(n^3)$ full refactorization.

B.3 Condition Number and Numerical Stability of Cholesky

Condition number. For $A \succ 0$, the condition number $\kappa(A) = \lambda_{\max}/\lambda_{\min}$. The relative error in solving $A\mathbf{x} = \mathbf{b}$ satisfies:

$$\frac{\|\hat{\mathbf{x}} - \mathbf{x}\|}{\|\mathbf{x}\|} \leq \kappa(A) \cdot \epsilon_{\text{mach}}.$$

For Cholesky specifically: the backward error bound on the Cholesky factor is $\|A - \hat{L}\hat{L}^\top\| \leq c_n\epsilon_{\text{mach}}\|A\|$. This means Cholesky is unconditionally backward stable for PD matrices. In contrast, LU without pivoting can have backward errors proportional to $\|L\|\|U\|$, which can be exponentially large.

Ill-conditioning warning. For large $\kappa(A)$:

• The log-det computation $2\sum\log L_{ii}$ involves some very small $L_{ii}$ (near-zero for near-singular dimensions) - those log terms dominate and may lose digits
• Linear solves $A^{-1}\mathbf{b}$ lose approximately $\log_{10}\kappa(A)$ digits of accuracy
• Cholesky itself remains stable; the error comes from the ill-conditioning of the problem, not the algorithm

Practical rule: If $\kappa(A) \approx 10^k$ and your floating-point arithmetic has $\approx 15$ decimal digits (double precision), you lose $k$ digits in the solution. For $k \geq 12$, the solution has no reliable digits.

B.4 Sparse Cholesky

For sparse PD matrices (many zero entries, arising in graph-structured models, FEM discretizations, and spatial statistics), general Cholesky is wasteful because it may "fill in" zeros and produce dense $L$.

The fill-in problem. Even if $A$ is sparse, $L$ may be dense. The fill-in pattern depends on the sparsity structure of $A$ and the ordering of variables.

Reordering algorithms. The approximate minimum degree (AMD) and nested dissection algorithms reorder the variables (permute $A \leftarrow PAP^\top$ for a permutation matrix $P$) to minimize fill-in. The standard library for sparse Cholesky is SuiteSparse (CHOLMOD).

For AI: In Gaussian Markov random fields (GMRFs), the precision matrix (inverse covariance) is sparse. Sparse Cholesky enables efficient inference in spatial models, graph neural networks with Gaussian priors, and structured prediction. GPyTorch's KeOps and PyTorch Sparse backends exploit sparsity for large-scale GP approximations.

---

Appendix C: Positive Definiteness in Complex Vector Spaces

C.1 Hermitian Positive Definite Matrices

Over $\mathbb{C}$, positive definiteness generalizes to Hermitian positive definite (HPD) matrices.

Definition. $A \in \mathbb{C}^{n \times n}$ is Hermitian positive definite if:

• $A = A^*$ (where $A^* = \bar{A}^\top$ is the conjugate transpose)
• $\mathbf{x}^* A \mathbf{x} > 0$ for all $\mathbf{x} \in \mathbb{C}^n \setminus \{0\}$

(Note: $\mathbf{x}^* A \mathbf{x}$ is always real for Hermitian $A$, so the positivity condition makes sense.)

Cholesky for HPD. Every HPD matrix $A$ has a unique Cholesky factorization $A = LL^*$ where $L$ is lower triangular with positive real diagonal entries.

For AI: Complex-valued neural networks (used in signal processing, quantum computing simulation) use HPD covariance matrices. Radar signal processing, MRI reconstruction, and quantum chemistry all work in $\mathbb{C}^n$ with HPD matrices.

---

Appendix D: Advanced Theory

D.1 The Polar Decomposition and Positive Definite Factors

Every invertible matrix $A \in \mathbb{R}^{n \times n}$ has a unique polar decomposition:

$$A = QP$$

where $Q$ is orthogonal ($Q^\top Q = I$) and $P = \sqrt{A^\top A} \succ 0$ is the unique PD symmetric positive definite "right polar factor." This is the matrix analogue of the polar form $z = e^{i\theta}|z|$ for complex numbers.

Existence and uniqueness. Using SVD $A = U\Sigma V^\top$:

$$A = U\Sigma V^\top = (UV^\top)(V\Sigma V^\top) = QP$$

where $Q = UV^\top$ (orthogonal) and $P = V\Sigma V^\top \succ 0$ (PSD square root of $A^\top A$).

Uniqueness of $P$: $P = (A^\top A)^{1/2}$, which is the unique PSD square root (Theorem 5.1). Given $P$, $Q = AP^{-1}$ is uniquely determined.

For AI: The polar decomposition appears in:

• Weight matrix orthogonalization: in Muon optimizer (Kosson et al., 2024), gradient updates are "orthogonalized" using the polar factor $Q = AP^{-1}$, which projects the weight update onto the manifold of orthogonal matrices, maintaining stable singular value spectrum during training
• Group equivariance: equivariant neural networks (E(n)-equivariant GNNs) use orthogonal matrices as symmetry group elements; the polar decomposition provides a differentiable projection onto $O(n)$

Computing the polar factor. Newton iteration: $X_{k+1} = \frac{1}{2}(X_k + X_k^{-\top})$ converges to $Q$ starting from $X_0 = A$. Quadratic convergence once close to $Q$. Each step costs one matrix inversion.

D.2 Positive Definite Functions and Bochner's Theorem

A function $k: \mathbb{R}^d \to \mathbb{R}$ is a positive definite function if for every $n$ and every set of points $\mathbf{x}_1,\ldots,\mathbf{x}_n \in \mathbb{R}^d$, the matrix $K_{ij} = k(\mathbf{x}_i - \mathbf{x}_j)$ is PSD.

Bochner's theorem. A continuous stationary kernel $k(\mathbf{x}-\mathbf{z})$ is a positive definite function on $\mathbb{R}^d$ if and only if it is the Fourier transform of a non-negative measure $\mu$ (spectral density):

$$k(\boldsymbol{\tau}) = \int_{\mathbb{R}^d} e^{i\boldsymbol{\omega}^\top\boldsymbol{\tau}} d\mu(\boldsymbol{\omega}).$$

Implications for kernel design:

• RBF kernel: $k(\boldsymbol{\tau}) = \exp(-\|\boldsymbol{\tau}\|^2/2\ell^2)$ is PD (its Fourier transform is a Gaussian, which is non-negative)
• Matern kernels: PD for all valid smoothness parameters $\nu > 0$
• Checking validity: A kernel is valid iff its Fourier transform is non-negative - this is the ultimate test

Random Fourier features (Rahimi & Recht, 2007). Bochner's theorem gives a randomized approximation to stationary PD kernels:

$$k(\mathbf{x}-\mathbf{z}) \approx z(\mathbf{x})^\top z(\mathbf{z}), \quad z(\mathbf{x}) = \sqrt{2/D}[\cos(\boldsymbol{\omega}_1^\top\mathbf{x}+b_1),\ldots,\cos(\boldsymbol{\omega}_D^\top\mathbf{x}+b_D)]$$

where $\boldsymbol{\omega}_i \sim \mu$ (the spectral distribution) and $b_i \sim \text{Uniform}[0,2\pi]$. This approximates the PD kernel as a finite inner product, enabling scalable kernel methods without forming the full $n \times n$ kernel matrix. Used in Performer transformer attention (FAVOR+, Choromanski et al., 2020).

D.3 The Determinant and Volume

Geometric interpretation of $\det A$ for $A \succ 0$.

For $A \succ 0$, the quadratic form $\mathbf{x}^\top A^{-1} \mathbf{x} \leq 1$ defines an ellipsoid $\mathcal{E}_A$ in $\mathbb{R}^n$. Its volume is:

$$\text{Vol}(\mathcal{E}_A) = \frac{\pi^{n/2}}{\Gamma(n/2+1)} \cdot (\det A)^{1/2}.$$

So $\det A$ measures the "volume" of the ellipsoid associated with $A$. Larger $\det A$ means a more "spread out" distribution.

In the Gaussian context: The normalizing constant of $\mathcal{N}(\boldsymbol{\mu}, \Sigma)$ is $(2\pi)^{n/2}(\det\Sigma)^{1/2}$. This is the volume of the "effective support" ellipsoid scaled by $(2\pi)^{n/2}$.

D-optimal experimental design. In Bayesian experimental design, the D-optimal criterion selects experiments $\mathbf{x}_1,\ldots,\mathbf{x}_n$ to maximize $\det(X^\top X + \sigma^{-2}I)$ - the determinant of the posterior precision matrix, which equals $1/\det(\Sigma_{\text{post}})$. Maximizing $\det(\text{posterior precision})$ minimizes the volume of the posterior confidence ellipsoid. This is an SDP with $X$ as the optimization variable.

Maximum entropy Gaussians. Among all distributions with mean $\boldsymbol{\mu}$ and covariance $\Sigma$, the maximum entropy distribution is $\mathcal{N}(\boldsymbol{\mu}, \Sigma)$, with entropy $H = \frac{n}{2}(1+\log(2\pi)) + \frac{1}{2}\log\det\Sigma$. Maximizing entropy subject to a trace constraint $\text{tr}(\Sigma) \leq c$ gives $\Sigma = (c/n)I$ (isotropic Gaussian) - the most "uninformed" Gaussian with bounded total variance.

D.4 Matrix Exponential and the PD Manifold

The set of PD matrices $\mathbb{S}_{++}^n$ is an open smooth manifold (in fact, a Riemannian manifold with the affine-invariant metric). This means concepts from differential geometry apply.

The matrix exponential. For any symmetric $S \in \mathbb{S}^n$, $e^S = Q e^\Lambda Q^\top \in \mathbb{S}_{++}^n$ (positive definite, since exponentials of real numbers are positive). The map $\exp: \mathbb{S}^n \to \mathbb{S}_{++}^n$ is a bijection - PD matrices are exactly exponentials of symmetric matrices.

Log-Euclidean metric. A simple Riemannian metric on $\mathbb{S}_{++}^n$ (used in diffusion tensor imaging) defines the "distance" between $A, B \succ 0$ as:

$$d_{\log}(A, B) = \|\log A - \log B\|_F.$$

This treats PD matrices as if they lived in a flat space via the log map. Computations (means, geodesics, interpolations) become Euclidean operations on $\log A$ and $\log B$.

Affine-invariant metric (Riemannian). The more geometrically natural metric:

$$d_{\text{AI}}(A, B) = \|\log(A^{-1/2}B A^{-1/2})\|_F = \left(\sum_{i=1}^n \log^2\lambda_i(A^{-1}B)\right)^{1/2}$$

is invariant under congruence $A \mapsto C^\top A C$ (affine-invariant), making it suitable for geometric statistics.

For AI: The SPD (symmetric positive definite) manifold appears in:

• Diffusion tensor MRI: each voxel has a $3 \times 3$ diffusion tensor $\in \mathbb{S}_{++}^3$; Riemannian means avoid non-PD artifacts
• Covariance-based classifiers: classifying neural activation covariance matrices using Riemannian geometry (SPDNet, 2017)
• Transformer positional encodings with structured covariances: attention with SPD constraints on the key-query metric tensors

D.5 Inequalities for Positive Definite Matrices

Several fundamental inequalities apply specifically to PD matrices:

Hadamard's inequality. For $A \succ 0$:

$$\det A \leq \prod_{i=1}^n A_{ii}.$$

Equality iff $A$ is diagonal. Interpretation: the volume of the ellipsoid defined by $A$ is at most the product of the diagonal entries (the volumes along coordinate axes).

Minkowski's inequality. For $A, B \succ 0$:

$$(\det(A+B))^{1/n} \geq (\det A)^{1/n} + (\det B)^{1/n}.$$

This is the matrix analogue of $\sqrt[n]{a+b} \geq \sqrt[n]{a} + \sqrt[n]{b}$ for scalars (Minkowski). It follows from the concavity of $A \mapsto (\det A)^{1/n}$ on $\mathbb{S}_{++}^n$.

Fan-Ky inequality. For $A \succ 0$ and $k \leq n$:

$$\sum_{i=1}^k \log\lambda_i(A+B) \geq \sum_{i=1}^k \log\lambda_i(A) + \sum_{i=1}^k \log\lambda_i(B).$$

(The sum of the $k$ largest log-eigenvalues of $A+B$ is at least the sum from $A$ plus the sum from $B$.)

Anderson's inequality. For $A, B, C \succ 0$ with $C = A + B$:

$$C^{-1} \preceq A^{-1} + B^{-1}.$$

Interpretation: combining two experiments (adding their information matrices) gives more information than either alone, but the combined precision matrix is "sandwiched" by the sum of individual inverse covariances.

For AI: Hadamard's inequality provides a bound on the log-det: $\log\det A \leq \sum_i \log A_{ii} = \log\prod A_{ii}$, useful when diagonal entries are cheaply computed. Minkowski's inequality is used in information-geometric proofs about combining Gaussian priors in Bayesian learning.

---

Appendix E: Extended Examples and Case Studies

E.1 Case Study: Covariance Estimation in High Dimensions

One of the most practically important applications of PSD matrix theory is estimating covariance matrices from data when the number of features $p$ is comparable to or larger than the number of samples $n$.

The problem. Given $n$ samples $\mathbf{x}_1,\ldots,\mathbf{x}_n \in \mathbb{R}^p$ (zero mean for simplicity), the sample covariance is:

$$\hat{S} = \frac{1}{n}\sum_{i=1}^n \mathbf{x}_i\mathbf{x}_i^\top = \frac{1}{n}X^\top X.$$

$\hat{S} \succeq 0$ always (it's a Gram matrix). But:

• If $n < p$: $\text{rank}(\hat{S}) \leq n < p$, so $\hat{S}$ is singular (not PD)
• Even for $n \geq p$: $\hat{S}$ can be ill-conditioned (extreme ratio of largest to smallest eigenvalue)

The Marchenko-Pastur law. For $n, p \to \infty$ with $n/p \to \gamma > 1$ and $\mathbf{x}_i \sim \mathcal{N}(0, I_p)$, the empirical eigenvalue distribution of $\hat{S}$ converges to the Marchenko-Pastur distribution:

$$f(\lambda) = \frac{\sqrt{(\lambda_+ - \lambda)(\lambda - \lambda_-)}}{2\pi\gamma\lambda}$$

for $\lambda \in [\lambda_-, \lambda_+]$, where $\lambda_\pm = (1 \pm 1/\sqrt{\gamma})^2$. This means even under the identity covariance, sample eigenvalues are spread over a wide range - the smallest eigenvalues of $\hat{S}$ are much smaller than 1 and the largest are much larger than 1.

Ledoit-Wolf shrinkage. The Ledoit-Wolf estimator regularizes the sample covariance toward a multiple of identity:

$$\hat{\Sigma}_{\text{LW}} = \alpha\hat{S} + (1-\alpha)\mu I$$

where $\mu = \text{tr}(\hat{S})/p$ (the average eigenvalue) and $\alpha$ is chosen to minimize mean squared error under the Frobenius norm. The result is always PD (since $\mu I \succ 0$ and $\hat{S} \succeq 0$).

Graphical lasso. Assuming the precision matrix (inverse covariance) $\Omega = \Sigma^{-1}$ is sparse (many zeros encode conditional independence), the graphical lasso estimates:

$$\hat{\Omega} = \arg\min_{\Omega \succ 0} \left[\text{tr}(\hat{S}\Omega) - \log\det\Omega + \lambda\|\Omega\|_1\right]$$

where $\|\cdot\|_1$ is the element-wise $\ell_1$ norm. The constraint $\Omega \succ 0$ is the PD cone constraint. The objective $\text{tr}(\hat{S}\Omega) - \log\det\Omega$ is the negative Gaussian log-likelihood (up to constants). Optimized by block coordinate descent (glasso algorithm) using the Schur complement for each block update.

For LLMs: Large language models trained with structured sparse attention patterns (local + global attention, as in Longformer and BigBird) can be understood as learning sparse precision matrices of the attention distributions.

E.2 Case Study: Gaussian Process Regression, Step by Step

Let's trace through a complete GP regression computation to see all the PD machinery in action.

Setup. Training data: $(x_1, y_1) = (0, 0.5)$, $(x_2, y_2) = (1, 1.2)$, $(x_3, y_3) = (2, 0.8)$. Test point: $x_* = 1.5$. Kernel: $k(x,z) = \exp(-(x-z)^2/2)$ (RBF with length-scale 1). Noise: $\sigma^2 = 0.1$.

Step 1: Compute the kernel matrix.

$$K_{ij} = \exp(-(x_i - x_j)^2/2).$$ $$K = \begin{pmatrix}1 & e^{-0.5} & e^{-2} \\ e^{-0.5} & 1 & e^{-0.5} \\ e^{-2} & e^{-0.5} & 1\end{pmatrix} \approx \begin{pmatrix}1 & 0.607 & 0.135 \\ 0.607 & 1 & 0.607 \\ 0.135 & 0.607 & 1\end{pmatrix}.$$

Step 2: Add noise and Cholesky factor.

$$K + \sigma^2 I = K + 0.1I \approx \begin{pmatrix}1.1 & 0.607 & 0.135 \\ 0.607 & 1.1 & 0.607 \\ 0.135 & 0.607 & 1.1\end{pmatrix}.$$

Compute $L = \text{chol}(K + 0.1I)$. (Details in theory.ipynb.)

Step 3: Predictive mean.

$\mathbf{k}_* = [k(x_*, x_i)]_i = [\exp(-2.25/2), \exp(-0.25/2), \exp(-0.25/2)] \approx [0.325, 0.882, 0.882]$.

$$\boldsymbol{\alpha} = (K + 0.1I)^{-1}\mathbf{y} = L^{-\top}L^{-1}\mathbf{y} \quad (\text{via Cholesky solve})$$ $$\mu_* = \mathbf{k}_*^\top\boldsymbol{\alpha}.$$

Step 4: Predictive variance.

$$\sigma_*^2 = k(x_*, x_*) - \mathbf{k}_*^\top(K+0.1I)^{-1}\mathbf{k}_* = 1 - \mathbf{v}^\top\mathbf{v}, \quad \mathbf{v} = L^{-1}\mathbf{k}_*.$$

The variance is the Schur complement of $(K + 0.1I)$ evaluated at $x_*$.

Interpretation. The Cholesky factorization appears in every step: solving for $\boldsymbol{\alpha}$, computing the predictive variance, and computing the log-marginal likelihood for hyperparameter optimization. The PD condition on $K + 0.1I$ guarantees $\sigma_*^2 \geq 0$ (the Schur complement is non-negative).

E.3 Case Study: The Multivariate Normal in Deep Learning

Parameterizing structured covariances. In deep generative models (VAEs, normalizing flows, diffusion models), we frequently need to parameterize a multivariate Gaussian $\mathcal{N}(\boldsymbol{\mu}, \Sigma)$ with $\Sigma$ learnable. The challenge: $\Sigma$ must be PD, and the parameterization must be differentiable.

Approach 1: Diagonal. $\Sigma = \text{diag}(\sigma_1^2,\ldots,\sigma_d^2)$ with $\sigma_i = \text{softplus}(s_i)$. Free parameters: $d$. Correlation structure: none. Used in standard VAE (Kingma & Welling, 2013).

Approach 2: Lower triangular Cholesky. $\Sigma = LL^\top$ where $L$ is lower triangular:

• Diagonal: $L_{ii} = \text{softplus}(\tilde{L}_{ii})$ ensures positivity
• Off-diagonal: $L_{ij}$ for $i > j$ unconstrained (any real number) Free parameters: $d(d+1)/2$. Full correlation structure. Used in full-covariance VAEs, normalizing flows.

Approach 3: Low-rank + diagonal. $\Sigma = FF^\top + D$ where $F \in \mathbb{R}^{d \times r}$ ($r \ll d$) and $D = \text{diag}(d_1,\ldots,d_d)$ with $d_i > 0$. Free parameters: $dr + d$. Captures $r$ "directions of correlation." Used in structured VBs, mean-field approximations.

The KL divergence between Gaussians. For $q = \mathcal{N}(\boldsymbol{\mu}_q, \Sigma_q)$ and $p = \mathcal{N}(\boldsymbol{\mu}_p, \Sigma_p)$:

$$D_{\text{KL}}(q\|p) = \frac{1}{2}\left[\text{tr}(\Sigma_p^{-1}\Sigma_q) + (\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)^\top\Sigma_p^{-1}(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q) - d + \log\frac{\det\Sigma_p}{\det\Sigma_q}\right].$$

For $p = \mathcal{N}(\mathbf{0}, I)$ (standard VAE prior):

$$D_{\text{KL}}(q\|\mathcal{N}(\mathbf{0},I)) = \frac{1}{2}\left[\text{tr}(\Sigma_q) + \boldsymbol{\mu}_q^\top\boldsymbol{\mu}_q - d - \log\det\Sigma_q\right].$$

For diagonal $\Sigma_q = \text{diag}(\sigma_i^2)$:

$$D_{\text{KL}} = \frac{1}{2}\sum_{i=1}^d \left[\sigma_i^2 + \mu_{q,i}^2 - 1 - \log\sigma_i^2\right].$$

For Cholesky $\Sigma_q = LL^\top$: $\log\det\Sigma_q = 2\sum_i \log L_{ii}$ and $\text{tr}(\Sigma_q) = \|L\|_F^2$.

All these computations - $\log\det\Sigma_q$, $\text{tr}(\Sigma_q)$, $\Sigma_q^{-1}$ - exploit the Cholesky and log-det tools developed in 4 and 7.

E.4 The Normal Equations and Ridge Regression

Least squares. Given $X \in \mathbb{R}^{n \times p}$ (data) and $\mathbf{y} \in \mathbb{R}^n$ (targets), ordinary least squares (OLS) minimizes:

$$\hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \|X\boldsymbol{\beta} - \mathbf{y}\|^2.$$

The normal equations: $X^\top X \hat{\boldsymbol{\beta}} = X^\top\mathbf{y}$. The matrix $X^\top X \succeq 0$ is always PSD (Gram matrix). If $X$ has full column rank ($n \geq p$, rank $p$), then $X^\top X \succ 0$ and the unique solution is $\hat{\boldsymbol{\beta}} = (X^\top X)^{-1}X^\top\mathbf{y}$.

Ridge regression. Add $\ell_2$ regularization: minimize $\|X\boldsymbol{\beta}-\mathbf{y}\|^2 + \lambda\|\boldsymbol{\beta}\|^2$.

Normal equations: $(X^\top X + \lambda I)\hat{\boldsymbol{\beta}} = X^\top\mathbf{y}$.

The matrix $A = X^\top X + \lambda I \succ 0$ for any $\lambda > 0$ (by Proposition 2.3.7). This is why ridge regression is always numerically stable: the PD matrix $A$ has all positive pivots, and Cholesky never fails.

Condition number improvement. $\kappa(X^\top X + \lambda I) = (\sigma_1^2 + \lambda)/(\sigma_p^2 + \lambda)$ where $\sigma_1 \geq \cdots \geq \sigma_p \geq 0$ are the singular values of $X$. For small $\lambda$, $\kappa \approx \sigma_1^2/\sigma_p^2 = \kappa(X)^2$. As $\lambda \to \infty$, $\kappa \to 1$ (well-conditioned). Choosing $\lambda$ balances bias (regularization) vs numerical stability.

Kernel ridge regression. For nonlinear regression using a PD kernel $k$, form the kernel matrix $K = k(\mathbf{X},\mathbf{X})$ and solve $(K + \lambda I)\boldsymbol{\alpha} = \mathbf{y}$. Predictions: $\hat{f}(\mathbf{x}_*) = \mathbf{k}_*^\top\boldsymbol{\alpha}$ where $\mathbf{k}_* = k(\mathbf{x}_*, \mathbf{X})$. This is equivalent to GP regression with the kernel $k$ and noise variance $\lambda$. The PD condition on $K + \lambda I$ is exactly the Cholesky solvability condition.

---

Appendix F: Quick Reference

F.1 Four Equivalent Characterizations of $A \succ 0$

POSITIVE DEFINITENESS: FOUR EQUIVALENT TESTS
========================================================================

  Let A \in \mathbb{R}^n^x^n be symmetric. The following are equivalent:

  DEFINITION   \forallx\neq0: x^TAx > 0            (quadratic form positive)
  SPECTRAL     \foralli: \lambda^i(A) > 0              (all eigenvalues positive)
  SYLVESTER    \forallk: \Delta_k = det(A[1:k,1:k]) > 0 (leading minors positive)
  CHOLESKY     \exists! lower triangular L with positive diagonal: A = LL^T

  COST:
  Definition  O(n) per test vector            <- manual/symbolic
  Spectral    O(n^3) eigendecomposition        <- most informative
  Sylvester   O(n^4) for all n determinants    <- manual/symbolic
  Cholesky    O(n^3/3) factorization           <- fastest in practice

========================================================================

F.2 Key Formulas

F.3 Notation Summary

Following the NOTATION_GUIDE.md conventions:

F.4 Python Cheatsheet

import numpy as np
import scipy.linalg as la

# === Cholesky factorization ===
L = np.linalg.cholesky(A)          # A = L @ L.T
# Or via scipy (can also compute upper factor):
L = la.cholesky(A, lower=True)

# === LDL^T factorization ===
L_ldl, D_ldl, _ = la.ldl(A)       # A = L @ D @ L.T

# === Log-determinant (numerically stable) ===
L = np.linalg.cholesky(A)
log_det = 2 * np.sum(np.log(np.diag(L)))

# === Solve A x = b via Cholesky ===
L = np.linalg.cholesky(A)
x = la.cho_solve((L, True), b)     # uses LAPACK dpotrs

# === PSD square root ===
eigvals, Q = np.linalg.eigh(A)     # symmetric eigendecomp
A_sqrt = Q @ np.diag(np.sqrt(eigvals)) @ Q.T

# === Woodbury: (A + U C V)^{-1} ===
# via np.linalg.solve on the smaller system

# === Reparameterization trick ===
L = np.linalg.cholesky(Sigma)      # Sigma = L @ L.T
eps = np.random.randn(n_samples, d)
z = mu + (L @ eps.T).T             # z ~ N(mu, Sigma)

# === Gram matrix ===
G = X @ X.T                        # G = X X^T, PSD

# === Check PD ===
try:
    np.linalg.cholesky(A)
    is_pd = True
except np.linalg.LinAlgError:
    is_pd = False

# === Check PSD ===
eigvals = np.linalg.eigvalsh(A)
is_psd = np.all(eigvals >= -1e-10)

---