Multivariate Gaussian Distribution

The multivariate Gaussian is the load-bearing distribution of quantitative physics: whenever several measured quantities are uncertain and their uncertainties are entangled, this is the law that describes them. Where the scalar normal needs a mean $\mu$ and a variance $\sigma^2$ , its $n$ -dimensional successor needs a mean vector $\boldsymbol{\mu}\in\mathbb{R}^n$ and a covariance matrix $\boldsymbol{\Sigma}\in\mathbb{R}^{n\times n}$ . The diagonal of $\boldsymbol{\Sigma}$ holds the individual variances; the off-diagonal entries record how pairs of variables co-vary, and it is exactly those off-diagonal numbers — the correlations — that the univariate picture cannot represent and that real experiments cannot avoid. Everything the distribution does is organized by a single scalar built from those two objects: the quadratic form $(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})$ . Its level sets are ellipsoids centred at $\boldsymbol{\mu}$ , tilted and stretched along the eigenvectors of $\boldsymbol{\Sigma}$ , with semi-axes scaling as the square roots of the eigenvalues. These are the "error ellipses" a physicist quotes for two correlated fit parameters, and the value of that quadratic form is the Mahalanobis distance — the honest, metric-aware count of how many standard deviations from the mean a point lies, once correlation is taken into account. What makes the family indispensable, though, is its algebraic closure. Marginalize out variables, condition on observed ones, or push the vector through any linear map, and the result is again Gaussian, with means and covariances given by tidy matrix formulas. Conditioning in particular reproduces linear regression plus a reduced Schur-complement covariance — it is the Bayesian update for jointly normal data, and iterating it is the Kalman filter. Add the multivariate central limit theorem, which makes this distribution the universal attractor of sums of independent vectors, and its maximum-entropy pedigree, and you see why Gaussian noise is the default model across physics, statistics, and machine learning.

Key takeaways

The multivariate Gaussian generalizes the scalar normal to $n$ correlated variables and is fixed entirely by a mean vector $\boldsymbol{\mu}$ and a symmetric positive-definite covariance matrix $\boldsymbol{\Sigma}$ .
Its density carries a quadratic form $-\tfrac12(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})$ in the exponent, whose level sets are concentric ellipsoids — the error ellipses of correlated measurements — with the value of the form being the squared Mahalanobis distance, $d_M^2\sim\chi^2_n$ .
Every marginal and every conditional of a Gaussian is Gaussian, and every affine image $A\mathbf{x}+\mathbf{b}$ is Gaussian with covariance $A\boldsymbol{\Sigma}A^{\top}$ — the whole of linear error propagation in one line.
Diagonalizing $\boldsymbol{\Sigma}$ rotates to principal axes on which the components become independent, with the eigenvalues as the variances along those axes; whitening turns any Gaussian into a standard isotropic one.
The characteristic function $\exp(i\boldsymbol{\mu}^{\top}\mathbf{t}-\tfrac12\mathbf{t}^{\top}\boldsymbol{\Sigma}\mathbf{t})$ has a quadratic exponent, so all cumulants above second order vanish and every higher moment is a sum over pairings — Isserlis' (Wick's) theorem.
Key idea: A multivariate Gaussian is completely determined by its first two moments $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ ; everything — marginals, conditionals, linear images, all higher moments — follows by linear algebra on these two objects.

Where this lives today

The multivariate Gaussian is quietly running much of quantitative science. In cosmology, the Planck satellite's constraints on the matter density $\Omega_m$ and the fluctuation amplitude $\sigma_8$ are published as tilted 2D error ellipses — literal level sets of a posterior Gaussian — and combining probes means intersecting ellipses. In navigation and control, the Kalman filter (equation (5) iterated) fuses noisy GPS, inertial, and radar data in real time aboard aircraft, spacecraft, and self-driving cars; it flew on Apollo and runs in every smartphone's location stack. In machine learning, Gaussian-process regression treats an unknown function as one enormous jointly Gaussian vector and conditions on observed points to predict the rest with calibrated error bars. In gravitational-wave astronomy, LIGO reports correlated posteriors on masses and spins as multivariate Gaussians whose Fisher-information matrix is the inverse covariance. And in lattice field theory, Wick's theorem (equation (7)) reduces every free-field correlator to a sum over pairings — the combinatorial engine behind Feynman diagrams.

The physical picture

Picture a fuzzy ellipsoidal cloud of points hovering around a centre $\boldsymbol{\mu}$ . Its shape is dictated entirely by $\boldsymbol{\Sigma}$ : a round ball when all variances are equal and nothing is correlated, a cigar tilted along a diagonal when two variables tend to rise and fall together. The eigenvectors of $\boldsymbol{\Sigma}$ are the principal axes of the cloud and the eigenvalues are the variances measured along them, so the geometry of the cloud and the spectral decomposition of the covariance are the same fact told two ways. The first figure makes this concrete with a real $2\times2$ matrix, its eigenvalues, its correlation, and a test point placed at a definite Mahalanobis distance.

Figure

