Training a Multilayer Neural Network With NumPy: Backpropagation From Scratch

In the previous post, we passed an input through multiple layers and got a prediction. The weights were random, though, so that prediction had not learned anything. Now we will teach a small network to change its own weights: calculate how far its prediction is from a target, send that error backward through the layers, and take a small step in the direction that reduces the error.

This is backpropagation combined with stochastic gradient descent. We will build both pieces with NumPy and ordinary Python loops, so we can see what each derivative does.

Create a small synthetic problem

We start with ten examples, each containing three numeric features. The values are divided by 10, bringing them into a smaller range. Smaller inputs help keep the initial weighted sums manageable when our network starts with random parameters.

raw_inputs = np.array([
    [3, 1, 7],
    [4, 2, 5],
    [1, 0, 3],
    [6, 4, 2],
    [8, 5, 1],
    [2, 3, 9],
    [7, 7, 4],
    [0, 4, 6],
    [9, 1, 8],
    [5, 0, 8],
], dtype=float)

inputs = raw_inputs / 10

Instead of hand-labeling every row, we use a simple rule to generate its target. First, a weighted sum uses the known weights [4, -2, 0] and bias -1. Then sigmoid squashes that result into the interval from 0 to 1.

TRUE_WEIGHTS = np.array([4.0, -2.0, 0.0])
TRUE_BIAS = -1.0

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

targets = sigmoid(inputs @ TRUE_WEIGHTS + TRUE_BIAS)

The @ operator performs a matrix-vector product. Because inputs has shape (10, 3) and TRUE_WEIGHTS has shape (3,), the result has one weighted sum per row. The third true weight is zero, so the third feature does not affect these targets. This rule is our teacher: training examples provide inputs and answers, but the network will not be given the teacher's weights.

Choose a network and a loss

Our network has layer sizes [3, 8, 8, 1]: three input features, two hidden layers of eight neurons each, and one output neuron. Each layer's weight matrix has shape (n_out, n_in), and each bias vector has shape (n_out,). A seeded random generator makes the starting parameters reproducible.

Hidden neurons use the ReLU activation max(0, z). The output neuron is linear, meaning it returns its weighted sum without an activation. We compare that prediction with its target using half the squared error:

loss = 0.5 * (target_value - predicted_value) ** 2

The factor of 0.5 is convenient because it cancels the 2 that appears when differentiating a squared value. The derivative of this loss with respect to a linear output is simply predicted_value - target_value.

Even though the targets were made with sigmoid and lie between 0 and 1, the network output here is linear and is not constrained to that interval. This example uses squared-error regression to approximate the teacher's outputs; if the model must produce probabilities, its output activation and loss should be chosen for that classification task.

Remember each layer's inputs

To send an error backward, we need the inputs that each neuron saw during the forward pass. The network resets input_values for each training example and saves each layer's input vector for each neuron. For the first layer, that vector is the original row; for later layers, it is the preceding layer's outputs.

def forward_pass(self, single_input):
    self.input_values = []
    last_output = []

    for i, layer_weights in enumerate(self.weights):
        outputs = []
        input_values = []
        input_value = single_input if i == 0 else last_output

        for j, neuron_weights in enumerate(layer_weights):
            input_values.append(input_value)
            z = np.dot(input_value, neuron_weights) + self.bias[i][j]
            outputs.append(self.activation_function(
                z, is_output_layer=(i == len(self.weights) - 1)
            ))

        last_output = outputs
        self.input_values.append(input_values)

    return last_output

This is the same layered calculation as before, with one extra job: cache the inputs needed later for the weight gradients. The output is a list with one value because this network has one output neuron.

Calculate the output error

Backpropagation begins at the output. For a neuron, its delta is the derivative of the loss with respect to that neuron's weighted sum, z. At the linear output, the activation derivative is 1, so the delta is just the prediction minus the target:

output_delta = predicted_value[0] - target_value

