Deep Learning from first principles in Python, R and Octave – Part 6

“Today you are You, that is truer than true. There is no one alive who is Youer than You.”
Dr. Seuss

“Explanations exist; they have existed for all time; there is always a well-known solution to every human problem — neat, plausible, and wrong.”
H L Mencken

Introduction

In this 6th instalment of ‘Deep Learning from first principles in Python, R and Octave-Part6’, I look at a couple of different initialization techniques used in Deep Learning, L2 regularization and the ‘dropout’ method. Specifically, I implement “He initialization” & “Xavier Initialization”. My earlier posts in this series of Deep Learning included

1. Part 1 – In the 1st part, I implemented logistic regression as a simple 2 layer Neural Network
2. Part 2 – In part 2, implemented the most basic of Neural Networks, with just 1 hidden layer, and any number of activation units in that hidden layer. The implementation was in vectorized Python, R and Octave
3. Part 3 -In part 3, I derive the equations and also implement a L-Layer Deep Learning network with either the relu, tanh or sigmoid activation function in Python, R and Octave. The output activation unit was a sigmoid function for logistic classification
4. Part 4 – This part looks at multi-class classification, and I derive the Jacobian of a Softmax function and implement a simple problem to perform multi-class classification.
5. Part 5 – In the 5th part, I extend the L-Layer Deep Learning network implemented in Part 3, to include the Softmax classification. I also use this L-layer implementation to classify MNIST handwritten digits with Python, R and Octave.

The code in Python, R and Octave are identical, and just take into account some of the minor idiosyncrasies of the individual language. In this post, I implement different initialization techniques (random, He, Xavier), L2 regularization and finally dropout. Hence my generic L-Layer Deep Learning network includes these additional enhancements for enabling/disabling initialization methods, regularization or dropout in the algorithm. It already included sigmoid & softmax output activation for binary and multi-class classification, besides allowing relu, tanh and sigmoid activation for hidden units.

This R Markdown file and the code for Python, R and Octave can be cloned/downloaded from Github at DeepLearning-Part6

1. Initialization techniques

The usual initialization technique is to generate Gaussian or uniform random numbers and multiply it by a small value like 0.01. Two techniques which are used to speed up convergence is the He initialization or Xavier. These initialization techniques enable gradient descent to converge faster.

1.1 a Default initialization – Python

This technique just initializes the weights to small random values based on Gaussian or uniform distribution

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn.linear_model
import pandas as pd
import sklearn
import sklearn.datasets
train_X, train_Y, test_X, test_Y = load_dataset()
# Set the layers dimensions
layersDimensions = [2,7,1]

