Why Machine Learning Needs Mathematics: Rewriting Learning as Loss Minimisation

0. Key points

• Supervised learning, once the data, the hypothesis class and the loss function have been fixed, can be rewritten as a single optimisation problem: find the parameters that minimise the empirical risk. Regression and classification both fit inside this frame.
• Arranging the data into a design matrix $X$ and a target vector $\boldsymbol{y}$ turns the minimisation of the squared loss into the normal equation $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$, and geometrically into the orthogonal projection of $\boldsymbol{y}$ onto the column space (Theorem 3.2). This is the role of linear algebra.
• The tool that carries out the minimisation is differentiation. For a convex function, “the gradient vanishes” and “the point is a global minimiser” are equivalent (Theorem 4.2), and both the range of learning rates for which gradient descent converges and the optimal learning rate are determined solely by the eigenvalues of $\frac{1}{n}X^{\mathsf{T}}X$ (Theorem 4.4).
• Probability answers the question “why the squared loss?”. Maximum likelihood estimation under Gaussian noise coincides exactly with least squares (Theorem 5.2), while assuming a Bernoulli distribution produces the cross-entropy (Example 5.3).
• Probability also explains why we may use the error measured on our finitely many samples in place of the error on unseen data (Proposition 5.4). This guarantee, however, breaks the moment the model is chosen after looking at the data. That is what overfitting really is.

1. Motivation: what a single line of a learning library hides

With a machine learning library, training a model takes one line. Hand over the data, call the fitting method, and a few seconds later you can make predictions. Seen from here, mathematics appears to have no part to play.

The trouble starts when that one line does not work. For instance, the following things happen.

1. Changing the learning rate from $0.001$ to $0.01$ makes the loss diverge instead of decrease.
2. Changing the units of a feature from centimetres to metres cuts the number of iterations needed for convergence by a factor of several dozen.
3. The error on the training data is almost zero, yet the model is hopeless on new data.
4. Regression uses the squared error, while classification uses a different quantity called the cross-entropy.

Trial and error against these problems has no end, but rewritten in mathematical language each of them is explained in a single short sentence. In the second half of this article we answer all four.

Looking back at the history, the problem of “fitting a model to data” is some two hundred years older than machine learning. Legendre published the method of least squares in 1805, as an appendix to a work on determining the orbits of comets. In his Theoria Motus Corporum Coelestium of 1809, Gauss showed that if observation errors are assumed to follow a normal distribution, then least squares yields the most likely estimate. The idea of deriving the loss function from a probabilistic model, still used verbatim in modern deep learning, was already present at that point.

This article has two goals. First, to formulate regression and classification, two apparently separate tasks, as one and the same optimisation problem, the minimisation of the empirical risk. Second, to see where and why linear algebra, calculus and probability each become necessary in solving that problem, by carrying a single example, linear regression, all the way through the computation.

2. Formulating the learning problem

We first fix the vocabulary. From now on vectors are written in bold, $\boldsymbol{x}$, and matrices as $X$; the standard inner product on $\mathbb{R}^n$ is $\langle\boldsymbol{a},\boldsymbol{b}\rangle=\sum_{i}a_ib_i$ and the associated norm is $\|\boldsymbol{a}\|=\sqrt{\langle\boldsymbol{a},\boldsymbol{a}\rangle}$.

Definition 2.1（The supervised learning setup）

Let $\mathcal{X}$ be an input space and $\mathcal{Y}$ an output space, both sets, and fix a probability distribution $P$ on $\mathcal{X}\times\mathcal{Y}$. The $n$ samples

$$D = \bigl((\boldsymbol{x}_1,y_1),\ldots,(\boldsymbol{x}_n,y_n)\bigr) \in (\mathcal{X}\times\mathcal{Y})^n$$

drawn independently from $P$, all with the same distribution, are called the training data. Fix a space $\mathcal{Y}'$ of predicted values; a map $f:\mathcal{X}\to\mathcal{Y}'$ is called a predictor, and a set $\mathcal{H}$ of predictors fixed in advance is called a hypothesis class. Finally, a function $\ell:\mathcal{Y}'\times\mathcal{Y}\to[0,\infty)$ is called a loss function, and $\ell(\hat{y},y)$ expresses the penalty for predicting $\hat{y}$ when the correct answer is $y$.

The case $\mathcal{Y}=\mathbb{R}$ is called regression, and the case where $\mathcal{Y}$ is a finite set is called classification.

The reason for introducing a space $\mathcal{Y}'$ of predicted values separate from $\mathcal{Y}$ is classification. Even in binary classification, where the labels are $\mathcal{Y}=\{0,1\}$, what the model outputs is normally the probability that the label is $1$, that is, a value in $\mathcal{Y}'=[0,1]$. This distinction earns its keep in Example 5.3.

Definition 2.2（Expected risk and empirical risk）

In the setting of Definition 2.1, for a predictor $f\in\mathcal{H}$ put

$$R(f) := \mathbb{E}_{(\boldsymbol{X},Y)\sim P}\bigl[\ell(f(\boldsymbol{X}),Y)\bigr], \qquad \hat{R}_n(f) := \frac{1}{n}\sum_{i=1}^{n}\ell\bigl(f(\boldsymbol{x}_i),y_i\bigr) .$$

These are called the expected risk (generalisation error) and the empirical risk (training error) of $f$. The learning principle of selecting the predictor that minimises $\hat{R}_n$ over $\mathcal{H}$ is called empirical risk minimisation.

The relation between these two quantities generates very nearly all of the difficulty in machine learning. What we truly want to make small is $R(f)$, but the distribution $P$ is unknown, so $R(f)$ cannot be computed. All we can compute is $\hat{R}_n(f)$, measured on the data at hand. How far this substitution is justified is treated in Proposition 5.4.

Once the hypothesis class is described by finitely many parameters, learning becomes a finite-dimensional optimisation problem. Parametrising $\mathcal{H}=\{f_{\boldsymbol{w}}\mid \boldsymbol{w}\in\mathbb{R}^d\}$ and setting

$$L(\boldsymbol{w}) := \hat{R}_n(f_{\boldsymbol{w}}),$$

learning is nothing other than the problem “find a minimiser of the function $L$ on $\mathbb{R}^d$”. The overall flow is as follows.

The basic pipeline of supervised learning, with the mathematics used at each stage

We now look at stages 2, 5 and 6 in turn. The subject matter is the simplest possible one, linear regression, and yet all three areas of mathematics show their faces there.

3. Linear algebra: arranging the data into a matrix

Take the inputs to be points of $\mathcal{X}=\mathbb{R}^d$, that is, tuples of $d$ real-valued features. Stacking the $n$ inputs $\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n$ vertically gives the matrix

$$X = \begin{pmatrix} \boldsymbol{x}_1^{\mathsf{T}} \\ \vdots \\ \boldsymbol{x}_n^{\mathsf{T}} \end{pmatrix} \in \mathbb{R}^{n\times d}, \qquad \boldsymbol{y} = \begin{pmatrix} y_1 \\ \vdots \\ y_n\end{pmatrix}\in\mathbb{R}^n$$

called the design matrix and the target vector. Rows correspond to data points and columns to features. When an intercept (bias) is needed, taking the first feature to be the constant $x_{i1}=1$ absorbs it into $\boldsymbol{w}$, so below we give the intercept no special treatment.

Taking the linear model $f_{\boldsymbol{w}}(\boldsymbol{x})=\langle\boldsymbol{w},\boldsymbol{x}\rangle$, all $n$ predictions are written together as a single matrix-vector product $X\boldsymbol{w}$. This is the first benefit of linear algebra as a language of representation. A loop over individual data points is replaced by one matrix operation, which not only shortens the description but also lets an implementation hand the work over to highly optimised matrix-product routines or to a GPU.

With the loss $\ell(\hat{y},y)=\tfrac{1}{2}(\hat{y}-y)^2$ (the factor $\tfrac12$ is merely a convenience so that the $2$ disappears on differentiation), the empirical risk becomes

$$L(\boldsymbol{w}) = \frac{1}{n}\sum_{i=1}^{n}\frac{1}{2}\bigl(\langle\boldsymbol{w},\boldsymbol{x}_i\rangle-y_i\bigr)^2 = \frac{1}{2n}\|X\boldsymbol{w}-\boldsymbol{y}\|^2 .$$