A positive delta means the prediction is too high; a negative delta means it is too low. The sign will determine which direction gradient descent changes the parameters.

Send the error through hidden layers

A hidden neuron affects the loss through the neurons in the next layer. Its downstream error is the sum of each next-layer delta times the weight connecting this neuron to that next neuron. Then we multiply by the derivative of its activation. ReLU's derivative is 1 when its weighted sum is positive and 0 otherwise.

downstream_error = 0.0
for k, next_neuron_weights in enumerate(self.weights[i + 1]):
    downstream_error += (
        deltas[i + 1][k] * next_neuron_weights[j]
    )

# Recompute this neuron's weighted sum from cached inputs and weights.
relu_derivative = 1.0 if z > 0 else 0.0
neuron_delta = downstream_error * relu_derivative

The loop processes hidden layers in reverse order, starting just before the output and moving toward the input. At each layer it computes all the deltas using the next layer's weights. This is the chain rule in code: a hidden neuron's influence is the combined influence of the neurons after it, scaled by the local activation derivative.

When z is zero or negative, ReLU's derivative is zero in this implementation, so that neuron receives no gradient for this example. This is why a hidden activation matters: it both transforms the signal in the forward direction and controls how the error flows backward.

Turn deltas into parameter updates

Once every layer's deltas have been calculated, we can update parameters. For a neuron's weight connected to input x_k, the gradient is delta * x_k. The bias gradient is just delta, since the bias is added directly to the weighted sum. Gradient descent subtracts the learning rate multiplied by each gradient.

gradient = neuron_delta * neuron_inputs[k]
self.weights[i][j][k] -= self.learning_rate * gradient
self.bias[i][j] -= self.learning_rate * neuron_delta

It is important that the code calculates all deltas before changing any weights. A hidden-layer delta depends on the weights from that layer to the next; updating those weights too early would mix old and new parameter values in one backward pass.

The implementation performs one update after each example. This is called stochastic gradient descent (often, for a small dataset, mini-batch gradient descent is used in practice). The order of examples can therefore affect the path the parameters take while learning.

Train and inspect the result

Training repeats the data pass for 1,000 epochs. An epoch is one visit to every example. We accumulate the loss to report an average, then call backward_pass for each prediction to update the model.

for epoch in range(1000):
    total_loss = 0.0

    for x_i, target_value in zip(inputs, targets):
        predicted_value = nn.forward_pass(x_i)
        total_loss += nn.loss(predicted_value[0], target_value)
        nn.backward_pass(predicted_value, target_value)

    if (epoch + 1) % 50 == 0:
        print(f"Epoch {epoch + 1:3d} | "
              f"Loss: {total_loss / len(inputs):.6f}")

With the supplied data, seed, architecture, and learning rate, the reported average loss falls from about 0.0170 after 50 epochs to about 0.000246 after 1,000 epochs. For example, one target of 0.9168 is predicted as about 0.9554, and a target of 0.1419 is predicted as about 0.1512. Small differences in the final values can occur if the code or training order is changed.

The loss goes down because the weights and biases are being adjusted to fit the examples. It does not prove the network has discovered the exact teacher rule or will perform well on new data. There are only ten examples, and the reported loss is measured on those same examples used for training. A separate validation set would be needed to check generalization.

What to improve next

This implementation is deliberately explicit, but it has constraints. It assumes a single output neuron because the output delta is initialized with one value. It handles one example at a time, stores the same layer input once per neuron, and recomputes weighted sums during backpropagation instead of caching them. Its standalone relu helper is also unused because ReLU is written directly inside activation_function.

Those choices keep the chain rule visible, but a more flexible version would cache every neuron's z, use NumPy arrays and matrix operations for whole layers, and support batches and multiple outputs. The central idea remains the same: make a prediction, measure its error, propagate that error backward through the network, and use the resulting gradients to improve the parameters.