8.1 Neural networks

A neural network stacks many simple units. Each one takes a weighted sum of its inputs and bends it. Together, and with enough data, they can learn almost any shape, and even their own measurements.

import warnings

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from sklearn import datasets
from sklearn.exceptions import ConvergenceWarning
from sklearn.model_selection import train_test_split
from sklearn.neural_network import MLPClassifier, MLPRegressor

from statly import style

style.use()
rng = np.random.default_rng(1)
warnings.filterwarnings(
    "ignore", category=ConvergenceWarning
)  # we stop training early on purpose

digits = datasets.load_digits()  # 1,797 handwritten digits, 8 × 8 pixels each
X_digits = digits.data / 16  # pixel darkness from 0 to 1
engel = sm.datasets.engel.load_pandas().data  # income and food spending, Belgian households
engel_income, engel_food = engel.income.to_numpy(), engel.foodexp.to_numpy()
X_train, X_test, y_train, y_test = train_test_split(
    X_digits, digits.target, test_size=0.3, random_state=1
)

8.1.1 A unit and a network

A single unit (or neuron) does two things. It takes a weighted sum of its inputs plus a constant, exactly like a regression. Then it passes the sum through an activation: a simple bend. The most common one, the ReLU, keeps positive values and turns negative ones into zero.

A network puts units side by side in a layer, feeds the inputs to all of them, and feeds their outputs on to the next layer. The last layer gives the prediction.

LAYERS = [("input", 3, style.BLUE), ("hidden", 4, style.GREEN), ("output", 1, style.CORAL)]


def layer_positions(count: int, x: float) -> list[tuple[float, float]]:
    """Evenly spaced positions for the units of one layer."""
    return [(x, (i + 1) / (count + 1)) for i in range(count)]


def plot_network(ax: Axes) -> None:
    """Circles for units, lines for weights, layer names underneath."""
    positions = [
        layer_positions(count, x) for x, (_, count, _) in zip((0.15, 0.5, 0.85), LAYERS)
    ]
    for left, right in zip(positions, positions[1:]):
        for a in left:
            for b in right:
                ax.plot([a[0], b[0]], [a[1], b[1]], color=style.INK, lw=0.8, alpha=0.5)
    for layer, (name, _, color) in zip(positions, LAYERS):
        for x, y in layer:
            ax.add_patch(plt.Circle((x, y), 0.045, color=color, alpha=0.85, zorder=3))
        ax.text(layer[0][0], 0.02, name, ha="center", fontsize=10, color=style.INK)
    ax.set(xlim=(0, 1), ylim=(0, 1), aspect="equal")
    ax.set_axis_off()


fig, ax = plt.subplots(figsize=(5.5, 3.2))
plot_network(ax)
plt.show()

Figure 8.1.1: A network with three inputs, one hidden layer of four units and one output. Every line is a weight to be learned: 3 × 4 + 4 × 1 = 16 weights, plus a constant for each unit.

Without the bend, a stack of weighted sums would just be one big weighted sum: a linear model. The activation is what lets a network draw curves. In symbols, a network with one hidden layer of K units computes, for an input $x = (x_1, \dots, x_p)$:

$$h_k = g\Big(w_{k0} + \sum_{j=1}^{p} w_{kj}\, x_j\Big), \quad k = 1, \dots, K, \qquad f(x) = \beta_0 + \sum_{k=1}^{K} \beta_k\, h_k$$

where g is the activation. The hidden values $h_k$ are new measurements the network makes for itself; the output is a linear regression on them. For classes, the output layer has one unit per class, and the softmax turns their scores $z_1, \dots, z_M$ into probabilities that add up to 1, $P(\text{class } m) = e^{z_m} / \sum_l e^{z_l}$, the multinomial logistic regression of Chapter 4.4.

