Neural Network From Scratch

Frameworks turn backpropagation into one line, which is also why it stays mysterious. Writing the forward and backward pass of a small network in plain NumPy shows how little is there: four matmuls, a mask, and a sum.

The idea

A hand-written network is a set of cached forward values plus four gradient rules.

• dlogits = (p − y) / B — sigmoid followed by binary cross-entropy collapses into this.
• dW₂ = aᵀ dlogits, db₂ = Σ dlogits — each weight's gradient is its input times the upstream gradient, summed over the batch.
• dh = (dlogits W₂ᵀ) · (h > 0) — back through the output layer, then times the ReLU's local derivative.
• dW₁ = Xᵀ dh, db₁ = Σ dh — the same rule, one layer deeper.

Each layer caches what backward needs: the linear layer its input, the ReLU its pre-activation, since the derivative is 1 where h > 0. Because the loss is a mean over the batch, every gradient carries a 1/B.

Worked example

Two noisy rings: 64 points at radius 1 and 64 at radius 2 with noise σ = 0.25, so the classes nearly touch. The model is 2 → 8 (ReLU) → 1 (sigmoid), trained full-batch with lr = 0.5.

• Loss at step 0: 1.0552
• Loss at step 200: 0.0805
• Final loss 0.0420, accuracy 0.9844

Accuracy stops short of 1.0 because the noisy rings overlap. The hand-derived gradient agrees with a central difference of the loss to eight digits: 0.0005647442 analytic versus 0.0005647442 numeric.

In code

import numpy as np

rng = np.random.default_rng(0)

def ring(radius, n=64):
    t = rng.uniform(0, 2 * np.pi, n)
    return np.c_[np.cos(t), np.sin(t)] * (radius + rng.normal(0, 0.25, (n, 1)))

X = np.vstack([ring(1.0), ring(2.0)])
y = np.concatenate([np.ones((64, 1)), np.zeros((64, 1))])
W1 = rng.normal(0, 1.0, (2, 8)); b1 = np.zeros(8)
W2 = rng.normal(0, 1.0, (8, 1)); b2 = np.zeros(1)

def forward(X):
    h = X @ W1 + b1                     # pre-activation, kept for the ReLU mask
    p = 1 / (1 + np.exp(-(np.maximum(h, 0) @ W2 + b2)))
    return h, p

def bce(p):
    return -np.mean(y * np.log(p + 1e-12) + (1 - y) * np.log(1 - p + 1e-12))

for step in range(2000):
    h, p = forward(X)
    loss = bce(p)
    dlogits = (p - y) / len(X)
    dW2 = np.maximum(h, 0).T @ dlogits; db2 = dlogits.sum(0)
    dh = (dlogits @ W2.T) * (h > 0)     # through the ReLU
    dW1 = X.T @ dh; db1 = dh.sum(0)
    for param, grad in ((W1, dW1), (b1, db1), (W2, dW2), (b2, db2)):
        param -= 0.5 * grad
    if step in (0, 200):
        print(step, round(float(loss), 4))          # 1.0552  0.0805

h, p = forward(X)
dW1 = X.T @ (((p - y) / len(X) @ W2.T) * (h > 0))
print(round(float(loss), 4), round(float(((p > 0.5) == (y > 0.5)).mean()), 4))   # 0.042 0.9844

e = 1e-6; W1[0, 0] += e; plus = bce(forward(X)[1])
W1[0, 0] -= 2 * e; minus = bce(forward(X)[1]); W1[0, 0] += e
print(dW1[0, 0], (plus - minus) / (2 * e))   # 0.0005647442 0.0005647442```

## Check yourself

1. Which line is the chain rule, and what local derivative does it multiply in?
2. Why is every gradient divided by `len(X)`?
3. Why must the ReLU mask come from the forward pass?

## Key takeaways

- Backprop is a short list of rules applied in reverse order, each reusing a cached forward value.
- Sigmoid plus BCE produces `dlogits = (p − y)/B`, which is why that pairing is everywhere.
- Check the derivation with a central difference; it costs three lines.