In the previous section, we covered a lot of the theory around neural networks, which can be a little bit overwhelming if you are new to this topic. Before we continue with the discussion of the algorithm for learning the weights of the MLP model, backpropagation, let's take a short break from the theory and see a neural network in action.
Note
The neural network theory can be quite complex, thus I want to recommend two additional resources, which cover some of the concepts that we discuss in this chapter in more detail:
- Chapter 6, Deep Feedforward Networks, Deep Learning, I. Goodfellow, Y. Bengio, and A. Courville, MIT Press, 2016. (Manuscripts freely accessible at http://www.deeplearningbook.org.)
- Pattern Recognition and Machine Learning, C. M. Bishop and others, Volume 1. Springer New York, 2006.
In this section, we will implement and train our first multilayer neural network to classify handwritten digits from the popular Mixed National Institute of Standards and Technology (MNIST) dataset that has been constructed by Yann LeCun and others, and serves as a popular benchmark dataset for machine learning algorithms (Gradient-Based Learning Applied to Document Recognition, Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, Proceedings of the IEEE, 86(11): 2278-2324, November 1998).
The MNIST dataset is publicly available at http://yann.lecun.com/exdb/mnist/ and consists of the following four parts:
- Training set images:
train-images-idx3-ubyte.gz(9.9 MB, 47 MB unzipped, and 60,000 samples) - Training set labels:
train-labels-idx1-ubyte.gz(29 KB, 60 KB unzipped, and 60,000 labels) - Test set images:
t10k-images-idx3-ubyte.gz(1.6 MB, 7.8 MB, unzipped and 10,000 samples) - Test set labels:
t10k-labels-idx1-ubyte.gz(5 KB, 10 KB unzipped, and 10,000 labels)
The MNIST dataset was constructed from two datasets of the US National Institute of Standards and Technology (NIST). The training set consists of handwritten digits from 250 different people, 50 percent high school students, and 50 percent employees from the Census Bureau. Note that the test set contains handwritten digits from different people following the same split. After downloading the files, I recommend that you unzip the files using the Unix/Linux gzip tool from the Terminal for efficiency, using the following command in your local MNIST download directory:
gzip *ubyte.gz -d
Alternatively, you could use your favorite unzipping tool if you are working with a machine running on Microsoft Windows. The images are stored in byte format, and we will read them into NumPy arrays that we will use to train and test our MLP implementation. In order to do that, we will define the following helper function:
import os
import struct
import numpy as np
def load_mnist(path, kind='train'):
"""Load MNIST data from `path`"""
labels_path = os.path.join(path,
'%s-labels-idx1-ubyte' % kind)
images_path = os.path.join(path,
'%s-images-idx3-ubyte' % kind)
with open(labels_path, 'rb') as lbpath:
magic, n = struct.unpack('>II',
lbpath.read(8))
labels = np.fromfile(lbpath,
dtype=np.uint8)
with open(images_path, 'rb') as imgpath:
magic, num, rows, cols = struct.unpack(">IIII",
imgpath.read(16))
images = np.fromfile(imgpath,
dtype=np.uint8).reshape(
len(labels), 784)
images = ((images / 255.) - .5) * 2
return images, labelsThe load_mnist function returns two arrays, the first being an n x m dimensional NumPy array (images), where n is the number of samples and m is the number of features (here, pixels). The training dataset consists of 60,000 training digits and the test set contains 10,000 samples, respectively. The images in the MNIST dataset consist of 28 x 28 pixels, and each pixel is represented by a gray scale intensity value. Here, we unroll the 28 x 28 pixels into one-dimensional row vectors, which represent the rows in our images array (784 per row or image). The second array (labels) returned by the load_mnist function contains the corresponding target variable, the class labels (integers 0-9) of the handwritten digits.
The way we read in the image might seem a little bit strange at first:
>>> magic, n = struct.unpack('>II', lbpath.read(8))
>>> labels = np.fromfile(lbpath, dtype=np.int8)To understand how those two lines of code work, let's take a look at the dataset description from the MNIST website:
[offset] [type] [value] [description] 0000 32 bit integer 0x00000801(2049) magic number (MSB first) 0004 32 bit integer 60000 number of items 0008 unsigned byte ?? label 0009 unsigned byte ?? label ........ xxxx unsigned byte ?? label
Using the two preceding lines of code, we first read in the magic number, which is a description of the file protocol as well as the number of items (n) from the file buffer before we read the following bytes into a NumPy array using the fromfile method. The fmt parameter value '>II' that we passed as an argument to struct.unpack can be composed into the two following parts:
>: This is big-endian—it defines the order in which a sequence of bytes is stored; if you are unfamiliar with the terms big-endian and little-endian, you can find an excellent article about Endianness on Wikipedia: https://en.wikipedia.org/wiki/EndiannessI: This is an unsigned integer
Finally, we also normalized the pixels values in MNIST to the range -1 to 1 (originally 0 to 255) via the following code line:
images = ((images / 255.) - .5) * 2
The reason behind this is that gradient-based optimization is much more stable under these conditions as discussed in Chapter 2, Training Simple Machine Learning Algorithms for Classification. Note that we scaled the images on a pixel-by-pixel basis, which is different from the feature scaling approach that we took in previous chapters. Previously, we derived scaling parameters from the training set and used these to scale each column in the training set and test set. However, when working with image pixels, centering them at zero and rescaling them to a [-1, 1] range is also common and usually works well in practice.
Note
Another recently developed trick to improve convergence in gradient-based optimization through input scaling is batch normalization, which is an advanced topic that we will not cover in this book. However, if you are interested in deep learning applications and research, I highly recommend that you read more about batch normalization in the excellent research article Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift by Sergey Ioffe and Christian Szegedy (2015, https://arxiv.org/abs/1502.03167).
By executing the following code, we will now load the 60,000 training instances as well as the 10,000 test samples from the local directory where we unzipped the MNIST dataset (in the following code snippet, it is assumed that the downloaded MNIST files were unzipped to the same directory in which this code was executed):
>>> X_train, y_train = load_mnist('', kind='train')
>>> print('Rows: %d, columns: %d'
... % (X_train.shape[0], X_train.shape[1]))
Rows: 60000, columns: 784
>>> X_test, y_test = load_mnist('', kind='t10k')
>>> print('Rows: %d, columns: %d'
... % (X_test.shape[0], X_test.shape[1]))
Rows: 10000, columns: 784To get an idea of how those images in MNIST look, let's visualize examples of the digits 0-9 after reshaping the 784-pixel vectors from our feature matrix into the original 28 × 28 image that we can plot via Matplotlib's imshow function:
>>> import matplotlib.pyplot as plt >>> fig, ax = plt.subplots(nrows=2, ncols=5, ... sharex=True, sharey=True) >>> ax = ax.flatten() >>> for i in range(10): ... img = X_train[y_train == i][0].reshape(28, 28) ... ax[i].imshow(img, cmap='Greys') >>> ax[0].set_xticks([]) >>> ax[0].set_yticks([]) >>> plt.tight_layout() >>> plt.show()
We should now see a plot of the 2 x 5 subfigures showing a representative image of each unique digit:
In addition, let's also plot multiple examples of the same digit to see how different the handwriting really is:
>>> fig, ax = plt.subplots(nrows=5, ... ncols=5, ... sharex=True, ... sharey=True) >>> ax = ax.flatten() >>> for i in range(25): ... img = X_train[y_train == 7][i].reshape(28, 28) ... ax[i].imshow(img, cmap='Greys') >>> ax[0].set_xticks([]) >>> ax[0].set_yticks([]) >>> plt.tight_layout() >>> plt.show()
After executing the code, we should now see the first 25 variants of the digit 7:
After we've gone through all the previous steps, it is a good idea to save the scaled images in a format that we can load more quickly into a new Python session to avoid the overhead of reading in and processing the data again. When we are working with NumPy arrays, an efficient yet most convenient method to save multidimensional arrays to disk is NumPy's savez function (the official documentation can be found here: https://docs.scipy.org/doc/numpy/reference/generated/numpy.savez.html).
In short, the savez function is analogous to Python's pickle module that we used in Chapter 9, Embedding a Machine Learning Model into a Web Application, but optimized for storing NumPy arrays. The savez function produces zipped archives of our data, producing .npz files that contain files in the .npy format; if you want to learn more about this format, you can find a nice explanation, including a discussion about advantages and disadvantages, in the NumPy documentation: https://docs.scipy.org/doc/numpy/neps/npy-format.html. Further, instead of using savez, we will use savez_compressed, which uses the same syntax as savez, but further compresses the output file down to substantially smaller file sizes (approximately 22 MB versus approximately 400 MB in this case). The following code snippet will save both the training and test datasets to the archive file 'mnist_scaled.npz':
>>> import numpy as np
>>> np.savez_compressed('mnist_scaled.npz',
... X_train=X_train,
... y_train=y_train,
... X_test=X_test,
... y_test=y_test)After we created the .npz files, we can load the preprocessed MNIST image arrays using NumPy's load function as follows:
>>> mnist = np.load('mnist_scaled.npz')The mnist variable now references to an object that can access the four data arrays as we provided them keyword arguments to the savez_compressed function, which are listed under the files attribute list of the mnist object:
>>> mnist.files ['X_train', 'y_train', 'X_test', 'y_test']
For instance, to load the training data into our current Python session, we will access the 'X_train' array as follows (similar to a Python dictionary):
>>> X_train = mnist['X_train']
Using a list comprehension, we can retrieve all four data arrays as follows:
>>> X_train, y_train, X_test, y_test = [mnist[f] for ... f in mnist.files]
Note that while the preceding np.savez_compressed and np.load examples are not essential for executing the code in this chapter, it serves as a demonstration of how to save and load NumPy arrays conveniently and efficiently.
In this subsection, we will now implement the code of an MLP with one input, one hidden, and one output layers to classify the images in the MNIST dataset. I have tried to keep the code as simple as possible. However, it may seem a little bit complicated at first, and I encourage you to download the sample code for this chapter from the Packt Publishing website or from GitHub (https://github.com/rasbt/python-machine-learning-book-2nd-edition) so that you can view this MLP implementation annotated with comments and syntax highlighting for better readability.
If you are not running the code from the accompanying Jupyter Notebook file or don't have access to the internet, I recommend that you copy the NeuralNetMLP code from this chapter into a Python script file in your current working directory, for example, neuralnet.py, which you can then import into your current Python session via the following command:
from neuralnet import NeuralNetMLP
The code will contain parts that we have not talked about yet, such as the backpropagation algorithm, but most of the code should look familiar to you based on the Adaline implementation in Chapter 2, Training Simple Machine Learning Algorithms for Classification, and the discussion of forward propagation in earlier sections.
Do not worry if not all of the code makes immediate sense to you; we will follow up on certain parts later in this chapter. However, going over the code at this stage can make it easier to follow the theory later.
The following is the implementation of a multilayer perceptron:
import numpy as np
import sys
class NeuralNetMLP(object):
""" Feedforward neural network / Multi-layer perceptron classifier.
Parameters
------------
n_hidden : int (default: 30)
Number of hidden units.
l2 : float (default: 0.)
Lambda value for L2-regularization.
No regularization if l2=0. (default)
epochs : int (default: 100)
Number of passes over the training set.
eta : float (default: 0.001)
Learning rate.
shuffle : bool (default: True)
Shuffles training data every epoch
if True to prevent circles.
minibatch_size : int (default: 1)
Number of training samples per minibatch.
seed : int (default: None)
Random seed for initializing weights and shuffling.
Attributes
-----------
eval_ : dict
Dictionary collecting the cost, training accuracy,
and validation accuracy for each epoch during training.
"""
def __init__(self, n_hidden=30,
l2=0., epochs=100, eta=0.001,
shuffle=True, minibatch_size=1, seed=None):
self.random = np.random.RandomState(seed)
self.n_hidden = n_hidden
self.l2 = l2
self.epochs = epochs
self.eta = eta
self.shuffle = shuffle
self.minibatch_size = minibatch_size
def _onehot(self, y, n_classes):
"""Encode labels into one-hot representation
Parameters
------------
y : array, shape = [n_samples]
Target values.
Returns
-----------
onehot : array, shape = (n_samples, n_labels)
"""
onehot = np.zeros((n_classes, y.shape[0]))
for idx, val in enumerate(y.astype(int)):
onehot[val, idx] = 1.
return onehot.T
def _sigmoid(self, z):
"""Compute logistic function (sigmoid)"""
return 1. / (1. + np.exp(-np.clip(z, -250, 250)))
def _forward(self, X):
"""Compute forward propagation step"""
# step 1: net input of hidden layer
# [n_samples, n_features] dot [n_features, n_hidden]
# -> [n_samples, n_hidden]
z_h = np.dot(X, self.w_h) + self.b_h
# step 2: activation of hidden layer
a_h = self._sigmoid(z_h)
# step 3: net input of output layer
# [n_samples, n_hidden] dot [n_hidden, n_classlabels]
# -> [n_samples, n_classlabels]
z_out = np.dot(a_h, self.w_out) + self.b_out
# step 4: activation output layer
a_out = self._sigmoid(z_out)
return z_h, a_h, z_out, a_out
def _compute_cost(self, y_enc, output):
"""Compute cost function.
Parameters
----------
y_enc : array, shape = (n_samples, n_labels)
one-hot encoded class labels.
output : array, shape = [n_samples, n_output_units]
Activation of the output layer (forward propagation)
Returns
---------
cost : float
Regularized cost
"""
L2_term = (self.l2 *
(np.sum(self.w_h ** 2.) +
np.sum(self.w_out ** 2.)))
term1 = -y_enc * (np.log(output))
term2 = (1. - y_enc) * np.log(1. - output)
cost = np.sum(term1 - term2) + L2_term
return cost
def predict(self, X):
"""Predict class labels
Parameters
-----------
X : array, shape = [n_samples, n_features]
Input layer with original features.
Returns:
----------
y_pred : array, shape = [n_samples]
Predicted class labels.
"""
z_h, a_h, z_out, a_out = self._forward(X)
y_pred = np.argmax(z_out, axis=1)
return y_pred
def fit(self, X_train, y_train, X_valid, y_valid):
""" Learn weights from training data.
Parameters
-----------
X_train : array, shape = [n_samples, n_features]
Input layer with original features.
y_train : array, shape = [n_samples]
Target class labels.
X_valid : array, shape = [n_samples, n_features]
Sample features for validation during training
y_valid : array, shape = [n_samples]
Sample labels for validation during training
Returns:
----------
self
"""
n_output = np.unique(y_train).shape[0] # no. of class
#labels
n_features = X_train.shape[1]
########################
# Weight initialization
########################
# weights for input -> hidden
self.b_h = np.zeros(self.n_hidden)
self.w_h = self.random.normal(loc=0.0, scale=0.1,
size=(n_features,
self.n_hidden))
# weights for hidden -> output
self.b_out = np.zeros(n_output)
self.w_out = self.random.normal(loc=0.0, scale=0.1,
size=(self.n_hidden,
n_output))
epoch_strlen = len(str(self.epochs)) # for progr. format.
self.eval_ = {'cost': [], 'train_acc': [], 'valid_acc':
[]}
y_train_enc = self._onehot(y_train, n_output)
# iterate over training epochs
for i in range(self.epochs):
# iterate over minibatches
indices = np.arange(X_train.shape[0])
if self.shuffle:
self.random.shuffle(indices)
for start_idx in range(0, indices.shape[0] -
self.minibatch_size +
1, self.minibatch_size):
batch_idx = indices[start_idx:start_idx +
self.minibatch_size]
# forward propagation
z_h, a_h, z_out, a_out =
self._forward(X_train[batch_idx])
##################
# Backpropagation
##################
# [n_samples, n_classlabels]
sigma_out = a_out - y_train_enc[batch_idx]
# [n_samples, n_hidden]
sigmoid_derivative_h = a_h * (1. - a_h)
# [n_samples, n_classlabels] dot [n_classlabels,
# n_hidden]
# -> [n_samples, n_hidden]
sigma_h = (np.dot(sigma_out, self.w_out.T) *
sigmoid_derivative_h)
# [n_features, n_samples] dot [n_samples,
# n_hidden]
# -> [n_features, n_hidden]
grad_w_h = np.dot(X_train[batch_idx].T, sigma_h)
grad_b_h = np.sum(sigma_h, axis=0)
# [n_hidden, n_samples] dot [n_samples,
# n_classlabels]
# -> [n_hidden, n_classlabels]
grad_w_out = np.dot(a_h.T, sigma_out)
grad_b_out = np.sum(sigma_out, axis=0)
# Regularization and weight updates
delta_w_h = (grad_w_h + self.l2*self.w_h)
delta_b_h = grad_b_h # bias is not regularized
self.w_h -= self.eta * delta_w_h
self.b_h -= self.eta * delta_b_h
delta_w_out = (grad_w_out + self.l2*self.w_out)
delta_b_out = grad_b_out # bias is not regularized
self.w_out -= self.eta * delta_w_out
self.b_out -= self.eta * delta_b_out
#############
# Evaluation
#############
# Evaluation after each epoch during training
z_h, a_h, z_out, a_out = self._forward(X_train)
cost = self._compute_cost(y_enc=y_train_enc,
output=a_out)
y_train_pred = self.predict(X_train)
y_valid_pred = self.predict(X_valid)
train_acc = ((np.sum(y_train ==
y_train_pred)).astype(np.float) /
X_train.shape[0])
valid_acc = ((np.sum(y_valid ==
y_valid_pred)).astype(np.float) /
X_valid.shape[0])
sys.stderr.write('
%0*d/%d | Cost: %.2f '
'| Train/Valid Acc.: %.2f%%/%.2f%% '
%
(epoch_strlen, i+1, self.epochs,
cost,
train_acc*100, valid_acc*100))
sys.stderr.flush()
self.eval_['cost'].append(cost)
self.eval_['train_acc'].append(train_acc)
self.eval_['valid_acc'].append(valid_acc)
return selfOnce you're done with executing this code, let's now initialize a new 784-100-10 MLP—a neural network with 784 input units (n_features), 100 hidden units (n_hidden), and 10 output units (n_output):
>>>nn = NeuralNetMLP(n_hidden=100, ... l2=0.01, ... epochs=200, ... eta=0.0005, ... minibatch_size=100, ... shuffle=True, ... seed=1)
If you read through the NeuralNetMLP code, you've probably already guessed what these parameters are for. Here, you find a short summary of these:
l2: This is the
parameter for L2 regularization to decrease the degree of overfitting.
epochs: This is the number of passes over the training set.eta: This is the learning rate
.
shuffle: This is for shuffling the training set prior to every epoch to prevent that the algorithm gets stuck in circles.seed: This is a random seed for shuffling and weight initialization.minibatch_size: This is the number of training samples in each mini-batch when splitting of the training data in each epoch for stochastic gradient descent. The gradient is computed for each mini-batch separately instead of the entire training data for faster learning.
Next, we train the MLP using 55,000 samples from the already shuffled MNIST training dataset and use the remaining 5,000 samples for validation during training. Note that training the neural network may take up to 5 minutes on standard desktop computer hardware.
As you may have noticed from the preceding code implementation, we implemented the fit method so that it takes four input arguments: training images, training labels, validation images, and validation labels. In neural network training, it is really useful to already compare training and validation accuracy during training, which helps us judge whether the network model performs well, given the architecture and hyperparameters.
In general, training (deep) neural networks is relatively expensive compared with the other models we discussed so far. Thus, we want to stop it early in certain circumstances and start over with different hyperparameter settings. Alternatively, if we find that it increasingly tends to overfit the training data (noticeable by an increasing gap between training and validation set performance), we may want to stop the training early as well.
Now, to start the training, we execute the following code:
>>> nn.fit(X_train=X_train[:55000], ... y_train=y_train[:55000], ... X_valid=X_train[55000:], ... y_valid=y_train[55000:]) 200/200 | Cost: 5065.78 | Train/Valid Acc.: 99.28%/97.98%
In our NeuralNetMLP implementation, we also defined an eval_ attribute that collects the cost, training, and validation accuracy for each epoch so that we can visualize the results using Matplotlib:
>>> import matplotlib.pyplot as plt
>>> plt.plot(range(nn.epochs), nn.eval_['cost'])
>>> plt.ylabel('Cost')
>>> plt.xlabel('Epochs')
>>> plt.show()The preceding code plots the cost over the 200 epochs, as shown in the following graph:
As we can see, the cost decreased substantially during the first 100 epochs and seems to slowly converge in the last 100 epochs. However, the small slope between epoch 175 and epoch 200 indicates that the cost would further decrease with a training over additional epochs.
Next, let's take a look at the training and validation accuracy:
>>> plt.plot(range(nn.epochs), nn.eval_['train_acc'],
... label='training')
>>> plt.plot(range(nn.epochs), nn.eval_['valid_acc'],
... label='validation', linestyle='--')
>>> plt.ylabel('Accuracy')
>>> plt.xlabel('Epochs')
>>> plt.legend()
>>> plt.show()The preceding code examples plot those accuracy values over the 200 training epochs, as shown in the following figure:
The plot reveals that the gap between training and validation accuracy increases the more epochs we train the network. At approximately the 50th epoch, the training and validation accuracy values are equal, and then, the network starts overfitting the training data.
Note that this example was chosen deliberately to illustrate the effect of overfitting and demonstrate why it is useful to compare the validation and training accuracy values during training. One way to decrease the effect of overfitting is to increase the regularization strength—for example, by setting l2=0.1. Another useful technique to tackle overfitting in neural networks, dropout, will be covered in Chapter 15, Classifying Images with Deep Convolutional Neural Networks.
Finally, let's evaluate the generalization performance of the model by calculating the prediction accuracy on the test set:
>>> y_test_pred = nn.predict(X_test)
>>> acc = (np.sum(y_test == y_test_pred)
... .astype(np.float) / X_test.shape[0])
>>> print('Training accuracy: %.2f%%' % (acc * 100))
Test accuracy: 97.54%Despite the slight overfitting on the training data, our relatively simple one-hidden layer neural network achieved a relatively good performance on the test dataset, similar to the validation set accuracy (97.98 percent).
To further fine-tune the model, we could change the number of hidden units, values of the regularization parameters, and the learning rate or use various other tricks that have been developed over the years but are beyond the scope of this book. In Chapter 14, Going Deeper – The Mechanics of TensorFlow, you will learn about a different neural network architecture that is known for its good performance on image datasets. Also, the chapter will introduce additional performance-enhancing tricks such as adaptive learning rates, momentum learning, and dropout.
Lastly, let's take a look at some of the images that our MLP struggles with:
>>> miscl_img = X_test[y_test != y_test_pred][:25]
>>> correct_lab = y_test[y_test != y_test_pred][:25]
>>> miscl_lab= y_test_pred[y_test != y_test_pred][:25]
>>> fig, ax = plt.subplots(nrows=5,
... ncols=5,
... sharex=True,
... sharey=True,)
>>> ax = ax.flatten()
>>> for i in range(25):
... img = miscl_img[i].reshape(28, 28)
... ax[i].imshow(img,
... cmap='Greys',
... interpolation='nearest')
... ax[i].set_title('%d) t: %d p: %d'
... % (i+1, correct_lab[i], miscl_lab[i]))
>>> ax[0].set_xticks([])
>>> ax[0].set_yticks([])
>>> plt.tight_layout()
>>> plt.show()We should now see a 5 x 5 subplot matrix where the first number in the subtitles indicates the plot index, the second number represents the true class label (t), and the third number stands for the predicted class label (p):
As we can see in the preceding figure, some of those images are even challenging for us humans to classify correctly. For example, the 6 in subplot 8 really looks like a carelessly drawn 0, and the 8 in subplot 23 could be a 9 due to the narrow lower part combined with the bold line.