Yet another introduction to Neural Networks

In this notebook, I will explain how to implement a neural network from scratch and use the version of MNIST dataset that is provided within Scikit-Learn for testing. I will specificallty illustrate the use of Python classes to define layers in the network as objects. Each layer object has forward and backward propagation methods which leads to compact, easily readable code. In writing this tutorial, I’ve had inspiration from Peter Roelants’ page.

## Imports

import numpy as np
from sklearn.datasets import load_digits

Data preparation

After loading the data, we divide it into three parts, training, validation and testing sets. The validation set is to be used to determine the hyperparameters (i.e. number and size of hidden layers and regularization parameter) and the testing set is a separate piece of data to assess the final model performance.

## Load the data and reshape images
digits = load_digits()
n_samples = len(digits['images'])
data = digits['images'].reshape((n_samples, -1)); targets = digits['target']

## Train-test splitting
from sklearn.model_selection import train_test_split
X, X_test, y, y_test = train_test_split(data, targets, test_size=0.33, random_state=111)

## Train and validation splitting
X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.40, random_state=123)

Activation function and classifier

We will use the sigmoid activation function (sigma(x) = 1/(1+e^{-x})) in the layers of the network. For the classifier, we will use the softmax functio, which results in class probabilities:

({rm softmax}_i({bf Y}) = frac{e^{Y_i}}{sum_{k=1}^K, e^{Y_k}}) where (Y_i) is vector from the output of the neural network for a given example, and (K) is the number of classes (10 in our case). Both functions are implemented below:

## Define the activation function, the Softmax classifier 
# Sigmoid function
def sigmoid(X):
    return 1.0 / (1.0 + np.exp(-X))

# Softmax function
def softmax(X):
    temp = np.max(X, axis = 0) # Determine the maximum of each column
    X = X - temp               # Subtract the max from each: does not change outcome
    return np.exp(X) / np.sum(np.exp(X), axis = 0)

In the above implementations, the assumption is that the arguments of both functions are cast in a matrix columnwise, so that each column represents one example: ({bf X} in {mathbb R}^{n_{rm feat} times N}) where (N) is the number of examples and (n_{rm feat}) is the number of features (i.e. number of pixels in our case). In softmax, we have subtracted the maximum of each column from each training example vector, an operation that does not change the outcome, for numerical stability.

Now that we have our data, the activation function and the classifier, we can construct the layers of the network.

Linear Update

First, we define the class LinearUpdate which performs the linear transformation of the (derived) features in the current layer:

({bf Y}^{(i)} = {bf W}^{T} cdot A^{(i)} + {bf b}) In this equation, ({bf A}^{(i)} in {mathbb R}^{n_{rm in}}) represents the current layer state (i.e features in case of input layer and derived features in case of hidden layers) with (n_{rm in}) being the number of neurons. (Y^{(i)} in {mathbb R}^{n_{in out}}) is the linear ouput of the current layer (which will later be the argument of an activation function) which is passed as the input to the next layer with (n_{rm out}) neurons. The supersript (i) refers to the training example being considered. Instead of using a for loop over the training examples, we can cast them in a matrix where each column is one training example vector (Y^{(i)}), i.e.

(Y^{(i)}_j rightarrow Y_{ji} : {bf Y} in {mathbb R}^{n_{rm out} times N}) where (N) is the number of training examples. Similary (A^{(i)}_j rightarrow A_{ji} : {bf A} in {mathbb R}^{n_{rm in} times N}).

Forward propagation

The above equation can be simply written in matrix notation as ({mathbf Y} = {mathbf W}^T cdot {bf A} + {bf b}). The weights ({bf W} in {mathbb R}^{n_{rm in} times n_{rm out}}) and biases ({bf b} in {mathbb R}^{n_{rm in}}) will be determined during training by the minimization of a loss function (L).

The LinearUpdate object is initialized with (n_{rm in}), (n_{rm out}), the weights and biases. The forward method below implements the above equation.

Back propagation

The backpropagation method simply takes the gradient of the loss function with respect to the state of the next layer ((nabla_{bf Y} L)) and computes the gradients with respect to the current state ((nabla_{bf A}L)), weights ((nabla_{bf W}L)) and biases ((nabla_{bf b}L)). The backpropagation rules can be derived by noting that:

(Y_{ji} = sum_{k=1}^{n_{rm in}}, left( W^T right)_{jk}, A_{ki} + b_j) First, the gradient of the loss function (L) with respect to the weights ({bf W}) can be written as:

(frac{partial L}{partial W_{ln}} = sum_{i=1}^N, sum_{j=1}^{n_{rm out}}, frac{partial L}{partial Y_{ji}}, frac{partial Y_{ji}}{partial W_{ln}} = sum_{i=1}^N, frac{partial L}{partial Y_{ni}}, A_{li}) We denote the gradient from the next layer (frac{partial L}{partial Y_{ji}}) as (nabla_{bf Y} L), and the above equation in matrix form becomes

(nabla_{bf W} L = {bf A} cdot (nabla_{bf Y} L)^T) Similarly, we can compute gradients with respect to ({bf A}) and ({bf b}), and following similar steps, we obtain

(( nabla_{bf b} L )_l = sum_{i=1}^N, (nabla_{bf Y} L)_{li}

nabla_{bf A} L = {bf W} cdot (nabla_{bf Y} L)) All of these three backpropagation rules are implemented below in the backward method of LinearUpdate.

## Define Linear Update
class LinearUpdate(object):
    def __init__ (self, n_in, n_out, W = None, b = None):
        # Initialize W randomly
        eps = np.sqrt(6.0) / np.sqrt(n_in + n_out)
        self.W = (np.random.uniform(-eps, eps, (n_in, n_out)) if W is None else W)
        
        # Initialize biases as zero
        self.b = np.zeros((n_out,1))
        
    def forward(self, A):
        "" Forward propagation method: A is current state, method results in next state linearly
            Y <- W.T * A + b
        ""
        return np.dot(self.W.T, A) + self.b
    
    def backward(self, A, gY):
        "" Backward propagation method: A is current state, gY is backpropagated derivative from 
            next layer.
        ""
        gW = np.dot(A, gY.T)       # dL/dW_ab = sum_i A_ai * (gY)_bi
        gA = np.dot(self.W, gY)    # dL/dA_ab = sum_j ( (W)_aj * gY_jb )  
        gb = np.sum(gY, axis=1)    # dL/db_l = sum_i (gY)_li
        
        return gA, gW, gb

Logistic Update

