Change of variables and the multivariate Gaussian

Before this lesson: lesson 04 (joint densities, the factorisation criterion, $\operatorname{Cov}(AX + b) = A\Sigma A^\top$) and lesson 03 (the normal and exponential distributions, moment generating functions). From other tracks: the change of variables in integrals (calculus 10), determinants (linear algebra 03) and the Cholesky factorisation (linear algebra 07). Time: about 85 minutes, plus the notebook.

---

Why this matters #

• It is the default model for noise and uncertainty. Affine maps, marginals and conditionals of a Gaussian are again Gaussian, with closed-form parameters.
• Gaussian conditioning under other names. The update step of the Kalman filter (weather forecasting, tracking), the prediction of a Gaussian process, and Bayesian linear regression are all one formula: the conditional of a joint Gaussian.
• Generative models. The VAE reparametrisation trick is the sampler $x = \mu + \sigma \odot z$. A diffusion model adds Gaussian noise with an affine map. A normalising flow trains with the change-of-variables formula.
• Numerics. In high dimension, $\det\Sigma$ underflows and $\Sigma^{-1}$ is inaccurate. Every serious implementation uses the Cholesky factor instead.

1. Predict first #

2. The picture #

Mass moves, density changes. Let $X$ be a random point in the plane with density $p_X$, and let $Y = g(X)$ for a smooth one-to-one map $g$. Take a tiny square of area $dA$ around a point $x$. It holds probability about $p_X(x)\,dA$. The map sends it to a small, slightly curved parallelogram around $y = g(x)$. Near $x$, the map looks like its Jacobian matrix $J_g(x)$ (the matrix of partial derivatives), so the image has area about $|\det J_g(x)|\,dA$ (calculus 10). The same mass now sits on that area:

Stretch the area, lower the density.

Interactive figure: loading…

From circles to ellipses. Let $z = (z_1, z_2)$ have independent $\mathcal N(0, 1)$ entries. Its density $\frac{1}{2\pi} e^{-(z_1^2 + z_2^2)/2}$ depends only on the length $\|z\|$, so its contour lines are circles and the cloud of samples is round. Now apply a matrix and a shift: $x = \mu + Lz$. Circles become ellipses, and by lesson 04

The contours of the density of $x$ are the ellipses $(x - \mu)^\top \Sigma^{-1} (x - \mu) = c^2$. The number $\sqrt{(x - \mu)^\top \Sigma^{-1}(x - \mu)}$ is the Mahalanobis distance of $x$ from $\mu$: the distance measured in standard deviations of the cloud, in each direction. On the ellipse with $c = 1$, every point is "one standard deviation" away, even when the Euclidean distances are very different.

Interactive figure: loading…

A marginal is a shadow; a conditional is a slice. Project the cloud onto the $x_1$ axis. The shadow is a 1-D Gaussian with mean $\mu_1$ and variance $\Sigma_{11}$: nothing to compute. Now fix $x_2$ and look at the density along that horizontal line, rescaled to integrate to 1. This slice is the conditional density of $x_1$ given $x_2$. It is again a 1-D Gaussian. With standard deviations $\sigma_1, \sigma_2$ and correlation $\rho$:

As $x_2$ moves, the centre of the slice moves on a straight line (the regression line of $x_1$ on $x_2$). The width of the slice stays the same. A strong correlation gives a narrow slice: knowing $x_2$ tells you a lot about $x_1$. Section 3 proves the general version.

Interactive figure: loading…

3. Definitions and main results #

Change of variables for a joint density #

Theorem. Let $X$ have density $p_X$ on an open set $U \subseteq \mathbb{R}^d$. Let $g: U \to V$ be a bijection onto an open set $V$, with $g$ and $g^{-1}$ continuously differentiable. Then $Y = g(X)$ has the density

where $J_{g^{-1}}(y)$ is the $d \times d$ matrix of partial derivatives of $g^{-1}$ at $y$.

Proof.

1. For any set $B \subseteq V$: $P(Y \in B) = P(X \in g^{-1}(B)) = \int_{g^{-1}(B)} p_X(x)\,dx$.
2. Substitute $x = g^{-1}(y)$ in this integral. The multivariable substitution rule says $dx = |\det J_{g^{-1}}(y)|\,dy$, so the integral equals $\int_B p_X(g^{-1}(y))\,|\det J_{g^{-1}}(y)|\,dy$.
3. This holds for every $B$, so the integrand is the density of $Y$.
4. For the second form: by the inverse function theorem $J_{g^{-1}}(y) = J_g(x)^{-1}$, so $|\det J_{g^{-1}}(y)| = 1/|\det J_g(x)|$. $\blacksquare$

It is enough that these conditions hold outside a set of probability 0. Two special cases:

• One dimension: $p_Y(y) = p_X(x)/|g'(x)|$, the rule of lesson 02.
• A linear map $Y = AX + b$ with invertible $A$: $p_Y(y) = p_X\big(A^{-1}(y - b)\big)/|\det A|$.

In log form, $\log p_Y(y) = \log p_X(x) - \log|\det J_g(x)|$. A normalising flow is a chain of invertible maps with cheap Jacobian determinants, trained by maximising exactly this log-density (generative models 04).

Box–Muller: two normals from two uniforms #

