Linear Learning

At the most fundamental, learning is the search for structure from noise. Given observations scattered in high dimensional space, mathematics must extract a structure that we can use to predict where in the space new 0.3

To understand learning, we begin with a fitting problem: for observations $(x_1, y_1), \ldots, (x_n, y_n)$, we seek a function $f$ such that $f(x_i) \approx y_i$. For linear models, this means finding parameters $\beta$ such that $X\beta \approx y$, for a matrix of features $X$ and a vector of observed responses $y$. In the space of our observations, we want to find the linear surface $X\beta$ that best fits the data, such that this linear surface will also predict new observations. We want to solve a system of linear equations. If the data is all signal and no noise, we will be able to find parameters $\beta$ such that $X\beta= y$. In other words, we will be able to find a linear law that fits the data precisely. This is called linear regression.

With real data, however, we typically have many more observations than parameters, and each observation carries some error. For instance, if we measured the extension of a spring under different applied forces, we would observe a linear trend of extension with force, but not a perfect linear law. There is no single line that would pass through all the points. Hooke’s law $F = kx$ predicts a perfectly linear relationship, in practice we see points scattered around a theoretical line due to shortcomings in experiment - it is impossible to conduct a perfect experiment. In cases such as these, we get an overdetermined system with no solution to $X\beta = y$. The vector of observations $y$ does not lie in the column space of matrix $X$, and noise prevents an exact learned law. It is impossible to find any linear surface that passes through every point.

Despite the imperfect data, however, Hooke was able to find his law $F = kx$. How, we must ask, did he do this? The answer is that he didn’t spot a perfect relationship, but spotted a linear trend. He implicitly fit a learned line of best fit to his data, finding the relationship that best explained and could best predict his noisy observations. This is the crux of linear regression. No computer tasked with solving $X\beta = y$ for Hooke’s data could have found Hooke’s law. We must change our question. Instead of seeking a perfect solution to $X\beta=y$ we should seek the best approximate solution to $X\beta=y$. We can find this using the principle of least squares, whereby we minimize $\lVert X\beta - y \rVert^{2}$. The algebra leads into geometry - this is the solution found by dropping a perpendicular from $y$ onto the space of possible predictions, finding the closest achievable target. Learning becomes the geometry of projections in high-dimensional spaces, where we find the orthogonal projection of $y$ onto the column space of $X$.

The Linear Equation

We begin at the most fundamental level. Consider a simple system of linear equations: $$\begin{align*} x_1 + x_2 + x_3 &= 2 \\ 2x_1 - x_2 + x_3 &= -1 \\ -x_1 + 2x_2 + 2x_3 &= 1 \end{align*}$$

Instead of repetitively writing $x_1$s, $x_2$s and $x_3$s in every line, we can more compactly write these equations in matrix form. We use a matrix $A \in \mathbb{R}^{3 \times 3}$ of coefficients and a vector $x \in \mathbb{R}^3$ of unknowns and $b \in \mathbb{R}^3$ of solutions. Then we can write: $$A = \begin{pmatrix} 1 & 1 & 1 \\ 2 & -1 & 1 \\ -1 & 2 & 2 \end{pmatrix}, \quad x = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}, \quad b = \begin{pmatrix} 2 \\ -1 \\ 1 \end{pmatrix}$$

Computing the matrix-vector product $Ax$ recovers exactly our original system of equations, with right hand sides $b$. The core problem of solving the system of equations, therefore, is to find a vector $x$ such that $$Ax = b.$$

The matrix $A$ is a linear map $A:\mathbb{R}^n\to\mathbb{R}^m$. $A$ transforms all of $\mathbb{R}^n$, sending every vector $v$ from that space to $Av$ in $\mathbb{R}^m$. Every row of the matrix $A$ (each equation in our system of equations) applies a linear operation onto the components of $x$ (such as $x_1 + x_2 + x_3$) to yield the components of the transformed vector. In our problem of linear equations we are given the matrix $A$ of coefficients and the vector $b$ of right hand sides. Our task then is to follow the transformation: to find the vector $x \in \mathbb{R}^n$ that is transformed by the linear map $A$ to give $b \in \mathbb{R}^m$.

The product $Ax$ reveals two different ways of analysing this problem. Firstly, we can see the product as a linear combination of the columns of $A$ with weights $x_1,x_2,x_3$, such that: $$\underbrace{\begin{pmatrix} 1 & 1 & 1\\ 2 & -1 & 1\\ -1 & 2 & 2 \end{pmatrix}}_{A} \underbrace{\begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix}}_{x} = x_1 \begin{pmatrix} 1\\ 2\\ -1 \end{pmatrix} + x_2 \begin{pmatrix} 1\\ -1\\ 2 \end{pmatrix} + x_3 \begin{pmatrix} 1\\ 1\\ 2 \end{pmatrix} = \underbrace{\begin{pmatrix} 2\\ -1\\ 1 \end{pmatrix}}_{b}.$$ When considering this system, the question we ask becomes: what combination of the columns of $A$ (the vectors $A_1$, $A_2$, $A_3$) yield the vector $b$? Our problem becomes finding a way to construct the vector $b$ from these columns. The columns of $A$ are the crux of its transformation, showing us where the standard basis vectors $e_1, e_2,e_3$ are transformed to by $A$.

These columns span a subspace consisting of all their linear combinations, called the column space of $A$. This column space is the output space of the transformation - every vector that is transformed by $A$ ends up in its column space. When we solve $Ax = b$, we find $x$, the coefficients telling us how to construct $b$ in terms of the basis of this column space.

The second way of analysing the problem is to look at the multiplication $Ax$ in terms of the rows of matrix $A$. Each entry of $Ax$ is the dot product of one row of $A$ with $x$, such that: $$Ax= \begin{pmatrix} (1,\,1,\,1)\cdot(x_1,\,x_2,\,x_3)\\[2pt] (2,\,-1,\,1)\cdot(x_1,\,x_2,\,x_3)\\[2pt] (-1,\,2,\,2)\cdot(x_1,\,x_2,\,x_3) \end{pmatrix} = \begin{pmatrix} x_1+x_2+x_3\\ 2x_1-x_2+x_3\\ -x_1+2x_2+2x_3 \end{pmatrix} = \begin{pmatrix} 2\\ -1\\ 1 \end{pmatrix}.$$

In this picture, every row of the matrix $A$ defines some constraint on the vector $x$. Consider, as an example, the first dot product

$$\begin{pmatrix}1&1&1\end{pmatrix} \cdot \begin{pmatrix}x_1\\ x_2\\ x_3\end{pmatrix} = x_1 + x_2 + x_3 = 2.$$

