Good Articles to learn how to implement a neural network 2

March 30, 2017 Author: david

Vectorization

This part will cover:

The previous tutorial described a very simple neural network with only one input, one hidden neuron and one output. This tutorial will describe a neural network that takes 2-dimensional input samples, projects them onto a 3-dimensional hidden layer, and classifies them with a 2-dimensional softmax output classfier, this softmax function is explained in intermezzo 2 . While we didn’t add the bias parameters to the previous 2 models, we will add them to this model. We will work out the vector computations needed to optimize and run this neural network in this tutorial. The network of this model is shown in the following figure:

The notebook starts out with importing the libraries we need:

In [1]:

# Python imports
import numpy as np  # Matrix and vector computation package
import sklearn.datasets # To generate the dataset
import matplotlib.pyplot as plt  # Plotting library
from matplotlib.colors import colorConverter, ListedColormap  # Some plotting functions
from mpl_toolkits.mplot3d import Axes3D  # 3D plots
from matplotlib import cm  # Colormaps
# Allow matplotlib to plot inside this notebook
%matplotlib inline
# Set the seed of the numpy random number generator so that the tutorial is reproducable
np.random.seed(seed=1)

# Python imports

import numpy as np # Matrix and vector computation package

import sklearn.datasets # To generate the dataset

import matplotlib.pyplot as plt # Plotting library

from matplotlib.colors import colorConverter, ListedColormap # Some plotting functions

from mpl_toolkits.mplot3d import Axes3D # 3D plots

from matplotlib import cm # Colormaps

# Allow matplotlib to plot inside this notebook

%matplotlib inline

# Set the seed of the numpy random number generator so that the tutorial is reproducable

np.random.seed(seed=1)

Define the dataset

In this example the target classes $t$

corresponding to the inputs $x$ will be generated from 2 class distributions: red ( $t = 0$ ) and blue ( $t = 1$

). Where the red class is a circular distribution that surrounds the distribution of the blue class. The dataset is generated by the scikit-learn make_circles method.

This results in a 2D dataset that is not linearly separable. The model from part 2 won’t be able to classify both classes correctly since it can learn only linear separators. By adding a hidden layer, the model will be able to train a non-linear separator.

The $N$

input samples with $2$ variables each are given as a $N \times 2$ matrix $X$

X = [\begin{matrix} x_{11} & x_{12} \\ ⋮ & ⋮ \\ x_{N 1} & x_{N 2} \end{matrix}]

Where $x_{i j}$

is the value of the $j$ -th variable of the $i$

-th input sample.

Since the softmax output function is used the targets for each input sample are given as a $N \times 2$

matrix $T$

T = [\begin{matrix} t_{11} & t_{12} \\ ⋮ & ⋮ \\ t_{N 1} & t_{N 2} \end{matrix}]

Where $t_{i j} = 1$

if and only if the $i$ -th input sample belongs to class $j$

. So blue points are labelled T = [0 1] and red points are labelled T = [1 0].

In [2]:

# Generate the dataset
X, t = sklearn.datasets.make_circles(n_samples=100, shuffle=False, factor=0.3, noise=0.1)
T = np.zeros((100,2)) # Define target matrix
T[t==1,1] = 1
T[t==0,0] = 1
# Separate the red and blue points for plotting
x_red = X[t==0]
x_blue = X[t==1]
print('shape of X: {}'.format(X.shape))
print('shape of T: {}'.format(T.shape))

# Generate the dataset

X, t = sklearn.datasets.make_circles(n_samples=100, shuffle=False, factor=0.3, noise=0.1)

T = np.zeros((100,2)) # Define target matrix

T[t==1,1] = 1

T[t==0,0] = 1

# Separate the red and blue points for plotting

x_red = X[t==0]

x_blue = X[t==1]

print('shape of X: {}'.format(X.shape))

print('shape of T: {}'.format(T.shape))

shape of X: (100, 2)
shape of T: (100, 2)

shape of X: (100, 2)

shape of T: (100, 2)

In [3]:

# Plot both classes on the x1, x2 plane
plt.plot(x_red[:,0], x_red[:,1], 'ro', label='class red')
plt.plot(x_blue[:,0], x_blue[:,1], 'bo', label='class blue')
plt.grid()
plt.legend(loc=1)
plt.xlabel('$x_1$', fontsize=15)
plt.ylabel('$x_2$', fontsize=15)
plt.axis([-1.5, 1.5, -1.5, 1.5])
plt.title('red vs blue classes in the input space')
plt.show()

# Plot both classes on the x1, x2 plane

plt.plot(x_red[:,0], x_red[:,1], 'ro', label='class red')

plt.plot(x_blue[:,0], x_blue[:,1], 'bo', label='class blue')

plt.grid()

plt.legend(loc=1)

plt.xlabel('$x_1$', fontsize=15)

plt.ylabel('$x_2$', fontsize=15)

plt.axis([-1.5, 1.5, -1.5, 1.5])

plt.title('red vs blue classes in the input space')

plt.show()

# Plot both classes on the x1, x2 plane
plt.plot(x_red[:,0], x_red[:,1], 'ro', label='class red')
plt.plot(x_blue[:,0], x_blue[:,1], 'bo', label='class blue')
plt.grid()
plt.legend(loc=1)
plt.xlabel('$x_1$', fontsize=15)
plt.ylabel('$x_2$', fontsize=15)
plt.axis([-1.5, 1.5, -1.5, 1.5])
plt.title('red vs blue classes in the input space')
plt.show()

# Plot both classes on the x1, x2 plane

plt.plot(x_red[:,0], x_red[:,1], 'ro', label='class red')

plt.plot(x_blue[:,0], x_blue[:,1], 'bo', label='class blue')

plt.grid()

plt.legend(loc=1)

plt.xlabel('$x_1$', fontsize=15)

plt.ylabel('$x_2$', fontsize=15)