Theorem. Let $Z_1, Z_2$ be i.i.d. $\mathcal N(0, 1)$, and write $(Z_1, Z_2) = (R\cos\Theta, R\sin\Theta)$ with $R > 0$ and $\Theta \in [0, 2\pi)$. Then $R$ and $\Theta$ are independent, $\Theta$ is uniform on $[0, 2\pi)$, and $R^2 \sim \text{Exp}(1/2)$ (rate $1/2$, mean 2).

Proof.

1. The joint density of $Z = (Z_1, Z_2)$ is $p_Z(z) = \frac{1}{2\pi} e^{-(z_1^2 + z_2^2)/2}$.
2. The polar map $(r, \theta) \mapsto (r\cos\theta, r\sin\theta)$ is a smooth bijection from $(0, \infty) \times (0, 2\pi)$ onto the plane without the half-line $\{(z_1, 0) : z_1 \ge 0\}$, which has probability 0. It is $g^{-1}$ in the theorem above, and its Jacobian determinant is

3. So $p_{R,\Theta}(r, \theta) = p_Z(r\cos\theta, r\sin\theta) \cdot r = \dfrac{1}{2\pi} \cdot r e^{-r^2/2}$ for $r > 0$, $\theta \in (0, 2\pi)$.
4. This is a function of $\theta$ times a function of $r$, on a product set. By the factorisation criterion of lesson 04, $R$ and $\Theta$ are independent, $\Theta$ has the constant density $1/(2\pi)$, and $R$ has density $re^{-r^2/2}$ (the Rayleigh distribution).
5. Let $S = R^2$. Then $P(S > s) = P(R > \sqrt s) = \int_{\sqrt s}^\infty r e^{-r^2/2}\,dr = e^{-s/2}$: the tail of $\text{Exp}(1/2)$. $\blacksquare$

The sampler runs the proof backwards. Take independent $U_1, U_2$ uniform on $(0, 1)$. Inverse-CDF sampling (lesson 02) gives $S = -2\log(1 - U_1) \sim \text{Exp}(1/2)$, and $\Theta = 2\pi U_2$ is uniform. Then

are two independent $\mathcal N(0, 1)$ numbers. (Many texts write $-2\log U_1$; it has the same distribution, because $1 - U_1$ is also uniform.) Libraries use faster methods today (NumPy uses the ziggurat algorithm), but the Jacobian argument is the same.

The multivariate Gaussian density #

Definition. For $\mu \in \mathbb{R}^d$ and a symmetric positive definite (PD) $d \times d$ matrix $\Sigma$, the Gaussian (normal) distribution $\mathcal N(\mu, \Sigma)$ has the density

$\mu$ is the mean and $\Sigma$ the covariance. The inverse $\Lambda = \Sigma^{-1}$ is the precision matrix, and $(x - \mu)^\top \Sigma^{-1} (x - \mu)$ is the squared Mahalanobis distance.

Theorem (from noise to density). Let $z \sim \mathcal N(0, I_d)$ ($d$ independent standard normals). Let $\Sigma = LL^\top$, where $L$ is lower triangular with a positive diagonal (the Cholesky factor), and let $x = \mu + Lz$. Then $x$ has the density above.

Proof.

1. The density of $z$ is a product of $d$ standard normal densities: $p_z(z) = (2\pi)^{-d/2} e^{-z^\top z/2}$.
2. The map $z \mapsto \mu + Lz$ is invertible. Its inverse is $x \mapsto L^{-1}(x - \mu)$, a linear map with Jacobian matrix $L^{-1}$.
3. By the change-of-variables theorem, $p_x(x) = p_z\big(L^{-1}(x - \mu)\big)\,|\det L^{-1}|$.
4. $\det\Sigma = \det L \cdot \det L^\top = (\det L)^2$, so $|\det L^{-1}| = 1/|\det L| = (\det\Sigma)^{-1/2}$.
5. Put $z = L^{-1}(x - \mu)$ into the exponent, and write $L^{-\top} = (L^{-1})^\top$: $z^\top z = (x - \mu)^\top L^{-\top} L^{-1} (x - \mu) = (x - \mu)^\top \Sigma^{-1} (x - \mu)$, because $\Sigma^{-1} = (LL^\top)^{-1} = L^{-\top} L^{-1}$. $\blacksquare$

The shape of the contours. With the eigendecomposition $\Sigma = Q\,\operatorname{diag}(\lambda_1, \dots, \lambda_d)\,Q^\top$, the contour $(x - \mu)^\top \Sigma^{-1} (x - \mu) = c^2$ is an ellipsoid with axes along the eigenvectors $q_k$ and half-lengths $c\sqrt{\lambda_k}$. In 2-D the probability inside the ellipse of Mahalanobis radius $c$ is $1 - e^{-c^2/2}$: about $39.3\%$ for $c = 1$ and $86.5\%$ for $c = 2$ (much less than the 68% and 95% of one dimension).

Jointly Gaussian, in general. A random vector $x$ is jointly Gaussian if every linear combination $a^\top x$ is a 1-D Gaussian (a constant counts, as $\mathcal N(c, 0)$). Then its moment generating function is

because $t^\top x \sim \mathcal N(t^\top\mu, t^\top\Sigma t)$ and $\mathbb{E}[e^W] = e^{m + s^2/2}$ for $W \sim \mathcal N(m, s^2)$ (lesson 03). The MGF depends only on $\mu$ and $\Sigma$, and it fixes the distribution, so a Gaussian vector is fully described by its mean and covariance. $\Sigma$ may be singular (only PSD); then there is no density on $\mathbb{R}^d$, and the mass lies on a lower-dimensional affine subspace. Below, $\Sigma$ is PD unless we say otherwise.

