The Gaussian Distribution

The previous two chapters built a language and a toolkit. Chapter 2 gave the rules for reasoning about uncertain quantities: distributions, the sum rule that marginalizes, and the product rule that conditions. Chapter 3 gave the matrices those rules turn into when there are many quantities at once: covariance matrices, their Cholesky factors, and block matrices. This chapter puts the two together in a single distribution, the Gaussian, which carries the rest of the book.

One distribution deserves a chapter because of what a Bayesian optimizer does with its beliefs. It keeps a belief about an unknown function, updates the belief after each evaluation, and asks it where to look next. Each step is an operation on a probability distribution, and for most distributions these operations are integrals with no closed form. For the Gaussian, each has an exact answer in a few lines of linear algebra: a linear map of a Gaussian is Gaussian, ignoring some of its coordinates leaves a Gaussian, and conditioning on some of its coordinates leaves a Gaussian. A Gaussian process (Section 7.3) is a Gaussian over function values, so these three closure properties are also what make Gaussian process regression and Bayesian optimization computable.

The chapter starts with one variable, where the formulas can be read at a glance, moves to many variables, where a covariance matrix gives the distribution its shape, and then derives the three operations in turn. The last of them, conditioning, is the chapter's main result; Chapter 8 applies it without change. A later section separates two ways of combining Gaussians that are easy to confuse, and prepares the Bayesian updating of Chapter 5.

4.1 One dimension #

Suppose a training run with a new learning rate has not happened yet. From experience with similar runs you expect a validation accuracy of about 0.82, rarely off by more than 0.06 in either direction. That belief has a center and a spread, treats deviations up and down alike, and makes large deviations rarer than small ones. The Gaussian distribution, also called the normal distribution, is the standard way to write such a belief with two numbers.

A random variable $X$ is Gaussian with mean $\mu$ and variance $\sigma^2 > 0$, written $X \sim \N(\mu, \sigma^2)$, when its density is

The semicolon separates the point $x$ where the density is evaluated from the parameters $\mu$ and $\sigma^2$. As Section 2.3 explained, a density is not a probability; probabilities are areas under it.

Read the formula from the inside out. The exponent $-(x - \mu)^2 / 2\sigma^2$ is a downward parabola with its peak at $x = \mu$. The exponential turns the parabola into a bell: it equals 1 at the peak and falls quickly as the parabola goes negative. The factor in front scales the bell so that its total area is one. Taking logarithms undoes the exponential:

A Gaussian is the exponential of a quadratic. The converse is the fact that does most of the work in this chapter: any density whose logarithm is a quadratic function of $x$, opening downward, is Gaussian, and its parameters can be read off the coefficients. If

To see it, complete the square: $-\tfrac12\alpha x^2 + \beta x = -\tfrac12\alpha(x - \beta/\alpha)^2 + \beta^2/2\alpha$, and the last term is a constant that the normalization absorbs. The coefficient $\alpha$, the reciprocal of the variance, is called the precision. To show that some distribution is Gaussian, we will repeatedly show that its log density is quadratic and then apply this recipe.

The two parameters mean what their names say. The mean is the expected value, $\E[X] = \mu$, and the variance is the expected squared deviation, $\Var[X] = \E[(X - \mu)^2] = \sigma^2$ (Section 2.6). The standard deviation $\sigma$, the square root of the variance, is in the units of $x$, so it is the number to think with. The accuracy belief above is about $\N(0.82, 0.03^2)$: a standard deviation of 0.03, with deviations of twice that size rare.

4.1.1 The standard normal #

Every Gaussian is a shifted and stretched copy of one reference distribution, the standard normal $Z \sim \N(0, 1)$. If $X \sim \N(\mu, \sigma^2)$, then $Z = (X - \mu)/\sigma$ is standard normal, and conversely $X = \mu + \sigma Z$. The value $z = (x - \mu)/\sigma$, the number of standard deviations by which $x$ lies above the mean, is called its z-score. The standard normal has its own symbols, used throughout the book:

The function $\phi$ is the density and $\Phi$ the cumulative distribution function: $\Phi(z)$ is the probability that $Z \le z$. There is no formula for $\Phi$ in elementary functions. Libraries compute it through the error function, $\Phi(z) = \tfrac12\left(1 + \operatorname{erf}(z/\sqrt{2})\right)$, to full floating-point precision.

Standardizing turns every probability question about $X$ into one about $\Phi$:

A few values are worth remembering. The interval $\mu \pm \sigma$ holds 68.3% of the probability, $\mu \pm 1.96\sigma$ holds 95%, $\mu \pm 2\sigma$ holds 95.4%, and $\mu \pm 3\sigma$ holds 99.7%. The tails fall off fast: a value more than five standard deviations from the mean, in either direction, has probability about $5.7 \times 10^{-7}$.

4.1.2 Why this distribution #

A reader may wonder why this particular bell curve, and not some other, is the default. There are three reasons, of different kinds.

The first comes from a theorem. The central limit theorem says that a sum of many independent random quantities with finite variance, none of which dominates, is approximately Gaussian once standardized, whatever the distribution of the individual terms (Blitzstein and Hwang, 2019). Measurement noise is often the sum of many small disturbances, such as thermal fluctuations, timing jitter, and rounding, and is then close to Gaussian. A classic demonstration: the sum of twelve independent uniform numbers on $[0, 1]$, minus 6, has mean 0 and variance 1, and its histogram is already hard to tell from $\phi$.

The second is a principle. Among all distributions on the real line with a given mean and variance, the Gaussian has the largest entropy, a measure of how spread out a distribution is that Section 6.1 makes precise (Cover and Thomas, 2006, ch. 12). If all we are willing to commit to is a center and a spread, the Gaussian is the choice that assumes nothing more.

The third is convenience, and it is why this book uses Gaussians even where the first two reasons do not apply: every operation in the rest of this chapter has a closed form. Convenience is a reason to choose a model, not evidence that the world is Gaussian. Real quantities can be skewed, bounded, or heavy-tailed, with extreme values far more common than the $5.7 \times 10^{-7}$ above suggests. A comparison between two options, the observation at the heart of Part IV, is not Gaussian at all; Chapter 17 shows how to approximate the resulting posterior by a Gaussian so that the machinery of this chapter still applies.

4.2 Many dimensions #

A Bayesian optimizer never holds a belief about a single number. It holds beliefs about the objective at many inputs at once, and those beliefs are linked: if the accuracy at a learning rate of 0.010 turns out high, the accuracy at 0.011 is probably high too. A list of separate one-dimensional Gaussians cannot express that link. We need a joint distribution over a vector of values that records how each pair of values moves together.

Start with two independent coordinates. If $x_1 \sim \N(0, \sigma_1^2)$ and $x_2 \sim \N(0, \sigma_2^2)$ are independent, their joint density is the product of the two densities (Section 2.7), and multiplying exponentials adds the exponents:

The density is constant wherever the exponent is, on the curves $x_1^2/\sigma_1^2 + x_2^2/\sigma_2^2 = r^2$. These are ellipses with axes along the coordinate directions, and circles when $\sigma_1 = \sigma_2$. To link the two coordinates, we allow the quadratic in the exponent a cross term $x_1 x_2$. The ellipses then tilt, so that a large $x_1$ makes a large $x_2$ more likely, or less likely, depending on the direction of the tilt. All of this is recorded in one matrix.

The matrix is the covariance matrix $\mSigma$ of Section 2.6.3. For a random vector $\vx = (x_1, \dots, x_d)^\T$ with mean vector $\vmu$, its entries are $\Sigma_{ij} = \Cov[x_i, x_j] = \E[(x_i - \mu_i)(x_j - \mu_j)]$, with the variances on the diagonal. It is symmetric and positive semidefinite (Section 3.3), and dividing an entry by the two standard deviations gives the correlation $\rho_{ij} = \Sigma_{ij} / \sqrt{\Sigma_{ii}\Sigma_{jj}}$.

