Machine learning

The question

You measured a rate constant at six temperatures and want the activation energy. You computed the adsorption energy of a molecule on twenty surfaces and want to predict it for a twenty-first. In both cases you have pairs $(x_i, y_i)$ and want a rule $f$ that turns a new $x$ into a good guess of $y$. That is supervised learning, and almost everything in this course is a version of it: the rule gets more flexible (a line, a neural network, an interatomic potential), but the recipe stays the same.

The recipe has three ingredients:

1. a model $f(x;\boldsymbol\theta)$, a family of functions indexed by parameters $\boldsymbol\theta$;
2. a loss $\ell(y, \hat y)$ that says how bad a prediction $\hat y$ is when the truth is $y$;
3. an optimiser that searches for the parameters with the smallest average loss on the data.

Minimising the average loss over a data set $\mathcal D = \{(x_i,y_i)\}_{i=1}^N$ is called empirical risk minimisation:

The word empirical matters. What we really want is a small error on data we have not seen, the expected loss over the whole population. We can only measure it on the sample we have, which is why the lesson on generalisation exists.

The simplest model: a straight line

Take $f(x; w, b) = wx + b$ and the squared loss $\ell = (y - \hat y)^2$. The objective is the sum of squared residuals

$L$ is a smooth bowl in $(w,b)$, so the minimum is where both partial derivatives vanish:

The first equation says the residuals sum to zero, so the line passes through the centroid $(\bar x, \bar y)$. Substituting $b = \bar y - w\bar x$ into the second gives the familiar result

Try three things in the lab. First, fit the line yourself with the sliders and notice how hard it is to beat least squares; no other line has a smaller mean squared error, by construction. Second, press Add an outlier and watch the least-squares line swing towards it. Third, switch to least absolute deviation with the outlier still there. The squared loss punishes a residual of 4 sixteen times more than a residual of 1, so a single bad point can dominate; the absolute loss punishes it only four times more.

The matrix form, which every later model reuses

Stack the inputs into a design matrix $\mathbf X$ whose rows are feature vectors, here $(1, x_i)$, and the targets into $\mathbf y$. Then $\hat{\mathbf y} = \mathbf X\boldsymbol\theta$ and

Setting the gradient $\nabla_{\boldsymbol\theta} L = -2\mathbf X^{\top}(\mathbf y - \mathbf X\boldsymbol\theta)$ to zero gives the normal equations

Geometrically, $\mathbf X\hat{\boldsymbol\theta}$ is the orthogonal projection of $\mathbf y$ onto the column space of $\mathbf X$: the residual vector is perpendicular to every feature column. Numerically, never form $(\mathbf X^\top\mathbf X)^{-1}$ explicitly; squaring $\mathbf X$ squares its condition number. Libraries solve the problem with a QR or singular-value decomposition instead (numpy.linalg.lstsq uses the SVD).

Where least squares comes from: maximum likelihood

Why squares and not absolute values or fourth powers? Suppose each measurement is the true line plus independent Gaussian noise of width $\sigma$:

The probability density of the whole data set, viewed as a function of the parameters, is the likelihood

Taking the negative logarithm,

and the only part that depends on $w$ and $b$ is the sum of squares. Least squares is maximum likelihood under Gaussian noise. Every other loss corresponds to a different noise model: the absolute loss is maximum likelihood for Laplace (double-exponential) noise, whose heavier tails make outliers less surprising, which is why that fit ignored the outlier.

The same view tells you what to do when the error bars differ from point to point. If point $i$ has standard deviation $\sigma_i$, the likelihood gives weighted least squares, $\sum_i (y_i-\hat y_i)^2/\sigma_i^2$: precise points count more. Fitting $\ln k$ instead of $k$ silently changes the weights, because a constant relative error in $k$ is a constant absolute error in $\ln k$. That is usually what you want for rate constants, and it is worth knowing you chose it.

How good is the fit?

Three numbers are worth reporting.

• The root-mean-square error, $\text{RMSE} = \sqrt{L/N}$, in the units of $y$. Compare it with the measurement uncertainty: a fit far better than your error bars is suspicious.
• The coefficient of determination, $R^2 = 1 - \sum_i (y_i - \hat y_i)^2 / \sum_i (y_i - \bar y)^2$, the fraction of the variance of $y$ the model explains. $R^2$ says nothing about whether the model is right: a straight line through curved data can still give $R^2 = 0.98$. Always look at the residuals; they should look like noise, with no trend.
• Uncertainties on the parameters. Under the Gaussian model the covariance of the estimate is $\operatorname{cov}(\hat{\boldsymbol\theta}) = \hat\sigma^2(\mathbf X^\top\mathbf X)^{-1}$ with $\hat\sigma^2 = L/(N-p)$ for $p$ parameters. An activation energy without an error bar is not a result.

In practice

The calculation in NumPy, here for an Arrhenius plot (run it locally, or in the Python lesson):

import numpy as np

T = np.array([300., 310., 320., 330., 340., 350.])           # K
k = np.array([1.2e-3, 2.4e-3, 4.6e-3, 8.5e-3, 1.5e-2, 2.6e-2])  # 1/s
R = 8.314462618                                              # J/(mol K)

X = np.column_stack([np.ones_like(T), 1.0 / T])   # features (1, 1/T)
y = np.log(k)
theta, res, rank, sv = np.linalg.lstsq(X, y, rcond=None)
lnA, slope = theta

N, p = X.shape
s2 = res[0] / (N - p)                         # residual variance
cov = s2 * np.linalg.inv(X.T @ X)             # acceptable for a tiny 2×2 problem
Ea = -slope * R / 1000                        # kJ/mol
dEa = np.sqrt(cov[1, 1]) * R / 1000
print(f"Ea = {Ea:.1f} ± {dEa:.1f} kJ/mol,  A = {np.exp(lnA):.2e} 1/s")

