The Analysis Behind Kernels

The three chapters before this one built a working model: a kernel says how the values of an unknown function move together (Chapter 7), conditioning turns it into predictions (Chapter 8), and the marginal likelihood picks its settings (Chapter 9). That is enough to run Bayesian optimization. It is not enough to say why it works. The guarantees of Chapter 13 are of two kinds. The Bayesian ones hold for a function drawn from the prior. The frequentist ones, and their counterparts for comparisons in Chapter 21 and Chapter 29, hold for one fixed function whose norm in a reproducing kernel Hilbert space is at most a number $B$. The rates of both depend on how fast the eigenvalues of the kernel decay. These statements treat a kernel not as a recipe for covariance matrices but as an object in its own right: a linear map on functions, with eigenvalues, a spectrum of frequencies, and a space of functions it calls simple.

This chapter builds that view from what the reader already has: inner products and eigenvectors (Chapter 3) and the information a Gaussian model gains from noisy data (Section 6.5). The plan is to treat a function as a very long vector and to carry each finite fact over. It ends where the results are used, with the growth of the maximum information gain and with a precise account of what a Gaussian process is as a random object. None of it is needed to use the methods of Part III; it is needed to read their guarantees.

10.1 Functions as vectors #

A vector $\vx \in \R^n$ is a list of $n$ numbers, which is the same thing as a function from the indices $\{1, \dots, n\}$ to $\R$: give it an index $i$, and it returns $x_i$. A function $f$ on $[0, 1]$ is the same kind of object with a continuum of indices. Section 7.4 used this reading when it drew a function as the vector of its values on a grid.

10.1.1 Inner products of functions #

Take the grid of midpoints $x_i = (i - \tfrac12)/n$ and the vectors $\vf$ and $\vg$ of two functions' values there. Their inner product $\vf^\T\vg$ (Equation (3.1)) grows with $n$, because it adds more terms. Divided by $n$ it settles, a Riemann sum turning into an integral as in Section 7.2.1:

The limit is the inner product of two functions, and length, $\lVert f \rVert = (\int_0^1 f^2\, \dd x)^{1/2}$, and orthogonality follow as in Section 3.1.1. A weighting density $q(x) > 0$ can make some inputs count more: $\langle f, g \rangle_q = \int f(x)\, g(x)\, q(x)\, \dd x$. This chapter uses the uniform weighting $q = 1$ on $[0, 1]$ unless it says otherwise.

10.1.2 Orthonormal bases: Fourier as the example #

An orthonormal basis gives vectors coordinates, and functions too. The best-known basis on $[0, 1]$ is Fourier's: the constant $1$ and the functions $\sqrt{2}\cos(2\pi jx)$ and $\sqrt{2}\sin(2\pi jx)$ for $j = 1, 2, \dots$, each of length 1 and orthogonal to the others. The coordinates of $f$ are its Fourier coefficients $a_j$ (cosine) and $b_j$ (sine), and its squared length is the sum of their squares, Parseval's identity.

Coefficients carry smoothness. For a differentiable $f$ with $f(0) = f(1)$, integrating by parts moves the derivative onto the basis function and brings out a factor $1/(2\pi j)$: the cosine coefficient of $f$ is $-1/(2\pi j)$ times the sine coefficient of $f'$, and the sine coefficient of $f$ is $1/(2\pi j)$ times the cosine coefficient of $f'$. Parseval's identity for $f'$ then gives

The energy of the slope, a natural measure of roughness, is a sum of squared coefficients with weights that grow with frequency. High frequencies are expensive and low ones cheap. The norm of Section 10.2 has this shape, with weights chosen by the kernel.

10.1.3 The kernel matrix as a view of an operator #

A matrix maps vectors to vectors. A kernel maps functions to functions:

a combination of the values of $g$, weighted by how strongly $x$ is correlated with each $x'$. On the grid the integral becomes $\frac1n\sum_j k(x_i, x_j)\, g(x_j)$, the matrix-vector product $\frac1n\mK\vg$. So $\mK/n$ is the operator $\mathcal{K}$ seen through $n$ points, and its eigenvalues approximate the operator's.