# Train a deep learning network with random initialization
parameters = L_Layer_DeepModel(train_X, train_Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid",learningRate = 0.6, num_iterations = 9000, initType="default", print_cost = True,figure="fig1.png")

# Clear the plot
plt.clf()
plt.close()

# Plot the decision boundary
plot_decision_boundary(lambda x: predict(parameters, x.T), train_X, train_Y,str(0.6),figure1="fig2.png")

1.1 b He initialization – Python

‘He’ initialization attributed to He et al, multiplies the random weights by
$\sqrt{\frac{2}{dimension\ of\ previous\ layer}}$

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn.linear_model
import pandas as pd
import sklearn
import sklearn.datasets

train_X, train_Y, test_X, test_Y = load_dataset()
# Set the layers dimensions
layersDimensions = [2,7,1]

# Train a deep learning network with He  initialization
parameters = L_Layer_DeepModel(train_X, train_Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate =0.6,    num_iterations = 10000,initType="He",print_cost = True,                           figure="fig3.png")

plt.clf()
plt.close()
# Plot the decision boundary
plot_decision_boundary(lambda x: predict(parameters, x.T), train_X, train_Y,str(0.6),figure1="fig4.png")

1.1 c Xavier initialization – Python

Xavier  initialization multiply the random weights by
$\sqrt{\frac{1}{dimension\ of\ previous\ layer}}$

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn.linear_model
import pandas as pd
import sklearn
import sklearn.datasets

train_X, train_Y, test_X, test_Y = load_dataset()
# Set the layers dimensions
layersDimensions = [2,7,1]

# Train a L layer Deep Learning network
parameters = L_Layer_DeepModel(train_X, train_Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid",
learningRate = 0.6,num_iterations = 10000, initType="Xavier",print_cost = True,
figure="fig5.png")

# Plot the decision boundary
plot_decision_boundary(lambda x: predict(parameters, x.T), train_X, train_Y,str(0.6),figure1="fig6.png")

1.2a Default initialization – R

source("DLfunctions61.R")
x <- z[,1:2]
y <- z[,3]
X <- t(x)
Y <- t(y)
#Set the layer dimensions
layersDimensions = c(2,11,1)
# Train a deep learning network
retvals = L_Layer_DeepModel(X, Y, layersDimensions,
hiddenActivationFunc='relu',
outputActivationFunc="sigmoid",
learningRate = 0.5,
numIterations = 9000,
initType="default",
print_cost = True)
#Plot the cost vs iterations
iterations <- seq(0,9000,1000)
costs=retvals$costs df=data.frame(iterations,costs) ggplot(df,aes(x=iterations,y=costs)) + geom_point() + geom_line(color="blue") + ggtitle("Costs vs iterations") + xlab("No of iterations") + ylab("Cost") # Plot the decision boundary plotDecisionBoundary(z,retvals,hiddenActivationFunc="relu",lr=0.5) 1.2b He initialization – R source("DLfunctions61.R") # Load the data z <- as.matrix(read.csv("circles.csv",header=FALSE)) x <- z[,1:2] y <- z[,3] X <- t(x) Y <- t(y) # Set the layer dimensions layersDimensions = c(2,11,1) # Train a deep learning network retvals = L_Layer_DeepModel(X, Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.5, numIterations = 9000, initType="He", print_cost = True) #Plot the cost vs iterations iterations <- seq(0,9000,1000) costs=retvals$costs
df=data.frame(iterations,costs)
ggplot(df,aes(x=iterations,y=costs)) + geom_point() + geom_line(color="blue") +
ggtitle("Costs vs iterations") + xlab("No of iterations") + ylab("Cost")

# Plot the decision boundary
plotDecisionBoundary(z,retvals,hiddenActivationFunc="relu",0.5,lr=0.5)

1.2c Xavier initialization – R

## Xav initialization
# Set the layer dimensions
layersDimensions = c(2,11,1)
# Train a deep learning network
retvals = L_Layer_DeepModel(X, Y, layersDimensions,
hiddenActivationFunc='relu',
outputActivationFunc="sigmoid",
learningRate = 0.5,
numIterations = 9000,
initType="Xav",
print_cost = True)
#Plot the cost vs iterations
iterations <- seq(0,9000,1000)
costs=retvals$costs df=data.frame(iterations,costs) ggplot(df,aes(x=iterations,y=costs)) + geom_point() + geom_line(color="blue") + ggtitle("Costs vs iterations") + xlab("No of iterations") + ylab("Cost") # Plot the decision boundary plotDecisionBoundary(z,retvals,hiddenActivationFunc="relu",0.5) 1.3a Default initialization – Octave source("DL61functions.m") # Read the data data=csvread("circles.csv"); X=data(:,1:2); Y=data(:,3); # Set the layer dimensions layersDimensions = [2 11 1]; #tanh=-0.5(ok), #relu=0.1 best! # Train a deep learning network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.5, lambd=0, keep_prob=1, numIterations = 10000, initType="default"); # Plot cost vs iterations plotCostVsIterations(10000,costs) #Plot decision boundary plotDecisionBoundary(data,weights, biases,keep_prob=1, hiddenActivationFunc="relu")  1.3b He initialization – Octave source("DL61functions.m") #Load data data=csvread("circles.csv"); X=data(:,1:2); Y=data(:,3); # Set the layer dimensions layersDimensions = [2 11 1]; #tanh=-0.5(ok), #relu=0.1 best! # Train a deep learning network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.5, lambd=0, keep_prob=1, numIterations = 8000, initType="He"); plotCostVsIterations(8000,costs) #Plot decision boundary plotDecisionBoundary(data,weights, biases,keep_prob=1,hiddenActivationFunc="relu")  1.3c Xavier initialization – Octave source("DL61functions.m") #Load data data=csvread("circles.csv"); X=data(:,1:2); Y=data(:,3); #Set layer dimensions layersDimensions = [2 11 1]; #tanh=-0.5(ok), #relu=0.1 best! # Train a deep learning network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.5, lambd=0, keep_prob=1, numIterations = 8000, initType="Xav"); plotCostVsIterations(8000,costs) plotDecisionBoundary(data,weights, biases,keep_prob=1,hiddenActivationFunc="relu")  2.1a Regularization : Circles data – Python The cross entropy cost for Logistic classification is given as $J = \frac{1}{m}\sum_{i=1}^{m}y^{i}log((a^{L})^{(i)}) - (1-y^{i})log((a^{L})^{(i)})$ The regularized L2 cost is given by $J = \frac{1}{m}\sum_{i=1}^{m}y^{i}log((a^{L})^{(i)}) - (1-y^{i})log((a^{L})^{(i)}) + \frac{\lambda}{2m}\sum \sum \sum W_{kj}^{l}$ Incidentally, a detailed implementation of L1 and L2 regularization in Machine Learning, besides other feature selection methods can be found in my book ‘Practical Machine Learning with R and Python: Machine Learning in stereo‘ available in Amazon. My book is compact and can be used as a quick reference for ML algorithms and associated metrics import numpy as np import matplotlib import matplotlib.pyplot as plt import sklearn.linear_model import pandas as pd import sklearn import sklearn.datasets exec(open("DLfunctions61.py").read()) #Load the data train_X, train_Y, test_X, test_Y = load_dataset() # Set the layers dimensions layersDimensions = [2,7,1] # Train a deep learning network parameters = L_Layer_DeepModel(train_X, train_Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid",learningRate = 0.6, lambd=0.1, num_iterations = 9000, initType="default", print_cost = True,figure="fig7.png") # Clear the plot plt.clf() plt.close() # Plot the decision boundary plot_decision_boundary(lambda x: predict(parameters, x.T), train_X, train_Y,str(0.6),figure1="fig8.png") plt.clf() plt.close() #Plot the decision boundary plot_decision_boundary(lambda x: predict(parameters, x.T,keep_prob=0.9), train_X, train_Y,str(2.2),"fig8.png",) 2.1 b Regularization: Spiral data – Python import numpy as np import matplotlib import matplotlib.pyplot as plt import sklearn.linear_model import pandas as pd import sklearn import sklearn.datasets exec(open("DLfunctions61.py").read()) N = 100 # number of points per class D = 2 # dimensionality K = 3 # number of classes X = np.zeros((N*K,D)) # data matrix (each row = single example) y = np.zeros(N*K, dtype='uint8') # class labels for j in range(K): ix = range(N*j,N*(j+1)) r = np.linspace(0.0,1,N) # radius t = np.linspace(j*4,(j+1)*4,N) + np.random.randn(N)*0.2 # theta X[ix] = np.c_[r*np.sin(t), r*np.cos(t)] y[ix] = j # Plot the data plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap=plt.cm.Spectral) plt.clf() plt.close() #Set layer dimensions layersDimensions = [2,100,3] y1=y.reshape(-1,1).T # Train a deep learning network parameters = L_Layer_DeepModel(X.T, y1, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 0.6,lambd=0.2, num_iterations = 5000, print_cost = True,figure="fig9.png") plt.clf() plt.close() W1=parameters['W1'] b1=parameters['b1'] W2=parameters['W2'] b2=parameters['b2'] plot_decision_boundary1(X, y1,W1,b1,W2,b2,figure2="fig10.png") 2.2a Regularization: Circles data – R source("DLfunctions61.R") #Load data df=read.csv("circles.csv",header=FALSE) z <- as.matrix(read.csv("circles.csv",header=FALSE)) x <- z[,1:2] y <- z[,3] X <- t(x) Y <- t(y) #Set layer dimensions layersDimensions = c(2,11,1) # Train a deep learning network retvals = L_Layer_DeepModel(X, Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.5, lambd=0.1, numIterations = 9000, initType="default", print_cost = True)  #Plot the cost vs iterations iterations <- seq(0,9000,1000) costs=retvals$costs
df=data.frame(iterations,costs)
ggplot(df,aes(x=iterations,y=costs)) + geom_point() + geom_line(color="blue") +
ggtitle("Costs vs iterations") + xlab("No of iterations") + ylab("Cost")

# Plot the decision boundary
plotDecisionBoundary(z,retvals,hiddenActivationFunc="relu",0.5)

2.2b Regularization:Spiral data – R

# Read the spiral dataset
source("DLfunctions61.R")

# Setup the data
X <- Z[,1:2]
y <- Z[,3]
X <- t(X)
Y <- t(y)

layersDimensions = c(2, 15, 6, 3)
# Train a deep learning network
retvals = L_Layer_DeepModel(X, Y, layersDimensions,
hiddenActivationFunc='relu',
outputActivationFunc="softmax",
learningRate = 5.1,
lambd=0.01,
numIterations = 9000,
print_cost = True)
parameters<-retvals$parameters plotDecisionBoundary1(Z,parameters) 2.3a Regularization: Circles data – Octave source("DL61functions.m") #Load data data=csvread("circles.csv"); X=data(:,1:2); Y=data(:,3); layersDimensions = [2 11 1]; #tanh=-0.5(ok), #relu=0.1 best! # Train a deep learning network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.5, lambd=0.2, keep_prob=1, numIterations = 8000, initType="default"); plotCostVsIterations(8000,costs) #Plot decision boundary plotDecisionBoundary(data,weights, biases,keep_prob=1,hiddenActivationFunc="relu")  2.3b Regularization:Spiral data 2 – Octave source("DL61functions.m") data=csvread("spiral.csv"); # Setup the data X=data(:,1:2); Y=data(:,3); layersDimensions = [2 100 3] # Train a deep learning network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 0.6, lambd=0.2, keep_prob=1, numIterations = 10000); plotCostVsIterations(10000,costs) #Plot decision boundary plotDecisionBoundary1(data,weights, biases,keep_prob=1,hiddenActivationFunc="relu")  3.1 a Dropout: Circles data – Python The ‘dropout’ regularization technique was used with great effectiveness, to prevent overfitting by Alex Krizhevsky, Ilya Sutskever and Prof Geoffrey E. Hinton in the Imagenet classification with Deep Convolutional Neural Networks The technique of dropout works by dropping a random set of activation units in each hidden layer, based on a ‘keep_prob’ criteria in the forward propagation cycle. Here is the code for Octave. A ‘dropoutMat’ is created for each layer which specifies which units to drop Note: The same ‘dropoutMat has to be used which computing the gradients in the backward propagation cycle. Hence the dropout matrices are stored in a cell array.  for l =1:L-1 ... D=rand(size(A)(1),size(A)(2)); D = (D < keep_prob) ; # Zero out some hidden units A= A .* D; # Divide by keep_prob to keep the expected value of A the same A = A ./ keep_prob; # Store D in a dropoutMat cell array dropoutMat{l}=D; ... endfor In the backward propagation cycle we have  for l =(L-1):-1:1 ... D = dropoutMat{l}; # Zero out the dAl based on same dropout matrix dAl= dAl .* D; # Divide by keep_prob to maintain the expected value dAl = dAl ./ keep_prob; ... endfor  import numpy as np import matplotlib import matplotlib.pyplot as plt import sklearn.linear_model import pandas as pd import sklearn import sklearn.datasets exec(open("DLfunctions61.py").read()) #Load the data train_X, train_Y, test_X, test_Y = load_dataset() # Set the layers dimensions layersDimensions = [2,7,1] # Train a deep learning network parameters = L_Layer_DeepModel(train_X, train_Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid",learningRate = 0.6, keep_prob=0.7, num_iterations = 9000, initType="default", print_cost = True,figure="fig11.png") # Clear the plot plt.clf() plt.close() # Plot the decision boundary plot_decision_boundary(lambda x: predict(parameters, x.T,keep_prob=0.7), train_X, train_Y,str(0.6),figure1="fig12.png")  3.1b Dropout: Spiral data – Python import numpy as np import matplotlib import matplotlib.pyplot as plt import sklearn.linear_model import pandas as pd import sklearn import sklearn.datasets exec(open("DLfunctions61.py").read()) # Create an input data set - Taken from CS231n Convolutional Neural networks, # http://cs231n.github.io/neural-networks-case-study/ N = 100 # number of points per class D = 2 # dimensionality K = 3 # number of classes X = np.zeros((N*K,D)) # data matrix (each row = single example) y = np.zeros(N*K, dtype='uint8') # class labels for j in range(K): ix = range(N*j,N*(j+1)) r = np.linspace(0.0,1,N) # radius t = np.linspace(j*4,(j+1)*4,N) + np.random.randn(N)*0.2 # theta X[ix] = np.c_[r*np.sin(t), r*np.cos(t)] y[ix] = j # Plot the data plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap=plt.cm.Spectral) plt.clf() plt.close() layersDimensions = [2,100,3] y1=y.reshape(-1,1).T # Train a deep learning network parameters = L_Layer_DeepModel(X.T, y1, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 0.6,keep_prob=0.5, num_iterations = 5000, print_cost = True,figure="fig13.png") plt.clf() plt.close() W1=parameters['W1'] b1=parameters['b1'] W2=parameters['W2'] b2=parameters['b2'] #Plot decision boundary plot_decision_boundary1(X, y1,W1,b1,W2,b2,figure2="fig14.png") 3.2a Dropout: Circles data – R source("DLfunctions61.R") #Load data df=read.csv("circles.csv",header=FALSE) z <- as.matrix(read.csv("circles.csv",header=FALSE)) x <- z[,1:2] y <- z[,3] X <- t(x) Y <- t(y) layersDimensions = c(2,11,1) # Train a deep learning network retvals = L_Layer_DeepModel(X, Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.5, keep_prob=0.8, numIterations = 9000, initType="default", print_cost = True) # Plot the decision boundary plotDecisionBoundary(z,retvals,keep_prob=0.6, hiddenActivationFunc="relu",0.5) 3.2b Dropout: Spiral data – R # Read the spiral dataset source("DLfunctions61.R") # Load data Z <- as.matrix(read.csv("spiral.csv",header=FALSE)) # Setup the data X <- Z[,1:2] y <- Z[,3] X <- t(X) Y <- t(y) layersDimensions = c(2, 15, 6, 3) # Train a deep learning network retvals = L_Layer_DeepModel(X, Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 5.1, keep_prob=0.9, numIterations = 9000, print_cost = True) parameters<-retvals$parameters
#Plot decision boundary
plotDecisionBoundary1(Z,parameters)

3.3a Dropout: Circles data – Octave

data=csvread("circles.csv");

X=data(:,1:2);
Y=data(:,3);
layersDimensions = [2 11  1]; #tanh=-0.5(ok), #relu=0.1 best!

# Train a deep learning network
[weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions,
hiddenActivationFunc='relu',
outputActivationFunc="sigmoid",
learningRate = 0.5,
lambd=0,
keep_prob=0.8,
numIterations = 10000,
initType="default");
plotCostVsIterations(10000,costs)
#Plot decision boundary
plotDecisionBoundary1(data,weights, biases,keep_prob=1, hiddenActivationFunc="relu")


3.3b Dropout  Spiral data – Octave

source("DL61functions.m")

# Setup the data
X=data(:,1:2);
Y=data(:,3);

layersDimensions = [numFeats numHidden  numOutput];
# Train a deep learning network
[weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions,
hiddenActivationFunc='relu',
outputActivationFunc="softmax",
learningRate = 0.1,
lambd=0,
keep_prob=0.8,
numIterations = 10000);

plotCostVsIterations(10000,costs)
#Plot decision boundary
plotDecisionBoundary1(data,weights, biases,keep_prob=1, hiddenActivationFunc="relu")


Note: The Python, R and Octave code can be cloned/downloaded from Github at DeepLearning-Part6
Conclusion
This post further enhances my earlier L-Layer generic implementation of a Deep Learning network to include options for initialization techniques, L2 regularization or dropout regularization

To see all posts click Index of posts

Deep Learning from first principles in Python, R and Octave – Part 5

Introduction

a. A robot may not injure a human being or, through inaction, allow a human being to come to harm.
b. A robot must obey orders given it by human beings except where such orders would conflict with the First Law.
c. A robot must protect its own existence as long as such protection does not conflict with the First or Second Law.

      Isaac Asimov's Three Laws of Robotics 

Any sufficiently advanced technology is indistinguishable from magic.

      Arthur C Clarke.   

In this 5th part on Deep Learning from first Principles in Python, R and Octave, I solve the MNIST data set of handwritten digits (shown below), from the basics. To do this, I construct a L-Layer, vectorized Deep Learning implementation in Python, R and Octave from scratch and classify the  MNIST data set. The MNIST training data set  contains 60000 handwritten digits from 0-9, and a test set of 10000 digits. MNIST, is a popular dataset for running Deep Learning tests, and has been rightfully termed as the ‘drosophila’ of Deep Learning, by none other than the venerable Prof Geoffrey Hinton.

The ‘Deep Learning from first principles in Python, R and Octave’ series, so far included  Part 1 , where I had implemented logistic regression as a simple Neural Network. Part 2 implemented the most elementary neural network with 1 hidden layer, but  with any number of activation units in that layer, and a sigmoid activation at the output layer.

This post, ‘Deep Learning from first principles in Python, R and Octave – Part 5’ largely builds upon Part3. in which I implemented a multi-layer Deep Learning network, with an arbitrary number of hidden layers and activation units per hidden layer and with the output layer was based on the sigmoid unit, for binary classification. In Part 4, I derive the Jacobian of a Softmax, the Cross entropy loss and the gradient equations for a multi-class Softmax classifier. I also  implement a simple Neural Network using Softmax classifications in Python, R and Octave.

In this post I combine Part 3 and Part 4 to to build a L-layer Deep Learning network, with arbitrary number of hidden layers and hidden units, which can do both binary (sigmoid) and multi-class (softmax) classification.

The generic, vectorized L-Layer Deep Learning Network implementations in Python, R and Octave can be cloned/downloaded from GitHub at DeepLearning-Part5. This implementation allows for arbitrary number of hidden layers and hidden layer units. The activation function at the hidden layers can be one of sigmoid, relu and tanh (will be adding leaky relu soon). The output activation can be used for binary classification with the ‘sigmoid’, or multi-class classification with ‘softmax’. Feel free to download and play around with the code!

I thought the exercise of combining the two parts(Part 3, & Part 4)  would be a breeze. But it was anything but. Incorporating a Softmax classifier into the generic L-Layer Deep Learning model was a challenge. Moreover I found that I could not use the gradient descent on 60,000 training samples as my laptop ran out of memory. So I had to implement Stochastic Gradient Descent (SGD) for Python, R and Octave. In addition, I had to also implement the numerically stable version of Softmax, as the softmax and its derivative would result in NaNs.

Numerically stable Softmax

The Softmax function $S_{j} =\frac{e^{Z_{j}}}{\sum_{i}^{k}e^{Z_{i}}}$ can be numerically unstable because of the division of large exponentials.  To handle this problem we have to implement stable Softmax function as below

$S_{j} =\frac{e^{Z_{j}}}{\sum_{i}^{k}e^{Z_{i}}}$
$S_{j} =\frac{e^{Z_{j}}}{\sum_{i}^{k}e^{Z_{i}}} = \frac{Ce^{Z_{j}}}{C\sum_{i}^{k}e^{Z_{i}}} = \frac{e^{Z_{j}+log(C)}}{\sum_{i}^{k}e^{Z_{i}+log(C)}}$
Therefore $S_{j} = \frac{e^{Z_{j}+ D}}{\sum_{i}^{k}e^{Z_{i}+ D}}$
Here ‘D’ can be anything. A common choice is
$D=-max(Z_{1},Z_{2},... Z_{k})$

Here is the stable Softmax implementation in Python

# A numerically stable Softmax implementation
def stableSoftmax(Z):
#Compute the softmax of vector x in a numerically stable way.
shiftZ = Z.T - np.max(Z.T,axis=1).reshape(-1,1)
exp_scores = np.exp(shiftZ)
# normalize them for each example
A = exp_scores / np.sum(exp_scores, axis=1, keepdims=True)
cache=Z
return A,cache


While trying to create a L-Layer generic Deep Learning network in the 3 languages, I found it useful to ensure that the model executed correctly on smaller datasets.  You can run into numerous problems while setting up the matrices, which becomes extremely difficult to debug. So in this post, I run the model on 2 smaller data for sets used in my earlier posts(Part 3 & Part4) , in each of the languages, before running the generic model on MNIST.

Here is a fair warning. if you think you can dive directly into Deep Learning, with just some basic knowledge of Machine Learning, you are bound to run into serious issues. Moreover, your knowledge will be incomplete. It is essential that you have a good grasp of Machine and Statistical Learning, the different algorithms, the measures and metrics for selecting the models etc.It would help to be conversant with all the ML models, ML concepts, validation techniques, classification measures  etc. Check out the internet/books for background.

Incidentally you could also check out my book, Practical Machine Learning in R and Python – Machine Learning in Stereo, available on Amazon. In this book, I implement different ML algorithms,  regularization techniques, model selection, ensemble methods, key measures and metrics. My book is compact, minimalist  and includes implementations of many ML algorithms in R and Python, and can be used as a handy reference.

1. Random dataset with Sigmoid activation – Python

This random data with 9 clusters, was used in my post Deep Learning from first principles in Python, R and Octave – Part 3 , and was used to test the complete L-layer Deep Learning network with Sigmoid activation.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import make_classification, make_blobs
exec(open("DLfunctions51.py").read()) # Cannot import in Rmd.
# Create a random data set with 9 centeres
X1, Y1 = make_blobs(n_samples = 400, n_features = 2, centers = 9,cluster_std = 1.3, random_state =4)

#Create 2 classes
Y1=Y1.reshape(400,1)
Y1 = Y1 % 2
X2=X1.T
Y2=Y1.T
# Set the dimensions of L -layer DL network
layersDimensions = [2, 9, 9,1] #  4-layer model
# Execute DL network with hidden activation=relu and sigmoid output function
parameters = L_Layer_DeepModel(X2, Y2, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid",learningRate = 0.3,num_iterations = 2500, print_cost = True)

2. Spiral dataset with Softmax activation – Python

The Spiral data was used in my post Deep Learning from first principles in Python, R and Octave – Part 4 and was used to test the complete L-layer Deep Learning network with multi-class Softmax activation at the output layer

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import make_classification, make_blobs

# Create an input data set - Taken from CS231n Convolutional Neural networks
# http://cs231n.github.io/neural-networks-case-study/
N = 100 # number of points per class
D = 2 # dimensionality
K = 3 # number of classes
X = np.zeros((N*K,D)) # data matrix (each row = single example)
y = np.zeros(N*K, dtype='uint8') # class labels
for j in range(K):
ix = range(N*j,N*(j+1))
t = np.linspace(j*4,(j+1)*4,N) + np.random.randn(N)*0.2 # theta
X[ix] = np.c_[r*np.sin(t), r*np.cos(t)]
y[ix] = j

X1=X.T
Y1=y.reshape(-1,1).T
numHidden=100 # No of hidden units in hidden layer
numFeats= 2 # dimensionality
numOutput = 3 # number of classes
# Set the dimensions of the layers
layersDimensions=[numFeats,numHidden,numOutput]
parameters = L_Layer_DeepModel(X1, Y1, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax",learningRate = 0.6,num_iterations = 9000, print_cost = True)
## Cost after iteration 0: 1.098759
## Cost after iteration 1000: 0.112666
## Cost after iteration 2000: 0.044351
## Cost after iteration 3000: 0.027491
## Cost after iteration 4000: 0.021898
## Cost after iteration 5000: 0.019181
## Cost after iteration 6000: 0.017832
## Cost after iteration 7000: 0.017452
## Cost after iteration 8000: 0.017161

3. MNIST dataset with Softmax activation – Python

In the code below, I execute Stochastic Gradient Descent on the MNIST training data of 60000. I used a mini-batch size of 1000. Python takes about 40 minutes to crunch the data. In addition I also compute the Confusion Matrix and other metrics like Accuracy, Precision and Recall for the MNIST data set. I get an accuracy of 0.93 on the MNIST test set. This accuracy can be improved by choosing more hidden layers or more hidden units and possibly also tweaking the learning rate and the number of epochs.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import math
from sklearn.datasets import make_classification, make_blobs
from sklearn.metrics import confusion_matrix
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
# Read the MNIST training and test sets
# Create labels and pixel arrays
lbls=[]
pxls=[]
print(len(training))
#for i in range(len(training)):
for i in range(60000):
l,p=training[i]
lbls.append(l)
pxls.append(p)
labels= np.array(lbls)
pixels=np.array(pxls)
y=labels.reshape(-1,1)
X=pixels.reshape(pixels.shape[0],-1)
X1=X.T
Y1=y.T
# Set the dimensions of the layers. The MNIST data is 28x28 pixels= 784
# Hence input layer is 784. For the 10 digits the Softmax classifier
# has to handle 10 outputs
layersDimensions=[784, 15,9,10] # Works very well,lr=0.01,mini_batch =1000, total=20000
np.random.seed(1)
costs = []
# Run Stochastic Gradient Descent with Learning Rate=0.01, mini batch size=1000
# number of epochs=3000
parameters = L_Layer_DeepModel_SGD(X1, Y1, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax",learningRate = 0.01 ,mini_batch_size =1000, num_epochs = 3000, print_cost = True)

# Compute the Confusion Matrix on Training set
# Compute the training accuracy, precision and recall
proba=predict_proba(parameters, X1,outputActivationFunc="softmax")
#A2, cache = forwardPropagationDeep(X1, parameters)
#proba=np.argmax(A2, axis=0).reshape(-1,1)
a=confusion_matrix(Y1.T,proba)
print(a)
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
print('Accuracy: {:.2f}'.format(accuracy_score(Y1.T, proba)))
print('Precision: {:.2f}'.format(precision_score(Y1.T, proba,average="micro")))
print('Recall: {:.2f}'.format(recall_score(Y1.T, proba,average="micro")))

lbls=[]
pxls=[]
print(len(test))
for i in range(10000):
l,p=test[i]
lbls.append(l)
pxls.append(p)
testLabels= np.array(lbls)
testPixels=np.array(pxls)
ytest=testLabels.reshape(-1,1)
Xtest=testPixels.reshape(testPixels.shape[0],-1)
X1test=Xtest.T
Y1test=ytest.T

# Compute the Confusion Matrix on Test set
# Compute the test accuracy, precision and recall
probaTest=predict_proba(parameters, X1test,outputActivationFunc="softmax")
#A2, cache = forwardPropagationDeep(X1, parameters)
#proba=np.argmax(A2, axis=0).reshape(-1,1)
a=confusion_matrix(Y1test.T,probaTest)
print(a)
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
print('Accuracy: {:.2f}'.format(accuracy_score(Y1test.T, probaTest)))
print('Precision: {:.2f}'.format(precision_score(Y1test.T, probaTest,average="micro")))
print('Recall: {:.2f}'.format(recall_score(Y1test.T, probaTest,average="micro")))

##1.  Confusion Matrix of Training set
0     1    2    3    4    5    6    7    8    9
## [[5854    0   19    2   10    7    0    1   24    6]
##  [   1 6659   30   10    5    3    0   14   20    0]
##  [  20   24 5805   18    6   11    2   32   37    3]
##  [   5    4  175 5783    1   27    1   58   60   17]
##  [   1   21    9    0 5780    0    5    2   12   12]
##  [  29    9   21  224    6 4824   18   17  245   28]
##  [   5    4   22    1   32   12 5799    0   43    0]
##  [   3   13  148  154   18    3    0 5883    4   39]
##  [  11   34   30   21   13   16    4    7 5703   12]
##  [  10    4    1   32  135   14    1   92  134 5526]]

##2. Accuracy, Precision, Recall of  Training set
## Accuracy: 0.96
## Precision: 0.96
## Recall: 0.96

##3. Confusion Matrix of Test set
0     1    2    3    4    5    6    7    8    9
## [[ 954    1    8    0    3    3    2    4    4    1]
##  [   0 1107    6    5    0    0    1    2   14    0]
##  [  11    7  957   10    5    0    5   20   16    1]
##  [   2    3   37  925    3   13    0    8   18    1]
##  [   2    6    1    1  944    0    7    3    4   14]
##  [  12    5    4   45    2  740   24    8   42   10]
##  [   8    4    4    2   16    9  903    0   12    0]
##  [   4   10   27   18    5    1    0  940    1   22]
##  [  11   13    6   13    9   10    7    2  900    3]
##  [   8    5    1    7   50    7    0   20   29  882]]
##4. Accuracy, Precision, Recall of  Training set
## Accuracy: 0.93
## Precision: 0.93
## Recall: 0.93

4. Random dataset with Sigmoid activation – R code

This is the random data set used in the Python code above which was saved as a CSV. The code is used to test a L -Layer DL network with Sigmoid Activation in R.

source("DLfunctions5.R")
# Read the random data set
x <- z[,1:2]
y <- z[,3]
X <- t(x)
Y <- t(y)
# Set the dimensions of the  layer
layersDimensions = c(2, 9, 9,1)

# Run Gradient Descent on the data set with relu hidden unit activation
# sigmoid activation unit in the output layer
retvals = L_Layer_DeepModel(X, Y, layersDimensions,
hiddenActivationFunc='relu',
outputActivationFunc="sigmoid",
learningRate = 0.3,
numIterations = 5000,
print_cost = True)
#Plot the cost vs iterations
iterations <- seq(0,5000,1000)
costs=retvals$costs df=data.frame(iterations,costs) ggplot(df,aes(x=iterations,y=costs)) + geom_point() + geom_line(color="blue") + ggtitle("Costs vs iterations") + xlab("Iterations") + ylab("Loss") 5. Spiral dataset with Softmax activation – R The spiral data set used in the Python code above, is reused to test multi-class classification with Softmax. source("DLfunctions5.R") Z <- as.matrix(read.csv("spiral.csv",header=FALSE)) # Setup the data X <- Z[,1:2] y <- Z[,3] X <- t(X) Y <- t(y) # Initialize number of features, number of hidden units in hidden layer and # number of classes numFeats<-2 # No features numHidden<-100 # No of hidden units numOutput<-3 # No of classes # Set the layer dimensions layersDimensions = c(numFeats,numHidden,numOutput) # Perform gradient descent with relu activation unit for hidden layer # and softmax activation in the output retvals = L_Layer_DeepModel(X, Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 0.5, numIterations = 9000, print_cost = True) #Plot cost vs iterations iterations <- seq(0,9000,1000) costs=retvals$costs
df=data.frame(iterations,costs)
ggplot(df,aes(x=iterations,y=costs)) + geom_point() + geom_line(color="blue") +
ggtitle("Costs vs iterations") + xlab("Iterations") + ylab("Costs")

6. MNIST dataset with Softmax activation – R

The code below executes a L – Layer Deep Learning network with Softmax output activation, to classify the 10 handwritten digits from MNIST with Stochastic Gradient Descent. The entire 60000 data set was used to train the data. R takes almost 8 hours to process this data set with a mini-batch size of 1000.  The use of ‘for’ loops is limited to iterating through epochs, mini batches and for creating the mini batches itself. All other code is vectorized. Yet, it seems to crawl. Most likely the use of ‘lists’ in R, to return multiple values is performance intensive. Some day, I will try to profile the code, and see where the issue is. However the code works!

Having said that, the Confusion Matrix in R dumps a lot of interesting statistics! There is a bunch of statistical measures for each class. For e.g. the Balanced Accuracy for the digits ‘6’ and ‘9’ is around 50%. Looks like, the classifier is confused by the fact that 6 is inverted 9 and vice-versa. The accuracy on the Test data set is just around 75%. I could have played around with the number of layers, number of hidden units, learning rates, epochs etc to get a much higher accuracy. But since each test took about 8+ hours, I may work on this, some other day!

source("DLfunctions5.R")
source("mnist.R")
show_digit(train$x[2,]) #Set the layer dimensions layersDimensions=c(784, 15,9, 10) # Works at 1500 x <- t(train$x)
X <- x[,1:60000]
y <-train$y y1 <- y[1:60000] y2 <- as.matrix(y1) Y=t(y2) # Execute Stochastic Gradient Descent on the entire training set # with Softmax activation retvalsSGD= L_Layer_DeepModel_SGD(X, Y, layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 0.01, mini_batch_size = 1000, num_epochs = 1, print_cost = True) # Compute the Confusion Matrix library(caret) library(e1071) predictions=predictProba(retvalsSGD[['parameters']], X,hiddenActivationFunc='relu', outputActivationFunc="softmax") confusionMatrix(predictions,Y) # Confusion Matrix on the Training set > confusionMatrix(predictions,Y) Confusion Matrix and Statistics Reference Prediction 0 1 2 3 4 5 6 7 8 9 0 5738 1 21 5 16 17 7 15 9 43 1 5 6632 21 24 25 3 2 33 13 392 2 12 32 5747 106 25 28 3 27 44 4779 3 0 27 12 5715 1 21 1 20 1 13 4 10 5 21 18 5677 9 17 30 15 166 5 142 21 96 136 93 5306 5884 43 60 413 6 0 0 0 0 0 0 0 0 0 0 7 6 9 13 13 3 4 0 6085 0 55 8 8 12 7 43 1 32 2 7 5703 69 9 2 3 20 71 1 1 2 5 6 19 Overall Statistics Accuracy : 0.777 95% CI : (0.7737, 0.7804) No Information Rate : 0.1124 P-Value [Acc > NIR] : < 2.2e-16 Kappa : 0.7524 Mcnemar's Test P-Value : NA Statistics by Class: Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5 Class: 6 Sensitivity 0.96877 0.9837 0.96459 0.93215 0.97176 0.97879 0.00000 Specificity 0.99752 0.9903 0.90644 0.99822 0.99463 0.87380 1.00000 Pos Pred Value 0.97718 0.9276 0.53198 0.98348 0.95124 0.43513 NaN Neg Pred Value 0.99658 0.9979 0.99571 0.99232 0.99695 0.99759 0.90137 Prevalence 0.09872 0.1124 0.09930 0.10218 0.09737 0.09035 0.09863 Detection Rate 0.09563 0.1105 0.09578 0.09525 0.09462 0.08843 0.00000 Detection Prevalence 0.09787 0.1192 0.18005 0.09685 0.09947 0.20323 0.00000 Balanced Accuracy 0.98314 0.9870 0.93551 0.96518 0.98319 0.92629 0.50000 Class: 7 Class: 8 Class: 9 Sensitivity 0.9713 0.97471 0.0031938 Specificity 0.9981 0.99666 0.9979464 Pos Pred Value 0.9834 0.96924 0.1461538 Neg Pred Value 0.9967 0.99727 0.9009521 Prevalence 0.1044 0.09752 0.0991500 Detection Rate 0.1014 0.09505 0.0003167 Detection Prevalence 0.1031 0.09807 0.0021667 Balanced Accuracy 0.9847 0.98568 0.5005701  # Confusion Matrix on the Training set xtest <- t(test$x) Xtest <- xtest[,1:10000] ytest <-test$y ytest1 <- ytest[1:10000] ytest2 <- as.matrix(ytest1) Ytest=t(ytest2)  Confusion Matrix and Statistics Reference Prediction 0 1 2 3 4 5 6 7 8 9 0 950 2 2 3 0 6 9 4 7 6 1 3 1110 4 2 9 0 3 12 5 74 2 2 6 965 21 9 14 5 16 12 789 3 1 2 9 908 2 16 0 21 2 6 4 0 1 9 5 938 1 8 6 8 39 5 19 5 25 35 20 835 929 8 54 67 6 0 0 0 0 0 0 0 0 0 0 7 4 4 7 10 2 4 0 952 5 6 8 1 5 8 14 2 16 2 3 876 21 9 0 0 3 12 0 0 2 6 5 1 Overall Statistics Accuracy : 0.7535 95% CI : (0.7449, 0.7619) No Information Rate : 0.1135 P-Value [Acc > NIR] : < 2.2e-16 Kappa : 0.7262 Mcnemar's Test P-Value : NA Statistics by Class: Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5 Class: 6 Sensitivity 0.9694 0.9780 0.9351 0.8990 0.9552 0.9361 0.0000 Specificity 0.9957 0.9874 0.9025 0.9934 0.9915 0.8724 1.0000 Pos Pred Value 0.9606 0.9083 0.5247 0.9390 0.9241 0.4181 NaN Neg Pred Value 0.9967 0.9972 0.9918 0.9887 0.9951 0.9929 0.9042 Prevalence 0.0980 0.1135 0.1032 0.1010 0.0982 0.0892 0.0958 Detection Rate 0.0950 0.1110 0.0965 0.0908 0.0938 0.0835 0.0000 Detection Prevalence 0.0989 0.1222 0.1839 0.0967 0.1015 0.1997 0.0000 Balanced Accuracy 0.9825 0.9827 0.9188 0.9462 0.9733 0.9043 0.5000 Class: 7 Class: 8 Class: 9 Sensitivity 0.9261 0.8994 0.0009911 Specificity 0.9953 0.9920 0.9968858 Pos Pred Value 0.9577 0.9241 0.0344828 Neg Pred Value 0.9916 0.9892 0.8989068 Prevalence 0.1028 0.0974 0.1009000 Detection Rate 0.0952 0.0876 0.0001000 Detection Prevalence 0.0994 0.0948 0.0029000 Balanced Accuracy 0.9607 0.9457 0.4989384  7. Random dataset with Sigmoid activation – Octave The Octave code below uses the random data set used by Python. The code below implements a L-Layer Deep Learning with Sigmoid Activation.  source("DL5functions.m") # Read the data data=csvread("data.csv"); X=data(:,1:2); Y=data(:,3); #Set the layer dimensions layersDimensions = [2 9 7 1]; #tanh=-0.5(ok), #relu=0.1 best! # Perform gradient descent [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="sigmoid", learningRate = 0.1, numIterations = 10000); # Plot cost vs iterations plotCostVsIterations(10000,costs);  8. Spiral dataset with Softmax activation – Octave The code below uses the spiral data set used by Python above. The code below implements a L-Layer Deep Learning with Softmax Activation. # Read the data data=csvread("spiral.csv"); # Setup the data X=data(:,1:2); Y=data(:,3); # Set the number of features, number of hidden units in hidden layer and number of classess numFeats=2; #No features numHidden=100; # No of hidden units numOutput=3; # No of classes # Set the layer dimensions layersDimensions = [numFeats numHidden numOutput]; #Perform gradient descent with softmax activation unit [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 0.1, numIterations = 10000);  9. MNIST dataset with Softmax activation – Octave The code below implements a L-Layer Deep Learning Network in Octave with Softmax output activation unit, for classifying the 10 handwritten digits in the MNIST dataset. Unfortunately, Octave can only index to around 10000 training at a time, and I was getting an error ‘error: out of memory or dimension too large for Octave’s index type error: called from…’, when I tried to create a batch size of 20000. So I had to come with a work around to create a batch size of 10000 (randomly) and then use a mini-batch of 1000 samples and execute Stochastic Gradient Descent. The performance was good. Octave takes about 15 minutes, on a batch size of 10000 and a mini batch of 1000. I thought if the performance was not good, I could iterate through these random batches and refining the gradients as follows # Pseudo code that could be used since Octave only allows 10K batches # at a time # Randomly create weights [weights biases] = initialize_weights() for i=1:k # Create a random permutation and create a random batch permutation = randperm(10000); X=trainX(permutation,:); Y=trainY(permutation,:); # Compute weights from SGD and update weights in the next batch update [weights biases costs]=L_Layer_DeepModel_SGD(X,Y,mini_bactch=1000,weights, biases,...); ... endfor # Load the MNIST data load('./mnist/mnist.txt.gz'); #Create a random permutatation from 60K permutation = randperm(10000); disp(length(permutation)); # Use this 10K as the batch X=trainX(permutation,:); Y=trainY(permutation,:); # Set layer dimensions layersDimensions=[784, 15, 9, 10]; # Run Stochastic Gradient descent with batch size=10K and mini_batch_size=1000 [weights biases costs]=L_Layer_DeepModel_SGD(X', Y', layersDimensions, hiddenActivationFunc='relu', outputActivationFunc="softmax", learningRate = 0.01, mini_batch_size = 2000, num_epochs = 5000);  9. Final thoughts Here are some of my final thoughts after working on Python, R and Octave in this series and in other projects 1. Python, with its highly optimized numpy library, is ideally suited for creating Deep Learning Models, which have a lot of matrix manipulations. Python is a real workhorse when it comes to Deep Learning computations. 2. R is somewhat clunky in comparison to its cousin Python in handling matrices or in returning multiple values. But R’s statistical libraries, dplyr, and ggplot are really superior to the Python peers. Also, I find R handles dataframes, much better than Python. 3. Octave is a no-nonsense,minimalist language which is very efficient in handling matrices. It is ideally suited for implementing Machine Learning and Deep Learning from scratch. But Octave has its problems and cannot handle large matrix sizes, and also lacks the statistical libaries of R and Python. They possibly exist in its sibling, Matlab Feel free to clone/download the code from GitHub at DeepLearning-Part5. Conclusion Building a Deep Learning Network from scratch is quite challenging, time-consuming but nevertheless an exciting task. While the statements in the different languages for manipulating matrices, summing up columns, finding columns which have ones don’t take more than a single statement, extreme care has to be taken to ensure that the statements work well for any dimension. The lessons learnt from creating L -Layer Deep Learning network are many and well worth it. Give it a try! Hasta la vista! I’ll be back, so stick around! Watch this space! To see all posts click Index of Posts Deep Learning from first principles in Python, R and Octave – Part 4 In this 4th post of my series on Deep Learning from first principles in Python, R and Octave – Part 4, I explore the details of creating a multi-class classifier using the Softmax activation unit in a neural network. The earlier posts in this series were 1. Deep Learning from first principles in Python, R and Octave – Part 1. In this post I implemented logistic regression as a simple Neural Network in vectorized Python, R and Octave 2. Deep Learning from first principles in Python, R and Octave – Part 2. This 2nd part implemented the most elementary neural network with 1 hidden layer and any number of activation units in the hidden layer with sigmoid activation at the output layer 3. Deep Learning from first principles in Python, R and Octave – Part 3. The 3rd implemented a multi-layer Deep Learning network with an arbitrary number if hidden layers and activation units per hidden layer. The output layer was for binary classification which was based on the sigmoid unit. This multi-layer deep network was implemented in vectorized Python, R and Octave This 4th part takes a swing at multi-class classification and uses the Softmax as the activation unit in the output layer. Inclusion of the Softmax activation unit in the activation layer requires us to compute the derivative of Softmax, or rather the “Jacobian” of the Softmax function, besides also computing the log loss for this Softmax activation during back propagation. Since the derivation of the Jacobian of a Softmax and the computation of the Cross Entropy/log loss is very involved, I have implemented a basic neural network with just 1 hidden layer with the Softmax activation at the output layer. I also perform multi-class classification based on the ‘spiral’ data set from CS231n Convolutional Neural Networks Stanford course, to test the performance and correctness of the implementations in Python, R and Octave. You can clone download the code for the Python, R and Octave implementations from Github at Deep Learning – Part 4 The Softmax function takes an N dimensional vector as input and generates a N dimensional vector as output. The Softmax function is given by $S_{j}= \frac{e_{j}}{\sum_{i}^{N}e_{k}}$ There is a probabilistic interpretation of the Softmax, since the sum of the Softmax values of a set of vectors will always add up to 1, given that each Softmax value is divided by the total of all values. As mentioned earlier, the Softmax takes a vector input and returns a vector of outputs. For e.g. the Softmax of a vector a=[1, 3, 6] is another vector S=[0.0063,0.0471,0.9464]. Notice that vector output is proportional to the input vector. Also, taking the derivative of a vector by another vector, is known as the Jacobian. By the way, The Matrix Calculus You Need For Deep Learning by Terence Parr and Jeremy Howard, is very good paper that distills all the main mathematical concepts for Deep Learning in one place. Let us take a simple 2 layered neural network with just 2 activation units in the hidden layer is shown below $Z_{1}^{1} =W_{11}^{1}x_{1} + W_{21}^{1}x_{2} + b_{1}^{1}$ $Z_{2}^{1} =W_{12}^{1}x_{1} + W_{22}^{1}x_{2} + b_{2}^{1}$ and $A_{1}^{1} = g'(Z_{1}^{1})$ $A_{2}^{1} = g'(Z_{2}^{1})$ where g'() is the activation unit in the hidden layer which can be a relu, sigmoid or a tanh function Note: The superscript denotes the layer. The above denotes the equation for layer 1 of the neural network. For layer 2 with the Softmax activation, the equations are $Z_{1}^{2} =W_{11}^{2}x_{1} + W_{21}^{2}x_{2} + b_{1}^{2}$ $Z_{2}^{2} =W_{12}^{2}x_{1} + W_{22}^{2}x_{2} + b_{2}^{2}$ and $A_{1}^{2} = S(A_{1}^{1})$ $A_{2}^{2} = S(A_{2}^{1})$ where S() is the Softmax activation function $S=\begin{pmatrix} S(Z_{1}^{2})\\ S(Z_{2}^{2}) \end{pmatrix}$ $S=\begin{pmatrix} \frac{e^{Z1}}{e^{Z1}+e^{Z2}}\\ \frac{e^{Z2}}{e^{Z1}+e^{Z2}} \end{pmatrix}$ The Jacobian of the softmax ‘S’ is given by $\begin{pmatrix} \frac {\partial S_{1}}{\partial Z_{1}} & \frac {\partial S_{1}}{\partial Z_{2}}\\ \frac {\partial S_{2}}{\partial Z_{1}} & \frac {\partial S_{2}}{\partial Z_{2}} \end{pmatrix}$ $\begin{pmatrix} \frac{\partial}{\partial Z_{1}} \frac {e^{Z1}}{e^{Z1}+ e^{Z2}} & \frac{\partial}{\partial Z_{2}} \frac {e^{Z1}}{e^{Z1}+ e^{Z2}}\\ \frac{\partial}{\partial Z_{1}} \frac {e^{Z2}}{e^{Z1}+ e^{Z2}} & \frac{\partial}{\partial Z_{2}} \frac {e^{Z2}}{e^{Z1}+ e^{Z2}} \end{pmatrix}$ – (A) Now the ‘division-rule’ of derivatives is as follows. If u and v are functions of x, then $\frac{d}{dx} \frac {u}{v} =\frac {vdu -udv}{v^{2}}$ Using this to compute each element of the above Jacobian matrix, we see that when i=j we have $\frac {\partial}{\partial Z1}\frac{e^{Z1}}{e^{Z1}+e^{Z2}} = \frac {\sum e^{Z1} - e^{Z1^{2}}}{\sum ^{2}}$ and when $i \neq j$ $\frac {\partial}{\partial Z1}\frac{e^{Z2}}{e^{Z1}+e^{Z2}} = \frac {0 - e^{z1}e^{Z2}}{\sum ^{2}}$ This is of the general form $\frac {\partial S_{j}}{\partial z_{i}} = S_{i}( 1-S_{j})$ when i=j and $\frac {\partial S_{j}}{\partial z_{i}} = -S_{i}S_{j}$ when $i \neq j$ Note: Since the Softmax essentially gives the probability the following notation is also used $\frac {\partial p_{j}}{\partial z_{i}} = p_{i}( 1-p_{j})$ when i=j and $\frac {\partial p_{j}}{\partial z_{i}} = -p_{i}p_{j} when i \neq j$ If you throw the “Kronecker delta” into the equation, then the above equations can be expressed even more concisely as $\frac {\partial p_{j}}{\partial z_{i}} = p_{i} (\delta_{ij} - p_{j})$ where $\delta_{ij} = 1$ when i=j and 0 when $i \neq j$ This reduces the Jacobian of the simple 2 output softmax vectors equation (A) as $\begin{pmatrix} p_{1}(1-p_{1}) & -p_{1}p_{2} \\ -p_{2}p_{1} & p_{2}(1-p_{2}) \end{pmatrix}$ The loss of Softmax is given by $L = -\sum y_{i} log(p_{i})$ For the 2 valued Softmax output this is $\frac {dL}{dp1} = -\frac {y_{1}}{p_{1}}$ $\frac {dL}{dp2} = -\frac {y_{2}}{p_{2}}$ Using the chain rule we can write $\frac {\partial L}{\partial w_{pq}} = \sum _{i}\frac {\partial L}{\partial p_{i}} \frac {\partial p_{i}}{\partial w_{pq}}$ (1) and $\frac {\partial p_{i}}{\partial w_{pq}} = \sum _{k}\frac {\partial p_{i}}{\partial z_{k}} \frac {\partial z_{k}}{\partial w_{pq}}$ (2) In expanded form this is $\frac {\partial L}{\partial w_{pq}} = \sum _{i}\frac {\partial L}{\partial p_{i}} \sum _{k}\frac {\partial p_{i}}{\partial z_{k}} \frac {\partial z_{k}}{\partial w_{pq}}$ Also $\frac {\partial L}{\partial Z_{i}} =\sum _{i} \frac {\partial L}{\partial p} \frac {\partial p}{\partial Z_{i}}$ Therefore $\frac {\partial L}{\partial Z_{1}} =\frac {\partial L}{\partial p_{1}} \frac {\partial p_{1}}{\partial Z_{1}} +\frac {\partial L}{\partial p_{2}} \frac {\partial p_{2}}{\partial Z_{1}}$ $\frac {\partial L}{\partial z_{1}}=-\frac {y1}{p1} p1(1-p1) - \frac {y2}{p2}*(-p_{2}p_{1})$ Since $\frac {\partial p_{j}}{\partial z_{i}} = p_{i}( 1-p_{j})$ when i=j and $\frac {\partial p_{j}}{\partial z_{i}} = -p_{i}p_{j}$ when $i \neq j$ which simplifies to $\frac {\partial L}{\partial Z_{1}} = -y_{1} + y_{1}p_{1} + y_{2}p_{1} =$ $p_{1}\sum (y_{1} + y_2) - y_{1}$ $\frac {\partial L}{\partial Z_{1}}= p_{1} - y_{1}$ Since $\sum_{i} y_{i} =1$ Similarly $\frac {\partial L}{\partial Z_{2}} =\frac {\partial L}{\partial p_{1}} \frac {\partial p_{1}}{\partial Z_{2}} +\frac {\partial L}{\partial p_{2}} \frac {\partial p_{2}}{\partial Z_{2}}$ $\frac {\partial L}{\partial z_{2}}=-\frac {y1}{p1}*(p_{1}p_{2}) - \frac {y2}{p2}*p_{2}(1-p_{2})$ $y_{1}p_{2} + y_{2}p_{2} - y_{2}$ $\frac {\partial L}{\partial Z_{2}} =p_{2}\sum (y_{1} + y_2) - y_{2}\\ = p_{2} - y_{2}$ In general this is of the form $\frac {\partial L}{\partial z_{i}} = p_{i} -y_{i}$ For e.g if the probabilities computed were p=[0.1, 0.7, 0.2] then this implies that the class with probability 0.7 is the likely class. This would imply that the ‘One hot encoding’ for yi would be yi=[0,1,0] therefore the gradient pi-yi = [0.1,-0.3,0.2] <strong>Note: Further, we could extend this derivation for a Softmax activation output that outputs 3 classes $S=\begin{pmatrix} \frac{e^{z1}}{e^{z1}+e^{z2}+e^{z3}}\\ \frac{e^{z2}}{e^{z1}+e^{z2}+e^{z3}} \\ \frac{e^{z3}}{e^{z1}+e^{z2}+e^{z3}} \end{pmatrix}$ We could derive $\frac {\partial L}{\partial z1}= \frac {\partial L}{\partial p_{1}} \frac {\partial p_{1}}{\partial z_{1}} +\frac {\partial L}{\partial p_{2}} \frac {\partial p_{2}}{\partial z_{1}} +\frac {\partial L}{\partial p_{3}} \frac {\partial p_{3}}{\partial z_{1}}$ which similarly reduces to $\frac {\partial L}{\partial z_{1}}=-\frac {y1}{p1} p1(1-p1) - \frac {y2}{p2}*(-p_{2}p_{1}) - \frac {y3}{p3}*(-p_{3}p_{1})$ $-y_{1}+ y_{1}p_{1} + y_{2}p_{1} + y_{3}p1 = p_{1}\sum (y_{1} + y_2 + y_3) - y_{1} = p_{1} - y_{1}$ Interestingly, despite the lengthy derivations the final result is simple and intuitive! As seen in my post ‘Deep Learning from first principles with Python, R and Octave – Part 3 the key equations for forward and backward propagation are Forward propagation equations layer 1 $Z_{1} = W_{1}X +b_{1}$ and $A_{1} = g(Z_{1})$ Forward propagation equations layer 1 $Z_{2} = W_{2}A_{1} +b_{2}$ and $A_{2} = S(Z_{2})$ Backward propagation equations layer 2 $\partial L/\partial W_{2} =\partial L/\partial Z_{2} *A_{1}$ $\partial L/\partial b_{2} =\partial L/\partial Z_{2}$ $\partial L/\partial A_{1} = \partial L/\partial Z_{2} * W_{2}$ Backward propagation equations layer 1 $\partial L/\partial W_{1} =\partial L/\partial Z_{1} *A_{0}$ $\partial L/\partial b_{1} =\partial L/\partial Z_{2}$ 2.0 Spiral data set As I mentioned earlier, I will be using the ‘spiral’ data from CS231n Convolutional Neural Networks to ensure that my vectorized implementations in Python, R and Octave are correct. Here is the ‘spiral’ data set. import numpy as np import matplotlib.pyplot as plt import os os.chdir("C:/junk/dl-4/dl-4") exec(open("././DLfunctions41.py").read()) # Create an input data set - Taken from CS231n Convolutional Neural networks # http://cs231n.github.io/neural-networks-case-study/ N = 100 # number of points per class D = 2 # dimensionality K = 3 # number of classes X = np.zeros((N*K,D)) # data matrix (each row = single example) y = np.zeros(N*K, dtype='uint8') # class labels for j in range(K): ix = range(N*j,N*(j+1)) r = np.linspace(0.0,1,N) # radius t = np.linspace(j*4,(j+1)*4,N) + np.random.randn(N)*0.2 # theta X[ix] = np.c_[r*np.sin(t), r*np.cos(t)] y[ix] = j # Plot the data plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap=plt.cm.Spectral) plt.savefig("fig1.png", bbox_inches='tight') The implementations of the vectorized Python, R and Octave code are shown diagrammatically below 2.1 Multi-class classification with Softmax – Python code A simple 2 layer Neural network with a single hidden layer , with 100 Relu activation units in the hidden layer and the Softmax activation unit in the output layer is used for multi-class classification. This Deep Learning Network, plots the non-linear boundary of the 3 classes as shown below import numpy as np import matplotlib.pyplot as plt import os os.chdir("C:/junk/dl-4/dl-4") exec(open("././DLfunctions41.py").read()) # Read the input data N = 100 # number of points per class D = 2 # dimensionality K = 3 # number of classes X = np.zeros((N*K,D)) # data matrix (each row = single example) y = np.zeros(N*K, dtype='uint8') # class labels for j in range(K): ix = range(N*j,N*(j+1)) r = np.linspace(0.0,1,N) # radius t = np.linspace(j*4,(j+1)*4,N) + np.random.randn(N)*0.2 # theta X[ix] = np.c_[r*np.sin(t), r*np.cos(t)] y[ix] = j # Set the number of features, hidden units in hidden layer and number of classess numHidden=100 # No of hidden units in hidden layer numFeats= 2 # dimensionality numOutput = 3 # number of classes # Initialize the model parameters=initializeModel(numFeats,numHidden,numOutput) W1= parameters['W1'] b1= parameters['b1'] W2= parameters['W2'] b2= parameters['b2'] # Set the learning rate learningRate=0.6 # Initialize losses losses=[] # Perform Gradient descent for i in range(10000): # Forward propagation through hidden layer with Relu units A1,cache1= layerActivationForward(X.T,W1,b1,'relu') # Forward propagation through output layer with Softmax A2,cache2 = layerActivationForward(A1,W2,b2,'softmax') # No of training examples numTraining = X.shape[0] # Compute log probs. Take the log prob of correct class based on output y correct_logprobs = -np.log(A2[range(numTraining),y]) # Conpute loss loss = np.sum(correct_logprobs)/numTraining # Print the loss if i % 1000 == 0: print("iteration %d: loss %f" % (i, loss)) losses.append(loss) dA=0 # Backward propagation through output layer with Softmax dA1,dW2,db2 = layerActivationBackward(dA, cache2, y, activationFunc='softmax') # Backward propagation through hidden layer with Relu unit dA0,dW1,db1 = layerActivationBackward(dA1.T, cache1, y, activationFunc='relu') #Update paramaters with the learning rate W1 += -learningRate * dW1 b1 += -learningRate * db1 W2 += -learningRate * dW2.T b2 += -learningRate * db2.T #Plot losses vs iterations i=np.arange(0,10000,1000) plt.plot(i,losses) plt.xlabel('Iterations') plt.ylabel('Loss') plt.title('Losses vs Iterations') plt.savefig("fig2.png", bbox="tight") #Compute the multi-class Confusion Matrix from sklearn.metrics import confusion_matrix from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score # We need to determine the predicted values from the learnt data # Forward propagation through hidden layer with Relu units A1,cache1= layerActivationForward(X.T,W1,b1,'relu') # Forward propagation through output layer with Softmax A2,cache2 = layerActivationForward(A1,W2,b2,'softmax') #Compute predicted values from weights and biases yhat=np.argmax(A2, axis=1) a=confusion_matrix(y.T,yhat.T) print("Multi-class Confusion Matrix") print(a) ## iteration 0: loss 1.098507 ## iteration 1000: loss 0.214611 ## iteration 2000: loss 0.043622 ## iteration 3000: loss 0.032525 ## iteration 4000: loss 0.025108 ## iteration 5000: loss 0.021365 ## iteration 6000: loss 0.019046 ## iteration 7000: loss 0.017475 ## iteration 8000: loss 0.016359 ## iteration 9000: loss 0.015703 ## Multi-class Confusion Matrix ## [[ 99 1 0] ## [ 0 100 0] ## [ 0 1 99]] To know more details about Confusion Matrix and other related measures check out my book My book ‘Practical Machine Learning with R and Python’ on Amazon 2.2 Multi-class classification with Softmax – R code The spiral data set created with Python was saved, and is used as the input with R code. The R Neural Network seems to perform much,much slower than both Python and Octave. Not sure why! Incidentally the computation of loss and the softmax derivative are identical for both R and Octave. yet R is much slower. To compute the softmax derivative I create matrices for the One Hot Encoded yi and then stack them before subtracting pi-yi. I am sure there is a more elegant and more efficient way to do this, much like Python. Any suggestions? library(ggplot2) library(dplyr) library(RColorBrewer) source("DLfunctions41.R") # Read the spiral dataset Z <- as.matrix(read.csv("spiral.csv",header=FALSE)) Z1=data.frame(Z) #Plot the dataset ggplot(Z1,aes(x=V1,y=V2,col=V3)) +geom_point() + scale_colour_gradientn(colours = brewer.pal(10, "Spectral")) # Setup the data X <- Z[,1:2] y <- Z[,3] X1 <- t(X) Y1 <- t(y) # Initialize number of features, number of hidden units in hidden layer and # number of classes numFeats<-2 # No features numHidden<-100 # No of hidden units numOutput<-3 # No of classes # Initialize model parameters <-initializeModel(numFeats, numHidden,numOutput) W1 <-parameters[['W1']] b1 <-parameters[['b1']] W2 <-parameters[['W2']] b2 <-parameters[['b2']] # Set the learning rate learningRate <- 0.5 # Initialize losses losses <- NULL # Perform gradient descent for(i in 0:9000){ # Forward propagation through hidden layer with Relu units retvals <- layerActivationForward(X1,W1,b1,'relu') A1 <- retvals[['A']] cache1 <- retvals[['cache']] forward_cache1 <- cache1[['forward_cache1']] activation_cache <- cache1[['activation_cache']] # Forward propagation through output layer with Softmax units retvals = layerActivationForward(A1,W2,b2,'softmax') A2 <- retvals[['A']] cache2 <- retvals[['cache']] forward_cache2 <- cache2[['forward_cache1']] activation_cache2 <- cache2[['activation_cache']] # No oftraining examples numTraining <- dim(X)[1] dA <-0 # Select the elements where the y values are 0, 1 or 2 and make a vector a=c(A2[y==0,1],A2[y==1,2],A2[y==2,3]) # Take log correct_probs = -log(a) # Compute loss loss= sum(correct_probs)/numTraining if(i %% 1000 == 0){ sprintf("iteration %d: loss %f",i, loss) print(loss) } # Backward propagation through output layer with Softmax units retvals = layerActivationBackward(dA, cache2, y, activationFunc='softmax') dA1 = retvals[['dA_prev']] dW2= retvals[['dW']] db2= retvals[['db']] # Backward propagation through hidden layer with Relu units retvals = layerActivationBackward(t(dA1), cache1, y, activationFunc='relu') dA0 = retvals[['dA_prev']] dW1= retvals[['dW']] db1= retvals[['db']] # Update parameters W1 <- W1 - learningRate * dW1 b1 <- b1 - learningRate * db1 W2 <- W2 - learningRate * t(dW2) b2 <- b2 - learningRate * t(db2) } ## [1] 1.212487 ## [1] 0.5740867 ## [1] 0.4048824 ## [1] 0.3561941 ## [1] 0.2509576 ## [1] 0.7351063 ## [1] 0.2066114 ## [1] 0.2065875 ## [1] 0.2151943 ## [1] 0.1318807 #Create iterations iterations <- seq(0,10) #df=data.frame(iterations,losses) ggplot(df,aes(x=iterations,y=losses)) + geom_point() + geom_line(color="blue") + ggtitle("Losses vs iterations") + xlab("Iterations") + ylab("Loss") plotDecisionBoundary(Z,W1,b1,W2,b2) Multi-class Confusion Matrix library(caret) library(e1071) # Forward propagation through hidden layer with Relu units retvals <- layerActivationForward(X1,W1,b1,'relu') A1 <- retvals[['A']] # Forward propagation through output layer with Softmax units retvals = layerActivationForward(A1,W2,b2,'softmax') A2 <- retvals[['A']] yhat <- apply(A2, 1,which.max) -1 Confusion Matrix and Statistics Reference Prediction 0 1 2 0 97 0 1 1 2 96 4 2 1 4 95 Overall Statistics Accuracy : 0.96 95% CI : (0.9312, 0.9792) No Information Rate : 0.3333 P-Value [Acc > NIR] : <2e-16 Kappa : 0.94 Mcnemar's Test P-Value : 0.5724 Statistics by Class: Class: 0 Class: 1 Class: 2 Sensitivity 0.9700 0.9600 0.9500 Specificity 0.9950 0.9700 0.9750 Pos Pred Value 0.9898 0.9412 0.9500 Neg Pred Value 0.9851 0.9798 0.9750 Prevalence 0.3333 0.3333 0.3333 Detection Rate 0.3233 0.3200 0.3167 Detection Prevalence 0.3267 0.3400 0.3333 Balanced Accuracy 0.9825 0.9650 0.9625  My book “Practical Machine Learning with R and Python” includes the implementation for many Machine Learning algorithms and associated metrics. Pick up your copy today! 2.3 Multi-class classification with Softmax – Octave code A 2 layer Neural network with the Softmax activation unit in the output layer is constructed in Octave. The same spiral data set is used for Octave also  source("DL41functions.m") # Read the spiral data data=csvread("spiral.csv"); # Setup the data X=data(:,1:2); Y=data(:,3); # Set the number of features, number of hidden units in hidden layer and number of classes numFeats=2; #No features numHidden=100; # No of hidden units numOutput=3; # No of classes # Initialize model [W1 b1 W2 b2] = initializeModel(numFeats,numHidden,numOutput); # Initialize losses losses=[] #Initialize learningRate learningRate=0.5; for k =1:10000 # Forward propagation through hidden layer with Relu units [A1,cache1 activation_cache1]= layerActivationForward(X',W1,b1,activationFunc ='relu'); # Forward propagation through output layer with Softmax units [A2,cache2 activation_cache2] = layerActivationForward(A1,W2,b2,activationFunc='softmax'); # No of training examples numTraining = size(X)(1); # Select rows where Y=0,1,and 2 and concatenate to a long vector a=[A2(Y==0,1) ;A2(Y==1,2) ;A2(Y==2,3)]; #Select the correct column for log prob correct_probs = -log(a); #Compute log loss loss= sum(correct_probs)/numTraining; if(mod(k,1000) == 0) disp(loss); losses=[losses loss]; endif dA=0; # Backward propagation through output layer with Softmax units [dA1 dW2 db2] = layerActivationBackward(dA, cache2, activation_cache2,Y,activationFunc='softmax'); # Backward propagation through hidden layer with Relu units [dA0,dW1,db1] = layerActivationBackward(dA1', cache1, activation_cache1, Y, activationFunc='relu'); #Update parameters W1 += -learningRate * dW1; b1 += -learningRate * db1; W2 += -learningRate * dW2'; b2 += -learningRate * db2'; endfor # Plot Losses vs Iterations iterations=0:1000:9000 plotCostVsIterations(iterations,losses) # Plot the decision boundary plotDecisionBoundary( X,Y,W1,b1,W2,b2) The code for the Python, R and Octave implementations can be downloaded from Github at Deep Learning – Part 4 Conclusion In this post I have implemented a 2 layer Neural Network with the Softmax classifier. In Part 3, I implemented a multi-layer Deep Learning Network. I intend to include the Softmax activation unit into the generalized multi-layer Deep Network along with the other activation units of sigmoid,tanh and relu. Stick around, I’ll be back!! Watch this space! To see all post click Index of posts Deep Learning from first principles in Python, R and Octave – Part 3 “Once upon a time, I, Chuang Tzu, dreamt I was a butterfly, fluttering hither and thither, to all intents and purposes a butterfly. I was conscious only of following my fancies as a butterfly, and was unconscious of my individuality as a man. Suddenly, I awoke, and there I lay, myself again. Now I do not know whether I was then a man dreaming I was a butterfly, or whether I am now a butterfly dreaming that I am a man.” from The Brain: The Story of you – David Eagleman “Thought is a great big vector of neural activity” Prof Geoffrey Hinton Introduction This is the third part in my series on Deep Learning from first principles in Python, R and Octave. In the first part Deep Learning from first principles in Python, R and Octave-Part 1, I implemented logistic regression as a 2 layer neural network. The 2nd part Deep Learning from first principles in Python, R and Octave-Part 2, dealt with the implementation of 3 layer Neural Networks with 1 hidden layer to perform classification tasks, where the 2 classes cannot be separated by a linear boundary. In this third part, I implement a multi-layer, Deep Learning (DL) network of arbitrary depth (any number of hidden layers) and arbitrary height (any number of activation units in each hidden layer). The implementations of these Deep Learning networks, in all the 3 parts, are based on vectorized versions in Python, R and Octave. The implementation in the 3rd part is for a L-layer Deep Netwwork, but without any regularization, early stopping, momentum or learning rate adaptation techniques. However even the barebones multi-layer DL, is a handful and has enough hyperparameters to fine-tune and adjust. The implementation of the vectorized L-layer Deep Learning network in Python, R and Octave were both exhausting, and exacting!! Keeping track of the indices, layer number and matrix dimensions required quite bit of focus. While the implementation was demanding, it was also very exciting to get the code to work. The trick was to be able to shift gears between the slight quirkiness between the languages. Here are some of challenges I faced 1. Python and Octave allow multiple return values to be unpacked in a single statement. With R, unpacking multiple return values from a list, requires the list returned, to be unpacked separately. I did see that there is a package gsubfn, which does this. I hope this feature becomes a base R feature. 2. Python and R allow dissimilar elements to be saved and returned from functions using dictionaries or lists respectively. However there is no real equivalent in Octave. The closest I got to this functionality in Octave, was the ‘cell array’. But the cell array can be accessed only by the index, and not with the key as in a Python dictionary or R list. This makes things just a bit more difficult in Octave. 3. Python and Octave include implicit broadcasting. In R, broadcasting is not implicit, but R has a nifty function, the sweep(), with which we can broadcast either by columns or by rows 4. The closest equivalent of Python’s dictionary, or R’s list, in Octave is the cell array. However I had to manage separate cell arrays for weights and biases and during gradient descent and separate gradients dW and dB 5. In Python the rank-1 numpy arrays can be annoying at times. This issue is not present in R and Octave. Though the number of lines of code for Deep Learning functions in Python, R and Octave are about ~350 apiece, they have been some of the most difficult code I have implemented. The current vectorized implementation supports the relu, sigmoid and tanh activation functions as of now. I will be adding other activation functions like the ‘leaky relu’, ‘softmax’ and others, to the implementation in the weeks to come. While testing with different hyper-parameters namely i) the number of hidden layers, ii) the number of activation units in each layer, iii) the activation function and iv) the number iterations, I found the L-layer Deep Learning Network to be very sensitive to these hyper-parameters. It is not easy to tune the parameters. Adding more hidden layers, or more units per layer, does not help and mostly results in gradient descent getting stuck in some local minima. It does take a fair amount of trial and error and very close observation on how the DL network performs for logical changes. We then can zero in on the most the optimal solution. Feel free to download/fork my code from Github DeepLearning-Part 3 and play around with the hyper-parameters for your own problems. Derivation of a Multi Layer Deep Learning Network Lets take a simple 3 layer Neural network with 3 hidden layers and an output layer In the forward propagation cycle the equations are $Z_{1} = W_{1}A_{0} +b_{1}$ and $A_{1} = g(Z_{1})$ $Z_{2} = W_{2}A_{1} +b_{2}$ and $A_{2} = g(Z_{2})$ $Z_{3} = W_{3}A_{2} +b_{3}$ and $A_{3} = g(Z_{3})$ The loss function is given by $L = -(ylogA3 + (1-y)log(1-A3))$ and $dL/dA3 = -(Y/A_{3} + (1-Y)/(1-A_{3}))$ For a binary classification the output activation function is the sigmoid function given by $A_{3} = 1/(1+ e^{-Z3})$. It can be shown that $dA_{3}/dZ_{3} = A_{3}(1-A_3)$ see equation 2 in Part 1 $\partial L/\partial Z_{3} = \partial L/\partial A_{3}* \partial A_{3}/\partial Z_{3} = A3-Y$ see equation (f) in Part 1 and since $\partial L/\partial A_{2} = \partial L/\partial Z_{3} * \partial Z_{3}/\partial A_{2} = (A_{3} -Y) * W_{3}$ because $\partial Z_{3}/\partial A_{2} = W_{3}$ -(1a) and $\partial L/\partial Z_{2} =\partial L/\partial A_{2} * \partial A_{2}/\partial Z_{2} = (A_{3} -Y) * W_{3} *g'(Z_{2})$ -(1b) $\partial L/\partial W_{2} = \partial L/\partial Z_{2} * A_{1}$ -(1c) since $\partial Z_{2}/\partial W_{2} = A_{1}$ and $\partial L/\partial b_{2} = \partial L/\partial Z_{2}$ -(1d) because $\partial Z_{2}/\partial b_{2} =1$ Also $\partial L/\partial A_{1} =\partial L/\partial Z_{2} * \partial Z_{2}/\partial A_{1} = \partial L/\partial Z_{2} * W_{2}$ – (2a) $\partial L/\partial Z_{1} =\partial L/\partial A_{1} * \partial A_{1}/\partial Z_{1} = \partial L/\partial A_{1} * W_{3} *g'(Z_{2})$ – (2b) $\partial L/\partial W_{1} = \partial L/\partial Z_{1} * A_{0}$ – (2c) $\partial L/\partial b_{1} = \partial L/\partial Z_{1}$ – (2d) Inspecting the above equations (1a – 1d & 2a-2d), our ‘Uber deep, bottomless’ brain can easily discern the pattern in these equations. The equation for any layer ‘l’ is of the form $Z_{l} = W_{l}A_{l-1} +b_{l}$ and $A_{l} = g(Z_{l})$ The equation for the backward propagation have the general form $\partial L/\partial A_{l} = \partial L/\partial Z_{l+1} * W^{l+1}$ $\partial L/\partial Z_{l}=\partial L/\partial A_{l} *g'(Z_{l})$ $\partial L/\partial W_{l} =\partial L/\partial Z_{l} *A^{l-1}$ $\partial L/\partial b_{l} =\partial L/\partial Z_{l}$ Some other important results The derivatives of the activation functions in the implemented Deep Learning network g(z) = sigmoid(z) = $1/(1+e^{-z})$ = a g’(z) = a(1-a) – See Part 1 g(z) = tanh(z) = a g’(z) = $1 - a^{2}$ g(z) = relu(z) = z when z>0 and 0 when z 0 and 0 when z <= 0 While it appears that there is a discontinuity for the derivative at 0 the small value at the discontinuity does not present a problem The implementation of the multi layer vectorized Deep Learning Network for Python, R and Octave is included below. For all these implementations, initially I create the size and configuration of the the Deep Learning network with the layer dimennsions So for example layersDimension Vector ‘V’ of length L indicating ‘L’ layers where V (in Python)= $[v_{0}, v_{1}, v_{2}$, … $v_{L-1}]$ V (in R)= $c(v_{1}, v_{2}, v_{3}$ , … $v_{L})$ V (in Octave)= [ $v_{1} v_{2} v_{3}$$v_{L}]$ In all of these implementations the first element is the number of input features to the Deep Learning network and the last element is always a ‘sigmoid’ activation function since all the problems deal with binary classification. The number of elements between the first and the last element are the number of hidden layers and the magnitude of each $v_{i}$ is the number of activation units in each hidden layer, which is specified while actually executing the Deep Learning network using the function L_Layer_DeepModel(), in all the implementations Python, R and Octave 1a. Classification with Multi layer Deep Learning Network – Relu activation(Python) In the code below a 4 layer Neural Network is trained to generate a non-linear boundary between the classes. In the code below the ‘Relu’ Activation function is used. The number of activation units in each layer is 9. The cost vs iterations is plotted in addition to the decision boundary. Further the accuracy, precision, recall and F1 score are also computed import os import numpy as np import matplotlib.pyplot as plt import matplotlib.colors import sklearn.linear_model from sklearn.model_selection import train_test_split from sklearn.datasets import make_classification, make_blobs from matplotlib.colors import ListedColormap import sklearn import sklearn.datasets #from DLfunctions import plot_decision_boundary execfile("./DLfunctions34.py") # os.chdir("C:\\software\\DeepLearning-Posts\\part3") # Create clusters of 2 classes X1, Y1 = make_blobs(n_samples = 400, n_features = 2, centers = 9, cluster_std = 1.3, random_state = 4) #Create 2 classes Y1=Y1.reshape(400,1) Y1 = Y1 % 2 X2=X1.T Y2=Y1.T # Set the dimensions of DL Network # Below we have # 2 - 2 input features # 9,9 - 2 hidden layers with 9 activation units per layer and # 1 - 1 sigmoid activation unit in the output layer as this is a binary classification # The activation in the hidden layer is the 'relu' specified in L_Layer_DeepModel layersDimensions = [2, 9, 9,1] # 4-layer model parameters = L_Layer_DeepModel(X2, Y2, layersDimensions,hiddenActivationFunc='relu', learning_rate = 0.3,num_iterations = 2500, fig="fig1.png") #Plot the decision boundary plot_decision_boundary(lambda x: predict(parameters, x.T), X2,Y2,str(0.3),"fig2.png") # Compute the confusion matrix yhat = predict(parameters,X2) from sklearn.metrics import confusion_matrix a=confusion_matrix(Y2.T,yhat.T) from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score print('Accuracy: {:.2f}'.format(accuracy_score(Y2.T, yhat.T))) print('Precision: {:.2f}'.format(precision_score(Y2.T, yhat.T))) print('Recall: {:.2f}'.format(recall_score(Y2.T, yhat.T))) print('F1: {:.2f}'.format(f1_score(Y2.T, yhat.T))) ## Accuracy: 0.90 ## Precision: 0.91 ## Recall: 0.87 ## F1: 0.89 For more details on metrics like Accuracy, Recall, Precision etc. used in classification take a look at my post Practical Machine Learning with R and Python – Part 2. More details about these and other metrics besides implementation of the most common machine learning algorithms are available in my book My book ‘Practical Machine Learning with R and Python’ on Amazon 1b. Classification with Multi layer Deep Learning Network – Relu activation(R) In the code below, binary classification is performed on the same data set as above using the Relu activation function. The DL network is same as above library(ggplot2) # Read the data z <- as.matrix(read.csv("data.csv",header=FALSE)) x <- z[,1:2] y <- z[,3] X1 <- t(x) Y1 <- t(y) # Set the dimensions of the Deep Learning network # No of input features =2, 2 hidden layers with 9 activation units and 1 output layer layersDimensions = c(2, 9, 9,1) # Execute the Deep Learning Neural Network retvals = L_Layer_DeepModel(X1, Y1, layersDimensions, hiddenActivationFunc='relu', learningRate = 0.3, numIterations = 5000, print_cost = True) library(ggplot2) source("DLfunctions33.R") # Get the computed costs costs <- retvals[['costs']] # Create a sequence of iterations numIterations=5000 iterations <- seq(0,numIterations,by=1000) df <-data.frame(iterations,costs) # Plot the Costs vs number of iterations ggplot(df,aes(x=iterations,y=costs)) + geom_point() +geom_line(color="blue") + xlab('No of iterations') + ylab('Cost') + ggtitle("Cost vs No of iterations") # Plot the decision boundary plotDecisionBoundary(z,retvals,hiddenActivationFunc="relu",0.3) library(caret) # Predict the output for the data values yhat <-predict(retvals$parameters,X1,hiddenActivationFunc="relu")
yhat[yhat==FALSE]=0
yhat[yhat==TRUE]=1
# Compute the confusion matrix
confusionMatrix(yhat,Y1)
## Confusion Matrix and Statistics
##
##           Reference
## Prediction   0   1
##          0 201  10
##          1  21 168
##
##                Accuracy : 0.9225
##                  95% CI : (0.8918, 0.9467)
##     No Information Rate : 0.555
##     P-Value [Acc > NIR] : < 2e-16
##
##                   Kappa : 0.8441
##  Mcnemar's Test P-Value : 0.07249
##
##             Sensitivity : 0.9054
##             Specificity : 0.9438
##          Pos Pred Value : 0.9526
##          Neg Pred Value : 0.8889
##              Prevalence : 0.5550
##          Detection Rate : 0.5025
##    Detection Prevalence : 0.5275
##       Balanced Accuracy : 0.9246
##
##        'Positive' Class : 0
## 

1c. Classification with Multi layer Deep Learning Network – Relu activation(Octave)

Included below is the code for performing classification. Incidentally Octave does not seem to have implemented the confusion matrix,  but confusionmat is available in Matlab. # Read the data data=csvread("data.csv"); X=data(:,1:2); Y=data(:,3); # Set layer dimensions layersDimensions = [2 9 7 1] #tanh=-0.5(ok), #relu=0.1 best! # Execute Deep Network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', learningRate = 0.1, numIterations = 10000); plotCostVsIterations(10000,costs); plotDecisionBoundary(data,weights, biases,hiddenActivationFunc="tanh") 

2a. Classification with Multi layer Deep Learning Network – Tanh activation(Python)

Below the Tanh activation function is used to perform the same classification. I found the Tanh activation required a simpler Neural Network of 3 layers.

# Tanh activation
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import sklearn.linear_model

from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification, make_blobs
from matplotlib.colors import ListedColormap
import sklearn
import sklearn.datasets

#from DLfunctions import plot_decision_boundary
os.chdir("C:\\software\\DeepLearning-Posts\\part3")
execfile("./DLfunctions34.py")
# Create the dataset
X1, Y1 = make_blobs(n_samples = 400, n_features = 2, centers = 9,
cluster_std = 1.3, random_state = 4)
#Create 2 classes
Y1=Y1.reshape(400,1)
Y1 = Y1 % 2
X2=X1.T
Y2=Y1.T
# Set the dimensions of the Neural Network
layersDimensions = [2, 4, 1] #  3-layer model
# Compute the DL network
parameters = L_Layer_DeepModel(X2, Y2, layersDimensions, hiddenActivationFunc='tanh', learning_rate = .5,num_iterations = 2500,fig="fig3.png")
#Plot the decision boundary
plot_decision_boundary(lambda x: predict(parameters, x.T), X2,Y2,str(0.5),"fig4.png")


2b. Classification with Multi layer Deep Learning Network – Tanh activation(R)

R performs better with a Tanh activation than the Relu as can be seen below

 #Set the dimensions of the Neural Network
layersDimensions = c(2, 9, 9,1)
library(ggplot2)
x <- z[,1:2]
y <- z[,3]
X1 <- t(x)
Y1 <- t(y)
# Execute the Deep Model
retvals = L_Layer_DeepModel(X1, Y1, layersDimensions,
hiddenActivationFunc='tanh',
learningRate = 0.3,
numIterations = 5000,
print_cost = True)
# Get the costs
costs <- retvals[['costs']]
iterations <- seq(0,numIterations,by=1000)
df <-data.frame(iterations,costs)
# Plot Cost vs number of iterations
ggplot(df,aes(x=iterations,y=costs)) + geom_point() +geom_line(color="blue") +
xlab('No of iterations') + ylab('Cost') + ggtitle("Cost vs No of iterations")

#Plot the decision boundary
plotDecisionBoundary(z,retvals,hiddenActivationFunc="tanh",0.3)

2c. Classification with Multi layer Deep Learning Network – Tanh activation(Octave)

The code below uses the   Tanh activation in the hidden layers for Octave # Read the data data=csvread("data.csv"); X=data(:,1:2); Y=data(:,3); # Set layer dimensions layersDimensions = [2 9 7 1] #tanh=-0.5(ok), #relu=0.1 best! # Execute Deep Network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='tanh', learningRate = 0.1, numIterations = 10000); plotCostVsIterations(10000,costs); plotDecisionBoundary(data,weights, biases,hiddenActivationFunc="tanh") 

3. Bernoulli’s Lemniscate

To make things  more interesting, I create a 2D figure of the Bernoulli’s lemniscate to perform non-linear classification. The Lemniscate is given by the equation
$(x^{2} + y^{2})^{2}$ = $2a^{2}*(x^{2}-y^{2})$

3a. Classifying a lemniscate with Deep Learning Network – Relu activation(Python)

import os
import numpy as np
import matplotlib.pyplot as plt
os.chdir("C:\\software\\DeepLearning-Posts\\part3")
execfile("./DLfunctions33.py")
x1=np.random.uniform(0,10,2000).reshape(2000,1)
x2=np.random.uniform(0,10,2000).reshape(2000,1)

X=np.append(x1,x2,axis=1)
X.shape

# Create a subset of values where squared is <0,4. Perform ravel() to flatten this vector
# Create the equation
# (x^{2} + y^{2})^2 - 2a^2*(x^{2}-y^{2}) <= 0
a=np.power(np.power(X[:,0]-5,2) + np.power(X[:,1]-5,2),2)
b=np.power(X[:,0]-5,2) - np.power(X[:,1]-5,2)
c= a - (b*np.power(4,2)) <=0
Y=c.reshape(2000,1)
# Create a scatter plot of the lemniscate
plt.scatter(X[:,0], X[:,1], c=Y, marker= 'o', s=15,cmap="viridis")
Z=np.append(X,Y,axis=1)
plt.savefig("fig50.png",bbox_inches='tight')
plt.clf()

# Set the data for classification
X2=X.T
Y2=Y.T
# These settings work the best
# Set the Deep Learning layer dimensions for a Relu activation
layersDimensions = [2,7,4,1]
#Execute the DL network
parameters = L_Layer_DeepModel(X2, Y2, layersDimensions, hiddenActivationFunc='relu', learning_rate = 0.5,num_iterations = 10000, fig="fig5.png")
#Plot the decision boundary
plot_decision_boundary(lambda x: predict(parameters, x.T), X2, Y2,str(2.2),"fig6.png")

# Compute the Confusion matrix
yhat = predict(parameters,X2)
from sklearn.metrics import confusion_matrix
a=confusion_matrix(Y2.T,yhat.T)
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
print('Accuracy: {:.2f}'.format(accuracy_score(Y2.T, yhat.T)))
print('Precision: {:.2f}'.format(precision_score(Y2.T, yhat.T)))
print('Recall: {:.2f}'.format(recall_score(Y2.T, yhat.T)))
print('F1: {:.2f}'.format(f1_score(Y2.T, yhat.T)))
## Accuracy: 0.93
## Precision: 0.77
## Recall: 0.76
## F1: 0.76

We could get better performance by tuning further. Do play around if you fork the code.
Note:: The lemniscate data is saved as a CSV and then read in R and also in Octave. I do this instead of recreating the lemniscate shape

3b. Classifying a lemniscate with Deep Learning Network – Relu activation(R code)

The R decision boundary for the Bernoulli’s lemniscate is shown below

Z <- as.matrix(read.csv("lemniscate.csv",header=FALSE))
Z1=data.frame(Z)
# Create a scatter plot of the lemniscate
ggplot(Z1,aes(x=V1,y=V2,col=V3)) +geom_point()
#Set the data for the DL network
X=Z[,1:2]
Y=Z[,3]

X1=t(X)
Y1=t(Y)

# Set the layer dimensions for the tanh activation function
layersDimensions = c(2,5,4,1)
# Execute the Deep Learning network with Tanh activation
retvals = L_Layer_DeepModel(X1, Y1, layersDimensions,
hiddenActivationFunc='tanh',
learningRate = 0.3,
numIterations = 20000, print_cost = True)
# Plot cost vs iteration
costs <- retvals[['costs']]
numIterations = 20000
iterations <- seq(0,numIterations,by=1000)
df <-data.frame(iterations,costs)
ggplot(df,aes(x=iterations,y=costs)) + geom_point() +geom_line(color="blue") +
xlab('No of iterations') + ylab('Cost') + ggtitle("Cost vs No of iterations")

#Plot the decision boundary
plotDecisionBoundary(Z,retvals,hiddenActivationFunc="tanh",0.3)

3c. Classifying a lemniscate with Deep Learning Network – Relu activation(Octave code)

Octave is used to generate the non-linear lemniscate boundary.
 # Read the data data=csvread("lemniscate.csv"); X=data(:,1:2); Y=data(:,3); # Set the dimensions of the layers layersDimensions = [2 9 7 1] # Compute the DL network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', learningRate = 0.20, numIterations = 10000); plotCostVsIterations(10000,costs); plotDecisionBoundary(data,weights, biases,hiddenActivationFunc="relu") 

4a. Binary Classification using MNIST – Python code

Finally I perform a simple classification using the MNIST handwritten digits, which according to Prof Geoffrey Hinton is “the Drosophila of Deep Learning”.

The Python code for reading the MNIST data is taken from Alex Kesling’s github link MNIST.

In the Python code below, I perform a simple binary classification between the handwritten digit ‘5’ and ‘not 5’ which is all other digits. I will perform the proper classification of all digits using the  Softmax classifier some time later.

import os
import numpy as np
import matplotlib.pyplot as plt
os.chdir("C:\\software\\DeepLearning-Posts\\part3")
execfile("./DLfunctions34.py")
lbls=[]
pxls=[]
print(len(training))

# Select the first 10000 training data and the labels
for i in range(10000):
l,p=training[i]
lbls.append(l)
pxls.append(p)
labels= np.array(lbls)
pixels=np.array(pxls)

#  Sey y=1  when labels == 5 and 0 otherwise
y=(labels==5).reshape(-1,1)
X=pixels.reshape(pixels.shape[0],-1)

# Create the necessary feature and target variable
X1=X.T
Y1=y.T

# Create the layer dimensions. The number of features are 28 x 28 = 784 since the 28 x 28
# pixels is flattened to single vector of length 784.
layersDimensions=[784, 15,9,7,1] # Works very well
parameters = L_Layer_DeepModel(X1, Y1, layersDimensions, hiddenActivationFunc='relu', learning_rate = 0.1,num_iterations = 1000, fig="fig7.png")

# Test data
lbls1=[]
pxls1=[]
for i in range(800):
l,p=test[i]
lbls1.append(l)
pxls1.append(p)

testLabels=np.array(lbls1)
testData=np.array(pxls1)

ytest=(testLabels==5).reshape(-1,1)
Xtest=testData.reshape(testData.shape[0],-1)
Xtest1=Xtest.T
Ytest1=ytest.T

yhat = predict(parameters,Xtest1)
from sklearn.metrics import confusion_matrix
a=confusion_matrix(Ytest1.T,yhat.T)
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score
print('Accuracy: {:.2f}'.format(accuracy_score(Ytest1.T, yhat.T)))
print('Precision: {:.2f}'.format(precision_score(Ytest1.T, yhat.T)))
print('Recall: {:.2f}'.format(recall_score(Ytest1.T, yhat.T)))
print('F1: {:.2f}'.format(f1_score(Ytest1.T, yhat.T)))

probs=predict_proba(parameters,Xtest1)
from sklearn.metrics import precision_recall_curve

precision, recall, thresholds = precision_recall_curve(Ytest1.T, probs.T)
closest_zero = np.argmin(np.abs(thresholds))
closest_zero_p = precision[closest_zero]
closest_zero_r = recall[closest_zero]
plt.xlim([0.0, 1.01])
plt.ylim([0.0, 1.01])
plt.plot(precision, recall, label='Precision-Recall Curve')
plt.plot(closest_zero_p, closest_zero_r, 'o', markersize = 12, fillstyle = 'none', c='r', mew=3)
plt.xlabel('Precision', fontsize=16)
plt.ylabel('Recall', fontsize=16)
plt.savefig("fig8.png",bbox_inches='tight')


## Accuracy: 0.99
## Precision: 0.96
## Recall: 0.89
## F1: 0.92

In addition to plotting the Cost vs Iterations, I also plot the Precision-Recall curve to show how the Precision and Recall, which are complementary to each other vary with respect to the other. To know more about Precision-Recall, please check my post Practical Machine Learning with R and Python – Part 4. You could also check out my book My book ‘Practical Machine Learning with R and Python’ on Amazon for details on the key metrics and algorithms for classification and regression problems. A physical copy of the book is much better than scrolling down a webpage. Personally, I tend to use my own book quite frequently to refer to R, Python constructs,  subsetting, machine Learning function calls and the necessary parameters etc. It is useless to commit any of this to memory, and a physical copy of a book is much easier to thumb through for the relevant code snippet. Pick up your copy today!

4b. Binary Classification using MNIST – R code

In the R code below the same binary classification of the digit ‘5’ and the ‘not 5’ is performed. The code to read and display the MNIST data is taken from Brendan O’ Connor’s github link at MNIST

source("mnist.R")
#show_digit(train$x[2,] layersDimensions=c(784, 7,7,3,1) # Works at 1500 x <- t(train$x)
# Choose only 5000 training data
x2 <- x[,1:5000]
y <-train$y # Set labels for all digits that are 'not 5' to 0 y[y!=5] <- 0 # Set labels of digit 5 as 1 y[y==5] <- 1 # Set the data y1 <- as.matrix(y) y2 <- t(y1) # Choose the 1st 5000 data y3 <- y2[,1:5000] #Execute the Deep Learning Model retvals = L_Layer_DeepModel(x2, y3, layersDimensions, hiddenActivationFunc='tanh', learningRate = 0.3, numIterations = 3000, print_cost = True) # Plot cost vs iteration costs <- retvals[['costs']] numIterations = 3000 iterations <- seq(0,numIterations,by=1000) df <-data.frame(iterations,costs) ggplot(df,aes(x=iterations,y=costs)) + geom_point() +geom_line(color="blue") + xlab('No of iterations') + ylab('Cost') + ggtitle("Cost vs No of iterations") # Compute probability scores scores <- computeScores(retvals$parameters, x2,hiddenActivationFunc='relu')
a=y3==1
b=y3==0

# Compute probabilities of class 0 and class 1
class1=scores[a]
class0=scores[b]

# Plot ROC curve
pr <-pr.curve(scores.class0=class1,
scores.class1=class0,
curve=T)

plot(pr)

The AUC curve hugs the top left corner and hence the performance of the classifier is quite good.

4c. Binary Classification using MNIST – Octave code

This code to load MNIST data was taken from Daniel E blog.
Precision recall curves are available in Matlab but are yet to be implemented in Octave’s statistics package.
 load('./mnist/mnist.txt.gz'); % load the dataset # Subset the 'not 5' digits a=(trainY != 5); # Subset '5' b=(trainY == 5); #make a copy of trainY #Set 'not 5' as 0 and '5' as 1 y=trainY; y(a)=0; y(b)=1; X=trainX(1:5000,:); Y=y(1:5000); # Set the dimensions of layer layersDimensions=[784, 7,7,3,1]; # Compute the DL network [weights biases costs]=L_Layer_DeepModel(X', Y', layersDimensions, hiddenActivationFunc='relu', learningRate = 0.1, numIterations = 5000); 

Conclusion

It was quite a challenge coding a Deep Learning Network in Python, R and Octave. The Deep Learning network implementation, in this post,is the base Deep Learning network, without any of the regularization methods included. Here are some key learning that I got while playing with different multi-layer networks on different problems

a. Deep Learning Networks come with many levers, the hyper-parameters,
– learning rate
– activation unit
– number of hidden layers
– number of units per hidden layer
– number of iterations while performing gradient descent
b. Deep Networks are very sensitive. A change in any of the hyper-parameter makes it perform very differently
c. Initially I thought adding more hidden layers, or more units per hidden layer will make the DL network better at learning. On the contrary, there is a performance degradation after the optimal DL configuration
d. At a sub-optimal number of hidden layers or number of hidden units, gradient descent seems to get stuck at a local minima
e. There were occasions when the cost came down, only to increase slowly as the number of iterations were increased. Probably early stopping would have helped.
f. I also did come across situations of ‘exploding/vanishing gradient’, cost went to Inf/-Inf. Here I would think inclusion of ‘momentum method’ would have helped

I intend to add the additional hyper-parameters of L1, L2 regularization, momentum method, early stopping etc. into the code in my future posts.
Feel free to fork/clone the code from Github Deep Learning – Part 3, and take the DL network apart and play around with it.

I will be continuing this series with more hyper-parameters to handle vanishing and exploding gradients, early stopping and regularization in the weeks to come. I also intend to add some more activation functions to this basic Multi-Layer Network.
Hang around, there are more exciting things to come.

Watch this space!

To see all posts see Index of posts

Deep Learning from first principles in Python, R and Octave – Part 2

“What does the world outside your head really ‘look’ like? Not only is there no color, there’s also no sound: the compression and expansion of air is picked up by the ears, and turned into electrical signals. The brain then presents these signals to us as mellifluous tones and swishes and clatters and jangles. Reality is also odorless: there’s no such thing as smell outside our brains. Molecules floating through the air bind to receptors in our nose and are interpreted as different smells by our brain. The real world is not full of rich sensory events; instead, our brains light up the world with their own sensuality.”
The Brain: The Story of You” by David Eagleman

The world is Maya, illusory. The ultimate reality, the Brahman, is all-pervading and all-permeating, which is colourless, odourless, tasteless, nameless and formless

1. Introduction

This post is a follow-up post to my earlier post Deep Learning from first principles in Python, R and Octave-Part 1. In the first part, I implemented Logistic Regression, in vectorized Python,R and Octave, with a wannabe Neural Network (a Neural Network with no hidden layers). In this second part, I implement a regular, but somewhat primitive Neural Network (a Neural Network with just 1 hidden layer). The 2nd part implements classification of manually created datasets, where the different clusters of the 2 classes are not linearly separable.

Neural Network perform really well in learning all sorts of non-linear boundaries between classes. Initially logistic regression is used perform the classification and the decision boundary is plotted. Vanilla logistic regression performs quite poorly. Using SVMs with a radial basis kernel would have performed much better in creating non-linear boundaries. To see R and Python implementations of SVMs take a look at my post Practical Machine Learning with R and Python – Part 4.

You could also check out my book on Amazon Practical Machine Learning with R and Python – Machine Learning in Stereo,  in which I implement several Machine Learning algorithms on regression and classification, along with other necessary metrics that are used in Machine Learning.

You can clone and fork this R Markdown file along with the vectorized implementations of the 3 layer Neural Network for Python, R and Octave from Github DeepLearning-Part2

2. The 3 layer Neural Network

A simple representation of a 3 layer Neural Network (NN) with 1 hidden layer is shown below.

In the above Neural Network, there are 2 input features at the input layer, 3 hidden units at the hidden layer and 1 output layer as it deals with binary classification. The activation unit at the hidden layer can be a tanh, sigmoid, relu etc. At the output layer the activation is a sigmoid to handle binary classification

# Superscript indicates layer 1
$z_{11} = w_{11}^{1}x_{1} + w_{21}^{1}x_{2} + b_{1}$
$z_{12} = w_{12}^{1}x_{1} + w_{22}^{1}x_{2} + b_{1}$
$z_{13} = w_{13}^{1}x_{1} + w_{23}^{1}x_{2} + b_{1}$

Also $a_{11} = tanh(z_{11})$
$a_{12} = tanh(z_{12})$
$a_{13} = tanh(z_{13})$

# Superscript indicates layer 2
$z_{21} = w_{11}^{2}a_{11} + w_{21}^{2}a_{12} + w_{31}^{2}a_{13} + b_{2}$
$a_{21} = sigmoid(z21)$

Hence
$Z1= \begin{pmatrix} z11\\ z12\\ z13 \end{pmatrix} =\begin{pmatrix} w_{11}^{1} & w_{21}^{1} \\ w_{12}^{1} & w_{22}^{1} \\ w_{13}^{1} & w_{23}^{1} \end{pmatrix} * \begin{pmatrix} x1\\ x2 \end{pmatrix} + b_{1}$
And
$A1= \begin{pmatrix} a11\\ a12\\ a13 \end{pmatrix} = \begin{pmatrix} tanh(z11)\\ tanh(z12)\\ tanh(z13) \end{pmatrix}$

Similarly
$Z2= z_{21} = \begin{pmatrix} w_{11}^{2} & w_{21}^{2} & w_{31}^{2} \end{pmatrix} *\begin{pmatrix} z_{11}\\ z_{12}\\ z_{13} \end{pmatrix} +b_{2}$
and $A2 = a_{21} = sigmoid(z_{21})$

These equations can be written as
$Z1 = W1 * X + b1$
$A1 = tanh(Z1)$
$Z2 = W2 * A1 + b2$
$A2 = sigmoid(Z2)$

I) Some important results (a memory refresher!)
$d/dx(e^{x}) = e^{x}$ and $d/dx(e^{-x}) = -e^{-x}$ -(a) and
$sinhx = (e^{x} - e^{-x})/2$ and $coshx = (e^{x} + e^{-x})/2$
Using (a) we can shown that $d/dx(sinhx) = coshx$ and $d/dx(coshx) = sinhx$ (b)
Now $d/dx(f(x)/g(x)) = (g(x)*d/dx(f(x)) - f(x)*d/dx(g(x)))/g(x)^{2}$ -(c)

Since $tanhx =z= sinhx/coshx$ and using (b) we get
$tanhx = (coshx*d/dx(sinhx) - coshx*d/dx(sinhx))/(1-sinhx^{2})$
Using the values of the derivatives of sinhx and coshx from (b) above we get
$d/dx(tanhx) = (coshx^{2} - sinhx{2})/coshx{2} = 1 - tanhx^{2}$
Since $tanhx =z$
$d/dx(tanhx) = 1 - tanhx^{2}= 1 - z^{2}$ -(d)

II) Derivatives
$L=-(Ylog(A2) + (1-Y)log(1-A2))$
$dL/dA2 = -(Y/A2 + (1-Y)/(1-A2))$
Since $A2 = sigmoid(Z2)$ therefore $dA2/dZ2 = A2(1-A2)$ see Part1
$Z2 = W2A1 +b2$
$dZ2/dW2 = A1$
$dZ2/db2 = 1$
$A1 = tanh(Z1)$ and $dA1/dZ1 = 1 - A1^{2}$
$Z1 = W1X + b1$
$dZ1/dW1 = X$
$dZ1/db1 = 1$

III) Back propagation
Using the derivatives from II) we can derive the following results using Chain Rule
$\partial L/\partial Z2 = \partial L/\partial A2 * \partial A2/\partial Z2$
$= -(Y/A2 + (1-Y)/(1-A2)) * A2(1-A2) = A2 - Y$
$\partial L/\partial W2 = \partial L/\partial A2 * \partial A2/\partial Z2 * \partial Z2/\partial W2$
$= (A2-Y) *A1$ -(A)
$\partial L/\partial b2 = \partial L/\partial A2 * \partial A2/\partial Z2 * \partial Z2/\partial b2 = (A2-Y)$ -(B)