Affine maps and marginals are Gaussian #

Theorem (affine maps). If $x \sim \mathcal N(\mu, \Sigma)$, $A \in \mathbb{R}^{m \times d}$ and $b \in \mathbb{R}^m$, then

Proof.

1. For $t \in \mathbb{R}^m$: $t^\top y = (A^\top t)^\top x + t^\top b$.
2. So $M_y(t) = e^{t^\top b}\,M_x(A^\top t) = \exp\!\big(t^\top (A\mu + b) + \tfrac12\,t^\top A\Sigma A^\top t\big)$.
3. This is the MGF of $\mathcal N(A\mu + b, A\Sigma A^\top)$, and the MGF fixes the distribution. $\blacksquare$

Corollary (marginals). Split $x = (x_a, x_b)$ into two blocks, and split $\mu$ and $\Sigma$ in the same way:

Then $x_a = \begin{pmatrix} I & 0 \end{pmatrix} x \sim \mathcal N(\mu_a, \Sigma_{aa})$: read off the blocks. The same works for any subset of coordinates: keep their entries of $\mu$ and their rows and columns of $\Sigma$.

Uncorrelated jointly Gaussian variables are independent #

Theorem. If $(x_a, x_b)$ is jointly Gaussian and $\Sigma_{ab} = 0$, then $x_a$ and $x_b$ are independent.

Proof (PD $\Sigma$).

1. $\Sigma = \begin{pmatrix} \Sigma_{aa} & 0 \\ 0 & \Sigma_{bb} \end{pmatrix}$ is block diagonal. So $\Sigma^{-1} = \begin{pmatrix} \Sigma_{aa}^{-1} & 0 \\ 0 & \Sigma_{bb}^{-1} \end{pmatrix}$ and $\det\Sigma = \det\Sigma_{aa} \det\Sigma_{bb}$.
2. The quadratic form splits into a part with $x_a$ and a part with $x_b$: $(x - \mu)^\top \Sigma^{-1} (x - \mu) = (x_a - \mu_a)^\top \Sigma_{aa}^{-1} (x_a - \mu_a) + (x_b - \mu_b)^\top \Sigma_{bb}^{-1} (x_b - \mu_b)$.
3. So $p(x) = \mathcal N(x_a; \mu_a, \Sigma_{aa}) \cdot \mathcal N(x_b; \mu_b, \Sigma_{bb})$, a function of $x_a$ times a function of $x_b$, on all of $\mathbb{R}^d$. By the factorisation criterion, $x_a$ and $x_b$ are independent. $\blacksquare$

The word "jointly" is essential. Gaussian marginals are not enough. Let $X \sim \mathcal N(0, 1)$, let $W = \pm 1$ with probability $1/2$ each, independent of $X$, and let $Y = WX$.

1. $Y \sim \mathcal N(0, 1)$, because the density of $X$ is symmetric: flipping the sign does not change it.
2. $\operatorname{Cov}(X, Y) = \mathbb{E}[WX^2] = \mathbb{E}[W]\,\mathbb{E}[X^2] = 0 \cdot 1 = 0$: uncorrelated.
3. But $X + Y = 0$ whenever $W = -1$, so $P(X + Y = 0) = 1/2$. A Gaussian has no point with positive probability (except the constant case), so $X + Y$ is not Gaussian, and $(X, Y)$ is not jointly Gaussian.
4. $X$ and $Y$ are dependent: $|Y| = |X|$ always.

Interactive figure: loading…

Cholesky in practice: sampling and a stable log-density #

Let $\Sigma = LL^\top$ be the Cholesky factorisation (linear algebra 07). It costs about $d^3/3$ floating-point operations, once.

• Sampling. Draw $z \sim \mathcal N(0, I_d)$ and set $x = \mu + Lz$. By the affine-map theorem, $x \sim \mathcal N(\mu, LL^\top) = \mathcal N(\mu, \Sigma)$. Each sample costs about $d^2$ operations. With $n$ samples as the rows of an $n \times d$ matrix $Z$, the samples are the rows of $\mathbf 1 \mu^\top + Z L^\top$. Any $B$ with $BB^\top = \Sigma$ works (for example $Q\,\operatorname{diag}(\sqrt{\lambda_k})$); Cholesky is the cheapest.
• Log-density. $\log p(x) = -\tfrac12\big[d\log(2\pi) + \log\det\Sigma + (x - \mu)^\top \Sigma^{-1}(x - \mu)\big]$. Two identities make it stable:
  1. The determinant of a triangular matrix is the product of its diagonal. So $\det\Sigma = (\det L)^2 = \prod_i L_{ii}^2$ and $\log\det\Sigma = 2\sum_i \log L_{ii}$.
  2. Solve $Lu = x - \mu$ by forward substitution (about $d^2$ operations). Then $(x - \mu)^\top \Sigma^{-1} (x - \mu) = (x - \mu)^\top L^{-\top} L^{-1} (x - \mu) = \|u\|^2$.
