Scientific Machine Learning:


Artificial Neural Networks (ANN)

An artificial neural network is the function approximator used throughout scientific machine learning. This notebook builds it up from the perceptron of the previous notebook to the multi-layer perceptron, explains how it is trained by backpropagation and gradient descent, and implements it in PyTorch on a regression problem and on image classification. The same construction, with the same training loop, will later be asked to satisfy differential equations rather than to fit labels.


Learning Objectives

After studying this notebook, you should be able to

  • state precisely what a perceptron computes and what it cannot do, and explain why a better optimizer or more data cannot repair that limitation;
  • explain how stacking perceptrons removes it, and read a hidden layer as a learned change of coordinates in which the classes become linearly separable;
  • build a small network in PyTorch with torch.nn.Sequential, train it with an optimizer and a loss function, and recover its weights and biases afterwards;
  • interpret a one-hidden-layer ReLU network as a sum of hinge functions, and predict how the fit improves as the number of hidden units is increased;
  • describe a general artificial neural network in terms of layers, weights, activation functions and a loss;
  • explain backpropagation as the chain rule applied recursively to a layered computation, and distinguish it from the optimization algorithm that consumes its gradients;
  • train a fully connected network on MNIST, choose the output layer and loss for a multi-class problem, and read the gap between training and test accuracy as the signature of overfitting.




1. Recall Perceptron

This notebook builds the artificial neural network from a single component that is already familiar. Before adding anything new, it is worth restating precisely what a perceptron is and precisely where it fails, because the entire architecture developed in this chapter is a response to that one failure.


The plan of the notebook is as follows. Chapter 1 recalls the perceptron and its limitation. Chapter 2 shows how stacking perceptrons removes that limitation, and interprets hidden layers as learned feature extractors. Chapters 3 and 4 examine the internal behavior of a trained network on a classification and a regression example. Chapter 5 assembles these observations into a general description of the artificial neural network. Chapter 6 addresses how such a network is trained, arriving at the backpropagation algorithm by way of recursion and dynamic programming. Chapter 7 applies everything to handwritten digit recognition.

In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('W5i9OA0bW-A', width = '560', height = '315', extras = ['referrerpolicy="strict-origin-when-cross-origin"'])
Out[ ]:

Perceptron



A perceptron models the decision boundary as a linear equation:


$$\omega_0 + \omega_1 x_1 + \omega_2 x_2 + \cdots + \omega_d x_d = 0$$


This means that a perceptron can only create a straight-line (hyperplane) to separate the data. It is instructive to conceptualize a single perceptron as representing a hyperplane in a geometric space.


This identification is worth holding onto, because it is the key to reading everything that follows:


$$ \text{one perceptron} \; = \; \text{one hyperplane} $$


The weights $\omega_1, \cdots, \omega_d$ determine the orientation of that hyperplane and $\omega_0$ determines its offset from the origin. Training a perceptron means moving a single flat surface into position, and nothing a perceptron does can produce a boundary that is not flat.


Once this is accepted, the question that organizes this entire notebook becomes almost obvious. If one perceptron gives one hyperplane, what does a collection of perceptrons give, and what happens when their outputs are fed into further perceptrons?

XOR Problem

Minsky-Papert Controversy on XOR

  • not linearly separable
  • limitation of perceptron

A perceptron is a simple linear classifier that separates data using a single hyperplane. However, there are some classification tasks that a perceptron cannot solve due to its linear nature. A classic example of this limitation is the XOR (exclusive OR) problem.


The XOR function is defined on four points, and its labels can be checked at a glance:


$x_1$ $x_2$ $y$
0 0 0
0 1 1
1 0 1
1 1 0

The two points labelled $1$ occupy one diagonal of the unit square and the two labelled $0$ occupy the other. Any straight line drawn in the plane leaves at least one point on the wrong side, so no choice of $\omega_0, \omega_1, \omega_2$ can classify all four correctly. The perceptron does not fail here because the training algorithm is inadequate; it fails because the model itself cannot express the required boundary.


This is worth emphasizing because it tells us where to look for a remedy. Improving the optimizer, collecting more data, or training for longer cannot help. Only enlarging the family of functions the model can represent can help, and there are exactly two ways to do that.








2. From Perceptron to Multi-Layer Perceptron (MLP)

To explain the concept of a multi-layer perceptron (MLP), we start with a simple perceptron and establish its connection to a single hyperplane. A perceptron functions as a linear classifier by learning a hyperplane that separates data into distinct classes. However, when the data is non-linearly separable, a single hyperplane is insufficient.


In such cases, two possible approaches arise:

  1. Using multiple hyperplanes (or multiple dividers): This can be achieved by stacking multiple layers of perceptrons, allowing the model to combine hyperplanes and approximate more complex decision boundaries.
  2. Learning non-linear boundaries (or curved surface): By introducing nonlinear activation functions and multiple layers, the model can learn non-linear transformations that capture complex patterns in the data.

This progression from a single hyperplane to non-linear boundaries forms the foundation of multi-layer perceptrons (MLPs), enabling them to handle non-linear separability.


The two approaches listed above are worth comparing carefully, because they sound different but converge on the same architecture.


The first says: keep the straight boundaries, but use several of them and combine the results. Two lines can carve the plane into four regions, which is already enough to isolate the XOR diagonals. Combining the outputs requires a second stage that takes the decisions of the first stage as its inputs, and that second stage is another perceptron. Multiple hyperplanes therefore imply multiple layers.


The second says: bend the boundary. A curved surface can separate what no straight one can, and curvature enters through nonlinear activation functions applied between layers.


These are not competing designs. A multi-layer perceptron does both at once: each hidden neuron contributes a hyperplane, and the nonlinear activation combines them into a boundary that is curved. The remainder of this chapter examines the same object from both points of view.


2.1. Perceptron for $h_{\omega}(x)$

  • Neurons compute the weighted sum of their inputs

  • A neuron is activated or fired when the sum $a$ is positive


$$ \begin{align*} a &= \omega_0 + \omega_1 x_1 + \omega_2 x_2 \\ \\ \hat{y} &= g(a) = \begin{cases} 1 & a > 0\\ 0 & \text{otherwise} \end{cases} \end{align*} $$



  • A step function (or sign function) is non-differentiable.
    • To address this limitation, we later replace the step function with differentiable non-linear activation functions such as the sigmoid, or tanh functions.

The reason this matters is worth spelling out. The step function has zero derivative everywhere except at the origin, where the derivative does not exist. Any gradient-based learning rule therefore receives no information from it: the gradient is either zero, which says nothing about which direction to move, or undefined.


For a single perceptron this was tolerable, because the perceptron update rule was derived geometrically rather than by differentiation. For a network of many layers it is fatal. Training a deep model requires propagating derivatives backwards through every layer, and a layer whose derivative is zero blocks that flow completely. Replacing the step function with a smooth activation is therefore not a cosmetic change. It is what makes multi-layer training possible at all.




  • A single layer is insufficient to solve the XOR problem.

2.2. Multi-Layer Perceptron (MLP) = Artificial Neural Networks (ANN)

  • Multi-neurons
    • A single perceptron is often insufficient for complex tasks, as it can only learn a single hyperplane to separate the data. When the data is not linearly separable, a single hyperplane cannot capture the intricate boundaries or patterns in the feature space, necessitating deeper or more complex architectures to model the underlying relationships effectively.
    • For example, if two perceptrons are stacked, it represents two hyperplanes.

Read the figures that follow with the correspondence of Chapter 1 in mind. Each neuron in the hidden layer draws its own hyperplane in the input space, so a hidden layer with three neurons draws three hyperplanes. The output neuron does not see the original inputs at all. It sees only the three numbers reporting where the input fell relative to each of those hyperplanes, and it draws a hyperplane in that new space of three numbers.


This is the whole mechanism, and it is worth stating in one sentence: a hidden layer changes the coordinates in which the problem is expressed, and the output layer solves the problem linearly in the new coordinates.



  • Switch to differentiable activation function (for example, the sigmoid function)


  • In a compact representation
    • Combining the summation and the sigmoid function forms another neuron.

This compaction is why network diagrams are drawn as circles connected by arrows rather than as explicit formulas. A single circle stands for two operations performed in sequence: an affine combination $z = \omega^T x + b$, followed by a nonlinear activation $a = g(z)$. Both are always present, and the distinction between them will matter in Chapter 6, since backpropagation differentiates them separately.



  • Multi-layer perceptron
    • Can be viewed as a sequence of feature extraction layers, where each layer transforms the input into a higher-level, more abstract representation. This process of hierarchical feature extraction allows the MLP to learn complex patterns from data.


Multi-Layer Perceptron as a Sequence of Feature Extraction

A multi-layer perceptron is not merely a collection of neurons but a structured sequence of feature extraction steps. Each layer extracts increasingly abstract and relevant features, allowing the model to handle complex, non-linear tasks that a single-layer perceptron cannot solve.


Intuition Behind Feature Extraction in MLP

  • The first layer extracts basic features from the input (e.g., edges or simple patterns in images).
  • The hidden layers combine these basic features to create more complex representations (e.g., shapes, contours).
  • The final layer maps the abstract features to the output (e.g., predicting a class label or regression value).

By stacking layers, the MLP gradually transforms simple input features into complex features that capture relationships in non-linear data.


Each layer can be thought of as performing a mapping from one feature space to another, progressively refining the representation to make the final decision easier.


This hierarchy of transformations allows the network to move from raw input data to task-specific representations.


A Change in Where the Engineering Happens

This idea deserves more emphasis than it usually receives, because it marks a genuine departure from classical modeling practice.


In traditional machine learning, a practitioner decides which features to compute from the raw data. For a vibration signal one might extract the RMS amplitude, the dominant frequency, and the kurtosis, and then fit a classifier to those three numbers. The quality of the result depends almost entirely on whether the chosen features happened to capture what mattered, and choosing them requires domain expertise that may not exist for a new problem.


A deep network removes that step. The features are parameters, and they are estimated from data by the same optimization that fits the final classifier. The engineer specifies the architecture and the objective; the network determines the representation.