$\partial L/\partial Z1 = \partial L/\partial A2 * \partial A2/\partial Z2 * \partial Z2/\partial A1 *\partial A1/\partial Z1 = (A2-Y) * W2 * (1-A1^{2})$
$\partial L/\partial W1 = \partial L/\partial A2 * \partial A2/\partial Z2 * \partial Z2/\partial A1 *\partial A1/\partial Z1 *\partial Z1/\partial W1$
$=(A2-Y) * W2 * (1-A1^{2}) * X$ -(C)
$\partial L/\partial b1 = \partial L/\partial A2 * \partial A2/\partial Z2 * \partial Z2/\partial A1 *dA1/dZ1 *dZ1/db1$
$= (A2-Y) * W2 * (1-A1^{2})$ -(D)

The key computations in the backward cycle are
$W1 = W1-learningRate * \partial L/\partial W1$ – From (C)
$b1 = b1-learningRate * \partial L/\partial b1$ – From (D)
$W2 = W2-learningRate * \partial L/\partial W2$ – From (A)
$b2 = b2-learningRate * \partial L/\partial b2$ – From (B)

The weights and biases (W1,b1,W2,b2) are updated for each iteration thus minimizing the loss/cost.

These derivations can be represented pictorially using the computation graph (from the book Deep Learning by Ian Goodfellow, Joshua Bengio and Aaron Courville)

