How Neural Networks Work: Forward Pass & Backprop in NumPy
This post builds a 2-layer network in plain NumPy and trains it on XOR: forward pass, binary cross-entropy loss, backpropagation by hand, and a numerical check that the gradients are right. It is the same cycle that runs when you call model.fit() in Keras or write a training loop in PyTorch, with every step written out.
pip install numpyThe Problem: Non-Linear Data
Consider the XOR problem:
import numpy as np
X = np.array([[0,0],[0,1],[1,0],[1,1]]) # inputs
y = np.array([[0],[1],[1],[0]]) # XOR outputs
No straight line separates the 1s from the 0s. Logistic regression on the raw inputs fails. You need a model that can bend the decision boundary, and that is what a hidden layer does.
What a Neuron Actually Computes
A single neuron takes a vector of inputs x, multiplies each by a weight, adds a bias, and passes the result through a non-linear activation function:
output = activation(W · x + b)That dot product W · x is just a weighted sum. The activation function is what makes it non-linear. Without it, stacking layers would still give you a linear model.
We'll use the sigmoid activation: σ(z) = 1 / (1 + e^(-z))
def sigmoid(z):
return 1 / (1 + np.exp(-z))
def sigmoid_deriv(z):
s = sigmoid(z)
return s * (1 - s)
Network Architecture
We'll build:
- Input layer: 2 neurons (x1, x2)
- Hidden layer: 4 neurons with sigmoid activation
- Output layer: 1 neuron with sigmoid activation (binary output)
np.random.seed(42)
# Weights and biases initialized randomly
W1 = np.random.randn(2, 4) * 0.1 # (input_size, hidden_size)
b1 = np.zeros((1, 4))
W2 = np.random.randn(4, 1) * 0.1 # (hidden_size, output_size)
b2 = np.zeros((1, 1))
Forward Pass
Data flows left to right through the network. Each layer applies the linear transformation then the activation:
def forward(X, W1, b1, W2, b2):
# Hidden layer
Z1 = X @ W1 + b1 # (4, 2) @ (2, 4) = (4, 4)
A1 = sigmoid(Z1) # apply activation elementwise
# Output layer
Z2 = A1 @ W2 + b2 # (4, 4) @ (4, 1) = (4, 1)
A2 = sigmoid(Z2) # final prediction (probability)
cache = (Z1, A1, Z2, A2)
return A2, cache
We store the intermediate values in cache because backprop needs them.
Loss Function
We use binary cross-entropy loss. It penalizes confident wrong predictions much harder than uncertain ones:
def compute_loss(y_true, y_pred):
m = y_true.shape[0]
# Clip to avoid log(0)
y_pred = np.clip(y_pred, 1e-9, 1 - 1e-9)
loss = -np.mean(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))
return loss
Backpropagation
Backprop is just the chain rule applied repeatedly from output to input. We compute how much each weight contributed to the loss and update it in the direction that reduces the loss.
def backward(X, y, W1, b1, W2, b2, cache, lr=0.1):
Z1, A1, Z2, A2 = cache
m = X.shape[0]
# Output layer gradients
dZ2 = A2 - y # per-example d(loss)/d(Z2); the mean adds the 1/m below
dW2 = A1.T @ dZ2 / m
db2 = np.sum(dZ2, axis=0, keepdims=True) / m
# Hidden layer gradients (chain rule through sigmoid)
dA1 = dZ2 @ W2.T
dZ1 = dA1 * sigmoid_deriv(Z1)
dW1 = X.T @ dZ1 / m
db1 = np.sum(dZ1, axis=0, keepdims=True) / m
# Gradient descent update
W2 -= lr * dW2
b2 -= lr * db2
W1 -= lr * dW1
b1 -= lr * db1
return W1, b1, W2, b2
Why is the output gradient just A2 - y? For one example, binary cross-entropy is L = -(y·log(a) + (1 - y)·log(1 - a)) with a = σ(z). Its derivative with respect to a is (a - y) / (a(1 - a)), and σ'(z) = a(1 - a). Multiply the two and the a(1 - a) terms cancel, leaving a - y. Because compute_loss averages over the m examples, every weight and bias gradient is divided by m.
For the hidden layer, the error flows back through W2 (dA1 = dZ2 @ W2.T) and is then scaled elementwise by the sigmoid derivative at Z1. Note that b1 and b2 must be passed in as arguments: if the function only reads them as globals, the line b2 -= lr * db2 raises UnboundLocalError, because augmented assignment makes b2 a local variable.
Training Loop
lr = 1.0
epochs = 10000
for epoch in range(epochs):
A2, cache = forward(X, W1, b1, W2, b2)
loss = compute_loss(y, A2)
W1, b1, W2, b2 = backward(X, y, W1, b1, W2, b2, cache, lr=lr)
if epoch % 1000 == 0:
print(f"Epoch {epoch:5d} | Loss: {loss:.4f}")
Output with np.random.seed(42):
Epoch 0 | Loss: 0.6934
Epoch 1000 | Loss: 0.6931
Epoch 2000 | Loss: 0.6931
Epoch 3000 | Loss: 0.6931
Epoch 4000 | Loss: 0.6931
Epoch 5000 | Loss: 0.6931
Epoch 6000 | Loss: 0.6931
Epoch 7000 | Loss: 0.6931
Epoch 8000 | Loss: 0.6928
Epoch 9000 | Loss: 0.0192
The long flat stretch is real, not a bug. 0.6931 is ln(2), the loss you get by predicting 0.5 for every input. With weights initialized at 0.1 scale, every hidden unit outputs close to 0.5 for all four inputs and the output weights are small, so the gradient that reaches W1 is tiny and the network sits near a saddle point for about 8,000 epochs. Once the hidden units start to differ, the loss drops quickly: in our run it went from 0.6928 at epoch 8000 to 0.0905 at epoch 8600 and 0.0192 at epoch 9000. The plateau length depends on the starting point. With np.random.seed(0) the drop came around epoch 3000 (loss 0.2264), and with * 1.0 instead of * 0.1 in the initialization the loss was already 0.0147 at epoch 1000.
Checking the Result
A2, _ = forward(X, W1, b1, W2, b2)
predictions = (A2 > 0.5).astype(int)
print("Input | Target | Predicted")
for xi, yi, pi in zip(X, y, predictions):
print(f" {xi} | {yi[0]} | {pi[0]}")
Output:
Input | Target | Predicted
[0 0] | 0 | 0
[0 1] | 1 | 1
[1 0] | 1 | 1
[1 1] | 0 | 0
The raw outputs after training were 0.0085, 0.995, 0.995 and 0.0044. The network learned XOR, a function no linear model can represent, by building a non-linear representation in the hidden layer.
Checking the Gradients Numerically
A hand-written backward pass is easy to get subtly wrong. The standard test is a finite-difference gradient check: nudge one weight up and down by a tiny eps, measure how the loss changes, and compare with the analytic gradient. Here it is for W1, the layer furthest from the loss:
def loss_at(W1, b1, W2, b2):
A2, _ = forward(X, W1, b1, W2, b2)
return compute_loss(y, A2)
# Random test point (not the trained weights)
rng = np.random.default_rng(0)
W1c, b1c = rng.normal(size=(2, 4)), rng.normal(size=(1, 4))
W2c, b2c = rng.normal(size=(4, 1)), rng.normal(size=(1, 1))
# Analytic gradient for W1, same formulas as backward()
_, (Z1, A1, Z2, A2) = forward(X, W1c, b1c, W2c, b2c)
m = X.shape[0]
dZ1 = ((A2 - y) @ W2c.T) * sigmoid_deriv(Z1)
dW1 = X.T @ dZ1 / m
# Numerical gradient with central differences
eps = 1e-6
num = np.zeros_like(W1c)
for i in np.ndindex(W1c.shape):
Wp, Wm = W1c.copy(), W1c.copy()
Wp[i] += eps
Wm[i] -= eps
num[i] = (loss_at(Wp, b1c, W2c, b2c) - loss_at(Wm, b1c, W2c, b2c)) / (2 * eps)
print(f"max abs difference: {np.abs(num - dW1).max():.1e}")
Output:
max abs difference: 2.1e-10A difference around 1e-10 means the formulas in backward() are correct. Running the same check on b1, W2 and b2 gave differences between 5e-11 and 2e-10. If you get something like 1e-2, there is a bug in a derivative or a transpose.
What This Means in Practice
When you call model.fit(X, y) in Keras or write a training loop in PyTorch, the same cycle runs under the hood, millions of times, with:
- Mini-batches instead of full-batch gradient descent
- Adam or RMSProp instead of vanilla SGD
- Dropout and batch norm for regularization
- Autograd computing the chain rule automatically
The plateau in the training log above is a small example of why this matters. The gradient that reaches W1 is multiplied by W2 and by σ'(Z1), and the sigmoid derivative is never larger than 0.25. Stack more sigmoid layers and that product shrinks at each one, which is the vanishing-gradient problem that pushed deep networks toward ReLU activations and careful initialization.
Frequently asked questions
Why does backpropagation use the chain rule?
Because the loss is a composition of functions: a linear transform, then an activation, layer after layer. The chain rule is how you compute how much each weight contributed to that loss, working backward from the output to the input.
What does the activation function actually do?
It introduces non-linearity. Without it, stacking any number of layers would still collapse to a single linear model, which is why a network with no activation function can't learn a non-linear boundary like XOR.
Why is XOR a good test problem for a neural network?
Because no straight line can separate its outputs, so logistic regression fails on it outright. Learning XOR proves the hidden layer is actually bending the decision boundary, not just memorizing a linear split.
How is this different from what PyTorch or TensorFlow do internally?
Mechanically, it's the same loop: forward pass, loss, backward pass, weight update. Production frameworks add mini-batches, adaptive optimizers like Adam, regularization, and autograd to compute the chain rule automatically, but the core mechanics shown here are what's running underneath.
Related articles
Where to next
New to PyTorch · step 2 of 7