Introduction
The previous post ended with a problem expressed in code. I built mlp(), a two-layer network that solves XOR when every weight is set by hand. Then I tried to train it with train(), the perceptron learning rule, and it failed. The learning rule computes error = target - output, but hidden neurons have no targets. The output neuron knows it got XOR wrong. The hidden neurons do not know what they should have produced. That is the credit assignment problem.
Rosenblatt spent the last decade of his life trying to solve it. He built multi-layer networks with random hidden projections and trainable output layers, but he never found a way to propagate learning into the hidden connections. Minsky and Papert proved the single-layer model was fundamentally limited. The field froze for over fifteen years.
The solution was hiding in calculus the entire time. Three changes turn a single-layer perceptron into a trainable deep network:
- Replace the step function with a smooth, differentiable activation
- Define a continuous measure of how wrong the output is (a loss function)
- Use the chain rule to propagate that error backward through every layer
Each change is simple on its own. Together, they solve the problem that killed neural network research and launch the field that eventually produces modern deep learning. Let me build each piece.
The Step Function Problem
The perceptron uses a step function: output 1 if z >= 0, output 0 otherwise. The output is binary. There is nothing in between.
For a single layer, this works fine. The perceptron learning rule does not need to know how wrong the output is, only that it is wrong. The update w = w + eta * (t - y) * x uses the raw error (which is -1, 0, or +1) as a correction signal. That is enough when every weight connects directly to the output.
For hidden layers, it is not enough. I need to answer a different question: "How much should I change this hidden weight to reduce the error at the output?" That question requires a derivative. If I change a hidden weight by a small amount, how much does the output change? For the step function, the answer is almost always zero. The output is flat (derivative zero) everywhere except at z = 0, where it jumps discontinuously (derivative undefined). A tiny change in any hidden weight produces either no change in the output or an infinite change. There is no gradient to follow.
Step function:
y = 1 if z >= 0, else 0
Derivative:
dy/dz = 0 everywhere except z = 0
dy/dz = undefined at z = 0
The fix: replace the step function with something smooth. Something whose output changes gradually as the input changes, so that the derivative exists everywhere and tells me which direction to push.
Activation Functions
Sigmoid: The Historical Fix
The function that made deep learning possible is the logistic sigmoid:
sigmoid(z) = 1 / (1 + exp(-z))
It looks like the step function with the hard edge smoothed out. For large positive z, the output approaches 1. For large negative z, it approaches 0. At z = 0, the output is exactly 0.5. The transition between 0 and 1 is smooth, and the derivative exists everywhere:
sigmoid'(z) = sigmoid(z) * (1 - sigmoid(z))
The derivative has a convenient property: it can be computed from the output alone. If a = sigmoid(z), then sigmoid'(z) = a * (1 - a). No need to store z. This matters for efficiency during backpropagation.
The maximum derivative is 0.25, at z = 0. For large |z|, the derivative approaches zero. The function saturates: when the neuron is confidently on or confidently off, the gradient vanishes. This will become a serious problem in deep networks, but for now it is enough that the derivative exists at all.
This is the activation function that Rumelhart, Hinton, and Williams used in 1986 to demonstrate backpropagation on multi-layer networks. It is the function that ended the AI winter.
ReLU: The Modern Default
Sigmoid's saturation problem motivated the search for alternatives. The breakthrough came from a deceptively simple function:
relu(z) = max(0, z)
Zero for negative inputs. Identity for positive inputs. The derivative is 0 or 1. No saturation for positive values, so gradients flow without shrinking. Nair and Hinton demonstrated in 2010 that ReLU significantly improved training of deep networks.
The tradeoff: neurons with negative pre-activations have zero gradient. If a neuron's weighted sum goes negative and stays there, it stops learning entirely. These "dead neurons" are the cost of ReLU's simplicity.
GELU: The Transformer's Choice
gelu(z) = z * Phi(z)
Where Phi(z) is the standard Gaussian cumulative distribution function. GELU smoothly gates the input by its own probability of being positive under a Gaussian distribution. For large positive values, GELU approaches the identity (like ReLU). For large negative values, it approaches zero (like ReLU). But the transition is smooth, with a small negative dip near z = -0.75. This is the activation used in GPT, BERT, and most modern transformers.
Figure 1: Four activation functions on the same axes. The step function (dashed gray) is the perceptron's original activation, binary with zero derivative. Sigmoid (cyan) was the first differentiable replacement, enabling backpropagation in 1986. ReLU (green) solved the vanishing gradient problem with a constant derivative of 1 for positive inputs. GELU (yellow) is the smooth variant used in modern transformers.
The MLP: Architecture and Forward Pass
With a differentiable activation function, I can revisit the multilayer perceptron from the previous post. The architecture is the same: input layer, one or more hidden layers, output layer. Each layer computes a weighted sum plus bias, then applies an activation function. The difference is that the activation is now sigmoid instead of step, and the output is a continuous value between 0 and 1 instead of a hard 0 or 1.
For a two-layer network (one hidden layer, one output layer), the forward pass is:
Hidden layer: z1 = W1 @ x + b1 a1 = sigmoid(z1)
Output layer: z2 = W2 @ a1 + b2 y = sigmoid(z2)
The complete function is a composition: f(x) = sigmoid(W2 @ sigmoid(W1 @ x + b1) + b2). Each layer transforms its input. The composition of transformations can represent nonlinear decision boundaries, exactly the kind the perceptron could not express.
How powerful is this? In 1989, George Cybenko proved the universal approximation theorem: a single hidden layer with enough neurons and a non-constant, bounded, continuous activation function (like sigmoid) can approximate any continuous function on a compact set to arbitrary accuracy. Kurt Hornik extended this in 1991 to show the result holds for a broader class of activations. The MLP is not just more powerful than a perceptron. It is a universal function approximator.
The catch is "enough neurons." A single hidden layer can represent anything in theory, but in practice it may need an astronomically large number of neurons. Deep networks, those with many layers, can often represent the same functions far more efficiently. Depth is not about what you can represent, it is about how efficiently you can represent it.
Figure 2: A 2-2-1 MLP for XOR. Two inputs feed into two hidden neurons (each applying sigmoid), which feed into one output neuron (also sigmoid). The weight matrices W1 (2x2) and W2 (1x2) are the learnable parameters. The forward pass flows left to right: inputs are transformed by the hidden layer into a new representation, then the output layer classifies that representation.
The Learning Problem: Loss Functions
Now I have a differentiable architecture that can represent nonlinear functions. The next question: how do I measure "wrong" precisely enough for calculus to work with?
The perceptron used error = target - output, which takes values in {-1, 0, 1}. With sigmoid outputs between 0 and 1, I need a continuous measure of distance between the output and the target.
Mean Squared Error
The simplest option: square the difference.
L = (target - output)^2
Derivative with respect to the output: dL/dy = -2(target - output). Simple and intuitive. But for classification with sigmoid, MSE has a problem. When the output is confidently wrong (sigmoid near 0 but target is 1, or near 1 but target is 0), sigmoid's derivative is tiny, so the combined gradient is small. The model learns slowly exactly when it is most wrong.
Cross-Entropy: The Right Loss for Classification
L = -[t * log(y) + (1-t) * log(1-y)]
Where t is the target (0 or 1) and y is the sigmoid output. When t = 1, the loss is -log(y): the closer y is to 1, the smaller the loss. When y approaches 0 (confidently wrong), the loss goes to infinity. The loss heavily penalizes confident wrong answers.
The derivative of cross-entropy with respect to the sigmoid output is (y - t) / (y * (1 - y)). Combined with sigmoid's derivative y * (1 - y), the compound gradient for the output layer simplifies to:
dL/dz = y - t
The y * (1 - y) terms cancel. There is no saturation. The gradient is simply the difference between the output and the target. This is why cross-entropy became the default loss for classification: the gradient signal stays strong regardless of how wrong the output is.
The Chain Rule: Why Calculus Solves Credit Assignment
This is the conceptual core of the post. Everything before this was setup. The credit assignment problem, the one that stumped Rosenblatt for a decade and froze the field, reduces to a single idea from calculus.
The question is: given a loss L computed at the output, how much is each weight in each layer responsible? Formally, I need dL/dw for every weight w in every layer.
Consider a two-layer network. The loss depends on the output y, which depends on the output pre-activation z2, which depends on the hidden activations a1, which depends on the hidden pre-activations z1, which depends on the first-layer weights W1. It is a chain of dependencies:
W1 -> z1 -> a1 -> z2 -> y -> L
The chain rule says: if L depends on y which depends on z2 which depends on a1 which depends on z1 which depends on W1, then:
dL/dW1 = dL/dy * dy/dz2 * dz2/da1 * da1/dz1 * dz1/dW1
Each factor in that product is a local derivative. dL/dy is the loss derivative. dy/dz2 is the activation derivative. dz2/da1 is just the output weight matrix. da1/dz1 is the hidden activation derivative. dz1/dW1 is just the input. Every factor can be computed from values already available during the forward pass.
The credit assignment problem was never about missing information. The forward pass computes everything needed. The chain rule tells you how to use it. The line that solves the problem is:
dL/da1 = W2^T @ dL/dz2
The output weight matrix W2, transposed, distributes the output error back to each hidden neuron in proportion to that neuron's connection strength. A hidden neuron connected to the output by a large positive weight gets a large share of the blame. A neuron connected by a small weight gets little blame. A neuron connected by a negative weight gets blame in the opposite direction.
This is exactly what Rosenblatt could not do. His learning rule needed a target for each neuron. The chain rule replaces targets with gradients. Hidden neurons do not need targets because the chain rule computes their contribution to the error directly from the network's own weights.
Figure 3: The computational graph of a two-layer network. Top row: the forward pass computes values from left to right. Bottom row: the backward pass propagates gradients from right to left by multiplying local derivatives at each step. The chain rule connects the output error to every weight in the network. The key step is W2^T, which distributes blame from the output to hidden neurons in proportion to their connection strength.
Backpropagation: The Algorithm
History
Paul Werbos described backpropagation in his 1974 PhD thesis Beyond Regression. The work went largely unnoticed. In 1986, David Rumelhart, Geoffrey Hinton, and Ronald Williams published Learning Representations by Back-Propagating Errors in Nature, demonstrating backpropagation on multi-layer networks and showing that hidden layers learn useful intermediate representations. This paper ended the AI winter and revived neural network research.
The algorithm has two phases: a forward pass that computes the output, and a backward pass that computes the gradient of the loss with respect to every weight.
Step by Step
1. FORWARD PASS
For each layer (input to output):
z = W @ a_prev + b (weighted sum)
a = sigmoid(z) (activation)
Store z and a for use in backward pass
2. COMPUTE LOSS
L = -[t * log(y) + (1-t) * log(1-y)] (cross-entropy)
3. BACKWARD PASS
Output layer:
delta = y - t (cross-entropy + sigmoid)
dL/dW = delta @ a_prev^T (weight gradient)
dL/db = delta (bias gradient)
For each hidden layer (output to input):
delta = (W_next^T @ delta_next) * sigmoid'(a) (propagate + activate)
dL/dW = delta @ a_prev^T (weight gradient)
dL/db = delta (bias gradient)
4. UPDATE
For every weight and bias:
w = w - eta * dL/dw
b = b - eta * dL/db
The key line is delta = (W_next^T @ delta_next) * sigmoid'(a). This propagates the error from the next layer back to the current layer. The transposed weight matrix distributes blame. The activation derivative gates it. This is the credit assignment step that Rosenblatt could never perform.
Worked Example
I will trace one forward and backward pass through a 2-2-1 network with specific weights. This is the XOR architecture from the perceptron post, but now with sigmoid activations and backpropagation.
Setup:
W1 = [[0.5, -0.3], b1 = [0.1, -0.2]
[0.2, 0.8]]
W2 = [[0.6, -0.4]] b2 = [0.3]
Input: x = [1, 0] Target: t = 1
Forward pass:
Hidden layer:
z1[0] = 0.5*1 + (-0.3)*0 + 0.1 = 0.6
z1[1] = 0.2*1 + 0.8*0 + (-0.2) = 0.0
a1[0] = sigmoid(0.6) = 0.6457
a1[1] = sigmoid(0.0) = 0.5000
Output layer:
z2 = 0.6*0.6457 + (-0.4)*0.5000 + 0.3 = 0.4874
y = sigmoid(0.4874) = 0.6195
Loss:
L = -[1*log(0.6195) + 0*log(0.3805)] = 0.4789
The network outputs 0.6195. The target is 1. The loss is 0.4789. Now I propagate the error backward.
Backward pass:
Output delta (cross-entropy + sigmoid simplification):
delta2 = y - t = 0.6195 - 1 = -0.3805
Output weight gradients:
dL/dW2[0] = delta2 * a1[0] = -0.3805 * 0.6457 = -0.2457
dL/dW2[1] = delta2 * a1[1] = -0.3805 * 0.5000 = -0.1903
dL/db2 = -0.3805
Now the key step. Propagate the error to the hidden layer:
Error distributed to hidden neurons:
dL/da1[0] = W2[0] * delta2 = 0.6 * (-0.3805) = -0.2283
dL/da1[1] = W2[1] * delta2 = (-0.4) * (-0.3805) = 0.1522
Hidden deltas (multiply by sigmoid derivative):
delta1[0] = -0.2283 * a1[0]*(1 - a1[0]) = -0.2283 * 0.2288 = -0.0522
delta1[1] = 0.1522 * a1[1]*(1 - a1[1]) = 0.1522 * 0.2500 = 0.0381
Hidden weight gradients:
dL/dW1[0][0] = delta1[0] * x[0] = -0.0522 * 1 = -0.0522
dL/dW1[0][1] = delta1[0] * x[1] = -0.0522 * 0 = 0.0000
dL/dW1[1][0] = delta1[1] * x[0] = 0.0381 * 1 = 0.0381
dL/dW1[1][1] = delta1[1] * x[1] = 0.0381 * 0 = 0.0000
Every weight in the network now has a gradient. The output weights get gradients directly from the output error. The hidden weights get gradients via the chain rule, with W2^T distributing blame and sigmoid' gating it. The weight update w = w - eta * dL/dw nudges every weight in the direction that reduces the loss. No hand-design. No per-neuron targets. Just the chain rule.
Figure 4: Backpropagation in action. Phase 1 (green/cyan): the forward pass computes values from inputs through hidden to output. Phase 2 (red): the backward pass propagates gradients from the loss back through the network. The output error y - t is distributed to hidden neurons via W2^T, then to input weights via the chain rule. Every weight in every layer gets a gradient and updates accordingly.
Gradient Descent
Backpropagation computes the gradient dL/dw for every weight. The gradient points in the direction of steepest increase of the loss. To reduce the loss, I move in the opposite direction:
w = w - eta * dL/dw
This is gradient descent. The learning rate eta controls step size. But there are several variants, and the choice matters.
Vanilla gradient descent computes the gradient over the entire training set, then takes one step. Stable but slow. For large datasets, computing the full gradient before each step is prohibitive.
Stochastic gradient descent (SGD) computes the gradient on a single random training example, then takes a step. Noisy but fast. Each step is cheap, and the noise actually helps escape shallow local minima. This is what makes training on millions of examples feasible.
Mini-batch SGD is the practical compromise: compute the gradient on a small batch of examples (typically 32 to 512), then step. It reduces noise while keeping updates frequent. This is the standard in practice.
Momentum adds inertia to the updates:
v = beta * v + dL/dw
w = w - eta * v
The velocity v accumulates a running average of past gradients. When gradients point in the same direction across multiple steps, the velocity builds up and the steps get larger. When gradients oscillate, the velocity dampens them. Like a ball rolling downhill: it accelerates on consistent slopes and resists zig-zagging.
Adam (Kingma and Ba, 2015) adapts the learning rate for each weight individually. It maintains running estimates of the first moment (mean) and second moment (uncentered variance) of the gradient, with bias correction for the early steps. Weights with consistently large gradients get smaller learning rates. Weights with small or noisy gradients get larger ones. Adam is the default optimizer for most modern deep learning.
Figure 5: A loss landscape with four optimizer paths. Vanilla gradient descent (gray) takes slow, steady steps. SGD (red, dashed) is noisy, zig-zagging toward the minimum. Momentum (cyan) follows a smooth curve with slight overshoot, then corrects. Adam (green) takes the most efficient path by adapting its learning rate per weight.
The Vanishing Gradient Problem
Backpropagation solved the credit assignment problem. But it introduced a new one.
The chain rule multiplies local derivatives along the path from output to input. With sigmoid activations, each local derivative is sigmoid'(z) = sigmoid(z) * (1 - sigmoid(z)), which has a maximum value of 0.25 at z = 0. In a 10-layer network, the gradient at the first layer passes through 10 sigmoid derivatives:
0.25 * 0.25 * 0.25 * ... (10 times) = 0.25^10 = 0.00000095
The gradient vanishes. The first layers receive a gradient that is a millionth of the original signal. They barely learn. This is why deep networks with sigmoid activations were so hard to train, even after backpropagation was discovered.
The opposite problem also exists: if weights are large, the chain of multiplications can grow exponentially. This is the exploding gradient problem. Gradients become enormous, weight updates become wild, and training diverges.
ReLU largely solved the vanishing gradient problem. For positive inputs, relu'(z) = 1. The gradient passes through without shrinking. A 10-layer network with ReLU activations and positive pre-activations has a gradient multiplier of 1^10 = 1. No vanishing. The cost is dead neurons: if a neuron's pre-activation goes negative, its gradient is zero and it stops learning entirely.
Weight Initialization
The other piece of the puzzle: how weights are initialized matters enormously for gradient flow. If initial weights are too large, activations saturate and gradients vanish (sigmoid) or explode. Too small, and signals die before reaching deeper layers.
Xavier/Glorot initialization (2010) sets weights to maintain variance across layers with sigmoid or tanh activations:
w ~ Normal(mean=0, std=sqrt(2 / (n_in + n_out)))
He initialization (2015) adjusts for ReLU, which zeros out half the distribution:
w ~ Normal(mean=0, std=sqrt(2 / n_in))
Both strategies keep the variance of activations and gradients roughly constant as signals pass through the network. Without them, deep networks fail to train. With them, networks of 10, 50, or 100 layers become feasible.
Figure 6: The vanishing gradient in a 5-layer sigmoid network. The gradient at the output is 1.0. After passing through each sigmoid derivative (max 0.25), it shrinks by at least 4x per layer. By layer 1, the gradient is 0.004, essentially zero. The first layers barely learn. ReLU fixes this: its derivative is 1 for positive inputs, so gradients pass through without shrinking.
Perceptron vs. MLP + Backpropagation
| Property | Perceptron (1958) | MLP + Backprop (1986) |
|---|---|---|
| Activation | Step function (binary) | Sigmoid, ReLU (continuous) |
| Learning rule | w += eta * (t-y) * x | w -= eta * dL/dw (chain rule) |
| Layers trainable | Output only | All layers |
| Decision boundary | Single hyperplane | Arbitrary nonlinear |
| Loss function | None (discrete error) | Cross-entropy, MSE |
| Credit assignment | Impossible for hidden layers | Solved by chain rule |
| XOR | Cannot learn | Learns automatically |
Backpropagation in Code
The code follows the same progression from the previous posts. neuron() became perceptron() became learn() and train(). Now I replace the step function with sigmoid, replace the perceptron update rule with backpropagation, and the MLP learns from data.
import math
import random
def sigmoid(z):
"""Logistic sigmoid: smooth, differentiable replacement for step."""
return 1 / (1 + math.exp(max(-500, min(500, -z))))
def forward(x, weights, biases):
"""Forward pass. Returns list of activations for every layer."""
activations = [list(x)]
for W, b in zip(weights, biases):
z = [sum(w * a for w, a in zip(row, activations[-1])) + bi
for row, bi in zip(W, b)]
activations.append([sigmoid(zi) for zi in z])
return activations
def backward(activations, weights, target):
"""Backward pass. Returns weight and bias gradients for every layer."""
output = activations[-1]
delta = [o - t for o, t in zip(output, target)]
w_grads, b_grads = [None] * len(weights), [None] * len(weights)
for layer in reversed(range(len(weights))):
a_prev = activations[layer]
w_grads[layer] = [[d * a for a in a_prev] for d in delta]
b_grads[layer] = list(delta)
if layer > 0:
a = activations[layer]
delta = [sum(weights[layer][k][j] * delta[k]
for k in range(len(delta)))
* a[j] * (1 - a[j])
for j in range(len(a))]
return w_grads, b_grads
def train_mlp(data, targets, sizes, eta=2.0, epochs=2000):
"""Train an MLP with backpropagation."""
random.seed(42)
weights, biases = [], []
for i in range(len(sizes) - 1):
s = (2 / (sizes[i] + sizes[i+1])) ** 0.5
weights.append([[random.gauss(0, s) for _ in range(sizes[i])]
for _ in range(sizes[i+1])])
biases.append([0.0] * sizes[i+1])
for epoch in range(epochs):
total_loss = 0
for x, t in zip(data, targets):
tv = [t] if isinstance(t, (int, float)) else t
acts = forward(x, weights, biases)
wg, bg = backward(acts, weights, tv)
for l in range(len(weights)):
for j in range(len(weights[l])):
for k in range(len(weights[l][j])):
weights[l][j][k] -= eta * wg[l][j][k]
biases[l][j] -= eta * bg[l][j]
y = max(1e-15, min(1-1e-15, acts[-1][0]))
total_loss += -(tv[0]*math.log(y) + (1-tv[0])*math.log(1-y))
if epoch == 0 or (epoch+1) % 500 == 0:
print(f"Epoch {epoch+1:4d} loss: {total_loss/len(data):.4f}")
return weights, biases
Training on XOR:
data = [[0,0], [0,1], [1,0], [1,1]]
targets = [0, 1, 1, 0]
weights, biases = train_mlp(data, targets, sizes=[2, 2, 1])
Epoch 1 loss: 1.1687
Epoch 500 loss: 0.0085
Epoch 1000 loss: 0.0018
Epoch 1500 loss: 0.0010
Epoch 2000 loss: 0.0007
for x, t in zip(data, targets):
y = forward(x, weights, biases)[-1][0]
print(f" {x} -> {y:.4f} (target: {t})")
[0, 0] -> 0.0006 (target: 0)
[0, 1] -> 0.9994 (target: 1)
[1, 0] -> 0.9994 (target: 1)
[1, 1] -> 0.0010 (target: 0)
XOR solved. Not hand-wired. Learned from data. The loss drops from 1.17 to 0.0007. The outputs are within a fraction of a percent of the targets. Compare this to the perceptron post, where train() ran for 1000 epochs on XOR and never converged, oscillating forever because a single hyperplane cannot separate the classes.
The difference is three functions: sigmoid() provides the smooth activation that makes calculus possible, backward() uses the chain rule to compute gradients for every weight in every layer, and train_mlp() uses those gradients to update all weights simultaneously. The credit assignment problem, the problem that killed Rosenblatt's research and froze neural networks for fifteen years, reduces to the chain rule applied to a composition of differentiable functions.
Key Takeaways
I started where the perceptron post left off: a multi-layer architecture that could represent XOR but could not learn it. The gap was the credit assignment problem. Hidden neurons have no targets, so the perceptron learning rule cannot update their weights.
The solution required three changes. First, replace the step function with sigmoid, a smooth activation whose derivative exists everywhere and can be computed from its own output. Second, define cross-entropy as the loss function, giving a continuous, gradient-friendly measure of error that does not saturate when combined with sigmoid. Third, apply the chain rule to propagate the error backward through every layer. The transposed weight matrix W2^T distributes blame from the output to hidden neurons in proportion to their connection strength. The activation derivative gates the signal. The result is a gradient for every weight in every layer.
Backpropagation came with its own challenges. Sigmoid's maximum derivative of 0.25 causes gradients to vanish exponentially in deep networks. ReLU fixed that with a constant derivative of 1 for positive inputs. Xavier and He initialization kept gradients flowing through many layers. Adam adapted learning rates per parameter. Each fix enabled deeper networks, which enabled more powerful representations.
The progression so far: McCulloch-Pitts gave us a neuron that computes but cannot learn. Rosenblatt's perceptron added learning but was limited to linear boundaries. Backpropagation solved the credit assignment problem and enabled multi-layer networks to learn arbitrary nonlinear functions. But there is a missing capability. These networks process fixed-size inputs. Every example has to be the same shape: a vector of numbers with a predetermined length. Real-world data, language especially, does not work that way. A sentence can be three words or three hundred. The order matters. The meaning of a word depends on which words came before it. How do you handle sequences? That question drives the next post.