USAAIO Lesson 23, from Week 8, fully worked. It builds the N(mu, Sigma) density term by term from the one-dimensional bell curve and proves from Cov(a^T x) why Sigma must be symmetric and PSD - and positive definite if it is to be inverted. It then covers the covariance ellipse from the eigendecomposition, derives and hand-computes the Mahalanobis distance, and covers the normalizing constant and the determinant. It derives the marginals and the bivariate conditional entry by entry, proves that a diagonal Sigma implies independence by factorizing the density, and gives the X^2 counterexample. It closes with Cholesky sampling, x = mu + Lz, proving Cov(mu + Lz) = Sigma line by line, verifying against 200k samples, and tying it to the VAE reparameterization trick. Every snippet runs standalone, and every number was produced by real execution. The lesson runs to 63 slides.
Subject: Machine Learning · 115 slides · code lesson
Open the interactive version of this deck · Homework for this lesson
Title
USAAIO · Lesson 23 · Week 8
The bell curve in many dimensions - and the distribution at the heart of the VAE. We build N(mu, Sigma) from the 1-D Gaussian term by term, prove why Sigma must be PSD, derive Mahalanobis distance, marginals, conditionals, and independence, then sample by Cholesky and prove x = mu + Lz has covariance Sigma.
Objectives
N(mu, Sigma) density from the 1-D bell curve and explain every factor - the quadratic form, the |Sigma|, and the (2pi) powerSigma must be symmetric PSD (and PD to be invertible) from Cov(a^T x) = a^T Sigma a >= 0x_2 | x_1 entry by entry, and prove diagonal Sigma => independent by factorizing the densityx = mu + Lz with L = chol(Sigma), prove Cov(x) = Sigma, recover it from 200k draws, and connect it to the VAE reparameterization trickWarm-up
Discussion prompt
Before we open Lesson 23: The Multivariate Gaussian: without looking back, what was the main idea of Affine Transforms, and what could you do by the end of it that you could not do before?
Hint: One sentence for the idea, one for the skill. If the second one is blank, that is the part to revisit.
Answer:
affine maps T(x)=Wx+b, why the composition of affine maps is affine, the resulting collapse of linear networks into a single layer, and why nonlinearity is essential — proven on XOR. Verify the collapse and watch a linear model fail XOR while an MLP solves it.
Section
Part 1 of 7 - the density
Concept
You already know the scalar bell curve: a random score x with mean mu and variance sigma^2 has density
\[ p(x) = \frac{1}{\sqrt{2\pi\,\sigma^2}}\,\exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \]
Two moving parts: a front constant that makes the area integrate to 1, and an exponent that is a squared, scaled distance from the center mu. Every piece of the multivariate version is one of these two, generalized.
Counterexample
Discussion prompt
You already know the scalar bell curve: a random score x with mean mu and variance sigma^2 has density
That is stated as though it always holds. Do one of two things: produce a case where it fails, or say precisely what rules such a case out. "It just does" is not on the menu.
Hint: Hunt at the extremes first — zero, one, negative, empty, equal. If every extreme survives, the reason they survive is the proof.
Concept
We carry ONE example through the whole lesson: a 2-D vector x = [x_1, x_2] - say a student's standardized reading and math scores - with center mu and covariance Sigma.
\[ \mu = \begin{bmatrix} 1 \\ -1 \end{bmatrix}, \qquad \Sigma = \begin{bmatrix} 2 & 0.8 \\ 0.8 & 1 \end{bmatrix} \]
| symbol | meaning | value |
|---|---|---|
| mu_1 | mean of x_1 | 1 |
| mu_2 | mean of x_2 | -1 |
| Sigma_11 | variance of x_1 | 2 |
| Sigma_22 | variance of x_2 | 1 |
| Sigma_12 = Sigma_21 | covariance of x_1, x_2 | 0.8 |
Comparison
Comparison matrix
From Our running example: two correlated features: refill the meaning column from what you know. The rest of the table is as it appeared.
| symbol | meaning | value |
|---|---|---|
| mu_1 | mean of x_1 | 1 |
| mu_2 | mean of x_2 | -1 |
| Sigma_11 | variance of x_1 | 2 |
| Sigma_22 | variance of x_2 | 1 |
| Sigma_12 = Sigma_21 | covariance of x_1, x_2 | 0.8 |
Concept
Sigma collects every pairwise covariance. Entry (i, j) is Cov(x_i, x_j) = E[(x_i - mu_i)(x_j - mu_j)]. The diagonal holds the variances; the off-diagonal holds the coupling.
\[ \Sigma = \mathbb{E}\big[(x-\mu)(x-\mu)^\top\big] = \begin{bmatrix} \operatorname{Var}(x_1) & \operatorname{Cov}(x_1,x_2) \\ \operatorname{Cov}(x_2,x_1) & \operatorname{Var}(x_2) \end{bmatrix} \]
Sigma is symmetric because Cov(x_i, x_j) = Cov(x_j, x_i)
Why: Swapping the two factors inside the expectation changes nothing, so Sigma_ij = Sigma_ji. Symmetry is not optional - it falls straight out of the definition, and we will lean on it constantly.
Analogy
Discussion prompt
Explain The covariance matrix, entry by entry by analogy to something with no Machine Learning in it at all — a queue, a recipe, a map, a bank balance, whatever fits. Then say where your analogy breaks.
Hint: An analogy that never breaks is not an analogy, it is the same idea wearing a hat. Find the seam — that is the part that is actually new.
Answer:
Sigma collects every pairwise covariance. Entry (i, j) is Cov(x_i, x_j) = E[(x_i - mu_i)(x_j - mu_j)]. The diagonal holds the variances; the off-diagonal holds the coupling.
Intuition
In 1-D the exponent divides by sigma^2 - it rescales the squared miss into variance units. A miss of 2 is a big deal if sigma^2 = 0.1 and trivial if sigma^2 = 100.
In many dimensions the single number 1/sigma^2 becomes the whole matrix Sigma^{-1}. It rescales the miss vector x - mu by variance and un-tilts the correlation, so the exponent measures distance the way the data actually spreads - not with a naive ruler.
Explain it
Discussion prompt
Explain What Sigma^{-1} is doing in the exponent to a student a year behind you. No notation, no jargon they have not met — and it still has to be true.
Hint: If your explanation needs a symbol they have never seen, you are describing the notation rather than the idea.
Answer:
In 1-D the exponent divides by sigma^2 - it rescales the squared miss into variance units. A miss of 2 is a big deal if sigma^2 = 0.1 and trivial if sigma^2 = 100.
Concept
Promote each 1-D piece to its matrix form: (x - mu)^2 / sigma^2 becomes the quadratic form (x - mu)^T Sigma^{-1} (x - mu), and sqrt(2 pi sigma^2) becomes (2 pi)^{d/2} |Sigma|^{1/2} for d dimensions.
\[ p(x) = \frac{1}{(2\pi)^{d/2}\,|\Sigma|^{1/2}}\,\exp\!\left(-\tfrac{1}{2}(x-\mu)^\top \Sigma^{-1}(x-\mu)\right) \]
| 1-D piece | generalizes to | role |
|---|---|---|
| (x - mu)^2 / sigma^2 | (x - mu)^T Sigma^{-1} (x - mu) | scaled squared distance |
| sqrt(2 pi) | (2 pi)^{d/2} | dimension count in the constant |
| sqrt(sigma^2) | |Sigma|^{1/2} | total spread (volume) |
Trade off
Comparison matrix
From The multivariate Gaussian density N(mu, Sigma): every row here is a choice with a cost. Fill the generalizes to column, then say which row you would actually pick and what you give up for it.
| 1-D piece | generalizes to | role |
|---|---|---|
| (x - mu)^2 / sigma^2 | (x - mu)^T Sigma^{-1} (x - mu) | scaled squared distance |
| sqrt(2 pi) | (2 pi)^{d/2} | dimension count in the constant |
| sqrt(sigma^2) | |Sigma|^{1/2} | total spread (volume) |
Concept
Look at what the formula demands: it uses Sigma^{-1} (so Sigma must be invertible) and |Sigma|^{1/2} (so |Sigma| must be positive, or the square root is imaginary).
It also needs the exponent to be <= 0 for every x, so the quadratic form (x - mu)^T Sigma^{-1} (x - mu) must never go negative. All three requirements point at one property of Sigma: positive definiteness. We prove that next.
Section
Part 2 of 7 - the constraint
Concept
Every covariance matrix is symmetric positive semidefinite (PSD): symmetric, and a^T Sigma a >= 0 for every vector a. For the density to be non-degenerate (invertible Sigma), we need positive definite (PD): a^T Sigma a > 0 for every a != 0.
\[ \Sigma = \Sigma^\top \quad\text{and}\quad a^\top \Sigma\, a \ge 0 \;\; \forall a \qquad (\text{PD: } > 0 \;\forall\, a \neq 0) \]
Ranking
Put in order
Put the moves of Prove a^T Sigma a >= 0, one line at a time into the order they have to happen.
Why: These are the moves of the worked example in the order it makes them, and each one is set up by the one before it. y is a single random number - the component of the centered vector along a.
Worked example
Pick any direction a and project the data onto it: define scalar y = a^T (x - mu)
Why: y is a single random number - the component of the centered vector along a. Its variance cannot be negative, and we will show that variance IS a^T Sigma a.
Write Var(y) = E[y^2] since E[y] = 0
Why: y is a^T times a centered vector, so E[y] = a^T E[x - mu] = 0. For a zero-mean quantity, variance is just the mean square.
\[ \operatorname{Var}(y) = \mathbb{E}[\,y^2\,] = \mathbb{E}\big[(a^\top(x-\mu))^2\big] \]
Expand the square as (a^T u)(u^T a) with u = x - mu
Why: A scalar equals itself transposed: a^T u = (a^T u)^T = u^T a. So y^2 = (a^T u)(u^T a) = a^T (u u^T) a.
\[ = \mathbb{E}\big[a^\top (x-\mu)(x-\mu)^\top a\big] \]
Pull the constant a outside the expectation
Why: a is not random, so E[a^T M a] = a^T E[M] a. The middle expectation is exactly Sigma by definition.
\[ = a^\top\,\mathbb{E}[(x-\mu)(x-\mu)^\top]\,a = a^\top \Sigma\, a \]
Conclude a^T Sigma a = Var(y) >= 0
Why: A real variance is never negative, and it equals a^T Sigma a for EVERY a. So Sigma is PSD. This is the whole reason covariance matrices are special - it is not an assumption, it is forced.
Notation
Annotate
From Prove a^T Sigma a >= 0, one line at a time — read this one piece at a time. What is each part doing?
On: \( \operatorname{Var}(y) = \mathbb{E}[\,y^2\,] = \mathbb{E}\big[(a^\top(x-\mu))^2\big] \)
Estimation
Predict first
PSD is guaranteed. For our density we need strictly PD - invertible. Confirm it three equivalent ways in one runnable block:
Commit before you compute: what does Our Sigma is PD: check the quadratic form and eigenvalues come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: Eigenvalues {0.5566, 2.4434} are both positive -> PD
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. A symmetric matrix is PD exactly when every eigenvalue is > 0.
Worked example
PSD is guaranteed. For our density we need strictly PD - invertible. Confirm it three equivalent ways in one runnable block:
import numpy as np
Sigma = np.array([[2., 0.8], [0.8, 1.]])
print(np.allclose(Sigma, Sigma.T)) # symmetric?
print(np.linalg.eigvalsh(Sigma).round(4)) # all eigenvalues > 0?
print(round(np.linalg.det(Sigma), 4)) # det > 0?
a = np.array([3., -1.])
print(round(a @ Sigma @ a, 4)) # a^T Sigma a > 0?Eigenvalues {0.5566, 2.4434} are both positive -> PD
Why: A symmetric matrix is PD exactly when every eigenvalue is > 0. Both are, det = 1.36 > 0, and the sample direction gives a^T Sigma a = 14.2 > 0. All three tests agree.
| test | value (verified) | verdict |
|---|---|---|
| symmetric | True | ok |
| eigenvalues | [0.5566, 2.4434] | both > 0 -> PD |
| det(Sigma) | 1.36 | > 0 |
| a^T Sigma a at a=[3,-1] | 14.2 | > 0 |
Discrimination
Sort into buckets
Sort these by verdict, from memory, without looking back at Our Sigma is PD: check the quadratic form and…. Telling them apart on the spot is the skill; the table is only where the answer happens to be written down.
Anomaly
Predict first
A student writes this, and it looks reasonable:
It is symmetric and has ones on the diagonal, so [[1, 2], [2, 1]] is a fine covariance matrix.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: Its eigenvalues are -1 and 3. Take a = [1, -1]: a^T Sigma a = 1 - 2 - 2 + 1 = -2 < 0.
Require every eigenvalue >= 0 (PSD), and > 0 (PD) before trusting it as a Sigma.
Why: Its eigenvalues are -1 and 3. Take a = [1, -1]: a^T Sigma a = 1 - 2 - 2 + 1 = -2 < 0. That would be a NEGATIVE variance along direction [1, -1] - impossible. The matrix is indefinite, not a covariance.
Trap
It is symmetric and has ones on the diagonal, so [[1, 2], [2, 1]] is a fine covariance matrix.
Treat [[1, 2], [2, 1]] as Sigma
Why: Its eigenvalues are -1 and 3. Take a = [1, -1]: a^T Sigma a = 1 - 2 - 2 + 1 = -2 < 0. That would be a NEGATIVE variance along direction [1, -1] - impossible. The matrix is indefinite, not a covariance.
Require every eigenvalue >= 0 (PSD), and > 0 (PD) before trusting it as a Sigma.
Check eigvalsh >= 0 (or attempt a Cholesky)
Why: Our [[2, 0.8], [0.8, 1]] has eigenvalues 0.5566 and 2.4434, both > 0. A failed np.linalg.cholesky is the fastest PD test - it raises the instant a matrix is not positive definite.
Concept
A PSD but not PD Sigma has a zero eigenvalue. Then |Sigma| = 0, Sigma^{-1} does not exist, and the standard density formula breaks - the distribution collapses onto a lower-dimensional line or plane (all its mass on a subspace).
That is still a valid Gaussian in a degenerate sense, but you cannot write its density with Sigma^{-1}. For everything below we assume Sigma is strictly PD, so |Sigma| > 0 and Sigma^{-1} exists.
Section
Part 3 of 7 - geometry & distance
Picture it
Figure (svg): Nested ellipses centered at a point mu, tilted along a diagonal direction, with a longer axis and a shorter axis drawn through the center.
Discussion prompt
Read the picture before the words. What is this showing, and what is the one thing it is built to make obvious? Commit to an answer, then read on.
Hint: Name the parts, then say what changes between them — and if nothing changes, say what is being held still.
Answer:
The density depends on x only through the exponent (x - mu)^T Sigma^{-1} (x - mu). Fixing that quadratic form to a constant c traces out all the equally-likely points - and that set is an ellipse centered at mu.
Intuition
The density depends on x only through the exponent (x - mu)^T Sigma^{-1} (x - mu). Fixing that quadratic form to a constant c traces out all the equally-likely points - and that set is an ellipse centered at mu.
Bigger c = further from center = lower density. So the density is a smooth hill over the plane whose contour lines are nested ellipses, tallest at mu.
Figure (svg): Nested ellipses centered at a point mu, tilted along a diagonal direction, with a longer axis and a shorter axis drawn through the center.
Concept
Sigma = Q Lambda Q^T (symmetric eigen-decomposition, Lesson 10). The eigenvectors (columns of Q) point along the ellipse's axes; the eigenvalues in Lambda are the variances along those axes.
A 1-standard-deviation contour reaches sqrt(lambda_k) along axis k. Large eigenvalue = long axis (high spread); small eigenvalue = short axis. Off-diagonal terms in Sigma are exactly what rotate the axes off the coordinate lines.
Missing information
Discussion prompt
Eigen-decompose our Sigma to get the axis directions and lengths. np.linalg.eigh is the symmetric solver - real eigenvalues, orthonormal eigenvectors, ascending order:
What do you need to know — or decide — before the first line can be written? List everything the problem has to hand you.
Hint: Anything you would have to invent to get started is a thing the problem must supply.
Answer:
trace = 2 + 1 = 3 and det = 1.36, so eigenvalues solve l^2 - (trace) l + det = 0. sqrt(3.56) = 1.8868, giving l = 0.5566 and 2.4434 - matching eigh exactly.
Worked example
Eigen-decompose our Sigma to get the axis directions and lengths. np.linalg.eigh is the symmetric solver - real eigenvalues, orthonormal eigenvectors, ascending order:
import numpy as np
Sigma = np.array([[2., 0.8], [0.8, 1.]])
vals, vecs = np.linalg.eigh(Sigma) # ascending; symmetric -> real
print(vals.round(4)) # axis variances
print(np.sqrt(vals).round(4)) # 1-sigma semi-axis lengths
print(vecs.round(4)) # columns = axis directionsChar. poly l^2 - 3l + 1.36 = 0 -> l = (3 +- sqrt(3.56))/2
Why: trace = 2 + 1 = 3 and det = 1.36, so eigenvalues solve l^2 - (trace) l + det = 0. sqrt(3.56) = 1.8868, giving l = 0.5566 and 2.4434 - matching eigh exactly.
| axis | eigenvalue (variance) | semi-axis sqrt(lambda) |
|---|---|---|
| short | 0.5566 | 0.7461 |
| long | 2.4434 | 1.5631 |
Reverse engineer
Discussion prompt
Work backwards. The example finished here:
Char. poly l^2 - 3l + 1.36 = 0 -> l = (3 +- sqrt(3.56))/2
What was it asked to do, and what must it have been given? Reconstruct the problem from its answer.
Hint: Every quantity in the result had to enter somewhere. Account for each one.
Answer:
Eigen-decompose our Sigma to get the axis directions and lengths. np.linalg.eigh is the symmetric solver - real eigenvalues, orthonormal eigenvectors, ascending order:
Concept
The quadratic form in the exponent has a name - the squared Mahalanobis distance from x to the center:
\[ d_M^2(x) = (x-\mu)^\top \Sigma^{-1}(x-\mu) \]
It measures how far x is from mu in standard-deviation units along the ellipse axes. Points on the same contour ellipse are all the same Mahalanobis distance from the center - equally surprising - even though their raw Euclidean distances differ.
Intuition
Euclidean distance treats every direction the same, as if the cloud were a perfect circle. But our cloud is a tilted, stretched ellipse.
A point 2 units out along the long axis is still deep inside the cloud - totally typical. The same 2 units across the short axis pokes out the side - a genuine outlier. Mahalanobis divides by the spread in each direction, so it calls the first point close and the second far. Euclidean cannot tell them apart.
Pattern
Predict first
The table runs: det(Sigma) | 1.36 · (Sigma^{-1})_11 | 0.7353 · (Sigma^{-1})_12 = _21 | -0.5882
In Invert Sigma by hand, given the rows so far: what is the next one — the row where entry is (Sigma^{-1})_22?
Correct: (Sigma^{-1})_22 | 1.4706
| entry | value (verified) |
|---|---|
| det(Sigma) | 1.36 |
| (Sigma^{-1})_11 | 0.7353 |
| (Sigma^{-1})_12 = _21 | -0.5882 |
| (Sigma^{-1})_22 | 1.4706 |
Why: The relationship between the columns, not the individual numbers, is what generates the next row. ad - bc with a=2, b=c=0.8, d=1.
Worked example
To compute Mahalanobis we need Sigma^{-1}. For a 2x2 [[a, b], [c, d]], the inverse is 1/det * [[d, -b], [-c, a]]. First the determinant:
det = 21 - 0.80.8 = 2 - 0.64 = 1.36
Why: ad - bc with a=2, b=c=0.8, d=1. Nonzero, so the inverse exists (PD confirmed again).
\[ \det \Sigma = 2\cdot 1 - 0.8\cdot 0.8 = 1.36 \]
Swap the diagonal, negate the off-diagonal, divide by 1.36
Why: The 2x2 inverse formula. Divide each entry of [[1, -0.8], [-0.8, 2]] by 1.36.
\[ \Sigma^{-1} = \frac{1}{1.36}\begin{bmatrix} 1 & -0.8 \\ -0.8 & 2 \end{bmatrix} = \begin{bmatrix} 0.7353 & -0.5882 \\ -0.5882 & 1.4706 \end{bmatrix} \]
| entry | value (verified) |
|---|---|
| det(Sigma) | 1.36 |
| (Sigma^{-1})_11 | 0.7353 |
| (Sigma^{-1})_12 = _21 | -0.5882 |
| (Sigma^{-1})_22 | 1.4706 |
Blank canvas
Draw it
Draw what Invert Sigma by hand just did — the shape of it, not the line-by-line working. One picture, labels only where you need them. Then check it against the steps: anything you could not draw is a step you followed rather than understood.
Fill the middle
Fill in the blanks
From Mahalanobis vs Euclidean to x = [3, 0] — one line has had its right-hand side removed. Put it back.
import numpy as np
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
x = np.array([3., 0.])
diff = x - mu
d2_M = diff @ np.linalg.inv(Sigma) @ diff
d2_E = (diff**2).sum()
print(round(d2_M, 4), round(d2_E, 4))
Why: diff is what everything below it consumes, so the wrong expression here fails later and somewhere else. Row 1: 12 + (-0.8)1 = 1.2. Row 2: (-0.8)2 + 21 = 0.4.
Worked example
Now measure the distance from x = [3, 0] to mu = [1, -1] both ways. The miss vector is x - mu = [2, 1].
Sigma^{-1}(x - mu) = (1/1.36)[[1,-0.8],[-0.8,2]] [2,1] = (1/1.36)[1.2, 0.4]
Why: Row 1: 12 + (-0.8)1 = 1.2. Row 2: (-0.8)2 + 21 = 0.4. Then divide by 1.36 -> [0.8824, 0.2941].
d_M^2 = (x-mu)^T [that] = 20.8824 + 10.2941 = 2.0588
Why: Equivalently (21.2 + 10.4)/1.36 = 2.8/1.36 = 2.0588. Verified in code below.
\[ d_M^2 = \frac{2\cdot 1.2 + 1\cdot 0.4}{1.36} = \frac{2.8}{1.36} = 2.0588 \]
import numpy as np
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
x = np.array([3., 0.])
diff = x - mu
d2_M = diff @ np.linalg.inv(Sigma) @ diff
d2_E = (diff**2).sum()
print(round(d2_M, 4), round(d2_E, 4))| distance | value | reading |
|---|---|---|
| Euclidean^2 | 5.0 | looks far |
| Mahalanobis^2 | 2.0588 | actually typical |
| ratio | 5.0 / 2.06 = 2.4x | Euclidean overstates by 2.4x |
Error analysis
Annotate
Walk the callouts on Mahalanobis vs Euclidean to x = [3, 0]. Each one is a place this is easy to get subtly wrong.
Concept
Set Sigma = I (unit variance, no correlation). Then Sigma^{-1} = I and d_M^2 = (x - mu)^T (x - mu) = ||x - mu||^2 - plain squared Euclidean distance.
So Euclidean is the special case of Mahalanobis for an isotropic (spherical) Gaussian. Mahalanobis is the honest generalization once the cloud stops being a perfect ball.
Anomaly
Predict first
A student writes this, and it looks reasonable:
To find anomalies in correlated data, rank points by raw distance from the mean, ||x - mu||.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: Ignores Sigma entirely. Our x = [3, 0] scores Euclidean^2 = 5.0 and looks like a far outlier, but along the high-variance tilt it is only 2.06 in Mahalanobis - a perfectly ordinary point.
Rank by Mahalanobis distance, which rescales by Sigma^{-1}.
Why: Ignores Sigma entirely. Our x = [3, 0] scores Euclidean^2 = 5.0 and looks like a far outlier, but along the high-variance tilt it is only 2.06 in Mahalanobis - a perfectly ordinary point. You would flag typical points and miss real ones.
Trap
To find anomalies in correlated data, rank points by raw distance from the mean, ||x - mu||.
Flag x by Euclidean ||x - mu||
Why: Ignores Sigma entirely. Our x = [3, 0] scores Euclidean^2 = 5.0 and looks like a far outlier, but along the high-variance tilt it is only 2.06 in Mahalanobis - a perfectly ordinary point. You would flag typical points and miss real ones.
Rank by Mahalanobis distance, which rescales by Sigma^{-1}.
Flag x by d_M^2 = (x-mu)^T Sigma^{-1} (x-mu)
Why: Distance in standard-deviation units along the ellipse axes - the correct notion of 'far' for correlated data. Under a Gaussian, d_M^2 follows a chi-square(d) law, so you can even set a principled threshold (e.g. 5.991 for 95% in 2-D).
Break the constraint
Discussion prompt
The rule this trap just fixed:
Rank by Mahalanobis distance, which rescales by Sigma^{-1}.
Now break it on purpose. Build a case that violates it and follow the consequences until something visibly fails. Where does the failure first show up — and would you have noticed it if you had not been looking?
Hint: The dangerous rules are the ones whose violation still produces an answer. If yours fails loudly, try to find one that fails quietly.
Answer:
Ignores Sigma entirely. Our x = [3, 0] scores Euclidean^2 = 5.0 and looks like a far outlier, but along the high-variance tilt it is only 2.06 in Mahalanobis - a perfectly ordinary point. You would flag typical points and miss real ones.
Section
Part 4 of 7 - determinant & PDF
Concept
The front constant makes the total probability integrate to 1. In 1-D that job needs 1/sqrt(2 pi sigma^2); the sigma^2 accounts for how wide the bump is.
In d dimensions the 'width' is a volume, and the volume of the covariance ellipse scales with sqrt(|Sigma|) - the product of the axis lengths sqrt(lambda_1) ... sqrt(lambda_d), since |Sigma| = product of eigenvalues. That is exactly the |Sigma|^{1/2} in the denominator.
Worked example
At x = mu the exponent is 0 (distance zero), so exp(0) = 1 and the density equals the bare front constant. For d = 2:
const = 1 / ((2 pi)^{d/2} |Sigma|^{1/2}) = 1 / (2 pi * sqrt(1.36))
Why: d/2 = 1, so (2 pi)^1 = 2 pi, and |Sigma|^{1/2} = sqrt(1.36). Multiply: 2 pi * 1.1662 = 7.328.
\[ p(\mu) = \frac{1}{2\pi\,\sqrt{1.36}} = \frac{1}{7.328} = 0.1365 \]
| quantity | value (verified) |
|---|---|
| det(Sigma) | 1.36 |
| sqrt(det) | 1.1662 |
| (2 pi)^{d/2} = 2 pi | 6.2832 |
| peak p(mu) | 0.1365 |
Estimation
Predict first
Combine everything: the constant 0.1365, and the exponent using the Mahalanobis d_M^2 = 2.0588 we already computed. Runnable, from scratch:
Commit before you compute: what does Evaluate the full density at x = [3, 0] come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: p(x) = 0.1365 * exp(-0.5 * 2.0588) = 0.1365 * 0.3573 = 0.0488
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. exp(-1.0294) = 0.3573; times the peak 0.1365 gives 0.0488.
Worked example
Combine everything: the constant 0.1365, and the exponent using the Mahalanobis d_M^2 = 2.0588 we already computed. Runnable, from scratch:
import numpy as np
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
d = len(mu)
x = np.array([3., 0.])
det = np.linalg.det(Sigma)
dM2 = (x - mu) @ np.linalg.inv(Sigma) @ (x - mu)
const = 1.0 / ((2*np.pi)**(d/2) * np.sqrt(det))
pdf = const * np.exp(-0.5 * dM2)
print(round(det, 4), round(dM2, 4))
print(round(const, 6), round(pdf, 6))p(x) = 0.1365 * exp(-0.5 * 2.0588) = 0.1365 * 0.3573 = 0.0488
Why: exp(-1.0294) = 0.3573; times the peak 0.1365 gives 0.0488. The point sits at about 36% of the peak height - down but not deep in the tail, consistent with its modest Mahalanobis distance.
| quantity | value (verified) |
|---|---|
| det(Sigma) | 1.36 |
| d_M^2 at [3,0] | 2.0588 |
| const (peak) | 0.136474 |
| p([3, 0]) | 0.048751 |
Pattern
Step through it
Step through Evaluate the full density at x = [3, 0] one row at a time. What is driving the change, and what would the row after the last one be?
Fill the middle
Fill in the blanks
From Cross-check against scipy — one line has had its right-hand side removed. Put it back.
import numpy as np
from scipy.stats import multivariate_normal
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
rv = multivariate_normal(mean=mu, cov=Sigma)
print(round(rv.pdf([3., 0.]), 6)) # our hand value
print(round(rv.pdf(mu), 6)) # peak, at the center
Why: Sigma is what everything below it consumes, so the wrong expression here fails later and somewhere else. Our term-by-term density is correct to 6 decimals against a battle-tested library.
Worked example
Never trust a hand PDF without an independent oracle. scipy.stats.multivariate_normal computes the same density - it must match our 0.048751.
import numpy as np
from scipy.stats import multivariate_normal
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
rv = multivariate_normal(mean=mu, cov=Sigma)
print(round(rv.pdf([3., 0.]), 6)) # our hand value
print(round(rv.pdf(mu), 6)) # peak, at the centerscipy prints 0.048751 and 0.136474 - exact match
Why: Our term-by-term density is correct to 6 decimals against a battle-tested library. The peak 0.136474 equals the front constant, confirming exp(0) = 1 at the center.
| point | our formula | scipy |
|---|---|---|
| x = [3, 0] | 0.048751 | 0.048751 |
| x = mu | 0.136474 | 0.136474 |
Section
Part 5 of 7 - closure properties
Concept
A defining superpower of the Gaussian: if x is jointly Gaussian, then every marginal (drop some coordinates) is Gaussian, and every conditional (fix some coordinates) is Gaussian. You never leave the family.
This is why Gaussians dominate ML: Kalman filters, Gaussian processes, and VAEs all rely on the fact that slicing and conditioning keep you Gaussian, so the math stays closed-form.
Fill the middle
Fill in the blanks
From The marginal of x_1: just read off Sigma — finish the line. Write what belongs on the right of the equals sign before you look.
x_1 \sim \mathcal\mathcal{N}(1,\ 2)(\mu_1,\ \Sigma____) = ___
Why: Producing the right-hand side unprompted is the difference between recognising this line and being able to use it. No conditioning, no correlation correction - the covariance term 0.8 does not enter a marginal.
Worked example
The marginal of a single coordinate is astonishingly simple - you ignore the other coordinates entirely and read the corresponding mean and variance straight out of mu and Sigma:
\[ x_1 \sim \mathcal{N}(\mu_1,\ \Sigma_{11}) = \mathcal{N}(1,\ 2) \]
Marginal mean = mu_1 = 1, marginal variance = Sigma_11 = 2
Why: No conditioning, no correlation correction - the covariance term 0.8 does not enter a marginal. Confirm on the sample: x_1 has mean 1.0018, variance 2.0067.
| quantity | theory | from 200k samples |
|---|---|---|
| mean of x_1 | 1 | 1.0018 |
| var of x_1 | 2 | 2.0067 |
Intuition
The marginal ignored x_2. But the conditional x_2 | x_1 = 3 asks: given we observed a specific x_1, what do we now believe about x_2? Because the two are correlated (Sigma_12 = 0.8 > 0), learning x_1 is high pulls our estimate of x_2 up too.
Two things change versus the marginal: the mean shifts toward the correlated evidence, and the variance shrinks because we now know more. Both effects are captured by the conditioning formulas next.
Worked example
For a bivariate Gaussian, conditioning x_2 on x_1 = v gives another Gaussian with a corrected mean and a reduced variance:
\[ \mathbb{E}[x_2 \mid x_1 = v] = \mu_2 + \frac{\Sigma_{21}}{\Sigma_{11}}\,(v - \mu_1) \]
\[ \operatorname{Var}[x_2 \mid x_1 = v] = \Sigma_{22} - \frac{\Sigma_{21}\,\Sigma_{12}}{\Sigma_{11}} \]
The mean gains a slope term; the variance loses a correlation term
Why: Sigma_21 / Sigma_11 is the regression slope of x_2 on x_1. The variance drop Sigma_21 Sigma_12 / Sigma_11 is the part of x_2's uncertainty explained by x_1 - it is never negative, so conditioning never increases variance.
Translation
\( \mathbb{E}[x_2 \mid x_1 = v] = \mu_2 + \frac{\Sigma_{21}}{\Sigma_{11}}\,(v - \mu_1) \)
Draw it
Translate both ways. First write the expression above as a sentence with no symbols in it at all. Then cover it, and write your sentence back as notation. If the two versions disagree, the disagreement is the thing to fix.
Pattern
Predict first
The table runs: conditional mean | -0.2 | -0.177 · conditional var | 0.68 | 0.684
In Condition x_2 on x_1 = 3, by hand then by sample, given the rows so far: what is the next one — the row where quantity is n in slice?
Correct: n in slice | - | 1934 samples
| quantity | by hand | empirical (x_1 approx 3) |
|---|---|---|
| conditional mean | -0.2 | -0.177 |
| conditional var | 0.68 | 0.684 |
| n in slice | - | 1934 samples |
Why: The relationship between the columns, not the individual numbers, is what generates the next row. Slope 0.8/2 = 0.4; the observed x_1 = 3 is 2 above its mean 1, so x_2's estimate rises by 0.4*2 = 0.8, from -1 up to -0.2.
Worked example
Mean: -1 + (0.8/2)(3 - 1) = -1 + 0.42 = -0.2
Why: Slope 0.8/2 = 0.4; the observed x_1 = 3 is 2 above its mean 1, so x_2's estimate rises by 0.4*2 = 0.8, from -1 up to -0.2.
\[ \mathbb{E}[x_2 \mid x_1 = 3] = -1 + 0.4\,(3-1) = -0.2 \]
Variance: 1 - (0.8*0.8)/2 = 1 - 0.32 = 0.68
Why: Knowing x_1 removes 0.32 of x_2's original variance of 1, leaving 0.68. Correlation buys certainty.
\[ \operatorname{Var}[x_2 \mid x_1 = 3] = 1 - \frac{0.64}{2} = 0.68 \]
import numpy as np
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
x1 = 3.0
cond_mean = mu[1] + Sigma[1,0]*(1/Sigma[0,0])*(x1 - mu[0])
cond_var = Sigma[1,1] - Sigma[1,0]*(1/Sigma[0,0])*Sigma[0,1]
print(round(cond_mean, 4), round(cond_var, 4))
L = np.linalg.cholesky(Sigma)
X = mu + np.random.default_rng(0).normal(size=(200000, 2)) @ L.T
m = np.abs(X[:, 0] - 3.0) < 0.05
print(round(X[m, 1].mean(), 3), round(X[m, 1].var(), 3))| quantity | by hand | empirical (x_1 approx 3) |
|---|---|---|
| conditional mean | -0.2 | -0.177 |
| conditional var | 0.68 | 0.684 |
| n in slice | - | 1934 samples |
Comparison
Comparison matrix
From Condition x_2 on x_1 = 3, by hand then by sample: refill the by hand column from what you know. The rest of the table is as it appeared.
| quantity | by hand | empirical (x_1 approx 3) |
|---|---|---|
| conditional mean | -0.2 | -0.177 |
| conditional var | 0.68 | 0.684 |
| n in slice | - | 1934 samples |
Concept
Special fact for Gaussians: if Sigma is diagonal (all off-diagonal covariances zero), the components are independent - not merely uncorrelated. For a Gaussian, zero covariance really is full independence.
The engine is factorization: when Sigma is diagonal, the joint density splits into a product of 1-D Gaussians, one per coordinate - and a density that factorizes is the definition of independence. We prove it next.
Step zero
Discussion prompt
Prove diagonal Sigma factorizes the density — before any calculation: what is the plan? Name the moves in order, in plain English, without doing the arithmetic.
Hint: It starts with: Take Sigma = diag(sigma_1^2, sigma_2^2); then Sigma^{-1} =…
Answer:
Worked example
Take Sigma = diag(sigma_1^2, sigma_2^2); then Sigma^{-1} = diag(1/sigma_1^2, 1/sigma_2^2)
Why: Inverting a diagonal matrix just inverts each diagonal entry. |Sigma| = sigma_1^2 sigma_2^2.
The quadratic form becomes a plain sum, no cross term
Why: With Sigma^{-1} diagonal, (x-mu)^T Sigma^{-1} (x-mu) = (x_1-mu_1)^2/sigma_1^2 + (x_2-mu_2)^2/sigma_2^2. No off-diagonal, so no x_1 x_2 coupling survives.
\[ (x-\mu)^\top \Sigma^{-1}(x-\mu) = \frac{(x_1-\mu_1)^2}{\sigma_1^2} + \frac{(x_2-\mu_2)^2}{\sigma_2^2} \]
exp of a sum = product of exps, and the constant splits too
Why: |Sigma|^{1/2} = sigma_1 sigma_2 and (2 pi)^{d/2} = (2 pi)^1 = sqrt(2 pi)^2, so the whole density separates into the two 1-D densities.
\[ p(x) = \underbrace{\frac{1}{\sqrt{2\pi}\,\sigma_1}e^{-\frac{(x_1-\mu_1)^2}{2\sigma_1^2}}}_{p(x_1)} \cdot \underbrace{\frac{1}{\sqrt{2\pi}\,\sigma_2}e^{-\frac{(x_2-\mu_2)^2}{2\sigma_2^2}}}_{p(x_2)} \]
import numpy as np
from scipy.stats import multivariate_normal, norm
Sd = np.array([[2., 0.], [0., 1.]]) # diagonal covariance
joint = multivariate_normal(mean=[0., 0.], cov=Sd).pdf([1.5, -0.5])
prod = norm.pdf(1.5, 0, np.sqrt(2.)) * norm.pdf(-0.5, 0, 1.)
print(round(joint, 8))
print(round(prod, 8))
print(np.isclose(joint, prod))| quantity | value (verified) |
|---|---|
| joint p(1.5, -0.5) | 0.05658843 |
| p(x_1) * p(x_2) | 0.05658843 |
| equal? | True -> independent |
Error analysis
Annotate
Walk the callouts on Prove diagonal Sigma factorizes the density. Each one is a place this is easy to get subtly wrong.
Anomaly
Predict first
A student writes this, and it looks reasonable:
Zero covariance means independent, for any distribution - so Cov(X, Y) = 0 always gives X independent of Y.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: Cov(X, X^2) = E[X^3] - E[X]E[X^2] = 0 - 0 = 0 (odd moments of a symmetric X vanish).
Uncorrelated implies independent only for jointly Gaussian vectors. That implication is special to the Gaussian.
Why: Cov(X, X^2) = E[X^3] - E[X]E[X^2] = 0 - 0 = 0 (odd moments of a symmetric X vanish). Verified: sampled Cov = 0.0006. Yet Y is a DETERMINISTIC function of X - knowing X fixes Y exactly. Total dependence, zero covariance.
Trap
Zero covariance means independent, for any distribution - so Cov(X, Y) = 0 always gives X independent of Y.
Take Y = X^2 with X ~ N(0, 1) and claim independence from Cov = 0
Why: Cov(X, X^2) = E[X^3] - E[X]E[X^2] = 0 - 0 = 0 (odd moments of a symmetric X vanish). Verified: sampled Cov = 0.0006. Yet Y is a DETERMINISTIC function of X - knowing X fixes Y exactly. Total dependence, zero covariance.
Uncorrelated implies independent only for jointly Gaussian vectors. That implication is special to the Gaussian.
Diagonal Sigma => independent, but ONLY inside the Gaussian family
Why: The density factorizes (previous slide) precisely because the pair is jointly Gaussian. Outside that family, zero correlation captures only linear dependence - X and X^2 are the classic warning.
Section
Part 6 of 7 - the reparameterization
Intuition
Your computer only knows how to draw standard normals z ~ N(0, I) - an isotropic round cloud centered at the origin, every direction identical.
We want a tilted, stretched, shifted ellipse. So we need a transform that takes the round cloud z and turns it into our Sigma-shaped cloud around mu. A single linear map does it: stretch/rotate by a matrix, then shift by mu.
Concept
The one rule we need: if z has covariance C and we form x = mu + A z, then x has covariance A C A^T. Shifting by mu moves the center but never changes the spread.
\[ \operatorname{Cov}(\mu + A z) = A\,\operatorname{Cov}(z)\,A^\top \]
With z ~ N(0, I) we have Cov(z) = I, so Cov(mu + A z) = A A^T. To hit target Sigma, we just need a matrix A with A A^T = Sigma - a square root of Sigma.
Concept
The Cholesky factorization (Lesson 13) of a PD matrix produces a lower-triangular L with Sigma = L L^T. That L is exactly the A we need.
\[ \Sigma = L L^\top \;\Longrightarrow\; x = \mu + L z,\quad z \sim \mathcal{N}(0, I) \]
And Gaussians stay Gaussian under linear maps, so x is not just any distribution with covariance Sigma - it is precisely N(mu, Sigma). Triangular L also makes the transform cheap and numerically clean.
Explain it
Discussion prompt
Explain Cholesky gives that square root to a student a year behind you. No notation, no jargon they have not met — and it still has to be true.
Hint: If your explanation needs a symbol they have never seen, you are describing the notation rather than the idea.
Answer:
The Cholesky factorization (Lesson 13) of a PD matrix produces a lower-triangular L with Sigma = L L^T. That L is exactly the A we need.
Ranking
Put in order
Put the moves of Prove Cov(mu + Lz) = Sigma, one line at a time into the order they have to happen.
Why: These are the moves of the worked example in the order it makes them, and each one is set up by the one before it. z is standard normal, so E[z] = 0 and E[z z^T] = I.
Worked example
Start: x = mu + L z, with Cov(z) = I
Why: z is standard normal, so E[z] = 0 and E[z z^T] = I. This is the only fact we need about z.
Center x: x - mu = L z (the mu cancels)
Why: E[x] = mu + L E[z] = mu, so the centered vector is exactly L z. Covariance only sees the centered vector.
\[ x - \mu = L z \]
Cov(x) = E[(x-mu)(x-mu)^T] = E[L z (L z)^T]
Why: Definition of covariance, substituting x - mu = L z. And (L z)^T = z^T L^T.
\[ \operatorname{Cov}(x) = \mathbb{E}[\,L z\, z^\top L^\top\,] \]
Pull the constants L and L^T outside the expectation
Why: L is not random, so E[L M L^T] = L E[M] L^T. The middle expectation E[z z^T] is I.
\[ = L\,\mathbb{E}[z z^\top]\,L^\top = L\,I\,L^\top = L L^\top \]
L L^T = Sigma by the Cholesky definition. Done.
Why: So Cov(mu + L z) = Sigma exactly - the construction is not an approximation, it is an identity. Any square root of Sigma would work; Cholesky is the cheap triangular choice.
\[ \operatorname{Cov}(x) = L L^\top = \Sigma \]
Blank canvas
Draw it
Draw what Prove Cov(mu + Lz) = Sigma, one line at a time just did — the shape of it, not the line-by-line working. One picture, labels only where you need them. Then check it against the steps: anything you could not draw is a step you followed rather than understood.
Ranking
Put in order
Put the moves of Compute L by hand into the order they have to happen.
Why: These are the moves of the worked example in the order it makes them, and each one is set up by the one before it. Top-left: L11^2 = Sigma_11 = 2, so L11 = sqrt(2).
Worked example
Cholesky of a 2x2 [[a, b], [b, c]] gives lower-triangular L = [[L11, 0], [L21, L22]] solved top-down: L11 = sqrt(a), L21 = b/L11, L22 = sqrt(c - L21^2).
L11 = sqrt(2) = 1.4142
Why: Top-left: L11^2 = Sigma_11 = 2, so L11 = sqrt(2).
L21 = 0.8 / sqrt(2) = 0.5657
Why: Off-diagonal: L21 * L11 = Sigma_21 = 0.8, so L21 = 0.8/1.4142.
L22 = sqrt(1 - 0.5657^2) = sqrt(1 - 0.32) = 0.8246
Why: Bottom-right: L21^2 + L22^2 = Sigma_22 = 1, so L22 = sqrt(1 - 0.32) = sqrt(0.68).
\[ L = \begin{bmatrix} 1.4142 & 0 \\ 0.5657 & 0.8246 \end{bmatrix} \]
import numpy as np
Sigma = np.array([[2., 0.8], [0.8, 1.]])
L = np.linalg.cholesky(Sigma) # lower-triangular
print(L.round(4))
print(np.allclose(L @ L.T, Sigma))| entry | hand value | np value |
|---|---|---|
| L11 | 1.4142 | 1.4142 |
| L21 | 0.5657 | 0.5657 |
| L22 | 0.8246 | 0.8246 |
| L L^T == Sigma | - | True |
Notation
Annotate
From Compute L by hand — read this one piece at a time. What is each part doing?
On: \( L = \begin{bmatrix} 1.4142 & 0 \\ 0.5657 & 0.8246 \end{bmatrix} \)
Missing information
Discussion prompt
Draw standard normals, transform with mu + z @ L.T (row-wise: each row is mu + L z), then check the empirical mean and covariance match. Seed 0 for reproducibility:
What do you need to know — or decide — before the first line can be written? List everything the problem has to hand you.
Hint: Anything you would have to invent to get started is a thing the problem must supply.
Answer:
Both recovered to ~2 decimals - the Cholesky transform reproduces mu AND Sigma, including the 0.8 correlation. With 200k draws the Monte-Carlo error is small; more samples tighten it further.
Worked example
Draw standard normals, transform with mu + z @ L.T (row-wise: each row is mu + L z), then check the empirical mean and covariance match. Seed 0 for reproducibility:
import numpy as np
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
L = np.linalg.cholesky(Sigma)
rng = np.random.default_rng(0)
z = rng.normal(size=(200000, 2)) # standard normals, N(0, I)
X = mu + z @ L.T # x = mu + L z, row-wise
print(X.mean(0).round(3))
print(np.cov(X.T).round(3))sample mean [1.002, -1.0]; sample cov [[2.007, 0.798], [0.798, 0.999]]
Why: Both recovered to ~2 decimals - the Cholesky transform reproduces mu AND Sigma, including the 0.8 correlation. With 200k draws the Monte-Carlo error is small; more samples tighten it further.
| quantity | target | from 200k samples |
|---|---|---|
| mean | [1, -1] | [1.002, -1.000] |
| Sigma_11 (var x_1) | 2.0 | 2.007 |
| Sigma_12 (cov) | 0.8 | 0.798 |
| Sigma_22 (var x_2) | 1.0 | 0.999 |
Anomaly
Predict first
A student writes this, and it looks reasonable:
To get covariance Sigma, scale the standard normals by Sigma directly: x = mu + Sigma z.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: Then Cov(x) = Sigma I Sigma^T = Sigma^2, NOT Sigma.
Scale by the square root L with L L^T = Sigma, not by Sigma itself.
Why: Then Cov(x) = Sigma I Sigma^T = Sigma^2, NOT Sigma. You have squared the covariance. For our Sigma, Sigma^2 = [[4.64, 2.4], [2.4, 1.64]] - the wrong spread and wrong correlation entirely.
Trap
To get covariance Sigma, scale the standard normals by Sigma directly: x = mu + Sigma z.
x = mu + Sigma z
Why: Then Cov(x) = Sigma I Sigma^T = Sigma^2, NOT Sigma. You have squared the covariance. For our Sigma, Sigma^2 = [[4.64, 2.4], [2.4, 1.64]] - the wrong spread and wrong correlation entirely.
Scale by the square root L with L L^T = Sigma, not by Sigma itself.
x = mu + L z, L = cholesky(Sigma)
Why: Cov(x) = L I L^T = L L^T = Sigma - exactly right, as we proved. The rule of thumb: to hit covariance Sigma you always need a matrix square root of Sigma, and Cholesky is the standard one.
Two truths and a lie
Sort into buckets
Some of these hold up and some are the exact mistakes this lesson is built to prevent. Sort them.
x with mean mu and variance sigma^2 has density; In 1-D the exponent divides by sigma^2 - it rescales the squared miss into variance units. A miss of 2 is a big deal if sigma^2 = 0.1 and trivial if sigma^2 = 100.; Bigger c = further from center = lower density. So the density is a smooth hill over the plane whose contour lines are nested ellipses, tallest at mu.[[1, 2], [2, 1]] is a fine covariance matrix.; To find anomalies in correlated data, rank points by raw distance from the mean, ||x - mu||.Section
Part 7 of 7 - why this matters
Concept
A variational autoencoder encodes each input into a Gaussian latent distribution: the encoder outputs a mean mu and a (usually diagonal) covariance Sigma per input, then the decoder reconstructs from a sample z of that Gaussian.
But sampling is random, and you cannot backpropagate through a raw random draw - the gradient has nowhere to flow. That is exactly the problem our Cholesky sampler solves.
Concept
Rewrite the draw as z = mu + L * eps with eps ~ N(0, I). Now the ONLY random node is eps; mu and L are deterministic outputs of the network.
\[ z = \mu + L\,\epsilon, \qquad \epsilon \sim \mathcal{N}(0, I) \]
Gradients flow straight through mu and L because they sit on a plain deterministic path - the randomness is pushed into the fixed eps. This is the identical x = mu + L z you just proved, wearing a training hat. For diagonal Sigma, L is just the per-dimension standard deviations.
Analogy
Discussion prompt
Explain The reparameterization trick by analogy to something with no Machine Learning in it at all — a queue, a recipe, a map, a bank balance, whatever fits. Then say where your analogy breaks.
Hint: An analogy that never breaks is not an analogy, it is the same idea wearing a hat. Find the seam — that is the part that is actually new.
Answer:
Rewrite the draw as z = mu + L * eps with eps ~ N(0, I). Now the ONLY random node is eps; mu and L are deterministic outputs of the network.
Estimation
Predict first
In a real VAE Sigma is diagonal, so L = diag(sigma) and the transform is an elementwise mu + sigma * eps - no matrix multiply. Demonstrate that it lands the right mean and per-dim std:
Commit before you compute: what does Diagonal-Sigma reparameterization (the VAE case) come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: sample mean [0.502, -2.0]; sample std [1.001, 0.3]
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. Elementwise mu + sigma * eps reproduces the target mean and standard deviations exactly - and because eps is the sole random node, autograd differentiates cleanly through mu and sigma.
Worked example
In a real VAE Sigma is diagonal, so L = diag(sigma) and the transform is an elementwise mu + sigma * eps - no matrix multiply. Demonstrate that it lands the right mean and per-dim std:
import numpy as np
mu = np.array([0.5, -2.0]) # encoder output: latent mean
sigma = np.array([1.0, 0.3]) # encoder output: per-dim std (diagonal L)
rng = np.random.default_rng(0)
eps = rng.normal(size=(500000, 2)) # noise, the ONLY random node
z = mu + sigma * eps # reparameterization: gradients flow to mu, sigma
print(z.mean(0).round(3))
print(z.std(0).round(3))sample mean [0.502, -2.0]; sample std [1.001, 0.3]
Why: Elementwise mu + sigma * eps reproduces the target mean and standard deviations exactly - and because eps is the sole random node, autograd differentiates cleanly through mu and sigma. This is the whole sampler inside a VAE forward pass.
| quantity | target | from samples |
|---|---|---|
| mean | [0.5, -2.0] | [0.502, -2.000] |
| std dim 0 | 1.0 | 1.001 |
| std dim 1 | 0.3 | 0.300 |
Reverse engineer
Discussion prompt
Work backwards. The example finished here:
sample mean [0.502, -2.0]; sample std [1.001, 0.3]
What was it asked to do, and what must it have been given? Reconstruct the problem from its answer.
Hint: Every quantity in the result had to enter somewhere. Account for each one.
Answer:
In a real VAE Sigma is diagonal, so L = diag(sigma) and the transform is an elementwise mu + sigma * eps - no matrix multiply. Demonstrate that it lands the right mean and per-dim std:
Constraint
Discussion prompt
Run The multivariate-Gaussian toolkit with this step confiscated:
Distance: Mahalanobis (x-mu)^T Sigma^{-1} (x-mu), not Euclidean - it rescales by the spread
Is it still possible? If it is, say what takes its place and what it costs you. If it is not, say exactly what that step was providing that nothing else does.
Hint: A step you can drop for free was never load-bearing. If you cannot drop it, name the thing that goes wrong the moment it is gone.
Answer:
N(mu, Sigma) = front constant 1/((2 pi)^{d/2} |Sigma|^{1/2}) times exp(-1/2 * quadratic form)Sigma symmetric PSD (PD to invert), because a^T Sigma a = Var(a^T x) >= 0Sigma = axes, sqrt(eigenvalues) = 1-sigma semi-axes(x-mu)^T Sigma^{-1} (x-mu), not Euclidean - it rescales by the spreadSigma directly; conditionals shift the mean and shrink the variance; diagonal Sigma => independent (Gaussian only)x = mu + L z, L = chol(Sigma) - proven Cov = L L^T = Sigma; the VAE reparameterization trickPattern
N(mu, Sigma) = front constant 1/((2 pi)^{d/2} |Sigma|^{1/2}) times exp(-1/2 * quadratic form)Sigma symmetric PSD (PD to invert), because a^T Sigma a = Var(a^T x) >= 0Sigma = axes, sqrt(eigenvalues) = 1-sigma semi-axes(x-mu)^T Sigma^{-1} (x-mu), not Euclidean - it rescales by the spreadSigma directly; conditionals shift the mean and shrink the variance; diagonal Sigma => independent (Gaussian only)x = mu + L z, L = chol(Sigma) - proven Cov = L L^T = Sigma; the VAE reparameterization trickEdge cases
Discussion prompt
The multivariate-Gaussian toolkit works on the cases you have just seen. Push it to the edge: what is the most degenerate input it still handles — empty, zero, one item, everything equal — and what is the first case where it stops being true? Name the case, not just "it breaks".
Hint: Try the smallest legal input, then the largest, then the one where two things collide. Methods are specified at their edges; the middle takes care of itself.
Answer:
N(mu, Sigma) = front constant 1/((2 pi)^{d/2} |Sigma|^{1/2}) times exp(-1/2 * quadratic form)Sigma symmetric PSD (PD to invert), because a^T Sigma a = Var(a^T x) >= 0Sigma = axes, sqrt(eigenvalues) = 1-sigma semi-axes(x-mu)^T Sigma^{-1} (x-mu), not Euclidean - it rescales by the spreadSigma directly; conditionals shift the mean and shrink the variance; diagonal Sigma => independent (Gaussian only)x = mu + L z, L = chol(Sigma) - proven Cov = L L^T = Sigma; the VAE reparameterization trickElimination
Eliminate the wrong options
For N(mu, Sigma) with the standard density formula (using Sigma^{-1} and |Sigma|^{1/2}) to be well-defined and non-degenerate, Sigma must be:
3 of these 4 are wrong. Strike them one at a time, and say what rules each one out before you strike the next. The survivor is the answer.
Survives elimination: A
Why: Sigma is a covariance so it is symmetric PSD (a^T Sigma a = Var(a^T x) >= 0). The formula uses Sigma^{-1} and |Sigma|^{1/2}, which need every eigenvalue strictly > 0 - positive DEFINITE - or the inverse and the square root of the determinant fail.
Check
What must Sigma be for the standard density to make sense?
Check your understanding
For N(mu, Sigma) with the standard density formula (using Sigma^{-1} and |Sigma|^{1/2}) to be well-defined and non-degenerate, Sigma must be:
Answer: A
Why: Sigma is a covariance so it is symmetric PSD (a^T Sigma a = Var(a^T x) >= 0). The formula uses Sigma^{-1} and |Sigma|^{1/2}, which need every eigenvalue strictly > 0 - positive DEFINITE - or the inverse and the square root of the determinant fail.
Prediction
Predict first
For a correlated Gaussian, Mahalanobis distance is preferred over Euclidean because it:
Answer it in your own words, now, with nothing to choose from. The options are on the next slide — and picking the right one off a list is an easier skill than producing it.
Correct: measures distance in standard-deviation units along the covariance axes
Why: Sigma^{-1} rescales each direction by its variance and un-tilts the correlation, so Mahalanobis reports how many standard deviations away a point is along the ellipse axes - the correct notion of 'far' for correlated data. Points on the same contour ellipse are equidistant.
Check
Why rescale by Sigma^{-1} instead of using raw distance?
Check your understanding
For a correlated Gaussian, Mahalanobis distance is preferred over Euclidean because it:
Answer: A
Why: Sigma^{-1} rescales each direction by its variance and un-tilts the correlation, so Mahalanobis reports how many standard deviations away a point is along the ellipse axes - the correct notion of 'far' for correlated data. Points on the same contour ellipse are equidistant.
Prediction
Predict first
Let X ~ N(0, 1) and Y = X^2. Which statement is correct?
Answer it in your own words, now, with nothing to choose from. The options are on the next slide — and picking the right one off a list is an easier skill than producing it.
Correct: Cov(X, Y) = 0, but X and Y are NOT independent
Why: Cov(X, X^2) = E[X^3] - E[X]E[X^2] = 0 - 0 = 0 because odd moments of a symmetric X vanish (verified: sampled 0.0006). Yet Y is a deterministic function of X, so they are totally dependent. Zero covariance captures only LINEAR dependence.
Check
Does zero covariance always mean independence?
Check your understanding
Let X ~ N(0, 1) and Y = X^2. Which statement is correct?
Answer: A
Why: Cov(X, X^2) = E[X^3] - E[X]E[X^2] = 0 - 0 = 0 because odd moments of a symmetric X vanish (verified: sampled 0.0006). Yet Y is a deterministic function of X, so they are totally dependent. Zero covariance captures only LINEAR dependence.
Elimination
Eliminate the wrong options
To sample x ~ N(mu, Sigma) from standard normals z ~ N(0, I), you compute:
3 of these 4 are wrong. Strike them one at a time, and say what rules each one out before you strike the next. The survivor is the answer.
Survives elimination: A
Why: With Sigma = L L^T, Cov(mu + L z) = L Cov(z) L^T = L I L^T = L L^T = Sigma. The Cholesky factor L is a square root of Sigma that reshapes the isotropic z into the correct correlated ellipse, then mu shifts it into place.
Check
How do you turn standard normals into correlated ones?
Check your understanding
To sample x ~ N(mu, Sigma) from standard normals z ~ N(0, I), you compute:
Answer: A
Why: With Sigma = L L^T, Cov(mu + L z) = L Cov(z) L^T = L I L^T = L L^T = Sigma. The Cholesky factor L is a square root of Sigma that reshapes the isotropic z into the correct correlated ellipse, then mu shifts it into place.
Section
Project
Concept
Build a function that samples N(mu, Sigma) by Cholesky and prove it recovers the mean and covariance you asked for. You have derived every piece - now assemble it.
| # | requirement | tool |
|---|---|---|
| 1 | L = Cholesky of Sigma; check L L^T = Sigma | np.linalg.cholesky |
| 2 | x = mu + z @ L.T for z ~ N(0, I) | rng.normal |
| 3 | Verify sample mean and covariance match | np.cov |
Build rules: type every line yourself, draw enough samples (>= 10^5) for the covariance to converge, and remember np.cov expects variables in rows (so transpose with X.T).
Counterexample
Discussion prompt
Build a function that samples N(mu, Sigma) by Cholesky and prove it recovers the mean and covariance you asked for. You have derived every piece - now assemble it.
That is stated as though it always holds. Do one of two things: produce a case where it fails, or say precisely what rules such a case out. "It just does" is not on the menu.
Hint: Hunt at the extremes first — zero, one, negative, empty, equal. If every extreme survives, the reason they survive is the proof.
Answer:
Build rules: type every line yourself, draw enough samples (>= 10^5) for the covariance to converge, and remember np.cov expects variables in rows (so transpose with X.T).
Worked example
Your turn: factor Sigma = [[2, 0.8], [0.8, 1]] and confirm L L^T = Sigma. Predict whether L is lower- or upper-triangular before you print it.
Hint: np.linalg.cholesky returns a lower-triangular L. Use np.allclose(L @ L.T, Sigma) to check the reconstruction.
import numpy as np
Sigma = np.array([[2., 0.8], [0.8, 1.]])
L = np.linalg.cholesky(Sigma)
print(L.round(4))
print(np.allclose(L @ L.T, Sigma))| check | value |
|---|---|
| L | [[1.4142, 0], [0.5657, 0.8246]] |
| lower-triangular? | yes (top-right is 0) |
| L L^T == Sigma | True |
Worked example
Your turn: draw 200k standard normals and transform them with mu + z @ L.T. Predict the sample mean before you print (it should be mu).
Hint: z = rng.normal(size=(200000, 2)), then X = mu + z @ L.T. Use np.random.default_rng(0) so your numbers match these exactly.
import numpy as np
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
L = np.linalg.cholesky(Sigma)
rng = np.random.default_rng(0)
z = rng.normal(size=(200000, 2))
X = mu + z @ L.T
print(X.mean(0).round(3))| quantity | value |
|---|---|
| sample mean | [1.002, -1.000] |
| target mu | [1, -1] |
Worked example
Your turn: estimate the sample covariance and compare to Sigma. Predict: will the 0.8 off-diagonal correlation come back?
Hint: np.cov(X.T) returns the 2x2 covariance (variables in rows). It should match Sigma to ~2 decimals - including the off-diagonal.
import numpy as np
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
L = np.linalg.cholesky(Sigma)
X = mu + np.random.default_rng(0).normal(size=(200000, 2)) @ L.T
print(np.cov(X.T).round(3))
print(Sigma)| value | |
|---|---|
| sample cov | [[2.007, 0.798], [0.798, 0.999]] |
| target Sigma | [[2.0, 0.8], [0.8, 1.0]] |
Trade off
Comparison matrix
From Milestone 3 - recover Sigma: every row here is a choice with a cost. Fill the value column, then say which row you would actually pick and what you give up for it.
| value | |
|---|---|
| sample cov | [[2.007, 0.798], [0.798, 0.999]] |
| target Sigma | [[2.0, 0.8], [0.8, 1.0]] |
Worked example
Wrap it in a reusable sample_gaussian, then sanity-check with the mean Mahalanobis distance - which should land near d = 2 (each sample averages about d in squared Mahalanobis units, a chi-square(d) fact).
import numpy as np
def sample_gaussian(mu, Sigma, n, seed=0):
L = np.linalg.cholesky(Sigma) # Sigma = L L^T
z = np.random.default_rng(seed).normal(size=(n, len(mu)))
return mu + z @ L.T # x = mu + L z
mu = np.array([1., -1.])
Sigma = np.array([[2., 0.8], [0.8, 1.]])
X = sample_gaussian(mu, Sigma, 200000)
print('mean', X.mean(0).round(3))
print('cov', np.cov(X.T).round(3))
Sinv = np.linalg.inv(Sigma)
dM2 = np.einsum('ij,jk,ik->i', X - mu, Sinv, X - mu)
print('mean Mahalanobis^2', round(dM2.mean(), 3))| output | value (verified) |
|---|---|
| mean | [1.002, -1.0] |
| cov | [[2.007, 0.798], [0.798, 0.999]] |
| mean Mahalanobis^2 | 2.006 (approx d = 2) |
If your samples recover both mu and Sigma and the mean Mahalanobis lands near 2 - you have implemented the exact sampler a VAE uses to generate its latents.
Comparison
Comparison matrix
From The full program: refill the value (verified) column from what you know. The rest of the table is as it appeared.
| output | value (verified) |
|---|---|
| mean | [1.002, -1.0] |
| cov | [[2.007, 0.798], [0.798, 0.999]] |
| mean Mahalanobis^2 | 2.006 (approx d = 2) |
Concept
Out loud, slides closed: explain (1) why Sigma must be PSD from a^T Sigma a = Var(a^T x), (2) what Mahalanobis distance corrects for that Euclidean misses, and (3) why x = mu + L z produces covariance exactly Sigma.
Stretch (homework): derive the marginal of x_1 by integrating out x_2, prove diagonal Sigma => independent by factorizing the density, and threshold outliers at d_M^2 = 5.991 (the 95% point of chi-square with 2 dof). This distribution powers GMMs (Week 24) and the VAE latent space (Week 44).
Connect it up
Draw it
One page, no notation unless you need it: draw how these connect — From one bell curve to many · Why Sigma is symmetric PSD · The ellipse and Mahalanobis · The normalizing constant · Marginals, conditionals, independence · Sampling by Cholesky. Put an arrow wherever one of them is what makes another possible, and label the arrow with why.
Recap
N(mu, Sigma) density term by term and explain the quadratic form, the |Sigma|^{1/2}, and the (2 pi)^{d/2}Sigma is symmetric PSD (PD to invert) from a^T Sigma a = Var(a^T x) >= 0d_M^2 = 2.06 at [3, 0])Sigma, derive the conditional x_2 | x_1 = 3 (mean -0.2, var 0.68), and prove diagonal Sigma => independentx = mu + L z, prove Cov = L L^T = Sigma, recover it from 200k draws, and connect it to the VAE reparameterization trick| idea | the one thing to remember |
|---|---|
| Sigma | covariance: symmetric PSD, PD to invert |
| distance | Mahalanobis (x-mu)^T Sigma^{-1} (x-mu), not Euclidean |
| geometry | eigenvectors = ellipse axes, sqrt(eigenvalues) = semi-axes |
| conditioning | mean shifts, variance shrinks; diagonal => independent (Gaussian only) |
| sampling | x = mu + L z, L = chol(Sigma); Cov = L L^T = Sigma |
Want this taught 1-on-1? Alexander tutors Machine Learning — $55/session, free consultation.