6 - The perceptron and multilayer neural networks, from scratch

In chapter 6 the Minermont team meets neural networks at AIA. Hazel starts with the simplest one, Rosenblatt's perceptron, and with its famous limitation, and then shows how stacking neurons into layers and training them with backpropagation overcomes it. In this notebook we build all of it with NumPy, without any deep learning library, as in Emma's "level 0" (section 6.8):

1. Import the libraries.
2. The perceptron and its learning rule.
3. The XOR problem.
4. A multilayer network that solves XOR, trained with backpropagation.
5. Activation functions.
6. Martha's complete example: one forward pass, one backward pass and one update (section 6.7).
7. Checking backpropagation with numerical derivatives.

Step 1: Import the libraries

---

import numpy as np
import matplotlib.pyplot as plt

np.set_printoptions(precision=3, suppress=True)

Step 2: The perceptron

---

A perceptron computes a weighted sum of its inputs plus a bias and passes it through a step function:

$$z = w \cdot x + b \qquad\qquad \hat{y} = \begin{cases} 1 & \text{if } z \ge 0 \\ 0 & \text{if } z < 0 \end{cases}$$

As Martha notices in the book, $z$ is the same formula as linear regression, and the step is an abrupt version of the sigmoid of logistic regression. Rosenblatt's learning rule goes through the examples one by one and, only when the prediction is wrong, moves the weights towards the correct answer:

$$w \leftarrow w + \eta\,(y - \hat{y})\,x \qquad\qquad b \leftarrow b + \eta\,(y - \hat{y})$$

We try it on the logical AND and OR functions, whose truth tables are our four training examples:

X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]], dtype=float)
targets = {'AND': np.array([0, 0, 0, 1]), 'OR': np.array([0, 1, 1, 1]), 'XOR': np.array([0, 1, 1, 0])}

def train_perceptron(X, y, learning_rate=0.1, epochs=20):
    w, b, errors = np.zeros(X.shape[1]), 0.0, []
    for _ in range(epochs):
        mistakes = 0
        for xi, yi in zip(X, y):
            y_hat = int(w @ xi + b >= 0)
            update = learning_rate * (yi - y_hat)       # zero when the prediction is right
            w, b = w + update * xi, b + update
            mistakes += int(update != 0)
        errors.append(mistakes)
    return w, b, errors

def predict_perceptron(X, w, b):
    return (X @ w + b >= 0).astype(int)

perceptrons = {}
for gate in ['AND', 'OR']:
    w, b, errors = train_perceptron(X, targets[gate])
    perceptrons[gate] = (w, b)
    print(f"{gate}: w = {w}, b = {b:.1f}, predictions {predict_perceptron(X, w, b)}, mistakes per epoch {errors[:6]}")

AND: w = [0.2 0.1], b = -0.2, predictions [0 0 0 1], mistakes per epoch [2, 3, 3, 0, 0, 0]
OR: w = [0.1 0.1], b = -0.1, predictions [0 1 1 1], mistakes per epoch [2, 2, 1, 0, 0, 0]

Both functions are learned in a few passes: from the fourth epoch on there are no more mistakes, and the weights define a line that separates the 1s from the 0s.

Step 3: The XOR problem

---

XOR is Hazel's light bulb with two switches: it turns on when exactly one switch is on, and off when both are in the same position. Let's train the same perceptron on it:

w_xor, b_xor, errors_xor = train_perceptron(X, targets['XOR'], epochs=20)
print(f"XOR: predictions {predict_perceptron(X, w_xor, b_xor)} (target {targets['XOR']})")
print(f"Mistakes per epoch: {errors_xor}")

fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
for ax, gate in zip(axes, ['AND', 'OR', 'XOR']):
    y = targets[gate]
    ax.scatter(*X[y == 0].T, marker='^', s=150, color='tab:red', label='0 (off)')
    ax.scatter(*X[y == 1].T, marker='o', s=150, color='tab:blue', label='1 (on)')
    if gate in perceptrons:
        w, b = perceptrons[gate]
        xs = np.linspace(-0.5, 1.5, 2)
        ax.plot(xs, -(w[0] * xs + b) / w[1], 'k--', label='perceptron boundary')
        ax.set_title(f'{gate}: one straight line is enough')
    else:
        ax.set_title('XOR: no straight line separates the classes')
    ax.set_xlim(-0.5, 1.5)
    ax.set_ylim(-0.5, 1.5)
    ax.set_xlabel('$x_1$ (switch 1)')
    ax.set_ylabel('$x_2$ (switch 2)')
    ax.legend(loc='upper right', fontsize=8)
plt.show()

