Sparse Kernel Machines

The kernel methods of Chapter 6 need the kernel evaluated between a new point and every training point. This chapter looks for solutions that depend on only a few training points: the support vector machine, which gets sparsity from a convex optimisation with a flat region in its loss, and the relevance vector machine, which gets sparser solutions and probabilistic outputs from a Bayesian prior with one precision per weight.

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

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

7.1 Maximum margin classifiers

pp. 326–331 Take the two-class problem with a linear model in a fixed feature space,

\[ y(\mathbf{x})=\mathbf{w}^\T\bm\phi(\mathbf{x})+b, \tag{7.1} \]

targets \(t_n\in\{-1,+1\}\), and classify a new point by the sign of \(y(\mathbf{x})\). Assume for now that the training set is linearly separable in feature space: some \((\mathbf{w},b)\) has \(t_ny(\mathbf{x}_n)\gt0\) for every \(n\). Usually many such hyperplanes exist, and the perceptron (Section 4.1.7) returns an arbitrary one that depends on its starting point and the order of the data. The support vector machine singles out one by a geometric criterion.

Proof

(a) Scaling multiplies numerator and denominator of the objective in (7.3) by \(\kappa\). Let \(\rho^\star\) be the optimal value of (7.3), over the \((\mathbf{w},b)\) with \(t_ny_n\gt0\) for all \(n\). Take any \((\mathbf{w},b)\) feasible for (7.5), i.e. with \(\min_nt_ny_n\ge1\). Then \(1/\|\mathbf{w}\|\)\({}\le\min_nt_ny_n/\|\mathbf{w}\|\)\({}\le\rho^\star\). Conversely if \((\mathbf{w}^\star,b^\star)\) attains \(\rho^\star\), put \(m=\min_nt_ny_n^\star\gt0\) and rescale by \(\kappa=1/m\): the new pair satisfies (7.5) with equality for the closest point, which is (7.4), and \(1/\|\kappa\mathbf{w}^\star\|=\rho^\star\). So the largest value of \(1/\|\mathbf{w}\|\) over (7.5) is \(\rho^\star\), and the two problems have the same solutions up to the scaling that sets the closest point to \(t_ny_n=1\). Maximising \(\|\mathbf{w}\|^{-1}\) is minimising \(\|\mathbf{w}\|^2\); the factor \(\frac12\) only simplifies the gradient.

(b) The feasible set of (7.5) is closed, convex and (by separability and rescaling) nonempty, and \(\frac12\|\mathbf{w}\|^2\) is strictly convex and coercive in \(\mathbf{w}\), so the optimal \(\mathbf{w}\) exists and is unique. Take an optimal \((\mathbf{w},b)\). Active constraints in both classes. Suppose no positive point is active, i.e. \(y_n\gt1\) for all \(t_n=+1\). Lower \(b\) by a small \(\delta\gt0\): the positive constraints stay strict, and the negative ones (\(-y_n\ge1\)) become looser, so now every constraint is strict. Then \((1-\varepsilon)\mathbf{w}\) with small \(\varepsilon\) is still feasible (finitely many strict inequalities) and has smaller norm, contradicting optimality. The same argument with \(b\) raised excludes the case of no active negative point. Uniqueness of \(b\). With an active positive point \(n\) we have \(\mathbf{w}^\T\bm\phi_n+b=1\), so \(b=1-\mathbf{w}^\T\bm\phi_n\) is determined by the unique \(\mathbf{w}\).

(c) If \(t_ny_n\ge\gamma\) then \(t_n(y_n/\gamma)\ge1\); so \((\mathbf{w},b)\) is feasible for the \(\gamma\)-problem iff \((\mathbf{w}/\gamma,b/\gamma)\) is feasible for (7.5). The objective of the \(\gamma\)-problem is \(\frac12\|\mathbf{w}\|^2\)\({}=\gamma^2\cdot\frac12\|\mathbf{w}/\gamma\|^2\), so the minimiser is \(\gamma\) times the minimiser of (7.5). Both give the same set \(\{y=0\}\) and the same classifier. ∎

The bias \(b\) seems to have dropped out of the objective (7.6); it is determined implicitly through the constraints. To solve the problem, attach a multiplier \(a_n\ge0\) to each constraint (toolkit: Lagrange multipliers), with a minus sign because we minimise over \((\mathbf{w},b)\):

\[ L(\mathbf{w},b,\mathbf{a})=\frac12\|\mathbf{w}\|^2-\sum_{n=1}^{N}a_n\big\{t_n(\mathbf{w}^\T\bm\phi(\mathbf{x}_n)+b)-1\big\}. \tag{7.7} \]

Proof

Stationarity. \(\nabla_{\mathbf{w}}L\)\({}=\mathbf{w}-\sum_na_nt_n\bm\phi_n\) (the gradient of \(\frac12\mathbf{w}^\T\mathbf{w}\) is \(\mathbf{w}\); the bracketed terms are linear in \(\mathbf{w}\) with gradient \(t_n\bm\phi_n\)), which vanishes at (7.8); \(\partial L/\partial b=-\sum_na_nt_n\), which vanishes at (7.9).

Why (7.12) appears. For fixed \(\mathbf{a}\ge\mathbf{0}\), \(L\) is a convex function of \((\mathbf{w},b)\), linear in \(b\) with slope \(-\sum_na_nt_n\). If that slope is nonzero, \(\inf_bL=-\infty\), so the dual function \(g(\mathbf{a})=\inf_{\mathbf{w},b}L\) is finite only on the set (7.12); there the infimum is attained at (7.8).

Substitution. Insert (7.8) and use (7.9) so that the \(b\) term vanishes:

\[ \tfrac12\|\mathbf{w}\|^2=\tfrac12\sum_{n,m}a_na_mt_nt_m\bm\phi_n^\T\bm\phi_m,\qquad \sum_na_nt_n\mathbf{w}^\T\bm\phi_n=\sum_{n,m}a_na_mt_nt_m\bm\phi_m^\T\bm\phi_n,\qquad b\sum_na_nt_n=0. \]

Therefore \(L=\frac12Q-Q-0+\sum_na_n\) with \(Q=\sum_{n,m}a_na_mt_nt_mk_{nm}\), i.e. (7.10). Substituting (7.8) into (7.1) gives (7.13).

Concavity and the role of positive definiteness. \(Q=\|\sum_na_nt_n\bm\phi_n\|^2\ge0\), i.e. the matrix \(t_nt_mk_{nm}\) is positive semidefinite whenever \(k\) is a valid kernel (Chapter 6); so \(\widetilde L\) is concave and any local maximum is global, and \(\widetilde L\le\sum_na_n\) is bounded above on each bounded set of \(\mathbf{a}\). For an indefinite \(k\), \(Q\) could be negative along some feasible direction, and the dual could be unbounded, which is the role of "positive definite" in the book's remark.

No duality gap. The primal is a convex quadratic problem with affine constraints, which is feasible by separability; for such problems strong duality holds (the Slater condition is only needed for nonlinear constraints), so the dual maximum equals the primal minimum. ∎

The primal has \(M\) variables (the dimension of feature space) and the dual has \(N\); with fixed \(M\lt N\) the dual looks worse, but it needs only the kernel, so \(\bm\phi\) may be infinite dimensional. Note also that (7.8) is a representer theorem for this problem: \(\mathbf{w}\) lies in the span of the training features.

Derivation of the KKT conditions

(7.14) is dual feasibility and (7.15) is primal feasibility (7.5). For (7.16) (complementary slackness) let \((\mathbf{w}^\star,b^\star)\) be primal optimal and \(\mathbf{a}^\star\) dual optimal. Strong duality (Result 7.2) gives \(\frac12\|\mathbf{w}^\star\|^2\)\({}=g(\mathbf{a}^\star)\). Since \(g(\mathbf{a}^\star)\) is the infimum of \(L(\cdot,\cdot,\mathbf{a}^\star)\),

\[ \tfrac12\|\mathbf{w}^\star\|^2=g(\mathbf{a}^\star)\le L(\mathbf{w}^\star,b^\star,\mathbf{a}^\star)=\tfrac12\|\mathbf{w}^\star\|^2-\sum_na_n^\star\{t_ny_n^\star-1\}\le\tfrac12\|\mathbf{w}^\star\|^2, \]

the last step because each \(a_n^\star\ge0\) and \(t_ny_n^\star-1\ge0\). Both ends are equal, so \(\sum_na_n^\star\{t_ny_n^\star-1\}=0\); a sum of nonnegative terms is zero only if every term is zero, which is (7.16). ∎

This is the central practical fact: once trained, only the support vectors are kept. Geometrically they are the points on the two margin hyperplanes (Figure 7.1), and the other points can be moved around freely, as long as they stay outside the margin, without changing the solution. The bias follows from any support vector.

Derivation

(7.17) is (7.16) with \(a_n\gt0\) and (7.13) (the sum over all \(m\) reduces to \(m\in\mathcal{S}\) because \(a_m=0\) otherwise). Multiply (7.17) by \(t_n\) and use \(t_n^2=1\): \(\sum_ma_mt_mk_{nm}+b=t_n\), i.e. \(b=t_n-\sum_{m\in\mathcal{S}}a_mt_mk_{nm}\) for each \(n\in\mathcal{S}\). Each is an exact solution; averaging them gives (7.18). The average is preferred because \(a_n\) and the kernel values are known only to the accuracy of the quadratic-programming solver, and the average reduces the error of a single noisy equation. ∎

The maximum margin classifier is also a regularised error minimisation,

\[ \sum_{n=1}^{N}E_\infty\big(y(\mathbf{x}_n)t_n-1\big)+\lambda\|\mathbf{w}\|^2,\tag{7.19} \]

where \(E_\infty(z)=0\) for \(z\ge0\) and \(+\infty\) otherwise. The infinite penalty enforces (7.5), and any \(\lambda\gt0\) gives the same minimiser, so its value is irrelevant.

Proof

\(\rho=1/\|\mathbf{w}\|\) is (7.4) with the closest point at distance \(t_ny_n/\|\mathbf{w}\|=1/\|\mathbf{w}\|\); so \(\rho^{-2}=\|\mathbf{w}\|^2\), which is (7.125). Next, by (7.8), then \(\bm\phi_n^\T\mathbf{w}=y_n-b\), then (7.9),

\[ \|\mathbf{w}\|^2=\sum_na_nt_n\mathbf{w}^\T\bm\phi_n=\sum_na_nt_n(y_n-b)=\sum_na_nt_ny_n-b\sum_na_nt_n=\sum_na_nt_ny_n. \]

By (7.16), \(a_nt_ny_n=a_n\) for every \(n\) (if \(a_n=0\) both sides vanish; otherwise \(t_ny_n=1\)). So \(\|\mathbf{w}\|^2=\sum_na_n\), which is (7.123). Finally the second term of (7.10) is \(\frac12\|\mathbf{w}\|^2\) by (7.8), so \(\widetilde L\)\({}=\sum_na_n-\frac12\|\mathbf{w}\|^2\)\({}=\|\mathbf{w}\|^2-\frac12\|\mathbf{w}\|^2\)\({}=\frac12\|\mathbf{w}\|^2\), i.e. \(\rho^{-2}=2\widetilde L\), (7.124). ∎