The covariance ellipse of Σ = (4, 2 ; 2, 3): its axes lie along the eigenvectors of Σ with semi-lengths √λ₁ and √λ₂, the eigenvalues are λ₁ ≈ 5.56 and λ₂ ≈ 1.44, and the tilt of ≈ 38° is set by the correlation ρ = 2/√12 ≈ 0.577; the point (3, 5) sits at Mahalanobis distance d_M² = 3.0, just outside the 68% ellipse where χ²₂ = 2.28 and inside the 95% ellipse where χ²₂ = 5.99

The inverse covariance $\boldsymbol{\Sigma}^{-1}$ — the precision or information matrix — supplies the natural metric on this cloud. Distances measured with it, the Mahalanobis distances, count standard deviations while discounting directions of large spread and accounting for correlation, so a point one error-ellipse out is exactly one Mahalanobis unit away in every direction. Because $\boldsymbol{\Sigma}$ is positive-definite it admits a whitening map $\mathbf{z}=\boldsymbol{\Sigma}^{-1/2}(\mathbf{x}-\boldsymbol{\mu})$ that turns any Gaussian into the standard isotropic one, so almost every computation reduces to that clean case and an affine change back. The second figure shows the most useful consequence of the geometry: conditioning on an observed variable is just slicing the ellipse, and the slice is again a narrower Gaussian.

Figure

Conditioning is Gaussian regression: fixing x₂ = 5 slices the joint ellipse along the dashed vertical, and the conditional law of x₁ is again Gaussian, centred on the regression line x₁ = 1 + (2/3)(x₂ − 2) = 3 with its variance sharpened from the marginal σ₁² = 4 to the Schur complement 4 − 2²/3 = 8/3 ≈ 2.67

The core equations

The density is the whole distribution in one line: a normalizing prefactor times the exponential of minus one-half the quadratic form, whose value is the squared Mahalanobis distance.

f(\mathbf{x})=\frac{1}{(2\pi)^{n/2}\,(\det\boldsymbol{\Sigma})^{1/2}}\exp\!\Big(-\tfrac12(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\Big) \qquad (1)

where: $\mathbf{x}$ — point in $\mathbb{R}^n$ , a vector of $n$ jointly random variables; $\boldsymbol{\mu}$ — mean vector, $\mu_i=\mathbb{E}[x_i]$ ; $\boldsymbol{\Sigma}$ — covariance matrix, symmetric positive-definite, $\Sigma_{ij}=\mathrm{Cov}(x_i,x_j)$ ; $\det\boldsymbol{\Sigma}$ — determinant, the squared volume factor of the error ellipsoid; $n$ — dimension. Reading it: the density is constant on ellipsoids of constant quadratic form and falls off as a Gaussian in the Mahalanobis distance from $\boldsymbol{\mu}$ . That quadratic form is the correct notion of distance in a correlated cloud, and under the model it is a random variable with a known law — which turns "how far out is this point?" into a calibrated probability.

d_M^2(\mathbf{x})=(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\ \sim\ \chi^2_n \qquad (2)

where: $\boldsymbol{\Sigma}^{-1}$ — precision (information) matrix, the metric; $d_M$ — Mahalanobis distance; $\chi^2_n$ — chi-squared law with $n$ degrees of freedom. Reading it: a level- $\alpha$ confidence region is the ellipsoid $d_M^2\le\chi^2_{n,\alpha}$ — the ellipse geometry and the $\chi^2$ table are one object. The characteristic function is the Fourier transform of the density; for a Gaussian it is again a Gaussian in $\mathbf{t}$ with a quadratic exponent — the single fact behind closure under linear maps, the central limit theorem, and Wick's theorem.

\varphi(\mathbf{t})=\mathbb{E}\big[e^{i\mathbf{t}^{\top}\mathbf{x}}\big]=\exp\!\Big(i\boldsymbol{\mu}^{\top}\mathbf{t}-\tfrac12\mathbf{t}^{\top}\boldsymbol{\Sigma}\mathbf{t}\Big) \qquad (3)

where: $\mathbf{t}$ — Fourier-conjugate vector in $\mathbb{R}^n$ ; $i$ — imaginary unit. Reading it: $\ln\varphi$ is quadratic, so the first cumulant is $\boldsymbol{\mu}$ , the second is $\boldsymbol{\Sigma}$ , and all cumulants of order $\ge3$ vanish — the algebraic definition of "Gaussian". Marginalizing is the easiest operation of all: to get the law of a sub-block, delete the other rows and columns of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ and read off what remains — no integration needed once you trust closure.

\begin{pmatrix}\mathbf{x}_a\\ \mathbf{x}_b\end{pmatrix}\sim\mathcal{N}\!\left(\begin{pmatrix}\boldsymbol{\mu}_a\\ \boldsymbol{\mu}_b\end{pmatrix},\begin{pmatrix}\boldsymbol{\Sigma}_{aa}&\boldsymbol{\Sigma}_{ab}\\ \boldsymbol{\Sigma}_{ba}&\boldsymbol{\Sigma}_{bb}\end{pmatrix}\right)\ \Rightarrow\ \mathbf{x}_a\sim\mathcal{N}(\boldsymbol{\mu}_a,\boldsymbol{\Sigma}_{aa}) \qquad (4)

where: $\mathbf{x}_a,\mathbf{x}_b$ — a partition of the components into two blocks; $\boldsymbol{\Sigma}_{aa}$ — the diagonal block for the retained variables. Reading it: marginalization ignores what you did not measure and keeps the corresponding sub-matrix untouched. Conditioning — updating $\mathbf{x}_a$ once $\mathbf{x}_b$ is observed — is what makes the Gaussian the workhorse of inference: the mean shifts along a regression line and the covariance shrinks to a Schur complement.

\mathbf{x}_a\mid\mathbf{x}_b\sim\mathcal{N}\!\Big(\boldsymbol{\mu}_a+\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}(\mathbf{x}_b-\boldsymbol{\mu}_b),\ \boldsymbol{\Sigma}_{aa}-\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}\boldsymbol{\Sigma}_{ba}\Big) \qquad (5)

