A Two-Layer Network from Scratch (Deep Learning Tutorial)

By the end of this lesson you will have written a complete neural network in about 40 lines of NumPy. It will learn a curved pattern that a straight line cannot, draw its own decision regions, and then read handwritten digits with high accuracy.

You already have every piece: layers from Layers & the Forward Pass, cross-entropy from Loss Functions, and the six backward lines from Backpropagation Step by Step. Now we join them.

The plan

Our network has two layers of weights:

Input of 2 numbers (a point's x and y position).
Hidden layer of 64 ReLU neurons.
Output layer of 3 neurons with softmax, one per class.

Training repeats four steps many times: forward pass, loss, backward pass, update. Think of a cricket bowler at practice. Bowl (forward), see where it lands (loss), work out what went wrong with the wrist and run-up (backward), adjust (update), and bowl again.

Step 1: a dataset with a twist

We need data that a straight line cannot separate. Imagine a rangoli drawn with three colours of powder, each colour swirling out from the centre as a spiral arm. Each point is a pinch of powder, and its class is its colour.

PythonRun it yourself

import numpy as np
import matplotlib.pyplot as plt

def make_spiral(n_per_class=100, classes=3, seed=0):
    rng = np.random.default_rng(seed)
    X, y = [], []
    for c in range(classes):
        r = np.linspace(0.05, 1, n_per_class)              # distance from centre
        angle = c * 2 * np.pi / classes + 4 * r + rng.normal(0, 0.3, n_per_class)
        X.append(np.column_stack([r * np.cos(angle), r * np.sin(angle)]))
        y.append(np.full(n_per_class, c))
    return np.vstack(X), np.concatenate(y)

X, y = make_spiral()
print("X shape:", X.shape, " y shape:", y.shape, " classes:", np.unique(y))

# a straight-line model first: softmax regression with no hidden layer
rng = np.random.default_rng(1)
W, b = rng.normal(0, 0.1, (2, 3)), np.zeros(3)
for _ in range(1000):
    Z = X @ W + b
    P = np.exp(Z - Z.max(axis=1, keepdims=True))
    P /= P.sum(axis=1, keepdims=True)
    dZ = (P - np.eye(3)[y]) / len(y)
    W -= 1.0 * X.T @ dZ
    b -= 1.0 * dZ.sum(axis=0)
print("straight-line model accuracy:", round(np.mean(P.argmax(axis=1) == y), 3))

plt.scatter(X[:, 0], X[:, 1], c=y, cmap="brg", s=15)
plt.title("Three-arm spiral")
plt.show()

Output

X shape: (300, 2)  y shape: (300,)  classes: [0 1 2]
straight-line model accuracy: 0.563

The straight-line model gets only about 56% right. Guessing would give 33%, so it learned something, but its boundaries are straight cuts through a swirl and each region mixes colours. We need a hidden layer.

Step 2: the network class

Here is the full network. Read the comments: each method is one stage of the training cycle. The starting weights are scaled by sqrt(2 / n_in); Weight Initialisation explains why.

PythonRun it yourself

import numpy as np

class TwoLayerNet:
    def __init__(self, n_in, n_hidden, n_out, seed=42):
        rng = np.random.default_rng(seed)
        self.W1 = rng.normal(0, 1, (n_in, n_hidden)) * np.sqrt(2 / n_in)
        self.b1 = np.zeros(n_hidden)
        self.W2 = rng.normal(0, 1, (n_hidden, n_out)) * np.sqrt(2 / n_hidden)
        self.b2 = np.zeros(n_out)

    def forward(self, X):
        self.X = X
        self.Z1 = X @ self.W1 + self.b1          # save values for the backward pass
        self.A1 = np.maximum(0, self.Z1)
        Z2 = self.A1 @ self.W2 + self.b2
        e = np.exp(Z2 - Z2.max(axis=1, keepdims=True))
        self.P = e / e.sum(axis=1, keepdims=True)
        return self.P

    def loss(self, y):
        return -np.log(self.P[np.arange(len(y)), y] + 1e-12).mean()

    def backward(self, y):
        n = len(y)
        dZ2 = self.P.copy()
        dZ2[np.arange(n), y] -= 1                  # probs - onehot
        dZ2 /= n
        self.dW2 = self.A1.T @ dZ2
        self.db2 = dZ2.sum(axis=0)
        dZ1 = (dZ2 @ self.W2.T) * (self.Z1 > 0)    # back through ReLU
        self.dW1 = self.X.T @ dZ1
        self.db1 = dZ1.sum(axis=0)

    def step(self, lr):
        self.W1 -= lr * self.dW1
        self.b1 -= lr * self.db1
        self.W2 -= lr * self.dW2
        self.b2 -= lr * self.db2

    def predict(self, X):
        return self.forward(X).argmax(axis=1)

net = TwoLayerNet(2, 64, 3)
X = np.array([[0.1, 0.2], [-0.5, 0.4]])
print("probabilities before training:\n", net.forward(X).round(3))
print("parameters:", sum(a.size for a in (net.W1, net.b1, net.W2, net.b2)))

Output

probabilities before training:
 [[0.287 0.33  0.383]
 [0.289 0.262 0.449]]
parameters: 387

Before training, the probabilities are close to one third each: the network has no opinion yet. It has 2×64 + 64 + 64×3 + 3 = 387 parameters to learn.

Notice one small trick in backward: instead of building a one-hot matrix, we copy the probabilities and subtract 1 at the true class. It gives exactly probs - onehot, using less memory.

Step 3: train and draw the decision regions

Now the training loop. Each pass over all 300 points is one epoch. We print the loss and accuracy now and then, and finally colour the whole plane by the network's prediction.

PythonRun it yourself

import numpy as np
import matplotlib.pyplot as plt

def make_spiral(n_per_class=100, classes=3, seed=0):
    rng = np.random.default_rng(seed)
    X, y = [], []
    for c in range(classes):
        r = np.linspace(0.05, 1, n_per_class)
        angle = c * 2 * np.pi / classes + 4 * r + rng.normal(0, 0.3, n_per_class)
        X.append(np.column_stack([r * np.cos(angle), r * np.sin(angle)]))
        y.append(np.full(n_per_class, c))
    return np.vstack(X), np.concatenate(y)

class TwoLayerNet:
    def __init__(self, n_in, n_hidden, n_out, seed=42):
        rng = np.random.default_rng(seed)
        self.W1 = rng.normal(0, 1, (n_in, n_hidden)) * np.sqrt(2 / n_in)
        self.b1 = np.zeros(n_hidden)
        self.W2 = rng.normal(0, 1, (n_hidden, n_out)) * np.sqrt(2 / n_hidden)
        self.b2 = np.zeros(n_out)

    def forward(self, X):
        self.X, self.Z1 = X, X @ self.W1 + self.b1
        self.A1 = np.maximum(0, self.Z1)
        Z2 = self.A1 @ self.W2 + self.b2
        e = np.exp(Z2 - Z2.max(axis=1, keepdims=True))
        self.P = e / e.sum(axis=1, keepdims=True)
        return self.P

    def loss(self, y):
        return -np.log(self.P[np.arange(len(y)), y] + 1e-12).mean()

    def backward(self, y):
        dZ2 = self.P.copy()
        dZ2[np.arange(len(y)), y] -= 1
        dZ2 /= len(y)
        self.dW2, self.db2 = self.A1.T @ dZ2, dZ2.sum(axis=0)
        dZ1 = (dZ2 @ self.W2.T) * (self.Z1 > 0)
        self.dW1, self.db1 = self.X.T @ dZ1, dZ1.sum(axis=0)

    def step(self, lr):
        for name in ("W1", "b1", "W2", "b2"):
            setattr(self, name, getattr(self, name) - lr * getattr(self, "d" + name))

    def predict(self, X):
        return self.forward(X).argmax(axis=1)

X, y = make_spiral()
net = TwoLayerNet(2, 64, 3)
for epoch in range(1, 1501):
    net.forward(X)                 # 1. forward
    loss = net.loss(y)             # 2. loss
    net.backward(y)                # 3. backward
    net.step(lr=1.0)               # 4. update
    if epoch in (1, 50, 200, 500, 1500):
        acc = np.mean(net.P.argmax(axis=1) == y)
        print(f"epoch {epoch:>4}: loss {loss:.3f}, accuracy {acc:.3f}")

xx, yy = np.meshgrid(np.linspace(-1.2, 1.2, 200), np.linspace(-1.2, 1.2, 200))
grid = np.column_stack([xx.ravel(), yy.ravel()])
plt.contourf(xx, yy, net.predict(grid).reshape(xx.shape), alpha=0.3, cmap="brg")
plt.scatter(X[:, 0], X[:, 1], c=y, cmap="brg", s=15)
plt.title("Decision regions learned by the network")
plt.show()

Output

epoch    1: loss 1.088, accuracy 0.293
epoch   50: loss 0.349, accuracy 0.913
epoch  200: loss 0.122, accuracy 0.973
epoch  500: loss 0.052, accuracy 0.997
epoch 1500: loss 0.014, accuracy 1.000

The loss falls and the accuracy climbs far above the straight-line model. The picture shows curved, swirling regions that follow the arms. The hidden layer has learned 64 new features, and the output layer combines them, just like the hand-made XOR network in the introduction, but learned automatically.

Step 4: the same network on handwritten digits

Now a real task. scikit-learn's digits dataset has 1,797 handwritten digits, each an 8 × 8 image, so 64 pixel inputs. We split it into training and test sets, scale the pixels to 0 to 1, and train the very same class with 64 inputs and 10 outputs.

PythonRun it yourself

import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

class TwoLayerNet:
    def __init__(self, n_in, n_hidden, n_out, seed=42):
        rng = np.random.default_rng(seed)
        self.W1 = rng.normal(0, 1, (n_in, n_hidden)) * np.sqrt(2 / n_in)
        self.b1 = np.zeros(n_hidden)
        self.W2 = rng.normal(0, 1, (n_hidden, n_out)) * np.sqrt(2 / n_hidden)
        self.b2 = np.zeros(n_out)

    def forward(self, X):
        self.X, self.Z1 = X, X @ self.W1 + self.b1
        self.A1 = np.maximum(0, self.Z1)
        Z2 = self.A1 @ self.W2 + self.b2
        e = np.exp(Z2 - Z2.max(axis=1, keepdims=True))
        self.P = e / e.sum(axis=1, keepdims=True)
        return self.P

    def backward(self, y):
        dZ2 = self.P.copy()
        dZ2[np.arange(len(y)), y] -= 1
        dZ2 /= len(y)
        self.dW2, self.db2 = self.A1.T @ dZ2, dZ2.sum(axis=0)
        dZ1 = (dZ2 @ self.W2.T) * (self.Z1 > 0)
        self.dW1, self.db1 = self.X.T @ dZ1, dZ1.sum(axis=0)

    def step(self, lr):
        for name in ("W1", "b1", "W2", "b2"):
            setattr(self, name, getattr(self, name) - lr * getattr(self, "d" + name))

X, y = load_digits(return_X_y=True)
X = X / 16.0                                     # pixels were 0 to 16
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=42, stratify=y)