z = np.linspace(-4, 4, 300)
ACTIVATIONS = [
    ("sigmoid: 1 / (1 + e^(−z))", 1 / (1 + np.exp(-z)), style.BLUE),
    ("tanh", np.tanh(z), style.GREEN),
    ("ReLU: max(0, z)", np.maximum(0, z), style.CORAL),
]

fig, axes = style.subplots(3, height=2.6, width=10.5)
for ax, (name, values, color) in zip(axes, ACTIVATIONS):
    ax.plot(z, values, color=color, lw=2.2)
    ax.axhline(0, color=style.INK, lw=0.6)
    ax.axvline(0, color=style.INK, lw=0.6)
    ax.set(title=name, xlabel="input z")
axes[0].set_ylabel("output g(z)")
fig.tight_layout()
plt.show()

Figure 8.1.2: Common activations. The sigmoid (the logistic S of Chapter 4.4) and tanh squash any input into a range; they were the classics. The ReLU, max(0, z), keeps positive values and zeroes the rest; it is the usual choice today because it is cheap and its slope does not fade for large inputs.

8.1.2 Hinges add up to a curve

With one input, each ReLU unit is a hinge: flat on one side, a straight slope on the other, with the bend wherever its weights put it. The output adds the hinges up, each with its own weight. A few hinges make a rough curve; many make a smooth one.

x = np.sort(rng.uniform(0, 1, 200))
truth = np.sin(2.2 * np.pi * x) * 0.9 + x
noisy = truth + rng.normal(0, 0.25, len(x))
GRID = np.linspace(0, 1, 500)


def fit_network(units: int, starts: int = 10) -> MLPRegressor:
    """A network with one hidden layer of ReLU units: the best of several random starts."""
    fits = [
        MLPRegressor(
            hidden_layer_sizes=(units,), solver="lbfgs", max_iter=20_000, random_state=seed
        ).fit(x[:, None], noisy)
        for seed in range(starts)
    ]
    scores = [net.score(x[:, None], noisy) for net in fits]
    return fits[int(np.argmax(scores))]


def plot_hinges(ax: Axes, net: MLPRegressor) -> None:
    """Each hidden unit's output times its weight on the output."""
    hidden = np.maximum(0, GRID[:, None] * net.coefs_[0] + net.intercepts_[0])
    for j, color in zip(range(hidden.shape[1]), style.PALETTE * 2):
        ax.plot(GRID, hidden[:, j] * net.coefs_[1][j, 0], color=color, lw=1.4)
    ax.axhline(0, color=style.INK, lw=0.6)
    ax.set(title="Six hinges", xlabel="x", ylabel="contribution")


def plot_fit(ax: Axes, net: MLPRegressor, title: str) -> None:
    """The data and the network's prediction."""
    ax.scatter(x, noisy, s=6, color=style.GREY, alpha=0.5)
    ax.plot(GRID, net.predict(GRID[:, None]), color=style.CORAL, lw=2)
    ax.set(title=title, xlabel="x", ylabel="y")


two, six = fit_network(2), fit_network(6)
fig, axes = style.subplots(3, height=3, width=10.5)
plot_fit(axes[0], two, "Two units")
plot_hinges(axes[1], six)
plot_fit(axes[2], six, "Six units: their sum")
fig.tight_layout()
plt.show()

Figure 8.1.3: The simulated curve from Chapter 1.3, fitted by networks with one hidden layer. Left: two units give two hinges, too few for the shape. Middle: with six units, each coloured line is one unit’s weighted hinge. Right: their sum (red) is a curve made of straight pieces that follows the shape.

This is why networks are called universal: with enough units, sums of hinges can follow any curve, in any number of dimensions. The same power lets them fit noise, so everything from Chapters 6.3 and 6.4 applies: hold data back, and limit the flexibility.

8.1.3 Learning, step by step

A network has thousands or millions of weights, and no formula sets them. Training needs two things: a loss that says how wrong the predictions are, and a way to reduce it.

