3 Supervised Learning, Prediction, and Evaluation

The most common loss for regression is squared loss,

\begin{align*} \lscr(y,f(x)) &= \frac{1}{2}(y-f(x))^2, & \lscr \in [0,\infty). \end{align*}

Clearly, this loss penalizes large prediction errors because then \((Y-f(X))^2\) is large. Likelihood maximization offers another view. Suppose, as an approximation, that conditional on the feature vectors, the labels \(Y_i\) are independent normal random variables with means \(f(X_i)\) and common variance \(v>0\):

\begin{align*} Y_i \mid X_i \sim \Norm{f(X_i),v}. \end{align*}

Then the likelihood of the observed labels in \(\Tset\) is

\begin{align*} L(f,v) &= \prod_{(x,y)\in\Tset} \frac{1}{\sqrt{2\pi v}} \exp\rb{-\frac{(y-f(x))^2}{2v}}. \end{align*}

Taking the negative logarithm¹⁶¹⁶Unless a base is shown explicitly, \(\log\) means the natural logarithm throughout these notes. to convert the product into a sum gives

\begin{align*} -\log L(f,v) &= \frac{n}{2}\log(2\pi v) + \frac{1}{2v}\sum_{(x,y)\in\Tset} (y-f(x))^2. \end{align*}

For fixed \(v>0\), maximizing the likelihood over \(f\) is therefore the same as minimizing the empirical squared loss.

To find the optimal rule for square loss, suppose we know \(F_{X,Y}\) and may choose any prediction function. As a first step, we decompose the risk in a way that will be useful throughout this section. It separates the part of the error that can be reduced by choosing a better prediction rule from the variation in \(Y\) that remains once \(X\) is known.

Proof

Write the prediction error \(Y-f(X)\) as \((Y-m(X))+(m(X)-f(X))\):

\begin{align*} \frac{1}{2}\E{(Y-f(X))^{2}} &= \frac{1}{2}\E{\rb{\rb{Y-m(X)} + \rb{m(X)-f(X)}}^{2}} \\ &= \frac{1}{2}\E{\rb{m(X)-f(X)}^{2}} + \E{\rb{m(X)-f(X)}\rb{Y-m(X)}}\\ &\quad + \frac{1}{2}\E{\rb{Y-m(X)}^{2}}. \end{align*}

By conditioning on \(X\) in the middle term,

\begin{align*} \E{\rb{m(X)-f(X)}\rb{Y-m(X)} \given X} = \rb{m(X)-f(X)} \E{Y-m(X) \given X}= 0, \end{align*}

because \(m(X) = \E{Y \given X}\).

With respect to the third term, conditioning on \(X\) yields

\begin{align*} \E{\rb{Y-m(X)}^{2}\given X} &= \E{Y^{2}\given X} - 2m(X)\E{Y\given X} + m(X)^{2}\\ &= \E{Y^{2}\given X} - m(X)^{2} &= \V{Y\mid X}. \end{align*}

The last equation follows from the definition of conditional variance

\begin{align*} \V{Y\mid X} &= \E{(Y-\E{Y\mid X})^{2}\mid X} = \E{Y^{2}\mid X} - m(X)^{2}. \end{align*}

Taking expectations completes the proof.

The Bayes regressor¹⁷¹⁷That is, the Bayes optimal predictor (3.1.4) for a regression problem. under squared loss now follows easily.

Fitting the best constant function requires minimizing \(\frac{1}{2}\E{(m(X) -c)^2}\) in (3.2.1). Expanding the square and using that \(\E{m(X)} = \E Y\), gives

\begin{align*} \frac{1}{2}\E{m(X)^{2}} +\frac{1}{2}c (c-2 \E Y). \end{align*}

This is minimal at \(c^* = \E Y\), and then the risk \(R(c^*)\) becomes

\begin{align*} c^{*} &= \E Y, & R(c^*) &= \frac{1}{2}\E{\V{Y\mid X}} + \frac{1}{2}\E{m(X)^{2}} - \frac{1}{2}(\E Y)^{2} = \frac{1}{2}\V Y, \tag{3.2.2} \end{align*}

where the last step follows from Eve’s law.

To obtain the empirical risk minimizer and minimal empirical risk, we replace \(F_{X,Y}\) by the probability measure induced by \(\Tset\) to get

\begin{align*} \hat c &= \frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}y, & \hat R_{\Tset}(\hat c) &= \frac{1}{2|\Tset|}\sum_{(x,y)\in\Tset}(y-\hat c)^2. \tag{3.2.3} \end{align*}

Fitting an affine prediction rule minimizes the excess-risk term in (3.2.1):

\begin{align*} \frac{1}{2}\E{\rb{m(X)-\beta(X)}^2}. \end{align*}

Setting the partial derivatives wrt \(\beta\) and \(\beta_0\) to zero gives the population normal equations:

\begin{align*} \beta(\E{X}) &= \E Y,\\ \E{\rb{m(X)-\beta(X)}X(j)} &= 0, \qquad j=1,\ldots,p. \end{align*}

Replacing \(F_{X,Y}\) by the normalized counting measure induced by \(\Tset\) and writing \(\hat\beta(x)=\hat\beta_0+\ip{\hat\beta,x}\) for the fitted rule gives

\begin{align*} \hat\beta(\hat x) &= \sum_{(x,y)\in \Tset} y, \\ \sum_{(x,y)\in\Tset}\rb{y-\hat\beta(x)}x(j) &= 0, \qquad j=1,\ldots,p. \end{align*}

These solve the ordinary least squares (OLS) problem:

\begin{align*} (\hat\beta_0,\hat\beta) &=\argmin_{\beta_0,\beta}\hat R_{\Tset}(\beta(\cdot)), \tag{3.2.4} \\ \hat R_{\Tset}(\beta(\cdot)) &=\frac{1}{2|\Tset|}\sum_{(x,y)\in\Tset} \rb{y-\beta(x)}^2. \tag{3.2.5} \end{align*}

Restricting the affine class to \(\beta=0\) recovers the constant class: the population equation gives \(\beta_0^*=\E{Y}\) and the sample equation gives \(\hat\beta_0=\hat c\) of (3.2.3).

