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)
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 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 output_layer = self.layers # 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, training_data network = Network(NETWORK_LAYER_SIZES, training_data, training_data) 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