Every piece has a one-dimensional counterpart. With $d = 1$ and $\mSigma = [\sigma^2]$ the formula is Equation (4.1). The exponent is again a quadratic, now a quadratic form in the vector $\vx - \vmu$, with the inverse covariance $\mSigma^{-1}$ in the role of $1/\sigma^2$. The inverse covariance is called the precision matrix, and it is the natural object in several derivations below. The determinant in the normalizer measures the volume over which the distribution spreads (Section 3.6) and plays the role of $\sigma$: a more spread-out distribution has a lower peak, so that the total probability stays one.

The recipe Equation (4.2) carries over unchanged. If a density satisfies

for a positive definite matrix $\bm{\Lambda}$ and a vector $\mathbf{h}$, then $p$ is Gaussian with precision $\bm{\Lambda}$, that is,

Completing the square works as before, and expanding the right side checks it:

The last term does not depend on $\vx$.

4.2.1 Shape #

The quantity in the exponent has a name. The Mahalanobis distance of $\vx$ from $\vmu$ is

the multivariate z-score. In one dimension it is $\lvert x - \mu\rvert / \sigma$, the number of standard deviations. In general it measures distance in units of the distribution's own spread, direction by direction: a point can be far from the mean in ordinary distance and still close in Mahalanobis distance, if it lies along a direction in which the distribution is wide. The density depends on $\vx$ only through $r$, so its contours are the sets where $r$ is constant: ellipses in two dimensions, ellipsoids in more.

Where do the ellipses point? Write the covariance in its eigendecomposition $\mSigma = \mathbf{U}\bm{\Lambda}_{\mathrm{e}}\mathbf{U}^\T$, with orthonormal eigenvectors $\mathbf{u}_i$ in the columns of $\mathbf{U}$ and eigenvalues $\lambda_i$ on the diagonal of $\bm{\Lambda}_{\mathrm{e}}$ (Section 3.4). In the rotated coordinates $\mathbf{y} = \mathbf{U}^\T(\vx - \vmu)$ the quadratic form becomes $\sum_i y_i^2/\lambda_i$, with no cross terms. The axes of every contour therefore point along the eigenvectors, and the ellipse at Mahalanobis distance $r$ has semi-axes of length $r\sqrt{\lambda_i}$. The eigenvalues are the variances along the principal directions, and their product is $\lvert\mSigma\rvert$. A multivariate Gaussian is an axis-aligned bell in some rotated coordinate system, and the eigenvectors say which one.

For two coordinates with standard deviations $\sigma_1, \sigma_2$ and correlation $\rho$, the covariance matrix and its determinant are

The figure below draws this case with mean zero.

Set $\rho$ to zero and the two standard deviations equal. The ellipses become circles. The distribution looks the same in every direction, and the eigenvectors could point anywhere.

Move $\rho$ toward 0.95. The ellipses narrow into a needle along the diagonal. With both standard deviations at 1 the eigenvalues are $1 + \rho$ and $1 - \rho$, so the first approaches 2 while the second and the determinant approach zero: the distribution concentrates near a line, and once $x_1$ is known, $x_2$ is nearly determined. At $\rho = \pm 1$ the covariance would be singular, the density of Equation (4.5) would not exist, and a Cholesky factorization would fail. This is the floating-point trouble that a small diagonal jitter repairs (Section 8.4).

Make $\rho$ negative. The ellipses tilt the other way: a large $x_1$ now goes with a small $x_2$.

Watch the two strips while you move $\rho$. They do not change. Section 4.4 explains why.

How much probability does each ellipse hold? Less than the one-dimensional numbers suggest. In two dimensions the probability inside the ellipse at Mahalanobis distance $r$ is $1 - e^{-r^2/2}$: 39% inside $r = 1$, 86% inside $r = 2$, and 99% inside $r = 3$. Enclosing 95% takes $r \approx 2.45$, not 1.96. The gap widens with the dimension.