Bias and variance are two forms of error that appear through the approximation of the Bayes predictor \(m\) by the empirical risk minimizer \(\hat f_{\Fset}\), (3.1.5). There are two sources of randomness to distinguish. The training set \(\Tset\) can itself be seen as a list of iid samples from \(F_{X,Y}\). Once \(\Tset\) is known and the rule \(\hat f_{\Fset}\) is fitted, a new observation \((X,Y)\), independent of \(\Tset\) and also with distribution \(F_{X,Y}\), offers the second source of randomness. Let us fix \(X =x\) first, and concentrate on the effect of the randomnes s of \(\Tset\). Write¹⁹¹⁹A subscript on an expectation indicates which random variable is averaged over; a vertical bar indicates conditioning.

\begin{equation*} \bar f(x)=\EE{\Tset}{\hat f_{\Fset}(x)} \tag{3.2.6} \end{equation*}

for the average prediction at \(x\) over the iid training sets. The bias \[ \bar f(x)-m(x) \] is the systematic discrepancy between the average fitted prediction and the conditional mean.²⁰²⁰Recall from Cor. 3.2.2 that \(m(x)\) is precisely the Bayes regressor. The variance

\begin{align*} \EE{\Tset}{\rb{\bar f(x) - \hat f_{\Fset}(x)}^{2}} \end{align*}

measures how the fitted prediction changes when we draw a new training sample.

Proof

We use (3.2.1). Conditional on \(\Tset\), the fitted rule is fixed, so that decomposition gives

\begin{align*} R(\hat f_{\Fset}) &= \frac{1}{2}\EE{X}{\V{Y\mid X}} + \frac{1}{2}\EE{X}{(m(X)-\hat f_{\Fset}(X))^2\mid\Tset}. \end{align*}

Averaging over training sets and exchanging the expectations in the second term yields

\begin{align*} \EE{\Tset}{R(\hat f_{\Fset})} &= \frac{1}{2}\EE{X}{\V{Y\mid X}} + \frac{1}{2}\EE{X}{\EE{\Tset}{(m(X)-\hat f_{\Fset}(X))^2\mid X}}. \end{align*}

Note that the exchange \(\EE{\Tset}{\EE{X}{\cdot}}\) to \(\EE{X}{\EE{\Tset}{\cdot}}\) is valid by independence and Tonelli’s theorem, since the squared deviation is nonnegative.

For fixed \(X\), expand the second term around \(\bar f(X)=\EE{\Tset}{\hat f_{\Fset}(X)\mid X}\):

\begin{align*} m(X)-\hat f_{\Fset}(X) &= (m(X)-\bar f(X))+(\bar f(X)-\hat f_{\Fset}(X)). \end{align*}

Squaring and taking the conditional expectation given \(X\) gives

\begin{align*} \EE{\Tset}{(m(X)-\hat f_{\Fset}(X))^2\mid X} &= (m(X)-\bar f(X))^2 \\ &\quad +2(m(X)-\bar f(X)) \EE{\Tset}{\bar f(X)-\hat f_{\Fset}(X)\mid X} \\ &\quad +\EE{\Tset}{(\bar f(X)-\hat f_{\Fset}(X))^2\mid X}. \end{align*}

Here \(m(X)\) and \(\bar f(X)\) are fixed conditional on \(X\). The cross term vanishes because

\begin{align*} \EE{\Tset}{\bar f(X)-\hat f_{\Fset}(X)\mid X} &= \bar f(X)-\EE{\Tset}{\hat f_{\Fset}(X)\mid X}=0. \end{align*}

The last term is the conditional variance of the fitted prediction, so

\begin{align*} \EE{\Tset}{(m(X)-\hat f_{\Fset}(X))^2\mid X} &= (m(X)-\bar f(X))^2 + \V{\hat f_{\Fset}(X)\mid X}. \end{align*}

Substituting this identity into the expression for \(\EE{\Tset}{R(\hat f_{\Fset})}\) above gives the claimed decomposition.

Feature selection affects both bias and variance. Suppose for this paragraph that the feature coordinates are uncorrelated, so that leaving one out does not disturb the estimates of the others. A feature with \(\beta_{j}=0\) has no contribution to the conditional mean, so estimating its coefficient adds variance without reducing bias. Omitting it therefore lowers the variance of the fitted rule without increasing bias. On the other hand, when we leave out a feature with \(\beta_{j}\neq0\) its contribution \(\beta_{j}x(j)\) to the residual

\begin{align*} y-\hat\beta_0-\sum_{k\neq j}\hat\beta_k x(k) \end{align*}

must be picked up by the other factors. This introduces bias and raises the variance of the remaining estimates.

Feature selection aims to retain useful features and omit the ones that add only estimation noise. The lasso in Section 3.5 is one method for selecting relevant features.

An alternative loss to squared loss is the absolute loss

\begin{align*} \lscr(y,f(x)) &= |y-f(x)|, & \lscr \in [0,\infty). \tag{3.2.8} \end{align*}

It is less sensitive to outliers than squared loss because it grows linearly with the magnitude of the error. Its Bayes regressor is the conditional median, as the following definition and theorem show.

To see why medians arise, consider the slopes of the graphs in Fig. 3.1--Fig. 3.4. Viewed as a function of \(c\), a term \(|y-c|\) has slope \(+1\) when \(y<c\) and slope \(-1\) when \(y>c\). Thus, to get the slope of the sum of such functions, when \(c\) is not equal to any of the \(y\), we can just count the number of labels \(y\) lying to the left and the right of \(c\). When changing \(c\) a tiny bit, all terms \(|y-c|\) with \(y>c\) contribute \(-1\) to the slope of the sum, and the ones with \(y<c\) contribute \(+1\). Clearly, any such sum is convex. When the number of labels is odd, the minimum of the sum is achieved at the median; otherwise it is flat between the two middle labels. For instance, for the two labels 1 and 3, the sum is constant between the two labels, so every \(c\in[1,3]\) is a minimizer. After adding the middle label 2, the sum has its unique minimum at \(c=2\).

⊕
Figure 3.1: The plots of \(|y_1-c|\) and \(|y_2-c|\), with \(y_1=1\) and \(y_2=3\).