• For a number, the loss is usually the squared error, $\sum_i (y_i - f(x_i))^2$, as in least squares.
• For classes, it is the cross-entropy, $-\sum_i \log \hat p_{i, y_i}$: minus the log of the probability the network gave to the true class. It is the negative log-likelihood of Chapters 2.3 and 4.4.

Gradient descent reduces the loss step by step. The gradient is the direction in which the loss rises fastest, one slope per weight. Step a little the other way and repeat:

$$w \leftarrow w - \rho\, \nabla L(w)$$

The step size ρ is the learning rate. Too small, and training crawls; too large, and it overshoots and can diverge. A loss with only two weights can be drawn as a map:

food_x = (
    engel_income - engel_income.mean()
) / engel_income.std() + 1.5  # off-centre: a tilted bowl
food_y = (engel_food - engel_food.mean()) / engel_food.std()
a_grid, b_grid = np.meshgrid(np.linspace(-3.2, 1.2, 200), np.linspace(-0.6, 2.4, 200))
loss_surface = (
    (food_y[None, None, :] - a_grid[..., None] - b_grid[..., None] * food_x[None, None, :]) ** 2
).mean(axis=2)
best_b, best_a = np.polyfit(food_x, food_y, 1)


def descend(
    rate: float, steps: int = 30, start: tuple[float, float] = (0.6, -0.3)
) -> np.ndarray:
    """The path of gradient descent on the mean squared error."""
    a, b = start
    path = [(a, b)]
    for _ in range(steps):
        residual = food_y - a - b * food_x
        a, b = a + rate * 2 * residual.mean(), b + rate * 2 * (residual * food_x).mean()
        path.append((a, b))
    return np.array(path)


def plot_descent(ax: Axes, rate: float, title: str) -> None:
    """Contours of the loss and one descent path."""
    ax.contour(
        a_grid,
        b_grid,
        loss_surface,
        levels=[0.3, 0.6, 1, 1.6, 2.5, 4, 6],
        colors=style.INK,
        linewidths=0.7,
    )
    path = descend(rate)
    ax.plot(path[:, 0], path[:, 1], "o-", ms=3, lw=1.2, color=style.CORAL)
    ax.plot(best_a, best_b, "x", ms=10, mew=2.5, color=style.BLUE)
    ax.set(title=title, xlabel="intercept", xlim=(-3.2, 1.2), ylim=(-0.6, 2.4))


fig, axes = style.subplots(3, height=3.2, width=10.5, sharey=True)
for ax, (rate, title) in zip(
    axes,
    ((0.02, "Learning rate 0.02"), (0.12, "Learning rate 0.12"), (0.26, "Learning rate 0.26")),
):
    plot_descent(ax, rate, title)
axes[0].set_ylabel("slope")
fig.tight_layout()
plt.show()

Figure 8.1.4: Gradient descent on the loss of a straight line (intercept and slope) fitted to the standardised food spending data, drawn as contours. The predictor is not centred, so the bowl is long, narrow and tilted, as loss surfaces usually are. Thirty steps from the same start with three learning rates. Left: too small, still far from the bottom (cross). Middle: about right, it drops into the valley at once and then follows it down. Right: too large, every step overshoots across the narrow valley and the zigzag grows: it diverges.

In a network the gradient has millions of parts, one per weight. Backpropagation computes all of them in one backward pass through the network, using the chain rule: the effect of a weight in an early layer on the loss is the product of the effects along the path from that weight to the output. And instead of using all the data for every step, stochastic gradient descent uses a small random minibatch of cases per step. Each step is noisier, but there are many more of them for the same work, and the noise even helps escape poor spots. One pass through all the training data is an epoch.

The 1,797 handwritten digits that come with scikit-learn are 8 × 8 pixel images. A network with 64 inputs (one per pixel), one hidden layer and 10 outputs (one per digit) learns to read them:

EPOCHS = 60
network = MLPClassifier(hidden_layer_sizes=(64,), learning_rate_init=0.002, random_state=1)
train_accuracy, test_accuracy = [], []
for _ in range(EPOCHS):
    network.partial_fit(X_train, y_train, classes=np.arange(10))
    train_accuracy.append(network.score(X_train, y_train))
    test_accuracy.append(network.score(X_test, y_test))


def plot_examples(ax: Axes, count: int = 24) -> None:
    """A grid of training images with their labels."""
    sheet = np.ones((4 * 9, 6 * 9))
    for i in range(count):
        row, col = divmod(i, 6)
        sheet[row * 9 : row * 9 + 8, col * 9 : col * 9 + 8] = 1 - X_train[i].reshape(8, 8)
    ax.imshow(sheet, cmap="gray", vmin=0, vmax=1)
    ax.set_axis_off()
    ax.set_title("Training digits")


def plot_curves(ax: Axes) -> None:
    """Accuracy per epoch on training and held-out digits."""
    epochs = np.arange(1, EPOCHS + 1)
    ax.plot(epochs, train_accuracy, color=style.GREY, label="training digits")
    ax.plot(epochs, test_accuracy, color=style.CORAL, label="held-out digits")
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8, loc="lower right")
    ax.set(title="Learning", xlabel="epoch", ylabel="accuracy", ylim=(0.5, 1.01))


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_examples(axes[0])
plot_curves(axes[1])
fig.tight_layout()
plt.show()

Figure 8.1.5: A network with 64 hidden units learning to read digits. Left: some training images. Right: accuracy on the training digits and on held-out digits, epoch by epoch. Both climb fast; then the training accuracy keeps creeping towards 100% while the held-out accuracy levels off just below it. That gap is the start of overfitting.

On digits it has never seen, the network is right about 98% of the time.

8.1.4 Keeping a network from overfitting

A network with many units can memorise its training data. Three guards are standard:

• Weight decay: add $\lambda \sum w^2$ to the loss, the ridge penalty of Chapter 6.4, which keeps weights small and the function smooth.
• Early stopping: keep an eye on the loss on held-out data, and stop when it stops falling.
• Dropout: at every training step, switch off a random share of the units (say half), so that no unit can rely on any other; at prediction time all units are used, with their outputs scaled to match.

curve_rng = np.random.default_rng(1)
curve_x = curve_rng.uniform(0, 1, 540)
curve_y = np.sin(2.2 * np.pi * curve_x) * 0.9 + curve_x + curve_rng.normal(0, 0.35, 540)
fit_x, held_x, fit_y, held_y = (
    curve_x[:40, None],
    curve_x[40:, None],
    curve_y[:40],
    curve_y[40:],
)
big = MLPRegressor(hidden_layer_sizes=(500,), learning_rate_init=0.01, alpha=0, random_state=0)
train_loss, held_loss = [], []
for _ in range(3_000):
    big.partial_fit(fit_x, fit_y)
    train_loss.append(np.mean((fit_y - big.predict(fit_x)) ** 2))
    held_loss.append(np.mean((held_y - big.predict(held_x)) ** 2))


def plot_dropout(ax: Axes) -> None:
    """A small network with some hidden units crossed out."""
    layers = [(3, 0.15), (6, 0.5), (1, 0.85)]
    positions = [[(x, (i + 1) / (count + 1)) for i in range(count)] for count, x in layers]
    dropped = {1, 4}
    for left, right, layer in (
        (positions[0], positions[1], 0),
        (positions[1], positions[2], 1),
    ):
        for i, a in enumerate(left):
            for j, b in enumerate(right):
                gone = (layer == 0 and j in dropped) or (layer == 1 and i in dropped)
                ax.plot(
                    [a[0], b[0]],
                    [a[1], b[1]],
                    color=style.INK,
                    lw=0.7,
                    alpha=0.08 if gone else 0.5,
                )
    for layer, units in enumerate(positions):
        for i, (x, y) in enumerate(units):
            gone = layer == 1 and i in dropped
            ax.add_patch(
                plt.Circle(
                    (x, y),
                    0.04,
                    color=style.GREY if gone else (style.BLUE, style.GREEN, style.CORAL)[layer],
                    alpha=0.3 if gone else 0.85,
                )
            )
            if gone:
                ax.plot([x - 0.04, x + 0.04], [y - 0.04, y + 0.04], color=style.CORAL, lw=2)
                ax.plot([x - 0.04, x + 0.04], [y + 0.04, y - 0.04], color=style.CORAL, lw=2)
    ax.set(xlim=(0, 1), ylim=(0, 1), aspect="equal", title="Dropout")
    ax.set_axis_off()