One more property is special to Gaussians. When $\mSigma$ is diagonal, the quadratic form has no cross terms and the density factorizes into a product of one-dimensional densities, so the coordinates are independent. For a Gaussian vector, uncorrelated therefore means independent. For other distributions it does not (Section 2.7), and Exercise 4.3 shows two variables that are each Gaussian and uncorrelated, yet dependent, because they are not jointly Gaussian.

4.3 Linear maps and sampling #

Two practical questions lead to the same result. First, if $f(\vx_1)$ and $f(\vx_2)$ are jointly Gaussian, what is the distribution of their difference, or of their average? Part IV needs the difference whenever a person compares two options. Second, a random number generator produces independent standard normal numbers. How do we turn them into a draw from $\N(\vmu, \mSigma)$ with an arbitrary covariance? Every picture of functions drawn from a model needs such draws, and so does Thompson sampling (Section 12.5), a rule that draws one plausible objective from the model and evaluates where that draw is highest.

Both answers follow from one closure property. Let $\vx \sim \N(\vmu, \mSigma)$ in $d$ dimensions, let $\mA$ be an $m \times d$ matrix, and let $\mathbf{c}$ be a vector in $\R^m$. Then

A linear map of a Gaussian is Gaussian. The mean is mapped like a point, and the covariance is sandwiched between the matrix and its transpose.

One caveat. The density of Definition 4.1 needs a positive definite covariance, and $\mA\mSigma\mA^\T$ is positive definite only when no row of $\mA$ is a combination of the others. When it is not, as when $\mA$ has more rows than columns, $\vy$ is still Gaussian in the sense used in step 6 below, but it is confined to a flat of lower dimension and has no density on $\R^m$. The formulas for its mean and covariance hold in both cases.

4.3.1 Sampling with the Cholesky factor #

Sampling runs the map in the useful direction. Let $\vz \sim \N(\mathbf{0}, \mI)$ be a vector of $d$ independent standard normal numbers, which every numerical library provides (the classic construction from uniform random numbers is the transform of Box and Muller (1958)). Choose any matrix $\mL$ with $\mL\mL^\T = \mSigma$ and set

By Equation (4.9), $\vx$ is Gaussian with mean $\vmu$ and covariance $\mL\mI\mL^\T = \mSigma$. The Cholesky factor of Section 3.5, lower triangular with a positive diagonal, is the usual choice of $\mL$: it exists for every positive definite $\mSigma$, it costs $O(d^3)$ operations to compute once, and after that each draw costs one triangular matrix-vector product (Rasmussen and Williams, 2006, app. A.2). Any other square root would do, such as $\mathbf{U}\bm{\Lambda}_{\mathrm{e}}^{1/2}$ from the eigendecomposition. It would map a given $\vz$ to a different point, but the distribution of the points would be the same.

In two dimensions the Cholesky factor of Equation (4.8) can be written down, and multiplying it by $\vz$ shows how each coordinate is built:

Multiplying out confirms $\mL\mL^\T = \mSigma$. The formula also says what correlation is, mechanically. The coordinate $x_2$ is built partly from the same random number $z_1$ that drives $x_1$ and partly from fresh randomness $z_2$; a fraction $\rho^2$ of its variance comes from the shared part. At $\rho = 0$ the two coordinates share nothing, and at $\rho = \pm 1$ they share everything.

Follow the arrows. Every gray point is moved by the same matrix. Because $\mL$ is lower triangular, $x_1 = \sigma_1 z_1$ depends on $z_1$ alone; with $\sigma_1 = 1$ the map leaves the first coordinate unchanged and every arrow is vertical. With a positive $\rho$, the map adds $\rho\sigma_2 z_1$ to the second coordinate, which pushes points on the left down and points on the right up, and it shrinks $z_2$ by the factor $\sqrt{1 - \rho^2}$. The round cloud becomes a tilted ellipse.

Compare the sample correlation with $\rho$. The readout computes the correlation of the blue points. It differs from $\rho$ by a few hundredths, and a fresh set of samples gives a different error: sampling noise, with a standard deviation of roughly $(1 - \rho^2)/\sqrt{n}$ for $n$ points.