⊕
Figure 3.2: The plot of \(|y_1-c|+|y_2-c|\), with \(y_1=1\) and \(y_2=3\).

⊕
Figure 3.3: The plots of \(|y_i-c|\), with \(y_1=1\), \(y_2=2\), and \(y_3=3\).

⊕
Figure 3.4: The plot of \(\sum_{i=1}^3|y_i-c|\), with \(y_1=1\), \(y_2=2\), and \(y_3=3\).

Let us apply this idea to understand the empirical risk of the constant prediction \(f(x)=c\) for a list of observations \(\Sset\):

\begin{align*} \hat R_{\Sset}(c)=\frac{1}{|\Sset|} \sum_{(x,y)\in\Sset}|y-c|. \end{align*}

Define

\begin{align*} \Lset(c)&=\{(x,y)\in\Sset:y\leq c\}, & \Rset(c)&=\Sset\setminus\Lset(c)=\{(x,y)\in\Sset:y>c\}. \end{align*}

It follows that for \(c\notin\set{y:(x,y)\in\Sset}\),

\begin{align*} \hat R'_{\Sset}(c) &=\frac{|\Lset(c)|-|\Rset(c)|}{|\Sset|} =2\frac{|\Lset(c)|}{|\Sset|}-1. \end{align*}

To specify this derivative everywhere on \(\R\), define it at \(c=y\) as²²²²It becomes a right-continuous function with this.

\begin{align*} \hat R'_{\Sset}(y) := \hat R'_{\Sset}(y+) = \lim_{c\downarrow y}\hat R'_{\Sset}(c). \end{align*}

The empirical risk is convex, so it is minimized where its derivative changes from negative to positive, possibly through an interval on which it is zero. At any such point at most half the labels lie strictly below a minimizer and at most half lie strictly above it: these are exactly the sample-median conditions.

The same counting argument, now applied to the conditional distribution of \(Y\) given \(X=x\), gives the Bayes predictor.

Proof

Conditioning on \(X\) gives \(R(f)=\E{\,\E{\abs{Y-f(X)}\given X}}\). Given \(X=x\), the prediction is the constant \(c=f(x)\), so we minimize \(\E{\abs{Y-c}\given X=x}\) over \(c\). Imagine increasingly large lists \(\Sset_n\) whose labels are sampled independently from the conditional distribution of \(Y\) given \(X=x\). By the law of large numbers,²³²³We omit the technical details needed to justify this passage to the limit. the empirical proportion \(|\Lset_n(c)|/|\Sset_n|\) converges to \(\P{Y\leq c\given X=x}\). Consequently, the empirical right slope \(2|\Lset_n(c)|/|\Sset_n|-1\) converges to \(2\P{Y\leq c\given X=x}-1\). The left limit of this limiting slope is \(2\P{Y<c\given X=x}-1\).

The conditional risk \(\E{|Y-c|\mid X}\) is convex and coercive,²⁴²⁴A function \(f\) is coercive if \(f(x)\to\infty\) when \(|x|\to\infty\). so it has a minimizer \(m\). At this minimizer,

\begin{align*} 2\P{Y<m\given X=x}-1 \leq 0 \leq 2\P{Y\leq m\given X=x}-1, \end{align*}

hence at \(m\),

\begin{align*} \P{Y<m\given X=x} \leq \frac 1 2 \leq \P{Y\leq m\given X=x}. \end{align*}

Thus, \(m\) satisfies the definition of the median; \(m=\med{Y\mid X=x}\). Choosing such a conditional median \(m(x)\) for for each \(x\) minimizes \(R(f)\).

Constant functions as predictors have the form \(f(x)=c\), \(c\in\R\). The population and empirical calculations above immediately give

\begin{align*} R(c)&=\E{\abs{Y-c}}, & c^*&=\med{Y}, \end{align*}

and the empirical risk minimizer is

\begin{align*} \hat m&=\med{y:(x,y)\in\Tset}, & \hat m&=\argmin_{c\in\R}\hat R_{\Tset}(c). \tag{3.2.9} \end{align*}

The attained cost has a simple form. Use the same three-way partition as Quicksort Alg. 2.3.2 with the median \(\hat m\) as the pivot:²⁵²⁵A median label is an optimal pivot for balancing Quicksort’s recursive subproblems: both \(\Lset\) and \(\Rset\) contain at most half the observations.

\begin{align*} \Lset&=\set{(x,y)\in\Tset:y<\hat m},& \Mset&=\set{(x,y)\in\Tset:y=\hat m},& \Rset&=\set{(x,y)\in\Tset:y>\hat m}. \end{align*}

The observations in \(\Mset\) contribute zero to the absolute loss. Resolving the remaining absolute values gives

\begin{align*} \hat R_{\Tset}(\hat m) &=\frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}\abs{y-\hat m}\\ &=\frac{1}{|\Tset|}\sum_{(x,y)\in\Lset}\rb{\hat m-y} +\frac{1}{|\Tset|}\sum_{(x,y)\in\Rset}\rb{y-\hat m}\\ &=\frac{1}{|\Tset|}\rb{ \sum_{(x,y)\in\Rset}y-\sum_{(x,y)\in\Lset}y +\hat m\rb{|\Lset|-|\Rset|}}. \end{align*}

If \(|\Lset|=|\Rset|\), the contributions involving \(\hat m\) cancel.

For affine functions as predictors the best-in-class predictor is not informative. However, the median regression

\begin{align*} \begin{split} (\hat\beta_0,\hat\beta) &=\argmin_{\beta_0,\beta}\hat R_{\Tset}(\beta(\cdot)),\\ \hat R_{\Tset}(\beta(\cdot)) &=\frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}\abs{y-\beta(x)} \end{split} \tag{3.2.10} \end{align*}

can be solved as an interesting linear program:

\begin{align*} \min_{\beta_0,\beta,t}\quad &\frac{1}{|\Tset|}\sum_{i=1}^{|\Tset|}t_i\\ \text{subject to}\quad &-t_i\le y_i-\beta(x_i)\le t_i,& i=1,\ldots,|\Tset|, \\ &t_i\ge 0,& i=1,\ldots,|\Tset|. \end{align*}