This prints Ea = 53.7 ± 0.1 kJ/mol, A = 2.65e+06 1/s. For nonlinear models, where no normal equations exist, scipy.optimize.curve_fit does the same job iteratively; the next lesson explains how such iterative optimisers work.

References and further reading

Why we need an iterative method

The straight line in the previous lesson had a closed-form answer because its loss is quadratic in the parameters. Most models do not: a nonlinear curve fit, a neural network or an interatomic potential has a loss $L(\boldsymbol\theta)$ with no formula for its minimum. What we can compute cheaply is the gradient $\nabla L(\boldsymbol\theta)$, the vector of partial derivatives that points in the direction of steepest increase. Gradient descent takes repeated small steps the other way:

where $\eta > 0$ is the learning rate (step size). The idea goes back to Cauchy in 1847, and with the improvements below it trains essentially every neural network in use today.

Why this direction? A first-order Taylor expansion around $\boldsymbol\theta$ gives $L(\boldsymbol\theta + \boldsymbol\delta) \approx L(\boldsymbol\theta) + \nabla L^{\top}\boldsymbol\delta$. Among all steps of a fixed length, the one that lowers $L$ most is the one anti-parallel to $\nabla L$ (by the Cauchy–Schwarz inequality). The expansion is only trustworthy for small $\boldsymbol\delta$, which is why the step size matters so much.

How large a step? The quadratic model

Near a minimum $\boldsymbol\theta^*$ every smooth loss looks like a quadratic bowl,

where $\mathbf H$ is the Hessian, the matrix of second derivatives. Write the error $\mathbf e_t = \boldsymbol\theta_t - \boldsymbol\theta^*$. Since $\nabla L = \mathbf H\mathbf e_t$, one gradient step gives

Decompose $\mathbf e$ along the eigenvectors of $\mathbf H$, with eigenvalues $\lambda_1 \ge \dots \ge \lambda_n > 0$ (the curvatures along each principal direction). Each component is multiplied by $(1-\eta\lambda_i)$ at every step. Two consequences follow, and you can see both in the lab.

1. Stability. Every factor must satisfy $|1-\eta\lambda_i| < 1$, so the iteration converges only if $$\eta < \frac{2}{\lambda_{\max}} .$$ The stretched bowl in the lab has curvatures 0.6 and 6, so it diverges for $\eta > 1/3$. Above that threshold each step overshoots the valley floor by more than it started with.
2. Speed. The slowest component shrinks by $|1-\eta\lambda_{\min}|$ per step. With the best fixed step, $\eta = 2/(\lambda_{\max}+\lambda_{\min})$, the error contracts by $$\frac{\kappa - 1}{\kappa + 1}, \qquad \kappa = \frac{\lambda_{\max}}{\lambda_{\min}},$$ per iteration. The condition number $\kappa$ measures how stretched the bowl is. For $\kappa = 10$ this is 0.82 per step; for $\kappa = 10^4$, common in real problems, it is 0.9998, and gradient descent crawls.

This is why the scale of input features matters. If one descriptor is in ångströms and another in electronvolts per atom, the loss surface is stretched along one axis, $\kappa$ is large, and training is slow. Standardising features (subtracting the mean and dividing by the standard deviation) is the cheapest optimisation trick there is.

Momentum

In a long, narrow valley, gradient descent zig-zags across the steep walls while making little progress along the floor. Polyak’s heavy-ball method adds a velocity that accumulates past gradients:

Components of the gradient that keep the same sign (along the valley floor) add up; components that flip sign from step to step (across the valley) cancel. With $\beta = 0.9$ the effective step along a consistent direction is up to $1/(1-\beta) = 10$ times larger. For quadratics, well-tuned momentum improves the convergence factor from about $1 - 1/\kappa$ to about $1 - 1/\sqrt\kappa$, a dramatic gain when $\kappa$ is large. Nesterov’s accelerated gradient, which evaluates the gradient at the look-ahead point $\boldsymbol\theta_t + \beta\mathbf v_t$, achieves this rate provably for convex problems.

Stochastic gradients

For a data set of $N$ points the gradient is a sum, $\nabla L = \frac1N\sum_i \nabla \ell_i$. When $N$ is a million training structures, computing it exactly for every step is wasteful. Stochastic gradient descent (SGD) estimates it from a random mini-batch $\mathcal B$ of, say, 32 points:

The estimate is unbiased but noisy. Robbins and Monro showed in 1951 that such a noisy iteration still converges if the learning rate decays so that $\sum_t \eta_t = \infty$ and $\sum_t \eta_t^2 < \infty$ (for example $\eta_t \propto 1/t$). In practice the learning rate follows a schedule: a short warm-up, then a cosine or step decay. A full pass through the data is called an epoch.

Adaptive methods: Adam

Different parameters often need different step sizes: weights feeding rare features see small, infrequent gradients. Adam (adaptive moment estimation) keeps running averages of the gradient and of its element-wise square,

corrects them for their initialisation at zero, $\hat{\mathbf m}_t = \mathbf m_t/(1-\beta_1^t)$ and $\hat{\mathbf s}_t = \mathbf s_t/(1-\beta_2^t)$, and steps by

Dividing by the root-mean-square gradient makes the step size roughly $\eta$ in every coordinate, whatever the gradient’s scale, which is a crude per-parameter correction for bad conditioning. The defaults $\beta_1 = 0.9$, $\beta_2 = 0.999$ work for most problems. In the lab, notice that Adam with $\eta = 0.1$ moves at a nearly constant speed even where the gradient is tiny, and how it can circle near the minimum. Most MLIP and GNN training scripts use Adam or its variant AdamW, which decouples weight decay from the adaptive step.