plt.axis([-1.5, 1.5, -1.5, 1.5])

plt.title('red vs blue classes in the input space')

plt.show()

Vectorization of backpropagation

1. Vectorization of the forward step

Compute activations of hidden layer

The 2-dimensional inputs $X$

are projected onto the 3 dimensions of the hidden layer $H$ by the weight matrix $W_{h}$ ( $w_{h i j}$ is the weight of the connection between input variable $i$ and hidden neuron activation $j$ ) and bias vector $b_{h}$

\begin{aligned} W_{h} = [\begin{matrix} w_{h 11} & w_{h 12} & w_{h 13} \\ w_{h 21} & w_{h 22} & w_{h 23} \end{matrix}] & b_{h} = [\begin{matrix} b_{h 1} & b_{h 2} & b_{h 3} \end{matrix}] \end{aligned}

following computation:

H = σ (X \cdot W_{h} + b_{h}) = \frac{1}{1 + e^{- (X \cdot W_{h} + b_{h})}} = [\begin{matrix} h_{11} & h_{12} & h_{13} \\ ⋮ & ⋮ & ⋮ \\ h_{N 1} & h_{N 2} & h_{N 3} \end{matrix}]

With $σ$

the logistic function, and with $H$ resulting in a $N \times 3$

matrix.

This computation flow is illustrated by the figure below. Each input $x_{i j}$

is multiplied with the weight parameters from the corresponding column $w_{h j 1}, w_{h j 2}, w_{h j 3}$ . The multiplication results for each row $k$ ( $x_{i j} * w_{h j k}$ ) are summed to form the $z_{i k}$ , after which these are transformed by the logistic function $σ$ to the $h_{i k}$ . Note that the bias $B_{h}$ can be incorporated in $W_{h}$ by adding a new variable to each input $x_{i}$ that is always $+ 1$

$b_{h}$

and $W_{h}$

are represented in the code below by bh and Wh respectively. The hidden layer activations are computed by the hidden_activations(X, Wh, bh) method.

Compute activations of output

To compute the output activations the hidden layer activations can be projected onto the 2-dimensional output layer. This is done by the $3 \times 2$

weight matrix $W_{o}$ ( $w_{o i j}$ is the weight of the connection between hidden layer neuron $i$ and output activation $j$ ) and $2 \times 1$ bias vector $b_{o}$

\begin{aligned} W_{o} = [\begin{matrix} w_{o 11} & w_{o 12} \\ w_{o 21} & w_{o 22} \\ w_{o 31} & w_{o 32} \end{matrix}] & b_{o} = [\begin{matrix} b_{o 1} & b_{o 2} \end{matrix}] \end{aligned}

following computation:

Y = ς (H \cdot W_{o} + b_{o}) = \frac{e^{Z_{o}}}{\sum_{d = 1}^{C} e^{z_{o d}}} = \frac{e^{H \cdot W_{o} + b_{o}}}{\sum_{d = 1}^{C} e^{H \cdot w_{o d} + b_{o d}}} = [\begin{matrix} y_{11} & y_{12} \\ ⋮ & ⋮ \\ y_{n 1} & y_{n 2} \end{matrix}]

With $ς$

the softmax function, $Y$ resulting in a $n \times 2$ matrix, $z_{o d}$ the $d$ -th column of the $Z_{o}$ matrix, $w_{o d}$ the $d$ -th column of the $W_{o}$ matrix, and $b_{o d}$ the $d$ -th element of the $b_{o}$

vector.

$b_{o}$

and $W_{o}$

are represented in the code below by bo and Wo respectively. The hidden layer activations are computed by the output_activations(H, Wo, bo) method.

In [4]:

# Define the logistic function
def logistic(z): 
    return 1 / (1 + np.exp(-z))
# Define the softmax function
def softmax(z): 
    return np.exp(z) / np.sum(np.exp(z), axis=1, keepdims=True)
# Function to compute the hidden activations
def hidden_activations(X, Wh, bh):
    return logistic(X.dot(Wh) + bh)
# Define output layer feedforward
def output_activations(H, Wo, bo):
    return softmax(H.dot(Wo) + bo)
# Define the neural network function
def nn(X, Wh, bh, Wo, bo): 
    return output_activations(hidden_activations(X, Wh, bh), Wo, bo)
# Define the neural network prediction function that only returns
#  1 or 0 depending on the predicted class
def nn_predict(X, Wh, bh, Wo, bo): 
    return np.around(nn(X, Wh, bh, Wo, bo))

# Define the logistic function

def logistic(z):

return 1 / (1 + np.exp(-z))

# Define the softmax function

def softmax(z):

return np.exp(z) / np.sum(np.exp(z), axis=1, keepdims=True)

# Function to compute the hidden activations

def hidden_activations(X, Wh, bh):

return logistic(X.dot(Wh) + bh)

# Define output layer feedforward

def output_activations(H, Wo, bo):

return softmax(H.dot(Wo) + bo)

# Define the neural network function

def nn(X, Wh, bh, Wo, bo):

return output_activations(hidden_activations(X, Wh, bh), Wo, bo)

# Define the neural network prediction function that only returns

# 1 or 0 depending on the predicted class

def nn_predict(X, Wh, bh, Wo, bo):

return np.around(nn(X, Wh, bh, Wo, bo))

# Define the logistic function
def logistic(z): 
    return 1 / (1 + np.exp(-z))
# Define the softmax function
def softmax(z): 
    return np.exp(z) / np.sum(np.exp(z), axis=1, keepdims=True)
# Function to compute the hidden activations
def hidden_activations(X, Wh, bh):
    return logistic(X.dot(Wh) + bh)
# Define output layer feedforward
def output_activations(H, Wo, bo):
    return softmax(H.dot(Wo) + bo)
# Define the neural network function
def nn(X, Wh, bh, Wo, bo): 
    return output_activations(hidden_activations(X, Wh, bh), Wo, bo)