• Why not det and inv. $\det\Sigma$ is a product of $d$ numbers, so it under- or overflows in float64 while $\log\det\Sigma$ is an ordinary number. An explicit $\Sigma^{-1}$ costs more and loses accuracy when $\Sigma$ is ill-conditioned; triangular solves are backward stable.
• A free check. $u = L^{-1}(x - \mu) \sim \mathcal N(0, I_d)$, so for real samples the squared Mahalanobis distance $\|u\|^2$ follows a $\chi^2_d$ distribution (mean $d$).
• When Cholesky fails: jitter. A covariance with eigenvalues close to 0 (for example a smooth Gaussian-process kernel on a dense grid) is PD in exact arithmetic, but rounding errors make it numerically indefinite, and Cholesky stops with an error. The standard fix is a small jitter: factor $\Sigma + \varepsilon I$, with $\varepsilon$ around $10^{-8}$ times the size of the variances.

Conditioning: the Schur complement #

For the split $\Sigma$ with $\Sigma_{bb}$ invertible, the Schur complement of $\Sigma_{bb}$ is

It is what remains of $\Sigma_{aa}$ after one step of block Gaussian elimination removes the $b$ block.

First, the 2-D case with numbers. Take $a = \{1\}$, $b = \{2\}$, $\Sigma = \begin{pmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}$. Then $K = \Sigma_{ab}\Sigma_{bb}^{-1} = \rho\sigma_1/\sigma_2$ and $S = \sigma_1^2 - (\rho\sigma_1\sigma_2)^2/\sigma_2^2 = \sigma_1^2(1 - \rho^2)$, which is the formula of section 2. With $\sigma_1 = \sigma_2 = 1$ and $\rho = 0.8$: $K = 0.8$ and $S = 1 - 0.64 = 0.36$.

Theorem (Gaussian conditioning). Let $x = (x_a, x_b) \sim \mathcal N(\mu, \Sigma)$ with $\Sigma_{bb}$ invertible, and let $K = \Sigma_{ab}\Sigma_{bb}^{-1}$. Then $x_a \mid x_b \sim \mathcal N(\mu_{a \mid b}, \Sigma_{a \mid b})$ with

In precision form, with $\Lambda = \Sigma^{-1}$ split into the same blocks:

Proof: remove the part of $x_a$ that $x_b$ explains

1. Let $z = x_a - Kx_b$. In matrix form $(z, x_b) = Tx$ with $T = \begin{pmatrix} I & -K \\ 0 & I \end{pmatrix}$, so by the affine-map theorem $(z, x_b)$ is jointly Gaussian.
2. By bilinearity, $\operatorname{Cov}(z, x_b) = \Sigma_{ab} - K\Sigma_{bb} = \Sigma_{ab} - \Sigma_{ab} = 0$. By the previous theorem, $z$ and $x_b$ are independent.
3. $\mathbb{E}[z] = \mu_a - K\mu_b$, and $\operatorname{Cov}(z) = \Sigma_{aa} - K\Sigma_{ba} - \Sigma_{ab}K^\top + K\Sigma_{bb}K^\top$. Each of the last three terms equals $\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}$ (use $\Sigma_{bb}^\top = \Sigma_{bb}$). So $\operatorname{Cov}(z) = S$, and $T\Sigma T^\top = \begin{pmatrix} S & 0 \\ 0 & \Sigma_{bb} \end{pmatrix}$. In particular $S$ is PD.
4. Write $x_a = z + Kx_b$. Given $x_b$, the term $Kx_b$ is a constant, and $z$ keeps its distribution because it is independent of $x_b$. So $x_a \mid x_b \sim \mathcal N(\mu_a - K\mu_b + Kx_b,\ S)$: the covariance form.
5. Precision form. $\det T = 1$, so $T$ is invertible, and step 3 gives $\Sigma = T^{-1}\begin{pmatrix} S & 0 \\ 0 & \Sigma_{bb} \end{pmatrix}T^{-\top}$. Invert both sides:

$$\Lambda = T^\top \begin{pmatrix} S^{-1} & 0 \\ 0 & \Sigma_{bb}^{-1} \end{pmatrix} T = \begin{pmatrix} S^{-1} & -S^{-1}K \\ -K^\top S^{-1} & \Sigma_{bb}^{-1} + K^\top S^{-1} K \end{pmatrix}.$$

So $\Lambda_{aa} = S^{-1}$ and $\Lambda_{ab} = -S^{-1}K$, that is, $S = \Lambda_{aa}^{-1}$ and $K = -\Lambda_{aa}^{-1}\Lambda_{ab}$. Put these into step 4. $\blacksquare$

Second route: completing the square

1. As a function of $x_a$, $\log p(x_a \mid x_b) = \log p(x_a, x_b) - \log p(x_b) = -\tfrac12 (x - \mu)^\top \Lambda (x - \mu) + \text{const}$.
2. Expand in blocks and keep only the terms with $x_a$: $-\tfrac12\,x_a^\top \Lambda_{aa} x_a + x_a^\top \big[\Lambda_{aa}\mu_a - \Lambda_{ab}(x_b - \mu_b)\big] + \text{const}$.
3. $\mathcal N(m, P^{-1})$ has log-density $-\tfrac12\,x_a^\top P x_a + x_a^\top P m + \text{const}$. Matching the quadratic terms gives $P = \Lambda_{aa}$; matching the linear terms gives $Pm = \Lambda_{aa}\mu_a - \Lambda_{ab}(x_b - \mu_b)$. This is the precision form.

This route is the one used for Bayesian posteriors (lesson 10) and the Laplace approximation.

