How to perform back propagation with different sized layers?

Question

I'm developing my first neural network, using the well known MNIST database of handwritten digit. I want the NN to be able to classify a number from 0 to 9 given an image.

My neural network consists of three layers: the input layer (784 neurons, each one for every pixel of the digit), a hidden layer of 30 neurons (it could also be 100 or 50, but I'm not too worried about hyperparameter tuning yet), and the output layer, 10 neurons, each one representing the activation for every digit. That gives to me two weight matrices: one of 30x724 and a second one of 10x30.

I know and understand the theory behind back propagation, optimization and the mathematical formulas behind that, that's not a problem as such. I can optimize the weights for the second matrix of weights, and the cost is indeed being reduced over time. But I'm not able to keep propagating that back because of the matrix structure.

Knowing that I have find the derivative of the cost w.r.t. the weights:

d(cost) / d(w) = d(cost) / d(f(z)) * d(f(z)) / d(z) * d(z) / d(w)

(Being f the activation function and z the dot product plus the bias of a neuron)

So I'm in the rightmost layer, with an output array of 10 elements. d(cost) / d(f(z)) is the subtraction of the observed an predicted values. I can multiply that by d(f(z)) / d(z), which is just f'(z) of the rightmost layer, also an unidimensional vector of 10 elements, having now d(cost) / d(z) calculated. Then, d(z)/d(w) is just the input to that layer, i.e. the output of the previous one, which is a vector of 30 elements. I figured that I can transpose d(cost) / d(z) so that T( d(cost) / d(z) ) * d(z) / d(w) gives me a matrix of (10, 30), which makes sense because it matches with the dimension of the rightmost weight matrix.

But then I get stuck. The dimension of d(cost) / d(f(z)) is (1, 10), for d(f(z)) / d(z) is (1, 30) and for d(z) / d(w) is (1, 784). I don't know how to come up with a result for this.

This is what I've coded so far. The incomplete part is the _propagate_back method. I'm not caring about the biases yet because I'm just stuck with the weights and first I want to figure this out.

import random
from typing import List, Tuple

import numpy as np
from matplotlib import pyplot as plt

import mnist_loader

np.random.seed(42)

NETWORK_LAYER_SIZES = [784, 30, 10]
LEARNING_RATE = 0.05
BATCH_SIZE = 20
NUMBER_OF_EPOCHS = 5000


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


def sigmoid_der(x):
    return sigmoid(x) * (1 - sigmoid(x))


class Layer:

    def __init__(self, input_size: int, output_size: int):
        self.weights = np.random.uniform(-1, 1, [output_size, input_size])
        self.biases = np.random.uniform(-1, 1, [output_size])
        self.z = np.zeros(output_size)
        self.a = np.zeros(output_size)
        self.dz = np.zeros(output_size)

    def feed_forward(self, input_data: np.ndarray):
        input_data_t = np.atleast_2d(input_data).T
        dot_product = self.weights.dot(input_data_t).T[0]
        self.z = dot_product + self.biases
        self.a = sigmoid(self.z)
        self.dz = sigmoid_der(self.z)


class Network:

    def __init__(self, layer_sizes: List[int], X_train: np.ndarray, y_train: np.ndarray):
        self.layers = [
            Layer(input_size, output_size)
            for input_size, output_size
            in zip(layer_sizes[0:], layer_sizes[1:])
        ]
        self.X_train = X_train
        self.y_train = y_train

    @property
    def predicted(self) -> np.ndarray:
        return self.layers[-1].a

    def _normalize_y(self, y: int) -> np.ndarray:
        output_layer_size = len(self.predicted)
        normalized_y = np.zeros(output_layer_size)
        normalized_y[y] = 1.

        return normalized_y

    def _calculate_cost(self, y_observed: np.ndarray) -> int:
        y_observed = self._normalize_y(y_observed)
        y_predicted = self.layers[-1].a

        squared_difference = (y_predicted - y_observed) ** 2

        return np.sum(squared_difference)

    def _get_training_batches(self, X_train: np.ndarray, y_train: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        train_batch_indexes = random.sample(range(len(X_train)), BATCH_SIZE)

        return X_train[train_batch_indexes], y_train[train_batch_indexes]

    def _feed_forward(self, input_data: np.ndarray):
        for layer in self.layers:
            layer.feed_forward(input_data)
            input_data = layer.a

    def _propagate_back(self, X: np.ndarray, y_observed: int):
        """
        der(cost) / der(weight) = der(cost) / der(predicted) * der(predicted) / der(z) * der(z) / der(weight)
        """
        y_observed = self._normalize_y(y_observed)
        d_cost_d_pred = self.predicted - y_observed

        hidden_layer = self.layers[0]
        output_layer = self.layers[1]

        # Output layer weights
        d_pred_d_z = output_layer.dz
        d_z_d_weight = hidden_layer.a  # Input to the current layer, i.e. the output from the previous one

        d_cost_d_z = d_cost_d_pred * d_pred_d_z
        d_cost_d_weight = np.atleast_2d(d_cost_d_z).T * np.atleast_2d(d_z_d_weight)

        output_layer.weights -= LEARNING_RATE * d_cost_d_weight

        # Hidden layer weights
        d_pred_d_z = hidden_layer.dz
        d_z_d_weight = X

        # ...

    def train(self, X_train: np.ndarray, y_train: np.ndarray):
        X_train_batch, y_train_batch = self._get_training_batches(X_train, y_train)
        cost_over_epoch = []

        for epoch_number in range(NUMBER_OF_EPOCHS):
            X_train_batch, y_train_batch = self._get_training_batches(X_train, y_train)

            cost = 0
            for X_sample, y_observed in zip(X_train_batch, y_train_batch):
                self._feed_forward(X_sample)
                cost += self._calculate_cost(y_observed)
                self._propagate_back(X_sample, y_observed)

            cost_over_epoch.append(cost / BATCH_SIZE)

        plt.plot(cost_over_epoch)
        plt.ylabel('Cost')
        plt.xlabel('Epoch')
        plt.savefig('cost_over_epoch.png')


training_data, validation_data, test_data = mnist_loader.load_data()
X_train, y_train = training_data[0], training_data[1]

network = Network(NETWORK_LAYER_SIZES, training_data[0], training_data[1])
network.train(X_train, y_train)

This is the code for mnist_loader, in case someone wanted to reproduce the example:

import pickle
import gzip


def load_data():
    f = gzip.open('data/mnist.pkl.gz', 'rb')
    training_data, validation_data, test_data = pickle.load(f, encoding='latin-1')
    f.close()

    return training_data, validation_data, test_data

Answer 1

I'm not familiar with python I'm afraid, but I'll give a go at presenting some of the maths...

In order to perform the back propagation using matrices, matrix transpose is used for different sized layers. Also note that when using matrices for this, we need to distinguish different types of matrix multiplication, namely matrix product, and the Hadamard product which operates differently (the latter of which is denoted by a circle with a dot in the centre).

Note the back propagation formulas:

(EQ1) \begin{equation*}\delta ^{l} = (w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l})\end{equation*}

(EQ2) \begin{equation*}\frac{\partial E}{\partial w}=a^{l-1}\delta ^{l}\end{equation*}

As you can see the transpose of the weight matrix is used to multiply against the delta of the layer.

As an example, consider a simple network with an input layer of 3 neurons, and an output layer of 2 neurons.

The unactivated feed forward output is given by...

\begin{equation*}z_{1} = h_{1} w_{1} + h_{2} w_{3} + h_{3} w_{5}\end{equation*} \begin{equation*}z_{2} = h_{1} w_{2} + h_{2} w_{4} + h_{3} w_{6}\end{equation*}

Which can be represented in matrix form as...

\begin{equation*} z = \begin{bmatrix}z_{1}\\z_{2}\end{bmatrix} = \begin{bmatrix}h_{1} w_{1} + h_{2} w_{3} + h_{3} w_{5}\\h_{1} w_{2} + h_{2} w_{4} + h_{3} w_{6}\end{bmatrix} \end{equation*}

(EQ3) \begin{equation*} = \begin{bmatrix}w_{1}&w_{3}&w_{5}\\w_{2}&w_{4}&w_{6}\end{bmatrix} \begin{bmatrix}h_{1}\\h_{2}\\h_{3}\end{bmatrix} \end{equation*}

Which allows us to forward propagate from a layer of 3 neurons, to a layer of 2 neurons (via matrix multiplication).

For the back propagation, we need the error of each neruon...

N.B Cost is a metric used to denote the error of the entire network and depends on your cost function. It is usually the mean of the sum of the errors of each nuerons, but of course depends on your cost function used.

e.g for MSE is... \begin{equation*} C_{MSE}=\frac{1}{N}\sum (o_{n}-t_{n})^{2} \end{equation*}

We are interested in the derivative of the error for each neuron (not the cost), which by the chain rule is...

\begin{equation*} \frac{\partial E}{\partial z} = \frac{\partial E}{\partial o} \frac{\partial o}{\partial z} \end{equation*}

Expressed in matrix form...

\begin{equation*} \frac{\partial E}{\partial z} = \begin{bmatrix} \frac{\partial E}{\partial z_{1}}\\ \frac{\partial E}{\partial z_{2}} \end{bmatrix} = \begin{bmatrix} \frac{\partial E}{\partial o_{1}} \frac{\partial o_{1}}{\partial z_{1}}\\ \frac{\partial E}{\partial o_{2}} \frac{\partial o_{2}}{\partial z_{2}} \end{bmatrix} \end{equation*}

(EQ4) \begin{equation*} = \begin{bmatrix} \frac{\partial E}{\partial o_{1}} \\ \frac{\partial E}{\partial o_{2}} \end{bmatrix} \odot \begin{bmatrix} \frac{\partial o_{1}}{\partial z_{1}} \\ \frac{\partial o_{2}}{\partial z_{2}} \end{bmatrix} \end{equation*}

Note the use of the hadamard product here. In fact, given these are simply vectors, vector dot product would work, but hadamard is used as this becomes important later when using this expression with the matrix equations as we want to distinguish the use of hadamard product as opposed to matrix product.

We start off with our first delta error, which for the first layer back propagation is...

\begin{equation*} \delta ^{L} = \frac{\partial E}{\partial o} \end{equation*}

And then we want to calculate the next delta error using the formula (EQ1)...

\begin{equation*} \delta ^{l} = (w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l}) = \frac{\partial E}{\partial h} \end{equation*}