net = TwoLayerNet(64, 64, 10)
for epoch in range(1, 501):
    net.forward(X_train)
    net.backward(y_train)
    net.step(lr=0.5)
    if epoch in (1, 50, 200, 500):
        train_acc = np.mean(net.P.argmax(axis=1) == y_train)
        test_acc = np.mean(net.forward(X_test).argmax(axis=1) == y_test)
        print(f"epoch {epoch:>3}: train accuracy {train_acc:.3f}, test accuracy {test_acc:.3f}")

Output

epoch   1: train accuracy 0.118, test accuracy 0.253
epoch  50: train accuracy 0.944, test accuracy 0.929
epoch 200: train accuracy 0.984, test accuracy 0.964
epoch 500: train accuracy 0.999, test accuracy 0.971

Your own network, written from scratch, reads digits it has never seen with high accuracy. The test accuracy is a little below the training accuracy, which is normal. Every deep learning library does the same four steps at its heart, only faster and with many more layer types.

Common mistakes

Not scaling inputs. Raw pixels from 0 to 16 (or 0 to 255 in photos) make the first updates huge. Divide to get 0 to 1.
Measuring accuracy only on training data. Always keep a test set the network never trains on.
Calling backward before forward. The backward pass uses values saved during the forward pass, so the order matters.
Starting all weights at zero. Every hidden neuron would compute the same thing and learn the same thing. Random starting weights break this symmetry.

