# Yet another introduction to Neural Networks

BlogDeep LearningIntro|Neural Networksposted by Burak Himmetoglu August 31, 2017 Burak Himmetoglu

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: classificationdef 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 probabilitiesreturn 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,

## Burak Himmetoglu

As an aspiring data scientist, I analyze large amounts of data, search for patterns, and solve interesting problems spanning a wide range of areas. I also work on applications of machine learning methods to predict electronic properties of molecules for discovering new compounds computationally. Currently, I am a staff member at UCSB as a computational physicist and High Performance Computing (HPC) specialist. I support to faculty and researchers in obtaining high performance computing resources, provide computational research consultation and manage supercomoputing allocations of UCSB. I help researchers port their codes, usually from Python, R and Matlab to C/C++ for high performance computing environments, I am enthusiastic about applying my skills to solve difficult business problems and develop new products using data science methods.

## 6 Reasons Your Business is Falling Behind the AI Race

BusinessAccelerate AIposted by Elizabeth Wallace Mar 22, 2019

## Watch: Unsupervised Feature Learning with Matrix Decomposition

ConferencesMachine LearningMachine Learningposted by ODSC Team Mar 22, 2019

## What Will the Next Generation of Data Scientists Look Like?

OpinionCareer Insightsposted by Daniel Gutierrez, ODSC Mar 21, 2019

- Learn Interpretability for Data Science 123 views | by Rajiv Shah | under Conferences, Featured Post
- Reviewing Amazon’s Machine Learning University – Is it Worth All of the Hype? 99 views | by Daniel Gutierrez, ODSC | under Featured Post, Modeling
- Introduction to R Shiny 61 views | by Alyssa Columbus | under Data Visualization, Featured Post, R, Tools & Languages