Set $\rho$ to zero with unequal standard deviations. $\mL$ is then diagonal and only stretches the cloud along the axes.

For a Gaussian process, $\vx$ holds the function values on a grid of inputs, one coordinate per grid point, and the cubic cost of the factorization is what limits posterior samples to grids of a few thousand points (Section 8.5). Grids run out quickly as the input dimension grows: 20 points per axis is 20 values on a line, 400 on a square, and 8,000 in a three-dimensional cube, whose covariance matrix has 64 million entries. Beyond two or three input dimensions, samples are drawn at a few thousand scattered candidate points instead of a grid.

The difference of two function values is the other question this section began with.

4.4 Marginalizing #

A Gaussian process describes infinitely many function values, but a computer holds finitely many. For that to make sense, the belief about a few values must not depend on which other values we chose to keep track of. In the terms of Section 2.4, we need the marginal distribution of a sub-vector, which the sum rule obtains by integrating out everything else. For most joint densities that integral is the hard part of a calculation. For a Gaussian it costs nothing.

Split the vector into two blocks: $\mathbf{a}$, the coordinates we keep, and $\mathbf{b}$, the coordinates we drop. Partition the mean and the covariance to match:

The diagonal blocks $\mSigma_{aa}$ and $\mSigma_{bb}$ hold the covariances within each block, and $\mSigma_{ab}$ holds the covariances between a coordinate of $\mathbf{a}$ and a coordinate of $\mathbf{b}$. The marginal distribution of $\mathbf{a}$ is

Marginalizing a Gaussian is reading off a sub-block. The proof is one line from the previous section: keeping $\mathbf{a}$ and dropping $\mathbf{b}$ is the linear map with matrix $[\mI \;\; \mathbf{0}]$, an identity block beside a block of zeros, and Equation (4.9) gives mean $[\mI \;\; \mathbf{0}]\vmu = \vmu_a$ and covariance $[\mI \;\; \mathbf{0}]\mSigma[\mI \;\; \mathbf{0}]^\T = \mSigma_{aa}$. The integral the sum rule calls for has been done once and for all; it can also be carried out directly by completing the square (Bishop, 2006, sec. 2.3.2).

Two consequences matter later. First, marginals are consistent: the distribution of $\mathbf{a}$ does not depend on how many other coordinates $\mathbf{b}$ contains, or which. A Gaussian process relies on this to be well defined, since its definition specifies only the joint distributions of finite sets of function values (Section 7.3). Second, the marginal discards the cross-covariances $\mSigma_{ab}$, and with them everything the two blocks say about each other. In Figure 4.1 the strips above and to the right of the plot are the two marginals, and moving the correlation slider rotates and squeezes the joint density without changing either strip. Many joint distributions share the same marginals.

The marginal of $x_2$ answers the question "what do I believe about $x_2$ if I ignore $x_1$?" The next section answers a different question: "what do I believe about $x_2$ once I know $x_1$?"

4.5 Conditioning #

This is the operation the rest of the book runs on. A Bayesian optimizer has evaluated the objective at a few inputs and wants its belief about the objective everywhere else. If the prior belief about the values at all these inputs is a joint Gaussian, the question becomes: once some coordinates of a Gaussian vector have been observed, what is the distribution of the others? Unlike the marginal, the answer must use what was observed.

4.5.1 A slice through the bell #

Look at two dimensions first. By the product rule, the conditional density of $x_2$ given $x_1 = a$ is

The numerator is the joint density along the vertical line $x_1 = a$: a slice through the bell. The denominator does not depend on $x_2$; it only rescales the slice so that its area is one. So the conditional density has the shape of the slice. Along the slice, the exponent of the joint density is a quadratic function of $x_2$, because fixing one variable of a quadratic in two variables leaves a quadratic in the other. The slice is therefore itself a bell, a Gaussian by Equation (4.2).

Working out that quadratic gives the slice's mean and variance. With the covariance Equation (4.8) and means $\mu_1, \mu_2$,