3. Manually create a data set that is not lineary separable

Initially I create a dataset with 2 classes which has around 9 clusters that cannot be separated by linear boundaries. Note: This data set is saved as data.csv and is used for the R and Octave Neural networks to see how they perform on the same dataset.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import sklearn.linear_model

from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification, make_blobs
from matplotlib.colors import ListedColormap
import sklearn
import sklearn.datasets

colors=['black','gold']
cmap = matplotlib.colors.ListedColormap(colors)
X, y = make_blobs(n_samples = 400, n_features = 2, centers = 7,
cluster_std = 1.3, random_state = 4)
#Create 2 classes
y=y.reshape(400,1)
y = y % 2
#Plot the figure
plt.figure()
plt.title('Non-linearly separable classes')
plt.scatter(X[:,0], X[:,1], c=y,
marker= 'o', s=50,cmap=cmap)
plt.savefig('fig1.png', bbox_inches='tight')

4. Logistic Regression

On the above created dataset, classification with logistic regression is performed, and the decision boundary is plotted. It can be seen that logistic regression performs quite poorly

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import sklearn.linear_model

from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification, make_blobs
from matplotlib.colors import ListedColormap
import sklearn
import sklearn.datasets

#from DLfunctions import plot_decision_boundary
execfile("./DLfunctions.py") # Since import does not work in Rmd!!!

