Kernel Methods

Most models so far learned a parameter vector and then threw the training data away. Kernel methods keep the data: a prediction is a weighted combination of similarities between the new input and the training inputs. This chapter shows how that form arises from ridge regression (the dual representation), which functions qualify as similarities (valid kernels), how kernels give radial basis function and Nadaraya-Watson regression, and how placing a Gaussian distribution directly on functions gives Gaussian processes for regression and classification, with hyperparameters learned by the marginal likelihood.

TrungTin Nguyen · Following Bishop (2006), Pattern Recognition and Machine Learning, Chapter 6 · Course home · Mathematical toolkit · Python code for the figures

Ch. 1 Introduction Ch. 2 Probability distributions Ch. 3 Linear regression Ch. 4 Linear classification Ch. 5 Neural networks Ch. 6 Kernel methods Ch. 7 Sparse kernel machines Ch. 8 Graphical models Ch. 9 Mixture models and EM Ch. 10 Approximate inference Ch. 11 Sampling methods Ch. 12 Continuous latent variables Ch. 13 Sequential data Ch. 14 Combining models

Memory-based methods and the kernel trick

pp. 291–292 A parametric model (Chapters 3 to 5) compresses the training set into a parameter vector and then ignores it. A memory-based method keeps the training points, or some of them, and uses them again at prediction time: a Parzen density estimate puts a small bump on every training point, and a nearest-neighbour classifier copies the label of the closest one. Such methods need a way to measure similarity between two inputs, and are cheap to "train" but costly to evaluate.

Many linear models can be rewritten so that predictions also have this form. The key object is the kernel

\[ k(\mathbf{x},\mathbf{x}')=\bm\phi(\mathbf{x})^\T\bm\phi(\mathbf{x}'), \tag{6.1} \]