This dot product establishes a constraint on $x$, requiring that its elements must sum to two. Any $x$ whose elements do not sum to two will not be allowed as a solution to the system. We can see each row of $A$ as a linear functional of $x$, taking an inner product with $x$ and acting as a linear map from $\mathbb{R}^3$ to the scalar field $\mathbb{R}$. Each dot product is a new constraint on $x$, requiring that the functional achieve a specific value. While the "column picture" looks at building the vector $b$ from components, this "row picture" tries to find the only vector $x$ that satisfies certain constraints.

Each constraint defines an $n-1$-dimensional subspace of $\mathbb{R}^n$, consisting of all vectors that satisfy the constraint. If we have $n$ independent rows in the matrix $X$, we will have $n$ such subspaces. Our solution $x$ must lie at the intersection of all these subspaces, where all constraints are simultaneously satisfied. The diagram below shows two lines, each line is the set of points that satisfies one of the constraints. The solution to the system lies in their intersection.

The column space $\mathrm{Col}(A)$, is the output space of the transformation. A solution exists if and only if the vector $b$ is inside this column space. Then it can be expressed as a part of the output space, in terms of its bases (the independent columns of $A$). For uniqueness of a solution $x$, there must be $n$ independent columns, such that $b$ has a unique representation in terms of those bases.

In the row picture, for there to be a unique solution, the intersection between the $n$ hyperplanes must be a single point. If the intersection is of higher dimension, there are an infinite number of solutions, since an infinite number of vectors satisfy all of the conditions. The row picture requires $n$ independent rows.

Linear Regression

Introduction to Least Squares

Consider the learning problem described earlier. We observe data points $(x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n)$ (where each $x_n$ is an input value and each $y_n$ is the observed output), and wish to find the line $y = \beta_0 + \beta_1 x$ that fits these points best. We can express this as a series of linear equations: $$\begin{align*} y_1 &= \beta_0 + \beta_1 x_1 \\ y_2 &= \beta_0 + \beta_1 x_2 \\ &\vdots \\ y_n &= \beta_0 + \beta_1 x_n \end{align*}$$

We must remember, however, that we will almost never find the exact values of the unknowns $\beta_0$ and $\beta_1$ that satisfy this set of linear equations. In this example we have $n$ equations and only two unknowns. As a consequence the column space of the coefficient matrix $X$, $\operatorname{Col}(X)$, lives in $\mathbb{R}^n$ but is at most two dimensional. Our matrix $X$ is very long but not wide.

In general, we will have many more observations than unknowns, since collecting data is always easier than adding new observation parameters. Our observed data vector $y = [y_1, y_2, \dots, y_n]^T$ lives in the full $n$-dimensional space $\mathbb{R}^n$ whereas $\operatorname{Col}(X)$ takes up only a very small portion of this space. Because $\mathbb{R}^n$ is so vast in comparison to $\operatorname{Col}(X)$, even small experimental errors are sufficient to prevent the data vector $y$ from falling into this narrow column space, meaning there is no exact solution $\beta$. A better formulation of the linear equations would be to write that we want to find a vector $\beta = [\beta_1, \beta_2]^T$ such that:

$$\begin{align*} y_1 &\approx \beta_0 + \beta_1 x_1 \\ y_2 &\approx \beta_0 + \beta_1 x_2 \\ &\vdots \\ y_n &\approx \beta_0 + \beta_1 x_n \end{align*}$$