Section 3.4.4 met one such matrix: the RBF kernel with lengthscale 0.1 on 100 evenly spaced inputs, with largest eigenvalues 23.9, 21.2, and 17.5. Divided by 100 they are 0.239, 0.212, and 0.175. Those inputs include both ends of $[0, 1]$. Placed instead at the 100 midpoints of Section 10.1.1, as in Figure 10.1 below, they give 0.241, 0.214, and 0.176, and finer grids of midpoints leave these digits unchanged. The rapid decay that broke the Cholesky factorization there belongs to the operator, not to the grid.

10.2 The space a kernel defines #

The frequentist regret theorems of Chapter 13, which bound how much an optimizer loses against the best value it could have found for one fixed function, need a class of functions large enough to contain realistic objectives and small enough that finitely many evaluations can pin a member down. A kernel defines one, and the construction starts from the weight-space view of Chapter 7.

10.2.1 Functions built from bumps #

Take a kernel that comes from features, $k(\vx, \vx') = \boldsymbol{\phi}(\vx)^\T\boldsymbol{\phi}(\vx')$, as in Equation (7.3) with prior weight covariance $\mI$. Each weight vector $\vw$ gives a function $\vw^\T\boldsymbol{\phi}(\vx)$, and the kernel bump $k(\cdot, \vx')$ is the one with weights $\boldsymbol{\phi}(\vx')$. A sum of bumps $f = \sum_i \alpha_i k(\cdot, \vx_i)$ therefore has weights $\sum_i \alpha_i \boldsymbol{\phi}(\vx_i)$, and for a second sum $g = \sum_j \beta_j k(\cdot, \vx'_j)$ the inner product of the two weight vectors is

The right side mentions no features, so it serves as a definition for any positive semidefinite kernel. It does not depend on how $f$ and $g$ are written as sums: grouping by $j$ gives $\sum_j \beta_j f(\vx'_j)$ and grouping by $i$ gives $\sum_i \alpha_i g(\vx_i)$, which depend only on the functions' values. The squared norm of $f$ is $\lVert f \rVert_k^2 = \bm{\alpha}^\T\mK\bm{\alpha} \ge 0$.

10.2.2 The reproducing property #

With $g = k(\cdot, \vx)$, a single bump with coefficient 1, the same grouping gives

Taking the inner product with a bump evaluates the function. This is the reproducing property. With $f = k(\cdot, \vx')$ it gives $\langle k(\cdot, \vx'), k(\cdot, \vx) \rangle_k = k(\vx, \vx')$: the bumps are features for the kernel. The Cauchy-Schwarz inequality, $\lvert\langle f, g\rangle_k\rvert \le \lVert f\rVert_k\lVert g\rVert_k$, applied to Equation (10.4) and to the difference of two bumps, gives two bounds (Chowdhury and Gopalan, 2017):

They say what the norm controls. Under a kernel with $k(\vx, \vx) \le 1$, a function of norm at most $B$ never exceeds $B$ in absolute value. And it cannot change quickly: for the RBF kernel the second square root is $\sqrt{2 - 2e^{-r^2/2\ell^2}}$ at distance $r$, which is 0.0999 at $r = 0.01$ when $\ell = 0.1$, close to $r/\ell$. Over a distance short compared with the lengthscale, such a function changes by at most about $B\,r/\ell$. Large values and fast changes cost norm.

The first bound also shows that only the zero function has norm zero, so $\lVert\cdot\rVert_k$ is a genuine length, and that sums of bumps converging in this norm converge at every input, so their limits are functions too. Adding those limits completes the construction.

For familiar kernels the space is recognizable. The line kernel $k(x, x') = xx'$ (Example 7.1 with $\sigma_0 = 0$, $\sigma_1 = 1$) has the lines $f(x) = wx$ as its RKHS, with norm $\lvert w \rvert$, the slope. For a kernel from finitely many features with weight covariance $\mI$, the norm of $f$ is the length of the shortest weight vector that produces it (Steinwart and Christmann, 2008, ch. 4). The RKHS of a Matérn kernel with smoothness $\nu$ on a bounded domain in $d$ dimensions holds the same functions as a Sobolev space, those whose derivatives up to order $\nu + d/2$ are square-integrable, with an equivalent norm, when $\nu + d/2$ is a whole number and the boundary is regular (Kanagawa et al., 2018, ex. 2.6). The RBF kernel's RKHS holds only functions whose Fourier transforms decay exponentially fast, so they are extremely smooth (Kanagawa et al., 2018, ex. 2.7).

10.2.3 The posterior mean lives in the space #

The posterior mean of Gaussian process regression is a weighted sum of bumps, one per observation (Equation (8.4)), so it lies in $\mathcal{H}_k$. It also solves a problem stated without probability.

That a penalized fit over an infinite-dimensional space has a solution with $n$ terms is the representer theorem, first stated for squared error by Kimeldorf and Wahba (1971); the match with the Gaussian process posterior mean goes back to Kimeldorf and Wahba (1970) (see also Kanagawa et al., 2018, prop. 3.6). With the data term a plain sum, as here, the penalty weight is the noise variance, the kernel ridge regression of Section 8.2.1. In weight space the match is expected: with $\vw \sim \N(\mathbf{0}, \mI)$, minus twice the log posterior is $\sigma_n^{-2}\sum_i(y_i - f(\vx_i))^2 + \lVert\vw\rVert^2$ up to a constant, and the shortest $\vw$ for $f$ has length $\lVert f\rVert_k$. The squared norm plays the part of minus twice the log prior.

The posterior mean of Figure 8.1 at its defaults (lengthscale 0.12, four observations, the largest 0.85 in absolute value) has $\lVert\mu\rVert_k^2 = \bm{\alpha}^\T\mK\bm{\alpha} = 1.39$, so Equation (10.5) caps it at 1.18, above its actual maximum: the norm is a guarantee, not a description.

10.2.4 What a norm bound assumes #

The frequentist regret theorems (Section 13.4.4, Section 21.3, Chapter 29), in which $f$ is one fixed function and the only randomness is the evaluation noise, assume that $f$ lies in $\mathcal{H}_k$ with a known bound on its norm. By the sections above, if $\lVert f\rVert_k \le B$ and $k(\vx, \vx) \le 1$, then $f$ is no larger than $B$ anywhere, changes by at most about $B$ per lengthscale, and puts little weight where the kernel says weight is expensive. The assumption depends on the kernel and its lengthscale, not only on $f$: a longer lengthscale makes fast changes more expensive, and under the RBF kernel a Gaussian bump of width $\ell/\sqrt{2} \approx 0.71\,\ell$ or less has infinite norm (Exercise 10.1).

Two conventions are in use, and they differ by a square. Srinivas et al. (2010) assume $\lVert f \rVert_k^2 \le B$, as Section 13.4.4 does; Chowdhury and Gopalan (2017) and the preference papers of Section 21.3 and Chapter 29 assume $\lVert f\rVert_k \le B$. In the second convention the bound enters the confidence width directly. For noise that is $R$-sub-Gaussian (tails no heavier than those of a Gaussian with standard deviation $R$), Chowdhury and Gopalan show that with probability at least $1 - \delta$, $\lvert\mu_{t-1}(\vx) - f(\vx)\rvert \le \beta_t^{1/2}\sigma_{t-1}(\vx)$ for all $\vx$ and all rounds $t \le T$, with $\beta_t^{1/2} = B + R\sqrt{2(\gamma_{t-1} + 1 + \log(1/\delta))}$ (Chowdhury and Gopalan, 2017, thm. 2); their posterior and their $\gamma_{t-1}$ use the noise variance $1 + 2/T$ in place of $\sigma_n^2$, and their $\beta_t$ is the square root of the book's. Chapter 13 explains where such widths come from. In practice $B$ is unknown, and the kernel that sets its units is fitted from the same data (Section 13.5.2).

10.2.5 Samples are rougher than the space #

The Bayesian theorem of Section 13.4.1 assumes instead that $f$ is a draw from $\GP(0, k)$. It is natural to guess that a typical draw has a moderate norm in $\mathcal{H}_k$. It has none at all.

The general statement is a zero-one law: a Gaussian process lies in a given RKHS with probability 0 or 1, and in the RKHS of its own kernel with probability 0 whenever that space is infinite-dimensional (Driscoll, 1973; Lukić and Beder, 2001; Kanagawa et al., 2018, thm. 4.9 and cor. 4.10). Read through $n$ inputs, a draw looks like a function of squared norm about $n$, and each new input adds about 1. The posterior mean is a finite sum of bumps, and averaging removes the roughness (Rasmussen and Williams, 2006, sec. 6.1).

Draws do lie in slightly larger spaces of rougher functions. For the RBF kernel the difference rarely matters (Kanagawa et al., 2018, cor. 4.13 and remark 4.13). For a Matérn kernel it does: the RKHS asks for Sobolev smoothness $\nu + d/2$, while draws have every order below $\nu$ and no more, rougher by $d/2$ (Kanagawa et al., 2018, cor. 4.15 and remarks 4.14 and 4.15).

So the two settings of Chapter 13 assume different things. The Bayesian theorem is about draws, which no bound $B$ covers; the frequentist theorems are about members of $\mathcal{H}_k$, a set to which the prior gives probability zero. Neither contains the other, as Srinivas et al. (2010) note. Choosing a Matérn 5/2 kernel in $d$ dimensions and then quoting a frequentist bound assumes smoothness $5/2 + d/2$, more than the draws of the same prior have (inference).

10.3 Mercer's theorem #

Bumps overlap and are not orthogonal, so they make awkward coordinates. For a symmetric matrix, Section 3.4.1 found better ones, the eigenvectors. The operator $\mathcal{K}$ is the continuous version of a symmetric matrix, and the same move works.

10.3.1 Eigenfunctions and the expansion of the kernel #

An eigenfunction of $\mathcal{K}$ is a function $\varphi$ that the operator only rescales, $\mathcal{K}\varphi = \lambda\varphi$, with eigenvalue $\lambda$. (This chapter writes $\varphi_i$ for eigenfunctions, to keep them apart from the features $\boldsymbol{\phi}$ of Chapter 7 and the normal density $\phi$.)

Equation (10.6) is the spectral theorem Theorem 3.1, $\mA = \sum_i \lambda_i\mathbf{u}_i\mathbf{u}_i^\T$, with functions in place of vectors. The conditions matter: where the weighting gives no weight, the expansion can fail (Kanagawa et al., 2018, remark 4.2). The eigenvalues and eigenfunctions depend on the weighting, while the kernel and its RKHS do not (Kanagawa et al., 2018, remarks 4.1 and 4.3). Three consequences follow.

The eigenvalues are a variance budget. Setting $\vx' = \vx$ in Equation (10.6) and integrating against $q$ gives, by orthonormality, $\sum_i \lambda_i = \int k(\vx, \vx)\, q(\vx)\, \dd\vx$. For a stationary kernel and a probability density $q$, the eigenvalues sum to $\sigma_f^2$, the prior variance, and say how it is shared among directions.

Every kernel is an inner product of features. With $\boldsymbol{\phi}(\vx) = (\sqrt{\lambda_1}\varphi_1(\vx), \sqrt{\lambda_2}\varphi_2(\vx), \dots)$, Equation (10.6) reads $k(\vx, \vx') = \boldsymbol{\phi}(\vx)^\T\boldsymbol{\phi}(\vx')$. This is the converse quoted in Section 7.2.2.

The RKHS norm in eigen-coordinates. Expand $f = \sum_i c_i\varphi_i$ with $c_i = \langle f, \varphi_i\rangle_q$. Then $f \in \mathcal{H}_k$ exactly when the following sum is finite, and

(Rasmussen and Williams, 2006, sec. 6.1; Kanagawa et al., 2018, thm. 4.2). The reproducing property checks it: by Equation (10.6) the bump $k(\cdot, \vx)$ has coefficients $\lambda_i\varphi_i(\vx)$, so $\langle f, k(\cdot,\vx)\rangle_k = \sum_i c_i\lambda_i\varphi_i(\vx)/\lambda_i = f(\vx)$. This is Equation (10.1) with weights $1/\lambda_i$ chosen by the kernel, and the continuous form of $\vf^\T\mK^{-1}\vf$ written in the eigenvectors of $\mK$ (Rasmussen and Williams, 2006, sec. 6.1). A direction with a small eigenvalue is expensive: a coefficient $c$ along it costs $c^2/\lambda_i$.

10.3.2 The eigen-expansion of a sample #

The same coordinates describe the prior. With independent standard normal numbers $z_1, z_2, \dots$, set

Then $\E[f(\vx)f(\vx')] = \sum_{i,j}\sqrt{\lambda_i\lambda_j}\,\E[z_iz_j]\,\varphi_i(\vx)\varphi_j(\vx') = \sum_i\lambda_i\varphi_i(\vx)\varphi_i(\vx') = k(\vx, \vx')$, since $\E[z_iz_j]$ is 1 for $i = j$ and 0 otherwise. So Equation (10.8) is a draw from $\GP(0, k)$ written in the eigenbasis, with independent coefficients of variance $\lambda_i$. This is the Karhunen-Loève expansion; the series converges in mean square, uniformly over the domain (Kanagawa et al., 2018, thm. 4.3; Berlinet and Thomas-Agnan, 2004, sec. 2.3). By orthonormality, keeping the first $m$ terms leaves an average squared error of

the eigenvalue mass left out.

The two readings of a small eigenvalue agree: the prior gives the direction little variance, $\lambda_i$, and the norm charges it heavily, $1/\lambda_i$. By Equation (10.7) the truncated draw has $\lVert f_m\rVert_k^2 = \sum_{i \le m} z_i^2$, about $m$, the picture of Section 10.2.5 in eigen-coordinates. This last step is intuition rather than proof, because the full series converges in mean square and not in the norm (Kanagawa et al., 2018, remark 4.9); Section 10.2.5 gave the proof.

10.3.3 Computing them #

Eigenfunctions rarely have a closed form. They are computed as Section 10.1.3 suggested: on a grid of $n$ points with uniform weighting, the eigenvalues of $\mK/n$ estimate the $\lambda_i$, and the eigenvectors times $\sqrt{n}$ estimate the $\varphi_i$ at the grid points. This is the Nyström method; it estimates the larger eigenvalues better than the smaller ones (Rasmussen and Williams, 2006, sec. 4.3.2).

Read the default. The RBF kernel with lengthscale 0.1 has largest eigenvalues 0.241, 0.214, and 0.176, the values of Section 10.1.3; all of its eigenvalues together sum to 1.00, the prior variance, and they reach the rounding floor by the 32nd. The first four eigenfunctions look like cosines with 0, 1, 2, and 3 sign changes.

Move the terms slider. One term holds 24% of the prior variance and ten hold 99.5%; the ten-term draw is hard to tell from the exact one. The squared norm of the truncated draw keeps growing, 8.7 with ten terms and 96.5 with all 100, though the draw barely changes.

Switch to Matérn 1/2. The eigenvalues fall on a straight line on log axes, a power law. Ten terms hold 80% of the variance and twenty about 90%. The truncated draw is smooth and the exact one jagged: the roughness lives in the many small eigenvalues.

Compare the slopes. The Matérn lines have slopes near $-2$, $-4$, and $-6$ for $\nu = 1/2$, $3/2$, and $5/2$; on a grid of 3,000 points, between the 20th and 60th eigenvalues, they are $-2.03$, $-4.00$, and $-5.89$. The eigenvalues decay like $i^{-(2\nu + 1)}$, in line with the result of Ritter and colleagues that a process with $r$ mean-square derivatives on $[0, 1]$ has eigenvalues decaying like $i^{-(2r + 2)}$ (Rasmussen and Williams, 2006, sec. 4.3).

Lengthen the lengthscale to 0.3. The first RBF eigenvalue rises to 0.590 and four terms hold 99.6% of the variance: a longer lengthscale concentrates the prior on fewer directions.

10.4 Bochner's theorem #

Mercer's eigenfunctions depend on the domain and the weighting, and must be computed. For a stationary kernel, which depends only on $\mathbf{r} = \vx - \vx'$ (Section 7.5.1), there is a description that needs neither: a list of frequencies and the variance each carries.

10.4.1 Kernels as spectra #

Start from one cosine with a random amplitude. If $f(x) = a\cos(\omega x) + b\sin(\omega x)$ with $a, b$ independent standard normals, then $\Cov[f(x), f(x')] = \cos\omega x\cos\omega x' + \sin\omega x\sin\omega x' = \cos\big(\omega(x - x')\big)$, a stationary kernel. Mixtures over frequencies give more, and Bochner's theorem says they give all.

For the kernels of this book $\Lambda$ has a density $s(\boldsymbol{\omega})$, the spectral density, symmetric in $\boldsymbol{\omega}$, and $k(\mathbf{r}) = \int s(\boldsymbol{\omega})\cos(\boldsymbol{\omega}^\T\mathbf{r})\,\dd\boldsymbol{\omega}$. Frequencies here are in radians per unit of input; Rasmussen and Williams use cycles, $\boldsymbol{\omega} = 2\pi\mathbf{s}$, which rescales the density but not its shape (Rasmussen and Williams, 2006, eq. 4.6). Since $\int s = k(\mathbf{0}) = \sigma_f^2$, $p = s/\sigma_f^2$ is a probability density, and

One direction of the theorem is short. For inputs $\vx_1, \dots, \vx_n$ and coefficients $a_i$, the identity $\cos(\theta - \theta') = \cos\theta\cos\theta' + \sin\theta\sin\theta'$ gives

If $s$ is positive at every frequency, as for the RBF and Matérn kernels, the integral is positive unless all $a_i$ are zero, because for distinct inputs the two sums cannot vanish together at every frequency. Every kernel matrix of distinct inputs is then invertible, which Section 10.2.5 used (Wendland, 2004, ch. 6). The converse, that every stationary kernel has such a spectrum, is the deep part.

10.4.2 The spectra of the RBF and Matérn kernels #

For the RBF kernel the frequencies are Gaussian, $p = \N(\mathbf{0}, \ell^{-2}\mI)$, so a typical frequency is about $1/\ell$ (Exercise 10.2). For the Matérn kernel they follow a Student-t distribution with $2\nu$ degrees of freedom and scale $1/\ell$,

which is eq. 4.15 of Rasmussen and Williams (2006) in radians. For $\nu = 1/2$ in one dimension it is the Cauchy distribution, $p(\omega) = (\ell/\pi)/(1 + \ell^2\omega^2)$.

The tails differ. The Gaussian falls faster than any power, while Equation (10.11) falls like $\lVert\boldsymbol{\omega}\rVert^{-(2\nu + d)}$, keeping a little variance at every high frequency, more for smaller $\nu$. The tails set the smoothness. Differentiating Equation (10.10) twice at $\mathbf{r} = 0$ in one dimension gives $-k''(0) = \sigma_f^2\int\omega^2 p(\omega)\,\dd\omega$, the variance of the mean-square derivative $f'(x)$ (Section 7.5.3). It is finite only if the tail falls faster than $\lvert\omega\rvert^{-3}$, which for the Matérn kernel means $2\nu + 1 > 3$, or $\nu > 1$. With higher even powers of $\omega$ in place of $\omega^2$ the same argument gives the rule of Table 9.1: a mean-square derivative of every whole order below $\nu$, and none of order $\nu$ or above. And $\sqrt{-k''(0)/k(0)}$, the root-mean-square frequency, is what Equation (7.6) counts: for Matérn 5/2 the Student-t with 5 degrees of freedom has variance $5/3$ in units of $1/\ell^2$, the $k''(0) = -5/(3\ell^2)$ of Exercise 7.4.

Mercer and Bochner describe the same kernel, and on a domain without ends they coincide. Join the ends of $[0, 1]$ into a circle and wrap the kernel around it, $k_\circ(r) = \sum_m k(r + m)$ over all integers $m$. Integrating $k_\circ(x - y)\cos(2\pi jy)$ over one turn equals integrating $k(x - y)\cos(2\pi jy)$ over the whole line: substituting $u = y - m$ turns the integral of the term $k(x - y + m)\cos(2\pi jy)$ over $[0, 1]$ into the integral of $k(x - u)\cos(2\pi ju)$ over $[-m, 1 - m]$, since the cosine repeats with period 1, and these intervals cover the line. With $r = x - y$, the cosine splits as $\cos(2\pi jx)\cos(2\pi jr) + \sin(2\pi jx)\sin(2\pi jr)$. The sine part integrates to zero because $k$ is even, and the cosine part gives $2\pi s(2\pi j)\cos(2\pi jx)$, because the spectral density is recovered from the kernel by $s(\omega) = \frac{1}{2\pi}\int k(r)\cos(\omega r)\,\dd r$ (Rasmussen and Williams, 2006, eq. 4.6). Sines work the same way. So the eigenfunctions on the circle are Fourier's, and the eigenvalues are $\lambda_j = 2\pi s(2\pi j)$: the spectral density sampled at the frequencies that fit around the circle. Eigenvalue decay is spectral decay. Since the $i$th eigenvalue sits near frequency $\pi i$, a Matérn spectrum falling like $\lvert\omega\rvert^{-(2\nu+1)}$ gives eigenvalues falling like $i^{-(2\nu + 1)}$, the slopes of Figure 10.1. The interval's ends change the values but not the tail: at lengthscale 0.1 the 80th Matérn 1/2 eigenvalue is $3.24 \times 10^{-4}$ on the interval and $3.16 \times 10^{-4}$ on the circle.

10.4.3 Random Fourier features #

Equation (10.10) turns a kernel into an expectation, which an average can estimate. Draw frequencies $\boldsymbol{\omega}_1, \dots, \boldsymbol{\omega}_M$ from $p$ and use the $2M$ features

By the cosine identity, $\mathbf{z}(\vx)^\T\mathbf{z}(\vx') = \frac{\sigma_f^2}{M}\sum_{j}\cos\boldsymbol{\omega}_j^\T(\vx - \vx')$, an average whose expectation is $k(\vx - \vx')$. These are the random Fourier features of Rahimi and Recht (2007). With $\sigma_f = 1$ each term lies in $[-1, 1]$ and has variance $\tfrac12(1 + k(2r)) - k(r)^2$, which tends to $\tfrac12$ at large distances (Exercise 10.3), so the error at one distance is typically about $1/\sqrt{2M}$. Accuracy $\varepsilon$ at all pairs of inputs in a bounded domain at once needs $M$ of order $(d/\varepsilon^2)\log(1/\varepsilon)$, with the domain's diameter and the spread of $p$ entering through the logarithm (Rahimi and Recht, 2007, claim 1).

A weighted sum $\mathbf{w}^\T\mathbf{z}(\vx)$ with $\mathbf{w} \sim \N(\mathbf{0}, \mI)$ is the linear model of Equation (7.1) with $2M$ features: an approximate draw from $\GP(0, k)$ that is a formula, costs $O(Md)$ per input, and can be maximized like any function. This is one way for Thompson sampling, which evaluates where a random draw from the posterior is largest (Section 12.5), to draw whole functions on a continuous domain, as Section 8.5 mentioned (Rahimi and Recht, 2007; Wilson et al., 2020).

Read the default. The Matérn 3/2 density falls like a power and is still near $10^{-3}$ of its peak at $\omega\ell = 10$, where the RBF density has vanished. With $M = 20$ the estimated kernel wanders around the true one, with a root-mean-square gap of 0.152 against $1/\sqrt{2M} = 0.158$.

Raise M to 200, then 2000. The gap falls to 0.058, then 0.023, against 0.050 and 0.016: roughly the $1/\sqrt{M}$ rate of an average, a factor of about three per tenfold increase in $M$.

Switch to Matérn 1/2 at M = 20. The heavy Cauchy tail puts a few of the twenty frequencies far out, and they show as a regular ripple, in the kernel estimate and in the feature draw, which is otherwise smooth. The exact draw is rough at every scale. A heavy-tailed spectrum spreads its roughness over many high frequencies, each with little variance, and twenty features sample only a few of them.

Switch to RBF. The frequencies all lie within a few multiples of $1/\ell$, and the feature draw has the character of the exact one at small $M$, because RBF draws are made of such frequencies. The kernel estimate is no more accurate than before: its gap at $M = 20$ is 0.161.

10.5 From eigenvalues to information gain #

The regret bounds of Chapter 13 (bounds on an optimizer's total shortfall from the best value) depend on the maximum information gain $\gamma_T$ (Definition 13.3), the most that $T$ noisy evaluations could reveal about $f$. Section 6.5 quoted its growth for the RBF and Matérn kernels. The rates come from the eigenvalues.

10.5.1 Counting resolved directions #

Write the kernel matrix of $T$ inputs through Equation (10.6) as $\mK_A = \boldsymbol{\Phi}\boldsymbol{\Lambda}\boldsymbol{\Phi}^\T$, where row $t$ of $\boldsymbol{\Phi}$ holds the eigenfunctions at $\vx_t$ and $\boldsymbol{\Lambda}$ holds the eigenvalues on its diagonal. The determinant lemma (Equation (B.7)) moves the information gain Equation (6.15) from the evaluations to the eigen-directions: $\det(\mI + \sigma_n^{-2}\boldsymbol{\Phi}\boldsymbol{\Lambda}\boldsymbol{\Phi}^\T) = \det(\mI + \sigma_n^{-2}\boldsymbol{\Lambda}^{1/2}\boldsymbol{\Phi}^\T\boldsymbol{\Phi}\boldsymbol{\Lambda}^{1/2})$. If the inputs are spread in proportion to $q$, then $\frac1T\sum_t\varphi_i(\vx_t)\varphi_j(\vx_t)$ approximates $\int\varphi_i\varphi_j\, q$, which is 1 or 0, as the Riemann sum of Section 10.1.1 did, so $\boldsymbol{\Phi}^\T\boldsymbol{\Phi} \approx T\mI$ and

Each eigen-direction is a separate experiment: its coefficient has prior variance $\lambda_i$, and $T$ spread-out evaluations measure it with noise variance about $\sigma_n^2/T$, which Equation (6.12) turns into the term above. A direction with $T\lambda_i \gg \sigma_n^2$ is resolved and adds about $\tfrac12\log(T\lambda_i/\sigma_n^2)$; one with $T\lambda_i \ll \sigma_n^2$ adds about $T\lambda_i/(2\sigma_n^2)$. The information is roughly the number of resolved directions times half a logarithm, plus $T/(2\sigma_n^2)$ times the unresolved eigenvalue mass.

Equation (10.13) is an approximation. Against the exact information of $T$ evenly spaced evaluations on $[0, 1]$ (lengthscale 0.1, $\sigma_n = 0.1$), it is within 0.25% at $T = 100$ for the RBF and the Matérn 5/2 and 3/2 kernels. For Matérn 1/2 it is 46% too high at $T = 100$, because it counts 143 resolved directions, more than 100 evaluations can resolve, and 7.4% too high at $T = 1000$.

10.5.2 A bound from the eigenvalues #

The count can be made rigorous for any design by cutting at $m$ directions and paying for the eigenvalue mass beyond them.

The tail $\delta_m$ is eigenvalue mass, not a failure probability like the $\delta$ of Section 10.2.4; the letter follows the paper. The bound on the eigenfunctions is an assumption, which the authors state holds for the kernels used in practice; we found no proof that it holds for the RBF and Matérn kernels on every domain and weighting. The proof needs only the tools of Section 6.3.

10.5.3 Polynomial and exponential decay #

The best $m$ balances the two terms of Equation (10.14), directions kept times $\log T$ against $T$ times the mass beyond them, and the balance depends only on how fast the eigenvalues fall.