colors=['black','gold']
cmap = matplotlib.colors.ListedColormap(colors)
X, y = make_blobs(n_samples = 400, n_features = 2, centers = 7,
cluster_std = 1.3, random_state = 4)
#Create 2 classes
y=y.reshape(400,1)
y = y % 2

# Train the logistic regression classifier
clf = sklearn.linear_model.LogisticRegressionCV();
clf.fit(X, y);

# Plot the decision boundary for logistic regression
plot_decision_boundary_n(lambda x: clf.predict(x), X.T, y.T,"fig2.png")


5. The 3 layer Neural Network in Python (vectorized)

The vectorized implementation is included below. Note that in the case of Python a learning rate of 0.5 and 3 hidden units performs very well.

## Random data set with 9 clusters
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn.linear_model
import pandas as pd

from sklearn.datasets import make_classification, make_blobs
execfile("./DLfunctions.py") # Since import does not work in Rmd!!!

X1, Y1 = make_blobs(n_samples = 400, n_features = 2, centers = 9,
cluster_std = 1.3, random_state = 4)
#Create 2 classes
Y1=Y1.reshape(400,1)
Y1 = Y1 % 2
X2=X1.T
Y2=Y1.T

parameters,costs = computeNN(X2, Y2, numHidden = 4, learningRate=0.5, numIterations = 10000)
plot_decision_boundary(lambda x: predict(parameters, x.T), X2, Y2,str(4),str(0.5),"fig3.png")
## Cost after iteration 0: 0.692669
## Cost after iteration 1000: 0.246650
## Cost after iteration 2000: 0.227801
## Cost after iteration 3000: 0.226809
## Cost after iteration 4000: 0.226518
## Cost after iteration 5000: 0.226331
## Cost after iteration 6000: 0.226194
## Cost after iteration 7000: 0.226085
## Cost after iteration 8000: 0.225994
## Cost after iteration 9000: 0.225915

6. The 3 layer Neural Network in R (vectorized)

For this the dataset created by Python is saved  to see how R performs on the same dataset. The vectorized implementation of a Neural Network was just a little more interesting as R does not have a similar package like ‘numpy’. While numpy handles broadcasting implicitly, in R I had to use the ‘sweep’ command to broadcast. The implementaion is included below. Note that since the initialization with random weights is slightly different, R performs best with a learning rate of 0.1 and with 6 hidden units

source("DLfunctions2_1.R")
z <- as.matrix(read.csv("data.csv",header=FALSE)) #
x <- z[,1:2]
y <- z[,3]
x1 <- t(x)
y1 <- t(y)
nn <-computeNN(x1, y1, 6, learningRate=0.1,numIterations=10000) # Good
## [1] 0.7075341
## [1] 0.2606695
## [1] 0.2198039
## [1] 0.2091238
## [1] 0.211146
## [1] 0.2108461
## [1] 0.2105351
## [1] 0.210211
## [1] 0.2099104
## [1] 0.2096437
## [1] 0.209409
plotDecisionBoundary(z,nn,6,0.1)

7.  The 3 layer Neural Network in Octave (vectorized)

This uses the same dataset that was generated using Python code.
source("DL-function2.m") data=csvread("data.csv"); X=data(:,1:2); Y=data(:,3); # Make sure that the model parameters are correct. Take the transpose of X & Y
#Perform gradient descent [W1,b1,W2,b2,costs]= computeNN(X', Y',4, learningRate=0.5, numIterations = 10000);

8a. Performance  for different learning rates (Python)

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn.linear_model
import pandas as pd

from sklearn.datasets import make_classification, make_blobs
execfile("./DLfunctions.py") # Since import does not work in Rmd!!!
# Create data
X1, Y1 = make_blobs(n_samples = 400, n_features = 2, centers = 9,
cluster_std = 1.3, random_state = 4)
#Create 2 classes
Y1=Y1.reshape(400,1)
Y1 = Y1 % 2
X2=X1.T
Y2=Y1.T
# Create a list of learning rates
learningRate=[0.5,1.2,3.0]
df=pd.DataFrame()
#Compute costs for each learning rate
for lr in learningRate:
parameters,costs = computeNN(X2, Y2, numHidden = 4, learningRate=lr, numIterations = 10000)
print(costs)
df1=pd.DataFrame(costs)
df=pd.concat([df,df1],axis=1)
#Set the iterations
iterations=[0,1000,2000,3000,4000,5000,6000,7000,8000,9000]
#Create data frame
#Set index
df1=df.set_index([iterations])
df1.columns=[0.5,1.2,3.0]
fig=df1.plot()
fig=plt.title("Cost vs No of Iterations for different learning rates")
plt.savefig('fig4.png', bbox_inches='tight')

8b. Performance  for different hidden units (Python)

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn.linear_model
import pandas as pd

from sklearn.datasets import make_classification, make_blobs
execfile("./DLfunctions.py") # Since import does not work in Rmd!!!
#Create data set
X1, Y1 = make_blobs(n_samples = 400, n_features = 2, centers = 9,
cluster_std = 1.3, random_state = 4)
#Create 2 classes
Y1=Y1.reshape(400,1)
Y1 = Y1 % 2
X2=X1.T
Y2=Y1.T
# Make a list of hidden unis
numHidden=[3,5,7]
df=pd.DataFrame()
#Compute costs for different hidden units
for numHid in numHidden:
parameters,costs = computeNN(X2, Y2, numHidden = numHid, learningRate=1.2, numIterations = 10000)
print(costs)
df1=pd.DataFrame(costs)
df=pd.concat([df,df1],axis=1)
#Set the iterations
iterations=[0,1000,2000,3000,4000,5000,6000,7000,8000,9000]
#Set index
df1=df.set_index([iterations])
df1.columns=[3,5,7]
#Plot
fig=df1.plot()
fig=plt.title("Cost vs No of Iterations for different no of hidden units")
plt.savefig('fig5.png', bbox_inches='tight')

9a. Performance  for different learning rates (R)

source("DLfunctions2_1.R")
x <- z[,1:2]
y <- z[,3]
x1 <- t(x)
y1 <- t(y)
#Loop through learning rates and compute costs
learningRate <-c(0.1,1.2,3.0)
df <- NULL
for(i in seq_along(learningRate)){
nn <-  computeNN(x1, y1, 6, learningRate=learningRate[i],numIterations=10000)
cost <- nn$costs df <- cbind(df,cost) }   #Create dataframe df <- data.frame(df) iterations=seq(0,10000,by=1000) df <- cbind(iterations,df) names(df) <- c("iterations","0.5","1.2","3.0") library(reshape2) df1 <- melt(df,id="iterations") # Melt the data #Plot ggplot(df1) + geom_line(aes(x=iterations,y=value,colour=variable),size=1) + xlab("Iterations") + ylab('Cost') + ggtitle("Cost vs No iterations for different learning rates") 9b. Performance for different hidden units (R) source("DLfunctions2_1.R") # Loop through Num hidden units numHidden <-c(4,6,9) df <- NULL for(i in seq_along(numHidden)){ nn <- computeNN(x1, y1, numHidden[i], learningRate=0.1,numIterations=10000) cost <- nn$costs
df <- cbind(df,cost)

}      
df <- data.frame(df)
iterations=seq(0,10000,by=1000)
df <- cbind(iterations,df)
names(df) <- c("iterations","4","6","9")
library(reshape2)
# Melt
df1 <- melt(df,id="iterations")
# Plot
ggplot(df1) + geom_line(aes(x=iterations,y=value,colour=variable),size=1)  +
xlab("Iterations") +
ylab('Cost') + ggtitle("Cost vs No iterations for  different number of hidden units")

10a. Performance of the Neural Network for different learning rates (Octave)

source("DL-function2.m") plotLRCostVsIterations() print -djph figa.jpg

10b. Performance of the Neural Network for different number of hidden units (Octave)

source("DL-function2.m") plotHiddenCostVsIterations() print -djph figa.jpg

11. Turning the heat on the Neural Network

In this 2nd part I create a a central region of positives and and the outside region as negatives. The points are generated using the equation of a circle (x – a)^{2} + (y -b) ^{2} = R^{2} . How does the 3 layer Neural Network perform on this?  Here’s a look! Note: The same dataset is also used for R and Octave Neural Network constructions

12. Manually creating a circular central region

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import sklearn.linear_model

from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification, make_blobs
from matplotlib.colors import ListedColormap
import sklearn
import sklearn.datasets

colors=['black','gold']
cmap = matplotlib.colors.ListedColormap(colors)
x1=np.random.uniform(0,10,800).reshape(800,1)
x2=np.random.uniform(0,10,800).reshape(800,1)
X=np.append(x1,x2,axis=1)
X.shape
# Create (x-a)^2 + (y-b)^2 = R^2
# Create a subset of values where squared is <0,4. Perform ravel() to flatten this vector
a=(np.power(X[:,0]-5,2) + np.power(X[:,1]-5,2) <= 6).ravel()
Y=a.reshape(800,1)