Exercise 4.2 obtains this from the general formula below. Each part has a reading. The observation enters through its z-score $(a - \mu_1)/\sigma_1$. The mean of $x_2$ moves away from $\mu_2$ by $\rho$ times that many of $x_2$'s own standard deviations: a value of $x_1$ one standard deviation above its mean predicts $x_2$ to lie $\rho$ standard deviations above its mean. The variance shrinks by the factor $1 - \rho^2$, the fraction of $x_2$'s variance that $x_1$ does not explain, and it does not depend on $a$ at all.

Drag the slice from left to right. The conditional density slides along the dashed line but keeps its width. The width depends on the correlation, never on the value observed.

Move $\rho$ toward $\pm 0.95$. The conditional density collapses, and the readout shows the fraction of variance removed, $\rho^2$, passing 90%. At $\rho = 0$ nothing is removed and the conditional equals the marginal: an observation uncorrelated with $x_2$ teaches nothing about it.

Compare the dashed line with the ellipses. The line of conditional means is not the long axis of the ellipses; it is flatter. It passes through the leftmost and rightmost points of every ellipse, because along a vertical slice the density is highest where the slice just touches an ellipse.

That flattening has a long history. Francis Galton noticed that the children of unusually tall parents were, on average, less unusual than their parents, and called the effect regression toward mediocrity (Galton, 1886); the statistical term "regression" descends from it. In Equation (4.14) with equal standard deviations, the predicted deviation of $x_2$ is $\rho$ times the observed deviation of $x_1$, closer to the mean whenever $\lvert\rho\rvert < 1$. Nothing pulls the children back. The shrinkage is what conditioning a correlated Gaussian does.

4.5.2 The general formula #

The same reasoning works in any number of dimensions, with blocks in place of numbers. Partition the vector into an observed block $\mathbf{a}$ and an unobserved block $\mathbf{b}$ as in Equation (4.12). Then

Compare the two-dimensional case. The matrix $\mSigma_{ab}^\T\mSigma_{aa}^{-1}$ plays the role of $\rho\sigma_2/\sigma_1 = \Sigma_{12}/\Sigma_{11}$, and the subtracted term $\mSigma_{ab}^\T\mSigma_{aa}^{-1}\mSigma_{ab}$ that of $\rho^2\sigma_2^2 = \Sigma_{12}^2/\Sigma_{11}$.

The dimensions in this formula are easy to misread, so consider a real-sized case. A machine learning model has six hyperparameters to tune, and 40 training runs have finished. In a Bayesian optimizer, $\mathbf{a}$ holds the 40 observed validation errors and $\mathbf{b}$ the unknown errors at, say, 1,000 untried configurations, so $\mSigma_{aa}$ is $40 \times 40$ and $\mSigma_{ab}$ is $40 \times 1000$. The number six appears nowhere: the Gaussian lives over function values, one coordinate per configuration, and the dimension of the input space enters only through the covariances that a kernel assigns to pairs of configurations (Chapter 9). The cost of conditioning grows with the number of evaluations, not with the number of hyperparameters, though in more dimensions more evaluations are needed to pin the function down (Chapter 30). This is the computation that Snoek et al. (2012) used to tune latent Dirichlet allocation, structured support vector machines, and convolutional networks, reaching or surpassing the settings chosen by human experts; Chapter 22 works through such a tuning problem on real data.

The derivation needs one fact from Section 3.7, restated in the notation of Equation (4.12). The Schur complement of the block $\mSigma_{aa}$ is

and the precision matrix $\bm{\Lambda} = \mSigma^{-1}$ has the blocks

Multiplying $\mSigma$ by this matrix gives the identity, block by block, which is how the formula is checked (Petersen and Pedersen, 2012). Only the bottom row is needed: $\bm{\Lambda}_{bb} = \mathbf{S}^{-1}$ and $\bm{\Lambda}_{ab}^\T = -\mathbf{S}^{-1}\mSigma_{ab}^\T\mSigma_{aa}^{-1}$.