At an optimum, \(t_i\) equals the absolute residual, so the objective is precisely (3.2.10).

Note that restricting to \(\beta=0\) recovers the marginal median at population level.

Zero-one loss compares samples based on class. Using (3.3.2),

\begin{align*} \lscr(r(y),f(x)) &= \1{y \neq d(f(x))} = \1{y \neq \argmax_{k\in\Kset}f_k(x)}, & \lscr \in \set{0,1}. \tag{3.3.3} \end{align*}

The next numerical example shows an important problem with zero-one loss. Take \(\Kset=\set{0,1}\) and an observed class \(y=1\), so that \(r(y)=r(1) = (0,1)\). Compare the probability predictions where the first is hesitant, the other confident:

\begin{align*} p^{A}&=(0.49,0.51), & p^{B}&=(0.01,0.99). \end{align*}

Both have maximizing class \(1\), since \(p_1^{A}>p_0^{A}\) and \(p_1^{B}>p_0^{B}\). Using (3.3.3),

\begin{align*} \lscr(r(1),p^{A}) &=\1{1\neq\argmax_{k\in\Kset}p_k^{A}}=\1{1\neq1}=0,\\ \lscr(r(1),p^{B}) &=\1{1\neq\argmax_{k\in\Kset}p_k^{B}}=\1{1\neq1}=0. \end{align*}

Thus zero-one loss gives the same value to a hesitant prediction and a confident prediction. Optimizing under this loss is difficult because small changes in the probabilities often leave the predicted label unchanged.

For any rule \(f:\Xset\to\Delta_K\), using (3.3.3) gives the conditional risk

\begin{align*} \E{\lscr(r(Y),f(X))\mid X} &= \E{\1{Y \neq \argmax_{k\in\Kset}f_k(X)}\mid X} \\ &= \P{Y \neq \argmax_{k\in\Kset}f_k(X)\mid X} \\ &= 1 - \P{Y = \argmax_{k\in\Kset}f_k(X)\mid X}. \end{align*}

Now take \(f(x) = p(x)\) with components

\begin{align*} f^{*}_{k}(x) = \P{Y=k \mid X=x} = p_{k}(x). \tag{3.3.4} \end{align*}

Then the expected loss is minimized when the predicted class is the most probable because

\begin{align*} \P{Y = \argmax_{k\in\Kset}f_k(X)\mid X} = \argmax_{k\in\Kset}\P{Y=k\mid X}. \end{align*}

Thus, \(f\) can serve as a Bayes predictor.

A constant predictor maps, via the decision map \(d\), any \(x\in \Xset\) to some constant \(k\in \Kset\). With this, the risk is easy:

\begin{align*} R(k) = \E{Y\neq k} = 1-\P{Y=k}. \end{align*}

Thus the best constant predictor is the majority class and the risk is the misclassification rate,

\begin{align*} c^{*} &= \argmax_{k} |\Tset_{k}|, & \hat R_{\Tset}(c^{*}) &= 1 - \max_{k} \hat p_{k}. \tag{3.3.5} \end{align*}

The Brier loss³²³²Wikipedia: Brier score.

\begin{align*} \lscr(r(y),f(x)) &= \sum_{k\in\Kset}(r_k(y)-f_k(x))^2, & \lscr \in [0,2], \tag{3.3.6} \end{align*}

compares the one-hot representation \(r(y)\) to the prediction \(f(x)\).

A numerical example shows the advantage of Brier loss over zero-one loss. For the example in zero loss,³³³³\(\Kset=\set{0,1}\), observed class \(y=1\), \(p^{A} = (0.49,0.51)\), \(p^{B} = (0.01,0.99)\). the Brier loss becomess \(p^{A}\),

\begin{align*} \lscr(r(1),p^{A}) &= (r_0(1)-p^{A}_{0})^{2} + (r_1(1)-p^{A}_{1})^{2} \\ &= (0-0.49)^2 + (1-0.51)^2 \\ &\approx 0.48. \\ \lscr(r(1),p^{B}) &= (0-p^{B}_{0})^{2} + (1-p^{B}_{1})^{2} \\ &= (0-0.01)^2 + (1-0.99)^2 \\ &= 0.0002. \end{align*}

Thus the second prediction receives a much smaller Brier loss.

For the Bayes predictor we split the risk in a similar way as Th. 3.2.1. As the \(k\)th component of \(r(Y)\) is \(\1{Y=k}\), the Brier loss becomes

\begin{align*} \E{\lscr(r(Y), f(X))} = \E{\sum_{k\in\Kset}\rb{\1{Y=k} - f_{k}(X)}^{2}} = \sum_{k\in\Kset} \E{\rb{\1{Y=k} - f_{k}(X)}^{2}}. \end{align*}

To minimize this, as \(f_k(X)\) is free, it suffices to minimize each terms separately. But observe that \(\E{\rb{\1{Y=k} - f_{k}(X)}^{2}}\) is a regression problem under square loss. We can therefore directly use the formulas of Th. 3.2.1 with \(\1{Y=k}\) in the role of \(Y\) and \(f_k(X)\) in the role of \(f(X)\). Thus, ³⁴³⁴Recalll \(\E{\1{X}} = p\) \(\implies\) \(\V{\1{X}} = p(1-p)\).

\begin{align*} \E{\1{Y=k}\mid X} &= p_{k}(X), & \V{\1{Y=k}\given X} &= p_{k}(X)\rb{1-p_{k}(X)}. \end{align*}

Summing the variances over \(k\) and using that the \(p_{k}(x)\) add up to one,

\begin{align*} \sum_{k\in\Kset} p_{k}(x)\rb{1-p_{k}(x)} = 1 - \sum_{k\in\Kset} p_{k}(x)^{2}. \end{align*}

Therefore, immediately from Th. 3.2.1:

From the above and Cor. 3.2.2 we obtain the following corollary.

A constant vector in \(\Delta_{K}\) has constant components \(c_k\). We can therefore use (3.2.2) with \(\1{Y=k}\) in the role of \(Y\). The optimal constant is \(\E{\1{Y=k}} = p_k\) and the variance is \(\V{\1{Y=k}} = p_k(1-p_k)\). Thus the best-in-class predictor and its risk are

