# Why Machine Learning Needs Mathematics: Rewriting Learning as Loss Minimisation

> Regression and classification both reduce to a single optimisation problem, empirical risk minimisation. We show where linear algebra, calculus and probability enter, with complete proofs for least squares and gradient descent.
> https://rikai.mugen-giken.com/en/computer-science/math-for-ml/why-math-for-ml

## 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 (<Ref to="thm-normal-equation" />). 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 (<Ref to="thm-convex-stationary" />), 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$ (<Ref to="thm-gd-quadratic" />).
- Probability answers the question "why the squared loss?". Maximum likelihood estimation under Gaussian noise coincides exactly with least squares (<Ref to="thm-mle-gaussian" />), while assuming a Bernoulli distribution produces the cross-entropy (<Ref to="ex-cross-entropy" />).
- Probability also explains why we may use the error measured on our finitely many samples in place of the error on unseen data (<Ref to="prop-erm-consistency" />). 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 id="def-supervised-setup" title="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.

</Definition>

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 <Ref to="ex-cross-entropy" />.

<Definition id="def-risks" title="Expected risk and empirical risk">

In the setting of <Ref to="def-supervised-setup" />, 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.

</Definition>

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 <Ref to="prop-erm-consistency" />.

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.

<Figure caption="The basic pipeline of supervised learning, with the mathematics used at each stage">
<Mermaid code={`flowchart TD
  A["1. Data: collect n pairs of input and correct answer"] --> B["2. Representation: arrange them into a design matrix X and target vector y — linear algebra"]
  B --> C["3. Hypothesis class: set up predictors f_w indexed by a parameter w"]
  C --> D["4. Loss: quantify the misfit as the empirical risk L(w) — probability and statistics"]
  D --> E["5. Optimisation: move w in a direction that decreases L — calculus"]
  E --> F["6. Evaluation: estimate the error R on unseen data — probability and statistics"]
  F --> C`} />
</Figure>

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 id="lem-expansion" title="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})$.

</Lemma>

<Proof of="lem-expansion">

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 <Ref to="mathematics/linear-algebra/inner-product-spaces#thm-cauchy-schwarz" text="Cauchy–Schwarz inequality" />, 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})$.

</Proof>

<Theorem id="thm-normal-equation" title="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.

</Theorem>

<Proof of="thm-normal-equation">

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 <Ref to="lem-expansion" /> 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}$.

</Proof>

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 id="lem-gram-kernel" title="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.

</Lemma>

<Proof of="lem-gram-kernel">

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 <Ref to="mathematics/linear-algebra/vector-spaces#thm-rank-nullity" text="rank–nullity theorem" /> 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](/en/mathematics/linear-algebra/vector-spaces).

</Proof>

<Corollary id="cor-unique-solution" title="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} .
$$

</Corollary>

<Proof of="cor-unique-solution">

By <Ref to="lem-gram-kernel" />, $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 <Ref to="thm-normal-equation" /> the set of solutions of the normal equation coincides with the set of minimisers of $L$, so this single point is the only minimiser.

</Proof>

<Remark id="rem-rank-deficient" title="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 <Ref to="lem-gram-kernel" /> 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 <Ref to="computer-science/math-for-ml/linear-regression#def-pseudoinverse" text="Moore–Penrose pseudoinverse" />.

</Remark>

The condition $\langle\boldsymbol{u},X^{\mathsf{T}}(X\boldsymbol{w}^{\star}-\boldsymbol{y})\rangle=0$ used in the proof of <Ref to="thm-normal-equation" /> 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 <Ref to="computer-science/math-for-ml/linear-regression#thm-projection" />.