where: $\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}$ — regression-coefficient matrix, the Kalman "gain"; $\boldsymbol{\Sigma}_{aa}-\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}\boldsymbol{\Sigma}_{ba}$ — Schur complement, the reduced covariance. Reading it: the posterior covariance never grows — observing data can only sharpen, never blur, a Gaussian estimate. Affine maps stay Gaussian, and the covariance transforms by sandwiching $\boldsymbol{\Sigma}$ between $A$ and its transpose. This one identity is the entire theory of linear error propagation, correlations included.

\mathbf{y}=A\mathbf{x}+\mathbf{b}\ \Rightarrow\ \mathbf{y}\sim\mathcal{N}\!\big(A\boldsymbol{\mu}+\mathbf{b},\,A\boldsymbol{\Sigma}A^{\top}\big) \qquad (6)

where: $A$ — an $m\times n$ matrix (a Jacobian, to first order, for nonlinear maps); $\mathbf{b}$ — constant shift; $A\boldsymbol{\Sigma}A^{\top}$ — propagated covariance. Reading it: "add errors in quadrature" is the special case of a row-vector $A$ with diagonal $\boldsymbol{\Sigma}$ ; the off-diagonal terms are the correction quadrature misses. Because the log-characteristic function is quadratic, every higher moment is a sum over pairings of the indices — Isserlis' theorem, known in physics as Wick's theorem, with the zero-mean fourth moment as the canonical case.

\mathbb{E}[x_i x_j x_k x_l]=\Sigma_{ij}\Sigma_{kl}+\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk} \qquad (7)

where: $\Sigma_{ij}$ — covariance entries; the three terms are the three pairings of four indices. Reading it: all odd moments vanish, and every even moment is a sum of products of covariances over pairings — the combinatorial backbone of Feynman diagrams. Finally, the differential entropy $H=\tfrac12\ln((2\pi e)^{n}\det\boldsymbol{\Sigma})$ depends on $\boldsymbol{\Sigma}$ only through its determinant — the log-volume of the error ellipsoid — so correlations shrink $\det\boldsymbol{\Sigma}$ and hence $H$ , and the Gaussian is the maximum-entropy law at fixed mean and covariance.

Derivation ANormalizing the multivariate Gaussian by whitening

The goal is to prove that the density (1) integrates to one over $\mathbb{R}^n$ , and to expose why the prefactor takes exactly its form — the $(2\pi)^{n/2}$ from $n$ one-dimensional Gaussian integrals, the $(\det\boldsymbol{\Sigma})^{1/2}$ from the eigenvalue product. The calculation is a change of variables that rotates to the principal axes, where the quadratic form falls into independent squares and the integral factorizes; that same rotation is the whitening map underlying almost every Gaussian computation.

Step 1. Centre the integral with $\mathbf{u}=\mathbf{x}-\boldsymbol{\mu}$ . Translation has unit Jacobian, so $\mathrm{d}^n x=\mathrm{d}^n u$ and the mean drops out, leaving $\int (2\pi)^{-n/2}(\det\boldsymbol{\Sigma})^{-1/2}\exp(-\tfrac12\mathbf{u}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{u})\,\mathrm{d}^n u$ . Why: the normalization cannot depend on where the cloud sits, only on its shape, so removing $\boldsymbol{\mu}$ costs nothing.

Step 2. Diagonalize the (symmetric positive-definite) covariance as $\boldsymbol{\Sigma}=Q\Lambda Q^{\top}$ with $Q$ orthogonal and $\Lambda=\mathrm{diag}(\lambda_1,\dots,\lambda_n)$ , all $\lambda_i>0$ . Then $\boldsymbol{\Sigma}^{-1}=Q\Lambda^{-1}Q^{\top}$ . Why: the spectral theorem guarantees a rotation to axes on which $\boldsymbol{\Sigma}$ is diagonal — the principal axes of the error ellipsoid — and positive-definiteness guarantees the eigenvalues are genuine positive variances.

Step 3. Rotate with $\mathbf{w}=Q^{\top}\mathbf{u}$ . An orthogonal change of variables has $|\det Q|=1$ , so the volume element is preserved, and the quadratic form becomes $\mathbf{u}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{u}=\sum_i w_i^2/\lambda_i$ . Why: rotating onto the eigenbasis is exactly what decouples the correlated variables into independent ones, which is the entire content of "whitening".

Step 4. Factorize. The determinant is $\det\boldsymbol{\Sigma}=\prod_i\lambda_i$ , and the integral splits into a product of $n$ one-dimensional Gaussians, each $\int_{-\infty}^{\infty} e^{-w_i^2/(2\lambda_i)}\,\mathrm{d}w_i=\sqrt{2\pi\lambda_i}$ . Why: once the exponent is a sum, the exponential is a product, and a product of integrals is trivial to evaluate term by term.