Zeros of the precision matrix. For $i \ne j$: $\Lambda_{ij} = 0$ if and only if $x_i$ and $x_j$ are conditionally independent given all other coordinates. Take $a = \{i, j\}$: the conditional covariance of $(x_i, x_j)$ is $\Lambda_{aa}^{-1}$, whose off-diagonal entry is $-\Lambda_{ij}/(\Lambda_{ii}\Lambda_{jj} - \Lambda_{ij}^2)$; it is 0 exactly when $\Lambda_{ij} = 0$. With $a = \{i\}$: $\operatorname{Var}(x_i \mid \text{rest}) = 1/\Lambda_{ii}$. A zero in $\Sigma$ means marginal independence; a zero in $\Lambda$ means conditional independence. Gaussian graphical models (graphical lasso) and Gaussian Markov random fields are built on the zeros of $\Lambda$.

Rule of thumb: marginals are easy in covariance form (read off $\Sigma_{aa}$); conditionals are easy in precision form (read off $\Lambda_{aa}$).

What the theorem says.

• The conditional mean is affine in $x_b$, and $K$ is the matrix of regression coefficients of $x_a$ on $x_b$. For Gaussians, the best predictor $\mathbb{E}[x_a \mid x_b]$ (lesson 04) is linear.
• The conditional covariance $S$ does not depend on the observed value $x_b$.
• $S \preceq \Sigma_{aa}$, because $\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}$ is PSD: for a Gaussian, conditioning never increases a variance, whatever you observe. (For other distributions, only the average $\mathbb{E}[\operatorname{Var}(X \mid Y)]$ is at most $\operatorname{Var}(X)$, by the law of total variance.)
• $\det\Sigma = \det\Sigma_{bb} \cdot \det S$, the determinant version of $p(x) = p(x_b)\,p(x_a \mid x_b)$.

Interactive figure: loading…

Answers to section 1. (1) $\mathbb{E}[x_1 \mid x_2 = 2] = 0.8 \cdot 2 = 1.6$, and $\operatorname{Var}(x_1 \mid x_2) = 1 - 0.8^2 = 0.36$. For $x_2 = -5$ the mean is $-4$, but the variance is the same $0.36$: the conditional covariance does not depend on the observed value. (2) $\det\Sigma = 0.1^{1000} = 10^{-1000}$ underflows, so det returns 0.0 (and its log would be $-\infty$). The true value is $\log\det\Sigma = 1000\log 0.1 \approx -2302.585$; Cholesky gives it as $2\sum_i \log\sqrt{0.1}$.

4. Worked example #

(a) Box–Muller by hand. Take $U_1 = 1 - e^{-1} \approx 0.632$ and $U_2 = 0.125$.

1. $S = -2\log(1 - U_1) = -2\log e^{-1} = 2$, so $R = \sqrt 2$.
2. $\Theta = 2\pi \cdot 0.125 = \pi/4$.
3. $Z_1 = \sqrt 2 \cos(\pi/4) = \sqrt 2 \cdot \tfrac{\sqrt 2}{2} = 1$ and $Z_2 = \sqrt 2 \sin(\pi/4) = 1$.

(b) Fusing two thermometers. The true temperature $\theta$ at a site has the prior $\mathcal N(20, 4)$ (in °C). Thermometer 1 reads $y_1 = \theta + e_1$ with error variance 1, and thermometer 2 reads $y_2 = \theta + e_2$ with error variance 4. The errors are independent of each other and of $\theta$. So $x = (\theta, y_1, y_2)$ is an affine map of the independent Gaussians $(\theta, e_1, e_2)$, hence jointly Gaussian, with

For example, $\operatorname{Cov}(y_1, y_2) = \operatorname{Var}(\theta) = 4$ (lesson 04, the two-thermometer example), and $\operatorname{Var}(y_2) = 4 + 4 = 8$.

Cholesky, column by column.

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

$L$ is the model itself: $\theta = 20 + 2z_1$, $y_1 = \theta + z_2$, $y_2 = \theta + 2z_3$. And $\det\Sigma = (2 \cdot 1 \cdot 2)^2 = 16$.

One sample by hand. $z = (0.5, -1, 1)$ gives $x = \mu + Lz = (20 + 1,\ 20 + 1 - 1,\ 20 + 1 + 2) = (21, 20, 23)$.

Log-density at $x = (21, 23, 17)$. Here $x - \mu = (1, 3, -3)$. Forward substitution in $Lu = x - \mu$: $u_1 = 1/2 = 0.5$, then $u_2 = (3 - 2 \cdot 0.5)/1 = 2$, then $u_3 = (-3 - 2 \cdot 0.5 - 0 \cdot 2)/2 = -2$. So $\|u\|^2 = 0.25 + 4 + 4 = 8.25$ and

$8.25$ is above the 95% quantile of $\chi^2_3$ (7.815), so this point is unusual. The reason: the thermometers disagree by 6, while $y_1 - y_2$ has variance $5 + 8 - 2 \cdot 4 = 5$. The gap is $6/\sqrt 5 \approx 2.7$ standard deviations.

Marginal. $(y_1, y_2) \sim \mathcal N\Big(\begin{pmatrix}20\\20\end{pmatrix}, \begin{pmatrix} 5 & 4 \\ 4 & 8 \end{pmatrix}\Big)$: the lower-right block.

Conditional, covariance form. Observe $y_1 = 23$ and $y_2 = 17$. Take $a = \{\theta\}$ and $b = \{y_1, y_2\}$: $\Sigma_{ab} = (4, 4)$, $\Sigma_{bb} = \begin{pmatrix} 5 & 4 \\ 4 & 8 \end{pmatrix}$, $\det\Sigma_{bb} = 40 - 16 = 24$, and $\Sigma_{bb}^{-1} = \frac{1}{24}\begin{pmatrix} 8 & -4 \\ -4 & 5 \end{pmatrix}$.