<Figure caption="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">
<svg viewBox="0 0 640 330" width="100%" role="img" aria-label="Orthogonal projection of the vector y onto the column space of the design matrix">
  <defs>
    <marker id="ls-head" viewBox="0 0 10 10" refX="9" refY="5" markerWidth="6" markerHeight="6" orient="auto-start-reverse">
      <path d="M 0 0 L 10 5 L 0 10 z" fill="currentColor" />
    </marker>
    <marker id="ls-head-accent" viewBox="0 0 10 10" refX="9" refY="5" markerWidth="6" markerHeight="6" orient="auto-start-reverse">
      <path d="M 0 0 L 10 5 L 0 10 z" fill="var(--sl-color-accent)" />
    </marker>
  </defs>
  <polygon points="70,250 300,200 570,250 340,300" fill="currentColor" fill-opacity="0.06" stroke="currentColor" stroke-opacity="0.45" stroke-width="1.5" />
  <text x="86" y="284" font-size="16" fill="currentColor">Im X</text>
  <text x="86" y="304" font-size="13" fill="currentColor" fill-opacity="0.75">all predictions the model can reach</text>
  <line x1="210" y1="252" x2="430" y2="240" stroke="currentColor" stroke-width="2" marker-end="url(#ls-head)" />
  <line x1="210" y1="252" x2="430" y2="80" stroke="var(--sl-color-accent)" stroke-width="2.5" marker-end="url(#ls-head-accent)" />
  <line x1="430" y1="80" x2="430" y2="240" stroke="currentColor" stroke-width="2" stroke-dasharray="6 5" marker-end="url(#ls-head)" />
  <path d="M 430 222 L 412 223 L 413 241" fill="none" stroke="currentColor" stroke-opacity="0.7" stroke-width="1.5" />
  <circle cx="210" cy="252" r="4" fill="currentColor" />
  <text x="186" y="276" font-size="15" fill="currentColor">0</text>
  <text x="437" y="70" font-size="16" fill="var(--sl-color-accent)">y (observed target)</text>
  <text x="322" y="274" font-size="16" fill="currentColor">Xw* (prediction)</text>
  <text x="444" y="162" font-size="15" fill="currentColor">residual r = Xw* - y</text>
</svg>
</Figure>

<Example id="ex-fit-three-points" title="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$, <Ref to="cor-unique-solution" /> 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>

<Example id="ex-collinear" title="Duplicated features destroy uniqueness of the solution">

For the same data, suppose we accidentally include both "$x$" and "$2x$" as features (for instance, entering a height as two columns, one in centimetres and one in metres). 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 <Ref to="lem-gram-kernel" /> 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 <Ref to="ex-fit-three-points" />. What changes is the number of solutions. The vector $\boldsymbol{w}=(\frac13,\frac32,0)^{\mathsf{T}}$ is one minimiser, but as stated in <Ref to="rem-rank-deficient" />, 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.

</Example>

## 4. Calculus: translating minimisation into the language of gradients

<Ref to="thm-normal-equation" /> 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 id="def-convex" title="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.

</Definition>

<Theorem id="thm-convex-stationary" title="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.

</Theorem>

<Proof of="thm-convex-stationary">

(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 <Ref to="def-convex" /> 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 (<Ref to="mathematics/calculus/mean-value-and-taylor#lem-fermat" text="Fermat's lemma" />) 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.

</Proof>

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 id="rem-convexity-of-l" title="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 <Ref to="def-convex" />. 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 <Ref to="thm-convex-stationary" />, solving the normal equation $\nabla L(\boldsymbol{w})=\boldsymbol{0}$ and minimising $L$ are equivalent. We proved <Ref to="thm-normal-equation" /> 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).

</Remark>

### 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 <Ref to="cor-unique-solution" />, 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 <Ref to="thm-convex-stationary" /> 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 id="thm-gd-quadratic" title="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 <Ref to="cor-unique-solution" />, 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$.

</Theorem>

<Proof of="thm-gd-quadratic">

**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 <Ref to="lem-gram-kernel" /> 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 <Ref to="lem-expansion" />, $\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 (<Ref to="thm-normal-equation" />), 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$ (<Ref to="mathematics/linear-algebra/spectral-theorem#cor-real-symmetric" />; for the full statement see [The spectral theorem](/mathematics/linear-algebra/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$).

</Proof>

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 id="ex-learning-rate" title="The condition number decides the number of iterations">

For the data of <Ref to="ex-fit-three-points" />, $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 <Ref to="thm-gd-quadratic" />, 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.

</Example>

<Aside type="tip">
Standardising features in practice is not about tidiness of appearance; it is about reducing the condition number $\kappa$ of <Ref to="thm-gd-quadratic" /> and thereby speeding up convergence. The choice of units (centimetres or metres) does not change the minimiser itself, but it does change the length of the route to it.
</Aside>

## 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 id="def-mle" title="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.

</Definition>

<Theorem id="thm-mle-gaussian" title="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$.

</Theorem>

<Proof of="thm-mle-gaussian">

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$.

</Proof>

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 (<Ref to="exr-laplace-loss" />).

<Example id="ex-cross-entropy" title="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 <Ref to="lem-expansion" />. 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 <Ref to="computer-science/math-for-ml/logistic-regression#thm-gradient" />; the details are treated in [Logistic regression](/en/computer-science/math-for-ml/logistic-regression).

</Example>

### 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 <Ref to="def-risks" />.

<Proposition id="prop-erm-consistency" title="For a fixed predictor, the empirical risk concentrates around the expected risk">

In the setting of <Ref to="def-supervised-setup" />, 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}$.

</Proposition>

<Proof of="prop-erm-consistency">

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 <Ref to="mathematics/probability/random-variables#cor-chebyshev" text="Chebyshev's inequality" /> 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](/en/mathematics/probability/random-variables) and [The law of large numbers and the central limit theorem](/en/mathematics/probability/limit-theorems).