Step 5. Collect the pieces: $\dfrac{1}{(2\pi)^{n/2}\prod_i\sqrt{\lambda_i}}\prod_i\sqrt{2\pi\lambda_i}=\dfrac{\prod_i\sqrt{2\pi\lambda_i}}{(2\pi)^{n/2}\prod_i\sqrt{\lambda_i}}=1$ . Why: the prefactor was reverse-engineered to be the reciprocal of the eigenvalue integrals, so the cancellation is exact and the density is normalized. *As a consistency check: each one-dimensional integral gives $\int e^{-w^2/(2\lambda)}\mathrm{d}w=\sqrt{2\pi\lambda}$ , and the eigenvalue product confirms $\det\boldsymbol{\Sigma}=\prod_i\lambda_i$ ; for the worked matrix, $\lambda_1\lambda_2=5.56\times1.44=8=\det\boldsymbol{\Sigma}$ .*

Go deeper · Derivation B · 1 minConditioning as Gaussian regression and the Schur complement

Derivation BConditioning as Gaussian regression and the Schur complement

The goal is to derive equation (5) — that conditioning $\mathbf{x}_a$ on an observed $\mathbf{x}_b$ yields a Gaussian whose mean rides a regression line and whose covariance is the Schur complement of $\boldsymbol{\Sigma}_{bb}$ . The strategy: write the joint exponent with the block inverse of $\boldsymbol{\Sigma}$ , treat $\mathbf{x}_b$ as fixed, and complete the square in $\mathbf{x}_a$ — the matrix on the quadratic term is the inverse conditional covariance, the linear term fixes the mean. This is exactly the calculation that becomes the Kalman-filter measurement update and Gaussian-process regression.

Step 1. Write the precision in blocks, $\boldsymbol{\Sigma}^{-1}=\begin{pmatrix}\Lambda_{aa}&\Lambda_{ab}\\ \Lambda_{ba}&\Lambda_{bb}\end{pmatrix}$ , where the block-inverse formula gives $\Lambda_{aa}=(\boldsymbol{\Sigma}_{aa}-\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}\boldsymbol{\Sigma}_{ba})^{-1}$ and $\Lambda_{ab}=-\Lambda_{aa}\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}$ . Why: the conditional law is governed by the precision, not the covariance, because the exponent of the joint density is written with $\boldsymbol{\Sigma}^{-1}$ , and the block $\Lambda_{aa}$ already contains the Schur complement.

Step 2. Freeze $\mathbf{x}_b$ and look at the joint exponent $-\tfrac12(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})$ purely as a function of $\mathbf{x}_a$ . It is a quadratic in $\mathbf{x}_a$ , so the conditional density $\propto\exp(\text{quadratic})$ is automatically Gaussian — closure is guaranteed before any algebra. Why: a Gaussian conditioned on part of itself cannot become anything but Gaussian, because restricting a quadratic form to an affine slice leaves a quadratic form.

Step 3. Complete the square. Expanding and grouping the $\mathbf{x}_a$ -dependence gives $-\tfrac12(\mathbf{x}_a-\boldsymbol{\mu}_a)^{\top}\Lambda_{aa}(\mathbf{x}_a-\boldsymbol{\mu}_a)-(\mathbf{x}_a-\boldsymbol{\mu}_a)^{\top}\Lambda_{ab}(\mathbf{x}_b-\boldsymbol{\mu}_b)+\text{const}$ , from which the coefficient of the quadratic term identifies the conditional covariance as $\Lambda_{aa}^{-1}=\boldsymbol{\Sigma}_{aa}-\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}\boldsymbol{\Sigma}_{ba}$ . Why: for any Gaussian the matrix multiplying the quadratic term is precisely the inverse covariance, so reading it off is the whole game.

Step 4. Read off the mean by locating the stationary point of the exponent, $\boldsymbol{\mu}_{a\mid b}=\boldsymbol{\mu}_a-\Lambda_{aa}^{-1}\Lambda_{ab}(\mathbf{x}_b-\boldsymbol{\mu}_b)=\boldsymbol{\mu}_a+\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}(\mathbf{x}_b-\boldsymbol{\mu}_b)$ , using $\Lambda_{aa}^{-1}\Lambda_{ab}=-\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}$ . Why: substituting the block-inverse relation collapses the abstract gain $\Lambda_{aa}^{-1}\Lambda_{ab}$ into the familiar regression matrix $\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}$ , exposing conditioning as least-squares regression. *As a numerical check: for the worked scalar blocks $\boldsymbol{\Sigma}_{aa}=4,\ \boldsymbol{\Sigma}_{ab}=2,\ \boldsymbol{\Sigma}_{bb}=3$ , the Schur complement is $4-2^2/3=8/3\approx2.67$ and the regression gain is $2/3$ , exactly the figures shown in the second figure.*

Go deeper · Derivation C · 1 minThe characteristic function, closure under linear maps, and Wick's theorem

Derivation CThe characteristic function, closure under linear maps, and Wick's theorem

