Linear Models for Regression

Regression predicts a continuous target from inputs. This chapter studies the simplest flexible family: models that are linear in their parameters but may be nonlinear in the inputs, through a fixed set of basis functions. Because the parameters enter linearly, everything is available in closed form: maximum likelihood is least squares (a projection), regularisation is a ridge or lasso penalty, the Bayesian posterior is Gaussian, the predictive distribution is Gaussian, and the marginal likelihood (the evidence) is an explicit determinant that can be maximised to set the regularisation strength from the training data alone. The chapter is also where the bias-variance trade-off and Bayesian model comparison first appear. Every step the book leaves to the reader is written out below.

TrungTin Nguyen · Following Bishop (2006), Pattern Recognition and Machine Learning, Chapter 3 · 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

3.1 Linear basis function models

pp. 137–139 The goal of regression is to predict a continuous target \(t\) from a \(D\)-dimensional input \(\mathbf{x}\), after seeing \(N\) training pairs. The simplest model is linear in the inputs, \(y(\mathbf{x},\mathbf{w})=w_0+w_1x_1+\dots+w_Dx_D\) (3.1), which is far too restrictive. The key idea of the chapter is to stay linear in the parameters while making the model nonlinear in the inputs, by feeding \(\mathbf{x}\) through fixed nonlinear functions \(\phi_j\).

Two remarks. First, "linear" means linear in \(\mathbf{w}\): the model can represent any function in the span of the \(\phi_j\), including highly nonlinear ones, and the learning problems below are convex. Second, the choice of basis is not unique in its effect on the parameters, but sometimes only in a trivial way. For instance, the logistic sigmoid and the hyperbolic tangent are related by \(\tanh a=2\sigma(2a)-1\), so a linear combination of sigmoids is a linear combination of \(\tanh\) functions with re-scaled weights and a shifted bias, and the two model families are the same set of functions (Exercise 3.1).

3.1.1 Maximum likelihood and least squares

pp. 140–143 Assume the target is a deterministic function plus Gaussian noise, \(t=y(\mathbf{x},\mathbf{w})+\epsilon\) (3.7), with \(\epsilon\sim\N(0,\beta^{-1})\), i.e.

\[ p(t\given\mathbf{x},\mathbf{w},\beta)=\N\big(t\given y(\mathbf{x},\mathbf{w}),\beta^{-1}\big). \tag{3.8} \]

Under squared loss the best prediction is the conditional mean, here \(\E[t\given\mathbf{x}]\)\({}=\int t\,p(t\given\mathbf{x})\,\d t\)\({}=y(\mathbf{x},\mathbf{w})\) (3.9). The Gaussian assumption makes the conditional distribution unimodal; mixtures of conditional Gaussians relax this (Section 14.5.1). For independent observations the likelihood (inputs are always conditioned on, never modelled) is

\[ p(\mathbf{t}\given\mathbf{X},\mathbf{w},\beta)=\prod_{n=1}^{N}\N\big(t_n\given\mathbf{w}^\T\bm\phi(\mathbf{x}_n),\beta^{-1}\big). \tag{3.10} \]

Proof

Log likelihood. Taking logs in (3.10) and using \(\ln\N(t\given\mu,\beta^{-1})\)\({}=\tfrac12\ln\beta-\tfrac12\ln2\pi\)\({}-\tfrac\beta2(t-\mu)^2\) with \(\mu=\mathbf{w}^\T\bm\phi_n\), summed over \(n\), gives (3.11).

Gradient. With \(r_n(\mathbf{w})=t_n-\mathbf{w}^\T\bm\phi_n\) we have \(\nabla_{\mathbf{w}}r_n=-\bm\phi_n\), hence \(\nabla_{\mathbf{w}}E_D\)\({}=\sum_n r_n\nabla_{\mathbf{w}}r_n\)\({}=-\sum_n r_n\bm\phi_n\), and \(\nabla\ln p\)\({}=-\beta\nabla E_D\)\({}=\beta\sum_n r_n\bm\phi_n\) (written as a row vector this is (3.13)).

Normal equations. Setting the gradient to zero: \(\sum_n t_n\bm\phi_n\)\({}=\big(\sum_n\bm\phi_n\bm\phi_n^\T\big)\mathbf{w}\) (this is (3.14) after transposing). Now \(\sum_n\bm\phi_n\bm\phi_n^\T=\bm\Phi^\T\bm\Phi\) (the rows of \(\bm\Phi\) are \(\bm\phi_n^\T\)) and \(\sum_nt_n\bm\phi_n=\bm\Phi^\T\mathbf{t}\), so \(\bm\Phi^\T\bm\Phi\mathbf{w}=\bm\Phi^\T\mathbf{t}\).

Global maximum. The Hessian of \(\ln p\) is \(-\beta\bm\Phi^\T\bm\Phi\), and \(\mathbf{v}^\T\bm\Phi^\T\bm\Phi\mathbf{v}\)\({}=\|\bm\Phi\mathbf{v}\|^2\)\({}\ge0\), so \(\ln p\) is concave and every stationary point is a global maximiser. The Hessian is negative definite (maximiser unique) exactly when \(\bm\Phi\mathbf{v}=0\Rightarrow\mathbf{v}=0\), i.e. the columns of \(\bm\Phi\) are independent; then \(\bm\Phi^\T\bm\Phi\) is invertible and (3.15) follows. If the columns are dependent, the normal equations still have solutions (they form an affine subspace); the pseudo-inverse then picks the one of minimum norm (Filling the gap 3.2). ∎

Proof