# Define the neural network prediction function that only returns
#  1 or 0 depending on the predicted class
def nn_predict(X, Wh, bh, Wo, bo): 
    return np.around(nn(X, Wh, bh, Wo, bo))

# Define the logistic function

def logistic(z):

return 1 / (1 + np.exp(-z))

# Define the softmax function

def softmax(z):

return np.exp(z) / np.sum(np.exp(z), axis=1, keepdims=True)

# Function to compute the hidden activations

def hidden_activations(X, Wh, bh):

return logistic(X.dot(Wh) + bh)

# Define output layer feedforward

def output_activations(H, Wo, bo):

return softmax(H.dot(Wo) + bo)

# Define the neural network function

def nn(X, Wh, bh, Wo, bo):

return output_activations(hidden_activations(X, Wh, bh), Wo, bo)

# Define the neural network prediction function that only returns

# 1 or 0 depending on the predicted class

def nn_predict(X, Wh, bh, Wo, bo):

return np.around(nn(X, Wh, bh, Wo, bo))

2. Vectorization of the backward step

Vectorization of the output layer backward step

Compute the error at the output

The output layer is a softmax layer with corresponding cross-entropy cost function, which are described in detail in intermezzo 2 . The cost function $ξ$

for $N$ samples and $C$

classes used in this model is defined as:

ξ (T, Y) = \sum_{i = 1}^{n} ξ (t_{i}, y_{i}) = - \sum_{i = 1}^{n} \sum_{i = c}^{C} t_{i c} \cdot l o g (y_{i c})

The error gradient $δ_{o}$

of this cost function at the softmax output layer is simply:

δ_{o} = \frac{\partial ξ}{\partial Z_{o}} = Y - T

With $Z_{o}$

a $n \times 2$ matrix of inputs to the softmax layer ( $Z_{o} = H \cdot W_{o} + b_{o}$ ), and $T$ a $n \times 2$ target matrix that corresponds to $Y$ . Note that $δ_{o}$ also results in a $n \times 2$

matrix.

$δ_{o}$

will be defined in the code as Eo and will be computed by the error_output(Y, T) method defined below.

Update the output layer weights

At the output the gradient $δ_{w_{o j}}$

over all $N$ samples is computed by $\partial ξ / \partial w_{o j}$

\frac{\partial ξ}{\partial w_{o j}} = \frac{\partial Z_{o}}{\partial w_{o j}} \frac{\partial Y}{\partial Z_{o}} \frac{\partial ξ}{\partial Y} = \frac{\partial Z_{o}}{\partial w_{o j i}} \frac{\partial ξ}{\partial Z_{o}} = \sum_{i = 1}^{N} h_{i j} (y_{i} - t_{i}) = \sum_{i = 1}^{N} h_{i j} δ_{o i}

With $w_{o j}$

the $j$ -th row of $W_{o}$ and thus a $1 \times 2$

vector. This formula can be written out as matrix operations with the parameters of the output layer:

\frac{\partial ξ}{\partial W_{o}} = H^{T} \cdot (Y - T) = H^{T} \cdot δ_{o}

The resulting gradient is a $3 \times 2$

Jacobian matrix :

J_{W_{o}} = \frac{\partial ξ}{\partial W_{o}} = [\begin{matrix} \frac{\partial ξ}{\partial w_{o 11}} & \frac{\partial ξ}{\partial w_{o 12}} \\ \frac{\partial ξ}{\partial w_{o 21}} & \frac{\partial ξ}{\partial w_{o 22}} \\ \frac{\partial ξ}{\partial w_{o 31}} & \frac{\partial ξ}{\partial w_{o 32}} \end{matrix}]

$J_{W_{o}}$

will be defined in the code as JWo and will be computed by the gradient_weight_out(H, Eo) method defined below.

Update the output layer bias

The bias $b_{o}$

can be updated in the same manner. The formula to compute the gradient $\partial ξ / \partial b_{o}$ for batch processing over all $N$

samples is:

\frac{\partial ξ}{\partial b_{o}} = \frac{\partial Z_{o}}{\partial b_{o}} \frac{\partial Y}{\partial Z_{o}} \frac{\partial ξ}{\partial Y} = \sum_{i = 1}^{N} 1 * (y_{i} - t_{i}) = \sum_{i = 1}^{N} δ_{o i}

The resulting gradient is a $2 \times 1$

Jacobian matrix:

J_{b_{o}} = \frac{\partial ξ}{\partial b_{o}} = [\begin{matrix} \frac{\partial ξ}{\partial b_{o 1}} & \frac{\partial ξ}{\partial b_{o 2}} \end{matrix}]

$J_{b_{o}}$

will be defined in the code as Jbo and will be computed by the gradient_bias_out(Eo) method defined below.

In [5]:

# Define the cost function
def cost(Y, T):
    return - np.multiply(T, np.log(Y)).sum()
# Define the error function at the output
def error_output(Y, T):
    return Y - T
# Define the gradient function for the weight parameters at the output layer
def gradient_weight_out(H, Eo): 
    return  H.T.dot(Eo)
# Define the gradient function for the bias parameters at the output layer
def gradient_bias_out(Eo): 
    return  np.sum(Eo, axis=0, keepdims=True)

# Define the cost function

def cost(Y, T):

return - np.multiply(T, np.log(Y)).sum()

# Define the error function at the output

def error_output(Y, T):

return Y - T

# Define the gradient function for the weight parameters at the output layer

def gradient_weight_out(H, Eo):

return H.T.dot(Eo)

# Define the gradient function for the bias parameters at the output layer

def gradient_bias_out(Eo):

return np.sum(Eo, axis=0, keepdims=True)

Vectorization of the hidden layer backward step

Compute the error at the hidden layer

The error gradient $δ_{h}$

of the cost function at the hidden layer is defined as:

δ_{h} = \frac{\partial ξ}{\partial Z_{h}} = \frac{\partial H}{\partial Z_{h}} \frac{\partial ξ}{\partial H} = \frac{\partial H}{\partial Z_{h}} \frac{\partial Z_{o}}{\partial H} \frac{\partial ξ}{\partial Z_{o}}

With $Z_{h}$

a $n \times 3$ matrix of inputs to the logistic functions in the hidden neurons ( $Z_{h} = X \cdot W_{h} + b_{h}$ ). Note that $δ_{h}$ will also result in a $N \times 3$

matrix.

Lets first derive the error gradient $δ_{h i j}$

for one sample $i$ at hidden neuron $j$ . The gradients that backpropagate from the previous layer via the weighted connections are summed for each origin $h_{i j}$

δ_{h i j} = \frac{\partial ξ}{\partial z_{h i j}} = \frac{\partial h_{i j}}{\partial z_{h i j}} \frac{\partial z_{o i}}{\partial h_{i j}} \frac{\partial ξ}{\partial z_{o i}} = h_{i j} (1 - h_{i j}) \sum_{k = 1}^{2} w_{o j k} (y_{i k} - t_{i k}) = h_{i j} (1 - h_{i j}) [δ_{o i} \cdot w_{o j}^{T}]

Where $w_{o j}$

is the $j$ -th row of $W_{o}$ , and thus a $1 \times 2$ vector and $δ_{o i}$ is a $1 \times 2$ vector. The full $N \times 3$ error matrix $δ_{h}$

can thus be calculated as:

δ_{h} = \frac{\partial ξ}{\partial Z_{h}} = H \circ (1 - H) \circ [δ_{o} \cdot W_{o}^{T}]

With $\circ$

the elementwise product .

$δ_{h}$

will be defined in the code as Eh and will be computed by the error_hidden(H, Wo, Eo) method defined below.

Update the hidden layer weights

At the hidden layer the gradient $\partial ξ / \partial w_{h j}$

over all $N$

samples can be computed by:

\frac{\partial ξ}{\partial w_{h j}} = \frac{\partial Z_{h}}{\partial w_{h j}} \frac{\partial H}{\partial Z_{h}} \frac{\partial ξ}{\partial H} = \frac{\partial Z_{h}}{\partial w_{h j}} \frac{\partial ξ}{\partial Z_{h}} = \sum_{i = 1}^{N} x_{i j} δ_{h i}

With $w_{h j}$

the $j$ -th row of $W_{h}$ and thus a $1 \times 3$

vector. We can write this formula out as matrix operations with the parameters of the hidden layer:

\frac{\partial ξ}{\partial W_{h}} = X^{T} \cdot δ_{h}

The resulting gradient is a $2 \times 3$

Jacobian matrix :

J_{W_{h}} = \frac{\partial ξ}{\partial W_{h}} = [\begin{matrix} \frac{\partial ξ}{\partial w_{h 11}} & \frac{\partial ξ}{\partial w_{h 12}} & \frac{\partial ξ}{\partial w_{h 13}} \\ \frac{\partial ξ}{\partial w_{h 21}} & \frac{\partial ξ}{\partial w_{h 22}} & \frac{\partial ξ}{\partial w_{h 23}} \end{matrix}]

$J_{W_{h}}$

will be defined in the code as JWh and will be computed by the gradient_weight_hidden(X, Eh) method defined below.

Update the hidden layer bias

The bias $b_{h}$

can be updated in the same manner. The formula to compute the gradient $\partial ξ / \partial b_{h}$ for batch processing over all $N$

samples is:

\frac{\partial ξ}{\partial b_{h}} = \frac{\partial Z_{h}}{\partial b_{h}} \frac{\partial H}{\partial Z_{h}} \frac{\partial ξ}{\partial H} = \sum_{j = 1}^{N} δ_{h j}

The resulting gradient is a $1 \times 3$

Jacobian matrix:

J_{b_{h}} = \frac{\partial ξ}{\partial b_{h}} = [\begin{matrix} \frac{\partial ξ}{\partial b_{h 1}} & \frac{\partial ξ}{\partial b_{h 2}} & \frac{\partial ξ}{\partial b_{h 3}} \end{matrix}]

$J_{b_{h}}$

will be defined in the code as Jbh and will be computed by the gradient_bias_hidden(Eh) method defined below.

In [6]:

# Define the error function at the hidden layer
def error_hidden(H, Wo, Eo):
    # H * (1-H) * (E . Wo^T)
    return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))
# Define the gradient function for the weight parameters at the hidden layer
def gradient_weight_hidden(X, Eh):
    return X.T.dot(Eh)
# Define the gradient function for the bias parameters at the output layer
def gradient_bias_hidden(Eh): 
    return  np.sum(Eh, axis=0, keepdims=True)

# Define the error function at the hidden layer

def error_hidden(H, Wo, Eo):

# H * (1-H) * (E . Wo^T)

return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))

# Define the gradient function for the weight parameters at the hidden layer

def gradient_weight_hidden(X, Eh):

return X.T.dot(Eh)

# Define the gradient function for the bias parameters at the output layer

def gradient_bias_hidden(Eh):

return np.sum(Eh, axis=0, keepdims=True)

# Define the error function at the hidden layer
def error_hidden(H, Wo, Eo):
    # H * (1-H) * (E . Wo^T)
    return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))
# Define the gradient function for the weight parameters at the hidden layer
def gradient_weight_hidden(X, Eh):
    return X.T.dot(Eh)
# Define the gradient function for the bias parameters at the output layer
def gradient_bias_hidden(Eh): 
    return  np.sum(Eh, axis=0, keepdims=True)

# Define the error function at the hidden layer

def error_hidden(H, Wo, Eo):

# H * (1-H) * (E . Wo^T)

return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))

# Define the gradient function for the weight parameters at the hidden layer