The derivation follows section 2.3.1 of Bishop (2006), and Section B.3 collects the result with related identities. Two remarks follow from it. The conditional covariance is the Schur complement Equation (4.16) itself. And step 5 shows that the conditional precision is the block $\bm{\Lambda}_{bb}$ of the joint precision: conditioning reads off a block of the precision matrix, just as marginalizing reads off a block of the covariance matrix.

The second half of the key idea is the reason a Gaussian process's uncertainty depends on where we evaluated and not on what we found (Section 8.1). A second derivation makes the Schur complement less mysterious.

A small numerical case shows the whole of Gaussian process regression in miniature.

4.5.3 Computing it #

In code, Equation (4.15) is evaluated with the Cholesky factor of $\mSigma_{aa}$ and triangular solves, never with an explicit inverse (Section 3.5). With $\mL\mL^\T = \mSigma_{aa}$, set $\mathbf{V} = \mL^{-1}\mSigma_{ab}$ and $\vw = \mL^{-1}(\mathbf{a} - \vmu_a)$. Since $\mSigma_{aa}^{-1} = \mL^{-\T}\mL^{-1}$, the mean offset is $\mathbf{V}^\T\vw$ and the subtracted covariance is $\mathbf{V}^\T\mathbf{V}$.

import numpy as np

def condition(mu, Sigma, ia, ib, a):
    """Mean and covariance of x[ib] given x[ia] = a, for x ~ N(mu, Sigma)."""
    Saa = Sigma[np.ix_(ia, ia)]
    Sab = Sigma[np.ix_(ia, ib)]
    Sbb = Sigma[np.ix_(ib, ib)]
    L = np.linalg.cholesky(Saa)
    V = np.linalg.solve(L, Sab)          # L^{-1} Sigma_ab
    w = np.linalg.solve(L, a - mu[ia])   # L^{-1} (a - mu_a)
    return mu[ib] + V.T @ w, Sbb - V.T @ V

When the observed coordinates nearly determine the others, the subtraction $\mSigma_{bb} - \mathbf{V}^\T\mathbf{V}$ cancels most of its digits, and the result can come out slightly asymmetric or with tiny negative eigenvalues. Symmetrizing it and adding a small jitter before factorizing it again, for example to draw samples, repairs the damage.

4.6 Sums and products #

Two more operations combine Gaussians, and they are easy to confuse because both take two Gaussians and return one. Adding two independent Gaussian random variables describes a quantity built from two uncertain parts, such as a function value plus measurement noise. Multiplying two Gaussian densities describes two independent pieces of evidence about one quantity, which is what Bayes' rule does with a Gaussian prior and a Gaussian likelihood. The first operation makes the uncertainty larger; the second makes it smaller.

4.6.1 Sums of independent Gaussians #

If $X \sim \N(\mu_1, \sigma_1^2)$ and $Y \sim \N(\mu_2, \sigma_2^2)$ are independent, the pair $(X, Y)$ is jointly Gaussian with a diagonal covariance, and $X + Y$ is the linear map with the single row $(1, 1)$. By Equation (4.9),

Variances add; standard deviations do not. Two independent errors with standard deviations 3 and 4 add up to an error with standard deviation 5, not 7. If $X$ and $Y$ are correlated, the same map gives the variance $\sigma_1^2 + \sigma_2^2 + 2\Cov[X, Y]$. Two uses recur in the book. An observation $y = f(\vx) + \varepsilon$ with independent noise $\varepsilon \sim \N(0, \sigma_n^2)$ has variance $\Var[f(\vx)] + \sigma_n^2$, the predictive variance of Section 8.3. And the average of $n$ independent measurements, each with variance $\sigma^2$, has variance $\sigma^2/n$, so its standard deviation falls like $1/\sqrt{n}$.

4.6.2 Products of Gaussian densities #

Now take two Gaussian densities over the same variable $x$ and multiply them pointwise. The result is not normalized, but its shape is Gaussian: the sum of two quadratic exponents is quadratic.

