The multivariate Gaussian

A joint distribution of $d$ variables is in general a function of $d$ arguments that no table can hold. One family is different. It is described completely by a mean vector and a covariance matrix, it stays in the family under linear maps, marginalization and conditioning, and all of these are a few lines of linear algebra. That family, the multivariate Gaussian, is the noise model of linear regression, the prior of Gaussian processes, the state of a Kalman filter and the latent space of most generative models. This note builds it from the one-dimensional normal and the covariance of the previous note.

Notation. $X \in \mathbb R^d$ is a random vector with mean $\mu = E[X]$ and covariance matrix $\Sigma$, $\Sigma_{ij} = \operatorname{Cov}(X_i, X_j)$, symmetric and positive definite. $X \sim N(\mu, \Sigma)$; $Z \sim N(0, I)$ is the standard normal vector. $\lvert \Sigma \rvert$ is the determinant, $A^\top$ the transpose. The note uses a little linear algebra: matrix products, inverses and determinants, the quadratic form $x^\top A x$, and the eigenvectors of a symmetric matrix.

The density#

$X \sim N(\mu, \Sigma)$ has the density

$$f(x) = \frac{1}{(2\pi)^{d/2} \lvert \Sigma \rvert^{1/2}} \exp\Big(-\frac12 (x - \mu)^\top \Sigma^{-1} (x - \mu)\Big) .$$

For $d = 1$ it is the normal density of the continuous distributions note, with $\Sigma = \sigma^2$. The quadratic form in the exponent, $\Delta^2(x) = (x - \mu)^\top \Sigma^{-1} (x - \mu)$, is the squared Mahalanobis distance: the distance from the mean measured in standard deviations along each direction. The density is constant on the sets where $\Delta$ is constant, which are ellipsoids.

If $\Sigma$ is diagonal, the quadratic form is a sum $\sum_i (x_i - \mu_i)^2 / \sigma_i^2$, the exponential of a sum is a product, and the density factorizes into one normal density per coordinate. The standard normal vector $Z \sim N(0, I)$ is simply $d$ independent $N(0, 1)$ coordinates.

Covariance as geometry#

A symmetric positive definite matrix has an orthonormal basis of eigenvectors $q_1, \dots, q_d$ with positive eigenvalues $\lambda_1, \dots, \lambda_d$: $\Sigma = Q \Lambda Q^\top$. In the coordinates of that basis the covariance is diagonal, so the density factorizes there, with variance $\lambda_i$ along $q_i$. The contours are ellipses with axes along the eigenvectors and half-lengths $k\sqrt{\lambda_i}$ at Mahalanobis distance $k$. For

$$\Sigma = \begin{pmatrix} 2 & 1.2 \\ 1.2 & 1 \end{pmatrix}, \qquad \lambda = 2.8,\ 0.2, \qquad \rho = \frac{1.2}{\sqrt{2 \cdot 1}} = 0.85 ,$$

the long axis points along $(0.83, 0.55)$:

Figure

A covariance matrix as an ellipse

The long axis carries $2.8 / 3 = 93\%$ of the total variance $\operatorname{tr} \Sigma = 3$, and projecting onto it gives the first principal component. The ellipses also correct a common intuition. In one dimension, one standard deviation around the mean holds 68% of the probability; in two dimensions the ellipse of Mahalanobis distance 1 holds only $1 - e^{-1/2} = 39.3\%$, and distance 2 holds $86.5\%$, because $\Delta^2$ is a sum of $d$ squared standard normals and has the chi-squared distribution with $d$ degrees of freedom.

Affine maps#

Linear maps keep Gaussians Gaussian. If $X \sim N(\mu, \Sigma)$, then for any matrix $B$ and vector $c$

$$BX + c \sim N(B\mu + c,\ B \Sigma B^\top) .$$

The mean follows from linearity of expectation and the covariance from bilinearity of covariance; that the result is again Gaussian follows, for an invertible square $B$, from the change-of-variables formula of the functions of random variables note in $d$ dimensions, and for any $B$ from the characteristic function; when $B \Sigma B^\top$ is singular, the result is a degenerate Gaussian concentrated on a subspace, without a density. Three consequences are used constantly:

• Sampling. $X = \mu + LZ$ with $LL^\top = \Sigma$ (a Cholesky factor) has distribution $N(\mu, \Sigma)$; the Gaussian distributions note of the generative models course works through it in PyTorch.
• Whitening. Conversely, $L^{-1}(X - \mu) \sim N(0, I)$: any Gaussian is a standard one in suitable coordinates.
• Linear combinations. $a^\top X \sim N(a^\top \mu, a^\top \Sigma a)$. For the example above, whose mean is 0, $X_1 + X_2 \sim N(0, 2 + 1 + 2 \cdot 1.2) = N(0, 5.4)$.