Finding a minimiser of this function is linear regression. The key is the following completely elementary identity.

Lemma 3.1（Quadratic expansion of the squared loss）

Let $X\in\mathbb{R}^{n\times d}$ and $\boldsymbol{y}\in\mathbb{R}^n$, and put $L(\boldsymbol{w})=\dfrac{1}{2n}\|X\boldsymbol{w}-\boldsymbol{y}\|^2$. For all $\boldsymbol{w},\boldsymbol{u}\in\mathbb{R}^d$,

$$L(\boldsymbol{w}+\boldsymbol{u}) = L(\boldsymbol{w}) + \frac{1}{n}\bigl\langle \boldsymbol{u},\, X^{\mathsf{T}}(X\boldsymbol{w}-\boldsymbol{y})\bigr\rangle + \frac{1}{2n}\|X\boldsymbol{u}\|^2 .$$

In particular $L$ is totally differentiable on $\mathbb{R}^d$, with gradient $\nabla L(\boldsymbol{w}) = \dfrac{1}{n}X^{\mathsf{T}}(X\boldsymbol{w}-\boldsymbol{y})$.

Proof（Lemma 3.1）

Put $\boldsymbol{r}:=X\boldsymbol{w}-\boldsymbol{y}$, so that $X(\boldsymbol{w}+\boldsymbol{u})-\boldsymbol{y}=\boldsymbol{r}+X\boldsymbol{u}$. Bilinearity and symmetry of the inner product give

$$\|\boldsymbol{r}+X\boldsymbol{u}\|^2 = \langle \boldsymbol{r}+X\boldsymbol{u},\,\boldsymbol{r}+X\boldsymbol{u}\rangle = \|\boldsymbol{r}\|^2 + 2\langle X\boldsymbol{u},\boldsymbol{r}\rangle + \|X\boldsymbol{u}\|^2 .$$

Applying the defining property of the transpose, $\langle X\boldsymbol{u},\boldsymbol{r}\rangle = \langle \boldsymbol{u}, X^{\mathsf{T}}\boldsymbol{r}\rangle$, to the second term and dividing throughout by $2n$ yields the stated identity.

We check differentiability. The second term is a linear form in $\boldsymbol{u}$. The third term is bounded, by the Cauchy–Schwarz inequality（Theorem 4.2）［内積空間とグラム・シュミット直交化］, as

$$\|X\boldsymbol{u}\|^2 = \sum_{i=1}^n \langle \boldsymbol{x}_i,\boldsymbol{u}\rangle^2 \le \Bigl(\sum_{i=1}^n\|\boldsymbol{x}_i\|^2\Bigr)\|\boldsymbol{u}\|^2 ,$$

and $C:=\sum_i\|\boldsymbol{x}_i\|^2$ is a constant independent of $\boldsymbol{u}$, so that $\frac{1}{2n}\|X\boldsymbol{u}\|^2 = O(\|\boldsymbol{u}\|^2) = o(\|\boldsymbol{u}\|)$ as $\|\boldsymbol{u}\|\to 0$. Hence, by the very definition of total differentiability, $L$ is differentiable and its gradient is the coefficient vector of the linear form, namely $\frac{1}{n}X^{\mathsf{T}}(X\boldsymbol{w}-\boldsymbol{y})$.

∎

Theorem 3.2（The normal equation）

Let $X\in\mathbb{R}^{n\times d}$ and $\boldsymbol{y}\in\mathbb{R}^n$, and put $L(\boldsymbol{w})=\dfrac{1}{2n}\|X\boldsymbol{w}-\boldsymbol{y}\|^2$. A point $\boldsymbol{w}^{\star}\in\mathbb{R}^d$ is a global minimiser of $L$ if and only if

$$X^{\mathsf{T}}X\boldsymbol{w}^{\star} = X^{\mathsf{T}}\boldsymbol{y} .$$

This equation is called the normal equation.

Proof（Theorem 3.2）

Put $\boldsymbol{g}:=X^{\mathsf{T}}X\boldsymbol{w}^{\star}-X^{\mathsf{T}}\boldsymbol{y} = X^{\mathsf{T}}(X\boldsymbol{w}^{\star}-\boldsymbol{y})$. Applying Lemma 3.1 with $\boldsymbol{w}=\boldsymbol{w}^{\star}$, we obtain, for every $\boldsymbol{u}\in\mathbb{R}^d$,

$$L(\boldsymbol{w}^{\star}+\boldsymbol{u}) - L(\boldsymbol{w}^{\star}) = \frac{1}{n}\langle\boldsymbol{u},\boldsymbol{g}\rangle + \frac{1}{2n}\|X\boldsymbol{u}\|^2 .$$

(Sufficiency.) If $\boldsymbol{g}=\boldsymbol{0}$, the first term on the right vanishes and the remaining $\frac{1}{2n}\|X\boldsymbol{u}\|^2$ is nonnegative. Hence $L(\boldsymbol{w}^{\star}+\boldsymbol{u})\ge L(\boldsymbol{w}^{\star})$ for every $\boldsymbol{u}$, that is, $\boldsymbol{w}^{\star}$ is a global minimiser.

(Necessity.) Suppose $\boldsymbol{w}^{\star}$ is a global minimiser. Fix $\boldsymbol{u}$ arbitrarily and replace $\boldsymbol{u}$ by $t\boldsymbol{u}$ ($t\in\mathbb{R}$) in the identity above; minimality gives

$$0 \le \frac{t}{n}\langle\boldsymbol{u},\boldsymbol{g}\rangle + \frac{t^2}{2n}\|X\boldsymbol{u}\|^2$$

for every $t$. Multiply both sides by $n$. For $t > 0$, dividing by $t$ gives $\langle\boldsymbol{u},\boldsymbol{g}\rangle \ge -\frac{t}{2}\|X\boldsymbol{u}\|^2$, and letting $t\to 0^{+}$ yields $\langle\boldsymbol{u},\boldsymbol{g}\rangle\ge 0$. For $t < 0$, dividing by $t$ reverses the inequality, giving $\langle\boldsymbol{u},\boldsymbol{g}\rangle \le -\frac{t}{2}\|X\boldsymbol{u}\|^2$, and letting $t\to 0^{-}$ yields $\langle\boldsymbol{u},\boldsymbol{g}\rangle\le 0$. Combining the two, $\langle\boldsymbol{u},\boldsymbol{g}\rangle=0$. Since $\boldsymbol{u}$ was arbitrary, taking $\boldsymbol{u}=\boldsymbol{g}$ gives $\|\boldsymbol{g}\|^2=0$, that is, $\boldsymbol{g}=\boldsymbol{0}$.

∎

The normal equation compresses $n$ equations (the number of data points) into $d$ equations (the number of parameters). The practical content of the theorem is that however many hundreds of millions of data points there are, the size of the linear system to be solved is determined by the number of features alone. We next investigate when that system has a unique solution.

Lemma 3.3（The kernel of the Gram matrix）

For every $X\in\mathbb{R}^{n\times d}$ we have $\ker(X^{\mathsf{T}}X)=\ker X$. Consequently $\operatorname{rank}(X^{\mathsf{T}}X)=\operatorname{rank}(X)$, and the $d\times d$ matrix $X^{\mathsf{T}}X$ is invertible if and only if the $d$ column vectors of $X$ are linearly independent.

Proof（Lemma 3.3）

If $X\boldsymbol{u}=\boldsymbol{0}$, multiplying on the left by $X^{\mathsf{T}}$ gives $X^{\mathsf{T}}X\boldsymbol{u}=\boldsymbol{0}$, so $\ker X\subseteq\ker(X^{\mathsf{T}}X)$. Conversely, suppose $X^{\mathsf{T}}X\boldsymbol{u}=\boldsymbol{0}$. Taking the inner product of both sides with $\boldsymbol{u}$,

$$0 = \langle \boldsymbol{u}, X^{\mathsf{T}}X\boldsymbol{u}\rangle = \langle X\boldsymbol{u}, X\boldsymbol{u}\rangle = \|X\boldsymbol{u}\|^2 ,$$

and only the zero vector has norm $0$, so $X\boldsymbol{u}=\boldsymbol{0}$. Hence the two kernels coincide.

As for the ranks, applying the rank–nullity theorem（Theorem 7.3）［Vector Spaces and Linear Maps］ to the two linear maps on $\mathbb{R}^d$ gives