def gradient_weight_hidden(X, Eh):

return X.T.dot(Eh)

# Define the gradient function for the bias parameters at the output layer

def gradient_bias_hidden(Eh):

return np.sum(Eh, axis=0, keepdims=True)

# Define the error function at the hidden layer
def error_hidden(H, Wo, Eo):
    # H * (1-H) * (E . Wo^T)
    return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))
# Define the gradient function for the weight parameters at the hidden layer
def gradient_weight_hidden(X, Eh):
    return X.T.dot(Eh)
# Define the gradient function for the bias parameters at the output layer
def gradient_bias_hidden(Eh): 
    return  np.sum(Eh, axis=0, keepdims=True)

# Define the error function at the hidden layer

def error_hidden(H, Wo, Eo):

# H * (1-H) * (E . Wo^T)

return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))

# Define the gradient function for the weight parameters at the hidden layer

def gradient_weight_hidden(X, Eh):

return X.T.dot(Eh)

# Define the gradient function for the bias parameters at the output layer

def gradient_bias_hidden(Eh):

return np.sum(Eh, axis=0, keepdims=True)

# Define the error function at the hidden layer
def error_hidden(H, Wo, Eo):
    # H * (1-H) * (E . Wo^T)
    return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))
# Define the gradient function for the weight parameters at the hidden layer
def gradient_weight_hidden(X, Eh):
    return X.T.dot(Eh)
# Define the gradient function for the bias parameters at the output layer
def gradient_bias_hidden(Eh): 
    return  np.sum(Eh, axis=0, keepdims=True)

# Define the error function at the hidden layer

def error_hidden(H, Wo, Eo):

# H * (1-H) * (E . Wo^T)

return np.multiply(np.multiply(H,(1 - H)), Eo.dot(Wo.T))

# Define the gradient function for the weight parameters at the hidden layer

def gradient_weight_hidden(X, Eh):

return X.T.dot(Eh)

# Define the gradient function for the bias parameters at the output layer

def gradient_bias_hidden(Eh):

return np.sum(Eh, axis=0, keepdims=True)

Gradient checking

Programming the computation of the backpropagation gradient is prone to bugs. This is why it is recommended to check the gradients of your models. Gradient checking is done by computing the numerical gradient of each parameter, and compare this value to the gradient found by backpropagation.

The numerical gradient $\partial ξ / \partial θ_{i}$

for a parameter $θ_{i}$

can be computed by:

\frac{\partial ξ}{\partial θ_{i}} = \frac{f (X, θ_{1}, \dots, θ_{i} + ϵ, \dots, θ_{m}) - f (X, θ_{1}, \dots, θ_{i} - ϵ, \dots, θ_{m})}{2 ϵ}

Where $f$

if the neural network function that takes the input $X$ and all parameters $θ$ , and $ϵ$ is the small change that is used to peturbate the parameter $θ_{i}$

The numerical gradient for each parameter should be close to the backpropagation gradient for that parameter.

In [7]:

# Initialize weights and biases
init_var = 1
# Initialize hidden layer parameters
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
# Initialize output layer parameters
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var
# Compute the gradients by backpropagation
# Compute the activations of the layers
H = hidden_activations(X, Wh, bh)
Y = output_activations(H, Wo, bo)
# Compute the gradients of the output layer
Eo = error_output(Y, T)
JWo = gradient_weight_out(H, Eo)
Jbo = gradient_bias_out(Eo)
# Compute the gradients of the hidden layer
Eh = error_hidden(H, Wo, Eo)
JWh = gradient_weight_hidden(X, Eh)
Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases

init_var = 1

# Initialize hidden layer parameters

bh = np.random.randn(1, 3) * init_var

Wh = np.random.randn(2, 3) * init_var

# Initialize output layer parameters

bo = np.random.randn(1, 2) * init_var

Wo = np.random.randn(3, 2) * init_var

# Compute the gradients by backpropagation

# Compute the activations of the layers

H = hidden_activations(X, Wh, bh)

Y = output_activations(H, Wo, bo)

# Compute the gradients of the output layer

Eo = error_output(Y, T)

JWo = gradient_weight_out(H, Eo)

Jbo = gradient_bias_out(Eo)

# Compute the gradients of the hidden layer

Eh = error_hidden(H, Wo, Eo)

JWh = gradient_weight_hidden(X, Eh)

Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases
init_var = 1
# Initialize hidden layer parameters
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
# Initialize output layer parameters
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var
# Compute the gradients by backpropagation
# Compute the activations of the layers
H = hidden_activations(X, Wh, bh)
Y = output_activations(H, Wo, bo)
# Compute the gradients of the output layer
Eo = error_output(Y, T)
JWo = gradient_weight_out(H, Eo)
Jbo = gradient_bias_out(Eo)
# Compute the gradients of the hidden layer
Eh = error_hidden(H, Wo, Eo)
JWh = gradient_weight_hidden(X, Eh)
Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases

init_var = 1

# Initialize hidden layer parameters

bh = np.random.randn(1, 3) * init_var

Wh = np.random.randn(2, 3) * init_var

# Initialize output layer parameters

bo = np.random.randn(1, 2) * init_var

Wo = np.random.randn(3, 2) * init_var

# Compute the gradients by backpropagation

# Compute the activations of the layers

H = hidden_activations(X, Wh, bh)

Y = output_activations(H, Wo, bo)

# Compute the gradients of the output layer

Eo = error_output(Y, T)

JWo = gradient_weight_out(H, Eo)

Jbo = gradient_bias_out(Eo)

# Compute the gradients of the hidden layer

Eh = error_hidden(H, Wo, Eo)