and $\beta$ is the best possible approximation. Just as with we did with earlier linear equations, we can represent our input data in a coefficient matrix $X \in \mathbb{R}^{n \times 2}$, and vectors $b \in \mathbb{R}^2$ and $y \in \mathbb{R}^2$. We consider then the product: $$X = \begin{pmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{pmatrix}, \qquad \beta = \begin{pmatrix} \beta_0\\ \beta_1 \end{pmatrix}$$

$$X\beta = \begin{pmatrix} \beta_0 + \beta_1 x_1\\ \beta_0 + \beta_1 x_2\\ \vdots\\ \beta_0 + \beta_1 x_n \end{pmatrix}$$

Since we want to find the best possible approximate solution, we need to formulate what it means to be the "best possible". To do this, we define a difference vector, often called the error $e$. This is the difference between the observed result $y$ and our predicted result $X\beta$. We write this as: $$y - X\beta = \begin{pmatrix} y_1 - (\beta_0 + \beta_1 x_1)\\ y_2 - (\beta_0 + \beta_1 x_2)\\ \vdots\\ y_n - (\beta_0 + \beta_1 x_n) \end{pmatrix}$$

Due to measurement error and the dimensions of the spaces, there will be no single line that passes perfectly through all points. Gaussian elimination, which relies on row operations to find a unique solution vector, will fail.

Consider four example measurements. They are each very close to a linear relationship, but each has small error such that it deviates slightly from the theoretical line: $$(0,\,2.10),\quad (1,\,2.60),\quad (2,\,3.20),\quad (3,\,3.70).$$ We can eagerly try to find a model $y=\beta_0+\beta_1 x$ by solving the linear system $X\beta=y$ exactly, where $$X=\begin{bmatrix} 1 & 0\\ 1 & 1\\ 1 & 2\\ 1 & 3 \end{bmatrix},\qquad \beta=\begin{bmatrix}\beta_0\\ \beta_1\end{bmatrix},\qquad y=\begin{bmatrix}2.10\\ 2.60\\ 3.20\\ 3.70\end{bmatrix}.$$ To do this we form an augmented matrix $[X y]$ and apply Gaussian elimination: $$\left[\!\begin{array}{cc|c} 1 & 0 & 2.10\\ 1 & 1 & 2.60\\ 1 & 2 & 3.20\\ 1 & 3 & 3.70 \end{array}\!\right] \;\xrightarrow{R_2\rightarrow R_2-R_1,\; R_3\rightarrow R_3-R_1,\; R_4\rightarrow R_4-R_1}\; \left[\!\begin{array}{cc|c} 1 & 0 & 2.10\\ 0 & 1 & 0.50\\ 0 & 2 & 1.10\\ 0 & 3 & 1.60 \end{array}\!\right]$$ $$\xrightarrow{R_3\rightarrow R_3-2R_2,\; R_4\rightarrow R_4-3R_2}\; \left[\!\begin{array}{cc|c} 1 & 0 & 2.10\\ 0 & 1 & 0.50\\ 0 & 0 & 0.10\\ 0 & 0 & 0.10 \end{array}\!\right].$$ The bottom rows assert that $0\cdot\beta_0+0\cdot\beta_1=0.10$, or that $0=0.10$, which is impossible. Gaussian elimination proves that the system is inconsistent and that the vector $y$ does not lie in the two-dimensional column space, $\operatorname{Col}(X)\subset\mathbb{R}^4$. There is no exact solution. The graphical representation of these data points makes this clear. When we plot the theoretical line of best fit (which this chapter will work out how to derive), we see that no point lies precisely on this line - there is no straight line passes through all of these observations.

The solution to this problem lies in minimizing the error between what we can achieve ($X\beta$, which must lie in the column space of $X$) and what we desire ($y$ which lies outside it). Since we cannot make $X\beta = y$ exactly, we seek $\beta$ that makes $X\beta$ as close as possible to $y$.

The principle of least squares, introduced independently by Carl Friedrich Gauss who applied it to astronomical problems in 1795, and by Adrien-Marie Legendre in 1805, suggests minimizing the sum of squared errors. Then we can write: $$S(\beta_0, \beta_1) = \sum_{i=1}^{n} \bigl(y_i - \beta_0 - \beta_1 x_i\bigr)^2 .$$ For each observation $(x_i, y_i)$ the term $\beta_0 + \beta_1x_i$ is our prediction, being the $y$ value our model assigned to $x_i$. The residual $y_i - (\beta_0 + \beta_1x_i)$ is the error. Since the squared Euclidean norm of any vector $v=(v_1,v_2,\ldots,v_n)$ is $$\|v\|^2 = v_1^2 + v_2^2 + \cdots + v_n^2 = \sum_{i=1}^{n} v_i^{2}.$$ We can write that $$\sum_{i=1}^{n} \bigl[y_i - (\beta_0 + \beta_1 x_i)\bigr]^2 = \bigl\|\,y - X\beta\,\bigr\|^{2}.$$

The sum of squared errors represents an elegant link between geometry and algebra. It gives both the squared numerical error between our prediction and the actual value, and the squared distance between the point $y$ (our observations) and the point $X\beta$ (our predictions) in $n$-dimensional space. Minimizing $S(\beta_0,\beta_1)$ means finding the values of $\beta_0$ and $\beta_1$ that make these two points as close as possible. Fundamentally, this minimisation problem represents finding the shortest path from our observations to the plane of possible predictions.

Beyond this, the choice of squared error is useful to us because of computational convenience. Consider, briefly, the unsquared error $\lvert y_i - \beta_0 - \beta_1 x_i\rvert$. As visible in the figure below, this function has a sharp corner at the point where the prediction equals the observation and we have zero error. At this point, the derivative jumps discontinuously from $-1$ to $+1$, meaning the function is not differentiable at $x = 0$.

If we consider, as a contrast, the parabolic graph of squared error we see that it is differentiable everywhere.

Since the function is differentiable everywhere, we can take partial derivatives with respect to $\beta_0$ and $\beta_1$, such that: $$\frac{\partial}{\partial \beta_0} S(\beta_0,\beta_1) = -2\sum_{i=1}^{n}\bigl(y_i - \beta_0 - \beta_1 x_i\bigr)$$

$$\frac{\partial}{\partial \beta_1} S(\beta_0,\beta_1) = -2\sum_{i=1}^{n} x_i\bigl(y_i - \beta_0 - \beta_1 x_i\bigr)$$

These derivatives exist everywhere and are continuous. Consider again the graph of the squared error parabola shown in the figure above. The minimum value of the parabola must come at a turning point, since that is where the function turns from moving down to moving up, and a parabola has one turning point only. This is obvious visually from the graph above. At a turning point, the partial derivatives are zero, so we can learn more about squared error’s turning point by setting the derivatives to zero:

: $$\frac{\partial S}{\partial \beta_0} = -2\sum_{i=1}^{n} \bigl(y_i - \beta_0 - \beta_1 x_i\bigr) = 0, \qquad \frac{\partial S}{\partial \beta_1} = -2\sum_{i=1}^{n} x_i \bigl(y_i - \beta_0 - \beta_1 x_i\bigr) = 0$$ Dividing by $-2$ yields: $$\sum_{i=1}^{n} \bigl(y_i - \beta_0 - \beta_1 x_i\bigr) = 0, \qquad \sum_{i=1}^{n} x_i \bigl(y_i - \beta_0 - \beta_1 x_i\bigr) = 0$$ Equivalently: $$\sum_{i=1}^{n} (y_i - \hat{y}_i) = 0, \qquad \sum_{i=1}^{n} x_i (y_i - \hat{y}_i) = 0$$ if we let $\hat{y}_i = \beta_0 + \beta_1 x_i$ be our predictions. Beautifully, these equations derived through relatively simple algebraic steps provide insight into what the ’best fit’ solution means. The first equation shows us that the sum of residuals must equal zero. This means geometrically that the residual errors balance out, such that we never systematically over- or under-predict. The line of best fit that we predict doesn’t consistently sit above or below the data points. The second equation shows us something more profound, demonstrating algebraically that the residuals are orthogonal to our feature vector $x$. The summed product represents the inner product between $x_i$ and $(y_i - \hat{y_i})$ $$\langle x, (y - \hat{y})\rangle = 0$$

A zero dot product is the condition for orthogonality. Since our design matrix $X$ has columns $[1,1, \dots, 1]^T$ and $[x_1, x_2, \dots, x_n]^T$, and $X\beta = \hat{y}$ it follows directly that:

$$X^T(y -X\beta)=0,$$ meaning that each column of $X$ has a zero dot product with $(y - X\beta)$, or equivalently that $(y - X\beta)$ is orthogonal to the whole column space $\operatorname{Col}(X)$ .

More generally, this shows us the geometric essence of least squares. The core idea is to find $\hat{\beta}$ such that $X\hat{\beta}$ is the point in $\operatorname{Col}(X)$ closest to $y$. As visible in the diagram below, the orthogonal projection of a vector $v$ onto a subspace $S$ is the unique point $p$ in $S$ such that the vector $(v - p)$ is orthogonal to $S$.

We can turn an algebraic problem into a geometric one. Since we derived that $(y - X\beta) \perp \operatorname{Col}(X)$ and $X\beta \in \operatorname{Col}(X)$, we can say that the point $\hat{y}=X\hat{\beta}$ (where $\hat{\beta}$ is the best fitted estimate), is the orthogonal projection of $y$ onto the column space. Then our residual vector $e=y-X\hat{\beta}$ represents the shortest path (shown as $y - p$ in the diagram above) from our observed data $y$ to its projection $X\beta$ in the column space. This is the irreducable error that is orthogonal to the model space.

To minimize $\lVert e\rVert$ the residual must be perpendicular to the column space. The geometric principle is that the shortest path from a point to subspace is the perpendicular. We see this in $2$-D, where the shortest distance from a point to a line is found by dropping the perpendicular. Projecting orthogonally onto a subspace in $n$-dimensions is the learned approximation with minimal error.

Fascinatingly, from what began as an optimisation problem using calculus has revealed the fundamental geometric principle of least squares and linear regression as a whole, that the best approximate solution is always the geometric projection. This is the solution that minimises the length of the orthogonal $e$, giving us the smallest distance between our prediction and the observed data. This is what corresponds to the minimum point on the quadratic from earlier.

The Normal Equations

Let us consider one more time the coefficient matrix $X$ $$X = \begin{bmatrix} 1 & x_1\\ 1 & x_2\\ \vdots & \vdots\\ 1 & x_n \end{bmatrix}$$ and the response vector $y$. We seek $\beta = [\beta_0, \beta_1]^{\top}$ such that $X\beta$ approximates $y$. The smallest possible error $y - X\beta$ must be the error directly orthogonal to the column space of $X$, as we demonstrated in the last section, such that: $$X^T(y - X\beta) = 0$$ Expanding this equation: $$X^Ty-X^TX\beta = 0$$

This represents the fundamental goal of least squares established at the beginning of this chapter. This is no longer an equation seeking a precise solution $\beta$ to $X\beta = y$, rather it is an equation that finds the $\beta$ that minimises $\|y - X\beta\|^2$. Since these equations are fundamentally based around the geometric base of orthogonality via $X^T(y - X\beta) = 0$, we call them the normal equations.

The matrix $X^TX$ is always $2\times2$, since $X$ is $n \times 2$, and the vector $X^Ty$ is $2 \times1$. Then $X^Ty-X^TX\beta = 0$ is a system of $2$ equations and $2$ unknowns. In general, if we have $X$ is $m\times n$, while before we would have had $n$ equations and $m$ unknowns for $n \gg m$, we now have a system of $n$ equations and $n$ unknowns. This is a square system that can be solved using standard techniques such as Gaussian elimination.

To demonstrate this, we can show algebraically what the system of equations will look like. For our specific design matrix $X$, we can compute $X^TX$. (I realise that most readers will prefer to gloss over the next few pages, those readers are excused and should feel no guilt in doing so. Please ensure, however, that you notice the highlights on pages 13 and 15.)

$$\begin{align*} X^TX &= \begin{bmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_n \end{bmatrix} \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{bmatrix} \\[0.5em] &= \begin{bmatrix} n & \sum_{i=1}^n x_i \\[0.3em] \sum_{i=1}^n x_i & \sum_{i=1}^n x_i^2 \end{bmatrix} \end{align*}$$

We compute $X^Ty$: $$\begin{align*} X^Ty &= \begin{bmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_n \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} \\[0.5em] &= \begin{bmatrix} \sum_{i=1}^n y_i \\[0.3em] \sum_{i=1}^n x_i y_i \end{bmatrix} \end{align*}$$

The normal equations $(X^TX)\beta = X^Ty$ therefore become: $$\begin{bmatrix} n & \sum_{i=1}^n x_i \\[0.3em] \sum_{i=1}^n x_i & \sum_{i=1}^n x_i^2 \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix} = \begin{bmatrix} \sum_{i=1}^n y_i \\[0.3em] \sum_{i=1}^n x_i y_i \end{bmatrix}$$

Expanding yields the system: $$\begin{align*} n\beta_0 + \left(\sum_{i=1}^n x_i\right)\beta_1 &= \sum_{i=1}^n y_i \\[0.3em] \left(\sum_{i=1}^n x_i\right)\beta_0 + \left(\sum_{i=1}^n x_i^2\right)\beta_1 &= \sum_{i=1}^n x_i y_i \end{align*}$$

This is our final system of equations, which has two equations and two unknowns $\beta_0$ and $\beta_1$. Since we have come this far already, the writer cannot resist going the extra step to solve for $\beta$:

We solve this system of equations using Gaussian elimination. To eliminate $\beta_0$ from the second equation, we multiply the first equation by $\frac{\sum_{i=1}^n x_i}{n}$:

$$\left(\sum_{i=1}^n x_i\right)\beta_0 + \frac{\left(\sum_{i=1}^n x_i\right)^2}{n}\beta_1 = \frac{\left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right)}{n}$$

