# How can I perform backpropagation directly in matrix form?

I had made a neural network library a few months ago, and I wasn't too familiar with matrices. So, instead of performing matrix dot products (between weights and inputs, then adding a bias matrix), I simply looped through each array element (for example, I looped through every weight), and as one would expect, it's slow.

I recently rewrote the library and employed matrix dot products (with numpy) for forward propagation, and it is way faster than the original library.

Original Forward Propagation Code:

x = 0
while x < layers: # Loop through layers
if x == 0: # If it is input layer
y = 0
while y < len(layervalues[x]): # Loop through neurons
layervalues[x][y] = sigmoid(np.sum(np.multiply(weight[x][y], input)) + bias[x][y]) # Multiplies inputs with weights and adds biases
y = y + 1
elif x == (layers - 1):
y = 0
while y < len(output): # Loop through neurons
output[y] = sigmoid(np.sum(np.multiply(weight[x][y], layervalues[(x - 1)])) + bias[x][y]) # Multiplies last hidden layer output with weights and adds biases
if pxl: # Used for debugging to prevent the output from passing through the sigmoid function
output[y] = np.sum(np.multiply(weight[x][y], layervalues[(x - 1)])) + bias[x][y]
y = y + 1
else:
y = 0
while y < len(layervalues[x]): # Loop through neurons
layervalues[x][y] = sigmoid(np.sum(np.multiply(weight[x][y], layervalues[(x - 1)])) + bias[x][y]) # Multiplies previous layer outputs with weights and adds biases
y = y + 1
x = x + 1


ip = np.asarray(ip)
ip = np.reshape(ip, (len(ip), 1))
for x in range(len(self.struct) - 1):
if x == 0:
food = np.dot(self.weights[x], ip)
else:
food = np.dot(self.weights[x], food)
if backprop:
np.reshape(self.logz[x], (1, len(self.logz[x])))
if backprop == True:
return food
else:
return np.reshape(food, (1, len(food))).tolist()[0]


Now that forward propagation works comparatively well, I'm trying to attempt the same for backpropagation. I need to do this because the backpropagation code takes an unacceptably long time to generate gradients for networks that have many neurons in their layers.

Here's my current backpropagation code:

netop = self.feedForward(ip, True)
output = np.asarray(output)
output = np.reshape(output, (len(output), 1))
delzdela = [None] * (len(self.struct) - 1)
weightgrad = [None] * (len(self.struct) - 1)
biasgrad = [None] * (len(self.struct) - 1)
for x in range(len(self.weights)): # Loop to calculate gradients of neurons of hidden layers
weightgrad[x] = np.zeros((self.struct[x + 1], self.struct[x]), dtype=np.float64)
biasgrad[x] = np.zeros((self.struct[x + 1], 1), dtype=np.float64)
for x in range(len(self.struct) - 1):
if x == 0:
delzdela[-(x + 1)] = np.dot(2, np.subtract(netop, output))
else:
delzdela[-(x + 1)] = np.zeros((self.struct[-(x + 1)], 1))
for y in range(self.struct[-(x + 1)]):
for z in range(self.struct[-x]):
delzdela[-(x + 1)][y] += derivative(self.logz[-x][z], self.actfunc[-x]) * delzdela[-x][z] * self.weights[-x][z][y]
for x in range(len(self.struct) - 1): # Calculating gradients of weights
for y in range(self.struct[x]):
for z in range(self.struct[x + 1]):
if x == 0:
weightgrad[x][z][y] = ip[y] * derivative(self.logz[x][z], self.actfunc[x]) * delzdela[x][z]
else:
weightgrad[x][z][y] = activation(self.logz[x - 1][z], self.actfunc[x - 1]) * derivative(self.logz[x][z], self.actfunc[x]) * delzdela[x][z]
for x in range(len(self.struct) - 1): # Calculating gradients of biases
for y in range(self.struct[x + 1]):
biasgrad[x][y] = derivative(self.logz[x][y], self.actfunc[x]) * delzdela[x][y]


This code takes a training example (ip) and its label (output) and returns the gradient for this single training example.

As you can see, I am looping through each index and dimension of the arrays, which is very slow.

My question is how I can directly use matrix operations to perform backpropagation, instead of looping through every index and performing operations on individual numbers.

You should avoid explicit for loops in Python, whenever possible. For that, you should use the power of broadcasting and vectorization of Python NumPy. You can refer to this short video:

For the vectorized application of backpropagation algorithm, have a look at (5:09), the algorithm at the right-hand side; you can apply this inside an (inevitable) for loop for each layer l:

Those videos are from Coursera course "Deep Learning". In the lab assigments, I became able to write forward and backpropagation algorithms by only 6 lines each by vectorizing the flows through the layers. The course instructor, Andrew Ng really forces the importance of this; Python NumPy can use parallellization structures of the CPUs and GPUs, you get both compitationally efficient (nearly about 300x faster in the course experiment) and you can write those in with few lines of codes, getting rid of the scary explicit for loops.

I wish I could share my own code here, yet it is prohibited by the Honor Code of Coursera, and Andrew Ng describes those basic concepts better than any way I can.

Here is how I did it:

public final Tensor3 backprop(final double learningRate, final Vector input, final Vector ideal) {
run(input);

final Matrix[] weightDeltaTensor = new Matrix[weightTensor.dimension];

Vector errorOutputDeriv = activations[layers - 1].elementOperation(ideal, MSE_DERIV);

for (int i = layers - 1; i >= 1; i--) {
final Vector sigmoidDeriv = activations[i].transform(FuncDerivPair.SIGMOID.partialDerivative);
final Matrix weightDeltas;

if (i == layers - 1) {
weightDeltas = errorInputDeriv.outerProduct(activations[i - 1]);
} else {
weightDeltas = errorInputDeriv.removeLastElement().outerProduct(activations[i - 1]);
}

weightDeltaTensor[i - 1] = weightDeltas;

final Matrix currentWeights = weightTensor.getLayer(i - 1);

if (i == layers - 1) {
errorOutputDeriv = currentWeights.transpose().multiply(errorInputDeriv);
} else {
errorOutputDeriv = currentWeights.transpose().multiply(errorInputDeriv.removeLastElement());
}
}

return new Tensor3(weightDeltaTensor);
}


It's in Java but you should be able to make sense of it. activations is a Vector[] and run(input) populates it with activation vectors for each layer, including the input vector. There is only one explicit for loop, and that is merely used to access each layer. weightTensor is really just an array of weight matrices, with biases in the final column of each matrix. This is the reason for if (i == layers - 1) and removeLastElement(), because the activation vectors for all but the output layer have an extra 1 appended to trigger the bias.

Vector.transform() takes a 1-argument function and returns a new vector with that function applied to each element. Similarly, Vector.elementOperation() takes a 2-argument function and another vector; it returns a new vector with the function applied element-wise to both vectors.

MSE_DERIV and FuncDerivPair.SIGMOID.partialDerivative can be replaced with functions of your choice. This code just assumes that every layer is using the sigmoid function, but you could extend this so that layers can have different activation functions.