The Linear Algebra of Uncertainty

Section 2.6.3 introduced a table: the covariances between every pair of uncertain quantities, the covariance matrix. A Bayesian optimizer keeps such a table for the objective's values at every input it has evaluated and every input it is considering, so after 50 evaluations the table has 50 rows at the least. Each update of its belief, each prediction, and each fit of the model's own settings is a computation with that table. This chapter is about those computations.

Linear algebra is a large subject, and this chapter takes only the part that the rest of the book uses. It treats a matrix as a map that moves the points of space, because that picture explains the properties that matter: which tables can be covariance matrices, what shape a covariance gives a cloud of uncertainty, how to solve the equations of a Gaussian process without ever inverting a matrix, and what a determinant measures. Two-dimensional pictures carry most of the intuition, but a few facts appear only from three dimensions on, so the chapter also turns a three-dimensional covariance in your hands.

The payoff comes in the next chapter, Chapter 4, where every result here becomes a statement about the Gaussian distribution, and in Chapter 8, where those statements become an algorithm.

3.1 Vectors #

A vector is an ordered list of numbers. Bold lowercase letters name vectors, and by convention a vector is a column:

The symbol $\R^d$ is the set of all such lists of $d$ real numbers, and $d$ is the vector's dimension. Writing $\vx^\T = (x_1, \dots, x_d)$, read "$\vx$ transpose", turns the column into a row.

Vectors appear in this book in two roles, and keeping them apart prevents a common confusion. In the first role a vector is a point in the space being searched: a configuration of hyperparameters, the settings of an exoskeleton, the parameters of a photo filter. Snoek et al. (2012), for instance, tuned nine settings of an image classifier's training procedure at once, so each configuration they tried was a vector in $\R^9$ (Section 9.2 returns to what they found). In the second role a vector holds the values of the objective at a list of inputs, $\vf = (f(\vx_1), \dots, f(\vx_n))^\T$. Its dimension is the number of inputs, not the number of hyperparameters, and it grows with every evaluation. The covariance matrices of a Gaussian process live in this second space.

Vectors add entry by entry, and multiplying a vector by a number multiplies every entry. Geometrically, $\vx + \vy$ places the arrow of $\vy$ at the tip of the arrow of $\vx$, and $2\vx$ is the same arrow twice as long.

3.1.1 Inner products, lengths, and angles #

The inner product, or dot product, of two vectors of the same dimension multiplies them entry by entry and adds the results:

The notation reads as a row times a column, which is the rule for multiplying matrices introduced in Section 3.2. From the inner product come length and distance. The norm of a vector, its length, is $\lVert\vx\rVert = \sqrt{\vx^\T\vx}$, Pythagoras's theorem in $d$ dimensions, and the distance between two points is the norm of their difference, $\lVert\vx - \vy\rVert$.

The inner product also measures angle. For any two vectors,

where $\theta$ is the angle between them. Two vectors are orthogonal, perpendicular, when their inner product is zero. Dividing by the two lengths gives $\cos\theta$, which is 1 for vectors pointing the same way, 0 for perpendicular ones, and $-1$ for opposite ones. The correlation of Section 2.6.3 is this cosine, computed for centered random variables instead of lists of numbers.