def plot_losses(ax: Axes) -> None:
    """Training and held-out loss by epoch, with the early-stopping point."""
    epochs = np.arange(1, len(train_loss) + 1)
    ax.plot(epochs, train_loss, color=style.GREY, label="training")
    ax.plot(epochs, held_loss, color=style.CORAL, label="held out")
    best = int(np.argmin(held_loss))
    ax.plot(best + 1, held_loss[best], "o", ms=8, color=style.CORAL)
    ax.axhline(0.35**2, color=style.INK, ls=":", lw=1)
    ax.legend(fontsize=8)
    ax.set(
        title="Early stopping", xlabel="epoch", ylabel="mean squared error", ylim=(0.04, 0.3)
    )


fig = plt.figure(figsize=(10, 3.2))
grid = fig.add_gridspec(1, 2, width_ratios=[1, 1.8])
plot_dropout(fig.add_subplot(grid[0]))
plot_losses(fig.add_subplot(grid[1]))
fig.tight_layout()
plt.show()

Figure 8.1.6: Left: dropout in one training step; the crossed-out units are switched off for this step only, and a different random set is dropped at the next. Right: a large network (500 hidden units) trained on only 40 noisy points from the curve of Chapter 1.3. The training loss (grey) keeps falling; the loss on held-out points (red) falls, then rises as the network starts to fit noise. Early stopping keeps the weights from the lowest point (dot).

The dotted line is the noise level: no model can beat it on new data. Training loss below it is a sure sign of memorising.

8.1.5 What the hidden units look for

Each hidden unit has one weight per pixel, so its weights can be drawn as an 8 × 8 image. Red pixels push the unit up, blue ones push it down. The patterns are clearest in a smaller network, with 16 hidden units and a penalty on large weights; it still reads 97% of the held-out digits. The units have learned their own measurements: strokes, loops and gaps that help tell digits apart.

small_network = MLPClassifier(
    hidden_layer_sizes=(16,), alpha=1.0, max_iter=2000, random_state=1
).fit(X_train, y_train)


def plot_weights(axes: np.ndarray, units: int = 16) -> None:
    """One small image per hidden unit: its weight on each pixel."""
    weights = small_network.coefs_[0]
    for ax, j in zip(axes, range(units)):
        w = weights[:, j].reshape(8, 8)
        limit = abs(w).max()
        ax.imshow(w, cmap="RdBu_r", vmin=-limit, vmax=limit)
        ax.set_axis_off()


fig, axes = plt.subplots(2, 8, figsize=(8, 2.2))
plot_weights(axes.ravel())
plt.show()

Figure 8.1.7: The pixel weights of the 16 hidden units of a smaller network. Each unit responds to its own pattern of strokes (red: ink here raises the unit, blue: ink here lowers it). Nobody told the network what to look for.

This is what sets networks apart from the other methods in this book. A tree or a regression works with the measurements it is given. A network with several layers builds new measurements in each layer from the ones below: edges from pixels, shapes from edges, objects from shapes. That is deep learning (Chapters 8.2 to 8.4), and it is why networks dominate for images, sound and text, where nobody can write the right measurements down by hand.

8.1.6 Interpolation and double descent