JWh = gradient_weight_hidden(X, Eh)

Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases
init_var = 1
# Initialize hidden layer parameters
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
# Initialize output layer parameters
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var
# Compute the gradients by backpropagation
# Compute the activations of the layers
H = hidden_activations(X, Wh, bh)
Y = output_activations(H, Wo, bo)
# Compute the gradients of the output layer
Eo = error_output(Y, T)
JWo = gradient_weight_out(H, Eo)
Jbo = gradient_bias_out(Eo)
# Compute the gradients of the hidden layer
Eh = error_hidden(H, Wo, Eo)
JWh = gradient_weight_hidden(X, Eh)
Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases

init_var = 1

# Initialize hidden layer parameters

bh = np.random.randn(1, 3) * init_var

Wh = np.random.randn(2, 3) * init_var

# Initialize output layer parameters

bo = np.random.randn(1, 2) * init_var

Wo = np.random.randn(3, 2) * init_var

# Compute the gradients by backpropagation

# Compute the activations of the layers

H = hidden_activations(X, Wh, bh)

Y = output_activations(H, Wo, bo)

# Compute the gradients of the output layer

Eo = error_output(Y, T)

JWo = gradient_weight_out(H, Eo)

Jbo = gradient_bias_out(Eo)

# Compute the gradients of the hidden layer

Eh = error_hidden(H, Wo, Eo)

JWh = gradient_weight_hidden(X, Eh)

Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases
init_var = 1
# Initialize hidden layer parameters
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
# Initialize output layer parameters
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var
# Compute the gradients by backpropagation
# Compute the activations of the layers
H = hidden_activations(X, Wh, bh)
Y = output_activations(H, Wo, bo)
# Compute the gradients of the output layer
Eo = error_output(Y, T)
JWo = gradient_weight_out(H, Eo)
Jbo = gradient_bias_out(Eo)
# Compute the gradients of the hidden layer
Eh = error_hidden(H, Wo, Eo)
JWh = gradient_weight_hidden(X, Eh)
Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases

init_var = 1

# Initialize hidden layer parameters

bh = np.random.randn(1, 3) * init_var

Wh = np.random.randn(2, 3) * init_var

# Initialize output layer parameters

bo = np.random.randn(1, 2) * init_var

Wo = np.random.randn(3, 2) * init_var

# Compute the gradients by backpropagation

# Compute the activations of the layers

H = hidden_activations(X, Wh, bh)

Y = output_activations(H, Wo, bo)

# Compute the gradients of the output layer

Eo = error_output(Y, T)

JWo = gradient_weight_out(H, Eo)

Jbo = gradient_bias_out(Eo)

# Compute the gradients of the hidden layer

Eh = error_hidden(H, Wo, Eo)

JWh = gradient_weight_hidden(X, Eh)

Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases
init_var = 1
# Initialize hidden layer parameters
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
# Initialize output layer parameters
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var
# Compute the gradients by backpropagation
# Compute the activations of the layers
H = hidden_activations(X, Wh, bh)
Y = output_activations(H, Wo, bo)
# Compute the gradients of the output layer
Eo = error_output(Y, T)
JWo = gradient_weight_out(H, Eo)
Jbo = gradient_bias_out(Eo)
# Compute the gradients of the hidden layer
Eh = error_hidden(H, Wo, Eo)
JWh = gradient_weight_hidden(X, Eh)
Jbh = gradient_bias_hidden(Eh)

# Initialize weights and biases

init_var = 1

# Initialize hidden layer parameters

bh = np.random.randn(1, 3) * init_var

Wh = np.random.randn(2, 3) * init_var

# Initialize output layer parameters

bo = np.random.randn(1, 2) * init_var

Wo = np.random.randn(3, 2) * init_var

# Compute the gradients by backpropagation

# Compute the activations of the layers

H = hidden_activations(X, Wh, bh)

Y = output_activations(H, Wo, bo)

# Compute the gradients of the output layer

Eo = error_output(Y, T)

JWo = gradient_weight_out(H, Eo)

Jbo = gradient_bias_out(Eo)

# Compute the gradients of the hidden layer

Eh = error_hidden(H, Wo, Eo)

JWh = gradient_weight_hidden(X, Eh)

Jbh = gradient_bias_hidden(Eh)

In [8]:

# Combine all parameter matrices in a list
params = [Wh, bh, Wo, bo]
# Combine all parameter gradients in a list
grad_params = [JWh, Jbh, JWo, Jbo]
# Set the small change to compute the numerical gradient
eps = 0.0001
# Check each parameter matrix
for p_idx in range(len(params)):
    # Check each parameter in each parameter matrix
    for row in range(params[p_idx].shape[0]):
        for col in range(params[p_idx].shape[1]):
            # Copy the parameter matrix and change the current parameter slightly
            p_matrix_min = params[p_idx].copy()
            p_matrix_min[row,col] -= eps
            p_matrix_plus = params[p_idx].copy()
            p_matrix_plus[row,col] += eps
            # Copy the parameter list, and change the updated parameter matrix
            params_min = params[:]
            params_min[p_idx] = p_matrix_min
            params_plus = params[:]
            params_plus[p_idx] =  p_matrix_plus
            # Compute the numerical gradient
            grad_num = (cost(nn(X, *params_plus), T)-cost(nn(X, *params_min), T))/(2*eps)
            # Raise error if the numerical grade is not close to the backprop gradient
            if not np.isclose(grad_num, grad_params[p_idx][row,col]):
                raise ValueError('Numerical gradient of {:.6f} is not close to the backpropagation gradient of {:.6f}!'.format(float(grad_num), float(grad_params[p_idx][row,col])))
print('No gradient errors found')

# Combine all parameter matrices in a list

params = [Wh, bh, Wo, bo]

# Combine all parameter gradients in a list

grad_params = [JWh, Jbh, JWo, Jbo]

# Set the small change to compute the numerical gradient

eps = 0.0001

# Check each parameter matrix

for p_idx in range(len(params)):

# Check each parameter in each parameter matrix

for row in range(params[p_idx].shape[0]):

for col in range(params[p_idx].shape[1]):