Second-order methods

If you can afford the Hessian, Newton’s method jumps straight to the minimum of the local quadratic model, $\boldsymbol\theta_{t+1} = \boldsymbol\theta_t - \mathbf H^{-1}\nabla L$, which removes the dependence on $\kappa$ entirely and converges quadratically near the minimum. For least-squares problems the Gauss–Newton approximation $\mathbf H \approx 2\mathbf J^\top\mathbf J$ (with $\mathbf J$ the Jacobian of the residuals) needs only first derivatives; adding a damping term that interpolates towards gradient descent gives Levenberg–Marquardt, the algorithm behind scipy.optimize.curve_fit and most curve-fitting software. L-BFGS builds a low-rank approximation of $\mathbf H^{-1}$ from recent gradients and is the standard for geometry optimisation of molecules and crystals. For networks with millions of parameters and noisy mini-batch gradients, first-order methods remain the practical choice.

Where the gradient comes from

Every method above needs $\nabla L$. For a model with a million parameters, finite differences would need a million extra loss evaluations per step. Reverse-mode automatic differentiation (backpropagation) computes the full gradient for about the cost of two or three evaluations of the loss, whatever the number of parameters. It is the subject of the lesson on neural networks in practice, and its forward-mode cousin reappears in differentiable molecular dynamics.

References and further reading

Fitting is not the goal

A model is useful only if it predicts data it has not seen. With enough parameters any data set can be fitted exactly: a polynomial of degree $N-1$ passes through $N$ points, noise included. Such a model has memorised the training set rather than learned the relationship behind it. The difference between error on the training data and error on new data is the generalisation gap, and controlling it is the central practical problem of machine learning.

In the lab, the dashed curve is the true function and each point carries Gaussian noise. Filled points are used for fitting; open points are held back for validation. Raise the polynomial degree from 0 and watch the two error curves on the right.

You should see three regimes:

• Underfitting (degree 0–2). The model is too rigid to follow the curve. Training and validation errors are both high and similar.
• A good fit (around degree 4–8 here). The validation error reaches its minimum, close to the noise level $\sigma$.
• Overfitting (high degree, $\lambda = 0$). The training error keeps falling towards zero, but the curve develops wild oscillations between the points and the validation error rises. Notice how large the weights $\lVert\mathbf w\rVert$ become: the model cancels huge positive and negative terms to thread through the noise.

The bias–variance decomposition

Why does a U-shape appear? Imagine repeating the whole experiment many times: new noise, new training set, new fit $\hat f$. At a fixed input $x$, with $y = f(x) + \varepsilon$ and noise variance $\sigma^2$, the expected squared error of the prediction splits exactly into three terms:

• Bias is systematic error: what the model gets wrong on average, because its family of functions cannot represent the truth. Rigid models have high bias.
• Variance is how much the fit changes from one training sample to the next. Flexible models fitted to few points have high variance; press New noisy sample at degree 11 and watch the curve jump about.
• Noise $\sigma^2$ is irreducible: no model can predict it. A validation RMSE well below $\sigma$ means something is wrong, usually leakage between training and test data.

Adding parameters lowers bias and raises variance. Adding data lowers variance without raising bias, which is why the overfitting at degree 11 largely disappears when you move the training-points slider to 60. More data is the most reliable regulariser.

Regularisation: preferring simple explanations

Instead of choosing the degree, keep a flexible model and penalise complexity. Ridge regression (Tikhonov regularisation) adds the squared size of the weights to the loss:

The extra $\lambda\mathbf I$ also makes the matrix well-conditioned, which is why ridge solutions are numerically stable even when features are nearly collinear. In the lab, keep degree 11 and slide $\lambda$ upwards: the oscillations disappear, the weights shrink, and the validation error drops well below its unregularised value. Too large a $\lambda$ flattens the model towards zero and underfits again; $\lambda$ is a hyperparameter that trades bias against variance continuously.

The Bayesian reading. Ridge regression is the most probable weight vector (the maximum a posteriori estimate) if, before seeing data, you believed each weight was drawn from a Gaussian $\mathcal N(0, \tau^2)$, with $\lambda = \sigma^2/\tau^2$. A penalty is a prior belief that the function is simple. The lasso uses $\lambda\lVert\mathbf w\rVert_1$ instead (a Laplace prior); it drives many weights exactly to zero and so selects features, which is useful when you suspect only a few descriptors matter. Neural networks use the same idea under the name weight decay, along with other implicit regularisers: early stopping, dropout and data augmentation.

Measuring generalisation honestly

Hold-out validation. Split the data into training, validation and test sets. Fit on the training set, choose hyperparameters (degree, $\lambda$, architecture, learning rate) on the validation set, and report performance on the test set once, at the very end. Using the test set to make choices turns it into a second validation set and makes the reported error optimistic.

$k$-fold cross-validation. With small data sets a single split is noisy. Divide the data into $k$ folds (5 or 10 is typical), train $k$ times, each time holding out one fold, and average the $k$ validation errors. Leave-one-out is the extreme case $k = N$.

In practice

With scikit-learn, a ridge model with a cross-validated penalty on standardised polynomial features takes a few lines:

import numpy as np
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import RidgeCV
from sklearn.model_selection import KFold, cross_val_score

rng = np.random.default_rng(0)
x = rng.uniform(0, 4, 40)[:, None]
y = 0.9*np.sin(2.6*x[:, 0])*np.exp(-0.25*x[:, 0]) + 0.15*x[:, 0] + 0.18*rng.normal(size=40)

model = make_pipeline(PolynomialFeatures(10), StandardScaler(),
                      RidgeCV(alphas=np.logspace(-8, 2, 41)))