XOR: predictions [1 1 0 0] (target [0 1 1 0])
Mistakes per epoch: [3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]

The perceptron never stops making mistakes: from the third epoch on it corrects its weights on all four examples in every epoch, and they go round in circles without settling. The plot shows why: the two "on" points are at opposite corners, and so are the two "off" points, so no straight line can leave each pair on its own side. It is the limitation that Minsky and Papert pointed out in 1969.

Step 4: A multilayer network for XOR

---

The solution is to add a hidden layer with a nonlinear activation. Our network has 2 inputs, 2 hidden neurons with the hyperbolic tangent, and 1 output neuron with the sigmoid, which gives the probability that the bulb is on:

$$a^{(1)} = \tanh\!\left(W^{(1)} x + b^{(1)}\right) \qquad\qquad \hat{y} = \sigma\!\left(W^{(2)} a^{(1)} + b^{(2)}\right)$$

The cost is the log loss (binary cross-entropy) of chapter 3. To train it we need the gradient of the cost with respect to every weight, and that is what backpropagation computes: the chain rule, applied from the output back to the input. For this network:

• Output error: $\delta^{(2)} = \hat{y} - y$ (the same simple expression as in logistic regression).
• Hidden error: $\delta^{(1)} = \left((W^{(2)})^\top \delta^{(2)}\right) \odot \left(1 - (a^{(1)})^2\right)$, since the derivative of $\tanh$ is $1 - \tanh^2$.
• Gradients: $\frac{\partial J}{\partial W^{(l)}} = \delta^{(l)}\,(a^{(l-1)})^\top$ and $\frac{\partial J}{\partial b^{(l)}} = \delta^{(l)}$, averaged over the examples.

We train it with gradient descent from several random starting weights:

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

def forward_xor(X, params):
    W1, b1, W2, b2 = params
    A1 = np.tanh(X @ W1.T + b1)                 # hidden layer (4 examples x 2 neurons)
    P = sigmoid(A1 @ W2.T + b2).ravel()         # output probability
    return A1, P

def train_xor_network(seed, hidden=2, learning_rate=0.5, epochs=5000):
    rng = np.random.default_rng(seed)
    params = [rng.normal(0, 1, (hidden, 2)), np.zeros(hidden), rng.normal(0, 1, (1, hidden)), np.zeros(1)]
    y = targets['XOR'].astype(float)
    losses = []
    for _ in range(epochs):
        W1, b1, W2, b2 = params
        A1, P = forward_xor(X, params)
        losses.append(-np.mean(y * np.log(P) + (1 - y) * np.log(1 - P)))
        d2 = (P - y)[:, None] / len(y)              # output error
        d1 = (d2 @ W2) * (1 - A1 ** 2)             # hidden error, back through W2 and tanh
        grads = [d1.T @ X, d1.sum(axis=0), d2.T @ A1, d2.sum(axis=0)]
        params = [p - learning_rate * g for p, g in zip(params, grads)]
    return params, losses

runs = {}
for seed in range(6):
    params, losses = train_xor_network(seed)
    _, P = forward_xor(X, params)
    runs[seed] = (params, losses)
    solved = np.array_equal((P >= 0.5).astype(int), targets['XOR'])
    print(f"Start {seed}: final cost {losses[-1]:.4f}, probabilities {P.round(3)} -> {'solved' if solved else 'stuck'}")

Start 0: final cost 0.3471, probabilities [0.    0.999 0.5   0.5  ] -> stuck
Start 1: final cost 0.3471, probabilities [0.001 0.5   0.999 0.5  ] -> stuck
Start 2: final cost 0.0013, probabilities [0.001 0.998 0.998 0.001] -> solved

Start 3: final cost 0.0013, probabilities [0.001 0.998 0.998 0.001] -> solved
Start 4: final cost 0.0013, probabilities [0.002 0.999 0.999 0.002] -> solved
Start 5: final cost 0.0013, probabilities [0.001 0.998 0.998 0.001] -> solved

Four of the six starts solve XOR, with probabilities very close to 0 and 1. The other two get stuck with a cost of 0.347, answering 0.5 for two of the examples: like K-Means in chapter 5, gradient descent can end up in a local minimum. It is one of the reasons why training neural networks requires care.

Let's look inside one of the networks that solved it. In the middle, the decision boundary in the original space is no longer a single straight line: the network classifies as "off" a diagonal band that contains (0, 0) and (1, 1), and as "on" the two regions on either side of it. On the right, the same four examples, but described by the two hidden neurons: the hidden layer has moved them to a new space where a single straight line does separate them.

params, losses = runs[2]
W1, b1, W2, b2 = params
A1, P = forward_xor(X, params)
y = targets['XOR']