cmap = matplotlib.colors.ListedColormap(colors)

plt.figure()
plt.title('Non-linearly separable classes')
plt.scatter(X[:,0], X[:,1], c=Y,
marker= 'o', s=15,cmap=cmap)
plt.savefig('fig6.png', bbox_inches='tight')

13a. Decision boundary with hidden units=4 and learning rate = 2.2 (Python)

With the above hyper parameters the decision boundary is triangular

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import sklearn.linear_model
execfile("./DLfunctions.py")
x1=np.random.uniform(0,10,800).reshape(800,1)
x2=np.random.uniform(0,10,800).reshape(800,1)
X=np.append(x1,x2,axis=1)
X.shape

# Create a subset of values where squared is <0,4. Perform ravel() to flatten this vector
a=(np.power(X[:,0]-5,2) + np.power(X[:,1]-5,2) <= 6).ravel()
Y=a.reshape(800,1)

X2=X.T
Y2=Y.T

parameters,costs = computeNN(X2, Y2, numHidden = 4, learningRate=2.2, numIterations = 10000)
plot_decision_boundary(lambda x: predict(parameters, x.T), X2, Y2,str(4),str(2.2),"fig7.png")

## Cost after iteration 0: 0.692836
## Cost after iteration 1000: 0.331052
## Cost after iteration 2000: 0.326428
## Cost after iteration 3000: 0.474887
## Cost after iteration 4000: 0.247989
## Cost after iteration 5000: 0.218009
## Cost after iteration 6000: 0.201034
## Cost after iteration 7000: 0.197030
## Cost after iteration 8000: 0.193507
## Cost after iteration 9000: 0.191949

13b. Decision boundary with hidden units=12 and learning rate = 2.2 (Python)

With the above hyper parameters the decision boundary is triangular

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import sklearn.linear_model
execfile("./DLfunctions.py")
x1=np.random.uniform(0,10,800).reshape(800,1)
x2=np.random.uniform(0,10,800).reshape(800,1)
X=np.append(x1,x2,axis=1)
X.shape

# Create a subset of values where squared is <0,4. Perform ravel() to flatten this vector
a=(np.power(X[:,0]-5,2) + np.power(X[:,1]-5,2) <= 6).ravel()
Y=a.reshape(800,1)

X2=X.T
Y2=Y.T

parameters,costs = computeNN(X2, Y2, numHidden = 12, learningRate=2.2, numIterations = 10000)
plot_decision_boundary(lambda x: predict(parameters, x.T), X2, Y2,str(12),str(2.2),"fig8.png")

## Cost after iteration 0: 0.693291
## Cost after iteration 1000: 0.383318
## Cost after iteration 2000: 0.298807
## Cost after iteration 3000: 0.251735
## Cost after iteration 4000: 0.177843
## Cost after iteration 5000: 0.130414
## Cost after iteration 6000: 0.152400
## Cost after iteration 7000: 0.065359
## Cost after iteration 8000: 0.050921
## Cost after iteration 9000: 0.039719

14a. Decision boundary with hidden units=9 and learning rate = 0.5 (R)

When the number of hidden units is 6 and the learning rate is 0,1, is also a triangular shape in R

source("DLfunctions2_1.R")
x <- z[,1:2]
y <- z[,3]
x1 <- t(x)
y1 <- t(y)
nn <-computeNN(x1, y1, 9, learningRate=0.5,numIterations=10000) # Triangular
## [1] 0.8398838
## [1] 0.3303621
## [1] 0.3127731
## [1] 0.3012791
## [1] 0.3305543
## [1] 0.3303964
## [1] 0.2334615
## [1] 0.1920771
## [1] 0.2341225
## [1] 0.2188118
## [1] 0.2082687
plotDecisionBoundary(z,nn,6,0.1)

14b. Decision boundary with hidden units=8 and learning rate = 0.1 (R)

source("DLfunctions2_1.R")
x <- z[,1:2]
y <- z[,3]
x1 <- t(x)
y1 <- t(y)
nn <-computeNN(x1, y1, 8, learningRate=0.1,numIterations=10000) # Hemisphere
## [1] 0.7273279
## [1] 0.3169335
## [1] 0.2378464
## [1] 0.1688635
## [1] 0.1368466
## [1] 0.120664
## [1] 0.111211
## [1] 0.1043362
## [1] 0.09800573
## [1] 0.09126161
## [1] 0.0840379
plotDecisionBoundary(z,nn,8,0.1)

15a. Decision boundary with hidden units=12 and learning rate = 1.5 (Octave)

source("DL-function2.m") data=csvread("data1.csv"); X=data(:,1:2); Y=data(:,3); # Make sure that the model parameters are correct. Take the transpose of X & Y [W1,b1,W2,b2,costs]= computeNN(X', Y',12, learningRate=1.5, numIterations = 10000); plotDecisionBoundary(data, W1,b1,W2,b2) print -djpg fige.jpg

Conclusion: This post implemented a 3 layer Neural Network to create non-linear boundaries while performing classification. Clearly the Neural Network performs very well when the number of hidden units and learning rate are varied.

To be continued…
Watch this space!!

To see all posts check Index of posts

Deep Learning from first principles in Python, R and Octave – Part 1

“You don’t perceive objects as they are. You perceive them as you are.”
“Your interpretation of physical objects has everything to do with the historical trajectory of your brain – and little to do with the objects themselves.”
“The brain generates its own reality, even before it receives information coming in from the eyes and the other senses. This is known as the internal model”

                          David Eagleman - The Brain: The Story of You

This is the first in the series of posts, I intend to write on Deep Learning. This post is inspired by the Deep Learning Specialization by Prof Andrew Ng on Coursera and Neural Networks for Machine Learning by Prof Geoffrey Hinton also on Coursera. In this post I implement Logistic regression with a 2 layer Neural Network i.e. a Neural Network that just has an input layer and an output layer and with no hidden layer.I am certain that any self-respecting Deep Learning/Neural Network would consider a Neural Network without hidden layers as no Neural Network at all!

This 2 layer network is implemented in Python, R and Octave languages. I have included Octave, into the mix, as Octave is a close cousin of Matlab. These implementations in Python, R and Octave are equivalent vectorized implementations. So, if you are familiar in any one of the languages, you should be able to look at the corresponding code in the other two. You can download this R Markdown file and Octave code from DeepLearning -Part 1

To start with, Logistic Regression is performed using sklearn’s logistic regression package for the cancer data set also from sklearn. This is shown below

1. Logistic Regression

import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_classification, make_blobs

from sklearn.metrics import confusion_matrix
from matplotlib.colors import ListedColormap
(X_cancer, y_cancer) = load_breast_cancer(return_X_y = True)
X_train, X_test, y_train, y_test = train_test_split(X_cancer, y_cancer,
random_state = 0)
# Call the Logisitic Regression function
clf = LogisticRegression().fit(X_train, y_train)
print('Accuracy of Logistic regression classifier on training set: {:.2f}'
.format(clf.score(X_train, y_train)))
print('Accuracy of Logistic regression classifier on test set: {:.2f}'
.format(clf.score(X_test, y_test)))
## Accuracy of Logistic regression classifier on training set: 0.96
## Accuracy of Logistic regression classifier on test set: 0.96

To check on other classification algorithms, check my post Practical Machine Learning with R and Python – Part 2. You could also take a look at my book, available on Amazon , in which I implement the most popular Machine Learning algorithms in both R and Python – My book ‘Practical Machine Learning with R and Python’ on Amazon

2. Logistic Regression as a 2 layer Neural Network

In the following section Logistic Regression is implemented as a 2 layer Neural Network in Python, R and Octave. The same cancer data set from sklearn will be used to train and test the Neural Network in Python, R and Octave. This can be represented diagrammatically as below

The cancer data set has 30 input features, and the target variable ‘output’ is either 0 or 1. Hence the sigmoid activation function will be used in the output layer for classification.

This simple 2 layer Neural Network is shown below
At the input layer there are 30 features and the corresponding weights of these inputs which are initialized to small random values.
$Z= w_{1}x_{1} +w_{2}x_{2} +..+ w_{30}x_{30} + b$
where ‘b’ is the bias term

The Activation function is the sigmoid function which is $a= 1/(1+e^{-z})$
The Loss, when the sigmoid function is used in the output layer, is given by
$L=-(ylog(a) + (1-y)log(1-a))$ (1)

Forward propagation

In forward propagation cycle of the Neural Network the output Z and the output of activation function, the sigmoid function, is first computed. Then using the output ‘y’ for the given features, the ‘Loss’ is computed using equation (1) above.

Backward propagation

The backward propagation cycle determines how the ‘Loss’ is impacted for small variations from the previous layers upto the input layer. In other words, backward propagation computes the changes in the weights at the input layer, which will minimize the loss. Several cycles of gradient descent are performed in the path of steepest descent to find the local minima. In other words the set of weights and biases, at the input layer, which will result in the lowest loss is computed by gradient descent. The weights at the input layer are decreased by a parameter known as the ‘learning rate’. Too big a ‘learning rate’ can overshoot the local minima, and too small a ‘learning rate’ can take a long time to reach the local minima. This is done for ‘m’ training examples.

Chain rule of differentiation
Let y=f(u)
and u=g(x) then
$\partial y/\partial x = \partial y/\partial u * \partial u/\partial x$

Derivative of sigmoid
$\sigma=1/(1+e^{-z})$
Let $x= 1 + e^{-z}$  then
$\sigma = 1/x$
$\partial \sigma/\partial x = -1/x^{2}$
$\partial x/\partial z = -e^{-z}$
Using the chain rule of differentiation we get
$\partial \sigma/\partial z = \partial \sigma/\partial x * \partial x/\partial z$
$=-1/(1+e^{-z})^{2}* -e^{-z} = e^{-z}/(1+e^{-z})^{2}$
Therefore $\partial \sigma/\partial z = \sigma(1-\sigma)$        -(2)

The 3 equations for the 2 layer Neural Network representation of Logistic Regression are
$L=-(y*log(a) + (1-y)*log(1-a))$      -(a)
$a=1/(1+e^{-z})$      -(b)
$Z= w_{1}x_{1} +w_{2}x_{2} +...+ w_{30}x_{30} +b$   -(c)

The back propagation step requires the computation of $dL/dw_{i}$ and $dL/db_{i}$. In the case of regression it would be $dE/dw_{i}$ and $dE/db_{i}$ where dE is the Mean Squared Error function.
Computing the derivatives for back propagation we have
$dL/da = -(y/a + (1-y)/(1-a))$          -(d)
because $d/dx(logx) = 1/x$
Also from equation (2) we get
$da/dZ = a (1-a)$                                  – (e)
By chain rule
$\partial L/\partial Z = \partial L/\partial a * \partial a/\partial Z$
therefore substituting the results of (d) & (e) we get
$\partial L/\partial Z = -(y/a + (1-y)/(1-a)) * a(1-a) = a-y$         (f)
Finally
$\partial L/\partial w_{i}= \partial L/\partial a * \partial a/\partial Z * \partial Z/\partial w_{i}$                                                           -(g)
$\partial Z/\partial w_{i} = x_{i}$            – (h)
and from (f) we have  $\partial L/\partial Z =a-y$
Therefore  (g) reduces to
$\partial L/\partial w_{i} = x_{i}* (a-y)$ -(i)
Also
$\partial L/\partial b = \partial L/\partial a * \partial a/\partial Z * \partial Z/\partial b$ -(j)
Since
$\partial Z/\partial b = 1$
$\partial L/\partial b = a-y$

The gradient computes the weights at the input layer and the corresponding bias by using the values
of $dw_{i}$ and $db$
$w_{i} := w_{i} -\alpha * dw_{i}$
$b := b -\alpha * db$
I found the computation graph representation in the book Deep Learning: Ian Goodfellow, Yoshua Bengio, Aaron Courville, very useful to visualize and also compute the backward propagation. For the 2 layer Neural Network of Logistic Regression the computation graph is shown below

3. Neural Network for Logistic Regression -Python code (vectorized)

import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

# Define the sigmoid function
def sigmoid(z):
a=1/(1+np.exp(-z))
return a

# Initialize
def initialize(dim):
w = np.zeros(dim).reshape(dim,1)
b = 0
return w

# Compute the loss
def computeLoss(numTraining,Y,A):
loss=-1/numTraining *np.sum(Y*np.log(A) + (1-Y)*(np.log(1-A)))
return(loss)

# Execute the forward propagation
def forwardPropagation(w,b,X,Y):
# Compute Z
Z=np.dot(w.T,X)+b
# Determine the number of training samples
numTraining=float(len(X))
# Compute the output of the sigmoid activation function
A=sigmoid(Z)
#Compute the loss
loss = computeLoss(numTraining,Y,A)
# Compute the gradients dZ, dw and db
dZ=A-Y
dw=1/numTraining*np.dot(X,dZ.T)
db=1/numTraining*np.sum(dZ)

# Return the results as a dictionary
"db": db}
loss = np.squeeze(loss)

def gradientDescent(w, b, X, Y, numIerations, learningRate):
losses=[]
idx =[]
# Iterate
for i in range(numIerations):
#Get the derivates
w = w-learningRate*dw
b = b-learningRate*db

# Store the loss
if i % 100 == 0:
idx.append(i)
losses.append(loss)
params = {"w": w,
"b": b}
"db": db}

# Predict the output for a training set
def predict(w,b,X):
size=X.shape[1]
yPredicted=np.zeros((1,size))
Z=np.dot(w.T,X)
# Compute the sigmoid
A=sigmoid(Z)
for i in range(A.shape[1]):
#If the value is > 0.5 then set as 1
if(A[0][i] > 0.5):
yPredicted[0][i]=1
else:
# Else set as 0
yPredicted[0][i]=0

return yPredicted

#Normalize the data
def normalize(x):
x_norm = None
x_norm = np.linalg.norm(x,axis=1,keepdims=True)
x= x/x_norm
return x

# Run the 2 layer Neural Network on the cancer data set

(X_cancer, y_cancer) = load_breast_cancer(return_X_y = True)
# Create train and test sets
X_train, X_test, y_train, y_test = train_test_split(X_cancer, y_cancer,
random_state = 0)
# Normalize the data for better performance
X_train1=normalize(X_train)

# Create weight vectors of zeros. The size is the number of features in the data set=30
w=np.zeros((X_train.shape[1],1))
#w=np.zeros((30,1))
b=0

#Normalize the training data so that gradient descent performs better
X_train1=normalize(X_train)
#Transpose X_train so that we have a matrix as (features, numSamples)
X_train2=X_train1.T

# Reshape to remove the rank 1 array and then transpose
y_train1=y_train.reshape(len(y_train),1)
y_train2=y_train1.T

# Run gradient descent for 4000 times and compute the weights
w = parameters["w"]
b = parameters["b"]

# Normalize X_test
X_test1=normalize(X_test)
#Transpose X_train so that we have a matrix as (features, numSamples)
X_test2=X_test1.T

#Reshape y_test
y_test1=y_test.reshape(len(y_test),1)
y_test2=y_test1.T

# Predict the values for
yPredictionTest = predict(w, b, X_test2)
yPredictionTrain = predict(w, b, X_train2)

# Print the accuracy
print("train accuracy: {} %".format(100 - np.mean(np.abs(yPredictionTrain - y_train2)) * 100))
print("test accuracy: {} %".format(100 - np.mean(np.abs(yPredictionTest - y_test)) * 100))

# Plot the Costs vs the number of iterations
fig1=plt.plot(idx,costs)
fig1=plt.title("Gradient descent-Cost vs No of iterations")
fig1=plt.xlabel("No of iterations")
fig1=plt.ylabel("Cost")
fig1.figure.savefig("fig1", bbox_inches='tight')
## train accuracy: 90.3755868545 %
## test accuracy: 89.5104895105 %

Note: It can be seen that the Accuracy on the training and test set is 90.37% and 89.51%. This is comparatively poorer than the 96% which the logistic regression of sklearn achieves! But this is mainly because of the absence of hidden layers which is the real power of neural networks.

4. Neural Network for Logistic Regression -R code (vectorized)

source("RFunctions-1.R")
# Define the sigmoid function
sigmoid <- function(z){
a <- 1/(1+ exp(-z))
a
}

# Compute the loss
computeLoss <- function(numTraining,Y,A){
loss <- -1/numTraining* sum(Y*log(A) + (1-Y)*log(1-A))
return(loss)
}

# Compute forward propagation
forwardPropagation <- function(w,b,X,Y){
# Compute Z
Z <- t(w) %*% X +b
#Set the number of samples
numTraining <- ncol(X)
# Compute the activation function
A=sigmoid(Z)

#Compute the loss
loss <- computeLoss(numTraining,Y,A)

# Compute the gradients dZ, dw and db
dZ<-A-Y
dw<-1/numTraining * X %*% t(dZ)
db<-1/numTraining*sum(dZ)

fwdProp <- list("loss" = loss, "dw" = dw, "db" = db)
return(fwdProp)
}