\begin{align*} c^{*} &=p, & c_{k} &= \E{\1{Y=k}}, & R(c^{*})&=1-\sum_{k\in\Kset}p_k^2. \end{align*}

Replacing the marginal probabilities by sample frequencies gives

\begin{align*} \hat c&=\hat p, & c_{k} &=\frac{|\Tset_{k}|}{|\Tset|}, & \hat R_{\Tset}(\hat c) &= 1 - \sum_{k\in\Kset} \hat p_{k}^{2}. \tag{3.3.7} \end{align*}

Affine scores as predictors need a transformation to give valid probability predictions, since an unconstrained affine rule does stay in the simplex \(\Delta_K\). The sigmoid and softmax transformations form one way to achieve the required mapping:

\begin{align*} \sm: \R &\to (0,1), & \sm(z) &= \frac{1}{1+e^{-z}}, \tag{3.3.8} \\ \sm: \R^{k} &\to (0,1)^{k}, & \sm_{k}(s) &= \frac{e^{s_{k}}}{\sum_{j\in\Kset} e^{s_{j}}}, & k&=2, 3, 4, \ldots. \tag{3.3.9} \end{align*}

Observe that for for binary classification, i.e., \(k=2\),

\begin{align*} \sigma_{0}(a,b) &= \frac{e^{a}}{e^{a}+e^{b}} = \frac{1}{1 + e^{b-a}}. \end{align*}

Thus, for \(k=2\), we can write the softmax in terms of the sigmoid function,

\begin{align*} \sigma_{0}(a,b) = \sigma(b-a) = 1- \sigma(a-b) = 1-\sigma_{1}(a,b). \end{align*}

Taking one intercept \(\beta_{0k}\) and one weight vector \(\beta_{k}\) per class, and writing \(\beta_{k}(x) = \beta_{0k} + \ip{\beta_{k}, x}\) for the affine score of class \(k\), write \(\beta(x)=(\beta_k(x))_{k\in\Kset}\) and let \(\theta\) collect these coefficients. The model class \(\Fset_{\mathrm{softmax}}=\{f_\theta\}\) consists of the rules

\begin{align*} f_{\theta}(x) &= \sm(\beta(x)) \in \Delta_K, & f_{\theta,k}(x) &= \sm_{k}(\beta(x)) = \PP{\theta}{Y=k\mid X=x}, \tag{3.3.10} \end{align*}

where the rhs is the probabilityunder the model with coefficients \(\theta\).

This model has \(K(p+1)\) parameters: for each of the \(K\) classes one intercept \(\beta_{0}\) and \(p\) weights for the \(\beta_{k}\). Not all of them can be recovered from data, because an extra factor \(e^{c_{0} + \ip{c, x}}\) cancels between the numerator and the denominator of (3.3.10). Only \((K-1)(p+1)\) parameters are therefore necessary. For \(\Kset=\set{0,1}\), we write \(\beta(x)=\beta_1(x)-\beta_0(x)\) and obtain

\begin{align*} f_{\theta,0}(x)=1-\sigma(\beta(x)),\qquad f_{\theta,1}(x)=\sigma(\beta(x)). \end{align*}

Under Brier loss, the training risk is

\begin{align*} \hat R_{\Tset}(f_\theta) =\frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}\sum_{k\in\Kset} \rb{r_k(y)-f_{\theta,k}(x)}^2. \end{align*}

Fitting means minimizing this expression over the coefficients \(\theta\). This cannot be solved algebra, we discuss the numerical optimization in Chapter 5.

The cross-entropy loss³⁵³⁵Wikipedia: Cross-entropy. arises from the conditional likelihood of the labels. Assuming the labels are independent conditional on their feature vectors, and taking the negative logarithm, the average log-likelihood under a candidate rule \(f\) gives

\begin{align*} \frac{1}{|\Tset|}\log \prod_{(x,y)\in\Tset}\prod_{k\in\Kset} f_k(x)^{r_{k}(y)} = \frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}\sum_{k\in\Kset} r_{k}(y)\log f_k(x) = - \frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}\lscr(r(y), f(x)), \end{align*}

where the cross-entropy loss is defined as

\begin{align*} \lscr(r(y),f(x)) &= -\sum_{k\in\Kset}r_k(y)\log f_k(x), & \lscr \in [0,\infty]. \tag{3.3.11} \end{align*}

A numerical example shows that cross-entropy loss³⁶³⁶Similar to Brier loss. assigns different losses to predictions with the same most likely class. For the example of the Brier loss, here

\begin{align*} \lscr(r(1),p^{A}) &= -\rb{0\cdot\log p^{A}_{0} + 1\cdot\log p^{A}_{1}} \approx 0.67, \\ \lscr(r(1),p^{B}) &= -\rb{0\cdot\log p^{B}_{0} + 1\cdot\log p^{B}_{1}} \approx 0.01. \end{align*}

Observe also that as \(0\log 0=0\), a class with zero observed weight, i.e., \(r_k(y)=0\), contributes zero, and by taking \(\log 0 = -\infty\), assigning probability zero to the observed class, i.e., \(r_k(y)=1\), gives infinite loss.

For the Bayes predictor we minimize population risk over all prediction functions.

Proof

Condition on \(X\). Then, as the \(k\)th component of \(r(Y)\) is \(\1{Y=k}\),

\begin{align*} \E{\lscr(r(Y), f(X)) \given X} &= \sum_{k\in\Kset} \E{\1{Y=k}\log\frac{1}{f_k(X)} \given X} \\ &= \sum_{k\in\Kset}p_k(X)\log\frac{1}{f_k(X)}. \end{align*}

Adding and subtracting \(\sum_{k}p_k(X)\log p_k(X)\) gives

\begin{align*} \sum_{k\in\Kset}p_k(X)\log\frac{p_k(X)}{f_k(X)} -\sum_{k\in\Kset}p_k(X)\log p_k(X). \end{align*}

Taking the expectation establishes the claim.