7.1.1 Overlapping class distributions

pp. 331–336 If the classes overlap in feature space, exact separation does not exist (and, if it did, would generalise poorly). Introduce one slack variable \(\xi_n\ge0\) per training point: \(\xi_n=0\) on or inside the correct margin boundary and \(\xi_n=|t_n-y(\mathbf{x}_n)|\) otherwise. A point on the decision surface has \(\xi_n=1\), a misclassified one \(\xi_n\gt1\) (Figure 7.3). The constraints (7.5) become

\[ t_ny(\mathbf{x}_n)\ge1-\xi_n,\qquad \xi_n\ge0.\tag{7.20} \]

We minimise

\[ C\sum_{n=1}^{N}\xi_n+\frac12\|\mathbf{w}\|^2,\tag{7.21} \]

where \(C\gt0\) trades margin against slack. Since a misclassified point has \(\xi_n\gt1\), \(\sum_n\xi_n\) bounds the number of training errors from above; \(C\to\infty\) recovers the hard-margin machine. Slack grows linearly with the distance past the margin, so the machine is still sensitive to outliers.

Derivation

Stationarity. \(\nabla_{\mathbf{w}}L\)\({}=\mathbf{w}-\sum_na_nt_n\bm\phi_n\)\({}=0\) gives (7.29); \(\partial L/\partial b=-\sum_na_nt_n=0\) gives (7.30); \(\partial L/\partial\xi_n=C-a_n-\mu_n=0\) gives (7.31).

The dual. Substitute (7.31) into (7.22): the slack terms combine as \(\xi_n(C-a_n-\mu_n)=0\), so all \(\xi\)-dependence cancels and \(L=\frac12\|\mathbf{w}\|^2-\sum_na_n(t_ny_n-1)\), exactly the separable-case expression. Then (7.29) and (7.30) reduce it to (7.32) as in Result 7.2.