Marginals and conditionals#

Split $X$ into two blocks $X_1$, $X_2$, with the mean and covariance split accordingly:

$$\mu = \begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix}, \qquad \Sigma = \begin{pmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{pmatrix} .$$

The marginal of $X_2$ is read off: $X_2 \sim N(\mu_2, \Sigma_{22})$, a special case of an affine map that selects coordinates. The conditional of $X_2$ given $X_1 = x_1$ is Gaussian as well, with

$$E[X_2 \mid X_1 = x_1] = \mu_2 + \Sigma_{21} \Sigma_{11}^{-1} (x_1 - \mu_1), \qquad \operatorname{Cov}(X_2 \mid X_1 = x_1) = \Sigma_{22} - \Sigma_{21} \Sigma_{11}^{-1} \Sigma_{12} .$$

For two scalar coordinates the blocks are numbers, and the formulas read $E[X_2 \mid X_1 = x_1] = \mu_2 + (\Sigma_{21} / \Sigma_{11})(x_1 - \mu_1)$ and $\operatorname{Var}(X_2 \mid X_1 = x_1) = \Sigma_{22} - \Sigma_{21}^2 / \Sigma_{11}$. For the example, $\Sigma_{21}/\Sigma_{11} = 1.2/2 = 0.6$, so given $X_1 = x_1$ the second coordinate is $N(0.6\, x_1,\ 1 - 1.2^2/2) = N(0.6\, x_1, 0.28)$:

Figure

Conditioning a Gaussian vector

Two features of the formulas matter far beyond this example. The conditional mean is a linear function of the observed value, so for jointly Gaussian data linear regression is not an approximation but the exact best predictor, the subject of the conditional expectation note. And the conditional covariance, the Schur complement $\Sigma_{22} - \Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}$, does not depend on the observed value: observing $X_1$ shrinks the uncertainty about $X_2$ by the same amount wherever $X_1$ lands. Gaussian process regression and the Kalman filter are these two formulas applied to large covariance matrices.

Uncorrelated means independent, for Gaussian vectors#

For a jointly Gaussian vector, zero covariance between two blocks makes $\Sigma$ block-diagonal, and the density factorizes: uncorrelated implies independent. This is a property of the joint distribution, not of the marginals. Let $X \sim N(0, 1)$ and let $Y = SX$, where $S = \pm 1$ is an independent fair sign. Then $Y \sim N(0, 1)$ as well, and $\operatorname{Cov}(X, Y) = E[S]\, E[X^2] = 0$, but $\lvert Y \rvert = \lvert X \rvert$: the two are strongly dependent. The pair is not jointly Gaussian: every linear combination of a Gaussian vector is normal, and $X + Y$ equals 0 with probability $1/2$, which no normal variable does.

High dimensions: the shell#

Optional depth: why high-dimensional Gaussian samples avoid the centre. The rest of the note does not depend on it.

The density of $N(0, I)$ is largest at the origin, and yet samples in high dimension are never near it. The squared length $\lVert Z \rVert^2 = Z_1^2 + \dots + Z_d^2$ has mean $d$ and variance $2d$, so the length concentrates around $\sqrt d$ with a spread that stays near $0.7$:

Figure

The length of a standard normal vector in d dimensions

In 100 dimensions the length is $10 \pm 0.7$, and a length below 5 has probability about $10^{-15}$. The density is high near the centre, but the centre has almost no volume; the volume, and with it the probability, is in a thin shell. Two practical consequences: interpolating linearly between two latent samples of a generative model passes through a region the model never saw during training, which is why spherical interpolation is often used instead; and in high dimensions, distances between independent Gaussian points are all about $\sqrt{2d}$, so nearest neighbours become less informative.

What changed#

	A general joint distribution	The multivariate Gaussian
Parameters	a function of $d$ arguments	$\mu$ and $\Sigma$: $d + d(d+1)/2$ numbers
Contours	anything	ellipsoids, axes from the eigenvectors of $\Sigma$
Linear maps, marginals	new distributions to compute	$N(B\mu + c, B\Sigma B^\top)$; read off blocks
Conditionals	renormalized slices	Gaussian, linear mean, Schur complement covariance
Uncorrelated	not enough for independence	equivalent to independence

Limitations#