g1, g2 = np.meshgrid(np.linspace(-0.5, 1.5, 200), np.linspace(-0.5, 1.5, 200))
_, grid_p = forward_xor(np.c_[g1.ravel(), g2.ravel()], params)

fig, axes = plt.subplots(1, 3, figsize=(17, 4.5))
axes[0].plot(losses)
axes[0].set_xlabel('Epoch')
axes[0].set_ylabel('Cost (log loss)')
axes[0].set_title('Training the XOR network')
axes[1].contourf(g1, g2, grid_p.reshape(g1.shape) >= 0.5, alpha=0.2, cmap='coolwarm_r')
axes[1].set_title('Original space: a non-linear boundary')
axes[2].set_title('Hidden layer: now a line separates them')
for ax, pts in [(axes[1], X), (axes[2], A1)]:
    ax.scatter(*pts[y == 0].T, marker='^', s=150, color='tab:red', label='0 (off)')
    ax.scatter(*pts[y == 1].T, marker='o', s=150, color='tab:blue', label='1 (on)')
    ax.legend(fontsize=8)
axes[1].set_xlabel('$x_1$')
axes[1].set_ylabel('$x_2$')
axes[2].set_xlabel('hidden neuron 1')
axes[2].set_ylabel('hidden neuron 2')
xs = np.linspace(-1.1, 1.1, 2)
axes[2].plot(xs, -(W2[0, 0] * xs + b2[0]) / W2[0, 1], 'k--')   # where the output neuron gives 0.5
axes[2].set_xlim(-1.1, 1.1)
axes[2].set_ylim(-1.1, 1.1)
plt.show()
print("Hidden-layer values of the four examples:")
for xi, ai in zip(X, A1):
    print(f"  {xi.astype(int)} -> {ai.round(3)}")

Hidden-layer values of the four examples:
  [0 0] -> [-0.947 -0.958]
  [0 1] -> [ 0.96 -1.  ]
  [1 0] -> [-1.    0.97]
  [1 1] -> [-0.958 -0.969]

In the hidden space the two "off" examples, (0, 0) and (1, 1), end up together near (−0.95, −0.96), while the two "on" examples move to opposite corners: (0, 1) goes to (0.96, −1) and (1, 0) to (−1, 0.97). One straight line now separates them. Each hidden neuron has learned a feature by itself, one of them "only switch 2 is on" and the other "only switch 1 is on", and the output neuron only has to combine them. That is what hidden layers do: build intermediate features that make the final decision simple.

Step 5: Activation functions

---

The nonlinearity is essential: without it, several layers would collapse into a single linear transformation. The book presents the most common activation functions. Their derivatives matter as much as the functions themselves, because backpropagation multiplies them layer after layer:

z = np.linspace(-5, 5, 400)
activations = {
    'step': (np.where(z >= 0, 1.0, 0.0), np.zeros_like(z)),
    'sigmoid': (sigmoid(z), sigmoid(z) * (1 - sigmoid(z))),
    'tanh': (np.tanh(z), 1 - np.tanh(z) ** 2),
    'ReLU': (np.maximum(0, z), (z > 0).astype(float)),
}
fig, axes = plt.subplots(1, 2, figsize=(13, 4))
for name, (f, df) in activations.items():
    axes[0].plot(z, f, label=name)
    axes[1].plot(z, df, label=name)
axes[0].set_title('Activation functions')
axes[1].set_title('Their derivatives')
for ax in axes:
    ax.set_xlabel('z')
    ax.grid(True)
    ax.legend()
plt.show()

print(f"Maximum slope of the sigmoid: {np.max(activations['sigmoid'][1]):.3f}")
print(f"Ten sigmoid layers can shrink a gradient by up to 0.25 ** 10 = {0.25 ** 10:.1e}")

Maximum slope of the sigmoid: 0.250
Ten sigmoid layers can shrink a gradient by up to 0.25 ** 10 = 9.5e-07

The step has zero derivative everywhere (and none at 0), so gradient descent cannot learn with it: that is why modern neurons use smooth functions. The sigmoid has a slope of at most 0.25, so a gradient that crosses many sigmoid layers can become tiny, the vanishing gradient problem the book describes. ReLU has slope 1 for every positive input, which is one of the reasons it made deep networks trainable.

Step 6: Martha's complete example

---

In her notes, Martha works through the small diagnostic network the team saw at AIA, with every number: two symptoms ($x_1 = 0.5$, $x_2 = 0.7$, the second one is cough), a hidden layer of 3 ReLU neurons and an output layer of 3 neurons with softmax, one for each diagnosis: cold, flu and pneumonia. The patient has flu. Softmax turns the output scores (logits) into probabilities, and the cost is the cross-entropy, which for a one-hot label is minus the logarithm of the probability of the correct class:

$$\hat{y}_i = \frac{e^{z_i}}{\sum_k e^{z_k}} \qquad\qquad \mathcal{L} = -\sum_i y_i \log \hat{y}_i$$

x = np.array([0.5, 0.7])
W1 = np.array([[0.8, 0.5], [0.3, 0.9], [1.2, 0.4]]);                b1 = np.array([0.1, -0.2, 0.0])
W2 = np.array([[1.1, 0.3, 0.5], [0.2, 1.0, 0.7], [0.6, 0.4, 1.3]]); b2 = np.array([0.0, 0.1, -0.1])
y_true = np.array([0, 1, 0])                                        # one-hot: flu
classes = ['cold', 'flu', 'pneumonia']

def softmax(z):
    e = np.exp(z - z.max())
    return e / e.sum()

def forward(x, W1, b1, W2, b2):
    z1 = W1 @ x + b1
    a1 = np.maximum(0, z1)          # ReLU
    z2 = W2 @ a1 + b2
    return z1, a1, z2, softmax(z2)

z1, a1, z2, y_hat = forward(x, W1, b1, W2, b2)
loss = -np.log(y_hat[1])
print("FORWARD PASS")
print(f"  hidden layer:  z1 = {z1},  a1 = {a1}")
print(f"  output logits: z2 = {z2}")
print(f"  probabilities: {dict((c, round(float(v), 4)) for c, v in zip(classes, y_hat))}")
print(f"  cross-entropy loss: {loss:.4f}")

FORWARD PASS
  hidden layer:  z1 = [0.85 0.58 0.88],  a1 = [0.85 0.58 0.88]
  output logits: z2 = [1.549 1.466 1.786]
  probabilities: {'cold': 0.3137, 'flu': 0.2887, 'pneumonia': 0.3976}
  cross-entropy loss: 1.2423

The probabilities are 0.314 for cold, 0.289 for flu and 0.398 for pneumonia: the network does not recognize the flu yet, and the loss is 1.242. In the book the same example shows 0.288 and a loss of 1.244, because the intermediate values are rounded to three decimals; the exact probability of flu is 0.2887.

Now the backward pass. With softmax and cross-entropy, the error signal of the output layer is simply $\delta^{(2)} = \hat{y} - y$. Each weight matrix gets the outer product of the error of its destination layer and the activations of its origin layer, and the error travels back to the hidden layer through $(W^{(2)})^\top$ and the derivative of ReLU:

$$\frac{\partial \mathcal{L}}{\partial W^{(2)}} = \delta^{(2)} (a^{(1)})^\top \qquad \delta^{(1)} = \left((W^{(2)})^\top \delta^{(2)}\right) \odot \text{ReLU}'(z^{(1)}) \qquad \frac{\partial \mathcal{L}}{\partial W^{(1)}} = \delta^{(1)} x^\top$$

def backward(x, y_true, z1, a1, y_hat, W2):
    d2 = y_hat - y_true                          # output error
    d1 = (W2.T @ d2) * (z1 > 0)                  # hidden error (ReLU derivative: 1 if z > 0)
    return {'W2': np.outer(d2, a1), 'b2': d2, 'W1': np.outer(d1, x), 'b1': d1}

grads = backward(x, y_true, z1, a1, y_hat, W2)
print("BACKWARD PASS")
for name, g in grads.items():
    print(f"  dL/d{name} =\n{g.round(3)}")
print(f"\n  The weight from h1 to 'flu', element (2, 1) of dL/dW2: {grads['W2'][1, 0]:.3f}")
print(f"  The weight from cough (x2) to h2, element (2, 2) of dL/dW1: {grads['W1'][1, 1]:.3f}")

BACKWARD PASS
  dL/dW2 =
[[ 0.267  0.182  0.276]
 [-0.605 -0.413 -0.626]
 [ 0.338  0.231  0.35 ]]
  dL/db2 =
[ 0.314 -0.711  0.398]
  dL/dW1 =
[[ 0.221  0.309]
 [-0.229 -0.321]
 [ 0.088  0.123]]
  dL/db1 =
[ 0.441 -0.458  0.176]

  The weight from h1 to 'flu', element (2, 1) of dL/dW2: -0.605
  The weight from cough (x2) to h2, element (2, 2) of dL/dW1: -0.321