Subtracting this from the second equation: $$\begin{align*} &\left(\sum_{i=1}^n x_i^2\right)\beta_1 - \frac{\left(\sum_{i=1}^n x_i\right)^2}{n}\beta_1 \\ &\quad = \sum_{i=1}^n x_i y_i - \frac{\left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right)}{n} \end{align*}$$

Factorising: $$\left[\sum_{i=1}^n x_i^2 - \frac{\left(\sum_{i=1}^n x_i\right)^2}{n}\right]\beta_1 = \sum_{i=1}^n x_i y_i - \frac{\left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right)}{n}$$

Multiplying both sides by $n$: $$\beta_1 = \frac{n\sum_{i=1}^n x_i y_i - \left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right)}{n\sum_{i=1}^n x_i^2 - \left(\sum_{i=1}^n x_i\right)^2}$$

Applying back substitution, substituting this result into the first equation: $$\beta_0 = \frac{\sum_{i=1}^n y_i - \left(\sum_{i=1}^n x_i\right)\beta_1}{n}$$

It turns out that there is a more elegant way in which one can rewrite these formulae for $\beta_0$ and $\beta_1$ in terms of the deviations from the means $\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$ and $\bar{y} = \frac{1}{n}\sum_{i=1}^n y_i$.

Considering the numerator of $\beta_1$, we can show that: $$\begin{align*} n\sum_{i=1}^n x_i y_i - \left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right) &= n\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}) \end{align*}$$

We now expand the right side: $$\begin{align*} \sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}) &= \sum_{i=1}^n (x_i y_i - x_i \bar{y} - \bar{x} y_i + \bar{x}\bar{y}) \\ &= \sum_{i=1}^n x_i y_i - \bar{y}\sum_{i=1}^n x_i - \bar{x}\sum_{i=1}^n y_i + n\bar{x}\bar{y} \\ &= \sum_{i=1}^n x_i y_i - \bar{y}(n\bar{x}) - \bar{x}(n\bar{y}) + n\bar{x}\bar{y} \\ &= \sum_{i=1}^n x_i y_i - n\bar{x}\bar{y} \end{align*}$$