an inner product in a feature space. It is symmetric, \(k(\mathbf{x},\mathbf{x}')\)\({}=k(\mathbf{x}',\mathbf{x})\). The simplest case is \(\bm\phi(\mathbf{x})=\mathbf{x}\), the linear kernel \(k=\mathbf{x}^\T\mathbf{x}'\). The kernel trick (kernel substitution): if an algorithm uses the inputs only through inner products, replace each inner product by a kernel and you obtain a nonlinear version of the algorithm, without ever computing \(\bm\phi\). Examples are kernel PCA (Section 12.3), kernel nearest neighbours (Exercise 6.3) and the kernel Fisher discriminant. Kernels were introduced to pattern recognition by Aizerman, Braverman and Rozonoer (1964) and brought back by Boser, Guyon and Vapnik (1992) in support vector machines (Chapter 7).

A kernel that depends only on the difference of its arguments, \(k(\mathbf{x}-\mathbf{x}')\), is stationary (unchanged by translating both inputs); one that depends only on \(\|\mathbf{x}-\mathbf{x}'\|\) is homogeneous or a radial basis function (Section 6.3).

6.1 Dual representations

pp. 293–294 Consider regularised least squares with a fixed feature map \(\bm\phi\) and \(\lambda\ge0\),

\[ J(\mathbf{w})=\frac12\sum_{n=1}^{N}\{\mathbf{w}^\T\bm\phi(\mathbf{x}_n)-t_n\}^2+\frac{\lambda}{2}\mathbf{w}^\T\mathbf{w}. \tag{6.2} \]

Other models have duals too. The perceptron's weight vector is a sum of the misclassified feature vectors, so it also has a dual form in which \(\bm\phi\) appears only through \(k\) (Exercise 6.2); and in Chapter 7 the dual is the main tool for support vector machines. Exercise 6.16 shows the general reason: whenever the objective depends on \(\mathbf{w}\) only through the values \(\mathbf{w}^\T\bm\phi(\mathbf{x}_n)\) and a monotone function of \(\mathbf{w}^\T\mathbf{w}\), the minimiser is a linear combination of the \(\bm\phi(\mathbf{x}_n)\) (the representer theorem, Kimeldorf and Wahba 1971). In Section 6.4 the same formula (6.9) reappears as the mean of a Gaussian process.

6.2 Constructing kernels

pp. 294–299 There are two ways to obtain a kernel: choose a feature map and compute the inner product, or write down a function \(k\) directly and make sure that it is an inner product of some feature map. The first route, with \(M\) basis functions \(\phi_i\), gives

\[ k(x,x')=\bm\phi(x)^\T\bm\phi(x')=\sum_{i=1}^{M}\phi_i(x)\phi_i(x'), \tag{6.10} \]

shown in Figure 6.1 for polynomial, Gaussian and logistic-sigmoid basis functions on a one-dimensional input space. For the second route, a simple example is

\[ k(\mathbf{x},\mathbf{z})=(\mathbf{x}^\T\mathbf{z})^2. \tag{6.11} \]

For \(\mathbf{x}=(x_1,x_2)\) expand the square:

\[ (\mathbf{x}^\T\mathbf{z})^2=(x_1z_1+x_2z_2)^2=x_1^2z_1^2+2x_1z_1x_2z_2+x_2^2z_2^2=(x_1^2,\sqrt2\,x_1x_2,x_2^2)(z_1^2,\sqrt2\,z_1z_2,z_2^2)^\T=\bm\phi(\mathbf{x})^\T\bm\phi(\mathbf{z}). \tag{6.12} \]

The feature map \(\bm\phi(\mathbf{x})\)\({}=(x_1^2,\sqrt2x_1x_2,x_2^2)^\T\) contains all second-order monomials with a specific weighting. In \(D\) dimensions, \((\mathbf{x}^\T\mathbf{z})^2\)\({}=\sum_{i,j}(x_ix_j)(z_iz_j)\), so the features are all the products \(x_ix_j\). More generally \((\mathbf{x}^\T\mathbf{z})^M\) contains all monomials of order exactly \(M\), and \((\mathbf{x}^\T\mathbf{z}+c)^M\), \(c>0\), all monomials up to order \(M\) (binomial theorem: \(\sum_j\binom{M}{j}c^{M-j}(\mathbf{x}^\T\mathbf{z})^j\), a sum of valid kernels with positive weights, see (6.13), (6.17), (6.18) below). For two images, \((\mathbf{x}^\T\mathbf{z})^M\) is a weighted sum of all products of \(M\) pixels of the first image with \(M\) pixels of the second.

Proof

(valid \(\Rightarrow\) PSD). Symmetry of \(k\) is symmetry of the inner product, and \(\mathbf{a}^\T\mathbf{K}\mathbf{a}\)\({}=\sum_{n,m}a_na_m\langle\bm\phi(\mathbf{x}_n),\bm\phi(\mathbf{x}_m)\rangle\)\({}=\big\|\sum_na_n\bm\phi(\mathbf{x}_n)\big\|^2\)\({}\ge0\).

(PSD \(\Rightarrow\) valid), finite set. For one fixed set of \(N\) points, \(\mathbf{K}=\mathbf{U}\bm\Lambda\mathbf{U}^\T\) with \(\bm\Lambda\ge0\) (spectral theorem), and the vectors \(\bm\phi_n\)\({}=\bm\Lambda^{1/2}\mathbf{U}^\T\mathbf{e}_n\in\R^N\) satisfy \(\bm\phi_n^\T\bm\phi_m\)\({}=(\mathbf{U}\bm\Lambda\mathbf{U}^\T)_{nm}\)\({}=K_{nm}\). (This is also why a Gram matrix with a negative eigenvalue cannot come from any feature map.)

(PSD \(\Rightarrow\) valid), whole space (the reproducing-kernel construction of Aronszajn, 1950). Let \(\mathcal{H}_0\) be the vector space of functions \(f=\sum_{i=1}^{p}\alpha_ik(\cdot,\mathbf{x}_i)\), and define \(\langle f,g\rangle\)\({}=\sum_{i,j}\alpha_i\beta_jk(\mathbf{x}_i,\mathbf{x}'_j)\) for \(g=\sum_j\beta_jk(\cdot,\mathbf{x}'_j)\). It depends only on the functions: \(\langle f,g\rangle=\sum_j\beta_jf(\mathbf{x}'_j)\), and symmetrically \(=\sum_i\alpha_ig(\mathbf{x}_i)\). It is bilinear and symmetric, and \(\langle f,f\rangle\)\({}=\bm\alpha^\T\mathbf{K}\bm\alpha\)\({}\ge0\) by hypothesis. It has the reproducing property \(\langle f,k(\cdot,\mathbf{x})\rangle\)\({}=f(\mathbf{x})\). Cauchy-Schwarz holds for any positive semidefinite bilinear form, so \(f(\mathbf{x})^2\)\({}=\langle f,k(\cdot,\mathbf{x})\rangle^2\)\({}\le\langle f,f\rangle k(\mathbf{x},\mathbf{x})\); hence \(\langle f,f\rangle=0\) forces \(f\equiv0\) and the form is a genuine inner product. Completing \(\mathcal{H}_0\) gives a Hilbert space \(\mathcal{H}\), and \(\bm\phi(\mathbf{x})=k(\cdot,\mathbf{x})\) satisfies \(\langle\bm\phi(\mathbf{x}),\bm\phi(\mathbf{x}')\rangle\)\({}=k(\mathbf{x},\mathbf{x}')\). This \(\mathcal{H}\) is the reproducing kernel Hilbert space of \(k\).∎

Remark. "Positive semidefinite" does not mean "entries \(\ge0\)": the matrix \(\begin{pmatrix}1&2\\2&1\end{pmatrix}\) has positive entries and eigenvalues \(3,-1\), while \(\begin{pmatrix}2&-1\\-1&2\end{pmatrix}\) has a negative entry and eigenvalues \(1,3\) (Exercise 6.4).

What is proved, and why it matters here

The proof (Mercer, 1909; see Rasmussen and Williams, 2006, Section 4.3) uses the spectral theorem for compact self-adjoint operators plus Dini's theorem for the uniform convergence; I do not repeat it. Two connections are worth recording.

Gram matrices. If every Gram matrix is PSD, approximate the double integral by a Riemann sum over a fine partition of \(\mathcal{X}\) (for continuous \(f\) and \(k\) the sums converge): each sum is \(\sum_{n,m}a_na_mk(\mathbf{x}_n,\mathbf{x}_m)\) with \(a_n=f(\mathbf{x}_n)\mu(\text{cell}_n)\), which is \(\ge0\), so the integral is \(\ge0\) for continuous \(f\), hence for all of \(L^2\) by density. Conversely, if the integral condition holds, choose \(f\) as a sharply peaked bump at each \(\mathbf{x}_n\) with weight \(a_n\) and let the bumps shrink; continuity of \(k\) gives \(\mathbf{a}^\T\mathbf{K}\mathbf{a}\ge0\).

Functions drawn from a Gaussian process. If \(\xi_i\sim\N(0,1)\) are independent, then \(y(\mathbf{x})\)\({}=\sum_i\sqrt{\lambda_i}\,\xi_i\,\psi_i(\mathbf{x})\) has \(\E[y(\mathbf{x})y(\mathbf{x}')]\)\({}=\sum_i\lambda_i\psi_i(\mathbf{x})\psi_i(\mathbf{x}')\)\({}=k(\mathbf{x},\mathbf{x}')\): the Gaussian process with kernel \(k\) is a Bayesian linear model with infinitely many basis functions \(\sqrt{\lambda_i}\psi_i\) and a standard normal prior on the weights (Section 6.4.1).

The Gram test is hard to apply directly, so kernels are built from simpler ones.

Proof (each rule, by the Gram criterion of Result 6.1)

Fix points \(\mathbf{x}_1,\dots,\mathbf{x}_N\), write \(\mathbf{K}_1,\mathbf{K}_2\) for the Gram matrices of \(k_1,k_2\) (both PSD), and let \(\mathbf{a}\in\R^N\) be arbitrary. Symmetry of each new function is obvious in every case.

(6.13). The Gram matrix is \(c\mathbf{K}_1\) and \(\mathbf{a}^\T(c\mathbf{K}_1)\mathbf{a}\)\({}=c\,\mathbf{a}^\T\mathbf{K}_1\mathbf{a}\)\({}\ge0\) for \(c\ge0\). (Feature map: \(\sqrt c\,\bm\phi\).)

(6.14). With \(\mathbf{D}\)\({}=\diag(f(\mathbf{x}_1),\dots,f(\mathbf{x}_N))\) the Gram matrix is \(\mathbf{D}\mathbf{K}_1\mathbf{D}\), and \(\mathbf{a}^\T\mathbf{D}\mathbf{K}_1\mathbf{D}\mathbf{a}\)\({}=(\mathbf{D}\mathbf{a})^\T\mathbf{K}_1(\mathbf{D}\mathbf{a})\)\({}\ge0\). (Feature map: \(f(\mathbf{x})\bm\phi(\mathbf{x})\).)

(6.17). \(\mathbf{a}^\T(\mathbf{K}_1+\mathbf{K}_2)\mathbf{a}\)\({}=\mathbf{a}^\T\mathbf{K}_1\mathbf{a}\)\({}+\mathbf{a}^\T\mathbf{K}_2\mathbf{a}\)\({}\ge0\). (Feature map: the concatenation \((\bm\phi_1,\bm\phi_2)\).)

(6.18). The Gram matrix is the entrywise (Hadamard) product \(\mathbf{K}_1\circ\mathbf{K}_2\). Write \(\mathbf{K}_1\)\({}=\sum_i\lambda_i\mathbf{u}_i\mathbf{u}_i^\T\) with \(\lambda_i\ge0\) (spectral theorem). Then \((\mathbf{K}_1\circ\mathbf{K}_2)_{nm}\)\({}=\sum_i\lambda_iu_{in}(\mathbf{K}_2)_{nm}u_{im}\), that is, \(\mathbf{K}_1\circ\mathbf{K}_2\)\({}=\sum_i\lambda_i\,\diag(\mathbf{u}_i)\mathbf{K}_2\diag(\mathbf{u}_i)\), a non-negative combination of matrices that are PSD by the argument for (6.14). This is Schur's product theorem. With explicit features \(k_1=\sum_m\phi_m\phi_m'\) and \(k_2=\sum_l\psi_l\psi_l'\) one can also write \(k_1k_2=\sum_{m,l}(\phi_m\psi_l)(\phi_m'\psi_l')\): the new features are all products \(\phi_m(\mathbf{x})\psi_l(\mathbf{x})\), of dimension \(M_1M_2\) (the tensor product).

(6.15). Let \(q(z)=\sum_{j=0}^{J}c_jz^j\) with all \(c_j\ge0\). The constant function \(1\) is a valid kernel (feature \(\phi\equiv1\)). By (6.18) and induction \(k_1^j\) is valid for each \(j\ge1\) (and \(k_1^0=1\)); by (6.13) \(c_jk_1^j\) is valid; by (6.17) the sum is valid.

(6.16). The partial sums \(s_J=\sum_{j=0}^{J}k_1^j/j!\) are of the form (6.15), hence valid, and \(s_J(\mathbf{x},\mathbf{x}')\)\({}\to\exp(k_1(\mathbf{x},\mathbf{x}'))\) for every pair. The Gram matrix of \(\exp(k_1)\) is therefore the entrywise limit of PSD matrices \(\mathbf{S}_J\), and \(\mathbf{a}^\T\mathbf{S}_J\mathbf{a}\ge0\) for all \(J\) gives, in the limit, \(\mathbf{a}^\T\mathbf{S}\mathbf{a}\ge0\) (the set of PSD matrices is closed).

(6.19). The Gram matrix of \(k\) at \(\mathbf{x}_1,\dots,\mathbf{x}_N\) is the Gram matrix of \(k_3\) at the points \(\bm\phi(\mathbf{x}_1),\dots,\bm\phi(\mathbf{x}_N)\in\R^M\), which is PSD (repeated points are allowed). (Feature map: \(\bm\psi\circ\bm\phi\), with \(\bm\psi\) a feature map of \(k_3\).)

(6.20). Write \(\mathbf{A}=\mathbf{L}^\T\mathbf{L}\) (for instance \(\mathbf{L}=\mathbf{A}^{1/2}\), from the spectral decomposition, which exists since the eigenvalues are \(\ge0\)). Then \(\mathbf{x}^\T\mathbf{A}\mathbf{x}'\)\({}=(\mathbf{L}\mathbf{x})^\T(\mathbf{L}\mathbf{x}')\), an inner product of the features \(\mathbf{L}\mathbf{x}\). (This is (6.19) with the linear kernel.)

(6.21), (6.22). The function \(\tilde k_a(\mathbf{x},\mathbf{x}')\)\({}=k_a(\mathbf{x}_a,\mathbf{x}_a')\) is (6.19) with \(\bm\phi(\mathbf{x})=\mathbf{x}_a\), hence valid on the whole space, and likewise \(\tilde k_b\). The sum (6.21) is (6.17) applied to \(\tilde k_a,\tilde k_b\), and the product (6.22) is (6.18). The parts \(\mathbf{x}_a,\mathbf{x}_b\) may share components.∎

Remark: what is not allowed. The difference \(k_1-k_2\) need not be valid (take \(k_2=2k_1\), then \(-k_1\) has a negative diagonal), and a general function \(g(k_1)\) of a kernel is valid only under conditions such as those of (6.15)–(6.16). Constants \(c<0\) are excluded in (6.13).

With these rules the polynomial kernels above are valid: \((\mathbf{x}^\T\mathbf{x}')^M\) by (6.18) repeatedly, \((\mathbf{x}^\T\mathbf{x}'+c)^M\) by (6.15) with \(q(z)=(z+c)^M\).

The Gaussian kernel

pp. 296–297; Exercise 6.11 The most used kernel is the "Gaussian"

\[ k(\mathbf{x},\mathbf{x}')=\exp\!\big(-\|\mathbf{x}-\mathbf{x}'\|^2/2\sigma^2\big). \tag{6.23} \]

It is not a probability density here, so the normalising constant is omitted.

Kernels for non-vector inputs

p. 297; Exercise 6.12 Kernels can be defined on graphs, sets, strings or documents. For subsets \(A_1,A_2\) of a fixed finite set \(D\),

\[ k(A_1,A_2)=2^{|A_1\cap A_2|}, \tag{6.27} \]

where \(|A|\) is the number of elements of \(A\) (the first printing says "subsets"; the errata corrects it). It is valid: take the feature space \(\R^{2^{|D|}}\) indexed by the subsets \(U\subseteq D\), with \(\phi_U(A)=1\) if \(U\subseteq A\) and \(0\) otherwise (6.95). Then \(\bm\phi(A_1)^\T\bm\phi(A_2)\) counts the subsets \(U\) with \(U\subseteq A_1\) and \(U\subseteq A_2\), that is \(U\subseteq A_1\cap A_2\), and a set with \(|A_1\cap A_2|\) elements has \(2^{|A_1\cap A_2|}\) subsets.

Kernels from probabilistic generative models

pp. 297–298, (6.28)–(6.31) Generative models handle missing data and variable-length sequences naturally; discriminative models usually classify better. One way to combine them is to let a generative model define the kernel and then use it in a discriminative method (Haussler, 1999). For a density \(p\),

\[ k(\mathbf{x},\mathbf{x}')=p(\mathbf{x})\,p(\mathbf{x}'), \tag{6.28} \]

is valid, an inner product in a one-dimensional feature space with \(\phi(\mathbf{x})=p(\mathbf{x})\): two inputs are similar if both are probable. Using (6.13) and (6.17), with weights \(p(i)\ge0\) summing to one,

\[ k(\mathbf{x},\mathbf{x}')=\sum_{i}p(\mathbf{x}\given i)\,p(\mathbf{x}'\given i)\,p(i), \tag{6.29} \]

whose feature map is \(\phi_i(\mathbf{x})\)\({}=\sqrt{p(i)}\,p(\mathbf{x}\given i)\). It is the probability that two independent draws \(\mathbf{x},\mathbf{x}'\) from the same component \(i\) (chosen with probability \(p(i)\)) take these values, so it is, up to a constant, a mixture whose components factorise over the pair, with \(i\) a latent variable. In the continuous limit,

\[ k(\mathbf{x},\mathbf{x}')=\int p(\mathbf{x}\given\mathbf{z})\,p(\mathbf{x}'\given\mathbf{z})\,p(\mathbf{z})\,\d\mathbf{z}, \tag{6.30} \]

the inner product in \(L^2(p(\mathbf{z})\d\mathbf{z})\) of the functions \(\mathbf{z}\mapsto p(\mathbf{x}\given\mathbf{z})\), hence valid. For sequences \(\mathbf{X}=(\mathbf{x}_1,\dots,\mathbf{x}_L)\) with a hidden Markov model (Section 13.2),

\[ k(\mathbf{X},\mathbf{X}')=\sum_{\mathbf{Z}}p(\mathbf{X}\given\mathbf{Z})\,p(\mathbf{X}'\given\mathbf{Z})\,p(\mathbf{Z}), \tag{6.31} \]

is (6.29) with the hidden state sequence \(\mathbf{Z}\) as the index: both observed sequences are generated by the same hidden path. It can be computed with a forward recursion on the product of the two emission models (a pair HMM), and extended to sequences of different length.

The Fisher kernel

pp. 298–299, (6.32)–(6.36); Exercises 6.13, 6.14 Given a parametric model \(p(\mathbf{x}\given\bm\theta)\), the Fisher score \(\mathbf{g}(\bm\theta,\mathbf{x})\)\({}=\nabla_{\bm\theta}\ln p(\mathbf{x}\given\bm\theta)\) (6.32) maps \(\mathbf{x}\) to a vector of the dimension of \(\bm\theta\): the direction in which the parameters should move to explain \(\mathbf{x}\) better. The Fisher kernel (Jaakkola and Haussler, 1999) is

\[ k(\mathbf{x},\mathbf{x}')=\mathbf{g}(\bm\theta,\mathbf{x})^\T\mathbf{F}^{-1}\mathbf{g}(\bm\theta,\mathbf{x}'),\qquad \mathbf{F}=\E_{\mathbf{x}}\big[\mathbf{g}(\bm\theta,\mathbf{x})\mathbf{g}(\bm\theta,\mathbf{x})^\T\big], \tag{6.33, 6.34} \]

with the expectation under \(p(\mathbf{x}\given\bm\theta)\).

Fisher kernels have been used for document retrieval (Hofmann, 2000) and, with HMMs, for biological sequences and speech.

The sigmoidal kernel is not valid

p. 299, (6.37) The sigmoidal kernel

\[ k(\mathbf{x},\mathbf{x}')=\tanh\big(a\,\mathbf{x}^\T\mathbf{x}'+b\big) \tag{6.37} \]

gives support-vector expansions a superficial resemblance to neural networks (Vapnik, 1995) but, as the book says, its Gram matrix is in general not positive semidefinite. Two explicit counterexamples:

6.3 Radial basis function networks

pp. 299–301 A radial basis function (RBF) depends only on the distance from a centre \(\bm\mu_j\): \(\phi_j(\mathbf{x})=h(\|\mathbf{x}-\bm\mu_j\|)\). Linear models in such basis functions (Chapter 3) are called RBF networks. Three different motivations lead to them.

Exact interpolation (Powell, 1987). Given noise-free data \(\{\mathbf{x}_n,t_n\}\), look for a function through every point, with one basis function centred on each input,

\[ f(\mathbf{x})=\sum_{n=1}^{N}w_n\,h(\|\mathbf{x}-\mathbf{x}_n\|). \tag{6.38} \]

Requiring \(f(\mathbf{x}_m)=t_m\) gives \(N\) linear equations \(\mathbf{H}\mathbf{w}=\mathbf{t}\), \(H_{mn}=h(\|\mathbf{x}_m-\mathbf{x}_n\|)\), for \(N\) unknowns, solvable when \(\mathbf{H}\) is invertible (for the Gaussian \(h\), \(\mathbf{H}\) is a Gram matrix of the Gaussian kernel, positive definite for distinct points). With noisy targets this over-fits: the fit passes through the noise (the \(\lambda\to0\) limit of (6.9) with \(\mathbf{K}=\mathbf{H}\)).

Regularisation theory (Poggio and Girosi, 1990). If the sum of squares is penalised by a differential operator, the optimal function is an expansion in the Green's functions of the operator, one per data point; for an isotropic operator they are radial, and the penalty means the fit no longer interpolates.

Noisy inputs (Webb, 1994; Bishop, 1995). If the input is corrupted, \(\mathbf{x}_n\to\mathbf{x}_n+\bm\xi\) with \(\bm\xi\sim\nu\), the expected sum of squares is

\[ E=\frac12\sum_{n=1}^{N}\int\{y(\mathbf{x}_n+\bm\xi)-t_n\}^2\,\nu(\bm\xi)\,\d\bm\xi. \tag{6.39} \]

The effect of normalisation is shown in Figure 6.2. Without it, the model output (a linear combination) decays to zero (or to the bias parameter) far from every centre; with it, the nearest centres dominate and the prediction approaches the target of the closest point. This is the Nadaraya-Watson model of the next subsection, which appears here as the solution for noisy inputs and there as the solution of kernel density estimation.

Fewer basis functions than data points. Because (6.38) and (6.40) have one basis function per data point, evaluating them is expensive. Broomhead and Lowe (1988), Moody and Darken (1989) and Poggio and Girosi (1990) keep the radial form but use \(Morthogonal least squares (Chen et al., 1991), which adds at each step the data point that most reduces the sum of squares; or cluster centres from K-means (Section 9.1), which no longer coincide with data points. The weights \(w_i\) are then found by ordinary least squares (Section 3.1.1).

6.3.1 Nadaraya-Watson model

pp. 301–303, (6.42)–(6.48) In Section 3.3.3 the prediction of linear regression was a linear combination of the training targets with the equivalent kernel (3.62). A different route to such a form starts from density estimation. Model the joint density of input and target by a Parzen estimator (Section 2.5.1) with one component \(f\) on every data point,

\[ p(\mathbf{x},t)=\frac1N\sum_{n=1}^{N}f(\mathbf{x}-\mathbf{x}_n,\,t-t_n), \tag{6.42} \]

and take as the regression function the conditional mean \(y(\mathbf{x})\)\({}=\E[t\given\mathbf{x}]\)\({}=\int tp(t\given\mathbf{x})\,\d t\)\({}=\int tp(\mathbf{x},t)\,\d t\big/\int p(\mathbf{x},t)\,\d t\).

The result is the Nadaraya-Watson model, or kernel regression (Nadaraya, 1964; Watson, 1964). For a localised \(g\) it gives most weight to the points near \(\mathbf{x}\). It is also the minimiser over a constant \(y\) of the locally weighted squares \(\sum_nk(\mathbf{x},\mathbf{x}_n)(y-t_n)^2\) (compare Gap 6.6): a local average. The width \(\sigma\) controls bias against variance: small \(\sigma\) follows the noise (nearly nearest-neighbour), large \(\sigma\) tends to the global mean.

More flexible components are possible, for instance different variances for \(\mathbf{x}\) and \(t\), or a Gaussian mixture fitted to \(p(t,\mathbf{x})\) by EM (Chapter 9; Ghahramani and Jordan, 1994) and then conditioned. The latter has fewer components than data points, so predictions are cheaper, but it is no longer a sum of kernels evaluated at the training points, and the cost has moved to training.

6.4 Gaussian processes

pp. 303–304 In Section 6.1 kernels arose from a non-probabilistic model. Here they arise in a Bayesian setting. The idea of a Gaussian process (GP) is to drop the parametric model \(y(\mathbf{x},\mathbf{w})\) and put a prior directly on the function \(y(\mathbf{x})\). This looks hard, since functions form an infinite-dimensional space, but only the function values at the training and test inputs ever matter, and those are a finite vector. A Gaussian process is a collection of random variables \(\{y(\mathbf{x}):\mathbf{x}\in\mathcal{X}\}\) such that every finite subcollection has a joint Gaussian distribution. Because a marginal of a Gaussian is Gaussian, the specifications for different finite sets are consistent with each other (the marginal of \(\N(\mathbf{0},\mathbf{K}_{N+1})\) over the first \(N\) variables is \(\N(\mathbf{0},\mathbf{K}_N)\), by the marginal formula in the toolkit), and the Kolmogorov extension theorem then gives a process on all of \(\mathcal{X}\).

The equivalents of GP regression have many names: kriging in geostatistics (Cressie, 1993); ARMA models, Kalman filters (Section 13.3) and RBF networks can all be viewed as Gaussian processes. Reviews: MacKay (1998), Williams (1999), Rasmussen and Williams (2006).

6.4.1 Linear regression revisited

pp. 304–306, (6.49)–(6.56) Start from a model one already knows: \(y(\mathbf{x})=\mathbf{w}^\T\bm\phi(\mathbf{x})\) with \(M\) fixed basis functions (6.49), and the prior

\[ p(\mathbf{w})=\N(\mathbf{w}\given\mathbf{0},\alpha^{-1}\mathbf{I}) \tag{6.50} \]

with precision \(\alpha\). Each draw of \(\mathbf{w}\) is a function. At the training inputs the function values form the vector

\[ \mathbf{y}=\bm\Phi\mathbf{w}, \tag{6.51} \]

with \(\bm\Phi_{nk}=\phi_k(\mathbf{x}_n)\). Since \(\mathbf{y}\) is a linear transformation of a Gaussian vector it is Gaussian (a linear map of a Gaussian: the characteristic function \(\E[e^{\mathrm{i}\mathbf{s}^\T\bm\Phi\mathbf{w}}]\)\({}=e^{-\frac{1}{2\alpha}\|\bm\Phi^\T\mathbf{s}\|^2}\) is again that of a Gaussian, see toolkit), so only its mean and covariance are needed.

In general, a stochastic process \(y(\mathbf{x})\) is specified by the joint distribution of any finite set of its values, in a consistent manner; for a Gaussian process the second-order statistics determine everything. With no prior knowledge of the mean we take it to be zero (by symmetry, as for the zero-mean weight prior), and the process is then specified by its covariance function, the kernel,

\[ \E[y(\mathbf{x}_n)y(\mathbf{x}_m)]=k(\mathbf{x}_n,\mathbf{x}_m). \tag{6.55} \]

(For two-dimensional inputs the process is also called a Gaussian random field.) The kernel can also be chosen directly. Figure 6.4 shows functions sampled from a GP with the Gaussian kernel (6.23) and with the exponential kernel

\[ k(x,x')=\exp(-\theta|x-x'|), \tag{6.56} \]

the covariance of the Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein, 1930).

6.4.2 Gaussian processes for regression

pp. 306–311, (6.57)–(6.68) Observed targets contain noise, \(t_n=y_n+\epsilon_n\) (6.57), with \(\epsilon_n\) independent for each \(n\) and Gaussian,

\[ p(t_n\given y_n)=\N(t_n\given y_n,\beta^{-1}), \tag{6.58} \] \[ p(\mathbf{t}\given\mathbf{y})=\N(\mathbf{t}\given\mathbf{y},\beta^{-1}\mathbf{I}_N), \tag{6.59} \]

and the GP prior supplies the marginal

\[ p(\mathbf{y})=\N(\mathbf{y}\given\mathbf{0},\mathbf{K}). \tag{6.60} \]

The kernel is chosen so that nearby inputs have strongly correlated function values. Integrating out \(\mathbf{y}\) gives the distribution of the targets,

\[ p(\mathbf{t})=\int p(\mathbf{t}\given\mathbf{y})p(\mathbf{y})\,\d\mathbf{y}=\N(\mathbf{t}\given\mathbf{0},\mathbf{C}),\qquad C(\mathbf{x}_n,\mathbf{x}_m)=k(\mathbf{x}_n,\mathbf{x}_m)+\beta^{-1}\delta_{nm}. \tag{6.61, 6.62} \]

This is the linear-Gaussian result (2.115) with the identity as the linear map (toolkit); it can also be seen directly, since \(\mathbf{t}=\mathbf{y}+\bm\epsilon\) is a sum of two independent Gaussian vectors with means \(\mathbf{0}\) and covariances \(\mathbf{K}\) and \(\beta^{-1}\mathbf{I}\), and covariances of independent summands add. A popular choice of kernel is the exponential of a quadratic form with constant and linear terms,

\[ k(\mathbf{x}_n,\mathbf{x}_m)=\theta_0\exp\!\Big\{-\frac{\theta_1}2\|\mathbf{x}_n-\mathbf{x}_m\|^2\Big\}+\theta_2+\theta_3\,\mathbf{x}_n^\T\mathbf{x}_m. \tag{6.63} \]

It is valid: the first term is \(\theta_0\) times a Gaussian kernel, the second a constant kernel, the third a linear kernel, combined by (6.13) and (6.17). Here \(\theta_0\) is the signal variance, \(\theta_1=\ell^{-2}\) is the inverse squared length scale, \(\theta_2\) allows a constant offset (a random bias of variance \(\theta_2\)), and \(\theta_3\) a random linear trend, the parametric linear model. Figure 6.5 shows samples for several values; Figure 6.6 shows data points drawn from (6.60)–(6.61).

The predictive distribution. Given the training targets \(\mathbf{t}_N=(t_1,\dots,t_N)^\T\) at \(\mathbf{x}_1,\dots,\mathbf{x}_N\), we want \(p(t_{N+1}\given\mathbf{t}_N)\) at a new input \(\mathbf{x}_{N+1}\) (inputs are suppressed in the notation). The joint distribution of the \(N+1\) targets is, by (6.61) applied to \(N+1\) points,

\[ p(\mathbf{t}_{N+1})=\N(\mathbf{t}_{N+1}\given\mathbf{0},\mathbf{C}_{N+1}),\qquad \mathbf{C}_{N+1}=\begin{pmatrix}\mathbf{C}_N&\mathbf{k}\\\mathbf{k}^\T&c\end{pmatrix}, \tag{6.64, 6.65} \]

where \(\mathbf{C}_N\) has entries (6.62), \(\mathbf{k}\) has entries \(k(\mathbf{x}_n,\mathbf{x}_{N+1})\), and \(c\)\({}=k(\mathbf{x}_{N+1},\mathbf{x}_{N+1})\)\({}+\beta^{-1}\).

Proof (two ways)

(a) From the conditional Gaussian formulas. For a zero-mean Gaussian vector \((\mathbf{x}_a,\mathbf{x}_b)\) with covariance blocks \(\bm\Sigma_{aa},\bm\Sigma_{ab},\bm\Sigma_{bb}\), the conditional is \(\N(\mathbf{x}_a\given\bm\Sigma_{ab}\bm\Sigma_{bb}^{-1}\mathbf{x}_b,\;\bm\Sigma_{aa}-\bm\Sigma_{ab}\bm\Sigma_{bb}^{-1}\bm\Sigma_{ba})\) (eqs (2.81), (2.82), proved in the toolkit). Take \(\mathbf{x}_a=t_{N+1}\), \(\mathbf{x}_b=\mathbf{t}\): \(\bm\Sigma_{aa}=c\), \(\bm\Sigma_{ab}=\mathbf{k}^\T\), \(\bm\Sigma_{bb}=\mathbf{C}_N\). This gives (6.66) and (6.67) at once. (The official solution orders the blocks the other way, with \(t_{N+1}\) first; the result is the same.)

(b) By completing the square. Let \(s\)\({}=(c-\mathbf{k}^\T\mathbf{C}_N^{-1}\mathbf{k})^{-1}\), the inverse Schur complement of \(\mathbf{C}_N\) in \(\mathbf{C}_{N+1}\). The block inverse is

\[ \mathbf{C}_{N+1}^{-1}=\begin{pmatrix}\mathbf{C}_N^{-1}+s\,\mathbf{C}_N^{-1}\mathbf{k}\mathbf{k}^\T\mathbf{C}_N^{-1}&-s\,\mathbf{C}_N^{-1}\mathbf{k}\\-s\,\mathbf{k}^\T\mathbf{C}_N^{-1}&s\end{pmatrix}. \]

As a function of \(\tau=t_{N+1}\) (with \(\mathbf{t}\) fixed) the exponent of the joint density is \(-\frac12(\mathbf{t}^\T,\tau)\mathbf{C}_{N+1}^{-1}(\mathbf{t}^\T,\tau)^\T\)\({}=-\frac12\{s\tau^2\)\({}-2s\tau\,\mathbf{k}^\T\mathbf{C}_N^{-1}\mathbf{t}\}\)\({}+\text{terms free of }\tau\)\({}=-\frac s2(\tau-\mathbf{k}^\T\mathbf{C}_N^{-1}\mathbf{t})^2\)\({}+\text{const}\). Since \(p(\tau\given\mathbf{t})\)\({}\propto p(\mathbf{t},\tau)\), it is Gaussian with precision \(s\) (variance \(1/s\), (6.67)) and mean \(\mathbf{k}^\T\mathbf{C}_N^{-1}\mathbf{t}\) (6.66).

The variance bound. \(\mathbf{C}_{N+1}\) is positive definite, hence so is its Schur complement \(c-\mathbf{k}^\T\mathbf{C}_N^{-1}\mathbf{k}>0\). Since \(\mathbf{k}^\T\mathbf{C}_N^{-1}\mathbf{k}\) is the variance of the part of the new function value explained by the data (a squared norm, \(\ge0\)), \(\sigma^2\le c\) (never worse than the prior), and \(\sigma^2-\beta^{-1}\) is a conditional variance of \(y_{N+1}\) (a Gaussian conditional covariance, \(\ge0\)).∎

Validity of the covariance. The only requirement on the kernel is that \(\mathbf{C}\) be positive definite. If \(\lambda_i\) are the eigenvalues of \(\mathbf{K}\) then those of \(\mathbf{C}=\mathbf{K}+\beta^{-1}\mathbf{I}\) are \(\lambda_i+\beta^{-1}\) (same eigenvectors), so a PSD Gram matrix (Section 6.2) suffices: a zero eigenvalue of \(\mathbf{K}\) still gives a positive eigenvalue of \(\mathbf{C}\) because \(\beta>0\). All the constructions of Section 6.2 can be used.

The mean is a kernel expansion. With \(a_n\) the \(n\)-th component of \(\mathbf{C}_N^{-1}\mathbf{t}\),

\[ m(\mathbf{x}_{N+1})=\sum_{n=1}^{N}a_n\,k(\mathbf{x}_n,\mathbf{x}_{N+1}), \tag{6.68} \]

so for a kernel depending on \(\|\mathbf{x}_n-\mathbf{x}_m\|\) only, this is an expansion in radial basis functions, one per data point. Compare (6.66) with (6.9): they are the same formula with \(\lambda=\beta^{-1}\). The posterior mean of a GP is kernel ridge regression with regularisation parameter equal to the noise variance; the GP adds the predictive variance (6.67). In Figure 6.8 the shaded band is \(\pm2\sigma\).

Kernels for inputs on a circle or a sphere

Sometimes the inputs live on a sphere \(S^{D-1}=\{\mathbf{x}:\|\mathbf{x}\|=1\}\) (directions, geographic positions) or a circle (angles, time of day). The rules of Section 6.2 give valid kernels there, and a tempting construction fails.

6.4.3 Learning the hyperparameters

pp. 311–312, (6.69), (6.70) Predictions depend on the covariance function, so rather than fixing it we use a parametric family and infer its parameters \(\bm\theta\) (length scales, signal and noise variances, ...), the hyperparameters of the GP. The simplest approach maximises the likelihood \(p(\mathbf{t}\given\bm\theta)\), the analogue of the type-2 maximum likelihood (evidence approximation) of Section 3.5. Because \(\mathbf{t}\sim\N(\mathbf{0},\mathbf{C}_N)\) by (6.61), the log likelihood is just the log of a Gaussian density:

\[ \ln p(\mathbf{t}\given\bm\theta)=-\frac12\ln|\mathbf{C}_N|-\frac12\mathbf{t}^\T\mathbf{C}_N^{-1}\mathbf{t}-\frac N2\ln(2\pi). \tag{6.69} \]

The first term penalises complexity (a large determinant means a broad, flexible prior over \(\mathbf{t}\)), the second measures fit; the balance is automatic (Occam's razor), with no cross-validation. With the finite basis kernel (6.54), (6.69) equals the log evidence (3.86) (matrix determinant lemma \(|\beta^{-1}\mathbf{I}_N\)\({}+\alpha^{-1}\bm\Phi\bm\Phi^\T|\)\({}=\beta^{-N}\alpha^{-M}|\alpha\mathbf{I}_M\)\({}+\beta\bm\Phi^\T\bm\Phi|\) and Woodbury for the quadratic form). To maximise it by gradient methods (conjugate gradients, Fletcher 1987; Nocedal and Wright, 1999) we need its gradient.

The log likelihood is generally not concave, so there can be several local maxima (for example "short length scale, small noise" versus "long length scale, large noise"); restarts help. It is easy to add a prior over \(\bm\theta\) and maximise the log posterior (then add \(\partial\ln p(\bm\theta)/\partial\theta_i\) to (6.70)). A fully Bayesian treatment would integrate over \(\bm\theta\), which is intractable and needs approximation (MCMC, Chapter 11). The model assumed that the noise variance \(\beta^{-1}\) is constant; for heteroscedastic problems a second GP can represent \(\ln\beta(\mathbf{x})\), the logarithm because a precision is positive (Goldberg et al., 1998).

6.4.4 Automatic relevance determination

pp. 312–313, (6.71), (6.72) Give each input dimension its own precision parameter \(\eta_i\). For \(\mathbf{x}=(x_1,x_2)\),

\[ k(\mathbf{x},\mathbf{x}')=\theta_0\exp\Big\{-\frac12\sum_{i=1}^{2}\eta_i(x_i-x_i')^2\Big\}. \tag{6.71} \]

This is a product over dimensions of Gaussian kernels with length scales \(\ell_i=\eta_i^{-1/2}\), valid by (6.22) and (6.13). When \(\eta_i\) is small the function varies slowly in the direction \(x_i\); in the limit \(\eta_i\to0\) the kernel (hence every sample function) does not depend on \(x_i\) at all (Figure 6.9). Maximising the marginal likelihood over the \(\eta_i\) (the gradient (6.70) with \(\partial C_{nm}/\partial\eta_i\)\({}=-\frac12(x_{ni}-x_{mi})^2\,k(\mathbf{x}_n,\mathbf{x}_m)\)) therefore detects irrelevant inputs: those whose \(\eta_i\) shrinks towards zero can be discarded. This is automatic relevance determination (ARD), originally formulated for neural networks (MacKay, 1994; Neal, 1996); the same mechanism in the weight-space view is in Section 7.2.2.

ARD is illustrated (Nabney, 2002) on a synthetic data set with three inputs (Figure 6.10): \(t\) is generated from a function of \(x_1\) (\(\sin(2\pi x_1)\) plus Gaussian noise), \(x_2\) is a noisy copy of \(x_1\), and \(x_3\) is independent noise. Maximising the marginal likelihood with scaled conjugate gradients, \(\eta_1\) converges to a relatively large value, \(\eta_2\) to a much smaller one and \(\eta_3\) to a very small value, indicating that \(x_3\) is irrelevant. The ARD idea is easily built into (6.63):

\[ k(\mathbf{x}_n,\mathbf{x}_m)=\theta_0\exp\Big\{-\frac12\sum_{i=1}^{D}\eta_i(x_{ni}-x_{mi})^2\Big\}+\theta_2+\theta_3\sum_{i=1}^{D}x_{ni}x_{mi}, \tag{6.72} \]

which has been found useful for a range of regression problems.

6.4.5 Gaussian processes for classification

pp. 313–315, (6.73)–(6.76) A GP produces values on the whole real line, but class probabilities must lie in \((0,1)\). As in logistic regression, put the GP on an activation function \(a(\mathbf{x})\) and squash it with the logistic sigmoid \(y=\sigma(a)\) (4.59). For two classes, \(t\in\{0,1\}\), the Bernoulli likelihood is

\[ p(t\given a)=\sigma(a)^t\big(1-\sigma(a)\big)^{1-t}. \tag{6.73} \]

A sample from a GP over \(a(\mathbf{x})\), transformed pointwise, is a non-Gaussian random function with values in \((0,1)\) (Figure 6.11).

Let the training inputs be \(\mathbf{x}_1,\dots,\mathbf{x}_N\) with labels \(\mathbf{t}_N=(t_1,\dots,t_N)^\T\) and a test point \(\mathbf{x}_{N+1}\). Put the GP prior on the vector \(\mathbf{a}_{N+1}\)\({}=(a(\mathbf{x}_1),\dots,a(\mathbf{x}_{N+1}))^\T\),

\[ p(\mathbf{a}_{N+1})=\N(\mathbf{a}_{N+1}\given\mathbf{0},\mathbf{C}_{N+1}),\qquad C(\mathbf{x}_n,\mathbf{x}_m)=k(\mathbf{x}_n,\mathbf{x}_m)+\nu\,\delta_{nm}. \tag{6.74, 6.75} \]

There is no noise term of the regression type, because the labels are assumed correctly assigned (the Bernoulli likelihood already contains the randomness). The small constant \(\nu\), typically fixed in advance (for instance \(10^{-6}\) up to \(10^{-2}\)), only keeps \(\mathbf{C}\) positive definite numerically; a larger \(\nu\) adds a small independent component to the latent function. The kernel \(k\) is any PSD kernel with parameters \(\bm\theta\). For two classes it suffices to predict

\[ p(t_{N+1}=1\given\mathbf{t}_N)=\int p(t_{N+1}=1\given a_{N+1})\,p(a_{N+1}\given\mathbf{t}_N)\,\d a_{N+1},\qquad p(t_{N+1}=1\given a_{N+1})=\sigma(a_{N+1}). \tag{6.76} \]

This integral is intractable. It can be approximated by sampling (Neal, 1997), or analytically if one has a Gaussian approximation to the posterior \(p(a_{N+1}\given\mathbf{t}_N)\), because the convolution of a sigmoid with a Gaussian has the approximation (4.153). The usual justification for a Gaussian posterior (central limit theorem) does not apply directly, since the number of variables grows with the data; but if the number of points in a fixed region increases, the uncertainty in \(a(\mathbf{x})\) there shrinks and the posterior tends to a Gaussian (Williams and Barber, 1998). Three approaches have been used: the local variational bound (10.144) of Section 10.6 (Gibbs and MacKay, 2000; extended to \(K>2\) by Gibbs, 1997), expectation propagation (Section 10.7; Opper and Winther, 2000; Minka, 2001; Seeger, 2003), which works well because the true posterior is unimodal, and the Laplace approximation, treated next.

6.4.6 Laplace approximation

pp. 315–319, (6.77)–(6.94) The Laplace approximation (Section 4.4) replaces a posterior by a Gaussian centred at its mode with covariance the inverse of the negative Hessian of the log density there. We need it for \(p(a_{N+1}\given\mathbf{t}_N)\). The route is: marginalise \(\mathbf{a}_N\) in two steps, first approximating \(p(\mathbf{a}_N\given\mathbf{t}_N)\), then using the GP conditional for the test point.

The prior is \(p(\mathbf{a}_N)\)\({}=\N(\mathbf{a}_N\given\mathbf{0},\mathbf{C}_N)\), and for independent labels the likelihood is, using \(\sigma(a)=e^a\sigma(-a)\) and \(1-\sigma(a)=\sigma(-a)\),

\[ p(\mathbf{t}_N\given\mathbf{a}_N)=\prod_{n=1}^{N}\sigma(a_n)^{t_n}\big(1-\sigma(a_n)\big)^{1-t_n}=\prod_{n=1}^{N}e^{a_nt_n}\sigma(-a_n). \tag{6.79} \]

Taylor-expand the log posterior up to a constant,

\[ \Psi(\mathbf{a}_N)=\ln p(\mathbf{a}_N)+\ln p(\mathbf{t}_N\given\mathbf{a}_N)=-\frac12\mathbf{a}_N^\T\mathbf{C}_N^{-1}\mathbf{a}_N-\frac N2\ln(2\pi)-\frac12\ln|\mathbf{C}_N|+\mathbf{t}_N^\T\mathbf{a}_N-\sum_{n=1}^{N}\ln(1+e^{a_n}). \tag{6.80} \]

Erratum. The book prints "+const" at the end of (6.80); it should be omitted, since the right side is already the exact sum \(\ln p(\mathbf{a}_N)\)\({}+\ln p(\mathbf{t}_N\given\mathbf{a}_N)\) (the "up to normalisation" refers to \(p(\mathbf{t}_N)\)). To see the last term: \(\ln\sigma(-a)=-\ln(1+e^{a})\).

At the mode we evaluate the Hessian, with \(\mathbf{W}_N\) computed from \(\mathbf{a}_N^\star\),

\[ \mathbf{H}=-\nabla\nabla\Psi(\mathbf{a}_N^\star)=\mathbf{W}_N+\mathbf{C}_N^{-1}, \tag{6.85} \]

which defines the Gaussian approximation to the posterior,

\[ q(\mathbf{a}_N)=\N\big(\mathbf{a}_N\given\mathbf{a}_N^\star,\mathbf{H}^{-1}\big). \tag{6.86} \]

By Woodbury, \(\mathbf{H}^{-1}\)\({}=\mathbf{C}\)\({}-\mathbf{C}(\mathbf{W}^{-1}+\mathbf{C})^{-1}\mathbf{C}\)\({}=\mathbf{C}(\mathbf{I}+\mathbf{W}\mathbf{C})^{-1}\): the posterior covariance is smaller than the prior covariance \(\mathbf{C}\) in the positive semidefinite order, as it should be.

Learning the covariance parameters

pp. 317–318, (6.89)–(6.94) The kernel parameters \(\bm\theta\) are found by maximising the likelihood, with a penalty if desired,

\[ p(\mathbf{t}_N\given\bm\theta)=\int p(\mathbf{t}_N\given\mathbf{a}_N)\,p(\mathbf{a}_N\given\bm\theta)\,\d\mathbf{a}_N. \tag{6.89} \]

Extension to \(K>2\) classes uses the softmax in place of the sigmoid (Williams and Barber, 1998). Figure 6.12 uses a synthetic two-class problem (Appendix A of the book).

6.4.7 Connection to neural networks

pp. 319–320 For a Bayesian neural network, a prior over the weights \(\mathbf{w}\) together with the network function induces a prior over functions. In the framework of maximum likelihood the number of hidden units \(M\) must be limited to avoid over-fitting, but from a Bayesian viewpoint there is no reason to restrict the number of parameters according to the size of the data set. Neal (1996) showed that for a broad class of priors the distribution over functions tends to a Gaussian process as \(M\to\infty\). The argument is the central limit theorem.

Proof

Write \(v_j=\sigma_v\zeta_j/\sqrt M\) with \(\zeta_j\) i.i.d., mean \(0\), variance \(1\), and let \(\mathbf{h}_j\)\({}=(h(\mathbf{x}_1;\mathbf{u}_j),\dots,h(\mathbf{x}_N;\mathbf{u}_j))^\T\in\R^N\). Then

\[ \big(f_M(\mathbf{x}_n)\big)_{n}=b\mathbf{1}+\sigma_v\frac{1}{\sqrt M}\sum_{j=1}^{M}\mathbf{z}_j,\qquad\mathbf{z}_j=\zeta_j\mathbf{h}_j. \]

The vectors \(\mathbf{z}_j\) are i.i.d. (they are functions of the independent pairs \((\zeta_j,\mathbf{u}_j)\)), with mean \(\E[\zeta_j]\E[\mathbf{h}_j]=\mathbf{0}\) and covariance \(\E[\zeta_j^2]\E[\mathbf{h}_j\mathbf{h}_j^\T]\)\({}=\bm\Sigma_h\), where \((\bm\Sigma_h)_{nm}\)\({}=\E_{\mathbf{u}}[h(\mathbf{x}_n;\mathbf{u})h(\mathbf{x}_m;\mathbf{u})]\), finite by the second-moment assumption and Cauchy-Schwarz. By the multivariate central limit theorem for i.i.d. vectors with finite covariance, \(M^{-1/2}\sum_j\mathbf{z}_j\)\({}\to\N(\mathbf{0},\bm\Sigma_h)\) in distribution. Independence of \(b\) from the rest, with \(b\mathbf{1}\)\({}\sim\N(\mathbf{0},\sigma_b^2\mathbf{1}\mathbf{1}^\T)\), adds the covariance \(\sigma_b^2\mathbf{1}\mathbf{1}^\T\), giving covariance \(\sigma_b^2\mathbf{1}\mathbf{1}^\T\)\({}+\sigma_v^2\bm\Sigma_h\)\({}=\mathbf{K}\).∎

Why the weights scale as \(1/\sqrt M\). Without the factor the variance of \(f_M\) would grow like \(M\) and no limit would exist; with it, each hidden unit contributes a vanishing share and the sum is Gaussian by the CLT. Other scalings (infinite-variance weight priors) give non-Gaussian limits (Neal, 1996).

Several outputs become independent. With \(K\) outputs \(f_k=b_k+\sum_jv_{jk}h(\mathbf{x};\mathbf{u}_j)\) and independent output weights, the same CLT gives a joint Gaussian limit whose cross-covariances \(\E[v_{jk}v_{jk'}]\) vanish for \(k\ne k'\); Gaussian and uncorrelated means independent. For finite \(M\) the outputs are dependent through the shared hidden units. In the limit this dependence disappears, so the "borrowing of statistical strength" between outputs, one of the merits of neural networks (hidden units are shaped by all the outputs), is lost in the Gaussian process limit; this is the remark in the book.

The kernel. Williams (1998) computed \(k\) in closed form for two activation functions. For the error function \(h(\mathbf{x};\mathbf{u})\)\({}=\operatorname{erf}(u_0+\sum_iu_ix_i)\) and a Gaussian prior \(\mathbf{u}\sim\N(\mathbf{0},\bm\Sigma)\), with the augmented input \(\tilde{\mathbf{x}}=(1,\mathbf{x}^\T)^\T\),

\[ \E_{\mathbf{u}}[h(\mathbf{x};\mathbf{u})h(\mathbf{x}';\mathbf{u})]=\frac2\pi\arcsin\!\Bigg(\frac{2\tilde{\mathbf{x}}^\T\bm\Sigma\tilde{\mathbf{x}}'}{\sqrt{(1+2\tilde{\mathbf{x}}^\T\bm\Sigma\tilde{\mathbf{x}})(1+2\tilde{\mathbf{x}}'^\T\bm\Sigma\tilde{\mathbf{x}}')}}\Bigg). \]

(I checked this against a Monte Carlo average with \(2\times10^6\) draws of \(\mathbf{u}\): the values agree to three decimals.) A second closed form exists for Gaussian-bump hidden units. These kernels are non-stationary: they cannot be written as a function of \(\mathbf{x}-\mathbf{x}'\), because the weight prior is centred at zero, which breaks translation invariance (as in the example above, where \(\tilde{\mathbf{x}}^\T\bm\Sigma\tilde{\mathbf{x}}'\) depends on both inputs).

Working directly with the covariance function marginalises over the weights implicitly. If the weight prior has hyperparameters, their values set the length scales of the function distribution (compare Figure 5.11 for a finite number of hidden units), and as in Section 6.4.3 they cannot be marginalised analytically but are handled by the techniques of Section 6.4.

Exercises

All 27 exercises of Chapter 6 (Bishop, 2006, pp. 320–323), restated in my own words, with complete solutions (hidden). Stars give the book's difficulty and "www" marks the exercises with official solutions. Many exercises are the derivations of the results above; their solutions summarise the argument and link to the full version.

1. Exercise 6.1 (★★, www). In the dual formulation of regularised least squares, show that the solution for the components \(a_n\) can be written as a linear combination of the elements of \(\bm\phi(\mathbf{x}_n)\). Calling the coefficients \(\mathbf{w}\), show that the dual of the dual is the original problem in \(\mathbf{w}\).

   Solution

   First part. From (6.4), \(a_n\)\({}=-\frac1\lambda\{\mathbf{w}^\T\bm\phi(\mathbf{x}_n)\)\({}-t_n\}\)\({}=\frac{t_n}\lambda\)\({}-\frac1\lambda\sum_iw_i\phi_i(\mathbf{x}_n)\), a linear combination of the elements \(\phi_i(\mathbf{x}_n)\) with coefficients \(-w_i/\lambda\), plus the offset \(t_n/\lambda\).

   A cleaner version (the official solution's route). The dual objective (6.7) can be written \(J(\mathbf{a})\)\({}=\frac12\|\mathbf{K}\mathbf{a}-\mathbf{t}\|^2\)\({}+\frac\lambda2\mathbf{a}^\T\mathbf{K}\mathbf{a}\). It does not change if \(\mathbf{a}\) is replaced by \(\mathbf{a}+\mathbf{z}\) with \(\mathbf{K}\mathbf{z}=\mathbf{0}\) (then \(\mathbf{K}(\mathbf{a}+\mathbf{z})\)\({}=\mathbf{K}\mathbf{a}\) and \((\mathbf{a}+\mathbf{z})^\T\mathbf{K}(\mathbf{a}+\mathbf{z})\)\({}=\mathbf{a}^\T\mathbf{K}\mathbf{a}\)). So the component of \(\mathbf{a}\) in the null space of \(\mathbf{K}\) is not determined by \(J\); remove the ambiguity by taking \(\mathbf{a}\) orthogonal to the null space, that is, in the range of \(\mathbf{K}=\bm\Phi\bm\Phi^\T\), which equals the column space of \(\bm\Phi\) (for a symmetric matrix the range is the orthogonal complement of the null space). Then \(\mathbf{a}=\bm\Phi\mathbf{u}\) for some \(\mathbf{u}\in\R^M\), that is \(a_n=\sum_iu_i\phi_i(\mathbf{x}_n)\), a linear combination of the elements of \(\bm\phi(\mathbf{x}_n)\) with coefficient vector \(\mathbf{u}\).

   Dual of the dual. Substitute \(\mathbf{a}=\bm\Phi\mathbf{u}\) into (6.7) and put \(\mathbf{w}=\bm\Phi^\T\bm\Phi\mathbf{u}\) (a bijection between \(\mathbf{u}\) and \(\mathbf{w}\) if \(\bm\Phi^\T\bm\Phi\) is invertible, i.e. \(\bm\Phi\) has full column rank). Then \(\mathbf{K}\mathbf{a}\)\({}=\bm\Phi\bm\Phi^\T\bm\Phi\mathbf{u}\)\({}=\bm\Phi\mathbf{w}\) and \(\mathbf{a}^\T\mathbf{K}\mathbf{a}\)\({}=\mathbf{u}^\T\bm\Phi^\T\bm\Phi\bm\Phi^\T\bm\Phi\mathbf{u}\)\({}=\mathbf{w}^\T\mathbf{w}\), so

   \[ J=\tfrac12\|\bm\Phi\mathbf{w}-\mathbf{t}\|^2+\tfrac\lambda2\mathbf{w}^\T\mathbf{w}, \]

   which is the original error function (6.2). (An alternative is Lagrange duality: minimise \(\frac12\|\bm\xi\|^2\)\({}+\frac\lambda2\|\mathbf{w}\|^2\) subject to \(\bm\xi=\bm\Phi\mathbf{w}-\mathbf{t}\); the Lagrangian dual in the multiplier \(\mathbf{a}=-\bm\xi/\lambda\) is, up to scaling, maximisation of \(\lambda\mathbf{a}^\T\mathbf{t}\)\({}-\frac\lambda2\mathbf{a}^\T(\mathbf{K}+\lambda\mathbf{I})\mathbf{a}\), whose maximiser is (6.8).) See also Gap 6.1.

2. Exercise 6.2 (★★). Derive a dual formulation of the perceptron learning algorithm: show that the learned weight vector is a linear combination of the \(t_n\bm\phi(\mathbf{x}_n)\), with \(t_n\in\{-1,+1\}\), find the algorithm and the prediction function in terms of the coefficients \(\alpha_n\), and show that \(\bm\phi\) appears only through \(k\).

   Solution

   The perceptron rule (4.55) updates, on each misclassified pattern \(n\), \(\mathbf{w}\leftarrow\mathbf{w}\)\({}+\eta\,t_n\bm\phi(\mathbf{x}_n)\) (take \(\eta=1\); another \(\eta\) rescales \(\mathbf{w}\) and leaves the classifier \(\operatorname{sign}(\mathbf{w}^\T\bm\phi)\) unchanged when \(\mathbf{w}\) starts at zero). Starting from \(\mathbf{w}=\mathbf{0}\), after any number of updates \(\mathbf{w}\)\({}=\sum_n\alpha_nt_n\bm\phi(\mathbf{x}_n)\), where \(\alpha_n\ge0\) is the number of times pattern \(n\) has been misclassified (times \(\eta\)). The prediction is

   \[ y(\mathbf{x})=\operatorname{sign}\big(\mathbf{w}^\T\bm\phi(\mathbf{x})\big)=\operatorname{sign}\Big(\sum_{n}\alpha_nt_n\,k(\mathbf{x}_n,\mathbf{x})\Big), \]

   and the algorithm becomes: initialise \(\bm\alpha=\mathbf{0}\); cycle through the patterns; for pattern \(m\), if \(t_m\sum_n\alpha_nt_nk(\mathbf{x}_n,\mathbf{x}_m)\)\({}\le0\) (misclassified), set \(\alpha_m\leftarrow\alpha_m+1\). The feature vector never appears except through \(k(\mathbf{x}_n,\mathbf{x}_m)\)\({}=\bm\phi(\mathbf{x}_n)^\T\bm\phi(\mathbf{x}_m)\), so any valid kernel can be used (the kernel perceptron). It converges for linearly separable data in the feature space (perceptron convergence theorem).

3. Exercise 6.3 (★). Write the nearest-neighbour classifier in terms of scalar products and use kernel substitution to get a nearest-neighbour classifier for a general kernel.

   Solution

   \(\|\mathbf{x}-\mathbf{x}_n\|^2\)\({}=\mathbf{x}^\T\mathbf{x}\)\({}-2\mathbf{x}^\T\mathbf{x}_n\)\({}+\mathbf{x}_n^\T\mathbf{x}_n\). Replacing each scalar product by a kernel,

   \[ d^2(\mathbf{x},\mathbf{x}_n)=k(\mathbf{x},\mathbf{x})-2k(\mathbf{x},\mathbf{x}_n)+k(\mathbf{x}_n,\mathbf{x}_n), \]

   the squared distance in feature space. Assign \(\mathbf{x}\) the class of \(\arg\min_nd^2(\mathbf{x},\mathbf{x}_n)\). The term \(k(\mathbf{x},\mathbf{x})\) is the same for all \(n\) and can be dropped; for kernels with \(k(\mathbf{x},\mathbf{x})=\) const (e.g. Gaussian) the rule reduces to picking the training point with the largest \(k(\mathbf{x},\mathbf{x}_n)\).

4. Exercise 6.4 (★). Find a \(2\times2\) matrix with positive eigenvalues but at least one negative element (the converse of the example of a positive matrix with a negative eigenvalue).

   Solution

   \(\mathbf{A}\)\({}=\begin{pmatrix}2&-1\\-1&2\end{pmatrix}\): the characteristic polynomial is \((2-\lambda)^2-1\), so \(\lambda=1,3>0\) (equivalently \(\mathbf{v}^\T\mathbf{A}\mathbf{v}\)\({}=v_1^2+v_2^2+(v_1-v_2)^2>0\)), and the off-diagonal entries are negative. In the other direction \(\begin{pmatrix}1&2\\2&1\end{pmatrix}\) has positive entries and eigenvalues \(3,-1\). Positive semidefiniteness is a property of the quadratic form, not of the signs of the entries.

5. Exercise 6.5 (★, www). Verify (6.13) and (6.14).

   Solution

   If \(k_1(\mathbf{x},\mathbf{x}')\)\({}=\bm\phi(\mathbf{x})^\T\bm\phi(\mathbf{x}')\), then \(ck_1\)\({}=\mathbf{u}(\mathbf{x})^\T\mathbf{u}(\mathbf{x}')\) with \(\mathbf{u}=\sqrt c\,\bm\phi\), and \(f(\mathbf{x})k_1(\mathbf{x},\mathbf{x}')f(\mathbf{x}')\)\({}=\mathbf{v}(\mathbf{x})^\T\mathbf{v}(\mathbf{x}')\) with \(\mathbf{v}=f\bm\phi\). Both are inner products of feature vectors, hence valid. In terms of Gram matrices: \(c\mathbf{K}_1\) and \(\mathbf{D}\mathbf{K}_1\mathbf{D}\), \(\mathbf{D}=\diag(f(\mathbf{x}_n))\), are PSD. See Result 6.3.

6. Exercise 6.6 (★). Verify (6.15) and (6.16).

   Solution

   (6.15): \(q(k_1)=\sum_jc_jk_1^j\) with \(c_j\ge0\). By (6.18) every power \(k_1^j\) is valid (the constant \(k_1^0=1\) is the inner product of the one-dimensional feature \(1\)), by (6.13) so is each \(c_jk_1^j\), and by (6.17) so is the finite sum. (6.16): \(\exp k_1=\lim_J\sum_{j\le J}k_1^j/j!\), a limit of kernels of the form (6.15); Gram matrices converge entrywise, and a limit of PSD matrices is PSD because \(\mathbf{a}^\T\mathbf{K}_J\mathbf{a}\ge0\) for all \(J\) passes to the limit. Full details in Result 6.3.

7. Exercise 6.7 (★, www). Verify (6.17) and (6.18).

   Solution

   (6.17): \(\mathbf{a}^\T(\mathbf{K}_1+\mathbf{K}_2)\mathbf{a}\)\({}=\mathbf{a}^\T\mathbf{K}_1\mathbf{a}\)\({}+\mathbf{a}^\T\mathbf{K}_2\mathbf{a}\)\({}\ge0\). (6.18): with feature vectors \(\bm\phi\in\R^M\), \(\bm\psi\in\R^L\) for \(k_1,k_2\), \(k_1k_2\)\({}=\sum_{m=1}^{M}\sum_{l=1}^{L}\phi_m(\mathbf{x})\phi_m(\mathbf{x}')\psi_l(\mathbf{x})\psi_l(\mathbf{x}')\)\({}=\sum_{m,l}[\phi_m\psi_l](\mathbf{x})[\phi_m\psi_l](\mathbf{x}')\), an inner product of the \(ML\)-vector of products (the official solution indexes this by \(k=(m-1)L+l\)). For infinite-dimensional feature spaces use the Hadamard (Schur) product argument of Result 6.3.

8. Exercise 6.8 (★). Verify (6.19) and (6.20).

   Solution

   (6.19): the Gram matrix of \(k_3(\bm\phi(\mathbf{x}),\bm\phi(\mathbf{x}'))\) at points \(\mathbf{x}_n\) is the Gram matrix of \(k_3\) at the points \(\bm\phi(\mathbf{x}_n)\), PSD by validity of \(k_3\). (6.20): for \(\mathbf{A}\) symmetric PSD, \(\mathbf{A}=\mathbf{L}^\T\mathbf{L}\) (e.g. \(\mathbf{L}\)\({}=\mathbf{A}^{1/2}\)\({}=\mathbf{U}\bm\Lambda^{1/2}\mathbf{U}^\T\)), so \(\mathbf{x}^\T\mathbf{A}\mathbf{x}'\)\({}=(\mathbf{L}\mathbf{x})^\T(\mathbf{L}\mathbf{x}')\). See Result 6.3.

9. Exercise 6.9 (★). Verify (6.21) and (6.22).

   Solution

   \(\tilde k_a(\mathbf{x},\mathbf{x}')\)\({}=k_a(\mathbf{x}_a,\mathbf{x}_a')\) is (6.19) with \(\bm\phi(\mathbf{x})=\mathbf{x}_a\), so valid on the whole space; likewise \(\tilde k_b\). The sum is then (6.17) and the product (6.18). The parts \(\mathbf{x}_a,\mathbf{x}_b\) may overlap. See Result 6.3.

10. Exercise 6.10 (★). Show that \(k(\mathbf{x},\mathbf{x}')\)\({}=f(\mathbf{x})f(\mathbf{x}')\) is an excellent kernel for learning \(f\), by showing that a linear learning machine on this kernel always finds a solution proportional to \(f(\mathbf{x})\).

    Solution

    The feature map is the scalar \(\phi(\mathbf{x})=f(\mathbf{x})\), so any linear machine \(y(\mathbf{x})=w\phi(\mathbf{x})=wf(\mathbf{x})\) is proportional to \(f\), for every \(w\). For instance, kernel ridge regression (6.9) with \(\mathbf{K}=\mathbf{f}\mathbf{f}^\T\) (\(\mathbf{f}=(f(\mathbf{x}_n))_n\)): since \((\mathbf{f}\mathbf{f}^\T+\lambda\mathbf{I})\mathbf{f}\)\({}=(\mathbf{f}^\T\mathbf{f}+\lambda)\mathbf{f}\), we get \((\mathbf{f}\mathbf{f}^\T+\lambda\mathbf{I})^{-1}\mathbf{f}\)\({}=\mathbf{f}/(\lambda+\mathbf{f}^\T\mathbf{f})\) and \(y(\mathbf{x})\)\({}=\mathbf{k}(\mathbf{x})^\T(\mathbf{K}+\lambda\mathbf{I})^{-1}\mathbf{t}\)\({}=f(\mathbf{x})\,\mathbf{f}^\T\mathbf{t}/(\lambda+\mathbf{f}^\T\mathbf{f})\), proportional to \(f\). The Gram matrix has rank one, the kernel encodes the answer in advance, and "learning" only fits one scale factor; the kernel is excellent only if you already know \(f\).

11. Exercise 6.11 (★). Using the expansion (6.25) and a power series for the middle factor, show that the Gaussian kernel is the inner product of an infinite-dimensional feature vector.

    Solution

    Done in Gap 6.3: \(\phi_{\bm\alpha}(\mathbf{x})\)\({}=e^{-\|\mathbf{x}\|^2/2\sigma^2}\mathbf{x}^{\bm\alpha}/(\sigma^{|\bm\alpha|}\sqrt{\bm\alpha!})\), indexed by all multi-indices; for \(D=1\), \(\phi_j(x)\)\({}=e^{-x^2/2\sigma^2}x^j/(\sigma^j\sqrt{j!})\) for \(j=0,1,\dots\), and \(\sum_j\phi_j(x)\phi_j(x')\)\({}=e^{-(x-x')^2/2\sigma^2}\).

12. Exercise 6.12 (★★, www). For subsets \(A\) of a fixed finite set \(D\), show that (6.27) is an inner product in a feature space of dimension \(2^{|D|}\) with \(\phi_U(A)=1\) if \(U\subseteq A\) and \(0\) otherwise (6.95).

    Solution

    \(\sum_{U\subseteq D}\phi_U(A_1)\phi_U(A_2)\)\({}=\#\{U:U\subseteq A_1\text{ and }U\subseteq A_2\}\)\({}=\#\{U:U\subseteq A_1\cap A_2\}\)\({}=2^{|A_1\cap A_2|}\), since a set with \(r\) elements has \(2^r\) subsets (including the empty set and itself). This is (6.27); the feature space has one coordinate per subset \(U\) of \(D\), \(2^{|D|}\) in all. (Note the errata: \(|A|\) is the number of elements.) The kernel counts the common sub-patterns of the two sets. See also Section 6.2.

13. Exercise 6.13 (★). Show that the Fisher kernel (6.33) is unchanged by an invertible differentiable reparametrisation \(\bm\theta\to\bm\psi(\bm\theta)\).

    Solution

    With \(\mathbf{J}=\partial\bm\theta/\partial\bm\psi\), \(\mathbf{g}_\psi=\mathbf{J}^\T\mathbf{g}_\theta\) and \(\mathbf{F}_\psi\)\({}=\mathbf{J}^\T\mathbf{F}_\theta\mathbf{J}\), so \(k_\psi\)\({}=\mathbf{g}_\theta^\T\mathbf{J}(\mathbf{J}^\T\mathbf{F}_\theta\mathbf{J})^{-1}\mathbf{J}^\T\mathbf{g}'_\theta\)\({}=\mathbf{g}_\theta^\T\mathbf{F}_\theta^{-1}\mathbf{g}'_\theta\)\({}=k_\theta\) because \(\mathbf{J}\) is square and invertible. See Gap 6.4.

14. Exercise 6.14 (★, www). Write down the Fisher kernel for \(p(\mathbf{x}\given\bm\mu)\)\({}=\N(\mathbf{x}\given\bm\mu,\mathbf{S})\) with \(\mathbf{S}\) fixed.

    Solution

    The score is \(\mathbf{g}(\bm\mu,\mathbf{x})\)\({}=\nabla_{\bm\mu}\ln\N(\mathbf{x}\given\bm\mu,\mathbf{S})\)\({}=\mathbf{S}^{-1}(\mathbf{x}-\bm\mu)\). The Fisher information is \(\mathbf{F}\)\({}=\E[\mathbf{g}\mathbf{g}^\T]\)\({}=\mathbf{S}^{-1}\E[(\mathbf{x}-\bm\mu)(\mathbf{x}-\bm\mu)^\T]\mathbf{S}^{-1}\)\({}=\mathbf{S}^{-1}\mathbf{S}\mathbf{S}^{-1}\)\({}=\mathbf{S}^{-1}\). Hence \(k(\mathbf{x},\mathbf{x}')\)\({}=(\mathbf{x}-\bm\mu)^\T\mathbf{S}^{-1}\mathbf{S}\,\mathbf{S}^{-1}(\mathbf{x}'-\bm\mu)\), that is

    \[ k(\mathbf{x},\mathbf{x}')=(\mathbf{x}-\bm\mu)^\T\mathbf{S}^{-1}(\mathbf{x}'-\bm\mu), \]

    the Mahalanobis inner product (its diagonal is the squared Mahalanobis distance to the mean), i.e. the linear kernel after centring and whitening.

15. Exercise 6.15 (★). By considering the determinant of a \(2\times2\) Gram matrix, show that a positive-definite kernel satisfies the Cauchy-Schwarz inequality (6.96).

    Solution

    For two points the Gram matrix is \(\begin{pmatrix}k(\mathbf{x}_1,\mathbf{x}_1)&k(\mathbf{x}_1,\mathbf{x}_2)\\k(\mathbf{x}_1,\mathbf{x}_2)&k(\mathbf{x}_2,\mathbf{x}_2)\end{pmatrix}\), symmetric PSD, so its eigenvalues are \(\ge0\) and so is their product, the determinant: \(k(\mathbf{x}_1,\mathbf{x}_1)k(\mathbf{x}_2,\mathbf{x}_2)\)\({}-k(\mathbf{x}_1,\mathbf{x}_2)^2\)\({}\ge0\), which is \(k(\mathbf{x}_1,\mathbf{x}_2)^2\)\({}\le k(\mathbf{x}_1,\mathbf{x}_1)k(\mathbf{x}_2,\mathbf{x}_2)\) (6.96). (It is the Cauchy-Schwarz inequality \(|\langle\bm\phi_1,\bm\phi_2\rangle|\)\({}\le\|\bm\phi_1\|\|\bm\phi_2\|\) in the feature space; the strict version holds for positive definite Gram matrices.) This is the quick test that rules out the sigmoidal kernel in Gap 6.5.

16. Exercise 6.16 (★★). Let \(J(\mathbf{w})\)\({}=f(\mathbf{w}^\T\bm\phi(\mathbf{x}_1),\dots,\mathbf{w}^\T\bm\phi(\mathbf{x}_N))\)\({}+g(\mathbf{w}^\T\mathbf{w})\) with \(g\) monotonically increasing. Writing \(\mathbf{w}\)\({}=\sum_n\alpha_n\bm\phi(\mathbf{x}_n)\)\({}+\mathbf{w}_\perp\) with \(\mathbf{w}_\perp^\T\bm\phi(\mathbf{x}_n)=0\) for all \(n\) (the errata add this condition), show that the minimiser is a linear combination of the \(\bm\phi(\mathbf{x}_n)\) (the representer theorem).

    Solution

    Any \(\mathbf{w}\) splits uniquely into its orthogonal projection onto \(\operatorname{span}\{\bm\phi(\mathbf{x}_n)\}\) and a component \(\mathbf{w}_\perp\) orthogonal to that span. Then \(\mathbf{w}^\T\bm\phi(\mathbf{x}_m)\)\({}=\sum_n\alpha_n\bm\phi(\mathbf{x}_n)^\T\bm\phi(\mathbf{x}_m)\): the first term \(f\) of \(J\) does not depend on \(\mathbf{w}_\perp\). The second satisfies \(\mathbf{w}^\T\mathbf{w}\)\({}=\|\sum_n\alpha_n\bm\phi(\mathbf{x}_n)\|^2\)\({}+\|\mathbf{w}_\perp\|^2\)\({}\ge\|\sum_n\alpha_n\bm\phi(\mathbf{x}_n)\|^2\) (Pythagoras), so \(g(\mathbf{w}^\T\mathbf{w})\)\({}\ge g(\|\sum_n\alpha_n\bm\phi(\mathbf{x}_n)\|^2)\) because \(g\) is increasing, with equality only if \(\mathbf{w}_\perp=\mathbf{0}\) when \(g\) is strictly increasing. Hence dropping \(\mathbf{w}_\perp\) never increases \(J\) (and strictly decreases it unless \(\mathbf{w}_\perp=\mathbf{0}\)), so a minimiser, and for strictly increasing \(g\) every minimiser, lies in the span: \(\mathbf{w}=\sum_n\alpha_n\bm\phi(\mathbf{x}_n)\). Then \(J\) depends on the data only through the Gram matrix, since \(\mathbf{w}^\T\bm\phi(\mathbf{x}_m)\)\({}=(\mathbf{K}\bm\alpha)_m\) and \(\mathbf{w}^\T\mathbf{w}\)\({}=\bm\alpha^\T\mathbf{K}\bm\alpha\). Ridge regression (6.2) is the case \(f=\frac12\sum(\cdot-t_n)^2\), \(g(z)=\lambda z/2\); SVMs (Chapter 7) are another. The result is due to Kimeldorf and Wahba (1971); the general form to Schölkopf, Herbrich and Smola (2001).

17. Exercise 6.17 (★★, www). Minimise (6.39) with respect to \(y(\mathbf{x})\) by the calculus of variations and show that the solution is (6.40) with basis functions (6.41).

    Solution

    See Gap 6.6. In first-principles form (the official solution): perturb \(y\to y+\epsilon\eta\), set the coefficient of \(\epsilon\) to zero, \(\sum_n\int\{y(\mathbf{x}_n+\bm\xi)\)\({}-t_n\}\eta(\mathbf{x}_n+\bm\xi)\nu(\bm\xi)\d\bm\xi\)\({}=0\) for all \(\eta\), and take \(\eta(\mathbf{x})=\delta(\mathbf{x}-\mathbf{z})\): this gives \(\sum_n\{y(\mathbf{z})\)\({}-t_n\}\nu(\mathbf{z}-\mathbf{x}_n)\)\({}=0\), hence (6.40)–(6.41). The same follows by changing variables to \(\mathbf{z}=\mathbf{x}_n+\bm\xi\) and minimising the integrand pointwise.

18. Exercise 6.18 (★). For a Nadaraya-Watson model with one input and one target and isotropic Gaussian components of covariance \(\sigma^2\mathbf{I}\), write \(p(t\given x)\), \(\E[t\given x]\) and \(\var[t\given x]\) in terms of \(k(x,x_n)\).

    Solution

    \(g(x)\)\({}=\int\N(x\given0,\sigma^2)\N(t\given0,\sigma^2)\d t\)\({}=\N(x\given0,\sigma^2)\), so \(k(x,x_n)\)\({}=\N(x-x_n\given0,\sigma^2)/\sum_m\N(x-x_m\given0,\sigma^2)\) and, from (6.48),

    \[ p(t\given x)=\sum_nk(x,x_n)\N(t\given t_n,\sigma^2),\quad\E[t\given x]=\sum_nk(x,x_n)t_n,\quad\var[t\given x]=\sigma^2+\sum_nk(x,x_n)t_n^2-\Big(\sum_nk(x,x_n)t_n\Big)^2. \]

    The variance is the mixture variance (law of total variance); see Gap 6.7.

19. Exercise 6.19 (★★). Targets are \(t_n=y(\mathbf{z}_n)+\)noise, but only \(\mathbf{x}_n=\mathbf{z}_n+\bm\xi_n\), with \(\bm\xi_n\sim g\), is observed. Use the error \(E\)\({}=\frac12\sum_n\int\{y(\mathbf{x}_n-\bm\xi_n)\)\({}-t_n\}^2g(\bm\xi_n)\d\bm\xi_n\) (6.99) and the calculus of variations to show that the optimal \(y\) is a Nadaraya-Watson regression of the form (6.45), (6.46).

    Solution

    Substitute \(\mathbf{z}=\mathbf{x}_n-\bm\xi_n\) in the \(n\)-th term: \(E\)\({}=\frac12\int\sum_ng(\mathbf{x}_n-\mathbf{z})\{y(\mathbf{z})\)\({}-t_n\}^2\d\mathbf{z}\). The integrand involves \(y\) only at the point \(\mathbf{z}\), so setting its derivative with respect to \(y(\mathbf{z})\) to zero (equivalently the functional derivative, as in Gap 6.6) gives \(\sum_ng(\mathbf{x}_n-\mathbf{z})\{y(\mathbf{z})\)\({}-t_n\}\)\({}=0\), hence

    \[ y(\mathbf{x})=\sum_nt_n\frac{g(\mathbf{x}_n-\mathbf{x})}{\sum_mg(\mathbf{x}_m-\mathbf{x})}=\sum_nk(\mathbf{x},\mathbf{x}_n)t_n. \]

    This is (6.45) with the kernel (6.46) when \(g\) is symmetric, \(g(-\bm\xi)=g(\bm\xi)\) (for instance any isotropic noise), so that \(g(\mathbf{x}_n-\mathbf{x})\)\({}=g(\mathbf{x}-\mathbf{x}_n)\). It is a minimum because the integrand is a convex quadratic in \(y(\mathbf{z})\).

20. Exercise 6.20 (★★, www). Verify (6.66) and (6.67).

    Solution

    See Result 6.4: identify \(t_{N+1}\) with \(\mathbf{x}_a\) and \(\mathbf{t}\) with \(\mathbf{x}_b\) in the conditional Gaussian formulas (2.81)–(2.82), with \(\bm\Sigma_{aa}=c\), \(\bm\Sigma_{ab}=\mathbf{k}^\T\), \(\bm\Sigma_{bb}=\mathbf{C}_N\), zero means; the second proof there completes the square with the block inverse.

21. Exercise 6.21 (★★, www). For a GP whose kernel is built from a fixed set of basis functions, show that the predictive distribution equals the Bayesian linear regression result (3.58), by comparing means (using (C.6)) and variances (using (C.7)).

    Solution

    See Gap 6.10: with \(\mathbf{C}_N\)\({}=\alpha^{-1}\bm\Phi\bm\Phi^\T\)\({}+\beta^{-1}\mathbf{I}\), the push-through identity (C.6) gives \(m(\mathbf{x}_{N+1})\)\({}=\beta\bm\phi^\T\mathbf{S}_N\bm\Phi^\T\mathbf{t}\)\({}=\bm\phi^\T\mathbf{m}_N\), and the Woodbury identity (C.7) gives \(\sigma^2\)\({}=\beta^{-1}+\bm\phi^\T\mathbf{S}_N\bm\phi\), both with \(\mathbf{S}_N^{-1}\)\({}=\alpha\mathbf{I}+\beta\bm\Phi^\T\bm\Phi\), as in (3.58).

22. Exercise 6.22 (★★). With \(N\) training inputs and \(L\) test inputs and a GP prior over \(t(\mathbf{x})\), derive the joint predictive distribution of \(t(\mathbf{x}_{N+1}),\dots,t(\mathbf{x}_{N+L})\) given the training targets, and show that the marginal for one test target is (6.66)–(6.67).

    Solution

    See Gap 6.11: the joint predictive is \(\N(\mathbf{K}_*^\T\mathbf{C}_N^{-1}\mathbf{t},\;\mathbf{K}_{**}+\beta^{-1}\mathbf{I}_L-\mathbf{K}_*^\T\mathbf{C}_N^{-1}\mathbf{K}_*)\), and the \(j\)-th marginal takes the \(j\)-th mean entry and \(j\)-th diagonal covariance entry.

23. Exercise 6.23 (★★, www). For a GP regression model with \(D\)-dimensional targets, write down the conditional distribution of \(\mathbf{t}_{N+1}\) given training inputs \(\mathbf{x}_1,\dots,\mathbf{x}_N\) (the errata fix the index) and targets \(\mathbf{t}_1,\dots,\mathbf{t}_N\).

    Solution

    Assuming the \(D\) targets are independent given the input, with a common kernel, \(p(\mathbf{t}_{N+1}\given\mathbf{T})\)\({}=\N(\mathbf{t}_{N+1}\given\mathbf{m}(\mathbf{x}_{N+1}),\sigma^2(\mathbf{x}_{N+1})\mathbf{I}_D)\), where \(\mathbf{T}\) is the \(N\times D\) matrix with rows \(\mathbf{t}_n^\T\), \(\mathbf{m}(\mathbf{x}_{N+1})^\T\)\({}=\mathbf{k}^\T\mathbf{C}_N^{-1}\mathbf{T}\), and \(\sigma^2\) is (6.67). \(\mathbf{C}_N\) depends on the inputs only, so it is shared. See Gap 6.11.

24. Exercise 6.24 (★). Show that a diagonal matrix \(\mathbf{W}\) with \(0

    Solution

    For \(\mathbf{v}\ne\mathbf{0}\): \(\mathbf{v}^\T\mathbf{W}\mathbf{v}\)\({}=\sum_iW_{ii}v_i^2>0\), since each \(W_{ii}>0\) and some \(v_i\ne0\). For PD \(\mathbf{A},\mathbf{B}\): \(\mathbf{v}^\T(\mathbf{A}+\mathbf{B})\mathbf{v}\)\({}=\mathbf{v}^\T\mathbf{A}\mathbf{v}\)\({}+\mathbf{v}^\T\mathbf{B}\mathbf{v}>0\). Application: \(\mathbf{W}_N+\mathbf{C}_N^{-1}\) is PD, hence \(\Psi\) is strictly concave (Gap 6.16).

25. Exercise 6.25 (★, www). Derive the update (6.83) from the Newton-Raphson formula (4.92).

    Solution

    \(\mathbf{a}^{\rm new}\)\({}=\mathbf{a}\)\({}+(\mathbf{C}^{-1}+\mathbf{W})^{-1}(\mathbf{t}-\bm\sigma-\mathbf{C}^{-1}\mathbf{a})\)\({}=(\mathbf{C}^{-1}+\mathbf{W})^{-1}[\mathbf{t}-\bm\sigma+\mathbf{W}\mathbf{a}]\)\({}=\mathbf{C}(\mathbf{I}+\mathbf{W}\mathbf{C})^{-1}[\mathbf{t}-\bm\sigma+\mathbf{W}\mathbf{a}]\), using \((\mathbf{C}^{-1}+\mathbf{W})^{-1}\)\({}=\mathbf{C}(\mathbf{I}+\mathbf{W}\mathbf{C})^{-1}\); full derivation in Gap 6.17.

26. Exercise 6.26 (★). Using (2.115), derive (6.87) and (6.88).

    Solution

    See Gap 6.18: apply (2.115) to the Gaussian \(q(\mathbf{a}_N)\)\({}=\N(\mathbf{a}^\star,\mathbf{H}^{-1})\) and the conditional (6.78) to get mean \(\mathbf{k}^\T\mathbf{C}^{-1}\mathbf{a}^\star\)\({}=\mathbf{k}^\T(\mathbf{t}-\bm\sigma)\) and variance \(c-\mathbf{k}^\T\mathbf{C}^{-1}\mathbf{k}\)\({}+\mathbf{k}^\T\mathbf{C}^{-1}\mathbf{H}^{-1}\mathbf{C}^{-1}\mathbf{k}\)\({}=c\)\({}-\mathbf{k}^\T(\mathbf{W}^{-1}+\mathbf{C})^{-1}\mathbf{k}\) by Woodbury.

27. Exercise 6.27 (★★★). Derive the Laplace approximation (6.90) to the log likelihood of the GP classifier, and the gradient terms (6.91), (6.92), (6.94).

    Solution

    (6.90): Gap 6.19 (Gaussian integral of \(e^{\Psi}\) around the mode). The gradient: Gap 6.20, where \(\nabla\Psi(\mathbf{a}^\star)=\mathbf{0}\) removes the dependence of \(\Psi(\mathbf{a}^\star)\) on the moving mode, the explicit parts of \(\Psi\) and \(-\frac12\ln|\mathbf{H}|\) combine into (6.91), the dependence of \(\mathbf{W}\) on \(\mathbf{a}^\star\) gives (6.92), and implicit differentiation of \(\mathbf{a}^\star\)\({}=\mathbf{C}(\mathbf{t}-\bm\sigma)\) gives (6.93)–(6.94). Two factor-order slips in the printed (6.91) and (6.94) are corrected there and verified by finite differences.

Demonstrations for this chapter

Each demonstration runs in the browser, with seeded random numbers so that a class sees the same pictures, and links back to these notes. The ▶ Demo boxes above open them already set to the case under discussion.

The Python code that produces every figure of this page, and checks the chapter's identities numerically, is code/ch06_kernel_methods.py (see the code index).

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 6, pp. 291–323 (and the official errata).
2. Aizerman, M. A., Braverman, E. M. and Rozonoer, L. I. (1964). Theoretical foundations of the potential function method in pattern recognition learning. Automation and Remote Control, 25, 821–837.
3. Aronszajn, N. (1950). Theory of reproducing kernels. Transactions of the American Mathematical Society, 68(3), 337–404. doi:10.1090/S0002-9947-1950-0051437-7.
4. Boser, B. E., Guyon, I. M. and Vapnik, V. N. (1992). A training algorithm for optimal margin classifiers. In Proceedings of the Fifth Annual Workshop on Computational Learning Theory, 144–152. doi:10.1145/130385.130401.
5. Broomhead, D. S. and Lowe, D. (1988). Multivariable functional interpolation and adaptive networks. Complex Systems, 2, 321–355.
6. Chen, S., Cowan, C. F. N. and Grant, P. M. (1991). Orthogonal least squares learning algorithm for radial basis function networks. IEEE Transactions on Neural Networks, 2(2), 302–309. doi:10.1109/72.80341.
7. Cressie, N. A. C. (1993). Statistics for Spatial Data, revised edition. Wiley, New York.
8. Ghahramani, Z. and Jordan, M. I. (1994). Supervised learning from incomplete data via an EM approach. In Advances in Neural Information Processing Systems, 6, 120–127.
9. Gibbs, M. N. and MacKay, D. J. C. (2000). Variational Gaussian process classifiers. IEEE Transactions on Neural Networks, 11(6), 1458–1464. doi:10.1109/72.883477.
10. Goldberg, P. W., Williams, C. K. I. and Bishop, C. M. (1998). Regression with input-dependent noise: a Gaussian process treatment. In Advances in Neural Information Processing Systems, 10, 493–499.
11. Haussler, D. (1999). Convolution kernels on discrete structures. Technical Report UCSC-CRL-99-10, University of California, Santa Cruz.
12. Jaakkola, T. S. and Haussler, D. (1999). Exploiting generative models in discriminative classifiers. In Advances in Neural Information Processing Systems, 11, 487–493.
13. Kimeldorf, G. and Wahba, G. (1971). Some results on Tchebycheffian spline functions. Journal of Mathematical Analysis and Applications, 33(1), 82–95. doi:10.1016/0022-247X(71)90184-3.
14. MacKay, D. J. C. (1998). Introduction to Gaussian processes. In C. M. Bishop (ed.), Neural Networks and Machine Learning, NATO ASI Series F, 168, 133–166. Springer.
15. Mercer, J. (1909). Functions of positive and negative type, and their connection with the theory of integral equations. Philosophical Transactions of the Royal Society A, 209, 415–446. doi:10.1098/rsta.1909.0016.
16. Minka, T. P. (2001). A family of algorithms for approximate Bayesian inference. PhD thesis, Massachusetts Institute of Technology.
17. Moody, J. and Darken, C. J. (1989). Fast learning in networks of locally-tuned processing units. Neural Computation, 1(2), 281–294. doi:10.1162/neco.1989.1.2.281.
18. Nadaraya, E. A. (1964). On estimating regression. Theory of Probability and its Applications, 9(1), 141–142. doi:10.1137/1109020.
19. Neal, R. M. (1996). Bayesian Learning for Neural Networks. Lecture Notes in Statistics 118, Springer, New York. doi:10.1007/978-1-4612-0745-0.
20. Neal, R. M. (1997). Monte Carlo implementation of Gaussian process models for Bayesian regression and classification. Technical Report 9702, Department of Statistics, University of Toronto.
21. Opper, M. and Winther, O. (2000). Gaussian processes for classification: mean-field algorithms. Neural Computation, 12(11), 2655–2684. doi:10.1162/089976600300014881.
22. Poggio, T. and Girosi, F. (1990). Networks for approximation and learning. Proceedings of the IEEE, 78(9), 1481–1497. doi:10.1109/5.58326.
23. Powell, M. J. D. (1987). Radial basis functions for multivariable interpolation: a review. In J. C. Mason and M. G. Cox (eds), Algorithms for Approximation, 143–167. Oxford University Press.
24. Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA. doi:10.7551/mitpress/3206.001.0001.
25. Schoenberg, I. J. (1942). Positive definite functions on spheres. Duke Mathematical Journal, 9(1), 96–108. doi:10.1215/S0012-7094-42-00908-6.
26. Schölkopf, B., Herbrich, R. and Smola, A. J. (2001). A generalized representer theorem. In Computational Learning Theory (COLT/EuroCOLT 2001), Lecture Notes in Computer Science 2111, 416–426. doi:10.1007/3-540-44581-1_27.
27. Schölkopf, B., Smola, A. J. and Müller, K.-R. (1998). Nonlinear component analysis as a kernel eigenvalue problem. Neural Computation, 10(5), 1299–1319. doi:10.1162/089976698300017467.
28. Shawe-Taylor, J. and Cristianini, N. (2004). Kernel Methods for Pattern Analysis. Cambridge University Press. doi:10.1017/CBO9780511809682.
29. Uhlenbeck, G. E. and Ornstein, L. S. (1930). On the theory of the Brownian motion. Physical Review, 36(5), 823–841. doi:10.1103/PhysRev.36.823.
30. Watson, G. S. (1964). Smooth regression analysis. Sankhyā: The Indian Journal of Statistics, Series A, 26(4), 359–372.
31. Williams, C. K. I. (1998). Computation with infinite neural networks. Neural Computation, 10(5), 1203–1216. doi:10.1162/089976698300017412.
32. Williams, C. K. I. and Barber, D. (1998). Bayesian classification with Gaussian processes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 20(12), 1342–1351. doi:10.1109/34.735807.

DOIs were checked against Crossref. The numerical statements made in the text (the Gram-matrix counterexamples to the sigmoidal and geodesic kernels, the Laplace evidence gradient (6.91)–(6.94) against finite differences, the predictive variance identity (6.88), Williams' arcsine kernel against Monte Carlo) come from short Python computations written for these notes; the figures come from the course code.