These are the gradients of Martha's notes. The error for flu is negative (−0.711): the network under-predicts flu, so the weights that feed it must grow. For example, the weight from $h_1$ to flu has a gradient of −0.605, and the weight from cough to $h_2$ one of −0.321, the two elements Martha checks by hand. The book shows −0.712, and a few third decimals differ by one unit, for the same rounding reason. All three hidden pre-activations are positive, so ReLU lets the whole error through.

Finally, one update with learning rate $\eta = 0.1$, $\theta \leftarrow \theta - \eta \frac{\partial \mathcal{L}}{\partial \theta}$, and a new forward pass to check that the network has moved in the right direction:

eta = 0.1
W1_new, b1_new = W1 - eta * grads['W1'], b1 - eta * grads['b1']
W2_new, b2_new = W2 - eta * grads['W2'], b2 - eta * grads['b2']
with np.printoptions(precision=4):
    print(f"W2 after the update:\n{W2_new}\nb2 after the update: {b2_new}")
    print(f"W1 after the update:\n{W1_new}\nb1 after the update: {b1_new}\n")

_, _, _, y_hat_new = forward(x, W1_new, b1_new, W2_new, b2_new)
print(f"Probabilities before: {dict((c, round(float(v), 3)) for c, v in zip(classes, y_hat))}")
print(f"Probabilities after:  {dict((c, round(float(v), 3)) for c, v in zip(classes, y_hat_new))}")
print(f"Loss: {loss:.3f} -> {-np.log(y_hat_new[1]):.3f}")

W2 after the update:
[[1.0733 0.2818 0.4724]
 [0.2605 1.0413 0.7626]
 [0.5662 0.3769 1.265 ]]
b2 after the update: [-0.0314  0.1711 -0.1398]
W1 after the update:
[[0.7779 0.4691]
 [0.3229 0.9321]
 [1.1912 0.3877]]
b1 after the update: [ 0.0559 -0.1542 -0.0176]

Probabilities before: {'cold': 0.314, 'flu': 0.289, 'pneumonia': 0.398}
Probabilities after:  {'cold': 0.274, 'flu': 0.378, 'pneumonia': 0.347}
Loss: 1.242 -> 0.972

After a single update the probability of flu rises from 0.289 to 0.378 and the loss falls from 1.242 to 0.972, exactly the direction we wanted (the book gets 0.379, working with rounded values). Training a network is nothing more than repeating this cycle, forward pass, backward pass and update, thousands of times over many examples.

Step 7: Checking backpropagation

---

How do we know our derivatives are right? A classic test is to compare them with numerical derivatives: nudge each weight by a tiny amount $\epsilon$ in both directions and measure how much the loss changes,

$$\frac{\partial \mathcal{L}}{\partial \theta} \approx \frac{\mathcal{L}(\theta + \epsilon) - \mathcal{L}(\theta - \epsilon)}{2\epsilon}$$

It is far too slow for training (two forward passes per weight), but it is a perfect check:

def loss_of(params):
    return -np.log(forward(x, *params)[3][1])

params = [W1, b1, W2, b2]
epsilon = 1e-6
largest_difference = 0.0
for k, name in enumerate(['W1', 'b1', 'W2', 'b2']):
    numeric = np.zeros_like(params[k])
    for idx in np.ndindex(params[k].shape):
        plus = [p.copy() for p in params]
        minus = [p.copy() for p in params]
        plus[k][idx] += epsilon
        minus[k][idx] -= epsilon
        numeric[idx] = (loss_of(plus) - loss_of(minus)) / (2 * epsilon)
    largest_difference = max(largest_difference, np.max(np.abs(numeric - grads[name])))
print(f"Largest difference between backpropagation and numerical derivatives: {largest_difference:.1e}")

Largest difference between backpropagation and numerical derivatives: 1.5e-10

The difference is of the order of $10^{-10}$, just rounding error: our backpropagation computes the correct derivatives. Deep learning libraries do all this bookkeeping automatically, as we will see in the next notebook.

Conclusion

---

• A perceptron is a weighted sum followed by a step. Its learning rule finds a separating line whenever one exists, as for AND and OR, but no single line can solve XOR.
• A hidden layer with a nonlinear activation transforms the inputs into a space where the classes become separable. That is the key idea of multilayer networks.
• Backpropagation is the chain rule applied from the output back to the input: an error signal $\delta$ per layer, and gradients that are outer products of errors and activations. We reproduced Martha's complete example and checked it against numerical derivatives.
• Training can get stuck, and the choice of activation matters: the sigmoid can make gradients vanish, which ReLU avoids.
• Writing all of this by hand is Emma's "level 0": ideal for understanding, and far too laborious for real networks. The next notebook climbs to the higher levels with TensorFlow and Keras, and builds the convolutional network for the dermatology images.