The first term of this decomposition of the risk is the relative entropy³⁷³⁷Wikipedia: Kullback-Leibler divergence. between \(p(X)\) and \(f(X)\), and the second is the conditional entropy³⁸³⁸Wikipedia: Conditional entropy. of \(Y\) given \(X\). The first is non-negative by Gibbs’ inequality³⁹³⁹Wikipedia: Gibbs’ inequality. and zero exactly when \(f(X)=p(X)\)

A constant vector as predictor has constant components \(c_k\). Substituting \(f_k(X)=c_k\) into the loss gives

\begin{align*} R(c) &= -\E{\sum_{k\in\Kset}r_k(Y)\log c_k} = -\sum_{k\in\Kset}\E{r_k(Y)}\log c_k &= -\sum_{k\in\Kset}p_k\log c_k \\ &= \sum_{k\in\Kset}p_k\log\frac{p_k}{c_k} -\sum_{k\in\Kset}p_k\log p_k, \end{align*}

since \(\E{r_k(Y)}=\P{Y=k}=p_k\). The first term in the last line is the relative entropy between \(p\) and \(c\), so by Gibbs’ inequality it is non-negative and vanishes precisely at \(c=p\). The optimal constant predictors are the same as in the Brier analysis: \(c^*=p\) and \(\hat c=\hat p\), with the empirical predictor given by (3.3.7). The corresponding risks are now the population and empirical entropies:

\begin{align*} R(c^*)&=-\sum_{k\in\Kset}p_k\log p_k, & \hat R_{\Tset}(\hat c)&=-\sum_{k\in\Kset}\hat p_k\log\hat p_k. \end{align*}

Affine scores as predictors use the same softmax model (3.3.10) as in the Brier analysis. Fitting this model with cross-entropy gives logistic regression. Fitting the coefficients requires numerical optimization, which we discuss in Chapter 5.

The absolute loss

\begin{align*} \lscr(r(y),f(x)) &= \sum_{k\in\Kset}|r_k(y)-f_k(x)|, & \lscr \in [0,2], \end{align*}

does not generally give the true class probabilities as its Bayes rule. Take \(\Kset=\set{0,1}\) and let \(f(x)=(1-u,u)\), with \(u\in[0,1]\). Here \(u\) is the predicted probability of class \(1\), while \(p_1(x)\) is its true conditional probability. Since \(r(0)=(1,0)\) and \(r(1)=(0,1)\), the two possible losses are

\begin{align*} \lscr(r(0),f(x)) &=|1-(1-u)|+|0-u|=2u,\\ \lscr(r(1),f(x)) &=|0-(1-u)|+|1-u|=2(1-u). \end{align*}

Weighting these losses by the true class probabilities gives

\begin{align*} \E{\lscr(r(Y),f(X))\mid X} &=2p_0(X)u+2p_1(X)(1-u)\\ &=2p_1(X)+2u\rb{1-2p_1(X)}. \end{align*}

This is linear in \(u\). If \(p_1(x)>1/2\), its coefficient of \(u\) is negative, so increasing \(u\) lowers the risk all the way to \(u=1\). If \(p_1(x)<1/2\), the coefficient is positive, so the minimum is at \(u=0\). If \(p_1(x)=1/2\), every \(u\in[0,1]\) gives the same risk.

Absolute loss can be used to learn class predictions, but its risk minimizer does not generally recover how plausible each class is. We prefer Brier loss and cross-entropy when we also want to learn these class probabilities.

If several training samples have exactly the feature vector \(x\), we can estimate \(m(x)\) by averaging their outcomes. For continuously distributed features, exact matches to a specified \(x\) generally do not occur. Instead, we can average outcomes at feature vectors near \(x\), using the idea that their conditional means should also be close.

Choose a number of neighbors \(k\), with \(1\leq k\leq n\), and measure closeness by Euclidean distance:

\begin{align*} \norm{x_i-x}_2 =\left(\sum_{j=1}^p (x_i(j)-x(j))^2\right)^{1/2}. \end{align*}

Let \(i_1(x),\ldots,i_n(x)\) be the training indices ordered by increasing distance to \(x\). Break equal-distance ties by the original sample index, so that the first \(k\) indices always identify exactly \(k\) samples, including repeated samples when present. The \(k\)-nearest-neighbor (KNN)⁴⁰⁴⁰We use lowercase \(k\) for the number of neighbors; uppercase \(K\) continues to denote the number of classes. prediction is

\begin{align*} \hat f_{k,\Tset}(x)=\frac{1}{k}\sum_{a=1}^k y_{i_a(x)}. \tag{3.4.1} \end{align*}

This is precisely the fitted constant of (3.2.3), but fitted to the local list \(\Tset_k(x)=[(x_{i_1(x)},y_{i_1(x)}),\ldots,(x_{i_k(x)},y_{i_k(x)})]\). The list changes with the prediction point. The resulting rule can therefore follow a nonlinear conditional mean even though each individual prediction is obtained by fitting a constant.

\begin{algorithm}
\caption{Predict by averaging the outcomes of exactly \(k\) nearest training samples.}
\begin{algorithmic}[1]
\State \textbf{Input:} a training list $\Tset$, a feature vector $x$, and an integer $1\leq k\leq|\Tset|$.
\State \textbf{Output:} a numerical prediction at $x$.
\Procedure{KnnPredict}{$\Tset,x,k$}
    \State Order the indices $i_1,\ldots,i_n$ by increasing $\norm{x_i-x}_2$, breaking ties by index.
    \Return $\pr{Average}([y_{i_1},\ldots,y_{i_k}])$
\EndProcedure
\end{algorithmic}

\end{algorithm}

Distances depend on the units of the features. A coordinate measured in euros can dominate one measured in kilometers even if both are useful. One common choice is to center each feature and divide it by its training standard deviation before calculating distances. The same transformation, with the same fitted means and standard deviations, is then applied to new inputs. For \(1\leq j\leq p\), let⁴¹⁴¹We apply the last transformation only when \(\sigma(j)>0\); constant features can be omitted.

\begin{align*} \mu(j)&=\frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}x(j),& \sigma(j)&=\sqrt{\frac{1}{|\Tset|}\sum_{(x,y)\in\Tset}(x(j)-\mu(j))^2},& \tilde x_i(j)&=\frac{x_i(j)-\mu(j)}{\sigma(j)}. \tag{3.4.2} \end{align*}