# Copy the parameter matrix and change the current parameter slightly

p_matrix_min = params[p_idx].copy()

p_matrix_min[row,col] -= eps

p_matrix_plus = params[p_idx].copy()

p_matrix_plus[row,col] += eps

# Copy the parameter list, and change the updated parameter matrix

params_min = params[:]

params_min[p_idx] = p_matrix_min

params_plus = params[:]

params_plus[p_idx] = p_matrix_plus

# Compute the numerical gradient

grad_num = (cost(nn(X, *params_plus), T)-cost(nn(X, *params_min), T))/(2*eps)

# Raise error if the numerical grade is not close to the backprop gradient

if not np.isclose(grad_num, grad_params[p_idx][row,col]):

raise ValueError('Numerical gradient of {:.6f} is not close to the backpropagation gradient of {:.6f}!'.format(float(grad_num), float(grad_params[p_idx][row,col])))

print('No gradient errors found')

# Combine all parameter matrices in a list
params = [Wh, bh, Wo, bo]
# Combine all parameter gradients in a list
grad_params = [JWh, Jbh, JWo, Jbo]
# Set the small change to compute the numerical gradient
eps = 0.0001
# Check each parameter matrix
for p_idx in range(len(params)):
    # Check each parameter in each parameter matrix
    for row in range(params[p_idx].shape[0]):
        for col in range(params[p_idx].shape[1]):
            # Copy the parameter matrix and change the current parameter slightly
            p_matrix_min = params[p_idx].copy()
            p_matrix_min[row,col] -= eps
            p_matrix_plus = params[p_idx].copy()
            p_matrix_plus[row,col] += eps
            # Copy the parameter list, and change the updated parameter matrix
            params_min = params[:]
            params_min[p_idx] = p_matrix_min
            params_plus = params[:]
            params_plus[p_idx] =  p_matrix_plus
            # Compute the numerical gradient
            grad_num = (cost(nn(X, *params_plus), T)-cost(nn(X, *params_min), T))/(2*eps)
            # Raise error if the numerical grade is not close to the backprop gradient
            if not np.isclose(grad_num, grad_params[p_idx][row,col]):
                raise ValueError('Numerical gradient of {:.6f} is not close to the backpropagation gradient of {:.6f}!'.format(float(grad_num), float(grad_params[p_idx][row,col])))
print('No gradient errors found')

# Combine all parameter matrices in a list

params = [Wh, bh, Wo, bo]

# Combine all parameter gradients in a list

grad_params = [JWh, Jbh, JWo, Jbo]

# Set the small change to compute the numerical gradient

eps = 0.0001

# Check each parameter matrix

for p_idx in range(len(params)):

# Check each parameter in each parameter matrix

for row in range(params[p_idx].shape[0]):

for col in range(params[p_idx].shape[1]):

# Copy the parameter matrix and change the current parameter slightly

p_matrix_min = params[p_idx].copy()

p_matrix_min[row,col] -= eps

p_matrix_plus = params[p_idx].copy()

p_matrix_plus[row,col] += eps

# Copy the parameter list, and change the updated parameter matrix

params_min = params[:]

params_min[p_idx] = p_matrix_min

params_plus = params[:]

params_plus[p_idx] = p_matrix_plus

# Compute the numerical gradient

grad_num = (cost(nn(X, *params_plus), T)-cost(nn(X, *params_min), T))/(2*eps)

# Raise error if the numerical grade is not close to the backprop gradient

if not np.isclose(grad_num, grad_params[p_idx][row,col]):

raise ValueError('Numerical gradient of {:.6f} is not close to the backpropagation gradient of {:.6f}!'.format(float(grad_num), float(grad_params[p_idx][row,col])))

print('No gradient errors found')

No gradient errors found

1 2	No gradient errors found

Backpropagation updates with momentum

In the previous examples we used simple gradient decent to optimize the parameters with respect to the cost function (which was convex in the first and second example). Multilayer neural networks with non-linear activation functions and many parameters are highly unlikely to have convex cost functions. Simple gradient decent is not the best method to find a global minimum of a non-convex cost function since it is a local optimization method that tends to converge to a local minimum .

To solve this issue this example uses a variant of gradient decent called momentum . The method of momentum can be visualized as a ball rolling down the surface of the cost function, this ball will increase in velocity while rolling downhill, and will decrease in velocity when rolling uphill. This momentum update is formulated as:

\begin{aligned} V (i + 1) & = λ V (i) - μ \frac{\partial ξ}{\partial θ (i)} \\ θ (i + 1) & = θ (i) + V (i + 1) \end{aligned}

Where $V (i)$

is the velocity of the parameters at iteration $i$ , $θ (i)$ is the location of the parameters at iteration $i$ , $\partial ξ / \partial θ (i)$ the gradient of the parameters with respect to the cost function at iteration $i$ , $λ$ how much the velocity decreases due to ‘resistance’ and $μ$

the learning rate (how much the gradient affects the velocity). This formula can be visualised as in the following illustration:

The velocities $V_{W_{h}}, V_{b_{h}}, V_{W_{o}}, V_{b_{o}},$

corresponding to the parameters $W_{h}, b_{h}, W_{o}, b_{o}$

are represented in the code by VWh , Vbh , VWo , Vbo . And kept in a list Vs . They are updated by the update_velocity(X, T, ls_of_params, Vs, momentum_term, learning_rate) method. After updating the velocities the parameters are updated via the update_params(ls_of_params, Vs) method.

In [9]:

# Define the update function to update the network parameters over 1 iteration
def backprop_gradients(X, T, Wh, bh, Wo, bo):
    # Compute the output of the network
    # Compute the activations of the layers
    H = hidden_activations(X, Wh, bh)
    Y = output_activations(H, Wo, bo)
    # Compute the gradients of the output layer
    Eo = error_output(Y, T)
    JWo = gradient_weight_out(H, Eo)
    Jbo = gradient_bias_out(Eo)
    # Compute the gradients of the hidden layer
    Eh = error_hidden(H, Wo, Eo)
    JWh = gradient_weight_hidden(X, Eh)
    Jbh = gradient_bias_hidden(Eh)
    return [JWh, Jbh, JWo, Jbo]