Since $\bar{x}\bar{y} = \frac{\sum_{i=1}^n x_i}{n} \cdot \frac{\sum_{i=1}^n y_i}{n} = \frac{(\sum_{i=1}^n x_i)(\sum_{i=1}^n y_i)}{n^2}$:

$$\begin{align*} n\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}) &= n\left[\sum_{i=1}^n x_i y_i - \frac{(\sum_{i=1}^n x_i)(\sum_{i=1}^n y_i)}{n^2}\right] \\ &= n\sum_{i=1}^n x_i y_i - \frac{(\sum_{i=1}^n x_i)(\sum_{i=1}^n y_i)}{n} \\ &= n\sum_{i=1}^n x_i y_i - \left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right) \end{align*}$$

The same applies to the denominator: $$\begin{align*} n\sum_{i=1}^n x_i^2 - \left(\sum_{i=1}^n x_i\right)^2 &= n\sum_{i=1}^n (x_i - \bar{x})^2 \end{align*}$$

Hence we write that: $$\begin{align*} \beta_1 &= \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2} \\[0.5em] \beta_0 &= \bar{y} - \beta_1 \bar{x} \end{align*}$$

A fundamental beauty of algebra is highlighted in the simplification from the complex sums of page 14 to the significantly simpler and more elegant sums above. Geometrically, the algebraic steps represented centering the data, eliminating the complicated cross-terms that absolute positioning brings. For instance, the complex numerator $$n\sum_{i=1}^n x_i y_i - \left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right)$$ transforms into the elegant centered form $$n\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y}),$$ where each data point is now measured relative to the sample means rather than the arbitrary origin. The complexity in the original equations arose because everything was measured relative to an origin that had no relevance to our data. Centering sets the origin of our coordinate system at the center of mass of our data, revealing the true geometric relationships.

The final formulae for $\beta_0$ and $\beta_1$ are profound because they reveal algebraically some fundamental insights about the nature of learning. This is because the formulae can be expressed in terms of variance and covariance. Variance is a statistical measure of how much a variable spreads out from its average. Since we measure deviation from the mean using $(x_i - \bar{x})^2$ (where the square eliminates the negatives) for some $x_i$, we can say that the variance of a variable $x$ is: $$\sum^n_i\frac{(x_i - \bar{x})^2}{n}$$

Covariance is a measure of of how two variables move together with respect to their means. If covariance is large and positive, then when $x$ is above its mean, $y$ tends to be above its mean too. And when $x$ is below its mean, $y$ tends to be below its mean too. They rise and fall together. We can write a formula for covariance, similarly as we did the formula for variance, as $$\frac{\sum_{i=1}^{n} \bigl(x_i - \bar{x}\bigr)\bigl(y_i - \bar{y}\bigr)} {n}$$

Now consider again our original equations for $\beta_0$ and $\beta_1$:

$$\begin{align*} \beta_1 &= \frac{\sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2} \\[0.5em] \beta_0 &= \bar{y} - \beta_1 \bar{x} \end{align*}$$

The denominator of our formula for $\beta_1$ is $n$ times the variance of $x$ and the numerator is $n$ times the covariance of $x$ and $y$! Our formula for $\beta_1$ is the ratio between the covariance of $x$ and $y$ and the variance of $x$. This means we can write: $$\beta_1 = \frac{\text{Cov}(x,y)}{\text{Var}(x)}$$

$\beta_1$ is the slope of the least squares solution. What we derived algebraically, that $\beta_1$ is covariance divided by variance gives us valuable information about the nature of the slope. Suppose our $x$ values were tightly clustered around their mean, with a low variance. Then even a small covariance between $x$ and $y$ would produce a steep slope, since we are dividing by a small variance. The $x$-$y$ slope measures how much $y$ changes with a small change in $x$. When $x$ has limited spread, even small changes in $y$ will result in large changes per unit $x$, since we are dividing by a large variance. If our $x$ values are widely scattered with high variance, we need a much larger covariance to produce the same slope.

This scaling builds a bridge between the algebra and geometry of least squares. Algebraically, we take the covariance variance ratio, but geometrically, this has the effect of ensuring that $\beta_1$ measures the predictive relationship between $x$ and $y$ and also accounts for the natural variability present in the input variable. The learned slope is calibrated to both the scale and spread of our data.

The Projection Matrix

To wrap up our study of linear regressions, we can now generalise our specific least squares solution to the full matrix solution. If we consider the normal equations: $$X^TX\beta = X^Ty$$ Since $(X^TX)^{-1}(X^TX) = I$, we can solve for the parameters $\beta$, such that: $$\beta = (X^TX)^{-1}X^Ty.$$

Matrices are so beautiful because of their versatility throughout mathematics. While they can be seen as just a simple table of numbers, they also represent transformations of space, systems of equations, the computational steps to solve the equations, geometric projections, and much more.

The matrix $(X^TX)^{-1}X^Ty$¹ accomplishes four things at once: it solves the normal equations $X^TX\beta = X^Ty$, minimises the sum of squared errors $\|y -X\beta\|^2$, computes the geometric projection of $y$ onto the column space of $X$, and finds the coefficients that extract the most information about the $x$-$y$ relationship. Here, geometry, algebra and statistics converge.

Let us conclude this section with the example from earlier. We had four measurements that were very close to a straight line, but each had a small experimental error: $$(0,\,2.10),\quad (1,\,2.60),\quad (2,\,3.20),\quad (3,\,3.70).$$ We wanted to fit the model $y=\beta_0+\beta_1 x$ by solving the linear system $X\beta=y$ exactly, but found that Gaussian elimination could not solve the system - $y$ was not in the column space of $X$. We now have the tools to find the best approximate solution, the line of best fit.

$$X^\top X=\begin{bmatrix}4&6\\6&14\end{bmatrix}, \qquad X^\top y=\begin{bmatrix}11.60\\20.10\end{bmatrix}.$$ To compute the inverse of $X^T X$, we find its determinant determinant $\det=4\cdot14-6\cdot6=56-36=20$, so $$(X^\top X)^{-1} =\frac{1}{20}\begin{bmatrix}14&-6\\-6&4\end{bmatrix}.$$ Therefore: $$\hat\beta = (X^\top X)^{-1}X^\top y =\frac{1}{20} \begin{bmatrix}14&-6\\-6&4\end{bmatrix} \begin{bmatrix}11.60\\20.10\end{bmatrix} =\frac{1}{20} \begin{bmatrix} 14(11.60)-6(20.10)\\ -6(11.60)+4(20.10) \end{bmatrix} =\frac{1}{20}\begin{bmatrix}41.8\\10.8\end{bmatrix} =\begin{bmatrix}2.09\\0.54\end{bmatrix}.$$ We have computed our linear regression, such that the line of best fit is therefore: $$\boxed{\ y\approx 2.09 + 0.54\,x\ }.$$