The goal is to evaluate $\varphi(\mathbf{t})=\mathbb{E}[e^{i\mathbf{t}^{\top}\mathbf{x}}]$ , obtaining the quadratic-exponent form (3), and to harvest its three consequences in one place: closure under any affine map, moments by differentiation, and Isserlis' pairing theorem. The computation is a single Gaussian integral done by completing the square with a complex shift, and everything else follows because the logarithm of the answer is quadratic — the structural signature that all cumulants above second order vanish.

Step 1. Combine the exponents of $e^{i\mathbf{t}^{\top}\mathbf{x}}$ and the density, shift to $\mathbf{u}=\mathbf{x}-\boldsymbol{\mu}$ , and complete the square: $i\mathbf{t}^{\top}\mathbf{u}-\tfrac12\mathbf{u}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{u}=-\tfrac12(\mathbf{u}-i\boldsymbol{\Sigma}\mathbf{t})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{u}-i\boldsymbol{\Sigma}\mathbf{t})-\tfrac12\mathbf{t}^{\top}\boldsymbol{\Sigma}\mathbf{t}$ . Why: the linear term $i\mathbf{t}^{\top}\mathbf{u}$ is absorbed by a shift of the integration variable into the complex plane, exactly as in the scalar Gaussian, leaving a pure quadratic that integrates to its normalization.

Step 2. The shifted Gaussian integrates to one against the density's normalizing factor, so $\varphi(\mathbf{t})=\exp(i\boldsymbol{\mu}^{\top}\mathbf{t}-\tfrac12\mathbf{t}^{\top}\boldsymbol{\Sigma}\mathbf{t})$ , and setting $\mathbf{t}=\mathbf{0}$ gives $\varphi(\mathbf{0})=1$ as a sanity check. Why: contour-shifting the Gaussian is legitimate because the integrand is entire and decays, so the complex shift changes nothing but the constant pulled outside.

Step 3. For an affine map, $\varphi_{A\mathbf{x}+\mathbf{b}}(\mathbf{t})=e^{i\mathbf{t}^{\top}\mathbf{b}}\varphi_{\mathbf{x}}(A^{\top}\mathbf{t})=\exp\!\big(i\mathbf{t}^{\top}(A\boldsymbol{\mu}+\mathbf{b})-\tfrac12\mathbf{t}^{\top}A\boldsymbol{\Sigma}A^{\top}\mathbf{t}\big)$ , which is again of Gaussian form with mean $A\boldsymbol{\mu}+\mathbf{b}$ and covariance $A\boldsymbol{\Sigma}A^{\top}$ . Why: because $\varphi$ of an affine image is just $\varphi$ evaluated at a transformed argument times a phase, the quadratic-exponent form is preserved, proving equation (6) without touching the density.

Step 4. Differentiate at the origin: $\mathbb{E}[x_j]=(-i)\partial_{t_j}\varphi|_0=\mu_j$ and $\mathbb{E}[x_jx_k]-\mu_j\mu_k=\Sigma_{jk}$ , while for zero mean the fact that $\ln\varphi$ is exactly quadratic forces all connected higher cumulants to zero, so higher moments are sums over pairings — Wick's theorem, e.g. $\mathbb{E}[x_ix_jx_kx_l]=\Sigma_{ij}\Sigma_{kl}+\Sigma_{ik}\Sigma_{jl}+\Sigma_{il}\Sigma_{jk}$ . Why: a moment is a derivative of $\varphi=e^{(\text{quadratic})}$ , and differentiating an exponential-of-quadratic generates precisely the sum over pairings of the covariance factors, which is the combinatorics of free-field contractions. *As consistency checks: $\varphi(\mathbf{0})=1$ ; the log-characteristic function is quadratic, the defining hallmark of Gaussianity; and the fourth-moment pairing gives $\mathbb{E}[x_a^2x_b^2]=\Sigma_{aa}\Sigma_{bb}+2\Sigma_{ab}^2=12+2\cdot4=20$ for the worked matrix.*

Limiting cases

The general $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ specializes in several instructive ways, each recovering a distribution a physicist already knows.

Case	Result	Meaning
Scalar $n=1$	$f(x)=\dfrac{1}{\sqrt{2\pi\sigma^2}}\exp\!\Big(-\dfrac{(x-\mu)^2}{2\sigma^2}\Big)$	With $\boldsymbol{\Sigma}=\sigma^2$ a $1\times1$ matrix the formula collapses to the ordinary normal distribution.
Standard isotropic	$\boldsymbol{\mu}=\mathbf{0},\ \boldsymbol{\Sigma}=\mathbf{I}_n\Rightarrow f(\mathbf{z})=(2\pi)^{-n/2}e^{-\tfrac12\mathbf{z}^{\top}\mathbf{z}}$	Components are i.i.d. $\mathcal{N}(0,1)$ ; level sets are spheres. Any Gaussian is an affine image of this one.
Diagonal $\boldsymbol{\Sigma}$	$\boldsymbol{\Sigma}=\mathrm{diag}(\sigma_1^2,\dots,\sigma_n^2)\Rightarrow f=\prod_i\dfrac{1}{\sqrt{2\pi\sigma_i^2}}e^{-(x_i-\mu_i)^2/2\sigma_i^2}$	Uncorrelated jointly Gaussian components are independent — the density factorizes.
Degenerate (singular)	$\det\boldsymbol{\Sigma}=0$	No density on $\mathbb{R}^n$ ; the mass lives on the affine subspace $\boldsymbol{\mu}+\mathrm{range}(\boldsymbol{\Sigma})$ — a constraint $\mathbf{a}^{\top}\mathbf{x}=c$ holds with probability one.