Explicitly, the equations are...

\begin{equation*}\frac{\partial E}{\partial h_{1}}=\frac{\partial E}{\partial z_{1}} w_{1} + \frac{\partial E}{\partial z_{2}} w_{2}\end{equation*} \begin{equation*}\frac{\partial E}{\partial h_{2}}=\frac{\partial E}{\partial z_{1}} w_{3} + \frac{\partial E}{\partial z_{2}} w_{4}\end{equation*} \begin{equation*}\frac{\partial E}{\partial h_{3}}=\frac{\partial E}{\partial z_{1}} w_{5} + \frac{\partial E}{\partial z_{2}} w_{6}\end{equation*}

Also given the transpose of the weight matrix which allows us to back propagate from a layer of 2 neurons to a layer of 3 neurons (via matrix multiplication) ...

\begin{equation*} \left (w \right )^{T} = \left ( \begin{bmatrix} w_{1} & w_{3} & w_{5}\\ w_{2} & w_{4} & w_{6} \end{bmatrix} \right )^{T} = \begin{bmatrix} w_{1} & w_{2}\\ w_{3} & w_{4}\\ w_{5} & w_{6} \end{bmatrix} \end{equation*}

So in a similar way to how we represented the forward pass (EQ3), this can be represented in matrix form...

\begin{equation*} \frac{\partial E}{\partial h} = \begin{bmatrix} \frac{\partial E}{\partial h_{1}}\\ \frac{\partial E}{\partial h_{2}}\\ \frac{\partial E}{\partial h_{3}} \end{bmatrix} = \begin{bmatrix} \frac{\partial E}{\partial z_{1}} w_{1} + \frac{\partial E}{\partial z_{2}} w_{2} \\ \frac{\partial E}{\partial z_{1}} w_{3} + \frac{\partial E}{\partial z_{2}} w_{4} \\ \frac{\partial E}{\partial z_{1}} w_{5} + \frac{\partial E}{\partial z_{2}} w_{6} \end{bmatrix} \end{equation*}