1. $K = \Sigma_{ab}\Sigma_{bb}^{-1} = \frac{1}{24}\,(4 \cdot 8 - 4 \cdot 4,\ -4 \cdot 4 + 4 \cdot 5) = \frac{1}{24}(16, 4) = (2/3,\ 1/6)$.
2. $\mu_{a \mid b} = 20 + \tfrac23 \cdot (23 - 20) + \tfrac16 \cdot (17 - 20) = 20 + 2 - 0.5 = 21.5$.
3. $\Sigma_{a \mid b} = 4 - K\Sigma_{ba} = 4 - \big(\tfrac23 \cdot 4 + \tfrac16 \cdot 4\big) = 4 - \tfrac{10}{3} = \tfrac23$: a standard deviation of $0.816$, down from $2$.
4. Check: $\det\Sigma_{bb} \cdot \Sigma_{a \mid b} = 24 \cdot \tfrac23 = 16 = \det\Sigma$. ✓

Conditional, precision form.

1. $\Sigma_{a \mid b} = 1/\Lambda_{\theta\theta} = 1/1.5 = 2/3$. ✓ And $1.5 = \tfrac14 + 1 + \tfrac14$: the prior precision plus the two thermometer precisions. Precisions add.
2. $\mu_{a \mid b} = 20 - \tfrac{1}{1.5}\big((-1) \cdot 3 + (-0.25) \cdot (-3)\big) = 20 + \tfrac{2.25}{1.5} = 21.5$. ✓ Equivalently, a precision-weighted average: $(20/4 + 23/1 + 17/4)/1.5 = 21.5$. The prior mean, $y_1$ and $y_2$ get the weights $1/6$, $2/3$ and $1/6$: their shares of the total precision.
3. $\Lambda_{y_1 y_2} = 0$: the two readings are conditionally independent given $\theta$ (their errors are independent), although $\operatorname{Cov}(y_1, y_2) = 4 \ne 0$.

5. In code #

SciPy and NumPy check the hand computations. The notebook for this lesson has these cells and the exercises.

import numpy as np
from scipy import stats

mu = np.array([20.0, 20.0, 20.0])                 # (θ, y1, y2)
Sigma = np.array([[4.0, 4.0, 4.0],
                  [4.0, 5.0, 4.0],
                  [4.0, 4.0, 8.0]])
L = np.linalg.cholesky(Sigma)
print(L)                                          # [[2. 0. 0.] [2. 1. 0.] [2. 0. 2.]]
print(round(2 * np.log(np.diag(L)).sum(), 4))     # 2.7726 = log det Σ = 4 log 2

x = np.array([21.0, 23.0, 17.0])
u = np.linalg.solve(L, x - mu)                    # forward substitution (use scipy.linalg.solve_triangular for speed)
logp = -0.5 * (3 * np.log(2 * np.pi) + 4 * np.log(2) + u @ u)
print(u, round(logp, 4))                          # [ 0.5  2.  -2. ] -8.2681
print(round(stats.multivariate_normal(mu, Sigma).logpdf(x), 4))   # -8.2681

# conditioning θ on (y1, y2) = (23, 17), covariance form
Sab, Sbb = Sigma[0, 1:], Sigma[1:, 1:]
K = np.linalg.solve(Sbb, Sab)                     # = Σ_ab Σ_bb^{-1} (Σ_bb is symmetric)
print(K.round(4), 20 + K @ (x[1:] - 20), round(4 - K @ Sab, 4))   # [0.6667 0.1667] 21.5 0.6667

# det underflows, log det does not
C = 0.1 * np.eye(1000)
print(np.linalg.det(C), round(np.linalg.slogdet(C)[1], 3))       # 0.0 -2302.585

# Gaussian marginals, uncorrelated, but not jointly Gaussian: Y = W X
rng = np.random.default_rng(1)
X1 = rng.standard_normal(200_000)
Y1 = rng.choice([-1.0, 1.0], size=X1.size) * X1
print(round(np.corrcoef(X1, Y1)[0, 1], 2), round(np.mean(X1 + Y1 == 0), 2))  # about 0.0 and 0.5

A least-squares fit on samples from the model finds $K$ and the Schur complement, because the conditional mean is affine:

rng = np.random.default_rng(0)
n = 1_000_000
theta = 20 + 2 * rng.standard_normal(n)           # θ ~ N(20, 4)
y1 = theta + rng.standard_normal(n)               # error variance 1
y2 = theta + 2 * rng.standard_normal(n)           # error variance 4
D = np.column_stack([np.ones(n), y1, y2])
coef, *_ = np.linalg.lstsq(D, theta, rcond=None)
print(coef[1:].round(3), (theta - D @ coef).var().round(3))      # about [0.667 0.167] 0.667

6. Where it shows up #