The cost of this convenience should also be recognized. The learned features are optimized for the training objective and nothing else, so they carry no guarantee of physical meaning and cannot generally be interpreted. This is precisely the loss of interpretability discussed in the SciML material, and it is the reason physical constraints are introduced later to restrict what the network is allowed to learn.


2.3. Hidden Layers as Kernel Learning

The XOR example can be solved by pre-processing the data to make the two populations linearly separable.


The construction is small enough to carry out by hand, and doing so makes the abstract claim concrete. Add a third coordinate $x_3 = x_1 x_2$ to each of the four XOR points:


$x_1$ $x_2$ $x_3 = x_1 x_2$ $y$
0 0 0 0
0 1 0 1
1 0 0 1
1 1 1 0

In this three-dimensional space the plane $x_1 + x_2 - 2x_3 = 0.5$ separates the two classes cleanly, as may be verified by substituting the four rows. The data that was not linearly separable in two dimensions has become linearly separable in three.


Nothing about the data changed. Only its representation changed, and the boundary that appears straight in the new coordinates corresponds to a curved boundary in the original ones.



In machine learning, defining the appropriate kernels for non-linear mapping is crucial for models like support vector machines (SVMs) and kernel-based methods. However, in deep learning, the kernel function is effectively represented by the hidden layers of the neural network, which are learned directly from the data rather than being predefined.

  • In traditional machine learning: Kernels (such as polynomial or RBF kernels) are manually chosen to project the input into a higher-dimensional space for better separation.
  • In deep learning: The hidden layers act as adaptive feature extractors, learning the non-linear transformations (analogous to kernels) automatically during training.

This flexibility allows deep learning models, such as multi-layer perceptrons (MLPs), to discover the best representations and transformations for the data without requiring explicit kernel design.


The comparison can be made precise. Both approaches compute a prediction of the form


$$ \hat{y} = \text{(linear model)} \circ \text{(nonlinear feature map)} , $$


and they differ only in where the feature map comes from.


Kernel methods Hidden layers
Feature map chosen by the user in advance learned from the data
Selection requires trial, error, and domain insight performed by gradient descent
Flexibility fixed once chosen adapts to the task
Interpretability the map is known explicitly the map is a set of learned weights

This is why the discussion of kernel selection was left open in the previous notebook. The question of which kernel to use for a given dataset is genuinely difficult, and deep learning does not answer it. It sidesteps it, by making the feature map itself a set of parameters to be estimated. What is gained is generality; what is given up is the explicit knowledge of what transformation is being applied.






3. Logistic Regression in a Form of Neural Network

Here, we will demonstrate a multi-layer perceptron (MLP) using a logistic regression example. In this demonstration, we will observe the hyperplanes and the features learned in the hidden layer at the end of the training process.

  • When we extend logistic regression to a multi-layer perceptron, the hidden layer transforms the input features through non-linear activations, creating intermediate feature spaces.
  • At the end of learning, we can visualize how the hidden layer features and multiple hyperplanes evolve to form more complex decision boundaries.

This approach will illustrate how the MLP builds upon logistic regression by stacking multiple transformations to handle non-linear separability.


It is worth pausing on the fact that logistic regression is already a neural network. It has an input layer, no hidden layer, a single output neuron, a sigmoid activation, and a cross-entropy loss. Everything derived for it in the previous notebook therefore transfers here unchanged.


What follows adds exactly one thing to that picture: a hidden layer between the input and the output. Watching what that single addition does, on data that logistic regression alone cannot handle, is the purpose of this chapter.


Let us start with logistic regression for a linearly separable case.


$$ \begin{align*} y^{(i)} &\in \{0, 1\}\\\\ y &= \sigma (\omega_0 + \omega_1 x_1 + \omega_2 x_2) \end{align*} $$


After training, $\omega_0 + \omega_1 x_1 + \omega_2 x_2 = 0$ will represent a linear classification boundary.



In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import torch

%matplotlib inline
In [ ]:
#training data gerneration
np.random.seed(0)
m = 1000
x1 = 8*np.random.rand(m, 1)
x2 = 7*np.random.rand(m, 1) - 4

g = 0.8*x1 + x2 - 3

C1 = np.where(g >= 0)[0]
C0 = np.where(g < 0)[0]
N = C1.shape[0]
M = C0.shape[0]
m = N + M

X1 = np.hstack([x1[C1], x2[C1]])
X0 = np.hstack([x1[C0], x2[C0]])

train_X = np.vstack([X1, X0])
train_y = np.vstack([np.ones([N,1]), np.zeros([M,1])])

train_X = np.asmatrix(train_X)
train_y = np.asmatrix(train_y)