$$\operatorname{rank}(X^{\mathsf{T}}X) = d - \dim\ker(X^{\mathsf{T}}X) = d - \dim\ker X = \operatorname{rank}(X) .$$

Finally, invertibility of $X^{\mathsf{T}}X$ is equivalent to $\ker(X^{\mathsf{T}}X)=\{\boldsymbol{0}\}$, which is equivalent to $\ker X=\{\boldsymbol{0}\}$; and $\ker X=\{\boldsymbol{0}\}$ says that $X\boldsymbol{u}=\boldsymbol{0}$ only for $\boldsymbol{u}=\boldsymbol{0}$, which is exactly linear independence of the columns of $X$. For the relation between linear maps and kernels see Vector spaces and linear transformations.

∎

Corollary 3.4（Closed form for the least-squares solution）

If the column vectors of $X\in\mathbb{R}^{n\times d}$ are linearly independent (that is, $\operatorname{rank}X=d$, which in particular forces $n\ge d$), then $L(\boldsymbol{w})=\frac{1}{2n}\|X\boldsymbol{w}-\boldsymbol{y}\|^2$ has exactly one minimiser, given by

$$\boldsymbol{w}^{\star} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} .$$

Proof（Corollary 3.4）

By Lemma 3.3, $X^{\mathsf{T}}X$ is invertible, so the normal equation $X^{\mathsf{T}}X\boldsymbol{w}=X^{\mathsf{T}}\boldsymbol{y}$ has the unique solution $(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$. By Theorem 3.2 the set of solutions of the normal equation coincides with the set of minimisers of $L$, so this single point is the only minimiser.

∎

Remark 3.5（When the columns are not linearly independent）

Minimisers exist even when $\operatorname{rank}X < d$. Indeed, from $X^{\mathsf{T}}X\boldsymbol{u}=X^{\mathsf{T}}(X\boldsymbol{u})$ we get $\operatorname{Im}(X^{\mathsf{T}}X)\subseteq\operatorname{Im}(X^{\mathsf{T}})$, while Lemma 3.3 gives $\dim\operatorname{Im}(X^{\mathsf{T}}X)=\operatorname{rank}X$, and equality of row rank and column rank gives $\dim\operatorname{Im}(X^{\mathsf{T}})=\operatorname{rank}(X^{\mathsf{T}})=\operatorname{rank}X$. An inclusion between subspaces of equal dimension is an equality, so $\operatorname{Im}(X^{\mathsf{T}}X)=\operatorname{Im}(X^{\mathsf{T}})$, and the right-hand side $X^{\mathsf{T}}\boldsymbol{y}$ always lies in the left-hand side. That is, the normal equation always has a solution.

The solution is not unique, however: given one solution $\boldsymbol{w}^{\star}$, the set of all solutions is the affine subspace $\boldsymbol{w}^{\star}+\ker X$. In practice uniqueness is imposed by ridge regularisation (minimising $L(\boldsymbol{w})+\lambda\|\boldsymbol{w}\|^2$) or by the Moore–Penrose pseudoinverse（Definition 8.6）［Linear Regression and Least Squares］.

The condition $\langle\boldsymbol{u},X^{\mathsf{T}}(X\boldsymbol{w}^{\star}-\boldsymbol{y})\rangle=0$ used in the proof of Theorem 3.2 reads $X^{\mathsf{T}}\boldsymbol{r}=\boldsymbol{0}$, that is, the residual vector $\boldsymbol{r}=X\boldsymbol{w}^{\star}-\boldsymbol{y}$ is orthogonal to every column of $X$. This is the geometric meaning of least squares. Equivalently, $X\boldsymbol{w}^{\star}$ is nothing other than the orthogonal projection of $\boldsymbol{y}$ onto the column space $\operatorname{Im}X$, a viewpoint treated head-on in Theorem 4.2［Linear Regression and Least Squares］.

The geometry of least squares. The prediction Xw* is the orthogonal projection of y onto the column space, and the residual is orthogonal to that space

Example 3.6（Fitting a line to three points, carried through to the end）

Take the data $(x_i,y_i)=(1,2),(2,3),(3,5)$ and the model $y=w_0+w_1x$. Placing the constant feature in the first column,

$$X=\begin{pmatrix}1&1\\1&2\\1&3\end{pmatrix},\qquad \boldsymbol{y}=\begin{pmatrix}2\\3\\5\end{pmatrix},\qquad X^{\mathsf{T}}X=\begin{pmatrix}3&6\\6&14\end{pmatrix},\qquad X^{\mathsf{T}}\boldsymbol{y}=\begin{pmatrix}10\\23\end{pmatrix}$$

(here $3=1+1+1$, $6=1+2+3$, $14=1+4+9$, $10=2+3+5$ and $23=1\cdot 2+2\cdot 3+3\cdot 5$). Since $\det(X^{\mathsf{T}}X)=3\cdot 14-6\cdot 6=6\ne 0$, Corollary 3.4 applies. In the normal equation

$$\begin{cases} 3w_0+6w_1=10\\ 6w_0+14w_1=23\end{cases}$$

multiplying the first equation by $2$ gives $6w_0+12w_1=20$; subtracting this from the second gives $2w_1=3$, that is, $w_1=\frac32$. Substituting back into the first gives $3w_0=10-9=1$, so $w_0=\frac13$. The fitted line is $y=\frac13+\frac32x$.

As a check, look at the residuals. The predicted values are $\frac13+\frac32=\frac{11}{6}$, $\frac13+3=\frac{10}{3}$ and $\frac13+\frac92=\frac{29}{6}$, so

$$\boldsymbol{r}=X\boldsymbol{w}^{\star}-\boldsymbol{y}=\Bigl(-\tfrac16,\ \tfrac13,\ -\tfrac16\Bigr)^{\mathsf{T}} .$$

The inner product with the first column (all $1$s) is $-\frac16+\frac13-\frac16=0$, and with the second column $(1,2,3)^{\mathsf{T}}$ it is $-\frac16+\frac23-\frac12=0$. The residual is indeed orthogonal to both columns. The minimum value is $L(\boldsymbol{w}^{\star})=\frac{1}{6}\bigl(\frac1{36}+\frac19+\frac1{36}\bigr)=\frac{1}{6}\cdot\frac16=\frac{1}{36}$.

Example 3.7（Duplicated features destroy uniqueness of the solution）

For the same data, suppose we accidentally include both "$x$" and "$2x$" as features (for instance, entering the same length twice, as two columns that differ only in the unit used). The design matrix is

$$X'=\begin{pmatrix}1&1&2\\1&2&4\\1&3&6\end{pmatrix} ,$$

and since the third column is twice the second, the columns are linearly dependent. Indeed $X'\boldsymbol{u}=\boldsymbol{0}$ for $\boldsymbol{u}=(0,2,-1)^{\mathsf{T}}$, so by Lemma 3.3 the matrix $X'^{\mathsf{T}}X'$ is not invertible.

The column space of $X'$ is the same as that of $X$ (the third column, being a constant multiple of the second, adds no new direction), so neither the minimum loss value nor the prediction vector differs from Example 3.6. What changes is the number of solutions. The vector $\boldsymbol{w}=(\frac13,\frac32,0)^{\mathsf{T}}$ is one minimiser, but as stated in Remark 3.5, so is every

$$\boldsymbol{w}(t)=\Bigl(\tfrac13,\ \tfrac32+2t,\ -t\Bigr)^{\mathsf{T}},\qquad t\in\mathbb{R} .$$

Indeed the prediction made by this $\boldsymbol{w}(t)$ is $\frac13+(\frac32+2t)x+(-t)(2x)=\frac13+\frac32x$, independent of $t$. The problem of multicollinearity, that reading meaning into the coefficient values themselves is dangerous, is in the language of linear algebra nothing more than the statement that the kernel of the design matrix is nontrivial.

4. Calculus: translating minimisation into the language of gradients

Theorem 3.2 relied heavily on the special form of the squared loss. Change the loss, or make the model nonlinear, and such closed forms are essentially never available. What is needed instead is a method that approaches a minimiser using only local information, namely derivatives. We first make clear what we should be aiming at.

Definition 4.1（Convex function）

A function $f:\mathbb{R}^d\to\mathbb{R}$ is convex if for all $\boldsymbol{v},\boldsymbol{w}\in\mathbb{R}^d$ and all $t\in[0,1]$,

$$f\bigl((1-t)\boldsymbol{w}+t\boldsymbol{v}\bigr) \le (1-t)f(\boldsymbol{w}) + t f(\boldsymbol{v}) .$$

That is, the segment joining any two points on the graph never dips below the graph.

Theorem 4.2（For convex functions, stationary points and global minimisers coincide）

Let $f:\mathbb{R}^d\to\mathbb{R}$ be totally differentiable.

1. If $f$ is convex, then $f(\boldsymbol{v}) \ge f(\boldsymbol{w}) + \langle \nabla f(\boldsymbol{w}),\,\boldsymbol{v}-\boldsymbol{w}\rangle$ for all $\boldsymbol{v},\boldsymbol{w}\in\mathbb{R}^d$.
2. If $f$ is convex, then $\boldsymbol{w}^{\star}$ is a global minimiser of $f$ if and only if $\nabla f(\boldsymbol{w}^{\star})=\boldsymbol{0}$. Moreover, the implication ”$\boldsymbol{w}^{\star}$ is a global minimiser $\implies \nabla f(\boldsymbol{w}^{\star})=\boldsymbol{0}$” holds without assuming convexity.

Proof（Theorem 4.2）

(1) Put $\boldsymbol{u}:=\boldsymbol{v}-\boldsymbol{w}$ and take $t\in(0,1]$. Since $(1-t)\boldsymbol{w}+t\boldsymbol{v}=\boldsymbol{w}+t\boldsymbol{u}$, the inequality of Definition 4.1 reads

$$f(\boldsymbol{w}+t\boldsymbol{u}) \le f(\boldsymbol{w}) + t\bigl(f(\boldsymbol{v})-f(\boldsymbol{w})\bigr) .$$

Moving $f(\boldsymbol{w})$ to the left and dividing by $t>0$,

$$\frac{f(\boldsymbol{w}+t\boldsymbol{u})-f(\boldsymbol{w})}{t} \le f(\boldsymbol{v})-f(\boldsymbol{w}) .$$

By total differentiability, $f(\boldsymbol{w}+t\boldsymbol{u}) = f(\boldsymbol{w}) + t\langle\nabla f(\boldsymbol{w}),\boldsymbol{u}\rangle + o(t)$, so the left-hand side converges to $\langle\nabla f(\boldsymbol{w}),\boldsymbol{u}\rangle$ as $t\to 0^{+}$. The right-hand side is a constant independent of $t$, so passing to the limit gives $\langle\nabla f(\boldsymbol{w}),\boldsymbol{v}-\boldsymbol{w}\rangle \le f(\boldsymbol{v})-f(\boldsymbol{w})$.

(2) First assume $\nabla f(\boldsymbol{w}^{\star})=\boldsymbol{0}$. Taking $\boldsymbol{w}=\boldsymbol{w}^{\star}$ in (1) gives $f(\boldsymbol{v})\ge f(\boldsymbol{w}^{\star})+0$ for every $\boldsymbol{v}$, so $\boldsymbol{w}^{\star}$ is a global minimiser.

Conversely, assume $\boldsymbol{w}^{\star}$ is a global minimiser. Take $\boldsymbol{u}\in\mathbb{R}^d$ arbitrarily and consider the one-variable function $g(t):=f(\boldsymbol{w}^{\star}+t\boldsymbol{u})$. By the chain rule $g$ is differentiable with $g'(0)=\langle\nabla f(\boldsymbol{w}^{\star}),\boldsymbol{u}\rangle$. Since $g$ attains its minimum at $t=0$, the necessary condition for an extremum at an interior point (Fermat's lemma（Lemma 2.5）［Mean Value Theorems and Taylor's Theorem］) gives $g'(0)=0$, that is, $\langle\nabla f(\boldsymbol{w}^{\star}),\boldsymbol{u}\rangle=0$. Taking $\boldsymbol{u}=\nabla f(\boldsymbol{w}^{\star})$ gives $\|\nabla f(\boldsymbol{w}^{\star})\|^2=0$, hence $\nabla f(\boldsymbol{w}^{\star})=\boldsymbol{0}$. Convexity was never used in this direction.

∎

The use of this theorem is plain. For a nonconvex function, “the gradient vanishes” is only a necessary condition, and one may have come to rest at a saddle point or a local minimiser. Convexity removes that worry: “making the gradient vanish” becomes exactly the same thing as “minimising”. The squared loss lies on the good side of this divide.

Remark 4.3（The squared loss is convex）

Put $h(\boldsymbol{z}):=\frac{1}{2n}\|\boldsymbol{z}\|^2$. For $t\in[0,1]$ the identity

$$(1-t)\|\boldsymbol{a}\|^2 + t\|\boldsymbol{b}\|^2 - \|(1-t)\boldsymbol{a}+t\boldsymbol{b}\|^2 = t(1-t)\|\boldsymbol{a}-\boldsymbol{b}\|^2 \ge 0$$

is verified by expansion alone (expanding the left-hand side gives $\bigl[(1-t)-(1-t)^2\bigr]\|\boldsymbol{a}\|^2+(t-t^2)\|\boldsymbol{b}\|^2-2t(1-t)\langle\boldsymbol{a},\boldsymbol{b}\rangle$, and $(1-t)-(1-t)^2=t-t^2=t(1-t)$), so $h$ is convex in the sense of Definition 4.1. Next, putting $\boldsymbol{w}_t:=(1-t)\boldsymbol{w}+t\boldsymbol{v}$,

$$X\boldsymbol{w}_t-\boldsymbol{y} = (1-t)X\boldsymbol{w}+tX\boldsymbol{v}-\boldsymbol{y} = (1-t)(X\boldsymbol{w}-\boldsymbol{y})+t(X\boldsymbol{v}-\boldsymbol{y})$$

(using $(1-t)(-\boldsymbol{y})+t(-\boldsymbol{y})=-\boldsymbol{y}$). That is, an affine map sends convex combinations to convex combinations. Hence $L(\boldsymbol{w}_t)=h(X\boldsymbol{w}_t-\boldsymbol{y})\le(1-t)h(X\boldsymbol{w}-\boldsymbol{y})+th(X\boldsymbol{v}-\boldsymbol{y})=(1-t)L(\boldsymbol{w})+tL(\boldsymbol{v})$, so $L$ too is convex.

Therefore, by part (2) of Theorem 4.2, solving the normal equation $\nabla L(\boldsymbol{w})=\boldsymbol{0}$ and minimising $L$ are equivalent. We proved Theorem 3.2 without using differentiation at all precisely to show that this fact is visible with algebra alone. One and the same fact has three faces: algebraic (the normal equation), geometric (orthogonal projection) and analytic (stationary point).

4.1. Why an iterative method rather than a closed form

Given the formula $\boldsymbol{w}^{\star}=(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$ of Corollary 3.4, it might seem that one need only evaluate it. In practice, iterative methods are used, for two reasons.

• Cost. Forming $X^{\mathsf{T}}X$ takes $O(nd^2)$ operations and solving the system takes $O(d^3)$. Once the number of features $d$ is of the order of $10^5$ or $10^6$, this is out of reach.
• Scope. For models such as neural networks, where $f_{\boldsymbol{w}}$ is nonlinear in $\boldsymbol{w}$, the equation $\nabla L(\boldsymbol{w})=\boldsymbol{0}$ is a nonlinear system and has in general no closed-form solution.

So we descend a little at a time, using only the gradient at the current point. With a constant learning rate $\eta>0$, the rule

$$\boldsymbol{w}_{k+1} = \boldsymbol{w}_{k} - \eta\,\nabla L(\boldsymbol{w}_{k}),\qquad k=0,1,2,\ldots$$

defines gradient descent. Part (1) of Theorem 4.2 can be read as saying that the gradient is the slope of a first-order approximation supporting the function from below, so moving in the opposite direction is a natural choice. But how large a step should we take? For the squared loss the answer can be written down completely.

Theorem 4.4（Convergence of gradient descent for the least-squares loss）

Let the column vectors of $X\in\mathbb{R}^{n\times d}$ be linearly independent, and put $L(\boldsymbol{w})=\frac{1}{2n}\|X\boldsymbol{w}-\boldsymbol{y}\|^2$ and $A:=\frac{1}{n}X^{\mathsf{T}}X$. Then $A$ is symmetric positive definite; list its eigenvalues with multiplicity as $0<\lambda_1\le\lambda_2\le\cdots\le\lambda_d$. Let $\boldsymbol{w}^{\star}$ be the unique minimiser of Corollary 3.4, let $\eta>0$ be a constant, and define the sequence by $\boldsymbol{w}_{k+1}=\boldsymbol{w}_k-\eta\nabla L(\boldsymbol{w}_k)$.

1. One has $\boldsymbol{w}_k\to\boldsymbol{w}^{\star}$ for every initial point $\boldsymbol{w}_0\in\mathbb{R}^d$ if and only if $\eta < 2/\lambda_d$.
2. If $\eta < 2/\lambda_d$, then setting $\rho(\eta):=\max\{|1-\eta\lambda_1|,\,|1-\eta\lambda_d|\}$ we have $\rho(\eta)<1$, and $\|\boldsymbol{w}_k-\boldsymbol{w}^{\star}\| \le \rho(\eta)^k\,\|\boldsymbol{w}_0-\boldsymbol{w}^{\star}\|$ for every $k\ge 0$.
3. The learning rate minimising $\rho(\eta)$ is $\eta^{\star}=\dfrac{2}{\lambda_1+\lambda_d}$, and then $\rho(\eta^{\star})=\dfrac{\kappa-1}{\kappa+1}$, where $\kappa:=\lambda_d/\lambda_1$ is the condition number of $A$.

Proof（Theorem 4.4）

Step 1 ($A$ is symmetric positive definite). From $(X^{\mathsf{T}}X)^{\mathsf{T}}=X^{\mathsf{T}}(X^{\mathsf{T}})^{\mathsf{T}}=X^{\mathsf{T}}X$, the matrix $A$ is symmetric. Moreover, for $\boldsymbol{u}\ne\boldsymbol{0}$ we have $\langle\boldsymbol{u},A\boldsymbol{u}\rangle=\frac{1}{n}\|X\boldsymbol{u}\|^2$, and linear independence of the columns together with Lemma 3.3 gives $X\boldsymbol{u}\ne\boldsymbol{0}$, so this is positive. The eigenvalues of a symmetric matrix are real, and positive definiteness makes them all positive.

Step 2 (recursion for the error). By Lemma 3.1, $\nabla L(\boldsymbol{w})=\frac1n X^{\mathsf{T}}(X\boldsymbol{w}-\boldsymbol{y})=A\boldsymbol{w}-\frac1n X^{\mathsf{T}}\boldsymbol{y}$. Since $\boldsymbol{w}^{\star}$ satisfies the normal equation (Theorem 3.2), we have $\frac1n X^{\mathsf{T}}\boldsymbol{y}=A\boldsymbol{w}^{\star}$, and therefore

$$\nabla L(\boldsymbol{w}) = A(\boldsymbol{w}-\boldsymbol{w}^{\star}) .$$

Writing the error as $\boldsymbol{e}_k:=\boldsymbol{w}_k-\boldsymbol{w}^{\star}$,

$$\boldsymbol{e}_{k+1} = \boldsymbol{w}_k - \eta A\boldsymbol{e}_k - \boldsymbol{w}^{\star} = (I-\eta A)\boldsymbol{e}_k,$$

hence $\boldsymbol{e}_k=(I-\eta A)^k\boldsymbol{e}_0$.

Step 3 (spectral decomposition). Since $A$ is real symmetric, the spectral theorem provides an orthonormal basis $\boldsymbol{q}_1,\ldots,\boldsymbol{q}_d$ of $\mathbb{R}^d$ with $A\boldsymbol{q}_j=\lambda_j\boldsymbol{q}_j$ (Corollary 4.3［スペクトル定理］; for the full statement see The spectral theorem). Expanding $\boldsymbol{e}_0=\sum_{j=1}^d c_j\boldsymbol{q}_j$ with $c_j=\langle\boldsymbol{e}_0,\boldsymbol{q}_j\rangle$ and using $(I-\eta A)\boldsymbol{q}_j=(1-\eta\lambda_j)\boldsymbol{q}_j$ repeatedly gives

$$\boldsymbol{e}_k = \sum_{j=1}^{d} c_j (1-\eta\lambda_j)^k \boldsymbol{q}_j, \qquad \|\boldsymbol{e}_k\|^2 = \sum_{j=1}^{d} c_j^2 (1-\eta\lambda_j)^{2k} .$$

The second identity uses orthonormality of the basis.

Step 4 (claim 1). If $|1-\eta\lambda_j|<1$ for every $j$, then each term of the finite sum above tends to $0$ as $k\to\infty$, so $\|\boldsymbol{e}_k\|\to0$. For $\eta>0$ and $\lambda_j>0$,

$$|1-\eta\lambda_j| < 1 \iff -1 < 1-\eta\lambda_j < 1 \iff 0 < \eta\lambda_j < 2 \iff \eta < 2/\lambda_j ,$$

and requiring this for every $j$ is equivalent to the single condition $\eta<2/\lambda_d$ on the largest eigenvalue. Conversely, if $\eta\ge 2/\lambda_d$ then $|1-\eta\lambda_d|\ge 1$. Taking the initial point $\boldsymbol{w}_0=\boldsymbol{w}^{\star}+\boldsymbol{q}_d$ gives $c_d=1$ and $c_j=0$ otherwise, so $\|\boldsymbol{e}_k\|=|1-\eta\lambda_d|^k\ge 1$, and there exists an initial point from which the iteration fails to converge.

Step 5 (claim 2). Since $\eta>0$, the quantity $1-\eta\lambda_j$ is decreasing in $\lambda_j$, so $1-\eta\lambda_d \le 1-\eta\lambda_j\le 1-\eta\lambda_1$. If a real number $x$ lies in a closed interval $[m,M]$ then $x\le M\le|M|$ and $-x\le -m\le |m|$, so $|x|\le\max\{|m|,|M|\}$. Hence $\max_j|1-\eta\lambda_j|=\rho(\eta)$ (the endpoints being attained at $j=1$ and $j=d$), and the formula of Step 3 gives

$$\|\boldsymbol{e}_k\|^2 \le \rho(\eta)^{2k}\sum_{j}c_j^2 = \rho(\eta)^{2k}\|\boldsymbol{e}_0\|^2 .$$

Taking square roots gives the claim. By Step 4, $\rho(\eta)<1$ when $\eta<2/\lambda_d$.

Step 6 (claim 3). We examine $\rho(\eta)=\max\{|1-\eta\lambda_1|,|1-\eta\lambda_d|\}$ for $\eta>0$. From $\lambda_1\le\lambda_d$ we get $1/\lambda_d\le 1/\lambda_1$, and there are three cases.

• For $0<\eta\le 1/\lambda_d$, both $1-\eta\lambda_1$ and $1-\eta\lambda_d$ are nonnegative, so $\rho(\eta)=1-\eta\lambda_1$, which is strictly decreasing in $\eta$.
• For $1/\lambda_d\le\eta\le 1/\lambda_1$, we have $\rho(\eta)=\max\{1-\eta\lambda_1,\ \eta\lambda_d-1\}$. The first expression decreases and the second increases, so the maximum is smallest where the two are equal. Solving $1-\eta\lambda_1=\eta\lambda_d-1$ gives $\eta=2/(\lambda_1+\lambda_d)$, and from $\lambda_1\le\lambda_d$ this value lies in the interval $[1/\lambda_d,1/\lambda_1]$ (indeed $2/(\lambda_1+\lambda_d)\ge 1/\lambda_d$ is equivalent to $2\lambda_d\ge\lambda_1+\lambda_d$, and $2/(\lambda_1+\lambda_d)\le 1/\lambda_1$ is equivalent to $2\lambda_1\le\lambda_1+\lambda_d$, both of which hold).
• For $\eta\ge 1/\lambda_1$, we have $\rho(\eta)=\max\{\eta\lambda_1-1,\ \eta\lambda_d-1\}=\eta\lambda_d-1$, strictly increasing.

Hence $\rho$ attains its minimum at $\eta^{\star}=2/(\lambda_1+\lambda_d)$, with value

$$\rho(\eta^{\star}) = 1-\frac{2\lambda_1}{\lambda_1+\lambda_d} = \frac{\lambda_d-\lambda_1}{\lambda_d+\lambda_1} = \frac{\kappa-1}{\kappa+1}$$

(the last step divides numerator and denominator by $\lambda_1$).

∎

This theorem answers both questions 1 and 2 raised at the outset. Raising the learning rate too far makes the iteration diverge because it crosses into $\eta\ge 2/\lambda_d$, and the speed of convergence is governed by the condition number $\kappa$ alone. When $\kappa$ is large, $\rho=(\kappa-1)/(\kappa+1)$ approaches $1$ and the number of iterations required grows roughly in proportion to $\kappa$. Let us see this in numbers.

Example 4.5（The condition number decides the number of iterations）

For the data of Example 3.6, $A=\frac13X^{\mathsf{T}}X=\begin{pmatrix}1&2\\2&14/3\end{pmatrix}$. The characteristic polynomial is $\lambda^2-\frac{17}{3}\lambda+\frac23=0$ (the trace is $1+\frac{14}{3}=\frac{17}{3}$ and the determinant is $\frac{14}{3}-4=\frac23$), so

$$\lambda = \frac{17\pm\sqrt{265}}{6},\qquad \lambda_1\approx 0.1202,\quad \lambda_2\approx 5.5465 .$$

By Theorem 4.4, the convergence condition is $\eta < 2/5.5465\approx 0.3606$ and the optimal learning rate is $\eta^{\star}=2/(17/3)=6/17\approx 0.3529$. The condition number is $\kappa\approx 46.1$, so $\rho\approx 45.1/47.1\approx 0.9576$. To shrink the error by a factor of $10^{-3}$,

$$\rho^k\le 10^{-3} \iff k \ge \frac{3\ln 10}{-\ln\rho} \approx \frac{6.908}{0.0433} \approx 159.4 ,$$

so $160$ iterations are needed. Moreover the gap between the upper limit $0.3606$ for convergence and the optimal value $0.3529$ is tiny, so a slightly greedy choice makes the iteration diverge.

Now centre the input, setting $\tilde{x}=x-2$ (the mean of $x$ is $2$). The second column of the design matrix becomes $(-1,0,1)^{\mathsf{T}}$, whose inner product with the first column is $0$, so

$$\tilde{X}^{\mathsf{T}}\tilde{X}=\begin{pmatrix}3&0\\0&2\end{pmatrix},\qquad A=\operatorname{diag}\Bigl(1,\ \tfrac23\Bigr),\qquad \kappa=\frac{1}{2/3}=1.5 .$$

Now $\eta^{\star}=2/(1+\frac23)=\frac65=1.2$ and $\rho=(1.5-1)/(1.5+1)=0.2$, so the number of iterations required is $k\ge 3\ln 10/\ln 5\approx 4.3$, that is, $5$.

	without centring	with centring
eigenvalues	$0.120,\ 5.547$	$0.667,\ 1.000$
condition number $\kappa$	$46.1$	$1.5$
learning rates that diverge	$\eta\ge 0.361$	$\eta\ge 2$
optimal learning rate $\eta^{\star}$	$0.353$	$1.2$
convergence rate $\rho$	$0.958$	$0.2$
iterations to shrink the error by $10^{-3}$	$160$	$5$

The fitted line is still $y=\frac13+\frac32x$; all we changed was the choice of coordinates. That alone cut the number of iterations by a factor of $32$. This is the answer to question 2.

Tip

Standardising features in practice is not about tidiness of appearance; it is about reducing the condition number $\kappa$ of Theorem 4.4 and thereby speeding up convergence. The choice of units (centimetres or metres) does not change the fitted function itself — the predictions and the minimum value of the loss are the same either way — but it does change the length of the route that takes us there.

5. Probability and statistics: where the loss comes from and what it guarantees

So far the loss function has been handed down from above. Why the square? Why not the absolute value? And why the cross-entropy for classification? Probability gives a principled answer to these questions.

Definition 5.1（Maximum likelihood estimation）

Given a probability model $p(\cdot\mid\boldsymbol{\theta})$ with parameter $\boldsymbol{\theta}$ and observed data $\boldsymbol{y}$, the map $\boldsymbol{\theta}\mapsto p(\boldsymbol{y}\mid\boldsymbol{\theta})$, regarded as a function of $\boldsymbol{\theta}$, is called the likelihood function. A maximiser $\hat{\boldsymbol{\theta}}$ of the likelihood is called a maximum likelihood estimator. Since the logarithm is strictly increasing, this is equivalent to maximising the log-likelihood $\log p(\boldsymbol{y}\mid\boldsymbol{\theta})$, or to minimising the negative log-likelihood.

Theorem 5.2（Under Gaussian noise, maximum likelihood coincides with least squares）

Let $\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n\in\mathbb{R}^d$ be fixed (nonrandom) inputs, let $\sigma>0$ be a known constant, and suppose the observations are generated by

$$y_i = \langle\boldsymbol{w},\boldsymbol{x}_i\rangle + \varepsilon_i \qquad (i=1,\ldots,n),$$

where $\varepsilon_1,\ldots,\varepsilon_n$ are independent with the normal distribution $\mathcal{N}(0,\sigma^2)$. Then the set of maximum likelihood estimators of $\boldsymbol{w}$ coincides with the set of minimisers of the squared loss $L(\boldsymbol{w})=\frac{1}{2n}\|X\boldsymbol{w}-\boldsymbol{y}\|^2$.

Proof（Theorem 5.2）

Since $\varepsilon_i\sim\mathcal{N}(0,\sigma^2)$, for fixed $\boldsymbol{w}$ the variable $y_i$ is normally distributed with mean $\langle\boldsymbol{w},\boldsymbol{x}_i\rangle$ and variance $\sigma^2$. The $\varepsilon_i$ are independent, so $y_1,\ldots,y_n$ are independent too and the joint density is a product.

$$p(\boldsymbol{y}\mid\boldsymbol{w}) = \prod_{i=1}^{n}\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{\bigl(y_i-\langle\boldsymbol{w},\boldsymbol{x}_i\rangle\bigr)^2}{2\sigma^2}\right)$$

Taking logarithms turns the product into a sum:

$$\log p(\boldsymbol{y}\mid\boldsymbol{w}) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{n}\bigl(y_i-\langle\boldsymbol{w},\boldsymbol{x}_i\rangle\bigr)^2 = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\|X\boldsymbol{w}-\boldsymbol{y}\|^2 .$$

The first term is a constant independent of $\boldsymbol{w}$, and the coefficient $\frac{1}{2\sigma^2}$ in the second is positive. Hence the set of $\boldsymbol{w}$ maximising $\log p$ and the set of $\boldsymbol{w}$ minimising $\|X\boldsymbol{w}-\boldsymbol{y}\|^2$ are exactly the same. Multiplying by the positive constant $\frac{1}{2n}$ does not change the minimisers, so this set coincides with the set of minimisers of $L$.

∎

In other words, the choice to “use the squared loss” is equivalent to the assumption that the errors are normally distributed, have the same spread at every data point, and are mutually independent. Choosing a loss function is implicitly choosing a probability model. If the assumptions do not match reality, the loss should be changed; for data contaminated by outliers, for instance, it is natural to assume a heavy-tailed distribution (Exercise 7.3).

Example 5.3（The cross-entropy comes from the Bernoulli distribution）

Consider binary classification. The labels are $y_i\in\{0,1\}$ and the model outputs the probability that the label is $1$:

$$p_i := \sigma(z_i),\qquad z_i := \langle\boldsymbol{w},\boldsymbol{x}_i\rangle,\qquad \sigma(z)=\frac{1}{1+e^{-z}}$$

(here $\sigma$ is the sigmoid function, whose range is the open interval $(0,1)$). Assuming the $y_i$ are independent with Bernoulli distribution $\mathrm{Be}(p_i)$, and using $y_i\in\{0,1\}$, the probability of a single point can be packed into the single expression $p_i^{y_i}(1-p_i)^{1-y_i}$ (which returns $p_i$ when $y_i=1$ and $1-p_i$ when $y_i=0$). Hence

$$-\frac{1}{n}\log p(\boldsymbol{y}\mid\boldsymbol{w}) = -\frac{1}{n}\sum_{i=1}^{n}\Bigl[y_i\log p_i + (1-y_i)\log(1-p_i)\Bigr] ,$$

and the right-hand side is precisely the cross-entropy loss. This is the answer to question 4. The difference between the squared loss and the cross-entropy is nothing more than a difference of probability model: whether the output is viewed as a real number or as a probability.

Let us also compute the gradient. From $\sigma'(z)=\sigma(z)(1-\sigma(z))$ we get $\frac{d}{dz}\log\sigma(z)=1-\sigma(z)$ and $\frac{d}{dz}\log(1-\sigma(z))=-\sigma(z)$, so for the loss at a single point, $\ell_i=-[y_i\log p_i+(1-y_i)\log(1-p_i)]$,

$$\frac{\partial \ell_i}{\partial z_i} = -\bigl[y_i(1-p_i) - (1-y_i)p_i\bigr] = -\bigl[y_i - y_ip_i - p_i + y_ip_i\bigr] = p_i - y_i .$$

By the chain rule, $\frac{\partial z_i}{\partial\boldsymbol{w}}=\boldsymbol{x}_i$, and therefore

$$\nabla \hat{R}_n(\boldsymbol{w}) = \frac{1}{n}\sum_{i=1}^n (p_i-y_i)\boldsymbol{x}_i = \frac{1}{n}X^{\mathsf{T}}(\boldsymbol{p}-\boldsymbol{y}),\qquad \boldsymbol{p}=(p_1,\ldots,p_n)^{\mathsf{T}} .$$

This has the same shape as $\nabla L(\boldsymbol{w})=\frac1n X^{\mathsf{T}}(X\boldsymbol{w}-\boldsymbol{y})$ from Lemma 3.1. The structure “apply the transpose of the design matrix to the residual” is common to linear regression and logistic regression. This form of the gradient is established again in Theorem 5.1［Logistic Regression］; the details are treated in Logistic regression.

5.1. Why the training error may be trusted, and why it may not

The other reason probability is needed is to justify the substitution of $\hat{R}_n$ for $R$ in Definition 2.2.

Proposition 5.4（For a fixed predictor, the empirical risk concentrates around the expected risk）

In the setting of Definition 2.1, let $(\boldsymbol{X}_1,Y_1),\ldots,(\boldsymbol{X}_n,Y_n)$ be drawn independently from the distribution $P$, all with the same distribution. Suppose the predictor $f$ is fixed in advance and does not depend on the training data, and that the random variables $Z_i:=\ell(f(\boldsymbol{X}_i),Y_i)$ have finite expectation $R(f)$ and finite variance $\sigma_{\ell}^2$. Then the following hold.

1. $\mathbb{E}\bigl[\hat{R}_n(f)\bigr]=R(f)$ (unbiasedness).
2. $\operatorname{Var}\bigl[\hat{R}_n(f)\bigr]=\sigma_{\ell}^2/n$.
3. For every $t>0$, $\Pr\bigl(|\hat{R}_n(f)-R(f)|\ge t\bigr) \le \dfrac{\sigma_{\ell}^2}{n t^2}$.

Proof（Proposition 5.4）

We have $\hat{R}_n(f)=\frac1n\sum_{i=1}^n Z_i$. The $Z_1,\ldots,Z_n$ are independent and identically distributed with $\mathbb{E}[Z_i]=\mathbb{E}[\ell(f(\boldsymbol{X}_i),Y_i)]=R(f)$ (since $f$ does not depend on the data, the distribution of $Z_i$ is determined by $P$ and $f$ alone).

1. By linearity of expectation, $\mathbb{E}[\hat{R}_n(f)]=\frac1n\sum_i\mathbb{E}[Z_i]=\frac1n\cdot nR(f)=R(f)$.
2. The variance of a sum of independent random variables is the sum of the variances, so $\operatorname{Var}[\sum_i Z_i]=n\sigma_{\ell}^2$; and since $\operatorname{Var}[cW]=c^2\operatorname{Var}[W]$, we get $\operatorname{Var}[\hat{R}_n(f)]=\frac{1}{n^2}\cdot n\sigma_{\ell}^2=\sigma_{\ell}^2/n$.
3. Applying Chebyshev's inequality（Corollary 6.2）［Random Variables and Expectation］ to the random variable $\hat{R}_n(f)$, with mean $R(f)$ and variance $\sigma_{\ell}^2/n$, gives the bound immediately.

For linearity of expectation, additivity of variance for independent random variables, and Chebyshev’s inequality, see Random variables and expectation and The law of large numbers and the central limit theorem.

∎

Claim 3 says that increasing $n$ brings the training error close to the generalisation error in the sense of probability. Since $n$ appears in the denominator, it also yields a quantitative guideline: to halve the accuracy $t$, take four times as many samples. This is the most naive mathematical backing for the rule of thumb that more data is better.

Remark 5.5（What breaks when the assumption that f is independent of the data is dropped）

Given $n$ distinct points $x_1,\ldots,x_n\in\mathbb{R}$, there always exists a polynomial of degree at most $n-1$ passing through all $n$ points (Lagrange interpolation; the coefficient matrix of the linear system determining the coefficients is a Vandermonde matrix, whose determinant is $\prod_{i<j}(x_j-x_i)\ne0$ when the points are distinct, so the matrix is invertible). Choosing this polynomial makes the empirical risk exactly $0$.

The expected risk of this predictor, however, is in general enormous, and it oscillates wildly between the training points. This looks like a contradiction with Proposition 5.4, but it is not. The polynomial was determined after looking at the training data, so the hypothesis of the proposition, that $f$ does not depend on the data, fails. When $f$ is chosen as a function of the data, neither the independence of the $Z_i$ nor the identity $\mathbb{E}[Z_i]=R(f)$ holds.

This is what overfitting really is, and it is the answer to question 3. As long as we select from within a hypothesis class $\mathcal{H}$, what we actually need is not a bound for each individual $f$ but a uniform bound on $\sup_{f\in\mathcal{H}}|\hat{R}_n(f)-R(f)|$. Controlling this quantity by the “size” of the hypothesis class is the subject of statistical learning theory (VC dimension, Rademacher complexity, and so on).

Caution

Accuracy measured on the training data becomes an optimistic figure the moment that same data is used to select the model. Setting aside validation data for tuning hyperparameters is a procedure for artificially restoring the hypothesis of Proposition 5.4.

6. The division of labour among the three areas, and what comes next

All four questions raised at the outset now have answers.

1. Raising the learning rate makes the iteration diverge. Because at $\eta\ge 2/\lambda_d$ the error along the direction of the largest eigenvalue is amplified (Theorem 4.4).
2. Changing the units speeds up convergence. Because the condition number $\kappa$ changes and the convergence rate $(\kappa-1)/(\kappa+1)$ improves (Example 4.5).
3. The training error is $0$ yet the model misses. Because unbiasedness of the empirical risk is lost the moment the model is chosen after looking at the data (Remark 5.5).
4. Regression and classification use different losses. Because the probability model placed on the output is normal in one case and Bernoulli in the other (Theorem 5.2, Example 5.3).

The division of labour among the three areas of mathematics can be summarised as follows.

Stage of learning	Mathematics chiefly used	Concrete instance seen here	Read on
Representing data and models	linear algebra (matrices, subspaces, orthogonal projection)	design matrix, normal equation, orthogonality of the residual	Linear regression and least squares
Compressing and diagnosing representations	linear algebra (eigenvalues, symmetric matrices)	the condition number $\kappa$ decides the iteration count	Principal component analysis
Minimising the loss	calculus (gradients, convexity, limits)	the bound $2/\lambda_d$ on the learning rate and its optimal value	Gradient descent
Computing gradients in deep models	calculus (chain rule)	the derivative $p-y$ through the sigmoid	Neural networks and backpropagation
Designing the loss	probability and statistics (likelihood)	squared loss from Gaussian noise, cross-entropy from Bernoulli	Logistic regression
Assessing uncertainty	probability and statistics (expectation, law of large numbers)	unbiasedness of the empirical risk and its variance $\sigma_{\ell}^2/n$	Probability and Bayesian statistics

The prerequisites are, for linear algebra, Vector spaces and linear transformations, Inner product spaces and Gram–Schmidt orthogonalisation and Eigenvalues and eigenvectors; for calculus, Differentiation of multivariable functions and partial derivatives and The mean value theorem and Taylor’s theorem; and for probability, Probability spaces and Kolmogorov’s axioms and Random variables and expectation. No deep excursion is required. All this article actually used was the relation between the inner product and the transpose, the rank–nullity theorem, the spectral theorem for symmetric matrices, the definition of total differentiability, linearity of expectation, and Chebyshev’s inequality.

7. Exercises

Exercise 7.1Easy

For the data $(x_i,y_i)=(0,1),(1,1),(2,4),(3,4)$, find the least-squares solution for the model $y=w_0+w_1x$ from the normal equation. Then verify that the resulting residual vector is orthogonal to both columns of the design matrix.

Solution

The design matrix and target vector are

$$X=\begin{pmatrix}1&0\\1&1\\1&2\\1&3\end{pmatrix},\qquad \boldsymbol{y}=\begin{pmatrix}1\\1\\4\\4\end{pmatrix} .$$

Then $X^{\mathsf{T}}X=\begin{pmatrix}4&6\\6&14\end{pmatrix}$ (the $4$ is the number of rows, $6=0+1+2+3$, $14=0+1+4+9$) and $X^{\mathsf{T}}\boldsymbol{y}=\begin{pmatrix}10\\21\end{pmatrix}$ ($10=1+1+4+4$, $21=0\cdot1+1\cdot1+2\cdot4+3\cdot4$). Since $\det(X^{\mathsf{T}}X)=56-36=20\ne0$, the solution is unique by Corollary 3.4. The normal equation is

$$\begin{cases}4w_0+6w_1=10\\ 6w_0+14w_1=21\end{cases}$$

Multiplying the first equation by $\frac32$ gives $6w_0+9w_1=15$; subtracting from the second gives $5w_1=6$, so $w_1=\frac65$. The first equation then gives $4w_0=10-6\cdot\frac65=10-\frac{36}{5}=\frac{14}{5}$, so $w_0=\frac{7}{10}$. The line is $y=0.7+1.2x$.

The predicted values are $0.7,\ 1.9,\ 3.1,\ 4.3$, so the residual is $\boldsymbol{r}=X\boldsymbol{w}^{\star}-\boldsymbol{y}=(-0.3,\ 0.9,\ -0.9,\ 0.3)^{\mathsf{T}}$. Its inner product with the first column $(1,1,1,1)^{\mathsf{T}}$ is $-0.3+0.9-0.9+0.3=0$, and with the second column $(0,1,2,3)^{\mathsf{T}}$ it is $0+0.9-1.8+0.9=0$. This is nothing but the componentwise form of $X^{\mathsf{T}}(X\boldsymbol{w}^{\star}-\boldsymbol{y})=\boldsymbol{0}$ from Theorem 3.2.

Exercise 7.2Standard

Let $L(\boldsymbol{w})=\frac12\bigl(w_1^2+100w_2^2\bigr)$. This is $L(\boldsymbol{w})=\frac12\langle\boldsymbol{w},A\boldsymbol{w}\rangle$ for $A=\operatorname{diag}(1,100)$, and the minimiser is the origin.

1. Find the range of learning rates $\eta$ for which gradient descent converges from every initial point.
2. Find the $\eta$ giving fastest convergence, and the corresponding convergence rate $\rho$.
3. Using the $\eta$ from part 2, find a number of iterations sufficient to guarantee $\|\boldsymbol{w}_k\|\le 10^{-3}\|\boldsymbol{w}_0\|$.

Solution

We have $\nabla L(\boldsymbol{w})=(w_1,100w_2)^{\mathsf{T}}=A\boldsymbol{w}$, and since $A$ is diagonal its eigenvalues are simply $\lambda_1=1$ and $\lambda_2=100$, with the standard basis as eigenvectors. Theorem 4.4 applies directly (with $\boldsymbol{w}^{\star}=\boldsymbol{0}$).

1. By claim 1, $0<\eta<2/\lambda_2=2/100=0.02$.
2. By claim 3, $\eta^{\star}=\dfrac{2}{1+100}=\dfrac{2}{101}\approx 0.0198$, and since $\kappa=100$ we get $\rho=\dfrac{100-1}{100+1}=\dfrac{99}{101}\approx 0.9802$.
3. By the bound in claim 2 it suffices that $\rho^k\le 10^{-3}$. Taking logarithms gives $k\ln\rho\le -3\ln 10$, and since $\ln\rho<0$,

$$k \ge \frac{3\ln 10}{\ln(101/99)} \approx \frac{6.9078}{0.0200} \approx 345.4 ,$$

so $346$ iterations suffice. A condition number of $100$, which is hardly extreme, already forces several hundred iterations. The value of putting features on a common scale shows itself here.

Exercise 7.3Standard

In the same setting as Theorem 5.2, suppose the noise $\varepsilon_i$ is independent with the Laplace density $p(\varepsilon)=\dfrac{1}{2b}\exp\bigl(-|\varepsilon|/b\bigr)$, where $b>0$ is known.

1. Derive which function maximum likelihood estimation of $\boldsymbol{w}$ amounts to minimising.
2. Explain, in terms of how the loss grows with the residual, why the resulting loss is less sensitive to outliers than the squared loss.

Solution

1. By independence the joint density is a product:

$$p(\boldsymbol{y}\mid\boldsymbol{w}) = \prod_{i=1}^{n}\frac{1}{2b}\exp\left(-\frac{\bigl|y_i-\langle\boldsymbol{w},\boldsymbol{x}_i\rangle\bigr|}{b}\right) .$$

Taking the negative logarithm,

$$-\log p(\boldsymbol{y}\mid\boldsymbol{w}) = n\log(2b) + \frac{1}{b}\sum_{i=1}^{n}\bigl|y_i-\langle\boldsymbol{w},\boldsymbol{x}_i\rangle\bigr| .$$

The first term is a constant independent of $\boldsymbol{w}$ and the coefficient $1/b$ of the second is positive, so by Definition 5.1 maximum likelihood is equivalent to minimising $\sum_i|y_i-\langle\boldsymbol{w},\boldsymbol{x}_i\rangle|$, that is, to the absolute-value loss (least absolute deviations regression).

2. For the residual $r_i=\langle\boldsymbol{w},\boldsymbol{x}_i\rangle-y_i$, the contribution of the squared loss is $\frac12 r_i^2$, its derivative with respect to the residual is $r_i$, and its contribution to the gradient is $r_i\boldsymbol{x}_i$. Moving one observation far away makes the contribution grow without bound in proportion to $|r_i|$, so that single point drags the solution around. The contribution of the absolute-value loss, by contrast, is $|r_i|$, its derivative at $r_i\ne0$ is $\operatorname{sign}(r_i)=\pm1$, and its contribution to the gradient is $\pm\boldsymbol{x}_i$, which is bounded. However far out a point lies, the size of its influence is capped.

In the language of probability models: the normal density decays rapidly as $e^{-r^2/2\sigma^2}$, so it judges large residuals to be all but impossible and works hard to exclude them, whereas the Laplace distribution has heavy tails, $e^{-|r|/b}$, and regards large residuals as reasonably possible. Note, though, that the absolute-value loss is not differentiable at $r_i=0$, so optimisation requires subgradient methods or a reformulation as a linear program.

References

• M. P. Deisenroth, A. A. Faisal, C. S. Ong, Mathematics for Machine Learning, Cambridge University Press, 2020 — Chapter 2 (linear algebra), Chapter 5 (vector calculus), Chapter 9 (linear regression). A full version released by the authors is available at mml-book.github.io.
• C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006 — Chapter 1 (decision theory and loss functions), Chapter 3 (linear regression and maximum likelihood), Chapter 4 (classification and logistic regression).
• I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016 — Part I (Chapter 2 linear algebra, Chapter 3 probability and information theory, Chapter 4 numerical computation, Chapter 5 machine learning basics). Available at deeplearningbook.org.
• S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004 — Chapter 3 (convex functions and the first-order condition), Chapter 9 (unconstrained minimisation and the convergence analysis of gradient descent). Available at web.stanford.edu/~boyd/cvxbook/.
• S. Shalev-Shwartz, S. Ben-David, Understanding Machine Learning: From Theory to Algorithms, Cambridge University Press, 2014 — Chapter 2 (empirical risk minimisation and overfitting), Chapter 4 (uniform convergence).
• S. M. Stigler, The History of Statistics: The Measurement of Uncertainty before 1900, Harvard University Press, 1986 — Part I treats the story of Legendre and Gauss around the method of least squares.