• Thin tails and one mode. Real data in high dimension are rarely Gaussian: images, text embeddings and returns have heavy tails, clusters and curved structure; mixtures, flows and diffusion models exist for that reason.
• The covariance matrix is expensive. It has $d(d+1)/2$ free entries, 500,500 for $d = 1{,}000$; estimating it needs many more samples than dimensions, and inverting it costs $O(d^3)$. Diagonal, low-rank or sparse-precision structures are used in practice.
• Marginal normality proves nothing about the joint. Each coordinate can be normal while the vector is not, as the sign example shows; formulas for conditionals and independence apply only when the vector is jointly Gaussian.

What you should remember#

• $N(\mu, \Sigma)$ has density $\propto \exp(-\frac12 (x - \mu)^\top \Sigma^{-1} (x - \mu))$; contours are ellipsoids of constant Mahalanobis distance.
• The eigenvectors of $\Sigma$ give the axes, the square roots of its eigenvalues the half-lengths; in 2D one Mahalanobis unit holds 39%, not 68%.
• $BX + c \sim N(B\mu + c, B\Sigma B^\top)$; $X = \mu + LZ$ with $LL^\top = \Sigma$ samples any Gaussian.
• Marginals are blocks of $\mu$ and $\Sigma$; conditionals are Gaussian with a mean linear in the observation and covariance $\Sigma_{22} - \Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12}$.
• For jointly Gaussian vectors, uncorrelated means independent.
• A standard normal vector in $d$ dimensions has length about $\sqrt d \pm 0.7$.

Code#

The example of the figures checked by sampling, with NumPy:

import numpy as np
from scipy import stats

S = np.array([[2.0, 1.2], [1.2, 1.0]])
w, V = np.linalg.eigh(S)
print(w.round(3), np.abs(V[:, 1]).round(3))                         # [0.2 2.8] [0.832 0.555]: the long axis
print(round(1.2 / np.sqrt(2), 2), [round(float(stats.chi2(2).cdf(k * k)), 3) for k in (1, 2)])   # 0.85 [0.393, 0.865]

rng = np.random.default_rng(0)
L = np.linalg.cholesky(S)                                            # sampling: X = L Z
X = rng.normal(size=(1_000_000, 2)) @ L.T
print(np.cov(X.T).round(2).tolist())                                 # [[2.0, 1.2], [1.2, 1.0]]
d2 = np.einsum('ij,jk,ik->i', X, np.linalg.inv(S), X)                # squared Mahalanobis distances
print(round(float((d2 <= 1).mean()), 3), round(float((d2 <= 4).mean()), 3))   # 0.393 0.865
print(round(float(X.sum(axis=1).var()), 2))                          # 5.39, against 2 + 1 + 2 * 1.2 = 5.4

# conditioning on X1 = 2: the samples with X1 close to 2
near = X[np.abs(X[:, 0] - 2) < 0.02, 1]
print(round(1.2 / 2 * 2, 2), round(1 - 1.2**2 / 2, 2), round(float(near.mean()), 2), round(float(near.var()), 2))   # 1.2 0.28 1.2 0.28

# normal marginals, zero correlation, and still dependent: Y = S X with a random sign S
x = rng.normal(size=1_000_000)
y = rng.choice([-1, 1], size=x.size) * x
print(round(float(np.corrcoef(x, y)[0, 1]), 3), round(float((x + y == 0).mean()), 3))   # -0.001 0.501

# the shell: the length of a standard normal vector in 100 dimensions
r = np.linalg.norm(rng.normal(size=(100_000, 100)), axis=1)
print(round(float(r.mean()), 2), round(float(r.std()), 2), f'{stats.chi(100).cdf(5):.1e}')   # 9.97 0.71 1.1e-15

Next#

For Gaussians the best prediction of $X_2$ from $X_1$ turned out to be linear, with a known variance left over. The question makes sense for any pair of random variables: what function of $X$ predicts $Y$ best, and how much error is left? The next note answers it with the conditional expectation $E[Y \mid X]$, the regression function that a model trained with squared error approximates; other losses target other features of the conditional distribution, such as its median or a quantile.

References#

• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Section 2.3 (the Gaussian distribution, conditional and marginal Gaussians).
• M. P. Deisenroth, A. A. Faisal, C. S. Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. Section 6.5 (Gaussian distribution).
• J. K. Blitzstein, J. Hwang. Introduction to Probability, 2nd ed. CRC Press, 2019. Section 7.5 (multivariate normal).
• T. White. Sampling Generative Networks. arXiv:1609.04468, 2016. Spherical interpolation in latent spaces.