plt.figure(figsize = (6, 4))
plt.plot(x1[C1], x2[C1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(x1[C0], x2[C0], 'bo', alpha = 0.4, label = 'C0')
plt.legend(loc = 1)
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.show()
No description has been provided for this image
In [ ]:
torch.manual_seed(0)

LogisticRegression = torch.nn.Sequential(
    torch.nn.Linear(2, 1),
    torch.nn.Sigmoid()
)
In [ ]:
optm = torch.optim.Adam(LogisticRegression.parameters(), lr = 0.1)
loss_fn = torch.nn.BCELoss()
In [ ]:
train_x = torch.tensor(np.asarray(train_X), dtype = torch.float32)
train_y_t = torch.tensor(np.asarray(train_y), dtype = torch.float32)

n_iter = 500

for i in range(n_iter):
    y_pred = LogisticRegression(train_x)
    loss = loss_fn(y_pred, train_y_t)

    optm.zero_grad()                    # gradients accumulate, so clear them
    loss.backward()
    optm.step()

    if (i + 1) % 100 == 0:
        print('iteration %4d   loss %.4f' % (i + 1, loss.item()))
iteration  100   loss 0.1104
iteration  200   loss 0.0773
iteration  300   loss 0.0623
iteration  400   loss 0.0534
iteration  500   loss 0.0473
In [ ]:
# .T puts the weights in the (input, output) layout used in the plots below
w = LogisticRegression[0].weight.detach().numpy().T
b = LogisticRegression[0].bias.detach().numpy()

print(w)
print('\n')
print(b)
[[3.5620832]
 [4.5382113]]


[-13.176879]
In [ ]:
x1p = np.arange(0, 8, 0.01).reshape(-1, 1)
x2p = - w[0,0]/w[1,0]*x1p - b[0]/w[1,0]

plt.figure(figsize = (6, 4))
plt.plot(x1[C1], x2[C1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(x1[C0], x2[C0], 'bo', alpha = 0.4, label = 'C0')
plt.plot(x1p, x2p, 'g', linewidth = 3, label = '')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.legend(loc = 1)
plt.show()
No description has been provided for this image

3.1. Looking at Parameters in Nonlinear Classification

Now, let us move on to a non-linearly separable case.


In this scenario, the dataset cannot be separated by a single linear hyperplane. Therefore, we use a multi-layer perceptron (MLP), which approximates complex non-linear decision boundaries by learning multiple hyperplanes and combining them using non-linear activation functions.


Let us see how the MLP performs in this case!


Before that, let us first mention the notation changes commonly used in neural network conventions. In neural network notation, it is common to represent: $\omega_0 \rightarrow b$

  • Weights as $\omega$
  • Bias as $b$ (previously denoted as $\omega_0$)

$$y = \sigma(\omega_0 + \omega_1 x_1 + \omega_2 x_2) \quad \longrightarrow \quad y = \sigma(b + \omega_1 x_1 + \omega_2 x_2)$$


The change is purely notational, but it reflects a practical reality. In the previous notebook the bias was folded into the weight vector by appending an artificial coordinate $x_0 = 1$, which kept the algebra compact. In a network with many layers this becomes awkward: every layer would need its own artificial coordinate inserted and then removed. Deep learning frameworks therefore keep the bias as a separate named quantity, and every layer carries a weight matrix $W$ and a bias vector $b$ as distinct objects. The layer computation is written


$$ z = W x + b, \qquad a = g(z) . $$


This is the form in which the code below, and essentially all deep learning code, is written.



In [ ]:
# training data gerneration
np.random.seed(0)

m = 1000
x1 = 10*np.random.rand(m, 1) - 5
x2 = 8*np.random.rand(m, 1) - 4

g = - 0.5*(x1-1)**2 + 2*x2 + 5

C1 = np.where(g >= 0)[0]
C0 = np.where(g < 0)[0]
N = C1.shape[0]
M = C0.shape[0]
m = N + M

X1 = np.hstack([x1[C1], x2[C1]])
X0 = np.hstack([x1[C0], x2[C0]])

train_X = np.vstack([X1, X0])
train_X = np.asmatrix(train_X)

train_y = np.vstack([np.ones([N,1]), np.zeros([M,1])])

plt.figure(figsize = (6, 4))
plt.plot(x1[C1], x2[C1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(x1[C0], x2[C0], 'bo', alpha = 0.4, label = 'C0')
plt.legend(loc = 1, fontsize = 15)
plt.xlabel(r'$x_1$', fontsize = 15)
plt.ylabel(r'$x_2$', fontsize = 15)
plt.xlim([-5, 5])
plt.ylim([-4, 4])
plt.show()
No description has been provided for this image

As illustrated in the figure, the data is non-linearly distributed. To approximate the non-linear decision boundary using two linear boundaries, a hidden layer with two neurons (plus a bias neuron) is intentionally added.


The word "intentionally" is important. Two hidden neurons were chosen because the boundary visible in the figure looks as though two straight cuts would suffice, and because a hidden layer of size two can be plotted directly, which makes the learned representation visible.


This is not how architectures are normally selected. In practice the number of hidden units is a hyperparameter chosen by experiment, and it is usually larger than strictly necessary. Here the smallest interesting value is chosen deliberately so that the internal state of the network can be examined rather than merely trusted.



In [ ]:
torch.manual_seed(0)

LogisticRegression = torch.nn.Sequential(
    torch.nn.Linear(2, 2),
    torch.nn.Sigmoid(),
    torch.nn.Linear(2, 1),
    torch.nn.Sigmoid()
)
In [ ]:
optm = torch.optim.Adam(LogisticRegression.parameters(), lr = 0.01)
loss_fn = torch.nn.BCELoss()
In [ ]:
train_x = torch.tensor(np.asarray(train_X), dtype = torch.float32)
train_y_t = torch.tensor(np.asarray(train_y), dtype = torch.float32)

n_iter = 1000

for i in range(n_iter):
    y_pred = LogisticRegression(train_x)
    loss = loss_fn(y_pred, train_y_t)

    optm.zero_grad()
    loss.backward()
    optm.step()

    if (i + 1) % 200 == 0:
        print('iteration %4d   loss %.4f' % (i + 1, loss.item()))
iteration  200   loss 0.4306
iteration  400   loss 0.2907
iteration  600   loss 0.1941
iteration  800   loss 0.1468
iteration 1000   loss 0.1223
In [ ]:
w1 = LogisticRegression[0].weight.detach().numpy().T
b1 = LogisticRegression[0].bias.detach().numpy()

w2 = LogisticRegression[2].weight.detach().numpy().T
b2 = LogisticRegression[2].bias.detach().numpy()
In [ ]:
H = train_X*w1 + b1
H = 1/(1 + np.exp(-H))

plt.figure(figsize = (6, 4))
plt.plot(H[0:N,0], H[0:N,1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(H[N:m,0], H[N:m,1], 'bo', alpha = 0.4, label = 'C0')
plt.xlabel('$z_1$', fontsize = 15)
plt.ylabel('$z_2$', fontsize = 15)
plt.legend(loc = 1, fontsize = 15)
plt.axis('equal')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.show()
WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
No description has been provided for this image

Here the features $z_1$ and $z_2$ learned in the hidden layer have two notable characteristics:

  • Bounded Output Range: Both features are restricted to the interval $(0,1)$ due to the use of the sigmoid activation function. This non-linear transformation ensures that the outputs remain within a normalized range, regardless of the input values.

  • Feature Redistribution: The hidden layer redistributes the input data into a new feature space where the transformed data becomes approximately linearly separable. This transformation allows the multi-layer perceptron to form a linear decision boundary in the higher-dimensional space that corresponds to a non-linear boundary in the original input space.


These learned features (or feature extractions) demonstrate how the hidden layer facilitates the approximation of complex patterns by creating new, informative representations of the input data.


Note what $z_1$ and $z_2$ actually are. Each hidden neuron computes $\sigma(b_j + \omega_{1j} x_1 + \omega_{2j} x_2)$, and the sigmoid argument is the signed distance from the input point to that neuron's hyperplane, up to a scale factor. The value of $z_j$ is therefore a smoothed indicator of which side of hyperplane $j$ the point lies on, close to $1$ well on one side, close to $0$ well on the other, and intermediate near the boundary itself.


The coordinates of the new representation are thus answers to the questions "which side of line 1?" and "which side of line 2?". The output layer then draws a straight boundary in the space of those two answers, which is what the next figure confirms.

In [ ]:
x1p = np.arange(0, 1, 0.01).reshape(-1, 1)
x2p = - w2[0,0]/w2[1,0]*x1p - b2[0]/w2[1,0]

plt.figure(figsize = (6, 4))
plt.plot(H[0:N,0], H[0:N,1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(H[N:m,0], H[N:m,1], 'bo', alpha = 0.4, label = 'C0')
plt.plot(x1p, x2p, 'k', linewidth = 3, label = '')
plt.xlabel('$z_1$', fontsize = 15)
plt.ylabel('$z_2$', fontsize = 15)
plt.legend(loc = 1, fontsize = 15)
plt.axis('equal')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.show()
WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
No description has been provided for this image

The result aligns with the linear boundary in the $z$-space, indicating that the transformed features $z_1$ and $z_2$ successfully reshape the input data into a feature space where a linear decision boundary can separate the classes. This transformation validates the role of the hidden layer in projecting non-linearly distributed data into a higher-dimensional space, making it linearly separable.

In [ ]:
x1p = np.arange(-5, 5, 0.01).reshape(-1, 1)
x2p = - w1[0,0]/w1[1,0]*x1p - b1[0]/w1[1,0]
x3p = - w1[0,1]/w1[1,1]*x1p - b1[1]/w1[1,1]

plt.figure(figsize = (6, 4))
plt.plot(x1[C1], x2[C1], 'ro', alpha = 0.4, label = 'C1')
plt.plot(x1[C0], x2[C0], 'bo', alpha = 0.4, label = 'C0')
plt.plot(x1p, x2p, 'k', linewidth = 3, label = '')
plt.plot(x1p, x3p, 'g', linewidth = 3, label = '')
plt.xlabel('$x_1$', fontsize = 15)
plt.ylabel('$x_2$', fontsize = 15)
plt.legend(loc = 1, fontsize = 15)
plt.axis('equal')
plt.xlim([-5, 5])
plt.ylim([-4, 4])
plt.show()
WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
No description has been provided for this image

After training, two linear boundaries in the $x$-space (input space) are learned, corresponding to the linear boundaries formed by the hidden layer neurons. These boundaries result in multiple lines that, when combined, approximate the non-linear classification boundary. This demonstrates how the multi-layer perceptron (MLP) constructs a complex decision boundary in the input space by learning and integrating multiple hyperplanes from the hidden layer.


Taking the three figures together completes the argument of this chapter.


  • In the input space, the hidden neurons place two straight lines. Each is an ordinary perceptron boundary, exactly as in Chapter 1.
  • Passing to $z$-coordinates transforms the data so that the two classes become separable by a single straight line.
  • The output neuron finds that line, which is an ordinary logistic regression in the transformed space.
  • Mapping the result back to the input space produces a curved boundary, formed by combining the two straight ones.

Every individual operation is linear-plus-sigmoid, and nothing more sophisticated occurs anywhere in the network. The nonlinearity of the final boundary arises entirely from composing these simple operations. This is the essential idea of deep learning, and here it can be seen in its entirety on a single page.




4. Regression in a Form of Neural Network

A multi-layer perceptron (MLP) can also be applied to nonlinear regression tasks. Unlike linear regression, which fits a straight line, MLPs can model complex, nonlinear relationships between input features and the output due to their use of multiple layers and nonlinear activation functions.


In this example, the Rectified Linear Unit (ReLU) will be used as the non-linear activation function. Therefore, let us first take a closer look at ReLU before proceeding to the example.


Rectified linear unit (ReLU activation function)


$$h(x) = \max(0, x) = \begin{cases} 0, &\text{if}\;\; x \leq 0 \\ x, &\text{if} \;\; x > 0 \end{cases}$$


Key Characteristics of ReLU:

  • Piecewise Linear: ReLU outputs zero for negative values of $x$ and passes the value of $x$ unchanged for positive values.
  • Sparsity: ReLU introduces sparsity by setting some neuron activations to zero, which can make the network more efficient.
  • Avoids Vanishing Gradients: Unlike sigmoid or tanh, ReLU does not squash values to a narrow range, thus avoiding the issue of vanishing gradients in deep networks.
  • Non-saturating: The gradient remains constant for positive values, making optimization faster.

Why ReLU Replaced the Sigmoid

The sigmoid was the natural choice in Chapter 3, since its output can be read as a probability. Inside the hidden layers of a deep network, however, that interpretation is not needed, and a different property becomes decisive.


Recall that $\sigma'(z) = \sigma(z)(1 - \sigma(z))$, whose maximum value is $0.25$ and which decays toward zero for large $\lvert z \rvert$. Backpropagation multiplies one such factor for every layer it passes through, so after ten layers the gradient has been scaled by at most $0.25^{10}$, roughly $10^{-6}$. The early layers of a deep sigmoid network therefore receive almost no gradient signal and effectively stop learning. This is the vanishing gradient problem, and it is the main reason deep networks were difficult to train before the mid-2000s.


The derivative of ReLU is exactly $1$ for every positive input. Multiplying many such factors together changes nothing, so the gradient reaches the early layers intact.


The Cost

ReLU has a corresponding weakness that should be known. Its derivative is exactly $0$ for negative inputs, so a neuron whose pre-activation is negative for every training example receives no gradient at all and can never recover. Such a unit is said to be dead, and a network can lose a substantial fraction of its capacity this way if the learning rate is too large. Variants such as leaky ReLU, which uses a small nonzero slope for negative inputs, exist to address this.


Note also that ReLU is not differentiable at the origin. As with the perceptron loss in the previous notebook, this causes no practical difficulty, since frameworks simply assign a value there by convention.

$$h(x)$$

In [ ]:
import numpy as np
import torch
import matplotlib.pyplot as plt

def relu(x):
    return np.maximum(0, x)

xp = np.linspace(-1.5, 1.5, 100)
yp = relu(xp)

plt.figure(figsize = (6, 4))
plt.plot(xp, yp, '--', color = 'red', alpha = 0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.grid(alpha = 0.3)
plt.show()
WARNING:matplotlib.axes._base:Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.
No description has been provided for this image

$$h(x-0.5)$$

In [ ]:
xp = np.linspace(-1.5, 1.5, 100)
yp = relu(xp - 0.5)

plt.figure(figsize = (6, 4))
plt.plot(xp, yp, '--', color = 'red', alpha = 0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.grid(alpha = 0.3)
plt.show()
WARNING:matplotlib.axes._base:Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.
No description has been provided for this image

$$h(-2x)$$

In [ ]:
xp = np.linspace(-1.5, 1.5, 100)
yp = relu(-2*xp)

plt.figure(figsize = (6, 4))
plt.plot(xp, yp, '--', color = 'red', alpha = 0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.grid(alpha = 0.3)
plt.show()
WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
No description has been provided for this image

$$- h(-2x-1)$$

In [ ]:
xp = np.linspace(-1.5, 1.5, 100)
yp = -relu(-2*xp - 1)

plt.figure(figsize = (6, 4))
plt.plot(xp, yp, '--', color = 'red', alpha = 0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.grid(alpha = 0.3)
plt.show()
WARNING:matplotlib.axes._base:Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.
No description has been provided for this image

The expression $\omega_2 h(\omega_1 x + b_1) + b_2 $, where $h(\cdot)$ represents the ReLU activation function, can describe a variety of transformations of the ReLU function, including horizontal and vertical shifts, as well as reflections across the $x$-axis or $y$-axis, depending on the values and signs of $\omega_1, \omega_2, b_1$, and $b_2$.


Let us consider the task of modeling a multi-layer perceptron (MLP) to approximate the given non-linear functions.


Before proceeding, it is worth confirming the claim by considering each parameter separately, since these four numbers are the entire vocabulary available to one hidden unit.


  • $b_1$ shifts the kink horizontally, since the elbow occurs where $\omega_1 x + b_1 = 0$, that is, at $x = -b_1/\omega_1$.
  • The sign of $\omega_1$ decides which side of the kink is flat and which side rises.
  • $\omega_2$ scales the slope of the rising part, and its sign reflects the function about the horizontal axis.
  • $b_2$ shifts the whole function vertically.

So a single hidden unit with ReLU produces one kink at a location it chooses, with a slope it chooses, oriented in one of four ways. The examples that follow simply ask what happens when several of these are added together.

In [ ]:
import numpy as np
import torch
import matplotlib.pyplot as plt

np.random.seed(0)
torch.manual_seed(0)

def function(x):
    return x**3 + 0.1*x**2 - x + 0.1

x = np.linspace(-1.5, 1.5, 1000)
y = function(x)

plt.figure(figsize = (6, 4))
plt.plot(x, y, '--', color = 'red', alpha = 0.5)
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha = 0.3)
plt.show()
No description has been provided for this image
In [ ]:
train_x = np.random.uniform(-1.5, 1.5, 1000)
train_y = function(train_x)

Nonlinear Regression as a Linear Combination of ReLU Functions

In an MLP with ReLU activation, the output of the model can be interpreted as a linear combination of ReLU-transformed features. Each hidden layer neuron applies a ReLU transformation to a weighted sum of the input features, and the final output is a linear combination of these transformed features.


Key Insight:


Nonlinear regression with MLPs can be thought of as creating a piecewise linear approximation of the true function, where the linear segments are defined by the ReLU activations in the hidden layer. The network learns how to place these segments and their slopes to best fit the target function.


This picture is worth taking literally, because it makes the behavior of the network in the examples below entirely predictable.


A network with $m$ ReLU units in one hidden layer computes a function with at most $m$ kinks, and is exactly linear between consecutive kinks. Training determines two things: where the kinks are placed, through the hidden weights and biases, and what the slope change is at each kink, through the output weights. Increasing $m$ increases the number of segments available and therefore the resolution of the approximation.


The comparison with classical approximation methods is instructive. Fitting a spline requires the knots to be specified in advance, and choosing them well requires knowing where the function has structure. The network places its own knots, moving them during training toward the regions where they reduce the loss most. This is the same trade seen throughout this notebook: the user gives up explicit control and receives adaptivity in return.


It is also the intuition behind the universal approximation theorem, which states that a single hidden layer with enough units can approximate any continuous function on a bounded interval to any desired accuracy. The theorem guarantees only that such a network exists. It says nothing about how many units are required, and nothing about whether gradient descent will find the right ones.

For a neural network with one hidden layer, the predicted output $ \hat{y} $ is:


$$ \hat{y} = \sum_{j=1}^{m} \omega_j^{(2)} h\left( \sum_{i=1}^{n} \omega_{ij}^{(1)} x_i + b_j^{(1)} \right) + b^{(2)} $$


$\quad$where:

  • $ \omega_{ij}^{(1)} $ and $ b_j^{(1)} $ are the weights and biases for the hidden layer.
  • $ \omega_j^{(2)} $ and $ b^{(2)} $ are the weights and bias for the output layer.

In the following examples, $b^{(2)}$, the bias term for the output layer, is omitted for simplicity. In other words, the vertical shift transformation is not applied.


Read the superscripts as layer indices rather than as exponents. The quantity $\omega_{ij}^{(1)}$ is the weight connecting input $i$ to hidden unit $j$ in the first layer, and $\omega_j^{(2)}$ is the weight from hidden unit $j$ to the output in the second layer. This convention is used consistently in the code that follows and in Chapter 6.


Note the structure of the expression: the sum over $j$ is a linear combination of the hidden activations, and the nonlinearity appears only inside $h(\cdot)$. This is the same form met throughout the previous notebook, a linear model applied to a feature map, except that here the feature map contains adjustable parameters.

Single Neuron with ReLU



$$\hat{y} = \omega^{(2)} \left( h \left( \omega^{(1)} x + b^{(1)} \right) \right)$$


With one hidden unit the model can produce a single kink and nothing more. Whatever the target function looks like, the best this network can do is one flat portion joined to one sloped portion. Comparing the fit obtained below with the target curve shows exactly how much of the shape is out of reach, and provides the baseline against which the next two examples should be read.

In [ ]:
torch.manual_seed(0)

model = torch.nn.Sequential(
    torch.nn.Linear(1, 1),
    torch.nn.ReLU(),
    torch.nn.Linear(1, 1, bias = False)
)

optm = torch.optim.Adam(model.parameters(), lr = 0.01)
loss_fn = torch.nn.MSELoss()

train_x = train_x.reshape(-1, 1)
train_y = train_y.reshape(-1, 1)

train_x_t = torch.tensor(train_x, dtype = torch.float32)
train_y_t = torch.tensor(train_y, dtype = torch.float32)

for i in range(500):
    y_pred = model(train_x_t)
    loss = loss_fn(y_pred, train_y_t)

    optm.zero_grad()
    loss.backward()
    optm.step()

print('final loss %.5f' % loss.item())
final loss 0.35221
In [ ]:
# the optimizer and the loss are set above, together with the model
In [ ]:
# the training loop above has already fitted the model
In [ ]:
weights = model[0].weight.detach().numpy().T
biases = model[0].bias.detach().numpy()
weights2 = model[2].weight.detach().numpy().T

print('Coefficients (Weights):', weights)
print('Intercepts (Biases):', biases)
print('Coefficients (Weights):', weights2)
Coefficients (Weights): [[-0.79910845]]
Intercepts (Biases): [-0.6336963]
Coefficients (Weights): [[-1.1237336]]
In [ ]:
real_x = np.linspace(-1.5, 1.5, 100)
real_y = function(real_x)
In [ ]:
def relu(x):
    return np.maximum(0, x)
In [ ]:
x1p = np.arange(-2, 2, 0.01).reshape(-1, 1)
x2p = weights2*(relu(weights*x1p + biases))

plt.figure(figsize = (6, 4))
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.plot(real_x, real_y, '--', color = 'red', alpha = 0.5)
plt.plot(x1p, x2p, 'c', linewidth = 3)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha = 0.3)
plt.show()
No description has been provided for this image

Two Neurons with ReLU



$$\hat{y} = \omega_1^{(2)} \left( h \left( \omega_1^{(1)} x + b_1^{(1)} \right) \right) + \omega_2^{(2)} \left( h \left( \omega_2^{(1)} x + b_2^{(1)} \right) \right)$$


Two units give two kinks and therefore up to three linear segments. Observe in the result where the network chooses to place them. They will not be evenly spaced, but concentrated where the target function curves most sharply, since that is where an additional kink reduces the loss the most.

In [ ]:
torch.manual_seed(0)

model_2 = torch.nn.Sequential(
    torch.nn.Linear(1, 2),
    torch.nn.ReLU(),
    torch.nn.Linear(2, 1, bias = False)
)

optm = torch.optim.Adam(model_2.parameters(), lr = 0.01)
loss_fn = torch.nn.MSELoss()

train_x = train_x.reshape(-1, 1)
train_y = train_y.reshape(-1, 1)

train_x_t = torch.tensor(train_x, dtype = torch.float32)
train_y_t = torch.tensor(train_y, dtype = torch.float32)

for i in range(3000):
    y_pred = model_2(train_x_t)
    loss = loss_fn(y_pred, train_y_t)

    optm.zero_grad()
    loss.backward()
    optm.step()

print('final loss %.5f' % loss.item())
final loss 0.15099
In [ ]:
weights = model_2[0].weight.detach().numpy().T
biases = model_2[0].bias.detach().numpy()
weights2 = model_2[2].weight.detach().numpy().T

print('Coefficients (Weights):', weights)
print('Intercepts (Biases):', biases)
print('Coefficients (Weights):', weights2)
Coefficients (Weights): [[-0.00748682  1.7940795 ]]
Intercepts (Biases): [-0.82304513 -1.7790197 ]
Coefficients (Weights): [[-0.27234524]
 [ 2.2294097 ]]
In [ ]:
x1p = np.arange(-2, 2, 0.01).reshape(-1, 1)
x2p = weights2[0]*relu(weights[0][0]*x1p + biases[0])
x3p = weights2[1]*relu(weights[0][1]*x1p + biases[1])

plt.figure(figsize = (6, 4))
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.plot(real_x, real_y, '--', color = 'red', alpha = 0.5)
plt.plot(x1p, x2p, 'b', linewidth = 3, alpha = 0.5)
plt.plot(x1p, x3p, 'k', linewidth = 3, alpha = 0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha = 0.3)
plt.show()
No description has been provided for this image
In [ ]:
x1p = np.arange(-2, 2, 0.01).reshape(-1, 1)
x2p = weights2[0]*relu(weights[0][0]*x1p + biases[0]) + weights2[1]*relu(weights[0][1]*x1p + biases[1])

plt.figure(figsize = (6, 4))
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.plot(real_x, real_y, '--', color = 'red', alpha = 0.5)
plt.plot(x1p, x2p, 'c', linewidth = 3)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha = 0.3)
plt.show()
No description has been provided for this image

Four Neurons with ReLU



$$\hat{y} = \omega_1^{(2)} \left( h \left( \omega_1^{(1)} x + b_1^{(1)} \right) \right) + \omega_2^{(2)} \left( h \left( \omega_2^{(1)} x + b_2^{(1)} \right) \right) + \omega_3^{(2)} \left( h \left( \omega_3^{(1)} x + b_3^{(1)} \right)\right) + \omega_4^{(2)} \left( h \left( \omega_4^{(1)} x + b_4^{(1)} \right) \right) $$

In [ ]:
torch.manual_seed(0)

model_4 = torch.nn.Sequential(
    torch.nn.Linear(1, 4),
    torch.nn.ReLU(),
    torch.nn.Linear(4, 1, bias = False)
)

optm = torch.optim.Adam(model_4.parameters(), lr = 0.01)
loss_fn = torch.nn.MSELoss()

train_x = train_x.reshape(-1, 1)
train_y = train_y.reshape(-1, 1)

train_x_t = torch.tensor(train_x, dtype = torch.float32)
train_y_t = torch.tensor(train_y, dtype = torch.float32)

for i in range(4000):
    y_pred = model_4(train_x_t)
    loss = loss_fn(y_pred, train_y_t)

    optm.zero_grad()
    loss.backward()
    optm.step()

print('final loss %.5f' % loss.item())
final loss 0.01411
In [ ]:
weights = model_4[0].weight.detach().numpy().T
biases = model_4[0].bias.detach().numpy()
weights2 = model_4[2].weight.detach().numpy().T
In [ ]:
x1p = np.arange(-2, 2, 0.01).reshape(-1, 1)
x2p = weights2[0]*relu(weights[0][0]*x1p + biases[0])
x3p = weights2[1]*relu(weights[0][1]*x1p + biases[1])
x4p = weights2[2]*relu(weights[0][2]*x1p + biases[2])
x5p = weights2[3]*relu(weights[0][3]*x1p + biases[3])

plt.figure(figsize = (6, 4))
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.plot(real_x, real_y, '--', color = 'red', alpha = 0.5)
plt.plot(x1p, x2p, 'b', linewidth = 3, alpha = 0.5)
plt.plot(x1p, x3p, 'k', linewidth = 3, alpha = 0.5)
plt.plot(x1p, x4p, 'g', linewidth = 3, alpha = 0.5)
plt.plot(x1p, x5p, 'y', linewidth = 3, alpha = 0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha = 0.3)
plt.show()
No description has been provided for this image
In [ ]:
x1p = np.arange(-2, 2, 0.01).reshape(-1, 1)
x2p = weights2[0]*relu(weights[0][0]*x1p + biases[0]) + weights2[1]*relu(weights[0][1]*x1p + biases[1]) + weights2[2]*relu(weights[0][2]*x1p + biases[2]) + weights2[3]*relu(weights[0][3]*x1p + biases[3])

plt.figure(figsize = (6, 4))
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.plot(real_x, real_y, '--', color = 'red', alpha = 0.5)
plt.plot(x1p, x2p, 'c', linewidth = 3)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha = 0.3)
plt.show()
No description has been provided for this image

100 Neurons with ReLU

In [ ]:
torch.manual_seed(0)

model_100 = torch.nn.Sequential(
    torch.nn.Linear(1, 100),
    torch.nn.ReLU(),
    torch.nn.Linear(100, 1, bias = False)
)

optm = torch.optim.Adam(model_100.parameters(), lr = 0.01)
loss_fn = torch.nn.MSELoss()

train_x = train_x.reshape(-1, 1)
train_y = train_y.reshape(-1, 1)

train_x_t = torch.tensor(train_x, dtype = torch.float32)
train_y_t = torch.tensor(train_y, dtype = torch.float32)

for i in range(1000):
    y_pred = model_100(train_x_t)
    loss = loss_fn(y_pred, train_y_t)

    optm.zero_grad()
    loss.backward()
    optm.step()

print('final loss %.5f' % loss.item())
final loss 0.00001
In [ ]:
weights = model_100[0].weight.detach().numpy().T
biases = model_100[0].bias.detach().numpy()
weights2 = model_100[2].weight.detach().numpy().T
In [ ]:
x1p = np.arange(-2, 2, 0.01).reshape(-1, 1)
x2p = 0

for i in range(100):
    x2p += weights2[i]*relu(weights[0][i]*x1p + biases[i])
In [ ]:
plt.figure(figsize = (6, 4))
plt.xlim([-2, 2])
plt.ylim([-1.7, 2.2])
plt.plot(real_x, real_y, '--', color = 'red', alpha = 0.5)
plt.plot(x1p, x2p, 'c', linewidth = 3)
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha = 0.3)
plt.show()
No description has been provided for this image

Key Points:

  • Each curve represents the output of a single ReLU neuron in response to the input feature $x$.
  • The activation values are zero for certain ranges of $x$, indicating that the neuron is "inactive" in those regions.
  • The piecewise linear segments created by the ReLU activations contribute to forming the overall non-linear regression curve.

This illustrates how the hidden layer's neurons, each applying ReLU activation, form different regions of the input space and collectively approximate the non-linear function.


Two observations are worth drawing from this figure.


  • Locality is learned, not built in. Nothing in the architecture assigns a region of the input to a particular neuron. Each unit is connected to the entire input, yet training produces units that are active only over part of the range. The division of labor emerges from the optimization.

  • The representation is redundant. With one hundred units approximating a smooth curve, many neurons contribute nearly identical or overlapping segments, and the same function could be represented by far fewer. Modern networks are routinely over-parameterized in this way. It is a practical observation that such networks are easier to train, even though the resulting model is larger than the task strictly requires.




5. Artificial Neural Networks

So far, we have conducted an in-depth examination of how an artificial neural network (ANN) operates at a micro level, focusing on the underlying mechanisms of individual neurons, layers, and their interactions. Having developed these foundational insights and intuitions, we are now in a position to broaden our perspective and analyze the ANN as a cohesive system, considering its overall architecture, functionality, and the role each component plays in contributing to the network's collective decision-making process.


Complex/Nonlinear Universal Function Approximator

  • ANNs are powerful universal function approximators that can model both linear and nonlinear relationships.

  • By stacking layers of neurons and using non-linear activation functions, ANNs can represent highly complex mappings, making them suitable for a wide range of tasks, such as image recognition, speech processing, and time-series predictions.


ANN Architecture

  • ANNs are typically organized as feedforward networks with layers that are fully connected. In these networks:

    • Each neuron in a layer receives inputs from all neurons in the previous layer.
    • The network propagates information forward, from input to output, without loops.
  • Each layer performs a linear transformation of the input.

  • Linear connections alone cannot capture complex relationships; this is where non-linear activation functions become essential.

  • Each neuron applies a non-linear activation function to the weighted sum of its inputs

  • These nonlinear neurons allow ANNs to learn non-linear decision boundaries, making them capable of solving complex, non-linearly separable problems.


Hidden Layers and Autonomous Feature Learning

  • Hidden layers are the intermediate layers between the input and output layers in a neural network.

  • Each hidden layer learns intermediate representations or features from the input data.

    • Shallow networks have fewer hidden layers and may struggle with highly complex patterns.

    • Deep networks with many hidden layers (Deep Neural Networks) can model highly intricate relationships but require large datasets and longer training times.

  • The neurons in the hidden layers create intermediate transformations that allow the network to construct hierarchical feature representations, leading to improved performance for complex tasks.

  • Unlike traditional machine learning models that rely on manual feature extraction, ANNs can automatically learn features from raw data.

    • Each layer in an ANN learns increasingly abstract features.

    • This autonomous feature learning makes ANNs robust across diverse domains and applications, eliminating the need for domain-specific feature engineering.


Notation for a Layer

It is useful to fix a compact notation before Chapter 6, since backpropagation is unreadable without one. A fully connected layer $l$ transforms its input vector $a^{(l-1)}$ into an output vector $a^{(l)}$ by


$$ z^{(l)} = W^{(l)} a^{(l-1)} + b^{(l)}, \qquad a^{(l)} = g\left(z^{(l)}\right) , $$


where $W^{(l)}$ is a matrix of weights, $b^{(l)}$ is a bias vector, and the activation $g$ is applied elementwise. The input to the network is $a^{(0)} = x$ and the prediction is the activation of the final layer. A network is completely specified by the number of layers, the size of each, the choice of $g$, and the values of all $W^{(l)}$ and $b^{(l)}$.


How Many Parameters?

Counting is a useful habit, and it explains why training is expensive. A layer with $n$ inputs and $m$ outputs has $nm$ weights and $m$ biases. For the MNIST network of Chapter 7, with $784$ inputs and a hidden layer of $100$ units followed by $10$ outputs, the count is


$$ (784 \times 100 + 100) + (100 \times 10 + 10) = 78{,}500 + 1{,}010 = 79{,}510 . $$


Nearly eighty thousand unknowns are estimated from sixty thousand training images. This is worth noticing, because it inverts the usual situation in engineering fitting, where one expects more equations than unknowns. The fact that such a model can still generalize is one of the genuinely surprising features of deep learning, and it is also why regularization and validation are indispensable rather than optional.


Depth Versus Width

The universal approximation theorem says that one sufficiently wide hidden layer suffices in principle, yet practice overwhelmingly favors depth. The reason is efficiency of representation rather than expressive power: functions built by composition, such as the hierarchy of edges to shapes to objects described above, are represented far more compactly by several narrow layers than by one very wide one. Depth reuses intermediate results, in much the same way that dynamic programming reuses subproblem solutions in Section 6.2.








6. ANN Learning

Now that we understand how an Artificial Neural Network (ANN) works, the next step is to determine the unknown parameters (the weights and biases) from the data. Once the structure of the ANN is designed, these parameters must be learned through training. In the following, we will discuss how this learning process takes place.


To begin with the conclusion, the learning process is performed using the backpropagation algorithm. Before diving into the details of backpropagation, it is essential to first review "recursive algorithms" and "dynamic programming", as they form the foundational concepts necessary for understanding backpropagation.


The detour through computer science may appear unmotivated, so it is worth stating the destination in advance.


Gradient descent requires the derivative of the loss with respect to every weight in the network. A network with eighty thousand parameters therefore requires eighty thousand partial derivatives at every single iteration. Computed naively, each of these is a separate application of the chain rule through the whole network, and training would be hopeless.


The saving comes from a structural observation: those eighty thousand derivatives are not independent. They share almost all of their intermediate factors, and if the shared quantities are computed once and stored, each individual derivative costs almost nothing. Recognizing shared subproblems and storing their solutions is exactly what dynamic programming does, and backpropagation is that idea applied to the chain rule. The next two sections develop the idea on small examples where it can be seen clearly.

In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('mGuXwSbantc', width = '560', height = '315', extras = ['referrerpolicy="strict-origin-when-cross-origin"'])
Out[ ]:

6.1. Recursive Algorithm

  • A fundamental concept in computer science:

    • Recursive algorithms play a pivotal role in solving problems by breaking them down into smaller, manageable subproblems.
  • Subproblem-based approach:

    • The solution to the main problem depends on solving smaller instances of the same problem (subproblems).
  • Self-referential function calls:

    • A recursive function calls itself with a modified input until reaching a base case.
    • This mechanism does not directly occur in the physical world, making recursion a unique and abstract computational concept.


In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('t4MSwiqfLaY', width = '560', height = '315', extras = ['referrerpolicy="strict-origin-when-cross-origin"'])
Out[ ]:
  • Factorial example

$$n ! = n \cdot (n-1) \cdots 2 \cdot 1$$

In [ ]:
n = 5

m = 1
for i in range(n):
    m = m*(i+1)

print(m)
120
In [ ]:
def fac(n):
    if n == 1:
        return 1
    else:
        return n*fac(n-1)
In [ ]:
# recursive

fac(5)
Out[ ]:
120

6.2. Dynamic Programming

Dynamic programming (DP) is a powerful algorithmic technique used to solve optimization problems by breaking them down into "overlapping" subproblems. It is particularly effective for problems with overlapping subproblems and optimal substructure, where solutions to smaller subproblems can be reused to build the solution to the overall problem.

  • Optimal Substructure

    • A problem exhibits optimal substructure if the optimal solution to the problem can be constructed from the optimal solutions of its subproblems.
    • Example: In the shortest path problem, the shortest path from point A to point C passing through point B consists of the shortest path from A to B combined with the shortest path from B to C.
  • Overlapping Subproblems

    • A problem has overlapping subproblems if the same subproblems are solved multiple times during the computation.
    • Instead of solving the same subproblem repeatedly, dynamic programming stores the results of solved subproblems in a lookup table (memoization) for future reuse.

Let us discuss dynamic programming using the Fibonacci sequence as an example:


The Fibonacci example is worth examining quantitatively, because the size of the saving is what makes the technique important rather than merely tidy.


Computing $F_n$ by direct recursion calls itself twice at every level, so the number of calls roughly doubles with each increase in $n$, giving a cost that grows like $2^n$. Computing $F_{40}$ this way requires on the order of a billion calls. Yet only forty distinct values are ever involved. The waste is entirely due to recomputing the same quantities: $F_{38}$ is computed twice, $F_{37}$ three times, $F_{36}$ five times, and so on.


Storing each value the first time it is computed reduces the cost from $2^n$ to $n$. The identical situation arises in a neural network, where the gradient with respect to an early weight requires quantities that the gradients of all later weights have already computed.


$$ \begin{align*} F_1 &= F_2 = 1 \\ F_n &= F_{n-1} + F_{n-2} \end{align*} $$

In [ ]:
# naive Fibonacci

def fib(n):
    if n <= 2:
        return 1
    else:
        return fib(n-1) + fib(n-2)
In [ ]:
fib(10)
Out[ ]:
55

The naive recursive method performs redundant calculations, repeatedly solving subproblems even when their solutions have already been computed. In contrast, dynamic programming allocates memory to store previously computed results. When a subproblem is encountered, DP first checks whether the result has already been calculated. If so, it retrieves the value from memory, avoiding redundant computations and significantly reducing the overall computational time. This approach exemplifies the trade-off where computational time is saved at the expense of additional memory usage, an illustration of the space-for-time principle.

In [ ]:
# Memorized DP Fibonacci

def mfib(n):
    global memo

    if memo[n-1] != 0:
        return memo[n-1]
    elif n <= 2:
        memo[n-1] = 1
        return memo[n-1]
    else:
        memo[n-1] = mfib(n-1) + mfib(n-2)
        return memo[n-1]
In [ ]:
import numpy as np

n = 10
memo = np.zeros(n)
mfib(n)
Out[ ]:
np.float64(55.0)

Dynamic Programming (DP) is much faster than a naive recursive algorithm due to its use of memoization to avoid redundant computations.

In [ ]:
n = 30
%timeit fib(30)
148 ms ± 35.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [ ]:
memo = np.zeros(n)
%timeit mfib(30)
351 ns ± 11.6 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

Dynamic Programming (DP) $\approx$ Recursion + Memoization

Dynamic programming can be thought of as an optimization technique that enhances recursion by storing the results of previously computed subproblems (memoization). This avoids redundant calculations and improves efficiency.

  • Recursion: Breaks a problem down into smaller subproblems and solves them recursively.
  • Memoization: Stores the solutions to subproblems in a cache (for example, a dictionary or an array) so that they can be reused instead of recomputed.

This combination allows DP to solve problems with overlapping subproblems and optimal substructure more efficiently.


6.3. Training Neural Networks

Now, let us return to our original problem: learning the weights and biases of a multi-layer perceptron (MLP).


  • Learning or estimating weights and biases of multi-layer perceptron from training data
    1. Forward Pass: Input $X$ passes through the network to generate $\hat Y$.
    2. Loss Calculation: The loss function compares $\hat Y$ and $Y$.
    3. Backward Pass: The optimizer computes gradients of the loss concerning the weights.
    4. Weight Update: Weights are adjusted based on the gradients.
    5. Repeat: The process continues until the loss converges or the training reaches a predefined number of iterations.

These five steps constitute one iteration of training, and every deep learning framework implements exactly this loop. Two terms attached to it should be defined now, since they appear in the code of Chapter 7.


  • A batch is the subset of training examples used in a single iteration. Using all $m$ examples every time is called batch gradient descent and is accurate but slow; using one example at a time is stochastic gradient descent and is fast but noisy. Practice uses mini-batches of perhaps $32$ to $256$ examples, which is a compromise that also happens to suit the parallel structure of a GPU.
  • An epoch is one complete pass through the training set. Training for ten epochs means every example has been visited ten times.

Note also that the noise introduced by mini-batching is not purely a nuisance. The stochastic fluctuation helps the optimizer escape shallow local minima, which matters because, unlike logistic regression, the loss surface of a neural network is not convex.



Loss Function

  • Measures error between target values and predictions

$$ \min_{\omega} \sum_{i=1}^{m}\ell\left( h_{\omega}\left(x^{(i)}\right),y^{(i)}\right)$$


  • Example
    • Squared loss (for regression):

      $$ \frac{1}{m} \sum_{i=1}^{m} \left(h_{\omega}\left(x^{(i)}\right) - y^{(i)}\right)^2 $$

    • Cross-entropy (for classification):

      $$ -\frac{1}{m}\sum_{i=1}^{m} \left[ y^{(i)}\log\left(h_{\omega}\left(x^{(i)}\right)\right) + \left(1-y^{(i)}\right)\log\left(1-h_{\omega}\left(x^{(i)}\right)\right) \right]$$

Learning

  • Learning weights and biases from data using gradient descent

$$\omega \leftarrow \omega - \alpha \nabla_{\omega} \ell \left( h_{\omega}\left(x^{(i)}\right), y^{(i)} \right)$$


  • $\frac{\partial \ell}{\partial \omega}$: too many computations are required for all $\omega$

  • Structural constraints of NN:

    • Composition of functions
    • Chain rule
    • Dynamic programming

These three lines are the entire argument of this section in compressed form, and they are worth unpacking.


A network is a composition of functions: the loss depends on the output, which depends on the last hidden layer, which depends on the one before it, and so on back to the input. Differentiating a composition requires the chain rule, which multiplies the derivatives of the individual stages. Because every weight in the network sits somewhere inside this same composition, the chain rule expressions for different weights share long common factors. Computing those shared factors once and reusing them is dynamic programming.


The name of the resulting algorithm follows from the order in which the shared factors become available. They are computed starting at the loss and moving backwards toward the input, which is why the method is called backpropagation.

Backpropagation

Backpropagation is a fundamental algorithm used for training neural networks by adjusting their weights to minimize the loss function. It operates by propagating errors backward through the network, from the output layer to the input layer, using the chain rule of differentiation.


  • Forward propagation

    • the initial information propagates up to the hidden units at each layer and finally produces output
  • Backpropagation

    • allows the information from the cost to flow backwards through the network in order to compute the gradients
  • Chain rule

    • Computing the derivative of the composition of functions

$$ \begin{align*} f(g(x))' &= f'(g(x))g'(x)\\\\ {dz \over dx} &= {dz \over dy} \cdot {dy \over dx}\\ {dz \over dw} &= \left({dz \over dy} \cdot {dy \over dx}\right) \cdot {dx \over dw}\\ {dz \over du} &= \left({dz \over dy} \cdot {dy \over dx} \cdot {dx \over dw}\right) \cdot {dw \over du} \end{align*} $$


  • Backpropagation
    • Update weights recursively with memory

A Worked Example

The mechanism is easiest to see on a network small enough to differentiate by hand. Consider a chain of three scalar layers with no bias terms,


$$ a_1 = g(\omega_1 x), \qquad a_2 = g(\omega_2 a_1), \qquad \hat{y} = \omega_3 a_2, \qquad \ell = (\hat{y} - y)^2 .$$


The three derivatives required for gradient descent are


$$ \begin{align*} \frac{\partial \ell}{\partial \omega_3} &= \underbrace{\frac{\partial \ell}{\partial \hat{y}}}_{\delta_3} \; a_2 \\\\ \frac{\partial \ell}{\partial \omega_2} &= \underbrace{\frac{\partial \ell}{\partial \hat{y}} \, \omega_3 \, g'(\omega_2 a_1)}_{\delta_2} \; a_1 \\\\ \frac{\partial \ell}{\partial \omega_1} &= \underbrace{\frac{\partial \ell}{\partial \hat{y}} \, \omega_3 \, g'(\omega_2 a_1) \, \omega_2 \, g'(\omega_1 x)}_{\delta_1} \; x \end{align*} $$


Look at the underbraced factors. Each one contains the entire preceding one: $\delta_2$ is $\delta_3$ multiplied by two new terms, and $\delta_1$ is $\delta_2$ multiplied by two more. Computing them independently would repeat the same multiplications over and over, which is precisely the redundancy that the naive Fibonacci recursion exhibited in Section 6.2.


Computing them in the order $\delta_3$, then $\delta_2$, then $\delta_1$, and storing each as it is produced, makes every step cost only two additional multiplications. That ordering runs from the output back toward the input, which is what gives backpropagation its name and its direction.


Two further observations generalize to networks of any size.


  • Every gradient has the form $\dfrac{\partial \ell}{\partial \omega} = \delta \cdot (\text{input to that weight})$. This is the same "error times input" structure met for the perceptron and for logistic regression in the previous notebook. Backpropagation does not introduce a new kind of update rule; it supplies the correct $\delta$ for a weight buried deep inside a composition.

  • The activations $a_1$ and $a_2$ appear in the gradients, so they must be retained from the forward pass. This is why training a network requires substantially more memory than merely evaluating one, and why memory rather than arithmetic is often the binding constraint on batch size.



Learning the weights and biases in an artificial neural network (ANN) can be viewed as an optimization process using gradient descent or its variants. However, a significant challenge arises due to the large number of parameters (or unknowns) that need to be estimated, as well as the extensive computations of gradients required during each iteration.


The uniqueness of this problem lies in the sequentially stacked structure of ANN layers, which can be interpreted as a composition of functions. As a result, the computation of gradients involves the application of the chain rule, leading to repeated calculations for intermediate layers.


By computing these gradients in a backward manner, previously computed results can be reused in subsequent steps. This is the fundamental concept behind backpropagation, making it analogous to dynamic programming. The term 'back' in 'backpropagation' highlights this backward computation process, where intermediate results are stored and revisited to avoid redundant computations.


The magnitude of the saving is worth stating. A naive approach that evaluates each of the $P$ partial derivatives separately costs roughly $P$ times the cost of one forward pass. Backpropagation obtains all $P$ of them in a single backward pass whose cost is comparable to one forward pass.


For the modest MNIST network of Chapter 7 this is the difference between roughly eighty thousand network evaluations per iteration and approximately two. No modern neural network could be trained without it, which is why a paper about an application of the chain rule is regarded as one of the foundational results of the field.


One clarification is worth making, since the terms are often confused. Backpropagation is not an optimization algorithm. It computes gradients. Gradient descent, or a variant such as Adam, is the optimization algorithm that consumes those gradients and updates the weights. The two are separate components that happen always to be used together.

Summary

  • Learning weights and biases from data using gradient descent $\approx$ backpropagation


  • It is not easy to numerically compute gradients in network in general.
    • The good news: people have already done all the "hardwork" of developing numerical solvers (or libraries)
    • There are a wide range of tools: PyTorch, TensorFlow, JAX

It should be emphasized that this convenience is not an approximation. The frameworks do not estimate derivatives by finite differences, which would require a step size and would introduce truncation error. They apply the chain rule symbolically to the recorded sequence of elementary operations, producing gradients that are exact to floating-point precision. This is automatic differentiation, and backpropagation is its application to the layered structure of a neural network.


The practical consequence for the remainder of this course is significant: once a model and a loss are written down, the derivatives are obtained without further effort. Attention shifts from deriving gradients to designing architectures and objectives.


6.4. Historical Notes

Backpropagation and Geoffrey Hinton

Backpropagation (short for "backward propagation of errors") is a core algorithm used to train artificial neural networks (ANNs). It calculates gradients of the loss function with respect to the network's weights by applying the chain rule of calculus in a systematic, layer-by-layer manner from the output layer back to the input layer. This enables the efficient updating of weights using optimization methods such as gradient descent.

While the mathematical foundations of backpropagation date back to the 1960s (notably by Werbos in 1974, who applied it to neural networks), it was largely unrecognized until Geoffrey Hinton, along with David Rumelhart and Ronald J. Williams, demonstrated its practical potential for deep learning in their seminal 1986 paper, "Learning Representations by Back-Propagating Errors".


Why Hinton's Contribution Was Pivotal

Revival of Neural Networks:

  • During the 1980s, symbolic AI dominated the field, and neural networks were seen as outdated and impractical. Hinton's work with backpropagation demonstrated that neural networks could learn complex representations and outperformed traditional AI approaches on several tasks.

Scalability and Deep Networks:

  • Hinton's research showed that backpropagation could handle multi-layer (deep) networks, allowing networks to automatically learn hierarchical representations of features.

The Legacy of Backpropagation:

  • The 1986 breakthrough laid the groundwork for modern deep learning.
  • Backpropagation became the standard training algorithm for neural networks.
  • Hinton continued to push the field forward, contributing to the resurgence of neural networks in the 2000s and co-developing groundbreaking architectures like deep belief networks (DBNs) and capsule networks.

The Broader Impact on Modern AI

Backpropagation remains essential in modern deep learning frameworks (such as TensorFlow and PyTorch), enabling the training of massive neural networks for applications in language models (e.g., GPT), computer vision, and beyond. Hinton's contributions earned him the nickname "the godfather of deep learning" and prestigious accolades, including the Turing Award (2018) and the Nobel Prize in Physics (2024).




7. ANN with MNIST

In [ ]:
from IPython.display import YouTubeVideo
YouTubeVideo('z-ZhKdQpF7I', width = '560', height = '315', extras = ['referrerpolicy="strict-origin-when-cross-origin"'])
Out[ ]:

7.1. What is MNIST?

From Wikipedia

  • The MNIST database (Mixed National Institute of Standards and Technology database) is a large database of handwritten digits that is commonly used for training various image processing systems. The database is also widely used for training and testing in the field of machine learning. It was created by "re-mixing" the samples from NIST's original datasets. The creators felt that since NIST's training dataset was taken from American Census Bureau employees, while the testing dataset was taken from American high school students, NIST's complete dataset was too hard.

  • MNIST (Mixed National Institute of Standards and Technology database) database

    • Handwritten digit database
    • $28 \times 28$ gray scaled image
    • (Flattened matrix into a vector of $28 \times 28 = 784$) $\rightarrow$ done by a Flatten layer inside the model




More here


We will use the MNIST dataset to create a multinomial classifier that can detect whether the image belongs to one of the classes: 0, 1, 2, 3, 4, 5, 6, 7, 8, or 9. In short, we are teaching a computer to recognize handwritten digits.

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import torch
from torchvision import datasets

%matplotlib inline

Let us download and load the dataset.

In [ ]:
train_set = datasets.MNIST(root = './data', train = True, download = True)
test_set = datasets.MNIST(root = './data', train = False, download = True)

train_x = train_set.data.numpy()/255.0
train_y = train_set.targets.numpy()

test_x = test_set.data.numpy()/255.0
test_y = test_set.targets.numpy()
100%|██████████| 9.91M/9.91M [00:00<00:00, 15.2MB/s]
100%|██████████| 28.9k/28.9k [00:00<00:00, 433kB/s]
100%|██████████| 1.65M/1.65M [00:00<00:00, 4.25MB/s]
100%|██████████| 4.54k/4.54k [00:00<00:00, 4.02MB/s]
In [ ]:
print ('The training data set is:\n')
print (train_x.shape)
print (train_y.shape)
The training data set is:

(60000, 28, 28)
(60000,)
In [ ]:
print ('The test data set is:')
print (test_x.shape)
print (test_y.shape)
The test data set is:
(10000, 28, 28)
(10000,)

Display a few random samples from it:

In [ ]:
# So now we have a 28x28 matrix, where each element is an intensity level from 0 to 1.

img = train_x[5]
img.shape
Out[ ]:
(28, 28)

Let us visualize what some of these images and their corresponding training labels look like.

In [ ]:
plt.figure(figsize = (4, 4))
plt.imshow(img, 'gray')
plt.xticks([])
plt.yticks([])
plt.show()
No description has been provided for this image
In [ ]:
print(train_y[5])
2

7.2. ANN with PyTorch

  • Feed a gray image to ANN

A fully connected network expects a vector, while an MNIST image is a $28 \times 28$ array, so each image is flattened into a vector of length $784$ by reading its pixels in order.


It is worth noticing what this discards. Flattening destroys the spatial relationship between neighboring pixels: two pixels that are vertically adjacent in the image end up $28$ positions apart in the vector, and the network has no way of knowing they were neighbors. The network must therefore learn from scratch, from data alone, a fact about the geometry of images that was available for free before flattening.


This is the specific weakness that convolutional networks remove, and it is the reason they outperform fully connected networks on images. The point is raised again at the end of the chapter.



  • Our network model
    • The hyperparameters of the network have not been optimized.
    • This example is intended solely to demonstrate how to implement the network in Python based on the provided structure.


  • Network training (learning)

    • Forward Pass: Input $X$ passes through the network to generate $\hat Y$.

    • Loss Calculation: The loss function compares $\hat Y$ and $Y$.

    • Backward Pass: The optimizer computes gradients of the loss concerning the weights.

    • Weight Update: Weights are adjusted based on the gradients.

    • Repeat: The process continues until the loss converges or the training reaches a predefined number of iterations.


$$\omega := \omega - \alpha \nabla_{\omega} \, \ell \left( h_{\omega} \left(x^{(i)}\right),\, y^{(i)}\right)$$

7.2.1. Import Library

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import torch
from torchvision import datasets

7.2.2. Load MNIST Data

  • Download MNIST data with torchvision.datasets
In [ ]:
train_set = datasets.MNIST(root = './data', train = True, download = True)
test_set = datasets.MNIST(root = './data', train = False, download = True)

train_x = train_set.data.numpy()/255.0
train_y = train_set.targets.numpy()

test_x = test_set.data.numpy()/255.0
test_y = test_set.targets.numpy()

7.2.3. Define an ANN Structure

  • Input size
  • Hidden layer size
  • The number of classes

Three design decisions are recorded here, and only one of them is genuinely free.


  • The input size is fixed by the data at $28 \times 28 = 784$.
  • The number of classes is fixed by the task at $10$, one output unit per digit.
  • The hidden layer size is a hyperparameter. It has no correct value determined by the problem, and it is chosen by experiment. Too few units and the network cannot represent the required function; too many and it fits the training data more closely than is warranted, while costing more to train.

The Output Layer

The output layer differs from the binary case of the previous notebook and deserves attention. With ten mutually exclusive classes, the network produces ten scores and converts them into a probability distribution using the softmax function


$$ \hat{y}_k = \frac{e^{z_k}}{\sum_{j=1}^{K} e^{z_j}} .$$


Exponentiating makes every value positive, and dividing by the sum forces the ten outputs to add to one, so they can be read as probabilities over the ten digits. Softmax is the multiclass generalization of the sigmoid, and for $K = 2$ the two coincide.


The labels are correspondingly represented as one-hot vectors: the digit $3$ becomes a vector of ten entries that is $1$ in position three and $0$ elsewhere. The cross-entropy loss defined in the next section then compares two probability distributions over the same ten outcomes.



In [ ]:
model = torch.nn.Sequential(
    torch.nn.Flatten(),
    torch.nn.Linear(28*28, 100),
    torch.nn.ReLU(),
    torch.nn.Linear(100, 10)          # logits; the softmax is inside the loss
)

7.2.4. Define Loss and Optimizer


Optimizer

  • This defines how the model is updated based on the data it sees and its loss function.
  • torch.optim.Adam: the most popular optimizer
  • Plain gradient descent uses a single learning rate for every parameter and every iteration. Adam maintains a running estimate of the recent magnitude of each parameter's gradient and scales the step accordingly, so that parameters with consistently small gradients still make progress. It also carries momentum, an averaging of recent gradients that damps oscillation. It is not guaranteed to be better than well-tuned gradient descent, but it is far less sensitive to the choice of learning rate, which is why it is the usual default.

Loss

  • In PyTorch the softmax is not placed in the model. torch.nn.CrossEntropyLoss takes the raw scores of the last layer, the logits, and applies the softmax and the logarithm together in one numerically stable step, exactly as BCEWithLogitsLoss did for the binary case in the previous notebook. The model therefore ends with a plain linear layer of ten units.
  • This defines how the accuracy of the model is measured during training. Training minimizes this function, which is what steers the model in the right direction.
  • Classification: Cross-entropy
    • The multiclass generalization of the logistic regression loss

$$ -\frac{1}{m}\sum_{i=1}^{m} \sum_{k=1}^{K} y_k^{(i)}\log\left(\hat{y}_k^{(i)}\right) $$


  • Here $K = 10$ is the number of classes, $y_k^{(i)}$ is $1$ if example $i$ belongs to class $k$ and $0$ otherwise, and $\hat{y}_k^{(i)}$ is the predicted probability of class $k$. For $K = 2$ this reduces to the binary cross-entropy used for logistic regression.

Metrics

  • Metrics are quantities monitored during training and testing. In this example we track the accuracy, that is, the fraction of images that are correctly classified.
In [ ]:
optm = torch.optim.Adam(model.parameters(), lr = 0.001)
loss_fn = torch.nn.CrossEntropyLoss()     # softmax and log in one stable step

7.2.5. Train Model



In [ ]:
from torch.utils.data import TensorDataset, DataLoader

# 1. Convert to Tensors
train_x_t = torch.tensor(train_x, dtype = torch.float32)
train_y_t = torch.tensor(train_y, dtype = torch.long)

# 2. Define Dataset and DataLoader (Standard Way)
dataset = TensorDataset(train_x_t, train_y_t)
train_loader = DataLoader(dataset, batch_size = 32, shuffle = True)

# 3. Model Training Loop
for epoch in range(5):
    running_loss = 0.0

    # DataLoader handles shuffling and batching automatically
    for batch_x, batch_y in train_loader:
        y_pred = model(batch_x)
        loss = loss_fn(y_pred, batch_y)

        optm.zero_grad()
        loss.backward()
        optm.step()

        running_loss += loss.item()

    # Tip: It is standard practice to output the average loss for the entire epoch rather than just the last batch's loss.
    epoch_loss = running_loss / len(train_loader)
    print('epoch %s   loss %.4f' % (epoch + 1, epoch_loss))
epoch 1   loss 0.3048
epoch 2   loss 0.1409
epoch 3   loss 0.0981
epoch 4   loss 0.0754
epoch 5   loss 0.0593

7.2.6. Test or Evaluate


The reason for holding out a test set should be stated explicitly, since it is the single most important discipline in applied machine learning.


A network with roughly eighty thousand parameters can achieve high accuracy on sixty thousand training images partly by learning genuine structure and partly by memorizing individual examples. Training accuracy therefore measures a mixture of the two and cannot distinguish them. The test set consists of images the optimizer never saw, so accuracy on it reflects only what generalizes.


A gap between the two numbers is the signature of overfitting. If training accuracy is high while test accuracy is much lower, the network has memorized more than it has learned, and the standard responses are to collect more data, reduce the capacity of the model, or add regularization such as weight decay or dropout.

In [ ]:
test_x_t = torch.tensor(test_x, dtype = torch.float32)
test_y_t = torch.tensor(test_y, dtype = torch.long)

# Key addition: Switch to evaluation mode
model.eval()

with torch.no_grad():
    y_pred = model(test_x_t)
    test_loss = loss_fn(y_pred, test_y_t).item()
    test_acc = (y_pred.argmax(dim = 1) == test_y_t).float().mean().item()

print('test loss %.4f   test accuracy %.4f' % (test_loss, test_acc))
test loss 0.0946   test accuracy 0.9695
In [ ]:
test_img = test_x[np.random.choice(test_x.shape[0], 1)]

model.eval()
with torch.no_grad():
    logits = model(torch.tensor(test_img, dtype = torch.float32))
    predict = torch.softmax(logits, dim = 1).numpy()

mypred = np.argmax(predict, axis = 1)

plt.figure(figsize = (8,4))
plt.subplot(1,2,1)
plt.imshow(test_img.reshape(28, 28), 'gray')
plt.axis('off')
plt.subplot(1,2,2)
plt.stem(predict[0])
plt.show()

print('Prediction : {}'.format(mypred[0]))
No description has been provided for this image
Prediction : 6

What is the highest accuracy you can achieve with this first fully connected model? Since the handwritten digit classification task is pretty straightforward, you may be wondering how we can do better...


$\Rightarrow$ Convolutional neural networks (CNNs) are particularly well suited for a variety of tasks in computer vision, and have achieved near-perfect accuracies on the MNIST dataset. We will build a CNN and ultimately output a probability distribution over the 10 digit classes (0-9) in the next lectures.


7.3. Historical Notes

The MNIST Dataset and Yann LeCun

The MNIST (Modified National Institute of Standards and Technology) dataset is one of the most famous benchmark datasets in machine learning, consisting of 60,000 training images and 10,000 test images of handwritten digits (0 to 9), each in a 28x28 grayscale pixel format. It is widely used for training and evaluating classification models, especially in the field of neural networks.


Yann LeCun's Role in the Creation of MNIST

Yann LeCun (Turing Award, 2018), a pioneer in the field of deep learning, created the MNIST dataset in collaboration with Corinna Cortes and Christopher J.C. Burges in 1998. LeCun, known for his groundbreaking work in neural networks, particularly convolutional neural networks (CNNs), developed MNIST as an easily accessible and standardized dataset to evaluate machine learning algorithms.

Before MNIST, many machine learning models struggled with inconsistent and poorly standardized datasets for real-world pattern recognition tasks. MNIST solved this by providing a well-curated, balanced, and preprocessed dataset that became the de facto benchmark for evaluating classification models.


Legacy of MNIST

Despite its simplicity, MNIST remains a crucial stepping stone in the history of artificial intelligence. It demonstrated the power of neural networks in pattern recognition and inspired the development of modern deep learning architectures. Yann LeCun's creation of MNIST not only provided a foundational benchmark for machine learning but also played a pivotal role in popularizing neural networks and shaping the trajectory of modern AI research.




Summary

Key takeaways

  • A perceptron produces a single linear boundary, so a problem that is not linearly separable is outside the family of functions it can represent. Enlarging the family — not training harder — is the only remedy.
  • Stacking perceptrons supplies that larger family. A hidden layer maps the inputs into new coordinates $z$, and the output layer draws a straight line in those coordinates; plotted in $z$, classes that were curved apart in $x$ become linearly separable. The hidden units are feature extractors, and their features are learned rather than designed.
  • With a ReLU activation, a network with one hidden layer of $m$ units computes a sum of $m$ hinge functions, which is a piecewise linear function with $m$ kinks. One unit gives a single kink, two give two segments, four begin to trace the curve, and a hundred are indistinguishable from it. This is the universal approximation property in its most concrete form.
  • A general artificial neural network is a composition of layers, each an affine map followed by a nonlinearity. Specifying a model means stating the architecture and the loss; everything else follows.
  • Backpropagation is the chain rule applied recursively to that composition, reusing the partial results of the layer above exactly as dynamic programming reuses subproblems. It computes gradients; it is not an optimization algorithm. Gradient descent, or a variant such as Adam, is what uses them.
  • In PyTorch a model is torch.nn.Sequential, and one training step is four lines: forward, zero_grad(), backward(), step(). Mini-batching makes each step cheap, and one pass through the data is an epoch; the noise it introduces also helps the optimizer escape shallow minima.
  • For multi-class classification the output layer has one unit per class and the loss is the cross-entropy. In PyTorch the softmax is not placed in the model: CrossEntropyLoss takes the logits and applies the softmax and logarithm together, for the same numerical reason that BCEWithLogitsLoss did in the binary case.
  • Flattening an image into a vector of $784$ pixels discards the fact that neighbouring pixels are neighbours, and the network has to learn that geometry from data. Test accuracy, measured on images the optimizer never saw, is what reports how much of the training accuracy generalizes.

Looking ahead. The weakness identified at the end — that a fully connected network is told nothing about the spatial structure of an image — is what convolutional networks repair, by sharing weights across positions and looking at small neighbourhoods at a time. The same idea of building known structure into the architecture, rather than hoping the network discovers it from data, is what physics-informed models do with governing equations.