Deep Learning

Linear regression methods find the best linear approximation $f(x) = x^T \beta$ by solving $\beta = (X^\top X)^{-1}X^\top y$. This fits the best possible linear surface to our data, such that $x^T \beta$ computes the surface height of at point $x$. Linear regression has one fundamental shortcoming: it confines us to the space of linear functions.

Consider the following problem² We want a linear model that can map the following input tuples to their respective outputs: $(0,0) \to 0$, $(0,1) \to 1$, $(1,0) \to 1$, $(1,1) \to 0$. This problem is called the XOR problem. XOR stands for "exclusive or", where the "or" part means the output is true when one of the inputs is $1$, and the exclusive part means that this excludes the case when both inputs are $1$. This differs slightly from the problems we saw earlier with linear regression because it is a classification problem. While least squares methods give us the best linear prediction $f(x)=x^\top\beta$ for a continuous target, here we want a linear classifier that will return either $1$ or $0$, True or False. To reuse our linear machinery, we establish a score: $$s(x)=w^T x + b,$$ and then classify by its sign, predicting True if $s(x)\geq0$ and False if $s(x)<0$. This exposes the geometric perspective of the linear classifier, where the set $\{x:\,s(x)=0\}$ is a hyperplane that divides space into two half-spaces. For all input points on one side of the hyperplane our model returns False; for all points on the other side it returns True. ³

A fundamental result of the 1960s was the discovery that no linear function $h(x) = w_1 x_1 + w_2 x_2 + b$ can separate the classes $\{0, 1\}$ from $(0,0) \to 0$, $(0,1) \to 1$, $(1,0) \to 1$, $(1,1) \to 0$. For this to happen, a linear model with weights $w_1$ and $w_2$ would have to satisfy: $$\begin{aligned} w_1\cdot 0 + w_2\cdot 0 + b &< 0,&& b < 0\\ w_1\cdot 0 + w_2\cdot 1 + b &> 0,&& w_2 + b > 0\\ w_1\cdot 1 + w_2\cdot 0 + b &> 0,&& w_1 + b > 0\\ w_1\cdot 1 + w_2\cdot 1 + b &< 0,&& w_1 + w_2 + b < 0 \end{aligned}$$ If the output of the model were positive we would assign one, if it were negative we would assign zero. A contradiction emerges when we sum the second and third inequalities: $$(w_2+b)+(w_1+b) > 0, \ \ \ w_1+w_2+2b>0,$$ since $b<0$ we have $-b>0$. Subtracting $b$ from both sides of the last inequality gives: $$w_1+w_2+b > -b > 0.$$ Then $w_1+w_2+b>0$, which contradicts $w_1+w_2+b<0$. This contradiction demonstrates that no linear model can separate the two classes, that is to say that no hyperplane can implement the XOR function on these four points. Figure 3.1 below shows visually why no decision boundary can divide the zeroes from the ones in XOR.

This demonstrates the fundamental principle of deep learning: to learn anything that contains nonlinearity, we need nonlinearity in our learning algorithm. XOR is fundamentally a nonlinear learning problem, since there is no way to draw a straight line separating the classes. Instead of a single linear transformation, we must compose multiple affine maps, each with nonlinearities. If we use a three layer network to learn our function, where each subsequent layer gets the outputs from the previous one, we can write this as $$h(x) = W_3\sigma\big(W_2\sigma(W_1 x + b_1) + b_2\big) + b_3$$ where $\sigma$ is some nonlinear activation function.

A famous result in the world of Deep Learning is Cybenko’s universal approximation theorem demonstrates that a network with a single hidden layer of sufficient width can approximate any continuous function on a bounded set to arbitrary precision. More precisely, for any continuous $f: K \to \mathbb{R}$ where $K \subset \mathbb{R}^n$ is a bounded set, $h(x)$ is our network, and any $\varepsilon > 0$, there exist weights such that: $$\big|f(x) - h(x)\big| < \varepsilon.$$

Moving from linear to composed models with nonlinearities means that we inevitably loose uniqueness of solution and must resort to iterative optimization with no guarantee of finding the optimum. The loss surface is filled with various saddle points and local minima.

What deep learning gains, however, is far better than this. Neural networks achieve generalisation. While linear regression worked only for linear relationships, neural networks have universal function approximation such that they can, in principle, learn to represent any continuous function to arbitrary precision. making them extraordinarily flexible learning machines. This is the power of deep learning with neural networks. It is the universal learning algorithm that can learn anything.

The Neuron

In 1837, Czech anatomist Jan Purkyně was examining thin slices of neural tissue, taken from the cerebellum. Purkyně discovered large, flask-shaped cells that were distinct from anything previously documented. The cells had a strange structure, with tree-like branches that extended in all directions. They were named the "Purkinje cells" and were postulated to represent the most basic nerve cells. Beyond this, however, the nervous system remained a mystery to scientists.

"A fundamental breakthrough came in 1873 when Italian physician Camillo Golgi, working in a home laboratory, discovered by accident that a solution of potassium dichromate and silver nitrate could stain individual neurons black. For the first time the full structure of these cells was visible. These cells were unlike anything anyone had ever seen before. Under the microscope emerged intricate cells with central bodies from which delicate, tree-like dendrites branched out on one side while long, thread-like axons extended from the other.

Golgi believed that these neurons were fused together in a continuous chain. Further study proved him false, instead they were found to be individual units that communicated across tiny gaps by electrical potential. These gaps were later called "synapses." Individual neurons were found to maintain a tiny voltage across their membranes, around -70 millivolts. When enough input signals arrived at the dendrites, the voltage would rise to a critical threshold of -55 millivolts which would trigger a rapid release of sodium ions. This would cause the voltage to spike to +40 millivolts, causing an electrical pulse, called the action potential, to travel down the axon to signal other neurons. The neuron has "fired".