# Perform one cycle of Gradient descent
gradientDescent <- function(w, b, X, Y, numIerations, learningRate){
losses <- NULL
idx <- NULL
# Loop through the number of iterations
for(i in 1:numIerations){
fwdProp <-forwardPropagation(w,b,X,Y)
#Get the derivatives
dw <- fwdProp$dw db <- fwdProp$db
w = w-learningRate*dw
b = b-learningRate*db
l <- fwdProp$loss # Stoe the loss if(i %% 100 == 0){ idx <- c(idx,i) losses <- c(losses,l) } } # Return the weights and losses gradDescnt <- list("w"=w,"b"=b,"dw"=dw,"db"=db,"losses"=losses,"idx"=idx) return(gradDescnt) } # Compute the predicted value for input predict <- function(w,b,X){ m=dim(X)[2] # Create a ector of 0's yPredicted=matrix(rep(0,m),nrow=1,ncol=m) Z <- t(w) %*% X +b # Compute sigmoid A=sigmoid(Z) for(i in 1:dim(A)[2]){ # If A > 0.5 set value as 1 if(A[1,i] > 0.5) yPredicted[1,i]=1 else # Else set as 0 yPredicted[1,i]=0 } return(yPredicted) } # Normalize the matrix normalize <- function(x){ #Create the norm of the matrix.Perform the Frobenius norm of the matrix n<-as.matrix(sqrt(rowSums(x^2))) #Sweep by rows by norm. Note '1' in the function which performing on every row normalized<-sweep(x, 1, n, FUN="/") return(normalized) } # Run the 2 layer Neural Network on the cancer data set # Read the data (from sklearn) cancer <- read.csv("cancer.csv") # Rename the target variable names(cancer) <- c(seq(1,30),"output") # Split as training and test sets train_idx <- trainTestSplit(cancer,trainPercent=75,seed=5) train <- cancer[train_idx, ] test <- cancer[-train_idx, ] # Set the features X_train <-train[,1:30] y_train <- train[,31] X_test <- test[,1:30] y_test <- test[,31] # Create a matrix of 0's with the number of features w <-matrix(rep(0,dim(X_train)[2])) b <-0 X_train1 <- normalize(X_train) X_train2=t(X_train1) # Reshape then transpose y_train1=as.matrix(y_train) y_train2=t(y_train1) # Perform gradient descent gradDescent= gradientDescent(w, b, X_train2, y_train2, numIerations=3000, learningRate=0.77) # Normalize X_test X_test1=normalize(X_test) #Transpose X_train so that we have a matrix as (features, numSamples) X_test2=t(X_test1) #Reshape y_test and take transpose y_test1=as.matrix(y_test) y_test2=t(y_test1) # Use the values of the weights generated from Gradient Descent yPredictionTest = predict(gradDescent$w, gradDescent$b, X_test2) yPredictionTrain = predict(gradDescent$w, gradDescent$b, X_train2) sprintf("Train accuracy: %f",(100 - mean(abs(yPredictionTrain - y_train2)) * 100)) ## [1] "Train accuracy: 90.845070" sprintf("test accuracy: %f",(100 - mean(abs(yPredictionTest - y_test)) * 100)) ## [1] "test accuracy: 87.323944" df <-data.frame(gradDescent$idx, gradDescent$losses) names(df) <- c("iterations","losses") ggplot(df,aes(x=iterations,y=losses)) + geom_point() + geom_line(col="blue") + ggtitle("Gradient Descent - Losses vs No of Iterations") + xlab("No of iterations") + ylab("Losses") 4. Neural Network for Logistic Regression -Octave code (vectorized)  1; # Define sigmoid function function a = sigmoid(z) a = 1 ./ (1+ exp(-z)); end # Compute the loss function loss=computeLoss(numtraining,Y,A) loss = -1/numtraining * sum((Y .* log(A)) + (1-Y) .* log(1-A)); end  # Perform forward propagation function [loss,dw,db,dZ] = forwardPropagation(w,b,X,Y) % Compute Z Z = w' * X + b; numtraining = size(X)(1,2); # Compute sigmoid A = sigmoid(Z);  #Compute loss. Note this is element wise product loss =computeLoss(numtraining,Y,A); # Compute the gradients dZ, dw and db dZ = A-Y; dw = 1/numtraining* X * dZ'; db =1/numtraining*sum(dZ); end  # Compute Gradient Descent function [w,b,dw,db,losses,index]=gradientDescent(w, b, X, Y, numIerations, learningRate) #Initialize losses and idx losses=[]; index=[]; # Loop through the number of iterations for i=1:numIerations, [loss,dw,db,dZ] = forwardPropagation(w,b,X,Y); # Perform Gradient descent w = w - learningRate*dw; b = b - learningRate*db; if(mod(i,100) ==0) # Append index and loss index = [index i]; losses = [losses loss]; endif end end  # Determine the predicted value for dataset function yPredicted = predict(w,b,X) m = size(X)(1,2); yPredicted=zeros(1,m); # Compute Z Z = w' * X + b; # Compute sigmoid A = sigmoid(Z); for i=1:size(X)(1,2), # Set predicted as 1 if A > 0,5 if(A(1,i) >= 0.5) yPredicted(1,i)=1; else yPredicted(1,i)=0; endif end end  # Normalize by dividing each value by the sum of squares function normalized = normalize(x) # Compute Frobenius norm. Square the elements, sum rows and then find square root a = sqrt(sum(x .^ 2,2)); # Perform element wise division normalized = x ./ a; end  # Split into train and test sets function [X_train,y_train,X_test,y_test] = trainTestSplit(dataset,trainPercent) # Create a random index ix = randperm(length(dataset)); # Split into training trainSize = floor(trainPercent/100 * length(dataset)); train=dataset(ix(1:trainSize),:); # And test test=dataset(ix(trainSize+1:length(dataset)),:); X_train = train(:,1:30); y_train = train(:,31); X_test = test(:,1:30); y_test = test(:,31); end  cancer=csvread("cancer.csv"); [X_train,y_train,X_test,y_test] = trainTestSplit(cancer,75); w=zeros(size(X_train)(1,2),1); b=0; X_train1=normalize(X_train); X_train2=X_train1'; y_train1=y_train'; [w1,b1,dw,db,losses,idx]=gradientDescent(w, b, X_train2, y_train1, numIerations=3000, learningRate=0.75); # Normalize X_test X_test1=normalize(X_test); #Transpose X_train so that we have a matrix as (features, numSamples) X_test2=X_test1'; y_test1=y_test'; # Use the values of the weights generated from Gradient Descent yPredictionTest = predict(w1, b1, X_test2); yPredictionTrain = predict(w1, b1, X_train2);   trainAccuracy=100-mean(abs(yPredictionTrain - y_train1))*100 testAccuracy=100- mean(abs(yPredictionTest - y_test1))*100 trainAccuracy = 90.845 testAccuracy = 89.510 graphics_toolkit('gnuplot') plot(idx,losses); title ('Gradient descent- Cost vs No of iterations'); xlabel ("No of iterations"); ylabel ("Cost"); Conclusion This post starts with a simple 2 layer Neural Network implementation of Logistic Regression. Clearly the performance of this simple Neural Network is comparatively poor to the highly optimized sklearn’s Logistic Regression. This is because the above neural network did not have any hidden layers. Deep Learning & Neural Networks achieve extraordinary performance because of the presence of deep hidden layers The Deep Learning journey has begun… Don’t miss the bus! Stay tuned for more interesting posts in Deep Learning!! To see all posts check Index of posts The 3rd paperback & kindle editions of my books on Cricket, now on Amazon The 3rd paperback & kindle edition of both my books on cricket is now available on Amazon a) Cricket analytics with cricketr, Third Edition. The paperback edition is$12.99 and the kindle edition is $4.99/Rs320. This book is based on my R package ‘cricketr‘, available on CRAN and uses ESPN Cricinfo Statsguru b) Beaten by sheer pace! Cricket analytics with yorkr, 3rd edition . The paperback is$12.99 and the kindle version is $6.99/Rs448. This is based on my R package ‘yorkr‘ on CRAN and uses data from Cricsheet Pick up your copies today!! Note: In the 3rd edition of the paperback book, the charts will be in black and white. If you would like the charts to be in color, please check out the 2nd edition of these books see More book, more cricket! 2nd edition of my books now on Amazon To see all posts see Index of posts My book ‘Practical Machine Learning with R and Python’ on Amazon My book ‘Practical Machine Learning with R and Python – Machine Learning in stereo’ is now available in both paperback ($9.99) and kindle ($6.97/Rs449) versions. In this book I implement some of the most common, but important Machine Learning algorithms in R and equivalent Python code. This is almost like listening to parallel channels of music in stereo! 1. Practical machine with R and Python – Machine Learning in Stereo (Paperback) 2. Practical machine with R and Python – Machine Learning in Stereo (Kindle) This book is ideal both for beginners and the experts in R and/or Python. Those starting their journey into datascience and ML will find the first 3 chapters useful, as they touch upon the most important programming constructs in R and Python and also deal with equivalent statements in R and Python. Those who are expert in either of the languages, R or Python, will find the equivalent code ideal for brushing up on the other language. And finally,those who are proficient in both languages, can use the R and Python implementations to internalize the ML algorithms better. Here is a look at the topics covered Table of Contents Essential R …………………………………….. 7 Essential Python for Datascience ……………….. 54 R vs Python ……………………………………. 77 Regression of a continuous variable ………………. 96 Classification and Cross Validation ……………….113 Regression techniques and regularization …………. 134 SVMs, Decision Trees and Validation curves …………175 Splines, GAMs, Random Forests and Boosting …………202 PCA, K-Means and Hierarchical Clustering …………. 234 Pick up your copy today!! Hope you have a great time learning as I did while implementing these algorithms! Practical Machine Learning with R and Python – Part 6 Introduction This is the final and concluding part of my series on ‘Practical Machine Learning with R and Python’. In this series I included the implementations of the most common Machine Learning algorithms in R and Python. The algorithms implemented were 1. Practical Machine Learning with R and Python – Part 1 In this initial post, I touch upon regression of a continuous target variable. Specifically I touch upon Univariate, Multivariate, Polynomial regression and KNN regression in both R and Python 2. Practical Machine Learning with R and Python – Part 2 In this post, I discuss Logistic Regression, KNN classification and Cross Validation error for both LOOCV and K-Fold in both R and Python 3. Practical Machine Learning with R and Python – Part 3 This 3rd part included feature selection in Machine Learning. Specifically I touch best fit, forward fit, backward fit, ridge(L2 regularization) & lasso (L1 regularization). The post includes equivalent code in R and Python. 4. Practical Machine Learning with R and Python – Part 4 In this part I discussed SVMs, Decision Trees, Validation, Precision-Recall, AUC and ROC curves 5. Practical Machine Learning with R and Python – Part 5 In this penultimate part, I touch upon B-splines, natural splines, smoothing spline, Generalized Additive Models(GAMs), Decision Trees, Random Forests and Gradient Boosted Treess. In this last part I cover Unsupervised Learning. Specifically I cover the implementations of Principal Component Analysis (PCA). K-Means and Heirarchical Clustering. You can download this R Markdown file from Github at MachineLearning-RandPython-Part6 The content of this post and much more is now available as a compact book on Amazon in both formats – as Paperback ($9.99) and a Kindle version($6.99/Rs449/). see ‘Practical Machine Learning with R and Python – Machine Learning in stereo 1.1a Principal Component Analysis (PCA) – R code Principal Component Analysis is used to reduce the dimensionality of the input. In the code below 8 x 8 pixel of handwritten digits is reduced into its principal components. Then a scatter plot of the first 2 principal components give a very good visial representation of the data library(dplyr) library(ggplot2) #Note: This example is adapted from an the example in the book Python Datascience handbook by # Jake VanderPlas (https://jakevdp.github.io/PythonDataScienceHandbook/05.09-principal-component-analysis.html) # Read the digits data (From sklearn datasets) digits= read.csv("digits.csv") # Create a digits classes target variable digitClasses <- factor(digits$X0.000000000000000000e.00.29)

#Invoke the Principal Componsent analysis on columns 1-64
digitsPCA=prcomp(digits[,1:64])

# Create a dataframe of PCA
df <- data.frame(digitsPCA$x) # Bind the digit classes df1 <- cbind(df,digitClasses) # Plot only the first 2 Principal components as a scatter plot. This plot uses only the # first 2 principal components ggplot(df1,aes(x=PC1,y=PC2,col=digitClasses)) + geom_point() + ggtitle("Top 2 Principal Components") 1.1 b Variance explained vs no principal components – R code In the code below the variance explained vs the number of principal components is plotted. It can be seen that with 20 Principal components almost 90% of the variance is explained by this reduced dimensional model. # Read the digits data (from sklearn datasets) digits= read.csv("digits.csv") # Digits target digitClasses <- factor(digits$X0.000000000000000000e.00.29)
digitsPCA=prcomp(digits[,1:64])

# Get the Standard Deviation
sd=digitsPCA$sdev # Compute the variance digitsVar=digitsPCA$sdev^2
#Compute the percent variance explained
percentVarExp=digitsVar/sum(digitsVar)

# Plot the percent variance exlained as a function of the  number of principal components
#plot(cumsum(percentVarExp), xlab="Principal Component",
#     ylab="Cumulative Proportion of Variance Explained",
#     main="Principal Components vs % Variance explained",ylim=c(0,1),type='l',lwd=2,
#       col="blue")

1.1c Principal Component Analysis (PCA) – Python code

import numpy as np
from sklearn.decomposition import PCA
from sklearn import decomposition
from sklearn import datasets
import matplotlib.pyplot as plt

# Select only the first 2 principal components
pca = PCA(2)  # project from 64 to 2 dimensions
#Compute the first 2 PCA
projected = pca.fit_transform(digits.data)

# Plot a scatter plot of the first 2 principal components
plt.scatter(projected[:, 0], projected[:, 1],
c=digits.target, edgecolor='none', alpha=0.5,
cmap=plt.cm.get_cmap('spectral', 10))
plt.xlabel('PCA 1')
plt.ylabel('PCA 2')
plt.colorbar();
plt.title("Top 2 Principal Components")
plt.savefig('fig1.png', bbox_inches='tight')

– Python code

import numpy as np
from sklearn.decomposition import PCA
from sklearn import decomposition
from sklearn import datasets
import matplotlib.pyplot as plt

# Select all 64 principal components
pca = PCA(64)  # project from 64 to 2 dimensions
projected = pca.fit_transform(digits.data)

# Obtain the explained variance for each principal component
varianceExp= pca.explained_variance_ratio_
# Compute the total sum of variance
totVarExp=np.cumsum(np.round(pca.explained_variance_ratio_, decimals=4)*100)

# Plot the variance explained as a function of the number of principal components
plt.plot(totVarExp)
plt.xlabel('No of principal components')
plt.ylabel('% variance explained')
plt.title('No of Principal Components vs Total Variance explained')
plt.savefig('fig2.png', bbox_inches='tight')

1.2a K-Means – R code

In the code first the scatter plot of the first 2 Principal Components of the handwritten digits is plotted as a scatter plot. Over this plot 10 centroids of the 10 different clusters corresponding the 10 diferent digits is plotted over the original scatter plot.

library(ggplot2)
# Create digit classes target variable
digitClasses <- factor(digits$X0.000000000000000000e.00.29) # Compute the Principal COmponents digitsPCA=prcomp(digits[,1:64]) # Create a data frame of Principal components and the digit classes df <- data.frame(digitsPCA$x)
df1 <- cbind(df,digitClasses)

# Pick only the first 2 principal components
a<- df[,1:2]
# Compute K Means of 10 clusters and allow for 1000 iterations
k<-kmeans(a,10,1000)

# Create a dataframe of the centroids of the clusters
df2<-data.frame(k$centers) #Plot the first 2 principal components with the K Means centroids ggplot(df1,aes(x=PC1,y=PC2,col=digitClasses)) + geom_point() + geom_point(data=df2,aes(x=PC1,y=PC2),col="black",size = 4) + ggtitle("Top 2 Principal Components with KMeans clustering")  1.2b K-Means – Python code The centroids of the 10 different handwritten digits is plotted over the scatter plot of the first 2 principal components. import numpy as np from sklearn.decomposition import PCA from sklearn import decomposition from sklearn import datasets import matplotlib.pyplot as plt from sklearn.datasets import load_digits from sklearn.cluster import KMeans digits = load_digits() # Select only the 1st 2 principal components pca = PCA(2) # project from 64 to 2 dimensions projected = pca.fit_transform(digits.data) # Create 10 different clusters kmeans = KMeans(n_clusters=10) # Compute the clusters kmeans.fit(projected) y_kmeans = kmeans.predict(projected) # Get the cluster centroids centers = kmeans.cluster_centers_ centers #Create a scatter plot of the first 2 principal components plt.scatter(projected[:, 0], projected[:, 1], c=digits.target, edgecolor='none', alpha=0.5, cmap=plt.cm.get_cmap('spectral', 10)) plt.xlabel('PCA 1') plt.ylabel('PCA 2') plt.colorbar(); # Overlay the centroids on the scatter plot plt.scatter(centers[:, 0], centers[:, 1], c='darkblue', s=100) plt.savefig('fig3.png', bbox_inches='tight') 1.3a Heirarchical clusters – R code Herirachical clusters is another type of unsupervised learning. It successively joins the closest pair of objects (points or clusters) in succession based on some ‘distance’ metric. In this type of clustering we do not have choose the number of centroids. We can cut the created dendrogram mat an appropriate height to get a desired and reasonable number of clusters These are the following ‘distance’ metrics used while combining successive objects • Ward • Complete • Single • Average • Centroid # Read the IRIS dataset iris <- datasets::iris iris2 <- iris[,-5] species <- iris[,5] #Compute the distance matrix d_iris <- dist(iris2) # Use the 'average' method to for the clsuters hc_iris <- hclust(d_iris, method = "average") # Plot the clusters plot(hc_iris) # Cut tree into 3 groups sub_grp <- cutree(hc_iris, k = 3) # Number of members in each cluster table(sub_grp) ## sub_grp ## 1 2 3 ## 50 64 36 # Draw rectangles around the clusters rect.hclust(hc_iris, k = 3, border = 2:5) 1.3a Heirarchical clusters – Python code from sklearn.datasets import load_iris import matplotlib.pyplot as plt from scipy.cluster.hierarchy import dendrogram, linkage # Load the IRIS data set iris = load_iris() # Generate the linkage matrix using the average method Z = linkage(iris.data, 'average') #Plot the dendrogram #dendrogram(Z) #plt.xlabel('Data') #plt.ylabel('Distance') #plt.suptitle('Samples clustering', fontweight='bold', fontsize=14); #plt.savefig('fig4.png', bbox_inches='tight') Conclusion This is the last and concluding part of my series on Practical Machine Learning with R and Python. These parallel implementations of R and Python can be used as a quick reference while working on a large project. A person who is adept in one of the languages R or Python, can quickly absorb code in the other language. Hope you find this series useful! More interesting things to come. Watch this space! References 1. Statistical Learning, Prof Trevor Hastie & Prof Robert Tibesherani, Online Stanford 2. Applied Machine Learning in Python Prof Kevyn-Collin Thomson, University Of Michigan, Coursera To see all posts see ‘Index of posts Practical Machine Learning with R and Python – Part 5 This is the 5th and probably penultimate part of my series on ‘Practical Machine Learning with R and Python’. The earlier parts of this series included 1. Practical Machine Learning with R and Python – Part 1 In this initial post, I touch upon univariate, multivariate, polynomial regression and KNN regression in R and Python 2.Practical Machine Learning with R and Python – Part 2 In this post, I discuss Logistic Regression, KNN classification and cross validation error for both LOOCV and K-Fold in both R and Python 3.Practical Machine Learning with R and Python – Part 3 This post covered ‘feature selection’ in Machine Learning. Specifically I touch best fit, forward fit, backward fit, ridge(L2 regularization) & lasso (L1 regularization). The post includes equivalent code in R and Python. 4.Practical Machine Learning with R and Python – Part 4 In this part I discussed SVMs, Decision Trees, validation, precision recall, and roc curves This post ‘Practical Machine Learning with R and Python – Part 5’ discusses regression with B-splines, natural splines, smoothing splines, generalized additive models (GAMS), bagging, random forest and boosting As with my previous posts in this series, this post is largely based on the following 2 MOOC courses 1. Statistical Learning, Prof Trevor Hastie & Prof Robert Tibesherani, Online Stanford 2. Applied Machine Learning in Python Prof Kevyn-Collin Thomson, University Of Michigan, Coursera You can download this R Markdown file and associated data files from Github at MachineLearning-RandPython-Part5 The content of this post and much more is now available as a compact book on Amazon in both formats – as Paperback ($9.99) and a Kindle version($6.99/Rs449/). see ‘Practical Machine Learning with R and Python – Machine Learning in stereo For this part I have used the data sets from UCI Machine Learning repository(Communities and Crime and Auto MPG) 1. Splines When performing regression (continuous or logistic) between a target variable and a feature (or a set of features), a single polynomial for the entire range of the data set usually does not perform a good fit.Rather we would need to provide we could fit regression curves for different section of the data set. There are several techniques which do this for e.g. piecewise-constant functions, piecewise-linear functions, piecewise-quadratic/cubic/4th order polynomial functions etc. One such set of functions are the cubic splines which fit cubic polynomials to successive sections of the dataset. The points where the cubic splines join, are called ‘knots’. Since each section has a different cubic spline, there could be discontinuities (or breaks) at these knots. To prevent these discontinuities ‘natural splines’ and ‘smoothing splines’ ensure that the seperate cubic functions have 2nd order continuity at these knots with the adjacent splines. 2nd order continuity implies that the value, 1st order derivative and 2nd order derivative at these knots are equal. A cubic spline with knots $\alpha_{k}$ , k=1,2,3,..K is a piece-wise cubic polynomial with continuous derivative up to order 2 at each knot. We can write $y_{i} = \beta_{0} +\beta_{1}b_{1}(x_{i}) +\beta_{2}b_{2}(x_{i}) + .. + \beta_{K+3}b_{K+3}(x_{i}) + \epsilon_{i}$. For each ($x{i},y{i}$), $b_{i}$ are called ‘basis’ functions, where $b_{1}(x_{i})=x_{i}$$b_{2}(x_{i})=x_{i}^2$, $b_{3}(x_{i})=x_{i}^3$, $b_{k+3}(x_{i})=(x_{i} -\alpha_{k})^3$ where k=1,2,3… K The 1st and 2nd derivatives of cubic splines are continuous at the knots. Hence splines provide a smooth continuous fit to the data by fitting different splines to different sections of the data 1.1a Fit a 4th degree polynomial – R code In the code below a non-linear function (a 4th order polynomial) is used to fit the data. Usually when we fit a single polynomial to the entire data set the tails of the fit tend to vary a lot particularly if there are fewer points at the ends. Splines help in reducing this variation at the extremities library(dplyr) library(ggplot2) source('RFunctions-1.R') # Read the data df=read.csv("auto_mpg.csv",stringsAsFactors = FALSE) # Data from UCI df1 <- as.data.frame(sapply(df,as.numeric)) #Select specific columns df2 <- df1 %>% dplyr::select(cylinder,displacement, horsepower,weight, acceleration, year,mpg) auto <- df2[complete.cases(df2),] # Fit a 4th degree polynomial fit=lm(mpg~poly(horsepower,4),data=auto) #Display a summary of fit summary(fit) ## ## Call: ## lm(formula = mpg ~ poly(horsepower, 4), data = auto) ## ## Residuals: ## Min 1Q Median 3Q Max ## -14.8820 -2.5802 -0.1682 2.2100 16.1434 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 23.4459 0.2209 106.161 <2e-16 *** ## poly(horsepower, 4)1 -120.1377 4.3727 -27.475 <2e-16 *** ## poly(horsepower, 4)2 44.0895 4.3727 10.083 <2e-16 *** ## poly(horsepower, 4)3 -3.9488 4.3727 -0.903 0.367 ## poly(horsepower, 4)4 -5.1878 4.3727 -1.186 0.236 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 4.373 on 387 degrees of freedom ## Multiple R-squared: 0.6893, Adjusted R-squared: 0.6861 ## F-statistic: 214.7 on 4 and 387 DF, p-value: < 2.2e-16 #Get the range of horsepower hp <- range(auto$horsepower)
#Create a sequence to be used for plotting
hpGrid <- seq(hp[1],hp[2],by=10)
#Predict for these values of horsepower. Set Standard error as TRUE
pred=predict(fit,newdata=list(horsepower=hpGrid),se=TRUE)
#Compute bands on either side that is 2xSE
seBands=cbind(pred$fit+2*pred$se.fit,pred$fit-2*pred$se.fit)
#Plot the fit with Standard Error bands
plot(auto$horsepower,auto$mpg,xlim=hp,cex=.5,col="black",xlab="Horsepower",
ylab="MPG", main="Polynomial of degree 4")
lines(hpGrid,pred$fit,lwd=2,col="blue") matlines(hpGrid,seBands,lwd=2,col="blue",lty=3) 1.1b Fit a 4th degree polynomial – Python code import numpy as np import pandas as pd import os import matplotlib.pyplot as plt from sklearn.model_selection import train_test_split from sklearn.preprocessing import PolynomialFeatures from sklearn.linear_model import LinearRegression #Read the auto data autoDF =pd.read_csv("auto_mpg.csv",encoding="ISO-8859-1") # Select columns autoDF1=autoDF[['mpg','cylinder','displacement','horsepower','weight','acceleration','year']] # Convert all columns to numeric autoDF2 = autoDF1.apply(pd.to_numeric, errors='coerce') #Drop NAs autoDF3=autoDF2.dropna() autoDF3.shape X=autoDF3[['horsepower']] y=autoDF3['mpg'] #Create a polynomial of degree 4 poly = PolynomialFeatures(degree=4) X_poly = poly.fit_transform(X) # Fit a polynomial regression line linreg = LinearRegression().fit(X_poly, y) # Create a range of values hpGrid = np.arange(np.min(X),np.max(X),10) hp=hpGrid.reshape(-1,1) # Transform to 4th degree poly = PolynomialFeatures(degree=4) hp_poly = poly.fit_transform(hp) #Create a scatter plot plt.scatter(X,y) # Fit the prediction ypred=linreg.predict(hp_poly) plt.title("Poylnomial of degree 4") fig2=plt.xlabel("Horsepower") fig2=plt.ylabel("MPG") # Draw the regression curve plt.plot(hp,ypred,c="red") plt.savefig('fig1.png', bbox_inches='tight') 1.1c Fit a B-Spline – R Code In the code below a B- Spline is fit to data. The B-spline requires the manual selection of knots #Splines library(splines) # Fit a B-spline to the data. Select knots at 60,75,100,150 fit=lm(mpg~bs(horsepower,df=6,knots=c(60,75,100,150)),data=auto) # Use the fitted regresion to predict pred=predict(fit,newdata=list(horsepower=hpGrid),se=T) # Create a scatter plot plot(auto$horsepower,auto$mpg,xlim=hp,cex=.5,col="black",xlab="Horsepower", ylab="MPG", main="B-Spline with 4 knots") #Draw lines with 2 Standard Errors on either side lines(hpGrid,pred$fit,lwd=2)
lines(hpGrid,pred$fit+2*pred$se,lty="dashed")
lines(hpGrid,pred$fit-2*pred$se,lty="dashed")
abline(v=c(60,75,100,150),lty=2,col="darkgreen")

1.1d Fit a Natural Spline – R Code

Here a ‘Natural Spline’ is used to fit .The Natural Spline extrapolates beyond the boundary knots and the ends of the function are much more constrained than a regular spline or a global polynomoial where the ends can wag a lot more. Natural splines do not require the explicit selection of knots

# There is no need to select the knots here. There is a smoothing parameter which
# can be specified by the degrees of freedom 'df' parameter. The natural spline

fit2=lm(mpg~ns(horsepower,df=4),data=auto)
pred=predict(fit2,newdata=list(horsepower=hpGrid),se=T)
plot(auto$horsepower,auto$mpg,xlim=hp,cex=.5,col="black",xlab="Horsepower",
ylab="MPG", main="Natural Splines")
lines(hpGrid,pred$fit,lwd=2) lines(hpGrid,pred$fit+2*pred$se,lty="dashed") lines(hpGrid,pred$fit-2*pred$se,lty="dashed") 1.1.e Fit a Smoothing Spline – R code Here a smoothing spline is used. Smoothing splines also do not require the explicit setting of knots. We can change the ‘degrees of freedom(df)’ paramater to get the best fit # Smoothing spline has a smoothing parameter, the degrees of freedom # This is too wiggly plot(auto$horsepower,auto$mpg,xlim=hp,cex=.5,col="black",xlab="Horsepower", ylab="MPG", main="Smoothing Splines") # Here df is set to 16. This has a lot of variance fit=smooth.spline(auto$horsepower,auto$mpg,df=16) lines(fit,col="red",lwd=2) # We can use Cross Validation to allow the spline to pick the value of this smpopothing paramter. We do not need to set the degrees of freedom 'df' fit=smooth.spline(auto$horsepower,auto$mpg,cv=TRUE) lines(fit,col="blue",lwd=2) 1.1e Splines – Python There isn’t as much treatment of splines in Python and SKLearn. I did find the LSQUnivariate, UnivariateSpline spline. The LSQUnivariate spline requires the explcit setting of knots import numpy as np import pandas as pd import os import matplotlib.pyplot as plt from scipy.interpolate import LSQUnivariateSpline autoDF =pd.read_csv("auto_mpg.csv",encoding="ISO-8859-1") autoDF.shape autoDF.columns autoDF1=autoDF[['mpg','cylinder','displacement','horsepower','weight','acceleration','year']] autoDF2 = autoDF1.apply(pd.to_numeric, errors='coerce') auto=autoDF2.dropna() auto=auto[['horsepower','mpg']].sort_values('horsepower') # Set the knots manually knots=[65,75,100,150] # Create an array for X & y X=np.array(auto['horsepower']) y=np.array(auto['mpg']) # Fit a LSQunivariate spline s = LSQUnivariateSpline(X,y,knots) #Plot the spline xs = np.linspace(40,230,1000) ys = s(xs) plt.scatter(X, y) plt.plot(xs, ys) plt.savefig('fig2.png', bbox_inches='tight')  1.2 Generalized Additiive models (GAMs) Generalized Additive Models (GAMs) is a really powerful ML tool. $y_{i} = \beta_{0} + f_{1}(x_{i1}) + f_{2}(x_{i2}) + .. +f_{p}(x_{ip}) + \epsilon_{i}$ In GAMs we use a different functions for each of the variables. GAMs give a much better fit since we can choose any function for the different sections 1.2a Generalized Additive Models (GAMs) – R Code The plot below show the smooth spline that is fit for each of the features horsepower, cylinder, displacement, year and acceleration. We can use any function for example loess, 4rd order polynomial etc. library(gam) # Fit a smoothing spline for horsepower, cyliner, displacement and acceleration gam=gam(mpg~s(horsepower,4)+s(cylinder,5)+s(displacement,4)+s(year,4)+s(acceleration,5),data=auto) # Display the summary of the fit. This give the significance of each of the paramwetr # Also an ANOVA is given for each combination of the features summary(gam) ## ## Call: gam(formula = mpg ~ s(horsepower, 4) + s(cylinder, 5) + s(displacement, ## 4) + s(year, 4) + s(acceleration, 5), data = auto) ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -8.3190 -1.4436 -0.0261 1.2279 12.0873 ## ## (Dispersion Parameter for gaussian family taken to be 6.9943) ## ## Null Deviance: 23818.99 on 391 degrees of freedom ## Residual Deviance: 2587.881 on 370 degrees of freedom ## AIC: 1898.282 ## ## Number of Local Scoring Iterations: 3 ## ## Anova for Parametric Effects ## Df Sum Sq Mean Sq F value Pr(>F) ## s(horsepower, 4) 1 15632.8 15632.8 2235.085 < 2.2e-16 *** ## s(cylinder, 5) 1 508.2 508.2 72.666 3.958e-16 *** ## s(displacement, 4) 1 374.3 374.3 53.514 1.606e-12 *** ## s(year, 4) 1 2263.2 2263.2 323.583 < 2.2e-16 *** ## s(acceleration, 5) 1 372.4 372.4 53.246 1.809e-12 *** ## Residuals 370 2587.9 7.0 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Anova for Nonparametric Effects ## Npar Df Npar F Pr(F) ## (Intercept) ## s(horsepower, 4) 3 13.825 1.453e-08 *** ## s(cylinder, 5) 3 17.668 9.712e-11 *** ## s(displacement, 4) 3 44.573 < 2.2e-16 *** ## s(year, 4) 3 23.364 7.183e-14 *** ## s(acceleration, 5) 4 3.848 0.004453 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 par(mfrow=c(2,3)) plot(gam,se=TRUE) 1.2b Generalized Additive Models (GAMs) – Python Code I did not find the equivalent of GAMs in SKlearn in Python. There was an early prototype (2012) in Github. Looks like it is still work in progress or has probably been abandoned. 1.3 Tree based Machine Learning Models Tree based Machine Learning are all based on the ‘bootstrapping’ technique. In bootstrapping given a sample of size N, we create datasets of size N by sampling this original dataset with replacement. Machine Learning models are built on the different bootstrapped samples and then averaged. Decision Trees as seen above have the tendency to overfit. There are several techniques that help to avoid this namely a) Bagging b) Random Forests c) Boosting Bagging, Random Forest and Gradient Boosting Bagging: Bagging, or Bootstrap Aggregation decreases the variance of predictions, by creating separate Decisiion Tree based ML models on the different samples and then averaging these ML models Random Forests: Bagging is a greedy algorithm and tries to produce splits based on all variables which try to minimize the error. However the different ML models have a high correlation. Random Forests remove this shortcoming, by using a variable and random set of features to split on. Hence the features chosen and the resulting trees are uncorrelated. When these ML models are averaged the performance is much better. Boosting: Gradient Boosted Decision Trees also use an ensemble of trees but they don’t build Machine Learning models with random set of features at each step. Rather small and simple trees are built. Successive trees try to minimize the error from the earlier trees. Out of Bag (OOB) Error: In Random Forest and Gradient Boosting for each bootstrap sample taken from the dataset, there will be samples left out. These are known as Out of Bag samples.Classification accuracy carried out on these OOB samples is known as OOB error 1.31a Decision Trees – R Code The code below creates a Decision tree with the cancer training data. The summary of the fit is output. Based on the ML model, the predict function is used on test data and a confusion matrix is output. # Read the cancer data library(tree) library(caret) library(e1071) cancer <- read.csv("cancer.csv",stringsAsFactors = FALSE) cancer <- cancer[,2:32] cancer$target <- as.factor(cancer$target) train_idx <- trainTestSplit(cancer,trainPercent=75,seed=5) train <- cancer[train_idx, ] test <- cancer[-train_idx, ] # Create Decision Tree cancerStatus=tree(target~.,train) summary(cancerStatus) ## ## Classification tree: ## tree(formula = target ~ ., data = train) ## Variables actually used in tree construction: ## [1] "worst.perimeter" "worst.concave.points" "area.error" ## [4] "worst.texture" "mean.texture" "mean.concave.points" ## Number of terminal nodes: 9 ## Residual mean deviance: 0.1218 = 50.8 / 417 ## Misclassification error rate: 0.02347 = 10 / 426 pred <- predict(cancerStatus,newdata=test,type="class") confusionMatrix(pred,test$target)
## Confusion Matrix and Statistics
##
##           Reference
## Prediction  0  1
##          0 49  7
##          1  8 78
##
##                Accuracy : 0.8944
##                  95% CI : (0.8318, 0.9397)
##     No Information Rate : 0.5986
##     P-Value [Acc > NIR] : 4.641e-15
##
##                   Kappa : 0.7795
##  Mcnemar's Test P-Value : 1
##
##             Sensitivity : 0.8596
##             Specificity : 0.9176
##          Pos Pred Value : 0.8750
##          Neg Pred Value : 0.9070
##              Prevalence : 0.4014
##          Detection Rate : 0.3451
##    Detection Prevalence : 0.3944
##       Balanced Accuracy : 0.8886
##
##        'Positive' Class : 0
## 
# Plot decision tree with labels
plot(cancerStatus)
text(cancerStatus,pretty=0)

1.31b Decision Trees – Cross Validation – R Code

We can also perform a Cross Validation on the data to identify the Decision Tree which will give the minimum deviance.

library(tree)
cancer <- cancer[,2:32]
cancer$target <- as.factor(cancer$target)
train_idx <- trainTestSplit(cancer,trainPercent=75,seed=5)
train <- cancer[train_idx, ]
test <- cancer[-train_idx, ]

# Create Decision Tree
cancerStatus=tree(target~.,train)

# Execute 10 fold cross validation
cvCancer=cv.tree(cancerStatus)
plot(cvCancer)

# Plot the
plot(cvCancer$size,cvCancer$dev,type='b')

prunedCancer=prune.tree(cancerStatus,best=4)
plot(prunedCancer)
text(prunedCancer,pretty=0)

pred <- predict(prunedCancer,newdata=test,type="class")
confusionMatrix(pred,test$target) ## Confusion Matrix and Statistics ## ## Reference ## Prediction 0 1 ## 0 50 7 ## 1 7 78 ## ## Accuracy : 0.9014 ## 95% CI : (0.8401, 0.945) ## No Information Rate : 0.5986 ## P-Value [Acc > NIR] : 7.988e-16 ## ## Kappa : 0.7948 ## Mcnemar's Test P-Value : 1 ## ## Sensitivity : 0.8772 ## Specificity : 0.9176 ## Pos Pred Value : 0.8772 ## Neg Pred Value : 0.9176 ## Prevalence : 0.4014 ## Detection Rate : 0.3521 ## Detection Prevalence : 0.4014 ## Balanced Accuracy : 0.8974 ## ## 'Positive' Class : 0 ##  1.31c Decision Trees – Python Code Below is the Python code for creating Decision Trees. The accuracy, precision, recall and F1 score is computed on the test data set. import numpy as np import pandas as pd import os import matplotlib.pyplot as plt from sklearn.metrics import confusion_matrix from sklearn import tree from sklearn.datasets import load_breast_cancer from sklearn.model_selection import train_test_split from sklearn.tree import DecisionTreeClassifier from sklearn.datasets import make_classification, make_blobs from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score import graphviz cancer = load_breast_cancer() (X_cancer, y_cancer) = load_breast_cancer(return_X_y = True) X_train, X_test, y_train, y_test = train_test_split(X_cancer, y_cancer, random_state = 0) clf = DecisionTreeClassifier().fit(X_train, y_train) print('Accuracy of Decision Tree classifier on training set: {:.2f}' .format(clf.score(X_train, y_train))) print('Accuracy of Decision Tree classifier on test set: {:.2f}' .format(clf.score(X_test, y_test))) y_predicted=clf.predict(X_test) confusion = confusion_matrix(y_test, y_predicted) print('Accuracy: {:.2f}'.format(accuracy_score(y_test, y_predicted))) print('Precision: {:.2f}'.format(precision_score(y_test, y_predicted))) print('Recall: {:.2f}'.format(recall_score(y_test, y_predicted))) print('F1: {:.2f}'.format(f1_score(y_test, y_predicted))) # Plot the Decision Tree clf = DecisionTreeClassifier(max_depth=2).fit(X_train, y_train) dot_data = tree.export_graphviz(clf, out_file=None, feature_names=cancer.feature_names, class_names=cancer.target_names, filled=True, rounded=True, special_characters=True) graph = graphviz.Source(dot_data) graph ## Accuracy of Decision Tree classifier on training set: 1.00 ## Accuracy of Decision Tree classifier on test set: 0.87 ## Accuracy: 0.87 ## Precision: 0.97 ## Recall: 0.82 ## F1: 0.89 1.31d Decision Trees – Cross Validation – Python Code In the code below 5-fold cross validation is performed for different depths of the tree and the accuracy is computed. The accuracy on the test set seems to plateau when the depth is 8. But it is seen to increase again from 10 to 12. More analysis needs to be done here  import numpy as np import pandas as pd import os import matplotlib.pyplot as plt from sklearn.datasets import load_breast_cancer from sklearn.tree import DecisionTreeClassifier (X_cancer, y_cancer) = load_breast_cancer(return_X_y = True) from sklearn.cross_validation import train_test_split, KFold def computeCVAccuracy(X,y,folds): accuracy=[] foldAcc=[] depth=[1,2,3,4,5,6,7,8,9,10,11,12] nK=len(X)/float(folds) xval_err=0 for i in depth: kf = KFold(len(X),n_folds=folds) for train_index, test_index in kf: X_train, X_test = X.iloc[train_index], X.iloc[test_index] y_train, y_test = y.iloc[train_index], y.iloc[test_index] clf = DecisionTreeClassifier(max_depth = i).fit(X_train, y_train) score=clf.score(X_test, y_test) accuracy.append(score) foldAcc.append(np.mean(accuracy)) return(foldAcc) cvAccuracy=computeCVAccuracy(pd.DataFrame(X_cancer),pd.DataFrame(y_cancer),folds=10) df1=pd.DataFrame(cvAccuracy) df1.columns=['cvAccuracy'] df=df1.reindex([1,2,3,4,5,6,7,8,9,10,11,12]) df.plot() plt.title("Decision Tree - 10-fold Cross Validation Accuracy vs Depth of tree") plt.xlabel("Depth of tree") plt.ylabel("Accuracy") plt.savefig('fig3.png', bbox_inches='tight') 1.4a Random Forest – R code A Random Forest is fit using the Boston data. The summary shows that 4 variables were randomly chosen at each split and the resulting ML model explains 88.72% of the test data. Also the variable importance is plotted. It can be seen that ‘rooms’ and ‘status’ are the most influential features in the model library(randomForest) df=read.csv("Boston.csv",stringsAsFactors = FALSE) # Data from MASS - SL # Select specific columns Boston <- df %>% dplyr::select("crimeRate","zone","indus","charles","nox","rooms","age", "distances","highways","tax","teacherRatio","color", "status","medianValue") # Fit a Random Forest on the Boston training data rfBoston=randomForest(medianValue~.,data=Boston) # Display the summatu of the fit. It can be seen that the MSE is 10.88 # and the percentage variance explained is 86.14%. About 4 variables were tried at each # #split for a maximum tree of 500. # The MSE and percent variance is on Out of Bag trees rfBoston ## ## Call: ## randomForest(formula = medianValue ~ ., data = Boston) ## Type of random forest: regression ## Number of trees: 500 ## No. of variables tried at each split: 4 ## ## Mean of squared residuals: 9.521672 ## % Var explained: 88.72 #List and plot the variable importances importance(rfBoston) ## IncNodePurity ## crimeRate 2602.1550 ## zone 258.8057 ## indus 2599.6635 ## charles 240.2879 ## nox 2748.8485 ## rooms 12011.6178 ## age 1083.3242 ## distances 2432.8962 ## highways 393.5599 ## tax 1348.6987 ## teacherRatio 2841.5151 ## color 731.4387 ## status 12735.4046 varImpPlot(rfBoston) 1.4b Random Forest-OOB and Cross Validation Error – R code The figure below shows the OOB error and the Cross Validation error vs the ‘mtry’. Here mtry indicates the number of random features that are chosen at each split. The lowest test error occurs when mtry = 8 library(randomForest) df=read.csv("Boston.csv",stringsAsFactors = FALSE) # Data from MASS - SL # Select specific columns Boston <- df %>% dplyr::select("crimeRate","zone","indus","charles","nox","rooms","age", "distances","highways","tax","teacherRatio","color", "status","medianValue") # Split as training and tst sets train_idx <- trainTestSplit(Boston,trainPercent=75,seed=5) train <- Boston[train_idx, ] test <- Boston[-train_idx, ] #Initialize OOD and testError oobError <- NULL testError <- NULL # In the code below the number of variables to consider at each split is increased # from 1 - 13(max features) and the OOB error and the MSE is computed for(i in 1:13){ fitRF=randomForest(medianValue~.,data=train,mtry=i,ntree=400) oobError[i] <-fitRF$mse[400]
pred <- predict(fitRF,newdata=test)
testError[i] <- mean((pred-test$medianValue)^2) } # We can see the OOB and Test Error. It can be seen that the Random Forest performs # best with the lowers MSE at mtry=6 matplot(1:13,cbind(testError,oobError),pch=19,col=c("red","blue"), type="b",xlab="mtry(no of varaibles at each split)", ylab="Mean Squared Error", main="Random Forest - OOB and Test Error") legend("topright",legend=c("OOB","Test"),pch=19,col=c("red","blue")) 1.4c Random Forest – Python code The python code for Random Forest Regression is shown below. The training and test score is computed. The variable importance shows that ‘rooms’ and ‘status’ are the most influential of the variables import numpy as np import pandas as pd import os import matplotlib.pyplot as plt from sklearn.model_selection import train_test_split from sklearn.ensemble import RandomForestRegressor df = pd.read_csv("Boston.csv",encoding = "ISO-8859-1") X=df[['crimeRate','zone', 'indus','charles','nox','rooms', 'age','distances','highways','tax', 'teacherRatio','color','status']] y=df['medianValue'] X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = 0) regr = RandomForestRegressor(max_depth=4, random_state=0) regr.fit(X_train, y_train) print('R-squared score (training): {:.3f}' .format(regr.score(X_train, y_train))) print('R-squared score (test): {:.3f}' .format(regr.score(X_test, y_test))) feature_names=['crimeRate','zone', 'indus','charles','nox','rooms', 'age','distances','highways','tax', 'teacherRatio','color','status'] print(regr.feature_importances_) plt.figure(figsize=(10,6),dpi=80) c_features=X_train.shape[1] plt.barh(np.arange(c_features),regr.feature_importances_) plt.xlabel("Feature importance") plt.ylabel("Feature name") plt.yticks(np.arange(c_features), feature_names) plt.tight_layout() plt.savefig('fig4.png', bbox_inches='tight')  ## R-squared score (training): 0.917 ## R-squared score (test): 0.734 ## [ 0.03437382 0. 0.00580335 0. 0.00731004 0.36461548 ## 0.00638577 0.03432173 0.0041244 0.01732328 0.01074148 0.0012638 ## 0.51373683] 1.4d Random Forest – Cross Validation and OOB Error – Python code As with R the ‘max_features’ determines the random number of features the random forest will use at each split. The plot shows that when max_features=8 the MSE is lowest import numpy as np import pandas as pd import os import matplotlib.pyplot as plt from sklearn.model_selection import train_test_split from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import cross_val_score df = pd.read_csv("Boston.csv",encoding = "ISO-8859-1") X=df[['crimeRate','zone', 'indus','charles','nox','rooms', 'age','distances','highways','tax', 'teacherRatio','color','status']] y=df['medianValue'] cvError=[] oobError=[] oobMSE=[] for i in range(1,13): regr = RandomForestRegressor(max_depth=4, n_estimators=400,max_features=i,oob_score=True,random_state=0) mse= np.mean(cross_val_score(regr, X, y, cv=5,scoring = 'neg_mean_squared_error')) # Since this is neg_mean_squared_error I have inverted the sign to get MSE cvError.append(-mse) # Fit on all data to compute OOB error regr.fit(X, y) # Record the OOB error for each max_features=i setting oob = 1 - regr.oob_score_ oobError.append(oob) # Get the Out of Bag prediction oobPred=regr.oob_prediction_ # Compute the Mean Squared Error between OOB Prediction and target mseOOB=np.mean(np.square(oobPred-y)) oobMSE.append(mseOOB) # Plot the CV Error and OOB Error # Set max_features maxFeatures=np.arange(1,13) cvError=pd.DataFrame(cvError,index=maxFeatures) oobMSE=pd.DataFrame(oobMSE,index=maxFeatures) #Plot fig8=df.plot() fig8=plt.title('Random forest - CV Error and OOB Error vs max_features') fig8.figure.savefig('fig8.png', bbox_inches='tight')  #Plot the OOB Error vs max_features plt.plot(range(1,13),oobError) fig2=plt.title("Random Forest - OOB Error vs max_features (variable no of features)") fig2=plt.xlabel("max_features (variable no of features)") fig2=plt.ylabel("OOB Error") fig2.figure.savefig('fig7.png', bbox_inches='tight')  1.5a Boosting – R code Here a Gradient Boosted ML Model is built with a n.trees=5000, with a learning rate of 0.01 and depth of 4. The feature importance plot also shows that rooms and status are the 2 most important features. The MSE vs the number of trees plateaus around 2000 trees library(gbm) # Perform gradient boosting on the Boston data set. The distribution is gaussian since we # doing MSE. The interaction depth specifies the number of splits boostBoston=gbm(medianValue~.,data=train,distribution="gaussian",n.trees=5000, shrinkage=0.01,interaction.depth=4) #The summary gives the variable importance. The 2 most significant variables are # number of rooms and lower status summary(boostBoston) ## var rel.inf ## rooms rooms 42.2267200 ## status status 27.3024671 ## distances distances 7.9447972 ## crimeRate crimeRate 5.0238827 ## nox nox 4.0616548 ## teacherRatio teacherRatio 3.1991999 ## age age 2.7909772 ## color color 2.3436295 ## tax tax 2.1386213 ## charles charles 1.3799109 ## highways highways 0.7644026 ## indus indus 0.7236082 ## zone zone 0.1001287 # The plots below show how each variable relates to the median value of the home. As # the number of roomd increase the median value increases and with increase in lower status # the median value decreases par(mfrow=c(1,2)) #Plot the relation between the top 2 features and the target plot(boostBoston,i="rooms") plot(boostBoston,i="status") # Create a sequence of trees between 100-5000 incremented by 50 nTrees=seq(100,5000,by=50) # Predict the values for the test data pred <- predict(boostBoston,newdata=test,n.trees=nTrees) # Compute the mean for each of the MSE for each of the number of trees boostError <- apply((pred-test$medianValue)^2,2,mean)
#Plot the MSE vs the number of trees
plot(nTrees,boostError,pch=19,col="blue",ylab="Mean Squared Error",
main="Boosting Test Error")

1.5b Cross Validation Boosting – R code

Included below is a cross validation error vs the learning rate. The lowest error is when learning rate = 0.09

cvError <- NULL
s <- c(.001,0.01,0.03,0.05,0.07,0.09,0.1)
for(i in seq_along(s)){
cvBoost=gbm(medianValue~.,data=train,distribution="gaussian",n.trees=5000,
shrinkage=s[i],interaction.depth=4,cv.folds=5)
cvError[i] <- mean(cvBoost\$cv.error)
}

# Create a data frame for plotting
a <- rbind(s,cvError)
b <- as.data.frame(t(a))
# It can be seen that a shrinkage parameter of 0,05 gives the lowes CV Error
ggplot(b,aes(s,cvError)) + geom_point() + geom_line(color="blue") +
xlab("Shrinkage") + ylab("Cross Validation Error") +
ggtitle("Gradient boosted trees - Cross Validation error vs Shrinkage")

1.5c Boosting – Python code

A gradient boost ML model in Python is created below. The Rsquared score is computed on the training and test data.

import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

X=df[['crimeRate','zone', 'indus','charles','nox','rooms', 'age','distances','highways','tax',
'teacherRatio','color','status']]
y=df['medianValue']

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = 0)