A numerical example demonstrates clearly the effect of \(k\). Consider the conditional mean

\begin{align*} m(x)=2-\frac{1}{2}x^2+3x+\sin(4x),\qquad -1\leq x\leq2. \tag{3.4.3} \end{align*}

Generate \(n=200\) independent samples according to

\begin{align*} X_i\sim\operatorname{Unif}(-1,2),\qquad Y_i=m(X_i)+\varepsilon_i,\qquad \varepsilon_i\sim\mathcal N(0,1.5^2), \end{align*}

with the noises independent of the inputs and of one another. Then \(\E{Y_i\given X_i=x}=m(x)\).⁴²⁴²Because this is a constructed example, we know \(m\); in an application we only observe the pairs \((x_i,y_i)\).

Fig. 3.5 compares \(k=200,100,50,10,5,1\) on the same sample. At \(k=200\) the rule is constant and cannot follow the trend or the oscillations. At \(k=1\), the fitted rule reproduces each training outcome, including its noise, which is clearly overfitting.

When the number of features \(p\) is large, e.g. \(p\geq 10\), the nearest neighbors of \(x\) are no longer close to \(x\). Two computations show this for samples that are uniformly distributed on \([0,1]^p\).

First, consider a sub-cube with sidelength⁴³⁴³Hence with volume \(a^p\). \(a\) contains a fraction \(s\) of the samples. Then, \(a = s^{1/p}\). For \(s=0.01\) and \(p=10\) this side length is \(0.01^{1/10}\approx 0.63\): a neighborhood that holds 1% of the samples spans 63% of the range of each feature.