Key points

A training loop repeats forward pass, loss, backward pass and update.
A two-layer network with ReLU learns curved boundaries that a straight-line model cannot.
Save forward values (X, Z1, A1, P) because the backward pass needs them.
The same small class can solve a real task such as handwritten digits.
Always judge a network by its accuracy on test data it has never seen.

Exercise

Train the spiral network from Step 3 with hidden layers of 4, 16 and 64 neurons (1,500 epochs, learning rate 1.0 each time). Print the final training accuracy for each size. How does the hidden layer size change what the network can learn?

Show solution

Python

import numpy as np

def make_spiral(n_per_class=100, classes=3, seed=0):
    rng = np.random.default_rng(seed)
    X, y = [], []
    for c in range(classes):
        r = np.linspace(0.05, 1, n_per_class)
        angle = c * 2 * np.pi / classes + 4 * r + rng.normal(0, 0.3, n_per_class)
        X.append(np.column_stack([r * np.cos(angle), r * np.sin(angle)]))
        y.append(np.full(n_per_class, c))
    return np.vstack(X), np.concatenate(y)

def train(n_hidden, epochs=1500, lr=1.0):
    X, y = make_spiral()
    rng = np.random.default_rng(42)
    W1 = rng.normal(0, 1, (2, n_hidden)) * np.sqrt(2 / 2)
    b1 = np.zeros(n_hidden)
    W2 = rng.normal(0, 1, (n_hidden, 3)) * np.sqrt(2 / n_hidden)
    b2 = np.zeros(3)
    for _ in range(epochs):
        Z1 = X @ W1 + b1
        A1 = np.maximum(0, Z1)
        Z2 = A1 @ W2 + b2
        P = np.exp(Z2 - Z2.max(axis=1, keepdims=True))
        P /= P.sum(axis=1, keepdims=True)
        dZ2 = P.copy()
        dZ2[np.arange(len(y)), y] -= 1
        dZ2 /= len(y)
        dZ1 = (dZ2 @ W2.T) * (Z1 > 0)
        W2 -= lr * A1.T @ dZ2
        b2 -= lr * dZ2.sum(axis=0)
        W1 -= lr * X.T @ dZ1
        b1 -= lr * dZ1.sum(axis=0)
    return np.mean(P.argmax(axis=1) == y)

for size in [4, 16, 64]:
    print(f"{size:>2} hidden neurons: training accuracy {train(size):.3f}")

Output

 4 hidden neurons: training accuracy 0.960
16 hidden neurons: training accuracy 1.000
64 hidden neurons: training accuracy 1.000

Even 4 hidden neurons reach 96%, but they cannot bend the boundaries enough to follow every point of the three arms. With 16 or 64 neurons there are more features to combine, and the network fits every training point. Bigger is not always better on new data, though, which is the topic of Overfitting & Regularisation.

Quick Quiz

1. What are the four steps repeated in every training iteration?

2. Why does the network save Z1 and A1 during the forward pass?

dW2 uses A1, and the ReLU step uses Z1 > 0.

3. Why should weights not all start at zero?

4. Why do we divide the digit pixels by 16?

The plan

Step 1: a dataset with a twist

Step 2: the network class

Step 3: train and draw the decision regions

Step 4: the same network on handwritten digits

Common mistakes

Key points

Exercise

Quick Quiz

Learn Deep Learning with a mentor