In $d$ dimensions the same steps, with Equation (4.6) in place of Equation (4.2), give

with the constant $Z = \N(\vmu_1;\, \vmu_2,\, \mSigma_1 + \mSigma_2)$ (Rasmussen and Williams, 2006, app. A.2).

Read the result in terms of precision. The precisions add, so the product is narrower than either factor. The mean is a weighted average of the two means, with weights proportional to the precisions, so the sharper density pulls harder. The constant $Z$ is the area under the raw product. It is large when the two densities agree and tiny when they put their mass in different places.

This is Bayes' rule for a Gaussian prior and a Gaussian measurement. Suppose a prior belief about a quantity $\theta$ is $\N(\mu_0, \sigma_0^2)$, and we observe $y = \theta + \varepsilon$ with noise $\varepsilon \sim \N(0, \sigma_n^2)$. As a function of $\theta$, the likelihood $\N(y;\, \theta, \sigma_n^2)$ is the density $\N(\theta;\, y, \sigma_n^2)$, since the formula is symmetric in $y$ and $\theta$. The posterior is proportional to prior times likelihood, so it is Gaussian with precision $1/\sigma_0^2 + 1/\sigma_n^2$ and a mean between the prior mean and the observation. The constant $Z = \N(y;\, \mu_0, \sigma_0^2 + \sigma_n^2)$ is the density of the observation under the prior, which Section 5.6 calls the model evidence. Notice that it is the distribution of the sum $\theta + \varepsilon$ from Equation (4.18): the two operations of this section meet. Section 5.1 develops this view, and Section 5.4 extends it from one number to a vector of weights.

Switch between the two operations with the same inputs. The sum is centered at $\mu_1 + \mu_2$ and is wider than either input. The product lies between $\mu_1$ and $\mu_2$ and is narrower than either.

In the product, make one input very wide. With $\sigma_1$ at 2.5 the product nearly coincides with the second input. A vague prior hardly changes what a precise measurement says.

Pull the two means apart. The normalized product keeps its width, because the precisions do not depend on the means, but the raw product sinks toward zero, and $Z$ with it. Two confident densities that disagree produce a confident compromise in a region where neither puts much mass. A small $Z$ is the warning sign that the prior and the measurement are in conflict.

4.7 Working with Gaussians in code #

Three habits avoid most numerical trouble.

Work with log densities. In many dimensions the density Equation (4.5) is a product of many small factors and underflows. For a standard Gaussian in $d = 1000$ dimensions, even the density at the mean, $(2\pi)^{-500}$, is about $10^{-399}$, below the smallest positive double-precision number. The log density is a sum of moderate numbers. Compute it from the Cholesky factor $\mL$ of $\mSigma$: with $\vw = \mL^{-1}(\vx - \vmu)$ the quadratic form is $\vw^\T\vw$, and $\log\lvert\mSigma\rvert = 2\sum_i \log L_{ii}$ (Section 3.6).

import numpy as np

def gaussian_logpdf(x, mu, Sigma):
    L = np.linalg.cholesky(Sigma)
    w = np.linalg.solve(L, x - mu)       # L^{-1} (x - mu)
    d = len(mu)
    return -0.5 * w @ w - np.log(np.diag(L)).sum() - 0.5 * d * np.log(2 * np.pi)

Never form $\mSigma^{-1}$. Every formula in this chapter that contains an inverse is evaluated with a Cholesky factorization and triangular solves, as in the conditioning code above, for the reasons of Section 3.5.1.

Check the parameterization. The notation $\N(\mu, \sigma^2)$ puts the variance second, but most libraries take the standard deviation: NumPy's random.normal(loc, scale), SciPy's stats.norm(loc, scale), and PyTorch's Normal(loc, scale) all expect $\sigma$. Passing a variance where a standard deviation is expected is a silent and common bug. Multivariate versions take the covariance matrix, or sometimes its Cholesky factor, as in PyTorch's MultivariateNormal(loc, scale_tril=L).

4.8 Exercises #