• Kalman filter and data assimilation. Forecast $x \sim \mathcal N(m, P)$, observation $y = Hx + e$ with $e \sim \mathcal N(0, R)$. Then $\operatorname{Cov}(x, y) = PH^\top$ and $\operatorname{Cov}(y) = HPH^\top + R$, and the conditioning theorem gives the update: the Kalman gain $K = PH^\top(HPH^\top + R)^{-1}$, the mean $m + K(y - Hm)$ and the covariance $P - KHP$, a Schur complement (ML for science 06).
• Gaussian processes. A GP prediction conditions the joint Gaussian of the training targets and the new function values. Implementations factor $k(X, X) + \sigma^2 I = LL^\top$ once, use triangular solves and jitter, and compute $\log\det$ as $2\sum_i \log L_{ii}$ (Rasmussen and Williams, Algorithm 2.1).
• VAEs and diffusion models. The VAE samples its latent code as $\mu_\phi(x) + \sigma_\phi(x) \odot \varepsilon$, which is $\mu + Lz$ with a diagonal $L$; gradients flow through $\mu$ and $\sigma$. The diffusion forward process $x_t = \sqrt{\bar\alpha_t}\,x_0 + \sqrt{1 - \bar\alpha_t}\,\varepsilon$ is an affine map of Gaussian noise, and the DDPM posterior $q(x_{t-1} \mid x_t, x_0)$ is Gaussian conditioning.
• Normalising flows maximise $\log p_X(g^{-1}(y)) + \log|\det J_{g^{-1}}(y)|$, the change-of-variables formula, with maps whose Jacobian determinant is cheap (triangular Jacobians) (generative models 04).
• Graphical models and random fields. The zeros of the precision matrix encode conditional independence: the graphical lasso estimates a sparse $\Lambda$, and Gaussian Markov random fields on large grids (weather, images) use a sparse precision with a sparse Cholesky factor.
• Anomaly detection. The Mahalanobis distance to the mean of a class, compared with a $\chi^2_d$ quantile, is a simple and strong out-of-distribution score, also on the features of a neural network.

7. Common mistakes #

• Multiplying by $|\det J_g|$ instead of dividing. Where the map stretches, the density goes down: $p_Y(y) = p_X(x)/|\det J_g(x)|$.
• "Gaussian marginals, so jointly Gaussian." False (the random-sign example). So "Gaussian and uncorrelated, so independent" is also false, unless the vector is jointly Gaussian.
• The wrong transpose in the sampler. $\mu + L^\top z$ has covariance $L^\top L \ne \Sigma$. With samples as rows, write mu + Z @ L.T.
• Mixing up the blocks of $\Sigma$ and $\Lambda$. $\Sigma_{aa}$ is the marginal covariance; $\Lambda_{aa}^{-1}$ is the conditional covariance. In general $(\Sigma^{-1})_{aa} \ne (\Sigma_{aa})^{-1}$.
• Computing det(Sigma) and inv(Sigma). They under- or overflow and lose accuracy. Use Cholesky; if it fails, add jitter rather than switching to a pseudo-inverse.
• Thinking the conditional variance depends on the observed value, or that conditioning can increase a Gaussian variance.

8. Check yourself #

1. $(x_1, x_2)$ is Gaussian with means 0, $\sigma_1 = 2$, $\sigma_2 = 3$ and $\rho = 0.5$. Find $\mathbb{E}[x_1 \mid x_2 = 3]$ and $\operatorname{Var}(x_1 \mid x_2)$.
   Answer

   $\Sigma_{12} = \rho\sigma_1\sigma_2 = 3$, so $K = 3/9 = 1/3$ and $\mathbb{E}[x_1 \mid x_2 = 3] = 1$. $\operatorname{Var}(x_1 \mid x_2) = 4 - 3^2/9 = 3$, which equals $\sigma_1^2(1 - \rho^2) = 4 \cdot 0.75$.
2. Find the Cholesky factor and the log-determinant of $\begin{pmatrix} 4 & 2 \\ 2 & 10 \end{pmatrix}$.
   Answer

   $L_{11} = 2$, $L_{21} = 2/2 = 1$, $L_{22} = \sqrt{10 - 1} = 3$. $\log\det = 2(\log 2 + \log 3) = 2\log 6 \approx 3.584$. Check: $\det = 40 - 4 = 36 = 6^2$.
3. $x \sim \mathcal N\Big(0, \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\Big)$. Give the joint distribution of $(x_1 + x_2,\ x_1 - x_2)$. Are the two independent?
   Answer

   With $A = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$: $A\Sigma A^\top = \operatorname{diag}(6, 2)$, so $\mathcal N(0, \operatorname{diag}(6, 2))$. They are jointly Gaussian (an affine map of a Gaussian) and uncorrelated, hence independent.
4. $X$ is uniform on the unit square and $Y = AX$ with $A = \begin{pmatrix} 2 & 1 \\ 0 & 3 \end{pmatrix}$. What is the density of $Y$?
   Answer

   $p_Y(y) = p_X(A^{-1}y)/|\det A| = 1/6$ on the parallelogram $A[0, 1]^2$ (area 6), and 0 elsewhere.
5. A Gaussian has the precision $\Lambda = \begin{pmatrix} 1 & 0 & 0.5 \\ 0 & 2 & 1 \\ 0.5 & 1 & 2 \end{pmatrix}$. Which pair is conditionally independent given the third variable? Is that pair marginally independent? What is $\operatorname{Var}(x_3 \mid x_1, x_2)$?
   Answer

   $x_1$ and $x_2$ are conditionally independent given $x_3$, because $\Lambda_{12} = 0$. They are not marginally independent: $\Sigma = \Lambda^{-1}$ has $\Sigma_{12} = 0.2 \ne 0$. $\operatorname{Var}(x_3 \mid x_1, x_2) = 1/\Lambda_{33} = 0.5$, against the marginal $\Sigma_{33} = 0.8$.