regr.fit(X_train, y_train)

print('R-squared score (training): {:.3f}'
.format(regr.score(X_train, y_train)))
print('R-squared score (test): {:.3f}'
.format(regr.score(X_test, y_test)))
## R-squared score (training): 0.983
## R-squared score (test): 0.821

1.5c Cross Validation Boosting – Python code

the cross validation error is computed as the learning rate is varied. The minimum CV eror occurs when lr = 0.04

import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score

X=df[['crimeRate','zone', 'indus','charles','nox','rooms', 'age','distances','highways','tax',
'teacherRatio','color','status']]
y=df['medianValue']

cvError=[]
learning_rate =[.001,0.01,0.03,0.05,0.07,0.09,0.1]
for lr in learning_rate:
mse= np.mean(cross_val_score(regr, X, y, cv=10,scoring = 'neg_mean_squared_error'))
# Since this is neg_mean_squared_error I have inverted the sign to get MSE
cvError.append(-mse)
learning_rate =[.001,0.01,0.03,0.05,0.07,0.09,0.1]
plt.plot(learning_rate,cvError)
plt.title("Gradient Boosting - 5-fold CV- Mean Squared Error vs max_features (variable no of features)")
plt.xlabel("max_features (variable no of features)")
plt.ylabel("Mean Squared Error")
plt.savefig('fig6.png', bbox_inches='tight')

Conclusion This post covered Splines and Tree based ML models like Bagging, Random Forest and Boosting. Stay tuned for further updates.

You may also like

To see all posts see Index of posts