Distance matters to Bayesian optimization because the most common kernels, the functions that decide how strongly the objective's values at two inputs are correlated (Chapter 9), depend on the inputs only through the distance between them. The RBF kernel, $k(\vx, \vx') = \exp\!\left(-\lVert\vx - \vx'\rVert^2 / 2\ell^2\right)$, says that nearby inputs have similar values and distant ones are unrelated, with the lengthscale $\ell$ setting what "nearby" means.

3.1.2 Distances in many dimensions #

Our intuition about distance comes from two and three dimensions, and it misleads in more. Scatter points uniformly at random in the unit cube $[0, 1]^d$. Each coordinate contributes on average $1/6$ to the squared distance between two points, so the typical distance grows like $\sqrt{d/6}$: about 0.58 in two dimensions, 1.0 in six, and 2.9 in fifty. The spread of the distances does not grow with it. In a simulation with 200 random points, the standard deviation of the pairwise distances stayed near 0.24 from two dimensions to a hundred. As a result, the ratio of each point's nearest neighbor distance to its farthest neighbor distance, which was 0.04 for a typical point in two dimensions, rose to 0.22 in six dimensions, 0.65 in fifty, and 0.74 in a hundred.

In many dimensions, then, every point is far from every other, and all of them are about equally far. A kernel with a lengthscale suited to a two-dimensional problem would declare every pair of points in a fifty-dimensional one unrelated, and the model would learn nothing from one evaluation about any other. This is one root of the difficulty of Bayesian optimization in many dimensions, and the remedies, such as lengthscales that grow with the dimension and a lengthscale for each input (Section 9.2), start from it (Chapter 30).

3.2 Matrices as transformations #

A matrix is a rectangular table of numbers, written with a bold capital. An $m \times n$ matrix $\mA$ has $m$ rows and $n$ columns, and $A_{ij}$ is the entry in row $i$ and column $j$. The table is the matrix's storage format. Its meaning is the map it defines.

3.2.1 The product of a matrix and a vector #

Multiplying an $m \times n$ matrix by a vector of dimension $n$ gives a vector of dimension $m$, whose $i$-th entry is the inner product of row $i$ with the vector:

The same product has a second reading, which is the one to picture. Write the columns of $\mA$ as $\mathbf{a}_1, \dots, \mathbf{a}_n$. Then

The vector $\mathbf{e}_1 = (1, 0, \dots, 0)^\T$ picks out the first column: $\mA\mathbf{e}_1 = \mathbf{a}_1$. So the columns of a matrix record where it sends the coordinate directions, and Equation (3.2) says that every other vector goes where its coordinates say: as much of $\mathbf{a}_1$ as $\vx$ had of $\mathbf{e}_1$, and so on.

A map of this kind is linear: it sends sums to sums and multiples to multiples, $\mA(\vx + \vy) = \mA\vx + \mA\vy$ and $\mA(c\vx) = c\,\mA\vx$. In the plane, a linear map sends straight lines to straight lines, keeps the origin in place, and turns a square grid into a grid of parallelograms. The unit circle, the set of vectors of length one, becomes an ellipse.

Some things to try:

• Drag the tip of $\mA\mathbf{e}_1$. Only the first column changes, and the whole grid follows: every point moves by its first coordinate times the change.
• Choose the Rotation, Shear, and Stretch presets. A rotation turns the circle without changing its shape; a stretch scales the two axes by different amounts; a shear slides one axis along the other and turns the square into a slanted parallelogram of the same area.
• Choose Singular. The two columns point the same way, the ellipse collapses to a line segment, and the whole plane is flattened onto a line. Every point on that line came from infinitely many points, so the map cannot be undone.

3.2.2 Composition, identity, and inverse #

Applying $\mathbf{B}$ and then $\mA$ is again a linear map, and its matrix is the product $\mA\mathbf{B}$, whose columns are $\mA$ applied to the columns of $\mathbf{B}$. Entry by entry, $(\mA\mathbf{B})_{ij} = \sum_k A_{ik} B_{kj}$, the inner product of row $i$ of $\mA$ with column $j$ of $\mathbf{B}$. Order matters: stretching and then rotating is not the same as rotating and then stretching, so in general $\mA\mathbf{B} \ne \mathbf{B}\mA$. Multiplying two $n \times n$ matrices this way takes $n^3$ multiplications, one for each combination of $i$, $j$, and $k$.

The identity matrix $\mI$, with ones on the diagonal and zeros elsewhere, leaves every vector where it is. A square matrix $\mA$ is invertible when another matrix $\mA^{-1}$ undoes it, $\mA^{-1}\mA = \mA\mA^{-1} = \mI$. A matrix is invertible exactly when it does not flatten space: when no nonzero vector is sent to zero, so no two points land on the same image. The Singular preset of Figure 3.1 is the opposite case. Inverting a product reverses the order, $(\mA\mathbf{B})^{-1} = \mathbf{B}^{-1}\mA^{-1}$, since the last map applied must be the first one undone.

3.2.3 Transpose #

The transpose $\mA^\T$ flips a matrix across its diagonal, so that $(\mA^\T)_{ij} = A_{ji}$ and an $m \times n$ matrix becomes $n \times m$. A column vector is an $n \times 1$ matrix, its transpose a $1 \times n$ row, and the inner product $\vx^\T\vy$ is the product of a row with a column. The reverse product $\vx\vy^\T$, a column times a row, is an $n \times n$ matrix called the outer product. Transposing a product also reverses the order, $(\mA\mathbf{B})^\T = \mathbf{B}^\T\mA^\T$.

The shapes in Chapter 8 follow these rules. With $n$ observed inputs and $m$ inputs to predict at, the kernel matrix $\mK$ among observed inputs is $n \times n$, the matrix $\mK_*$ between observed and new inputs is $n \times m$, and the posterior mean $\mK_*^\T\mK^{-1}\vy$ multiplies an $m \times n$ matrix, an $n \times n$ matrix, and an $n$-vector to give one prediction for each of the $m$ new inputs.

3.3 Symmetric and positive definite matrices #

Not every square table of numbers can be a covariance matrix. This section finds the condition, which turns out to be geometric: a covariance matrix must not let any direction have negative variance.

3.3.1 Quadratic forms #

A matrix is symmetric when it equals its transpose, $A_{ij} = A_{ji}$. A covariance matrix is symmetric because $\Cov[x_i, x_j] = \Cov[x_j, x_i]$, and so is a kernel matrix, because $k(\vx, \vx') = k(\vx', \vx)$.

For a symmetric matrix $\mA$, the expression $\vw^\T\mA\vw$ is a single number for each vector $\vw$, called a quadratic form. In two dimensions, with $\mA = \begin{bmatrix} a & b \\ b & d \end{bmatrix}$,

a quadratic polynomial in the entries of $\vw$. Quadratic forms appear whenever a covariance matrix meets a weighted sum, because of the following identity.

3.3.2 Every direction has nonnegative variance #

A variance cannot be negative. By the derivation, a covariance matrix must therefore satisfy $\vw^\T\mSigma\vw \ge 0$ for every $\vw$. This property has a name.

Read with $\mA = \mSigma$, positive definiteness says that every direction in the space has positive variance: no weighted combination of the variables is known exactly. A covariance matrix that is only semidefinite has some combination with zero variance, a combination that is fixed. That happens in a Gaussian process when two observed inputs coincide: the values $f(\vx)$ and $f(\vx')$ are then the same random variable, and $f(\vx) - f(\vx')$ is exactly zero. Every result below assumes the strict version, which is why implementations add a little noise or "jitter" to the diagonal (Section 8.4).

In two dimensions the condition is easy to check. A covariance matrix with standard deviations $\sigma_1, \sigma_2$ and correlation $\rho$ is positive definite exactly when $\sigma_1, \sigma_2 > 0$ and $-1 < \rho < 1$, and every correlation in that range is possible. From three dimensions on, this is no longer enough.

3.3.3 Three correlations that cannot coexist #

Suppose three quantities each have variance 1, the first is correlated $0.8$ with the second and $0.8$ with the third, and the second and third are uncorrelated:

Each correlation is between $-1$ and $1$, and each pair on its own is a valid two-dimensional covariance. Yet the weights $\vw = (-2, 1, 1)^\T$ give

a negative variance for the quantity $x_2 + x_3 - 2x_1$. No three random variables have these correlations. Intuitively, if $x_1$ moves closely with $x_2$ and closely with $x_3$, then $x_2$ and $x_3$ must move together to some extent; they cannot be unrelated. Exercise 3.2 finds how strongly they must be correlated.

This is not a curiosity. In finance, tables of correlations between many assets can fail to be positive semidefinite, and the problem of repairing such a table into the nearest valid correlation matrix has a literature of its own (Higham, 2002). For a Gaussian process it is the reason a kernel cannot be an arbitrary similarity score. A natural-looking rule such as "two inputs are fully correlated if they are closer than 1 and uncorrelated otherwise" fails on the inputs $0$, $0.6$, and $1.2$: the first two are fully correlated, the last two are fully correlated, and the first and last are uncorrelated, an impossible triple of the same kind. The kernels of Chapter 9 are built so that this can never happen: a valid kernel is one whose matrices are positive semidefinite for every choice of inputs, and no other function is (Section 7.2).

Figure 3.2 in the next section lets you try to build such a matrix and see what goes wrong.

3.4 Eigenvectors: the axes a matrix stretches #

The entries of a matrix say little at a glance about what it does. In Figure 3.1, most vectors $\vx$ are both turned and stretched by $\mA$. Some directions, though, are only stretched: the image $\mA\vx$ points along $\vx$ itself. If we could find those directions, we could describe the map as a set of independent stretches.

An eigenvector is a direction the matrix does not turn; its eigenvalue is the factor by which the direction is stretched, and a negative eigenvalue reverses it. In Figure 3.1, drag $\vx$ around the circle until $\vx$ and $\mA\vx$ line up, and the readout reports the eigenvalue. For the default matrix this happens near $18°$, with eigenvalue 1.4, and near $135°$, with eigenvalue 0.6. Turn on Show eigenvectors to see both directions. The Rotation preset has no real eigenvectors at all, because it turns every direction.

3.4.1 Symmetric matrices #

For a general matrix the eigenvectors need not be perpendicular, as the default matrix shows. Symmetric matrices, the only kind a covariance can be, are much better behaved.

The proof that the eigenvalues are real takes more space than it deserves here (Strang, 2016), but the perpendicularity is short.

The decomposition reads from right to left as a recipe. $\mathbf{U}^\T$ turns space so that the eigenvectors line up with the coordinate axes, $\bm{\Lambda}$ stretches each axis by its eigenvalue, and $\mathbf{U}$ turns space back. A symmetric matrix is a set of perpendicular stretches, nothing more. Choose the Symmetric preset in Figure 3.1 and turn on the eigenvectors: they are perpendicular, and they are the axes of the ellipse.

3.4.2 Eigenvalues and positive definiteness #

The spectral theorem turns the definition of positive definiteness into a condition on the eigenvalues.

For a covariance matrix the derivation says more. Take $\vw$ of length one, a direction. Its coordinates along the eigenvectors then satisfy $\sum_i c_i^2 = 1$, so the variance in direction $\vw$ is a weighted average of the eigenvalues, with weights $c_i^2$. The eigenvalues are the variances along the eigenvectors, and every other direction's variance lies between the smallest and the largest. The eigenvector with the largest eigenvalue is the direction in which the uncertainty is widest. The shape of the uncertainty is an ellipsoid whose axes point along the eigenvectors, with half-lengths $\sqrt{\lambda_i}$ (the standard deviations along the axes), which Section 4.2.1 derives for the Gaussian.

3.4.3 A covariance in three dimensions #

In two dimensions the ellipse summarizes everything. In three, two things appear that the plane cannot show. Figure 3.2 draws the ellipsoid of a three-dimensional covariance, $\{\vx : \vx^\T\mSigma^{-1}\vx = 1\}$, the set of points one standard deviation from the center in the sense that Section 4.2.1 makes precise, together with its three axes and its shadows on the walls of the cube. Beside it are the three pairwise views, the ellipses of the $2 \times 2$ blocks of $\mSigma$ that involve two of the three coordinates.

Some things to try:

• Turn the cube. The ellipsoid has three perpendicular axes, the eigenvectors, and from most angles none of them is aligned with a coordinate axis. The readout lists the eigenvalues; their square roots are the axes' half-lengths.
• Compare the shadows with the pairwise views. They are the same ellipses. The shadow of a covariance's ellipsoid on the plane of two coordinates is the ellipse of the corresponding $2 \times 2$ block, which is why ignoring a variable of a Gaussian amounts to deleting its row and column (Section 4.4).
• Press An impossible triple. The pairwise views show three valid ellipses, with correlations 0.8, 0.8, and 0. The ellipsoid vanishes, and the red dashed line marks the direction whose variance would be negative, the eigenvector of the eigenvalue $-0.13$. Now raise $\rho_{23}$ slowly: the ellipsoid reappears once $\rho_{23}$ passes 0.28 (at 0.29 on the slider), as a flat disk that thickens as $\rho_{23}$ grows.
• Press A chain. Here $x_1$ is linked to $x_2$, and $x_2$ to $x_3$, with correlation 0.8 each, and $\rho_{13} = 0.64$ is the product of the two. The ellipsoid becomes a cigar along the diagonal, with one eigenvalue much larger than the others: most of the uncertainty lies along a single direction.

The first lesson is that pairwise checks do not certify a covariance matrix; its validity is a property of the whole table. The second is that a covariance can be nearly flat in a direction that no single variable reveals: none of the pairwise views of a nearly singular matrix need look degenerate.

3.4.4 The eigenvalues of a kernel matrix #

That second lesson is the everyday situation of a Gaussian process. The kernel matrix of 100 inputs spaced evenly in $[0, 1]$, under an RBF kernel with lengthscale 0.1, is a $100 \times 100$ positive definite matrix in exact arithmetic. Its largest eigenvalues are 23.9, 21.2, and 17.5, but they decay so quickly that only 28 of the 100 exceed $10^{-10}$. In double-precision arithmetic the smallest ones are lost in rounding error; in our run with NumPy the smallest computed eigenvalue even came out slightly negative, $-4 \times 10^{-15}$. The function values at nearby inputs are so strongly correlated that most directions in the 100-dimensional space have almost no variance. In the same run, NumPy's Cholesky factorization of this matrix failed; adding $10^{-6}$ to the diagonal, the jitter of Section 8.4, makes it succeed. A smooth prior over functions has, in effect, far fewer than 100 independent directions of uncertainty, which is part of why a Gaussian process can learn a smooth function from a few evaluations.

The ratio of the largest to the smallest eigenvalue, the condition number, measures how close to singular a matrix is. Double-precision numbers carry about 16 significant decimal digits, and solving a system with a matrix of condition number $10^k$ can lose about $k$ of them (Golub and Van Loan, 2013). Kernel matrices of smooth kernels on closely spaced inputs routinely reach condition numbers near $10^{16}$, where no digits remain, so jitter or observation noise is a numerical necessity, not a modeling choice.

Eigenvalues of exactly zero carry information too. In Section 18.5, the eigenvalues of a matrix built from a set of pairwise comparisons count the directions in which the comparisons say nothing at all about a person's preferences.

3.5 Solving systems without inverting #

The formulas of Gaussian process regression are full of inverse matrices: the posterior mean $\vk(\vx)^\T\mK^{-1}\vy$ and the variance $k(\vx, \vx) - \vk(\vx)^\T\mK^{-1}\vk(\vx)$ of Equation (8.3). Read literally, they say: invert $\mK$, then multiply. Numerical practice says never to do that. This section explains why, and what to do instead.

3.5.1 Why not invert #

A formula containing $\mK^{-1}\vy$ never needs $\mK^{-1}$ itself. It needs the vector $\bm{\alpha}$ that solves the linear system $\mK\bm{\alpha} = \vy$. Two reasons favor solving over inverting. Forming the inverse costs more, several times as much as the factorization used below. And the computed inverse is less accurate: multiplying by a rounded inverse loses more digits than solving directly with a factorization (Golub and Van Loan, 2013). When the matrix is badly conditioned, as kernel matrices often are, the difference decides whether the answer is usable.

3.5.2 Triangular systems #

Some systems are easy. A matrix is lower triangular when every entry above the diagonal is zero. The system $\mL\vz = \mathbf{b}$ with such a matrix can be solved one entry at a time, top to bottom: the first equation involves only $z_1$, the second only $z_1$ and $z_2$, and so on.

This forward substitution costs about $n^2$ operations, against $n^3$ for general methods. An upper triangular system, with zeros below the diagonal, is solved the same way from the bottom up, by back substitution. The whole strategy for symmetric positive definite systems is to reduce them to two triangular ones.

3.5.3 The Cholesky factorization #

The factor is a square root of the matrix. For a $1 \times 1$ matrix $[a]$ it is $[\sqrt{a}]$, which exists exactly when $a > 0$. The $2 \times 2$ case shows how the general algorithm works and why it needs positive definiteness.

The quantity $d - b^2/a$ in step 4 is a first glimpse of the Schur complement of Section 3.7. For a covariance matrix with $a = \sigma_1^2$, $b = \rho\sigma_1\sigma_2$, and $d = \sigma_2^2$, it equals $\sigma_2^2(1 - \rho^2)$: the variance of the second variable that remains once the first is known, as Section 4.5 will show.

The general algorithm fills in $\mL$ one column at a time, in the same way.

Step 2 makes the algorithm a test as well as a factorization: it succeeds exactly when the matrix is positive definite, and it is the cheapest such test in practice. In Figure 3.2, the impossible triple fails at the third column, where $s$ comes out negative. The algorithm costs about $n^3/3$ floating-point operations (Golub and Van Loan, 2013).

3.5.4 Solving with the factor #

With $\mA = \mL\mL^\T$, the system $\mA\vx = \mathbf{b}$ splits into two triangular systems. First solve $\mL\vz = \mathbf{b}$ by forward substitution, then $\mL^\T\vx = \vz$ by back substitution. Then $\mA\vx = \mL(\mL^\T\vx) = \mL\vz = \mathbf{b}$, as required. Once $\mL$ is known, each new right-hand side costs only $O(n^2)$. This is Algorithm 8.1: one factorization of the kernel matrix, then triangular solves for the weights $\bm{\alpha}$ and for each predictive variance.

import numpy as np
from scipy.linalg import cho_factor, cho_solve, solve_triangular

A = np.array([[4.0, 2.0], [2.0, 3.0]])
b = np.array([2.0, 4.0])

L = np.linalg.cholesky(A)                  # lower triangular, L @ L.T == A
z = solve_triangular(L, b, lower=True)     # forward substitution
x = solve_triangular(L.T, z, lower=False)  # back substitution
print(x)                                   # [-0.25  1.5 ]

c = cho_factor(A)                          # the same, packaged
print(cho_solve(c, b))                     # [-0.25  1.5 ]

3.5.5 The cost, in practice #

The cubic cost decides how many observations an exact Gaussian process can handle. Table 3.1 gives timings measured with NumPy in double precision on the laptop used to prepare this chapter (an Apple M5 Pro); other machines will differ, but the growth rate will not.

Doubling $n$ multiplies the time by roughly eight, as $n^3$ predicts, and the inverse costs about five times the factorization. A Bayesian optimization run rarely has more than a few hundred observations, so a single factorization is cheap. Fitting the kernel's own settings, such as its lengthscale, is what multiplies the cost, because every step of that fit refactorizes the matrix (Section 9.3). For larger data sets, the library GPyTorch replaces the factorization by conjugate gradients, an iterative method that needs only products of the kernel matrix with vectors, and runs them on graphics processors; its authors report that this reduces the asymptotic cost of exact inference from $O(n^3)$ to $O(n^2)$ (Gardner et al., 2018).

3.5.6 The square root that makes samples #

The Cholesky factor is also a recipe for building correlated randomness out of independent randomness. If $\vz$ has independent coordinates with variance 1, then $\mL\vz$ has covariance $\mL\mL^\T = \mSigma$, a consequence of the rule $\Cov[\mA\vz] = \mA\Cov[\vz]\mA^\T$ that Section 4.3 derives. Turn on Show samples in Figure 3.2 to see 160 such points: each is $\mL\vz$ for a draw $\vz$ of three independent random numbers with mean 0 and variance 1 (standard normal numbers, the bell curve of Section 4.1.1), and together they fill the ellipsoid's shape. Because $\mL$ is lower triangular, the first coordinate uses only $z_1$, the second $z_1$ and $z_2$, and the third all three, so each new coordinate is built from the randomness of the earlier ones plus a fresh part of its own. This is how the posterior samples of Section 8.5 are drawn.

3.6 Determinants as volume #

The density of a multivariate Gaussian (Section 4.2) and the marginal likelihood used to fit kernels (Section 9.3) both contain a determinant. Its meaning is geometric.

The determinant of a square matrix, $\det\mA$ or $\lvert\mA\rvert$, is the factor by which the map $\mA$ multiplies areas in two dimensions, volumes in three, and $d$-dimensional volumes in general. Its sign records whether the map flips orientation, as a mirror does. For a $2 \times 2$ matrix,

the signed area of the parallelogram spanned by the two columns. This is the yellow parallelogram in Figure 3.1, the image of the unit square; the Singular preset flattens it to zero area, and the Reflection preset turns the determinant negative.

Three properties follow from the picture. Applying two maps multiplies their volume factors, so $\det(\mA\mathbf{B}) = \det\mA \cdot \det\mathbf{B}$, and in particular $\det(\mA^{-1}) = 1/\det\mA$. A matrix is invertible exactly when its determinant is nonzero, when it does not flatten space. And for a symmetric matrix $\mA = \mathbf{U}\bm{\Lambda}\mathbf{U}^\T$, the rotations $\mathbf{U}$ and $\mathbf{U}^\T$ preserve volume while $\bm{\Lambda}$ stretches axis $i$ by $\lambda_i$, so

For a covariance matrix the determinant is therefore a single number for the overall size of the uncertainty, sometimes called the generalized variance. The ellipsoid of Figure 3.2 has half-axes $\sqrt{\lambda_i}$ and so volume $\frac{4}{3}\pi\sqrt{\lambda_1\lambda_2\lambda_3} = \frac{4}{3}\pi \sqrt{\det\mSigma}$. Correlation shrinks it: for three variables with unit variance and correlations $\rho_{12}, \rho_{13}, \rho_{23}$,

which falls from 1 for independent variables toward 0 as the ellipsoid flattens, and becomes negative for the impossible triple ($1 - 0.64 - 0.64 = -0.28$). In Section 6.5, a determinant of this kind measures how much a set of evaluations can teach about the objective, and the evaluations that make it largest are the most informative.

3.6.1 The log-determinant from the Cholesky factor #

The determinant of a triangular matrix is the product of its diagonal entries. (In the $2 \times 2$ formula, $ad - bc$ with $b = 0$ is $ad$. In general, a triangular matrix stretches coordinate axis $i$ by its $i$-th diagonal entry and then shears, and a shear moves points parallel to the other axes without changing any volume, as the Shear preset of Figure 3.1 shows.) With $\mA = \mL\mL^\T$, the product rule for determinants gives

The logarithm is not decoration. A kernel matrix with many tiny eigenvalues has a determinant far below the smallest positive double-precision number, about $5 \times 10^{-324}$. For the 100-input example of Section 3.4.4, with the jitter that makes it factorizable, the product of the computed eigenvalues underflows to zero, while Equation (3.6) gives $\log\det\mA \approx -1139$, a moderate number that costs nothing once the factor is known. Every implementation of the Gaussian log density and of the log marginal likelihood computes the determinant this way (Section 4.7).

3.7 Block matrices and the Schur complement #

A Gaussian process works with two groups of variables at once: the function values already observed and those to be predicted. Its matrices are therefore naturally split into blocks, as in Equation (8.1). This section develops the algebra of such blocks and ends with the one formula Section 4.5 needs.

3.7.1 Partitioned matrices #

A symmetric matrix partitioned into blocks looks like

where $\mA$ is $p \times p$, $\mathbf{D}$ is $q \times q$, and $\mathbf{B}$ is $p \times q$. Blocks multiply like numbers, with the rule that the order of each product is kept, since matrices do not commute. For a covariance, $\mA$ and $\mathbf{D}$ are the covariances within each group and $\mathbf{B}$ holds the covariances across groups.

3.7.2 Elimination by blocks #

Solving two equations in two unknowns by hand, we use the first equation to eliminate the first unknown from the second. The same step with blocks produces the central object of this section.

For a $2 \times 2$ matrix with entries $a, b, d$, the Schur complement of $a$ is the number $d - b^2/a$ from the Cholesky derivation. In general it is what remains of $\mathbf{D}$ after the part explained by the first group is removed.

Three consequences follow, each used later in the book.

Determinants. $\mathbf{E}$ is block triangular with identity blocks on the diagonal, so its determinant is 1 and $\det\mathbf{M} = \det\mA \cdot \det\mathbf{S}$.

Positive definiteness. For any vector $\vw$, set $\mathbf{v} = \mathbf{E}^{-\T}\vw$; then $\vw^\T\mathbf{M}\vw = \mathbf{v}_1^\T\mA\mathbf{v}_1 + \mathbf{v}_2^\T\mathbf{S}\mathbf{v}_2$, where $\mathbf{v}_1$ and $\mathbf{v}_2$ are the two parts of $\mathbf{v}$. So $\mathbf{M}$ is positive definite exactly when both $\mA$ and its Schur complement $\mathbf{S}$ are. This is the precise version of the impossible triple: once $x_1$ is accounted for, what is left of $x_2$ and $x_3$ must still be a valid covariance (Exercise 3.2).

The inverse. Inverting step 3 block by block, using $(\mathbf{E}^{-1})^{-1} = \mathbf{E}$, gives $\mathbf{M}^{-1} = \mathbf{E}^\T \begin{bmatrix} \mA^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{S}^{-1} \end{bmatrix} \mathbf{E}$, and multiplying out,

The bottom-right block of the inverse is the inverse of the Schur complement. Section 4.5 uses this, with $\mA$ the covariance of the observed values and $\mathbf{D}$ that of the unobserved ones, to show that the covariance after conditioning is the Schur complement: $\mathbf{D} - \mathbf{B}^\T\mA^{-1}\mathbf{B}$, the prior covariance minus the part the observations explain. Section 8.1 is the same formula with kernel matrices in the blocks, and Section B.1 and Section B.2 collect it with related identities (Petersen and Pedersen, 2012).

3.7.3 Adding one observation #

The block view also says how to update a Cholesky factor cheaply, which a Bayesian optimizer needs after every evaluation. Suppose $\mK = \mL\mL^\T$ for the $n$ inputs observed so far, and a new input arrives with covariances $\vk$ to the old inputs and variance $\kappa$. The new kernel matrix and its factor have the block forms

Matching blocks, $\mL\mathbf{l} = \vk$, one forward substitution costing $O(n^2)$, and $l_{\ast}^2 = \kappa - \mathbf{l}^\T\mathbf{l}$. Since $\mathbf{l}^\T\mathbf{l} = \vk^\T\mL^{-\T}\mL^{-1}\vk = \vk^\T\mK^{-1}\vk$, the new diagonal entry squared is $\kappa - \vk^\T\mK^{-1}\vk$, the Schur complement of $\mK$, which is also the Gaussian process's posterior variance at the new input (Equation (8.3)). Growing the factor by one row costs $O(n^2)$ instead of the $O(n^3)$ of starting over, as long as the kernel's settings stay fixed. And if the new input repeats an old one without observation noise, $l_{\ast}$ is zero and the factorization breaks: the matrix has become singular, for the reason given in Section 3.3.2.

3.8 Exercises #

Further reading #

• Strang (2016) is a patient introduction to matrices as maps, eigenvalues, and positive definite matrices, with the geometric emphasis of this chapter.
• Golub and Van Loan (2013) is the standard reference for how these computations are done in floating point: triangular systems, the Cholesky factorization, operation counts, and conditioning.
• Rasmussen and Williams (2006), appendix A, collects the matrix identities and the Cholesky-based computations that Gaussian process regression uses.
• Petersen and Pedersen (2012) is a compact catalog of identities, including the block inverse and determinant formulas, useful for checking a derivation.
• Sanderson (2016) is an animated video series that shows matrices acting on the plane and in space, the picture behind Figure 3.1.
• Higham (2002) treats the problem of repairing an invalid correlation matrix, the practical side of Section 3.3.3.