def update_velocity(X, T, ls_of_params, Vs, momentum_term, learning_rate):
    # ls_of_params = [Wh, bh, Wo, bo]
    # Js = [JWh, Jbh, JWo, Jbo]
    Js = backprop_gradients(X, T, *ls_of_params)
    return [momentum_term * V - learning_rate * J for V,J in zip(Vs, Js)]
def update_params(ls_of_params, Vs):
    # ls_of_params = [Wh, bh, Wo, bo]
    # Vs = [VWh, Vbh, VWo, Vbo]
    return [P + V for P,V in zip(ls_of_params, Vs)]

# Define the update function to update the network parameters over 1 iteration

def backprop_gradients(X, T, Wh, bh, Wo, bo):

# Compute the output of the network

# Compute the activations of the layers

H = hidden_activations(X, Wh, bh)

Y = output_activations(H, Wo, bo)

# Compute the gradients of the output layer

Eo = error_output(Y, T)

JWo = gradient_weight_out(H, Eo)

Jbo = gradient_bias_out(Eo)

# Compute the gradients of the hidden layer

Eh = error_hidden(H, Wo, Eo)

JWh = gradient_weight_hidden(X, Eh)

Jbh = gradient_bias_hidden(Eh)

return [JWh, Jbh, JWo, Jbo]

def update_velocity(X, T, ls_of_params, Vs, momentum_term, learning_rate):

# ls_of_params = [Wh, bh, Wo, bo]

# Js = [JWh, Jbh, JWo, Jbo]

Js = backprop_gradients(X, T, *ls_of_params)

return [momentum_term * V - learning_rate * J for V,J in zip(Vs, Js)]

def update_params(ls_of_params, Vs):

# ls_of_params = [Wh, bh, Wo, bo]

# Vs = [VWh, Vbh, VWo, Vbo]

return [P + V for P,V in zip(ls_of_params, Vs)]

Training the neural network model with this momentum method over $300$

iterations on the dataset $X$ and $T$

defined in the beginning results in the training convergence plotted in the figure below. This smooth convergence would almost never result from simple gradient descent. You can try it out yourself and implement gradient descent optimization for this model. You will notice that if you run your gradient descent a couple of times from different starting points that it will almost always converge to a sub-optimal solution. Note that momentum is less sensitive to the learning rate than gradient descent, but the learning rate hyperparameter still needs to be tuned separately. In this tutorial, this learning rate was tuned by experimentation.

In [10]:

# Run backpropagation
# Initialize weights and biases
init_var = 0.1
# Initialize hidden layer parameters
bh = np.random.randn(1, 3) * init_var
Wh = np.random.randn(2, 3) * init_var
# Initialize output layer parameters
bo = np.random.randn(1, 2) * init_var
Wo = np.random.randn(3, 2) * init_var
# Parameters are already initilized randomly with the gradient checking
# Set the learning rate
learning_rate = 0.02
momentum_term = 0.9
# define the velocities Vs = [VWh, Vbh, VWo, Vbo]
Vs = [np.zeros_like(M) for M in [Wh, bh, Wo, bo]]
# Start the gradient descent updates and plot the iterations
nb_of_iterations = 300  # number of gradient descent updates
lr_update = learning_rate / nb_of_iterations # learning rate update rule
ls_costs = [cost(nn(X, Wh, bh, Wo, bo), T)]  # list of cost over the iterations
for i in range(nb_of_iterations):
    # Update the velocities and the parameters
    Vs = update_velocity(X, T, [Wh, bh, Wo, bo], Vs, momentum_term, learning_rate)
    Wh, bh, Wo, bo = update_params([Wh, bh, Wo, bo], Vs)
    ls_costs.append(cost(nn(X, Wh, bh, Wo, bo), T))

# Run backpropagation

# Initialize weights and biases

init_var = 0.1

# Initialize hidden layer parameters

bh = np.random.randn(1, 3) * init_var

Wh = np.random.randn(2, 3) * init_var

# Initialize output layer parameters

bo = np.random.randn(1, 2) * init_var

Wo = np.random.randn(3, 2) * init_var

# Parameters are already initilized randomly with the gradient checking

# Set the learning rate

learning_rate = 0.02

momentum_term = 0.9

# define the velocities Vs = [VWh, Vbh, VWo, Vbo]

Vs = [np.zeros_like(M) for M in [Wh, bh, Wo, bo]]

# Start the gradient descent updates and plot the iterations

nb_of_iterations = 300 # number of gradient descent updates

lr_update = learning_rate / nb_of_iterations # learning rate update rule

ls_costs = [cost(nn(X, Wh, bh, Wo, bo), T)] # list of cost over the iterations

for i in range(nb_of_iterations):

# Update the velocities and the parameters

Vs = update_velocity(X, T, [Wh, bh, Wo, bo], Vs, momentum_term, learning_rate)

Wh, bh, Wo, bo = update_params([Wh, bh, Wo, bo], Vs)

ls_costs.append(cost(nn(X, Wh, bh, Wo, bo), T))

In [11]:

# Plot the cost over the iterations
plt.plot(ls_costs, 'b-')
plt.xlabel('iteration')
plt.ylabel('$\\xi$', fontsize=15)
plt.title('Decrease of cost over backprop iteration')
plt.grid()
plt.show()

# Plot the cost over the iterations

plt.plot(ls_costs, 'b-')

plt.xlabel('iteration')

plt.ylabel('$\\xi$', fontsize=15)

plt.title('Decrease of cost over backprop iteration')

plt.grid()

plt.show()