Common trapUncorrelated does not mean independent — except for jointly Gaussian variables

The slogan "zero covariance implies independence" is false in general and true in exactly the case at hand — confusing the two is the most common Gaussian blunder. Two marginally normal variables can be uncorrelated yet strongly dependent: take $X\sim\mathcal{N}(0,1)$ and $Y=SX$ with $S=\pm1$ a random sign, so $\mathrm{Cov}(X,Y)=0$ while $|X|=|Y|$ always; the pair is not jointly Gaussian. For a genuine multivariate Gaussian, a zero off-diagonal entry makes that block of $\boldsymbol{\Sigma}$ diagonal, the density factorizes, and independence really does follow. A second trap: formula (1) requires $\boldsymbol{\Sigma}$ positive-definite, not merely positive-semi-definite. A singular $\boldsymbol{\Sigma}$ has no inverse and no usable determinant, so there is no density on $\mathbb{R}^n$ — the distribution is a degenerate Gaussian on a lower-dimensional subspace, handled with the pseudo-inverse or by dropping the redundant directions.

Worked examples

1 · Error ellipse, eigenstructure, and Mahalanobis distance (fully worked)

A two-parameter fit returns $\boldsymbol{\mu}=(1,2)$ and covariance $\boldsymbol{\Sigma}=\begin{pmatrix}4&2\\2&3\end{pmatrix}$ . Find the correlation coefficient, the principal axes of the error ellipse, and decide whether the point $(3,5)$ lies inside the 68% or 95% confidence region.

Step 1. Correlation. By definition $\rho=\Sigma_{12}/\sqrt{\Sigma_{11}\Sigma_{22}}=2/\sqrt{4\cdot3}=2/\sqrt{12}=\mathbf{0.577}$ . A positive $\rho$ tilts the ellipse up to the right; its magnitude below one confirms $\boldsymbol{\Sigma}$ is non-degenerate.

Step 2. Eigenstructure. Solve $\det(\boldsymbol{\Sigma}-\lambda\mathbf{I})=0$ , i.e. $\lambda^2-7\lambda+8=0$ , giving $\lambda=\tfrac{7\pm\sqrt{17}}{2}$ , so $\lambda_1=\mathbf{5.56}$ and $\lambda_2=\mathbf{1.44}$ . As a check, $\lambda_1+\lambda_2=7=\mathrm{tr}\,\boldsymbol{\Sigma}$ and $\lambda_1\lambda_2=8=\det\boldsymbol{\Sigma}$ . The semi-axes of the one-sigma ellipse are $\sqrt{\lambda_1}=\mathbf{2.36}$ and $\sqrt{\lambda_2}=\mathbf{1.20}$ , and the major axis tilts at $\approx\mathbf{38^\circ}$ (the eigenvector direction $(0.79,0.62)$ ).

Step 3. Mahalanobis distance of $(3,5)$ . The inverse is $\boldsymbol{\Sigma}^{-1}=\tfrac18\begin{pmatrix}3&-2\\-2&4\end{pmatrix}$ , and with $\mathbf{x}-\boldsymbol{\mu}=(2,3)$ we get $d_M^2=\tfrac18\,(2,3)\begin{pmatrix}3&-2\\-2&4\end{pmatrix}\begin{pmatrix}2\\3\end{pmatrix}=\tfrac18(0\cdot2+8\cdot3)=\mathbf{3.0}$ .

Step 4. Verdict. For $n=2$ the thresholds are $\chi^2_{2,0.68}=2.28$ and $\chi^2_{2,0.95}=5.99$ . Since $2.28<3.0<5.99$ , the point lies outside the 68% ellipse but inside the 95% ellipse — at the $\Pr(\chi^2_2\le3.0)=\mathbf{77.7\%}$ contour. Correlation matters: naively rescaling each coordinate by its own $\sigma$ (ignoring $\rho$ ) gives $(2/2)^2+(3/\sqrt3)^2=4\neq3.0$ , over-stating the distance by double-counting the shared part.

2 · Conditioning as a measurement update

With the same $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ , a noiseless measurement pins $x_2=5$ . Update the estimate of $x_1$ .

Step 1. The gain (regression coefficient) is $\Sigma_{12}\Sigma_{22}^{-1}=2/3$ , so the conditional mean is $\mu_{1\mid2}=\mu_1+\tfrac23(x_2-\mu_2)=1+\tfrac23(5-2)=\mathbf{3}$ . The estimate of $x_1$ moves up because $x_1$ and $x_2$ are positively correlated and $x_2$ came in above its mean.

Step 2. The conditional variance is the Schur complement $\Sigma_{11}-\Sigma_{12}^2/\Sigma_{22}=4-4/3=\mathbf{8/3\approx2.67}$ , down from the prior variance $4$ . Observing a correlated quantity has sharpened our knowledge of $x_1$ by a factor $\sigma^2_{\text{prior}}/\sigma^2_{\text{post}}=4/(8/3)=1.5$ .