The maximisation over \((\mathbf{w},\beta)\) is done in two steps because \(\mathbf{w}_{\rm ML}\) does not depend on \(\beta\) (Result 3.1). For fixed \(\mathbf{w}=\mathbf{w}_{\rm ML}\), \(f(\beta)=\tfrac N2\ln\beta-\beta E_D\) has \(f'(\beta)=\tfrac N{2\beta}-E_D\) and \(f''(\beta)=-\tfrac N{2\beta^2}<0\), so \(f\) is concave with its unique maximum at \(\beta=N/(2E_D)\), which is (3.21).

Bias. Let \(\mathbf{t}=\bm\Phi\mathbf{w}^\star+\bm\epsilon\) with \(\bm\epsilon\sim\N(0,\beta_0^{-1}\mathbf{I}_N)\). The fitted vector is \(\bm\Phi\mathbf{w}_{\rm ML}=\mathbf{P}\mathbf{t}\), where \(\mathbf{P}\)\({}=\bm\Phi(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\) is the projection of Result 3.3. Since \(\mathbf{P}\bm\Phi=\bm\Phi\), the residual is \(\mathbf{t}-\mathbf{P}\mathbf{t}\)\({}=(\mathbf{I}-\mathbf{P})\bm\epsilon\). Using \(\mathbf{I}-\mathbf{P}\) symmetric and idempotent, \(\|(\mathbf{I}-\mathbf{P})\bm\epsilon\|^2\)\({}=\bm\epsilon^\T(\mathbf{I}-\mathbf{P})\bm\epsilon\), with expectation \(\beta_0^{-1}\Tr(\mathbf{I}-\mathbf{P})\)\({}=\beta_0^{-1}(N-\Tr\mathbf{P})\)\({}=\beta_0^{-1}(N-M)\) since \(\Tr\mathbf{P}\)\({}=\Tr[(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\bm\Phi]\)\({}=M\) (cyclic property of the trace, toolkit). Dividing by \(N\) gives the claim. Intuitively, fitting \(M\) weights uses up \(M\) degrees of freedom of the noise. ∎

3.1.2 Geometry of least squares

pp. 143–144 Think of \(\R^N\) with one axis per training point. The target vector \(\mathbf{t}=(t_1,\dots,t_N)^\T\) is a point of this space, and so is each basis function evaluated at the training inputs, \(\bm\varphi_j\)\({}=(\phi_j(\mathbf{x}_1),\dots,\phi_j(\mathbf{x}_N))^\T\), the \(j\)th column of \(\bm\Phi\) (the \(n\)th row is \(\bm\phi_n^\T\)). When \(M

Proof

Symmetry: \(\mathbf{P}^\T\)\({}=\bm\Phi[(\bm\Phi^\T\bm\Phi)^{-1}]^\T\bm\Phi^\T\)\({}=\mathbf{P}\) because \(\bm\Phi^\T\bm\Phi\) is symmetric. Idempotence: \(\mathbf{P}^2\)\({}=\bm\Phi(\bm\Phi^\T\bm\Phi)^{-1}[\bm\Phi^\T\bm\Phi](\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\)\({}=\mathbf{P}\). Range: \(\mathbf{P}\mathbf{v}\)\({}=\bm\Phi[(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\mathbf{v}]\) is a combination of the columns. Orthogonality: for any \(\mathbf{u}=\bm\Phi\mathbf{c}\in\mathcal S\), \(\mathbf{u}^\T(\mathbf{v}-\mathbf{P}\mathbf{v})\)\({}=\mathbf{c}^\T(\bm\Phi^\T\mathbf{v}-\bm\Phi^\T\bm\Phi(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\mathbf{v})\)\({}=0\). Finally \(\bm\Phi\mathbf{w}_{\rm ML}\)\({}=\bm\Phi(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\mathbf{t}\)\({}=\mathbf{P}\mathbf{t}\), and Pythagoras applies to the orthogonal decomposition \(\mathbf{t}\)\({}=\mathbf{P}\mathbf{t}\)\({}+(\mathbf{I}-\mathbf{P})\mathbf{t}\). That the projection is the closest point of \(\mathcal S\): for any \(\mathbf{u}\in\mathcal S\), \(\|\mathbf{t}-\mathbf{u}\|^2\)\({}=\|\mathbf{t}-\mathbf{y}\|^2+\|\mathbf{y}\)\({}-\mathbf{u}\|^2\)\({}\ge\|\mathbf{t}-\mathbf{y}\|^2\) (the cross term vanishes since \(\mathbf{t}\)\({}-\mathbf{y}\perp\mathcal S\ni\mathbf{y}\)\({}-\mathbf{u}\)). ∎

3.1.3 Sequential learning

pp. 143–144 For large or streaming data sets we prefer to update the weights one observation at a time. If the error is a sum over data points, \(E=\sum_nE_n\), stochastic (sequential) gradient descent updates after presenting pattern \(n\) by \(\mathbf{w}^{(\tau+1)}\)\({}=\mathbf{w}^{(\tau)}-\eta\nabla E_n\) (3.22), with learning rate \(\eta>0\).

Proof

(a) \(\nabla E_n\)\({}=-(t_n-\mathbf{w}^\T\bm\phi_n)\bm\phi_n\) gives (3.23). After the update the new residual on point \(n\) is \(t_n-\mathbf{w}'^\T\bm\phi_n\)\({}=(t_n-\mathbf{w}^\T\bm\phi_n)\)\({}-\eta(t_n-\mathbf{w}^\T\bm\phi_n)\|\bm\phi_n\|^2\)\({}=(1-\eta\|\bm\phi_n\|^2)(t_n-\mathbf{w}^\T\bm\phi_n)\). So \(\eta<1/\|\bm\phi_n\|^2\) makes the update shrink the error on the current point without overshooting.

(b) \(\nabla E_D\)\({}=\bm\Phi^\T\bm\Phi\mathbf{w}\)\({}-\bm\Phi^\T\mathbf{t}\)\({}=\bm\Phi^\T\bm\Phi(\mathbf{w}-\mathbf{w}_{\rm ML})\), so \(\mathbf{w}^{(\tau+1)}-\mathbf{w}_{\rm ML}\)\({}=(\mathbf{I}-\eta\bm\Phi^\T\bm\Phi)(\mathbf{w}^{(\tau)}-\mathbf{w}_{\rm ML})\). In the eigenbasis of \(\bm\Phi^\T\bm\Phi\) (eigenvalues \(d_i\ge0\)) each component is multiplied by \(1-\eta d_i\) per step. All components vanish iff \(|1-\eta d_i|<1\) for all \(i\), i.e. \(0<\eta<2/d_i\), i.e. \(\eta<2/d_{\max}\). The slowest component has factor \(1-\eta d_{\min}\); with \(\eta\) near \(2/d_{\max}\) the contraction per step is only about \(1-2d_{\min}/d_{\max}\).

(c) At \(\mathbf{w}_{\rm ML}\) the batch gradient vanishes but the individual \(\nabla E_n=-r_n\bm\phi_n\) do not (the residuals \(r_n\) are non-zero), so each step moves \(\mathbf{w}\) away by an amount proportional to \(\eta\); the weights therefore fluctuate with a spread that shrinks with \(\eta\). Decreasing step sizes with the stated conditions are the classical stochastic-approximation condition (Robbins and Monro, 1951; Section 2.3.5 of the book). ∎

The exact Bayesian counterpart of LMS, in which a matrix gain that shrinks as data accumulate replaces \(\eta\) and the batch answer is reproduced exactly after one pass, is the recursive update of Result 3.12.

3.1.4 Regularized least squares

pp. 144–147 To control over-fitting, add a penalty to the data error: minimise \(E_D(\mathbf{w})+\lambda E_W(\mathbf{w})\) (3.24), where \(\lambda\ge0\) sets the relative weight. The simplest penalty is \(E_W=\tfrac12\mathbf{w}^\T\mathbf{w}\) (3.25), known as weight decay in machine learning (because in sequential learning it pulls weights towards zero unless the data resist) and as ridge regression in statistics (Hoerl and Kennard, 1970). The total error is

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

Erratum. The book's (3.26) names the sum of squares \(E\) on its left-hand side; it is \(E_D\).

Proof

The gradient of (3.27) is \(-\bm\Phi^\T(\mathbf{t}-\bm\Phi\mathbf{w})\)\({}+\lambda\mathbf{w}\) and the Hessian is \(\bm\Phi^\T\bm\Phi+\lambda\mathbf{I}\), which is positive definite for \(\lambda>0\) (its eigenvalues are \(d_i+\lambda\ge\lambda\)). So the function is strictly convex, any stationary point is the unique global minimiser, and \((\bm\Phi^\T\bm\Phi+\lambda\mathbf{I})\mathbf{w}\)\({}=\bm\Phi^\T\mathbf{t}\) gives (3.28). In the eigenbasis, \(\mathbf{v}_i^\T\mathbf{w}\)\({}=\frac1{d_i+\lambda}\mathbf{v}_i^\T\bm\Phi^\T\mathbf{t}\) while the least-squares component is \(\frac1{d_i}\mathbf{v}_i^\T\bm\Phi^\T\mathbf{t}\) (when \(d_i>0\)); the ratio is \(d_i/(d_i+\lambda)\). ∎

Weight decay in a gradient step: \(\mathbf{w}\leftarrow\mathbf{w}\)\({}-\eta[\lambda\mathbf{w}-\bm\Phi^\T(\mathbf{t}-\bm\Phi\mathbf{w})]\)\({}=(1-\eta\lambda)\mathbf{w}\)\({}+\eta\bm\Phi^\T(\mathbf{t}-\bm\Phi\mathbf{w})\): the weights are multiplied by \(1-\eta\lambda<1\) before the data term acts. A further reading of the penalty, as the effect of adding noise to the inputs, is Exercise 3.4.

A more general penalty uses \(\tfrac\lambda2\sum_j|w_j|^q\) (3.29), with \(q=2\) giving (3.27). Figure 3.3 shows the contours of the penalty. The case \(q=1\) is the lasso (Tibshirani, 1996).

Proof

Penalised implies constrained. Write \(R(\mathbf{w})=\sum_j|w_j|^q\). Suppose some \(\mathbf{w}\) with \(R(\mathbf{w})\le\eta(\lambda)\) had \(E_D(\mathbf{w})\({}+\tfrac\lambda2R(\mathbf{w})\({}+\tfrac\lambda2\eta(\lambda)\)\({}=E_D(\mathbf{w}^\star)\)\({}+\tfrac\lambda2R(\mathbf{w}^\star)\), contradicting that \(\mathbf{w}^\star\) minimises (3.29). (No convexity is needed in this direction.)

Lagrangian. The Lagrangian of "minimise \(E_D\) s.t. \(\tfrac12(R(\mathbf{w})-\eta)\le0\)" is \(E_D(\mathbf{w})\)\({}+\tfrac\lambda2(R(\mathbf{w})-\eta)\), which as a function of \(\mathbf{w}\) is (3.29) up to the constant \(-\lambda\eta/2\). By the KKT conditions (toolkit), at a solution with an active constraint \(\lambda\ge0\) and \(R(\mathbf{w}^\star)=\eta\), which is the formula for \(\eta(\lambda)\). For convex \(E_D\) and \(R\) (\(q\ge1\)) with \(\eta>0\) (so \(\mathbf{w}=0\) is strictly feasible, Slater's condition), the KKT conditions are also sufficient and a multiplier \(\lambda\) exists, which gives the converse.

Monotonicity. Let \(\lambda_1<\lambda_2\) with minimisers \(\mathbf{w}_1,\mathbf{w}_2\), \(R_i=R(\mathbf{w}_i)\), \(D_i=E_D(\mathbf{w}_i)\). Optimality of each for its own objective gives \(D_1+\tfrac{\lambda_1}2R_1\)\({}\le D_2+\tfrac{\lambda_1}2R_2\) and \(D_2+\tfrac{\lambda_2}2R_2\)\({}\le D_1+\tfrac{\lambda_2}2R_1\). Adding, \(\lambda_1R_1+\lambda_2R_2\)\({}\le\lambda_1R_2+\lambda_2R_1\), i.e. \((\lambda_2-\lambda_1)(R_1-R_2)\ge0\), so \(R_1\ge R_2\). The multiplier is the price of the budget: by the envelope theorem \(\lambda/2=-\d E_D^\star/\d\eta\). ∎

Regularisation lets complex models be trained on small data sets without severe over-fitting, but the problem of choosing model complexity is moved to choosing \(\lambda\) (or \(q\)). Sections 3.2 to 3.5 are about exactly this. The book then fixes the quadratic regulariser for the rest of the chapter.

3.1.5 Multiple outputs

pp. 146–147 For \(K>1\) targets \(\mathbf{t}\in\R^K\) one can share the same basis functions: \(\mathbf{y}(\mathbf{x},\mathbf{W})\)\({}=\mathbf{W}^\T\bm\phi(\mathbf{x})\) (3.31), with \(\mathbf{W}\) an \(M\times K\) matrix, and an isotropic Gaussian noise model \(p(\mathbf{t}\given\mathbf{x},\mathbf{W},\beta)\)\({}=\N(\mathbf{t}\given\mathbf{W}^\T\bm\phi(\mathbf{x}),\beta^{-1}\mathbf{I})\) (3.32). Stack the targets in an \(N\times K\) matrix \(\mathbf{T}\) whose \(n\)th row is \(\mathbf{t}_n^\T\). Then

\[ \ln p(\mathbf{T}\given\mathbf{X},\mathbf{W},\beta)=\frac{NK}2\ln\frac\beta{2\pi}-\frac\beta2\sum_{n=1}^{N}\big\|\mathbf{t}_n-\mathbf{W}^\T\bm\phi(\mathbf{x}_n)\big\|^2. \tag{3.33} \]

Proof

The squared norm splits over output components: \(\sum_n\|\mathbf{t}_n-\mathbf{W}^\T\bm\phi_n\|^2\)\({}=\sum_{k=1}^{K}\sum_n(t_{nk}-\mathbf{w}_k^\T\bm\phi_n)^2\)\({}=\sum_k\|\mathbf{t}_k-\bm\Phi\mathbf{w}_k\|^2\). So (3.33) is the sum of \(K\) separate single-output log likelihoods in \(\mathbf{w}_1,\dots,\mathbf{w}_K\) which share no parameters; each is maximised by Result 3.1, giving (3.35); stacking the columns gives (3.34). (In matrix form: \(\|\mathbf{T}-\bm\Phi\mathbf{W}\|_F^2\) has gradient \(-2\bm\Phi^\T(\mathbf{T}-\bm\Phi\mathbf{W})\) by matrix calculus, zero at (3.34).) The maximum over \(\beta\) is \(1/\beta_{\rm ML}\)\({}=\frac1{NK}\sum_n\|\mathbf{t}_n\)\({}-\mathbf{W}_{\rm ML}^\T\bm\phi_n\|^2\), by the same argument as Result 3.2 with \(NK\) in place of \(N\). ∎

From now on the book returns to a single target \(t\).

3.2 The bias-variance decomposition

pp. 147–152 Maximum likelihood over-fits flexible models, and limiting the number of basis functions limits flexibility; a regulariser seems to dissolve the dilemma but moves it to the choice of \(\lambda\). Minimising the regularised error over \(\lambda\) as well as \(\mathbf{w}\) is not an option, since it returns \(\lambda=0\). Before the Bayesian answer, the book gives the frequentist view of model complexity: the bias-variance trade-off.

Let \(h(\mathbf{x})\)\({}=\E[t\given\mathbf{x}]\)\({}=\int t\,p(t\given\mathbf{x})\,\d t\) (3.36) be the optimal prediction under squared loss (Section 1.5.5). For any prediction function \(y(\mathbf{x})\) the expected squared loss is

\[ \E[L]=\int\{y(\mathbf{x})-h(\mathbf{x})\}^2p(\mathbf{x})\,\d\mathbf{x}+\iint\{h(\mathbf{x})-t\}^2p(\mathbf{x},t)\,\d\mathbf{x}\,\d t. \tag{3.37} \]

Erratum. The book prints a single integral sign on the second term of (3.37), and likewise in (3.44); both are double integrals (over \(\mathbf{x}\) and \(t\)), as written here.

The first term of (3.37) depends on the particular training set \(\mathcal D\) through \(y(\mathbf{x};\mathcal D)\). Our uncertainty about \(\mathcal D\) is handled by averaging over data sets drawn from \(p(\mathcal D)\), written \(\E_{\mathcal D}[\cdot]\).

Proof

Fix \(\mathbf{x}\) and abbreviate \(y=y(\mathbf{x};\mathcal D)\), \(\bar y=\E_{\mathcal D}[y]\), \(h=h(\mathbf{x})\). Add and subtract \(\bar y\):

\[ \{y-h\}^2=\{y-\bar y\}^2+\{\bar y-h\}^2+2\{y-\bar y\}\{\bar y-h\}. \tag{3.39} \]

Take \(\E_{\mathcal D}\). The term \(\{\bar y-h\}^2\) is a constant with respect to \(\mathcal D\). The cross term is \(2\{\bar y-h\}\,\E_{\mathcal D}[y-\bar y]\)\({}=2\{\bar y-h\}(\bar y-\bar y)\)\({}=0\). This proves (3.40). Now average (3.37) over \(\mathcal D\) (the noise term does not involve \(\mathcal D\) as \(h\) and \(t\) refer to a fresh test point) and substitute (3.40) in its first term: integrating the two terms of (3.40) against \(p(\mathbf{x})\) gives (3.42) and (3.43), and the second term of (3.37) is (3.44). ∎

The trade-off: very flexible models have low bias and high variance, rigid models high bias and low variance, and the best predictor balances the two (the noise is a floor). Intuition from the polynomial example: a straight line fitted to a sine curve is always about the same wrong line (small variance, large bias), whereas a ninth-order polynomial follows each noisy data set, so the average over data sets is close to the truth but individual fits vary enormously.

The decomposition for ridge regression, in closed form

For the ridge estimator the bias and variance can be computed exactly, which explains the shapes in Figure 3.6.

Proof

Let \(\mathbf{B}\)\({}=(\bm\Phi^\T\bm\Phi+\lambda\mathbf{I})^{-1}\). Then \(\hat{\mathbf{w}}_\lambda\)\({}=\mathbf{B}\bm\Phi^\T\bm\Phi\mathbf{w}^\star\)\({}+\mathbf{B}\bm\Phi^\T\bm\epsilon\). The mean is \(\mathbf{B}\bm\Phi^\T\bm\Phi\mathbf{w}^\star\)\({}=\mathbf{B}(\bm\Phi^\T\bm\Phi+\lambda\mathbf{I}-\lambda\mathbf{I})\mathbf{w}^\star\)\({}=\mathbf{w}^\star\)\({}-\lambda\mathbf{B}\mathbf{w}^\star\), so the bias vector is \(-\lambda\mathbf{B}\mathbf{w}^\star\), with \(\mathbf{v}_i\)-component \(-\lambda w^\star_i/(d_i+\lambda)\). The covariance is \(\sigma^2\mathbf{B}\bm\Phi^\T\bm\Phi\mathbf{B}\)\({}=\sigma^2\mathbf{V}\operatorname{diag}\big(d_i/(d_i+\lambda)^2\big)\mathbf{V}^\T\), because \(\mathbf{B}\) and \(\bm\Phi^\T\bm\Phi\) share eigenvectors; this gives the variance along \(\mathbf{v}_i\). For the prediction error, \(\|\bm\Phi\mathbf{z}\|^2\)\({}=\sum_id_i(\mathbf{v}_i^\T\mathbf{z})^2\) and \(\E(\mathbf{v}_i^\T(\hat{\mathbf{w}}-\mathbf{w}^\star))^2\)\({}=\text{bias}_i^2+\text{variance}_i\). Monotonicity: \(\lambda/(\lambda+d_i)\) increases and \(d_i/(\lambda+d_i)^2\) decreases in \(\lambda\). ∎

The link with the Bayesian view

The frequentist decomposition needs many data sets, which in practice are not available (if they were, one would pool them into one big set and over-fit less). Bayesian inference replaces the average over hypothetical data sets by an average over parameters given the single observed set, and this section's trade-off survives in three precise ways.

Proof

(b) Average Result 3.9 over \(w_i^\star\) with \(\E w_i^{\star2}=1/\alpha\) (the components of \(\mathbf{w}^\star\) in any orthonormal basis are independent \(\N(0,\alpha^{-1})\)). Let \(f(\lambda)\)\({}=(\lambda^2/\alpha+\sigma^2d)/(\lambda+d)^2\). Then \[ f'(\lambda)=\frac{2(\lambda/\alpha)(\lambda+d)-2(\lambda^2/\alpha+\sigma^2d)}{(\lambda+d)^3}=\frac{2d}{\alpha}\,\frac{\lambda-\sigma^2\alpha}{(\lambda+d)^3}, \] which is negative for \(\lambda<\sigma^2\alpha\) and positive after, so the minimum is at \(\lambda=\sigma^2\alpha\) whatever \(d\) is. There \(f\)\({}=(\sigma^4\alpha+\sigma^2d)/(\sigma^2\alpha+d)^2\)\({}=\sigma^2/(\sigma^2\alpha+d)\)\({}=1/(\alpha+\beta d)\) with \(\sigma^2=1/\beta\). Because the optimal \(\lambda\) is the same for each \(i\), it also minimises the weighted sum in the prediction error. Equivalently, the sum over \(i\) is \(\Tr\mathbf{S}_N=\E\|\mathbf{w}-\mathbf{m}_N\|^2\), the Bayes risk of the posterior mean, which is minimal among all estimators.

(c) Given the data, \(t=\mathbf{w}^\T\bm\phi+\epsilon\) with \(\mathbf{w}\sim\N(\mathbf{m}_N,\mathbf{S}_N)\) and \(\epsilon\) independent. The law of total variance gives \(\var[t]\)\({}=\E_{\mathbf{w}}[\var(t\given\mathbf{w})]\)\({}+\var_{\mathbf{w}}[\E(t\given\mathbf{w})]\)\({}=\beta^{-1}+\bm\phi^\T\mathbf{S}_N\bm\phi\), which is (3.59). ∎

So the Bayesian treatment does not escape the trade-off, it resolves it: the prior says how large the true weights are expected to be, and that determines how much bias to accept in exchange for lower variance. When the prior is not known, Section 3.5 learns \(\alpha\) and \(\beta\) from the data, and so \(\lambda=\alpha/\beta\), without a validation set.

3.3 Bayesian linear regression

pp. 152–161 Maximum likelihood always prefers the more complex model, and hold-out data are wasteful. A Bayesian treatment places a prior on the weights, integrates over them, and thereby avoids over-fitting and gives a way to choose complexity from the training data alone. For simplicity there is one target; the extension to several follows Section 3.1.5.

3.3.1 Parameter distribution

pp. 152–156 Treat the noise precision \(\beta\) as a known constant for now. The likelihood (3.10) is the exponential of a quadratic function of \(\mathbf{w}\), so the conjugate prior is Gaussian: \(p(\mathbf{w})\)\({}=\N(\mathbf{w}\given\mathbf{m}_0,\mathbf{S}_0)\) (3.48).

Proof (completing the square)

By Bayes' theorem \(p(\mathbf{w}\given\mathbf{t})\)\({}\propto p(\mathbf{t}\given\mathbf{w})p(\mathbf{w})\). Keep only terms involving \(\mathbf{w}\):

\[ \ln p(\mathbf{w}\given\mathbf{t})=-\frac\beta2(\mathbf{t}-\bm\Phi\mathbf{w})^\T(\mathbf{t}-\bm\Phi\mathbf{w})-\frac12(\mathbf{w}-\mathbf{m}_0)^\T\mathbf{S}_0^{-1}(\mathbf{w}-\mathbf{m}_0)+\const. \]

Expand both quadratics and collect by powers of \(\mathbf{w}\):

\[ \ln p(\mathbf{w}\given\mathbf{t})=-\frac12\Big[\mathbf{w}^\T(\mathbf{S}_0^{-1}+\beta\bm\Phi^\T\bm\Phi)\mathbf{w}-2\mathbf{w}^\T(\mathbf{S}_0^{-1}\mathbf{m}_0+\beta\bm\Phi^\T\mathbf{t})\Big]+\const. \]

With \(\mathbf{S}_N^{-1}:\)\({}=\mathbf{S}_0^{-1}+\beta\bm\Phi^\T\bm\Phi\) (symmetric positive definite) and \(\mathbf{m}_N:\)\({}=\mathbf{S}_N(\mathbf{S}_0^{-1}\mathbf{m}_0+\beta\bm\Phi^\T\mathbf{t})\), the identity \((\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)\)\({}=\mathbf{w}^\T\mathbf{S}_N^{-1}\mathbf{w}\)\({}-2\mathbf{w}^\T\mathbf{S}_N^{-1}\mathbf{m}_N\)\({}+\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N\) shows the bracket equals \((\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)\) minus a term free of \(\mathbf{w}\). Hence \(\ln p\)\({}=-\frac12(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)\)\({}+\const\), the kernel of \(\N(\mathbf{m}_N,\mathbf{S}_N)\); the constant is fixed by normalisation (toolkit). The same result follows from the linear-Gaussian formula (toolkit, book (2.116)) with prior \(\N(\mathbf{m}_0,\mathbf{S}_0)\), matrix \(\bm\Phi\) and noise precision \(\beta\).

For (3.55): set \(\mathbf{m}_0\)\({}=0,\ \mathbf{S}_0\)\({}=\alpha^{-1}\mathbf{I}\) in the display above. Maximising it is minimising \(\frac12\sum_n\{\cdot\}^2\)\({}+\frac{\alpha/\beta}2\mathbf{w}^\T\mathbf{w}\), i.e. (3.27) with \(\lambda=\alpha/\beta\); by Result 3.5, \(\mathbf{m}_N\)\({}=(\lambda\mathbf{I}+\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\mathbf{t}\), which is also \(\beta\mathbf{S}_N\bm\Phi^\T\mathbf{t}\) since \(\beta\mathbf{S}_N\)\({}=(\frac\alpha\beta\mathbf{I}+\bm\Phi^\T\bm\Phi)^{-1}\). The limit \(\alpha\to0\) gives (3.15). ∎

A reading worth keeping: the Bayesian posterior mean is ridge regression with a specific, interpretable \(\lambda=\alpha/\beta\) (the ratio of the noise variance to the prior variance), and in addition the posterior supplies a covariance \(\mathbf{S}_N\), which is what the predictive distribution below uses.

Proof

Treat \(\N(\mathbf{w}\given\mathbf{m}_N,\mathbf{S}_N)\) as the prior and multiply by the likelihood of the new point, \(\propto\exp\{-\frac\beta2(t_{N+1}-\mathbf{w}^\T\bm\phi)^2\}\). The exponent, up to the factor \(-\frac12\), is \((\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)\)\({}+\beta(t_{N+1}-\mathbf{w}^\T\bm\phi)^2\)\({}=\mathbf{w}^\T(\mathbf{S}_N^{-1}+\beta\bm\phi\bm\phi^\T)\mathbf{w}\)\({}-2\mathbf{w}^\T(\mathbf{S}_N^{-1}\mathbf{m}_N+\beta\bm\phi\,t_{N+1})\)\({}+\const\), which has the same form as in the proof of Result 3.11, giving the first pair of formulas (Exercise 3.8). The linear-Gaussian route (Exercise 3.9) reads off the same formulas from (2.116) with \(\bm\Lambda=\mathbf{S}_N^{-1}\), \(\mathbf{A}=\bm\phi^\T\), \(\mathbf{L}=\beta\).

Gain form. Sherman–Morrison (Woodbury with rank one): \((\mathbf{S}_N^{-1}+\beta\bm\phi\bm\phi^\T)^{-1}\)\({}=\mathbf{S}_N\)\({}-\dfrac{\beta\mathbf{S}_N\bm\phi\bm\phi^\T\mathbf{S}_N}{1+\beta\bm\phi^\T\mathbf{S}_N\bm\phi}\)\({}=\mathbf{S}_N-\mathbf{k}\bm\phi^\T\mathbf{S}_N\), with \(\mathbf{k}\) as above (multiply numerator and denominator by \(\beta^{-1}\)). Also \(\mathbf{S}_{N+1}\bm\phi\)\({}=\mathbf{S}_N\bm\phi\)\({}-\dfrac{\beta\mathbf{S}_N\bm\phi(\bm\phi^\T\mathbf{S}_N\bm\phi)}{1+\beta\bm\phi^\T\mathbf{S}_N\bm\phi}\)\({}=\dfrac{\mathbf{S}_N\bm\phi}{1+\beta\bm\phi^\T\mathbf{S}_N\bm\phi}\)\({}=\beta^{-1}\mathbf{k}\), so \(\mathbf{k}=\beta\mathbf{S}_{N+1}\bm\phi\). Using \(\mathbf{S}_N^{-1}\)\({}=\mathbf{S}_{N+1}^{-1}-\beta\bm\phi\bm\phi^\T\) in the mean, \[ \mathbf{m}_{N+1}=\mathbf{S}_{N+1}\mathbf{S}_{N+1}^{-1}\mathbf{m}_N-\beta\mathbf{S}_{N+1}\bm\phi\bm\phi^\T\mathbf{m}_N+\beta\mathbf{S}_{N+1}\bm\phi\,t_{N+1}=\mathbf{m}_N+\beta\mathbf{S}_{N+1}\bm\phi\,(t_{N+1}-\bm\phi^\T\mathbf{m}_N). \]

Order independence: after any ordering the precision is \(\mathbf{S}_0^{-1}\)\({}+\beta\sum_n\bm\phi_n\bm\phi_n^\T\) and the linear term is \(\mathbf{S}_0^{-1}\mathbf{m}_0\)\({}+\beta\sum_nt_n\bm\phi_n\), both symmetric functions of the data. ∎

A worked example: fitting a straight line

To see sequential learning, take \(y(x,\mathbf{w})=w_0+w_1x\), so the two-dimensional posterior can be drawn. Data come from \(f(x,\mathbf{a})=a_0+a_1x\) with \(a_0=-0.3\), \(a_1=0.5\): inputs \(x_n\sim\mathcal U(-1,1)\), noise standard deviation \(0.2\), so the true \(\beta=(1/0.2)^2=25\), which is also used in the model, and \(\alpha=2.0\). Figure 3.7 shows, row by row: the prior and six lines sampled from it; after one data point, the likelihood of that point (a ridge in \(\mathbf{w}\)-space, the soft constraint that the line pass near it), the posterior (the product with the previous row) and six lines sampled from it, which now pass close to the point; after two points, a compact posterior (two points determine a line); after 20 points, a sharp posterior around the true values (white cross). The key lesson, which Result 3.12 makes exact: the posterior after the second point is the same as the one obtained by combining the original prior with the likelihood of both points at once.

3.3.2 Predictive distribution

pp. 156–158 One rarely wants \(\mathbf{w}\) itself, but predictions at new inputs. The predictive distribution averages the likelihood of a new target over the posterior,

\[ p(t\given\mathbf{t},\alpha,\beta)=\int p(t\given\mathbf{w},\beta)\,p(\mathbf{w}\given\mathbf{t},\alpha,\beta)\,\d\mathbf{w}. \tag{3.57} \]

Proof (two ways)

(i) Constructive. Under the posterior, \(\mathbf{w}\sim\N(\mathbf{m}_N,\mathbf{S}_N)\), and a new target is \(t=\mathbf{w}^\T\bm\phi+\epsilon\) with \(\epsilon\sim\N(0,\beta^{-1})\) independent of \(\mathbf{w}\) (given \(\mathbf{w}\), the new noise is independent of the old). Then \(t\) is a linear function of the jointly Gaussian pair \((\mathbf{w},\epsilon)\), hence Gaussian, with mean \(\bm\phi^\T\mathbf{m}_N\) and variance \(\bm\phi^\T\mathbf{S}_N\bm\phi+\beta^{-1}\) (variances of independent summands add). This is (3.58) and (3.59).

(ii) By the book's route. Write the integrand as \(\N(t\given\bm\phi^\T\mathbf{w},\beta^{-1})\N(\mathbf{w}\given\mathbf{m}_N,\mathbf{S}_N)\). This is a linear-Gaussian model, \(p(\mathbf{w})=\N(\bm\mu,\bm\Lambda^{-1})\), \(p(t\given\mathbf{w})\)\({}=\N(\mathbf{A}\mathbf{w}+\mathbf{b},\mathbf{L}^{-1})\) with \(\mathbf{A}=\bm\phi^\T\), \(\mathbf{b}=0\), \(\mathbf{L}^{-1}=\beta^{-1}\), and the marginal is \(\N(\mathbf{A}\bm\mu+\mathbf{b},\mathbf{L}^{-1}+\mathbf{A}\bm\Lambda^{-1}\mathbf{A}^\T)\) (book (2.115), toolkit), i.e. mean \(\bm\phi^\T\mathbf{m}_N\) and variance \(\beta^{-1}+\bm\phi^\T\mathbf{S}_N\bm\phi\). Since \(\mathbf{S}_N\) is positive definite, \(\bm\phi^\T\mathbf{S}_N\bm\phi\ge0\). ∎

Proof

By Result 3.12, \(\mathbf{S}_{N+1}\)\({}=\mathbf{S}_N\)\({}-\dfrac{\beta\mathbf{S}_N\bm\phi\bm\phi^\T\mathbf{S}_N}{1+\beta\bm\phi^\T\mathbf{S}_N\bm\phi}\); this is the identity (3.110), \((\mathbf{M}+\mathbf{v}\mathbf{v}^\T)^{-1}\)\({}=\mathbf{M}^{-1}\)\({}-\dfrac{(\mathbf{M}^{-1}\mathbf{v})(\mathbf{v}^\T\mathbf{M}^{-1})}{1+\mathbf{v}^\T\mathbf{M}^{-1}\mathbf{v}}\), with \(\mathbf{M}=\mathbf{S}_N^{-1}\) and \(\mathbf{v}=\sqrt\beta\bm\phi\). Sandwiching between \(\bm\phi(\mathbf{x})^\T\) and \(\bm\phi(\mathbf{x})\) (call it \(\bm\psi\) to avoid a clash): \(\bm\psi^\T\mathbf{S}_{N+1}\bm\psi\)\({}=\bm\psi^\T\mathbf{S}_N\bm\psi\)\({}-\dfrac{\beta(\bm\psi^\T\mathbf{S}_N\bm\phi)^2}{1+\beta\bm\phi^\T\mathbf{S}_N\bm\phi}\). The subtracted term is a non-negative number over a number \(\ge1\) (as \(\mathbf{S}_N\succ0\)), which is the inequality. For the limit: \(\mathbf{S}_N^{-1}\)\({}=\alpha\mathbf{I}\)\({}+\beta\bm\Phi^\T\bm\Phi\succeq(\alpha+\beta\lambda_{\min})\mathbf{I}\), hence \(\mathbf{S}_N\preceq(\alpha+\beta\lambda_{\min})^{-1}\mathbf{I}\) and \(\bm\phi^\T\mathbf{S}_N\bm\phi\)\({}\le\|\bm\phi\|^2/(\alpha+\beta\lambda_{\min})\). ∎

The plots of Figure 3.8 show only the pointwise predictive variance. The joint uncertainty across inputs is visible by drawing weight vectors from the posterior and plotting the functions \(y(x,\mathbf{w})\) (Figure 3.9); it is governed by the covariance (3.63) of Section 3.3.3.

3.3.3 Equivalent kernel

pp. 159–161 Substituting the posterior mean (3.53) into \(y(\mathbf{x},\mathbf{m}_N)\)\({}=\mathbf{m}_N^\T\bm\phi(\mathbf{x})\) shows that the prediction is a linear combination of the training targets:

\[ y(\mathbf{x},\mathbf{m}_N)=\mathbf{m}_N^\T\bm\phi(\mathbf{x})=\beta\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\Phi^\T\mathbf{t}=\sum_{n=1}^{N}\beta\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\phi(\mathbf{x}_n)\,t_n. \tag{3.60} \]

(The scalar \(\mathbf{m}_N^\T\bm\phi\) equals its transpose \(\bm\phi^\T\mathbf{m}_N\), and \(\bm\Phi^\T\mathbf{t}=\sum_n\bm\phi_nt_n\).) Writing

\[ y(\mathbf{x},\mathbf{m}_N)=\sum_{n=1}^{N}k(\mathbf{x},\mathbf{x}_n)\,t_n,\qquad k(\mathbf{x},\mathbf{x}')=\beta\,\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\phi(\mathbf{x}'), \tag{3.61, 3.62} \]

\(k\) is the equivalent kernel or smoother matrix, and predictors of this form are linear smoothers. Note that \(k\) depends on the training inputs through \(\mathbf{S}_N\) (but not on the targets).

Proof

(a) \(y(\mathbf{x})=\bm\phi(\mathbf{x})^\T\mathbf{w}\) with \(\cov\mathbf{w}=\mathbf{S}_N\), so \(\cov[\bm\phi(\mathbf{x})^\T\mathbf{w},\mathbf{w}^\T\bm\phi(\mathbf{x}')]\)\({}=\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\phi(\mathbf{x}')\). (b) \(\mathbf{S}_N\) is symmetric positive definite, so it has a symmetric square root \(\mathbf{S}_N^{1/2}\) (from its eigendecomposition) and \(\beta\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\phi(\mathbf{z})\)\({}=(\beta^{1/2}\mathbf{S}_N^{1/2}\bm\phi(\mathbf{x}))^\T(\beta^{1/2}\mathbf{S}_N^{1/2}\bm\phi(\mathbf{z}))\); a Gram matrix \(\sum_{ab}c_ac_bk(\mathbf{x}_a,\mathbf{x}_b)\)\({}=\|\sum_ac_a\bm\psi(\mathbf{x}_a)\|^2\)\({}\ge0\).

(c) Direct proof. For \(\alpha=0\), \(\mathbf{S}_N=(\beta\bm\Phi^\T\bm\Phi)^{-1}\), so \(\sum_nk(\mathbf{x},\mathbf{x}_n)\)\({}=\bm\phi(\mathbf{x})^\T(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\mathbf{1}\), where \(\mathbf{1}=(1,\dots,1)^\T\). Since one basis function is constant, say \(\phi_0\equiv1\), the vector \(\mathbf{1}\) is the first column of \(\bm\Phi\), so \(\bm\Phi^\T\mathbf{1}\)\({}=\bm\Phi^\T\bm\Phi\mathbf{e}_0\) and the sum equals \(\bm\phi(\mathbf{x})^\T\mathbf{e}_0\)\({}=\phi_0(\mathbf{x})\)\({}=1\). (This is the book's informal argument that least squares fits constant targets exactly.) For \(\alpha>0\), \(\sum_nk\)\({}=\bm\phi^\T(\frac\alpha\beta\mathbf{I}+\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\bm\Phi\mathbf{e}_0\)\({}=\bm\phi^\T\mathbf{e}_0\)\({}-\frac\alpha\beta\bm\phi^\T(\frac\alpha\beta\mathbf{I}+\bm\Phi^\T\bm\Phi)^{-1}\mathbf{e}_0\), which differs from 1 by a term of order \(\alpha/\beta\), vanishing as \(\alpha\to0\) or as \(N\to\infty\). A fully worked version using orthonormalised basis functions is Exercise 3.14.

(d) On the training inputs \(\mathbf{K}\)\({}=\beta\bm\Phi\mathbf{S}_N\bm\Phi^\T\)\({}=\bm\Phi(\frac\alpha\beta\mathbf{I}+\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\); with the SVD of Gap 3.2 this is \(\mathbf{U}\operatorname{diag}\big(\frac{s_i^2}{s_i^2+\alpha/\beta}\big)\mathbf{U}^\T\), whose eigenvalues are in \([0,1)\). Its trace, by cyclicity, is \(\Tr[\beta\mathbf{S}_N\bm\Phi^\T\bm\Phi]\)\({}=\Tr[\mathbf{S}_N(\mathbf{S}_N^{-1}-\alpha\mathbf{I})]\)\({}=M-\alpha\Tr\mathbf{S}_N\)\({}=\sum_i\big(1-\frac\alpha{\lambda_i+\alpha}\big)\)\({}=\sum_i\frac{\lambda_i}{\lambda_i+\alpha}\) where \(\lambda_i=\beta d_i\) are the eigenvalues of \(\beta\bm\Phi^\T\bm\Phi\). ∎

This suggests an alternative way to do regression: instead of choosing basis functions (which implicitly fix a kernel), choose a localised kernel directly and predict with \(\sum_nk(\mathbf{x},\mathbf{x}_n)t_n\). That is the programme of kernel methods (Chapter 6) and Gaussian processes (Section 6.4): by (a) and (b), the posterior distribution of the function \(y(\mathbf{x})=\mathbf{w}^\T\bm\phi(\mathbf{x})\) in Bayesian linear regression is itself a Gaussian process, with mean function \(y(\mathbf{x},\mathbf{m}_N)\) and covariance function \(\beta^{-1}k\). Erratum. The book says "effective kernel" once at the start of the paragraph on p. 160; the term is "equivalent kernel".

3.4 Bayesian model comparison

pp. 161–165 Over-fitting is a symptom of maximum likelihood, not of Bayesian inference: if parameters are integrated out instead of point-estimated, models can be compared on the training data alone, with no validation set, and several complexity parameters can be tuned at once (the relevance vector machine of Chapter 7 has one per data point). Suppose we have \(L\) candidate models \(\{\mathcal M_i\}\), each a probability distribution over the observed data \(\mathcal D\) (for curve fitting, over the targets \(\mathbf{t}\) given known inputs). We are uncertain which generated the data, expressed by a prior \(p(\mathcal M_i)\). Given \(\mathcal D\), Bayes' theorem gives

\[ p(\mathcal M_i\given\mathcal D)\propto p(\mathcal M_i)\,p(\mathcal D\given\mathcal M_i). \tag{3.66} \]

With equal priors the interesting term is the model evidence or marginal likelihood \(p(\mathcal D\given\mathcal M_i)\), the likelihood as a function of the model with the parameters marginalised out. The ratio \(p(\mathcal D\given\mathcal M_i)/p(\mathcal D\given\mathcal M_j)\) is a Bayes factor (Kass and Raftery, 1995). Predictions follow from the sum and product rules as a model average,

\[ p(t\given\mathbf{x},\mathcal D)=\sum_{i=1}^{L}p(t\given\mathbf{x},\mathcal M_i,\mathcal D)\,p(\mathcal M_i\given\mathcal D), \tag{3.67} \]

a mixture of the individual predictive distributions. If two models are equally probable and one predicts a narrow density near \(t=a\), the other near \(t=b\), the average is bimodal with modes at \(a\) and \(b\), not a single mode at \((a+b)/2\). Model selection is the cheap approximation that keeps only the most probable model. For a model with parameters \(\mathbf{w}\),

\[ p(\mathcal D\given\mathcal M_i)=\int p(\mathcal D\given\mathbf{w},\mathcal M_i)\,p(\mathbf{w}\given\mathcal M_i)\,\d\mathbf{w}, \tag{3.68} \]

which is the probability of generating \(\mathcal D\) from a model whose parameters are drawn at random from the prior; it is also exactly the normaliser in Bayes' theorem for the parameters,

\[ p(\mathbf{w}\given\mathcal D,\mathcal M_i)=\frac{p(\mathcal D\given\mathbf{w},\mathcal M_i)\,p(\mathbf{w}\given\mathcal M_i)}{p(\mathcal D\given\mathcal M_i)}. \tag{3.69} \]

Why the evidence penalises complexity

Consider one parameter \(w\) and drop \(\mathcal M_i\) from the notation. If the posterior is sharply peaked at \(w_{\rm MAP}\) with width \(\Delta w_{\rm posterior}\) and the prior is flat of width \(\Delta w_{\rm prior}\), so \(p(w)=1/\Delta w_{\rm prior}\), then approximating the integral by the height times the width of the peak,

\[ p(\mathcal D)=\int p(\mathcal D\given w)p(w)\,\d w\simeq p(\mathcal D\given w_{\rm MAP})\,\frac{\Delta w_{\rm posterior}}{\Delta w_{\rm prior}}, \tag{3.70} \] \[ \ln p(\mathcal D)\simeq\ln p(\mathcal D\given w_{\rm MAP})+\ln\frac{\Delta w_{\rm posterior}}{\Delta w_{\rm prior}}. \tag{3.71} \]

The first term is the fit to the data at the best parameter; the second is negative (\(\Delta w_{\rm posterior}<\Delta w_{\rm prior}\)) and is larger in magnitude the more finely the parameter must be tuned. With \(M\) parameters, each with the same ratio,

\[ \ln p(\mathcal D)\simeq\ln p(\mathcal D\given\mathbf{w}_{\rm MAP})+M\ln\frac{\Delta w_{\rm posterior}}{\Delta w_{\rm prior}}, \tag{3.72} \]

so the penalty grows linearly with \(M\). Adding parameters raises the first term (better fit) and lowers the second; the maximum evidence is the best compromise.

Proof

By Bayes' theorem (3.69) rearranged at \(\mathbf{w}=\mathbf{m}_N\), \(p(\mathbf{t}\given\alpha,\beta)\)\({}=p(\mathbf{t}\given\mathbf{m}_N,\beta)\,p(\mathbf{m}_N\given\alpha)/p(\mathbf{m}_N\given\mathbf{t})\). The posterior density at its mean is \(\N(\mathbf{m}_N\given\mathbf{m}_N,\mathbf{S}_N)\)\({}=(2\pi)^{-M/2}|\mathbf{S}_N|^{-1/2}\)\({}=(2\pi)^{-M/2}|\mathbf{A}|^{1/2}\) with \(\mathbf{A}=\mathbf{S}_N^{-1}\), and the prior density is \((\alpha/2\pi)^{M/2}e^{-\alpha\mathbf{m}_N^\T\mathbf{m}_N/2}\). Hence \(\ln p(\mathbf{t}\given\alpha,\beta)\)\({}=\ln p(\mathbf{t}\given\mathbf{m}_N,\beta)\)\({}+\frac M2\ln\alpha\)\({}-\frac\alpha2\mathbf{m}_N^\T\mathbf{m}_N\)\({}-\frac12\ln|\mathbf{A}|\). Since \(\mathbf{A}\)\({}=\alpha\mathbf{I}+\beta\bm\Phi^\T\bm\Phi\) has eigenvalues \(\alpha+\lambda_i\), \(\frac M2\ln\alpha-\frac12\ln|\mathbf{A}|\)\({}=-\frac12\sum_i\ln\frac{\alpha+\lambda_i}\alpha\). The posterior variance along \(\mathbf{u}_i\) is \(1/(\alpha+\lambda_i)\) and the prior variance \(1/\alpha\), whose ratio is \((1+\lambda_i/\alpha)^{-1}\). Substituting \(\ln p(\mathbf{t}\given\mathbf{m}_N,\beta)\)\({}=\frac N2\ln\frac\beta{2\pi}\)\({}-\beta E_D(\mathbf{m}_N)\) reproduces (3.86) of Result 3.18, which verifies the claim. The chain-rule identity is the product rule applied to \(p(t_1,\dots,t_N\given\mathcal M)\), each factor being a predictive distribution of the kind (3.57) (for the linear model, the Gaussian (3.58) with the data so far). ∎

Why intermediate complexity can win (Figure 3.13)

Put a one-dimensional axis for "the space of possible data sets" and run three models of increasing complexity generatively: pick parameters from the prior, then sample data. A simple model (a first-order polynomial) produces similar data sets every time, so its distribution \(p(\mathcal D\given\mathcal M_1)\) is concentrated in a small region. A complex model (ninth order) produces a great variety, so \(p(\mathcal D\given\mathcal M_3)\) is spread over a large region. Because each distribution integrates to one, the broad one is low everywhere. A data set \(\mathcal D_0\) in the region where the simple model is too rigid but the complex one is too diffuse has the largest evidence under the model of intermediate complexity \(\mathcal M_2\). Complexity is penalised not by counting parameters but by the normalisation of \(p(\mathcal D\given\mathcal M)\). (The book's Figure 3.13 is a schematic and is not reproduced numerically here; Figure 3.14 below is its numerical counterpart.)

Proof

The left side of (3.73) is the definition of \(\KL(p_1\|p_2)\) with \(p_i=p(\cdot\given\mathcal M_i)\). Non-negativity: let \(f=-\ln\) (strictly convex). By Jensen's inequality (toolkit, KL), \(\KL(p_1\|p_2)\)\({}=\E_{p_1}[-\ln(p_2/p_1)]\)\({}\ge-\ln\E_{p_1}[p_2/p_1]\)\({}=-\ln\int p_2\,\d\mathcal D\)\({}\ge-\ln1\)\({}=0\), the last inequality because \(\int_{\{p_1>0\}}p_2\le1\); equality in Jensen needs \(p_2/p_1\) constant \(p_1\)-a.e., so \(p_2=p_1\). For the second claim, \(\E_{p_1}[p_2/p_1]=\int_{\{p_1>0\}}p_2=1\) under the support condition. The contrast between the two statements is itself instructive: the Bayes factor \(p_2/p_1\) is 1 on average, yet its logarithm is negative on average (Jensen), so the wrong model usually loses and occasionally wins big. ∎

This relies on the true distribution being among the models considered (the "closed-world" assumption). Two further cautions from the book: the evidence is sensitive to the prior's details (Figure 3.12, e.g. the tails), and it is undefined for an improper prior (an arbitrary scale factor, or, in the limit of a Gaussian prior with infinite variance, the evidence tends to zero as (3.70) shows); one can still compare an evidence ratio of two models first and take the limit afterwards. Hence, in applications, keep an independent test set to assess the final system.

3.5 The evidence approximation

pp. 165–172 A fully Bayesian treatment would also put hyperpriors on \(\alpha\) and \(\beta\) and predict by

\[ p(t\given\mathbf{t})=\iiint p(t\given\mathbf{w},\beta)\,p(\mathbf{w}\given\mathbf{t},\alpha,\beta)\,p(\alpha,\beta\given\mathbf{t})\,\d\mathbf{w}\,\d\alpha\,\d\beta. \tag{3.74} \]

One can integrate analytically over \(\mathbf{w}\) or over the hyperparameters, but not over everything. The evidence approximation (also called empirical Bayes, type-2 maximum likelihood or generalised maximum likelihood; Gull, 1989; MacKay, 1992) assumes that \(p(\alpha,\beta\given\mathbf{t})\) is sharply peaked at \(\hat\alpha,\hat\beta\), so that

\[ p(t\given\mathbf{t})\simeq p(t\given\mathbf{t},\hat\alpha,\hat\beta)=\int p(t\given\mathbf{w},\hat\beta)\,p(\mathbf{w}\given\mathbf{t},\hat\alpha,\hat\beta)\,\d\mathbf{w}. \tag{3.75} \]

By Bayes' theorem \(p(\alpha,\beta\given\mathbf{t})\)\({}\propto p(\mathbf{t}\given\alpha,\beta)\,p(\alpha,\beta)\) (3.76), so if the hyperprior is relatively flat, \(\hat\alpha,\hat\beta\) maximise the marginal likelihood \(p(\mathbf{t}\given\alpha,\beta)\), obtained by first integrating out \(\mathbf{w}\). This sets the hyperparameters, and hence the effective regularisation \(\lambda=\alpha/\beta\), from the training data alone. (With Gamma hyperpriors, the integral over \(\alpha,\beta\) gives a Student t over \(\mathbf{w}\); the remaining integral over \(\mathbf{w}\) could then be Laplace-approximated, but its skewed integrand makes that approximation poorer than the evidence procedure (MacKay, 1999).) The maximisation can be done either by differentiating the evidence analytically, as below, or by EM (Section 9.3.4), which reaches the same solution.

3.5.1 Evaluation of the evidence function

pp. 166–168 The marginal likelihood integrates the weights out,

\[ p(\mathbf{t}\given\alpha,\beta)=\int p(\mathbf{t}\given\mathbf{w},\beta)\,p(\mathbf{w}\given\alpha)\,\d\mathbf{w}. \tag{3.77} \]

Proof

(3.78). From (3.11), (3.12): \(p(\mathbf{t}\given\mathbf{w},\beta)\)\({}=(\beta/2\pi)^{N/2}\exp\{-\beta E_D(\mathbf{w})\}\), and from (3.52): \(p(\mathbf{w}\given\alpha)\)\({}=(\alpha/2\pi)^{M/2}\exp\{-\alpha E_W(\mathbf{w})\}\). Multiply and integrate; the exponents add to \(-E(\mathbf{w})\) (Exercise 3.17).

(3.80)–(3.84). Expand: \(E(\mathbf{w})\)\({}=\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\beta\mathbf{t}^\T\bm\Phi\mathbf{w}\)\({}+\frac12\mathbf{w}^\T(\beta\bm\Phi^\T\bm\Phi+\alpha\mathbf{I})\mathbf{w}\)\({}=\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\beta\mathbf{t}^\T\bm\Phi\mathbf{w}\)\({}+\frac12\mathbf{w}^\T\mathbf{A}\mathbf{w}\). The Hessian is \(\mathbf{A}\) (matrix calculus), positive definite for \(\alpha>0\). Since \(\nabla E\)\({}=\mathbf{A}\mathbf{w}\)\({}-\beta\bm\Phi^\T\mathbf{t}\), the unique minimiser is \(\mathbf{m}_N\)\({}=\beta\mathbf{A}^{-1}\bm\Phi^\T\mathbf{t}\), and a quadratic equals its minimum plus half the Hessian form in the displacement: \(E(\mathbf{w})\)\({}=E(\mathbf{m}_N)\)\({}+\frac12(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{A}(\mathbf{w}-\mathbf{m}_N)\), which is (3.80). To get the value (3.82), evaluate \(E\) at \(\mathbf{m}_N\): this is directly \(\frac\beta2\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2\)\({}+\frac\alpha2\mathbf{m}_N^\T\mathbf{m}_N\) by the definition (3.79). Equivalently, using \(\mathbf{A}\mathbf{m}_N=\beta\bm\Phi^\T\mathbf{t}\): \(E(\mathbf{m}_N)\)\({}=\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\beta\mathbf{t}^\T\bm\Phi\mathbf{m}_N\)\({}+\frac12\mathbf{m}_N^\T\mathbf{A}\mathbf{m}_N\)\({}=\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\frac12\mathbf{m}_N^\T\mathbf{A}\mathbf{m}_N\) (Exercise 3.18). That \(\mathbf{m}_N\) equals the mean (3.53) follows from \(\mathbf{A}=\mathbf{S}_N^{-1}\).

(3.85), (3.86). The remaining integral is a Gaussian integral, \(\int\exp\{-\frac12(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{A}(\mathbf{w}-\mathbf{m}_N)\}\d\mathbf{w}\)\({}=(2\pi)^{M/2}|\mathbf{A}|^{-1/2}\), the reciprocal of the normaliser of \(\N(\mathbf{m}_N,\mathbf{A}^{-1})\). Taking logs of (3.78) with this value: \(\frac N2\ln\beta-\frac N2\ln2\pi\)\({}+\frac M2\ln\alpha-\frac M2\ln2\pi\)\({}-E(\mathbf{m}_N)+\frac M2\ln2\pi\)\({}-\frac12\ln|\mathbf{A}|\); the two \(\frac M2\ln2\pi\) terms cancel, giving (3.86). ∎

3.5.2 Maximizing the evidence function

pp. 168–170 Define the eigenvalues of the data term by

\[ \big(\beta\bm\Phi^\T\bm\Phi\big)\mathbf{u}_i=\lambda_i\mathbf{u}_i. \tag{3.87} \]

The eigenvectors of \(\mathbf{A}\)\({}=\alpha\mathbf{I}+\beta\bm\Phi^\T\bm\Phi\) are the same \(\mathbf{u}_i\), with eigenvalues \(\alpha+\lambda_i\), because adding \(\alpha\mathbf{I}\) shifts every eigenvalue by \(\alpha\).

Proof

The \(\alpha\) equation. Differentiate the four \(\alpha\)-dependent terms of (3.86). First, \(\frac{\d}{\d\alpha}\frac M2\ln\alpha\)\({}=\frac M{2\alpha}\). Second, by Gap 3.9, \(-\frac{\d}{\d\alpha}E(\mathbf{m}_N)\)\({}=-\frac12\mathbf{m}_N^\T\mathbf{m}_N\). Third, since \(|\mathbf{A}|=\prod_i(\lambda_i+\alpha)\) (the determinant is the product of the eigenvalues) and \(\lambda_i=\beta d_i\) does not depend on \(\alpha\),

\[ \frac{\d}{\d\alpha}\ln|\mathbf{A}|=\frac{\d}{\d\alpha}\sum_i\ln(\lambda_i+\alpha)=\sum_i\frac1{\lambda_i+\alpha}. \tag{3.88} \]

Therefore \(\dfrac{\d\ln p}{\d\alpha}\)\({}=\dfrac M{2\alpha}\)\({}-\dfrac12\mathbf{m}_N^\T\mathbf{m}_N\)\({}-\dfrac12\sum_i\dfrac1{\lambda_i+\alpha}\) (3.89). Multiplying the stationarity condition by \(2\alpha\), \(0\)\({}=M-\alpha\mathbf{m}_N^\T\mathbf{m}_N\)\({}-\alpha\sum_i\frac1{\lambda_i+\alpha}\), i.e. \(\alpha\mathbf{m}_N^\T\mathbf{m}_N\)\({}=M-\alpha\sum_i\frac1{\lambda_i+\alpha}\)\({}=\gamma\) (3.90), where we used that the sum has \(M\) terms to write \(M\)\({}=\sum_i\frac{\lambda_i+\alpha}{\lambda_i+\alpha}\), so \(M-\alpha\sum_i\frac1{\lambda_i+\alpha}\)\({}=\sum_i\frac{\lambda_i}{\lambda_i+\alpha}\)\({}=\gamma\) (3.91). This is (3.92).

The \(\beta\) equation. The eigenvalues are proportional to \(\beta\), so \(\d\lambda_i/\d\beta=\lambda_i/\beta\), and

\[ \frac{\d}{\d\beta}\ln|\mathbf{A}|=\sum_i\frac{1}{\lambda_i+\alpha}\frac{\lambda_i}\beta=\frac\gamma\beta. \tag{3.93} \]

With \(\frac{\d}{\d\beta}\frac N2\ln\beta\)\({}=\frac N{2\beta}\) and \(\frac{\d E(\mathbf{m}_N)}{\d\beta}\)\({}=\frac12\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2\) (Gap 3.9):

\[ 0=\frac N{2\beta}-\frac12\sum_{n=1}^N\big(t_n-\mathbf{m}_N^\T\bm\phi_n\big)^2-\frac\gamma{2\beta}. \tag{3.94} \]

Multiplying by \(2\beta\) and solving for \(1/\beta\): \(N-\gamma=\beta\sum_n(\cdot)^2\), which is (3.95).

Consistency check. Substituting (3.92) and (3.95) into (3.82): \(E(\mathbf{m}_N)\)\({}=\frac\beta2\|\mathbf{t}\)\({}-\bm\Phi\mathbf{m}_N\|^2\)\({}+\frac\alpha2\mathbf{m}_N^\T\mathbf{m}_N\)\({}=\frac{N-\gamma}2+\frac\gamma2\)\({}=\frac N2\). ∎

3.5.3 Effective number of parameters

pp. 170–172 Rotate the axes of weight space to the eigenvectors \(\mathbf{u}_i\) of (3.87). Likelihood contours become axis-aligned ellipses, the eigenvalue \(\lambda_i\) is the curvature of the (scaled) log likelihood in direction \(i\) (small curvature = elongated contour), and the prior contours are circles.

Proof

(a) Each ratio \(\lambda_i/(\lambda_i+\alpha)\in[0,1)\) for \(\lambda_i\ge0\), so \(0\le\gamma\le M\); for \(\lambda_i>0\) it decreases strictly in \(\alpha\), from 1 to 0, and for \(\lambda_i=0\) it is identically 0. (b) \(\mathbf{m}_N\)\({}=\beta\mathbf{A}^{-1}\bm\Phi^\T\mathbf{t}\)\({}=\mathbf{A}^{-1}(\beta\bm\Phi^\T\bm\Phi)\mathbf{w}_{\rm ML}\) since \(\bm\Phi^\T\mathbf{t}\)\({}=\bm\Phi^\T\bm\Phi\mathbf{w}_{\rm ML}\); apply \(\mathbf{u}_i^\T\) using \(\mathbf{A}^{-1}\mathbf{u}_i\)\({}=\mathbf{u}_i/(\alpha+\lambda_i)\) and \(\beta\bm\Phi^\T\bm\Phi\mathbf{u}_i\)\({}=\lambda_i\mathbf{u}_i\). So parameters with \(\lambda_i\gg\alpha\) (well determined by the data) keep their ML values, those with \(\lambda_i\ll\alpha\) are set near zero by the prior (Figure 3.15). (c) The posterior covariance \(\mathbf{S}_N=\mathbf{A}^{-1}\) has eigenvalue \(1/(\alpha+\lambda_i)\), the prior \(1/\alpha\). (d) was shown in Result 3.15(d). (e) With \(\alpha\to0\), each \(\lambda_i/(\lambda_i+\alpha)\to1\), and (3.95) has \(N-M\) in the denominator. ∎

Why \(N-\gamma\) in (3.95). Compare (3.95) with the maximum likelihood (3.21): the divisor \(N\) is replaced by \(N-\gamma\). The univariate analogue: the ML variance \(\sigma^2_{\rm ML}\)\({}=\frac1N\sum_n(x_n-\mu_{\rm ML})^2\) (3.96) is biased because \(\mu_{\rm ML}\) has absorbed part of the noise, using up one degree of freedom, and the unbiased estimate divides by \(N-1\) (3.97). In the regression model the mean function has \(M\) parameters, but only \(\gamma\) of them are determined by the data (the other \(M-\gamma\) are held near zero by the prior and absorb little noise), so the correction is by \(\gamma\), not \(M\). Erratum. The first sentence after (3.97) in the book is to be omitted in corrected printings.

Proof

The matrix \(\bm\Phi^\T\bm\Phi=\sum_n\bm\phi_n\bm\phi_n^\T\) is a sum over data points; if the inputs are i.i.d. and \(\E[\bm\phi\bm\phi^\T]\) is positive definite, its eigenvalues are of order \(N\), so \(\lambda_i=\beta d_i\) grows without bound with \(N\) while \(\alpha\) stays fixed. Then \(\lambda_i/(\lambda_i+\alpha)=1-O(\alpha/N)\) for every \(i\), \(\gamma\)\({}=M-\alpha\sum_i(\lambda_i+\alpha)^{-1}\)\({}=M-O(1/N)\)\({}\to M\). Put \(\gamma=M\) in (3.92) and (3.95) and \(N-\gamma\simeq N\) (since \(M\ll N\)): (3.98) and (3.99), where \(2E_W=\mathbf{m}_N^\T\mathbf{m}_N\) and \(2E_D=\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2\). ∎

3.6 Limitations of fixed basis functions

pp. 172–173 Models that are linear combinations of fixed nonlinear basis functions have closed-form least squares and a tractable Bayesian treatment, and with suitable basis functions can approximate any reasonable mapping. They also have a basic weakness: the basis functions are chosen before the data are seen, which is a form of the curse of dimensionality (Section 1.4). The number of basis functions must grow quickly, often exponentially, with the input dimension \(D\).

Two properties of real data help. First, the inputs typically lie near a lower-dimensional nonlinear manifold (strong correlations among inputs; images of handwritten digits in Chapter 12), so local basis functions need only be placed in regions that contain data. This is the idea of radial basis function networks, support vector machines and relevance vector machines (Chapters 6 and 7). Neural networks, with adaptive sigmoidal basis functions, can instead adapt the directions in input space to which the basis functions respond so that they follow the manifold (Chapter 5). Second, the targets may depend on only a few directions within the manifold, which neural networks can also exploit. Together with the over-confident extrapolation noted in Section 3.3.2, these are the reasons the course turns next to classification (Chapter 4) and then to kernels, sparsity and adaptive bases.

Exercises

All 24 exercises of Chapter 3 (Bishop, 2006, pp. 173–177), 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 then give the argument in full or point to the exact result that contains it.

1. Exercise 3.1 (★, www). Show that \(\tanh a=2\sigma(2a)-1\) for the logistic sigmoid \(\sigma\) of (3.6), and hence that \(y(x,\mathbf{w})\)\({}=w_0\)\({}+\sum_{j=1}^Mw_j\sigma\big(\frac{x-\mu_j}s\big)\) (3.101) equals \(y(x,\mathbf{u})\)\({}=u_0\)\({}+\sum_{j=1}^Mu_j\tanh\big(\frac{x-\mu_j}{2s}\big)\) (3.102). Find \(\mathbf{u}\) in terms of \(\mathbf{w}\). Erratum: the book's (3.102) lacks the factor 2 in the denominator of the argument of tanh, and its last line should read \(u_0,w_0\) where it says \(u_1,w_1\).

   Solution

   \(2\sigma(2a)-1\)\({}=\dfrac2{1+e^{-2a}}-1\)\({}=\dfrac{1-e^{-2a}}{1+e^{-2a}}\)\({}=\dfrac{e^a-e^{-a}}{e^a+e^{-a}}\)\({}=\tanh a\) (multiply numerator and denominator by \(e^a\)). Hence \(\sigma(2a)=\tfrac12(1+\tanh a)\). Put \(a_j=\frac{x-\mu_j}{2s}\), so \(\sigma(\frac{x-\mu_j}s)\)\({}=\sigma(2a_j)\)\({}=\frac12+\frac12\tanh a_j\), and \[ y(x,\mathbf{w})=w_0+\sum_jw_j\Big(\frac12+\frac12\tanh a_j\Big)=\Big(w_0+\frac12\sum_jw_j\Big)+\sum_j\frac{w_j}2\tanh\Big(\frac{x-\mu_j}{2s}\Big). \] So \(u_j=w_j/2\) for \(j\ge1\) and \(u_0=w_0+\frac12\sum_{j=1}^Mw_j\). The map \(\mathbf{w}\mapsto\mathbf{u}\) is invertible (\(w_j=2u_j\), \(w_0=u_0-\sum_ju_j\)), so the two families represent exactly the same set of functions.

2. Exercise 3.2 (★★). Show that \(\bm\Phi(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\) (3.103) projects any vector onto the space spanned by the columns of \(\bm\Phi\), and use this to show that the least-squares solution (3.15) is the orthogonal projection of \(\mathbf{t}\) onto the subspace \(\mathcal S\) of Figure 3.2.

   Solution

   This is Result 3.3: \(\mathbf{P}\) is symmetric and idempotent, its range is \(\operatorname{col}(\bm\Phi)\), and \(\mathbf{v}-\mathbf{P}\mathbf{v}\) is orthogonal to every column. Since \(\bm\Phi\mathbf{w}_{\rm ML}\)\({}=\bm\Phi(\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\mathbf{t}\)\({}=\mathbf{P}\mathbf{t}\), the fitted vector is the orthogonal projection of \(\mathbf{t}\) onto \(\mathcal S\), the nearest point of \(\mathcal S\) to \(\mathbf{t}\).

3. Exercise 3.3 (★). Each data point has a weight \(r_n>0\) and the error is \(E_D(\mathbf{w})\)\({}=\frac12\sum_nr_n\{t_n\)\({}-\mathbf{w}^\T\bm\phi(\mathbf{x}_n)\}^2\) (3.104). Find the minimiser, and give two interpretations: (i) data-dependent noise variance, (ii) replicated data points.

   Solution

   With \(\mathbf{R}=\diag(r_1,\dots,r_N)\), \(E_D\)\({}=\frac12(\mathbf{t}-\bm\Phi\mathbf{w})^\T\mathbf{R}(\mathbf{t}-\bm\Phi\mathbf{w})\), with gradient \(-\bm\Phi^\T\mathbf{R}(\mathbf{t}-\bm\Phi\mathbf{w})\). Setting it to zero,

   \[ \mathbf{w}^\star=(\bm\Phi^\T\mathbf{R}\bm\Phi)^{-1}\bm\Phi^\T\mathbf{R}\mathbf{t}, \]

   the minimiser since the Hessian \(\bm\Phi^\T\mathbf{R}\bm\Phi\) is positive definite (for independent columns). (i) If \(t_n\)\({}\sim\N(\mathbf{w}^\T\bm\phi_n,(\beta r_n)^{-1})\), i.e. point \(n\) has its own noise precision \(\beta r_n\), the log likelihood is \(\sum_n\frac12\ln(\beta r_n)\)\({}-\beta E_D(\mathbf{w})+\const\), so its maximiser in \(\mathbf{w}\) is \(\mathbf{w}^\star\): precisely measured points (large \(r_n\)) count more. (ii) For integer \(r_n\), the data set in which point \(n\) is repeated \(r_n\) times has plain sum of squares \(\frac12\sum_nr_n(t_n-\mathbf{w}^\T\bm\phi_n)^2\), the same function.

4. Exercise 3.4 (★, www). For \(y(\mathbf{x},\mathbf{w})=w_0+\sum_{i=1}^Dw_ix_i\) and the sum of squares (3.106), add independent noise \(\epsilon_i\) (mean 0, variance \(\sigma^2\)) to each input. Show that the error averaged over the noise equals the noise-free error plus a weight-decay term that omits \(w_0\).

   Solution

   With noisy inputs \(\tilde x_{ni}=x_{ni}+\epsilon_{ni}\), the model output is \(\tilde y_n\)\({}=w_0+\sum_iw_i(x_{ni}+\epsilon_{ni})\)\({}=y_n+\sum_iw_i\epsilon_{ni}\). Then

   \[ \tilde E=\frac12\sum_n(\tilde y_n-t_n)^2=\frac12\sum_n\Big[(y_n-t_n)^2+2(y_n-t_n)\sum_iw_i\epsilon_{ni}+\Big(\sum_iw_i\epsilon_{ni}\Big)^2\Big]. \]

   Take expectations over the \(\epsilon\): the middle term has zero mean since \(\E\epsilon_{ni}=0\) and \(y_n,t_n\) are not random; for the last, \(\E[(\sum_iw_i\epsilon_{ni})^2]\)\({}=\sum_{i,j}w_iw_j\E[\epsilon_{ni}\epsilon_{nj}]\)\({}=\sigma^2\sum_iw_i^2\) using \(\E[\epsilon_i\epsilon_j]=\delta_{ij}\sigma^2\). Hence

   \[ \E[\tilde E]=E_D+\frac{N\sigma^2}2\sum_{i=1}^Dw_i^2 . \]

   (The book's solution writes the penalty as \(\frac12\sum_iw_i^2\sigma^2\) per data point, i.e. summed over the \(N\) points gives the factor \(N\) above.) This is the sum of squares for noise-free inputs plus weight decay with \(\lambda=N\sigma^2\), and \(w_0\) is absent because it multiplies no noisy input. Training with input noise is therefore equivalent to Tikhonov regularisation (Bishop, 1995).

5. Exercise 3.5 (★, www). Using Lagrange multipliers, show that minimising the regularised error (3.29) is equivalent to minimising the sum of squares (3.12) subject to (3.30), and discuss the relation between \(\eta\) and \(\lambda\).

   Solution

   Write the constraint as \(\frac12(\sum_j|w_j|^q-\eta)\le0\). The Lagrangian is \(L(\mathbf{w},\lambda)\)\({}=E_D(\mathbf{w})\)\({}+\frac\lambda2(\sum_j|w_j|^q-\eta)\), which in its dependence on \(\mathbf{w}\) is (3.29) (plus the constant \(-\lambda\eta/2\)). For fixed \(\lambda\ge0\), minimising \(L\) over \(\mathbf{w}\) is therefore minimising (3.29); complementary slackness fixes \(\eta=\sum_j|w_j^\star(\lambda)|^q\) when \(\lambda>0\). \(\eta(\lambda)\) is non-increasing in \(\lambda\), and in the convex case \(q\ge1\) every binding \(\eta\) corresponds to some \(\lambda\). Full argument, including the monotonicity and the envelope relation \(\lambda/2=-\d E_D^\star/\d\eta\): Result 3.6.

6. Exercise 3.6 (★, www). For the multivariate model \(p(\mathbf{t}\given\mathbf{W},\bm\Sigma)\)\({}=\N(\mathbf{t}\given\mathbf{W}^\T\bm\phi(\mathbf{x}),\bm\Sigma)\) (3.107, 3.108), show that each column of \(\mathbf{W}_{\rm ML}\) is given by (3.15), independently of \(\bm\Sigma\), and that \(\bm\Sigma_{\rm ML}\)\({}=\frac1N\sum_n(\mathbf{t}_n-\mathbf{W}_{\rm ML}^\T\bm\phi_n)(\mathbf{t}_n-\mathbf{W}_{\rm ML}^\T\bm\phi_n)^\T\) (3.109).

   Solution

   With residuals \(\mathbf{r}_n=\mathbf{t}_n-\mathbf{W}^\T\bm\phi_n\), \(\ln L\)\({}=-\frac N2\ln|\bm\Sigma|\)\({}-\frac12\sum_n\mathbf{r}_n^\T\bm\Sigma^{-1}\mathbf{r}_n\)\({}+\const\). Since \(\partial\mathbf{r}_n/\partial\mathbf{W}\) gives \(\nabla_{\mathbf{W}}(\frac12\mathbf{r}^\T\bm\Sigma^{-1}\mathbf{r})\)\({}=-\bm\phi\,\mathbf{r}^\T\bm\Sigma^{-1}\) (an \(M\times K\) matrix; matrix calculus), the stationarity condition is \(\sum_n\bm\phi_n\mathbf{r}_n^\T\bm\Sigma^{-1}=0\). Multiply on the right by \(\bm\Sigma\) (invertible): \(\sum_n\bm\phi_n\mathbf{r}_n^\T\)\({}=\bm\Phi^\T(\mathbf{T}-\bm\Phi\mathbf{W})\)\({}=0\), so \(\bm\Phi^\T\bm\Phi\mathbf{W}=\bm\Phi^\T\mathbf{T}\), i.e. \(\mathbf{W}_{\rm ML}\) of (3.34), whose \(k\)th column is (3.15) for \(\mathbf{t}_k\). \(\bm\Sigma\) has dropped out. For \(\bm\Sigma\): with \(\mathbf{W}=\mathbf{W}_{\rm ML}\) and \(\bm\Lambda=\bm\Sigma^{-1}\), the log likelihood is \(\frac N2\ln|\bm\Lambda|\)\({}-\frac12\Tr(\bm\Lambda\sum_n\mathbf{r}_n\mathbf{r}_n^\T)\); its derivative with respect to \(\bm\Lambda\) is \(\frac N2\bm\Lambda^{-1}\)\({}-\frac12\sum_n\mathbf{r}_n\mathbf{r}_n^\T\), zero at \(\bm\Sigma\)\({}=\frac1N\sum_n\mathbf{r}_n\mathbf{r}_n^\T\), which is (3.109) (the standard result for the covariance of a Gaussian, Section 2.3.4). The maximisation is joint, and since \(\mathbf{W}_{\rm ML}\) does not depend on \(\bm\Sigma\) it is the same \(\mathbf{W}_{\rm ML}\) that appears. See also Result 3.7.

7. Exercise 3.7 (★). By completing the square verify (3.49) for the posterior with \(\mathbf{m}_N\), \(\mathbf{S}_N\) of (3.50), (3.51).

   Solution

   This is the proof of Result 3.11: the log of likelihood times prior is quadratic in \(\mathbf{w}\) with precision \(\mathbf{S}_0^{-1}+\beta\bm\Phi^\T\bm\Phi\) and linear coefficient \(\mathbf{S}_0^{-1}\mathbf{m}_0\)\({}+\beta\bm\Phi^\T\mathbf{t}\); completing the square identifies the Gaussian with these mean and covariance.

8. Exercise 3.8 (★★, www). Having observed \(N\) points, treat the posterior (3.49) as the prior for a new point \((\mathbf{x}_{N+1},t_{N+1})\) and show by completing the square that the new posterior is (3.49) with \(\mathbf{S}_{N+1}\), \(\mathbf{m}_{N+1}\).

   Solution

   Multiplying \(\N(\mathbf{w}\given\mathbf{m}_N,\mathbf{S}_N)\) by \(\exp\{-\frac\beta2(t_{N+1}-\mathbf{w}^\T\bm\phi_{N+1})^2\}\), the exponent (times \(-2\)) is \(\mathbf{w}^\T(\mathbf{S}_N^{-1}+\beta\bm\phi_{N+1}\bm\phi_{N+1}^\T)\mathbf{w}\)\({}-2\mathbf{w}^\T(\mathbf{S}_N^{-1}\mathbf{m}_N+\beta\bm\phi_{N+1}t_{N+1})\)\({}+\const\). Reading off: \(\mathbf{S}_{N+1}^{-1}\)\({}=\mathbf{S}_N^{-1}\)\({}+\beta\bm\phi_{N+1}\bm\phi_{N+1}^\T\) and \(\mathbf{m}_{N+1}\)\({}=\mathbf{S}_{N+1}(\mathbf{S}_N^{-1}\mathbf{m}_N+\beta\bm\phi_{N+1}t_{N+1})\). This is Result 3.12, which also gives the equivalent gain form and shows that the result equals the batch posterior.

9. Exercise 3.9 (★★). Repeat Exercise 3.8 using the linear-Gaussian result (2.116).

   Solution

   Take the marginal \(p(\mathbf{w})\)\({}=\N(\mathbf{w}\given\bm\mu,\bm\Lambda^{-1})\) with \(\bm\mu=\mathbf{m}_N\), \(\bm\Lambda=\mathbf{S}_N^{-1}\), and the conditional \(p(t\given\mathbf{w})\)\({}=\N(t\given\mathbf{A}\mathbf{w}+\mathbf{b},\mathbf{L}^{-1})\) with \(\mathbf{A}=\bm\phi_{N+1}^\T\), \(\mathbf{b}=0\), \(\mathbf{L}^{-1}=\beta^{-1}\). Book (2.116) gives the posterior \(\N(\mathbf{w}\given\mathbf{S}\{\mathbf{A}^\T\mathbf{L}(t-\mathbf{b})+\bm\Lambda\bm\mu\},\mathbf{S})\) with \(\mathbf{S}\)\({}=(\bm\Lambda+\mathbf{A}^\T\mathbf{L}\mathbf{A})^{-1}\). Substituting: \(\mathbf{S}\)\({}=(\mathbf{S}_N^{-1}+\beta\bm\phi_{N+1}\bm\phi_{N+1}^\T)^{-1}\)\({}=\mathbf{S}_{N+1}\) and the mean is \(\mathbf{S}_{N+1}(\beta\bm\phi_{N+1}t_{N+1}+\mathbf{S}_N^{-1}\mathbf{m}_N)\)\({}=\mathbf{m}_{N+1}\). (toolkit.)

10. Exercise 3.10 (★★, www). Using (2.115) to evaluate the integral (3.57), verify the predictive distribution (3.58) with variance (3.59).

    Solution

    By (3.3), (3.8) and (3.49) the integrand is \(\N(t\given\bm\phi(\mathbf{x})^\T\mathbf{w},\beta^{-1})\,\N(\mathbf{w}\given\mathbf{m}_N,\mathbf{S}_N)\). Matching the first factor to the book's (2.114) \(p(\mathbf{y}\given\mathbf{x})\)\({}=\N(\mathbf{y}\given\mathbf{A}\mathbf{x}+\mathbf{b},\mathbf{L}^{-1})\) with \(\mathbf{A}=\bm\phi(\mathbf{x})^\T\) and the second to (2.113) \(p(\mathbf{x})=\N(\bm\mu,\bm\Lambda^{-1})\), (2.115) gives \(\N(t\given\bm\phi^\T\mathbf{m}_N,\beta^{-1}+\bm\phi^\T\mathbf{S}_N\bm\phi)\). A second, direct proof is given in Result 3.13.

11. Exercise 3.11 (★★). Use \((\mathbf{M}+\mathbf{v}\mathbf{v}^\T)^{-1}\)\({}=\mathbf{M}^{-1}\)\({}-\frac{(\mathbf{M}^{-1}\mathbf{v})(\mathbf{v}^\T\mathbf{M}^{-1})}{1+\mathbf{v}^\T\mathbf{M}^{-1}\mathbf{v}}\) (3.110) to show \(\sigma_{N+1}^2(\mathbf{x})\)\({}\le\sigma_N^2(\mathbf{x})\) (3.111).

    Solution

    Take \(\mathbf{M}=\mathbf{S}_N^{-1}\), \(\mathbf{v}=\sqrt\beta\bm\phi_{N+1}\); then \(\mathbf{M}+\mathbf{v}\mathbf{v}^\T\)\({}=\mathbf{S}_{N+1}^{-1}\) and \(\mathbf{S}_{N+1}\)\({}=\mathbf{S}_N\)\({}-\frac{\beta\mathbf{S}_N\bm\phi\bm\phi^\T\mathbf{S}_N}{1+\beta\bm\phi^\T\mathbf{S}_N\bm\phi}\). Hence \(\sigma_N^2(\mathbf{x})\)\({}-\sigma_{N+1}^2(\mathbf{x})\)\({}=\frac{\beta(\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\phi_{N+1})^2}{1+\beta\bm\phi_{N+1}^\T\mathbf{S}_N\bm\phi_{N+1}}\)\({}\ge0\), a square over a number at least one (as \(\mathbf{S}_N\) is positive definite). See Result 3.14.

12. Exercise 3.12 (★★). With both \(\mathbf{w}\) and \(\beta\) unknown, the conjugate prior is \(p(\mathbf{w},\beta)\)\({}=\N(\mathbf{w}\given\mathbf{m}_0,\beta^{-1}\mathbf{S}_0)\Gam(\beta\given a_0,b_0)\) (3.112). Show that the posterior is \(\N(\mathbf{w}\given\mathbf{m}_N,\beta^{-1}\mathbf{S}_N)\Gam(\beta\given a_N,b_N)\) (3.113) and find the parameters.

    Solution

    Multiply the likelihood \((\beta/2\pi)^{N/2}\exp\{-\frac\beta2\|\mathbf{t}\)\({}-\bm\Phi\mathbf{w}\|^2\}\) by the prior

    \[ \Big(\frac\beta{2\pi}\Big)^{M/2}|\mathbf{S}_0|^{-1/2}\exp\Big\{-\frac\beta2(\mathbf{w}-\mathbf{m}_0)^\T\mathbf{S}_0^{-1}(\mathbf{w}-\mathbf{m}_0)\Big\}\cdot\frac{b_0^{a_0}}{\Gamma(a_0)}\beta^{a_0-1}e^{-b_0\beta}. \]

    The exponent in \(\mathbf{w}\) is \(-\frac\beta2\) times a quadratic with precision \(\mathbf{S}_0^{-1}+\bm\Phi^\T\bm\Phi\) (no factor \(\beta\), since \(\beta\) is pulled out) and linear coefficient \(\mathbf{S}_0^{-1}\mathbf{m}_0\)\({}+\bm\Phi^\T\mathbf{t}\). Completing the square (toolkit): \[ \mathbf{S}_N^{-1}=\mathbf{S}_0^{-1}+\bm\Phi^\T\bm\Phi,\qquad \mathbf{m}_N=\mathbf{S}_N(\mathbf{S}_0^{-1}\mathbf{m}_0+\bm\Phi^\T\mathbf{t}), \] \[ -\frac\beta2\Big[(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)+\mathbf{t}^\T\mathbf{t}+\mathbf{m}_0^\T\mathbf{S}_0^{-1}\mathbf{m}_0-\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N\Big]. \] The first piece, with the factor \(\beta^{M/2}\) present in the product, is the kernel of \(\N(\mathbf{w}\given\mathbf{m}_N,\beta^{-1}\mathbf{S}_N)\). What is left is a function of \(\beta\) alone, \[ \beta^{a_0-1+N/2}\exp\Big\{-\beta\Big(b_0+\frac12\big[\mathbf{t}^\T\mathbf{t}+\mathbf{m}_0^\T\mathbf{S}_0^{-1}\mathbf{m}_0-\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N\big]\Big)\Big\}, \] a Gamma kernel with \[ a_N=a_0+\frac N2,\qquad b_N=b_0+\frac12\Big(\mathbf{m}_0^\T\mathbf{S}_0^{-1}\mathbf{m}_0-\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N+\sum_{n=1}^Nt_n^2\Big). \] An equivalent form is \[ b_N=b_0+\frac12\Big[\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2+(\mathbf{m}_N-\mathbf{m}_0)^\T\mathbf{S}_0^{-1}(\mathbf{m}_N-\mathbf{m}_0)\Big]: \] the residual at the posterior mean plus the distance from the prior mean, so \(b_N\) accumulates half the "unexplained" squared error. (Both forms agree because the quadratic form in \(\mathbf{w}\) equals its value at \(\mathbf{m}_N\) plus the displacement term.) Hence the posterior has the same form: the prior is conjugate.

13. Exercise 3.13 (★★). Show that the predictive distribution for the model of Exercise 3.12 is a Student t, \(p(t\given\mathbf{x},\mathbf{t})\)\({}=\St(t\given\mu,\lambda,\nu)\) (3.114), and find \(\mu,\lambda,\nu\).

    Solution

    Given \(\beta\), the argument of Result 3.13 applies with covariances scaled by \(\beta^{-1}\): \(t=\mathbf{w}^\T\bm\phi+\epsilon\) with \(\mathbf{w}\)\({}\sim\N(\mathbf{m}_N,\beta^{-1}\mathbf{S}_N)\) and \(\epsilon\sim\N(0,\beta^{-1})\), so \(p(t\given\mathbf{x},\mathbf{t},\beta)\)\({}=\N\big(t\given\mathbf{m}_N^\T\bm\phi,\ \beta^{-1}(1+\bm\phi^\T\mathbf{S}_N\bm\phi)\big)\). Now integrate over \(\beta\sim\Gam(a_N,b_N)\): a Gaussian whose precision is \(\beta\Lambda\), with \(\Lambda=1/(1+\bm\phi^\T\mathbf{S}_N\bm\phi)\), mixed over a Gamma-distributed \(\beta\), is a Student t (Section 2.3.7: \(\int\N(x\given\mu,(\beta\Lambda)^{-1})\Gam(\beta\given a,b)\d\beta\)\({}=\St(x\given\mu,\Lambda a/b,2a)\)). Therefore \[ \mu=\mathbf{m}_N^\T\bm\phi(\mathbf{x}),\qquad \lambda=\frac{a_N}{b_N}\,\frac1{1+\bm\phi(\mathbf{x})^\T\mathbf{S}_N\bm\phi(\mathbf{x})},\qquad \nu=2a_N=2a_0+N. \] The degrees of freedom grow with \(N\), so the tails become Gaussian as data accumulate; \(b_N/a_N\) plays the role of the noise variance estimate.

14. Exercise 3.14 (★★). Let the basis functions be linearly independent, \(N>M\), and one of them constant. Take new basis functions \(\psi_j\) spanning the same space and orthonormal on the data, \(\sum_n\psi_j(\mathbf{x}_n)\psi_k(\mathbf{x}_n)\)\({}=I_{jk}\) (3.115). Show that for \(\alpha=0\) the equivalent kernel is \(k(\mathbf{x},\mathbf{x}')\)\({}=\bm\psi(\mathbf{x})^\T\bm\psi(\mathbf{x}')\), and deduce the summation constraint (3.116), \(\sum_nk(\mathbf{x},\mathbf{x}_n)=1\). Erratum: the book's range \(\{\psi_1,\dots,\psi_M\}\) in the exercise should be \(\{\psi_0,\dots,\psi_{M-1}\}\).

    Solution

    Let \(\bm\psi(\mathbf{x})\)\({}=\mathbf{R}\bm\phi(\mathbf{x})\) for an invertible \(M\times M\) matrix \(\mathbf{R}\) (a change of basis within the same span; for example Gram–Schmidt starting from the constant function, or \(\mathbf{R}=(\bm\Phi^\T\bm\Phi)^{-1/2}\)), and \(\bm\Psi=\bm\Phi\mathbf{R}^\T\). Orthonormality (3.115) says \(\bm\Psi^\T\bm\Psi\)\({}=\mathbf{R}\bm\Phi^\T\bm\Phi\mathbf{R}^\T\)\({}=\mathbf{I}\), so \((\bm\Phi^\T\bm\Phi)^{-1}=\mathbf{R}^\T\mathbf{R}\). For \(\alpha=0\), \(\mathbf{S}_N=(\beta\bm\Phi^\T\bm\Phi)^{-1}\) and \[ k(\mathbf{x},\mathbf{x}')=\beta\bm\phi^\T(\beta\bm\Phi^\T\bm\Phi)^{-1}\bm\phi'=\bm\phi^\T\mathbf{R}^\T\mathbf{R}\bm\phi'=\bm\psi(\mathbf{x})^\T\bm\psi(\mathbf{x}'). \] Summation. \(\sum_nk(\mathbf{x},\mathbf{x}_n)\)\({}=\sum_j\psi_j(\mathbf{x})\sum_n\psi_j(\mathbf{x}_n)\). The span contains the constant function, so we may take \(\psi_0\) constant, \(\psi_0\equiv c\); normalisation \(\sum_n\psi_0(\mathbf{x}_n)^2=Nc^2=1\) forces \(c=1/\sqrt N\). (The book's "\(\psi_0=1\)" is not compatible with (3.115) as printed, which requires unit norm; it holds if the normalisation is rescaled. The conclusion is unchanged.) By orthogonality to \(\psi_0\), for each \(j\): \(\sum_n\psi_j(\mathbf{x}_n)\)\({}=\frac1c\sum_n\psi_j(\mathbf{x}_n)\psi_0(\mathbf{x}_n)\)\({}=\delta_{j0}/c\). Hence \(\sum_nk(\mathbf{x},\mathbf{x}_n)\)\({}=\psi_0(\mathbf{x})/c\)\({}=1\). A direct proof is in Result 3.15(c).

15. Exercise 3.15 (★, www). With \(\alpha\) and \(\beta\) set by the evidence framework, show that \(E(\mathbf{m}_N)\) of (3.82) satisfies \(2E(\mathbf{m}_N)=N\).

    Solution

    At the stationary values, (3.92) gives \(\alpha\mathbf{m}_N^\T\mathbf{m}_N=\gamma\) and (3.95) gives \(\beta\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2\)\({}=N-\gamma\). Then \(2E(\mathbf{m}_N)\)\({}=\beta\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2\)\({}+\alpha\mathbf{m}_N^\T\mathbf{m}_N\)\({}=(N-\gamma)+\gamma\)\({}=N\).

16. Exercise 3.16 (★★). Derive the log evidence (3.86) by using (2.115) to evaluate the integral (3.77) directly.

    Solution

    See Filling the gap 3.8: (2.115) gives \(\mathbf{t}\)\({}\sim\N(0,\beta^{-1}\mathbf{I}+\alpha^{-1}\bm\Phi\bm\Phi^\T)\); the determinant lemma and the Woodbury identity turn its log density into (3.86).

17. Exercise 3.17 (★). Show that the evidence can be written in the form (3.78) with \(E(\mathbf{w})\) of (3.79).

    Solution

    \(p(\mathbf{t}\given\mathbf{w},\beta)\)\({}=(\beta/2\pi)^{N/2}e^{-\beta E_D(\mathbf{w})}\) by (3.10)–(3.12) and \(p(\mathbf{w}\given\alpha)\)\({}=(\alpha/2\pi)^{M/2}e^{-\alpha\mathbf{w}^\T\mathbf{w}/2}\) by (3.52). The product is \((\beta/2\pi)^{N/2}(\alpha/2\pi)^{M/2}\exp\{-\beta E_D\)\({}-\alpha E_W\}\), with \(E_W=\frac12\mathbf{w}^\T\mathbf{w}\); integrating over \(\mathbf{w}\) gives (3.78).

18. Exercise 3.18 (★★, www). By completing the square over \(\mathbf{w}\), show that \(E(\mathbf{w})\) of (3.79) equals (3.80).

    Solution

    Expand \(E(\mathbf{w})\)\({}=\frac12[\beta\mathbf{t}^\T\mathbf{t}-2\beta\mathbf{t}^\T\bm\Phi\mathbf{w}+\mathbf{w}^\T\mathbf{A}\mathbf{w}]\) with \(\mathbf{A}\)\({}=\alpha\mathbf{I}+\beta\bm\Phi^\T\bm\Phi\). Since \(\beta\bm\Phi^\T\mathbf{t}=\mathbf{A}\mathbf{m}_N\) (definition (3.84)), \(\beta\mathbf{t}^\T\bm\Phi\mathbf{w}\)\({}=\mathbf{m}_N^\T\mathbf{A}\mathbf{w}\); adding and subtracting \(\mathbf{m}_N^\T\mathbf{A}\mathbf{m}_N\): \[ E(\mathbf{w})=\tfrac12\beta\mathbf{t}^\T\mathbf{t}-\tfrac12\mathbf{m}_N^\T\mathbf{A}\mathbf{m}_N+\tfrac12(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{A}(\mathbf{w}-\mathbf{m}_N). \] The first two terms equal \(\frac\beta2\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2\)\({}+\frac\alpha2\mathbf{m}_N^\T\mathbf{m}_N\)\({}=E(\mathbf{m}_N)\): expand \(\frac\beta2(\mathbf{t}-\bm\Phi\mathbf{m}_N)^\T(\mathbf{t}-\bm\Phi\mathbf{m}_N)\)\({}+\frac\alpha2\mathbf{m}_N^\T\mathbf{m}_N\)\({}=\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\beta\mathbf{t}^\T\bm\Phi\mathbf{m}_N\)\({}+\frac12\mathbf{m}_N^\T\mathbf{A}\mathbf{m}_N\), and \(\beta\mathbf{t}^\T\bm\Phi\mathbf{m}_N\)\({}=\mathbf{m}_N^\T\mathbf{A}\mathbf{m}_N\), giving \(\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\frac12\mathbf{m}_N^\T\mathbf{A}\mathbf{m}_N\). So (3.80) and (3.82). See Result 3.18.

19. Exercise 3.19 (★★). Show that integrating over \(\mathbf{w}\) gives (3.85), hence the log marginal likelihood (3.86).

    Solution

    By (3.80), \[ \int e^{-E(\mathbf{w})}\d\mathbf{w}=e^{-E(\mathbf{m}_N)}\int\exp\Big\{-\frac12(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{A}(\mathbf{w}-\mathbf{m}_N)\Big\}\d\mathbf{w}, \] and the last integral is the reciprocal of the normalisation constant \((2\pi)^{-M/2}|\mathbf{A}|^{1/2}\) of \(\N(\mathbf{m}_N,\mathbf{A}^{-1})\), i.e. \((2\pi)^{M/2}|\mathbf{A}|^{-1/2}\): (3.85). Insert in (3.78) and take logs: \(\frac N2\ln\frac\beta{2\pi}\)\({}+\frac M2\ln\frac\alpha{2\pi}-E(\mathbf{m}_N)\)\({}+\frac M2\ln2\pi-\frac12\ln|\mathbf{A}|\), which is (3.86).

20. Exercise 3.20 (★★, www). Verify all steps leading from the maximisation of (3.86) with respect to \(\alpha\) to the re-estimation equation (3.92). Erratum: the book's wording "Starting from (3.86) verify ..." should read "Verify ...".

    Solution

    Only the first, third and fourth terms of (3.86) depend on \(\alpha\) (for fixed \(\beta\)). Eigenvalues of \(\mathbf{A}\) are \(\alpha+\lambda_i\) (the eigenvectors \(\mathbf{u}_i\) of \(\beta\bm\Phi^\T\bm\Phi\) are orthonormal and \(\alpha\mathbf{I}\) shifts them), giving (3.88). The derivative of \(E(\mathbf{m}_N)\) is exactly \(\frac12\mathbf{m}_N^\T\mathbf{m}_N\) (Gap 3.9; this justifies the book's "neglected" derivative). Combining gives (3.89), and multiplying by \(2\alpha\) and using \(M\)\({}=\sum_i\frac{\alpha+\lambda_i}{\alpha+\lambda_i}\) gives (3.90), (3.91), (3.92). Every step is written out in Result 3.19.

21. Exercise 3.21 (★★). Prove \(\frac{\d}{\d\alpha}\ln|\mathbf{A}|\)\({}=\Tr(\mathbf{A}^{-1}\frac{\d\mathbf{A}}{\d\alpha})\) (3.117) using the eigenvalue expansion of a real symmetric matrix, and use it to derive (3.92) from (3.86).

    Solution

    The proof is in Filling the gap 3.10 (differentiate \(\mathbf{A}\mathbf{e}_i=\mu_i\mathbf{e}_i\), obtain \(\mu_i'=\mathbf{e}_i^\T\mathbf{A}'\mathbf{e}_i\), and sum \(\mu_i'/\mu_i\)). With \(\mathbf{A}'=\mathbf{I}\), \(\Tr\mathbf{A}^{-1}=\sum_i(\alpha+\lambda_i)^{-1}\) and (3.89) follows as before; the rest is Result 3.19.

22. Exercise 3.22 (★★). Verify all steps leading from the maximisation of (3.86) with respect to \(\beta\) to (3.95).

    Solution

    The \(\beta\)-dependent terms of (3.86) are \(\frac N2\ln\beta\), \(-E(\mathbf{m}_N)\) and \(-\frac12\ln|\mathbf{A}|\). Their derivatives are \(\frac N{2\beta}\), \(-\frac12\|\mathbf{t}-\bm\Phi\mathbf{m}_N\|^2\) (envelope theorem, Gap 3.9) and \(-\frac\gamma{2\beta}\) (since \(\d\lambda_i/\d\beta=\lambda_i/\beta\), (3.93)). Setting the sum to zero is (3.94); rearranged it is (3.95). See Result 3.19.

23. Exercise 3.23 (★★, www). For the model of Exercise 3.12, show that the evidence is \[ p(\mathbf{t})=\frac1{(2\pi)^{N/2}}\frac{b_0^{a_0}}{b_N^{a_N}}\frac{\Gamma(a_N)}{\Gamma(a_0)}\frac{|\mathbf{S}_N|^{1/2}}{|\mathbf{S}_0|^{1/2}} \tag{3.118} \] by first marginalising over \(\mathbf{w}\) and then over \(\beta\).

    Solution

    Use the notation of Exercise 3.12 and write \(2(b_N-b_0)\)\({}=\mathbf{t}^\T\mathbf{t}\)\({}+\mathbf{m}_0^\T\mathbf{S}_0^{-1}\mathbf{m}_0\)\({}-\mathbf{m}_N^\T\mathbf{S}_N^{-1}\mathbf{m}_N\). The integrand over \(\mathbf{w}\) is \[ \Big(\frac\beta{2\pi}\Big)^{(N+M)/2}|\mathbf{S}_0|^{-1/2}\exp\Big\{-\frac\beta2\big[(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)+2(b_N-b_0)\big]\Big\}. \] The Gaussian integral over \(\mathbf{w}\) is \((2\pi/\beta)^{M/2}|\mathbf{S}_N|^{1/2}\), leaving \((\beta/2\pi)^{N/2}|\mathbf{S}_N|^{1/2}|\mathbf{S}_0|^{-1/2}e^{-\beta(b_N-b_0)}\). Multiply by the Gamma prior \(\frac{b_0^{a_0}}{\Gamma(a_0)}\beta^{a_0-1}e^{-b_0\beta}\) and integrate over \(\beta\): the exponent combines to \(e^{-b_N\beta}\) and the power of \(\beta\) is \(\beta^{a_0+N/2-1}=\beta^{a_N-1}\), with \(\int_0^\infty\beta^{a_N-1}e^{-b_N\beta}\d\beta\)\({}=\Gamma(a_N)/b_N^{a_N}\). Collecting constants gives (3.118).

24. Exercise 3.24 (★★). Repeat Exercise 3.23 using Bayes' theorem in the form \(p(\mathbf{t})\)\({}=p(\mathbf{t}\given\mathbf{w},\beta)p(\mathbf{w},\beta)/p(\mathbf{w},\beta\given\mathbf{t})\) (3.119).

    Solution

    The identity holds for every \((\mathbf{w},\beta)\); substitute the three normalised densities. Powers of \(\beta\): the numerator has \(\beta^{N/2}\) (likelihood), \(\beta^{M/2}\) (Gaussian prior), \(\beta^{a_0-1}\) (Gamma prior); the denominator has \(\beta^{M/2}\) (Gaussian posterior) and \(\beta^{a_N-1}\) (Gamma posterior). Their ratio is \(\beta^{N/2+a_0-a_N}=\beta^0\). Powers of \(2\pi\): \((2\pi)^{-N/2-M/2}\) over \((2\pi)^{-M/2}\) gives \((2\pi)^{-N/2}\). Determinants and Gamma constants: \(\frac{|\mathbf{S}_0|^{-1/2}}{|\mathbf{S}_N|^{-1/2}}\cdot\frac{b_0^{a_0}/\Gamma(a_0)}{b_N^{a_N}/\Gamma(a_N)}\). The exponents: the numerator has (Exercise 3.12) \[ -\frac\beta2\Big[\|\mathbf{t}-\bm\Phi\mathbf{w}\|^2+(\mathbf{w}-\mathbf{m}_0)^\T\mathbf{S}_0^{-1}(\mathbf{w}-\mathbf{m}_0)\Big]-b_0\beta=-\frac\beta2\Big[(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)+2(b_N-b_0)\Big]-b_0\beta, \] the denominator \(-\frac\beta2(\mathbf{w}-\mathbf{m}_N)^\T\mathbf{S}_N^{-1}(\mathbf{w}-\mathbf{m}_N)\)\({}-b_N\beta\). The \(\mathbf{w}\)-dependence cancels and so does the rest: \(-\beta(b_N-b_0)-b_0\beta+b_N\beta=0\). What remains is (3.118), as it must, because the left side does not depend on \((\mathbf{w},\beta)\).

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/ch03_linear_regression.py (see the code index).

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 3, pp. 137–177 (and the official errata). doi:10.1007/978-0-387-45528-0.
2. Bishop, C. M. (1995). Training with noise is equivalent to Tikhonov regularization. Neural Computation, 7(1), 108–116. doi:10.1162/neco.1995.7.1.108.
3. Geman, S., Bienenstock, E. and Doursat, R. (1992). Neural networks and the bias/variance dilemma. Neural Computation, 4(1), 1–58. doi:10.1162/neco.1992.4.1.1.
4. Gull, S. F. (1989). Developments in maximum entropy data analysis. In J. Skilling (ed.), Maximum Entropy and Bayesian Methods, 53–71. Kluwer, Dordrecht.
5. Hoerl, A. E. and Kennard, R. W. (1970). Ridge regression: biased estimation for nonorthogonal problems. Technometrics, 12(1), 55–67. doi:10.1080/00401706.1970.10488634.
6. Kass, R. E. and Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773–795. doi:10.1080/01621459.1995.10476572.
7. MacKay, D. J. C. (1992). Bayesian interpolation. Neural Computation, 4(3), 415–447. doi:10.1162/neco.1992.4.3.415.
8. MacKay, D. J. C. (1999). Comparison of approximate methods for handling hyperparameters. Neural Computation, 11(5), 1035–1068. doi:10.1162/089976699300016331.
9. Penrose, R. (1955). A generalized inverse for matrices. Mathematical Proceedings of the Cambridge Philosophical Society, 51(3), 406–413. doi:10.1017/S0305004100030401.
10. Qazaz, C. S., Williams, C. K. I. and Bishop, C. M. (1997). An upper bound on the Bayesian error bars for generalized linear regression. In S. W. Ellacott, J. C. Mason and I. J. Anderson (eds), Mathematics of Neural Networks: Models, Algorithms and Applications, 295–299. Kluwer. doi:10.1007/978-1-4615-6099-9_51.
11. Robbins, H. and Monro, S. (1951). A stochastic approximation method. The Annals of Mathematical Statistics, 22(3), 400–407. doi:10.1214/aoms/1177729586.
12. Sherman, J. and Morrison, W. J. (1950). Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. The Annals of Mathematical Statistics, 21(1), 124–127. doi:10.1214/aoms/1177729893.
13. Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B, 58(1), 267–288. doi:10.1111/j.2517-6161.1996.tb02080.x.

DOIs were checked against Crossref. The numerical statements made in the text (the evidence formula against the direct Gaussian marginal, the exactness of the envelope-theorem derivatives, the fixed point with \(2E(\mathbf{m}_N)=N\), the optimal ridge parameter averaged over the prior, the evidence (3.118) against numerical integration over \(\beta\), the summation constraint of the equivalent kernel) come from short Python computations written for these notes; the figures come from the course code.