Chapter 1.3 drew the U: as flexibility grows, test error first falls, then rises. Huge networks seem to break that rule: they fit their training data perfectly and still predict well. A simple experiment shows the curious shape behind this. Take 100 training digits, compute p random ReLU features of their pixels (a stand-in for a network’s hidden units, with random, untrained weights), and fit a least squares classifier on them, for growing p. When there are more features than images, many fits pass through every training image exactly; take the one with the smallest coefficients.

N_TRAIN = 100
FEATURE_COUNTS = np.unique(np.geomspace(5, 3000, 36).astype(int))
pixels = digits.data / 16
labels = digits.target


def error_rate(p: int, repetitions: int = 20) -> float:
    """Median held-out error of the minimum-norm least squares fit on p random ReLU features."""
    errors = []
    for r in range(repetitions):
        local = np.random.default_rng(r)
        order = local.permutation(len(labels))
        train, test = order[:N_TRAIN], order[N_TRAIN : N_TRAIN + 500]
        weights, offsets = local.normal(0, 1, (64, p)) / 8, local.normal(0, 0.5, p)

        def features(rows: np.ndarray) -> np.ndarray:
            return np.maximum(0, pixels[rows] @ weights + offsets)

        coef = (
            np.linalg.pinv(features(train)) @ np.eye(10)[labels[train]]
        )  # smallest coefficients that fit
        errors.append(np.mean((features(test) @ coef).argmax(axis=1) != labels[test]))
    return float(np.median(errors))


errors = [error_rate(p) for p in FEATURE_COUNTS]

fig, ax = plt.subplots(figsize=(8, 3.3))
ax.plot(FEATURE_COUNTS, errors, "o-", ms=3, color=style.CORAL)
ax.axvline(N_TRAIN, color=style.INK, ls="--", lw=1)
ax.text(N_TRAIN * 1.1, max(errors) * 0.92, "features =\ntraining images", fontsize=9)
ax.set_xscale("log")
ax.yaxis.set_major_formatter(style.PERCENT)
ax.set(xlabel="random features p (log scale)", ylabel="error rate on new digits")
plt.show()

Figure 8.1.8: Error rate on 500 held-out digits of a least squares classifier on p random ReLU features, trained on 100 digits, median over 20 repetitions. Up to about 50 features the familiar U begins; near p = 100, where the fit can just pass through every training image, the error peaks; with many more features than images it falls again, below the first minimum: double descent.

The peak sits where the model has just enough flexibility to pass through every training image, and has to contort itself to do so. With far more flexibility than that, there are many perfect fits, and the one with the smallest coefficients is gentle. Large networks trained by gradient descent seem to find such gentle solutions on their own. The bias and variance trade-off still holds; it is the meaning of “flexibility” that is more subtle than the number of parameters.

8.1.7 When to use a network

• Images, sound, text: yes. Special layers that reuse the same weights across an image (convolutional networks) or across a sentence (transformers) make them far better than anything else here.
• Tables of a few hundred or thousand rows: usually not. Regression, GAMs and boosted trees are as accurate, faster, and easier to check.
• Always: hold data back, and compare with a simple model.

Networks are also hard to read. Methods that explain a single prediction exist, but a network is never as transparent as a regression table.

8.1.8 In one picture

• A unit takes a weighted sum and bends it with an activation; a network stacks layers of units.
• With one input, ReLU units are hinges, and their sum can follow any curve.
• $h_k = g(w_{k0} + \sum_j w_{kj} x_j)$, $f(x) = \beta_0 + \sum_k \beta_k h_k$; softmax for classes.
• Training minimises a loss (squared error, cross-entropy) by gradient descent, $w \leftarrow w - \rho \nabla L$; backpropagation computes the gradient, minibatches make it cheap. The learning rate ρ must be neither too small nor too large.
• Guard against overfitting with held-out data, early stopping, weight decay and dropout.
• Past the point where a model fits every training point, test error can fall again: double descent.
• Deep networks learn their own measurements, which is why they win on images, sound and text, and rarely on small tables.

Next: networks that see, classifying images.