Step 3. So $x_1\mid x_2{=}5\sim\mathcal{N}(3,\,8/3)$ . Iterate this over a stream of noisy measurements and it is the Kalman filter; make the index a continuous input and it is Gaussian-process regression. The structural point: the posterior variance $8/3$ is independent of the measured value $5$ — only the mean moves — a hallmark of linear-Gaussian inference.

3 · Correlated error propagation and a Wick moment

Two zero-mean detector readings have covariance $\boldsymbol{\Sigma}=\begin{pmatrix}4&2\\2&3\end{pmatrix}$ . Find the variance of the sum $y_+=x_1+x_2$ and the difference $y_-=x_1-x_2$ , and evaluate $\mathbb{E}[x_1^2x_2^2]$ .

Step 1. Sum. With $A=(1,\ 1)$ , equation (6) gives $\mathrm{Var}(y_+)=A\boldsymbol{\Sigma}A^{\top}=\Sigma_{11}+\Sigma_{22}+2\Sigma_{12}=4+3+4=\mathbf{11}$ . The naive "add in quadrature" answer $\Sigma_{11}+\Sigma_{22}=7$ is too small: positive correlation makes the readings fluctuate together, amplifying the sum.

Step 2. Difference. With $A=(1,\ -1)$ , $\mathrm{Var}(y_-)=\Sigma_{11}+\Sigma_{22}-2\Sigma_{12}=4+3-4=\mathbf{3}$ . The same correlation that inflated the sum cancels in the difference — which is exactly why differential (common-mode-rejecting) measurements are so powerful when the noise is correlated.

Step 3. Wick. For zero-mean jointly Gaussian variables, equation (7) with $i=j=1,\ k=l=2$ gives $\mathbb{E}[x_1^2x_2^2]=\Sigma_{11}\Sigma_{22}+2\Sigma_{12}^2=4\cdot3+2\cdot4=\mathbf{20}$ . The " $+2\Sigma_{12}^2$ " is the correlation contribution; for independent readings it would vanish and leave just $\Sigma_{11}\Sigma_{22}=12$ .

Physics of the portrait

The multivariate Gaussian has no single inventor — it accreted across a century, from Gauss's law of errors, through Bravais's correlated ellipses, to Galton's discovery that regression and correlation are geometric features of one joint law. The most human chapter belongs to Francis Galton, who saw the ellipse before he could derive it.

Origins · 1 minFrancis Galton and how this was discovered

Francis Galton

1822–1911 · the bivariate normal, regression, and correlation

Francis Galton — polymath, geographer, and cousin of Charles Darwin — came to the bivariate normal through data, not mathematics. Tabulating the heights of parents and their grown children in the 1870s and 1880s, he connected cells of equal joint frequency and was startled to find they traced concentric, tilted ellipses: the level sets of a two-dimensional Gaussian. The locus of mean child-height for each parent-height fell on a straight line of slope less than one — tall parents had tall but, on average, less-extreme children — and he named this pull toward the mean regression. Unable to prove the ellipses exact, he handed the picture to the Cambridge mathematician J. Hamilton Dickson, who in 1886 recovered them from the bivariate normal density, confirming that the "regression line" was precisely the conditional-mean formula $\mu_{a}+\rho(\sigma_a/\sigma_b)(x_b-\mu_b)$ . In 1888 Galton distilled the tilt into a single index of co-relation — the correlation coefficient $\rho$ — and modern multivariate statistics was born. Karl Pearson made $\rho$ rigorous; Auguste Bravais had already written the correlated $n$ -dimensional density and its error ellipsoid in 1846; Leon Isserlis proved the moment-pairing theorem in 1918 that physicists rediscovered as Wick's theorem; and Prasanta Chandra Mahalanobis defined his correlation-aware distance in 1936. Galton's gift was to see the geometry — ellipse, regression line, correlation — in a cloud of ordinary measurements, and to insist those pictures were laws.

History

1809

— Gauss frames the normal law of errors and least squares in Theoria Motus

1846

— Bravais derives the correlated multivariate normal and its elliptical error contours

1885

— Galton discovers regression and the ellipses of equal frequency in his height data

1918

— Isserlis proves the moment-pairing theorem later rediscovered as Wick's theorem

1936

— Mahalanobis introduces his correlation-aware generalized distance

Going further

Texts. Stone & Goldbart, Mathematics for Physics, probability and statistics chapter · Arfken & Weber, Mathematical Methods for Physicists, probability chapter · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, statistics chapters · C. M. Bishop, Pattern Recognition and Machine Learning, §2.3 (Gaussians) · A. Papoulis & S. U. Pillai, Probability, Random Variables and Stochastic Processes · C. E. Rasmussen & C. K. I. Williams, Gaussian Processes for Machine Learning, Ch. 2. The key numerical anchors used throughout this article are: $\det\boldsymbol{\Sigma}=\prod_i\lambda_i$ (here $5.56\times1.44=8$ ); characteristic function $\exp(i\boldsymbol{\mu}^{\top}\mathbf{t}-\tfrac12\mathbf{t}^{\top}\boldsymbol{\Sigma}\mathbf{t})$ ; chi-squared thresholds $\chi^2_{2,0.68}=2.28$ , $\chi^2_{2,0.95}=5.99$ , and $\Pr(\chi^2_2\le3.0)=0.777$ ; worked figures $\rho=0.577$ , $\lambda=\{5.56,1.44\}$ , $d_M^2(3,5)=3.0$ , conditional $\mathcal{N}(3,8/3)$ , $\mathrm{Var}(x_1\pm x_2)=\{11,3\}$ , Wick moment $\mathbb{E}[x_1^2x_2^2]=20$ , and differential entropy $\tfrac12\ln((2\pi e)^n\det\boldsymbol{\Sigma})$ . All results are standard textbook identities; no nonstandard constants are introduced.

