# What's the bug in my implementation/understanding of backpropagation?

For learning purposes, I'm trying to implement a simple neural network with only linear layers followed by logistic activation. As far as I understand, the backpropagation algorithm exploits the structure of the chain rule in order to efficiently decompose the derivative of the loss (with respect to the weights and inputs) into factors corresponding to each layer.

So, I've written some training code corresponding to the following steps:

1. Initialize weights for each layer to random values centered around 0
2. Forward propagate a training example through the network, i.e.:

• For each layer in forward order, compute the output of the layer:

y = logistic(w * [1 | x])

where w is the matrix of weights of which the first column is the bias, y is the output of the layer, x is the input of the layer (i.e. the output of the previous layer or the training example for the first layer)

• Store the output of the last layer
3. Backward propagate error gradients, i.e.:
• For each layer in backward order:
• compute the gradient of the error in the layer's output with respect to its input
• compute the same error gradient with respect to the weights of the layer, which, by the chain rule, is the product of the gradient of the error w.r.t. the output and the gradient of the output w.r.t the weights
• adjust weights in the direction of this (second) antigradient (this is effectively stochastic gradient descent). The input for this procedure for the last layer is the gradient of the loss w.r.t. the output of the last layer; for the rest of the layers, it's the gradient w.r.t. the input of the next layer.
4. Repeat for another randomly-selected training example several times until convergence.

The code for the neural network is as follows. First, some helpers:

from scipy.special import erf
import numpy as np

# nonlinearity
def logistic(x):
return 0.5 * (1 + erf(x))

# Computes gradient (possibly, the Jacobian) of y with respect to xs
dx = 0.0001
return np.stack(
(y(xs + dxs) - y(xs - dxs)) / (2 * dx)
for dxs in dx * np.eye(len(xs.T))
).T

# Helper for Layer in the forward phase
def forward_pure(weights, inputs):
return logistic(weights @ inputs)

def loss(actual, predicted):
return np.sum((predicted - actual) ** 2)

# Gradient of loss w.r.t. predictions
return grad(lambda xs: loss(actual, xs), predicted)


Then, the class for a single layer:

class Layer:
def __init__(self, input_size, output_size, learning_rate=1):
# weights column 0 is bias
self.weights = np.random.normal(size=(output_size, input_size + 1))
self.learning_rate = learning_rate

def forward(self, inputs):
# Concatenate leading 1 to inputs for applying the bias too
self.inputs = np.array([1, *inputs])
self.outputs = forward_pure(self.weights, self.inputs)
return self.outputs

# Derivative of output w.r.t. weights
self.dout_dw = grad(lambda w: forward_pure(w, self.inputs), self.weights)
# Derivative of output w.r.t. inputs
# (input 0 is constant, so differentiating wrt. it is useless)
self.dout_din = grad(lambda i: forward_pure(self.weights, i), self.inputs)[:, 1:]

# We note that according to the chain rule:
# d[Error]/d[weights] = d[Error]/d[outputs] @ d[outputs]/d[weights]
# self.derror_dw = error_grad @ self.dout_dw
self.weights -= self.learning_rate * self.derror_dw

return self.dout_din


Finally, the class for a full multi-layer network:

class NN:
def __init__(self, layers):
self.layers = layers

def forward(self, inputs):
for layer in self.layers:
inputs = layer.forward(inputs)

self.outputs = inputs
return inputs

def backward(self, actual):
# derivative of error w.r.t. output

for layer in reversed(self.layers):
de_do = layer.backward(de_do)

return de_do


Although I'm using linear layers only, each followed by the logistic activation, I'm performing numeric differentiation instead of analytically computing the derivatives because I'd like to abstract away the details of the layers, I'm trying to understand backpropagation in general.

The problem is that I've got an interesting bug. When my output is a single value (i.e. each training example comes with a label that is a one-element array), then the network learns perfectly fine. However, when the output is a vector with more than one element, then it converges to the wrong values. For example, with a single-layer network, I can successfully predict the AND and the OR functions, however, trying to do both at the same time fails:

The training code:

def train(nn, num_epochs, Xs, ys):
for i in range(num_epochs):
samples = list(zip(Xs, ys))
np.random.shuffle(samples)

for X, y in samples:
pred = nn.forward(X)
err = nn.backward(y)

for X, y in zip(Xs, ys):
pred = nn.forward(X)
print(X, y, np.round(pred, 3))


And its invocation:

print("Learning logical AND")
train(
nn = NN(layers=[
Layer(input_size=2, output_size=1),
]),
num_epochs=10000,
Xs = np.array([
[0, 0],
[0, 1],
[1, 0],
[1, 1],
]),
ys = np.array([
[0],
[0],
[0],
[1],
])
)
print("Learning logical OR")
train(
nn = NN(layers=[
Layer(input_size=2, output_size=1),
]),
num_epochs=10000,
Xs = np.array([
[0, 0],
[0, 1],
[1, 0],
[1, 1],
]),
ys = np.array([
[0],
[1],
[1],
[1],
])
)
print("Learning both at the same time")
train(
nn = NN(layers=[
Layer(input_size=2, output_size=2),
]),
num_epochs=10000,
Xs = np.array([
[0, 0],
[0, 1],
[1, 0],
[1, 1],
]),
ys = np.array([
[0, 0],
[0, 1],
[0, 1],
[1, 1],
])
)


The output is:

Learning logical AND
[0 0] [0] [0.]
[0 1] [0] [0.003]
[1 0] [0] [0.003]
[1 1] [1] [0.997]
Learning logical OR
[0 0] [0] [0.003]
[0 1] [1] [0.998]
[1 0] [1] [0.998]
[1 1] [1] [1.]
Learning both at the same time
[0 0] [0 0] [0. 0.]
[0 1] [0 1] [0.127 0.279]
[1 0] [0 1] [1.    0.996]
[1 1] [1 1] [1. 1.]


For your reference, the full code can be found here so you don't have to copy-paste together the above snippets.