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