Prerequisites. Gaussian (Normal) Distribution · Joint, Marginal and Conditional Distributions · Expectation, Variance and Moments · Moment-Generating and Characteristic Functions · Cumulative Distribution Functions and Change of Variables. Leads to. Central Limit Theorem · Conditional Probability and Bayes' Theorem · Laws of Large Numbers and Convergence of Random Variables.

Ask the tutor: Hand me any $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ and I will draw its error ellipse, condition it on a measurement, or propagate it through your favourite linear map, step by step.

Why is the multivariate Gaussian so much more useful than just stacking independent normals?

Because real quantities are correlated, and the off-diagonal entries of $\boldsymbol{\Sigma}$ are the only place that information can live. Stacking $n$ independent normals forces a diagonal covariance, which cannot represent that two fit parameters trade off or that two detectors share common-mode noise. The multivariate Gaussian keeps the full $\boldsymbol{\Sigma}$ and stays closed under marginalizing, conditioning, and linear maps even when $\boldsymbol{\Sigma}$ is dense — which is what lets you propagate correlated errors, fuse correlated measurements, and read off regression lines by matrix algebra rather than by integration.

What exactly is the Mahalanobis distance, and why not just use ordinary distance?

It is the distance $\sqrt{(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})}$ measured with the precision matrix as its metric. Ordinary Euclidean distance treats all directions alike, so it calls a point "far" even along a direction of huge natural spread. The Mahalanobis distance rescales each direction by its standard deviation and removes the shared, correlated part, so one unit always means "one error-ellipse out". That makes it the right tool for outlier detection and confidence regions: its square follows a $\chi^2_n$ law, converting geometry directly into a probability.

Why does conditioning always shrink the covariance?

Because information can only help. The conditional covariance is the Schur complement $\boldsymbol{\Sigma}_{aa}-\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}\boldsymbol{\Sigma}_{ba}$ , and the subtracted term is positive-semi-definite, so the posterior covariance is never larger than the prior $\boldsymbol{\Sigma}_{aa}$ . Observing a variable correlated with $\mathbf{x}_a$ tightens the estimate; if the two are uncorrelated ( $\boldsymbol{\Sigma}_{ab}=0$ ) the subtracted term vanishes and you learn nothing. Note the shrinkage is by a fixed matrix, independent of the value observed — a peculiarly linear-Gaussian feature.

Can I have a Gaussian with a singular covariance matrix?

Yes, but it has no density on $\mathbb{R}^n$ . A singular (positive-semi-definite but not definite) $\boldsymbol{\Sigma}$ describes a degenerate Gaussian whose mass is confined to the affine subspace $\boldsymbol{\mu}+\mathrm{range}(\boldsymbol{\Sigma})$ — three variables constrained by $x_1+x_2+x_3=0$ live on a plane in $\mathbb{R}^3$ . The density formula breaks because $\boldsymbol{\Sigma}^{-1}$ and $(\det\boldsymbol{\Sigma})^{-1/2}$ do not exist; you handle it with the Moore–Penrose pseudo-inverse, or by dropping the redundant directions. The characteristic function $\exp(i\boldsymbol{\mu}^{\top}\mathbf{t}-\tfrac12\mathbf{t}^{\top}\boldsymbol{\Sigma}\mathbf{t})$ , by contrast, stays perfectly well defined even when $\boldsymbol{\Sigma}$ is singular.

What is Wick's theorem really saying, and why do physicists love it?

It says that for zero-mean jointly Gaussian variables, every higher moment is a sum of products of covariances over all ways of pairing the factors — nothing new appears beyond $\boldsymbol{\Sigma}$ . Odd moments vanish (an odd number of factors cannot be fully paired), and $\mathbb{E}[x_ix_jx_kx_l]$ is the three-term sum of pairings. Physicists love it because a free quantum field has exactly this structure: every correlator reduces to a sum over pairings of two-point functions (propagators), drawn as the lines of Feynman diagrams. Perturbation theory around a Gaussian then becomes bookkeeping of contractions.

How is the Kalman filter just the conditioning rule?

A Kalman filter alternates two Gaussian operations. The predict step pushes the current estimate through the linear dynamics — an affine map, so by equation (6) it stays Gaussian with covariance $A\boldsymbol{\Sigma}A^{\top}$ plus process noise. The update step folds in a noisy measurement by conditioning the joint state-plus-measurement Gaussian on the observed value — exactly equation (5), with $\boldsymbol{\Sigma}_{ab}\boldsymbol{\Sigma}_{bb}^{-1}$ as the Kalman gain and the Schur complement as the shrunken posterior. Iterating predict-then-condition over a data stream is the entire algorithm — no extra machinery, only closure under affine maps and conditioning.