Second, let \(X\) and \(X'\) be independent uniform points in \([0,1]^p\). Then

\begin{align*} \E{\norm{X-X'}_2^2} &= \sum_{j=1}^{p}\E{(X(j)-X'(j))^2} = \frac{p}{6}, \tag{3.4.4} \\ \V{\norm{X-X'}_2^2} &= \sum_{j=1}^{p}\V{(X(j)-X'(j))^2} = \frac{7p}{180}, \tag{3.4.5} \end{align*}

where the second line uses that the coordinate differences are independent. The squared coefficient of variation⁴⁴⁴⁴The variance divided by the squared mean. is

\begin{align*} \frac{\V{\norm{X-X'}_2^2}}{\E{\norm{X-X'}_2^2}^{\,2}} = \frac{7}{5p}. \tag{3.4.6} \end{align*}

Thus, as \(p\) grows, the variation of squared distances relative to their mean becomes small: all points lie roughly equally far away. The \(k\) nearest samples are then barely closer to \(x\) than the other samples, so their average says little more about \(m(x)\) than an average over arbitrary samples. Yet using all \(n\) observations is also undesirable, since it gives the global constant predictor.

This is the curse of dimensionality: local averaging becomes difficult in high dimension. An alternative restricts the rule to a class \(\Fset\) with a fixed number of parameters. All \(n\) samples then contribute to fitting each parameter, not only the \(k\) samples nearest to \(x\).

Polynomial regression is one such class. For a nonnegative integer degree \(q\), define

\begin{align*} \Fset_q =\left\{f_\beta:f_\beta(x)=\beta_0+\sum_{j=1}^q\beta_jx^j,\ \beta\in\R^{q+1}\right\}. \tag{3.4.7} \end{align*}

The cases \(q=0\) and \(q=1\) recover the constant and affine classes for one numerical input. For \(q>1\), the rule is generally nonlinear in \(x\), but remains affine in the transformed features \((x,x^2,\ldots,x^q)\) and linear in the coefficients. Fitting therefore uses the same least-squares method as before:

\begin{align*} \hat\beta=\argmin_{\beta\in\R^{q+1}} \frac{1}{|\Tset|}\sum_{(x,y)\in\Tset} \left(y-\beta_0-\sum_{j=1}^q\beta_jx^j\right)^2. \end{align*}

The coefficients \(\beta\) are parameters: least squares fits them for a given degree. The degree \(q\) is a hyperparameter: it fixes the class \(\Fset_q\) in which we fit.

Fig. 3.6 compares degrees \(q=1,4,10,30\) on the same observations used for KNN. An affine rule, \(q=1\), cannot follow the oscillations in \(m\). A polynomial of degree \(q\) has at most \(q-1\) turning points, so a higher degree can improve the approximation of \(m\), but it also lets the fit respond more strongly to the noise in the sample. The \(q=10\) fit follows \(m\) closely; the \(q=30\) fit adds oscillations, especially near the endpoints.

Choosing and assessing a prediction rule takes three steps. First, we select a function class \(\Fset\) over which we fit. A function in \(\Fset\) is specified by hyperparameters and parameters. For instance, for the union of the polynomial classes (3.4.7), the degree \(q\) is a hyperparameter and the coefficients \(\beta_0,\ldots,\beta_q\) are the parameters. Here we discuss how to tune the hyperparameters and the coeficients and assess how well the final result would do on a new observation \(X\). To tune a hyperparameter, we fix in advance a finite set of values to compare, for instance the degrees \(\Gamma = \set{0,1,\ldots,30}\). Each value in this set is a candidate, and the rule fitted at a candidate is a candidate rule.

For independent observations from a common population, we can randomly divide the available sample into three disjoint lists:

• The training list \(\Tset\): used to fit candidate rules.
• The validation list \(\Vset\): used to select a candidate.
• The test list \(\Eset\): kept aside to assess the final rule after selection and fitting are complete.

Splitting the samples into these subsets affects both fitting accuracy and evaluation precision: more training data can improve a fitted rule, while more evaluation data make its measured risk more precise. There is no single split proportion that is appropriate for every sample size and prediction problem.

Let \(\Gamma\) be a finite collection of candidate settings and let \(\hat f_{\gamma,\Tset}\) be the rule fitted with setting \(\gamma\). Now we use the validation set \(\Vset\) to find the best \(\gamma\). Writing the risk on the validation set as \(\hat R_\Vset(\hat f_{\gamma, \Tset})\),⁴⁸⁴⁸Note: each candidate is fitted without using the validation observations. the best candidate is

\begin{align*} \hat\gamma=\argmin_{\gamma\in\Gamma} \hat R_{\Vset}(\hat f_{\gamma,\Tset}). \tag{3.6.1} \end{align*}

After choosing \(\hat\gamma\), we can refit on the concatenated list \( \Dset=\Tset+\Vset. \) The final rule is then \(\hat f_{\hat\gamma,\Dset}\).

Finally, we use the test list to assess the final rule. Its test risk

\begin{align*} \hat R_{\Eset}(\hat f_{\hat\gamma,\Dset}) \end{align*}

is the average loss of \(\hat f_{\hat\gamma,\Dset}\) over the test samples. The test samples are used neither to fit the parameters nor to select \(\hat\gamma\). For the final rule they are therefore new observations, and the test risk estimates its risk \(R(\hat f_{\hat\gamma,\Dset})\) on a new observation \((X,Y)\). Since the test risk is an average of losses on independent samples, a larger test list gives a more precise estimate.

Preprocessing is part of fitting. Means and standard deviations used to scale inputs must be estimated on the training part, then applied unchanged to the validation part. Feature selection, missing-value imputation, and other transformations estimated from data must follow the same procedure. Allowing the validation or test data to influence fitting is called data leakage.

In the procedure above the candidate rules are fitted on \(\Tset\) only, and \(\Vset\) is never used for this purpose. Cross-validation divides \(\Dset=\Tset+\Vset\) instead into several parts. Each part serves once as validation list, while the rules are fitted on the other parts. Every sample in \(\Dset\) is then used for fitting and for validation, but never both in the same fit.

Let \(N=|\Dset|\), and divide its sample indices into \(F\) disjoint non-empty folds \(\Iset_1,\ldots,\Iset_F\), usually of nearly equal size. For fold \(b\), write

\begin{align*} \Dset_b &=[(x_i,y_i):i\in\Iset_b],& \Dset_{-b} &= \Dset\setminus \Dset_{b} = [(x_i,y_i):i\notin\Iset_b]. \end{align*}

In \(F\)-fold cross-validation, we fit a candiate \(\gamma\) on \(\Dset_{-b}\) and we evaluate its risk on \(\Dset_b\).⁴⁹⁴⁹Use the same folds for every candidate, so their errors are compared on the same validation observations. Writing \(\hat f_{\gamma,\Dset_{-b}}\) for such a fit, the cross-validation criterion for candidate \(\gamma\) is

\begin{align*} \mathrm{CV}(\gamma) &=\frac{1}{|\Dset|}\sum_{b=1}^F\sum_{(x,y)\in\Dset_b} \lscr(r(y),\hat f_{\gamma,\Dset_{-b}}(x)) \\ &=\sum_{b=1}^F\frac{|\Dset_b|}{|\Dset|} \hat R_{\Dset_b}(\hat f_{\gamma,\Dset_{-b}}). \tag{3.6.2} \end{align*}

Every sample contributes one loss, from a rule fitted without that sample’s fold. The weights account for unequal fold sizes. Now select

\begin{align*} \hat\gamma=\argmin_{\gamma\in\Gamma}\mathrm{CV}(\gamma). \end{align*}

Again, once \(\hat \gamma\) is computed, fit the model on \(\Dset\) for this candidate.

A special case is when \(F=N\); this is called leave-one-out cross-validation: each fit omits just one sample.

\begin{algorithm}
\caption{Select a setting by cross-validation and refit on all development data.}
\begin{algorithmic}[1]
\State \textbf{Input:} a development list $\Dset$, candidate settings $\Gamma$, and $2\leq F\leq|\Dset|$.
\State \textbf{Output:} a selected setting and the final fitted rule.
\Procedure{CrossValidate}{$\Dset,\Gamma,F$}
    \State Divide the indices of $\Dset$ into folds $\Iset_1,\ldots,\Iset_F$.
    \For{$\gamma\in\Gamma$}
        \State $\mathrm{CV}(\gamma)\gets0$
        \For{$b=1,\ldots,F$}
            \State Fit all preprocessing and the rule at $\gamma$ on $\Dset_{-b}$.
            \State $\mathrm{CV}(\gamma)\gets\mathrm{CV}(\gamma)+
            \dfrac{|\Dset_b|}{|\Dset|}\cdot\hat R_{\Dset_b}(\hat f_{\gamma,\Dset_{-b}})$
        \EndFor
    \EndFor
    \State $\hat\gamma\gets$ a minimizer of $\mathrm{CV}(\gamma)$ over $\Gamma$
    \State Refit preprocessing and the rule at $\hat\gamma$ on all of $\Dset$.
    \Return $(\hat\gamma,\hat f_{\hat\gamma,\Dset})$
\EndProcedure
\end{algorithmic}

\end{algorithm}

The number of folds \(F\) affects cross-validation in three ways.

First, with equal folds, each fold rule is fitted on \(N(1-1/F)\) samples, while the final rule \(\hat f_{\hat\gamma,\Dset}\) is fitted on all \(N\). A rule fitted on fewer samples has a larger estimation error, so with cross-validation we obtain usually less good predictors. Hence \(\mathrm{CV}(\gamma)\) tends to be somewhat larger than the risk of the rule fitted at \(\gamma\) on all of \(\Dset\).

Second, cross-validation fits every candidate \(F\) times, so \(F|\Gamma|\) fits in total, plus the final refit. Hence, CV increases the compuational load.

Third, the random variation of \(\mathrm{CV}(\gamma)\). Whatever \(F\), each sample of \(\Dset\) is validated exactly once, so \(\mathrm{CV}(\gamma)\) is always an average of \(N\) losses. What changes with \(F\) is how much the fold rules have in common. For large \(F\), any two fitting lists \(\Dset_{-b}\) share all but \(2N/F\) samples, so the fold rules are nearly the same rule. If that rule happens to be poor, all its losses are large together, and averaging them removes less of this chance variation than averaging over rules fitted on different samples. Taking \(F=5\) or \(F=10\) keeps the fold training size at 80 or 90 percent of \(N\), with 5 or 10 fits per candidate.