</Proof>

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 id="rem-overfitting" title="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 <Ref to="prop-erm-consistency" />, 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).

</Remark>

<Aside type="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 <Ref to="prop-erm-consistency" />.
</Aside>

## 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 (<Ref to="thm-gd-quadratic" />).
2. Changing the units speeds up convergence. Because the condition number $\kappa$ changes and the convergence rate $(\kappa-1)/(\kappa+1)$ improves (<Ref to="ex-learning-rate" />).
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 (<Ref to="rem-overfitting" />).
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 (<Ref to="thm-mle-gaussian" />, <Ref to="ex-cross-entropy" />).

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](/en/computer-science/math-for-ml/linear-regression) |
| Compressing and diagnosing representations | linear algebra (eigenvalues, symmetric matrices) | the condition number $\kappa$ decides the iteration count | [Principal component analysis](/computer-science/math-for-ml/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](/computer-science/math-for-ml/gradient-descent) |
| Computing gradients in deep models | calculus (chain rule) | the derivative $p-y$ through the sigmoid | [Neural networks and backpropagation](/computer-science/math-for-ml/backpropagation) |
| Designing the loss | probability and statistics (likelihood) | squared loss from Gaussian noise, cross-entropy from Bernoulli | [Logistic regression](/en/computer-science/math-for-ml/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](/computer-science/math-for-ml/bayesian-statistics) |

The prerequisites are, for linear algebra, [Vector spaces and linear transformations](/en/mathematics/linear-algebra/vector-spaces), [Inner product spaces and Gram–Schmidt orthogonalisation](/mathematics/linear-algebra/inner-product-spaces) and [Eigenvalues and eigenvectors](/mathematics/linear-algebra/eigenvalues); for calculus, [Differentiation of multivariable functions and partial derivatives](/mathematics/calculus/multivariable-differentiation) and [The mean value theorem and Taylor's theorem](/en/mathematics/calculus/mean-value-and-taylor); and for probability, [Probability spaces and Kolmogorov's axioms](/en/mathematics/probability/probability-spaces) and [Random variables and expectation](/en/mathematics/probability/random-variables). 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 id="exr-normal-equation-4points" difficulty="Easy">

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 <Ref to="cor-unique-solution" />. 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 <Ref to="thm-normal-equation" />.

</Solution>
</Exercise>

<Exercise id="exr-learning-rate" difficulty="Standard">

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. <Ref to="thm-gd-quadratic" /> 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.

</Solution>
</Exercise>

<Exercise id="exr-laplace-loss" difficulty="Standard">

In the same setting as <Ref to="thm-mle-gaussian" />, 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 <Ref to="def-mle" /> 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.

</Solution>
</Exercise>

## 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](https://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](https://www.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/](https://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.