cv = KFold(n_splits=5, shuffle=True, random_state=0)
scores = cross_val_score(model, x, y, cv=cv, scoring="neg_root_mean_squared_error")
print(f"5-fold RMSE: {-scores.mean():.3f} ± {scores.std():.3f}")

RidgeCV chooses $\lambda$ by an efficient leave-one-out estimate inside each training fold, so the outer cross-validation still measures an honest error. For grouped data, replace KFold with GroupKFold and pass a group label, such as the chemical system or the trajectory each sample came from.

References and further reading

Why uncertainty is not optional

A model predicts that a new alloy has a formation energy of −0.21 eV/atom. Should you synthesise it? That depends entirely on whether the prediction is good to ±0.01 or ±0.3 eV/atom, and on whether the alloy resembles anything in the training data. Scientific uses of machine learning, such as choosing the next experiment, flagging an MD configuration the potential has never seen, or deciding when to stop a screening campaign, all need the model to say how sure it is.

It helps to separate two kinds of uncertainty.

• Aleatoric uncertainty is noise in the data themselves: measurement error, or the variance of a property across nominally identical samples. More data does not remove it.
• Epistemic uncertainty is the model’s ignorance: it grows away from the training data and shrinks as data are added. This is the kind that tells you where to measure next.

Gaussian processes (GPs) model both, in closed form, and are the reference against which other uncertainty methods are judged.

A distribution over functions