(EQ5) \begin{equation*} = \begin{bmatrix} w_{1} & w_{2}\\ w_{3} & w_{4}\\ w_{5} & w_{6} \end{bmatrix} \begin{bmatrix} \frac{\partial E}{\partial z_{1}} \\ \frac{\partial E}{\partial z_{2}} \end{bmatrix} \end{equation*}

And then plugging the hadamard version of the delta (EQ4) into this, we get...

\begin{equation*} \begin{bmatrix} w_{1} & w_{2}\\ w_{3} & w_{4}\\ w_{5} & w_{6} \end{bmatrix} \begin{bmatrix} \frac{\partial E}{\partial o_{1}} \\ \frac{\partial E}{\partial o_{2}} \end{bmatrix} \odot \begin{bmatrix} \frac{\partial o_{1}}{\partial z_{1}} \\ \frac{\partial o_{2}}{\partial z_{2}} \end{bmatrix} \end{equation*}

aka (EQ1) ...

\begin{equation*}(w^{l+1})^{T} \delta ^{l+1} \odot \sigma {}' (z ^{l})\end{equation*}

And thus we have back propagated from a layer of 2 neurons, to a layer of 3 neurons (via matrix multiplication) thanks to transpose.

For completeness... the other aspect of the back propagation to uses matrices, is the delta weight matrix....

