Linear Models for Classification

Classification assigns an input to one of a finite number of classes. This chapter studies models whose decision boundaries are hyperplanes in a (possibly fixed non-linear) feature space. It first fits such boundaries directly (least squares, Fisher's discriminant, the perceptron), then derives them from generative class models (a Gaussian or exponential-family density per class, which makes the posterior a logistic sigmoid or softmax of a linear function), then fits the posterior directly by logistic regression with an iterative reweighted least squares algorithm. The last two sections introduce the Laplace approximation and use it to make logistic regression Bayesian. Every step the book leaves to the reader is written out below.

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

4.1 Discriminant functions

pp. 181–196 A discriminant is a function that maps an input straight to a class label, with no probabilities in between. The word "linear" says that the surfaces separating the classes are hyperplanes; the chapter starts with the geometry of such surfaces and then meets three ways of choosing them: least squares, Fisher's criterion and the perceptron.

4.1.1 Two classes

pp. 181–182 The simplest discriminant is an affine function of the input,

\[ y(\mathbf{x})=\mathbf{w}^\T\mathbf{x}+w_0, \tag{4.4} \]

with the rule "assign \(\mathbf{x}\) to \(\mathcal{C}_1\) if \(y(\mathbf{x})\ge0\), otherwise to \(\mathcal{C}_2\)". The weight vector \(\mathbf{w}\) is a normal to the decision surface \(y(\mathbf{x})=0\) and \(w_0\) is the bias (not the bias of statistics; \(-w_0\) is a threshold).

4.1.2 Multiple classes

pp. 182–184 With \(K>2\) classes one is tempted to combine two-class discriminants. One-versus-the-rest trains \(K-1\) classifiers, the \(k\)th separating \(\mathcal{C}_k\) from the union of the others; one-versus-one trains \(K(K-1)/2\) classifiers, one per pair, and takes a majority vote. Both leave ambiguous regions (Figure 4.2 of the book, a schematic with three classes in the plane):

• One-versus-the-rest. Take two classifiers, "\(\mathcal{C}_1\) or not" and "\(\mathcal{C}_2\) or not". A point on the \(\mathcal{C}_1\) side of the first and the \(\mathcal{C}_2\) side of the second is claimed by both; a point on the "not" side of both is claimed by neither.
• One-versus-one. With three classes the three pairwise lines can form a triangle in which the votes are cyclic, \(\mathcal{C}_1\) beats \(\mathcal{C}_2\) beats \(\mathcal{C}_3\) beats \(\mathcal{C}_1\), so there is no majority.

The remedy is a single \(K\)-class discriminant of \(K\) linear functions,

\[ y_k(\mathbf{x})=\mathbf{w}_k^\T\mathbf{x}+w_{k0}, \tag{4.9} \]

with the rule: assign \(\mathbf{x}\) to \(\mathcal{C}_k\) if \(y_k(\mathbf{x})>y_j(\mathbf{x})\) for all \(j\ne k\). Every point has one largest value (ties occur only on the boundaries), so there are no ambiguous regions. The boundary between \(\mathcal{C}_k\) and \(\mathcal{C}_j\) is where the two functions are equal, which is again a hyperplane,

\[ (\mathbf{w}_k-\mathbf{w}_j)^\T\mathbf{x}+(w_{k0}-w_{j0})=0, \tag{4.10} \]

with the geometry of Result 4.1 (normal \(\mathbf{w}_k-\mathbf{w}_j\)).

For \(K=2\) the one-function form (4.4) and the two-function form (4.9) are equivalent: take \(y=y_1-y_2\), \(\mathbf{w}=\mathbf{w}_1-\mathbf{w}_2\), \(w_0=w_{10}-w_{20}\).

4.1.3 Least squares for classification

pp. 184–186 Chapter 3 gave a closed-form solution for linear models under a sum-of-squares error, so it is tempting to reuse it. Code the class of \(\mathbf{x}_n\) by the 1-of-\(K\) target \(\mathbf{t}_n\) and give each class its own linear model, collected as

\[ \mathbf{y}(\mathbf{x})=\widetilde{\mathbf{W}}^\T\widetilde{\mathbf{x}}, \tag{4.14} \]

where the \(k\)th column of \(\widetilde{\mathbf{W}}\) is \(\tilde{\mathbf{w}}_k=(w_{k0},\mathbf{w}_k^\T)^\T\) (so (4.14) stacks the \(K\) functions (4.13)), and assign \(\mathbf{x}\) to the class with the largest \(y_k\). One motivation is that least squares estimates the conditional mean \(\mathbb{E}[\mathbf{t}\given\mathbf{x}]\), and for a 1-of-\(K\) code that is the vector of posterior class probabilities.

The sum-to-one property is not enough for probabilities: nothing keeps the outputs inside \((0,1)\), and a linear model must leave that interval somewhere.

Where least squares fails

• Not robust. A point far on the correct side of the boundary is "too correct": its residual is large, so it drags the fit (Figure 4.4, right). Section 7.1.2 meets error functions without this flaw.
• Masking. With three or more classes whole classes can disappear (Figure 4.5, left). A closed-form example: put the same number of points at the three positions \(x=-2,\,0,\,2\) on a line, one class per position. Least squares on each 1-of-3 column gives (computed exactly, one point per class) \(y_1=\tfrac13-\tfrac14x\), \(y_2=\tfrac13\), \(y_3=\tfrac13+\tfrac14x\): the middle class has a flat line that never strictly exceeds the maximum of the other two (all three are equal at \(x=0\)), so the middle class owns no region at all, although it is trivially separable by two thresholds.
• Wrong noise model. Least squares is maximum likelihood for Gaussian conditional targets (Chapter 3), and binary 1-of-\(K\) vectors are anything but Gaussian. The probabilistic models below use a likelihood that matches the data.

4.1.4 Fisher's linear discriminant

pp. 186–189 A linear classifier can be seen as a dimension reduction: project \(\mathbf{x}\in\R^D\) to one number,

\[ y=\mathbf{w}^\T\mathbf{x}, \tag{4.20} \]

and classify by a threshold on \(y\) (the earlier rule with \(y_0=-w_0\)). Let \(N_1,N_2\) be the class sizes, with class means

\[ \mathbf{m}_1=\frac1{N_1}\sum_{n\in\mathcal{C}_1}\mathbf{x}_n,\qquad \mathbf{m}_2=\frac1{N_2}\sum_{n\in\mathcal{C}_2}\mathbf{x}_n. \tag{4.21} \]

The first idea is to separate the projected means: maximise \(m_2-m_1=\mathbf{w}^\T(\mathbf{m}_2-\mathbf{m}_1)\) (4.22), with \(m_k=\mathbf{w}^\T\mathbf{m}_k\) (4.23). This is unbounded in \(\|\mathbf{w}\|\), so constrain \(\mathbf{w}^\T\mathbf{w}=1\).

This is poor when the class covariances are strongly non-diagonal: the classes can be well separated in \(\R^D\) yet overlap on the line joining their means (Figure 4.6, left). Fisher's idea is to ask for a large separation of the projected means relative to the projected within-class spread. With \(y_n=\mathbf{w}^\T\mathbf{x}_n\), define

\[ s_k^2=\sum_{n\in\mathcal{C}_k}(y_n-m_k)^2, \tag{4.24} \]

Rewriting (Exercise 4.5). The numerator is \((m_2-m_1)^2\)\({}=(\mathbf{w}^\T(\mathbf{m}_2-\mathbf{m}_1))^2\)\({}=\mathbf{w}^\T(\mathbf{m}_2-\mathbf{m}_1)(\mathbf{m}_2-\mathbf{m}_1)^\T\mathbf{w}\). In the denominator \(s_k^2\)\({}=\sum_{n\in\mathcal{C}_k}(\mathbf{w}^\T\mathbf{x}_n-\mathbf{w}^\T\mathbf{m}_k)^2\)\({}=\mathbf{w}^\T\big[\sum_{n\in\mathcal{C}_k}(\mathbf{x}_n-\mathbf{m}_k)(\mathbf{x}_n-\mathbf{m}_k)^\T\big]\mathbf{w}\); summing over \(k\) gives \(\mathbf{w}^\T\mathbf{S}_\mathrm{W}\mathbf{w}\).

If \(\mathbf{S}_\mathrm{W}\propto\mathbf{I}\) the direction reduces to \(\mathbf{m}_2-\mathbf{m}_1\). The result is called Fisher's linear discriminant, though strictly it only chooses a direction. To get a classifier, model the projected classes \(p(y\given\mathcal{C}_k)\) (for example as Gaussians, plausible by the central limit theorem since \(y\) is a sum of random variables) and choose the threshold \(y_0\) by decision theory (Chapter 1). If both projected classes are Gaussians with the same variance \(\sigma^2\) and means \(\mu_1,\mu_2\), the Bayes threshold satisfies \(y_0\)\({}=\tfrac12(\mu_1+\mu_2)\)\({}+\sigma^2\ln\{p(\mathcal{C}_2)/p(\mathcal{C}_1)\}/(\mu_1-\mu_2)\).

Fisher's direction is also the discriminant (4.66) that Gaussian class models with a shared covariance produce, with \(\mathbf{S}_\mathrm{W}/N\) estimating \(\bm\Sigma\) (Section 4.2.1).

4.1.5 Relation to least squares

pp. 189–190 For two classes Fisher's direction is a least-squares solution for a particular coding of the targets. Use targets \(t_n=N/N_1\) for \(\mathcal{C}_1\) and \(t_n=-N/N_2\) for \(\mathcal{C}_2\) (roughly the reciprocal of the class prior, with a sign), and minimise

\[ E=\frac12\sum_{n=1}^N(\mathbf{w}^\T\mathbf{x}_n+w_0-t_n)^2. \tag{4.31} \]

4.1.6 Fisher's discriminant for multiple classes

pp. 191–192 For \(K>2\) classes and \(D>K\) take \(D'>1\) linear features \(\mathbf{y}=\mathbf{W}^\T\mathbf{x}\) (4.39), the columns of \(\mathbf{W}\) being the weight vectors (no biases). The within-class scatter generalises (4.28):

\[ \mathbf{S}_\mathrm{W}=\sum_{k=1}^K\mathbf{S}_k,\quad \mathbf{S}_k=\sum_{n\in\mathcal{C}_k}(\mathbf{x}_n-\mathbf{m}_k)(\mathbf{x}_n-\mathbf{m}_k)^\T,\quad \mathbf{m}_k=\frac1{N_k}\sum_{n\in\mathcal{C}_k}\mathbf{x}_n, \tag{4.40–4.42} \]

the total scatter is \(\mathbf{S}_\mathrm{T}\)\({}=\sum_{n=1}^N(\mathbf{x}_n-\mathbf{m})(\mathbf{x}_n-\mathbf{m})^\T\) with \(\mathbf{m}\)\({}=\frac1N\sum_n\mathbf{x}_n\)\({}=\frac1N\sum_kN_k\mathbf{m}_k\) (4.43, 4.44), and it splits as

\[ \mathbf{S}_\mathrm{T}=\mathbf{S}_\mathrm{W}+\mathbf{S}_\mathrm{B},\qquad \mathbf{S}_\mathrm{B}=\sum_{k=1}^KN_k(\mathbf{m}_k-\mathbf{m})(\mathbf{m}_k-\mathbf{m})^\T. \tag{4.45, 4.46} \]

Proof of the split (the book states it): as in the two-class proof, \(\sum_{n\in\mathcal{C}_k}(\mathbf{x}_n-\mathbf{m})(\mathbf{x}_n-\mathbf{m})^\T\)\({}=\mathbf{S}_k\)\({}+N_k(\mathbf{m}_k-\mathbf{m})(\mathbf{m}_k-\mathbf{m})^\T\) (the cross terms contain \(\sum_{n\in\mathcal{C}_k}(\mathbf{x}_n-\mathbf{m}_k)\)\({}=\mathbf{0}\)); sum over \(k\). For \(K=2\), \(\mathbf{S}_\mathrm{B}\) of (4.46) equals \(\frac{N_1N_2}N\) times the two-class \(\mathbf{S}_\mathrm{B}\) of (4.27) (see the proof of Result 4.5), which changes nothing since only directions matter.

In the projected space the same matrices are \(\mathbf{s}_\mathrm{W}\)\({}=\mathbf{W}^\T\mathbf{S}_\mathrm{W}\mathbf{W}\) and \(\mathbf{s}_\mathrm{B}\)\({}=\mathbf{W}^\T\mathbf{S}_\mathrm{B}\mathbf{W}\) (4.47, 4.48, with the class means (4.49) of the projected data). A scalar that is large when the between-class scatter is large relative to the within-class scatter is

\[ J(\mathbf{W})=\operatorname{Tr}\{\mathbf{s}_\mathrm{W}^{-1}\mathbf{s}_\mathrm{B}\}=\operatorname{Tr}\{(\mathbf{W}^\T\mathbf{S}_\mathrm{W}\mathbf{W})^{-1}(\mathbf{W}^\T\mathbf{S}_\mathrm{B}\mathbf{W})\}. \tag{4.50, 4.51} \]

Erratum. The book's (4.47), (4.48), (4.50) print \(\mathbf{s}_\mathrm{W},\mathbf{s}_\mathrm{B}\) with \(s\) for the projected matrices where capital \(\mathbf{S}\) is meant, and (4.51) prints \((\mathbf{W}\mathbf{S}_\mathrm{W}\mathbf{W}^\T)^{-1}(\mathbf{W}\mathbf{S}_\mathrm{B}\mathbf{W}^\T)\); with \(\mathbf{y}=\mathbf{W}^\T\mathbf{x}\) the correct form is the one above (errata, p. 192).

For \(K=2\) there is a single non-zero eigenvalue, with eigenvector \(\mathbf{S}_\mathrm{W}^{-1}(\mathbf{m}_2-\mathbf{m}_1)\), recovering Result 4.4. The projected features are then fed to any classifier (Fukunaga, 1990).

4.1.7 The perceptron algorithm

pp. 192–196 Rosenblatt's perceptron is a two-class model: a fixed non-linear feature vector \(\bm\phi(\mathbf{x})\) (including \(\phi_0=1\)) is passed through a step function,

\[ y(\mathbf{x})=f(\mathbf{w}^\T\bm\phi(\mathbf{x})),\qquad f(a)=\begin{cases}+1,&a\ge0\\-1,&a<0.\end{cases} \tag{4.52, 4.53} \]

Targets are \(t=+1\) for \(\mathcal{C}_1\) and \(t=-1\) for \(\mathcal{C}_2\). We want \(\mathbf{w}^\T\bm\phi_nt_n>0\) for every \(n\). The number of misclassified points is a bad error function: it is piecewise constant in \(\mathbf{w}\), with zero gradient almost everywhere. The perceptron criterion is instead

\[ E_\mathrm{P}(\mathbf{w})=-\sum_{n\in\mathcal{M}}\mathbf{w}^\T\bm\phi_nt_n, \tag{4.54} \]

where \(\mathcal{M}\) is the set of currently misclassified points and \(\bm\phi_n=\bm\phi(\mathbf{x}_n)\) (erratum, p. 194).

One update reduces the error contribution of the point it was made on: using \(\eta=1\),

\[ -\mathbf{w}^{(\tau+1)\T}\bm\phi_nt_n=-\mathbf{w}^{(\tau)\T}\bm\phi_nt_n-(\bm\phi_nt_n)^\T\bm\phi_nt_n<-\mathbf{w}^{(\tau)\T}\bm\phi_nt_n, \tag{4.56} \]

because \(\|\bm\phi_nt_n\|^2=\|\bm\phi_n\|^2>0\). But that does not make \(E_\mathrm{P}\) decrease: other points may become misclassified, so the criterion can go up at some steps. The following theorem is why the algorithm nevertheless works. The book only cites it.

The remaining limitations the book lists are that the perceptron has no probabilistic output, does not generalise readily to \(K>2\) classes, and, like every model of this chapter, uses fixed basis functions. Figure 4.8 of the book, a photograph of Rosenblatt's Mark 1 hardware, is not reproduced here. Widrow and Hoff's adaline (1960) had the same form but was trained by least squares on the activation.

4.2 Probabilistic generative models

pp. 196–203 Now model how the data arise: a class-conditional density \(p(\mathbf{x}\given\mathcal{C}_k)\) and a prior \(p(\mathcal{C}_k)\), and obtain the posterior by Bayes' theorem. The question is which assumptions give a linear classifier.

The name "softmax" is explained by the case \(a_k\gg a_j\) for all \(j\ne k\), where \(p(\mathcal{C}_k\given\mathbf{x})\approx1\): it is a smoothed arg max. Rewriting Bayes' theorem as a sigmoid is only useful if \(a(\mathbf{x})\) has a simple form; the following sections show when it is linear.

4.2.1 Continuous inputs

pp. 198–200 Let every class density be Gaussian with the same covariance \(\bm\Sigma\):

\[ p(\mathbf{x}\given\mathcal{C}_k)=\frac1{(2\pi)^{D/2}|\bm\Sigma|^{1/2}}\exp\Big\{-\tfrac12(\mathbf{x}-\bm\mu_k)^\T\bm\Sigma^{-1}(\mathbf{x}-\bm\mu_k)\Big\}. \tag{4.64} \]

Erratum. The book's (4.63) lacks a pair of parentheses; the correct form is \(a_k\)\({}=\ln\big(p(\mathbf{x}\given\mathcal{C}_k)p(\mathcal{C}_k)\big)\), used above.

Geometry. Since \(\bm\mu_1^\T\bm\Sigma^{-1}\bm\mu_1\)\({}-\bm\mu_2^\T\bm\Sigma^{-1}\bm\mu_2\)\({}=\mathbf{w}^\T(\bm\mu_1+\bm\mu_2)\), (4.67) can be written \(w_0\)\({}=-\tfrac12\mathbf{w}^\T(\bm\mu_1+\bm\mu_2)\)\({}+\ln\frac{p(\mathcal{C}_1)}{p(\mathcal{C}_2)}\), so

\[ a(\mathbf{x})=\mathbf{w}^\T\big(\mathbf{x}-\tfrac12(\bm\mu_1+\bm\mu_2)\big)+\ln\frac{p(\mathcal{C}_1)}{p(\mathcal{C}_2)}. \]

The boundary \(a=0\) is perpendicular to \(\bm\Sigma^{-1}(\bm\mu_1-\bm\mu_2)\) and passes through \(\tfrac12(\bm\mu_1+\bm\mu_2)\)\({}-\ln\frac{p(\mathcal{C}_1)}{p(\mathcal{C}_2)}\,\frac{\bm\mu_1-\bm\mu_2}{\Delta^2}\), where \(\Delta^2\)\({}=(\bm\mu_1-\bm\mu_2)^\T\bm\Sigma^{-1}(\bm\mu_1-\bm\mu_2)\) is the squared Mahalanobis distance between the means. The priors enter only through the bias, so changing them shifts the boundary and the contours of constant posterior in parallel, without rotating them.

Quadratic discriminant. If each class has its own \(\bm\Sigma_k\), the terms \(-\tfrac12\mathbf{x}^\T\bm\Sigma_k^{-1}\mathbf{x}\) no longer cancel and \(a_k(\mathbf{x})\)\({}=-\tfrac12\mathbf{x}^\T\bm\Sigma_k^{-1}\mathbf{x}\)\({}+\mathbf{x}^\T\bm\Sigma_k^{-1}\bm\mu_k\)\({}-\tfrac12\bm\mu_k^\T\bm\Sigma_k^{-1}\bm\mu_k\)\({}-\tfrac12\ln|\bm\Sigma_k|+\ln p(\mathcal{C}_k)\) is quadratic, giving quadratic decision boundaries (Figure 4.11; the regions are then no longer convex in general).

4.2.2 Maximum likelihood solution

pp. 200–202 Given labelled data \(\{\mathbf{x}_n,t_n\}\), fit the parameters by maximum likelihood. Two classes: \(t_n=1\) for \(\mathcal{C}_1\), \(t_n=0\) for \(\mathcal{C}_2\), \(p(\mathcal{C}_1)=\pi\). A point of class \(\mathcal{C}_1\) contributes \(p(\mathbf{x}_n,\mathcal{C}_1)\)\({}=\pi\N(\mathbf{x}_n\given\bm\mu_1,\bm\Sigma)\), one of \(\mathcal{C}_2\) contributes \((1-\pi)\N(\mathbf{x}_n\given\bm\mu_2,\bm\Sigma)\), so (erratum: the left side is the joint density of labels and inputs)

\[ p(\mathbf{t},\mathbf{X}\given\pi,\bm\mu_1,\bm\mu_2,\bm\Sigma)=\prod_{n=1}^N\big[\pi\N(\mathbf{x}_n\given\bm\mu_1,\bm\Sigma)\big]^{t_n}\big[(1-\pi)\N(\mathbf{x}_n\given\bm\mu_2,\bm\Sigma)\big]^{1-t_n}. \tag{4.71} \]

Fitting Gaussians is not robust to outliers, because maximum likelihood for a Gaussian is not (Section 2.3.7).

4.2.3 Discrete features

p. 202 For binary features \(x_i\in\{0,1\}\), \(i=1,\dots,D\), a general class-conditional distribution needs \(2^D-1\) numbers per class. The naive Bayes assumption (features independent given the class, Section 8.2.2) uses \(D\):

\[ p(\mathbf{x}\given\mathcal{C}_k)=\prod_{i=1}^D\mu_{ki}^{x_i}(1-\mu_{ki})^{1-x_i}. \tag{4.81} \]

Substituting into (4.63),

\[ a_k(\mathbf{x})=\sum_{i=1}^D\{x_i\ln\mu_{ki}+(1-x_i)\ln(1-\mu_{ki})\}+\ln p(\mathcal{C}_k), \tag{4.82} \]

which is linear in \(\mathbf{x}\): slope \(\ln\frac{\mu_{ki}}{1-\mu_{ki}}\) on \(x_i\), intercept \(\sum_i\ln(1-\mu_{ki})+\ln p(\mathcal{C}_k)\). For two classes, \(a\)\({}=\sum_ix_i\ln\frac{\mu_{1i}(1-\mu_{2i})}{\mu_{2i}(1-\mu_{1i})}\)\({}+\sum_i\ln\frac{1-\mu_{1i}}{1-\mu_{2i}}\)\({}+\ln\frac{p(\mathcal{C}_1)}{p(\mathcal{C}_2)}\) in (4.57). The ML estimate is \(\mu_{ki}\)\({}=\frac1{N_k}\sum_{n\in\mathcal{C}_k}x_{ni}\) (the fraction of class-\(k\) points with feature \(i\) on; a Bernoulli ML estimate per feature and class).

4.2.4 Exponential family

pp. 202–203 Both cases are instances of one result. Take class-conditional densities in the exponential family (2.194) with \(\mathbf{u}(\mathbf{x})=\mathbf{x}\) and a scaling parameter \(s\) common to all classes (2.236):

\[ p(\mathbf{x}\given\bm\lambda_k,s)=\frac1sh\Big(\frac1s\mathbf{x}\Big)g(\bm\lambda_k)\exp\Big\{\frac1s\bm\lambda_k^\T\mathbf{x}\Big\}. \tag{4.83, 4.84} \]

4.3 Probabilistic discriminative models

pp. 203–213 Section 4.2 fitted class densities and then inverted them. The alternative is to write the generalised linear model for \(p(\mathcal{C}_k\given\bm\phi)\) down directly and maximise the conditional likelihood. This is discriminative training. It models less (no density of \(\mathbf{x}\)), typically needs fewer parameters, and can predict better when the class-density assumptions are wrong. The price is that the data cannot be generated from the model. For \(\bm\phi\in\R^M\), logistic regression has \(M\) parameters; fitting Gaussians with a shared covariance instead needs \(2M\) (means) \(+\,M(M+1)/2\) (covariance) \(+\,1\) (prior) \(=M(M+5)/2+1\), quadratic in \(M\).

4.3.1 Fixed basis functions

pp. 204–205 Everything so far works equally on \(\bm\phi(\mathbf{x})\) in place of \(\mathbf{x}\): decision boundaries are linear in feature space, non-linear in input space, and classes inseparable in \(\mathbf{x}\) can become separable in \(\bm\phi\) (Figure 4.12). One basis function is conventionally \(\phi_0=1\), carrying the bias. A fixed \(\bm\phi\) cannot remove class overlap, because \(p(\mathcal{C}_k\given\mathbf{x})\) already lies strictly between 0 and 1 there and the posterior is a function of \(\mathbf{x}\) only; it can even create overlap (a many-to-one map identifies inputs of different classes) or make the posterior easier to model. The limitation that the basis is fixed motivates the adaptive models of later chapters (Section 3.6).

4.3.2 Logistic regression

pp. 205–206 For two classes the model is

\[ p(\mathcal{C}_1\given\bm\phi)=y(\bm\phi)=\sigma(\mathbf{w}^\T\bm\phi),\qquad p(\mathcal{C}_2\given\bm\phi)=1-y(\bm\phi). \tag{4.87} \]

(Despite the name, this is a classification model.) For targets \(t_n\in\{0,1\}\), \(y_n=\sigma(a_n)\), \(a_n=\mathbf{w}^\T\bm\phi_n\), the likelihood of independent Bernoulli observations is

\[ p(\mathbf{t}\given\mathbf{w})=\prod_{n=1}^Ny_n^{t_n}\{1-y_n\}^{1-t_n}, \tag{4.89} \]

and the negative log likelihood is the cross-entropy error

\[ E(\mathbf{w})=-\ln p(\mathbf{t}\given\mathbf{w})=-\sum_{n=1}^N\{t_n\ln y_n+(1-t_n)\ln(1-y_n)\}. \tag{4.90} \]

The factor \(\sigma'\) cancels and the contribution of point \(n\) is its "error" \(y_n-t_n\) times \(\bm\phi_n\), the same form as the gradient of the sum-of-squares error for linear regression (3.13). A sequential algorithm updates with one term at a time. A numerically safe form of the error is \(E=\sum_n\{\ln(1+e^{a_n})-t_na_n\}\), obtained by substituting \(y_n=\sigma(a_n)\) in (4.90) (use \(\operatorname{logaddexp}\)).

4.3.3 Iterative reweighted least squares

pp. 207–208 The logistic-regression error has no closed-form minimiser, but it is convex and close to quadratic, so Newton-Raphson works well. For minimising \(E(\mathbf{w})\) it updates

\[ \mathbf{w}^{(\mathrm{new})}=\mathbf{w}^{(\mathrm{old})}-\mathbf{H}^{-1}\nabla E(\mathbf{w}), \tag{4.92} \]

with \(\mathbf{H}=\nabla\nabla E\) the Hessian.

Convergence. Near the minimum the convergence is quadratic (standard Newton theory: the error is squared at each step). Far from it, Newton's step can overshoot because \(E\) is not quadratic, so robust implementations add a line search or damping. Because the logistic model uses the canonical link (Section 4.3.6), the Hessian does not depend on the targets, so IRLS is the same as Fisher scoring. On separable data (Result 4.13) \(R_{nn}\to0\) and \(\bm\Phi^\T\mathbf{R}\bm\Phi\) becomes nearly singular: IRLS does not converge but drives \(\|\mathbf{w}\|\) to infinity, with a surprisingly regular pattern. In the one-point example \(\phi=1\), \(t=1\), \(E(w)=\ln(1+e^{-w})\), the Newton step is \(-E'/E''=1/\{\sigma(w)\}\), which tends to 1 as \(w\) grows: the weight increases by about one per iteration (computed: \(w\)\({}=2,\,3.14,\,4.18,\,5.19,\,6.20,\,7.20,\,8.20,\,9.20\) for iterations 1 to 8 from \(w=0\)) while the error falls like \(e^{-w}\). In the demo, on separable data \(\|\mathbf{w}\|\) is 1.9, 5.3 and 10.7 after 5, 15 and 30 IRLS steps while \(\ln L\to0\).

4.3.4 Multiclass logistic regression

pp. 209–210 For \(K\) classes use the softmax of linear activations,

\[ p(\mathcal{C}_k\given\bm\phi)=y_k(\bm\phi)=\frac{\exp(a_k)}{\sum_j\exp(a_j)},\qquad a_k=\mathbf{w}_k^\T\bm\phi. \tag{4.104, 4.105} \]

With 1-of-\(K\) targets \(\mathbf{t}_n\) (matrix \(\mathbf{T}\), entries \(t_{nk}\)) and \(y_{nk}=y_k(\bm\phi_n)\), the likelihood and cross-entropy error are

\[ p(\mathbf{T}\given\mathbf{w}_1,\dots,\mathbf{w}_K)=\prod_{n=1}^N\prod_{k=1}^Ky_{nk}^{t_{nk}},\qquad E(\mathbf{w}_1,\dots,\mathbf{w}_K)=-\sum_{n=1}^N\sum_{k=1}^Kt_{nk}\ln y_{nk}. \tag{4.107, 4.108} \]

The gradient has the same "error times feature" form as linear regression and two-class logistic regression, an instance of the canonical-link result (Result 4.17). The Newton-Raphson (IRLS) update uses the block Hessian (4.110) and gradient (4.109); a practical implementation fixes \(\mathbf{w}_K=\mathbf{0}\) (or adds a ridge) so that the \(M(K-1)\times M(K-1)\) system is non-singular (Bishop and Nabney, 2008). A Figure 4.5-style comparison shows the effect: logistic regression separates the three classes on which least squares masks one.

4.3.5 Probit regression

pp. 210–212 Not every class-conditional model gives a logistic posterior (Gaussian mixtures do not), so other link functions are worth considering within the generalised linear model \(p(t=1\given a)=f(a)\), \(a=\mathbf{w}^\T\bm\phi\) (4.111). One motivation is a noisy threshold model: set \(t_n=1\) if \(a_n\ge\theta\) and \(t_n=0\) otherwise (4.112), with \(\theta\) drawn at random from a density \(p(\theta)\). Then

\[ p(t=1\given a)=P(\theta\le a)=f(a)=\int_{-\infty}^ap(\theta)\,\d\theta, \tag{4.113} \]

the cumulative distribution function of \(\theta\) (Figure 4.13: the density is the slope of the cdf, the cdf is the area under the density). For a standard normal \(\theta\),

\[ \Phi(a)=\int_{-\infty}^a\N(\theta\given0,1)\,\d\theta, \tag{4.114} \]

the inverse probit function (the probit is its inverse; errata). A \(\theta\) with another mean and variance is absorbed in the bias and in the scale of \(\mathbf{w}\), so nothing is lost. Many libraries provide instead the error function, with the corrected definition

\[ \operatorname{erf}(a)=\frac2{\sqrt\pi}\int_0^a\exp(-\theta^2)\,\d\theta,\qquad \Phi(a)=\frac12\Big\{1+\operatorname{erf}\Big(\frac a{\sqrt2}\Big)\Big\}. \tag{4.115, 4.116} \]

Erratum. The book's (4.115) has a factor \(\tfrac12\) in the exponent (\(\exp(-\theta^2/2)\)); the correct integrand is \(\exp(-\theta^2)\), as above. (4.116) is the corrected relation, with \(\Phi\) not bold.

Logit versus probit. The two fits are usually similar. Since \(\sigma(a)\approx\Phi(\lambda a)\) with \(\lambda^2=\pi/8\) (Figure 4.9), logistic weights are expected to be about \(1/\lambda\approx1.6\) times the probit weights (in our demo fits the ratio \(|\mathbf{w}_\mathrm{probit}|/|\mathbf{w}_\mathrm{logistic}|\) is about 0.59, against \(1/\lambda=1.596\) for the exact reciprocal 0.63). They differ in the tails: \(1-\sigma(a)\sim e^{-a}\) but \(1-\Phi(a)\sim e^{-a^2/2}/a\), so a point far on the wrong side of the boundary (an outlier or a mislabel) costs only linearly in the logistic error but quadratically in the probit error, which makes probit regression more sensitive to such points. Both assume correct labels. Mislabelling with flip probability \(\epsilon\) is modelled by

\[ p(t\given\mathbf{x})=(1-\epsilon)\sigma(\mathbf{x})+\epsilon\{1-\sigma(\mathbf{x})\}=\epsilon+(1-2\epsilon)\sigma(\mathbf{x}), \tag{4.117} \]

with \(\epsilon\) fixed or inferred (Opper and Winther, 2000). The label is \(t=1\) either if it is truly 1 and not flipped, or truly 0 and flipped; the second form follows by expanding. The posterior is then squeezed into \([\epsilon,1-\epsilon]\), so a point far from the boundary has a bounded cost \(-\ln\epsilon\).

4.3.6 Canonical link functions

pp. 212–213 The simple gradient "(prediction minus target) times feature" appeared for linear regression (3.12), the sigmoid with cross-entropy (4.90) and the softmax with its cross-entropy (4.108). It is a general property of the exponential family with the right link. Let the target have the conditional density (cf. (4.84), here applied to \(t\), not to the input)

\[ p(t\given\eta,s)=\frac1sh\Big(\frac ts\Big)g(\eta)\exp\Big\{\frac{\eta t}s\Big\}. \tag{4.118} \]

So \(y\) and \(\eta\) are linked, \(\eta=\psi(y)\). A generalised linear model (Nelder and Wedderburn, 1972) has \(y=f(\mathbf{w}^\T\bm\phi)\) (4.120): \(f\) is the activation function, \(f^{-1}\) the link function.

Examples. (i) Gaussian target with precision \(\beta\): \(s=\beta^{-1}\), \(\eta=\mu\), \(g(\eta)=e^{-\beta\eta^2/2}\), \(y=-s\,\d\ln g/\d\eta=\eta\), so the canonical link is the identity (linear regression). (ii) Bernoulli: \(s=1\), \(\eta=\ln\frac y{1-y}\), \(g=1-y=\sigma(-\eta)\), \(y=\sigma(\eta)\): the logit link, hence logistic regression. (iii) Categorical: the softmax. (iv) Poisson: \(p(t)=\lambda^te^{-\lambda}/t!\) has \(\eta=\ln\lambda\), \(g(\eta)=\exp(-e^\eta)\), \(y=\lambda\), so the canonical link is the log, \(\lambda=\exp(\mathbf{w}^\T\bm\phi)\), and again \(\nabla E=\sum_n(y_n-t_n)\bm\phi_n\). The probit link is not canonical for the Bernoulli distribution, as the gradient above shows.

For a canonical link the Hessian is \(\frac1s\sum_nf'(a_n)\bm\phi_n\bm\phi_n^\T\) (differentiate (4.124) once, with \(\d y_n/\d a_n=f'(a_n)\)), which does not contain the targets: Newton-Raphson then equals Fisher scoring, and for the Bernoulli case it is (4.97).

4.4 The Laplace approximation

pp. 213–217 The posterior of logistic regression is not Gaussian, so the Bayesian treatment of Section 3.3 cannot be copied. The Laplace approximation fits a Gaussian to a density by matching its mode and its curvature there. Let \(p(z)=f(z)/Z\) be a density of one continuous variable with unknown normaliser \(Z=\int f(z)\,\d z\) (4.125).

Example (Figure 4.14). For \(f(z)=e^{-z^2/2}\sigma(20z+4)\), \(\ln f=-z^2/2+\ln\sigma(20z+4)\); the mode solves \(z=20\{1-\sigma(20z+4)\}\), and \(A=1+400\,\sigma(20z_0+4)\{1-\sigma(20z_0+4)\}\). Numerically (my computation) \(z_0=0.0775\) and \(A=2.54\); the exact normaliser is \(Z=\int f\,\d z=1.451\) while the Laplace estimate (4.135) is \(f(z_0)\sqrt{2\pi/A}=1.561\), about 8% high. The true density is skewed (heavier on the side where the sigmoid is flat), so a symmetric Gaussian cannot match it; the right panel (negative logs) shows the quadratic against the true curve.

Several variables

For \(\mathbf{z}\in\R^M\) and \(p(\mathbf{z})=f(\mathbf{z})/Z\), at a stationary point \(\mathbf{z}_0\) (\(\nabla f=\mathbf{0}\)) the second-order expansion is

\[ \ln f(\mathbf{z})\simeq\ln f(\mathbf{z}_0)-\tfrac12(\mathbf{z}-\mathbf{z}_0)^\T\mathbf{A}(\mathbf{z}-\mathbf{z}_0),\qquad \mathbf{A}=-\nabla\nabla\ln f(\mathbf{z})\big|_{\mathbf{z}=\mathbf{z}_0}, \tag{4.131, 4.132} \]

the linear term vanishing because \(\nabla\ln f=\nabla f/f=\mathbf{0}\), and the Hessian being symmetric. Exponentiating gives (4.133). Comparing with the Gaussian density (2.43), the normalised approximation is

\[ q(\mathbf{z})=\frac{|\mathbf{A}|^{1/2}}{(2\pi)^{M/2}}\exp\Big\{-\tfrac12(\mathbf{z}-\mathbf{z}_0)^\T\mathbf{A}(\mathbf{z}-\mathbf{z}_0)\Big\}=\N(\mathbf{z}\given\mathbf{z}_0,\mathbf{A}^{-1}). \tag{4.134} \]

It is a proper density iff \(\mathbf{A}\succ0\), i.e. the stationary point is a local maximum, not a minimum or saddle.

Properties and limits

• Easy to apply: it needs a mode (numerical optimisation) and one Hessian, not \(Z\). Multimodal densities give a different approximation per mode.
• Large \(N\): as the data accumulate the posterior concentrates and becomes nearly Gaussian (central limit behaviour), so the approximation improves.
• Real variables only, and not invariant to re-parametrisation. Example: a Gamma-shaped \(p(\tau)\propto\tau^{a-1}e^{-b\tau}\) has mode \((a-1)/b\) in \(\tau\), but in \(u=\ln\tau\) the density is \(\propto e^{au-be^u}\) with mode \(\ln(a/b)\) and curvature \(a\); a Laplace approximation in \(u\) is a log-normal in \(\tau\) with median \(a/b\), different from the Gaussian in \(\tau\). One chooses a parametrisation in which the posterior is closer to Gaussian (for a positive scale, its log).
• Local: it uses only the neighbourhood of one point and can miss global properties (the skew in Figure 4.14, mass in other modes). Chapter 10 gives global alternatives.

4.4.1 Model comparison and BIC

pp. 216–217 The same expansion approximates the normaliser:

\[ Z=\int f(\mathbf{z})\,\d\mathbf{z}\simeq f(\mathbf{z}_0)\int\exp\Big\{-\tfrac12(\mathbf{z}-\mathbf{z}_0)^\T\mathbf{A}(\mathbf{z}-\mathbf{z}_0)\Big\}\d\mathbf{z}=f(\mathbf{z}_0)\frac{(2\pi)^{M/2}}{|\mathbf{A}|^{1/2}}, \tag{4.135} \]

because the integrand is an unnormalised Gaussian whose integral is the reciprocal of the constant in (4.134). Apply this to the model evidence (Section 3.4) of a model with parameters \(\bm\theta\in\R^M\),

\[ p(\mathcal{D})=\int p(\mathcal{D}\given\bm\theta)p(\bm\theta)\,\d\bm\theta, \tag{4.136} \]

identifying \(f(\bm\theta)\)\({}=p(\mathcal{D}\given\bm\theta)p(\bm\theta)\) (so that \(f/Z\) is the posterior) and \(Z=p(\mathcal{D})\).

(4.139) is the Bayesian information criterion (Schwarz, 1978). Compared with the Akaike criterion (1.73), \(\ln p(\mathcal{D}\given\mathbf{w}_\mathrm{ML})-M\), its penalty \(\frac12M\ln N\) is heavier once \(\ln N>2\), i.e. \(N>e^2\approx7.4\). Both are cheap but crude: the full-rank assumption fails when some parameters are poorly determined by the data (the "well-determined parameters" of Section 3.5.3), and then (4.137) with the actual \(\mathbf{A}\) is more accurate, as Section 5.7 illustrates for neural networks. Note that (4.139) replaces the posterior mode by the log likelihood at \(\bm\theta_\mathrm{MAP}\); at large \(N\) the maximum-likelihood value is equivalent.

4.5 Bayesian logistic regression

pp. 217–220 Exact Bayesian inference for logistic regression is intractable: the posterior needs the normaliser of a prior times a product of \(N\) sigmoids, and the predictive distribution integrates a sigmoid against that posterior. The Laplace approximation of Section 4.4 handles the first, and the probit trick the second (Spiegelhalter and Lauritzen, 1990; MacKay, 1992).

4.5.1 Laplace approximation

pp. 217–218 Take a Gaussian prior \(p(\mathbf{w})\)\({}=\N(\mathbf{w}\given\mathbf{m}_0,\mathbf{S}_0)\) (4.140). The posterior is \(p(\mathbf{w}\given\mathbf{t})\)\({}\propto p(\mathbf{w})p(\mathbf{t}\given\mathbf{w})\) (4.141), so with the likelihood (4.89),

\[ \ln p(\mathbf{w}\given\mathbf{t})=-\tfrac12(\mathbf{w}-\mathbf{m}_0)^\T\mathbf{S}_0^{-1}(\mathbf{w}-\mathbf{m}_0)+\sum_{n=1}^N\{t_n\ln y_n+(1-t_n)\ln(1-y_n)\}+\mathrm{const}, \tag{4.142} \]

with \(y_n=\sigma(\mathbf{w}^\T\bm\phi_n)\).

4.5.2 Predictive distribution

pp. 218–220 The predictive probability of \(\mathcal{C}_1\) for a new feature vector marginalises over \(q(\mathbf{w})\):

\[ p(\mathcal{C}_1\given\bm\phi,\mathbf{t})=\int p(\mathcal{C}_1\given\bm\phi,\mathbf{w})p(\mathbf{w}\given\mathbf{t})\,\d\mathbf{w}\simeq\int\sigma(\mathbf{w}^\T\bm\phi)\,q(\mathbf{w})\,\d\mathbf{w}. \tag{4.145} \]

The integrand depends on \(\mathbf{w}\) only through \(a=\mathbf{w}^\T\bm\phi\), so \(\sigma(\mathbf{w}^\T\bm\phi)\)\({}=\int\delta(a-\mathbf{w}^\T\bm\phi)\sigma(a)\,\d a\) (4.146) and

\[ \int\sigma(\mathbf{w}^\T\bm\phi)q(\mathbf{w})\,\d\mathbf{w}=\int\sigma(a)p(a)\,\d a,\qquad p(a)=\int\delta(a-\mathbf{w}^\T\bm\phi)q(\mathbf{w})\,\d\mathbf{w}. \tag{4.147, 4.148} \]

The integral (4.151), a sigmoid convolved with a Gaussian, has no closed form. The approximation uses the near-equality of \(\sigma\) and the inverse probit \(\Phi\) once the horizontal axis is rescaled: \(\sigma(a)\approx\Phi(\lambda a)\).

Consequences. The decision boundary \(p(\mathcal{C}_1\given\bm\phi,\mathbf{t})\)\({}=\tfrac12\) is where \(\kappa\mu_a=0\), i.e. \(\mu_a=0\): the same as the MAP boundary \(\mathbf{w}_\mathrm{MAP}^\T\bm\phi=0\). With equal priors and the minimum-misclassification-rate criterion, marginalisation changes nothing. For other criteria (a loss matrix, a reject option) it matters, because \(\kappa<1\) pulls the probability towards \(\tfrac12\) by an amount that grows with the posterior uncertainty \(\sigma_a^2=\bm\phi^\T\mathbf{S}_N\bm\phi\), which is large far from the training data. In the \(\bm\phi\)-plane the contours of constant \(p(\mathcal{C}_1\given\bm\phi,\mathbf{t})\) are therefore no longer parallel to the MAP boundary: they splay out away from the data (compare Figure 10.13). The approximation also applies only to two classes.

Exercises

All 26 exercises of Chapter 4 (Bishop, 2006, pp. 220–224), 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 summarise the argument and link to the full version.

1. Exercise 4.1 (★★). The convex hull of points \(\{\mathbf{x}_n\}\) is the set of all \(\sum_n\alpha_n\mathbf{x}_n\) with \(\alpha_n\ge0\), \(\sum_n\alpha_n=1\) (4.156). Two point sets \(\{\mathbf{x}_n\}\), \(\{\mathbf{y}_m\}\) are linearly separable if there are \(\hat{\mathbf{w}},w_0\) with \(\hat{\mathbf{w}}^\T\mathbf{x}_n+w_0>0\) for all \(n\) and \(\hat{\mathbf{w}}^\T\mathbf{y}_m+w_0<0\) for all \(m\). Show that if the convex hulls intersect the sets are not separable, and conversely that if they are separable the hulls do not intersect.

   Solution

   Separable \(\Rightarrow\) hulls disjoint. Let \(\mathbf{x}=\sum_n\alpha_n\mathbf{x}_n\) be in the first hull. Then \(\hat{\mathbf{w}}^\T\mathbf{x}+w_0\)\({}=\sum_n\alpha_n(\hat{\mathbf{w}}^\T\mathbf{x}_n+w_0)>0\) (using \(\sum\alpha_n=1\) and at least one \(\alpha_n>0\)). Likewise every point of the second hull has \(\hat{\mathbf{w}}^\T\mathbf{y}+w_0<0\). No point can satisfy both, so the hulls do not meet. The contrapositive is the first claim: intersecting hulls cannot be separable.

   Hulls disjoint \(\Rightarrow\) separable. The two hulls \(H_x,H_y\) are compact (images of the simplex under linear maps) and convex. Choose \(\mathbf{p}\in H_x\), \(\mathbf{q}\in H_y\) at minimum distance (they exist by compactness, and \(\mathbf{p}\ne\mathbf{q}\) since the hulls are disjoint). Put \(\hat{\mathbf{w}}=\mathbf{p}-\mathbf{q}\), \(w_0\)\({}=-\hat{\mathbf{w}}^\T(\mathbf{p}+\mathbf{q})/2\). For \(\mathbf{x}\in H_x\) and \(s\in(0,1]\) the point \(\mathbf{p}+s(\mathbf{x}-\mathbf{p})\) lies in \(H_x\), so \(\|\mathbf{p}+s(\mathbf{x}-\mathbf{p})\)\({}-\mathbf{q}\|^2\)\({}\ge\|\mathbf{p}-\mathbf{q}\|^2\), i.e. \(2s(\mathbf{x}-\mathbf{p})^\T(\mathbf{p}-\mathbf{q})\)\({}+s^2\|\mathbf{x}-\mathbf{p}\|^2\)\({}\ge0\). Divide by \(s\) and let \(s\to0\): \(\hat{\mathbf{w}}^\T(\mathbf{x}-\mathbf{p})\ge0\). Hence \(\hat{\mathbf{w}}^\T\mathbf{x}+w_0\)\({}\ge\hat{\mathbf{w}}^\T\mathbf{p}+w_0\)\({}=\tfrac12\|\hat{\mathbf{w}}\|^2>0\). The same argument with the roles exchanged gives \(\hat{\mathbf{w}}^\T\mathbf{y}+w_0\)\({}\le-\tfrac12\|\hat{\mathbf{w}}\|^2<0\) for \(\mathbf{y}\in H_y\). In particular the finite point sets are separated. ∎

2. Exercise 4.2 (★★, www). Suppose all training targets satisfy a linear constraint \(\mathbf{a}^\T\mathbf{t}_n+b=0\) (4.157). Show that the least-squares prediction (4.17) satisfies \(\mathbf{a}^\T\mathbf{y}(\mathbf{x})+b=0\) (4.158), assuming a constant basis function \(\phi_0=1\).

   Solution

   This is Result 4.3. In short: the constraint reads \(\mathbf{T}\mathbf{a}=-b\mathbf{1}\); then \(\mathbf{a}^\T\mathbf{y}(\mathbf{x})\)\({}=(\mathbf{T}\mathbf{a})^\T\widetilde{\mathbf{X}}(\widetilde{\mathbf{X}}^\T\widetilde{\mathbf{X}})^{-1}\tilde{\mathbf{x}}\)\({}=-b\,\mathbf{1}^\T\widetilde{\mathbf{X}}(\widetilde{\mathbf{X}}^\T\widetilde{\mathbf{X}})^{-1}\tilde{\mathbf{x}}\), and because \(\mathbf{1}=\widetilde{\mathbf{X}}\mathbf{e}_0\) (the first column of \(\widetilde{\mathbf{X}}\)), \(\mathbf{1}^\T\widetilde{\mathbf{X}}(\widetilde{\mathbf{X}}^\T\widetilde{\mathbf{X}})^{-1}\)\({}=\mathbf{e}_0^\T\) and \(\mathbf{e}_0^\T\tilde{\mathbf{x}}=1\). So \(\mathbf{a}^\T\mathbf{y}=-b\).

3. Exercise 4.3 (★★). Extend Exercise 4.2 to several simultaneous linear constraints.

   Solution

   Suppose \(\mathbf{a}_i^\T\mathbf{t}_n+b_i=0\) for \(i=1,\dots,I\) and all \(n\). The proof of Result 4.3 uses only one constraint at a time and the fixed matrix \(\widetilde{\mathbf{X}}(\widetilde{\mathbf{X}}^\T\widetilde{\mathbf{X}})^{-1}\tilde{\mathbf{x}}\), which does not depend on the constraint. Applying it for each \(i\) gives \(\mathbf{a}_i^\T\mathbf{y}(\mathbf{x})+b_i=0\) for all \(i\). In matrix form: if \(\mathbf{A}\mathbf{t}_n+\mathbf{b}=\mathbf{0}\) for all \(n\) then \(\mathbf{A}\mathbf{y}(\mathbf{x})+\mathbf{b}\)\({}=\mathbf{0}\) for all \(\mathbf{x}\), that is, the prediction lies in the same affine subspace as the targets.

4. Exercise 4.4 (★, www). Show that maximising the class separation (4.22) with respect to \(\mathbf{w}\), under \(\mathbf{w}^\T\mathbf{w}=1\) with a Lagrange multiplier, gives \(\mathbf{w}\propto\mathbf{m}_2-\mathbf{m}_1\). (The errata change the pointer from (4.23) to (4.22).)

   Solution

   Stationarity of \(\mathbf{w}^\T(\mathbf{m}_2-\mathbf{m}_1)\)\({}+\lambda(\mathbf{w}^\T\mathbf{w}-1)\) gives \(\mathbf{m}_2-\mathbf{m}_1+2\lambda\mathbf{w}\)\({}=\mathbf{0}\), so \(\mathbf{w}\propto\mathbf{m}_2-\mathbf{m}_1\); details in Filling the gap 4.2.

5. Exercise 4.5 (★). Using (4.20), (4.23), (4.24), show that the Fisher criterion (4.25) can be written as (4.26).

   Solution

   \((m_2-m_1)^2\)\({}=\{\mathbf{w}^\T(\mathbf{m}_2-\mathbf{m}_1)\}^2\)\({}=\mathbf{w}^\T\mathbf{S}_\mathrm{B}\mathbf{w}\) and \(s_k^2\)\({}=\sum_{n\in\mathcal{C}_k}\{\mathbf{w}^\T(\mathbf{x}_n-\mathbf{m}_k)\}^2\)\({}=\mathbf{w}^\T\big[\sum_{n\in\mathcal{C}_k}(\mathbf{x}_n-\mathbf{m}_k)(\mathbf{x}_n-\mathbf{m}_k)^\T\big]\mathbf{w}\); summing over the two classes gives \(\mathbf{w}^\T\mathbf{S}_\mathrm{W}\mathbf{w}\) (see the text after Definition 4.2).

6. Exercise 4.6 (★). With the definitions (4.27), (4.28), the results (4.34), (4.36) and the target coding \(N/N_1\), \(-N/N_2\), show that the least-squares solution (4.33) can be written in the form (4.37).

   Solution

   Substitute \(w_0=-\mathbf{w}^\T\mathbf{m}\) in (4.33): \(\sum_n\mathbf{x}_n(\mathbf{x}_n-\mathbf{m})^\T\mathbf{w}\)\({}=\sum_nt_n\mathbf{x}_n\)\({}=N(\mathbf{m}_1-\mathbf{m}_2)\). The left matrix is the total scatter, which splits as \(\mathbf{S}_\mathrm{W}\)\({}+\frac{N_1N_2}N\mathbf{S}_\mathrm{B}\) (the algebra with \(\mathbf{m}_1-\mathbf{m}\)\({}=\frac{N_2}N(\mathbf{m}_1-\mathbf{m}_2)\), \(\mathbf{m}_2-\mathbf{m}\)\({}=\frac{N_1}N(\mathbf{m}_2-\mathbf{m}_1)\)), giving (4.37). Full proof, and the sign of the proportionality constant, in Result 4.5.

7. Exercise 4.7 (★, www). Show that \(\sigma(-a)=1-\sigma(a)\) and that the inverse of \(\sigma\) is \(\sigma^{-1}(y)=\ln\{y/(1-y)\}\).

   Solution

   \(\sigma(-a)\)\({}=\frac1{1+e^a}\)\({}=\frac{e^{-a}}{e^{-a}+1}\)\({}=1-\frac1{1+e^{-a}}\). If \(y=\frac1{1+e^{-a}}\) then \(e^{-a}=\frac{1-y}y\), so \(a=\ln\frac y{1-y}\). See Result 4.8(b).

8. Exercise 4.8 (★). Using (4.57), (4.58), derive (4.65) for the two-class Gaussian generative model and verify (4.66), (4.67).

   Solution

   Write \(p(\mathcal{C}_1\given\mathbf{x})=\sigma(a)\) with \(a\)\({}=\ln\frac{p(\mathbf{x}\given\mathcal{C}_1)p(\mathcal{C}_1)}{p(\mathbf{x}\given\mathcal{C}_2)p(\mathcal{C}_2)}\). With the shared-covariance Gaussians the normalisers cancel, and expanding the two quadratic forms the terms \(\mp\frac12\mathbf{x}^\T\bm\Sigma^{-1}\mathbf{x}\) cancel, leaving \(a\)\({}=\mathbf{x}^\T\bm\Sigma^{-1}(\bm\mu_1-\bm\mu_2)\)\({}-\frac12\bm\mu_1^\T\bm\Sigma^{-1}\bm\mu_1\)\({}+\frac12\bm\mu_2^\T\bm\Sigma^{-1}\bm\mu_2\)\({}+\ln\frac{p(\mathcal{C}_1)}{p(\mathcal{C}_2)}\). Full details in Result 4.9.

9. Exercise 4.9 (★, www). A generative model for \(K\) classes has priors \(p(\mathcal{C}_k)=\pi_k\) and arbitrary class-conditional densities \(p(\bm\phi\given\mathcal{C}_k)\); the data are i.i.d. with 1-of-\(K\) targets \(t_{nj}=I_{jk}\) if \(n\in\mathcal{C}_k\). Show that the maximum-likelihood priors are \(\pi_k=N_k/N\) (4.159).

   Solution

   The joint probability of \((\bm\phi_n,\mathbf{t}_n)\) is \(\prod_k\{\pi_kp(\bm\phi_n\given\mathcal{C}_k)\}^{t_{nk}}\), so the log likelihood is \(\sum_n\sum_kt_{nk}\ln\pi_k\)\({}+\sum_n\sum_kt_{nk}\ln p(\bm\phi_n\given\mathcal{C}_k)\). Only the first term involves \(\bm\pi\): maximise \(\sum_kN_k\ln\pi_k\) subject to \(\sum_k\pi_k=1\). The Lagrangian condition \(N_k/\pi_k+\lambda=0\) gives \(\pi_k=-N_k/\lambda\), and summing over \(k\) gives \(\lambda=-N\), hence \(\pi_k=N_k/N\). It is a maximum since \(\sum N_k\ln\pi_k\) is concave. (Result 4.10.)

10. Exercise 4.10 (★★). In the model of Exercise 4.9 let \(p(\bm\phi\given\mathcal{C}_k)\)\({}=\N(\bm\phi\given\bm\mu_k,\bm\Sigma)\). Show that the maximum-likelihood means are \(\bm\mu_k=\frac1{N_k}\sum_nt_{nk}\bm\phi_n\) (4.161) and the shared covariance is \(\bm\Sigma=\sum_k\frac{N_k}N\mathbf{S}_k\) with \(\mathbf{S}_k\)\({}=\frac1{N_k}\sum_nt_{nk}(\bm\phi_n-\bm\mu_k)(\bm\phi_n-\bm\mu_k)^\T\) (4.162, 4.163).

    Solution

    The terms in \(\bm\mu_k\) are \(-\frac12\sum_nt_{nk}(\bm\phi_n-\bm\mu_k)^\T\bm\Sigma^{-1}(\bm\phi_n-\bm\mu_k)\); setting the gradient \(\bm\Sigma^{-1}\sum_nt_{nk}(\bm\phi_n-\bm\mu_k)\) to zero gives \(\bm\mu_k=\frac1{N_k}\sum_nt_{nk}\bm\phi_n\). For \(\bm\Sigma\), the log likelihood part is \(-\frac N2\ln|\bm\Sigma|\)\({}-\frac N2\operatorname{Tr}(\bm\Sigma^{-1}\mathbf{S})\) with \(\mathbf{S}\)\({}=\frac1N\sum_{n,k}t_{nk}(\bm\phi_n-\bm\mu_k)(\bm\phi_n-\bm\mu_k)^\T\)\({}=\sum_k\frac{N_k}N\mathbf{S}_k\); its maximiser is \(\bm\Sigma=\mathbf{S}\) (derivative \(-\frac N2\bm\Sigma^{-1}\)\({}+\frac N2\bm\Sigma^{-1}\mathbf{S}\bm\Sigma^{-1}\)). The covariance is thus the prior-weighted average of the class covariances. (Result 4.10.)

11. Exercise 4.11 (★★). A \(K\)-class problem has a feature vector \(\bm\phi\) of \(M\) components, each with \(L\) discrete states coded 1-of-\(L\); the components are independent given the class. Show that the \(a_k\) of (4.63) are linear functions of the components of \(\bm\phi\).

    Solution

    \(p(\bm\phi\given\mathcal{C}_k)\)\({}=\prod_{i=1}^M\prod_{l=1}^L\mu_{kil}^{\phi_{il}}\), so \(a_k\)\({}=\sum_{i,l}\phi_{il}\ln\mu_{kil}\)\({}+\ln p(\mathcal{C}_k)\), linear in the \(\phi_{il}\). See Filling the gap 4.4, which also gives the ML estimates.

12. Exercise 4.12 (★, www). Verify \(\d\sigma/\d a=\sigma(1-\sigma)\) (4.88).

    Solution

    \(\dfrac{\d}{\d a}(1+e^{-a})^{-1}\)\({}=\dfrac{e^{-a}}{(1+e^{-a})^2}\)\({}=\dfrac1{1+e^{-a}}\cdot\dfrac{e^{-a}}{1+e^{-a}}\)\({}=\sigma(a)(1-\sigma(a))\), since \(1-\sigma=\frac{e^{-a}}{1+e^{-a}}\).

13. Exercise 4.13 (★, www). Use (4.88) to show that the derivative of the cross-entropy (4.90) is (4.91).

    Solution

    See Result 4.12: \(\nabla E\)\({}=-\sum_n\{\frac{t_n}{y_n}\)\({}-\frac{1-t_n}{1-y_n}\}y_n(1-y_n)\bm\phi_n\)\({}=\sum_n(y_n-t_n)\bm\phi_n\).

14. Exercise 4.14 (★). Show that for linearly separable data the maximum-likelihood solution of logistic regression is obtained by finding a separating \(\mathbf{w}\) and taking its magnitude to infinity.

    Solution

    Result 4.13: the error is strictly positive for every finite \(\mathbf{w}\), but along \(c\hat{\mathbf{w}}\), with \(\hat{\mathbf{w}}\) separating, each \(\sigma(c\hat a_n)\) tends to 1 for \(t_n=1\) and to 0 for \(t_n=0\), so \(E(c\hat{\mathbf{w}})\to0\) as \(c\to\infty\). The infimum 0 is therefore approached only at infinity. Conversely, every sequence along which \(E\to0\) must have all \(y_n\to t_n\), i.e. eventually has \(\mathbf{w}^\T\bm\phi_n\) of the correct sign for every \(n\), so the iterates eventually separate the data and \(\|\mathbf{w}\|\to\infty\).

15. Exercise 4.15 (★★). Show that the logistic-regression Hessian \(\mathbf{H}=\bm\Phi^\T\mathbf{R}\bm\Phi\) (4.97), \(R_{nn}=y_n(1-y_n)\), is positive definite, and hence that \(E\) is convex with a unique minimum. (The errata correct "concave" to "convex".)

    Solution

    For any \(\mathbf{u}\), \(\mathbf{u}^\T\mathbf{H}\mathbf{u}\)\({}=\sum_ny_n(1-y_n)(\mathbf{u}^\T\bm\phi_n)^2\)\({}\ge0\), strictly positive for \(\mathbf{u}\ne\mathbf{0}\) when \(\bm\Phi\) has full column rank (otherwise zero for \(\mathbf{u}\) in its null space). So \(E\) is convex, strictly convex when \(\operatorname{rank}\bm\Phi=M\), and a stationary point is then the unique minimum; for separable data none exists. See Result 4.14.

16. Exercise 4.16 (★). Each observation \(\mathbf{x}_n\) has an uncertain label: instead of \(t_n\in\{0,1\}\) we are given \(\pi_n=P(t_n=1)\) (errata: \(t\) is \(t_n\)). For a model \(p(t=1\given\bm\phi)=y(\bm\phi)\) write down the appropriate log likelihood.

    Solution

    If the true label of point \(n\) were known, its log likelihood contribution would be \(t_n\ln y_n+(1-t_n)\ln(1-y_n)\). The label is a Bernoulli variable with \(P(t_n=1)=\pi_n\), so the expected log likelihood is

    \[ \sum_{n=1}^N\{\pi_n\ln y(\bm\phi_n)+(1-\pi_n)\ln(1-y(\bm\phi_n))\}. \]

    This is the cross-entropy (4.90) with the soft target \(\pi_n\) in place of \(t_n\); maximising it minimises \(\sum_n\mathrm{KL}(\mathrm{Bern}(\pi_n)\|\mathrm{Bern}(y_n))\), and it is exactly the expected complete-data log likelihood (the E-step quantity) when the label is a latent variable with posterior \(\pi_n\). The gradient is \(\sum_n(y_n-\pi_n)\bm\phi_n\), so IRLS applies unchanged with \(t_n\to\pi_n\).

17. Exercise 4.17 (★, www). Show that the derivatives of the softmax (4.104) with respect to the activations are (4.106).

    Solution

    Result 4.15(a): \(\partial y_k/\partial a_j\)\({}=y_kI_{kj}-y_ky_j\)\({}=y_k(I_{kj}-y_j)\).

18. Exercise 4.18 (★). Using (4.106) (errata: the book says (4.91)), show that the gradients of the cross-entropy error (4.108) are (4.109).

    Solution

    Result 4.15(b): \(\partial E/\partial a_{nj}\)\({}=-\sum_kt_{nk}(I_{kj}-y_{nj})\)\({}=y_{nj}-t_{nj}\) because \(\sum_kt_{nk}=1\); multiply by \(\nabla_{\mathbf{w}_j}a_{nj}=\bm\phi_n\) and sum over \(n\).

19. Exercise 4.19 (★, www). Write down the gradient and the Hessian of the log likelihood of the probit regression model of Section 4.3.5, the quantities needed by IRLS.

    Solution

    With \(y_n=\Phi(\mathbf{w}^\T\bm\phi_n)\), \(a_n=\mathbf{w}^\T\bm\phi_n\), \(\varphi=\N(\cdot\given0,1)\) and \(E=-\ln p(\mathbf{t}\given\mathbf{w})\):

    \[ \nabla E=\sum_n\frac{y_n-t_n}{y_n(1-y_n)}\varphi(a_n)\bm\phi_n,\qquad \mathbf{H}=\sum_n\Big\{\frac{\varphi^2}{y(1-y)}-(y-t)\Big[\frac{a\varphi}{y(1-y)}+\frac{(1-2y)\varphi^2}{y^2(1-y)^2}\Big]\Big\}_{n}\bm\phi_n\bm\phi_n^\T. \]

    Derivation in Filling the gap 4.7 (checked against finite differences). The log likelihood has the opposite sign: gradient \(-\nabla E\), Hessian \(-\mathbf{H}\).

20. Exercise 4.20 (★★). Show that the Hessian of multiclass logistic regression (4.110) is positive semi-definite, by considering \(\mathbf{u}^\T\mathbf{H}\mathbf{u}\) for an arbitrary \(\mathbf{u}\in\R^{MK}\) and applying Jensen's inequality.

    Solution

    Splitting \(\mathbf{u}=(\mathbf{u}_1,\dots,\mathbf{u}_K)\) and \(b_{nk}=\mathbf{u}_k^\T\bm\phi_n\), \(\mathbf{u}^\T\mathbf{H}\mathbf{u}\)\({}=\sum_n\{\sum_ky_{nk}b_{nk}^2\)\({}-(\sum_ky_{nk}b_{nk})^2\}\). For each \(n\) the bracket is \(\mathbb{E}[b^2]-\mathbb{E}[b]^2\ge0\) under the probability vector \((y_{nk})_k\), which is Jensen's inequality for \(x\mapsto x^2\). It vanishes for \(\mathbf{u}_1=\dots=\mathbf{u}_K\), so \(\mathbf{H}\) is semi-definite but singular (Result 4.15(d)).

21. Exercise 4.21 (★). Show that the inverse probit function (4.114) and the erf function (4.115) are related by (4.116).

    Solution

    With the corrected \(\operatorname{erf}(a)\)\({}=\frac2{\sqrt\pi}\int_0^ae^{-\theta^2}\d\theta\): \(\Phi(a)\)\({}=\frac12\)\({}+\int_0^a\frac{e^{-\theta^2/2}}{\sqrt{2\pi}}\d\theta\), and \(\theta=\sqrt2u\) turns the integral into \(\frac1{\sqrt\pi}\int_0^{a/\sqrt2}e^{-u^2}\d u\)\({}=\frac12\operatorname{erf}(a/\sqrt2)\). Hence \(\Phi(a)\)\({}=\frac12\{1+\operatorname{erf}(a/\sqrt2)\}\) (Filling the gap 4.7).

22. Exercise 4.22 (★). Using (4.135), derive the log evidence (4.137).

    Solution

    Apply (4.135) to \(f(\bm\theta)\)\({}=p(\mathcal{D}\given\bm\theta)p(\bm\theta)\), \(Z=p(\mathcal{D})\), and take logs: \(\ln p(\mathcal{D})\)\({}\simeq\ln p(\mathcal{D}\given\bm\theta_\mathrm{MAP})\)\({}+\ln p(\bm\theta_\mathrm{MAP})+\frac M2\ln2\pi\)\({}-\frac12\ln|\mathbf{A}|\) (Result 4.18).

23. Exercise 4.23 (★★, www). Starting from the Laplace evidence (4.137), show that for a Gaussian prior \(\N(\bm\theta\given\mathbf{m},\mathbf{V}_0)\) it takes the form \(\ln p(\mathcal{D})\)\({}\simeq\ln p(\mathcal{D}\given\bm\theta_\mathrm{MAP})\)\({}-\frac12(\bm\theta_\mathrm{MAP}-\mathbf{m})^\T\mathbf{V}_0^{-1}(\bm\theta_\mathrm{MAP}-\mathbf{m})\)\({}-\frac12\ln|\mathbf{H}|+\mathrm{const}\), \(\mathbf{H}\) being the Hessian of the negative log likelihood at \(\bm\theta_\mathrm{MAP}\) (errata: "negative"); then derive the BIC (4.139) for a broad prior and i.i.d. data.

    Solution

    See Result 4.19. Briefly: \(\mathbf{A}=\mathbf{H}+\mathbf{V}_0^{-1}\) and \(\ln p(\bm\theta_\mathrm{MAP})\)\({}=-\frac12(\cdot)^\T\mathbf{V}_0^{-1}(\cdot)\)\({}-\frac12\ln|\mathbf{V}_0|-\frac M2\ln2\pi\); so (4.137) becomes \(\ln p(\mathcal{D}\given\bm\theta_\mathrm{MAP})\)\({}-\frac12(\cdot)^\T\mathbf{V}_0^{-1}(\cdot)\)\({}-\frac12\ln|\mathbf{V}_0|\)\({}-\frac12\ln|\mathbf{H}\)\({}+\mathbf{V}_0^{-1}|\). Dropping the \(\mathbf{V}_0^{-1}\) terms for a broad prior leaves the stated form (with \(-\frac12\ln|\mathbf{V}_0|\) in the constant). For i.i.d. data \(\mathbf{H}=N\bar{\mathbf{H}}\), \(\ln|\mathbf{H}|=M\ln N+\ln|\bar{\mathbf{H}}|\), and keeping the term that grows with \(N\) gives \(-\frac12M\ln N\).

24. Exercise 4.24 (★★). Use the results of Section 2.3.2 to derive the marginalisation (4.151) over a Gaussian posterior.

    Solution

    Under \(q(\mathbf{w})\)\({}=\N(\mathbf{w}_\mathrm{MAP},\mathbf{S}_N)\), the scalar \(a=\bm\phi^\T\mathbf{w}\) is a linear function of a Gaussian vector, hence Gaussian (Section 2.3.2; the marginal of the linear-Gaussian model \(a\given\mathbf{w}=\bm\phi^\T\mathbf{w}\), noise-free) with mean \(\bm\phi^\T\mathbf{w}_\mathrm{MAP}\) and variance \(\bm\phi^\T\mathbf{S}_N\bm\phi\). Since \(\sigma(\mathbf{w}^\T\bm\phi)\) depends on \(\mathbf{w}\) only through \(a\), \(\mathbb{E}_q[\sigma(\mathbf{w}^\T\bm\phi)]\)\({}=\int\sigma(a)\N(a\given\mu_a,\sigma_a^2)\d a\). (Result 4.21.)

25. Exercise 4.25 (★★). Approximate \(\sigma(a)\) by a scaled inverse probit \(\Phi(\lambda a)\). Show that choosing \(\lambda\) so that the derivatives agree at \(a=0\) gives \(\lambda^2=\pi/8\).

    Solution

    \(\sigma'(0)=\frac14\) and \(\frac{\d}{\d a}\Phi(\lambda a)|_0\)\({}=\lambda/\sqrt{2\pi}\). Setting them equal gives \(\lambda=\sqrt{2\pi}/4\), \(\lambda^2=\pi/8\) (Result 4.22). Numerically \(\lambda=0.6267\).

26. Exercise 4.26 (★★). Prove the convolution formula (4.152) by differentiating both sides with respect to \(\mu\), integrating, and showing that the constant of integration vanishes.

    Solution

    See Result 4.23: with \(a=\mu+\sigma z\) the derivative of the left side is a Gaussian integral over \(z\) which, by completing the square with \(c=1+\lambda^2\sigma^2\), equals \(\frac\lambda{\sqrt c}\N(\frac{\lambda\mu}{\sqrt c}\given0,1)\); this is the derivative of the right side. Both sides tend to 0 as \(\mu\to-\infty\), so the constant is zero.

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

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 4, pp. 179–224 (and the official errata).
2. Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. doi:10.1109/TAC.1974.1100705.
3. Block, H. D. (1962). The perceptron: a model for brain functioning. I. Reviews of Modern Physics, 34(1), 123–135. doi:10.1103/RevModPhys.34.123.
4. Duda, R. O. and Hart, P. E. (1973). Pattern Classification and Scene Analysis. Wiley, New York.
5. Fan, K. (1949). On a theorem of Weyl concerning eigenvalues of linear transformations I. Proceedings of the National Academy of Sciences, 35(11), 652–655. doi:10.1073/pnas.35.11.652. (The trace inequality used in the proof of Result 4.6.)
6. Fisher, R. A. (1936). The use of multiple measurements in taxonomic problems. Annals of Eugenics, 7(2), 179–188. doi:10.1111/j.1469-1809.1936.tb02137.x.
7. Fukunaga, K. (1990). Introduction to Statistical Pattern Recognition, 2nd edition. Academic Press, San Diego.
8. MacKay, D. J. C. (1992). The evidence framework applied to classification networks. Neural Computation, 4(5), 720–736. doi:10.1162/neco.1992.4.5.720.
9. Minsky, M. L. and Papert, S. A. (1969). Perceptrons. MIT Press, Cambridge, MA.
10. Nelder, J. A. and Wedderburn, R. W. M. (1972). Generalized linear models. Journal of the Royal Statistical Society: Series A, 135(3), 370–384. doi:10.2307/2344614.
11. Novikoff, A. B. J. (1962). On convergence proofs on perceptrons. In Proceedings of the Symposium on the Mathematical Theory of Automata, Vol. 12, 615–622. Polytechnic Institute of Brooklyn.
12. Opper, M. and Winther, O. (2000). Gaussian processes for classification: mean-field algorithms. Neural Computation, 12(11), 2655–2684. doi:10.1162/089976600300014881.
13. Rosenblatt, F. (1958). The perceptron: a probabilistic model for information storage and organization in the brain. Psychological Review, 65(6), 386–408. doi:10.1037/h0042519.
14. Rubin, D. B. (1983). Iteratively reweighted least squares. In Encyclopedia of Statistical Sciences, Vol. 4, 272–275. Wiley.
15. Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. doi:10.1214/aos/1176344136.
16. Spiegelhalter, D. J. and Lauritzen, S. L. (1990). Sequential updating of conditional probabilities on directed graphical structures. Networks, 20(5), 579–605. doi:10.1002/net.3230200507.
17. Widrow, B. and Hoff, M. E. (1960). Adaptive switching circuits. In IRE WESCON Convention Record, Part 4, 96–104.

DOIs were checked against Crossref. The numerical statements made in the text (the equivalence of the least-squares and Fisher directions including the constant \(\kappa\), the masking example, the probit gradient and Hessian against finite differences, the multiclass Hessian eigenvalues, the Laplace approximation of Figure 4.14, the accuracy of the probit approximation (4.153) and the Newton iterates on separable data) come from short Python computations written for these notes; the figures come from the course code.