Now we implement the object which takes the outcome of the linear update, and transforms it with the actiovation function. Our choice for the activation is the sigmoid function defined above. The activation function takes the input into the layer (let’s denote by (Z) and generates an output (sigma(Z)) which will be passed to the next layer. Specifically, we will have:

({bf Z} = {bf W}^T cdot {bf A} + {bf b}

{bf Y} = sigma({bf Z})) i.e., the linear output ({bf Z}) will be transformed by the activation function and then passed to the next layer.

The forward and backward propagation methods are easily obtained as follows:

({bf Y} = sigma(Z), quad
frac{partial L}{partial {bf Z}} = frac{partial L}{partial {bf Y}} cdot frac{partial {bf Y}}{partial {bf Z}} = {bf Y} odot (1-{bf Y}) odot (nabla_{bf Y}L)) where ({bf Y}) is the input into the next layer and (nabla_{bf Y}L) is the gradient with respect to the input of the next layer and (odot) denote elementwise multiplication.

## Define Logistic Update
class LogisticUpdate(object):
    def forward(self, Z):
        "" Sigmoid activation: ""
        return sigmoid(Z) #  Y = sigmoid(Z)
    
    def backward(self, Y, grad_out):
        "" Backward propagation: ""
        return np.multiply(Y * (1 - Y), grad_out)   # dL/dZ = dL/dY * Y * (1-Y)

Classifier

The output layer will be the classifier which deserves its own class. The forward method is still the sigmoid function, which will output class probabilities. The backward method implements the gradient of the loss function with respect to the outputs of the network. Finally, the crossEntropy method defines the loss function:

(L = -frac{1}{N}, sum_{i=1}^{N}, sum_{k=1}^K, log(P_{ki}), I(T_i = k)) where (P_{ki}) is the calculated probability of training example ((i)) being of class (k), and the function (I) is (1) when the target (i.e. the actual value of the digit) is of class (k), and zero otherwise. Evaluating the derivative of the loss function with respect to the ouput layer state ({bf Y}) results in

(frac{partial L}{partial {bf Y}} = frac{1}{N}, ({bf P} – {bf I})) These methods are implemented below.

# Define Softmax Classifier layer
class SoftmaxClassifier(object):
    def forward (self, A): 
        "" Given state A, produces output probs P by softmax ""
        return softmax(A)
    
    def backward(self, P, T):
        "" Given output probs P and targets T, produces output gradient ""
        expansionMatrix = np.eye(P.shape[0], dtype = int)
        expandedTarget = expansionMatrix[:, T]
        
        return P - expandedTarget # No division by number of samples yet.
    
    def crossEntropy(self, P, T):
        "" Computes cross entropy ""
        expansionMatrix = np.eye(P.shape[0], dtype = int)
        expandedTarget = expansionMatrix[:, T]
        
        CE = -np.sum(expandedTarget * np.log(P + 1e-30))/P.shape[1]
        return CE

Neurons

Now that we have defined how states in a layer is forward propagated (first linearly, then by the activation function), we can use the LinearUpdate and LogisticUpdate classes to define the Layer class. The layer class first linearly transforms the current state vectors ({bf A}) and then feeds them into the activation layer to yield the input to the next layer ({bf Y}) by the forward method. In the backward method, the incoming gradient is first backpropagated through the logistic update, then by the linear update to yield the gradients with respect to curent layer states, weights and biases.

Just for a sanity check, we can explicitly evaluate the backpropagation for ({bf A}) as an example as follows:

which are is implemented in two steps; first ({bf Y} odot (1-{bf Y}) odot (nabla_{bf Y}L)) as the backward method of the LogisticUpdate which is fed into the backward method of LinearUpdate to yield its dot product with ({bf W}).

## Define Layes (combine linear and logistic updates)
class Layer(object):
    def __init__ (self, n_in, n_out):
        self.linear = LinearUpdate(n_in, n_out)
        self.logistic = LogisticUpdate()
        
    def forward(self, A):
        "" Forward propagation method ""
        return self.logistic.forward(self.linear.forward(A))
    
    def backward(self, A, Y, grad_out):
        "" Backward propagation method ""
        # First the derivative of the logistic State
        gZ = self.logistic.backward(Y, grad_out)
        # Then, gZ becomes gY for the linear state
        gA, gW, gb = self.linear.backward(A, gZ) 
        
        return gA, gW, gb

Construct the network

Now we have the Layers and the Classifier objects, we can construct the neural network. We also add a regularization term to the loss function, which will be determined by the validation set performance below. The regularization term is given by

(L_{rm reg} = frac{lambda}{2 N}, sum_{l=1}^{N_{rm layers}}, vert {bf W}_l vert^2) which is a sum over all the layers and penalizes the weights in all the layers using their square norm defined as

(vert {bf W}_l vert^2 = sum_{i,j}, (W_l)_{ij}^2)

## Construct a two hidden layer network

class nnet(object):

    def __init__ (self, nInput, nHidden1, nHidden2, K):

        "" Initiate method for the net ""

        self.inputLayer = Layer(nInput, nHidden1)     # Input layer

        self.hiddenLayer1 = Layer(nHidden1, nHidden2) # 1st hidden layers

        self.hiddenLayer2 = Layer(nHidden2, K)        # 2nd hidden layer

        self.classifier = SoftmaxClassifier()         # Output: classification
    def forward(self, input_train):

        "" Perform forward propagation through all layers ""

        A1 = input_train.T                 # Initial data

        A2 = self.inputLayer.forward(A1)   # inp -> hid1 

        A3 = self.hiddenLayer1.forward(A2) # hid1 -> hid2

        A4 = self.hiddenLayer2.forward(A3) # hid2 -> out

        P = self.classifier.forward(A4)    # output probabilities
        return A1, A2, A3, A4, P
    def backward(self, P, A4, A3, A2, A1, T, lam = 0.0):

        "" Back propagation method through all layers ""

        grad_out = self.classifier.backward(P, T)                    # output grads

        gA3, gW3, gb3 = self.hiddenLayer2.backward(A3, A4, grad_out) # hid2 grads

        gA2, gW2,