Although the biology of the nervous system is extraordinarily complex, neurons themselves are not too difficult to model mathematically. The neuron receives as input a series of voltages $[x_1, x_2, \dots x_n]^T$ via the synapses. Each synapse has an efficacy (a measure of how well it transmits the charge). We can model this by assigning each connection a "synaptic strength". Each neuron weights its inputs by multiplying by a weight vector $[w_1, w_2, \dots w_n]^T$, where $w_i$ represents the synaptic strength of the connection with the $i$th neuron. The input charges combine in the cell body, so in our artificial neuron the cell body is the sum of all $x_iw_i$. This is equivalent to taking the dot product between the vectors: $$(x_1 \ \ x_2 \ \ \dots \ \ x_n) \begin{pmatrix} w_1\\ w_2\\ w_3\\ \vdots\\ w_n \end{pmatrix} = \sum_iw_ix_i$$

Biologically, every neuron has a different baseline activation level, such that some neurons are naturally more or less ready to fire. We model this in an artificial neuron by adding a bias $b$ to our sum, such that we have: $$\sum_iw_ix_i +b.$$ With a large positive bias, the neuron fires with less input. Large negative bias raises the neuron’s natural inhibition and a greater input is needed to fire.

The neuron obeys a very binary firing pattern - either its threshold is reached and it fires, or its threshold isn’t reached, it doesn’t fire and the charge is dissipated. We model this through use of an activation function. This is a nonlinear function that determines the neuron’s final output to subsequent neurons based on the its aggregated input $w^Tx + b$. A neuron with weights $w$, bias $b$ and nonlinearity $\varphi$, given an $m$-dimensional vector $x$ as input produces as an output $${\displaystyle \varphi \left(\sum _{j=0}^{m}w_{kj}x_{j}\right)},$$

Figure 3.3 shows a schematic diagram of this process.

As we saw in the example of XOR, nonlinear activation functions are very important for neural networks. Linearity means the output for a sum of inputs equals the sum of the outputs. Most things don’t add - for instance images or patterns do not combine into new useful features by simple addition. They cannot be learned through linear methods alone.

Without nonlinearities, even the deepest networks would collapse to a single linear transformation. Consider a network with $L$ layers and no nonlinear activation function (such that each layer is just a matrix $W$): $$y = W^{(L)} W^{(L-1)} \cdots W^{(2)} W^{(1)}\, x = Wx.$$ The product of two linear transformations is itself linear. No matter how many we stack, we can still learn linear functions only, and the whole premise of deep learning pointless.

To model the complex nonlinearity found in neurons, we use simplified continuous functions that are enough to capture the nonlinearity but are not overly complex. One such function often used is the hyperbolic tangent function: $$\tanh(z)=\frac{e^{z}-e^{-z}}{e^{z}+e^{-z}}=\frac{e^{2z}-1}{e^{2z}+1}$$ $\tanh$ is continuous and differentiable everywhere. We can compute its derivative from the earlier expression: $$\frac{d}{dz}\tanh z = \frac{(2e^{2z})(e^{2z}+1)-(e^{2z}-1)(2e^{2z})}{(e^{2z}+1)^{2}} = \frac{4e^{2z}}{(e^{2z}+1)^{2}} = 1-\left(\frac{e^{2z}-1}{e^{2z}+1}\right)^{2} = 1-\tanh^{2} z.$$

Consider the graph of $\tanh$.

This graph highlights why $\tanh$ is often used to model the nonlinearities in neurons. Firstly, for $\lvert z\rvert > 2$, $\tanh(z)\approx\pm1$. The artificial neuron using $\tanh$ has a steep threshold such that it is almost always either "on" ($+1$) or "off" ($-1$), just like the binary nature of action potentials in neurons.

The graph of $\tanh$ also has a near linear region for small $z$. Near $z=0$, $\tanh(z)\approx z$, the derivative $\tanh'(0)=1$. This has the convenient effect that small signals pass through unchanged. This is valuable to us because it means the neuron doesn’t muffle small signals which are often very important for later layers. At its core, $\tanh$ performs a smooth squashing of space such that an infinite input range is compressed to finite outputs, but order is preserved, if $a < b$, then $\tanh(a) < \tanh (b)$.

In 2010, a new activation function, called the Rectified Linear Unit or ReLU for short, became popular. ReLU is defined as: $$\operatorname{ReLU}(z)=\max(0,z)= \begin{cases} z & \text{if } z>0,\\[4pt] 0 & \text{if } z\le 0. \end{cases}$$

ReLU is a very simple function, outputting the input value for positive numbers and zero for negative numbers. Despite its simplicity, ReLU has several remarkable properties that make it highly effective for neural networks. For any given input, many ReLU neurons will output exactly zero, meaning only a subset of neurons remain active. This creates a sparse activation pattern that helps the network focus on the most relevant features for each input. When a ReLU neuron outputs zero, it effectively "turns off" that pathway, allowing the network to ignore irrelevant information and concentrate computational resources on useful signals. This selective activation mirrors biological neural networks, where only a small percentage of neurons (approximately 1-4%) are active at any given moment, creating efficient and focused information processing.

As we can see from its graph, ReLU has a simple derivative:

$$\operatorname{ReLU}'(z)= \begin{cases} 1 & \text{if } z>0,\\[4pt] 0 & \text{if } z\le 0. \end{cases}$$

This means that for active neurons ($z>0$), gradients flow through unchanged during backpropogation. Computing this derivative is also very computationally efficient. ReLU helps to mitigate the gradient vanishing problem that occurs when training deep networks, when during backpropagation through many layers gradients can become progressively smaller due to multiplication of small derivatives. Since ReLU’s derivative is one for active neurons, gradients can flow through these neurons without being dampened allowing for the training of deeper networks.

The Neural Network

A single artificial neuron computes: $$h = \sigma(w^T x + b)$$ where: $x\in\mathbb{R}^n$ is the input vector, $w\in\mathbb{R}^n$ is a weight vector of learnable parameters, $b\in\mathbb{R}$ is a learnable bias term and $\sigma$ is an activation function. Then $h\in\mathbb{R}$ is the output.

A layer of neurons is a parallel collection of neurons that operate in parallel. Each neuron in the layer receives the same input vector $x$ from the previous layer, where the $i$th component of $x$ is the output from the $i$th neuron in the previous layer. Each neuron then passes its output to every neuron in the next layer. The strength of each connection is determined by the weights in the receiving neurons’ weight vectors. If we have $m$ neurons in one layer, we can write the entire layer’s computation concisely with matrix notation: $$h = \sigma(Wx + b).$$

If there are $m$ neurons in the current layer and $n$ neurons in the previous layer, then $W\in \mathbb{R}^{m\times n}$ is the layer’s weight matrix, where each row in $W$ represents the weights for one neuron. The element $W_{ij}$ represents the connection strength between the $i$th neuron in the layer the $j$th neuron in the previous layer. $b \in \mathbb{R}^m$is the bias vector. $\sigma$ is the nonlinearity. Note that $\sigma$ is applied element-wise, such that $\sigma(v)=[\sigma(v_1),\ldots,\sigma(v_m)]^\top$. Finally $\mathbf{h}\in\mathbb{R}^m$ is the layer’s output vector, and is passed on as input to the subsequent layer.