A Gaussian process is a probability distribution over functions $f(x)$ such that the values at any finite set of inputs are jointly Gaussian. It is specified by a mean function (often zero) and a covariance function, or kernel, $k(x, x')$, which says how strongly the function values at two inputs are correlated:

The most common choice is the squared-exponential (RBF) kernel

It encodes a belief about smoothness: points closer than about one length-scale $\ell$ have nearly equal function values, while points much further apart are nearly independent. $\sigma_f$ sets the typical amplitude. Press Clear, then Draw 3 functions in the lab: before any data, these are samples from the prior. Change $\ell$ and draw again to see what the kernel assumes.

Conditioning on data

Suppose we have noisy observations $y_i = f(x_i) + \varepsilon_i$, $\varepsilon_i \sim \mathcal N(0,\sigma_n^2)$, collected in a vector $\mathbf y$. By the definition of a GP, the observations and the unknown value $f_* = f(x_*)$ at a new input are jointly Gaussian:

where $K_{ij} = k(x_i,x_j)$, $(\mathbf k_*)_i = k(x_i, x_*)$ and $k_{**} = k(x_*,x_*)$. Conditioning a multivariate Gaussian on some of its components is a standard identity, and it gives the posterior predictive distribution in closed form:

Read these two lines carefully, because they explain everything the lab shows.

• The mean $\mu_*$ is a weighted sum of kernel functions centred on the data points, $\mu_* = \sum_i \alpha_i\,k(x_i, x_*)$ with $\boldsymbol\alpha = (\mathbf K + \sigma_n^2\mathbf I)^{-1}\mathbf y$. Far from all data the kernel vanishes and the mean returns to the prior mean, zero.
• The variance starts at the prior variance $k_{**} = \sigma_f^2$ and is reduced by an amount that depends only on where the data are, not on their values. Uncertainty shrinks near observations and grows back over a distance of about $\ell$.
• Near an observation the variance does not go to zero but to about $\sigma_n^2$: the model knows the measurement is noisy.

In code, $(\mathbf K + \sigma_n^2\mathbf I)$ is factorised once by Cholesky decomposition, $\mathbf L\mathbf L^\top$, and both formulas follow by triangular solves. The cost is $\mathcal O(N^3)$ for $N$ training points, comfortable up to a few thousand points and the main limitation of exact GPs. Sparse approximations that summarise the data with $M \ll N$ inducing points reduce it to $\mathcal O(NM^2)$.

Choosing the kernel’s hyperparameters

$\ell$, $\sigma_f$ and $\sigma_n$ are not fitted by the formulas above; they are chosen by maximising the log marginal likelihood, the probability of the observed data under the GP prior with all possible functions integrated out:

The first term rewards fitting the data; the second penalises kernels that could explain too many different data sets (a very short length-scale can fit anything). The balance is an automatic Occam’s razor, and it is why GPs rarely overfit badly. Press Fit hyperparameters in the lab after adding six or eight points; the readout shows the log marginal likelihood for whatever settings you choose by hand, and you can try to beat the optimiser.

In practice

scikit-learn’s implementation maximises the marginal likelihood for you:

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernel

rng = np.random.default_rng(1)
X = rng.uniform(0, 10, 12)[:, None]
y = np.sin(1.1 * X[:, 0]) + 0.35 * np.cos(2.7 * X[:, 0]) + 0.08 * X[:, 0] + 0.05 * rng.normal(size=12)

kernel = ConstantKernel(1.0) * RBF(length_scale=1.0) + WhiteKernel(noise_level=0.01)
gp = GaussianProcessRegressor(kernel, normalize_y=True, n_restarts_optimizer=5, random_state=0)
gp.fit(X, y)                                   # maximises the log marginal likelihood
print(gp.kernel_)                              # fitted hyperparameters

Xq = np.array([[2.0], [5.0], [12.0]])
mu, sd = gp.predict(Xq, return_std=True)
for x, m, s in zip(Xq[:, 0], mu, sd):
    print(f"x = {x:5.1f}:  {m:6.3f} ± {2*s:.3f}")

The output shows the behaviour you want: ±0.13 at $x = 5$, inside the data; ±0.58 at $x = 2$, where samples are sparse; and ±1.3 at $x = 12$, outside the training range, where the model honestly falls back to its prior.

1.08**2 * RBF(length_scale=0.846) + WhiteKernel(noise_level=0.00489)
x =   2.0:   0.705 ± 0.575
x =   5.0:  -0.152 ± 0.127
x =  12.0:   0.208 ± 1.295

Uncertainty for neural networks

Neural networks do not provide uncertainties by default, and a single network can be confidently wrong far from its training data. The practical options are:

• Deep ensembles. Train $M$ networks (typically 5) from different random initialisations and use the spread of their predictions as the uncertainty. Simple, robust and the usual choice for machine-learned interatomic potentials, where the disagreement between ensemble members flags configurations that need a new DFT calculation.
• Monte Carlo dropout. Keep dropout active at prediction time and average many stochastic forward passes; cheaper, but usually less reliable than ensembles.
• Last-layer and Laplace approximations. Treat only the final layer (or a local Gaussian around the trained weights) probabilistically.
• Conformal prediction. A distribution-free wrapper: use the errors on a held-out calibration set to turn any model’s predictions into intervals with a guaranteed coverage, for example 90% of future values fall inside. It assumes only that calibration and future data are exchangeable, so it cannot rescue extrapolation.

Whatever the method, check the calibration: on held-out data, about 95% of true values should fall inside the predicted 95% intervals. Over-confident uncertainties are worse than none, because they are trusted.

References and further reading

Where the idea came from

In 1943 the neurophysiologist Warren McCulloch and the logician Walter Pitts asked a bold question: could a brain cell be described as a tiny logical device? A biological neuron receives signals from thousands of others through synapses, some exciting it and some inhibiting it; if the combined stimulation crosses a threshold, it fires a pulse down its axon. McCulloch and Pitts wrote this as arithmetic: add up the inputs, each multiplied by a strength, and output 1 if the total exceeds a threshold, 0 otherwise. In 1958 Frank Rosenblatt built a machine, the perceptron, that could learn those strengths from examples.

Modern artificial neurons are only loosely inspired by biology, but the arithmetic has hardly changed. Everything in this volume, up to graph networks and transformers, is built from millions of copies of the unit in this chapter.

What a neuron computes

A neuron takes several numbers in, $x_1, x_2, \dots, x_n$, and gives one number out. It does two things.

A positive weight means “more of this input makes me fire more”; a negative weight means “more of this input makes me fire less”; a weight near zero means “I ignore this input”. The bias shifts the threshold: a large positive bias makes the neuron fire even with weak input.

Seeing a neuron as geometry

With two inputs, every possible input is a point in a plane, and the neuron gives each point a value: the colours in the lab. Three facts are visible there and are worth understanding exactly.

The boundary is a straight line. The output is halfway (for the sigmoid, $y = 0.5$) where $z = 0$, that is where $w_1x_1 + w_2x_2 + b = 0$. That equation describes a straight line; with three inputs it is a flat plane, and in $n$ dimensions a hyperplane. Whatever activation you choose, a single neuron divides input space with a flat cut.

The weights point across the line. The weight vector $\mathbf w = (w_1, w_2)$ is perpendicular to the boundary and points towards higher output. Moving the input in the direction of $\mathbf w$ raises $z$ fastest; moving along the boundary does not change it at all. Rotating $\mathbf w$ in the lab rotates the boundary.

The bias slides the line; the size of the weights sets its sharpness. The distance from the origin to the boundary is $|b|/\lVert\mathbf w\rVert$. Doubling both weights and the bias leaves the boundary where it is but makes the transition from blue to orange twice as abrupt.

The activation function

The activation turns the weighted sum into a response. The common choices, all in the lab:

Why not leave the activation out? Then $y = \mathbf w\cdot\mathbf x + b$ is just linear regression (the chapter on fitting), and, as later chapters show, stacking any number of linear neurons still gives only a linear function. The bend in the activation is what lets networks of neurons represent curved relationships. Why smooth curves rather than the step? Because learning, the subject of the next chapter, works by nudging the weights in the direction that improves the output, and a step function gives no hint of which direction that is: its slope is zero everywhere except at the cliff.

A neuron as a logic gate, and its limit

With suitable weights a single step-function neuron can compute simple logic. With inputs 0 or 1:

• AND: $w_1 = w_2 = 1$, $b = -1.5$. Only $(1, 1)$ gives $z > 0$.
• OR: $w_1 = w_2 = 1$, $b = -0.5$. Any single 1 suffices.
• NOT (one input): $w = -1$, $b = 0.5$.

But no choice of weights computes XOR (output 1 when exactly one input is 1). Plot the four corners: $(0,1)$ and $(1,0)$ must be on one side of the line, $(0,0)$ and $(1,1)$ on the other, and no straight line does that. Minsky and Papert’s 1969 book Perceptrons made this limitation famous, and it contributed to a long decline of interest in neural networks. The way out, layers of neurons, is the subject of the chapter after next, where you will watch a small network solve XOR by bending the plane.

A neuron in science

A single neuron with a sigmoid output is exactly logistic regression, one of the most used models in science: will this compound be stable or not, is this spectrum from phase A or B. Its weights are directly interpretable as the influence of each feature. A neuron with no activation is linear regression. Everything more powerful is built by connecting many of them.

References and further reading

From choosing weights to learning them

In the previous chapter you set the weights by hand. With two inputs that is a game; with two hundred it is impossible. Learning means letting the data choose the weights. The recipe has three parts, and every network in this volume uses the same three:

1. a loss function that turns “how wrong are we?” into a single number;
2. its gradient, which says how the loss changes when each weight is nudged;
3. a rule that moves the weights a little way downhill, repeated many times.

Measuring error: from probability to cross-entropy

Suppose each training example is a point $\mathbf x_i$ with a label $y_i$ that is 1 (orange class) or 0 (blue class). A sigmoid neuron outputs a number between 0 and 1, and the natural reading is a probability: $\hat y_i = \sigma(\mathbf w\cdot\mathbf x_i + b)$ is the neuron’s belief that example $i$ is orange.

How good is that belief? A principled answer is the probability the model assigns to the labels that actually occurred, the likelihood. For one example it is $\hat y_i$ if the label is 1 and $1 - \hat y_i$ if it is 0, which can be written in one expression as $\hat y_i^{\,y_i}(1-\hat y_i)^{1-y_i}$. Multiplying over all examples and taking the negative logarithm (so that products become sums and maximising becomes minimising) gives the cross-entropy loss:

Look at what one term does. If the label is 1 and the neuron says $\hat y = 0.99$, the cost is $-\ln 0.99 \approx 0.01$. If it says $0.5$, the cost is $0.69$. If it confidently says $0.01$, the cost is $4.6$, and it grows without limit as the prediction approaches 0. Cross-entropy punishes confident mistakes severely, which is exactly what we want from a classifier.

The gradient, in one line

To go downhill we need the slope of $\mathcal L$ with respect to every weight. Here the sigmoid and the logarithm cancel in a way that is almost too neat. Use the chain rule on one example, with $z = \mathbf w\cdot\mathbf x + b$:

The meaning is intuitive. If the neuron says 0.9 for an example that is really 0, the error $\hat y - y = 0.9$ is large and positive, so every weight whose input was positive is pushed down. If the prediction is right, the error is near zero and nothing changes. Learning is driven by mistakes.

Walking downhill

Gradient descent repeats one update: compute the average gradient over the data and step against it,

where $\eta$, the learning rate, sets the step size.

The right-hand panel is the most important picture in machine learning. Every pair of weights is a point on the loss landscape; training is a walk downhill on it. Watch the path: it starts steep and fast, then slows down as the slope flattens near the bottom. With a large learning rate the walk overshoots and zig-zags. For a single sigmoid neuron with cross-entropy, the landscape is convex (bowl-shaped, with no false valleys), so gradient descent always reaches the best weights. That stops being true once neurons are stacked into networks, which is why the optimisation chapter is needed later.

Two things that go wrong already

• Perfectly separable data. If a straight line separates the classes exactly, the loss can always be lowered by making the weights larger and the boundary sharper, so they grow forever. In the lab, drag the classes apart and watch $\lVert\mathbf w\rVert$ keep increasing. The cure is regularisation, a penalty on large weights, discussed in the generalisation chapter.
• Unscaled inputs. If one input is in nanometres and another in kelvin, the landscape becomes a long thin valley and gradient descent zig-zags across it. Standardising inputs (subtract the mean, divide by the standard deviation) makes the bowl round.

The perceptron rule

Rosenblatt’s original perceptron used a step function and an even simpler update: for each misclassified example, add $\eta\,\mathbf x_i$ to the weights if it should have fired, and subtract it if it should not. He proved that if the classes can be separated by a line at all, this rule finds one in a finite number of steps (the perceptron convergence theorem). Logistic regression replaces the hard step by a smooth sigmoid so that every example, not only misclassified ones, contributes a gradient, and the output becomes a calibrated probability.

References and further reading

Many neurons side by side

A single neuron draws one straight cut (chapter on the neuron). Put several neurons side by side, all reading the same inputs, and you have a layer: each neuron computes its own weighted sum and activation, so a layer of $m$ neurons turns an input vector into an output vector of $m$ numbers. In matrix form, a whole layer is one line:

A network feeds the output of one layer into the next, $\mathbf h^{(1)} = \sigma(\mathbf W^{(1)}\mathbf x + \mathbf b^{(1)})$, $\mathbf h^{(2)} = \sigma(\mathbf W^{(2)}\mathbf h^{(1)} + \mathbf b^{(2)})$, and so on, and a final neuron reads the last representation. The layers between input and output are hidden layers: nobody tells them what to compute; they discover it during training.

A layer is a transformation of space

There are two ways to look at a network. The first, from the last chapter, asks what boundary each neuron draws in input space. The second, made famous by Christopher Olah and by 3Blue1Brown’s videos, is more revealing: a layer moves every point of space somewhere else. Think of the input plane as a rubber sheet with a grid drawn on it.

1. The matrix $\mathbf W$ is a linear transformation: it rotates, stretches, shears or reflects the sheet. Straight grid lines stay straight and parallel, and the origin stays put.
2. Adding $\mathbf b$ translates the whole sheet.
3. The activation $\tanh$ bends it: it squashes each coordinate into $(-1, 1)$, pulling far-away regions in towards the edges of a square while leaving the middle almost unchanged. Grid lines now curve.

The final neuron then draws a single straight line through the transformed sheet. The whole network succeeds if the hidden layers have deformed the sheet so that the two classes end up on opposite sides of some straight line.

Train the network on two moons, then play the transformation. The first linear map tilts and stretches the plane; the first tanh squeezes it so that the tips of the moons fold inwards; the second pair finishes the job, and in the final space a straight line separates the classes cleanly. Nothing about the classifier is curved. The curved boundary you would see in input space is a straight line in a bent space.

XOR, impossible for one neuron, falls to the same trick: the hidden layer folds the plane so that opposite corners land on the same side.

What layers cannot do: a lesson from topology

Now try the circle inside a ring with two neurons per layer. However long you train, the accuracy stalls. This is not a failure of the optimiser. Each layer here maps the plane continuously into the plane (the linear map and tanh are both continuous, and with two neurons the output is two-dimensional, and for these invertible weights, one-to-one). A continuous, one-to-one map of the plane can stretch and fold the rubber sheet but never tear it, and there is no way to get the inner disc outside the ring without tearing. The two classes stay linked, so no straight line in the final space can separate them.

The fix is a third neuron, a third dimension: with room to move upwards, a layer can lift the centre of the sheet out of the plane like a tent, and a flat plane can then slice between the tent’s peak and the ring. The general lesson, which Olah set out with pictures, is that width gives a network room to untangle; depth gives it more folds. Real networks use tens to thousands of neurons per layer for exactly this reason.

How expressive are layered networks?

The universal approximation theorem (Cybenko 1989; Hornik, Stinchcombe and White 1989) says that a network with one hidden layer of enough neurons and a non-polynomial activation can approximate any continuous function on a bounded region as closely as you like. The intuition is visible in the practical neural-network chapter: pairs of sigmoid neurons make “bumps”, and enough bumps can trace any curve.

The theorem says nothing about how many neurons are needed, and the answer can be astronomically many for one layer. Depth is what makes networks practical. A function made of repeated simple pieces, such as a sum over many atoms, or the layered structure of images, can need exponentially fewer neurons when represented by several layers each building on the last. Early layers learn simple features; later layers combine them into complex ones.

Why the nonlinearity is essential

Remove the activations and every layer is $\mathbf h = \mathbf W\mathbf x + \mathbf b$. Composing two such layers gives

a single linear layer again. A hundred linear layers are no more powerful than one: they can only rotate and stretch the sheet, never bend it. The bend is the whole point.

References and further reading

The credit-assignment problem

For one neuron, the gradient was a single neat formula. A network has layers of neurons, and a weight deep inside it affects the loss only through everything downstream of it. If the prediction is wrong, how much is each of a million weights to blame, and in which direction should each move? This is the credit-assignment problem, and its solution, backpropagation, is the algorithm that made deep learning possible.

The solution rests on one idea from school calculus, the chain rule, organised cleverly.

Computations as graphs

Write any calculation as a computational graph: a chain (or more generally a network) of simple operations, each node computing one value from the values flowing into it. Here is the smallest interesting network, one input, one hidden sigmoid neuron and one output neuron, with a squared-error loss, as a graph of seven nodes:

Each node knows two things about itself: how to compute its output from its input (the forward rule) and how its output changes when its input changes (its local derivative). That is all backpropagation needs.

The chain rule as a relay

The chain rule says that if $L$ depends on $h$, which depends on $z$, then $\partial L/\partial z = (\partial L/\partial h)(\partial h/\partial z)$. Applied along the whole chain, the gradient with respect to the first weight is a product of local derivatives:

Backpropagation evaluates this product from the right. Call $\bar v = \partial L/\partial v$ the gradient arriving at node $v$. Start with $\bar L = 1$ and pass it leftwards, each node multiplying by its local derivative:

The parameter gradients fall out along the way: $\bar b_2 = \bar{\hat y}$, $\bar w_2 = \bar m\, h$, $\bar b_1 = \bar z$. A worked example with the lab’s starting values ($x = 1.5$, $y = 0.2$, $w_1 = 0.8$, $b_1 = -0.4$, $w_2 = 1.3$, $b_2 = 0.1$):

so $\partial L/\partial w_1 = 0.4432 \times 1.5 = 0.665$, which the lab confirms against a finite difference.

Why it is cheap: reverse mode

Notice what the backward pass did. Each node was visited once going forward and once coming back, and every parameter’s gradient was produced along the way. The cost of all the gradients together is a small multiple (typically two to three times) of the cost of computing the loss once, whether the network has seven parameters or seventy billion.

Compare the obvious alternative: nudge each weight in turn and re-run the network to see how the loss changes (finite differences). That needs one extra forward pass per parameter, a million passes for a million weights. Backpropagation is a special case of reverse-mode automatic differentiation: for a function with many inputs (the weights) and one output (the loss), sweeping derivatives backwards from the output is the efficient direction. Forward-mode differentiation, which carries derivatives along with the values from the inputs, is efficient in the opposite situation, few inputs and many outputs; you will meet it in differentiable molecular dynamics.

In a real network each node is a whole layer, the values are vectors and the local derivatives are matrices, but nothing else changes. For a layer $\mathbf h = \sigma(\mathbf W\mathbf x + \mathbf b)$ the backward rules are:

Vanishing and exploding gradients

The backward pass multiplies many factors together, one per layer. If they are typically smaller than 1, the gradient shrinks exponentially with depth and early layers stop learning (vanishing gradients); if larger than 1, it grows exponentially (exploding gradients). The sigmoid’s derivative is at most 0.25, so deep sigmoid networks suffered badly. The cures are now standard: activations whose slope does not shrink (ReLU and its relatives), careful initialisation of the weights (the practical chapter), normalisation layers, and residual connections $\mathbf h^{(l+1)} = \mathbf h^{(l)} + f(\mathbf h^{(l)})$, which give the gradient a direct path backwards through every layer.

A short history

Reverse-mode differentiation was described by Seppo Linnainmaa in his 1970 master’s thesis, applied to neural networks by Paul Werbos in his 1974 doctoral thesis, and popularised by David Rumelhart, Geoffrey Hinton and Ronald Williams in a 1986 Nature paper, which showed that hidden layers trained this way learn useful internal representations. Today every deep-learning framework (PyTorch, JAX, TensorFlow) builds the computational graph as your code runs and performs the backward pass automatically when you call loss.backward().

References and further reading

From fixed features to learned features

A linear model predicts with fixed features: $\hat y = \sum_k w_k\,\phi_k(x)$. Someone has to choose the $\phi_k$, and for a complex relationship (the energy of a configuration of 200 atoms, say) nobody knows the right ones. A neural network keeps the same form but makes the features themselves adjustable:

Each term is a neuron: it takes a weighted sum of its inputs, adds a bias, and passes the result through a fixed nonlinear activation function $\sigma$. The weights $w_k, b_k$ decide where and how sharply each neuron bends; the output weights $v_k$ decide how much each one contributes. Everything is fitted together by gradient descent. This is a network with one hidden layer of $H$ neurons, a multilayer perceptron (MLP).

With one tanh neuron the network can only draw one smooth step. Two neurons with opposite signs can draw a bump. Many bumps placed side by side can trace any curve, which is the intuition behind the theorem below. Watch the faint curves during training: the neurons move, stretch and flip until their sum matches the data.

Layers, in matrix form

For a vector input $\mathbf x \in \mathbb R^{d}$, a layer with $H$ neurons computes

with $\sigma$ applied element-wise. A deep network stacks $L$ such layers, each feeding the next:

Without the nonlinearity the whole stack would collapse into a single linear map, $\mathbf W^{(L)}\cdots\mathbf W^{(1)}\mathbf x$, no more expressive than linear regression. The activation is what makes depth useful.

Common activations.

Why networks can represent almost anything

The universal approximation theorem (Cybenko 1989; Hornik, Stinchcombe and White 1989) states that a network with a single hidden layer and a non-polynomial activation can approximate any continuous function on a bounded domain to any accuracy, given enough neurons. It is an existence result: it says nothing about how many neurons are needed, or whether gradient descent will find the weights. In practice, depth is far more efficient than width: functions with repeated structure, such as sums over many atoms, can need exponentially fewer neurons when represented by several layers.

Backpropagation: the gradient for the price of one extra pass

Training needs $\partial L/\partial\theta$ for every weight. Consider a single data point and squared loss $L = (\hat y - y)^2$. The chain rule, applied from the output backwards, gives a recursion. Define the error signal of layer $l$ as $\boldsymbol\delta^{(l)} = \partial L/\partial\mathbf z^{(l)}$. Then

and the gradients with respect to the parameters of each layer are outer products of what flowed in and what flows back:

The cost of the backward pass is about the same as the forward pass: each weight is touched once going forward and once coming back. So the gradient with respect to all parameters, whether there are ten or ten billion, costs roughly two to three forward evaluations. This is what made deep learning possible. The algorithm was popularised for neural networks by Rumelhart, Hinton and Williams (1986). It is a special case of reverse-mode automatic differentiation, which modern frameworks (PyTorch, JAX, TensorFlow) apply automatically to any program built from differentiable operations.

Here is the whole algorithm in NumPy, a one-hidden-layer network fitted to $\sin 3x$:

import numpy as np

rng = np.random.default_rng(0)
x = np.linspace(-1, 1, 64)[:, None]               # inputs, shape (N, 1)
y = np.sin(3 * x)                                  # targets

H = 16                                             # hidden neurons
W1 = rng.normal(0, 1.0, (1, H)); b1 = np.zeros(H)
W2 = rng.normal(0, 1/np.sqrt(H), (H, 1)); b2 = np.zeros(1)
lr = 0.05

for step in range(5001):
    # forward pass
    z = x @ W1 + b1                    # (N, H) pre-activations
    a = np.tanh(z)                     # (N, H) activations
    yhat = a @ W2 + b2                 # (N, 1) output
    loss = np.mean((yhat - y) ** 2)
    # backward pass: chain rule, from the loss back to every weight
    d_yhat = 2 * (yhat - y) / len(x)   # dL/dyhat
    dW2 = a.T @ d_yhat;  db2 = d_yhat.sum(0)
    d_a = d_yhat @ W2.T                # dL/da
    d_z = d_a * (1 - a ** 2)           # tanh'(z) = 1 - tanh(z)^2
    dW1 = x.T @ d_z;     db1 = d_z.sum(0)
    # plain gradient-descent update
    W1 -= lr * dW1; b1 -= lr * db1; W2 -= lr * dW2; b2 -= lr * db2
    if step % 1000 == 0:
        print(f"step {step:5d}   loss {loss:.5f}")

The loss falls from 0.25 to 0.0044 in 5000 steps of plain gradient descent. The same network in PyTorch, where loss.backward() does the backward pass for you and Adam replaces the update, reaches $4\times10^{-6}$ in 2000 steps:

import torch

torch.manual_seed(0)
x = torch.linspace(-1, 1, 64).unsqueeze(1)
y = torch.sin(3 * x)

model = torch.nn.Sequential(torch.nn.Linear(1, 16), torch.nn.Tanh(), torch.nn.Linear(16, 1))
opt = torch.optim.Adam(model.parameters(), lr=1e-2)

for step in range(2001):
    opt.zero_grad()
    loss = torch.mean((model(x) - y) ** 2)
    loss.backward()            # reverse-mode autodiff fills .grad for every parameter
    opt.step()
    if step % 500 == 0:
        print(f"step {step:5d}   loss {loss.item():.2e}")

Initialisation and training in practice

• Initial weights must be random, or every neuron in a layer receives the same gradient and they stay identical forever. Their scale matters too: too large and tanh neurons saturate (their gradient vanishes); too small and the signal dies out through the layers. Glorot (Xavier) initialisation uses variance $1/d_{\text{in}}$ or $2/(d_{\text{in}}+d_{\text{out}})$ for tanh; He initialisation uses $2/d_{\text{in}}$ for ReLU. Press New initial weights in the lab: different starting points reach different, comparably good fits, a reminder that the loss surface of a network has many minima.
• Standardise inputs and outputs, for the conditioning reasons in the gradient-descent lesson.
• Monitor a validation set and stop or decay the learning rate when it stops improving (the lesson on generalisation).