6. You have the Cholesky factor $R$ of the precision, $\Lambda = RR^\top$ (common for sparse models). How do you sample $\mathcal N(\mu, \Lambda^{-1})$ without forming $\Sigma$?
   Answer

   Draw $z \sim \mathcal N(0, I)$, solve the triangular system $R^\top v = z$ by back substitution, and set $x = \mu + v$. Then $\operatorname{Cov}(v) = R^{-\top}R^{-1} = (RR^\top)^{-1} = \Lambda^{-1} = \Sigma$.

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

1. State the change-of-variables formula for an invertible smooth map $Y = g(X)$ in $d$ dimensions.

   Show answer

2. Why is the density divided by $|\det J_g|$ and not multiplied?

   Show answer

3. In Box–Muller, what are the distributions of $\Theta$ and $R^2$ for two i.i.d. standard normals?

   Show answer

4. Write the density of $\mathcal N(\mu, \Sigma)$ in $d$ dimensions.

   Show answer

5. How do you sample $\mathcal N(\mu, \Sigma)$, and why does it work?

   Show answer

6. How do you compute $\log\det\Sigma$ from the Cholesky factor $L$?

   Show answer

7. What is the marginal of $x_a$ when $(x_a, x_b) \sim \mathcal N(\mu, \Sigma)$?

   Show answer

8. Give the conditional mean and covariance of $x_a$ given $x_b$ for a joint Gaussian.

   Show answer

9. For a 2-D Gaussian with $\rho = 0.8$ and unit variances, what is $\operatorname{Var}(x_1 \mid x_2)$?

   Show answer

10. When do zero covariance and independence mean the same thing?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "How do you sample from $\mathcal N(\mu, \Sigma)$?" — Cholesky $\Sigma = LL^\top$ and $x = \mu + Lz$, since $\operatorname{Cov}(Lz) = LL^\top$. Cost: $d^3/3$ once, $d^2$ per sample. If $\Sigma$ is only PSD: the eigendecomposition $Q\,\operatorname{diag}(\sqrt{\lambda_k})$, a pivoted Cholesky, or jitter.
• "Derive the conditional of a bivariate Gaussian." — Mean $\mu_1 + \rho\frac{\sigma_1}{\sigma_2}(x_2 - \mu_2)$, variance $\sigma_1^2(1 - \rho^2)$; the proof uses $z = x_1 - Kx_2$, which is uncorrelated with $x_2$ and therefore independent of it.
• "Are uncorrelated Gaussian variables independent?" — Yes if they are jointly Gaussian; the random-sign example shows that Gaussian marginals are not enough.
• "Why compute a Gaussian log-likelihood with Cholesky?" — $\det$ under- or overflows and an inverse is slow and inaccurate; $2\sum_i \log L_{ii}$ and one triangular solve are stable.
• "What does a zero in the precision matrix mean? Where does the Kalman gain come from?" — Conditional independence given all other variables; the gain is $K = \Sigma_{xy}\Sigma_{yy}^{-1}$ from Gaussian conditioning.
• "How does a normalising flow compute $\log p(x)$?" — Map $x$ to the base variable $z$, then $\log p(x) = \log p_Z(z) + \log|\det J_{g^{-1}}(x)|$, with a Jacobian that is triangular so the determinant is cheap.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. A Box–Muller sampler, tested with a Kolmogorov–Smirnov test against $\mathcal N(0, 1)$.
2. mvn_sample with the Cholesky factor; check the sample mean and covariance.
3. mvn_logpdf with Cholesky, compared with scipy.stats.multivariate_normal, also for $d = 500$.
4. mvn_condition with the Schur complement, compared with a least-squares regression on samples.
5. Jitter: the smallest jitter that lets Cholesky factor a smooth kernel matrix on a dense grid.

11. Watch and read #

Lecture 22: Transformations and Convolutions | Statistics 110Harvard University · Joe Blitzstein, Statistics 110 · 48 min

The change-of-variables formula with the Jacobian, in one and in two dimensions.

Lecture 30: Chi-Square, Student-t, Multivariate Normal | Statistics 110Harvard University · Joe Blitzstein, Statistics 110 · 47 min

The second half defines the multivariate normal through linear combinations and shows why uncorrelated implies independent for it.

A pretty reason why Gaussian + Gaussian = Gaussian3Blue1Brown · 13 min

The rotational symmetry of the 2-D standard Gaussian, the same fact that makes Box–Muller work.

Multivariate Normal (Gaussian) Distribution ExplainedDataMListic · 7 min

A short visual summary of the density, the role of the covariance matrix and the ellipse contours.

• Christopher Bishop, Pattern Recognition and Machine Learning (free PDF from Microsoft Research): §2.3 on the geometry of the Gaussian, §2.3.1 conditional Gaussians (completing the square and the Schur complement), §2.3.2 marginals, §2.3.3 linear Gaussian models.
• Carl Rasmussen and Christopher Williams, Gaussian Processes for Machine Learning, free chapters: Appendix A.2 (Gaussian identities and sampling with Cholesky), A.3 (inverse of a partitioned matrix), A.4 (Cholesky and $\log\det$).
• Kevin Murphy, Probabilistic Machine Learning: An Introduction, free PDF: §3.2 (multivariate Gaussian, Mahalanobis distance, marginals and conditionals) and §3.3 (linear Gaussian systems).
• Joseph Blitzstein and Jessica Hwang, Introduction to Probability, 2nd edition: §8.1 (change of variables, with Box–Muller as an example).