The box. \(a_n\ge0\) is (7.23); \(\mu_n\ge0\) and (7.31) give \(a_n=C-\mu_n\le C\). (The dual is a maximisation of (7.32); the book's sentence above (7.33) says "minimize", which the errata correct.)

The cases. If \(a_n\gt0\), (7.25) forces \(t_ny_n=1-\xi_n\) (7.35). If \(a_n\lt C\), (7.31) gives \(\mu_n=C-a_n\gt0\), and (7.28) forces \(\xi_n=0\), so \(t_ny_n=1\). If \(a_n=C\), then \(\mu_n=0\) and (7.28) puts no restriction on \(\xi_n\): any \(\xi_n\ge0\) is allowed, subject to (7.35). If \(a_n=0\), then \(\mu_n=C\gt0\) so \(\xi_n=0\) and (7.24) gives \(t_ny_n\ge1\).

Bias. For \(n\in\mathcal{M}\), \(\xi_n=0\) and \(a_n\gt0\) so \(t_ny_n=1\), which is (7.36); multiplying by \(t_n\) and averaging as in Result 7.4 gives (7.37) (only unbounded multipliers can be used, since for \(a_n=C\) the slack is unknown). ∎

The parameter \(C\) is like an inverse regularisation constant. The ν-SVM replaces it by a parameter with a direct meaning.

Derivation

Lagrangian. With multipliers \(a_n\), \(\mu_n\), \(\delta\ge0\),

\[ L=\tfrac12\|\mathbf{w}\|^2-\nu\rho+\tfrac1N\sum_n\xi_n-\sum_na_n\{t_ny_n-\rho+\xi_n\}-\sum_n\mu_n\xi_n-\delta\rho. \]

Stationarity: \(\mathbf{w}=\sum_na_nt_n\bm\phi_n\); \(\sum_na_nt_n=0\); \(a_n+\mu_n=1/N\), so \(0\le a_n\le1/N\); and \(\partial L/\partial\rho=-\nu+\sum_na_n-\delta=0\), so \(\sum_na_n=\nu+\delta\ge\nu\), which is (7.41).

Dual. Substituting, the \(\xi\) terms cancel by \(a_n+\mu_n=1/N\), the \(\rho\) terms by the last equation, and the \(b\) term by the balance condition, leaving \(\frac12\|\mathbf{w}\|^2\)\({}-\sum_na_nt_n\mathbf{w}^\T\bm\phi_n\)\({}=-\frac12\|\mathbf{w}\|^2\), which is (7.38).

The bounds. If \(\rho\gt0\) then \(\delta=0\) (complementary slackness for \(\delta\rho=0\)), so \(\sum_na_n=\nu\). A point with \(\xi_n\gt0\) has \(\mu_n=0\) (since \(\mu_n\xi_n=0\)), hence \(a_n=1/N\); thus \(\#\{\xi_n\gt0\}/N\le\sum_na_n=\nu\). Every \(a_n\) is at most \(1/N\), so \(\nu=\sum_na_n\le\#\{a_n\gt0\}/N\): at least a fraction \(\nu\) are support vectors.

Link to \(C\). Rescale \(\tilde{\mathbf{w}}=\mathbf{w}/\rho\), \(\tilde b=b/\rho\), \(\tilde\xi_n=\xi_n/\rho\), so that \(t_n\tilde y_n\ge1-\tilde\xi_n\). The objective becomes \(\rho^2\big[\frac12\|\tilde{\mathbf{w}}\|^2+\frac1{N\rho}\sum_n\tilde\xi_n\big]\)\({}-\nu\rho\), whose bracket is the C-SVM objective (7.21) with \(C=1/(N\rho)\); the multipliers are rescaled by \(a_n^{(C)}=a_n^{(\nu)}/\rho\), so the box \(a_n^{(\nu)}\le1/N\) becomes \(a_n^{(C)}\le1/(N\rho)=C\). ∎

Two further remarks of the book. (i) The SVM outputs a decision, not a probability. Platt scaling fits a sigmoid to the machine's outputs,

\[ p(t=1\given\mathbf{x})=\sigma(Ay(\mathbf{x})+B),\tag{7.43} \]

by minimising the cross-entropy on data held out from the SVM training (otherwise the outputs are over-confident because the training points sit at \(|y|\ge1\) by construction); this assumes \(y\) is a log-odds, which SVM training does not encourage, so the fit can be a poor approximation of the posterior (Tipping, 2001). Veropoulos et al. (1999) control the false positive/negative trade-off by using different box sizes \(C_+\), \(C_-\) for the two classes. (ii) Tong and Koller (2000) show that a hybrid of generative and discriminative models gives the maximum margin hyperplane in the limit where the Parzen width \(\sigma^2\to0\), because the hyperplane is then dominated by the nearest points.

7.1.2 Relation to logistic regression

pp. 336–338 The soft-margin machine is also a regularised error minimisation, which makes the comparison with logistic regression (Section 4.3.2) transparent.

Proof

Hinge. For fixed \((\mathbf{w},b)\), the constraints (7.20) say \(\xi_n\ge\max(0,1-t_ny_n)\) and \(\xi_n\) enters (7.21) with positive coefficient \(C\), so the best \(\xi_n\) is the smallest allowed one: \(\xi_n=[1-y_nt_n]_+\) (zero when \(y_nt_n\ge1\), and \(1-y_nt_n\) otherwise). Substituting, (7.21) becomes \(C\sum_n[1-y_nt_n]_++\frac12\|\mathbf{w}\|^2\)\({}=C\big\{\sum_nE_{\rm SV}(y_nt_n)\)\({}+\frac1{2C}\|\mathbf{w}\|^2\big\}\); a positive multiple does not change the minimiser, giving (7.44).

Logistic. If \(p(t=1\given y)=\sigma(y)\) then \(p(t=-1\given y)=1-\sigma(y)=\sigma(-y)\) (symmetry of the sigmoid), so both cases are \(\sigma(yt)\). The negative log likelihood of independent data is \(-\sum_n\ln\sigma(y_nt_n)\)\({}=\sum_n\ln(1+e^{-y_nt_n})\) because \(\sigma(z)=1/(1+e^{-z})\). Adding the quadratic regulariser gives (7.47). ∎

Both are continuous (convex) upper-bound surrogates of the 0/1 error. The key difference: the hinge is exactly zero for \(yt\ge1\) (the flat region), so correctly classified points beyond the margin have zero loss and zero gradient and drop out of the solution (\(a_n=0\)); the logistic loss is never zero, so every point keeps a nonzero weight \(\sigma(-y_nt_n)\) in the gradient and the solution is not sparse. The squared error \((yt-1)^2\) punishes points that are far on the correct side, which makes no sense for classification.

7.1.3 Multiclass SVMs

pp. 338–339 The SVM is a two-class machine. The common recipes for \(K\) classes:

• One-versus-the-rest (Vapnik, 1998): \(K\) machines \(y_k\), the \(k\)th trained with class \(k\) as positive and the rest as negative; predict \(\arg\max_ky_k(\mathbf{x})\) (7.49). Problems: the \(y_k\) are trained on different tasks and need not be on a common scale; the training sets are imbalanced (with ten equal classes, 90% negatives); the plain sign rule can assign an input to several classes or none (Figure 4.2). Lee et al. (2001) set the negative target to \(-1/(K-1)\) to restore a symmetry.
• All-in-one (Weston and Watkins, 1999): one optimisation of size \((K-1)N\) instead of \(K\) problems of size \(N\); the cost goes from \(O(KN^2)\) to \(O(K^2N^2)\).
• One-versus-one: \(K(K-1)/2\) machines, one per pair, and a vote. Each machine sees only two classes' data, so the total training cost can be small, but for large \(K\) it needs many machines and ties are possible. The DAGSVM (Platt et al., 2000) arranges the pairwise machines in a directed acyclic graph so that a test point needs only \(K-1\) evaluations.
• Error-correcting output codes (Dietterich and Bakiri, 1995; Allwein et al., 2000): each class gets a binary code word, one binary machine per code bit, and a new point is assigned to the nearest code word. One-versus-one voting is a special case; the redundancy corrects some wrong outputs.

In practice one-versus-the-rest is most used, in spite of its ad hoc nature. The RVM (Section 7.2.3) has a principled multiclass form through the softmax.

7.1.4 SVMs for regression

pp. 339–343 Regularised least squares minimises \(\frac12\sum_n(y_n-t_n)^2\)\({}+\frac\lambda2\|\mathbf{w}\|^2\) (7.50), and every point contributes. To get sparsity replace the squared error by the ε-insensitive error (Figure 7.6), which is zero inside a tube of half-width \(\epsilon\gt0\) and grows linearly outside:

\[ E_\epsilon(y(\mathbf{x})-t)=\begin{cases}0,&|y(\mathbf{x})-t|\lt\epsilon;\\ |y(\mathbf{x})-t|-\epsilon,&\text{otherwise,}\end{cases}\tag{7.51} \]

and minimise, with the convention that the factor \(C\) multiplies the error,

\[ C\sum_{n=1}^{N}E_\epsilon(y(\mathbf{x}_n)-t_n)+\frac12\|\mathbf{w}\|^2.\tag{7.52} \]

As before, introduce two slack variables per point: \(\xi_n\ge0\) if \(t_n\) lies above the tube and \(\hat\xi_n\ge0\) if below (Figure 7.7), with

\[ t_n\le y(\mathbf{x}_n)+\epsilon+\xi_n,\tag{7.53} \] \[ t_n\ge y(\mathbf{x}_n)-\epsilon-\hat\xi_n,\tag{7.54} \]

and minimise

\[ C\sum_{n=1}^{N}(\xi_n+\hat\xi_n)+\frac12\|\mathbf{w}\|^2.\tag{7.55} \]

At the optimum \(\xi_n+\hat\xi_n=E_\epsilon(y_n-t_n)\) (the same argument as for the hinge: each slack takes its smallest allowed value, and a point cannot be both above and below the tube), so (7.55) is (7.52).

Derivation

Stationarity. Substitute \(y_n=\mathbf{w}^\T\bm\phi_n+b\). \(\nabla_{\mathbf{w}}L\)\({}=\mathbf{w}-\sum_na_n\bm\phi_n\)\({}+\sum_n\hat a_n\bm\phi_n\), giving (7.57). \(\partial L/\partial b=-\sum_na_n+\sum_n\hat a_n\), giving (7.58). \(\partial L/\partial\xi_n=C-\mu_n-a_n\), giving (7.59), and likewise (7.60).

The dual (Exercise 7.7). The \(\xi_n\) terms of (7.56) are \(\xi_n(C-\mu_n-a_n)=0\) and similarly for \(\hat\xi_n\), so they vanish. What remains is

\[ L=\tfrac12\|\mathbf{w}\|^2-\epsilon\sum_n(a_n+\hat a_n)+\sum_n(a_n-\hat a_n)t_n-\sum_n(a_n-\hat a_n)y_n. \]

Now \(\sum_n(a_n-\hat a_n)y_n\)\({}=\mathbf{w}^\T\sum_n(a_n-\hat a_n)\bm\phi_n\)\({}+b\sum_n(a_n-\hat a_n)\)\({}=\|\mathbf{w}\|^2+0\) by (7.57), (7.58). So \(L\)\({}=-\frac12\|\mathbf{w}\|^2\)\({}-\epsilon\sum_n(a_n+\hat a_n)\)\({}+\sum_n(a_n-\hat a_n)t_n\), and (7.57) gives \(\|\mathbf{w}\|^2\)\({}=\sum_{n,m}(a_n-\hat a_n)(a_m-\hat a_m)k_{nm}\): this is (7.61).

Constraints. \(a_n\ge0\) and \(\mu_n\ge0\) with (7.59) give \(0\le a_n\le C\), likewise (7.63); (7.58) stays. Substituting (7.57) in (7.1) gives (7.64).

KKT consequences. (7.65)–(7.66) are complementary slackness for (7.53), (7.54); (7.67)–(7.68) are \(\mu_n\xi_n=0\), \(\hat\mu_n\hat\xi_n=0\) rewritten with (7.59), (7.60). (i) If \(\xi_n\gt0\), (7.67) forces \(C-a_n=0\). (ii) If both \(a_n\gt0\) and \(\hat a_n\gt0\), (7.65) and (7.66) give \(\epsilon+\xi_n+y_n-t_n=0\) and \(\epsilon+\hat\xi_n-y_n+t_n=0\); adding gives \(2\epsilon+\xi_n+\hat\xi_n=0\), impossible since \(\epsilon\gt0\) and the slacks are \(\ge0\). (iii) \(a_n\ne0\) needs \(\epsilon+\xi_n+y_n-t_n=0\), i.e. \(t_n=y_n+\epsilon+\xi_n\ge y_n+\epsilon\): on or above the upper edge of the tube; likewise for \(\hat a_n\) below the lower edge. For \(|t_n-y_n|\lt\epsilon\) neither equation can hold with nonnegative slack, so both multipliers vanish. (iv) For \(0\lt a_n\lt C\), (7.67) gives \(\xi_n=0\), then (7.65) gives \(\epsilon+y_n-t_n=0\), and (7.1) with (7.57) gives (7.69). ∎

Derivation

Add a multiplier \(\delta\ge0\) for \(\epsilon\ge0\). The terms in \(\epsilon\) in the Lagrangian are \(C\nu\epsilon-\epsilon\sum_n(a_n+\hat a_n)\)\({}-\delta\epsilon\); stationarity in \(\epsilon\) gives \(\sum_n(a_n+\hat a_n)=C\nu-\delta\le C\nu\), which is (7.74), and removes \(\epsilon\) from the dual so that (7.61) becomes (7.70). The slack coefficient is \(C/N\), so (7.59), (7.60) become \(a_n+\mu_n=C/N\) etc., which gives (7.71), (7.72). Counts. A point outside the tube has \(\xi_n\gt0\) or \(\hat\xi_n\gt0\), hence \(a_n=C/N\) or \(\hat a_n=C/N\) as in (i) of Result 7.9; each such point contributes \(C/N\) to \(\sum_n(a_n+\hat a_n)\le\nu C\), so there are at most \(\nu N\) of them. If \(\epsilon\gt0\) then \(\delta=0\) and \(\sum_n(a_n+\hat a_n)=\nu C\); each support vector contributes at most \(C/N\) (since at most one of \(a_n,\hat a_n\) is nonzero), so at least \(\nu N\) support vectors are needed. ∎

7.1.5 Computational learning theory

pp. 344–345 SVMs were originally motivated by statistical (computational) learning theory. In the PAC (probably approximately correct) framework of Valiant (1984), the data \(\mathcal{D}\) of size \(N\) are drawn from \(p(\mathbf{x},t)\) with labels from an unknown function \(t=g(\mathbf{x})\), and a learned function \(f(\mathbf{x};\mathcal{D})\) from a class \(\mathcal{F}\) has good generalisation if

\[ \mathbb{E}_{\mathbf{x},t}\big[\mathbb{I}(f(\mathbf{x};\mathcal{D})\ne t)\big]\lt\epsilon\tag{7.75} \]

with probability at least \(1-\delta\) over \(\mathcal{D}\). The theory asks how large \(N\) must be as a function of \(\epsilon\), \(\delta\) and the complexity of \(\mathcal{F}\).

Proof

A fixed function with true error at least \(\epsilon\) classifies one random example correctly with probability at most \(1-\epsilon\), hence all \(N\) independent examples with probability at most \((1-\epsilon)^N\le e^{-\epsilon N}\). By the union bound, the probability that some function in \(\mathcal{F}\) with error at least \(\epsilon\) is consistent with the data is at most \(|\mathcal{F}|e^{-\epsilon N}\). Requiring this to be at most \(\delta\) gives the condition. ∎

For infinite classes \(\ln|\mathcal{F}|\) is replaced by the Vapnik-Chervonenkis (VC) dimension (Vapnik and Chervonenkis, 1971), the size of the largest set of points that \(\mathcal{F}\) can label in all possible ways ("shatter"). For hyperplanes in \(\mathbb{R}^D\) it is \(D+1\); for large-margin hyperplanes, Vapnik's bound replaces the dimension of feature space by \(\min(D,\lceil R^2/\rho^2\rceil)\) (up to a constant), where \(R\) is the radius of a ball containing the data and \(\rho\) the margin. This is the learning-theory reason to maximise the margin: a large margin lowers the effective capacity, whatever the dimension of the feature space. I state this bound without proof. PAC bounds are worst case: they hold for every distribution and every function in the class, hence are very conservative and have found few practical applications. The PAC-Bayesian framework (McAllester, 2003) puts a distribution over \(\mathcal{F}\), like a prior, and gives tighter but still conservative bounds.

7.2 Relevance vector machines

pp. 345–356 The SVM gives decisions, not posterior probabilities; its extension to \(K\gt2\) classes is awkward; \(C\) and \(\epsilon\) (or \(\nu\)) need cross-validation; and the kernels must be positive definite and centred on the data. The relevance vector machine (Tipping, 2001) is a Bayesian linear model that keeps the sparse, kernel-expansion form but removes these limits, and typically uses far fewer basis functions than the SVM. As in the book, regression comes first.

7.2.1 RVM for regression

pp. 345–349 The model is the linear regression of Chapter 3,

\[ p(t\given\mathbf{x},\mathbf{w},\beta)=\N(t\given y(\mathbf{x}),\beta^{-1}),\qquad y(\mathbf{x})=\sum_{i=1}^{M}w_i\phi_i(\mathbf{x})=\mathbf{w}^\T\bm\phi(\mathbf{x}),\tag{7.76, 7.77} \]

with noise precision \(\beta=\sigma^{-2}\). For the RVM proper the basis functions are kernels centred on the training points,

\[ y(\mathbf{x})=\sum_{n=1}^{N}w_nk(\mathbf{x},\mathbf{x}_n)+b,\tag{7.78} \]

so \(M=N+1\) and \(\bm\Phi=\mathbf{K}\) (the \((N+1)\times(N+1)\) kernel matrix with a column of ones for the bias). The analysis below holds for arbitrary basis functions: no positive definiteness is needed and the basis functions need not be tied to the data. With a data matrix \(\mathbf{X}\) and targets \(\mathbf{t}\) the likelihood is

\[ p(\mathbf{t}\given\mathbf{X},\mathbf{w},\beta)=\prod_{n=1}^{N}p(t_n\given\mathbf{x}_n,\mathbf{w},\beta).\tag{7.79} \]

Erratum. The book prints \(p(t_n\given\mathbf{x}_n,\mathbf{w},\beta^{-1})\) on the right of (7.79); the argument is the precision \(\beta\), as in (7.76).

The key difference from Chapter 3 is the ARD prior: one precision per weight,

\[ p(\mathbf{w}\given\bm\alpha)=\prod_{i=1}^{M}\N(w_i\given0,\alpha_i^{-1}).\tag{7.80} \]

Maximising the evidence over \(\bm\alpha\) will send many \(\alpha_i\) to infinity; the corresponding weight then has a posterior concentrated at zero and its basis function is pruned.

Proof

By Bayes, \(\ln p(\mathbf{w}\given\mathbf{t})\)\({}=\ln p(\mathbf{t}\given\mathbf{w})\)\({}+\ln p(\mathbf{w})+\const\), and as a function of \(\mathbf{w}\)

\[ -\frac\beta2\|\mathbf{t}-\bm\Phi\mathbf{w}\|^2-\frac12\mathbf{w}^\T\mathbf{A}\mathbf{w}=-\frac12\mathbf{w}^\T(\mathbf{A}+\beta\bm\Phi^\T\bm\Phi)\mathbf{w}+\beta\mathbf{w}^\T\bm\Phi^\T\mathbf{t}+\const. \]

Compare with the Gaussian form \(-\frac12\mathbf{w}^\T\bm\Sigma^{-1}\mathbf{w}\)\({}+\mathbf{w}^\T\bm\Sigma^{-1}\mathbf{m}+\const\) (completing the square): the quadratic coefficient gives \(\bm\Sigma^{-1}=\mathbf{A}+\beta\bm\Phi^\T\bm\Phi\), (7.83), and the linear one \(\bm\Sigma^{-1}\mathbf{m}\)\({}=\beta\bm\Phi^\T\mathbf{t}\), (7.82). The matrix \(\mathbf{A}+\beta\bm\Phi^\T\bm\Phi\) is positive definite when all \(\alpha_i\gt0\), so the Gaussian is proper. This is also the linear-Gaussian result of the toolkit, and it is (3.49) with \(\alpha\mathbf{I}\) replaced by \(\mathbf{A}\). ∎

The hyperparameters are set by type-2 maximum likelihood (the evidence approximation, Section 3.5): maximise the marginal likelihood

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

Proof (two routes)

Route 1: linear-Gaussian rule (Exercise 7.11). Write \(\mathbf{t}=\bm\Phi\mathbf{w}+\bm\varepsilon\) with \(\mathbf{w}\sim\N(\mathbf{0},\mathbf{A}^{-1})\), \(\bm\varepsilon\)\({}\sim\N(\mathbf{0},\beta^{-1}\mathbf{I})\) independent. This is the setting of (2.113)–(2.115) (toolkit): the marginal of \(\mathbf{t}\) is Gaussian with mean \(\bm\Phi\cdot\mathbf{0}=\mathbf{0}\) and covariance \(\beta^{-1}\mathbf{I}\)\({}+\bm\Phi\mathbf{A}^{-1}\bm\Phi^\T\)\({}=\mathbf{C}\). That is (7.85), (7.86).

Route 2: complete the square (Exercise 7.10). The integrand is \(\big(\tfrac\beta{2\pi}\big)^{N/2}\prod_i\big(\tfrac{\alpha_i}{2\pi}\big)^{1/2}\exp\{-E(\mathbf{w})\}\) with \(E(\mathbf{w})\)\({}=\frac\beta2\|\mathbf{t}-\bm\Phi\mathbf{w}\|^2\)\({}+\frac12\mathbf{w}^\T\mathbf{A}\mathbf{w}\). Expanding and using (7.82), (7.83),

\[ E(\mathbf{w})=\tfrac12(\mathbf{w}-\mathbf{m})^\T\bm\Sigma^{-1}(\mathbf{w}-\mathbf{m})+E(\mathbf{m}),\qquad E(\mathbf{m})=\tfrac\beta2\|\mathbf{t}-\bm\Phi\mathbf{m}\|^2+\tfrac12\mathbf{m}^\T\mathbf{A}\mathbf{m}. \]

(Check: the quadratic terms match by (7.83); the linear term of \(E\) is \(-\beta\mathbf{w}^\T\bm\Phi^\T\mathbf{t}\), that of the right-hand side is \(-\mathbf{w}^\T\bm\Sigma^{-1}\mathbf{m}\), equal by (7.82); the constant of \(E\) is \(\frac\beta2\mathbf{t}^\T\mathbf{t}\), and \(E(\mathbf{m})\)\({}=\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\beta\mathbf{m}^\T\bm\Phi^\T\mathbf{t}\)\({}+\frac12\mathbf{m}^\T\bm\Sigma^{-1}\mathbf{m}\)\({}=\frac\beta2\mathbf{t}^\T\mathbf{t}\)\({}-\frac12\mathbf{m}^\T\bm\Sigma^{-1}\mathbf{m}\), using \(\bm\Sigma^{-1}\mathbf{m}\)\({}=\beta\bm\Phi^\T\mathbf{t}\) twice.) The integral over \(\mathbf{w}\) is a Gaussian integral, \((2\pi)^{M/2}|\bm\Sigma|^{1/2}\); with \(\prod_i\alpha_i=|\mathbf{A}|\) this gives \((*)\).

The two forms agree. The determinant lemma gives \(|\mathbf{C}|\)\({}=|\beta^{-1}\mathbf{I}|\,|\mathbf{A}^{-1}|\,|\mathbf{A}\)\({}+\beta\bm\Phi^\T\bm\Phi|\)\({}=\beta^{-N}|\mathbf{A}|^{-1}|\bm\Sigma|^{-1}\), so \(-\frac12\ln|\mathbf{C}|\)\({}=\frac N2\ln\beta+\frac12\ln|\mathbf{A}|\)\({}+\frac12\ln|\bm\Sigma|\). The Woodbury identity gives \(\mathbf{C}^{-1}\)\({}=\beta\mathbf{I}\)\({}-\beta^2\bm\Phi\bm\Sigma\bm\Phi^\T\), so \(\mathbf{t}^\T\mathbf{C}^{-1}\mathbf{t}\)\({}=\beta\mathbf{t}^\T\mathbf{t}\)\({}-\mathbf{m}^\T\bm\Sigma^{-1}\mathbf{m}\)\({}=2E(\mathbf{m})\). Substituting both into (7.85) gives \((*)\). ∎

Derivation

Differentiate the form \((*)\) of Result 7.13. The term \(E(\mathbf{m})\) depends on the hyperparameters directly and through \(\mathbf{m}\); but \(\mathbf{m}\) minimises \(E(\mathbf{w})\) over \(\mathbf{w}\), so \(\partial E/\partial\mathbf{w}=\mathbf{0}\) at \(\mathbf{m}\), and the indirect dependence drops out: only the explicit derivatives matter.

For \(\alpha_i\). \(\partial\frac12\ln|\mathbf{A}|/\partial\alpha_i\)\({}=\frac1{2\alpha_i}\); \(\partial\frac12\mathbf{m}^\T\mathbf{A}\mathbf{m}/\partial\alpha_i|_{\mathbf{m}}\)\({}=\frac12m_i^2\); and \(\ln|\bm\Sigma|\)\({}=-\ln|\mathbf{A}+\beta\bm\Phi^\T\bm\Phi|\), so by \(\partial\ln|\mathbf{X}|\)\({}=\Tr(\mathbf{X}^{-1}\partial\mathbf{X})\) (matrix calculus), \(\partial\frac12\ln|\bm\Sigma|/\partial\alpha_i\)\({}=-\frac12\Tr(\bm\Sigma\,\mathbf{e}_i\mathbf{e}_i^\T)\)\({}=-\frac12\Sigma_{ii}\). Hence

\[ \frac{\partial\ln p}{\partial\alpha_i}=\frac1{2\alpha_i}-\frac12\Sigma_{ii}-\frac12m_i^2=0\iff\alpha_im_i^2=1-\alpha_i\Sigma_{ii}=\gamma_i, \]

after multiplying by \(2\alpha_i\); this is (7.87).

For \(\beta\). \(\partial\frac N2\ln\beta/\partial\beta\)\({}=\frac N{2\beta}\); \(\partial(-\frac\beta2\|\mathbf{t}-\bm\Phi\mathbf{m}\|^2)/\partial\beta|_{\mathbf{m}}\)\({}=-\frac12\|\mathbf{t}-\bm\Phi\mathbf{m}\|^2\); \(\partial\frac12\ln|\bm\Sigma|/\partial\beta\)\({}=-\frac12\Tr(\bm\Sigma\bm\Phi^\T\bm\Phi)\). Since \(\beta\bm\Phi^\T\bm\Phi=\bm\Sigma^{-1}-\mathbf{A}\), \(\beta\Tr(\bm\Sigma\bm\Phi^\T\bm\Phi)\)\({}=\Tr(\mathbf{I}-\bm\Sigma\mathbf{A})\)\({}=M-\sum_i\alpha_i\Sigma_{ii}\)\({}=\sum_i\gamma_i\). Setting the derivative to zero and multiplying by \(2\beta\): \(N-\beta\|\mathbf{t}-\bm\Phi\mathbf{m}\|^2\)\({}-\sum_i\gamma_i\)\({}=0\), which is (7.88).

Range of \(\gamma_i\). \(\bm\Sigma^{-1}\)\({}=\mathbf{A}\)\({}+\beta\bm\Phi^\T\bm\Phi\succeq\mathbf{A}\), so \(\bm\Sigma\preceq\mathbf{A}^{-1}\) and \(\Sigma_{ii}\le1/\alpha_i\), i.e. \(\gamma_i\ge0\); and \(\Sigma_{ii}\gt0\) gives \(\gamma_i\lt1\). Interpretation: \(\gamma_i\) measures how well \(w_i\) is determined by the data (prior variance \(1/\alpha_i\) against posterior variance \(\Sigma_{ii}\)); \(\sum_i\gamma_i\) is the effective number of well-determined parameters (Section 3.5.3). In (7.87) the new precision is large when the posterior mean \(m_i\) is small compared to what the data determine, which pushes \(\alpha_i\) up further: a runaway for irrelevant weights. ∎

Proof

Given \(\mathbf{x}\), \(t=\mathbf{w}^\T\bm\phi(\mathbf{x})+\varepsilon\) with \(\varepsilon\sim\N(0,(\beta^\star)^{-1})\) independent of \(\mathbf{w}\sim\N(\mathbf{m},\bm\Sigma)\). A linear function of Gaussians is Gaussian (toolkit, (2.115)): mean \(\mathbf{m}^\T\bm\phi\) and variance \(\bm\phi^\T\bm\Sigma\bm\phi+(\beta^\star)^{-1}\), as variances of independent terms add. ∎

The first term of (7.91) is noise, the second is uncertainty in \(\mathbf{w}\). With localised basis functions centred on data (the RVM), \(\bm\phi(\mathbf{x})\to\mathbf{0}\) far from the data, so \(\sigma^2(\mathbf{x})\to(\beta^\star)^{-1}\): the model becomes over-confident away from the data, unlike a Gaussian process, whose variance returns to the prior variance (Rasmussen and Quiñonero-Candela, 2005). Training costs \(O(M^3)\) for the \(M\times M\) inversion (\(M=N+1\) at the start), more than SVM training (about \(O(N^2)\)), but one RVM run replaces the cross-validation over \(C,\epsilon\) needed by the SVM, and typically leaves an order of magnitude fewer basis functions, so prediction is much cheaper.

7.2.2 Analysis of sparsity

pp. 349–353 Why does maximising the evidence prune basis functions? First a geometric view, then an exact analysis that also yields a faster algorithm.

Proof

\(\mathbf{C}\) has eigenvalue \(v=\beta^{-1}+\|\bm\varphi\|^2/\alpha\) along \(\bm\varphi\) and \(\beta^{-1}\) in the \(N-1\) orthogonal directions. Hence \(\ln|\mathbf{C}|=(N-1)\ln\beta^{-1}+\ln v\) and \(\mathbf{t}^\T\mathbf{C}^{-1}\mathbf{t}\)\({}=\beta\|\mathbf{t}_\perp\|^2+t_\parallel^2/v\). Only the part \(f(v)=\ln v+t_\parallel^2/v\) depends on \(\alpha\), and the log marginal likelihood is \(-\frac12f(v)+\const\). Now \(f'(v)=1/v-t_\parallel^2/v^2\) is negative for \(v\lt t_\parallel^2\) and positive for \(v\gt t_\parallel^2\), so \(f\) has its minimum at \(v=t_\parallel^2\). The admissible values are \(v\in(\beta^{-1},\infty)\) (as \(\alpha\) runs over \((\infty,0)\), \(\alpha=\infty\) being \(v=\beta^{-1}\)). If \(t_\parallel^2\gt\beta^{-1}\), the unconstrained minimiser is admissible, and \(\|\bm\varphi\|^2/\alpha=t_\parallel^2-\beta^{-1}\) gives \(\alpha^\star\). If \(t_\parallel^2\le\beta^{-1}\), \(f\) is increasing on the whole admissible range, so the best value is the endpoint \(v=\beta^{-1}\), i.e. \(\alpha=\infty\). ∎

For several basis vectors the same intuition holds. To make it exact, isolate the contribution of one hyperparameter. In (7.87) the right-hand side depends on \(\alpha_i\) itself (through \(m_i\) and \(\Sigma_{ii}\)), so it is only an implicit equation; instead we make the dependence on \(\alpha_i\) explicit (Faul and Tipping, 2002; Tipping and Faul, 2003). Write \(\bm\varphi_i\) for the \(i\)th column of \(\bm\Phi\) (the \(N\)-vector \((\phi_i(\mathbf{x}_1),\dots,\phi_i(\mathbf{x}_N))^\T\)), as opposed to \(\bm\phi_n\), the \(n\)th row.

Proof

(7.93): in (7.86), \(\bm\Phi\mathbf{A}^{-1}\bm\Phi^\T\)\({}=\sum_j\alpha_j^{-1}\bm\varphi_j\bm\varphi_j^\T\). (7.94) is the matrix determinant lemma \(|\mathbf{M}+\mathbf{u}\mathbf{v}^\T|\)\({}=|\mathbf{M}|(1+\mathbf{v}^\T\mathbf{M}^{-1}\mathbf{u})\) with \(\mathbf{M}=\mathbf{C}_{-i}\), \(\mathbf{u}=\alpha_i^{-1}\bm\varphi_i\), \(\mathbf{v}=\bm\varphi_i\) (the Sherman-Morrison/Woodbury case with a rank-one update, toolkit). (7.95) is the Sherman-Morrison formula \((\mathbf{M}+\mathbf{u}\mathbf{v}^\T)^{-1}\)\({}=\mathbf{M}^{-1}\)\({}-\mathbf{M}^{-1}\mathbf{u}\mathbf{v}^\T\mathbf{M}^{-1}/(1+\mathbf{v}^\T\mathbf{M}^{-1}\mathbf{u})\); here \(\alpha_i^{-1}/(1+\alpha_i^{-1}s_i)\)\({}=1/(\alpha_i+s_i)\).

Now take logarithms. From (7.94), with \(1+s_i/\alpha_i=(\alpha_i+s_i)/\alpha_i\),

\[ \ln|\mathbf{C}|=\ln|\mathbf{C}_{-i}|+\ln(\alpha_i+s_i)-\ln\alpha_i. \]

From (7.95), \(\mathbf{t}^\T\mathbf{C}^{-1}\mathbf{t}\)\({}=\mathbf{t}^\T\mathbf{C}_{-i}^{-1}\mathbf{t}\)\({}-q_i^2/(\alpha_i+s_i)\). Insert both into (7.85): the \(\alpha_i\)-independent parts form \(\mathcal{L}(\bm\alpha_{-i})\) and the rest is \(-\frac12\{\ln(\alpha_i+s_i)-\ln\alpha_i\)\({}-q_i^2/(\alpha_i+s_i)\}\)\({}=\lambda(\alpha_i)\). ∎

Interpretation: \(s_i\) measures how much \(\bm\varphi_i\) overlaps with the basis vectors already in the model (in the metric \(\mathbf{C}_{-i}^{-1}\)); \(q_i\) measures how well \(\bm\varphi_i\) is aligned with the error between \(\mathbf{t}\) and the prediction \(\mathbf{y}_{-i}\) of the model without \(\bm\varphi_i\).

Proof

(7.100). \(\lambda'\)\({}=\frac12\big[\frac1{\alpha}-\frac1{\alpha+s}-\frac{q^2}{(\alpha+s)^2}\big]\). Put over a common denominator: \(\frac1\alpha-\frac1{\alpha+s}\)\({}=\frac{s}{\alpha(\alpha+s)}\), so \(\lambda'\)\({}=\frac12\cdot\frac{s(\alpha+s)-\alpha q^2}{\alpha(\alpha+s)^2}\)\({}=\frac{s^2/\alpha+s-q^2}{2(\alpha+s)^2}\), which is (7.100).

Signs. The denominator is positive, so the sign of \(\lambda'\) is that of \(h(\alpha)=s^2/\alpha-(q^2-s)\), a strictly decreasing function of \(\alpha\gt0\) going from \(+\infty\) at \(0^+\) to \(-(q^2-s)\) at \(\infty\). (a) If \(q^2\gt s\), \(h\) changes sign exactly once, at \(s^2/\alpha=q^2-s\), i.e. (7.101). \(\lambda'\gt0\) before and \(\lt0\) after, so this is the global maximum. (This is the "second derivative" check of Exercise 7.16: at the zero of \(h\), \(\lambda''\)\({}=h'(\alpha)/(2(\alpha+s)^2)\)\({}=-s^2/(2\alpha^2(\alpha+s)^2)\)\({}\lt0\).) (b) If \(q^2\le s\), \(h(\alpha)\gt s^2/\alpha\gt0\) always, so \(\lambda\) increases and its supremum is the limit \(\alpha\to\infty\), where \(\lambda\)\({}\to\frac12[\ln\frac\alpha{\alpha+s}+\frac{q^2}{\alpha+s}]\)\({}\to0\).

(c) At (7.101), \(\alpha+s=sq^2/(q^2-s)\), so \(\alpha/(\alpha+s)=s/q^2\) and \(q^2/(\alpha+s)=(q^2-s)/s\). Therefore \(\lambda(\alpha^\star)\)\({}=\frac12\{\ln\frac s{q^2}+\frac{q^2}s-1\}\)\({}=\frac12(x-1-\ln x)\), which is positive for \(x\gt1\) (since \(\ln x\lt x-1\)). ∎

Compare with Result 7.16: for one basis vector, \(\mathbf{C}_{-i}=\beta^{-1}\mathbf{I}\), so \(s=\beta\|\bm\varphi\|^2\), \(q=\beta\bm\varphi^\T\mathbf{t}\), and the condition \(q^2\gt s\) is \(\beta(\bm\varphi^\T\mathbf{t})^2\)\({}\gt\|\bm\varphi\|^2\), i.e. \(\beta t_\parallel^2\gt1\); and (7.101) becomes \(\alpha^\star\)\({}=\|\bm\varphi\|^2/(t_\parallel^2-\beta^{-1})\). The two analyses agree.

Proof

(7.106), (7.107). By the Woodbury identity applied to (7.86), \(\mathbf{C}^{-1}\)\({}=\beta\mathbf{I}\)\({}-\beta^2\bm\Phi(\mathbf{A}+\beta\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\)\({}=\beta\mathbf{I}\)\({}-\beta^2\bm\Phi\bm\Sigma\bm\Phi^\T\) by (7.83). Multiply on the left by \(\bm\varphi_i^\T\) and on the right by \(\mathbf{t}\) or \(\bm\varphi_i\). This costs \(O(M^2)\) per \(i\) given \(\bm\Sigma\), instead of an \(N\times N\) inverse.

(7.104), (7.105). Multiply (7.95) by \(\bm\varphi_i^\T\) on the left and \(\mathbf{t}\) (or \(\bm\varphi_i\)) on the right. Using the definitions (7.98), (7.99):

\[ Q_i=q_i-\frac{s_iq_i}{\alpha_i+s_i}=\frac{\alpha_iq_i}{\alpha_i+s_i},\qquad S_i=s_i-\frac{s_i^2}{\alpha_i+s_i}=\frac{\alpha_is_i}{\alpha_i+s_i}. \]

Solve for \(s_i\): \(S_i(\alpha_i+s_i)\)\({}=\alpha_is_i\)\({}\Rightarrow S_i\alpha_i\)\({}=s_i(\alpha_i-S_i)\)\({}\Rightarrow s_i\)\({}=\alpha_iS_i/(\alpha_i-S_i)\). Then \(q_i\)\({}=Q_i(\alpha_i+s_i)/\alpha_i\)\({}=Q_i\big(1+\tfrac{S_i}{\alpha_i-S_i}\big)\)\({}=\alpha_iQ_i/(\alpha_i-S_i)\). For \(\alpha_i\to\infty\) the factors tend to 1. ∎

7.2.3 RVM for classification

pp. 353–356 For two classes with \(t\in\{0,1\}\) (the book's convention for the RVM) the model is logistic regression with the ARD prior (7.80):

\[ y(\mathbf{x},\mathbf{w})=\sigma\big(\mathbf{w}^\T\bm\phi(\mathbf{x})\big).\tag{7.108} \]

The integral over \(\mathbf{w}\) is no longer analytic, so use the Laplace approximation (Section 4.4) at the current \(\bm\alpha\), maximise the resulting approximate marginal likelihood over \(\bm\alpha\), and repeat. The log posterior of \(\mathbf{w}\) is

\[ \ln p(\mathbf{w}\given\mathbf{t},\bm\alpha)=\ln\{p(\mathbf{t}\given\mathbf{w})p(\mathbf{w}\given\bm\alpha)\}-\ln p(\mathbf{t}\given\bm\alpha)=\sum_{n=1}^{N}\{t_n\ln y_n+(1-t_n)\ln(1-y_n)\}-\frac12\mathbf{w}^\T\mathbf{A}\mathbf{w}+\const.\tag{7.109} \]

Proof

With \(a_n=\mathbf{w}^\T\bm\phi_n\) and \(y_n=\sigma(a_n)\), \(\sigma'(a)=\sigma(1-\sigma)\) (4.88). The data term has \(\frac{\partial}{\partial a_n}\{t_n\ln y_n\)\({}+(1-t_n)\ln(1-y_n)\}\)\({}=\frac{t_n}{y_n}y_n(1-y_n)\)\({}-\frac{1-t_n}{1-y_n}y_n(1-y_n)\)\({}=t_n(1-y_n)-(1-t_n)y_n\)\({}=t_n-y_n\). By the chain rule (\(\partial a_n/\partial\mathbf{w}=\bm\phi_n\)) the gradient is \(\sum_n(t_n-y_n)\bm\phi_n\)\({}=\bm\Phi^\T(\mathbf{t}-\mathbf{y})\); the prior term \(-\frac12\mathbf{w}^\T\mathbf{A}\mathbf{w}\) has gradient \(-\mathbf{A}\mathbf{w}\). This is (7.110). Differentiating again, \(\partial(t_n-y_n)/\partial\mathbf{w}\)\({}=-y_n(1-y_n)\bm\phi_n^\T\), so the data term contributes \(-\sum_nb_n\bm\phi_n\bm\phi_n^\T\)\({}=-\bm\Phi^\T\mathbf{B}\bm\Phi\), and the prior \(-\mathbf{A}\): (7.111). \(\mathbf{B}\succ0\) because \(0\lt y_n\lt1\), and \(\mathbf{A}\succ0\), so the Hessian is negative definite. Setting (7.110) to zero gives (7.112); the Laplace approximation takes the covariance to be the inverse of the negative Hessian, (7.113). ∎

The Laplace approximation (4.135) of the marginal likelihood is

\[ p(\mathbf{t}\given\bm\alpha)=\int p(\mathbf{t}\given\mathbf{w})p(\mathbf{w}\given\bm\alpha)\,\d\mathbf{w}\simeq p(\mathbf{t}\given\mathbf{w}^\star)\,p(\mathbf{w}^\star\given\bm\alpha)\,(2\pi)^{M/2}|\bm\Sigma|^{1/2}.\tag{7.114} \]

Derivation

(7.115). Take the log of (7.114): \(\ln p(\mathbf{t}\given\bm\alpha)\)\({}\simeq\ln p(\mathbf{t}\given\mathbf{w}^\star)\)\({}+\ln p(\mathbf{w}^\star\given\bm\alpha)\)\({}+\frac M2\ln2\pi+\frac12\ln|\bm\Sigma|\), with \(\ln p(\mathbf{w}^\star\given\bm\alpha)\)\({}=\frac12\sum_j\ln\alpha_j\)\({}-\frac12\mathbf{w}^{\star\T}\mathbf{A}\mathbf{w}^\star\)\({}-\frac M2\ln2\pi\). Differentiate in \(\alpha_i\). The mode \(\mathbf{w}^\star\) depends on \(\alpha_i\), but the function \(\ln p(\mathbf{t}\given\mathbf{w})\)\({}+\ln p(\mathbf{w}\given\bm\alpha)\) is stationary in \(\mathbf{w}\) at \(\mathbf{w}^\star\) (that is (7.112)), so that indirect dependence drops out of the first two terms. This leaves \(\frac1{2\alpha_i}-\frac12(w_i^\star)^2\) from the prior. For the last term, \(\ln|\bm\Sigma|\)\({}=-\ln|\bm\Phi^\T\mathbf{B}\bm\Phi+\mathbf{A}|\), so with \(\mathbf{B}\) held fixed (the book ignores the dependence of \(\mathbf{B}\) on \(\mathbf{w}^\star\), a standard approximation) \(\partial\frac12\ln|\bm\Sigma|/\partial\alpha_i\)\({}=-\frac12\Sigma_{ii}\). Adding the three pieces gives (7.115).

(7.116). Multiply (7.115) by \(2\alpha_i\): \(-\alpha_i(w_i^\star)^2+1-\alpha_i\Sigma_{ii}=0\), i.e. \(\alpha_i(w_i^\star)^2=\gamma_i\).

(7.117)–(7.119). Consider the Gaussian regression problem \(\hat{\mathbf{t}}\)\({}=\bm\Phi\mathbf{w}+\bm\varepsilon\), \(\bm\varepsilon\sim\N(\mathbf{0},\mathbf{B}^{-1})\), \(\mathbf{w}\sim\N(\mathbf{0},\mathbf{A}^{-1})\). Its posterior has precision \(\mathbf{A}+\bm\Phi^\T\mathbf{B}\bm\Phi\)\({}=\bm\Sigma^{-1}\) (Result 7.12 with \(\beta\mathbf{I}\to\mathbf{B}\)) and mean \(\bm\Sigma\bm\Phi^\T\mathbf{B}\hat{\mathbf{t}}\). By (7.117), \(\bm\Phi^\T\mathbf{B}\hat{\mathbf{t}}\)\({}=\bm\Phi^\T\mathbf{B}\bm\Phi\mathbf{w}^\star\)\({}+\bm\Phi^\T(\mathbf{t}-\mathbf{y})\), and by (7.112) \(\bm\Phi^\T(\mathbf{t}-\mathbf{y})\)\({}=\mathbf{A}\mathbf{w}^\star\); hence the mean is \(\bm\Sigma(\bm\Phi^\T\mathbf{B}\bm\Phi+\mathbf{A})\mathbf{w}^\star\)\({}=\mathbf{w}^\star\). So this regression problem has exactly the Laplace posterior \(\N(\mathbf{w}^\star,\bm\Sigma)\). Its log marginal likelihood is, by Result 7.13 with \(\beta\mathbf{I}\to\mathbf{B}\), the expression (7.118) with \(\mathbf{C}\) as in (7.119), and equals \(\frac12\ln|\mathbf{B}|+\frac12\ln|\mathbf{A}|\)\({}+\frac12\ln|\bm\Sigma|\)\({}-\frac12(\hat{\mathbf{t}}-\bm\Phi\mathbf{w}^\star)^\T\mathbf{B}(\hat{\mathbf{t}}-\bm\Phi\mathbf{w}^\star)\)\({}-\frac12\mathbf{w}^{\star\T}\mathbf{A}\mathbf{w}^\star\)\({}+\const\). The \(\bm\alpha\)-gradient of that expression (again the indirect dependence through \(\mathbf{w}^\star\) vanishes, as \(\partial/\partial\mathbf{w}\) of the last two terms is \(\bm\Phi^\T\mathbf{B}(\hat{\mathbf{t}}-\bm\Phi\mathbf{w}^\star)\)\({}-\mathbf{A}\mathbf{w}^\star\)\({}=\bm\Phi^\T(\mathbf{t}-\mathbf{y})\)\({}-\mathbf{A}\mathbf{w}^\star\)\({}=0\)) is again \(\frac1{2\alpha_i}-\frac12\Sigma_{ii}\)\({}-\frac12(w_i^\star)^2\), the same as (7.115). The two approximate evidences therefore have the same stationary points in \(\bm\alpha\). ∎

Why do the relevance vectors avoid the boundary? A basis function centred on a point near the decision boundary has a vector \(\bm\varphi_i\) that is poorly aligned with the target vector (the kernel is large for points of both classes), so \(q_i^2\) is small compared with \(s_i\) and it is pruned (Result 7.18). An SVM, in contrast, keeps exactly the points at the boundary.

More than two classes. Use \(K\) linear models \(a_k=\mathbf{w}_k^\T\mathbf{x}\) (7.120), a softmax (7.121), \(y_k=\exp(a_k)/\sum_j\exp(a_j)\), and the log likelihood with 1-of-\(K\) targets \(t_{nk}\),

\[ \ln p(\mathbf{T}\given\mathbf{w}_1,\dots,\mathbf{w}_K)=\sum_{n=1}^{N}\sum_{k=1}^{K}t_{nk}\ln y_{nk}\tag{7.122} \]

(The book prints the product, i.e. \(p\) and not \(\ln p\); the logarithm is the sum.) The Laplace approximation and hyperparameter optimisation proceed as before using IRLS (Tipping, 2001). This is a more principled multiclass method than the pairwise SVM schemes and gives probabilities, but the Hessian is \(MK\times MK\), so training is a factor \(K^3\) more costly than the two-class RVM.

Exercises

All 19 exercises of Chapter 7 (Bishop, 2006, pp. 357–358), restated in my own words, with complete solutions (hidden). Stars give the book's difficulty and "www" marks the exercises with official solutions. Most exercises are derivations used above; their solutions give the argument and link to the full version where it already appears. Solution manuals contain slips: the closest-mean remark in Exercise 7.1 is one.

1. Exercise 7.1 (★★, www). Model the inputs of each class (labels \(t_n\in\{-1,1\}\)) with a Parzen density estimator with kernel \(k(\mathbf{x},\mathbf{x}')\). Write the minimum misclassification-rate rule for equal priors. Show that for \(k(\mathbf{x},\mathbf{x}')\)\({}=\mathbf{x}^\T\mathbf{x}'\) it reduces to comparing inner products with the class means, and that for \(k=\bm\phi(\mathbf{x})^\T\bm\phi(\mathbf{x}')\) the same holds in feature space.

   Solution

   Let \(N_\pm\) be the class sizes. The Parzen estimates are \(p(\mathbf{x}\given t=\pm1)\)\({}=\frac1{N_\pm}\sum_{n:t_n=\pm1}\frac1{Z_k}k(\mathbf{x},\mathbf{x}_n)\), with \(Z_k\) the normalising constant of the kernel (the same for both classes). With equal priors, \(p(t\given\mathbf{x})\)\({}\propto p(\mathbf{x}\given t)\), and the misclassification rate is minimised by choosing the class with the larger posterior, so

   \[ \hat t(\mathbf{x})=\operatorname{sign}\Big(\frac1{N_+}\sum_{n:t_n=+1}k(\mathbf{x},\mathbf{x}_n)-\frac1{N_-}\sum_{n:t_n=-1}k(\mathbf{x},\mathbf{x}_n)\Big)=\operatorname{sign}\Big(\sum_n\frac{t_n}{N_{t_n}}k(\mathbf{x},\mathbf{x}_n)\Big). \]

   For \(k(\mathbf{x},\mathbf{x}')\)\({}=\mathbf{x}^\T\mathbf{x}'\) (improper as a density, but the comparison is still meaningful), \(\frac1{N_\pm}\sum_{\pm}\mathbf{x}^\T\mathbf{x}_n\)\({}=\mathbf{x}^\T\bar{\mathbf{x}}_\pm\) with the class mean \(\bar{\mathbf{x}}_\pm\), so \(\hat t\)\({}=\operatorname{sign}\{\mathbf{x}^\T(\bar{\mathbf{x}}_+-\bar{\mathbf{x}}_-)\}\): assign to the class whose mean has the larger inner product with \(\mathbf{x}\). The same algebra with \(k=\bm\phi^\T\bm\phi\) gives \(\hat t\)\({}=\operatorname{sign}\{\bm\phi(\mathbf{x})^\T(\bar{\bm\phi}_+-\bar{\bm\phi}_-)\}\), with the means of the feature vectors.

   Caveat. This is "closest mean" in the sense of the larger inner product. Exactly closest in Euclidean distance would be \(\|\mathbf{x}-\bar{\mathbf{x}}_+\|^2\)\({}\lt\|\mathbf{x}-\bar{\mathbf{x}}_-\|^2\)\({}\iff\mathbf{x}^\T(\bar{\mathbf{x}}_+-\bar{\mathbf{x}}_-)\)\({}\gt\frac12(\|\bar{\mathbf{x}}_+\|^2-\|\bar{\mathbf{x}}_-\|^2)\), which differs by a threshold unless the two mean vectors have equal norm. Using a normalised kernel such as the Gaussian removes the issue in the Parzen setting.

2. Exercise 7.2 (★). Show that replacing the 1 in the constraint (7.5) by any \(\gamma\gt0\) leaves the maximum margin hyperplane unchanged.

   Solution

   Result 7.1(c): feasible set of the \(\gamma\)-problem is \(\gamma\) times that of (7.5), the objective scales by \(\gamma^2\), so the minimiser is \(\gamma(\mathbf{w},b)\), the same hyperplane \(\{y=0\}\) with margin \(\gamma/\|\gamma\mathbf{w}\|=1/\|\mathbf{w}\|\) unchanged. (This is the same as the freedom \(\kappa\) in (7.4).)

3. Exercise 7.3 (★★). Show that, whatever the dimension of the data space, a data set of just two points, one from each class, determines the maximum margin hyperplane.

   Solution

   Let the feature vectors be \(\bm\phi_+\ne\bm\phi_-\), \(d=\|\bm\phi_+-\bm\phi_-\|\). The balance constraint (7.12) is \(a_+-a_-=0\), so \(a_+=a_-=a\ge0\). Since \(t_+t_-=-1\), the dual (7.10) is

   \[ \widetilde L(a)=2a-\tfrac12a^2(k_{++}+k_{--}-2k_{+-})=2a-\tfrac12a^2d^2, \]

   maximised at \(a=2/d^2\gt0\). By (7.8), \(\mathbf{w}\)\({}=a(\bm\phi_+-\bm\phi_-)\)\({}=2(\bm\phi_+-\bm\phi_-)/d^2\), so \(\|\mathbf{w}\|=2/d\), the margin is \(\rho=d/2\), and both points are support vectors. The bias comes from \(t_+y_+=1\): \(b=1-\mathbf{w}^\T\bm\phi_+\). The midpoint \(\bar{\bm\phi}=(\bm\phi_++\bm\phi_-)/2\) has \(y(\bar{\bm\phi})\)\({}=\mathbf{w}^\T\bar{\bm\phi}+1\)\({}-\mathbf{w}^\T\bm\phi_+\)\({}=\mathbf{w}^\T(\bm\phi_--\bm\phi_+)/2+1\)\({}=-\frac{2}{d^2}\cdot\frac{d^2}2+1\)\({}=0\). Hence the hyperplane is the perpendicular bisector of the segment joining the two points (its normal \(\mathbf{w}\) is parallel to \(\bm\phi_+-\bm\phi_-\)). Geometrically: any separating hyperplane has distances \(r_+\) and \(r_-\) to the two points on opposite sides with \(r_++r_-\le d\) (the segment crosses it), so the margin is at most \(d/2\), with equality only for the bisector. Only \(d\) enters, so the dimension of the space plays no role. (In input space with a nonlinear kernel, this bisector is a curved surface.)

4. Exercise 7.4 (★★, www). Show that the margin \(\rho\) of the maximum margin hyperplane satisfies \(\rho^{-2}=\sum_na_n\), with \(\{a_n\}\) maximising (7.10) under (7.11), (7.12).

   Solution

   See Result 7.5. In short: \(\rho=1/\|\mathbf{w}\|\); by (7.8), (7.9) and \(\bm\phi_n^\T\mathbf{w}=y_n-b\), \(\|\mathbf{w}\|^2=\sum_na_nt_ny_n\); by (7.16) \(a_nt_ny_n=a_n\); so \(\rho^{-2}=\|\mathbf{w}\|^2=\sum_na_n\).

5. Exercise 7.5 (★★). Show that also \(\rho^{-2}=2\widetilde L(\mathbf{a})\) and \(\rho^{-2}=\|\mathbf{w}\|^2\).

   Solution

   See Result 7.5. The second is the definition of the margin in canonical form, (7.4). For the first, (7.8) turns the quadratic term of (7.10) into \(\frac12\|\mathbf{w}\|^2\), so \(\widetilde L=\sum_na_n-\frac12\|\mathbf{w}\|^2\), and with \(\sum_na_n=\|\mathbf{w}\|^2\) from Exercise 7.4, \(\widetilde L\)\({}=\frac12\|\mathbf{w}\|^2\)\({}=\frac1{2\rho^2}\). (Consistency: at the optimum the primal value \(\frac12\|\mathbf{w}\|^2\) equals the dual value, as strong duality says.)

6. Exercise 7.6 (★). For logistic regression with \(t\in\{-1,1\}\), \(p(t=1\given y)=\sigma(y)\), show that the regularised negative log likelihood is (7.47).

   Solution

   See Result 7.8. \(p(t=-1\given y)=1-\sigma(y)=\sigma(-y)\), so \(p(t\given y)=\sigma(yt)\) for both values of \(t\). The negative log of the product over independent observations is \(-\sum_n\ln\sigma(y_nt_n)\)\({}=\sum_n\ln(1+e^{-y_nt_n})\), and adding the regulariser \(\lambda\|\mathbf{w}\|^2\) (the negative log of a Gaussian prior up to constants) gives (7.47), (7.48).

7. Exercise 7.7 (★). From the Lagrangian (7.56) of SVM regression, set the derivatives with respect to \(\mathbf{w},b,\xi_n,\hat\xi_n\) to zero and back-substitute to obtain the dual (7.61).

   Solution

   Done in full in Result 7.9: the stationarity conditions are (7.57)–(7.60); the slack terms cancel because of (7.59), (7.60); the remaining terms are \(\frac12\|\mathbf{w}\|^2\)\({}-\epsilon\sum(a_n+\hat a_n)\)\({}+\sum(a_n-\hat a_n)t_n-\sum(a_n-\hat a_n)y_n\), and \(\sum(a_n-\hat a_n)y_n=\|\mathbf{w}\|^2\) by (7.57), (7.58), which leaves (7.61).

8. Exercise 7.8 (★, www). In SVM regression, show that every point with \(\xi_n\gt0\) has \(a_n=C\), and every point with \(\hat\xi_n\gt0\) has \(\hat a_n=C\).

   Solution

   The KKT conditions for the constraints \(\xi_n\ge0\) are \(\mu_n\ge0\) and \(\mu_n\xi_n=0\). If \(\xi_n\gt0\), then \(\mu_n=0\); stationarity in \(\xi_n\), (7.59), gives \(a_n=C-\mu_n=C\). The same with hats and (7.60) gives \(\hat a_n=C\). Equivalently (7.67), (7.68).

9. Exercise 7.9 (★). Verify the posterior mean and covariance (7.82), (7.83) of the regression RVM.

   Solution

   Result 7.12: the log posterior is \(-\frac12\mathbf{w}^\T(\mathbf{A}+\beta\bm\Phi^\T\bm\Phi)\mathbf{w}\)\({}+\beta\mathbf{w}^\T\bm\Phi^\T\mathbf{t}+\const\); matching to \(-\frac12(\mathbf{w}-\mathbf{m})^\T\bm\Sigma^{-1}(\mathbf{w}-\mathbf{m})\) gives \(\bm\Sigma^{-1}=\mathbf{A}+\beta\bm\Phi^\T\bm\Phi\) and \(\bm\Sigma^{-1}\mathbf{m}\)\({}=\beta\bm\Phi^\T\mathbf{t}\).

10. Exercise 7.10 (★★, www). Derive the marginal likelihood (7.85) by performing the Gaussian integral over \(\mathbf{w}\) in (7.84) with completing the square.

    Solution

    Route 2 of Result 7.13: the integrand is a constant times \(\exp\{-\frac12(\mathbf{w}-\mathbf{m})^\T\bm\Sigma^{-1}(\mathbf{w}-\mathbf{m})\)\({}-E(\mathbf{m})\}\); integrating gives \((2\pi)^{M/2}|\bm\Sigma|^{1/2}\) and the form \((*)\), which equals (7.85) by the determinant lemma and the Woodbury identity (shown there).

11. Exercise 7.11 (★★). Repeat Exercise 7.10 using the general result (2.115).

    Solution

    Route 1 of Result 7.13. With prior \(\mathbf{w}\)\({}\sim\N(\bm\mu=\mathbf{0},\bm\Lambda^{-1}=\mathbf{A}^{-1})\) and conditional \(\mathbf{t}\given\mathbf{w}\)\({}\sim\N(\bm\Phi\mathbf{w},\mathbf{L}^{-1}=\beta^{-1}\mathbf{I})\), (2.115) says \(p(\mathbf{t})\)\({}=\N(\bm\Phi\bm\mu,\mathbf{L}^{-1}+\bm\Phi\bm\Lambda^{-1}\bm\Phi^\T)\)\({}=\N(\mathbf{0},\beta^{-1}\mathbf{I}+\bm\Phi\mathbf{A}^{-1}\bm\Phi^\T)\)\({}=\N(\mathbf{0},\mathbf{C})\).

12. Exercise 7.12 (★★, www). Show that direct maximisation of (7.85) gives the re-estimation equations (7.87), (7.88) with \(\gamma_i\) as in (7.89).

    Solution

    Result 7.14: differentiate \((*)\); the dependence through \(\mathbf{m}\) vanishes because \(\mathbf{m}\) minimises \(E\). The derivative in \(\alpha_i\) is \(\frac1{2\alpha_i}-\frac12\Sigma_{ii}\)\({}-\frac12m_i^2\); in \(\beta\), \(\frac N{2\beta}-\frac12\|\mathbf{t}\)\({}-\bm\Phi\mathbf{m}\|^2\)\({}-\frac1{2\beta}\sum_i\gamma_i\), using \(\beta\Tr(\bm\Sigma\bm\Phi^\T\bm\Phi)\)\({}=\sum_i\gamma_i\).

13. Exercise 7.13 (★★). Add Gamma hyperpriors \(p(\alpha_i)=\Gam(\alpha_i\given a,b)\), \(p(\beta)=\Gam(\beta\given c,d)\) and obtain the re-estimation formulae by maximising the posterior \(p(\mathbf{t},\bm\alpha,\beta\given\mathbf{X})\).

    Solution

    With \(\Gam(\lambda\given a,b)\)\({}\propto\lambda^{a-1}e^{-b\lambda}\) (toolkit), the objective is \(\ln p(\mathbf{t}\given\bm\alpha,\beta)\)\({}+\sum_i\{(a-1)\ln\alpha_i-b\alpha_i\}\)\({}+(c-1)\ln\beta-d\beta+\const\). The extra derivatives are \((a-1)/\alpha_i-b\) and \((c-1)/\beta-d\). Add them to the derivatives of Result 7.14.

    For \(\alpha_i\): \(\frac1{2\alpha_i}-\frac12\Sigma_{ii}\)\({}-\frac12m_i^2+\frac{a-1}{\alpha_i}-b\)\({}=0\). Multiply by \(2\alpha_i\): \(1-\alpha_i\Sigma_{ii}-\alpha_im_i^2+2(a-1)\)\({}-2b\alpha_i\)\({}=0\), i.e. \(\alpha_i(m_i^2+2b)=\gamma_i+2a-2\):

    \[ \alpha_i^{\rm new}=\frac{\gamma_i+2a-2}{m_i^2+2b}. \]

    For \(\beta\): \(\frac N{2\beta}-\frac12\|\mathbf{t}\)\({}-\bm\Phi\mathbf{m}\|^2\)\({}-\frac{\sum_i\gamma_i}{2\beta}+\frac{c-1}\beta\)\({}-d\)\({}=0\). Multiply by \(2\beta\): \(N-\sum_i\gamma_i-\beta\|\mathbf{t}\)\({}-\bm\Phi\mathbf{m}\|^2+2(c-1)-2d\beta\)\({}=0\):

    \[ \beta^{\rm new}=\frac{N-\sum_i\gamma_i+2c-2}{\|\mathbf{t}-\bm\Phi\mathbf{m}\|^2+2d}. \]

    Checks: \(a=1,b=0,c=1,d=0\) (flat hyperpriors) recovers (7.87), (7.88). The Gamma shape \(a\lt1\) with small rate (a typical "non-informative" choice, \(a,b\to0\)) lowers the numerator by about 2 and makes pruning more aggressive (this is a statement about the mode of the density in \(\alpha\), which is not invariant under reparametrisation).

14. Exercise 7.14 (★★). Derive the predictive distribution (7.90) of the regression RVM and its variance (7.91).

    Solution

    Result 7.15: \(t=\mathbf{w}^\T\bm\phi+\varepsilon\), with \(\mathbf{w}\sim\N(\mathbf{m},\bm\Sigma)\) and independent noise \(\varepsilon\sim\N(0,(\beta^\star)^{-1})\); the sum is Gaussian with mean \(\mathbf{m}^\T\bm\phi\) and variance \((\beta^\star)^{-1}+\bm\phi^\T\bm\Sigma\bm\phi\) (the linear-Gaussian rule, (2.115)).

15. Exercise 7.15 (★★, www). Using (7.94), (7.95), show that (7.85) can be written as (7.96) with \(\lambda\) in (7.97), \(s_i\) and \(q_i\) in (7.98), (7.99).

    Solution

    Result 7.17: \(\ln|\mathbf{C}|\)\({}=\ln|\mathbf{C}_{-i}|+\ln(\alpha_i+s_i)\)\({}-\ln\alpha_i\) and \(\mathbf{t}^\T\mathbf{C}^{-1}\mathbf{t}\)\({}=\mathbf{t}^\T\mathbf{C}_{-i}^{-1}\mathbf{t}\)\({}-q_i^2/(\alpha_i+s_i)\); insert in (7.85).

16. Exercise 7.16 (★). Show by the second derivative that the stationary point (7.101) of the log marginal likelihood is a maximum.

    Solution

    Result 7.18. Write \(\lambda'=h(\alpha)/(2(\alpha+s)^2)\) with \(h=s^2/\alpha-(q^2-s)\). At a zero of \(h\), the second term of the product rule vanishes and \(\lambda''\)\({}=h'(\alpha)/(2(\alpha+s)^2)\)\({}=-s^2/(2\alpha^2(\alpha+s)^2)\)\({}\lt0\). In fact \(h\) is strictly decreasing, which also proves uniqueness (and global maximality) of the stationary point when \(q^2\gt s\).

17. Exercise 7.17 (★★). Using (7.83), (7.86) and the Woodbury identity (C.7), show that \(Q_i\) and \(S_i\) of (7.102), (7.103) take the form (7.106), (7.107).

    Solution

    Result 7.19: Woodbury applied to \(\mathbf{C}\)\({}=\beta^{-1}\mathbf{I}\)\({}+\bm\Phi\mathbf{A}^{-1}\bm\Phi^\T\) gives \(\mathbf{C}^{-1}\)\({}=\beta\mathbf{I}\)\({}-\beta^2\bm\Phi(\mathbf{A}+\beta\bm\Phi^\T\bm\Phi)^{-1}\bm\Phi^\T\)\({}=\beta\mathbf{I}\)\({}-\beta^2\bm\Phi\bm\Sigma\bm\Phi^\T\). Then \(Q_i=\bm\varphi_i^\T\mathbf{C}^{-1}\mathbf{t}\) and \(S_i=\bm\varphi_i^\T\mathbf{C}^{-1}\bm\varphi_i\) give (7.106), (7.107).

18. Exercise 7.18 (★, www). Show that the gradient and Hessian of the log posterior (7.109) of the classification RVM are (7.110), (7.111).

    Solution

    Result 7.20: \(\partial\ell_n/\partial a_n=t_n-y_n\) for \(\ell_n=t_n\ln y_n+(1-t_n)\ln(1-y_n)\) with \(y_n=\sigma(a_n)\) and \(\sigma'=\sigma(1-\sigma)\); \(\partial a_n/\partial\mathbf{w}=\bm\phi_n\); the prior contributes \(-\mathbf{A}\mathbf{w}\) and, on differentiating again, \(-\mathbf{A}\). The second derivative of the data term is \(-\sum_ny_n(1-y_n)\bm\phi_n\bm\phi_n^\T\)\({}=-\bm\Phi^\T\mathbf{B}\bm\Phi\).

19. Exercise 7.19 (★★). Show that maximising the approximate marginal likelihood (7.114) of the classification RVM gives the re-estimation (7.116).

    Solution

    Result 7.21. Differentiate \(\ln p(\mathbf{t}\given\mathbf{w}^\star)\)\({}+\ln p(\mathbf{w}^\star\given\bm\alpha)\)\({}+\frac12\ln|\bm\Sigma|\) in \(\alpha_i\): the dependence through \(\mathbf{w}^\star\) vanishes (stationarity), the prior gives \(\frac1{2\alpha_i}-\frac12(w_i^\star)^2\), and \(\frac12\ln|\bm\Sigma|\) gives \(-\frac12\Sigma_{ii}\) when \(\mathbf{B}\) is held fixed, which is (7.115); then \(\alpha_i(w_i^\star)^2\)\({}=1-\alpha_i\Sigma_{ii}\)\({}=\gamma_i\). Note that holding \(\mathbf{B}\) fixed is an approximation (the book's), not an exact statement.

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

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 7, pp. 325–358.
2. Allwein, E. L., Schapire, R. E. and Singer, Y. (2000). Reducing multiclass to binary: a unifying approach for margin classifiers. Journal of Machine Learning Research, 1, 113–141.
3. 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.
4. Burges, C. J. C. (1998). A tutorial on support vector machines for pattern recognition. Data Mining and Knowledge Discovery, 2(2), 121–167. doi:10.1023/A:1009715923555.
5. Cortes, C. and Vapnik, V. (1995). Support-vector networks. Machine Learning, 20(3), 273–297. doi:10.1007/BF00994018.
6. Dietterich, T. G. and Bakiri, G. (1995). Solving multiclass learning problems via error-correcting output codes. Journal of Artificial Intelligence Research, 2, 263–286. doi:10.1613/jair.105.
7. Faul, A. C. and Tipping, M. E. (2002). Analysis of sparse Bayesian learning. In Advances in Neural Information Processing Systems 14, 383–389. MIT Press.
8. Lee, Y., Lin, Y. and Wahba, G. (2004). Multicategory support vector machines: theory and application to the classification of microarray data and satellite radiance data. Journal of the American Statistical Association, 99(465), 67–81. doi:10.1198/016214504000000098. (The book cites the 2001 technical report.)
9. MacKay, D. J. C. (1992). Bayesian interpolation. Neural Computation, 4(3), 415–447. doi:10.1162/neco.1992.4.3.415. (The evidence approximation behind (7.87)–(7.88).)
10. McAllester, D. A. (2003). PAC-Bayesian stochastic model selection. Machine Learning, 51(1), 5–21. doi:10.1023/A:1021840411064.
11. Osuna, E., Freund, R. and Girosi, F. (1996). An improved training algorithm for support vector machines. In Neural Networks for Signal Processing VII, 276–285.
12. Platt, J. C. (1999). Fast training of support vector machines using sequential minimal optimization. In Advances in Kernel Methods: Support Vector Learning, 185–208. MIT Press.
13. Platt, J. C. (2000). Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. In Advances in Large Margin Classifiers, 61–73. MIT Press.
14. Platt, J. C., Cristianini, N. and Shawe-Taylor, J. (2000). Large margin DAGs for multiclass classification. In Advances in Neural Information Processing Systems 12, 547–553. MIT Press.
15. Rasmussen, C. E. and Quiñonero-Candela, J. (2005). Healing the relevance vector machine through augmentation. In Proceedings of the 22nd International Conference on Machine Learning, 689–696.
16. Schölkopf, B., Smola, A. J., Williamson, R. C. and Bartlett, P. L. (2000). New support vector algorithms. Neural Computation, 12(5), 1207–1245. doi:10.1162/089976600300015565.
17. Schölkopf, B., Platt, J. C., Shawe-Taylor, J., Smola, A. J. and Williamson, R. C. (2001). Estimating the support of a high-dimensional distribution. Neural Computation, 13(7), 1443–1471. doi:10.1162/089976601750264965.
18. Smola, A. J. and Schölkopf, B. (2004). A tutorial on support vector regression. Statistics and Computing, 14(3), 199–222. doi:10.1023/B:STCO.0000035301.49549.88.
19. Tax, D. M. J. and Duin, R. P. W. (1999). Support vector domain description. Pattern Recognition Letters, 20(11–13), 1191–1199. doi:10.1016/S0167-8655(99)00087-2.
20. Tipping, M. E. (2001). Sparse Bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 1, 211–244.
21. Tipping, M. E. and Faul, A. C. (2003). Fast marginal likelihood maximisation for sparse Bayesian models. In Proceedings of the Ninth International Workshop on Artificial Intelligence and Statistics.
22. Tong, S. and Koller, D. (2000). Restricted Bayes optimal classifiers. In Proceedings of the 17th National Conference on Artificial Intelligence, 658–664.
23. Valiant, L. G. (1984). A theory of the learnable. Communications of the ACM, 27(11), 1134–1142. doi:10.1145/1968.1972.
24. Vapnik, V. N. (1995). The Nature of Statistical Learning Theory. Springer, New York. (2nd edition 2000.)
25. Vapnik, V. N. (1998). Statistical Learning Theory. Wiley, New York.
26. Vapnik, V. N. and Chervonenkis, A. Ya. (1971). On the uniform convergence of relative frequencies of events to their probabilities. Theory of Probability and its Applications, 16(2), 264–280. doi:10.1137/1116025.
27. Veropoulos, K., Campbell, C. and Cristianini, N. (1999). Controlling the sensitivity of support vector machines. In Proceedings of the International Joint Conference on Artificial Intelligence, 55–60.
28. Weston, J. and Watkins, C. (1999). Support vector machines for multi-class pattern recognition. In Proceedings of the Seventh European Symposium on Artificial Neural Networks, 219–224.

DOIs were checked against Crossref; references without a DOI have none that I could verify. The numerical statements made in the text (the RVM marginal likelihood in its two forms, the stationarity conditions for \(\alpha_i\) and \(\beta\), the sparsity decomposition (7.96)–(7.99), the formulas (7.104)–(7.107) and the Laplace-approximation equivalence) were checked on random examples with numpy.