\begin{equation*} \frac{\partial E}{\partial w} = \begin{bmatrix} \frac{\partial E}{\partial w_{1}} & \frac{\partial E}{\partial w_{3}} & \frac{\partial E}{\partial w_{5}}\\ \frac{\partial E}{\partial w_{2}} & \frac{\partial E}{\partial w_{4}} & \frac{\partial E}{\partial w_{6}} \end{bmatrix} \end{equation*}

As mentioned before, you need to cache weight matrix and the activated layer output of a forward pass of the network.

In a similar vein to (EQ3)...

We explcitily have the equations...

\begin{equation*}\frac{\partial E}{\partial w_{1}}=\frac{\partial E}{\partial z_{1}}h_{1}\end{equation*} \begin{equation*}\frac{\partial E}{\partial w_{2}}=\frac{\partial E}{\partial z_{2}}h_{1}\end{equation*} \begin{equation*}\frac{\partial E}{\partial w_{3}}=\frac{\partial E}{\partial z_{1}}h_{2}\end{equation*} \begin{equation*}\frac{\partial E}{\partial w_{4}}=\frac{\partial E}{\partial z_{2}}h_{2}\end{equation*} \begin{equation*}\frac{\partial E}{\partial w_{5}}=\frac{\partial E}{\partial z_{1}}h_{3}\end{equation*} \begin{equation*}\frac{\partial E}{\partial w_{6}}=\frac{\partial E}{\partial z_{2}}h_{3}\end{equation*}

Note the use of h1, h2 and h3, which are the activated outputs of the previous layer.. (or in the case of our example these are the inputs).

Which we represent in matrix form...

\begin{equation*} \frac{\partial E}{\partial w} = \begin{bmatrix} \frac{\partial E}{\partial z_{1}}h_{1} & \frac{\partial E}{\partial z_{1}}h_{2} & \frac{\partial E}{\partial z_{1}}h_{3}\\ \frac{\partial E}{\partial z_{2}}h_{1} & \frac{\partial E}{\partial z_{2}}h_{2} & \frac{\partial E}{\partial z_{2}}h_{3} \end{bmatrix} \end{equation*} (EQ6) \begin{equation*} = \begin{bmatrix} \frac{\partial E}{\partial z_{1}} \\ \frac{\partial E}{\partial z_{2}} \end{bmatrix} \begin{bmatrix} h_{1} & h_{2} & h_{3} \end{bmatrix} \end{equation*}

Which just so happens to be (EQ2) ...

\begin{equation*}\frac{\partial E}{\partial w}=\delta ^{l}a^{l-1}\end{equation*}

:)

Since the delta weight matrix and the original weight matrix are the same dimensions, it is trivia to apply the learning rate...

\begin{equation*} w = \begin{bmatrix} w_{1} - \alpha \frac{\partial E}{\partial w_{1}} & w_{3} - \alpha \frac{\partial E}{\partial w_{3}} & w_{5} - \alpha \frac{\partial E}{\partial w_{5}}\\ w_{2} - \alpha \frac{\partial E}{\partial w_{2}} & w_{4} - \alpha \frac{\partial E}{\partial w_{4}} & w_{6} - \alpha \frac{\partial E}{\partial w_{6}} \end{bmatrix} \end{equation*}

\begin{equation*} = \begin{bmatrix}w_{1}&w_{3}&w_{5}\\w_{2}&w_{4}&w_{6}\end{bmatrix} - \alpha \begin{bmatrix} \frac{\partial E}{\partial w_{1}} & \frac{\partial E}{\partial w_{3}} & \frac{\partial E}{\partial w_{5}}\\ \frac{\partial E}{\partial w_{2}} & \frac{\partial E}{\partial w_{4}} & \frac{\partial E}{\partial w_{6}} \end{bmatrix} \end{equation*}

\begin{equation*} = w - \alpha \frac{\partial E}{\partial w}\end{equation*}

I'll omit the equations for the bias, as these are easy to put in and don't require any other equations other than the ones presented here. So in summary, the equations you want are (EQ3), (EQ5) and (EQ6) to use matrices, and the ability to change move between layers of differing number of neurons is the matrix transpose. Hope this helps, and let me know if you want me to expand on anything.

It might be also important to note that you appear to be using MSE as the cost. This is perhaps not optimal for the MNIST data set which is a clasification. A better cost function to use would be cross entropy.

Happy propagating!