Each layer performs two operations. Firstly, it takes a weighted combination of the inputs using the synaptic strengths and adds the bias term. It computes an affine transformation, to yield $z = Wx + b$. Secondly, it applies the nonlinear activation function element-wise to produce the output $h = \sigma(z)$. .

The linear part $Wx + b$ applies some transformation on the input space. The nonlinear part ensures that every neuron adds something new to the network, preventing the network from collapsing to a single linear transformation. The standard notation for neural nets is to denote the $l$-th layer’s weights as $W^{(l)}$ and its biases as $b^{(l)}$. Then we can write a neural network as a composed function $$f(x) = (f_L\circ f_{L-1}\circ\cdots\circ f_2\circ f_1)(x),$$ where $f_\ell(x)=\sigma_\ell(W^{(\ell)}x+b^{(\ell)})$.

We can use neural nets to correctly solve the XOR problem. As inputs we provide:

$$x=\begin{bmatrix}x_1\\[2pt]x_2\end{bmatrix}\in \left\{ \begin{bmatrix}0\\[2pt]0\end{bmatrix}, \begin{bmatrix}0\\[2pt]1\end{bmatrix}, \begin{bmatrix}1\\[2pt]0\end{bmatrix}, \begin{bmatrix}1\\[2pt]1\end{bmatrix} \right\}, \qquad \operatorname{XOR}(x_1,x_2)\in\{0,1,1,0\}.$$

If $W^{(n)}$ represents the weight matrix of the $n$th layer, and $b^n$ represents the bias of the $n$th layer, we describe our network $$h=\operatorname{ReLU}\!\big(W^{(1)}x+b^{(1)}\big),\qquad \hat y=W^{(2)}h+b^{(2)},$$

We set our weights to: $$W^{(1)}=\begin{bmatrix}1&1\\[2pt]1&1\end{bmatrix},\quad b^{(1)}=\begin{bmatrix}0\\[2pt]-1\end{bmatrix},\quad W^{(2)}=\begin{bmatrix}1&-2\end{bmatrix},\quad b^{(2)}=0.$$

The first layer’s weight matrix $W^{(1)}$ has identical rows $[1 \ \ 1]$. This means that both hidden neurons in the network compute the same linear combination $x_1$ + $x_2$. Their different biases create distinct detection mechanisms. The first neuron computes $h_1 = \operatorname{ReLU}(x_1 + x_2)$, which activates whenever one neuron is positive, and the second neuron computes $h_2 = \operatorname{ReLU}(x_1 + x_2 -1)$, which only activates when both neurons are positive (since $x_1 + x_2 \geq 1$ is required for positive argument to $\operatorname{ReLU}$).

The second layer the detectors from the previous layer with weights $W^{(2)} = [1, \ -2]$ to compute $\hat{y} = h_1 - 2h_2$ This implements the logical expression “at least one input active” minus “both inputs active”. The coefficient $-2$ (rather than $-1$) is crucial because when both inputs are active $$h_1=\operatorname{ReLU}(2)=2$$ rather than $1$, so we need $$\hat{y}=2-2\cdot 1=0$$ to correctly output zero for the $(1,1)$ case. For the cases $(1,0)$ and $(0,1)$ we have $$h_1=1,\qquad h_2=0,$$ yielding $$\hat{y}=1-2\cdot 0=1.$$ For $(0,0)$ both hidden units output zero, giving $\hat{y}=0$. Thus the network converts the linearly inseparable XOR pattern into a linearly separable one in the two-dimensional hidden space, where the final layer needs only a simple linear combination to distinguish the desired outputs.

Let $s=x_1+x_2$. Then $$h=\begin{bmatrix}\operatorname{ReLU}(s)\\[2pt]\operatorname{ReLU}(s-1)\end{bmatrix}, \qquad \hat y=\operatorname{ReLU}(s)-2\,\operatorname{ReLU}(s-1).$$

We verify on the set $\{0,1\}^2$: $$\begin{array}{c|c|c|c} (x_1,x_2) & s=x_1+x_2 & h=\big[\operatorname{ReLU}(s),~\operatorname{ReLU}(s-1)\big]^{\mathsf T} & \hat y \\ \hline (0,0) & 0 & [0,0]^{\mathsf T} & 0\\ (1,0) & 1 & [1,0]^{\mathsf T} & 1\\ (0,1) & 1 & [1,0]^{\mathsf T} & 1\\ (1,1) & 2 & [2,1]^{\mathsf T} & 0 \end{array}$$

Backpropagation

Modern artificial intelligence models, which are all based on neural networks have huge amounts of parameters. The parameters of a neural network are its learnable parts, consisting of its weights in the matrices $W$ and its biases $b$. GPT-5 has approximately 600 billion parameters. As of September 2025, Meta’s Llama-4 "Behemoth" has the greatest number of parameters of any model, with a total of 2 billion parameters. Each parameter is a continuous value that is typically represented as a 32-bit floating-point number. If we consider GPT-5 again, the space of possible parameter configurations is therefore $\mathbb{R}^{600000000000}$.

We could not possibly search this space for the best parameter configuration for our learning problem, since we have finite time and finite compute. GPT-5 was trained over a period of a few months using approximately $10^{26}$ FLOPS.

To do this, we need a good algorithm to train our model, updating the parameters (weights of the neural network) until we find the best parameters solve our problem or learn from our data. The algorithm we use to do this is called backpropagation, a remarkably simple application of the chain rule.

To measure how far away our prediction is from our expected result, we define a loss function $L(\theta)$. To update the parameters in our favour, we need a way that is better than brute-force to measure how much changing a specific parameter changes the final output. That is to say for a neural network with parameters $\theta=(\theta_1,\theta_2,\dots,\theta_n)$ and loss function $L(\theta)$, we need the gradient: $$\nabla L = \begin{bmatrix} \frac{\partial L}{\partial \theta_1} \\[4pt] \frac{\partial L}{\partial \theta_2} \\[4pt] \vdots \\[4pt] \frac{\partial L}{\partial \theta_n} \end{bmatrix}.$$ Each component $\partial L/\partial\theta_i$ represents the change in the loss $\varepsilon\cdot\partial L/\partial\theta_i$ if parameter $\theta_i$ is incremented by $\varepsilon$.