Linux, OpenBlas and Numpy matrix multiplications – avoid using all processor cores

Recently, I tested the propagation methods of a small Python3/Numpy class for a multilayer perceptron [MLP]. I unexpectedly ran into a performance problem with OpenBlas.

The problem had to do with the required vectorized matrix operations for forward propagation - in my case through an artificial neural network [ANN] with 4 layers. In a first approach I used 784, 100, 50, 10 neurons in 4 consecutive layers of the MLP. The weight matrices had corresponding dimensions.

The performance problem was caused by extensive multi-threading; it showed a strong dependency on mini-batch sizes and on basic matrix dimensions related to the neuron numbers per layer:

  • For the given relatively small number of neurons per layer and for mini-batch sizes above a critical value (N > 255) OpenBlas suddenly occupied all processor cores with a 100% work load. This had a disastrous impact on performance.
  • For neuron numbers as 784, 300, 140, 10 OpenBlas used all processor cores with a 100% work load right from the beginning, i.e. even for small batch sizes. With a seemingly bad performance over all batch sizes - but decreasing somewhat with large batch sizes.

This problem has been discussed elsewhere with respect to the matrix dimensions relevant for the core multiplication and summation operations - i.e. the neuron numbers per layer. However, the vectorizing aspect of matrix multiplications is interesting, too:

One can imagine that splitting the operations for multiple independent test samples is in principle ideal for multi-threading. So, using as many processor cores as possible (in my case 8) does not look like a wrong decision of OpenBlas at first.

Then I noticed that for mini-batch sizes "N" below a certain number (N < 250) the system only seemed to use up to 3-4 cores; so there remained plenty of CPU capacity left for other tasks. Performance for N < 250 was better by at least a factor of 2 compared to a situation with an only slightly bigger batch size (N ≥ 260). I got the impression that OpenBLAS under certain conditions just decides to use as many threads as possible - which no good outcome.

In the last years I sometimes had to work with optimizing multi-threaded database operations on Linux systems. I often got the impression that you have to be careful and leave some CPU resources free for other tasks and to avoid heavy context switching. In addition bottlenecks appeared due to the concurrent access of may processes to the CPU cache. (RAM limitations were an additional factor; but this should not be the case for my Python program.) Furthermore, one should not forget that Python/Numpy experiments on Jupyter notebooks require additional resources to handle the web page output and page update on the browser. And Linux itself also requires some free resources.

So, I wanted to find out whether reducing the number of threads - or available cores - for Numpy and OpenBlas would be helpful in the sense of an overall gain in performance.

All data shown below were gathered on a desktop system with some background activity due to several open browsers, clementine and pulse-audio as active audio components, an open mail client (kontact), an open LXC container, open Eclipse with Pydev and open ssh connections. Program tests were performed with the help of Jupyter notebooks. Typical background CPU consumption looks like this on Ksysguard:

Most of the consumption is due to audio. Small spikes on one CPU core due to the investigation of incoming mails were possible - but always below 20%.

Basics

One of the core ingredients to get an ANN running are matrix operations. More precisely: multiplications of 2-dim Numpy matrices (weight matrices) with input vectors. The dimensions of the weight matrices reflect the node-numbers of consecutive ANN-layers. The dimension of the input vector depends on the node number of the lower of two neighbor layers.

However, we perform such matrix operations NOT sequentially sample for sample of a collection of training data - we do it vectorized for so called mini-batches consisting of between 50 and 600000 individual samples of training data. Instead of operating with a matrix on just one feature vector of one training sample we use matrix multiplications whereby the second matrix often comprises many vectors of data samples.

I have described such multiplications already in a previous blog article; see Numpy matrix multiplication for layers of simple feed forward ANNs.

In the most simple case of an MLP with e.g.

  • an input layer of 784 nodes (suitable for the MNIST dataset),
  • one hidden layer with 100 nodes,
  • another hidden layer with 50 nodes
  • and an output layer of 10 nodes (fitting again the MNIST dataset)

and "mini"-batches of different sizes (between 20 and 20000). An input vector to the first hidden layer has a dimension of 100, so the weight matrix creating this input vector from the "output" of the MLP's input layer has a shape of 784x100. Multiplication and summation in this case is done over the dimension covering 784 features. When we work with mini-batches we want to do these operations in parallel for as many elements of a mini-batch as possible.

All in all we have to perform 3 matrix operations

(784x100) matrix on (784)-vector, (100x50) matrix on (100)-vector, (50x10) matrix on (50) vector

on our example ANN with 4 layers. However, we collect the data for N mini-batch samples in an array. This leads to Numpy matrix multiplications of the kind

(784x100) matrix on an (784, N)-array, (100x50) matrix on an (100, N)-array, (50x10) matrix on an (50, N)-array.

Thus, we deal with matrix multiplications of two 2-dim matrices. Linear algebra libraries should optimize such operations for different kinds of processors.

The reaction of OpenBlas to an MLP with 4 layers comprising 784, 100, 50, 10 nodes

On my Linux system Python/Numpy use the openblas-library. This is confirmed by the output of command "np.__config__.show()":

openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]

and by

(ml1) myself@mytux:/projekte/GIT/ai/ml1/lib64/python3.6/site-packages/numpy/core> ldd  _multiarray_umath.cpython-36m-x86_64-linux-gnu.so
        linux-vdso.so.1 (0x00007ffe8bddf000)
        libopenblasp-r0-2ecf47d5.3.7.dev.so => /projekte/GIT/ai/ml1/lib/python3.6/site-packages/numpy/core/./../.libs/libopenblasp-r0-2ecf47d5.3.7.dev.so (0x00007fdd9d15f000)
        libm.so.6 => /lib64/libm.so.6 (0x00007fdd9ce27000)
        libpthread.so.0 => /lib64/libpthread.so.0 (0x00007fdd9cc09000)
        libc.so.6 => /lib64/libc.so.6 (0x00007fdd9c84f000)
        /lib64/ld-linux-x86-64.so.2 (0x00007fdd9f4e8000)
        libgfortran-ed201abd.so.3.0.0 => /projekte/GIT/ai/ml1/lib/python3.6/site-packages/numpy/core/./../.libs/libgfortran-ed201abd.so.3.0.0 (0x00007fdd9c555000)

In all tests discussed below I performed a series of calculations for different batch sizes

N = 50, 100, 200, 250, 260, 500, 2000, 10000, 20000

and repeated the full forward propagation 30 times (corresponding to 30 epochs in a full training series - but here without cost calculation and weight adjustment. I just did forward propagation.)

In a first experiment, I did not artificially limit the number of cores to be used. Measured response times in seconds are indicated in the following plot:

Runtime for a free number of cores to use and different batch-sizes N

We see that something dramatic happens between a batch size of 250 and 260. Below you see the plots for CPU core consumption for N=50, N=200, N=250, N=260 and N=2000.

N=50:

N=200:

N=250:

N=260:

N=2000:

The plots indicate that everything goes well up to N=250. Up to this point around 4 cores are used - leaving 4 cores relatively free. After N=260 OpenBlas decides to use all 8 cores with a load of 100% - and performance suffers by more than a factor of 2.

This result support the idea to look for an optimum of the number of cores "C" to use.

The reaction of OpenBlas to an MLP with layers comprising 784, 300, 140, 10 nodes

For a MLP with neuron numbers (784, 300, 140, 10) I got the red curve for response time in the plot below. The second curve shows what performance is possible with just using 4 cores:

Note the significantly higher response times. We also see again that something strange happens at the change of the batch-size from 250 to 260.

The 100% CPU consumption even for a batch-size of only 50 is shown below:

Though different from the first test case also these plots indicate that - somewhat paradoxically - reducing the number of CPU cores available to OpenBlas could have a performance enhancing effect.

Limiting the number of available cores to OpenBlas

A bit of Internet research shows that one can limit the number of cores to use by OpenBlas e.g. via an environment variable for the shell, in which we start a Jupyter notebook. The relevant command to limit the number of cores "C" to 3 is :

export OPENBLAS_NUM_THREADS=3

Below you find plots for the response times required for the batch sizes N listed above and core numbers of

C=1, C=2, C=3, C=4, C=5, C=6, C=7, C=8 :

For C=5 I did 2 different runs; the different results for C=5 show that the system reacts rather sensitively. It changes its behavior for larger core number drastically.

We also find an overall minimum of the response time:
The overall optimum occurs for 400 < N < 500 for C=1, 2, 3, 4 - with the minimum region being broadest for C=3. The absolute minimum is reached on my CPU for C=4.

We understand from the plots above that the number of cores to use become hyper-parameters for the tuning of the performance of ANNs - at least as long as a standard multicore-CPU is used.

CPU-consumption

CPU consumption for N=50 and C=2 looks like:

For comparison see the CPU consumption for N=20000 and C=4:

CPU consumption for N=20000 and C=6:

We see that between C=5 and C=6 CPU resources get heavily consumed; there are almost no reserves left in the Linux system for C ≥ 6.

Dependency on the size of the weight-matrices and the node numbers

For a full view on the situation I also looked at the response time variation with node numbers for a given number of CPU cores.

For C=4 and node number cases

  • 784, 300, 140, 10
  • 784, 200, 100, 10
  • 784, 100, 50, 10
  • 784, 50, 20, 10

I got the following results:

There is some broad variation with the weight-matrix size; the bigger the weight-matrix the longer the calculation time. This is, of course, to be expected. Note that the variation with the batch-size number is relatively smooth - with an optimum around 400.

Now, look at the same plot for C=6:

Note that the response time is significantly bigger in all cases compared to the previous situation with C=4. In cases of a large matrix by around 36% for N=2000. Also the variation with batch-size is more pronounced.

Still, even with 6 cores you do not get factors between 1.4 and 2.0 as compared to the case of C=8 (see above)!

Conclusion

As I do not know what the authors of OpenBlas are doing exactly, I refrain from technically understanding and interpreting the causes of the data shown above.

However, some consequences seem to be clear:

  • It is a bad idea to provide all CPU cores to OpenBlas - which unfortunately is the default.
  • The data above indicate that using only 4 out of 8 core threads is reasonable to get an optimum performance for vectorized matrix multiplications on my CPU.
  • Not leaving at least 2 CPU cores free for other tasks is punished by significant performance losses.
  • When leaving the decision of how many cores to use to OpenBlas a critical batch-size may exist for which the performance suddenly breaks down due to heavy multi-threading.

Whenever you deal with ANN or MLP simulations on a standard CPU (not GPU!) you should absolutely care about how many cores and related threads you want to offer to OpenBlas. As far as I understood from some Internet articles the number of cores to be used can be not only be controlled by Linux (shell) environment variables but also by os-commands in a Python program. You should perform tests to find optimum values for your CPU.

Links

stackoverflow: numpy-suddenly-uses-all-cpus

stackoverflow: run-openblas-on-multicore

stackoverflow: multiprocessing-pool-makes-numpy-matrix-multiplication-slower

scicomp: why-isnt-my-matrix-vector-multiplication-scaling/1729

Setting the number of threads via Python
stackoverflow:
set-max-number-of-threads-at-runtime-on-numpy-openblas

codereview.stackexchange: better-way-to-set-number-of-threads-used-by-numpy

 

A simple program for an ANN to cover the Mnist dataset – I

For beginners both in Python and Machine Learning [ML] the threshold to do some real programming and create your own ANN seems to be relatively high. Well, some readers might say: Why program an ANN by yourself at a basic Python level at all when Keras and TensorFlow [TF] are available? Answer: For learning! And eventually to be able to do some things TF has not been made for. And as readers of this blog will see soon, I have some ideas along this line ...

So, I thought, just let us set up a small Python3 and Numpy based program to create a simple ANN and train it for the MNIST dataset.

I take a shortcut and assume that readers are already acquainted with the following topics:

  1. what simple neural networks with a hidden layer look like,
  2. what a cost function and the gradient descent method is,
  3. what logistic regression is and what the cost function J for it looks like,
  4. why the back propagation of deviations from known values gives you the required partial derivatives for typical cost functions,
  5. what a mini-batch approach to optimization is.

These basics are well documented in literature; e.g. in the books of Geron and Rashka (see the references at the end of this article). However, I shall briefly comment on these topics whilst building the code.

We need a relatively well defined first objective for the usage of our ANN. We shall concentrate on classification tasks. As a first example we shall use the conventional MNIST data set. The MNIST data set consists of 28x28 px images of handwritten numbers. It is a standard data set used in many elementary courses on ML. The challenge for the ANN is that it should be able to recognize hand-written digits from a digitized gray-color image after some training.

Note that this task does NOT require the use of an ANN. "Stochastic Gradient Descent"-approaches to determine (linear) separation surfaces in combination with a One-versus-All strategy for multi-category-classification may be sufficient. See chapters 3 to 5 in the book of Geron for more information.

Regarding the build up of the ANN program, I basically follow an approach described by Raschka in his book. However, at multiple points I take the freedom to organize the code differently and comment in my own way ... I am only a beginner in Python; I hope my insights are helpful for others in the same situation.

In any case you should make yourself familiar with numpy arrays and their "shapes". I assume that you understand the multidimensional structure of Numpy arrays ....

Wording

To avoid confusion, I use the following wording:

Category: Each input data element is associated with a category to which it belongs. The classification algorithm (here: the ANN) may achieve a distinction between the association of input data with multiple categories. It should - after some training - detect (non-linear) separation interfaces for categories in a multidimensional feature space. In the case of MNIST we speak about ten categories corresponding to 10 digits, including zero.

Label: A category may be described by a label. Training data may provide a so called "target label array" _y_train for all input data. We must be prepared to transform target labels for input data into a usable form for an ANN, i.e. into a vectorized form which selects a specific category out of many. This process is called "label encoding".

Input data set: A complete set of input data. Such a set consists of individual "elements" or "records". The MNIST input set of training data has 60000 records - which we provide via an array _X_train.

Output data set: A complete set of output data after propagation through the ANN for the input data set. The number of elements or records of the output data sets is equal to the number of records in the input data set. The output set will be represented by an array "_ay_a_Out".

A data record of the input data set: One distinct element of the input data set (and its array). Note that such an element itself may be a multidimensional array covering all features in a distinct form. Such arrays represents a so called "tensor".

A data record of the output data set: One distinct element of the output data set. Note that such an element itself may be an array covering all possible categories in a distinct form.

A simple neural network - layers, nodes, weights

A neural network is composed of a series of sequential layers with nodes. All nodes of a specific layer can be connected with all nodes of neighboring layers. The simplified sketch below displays an ANN with just three layers - an input layer, a "hidden" middle layer and an output layer. Note that there can be (many) more hidden layers than just one.

Input layer and its number of nodes
To feed input data into the ANN we need an input layer with input nodes. How many? Well, this depends on the number of features your data set represents. In the MNIST case an 28x28 px image with simple gray values (integer number between 0 and 256) represents 28x28 = 768 different "features" - i.e. numbers for color values. We need as many nodes to represent the full image information by the input layer.

For other input data the number and structure of features may be different; the number of input nodes must then be adjusted accordingly. The number of nodes should, therefore, be a parameter or be derived from information on the type of input data. How complicated and structured features are mapped to a one dimensional input vector is a question one should think about carefully. (Most people today deal with a time set of data as just a special form of a feature - I regard this as questionable in some cases, but this is beyond this article series ...)

Output layer and its number of nodes
We shall use our ANN for classification tasks in the beginning. We, therefore, assume that the output of the ANN should distinguish between "strong>NC" different categories an input data set can belong to. In case of the MNIST dataset we can distinguish between 10 different digits. Thus an output layer must in this case at least comprise 10 different nodes. To be able to cover other data sets with a different number of categories the number of output nodes must be a parameter, too.

How we indicate the belonging to a category numerically - by a probability number between "0" and "1" or just a "1" at the right category and zeros otherwise - can be a matter of discussion. It is also a question of the cost function we wish to use. We will come back to this point later.

The numbers of "hidden layers" and their nodes
We want these numbers to be parameters, too. For simple datasets we do not need big networks, but we want to able to play around a bit with 1 up to 3 layers. (For an ANN to recognize hand written digits an input layer "layer 1" and only one hidden layer "layer 2" before the output layer"layer 3" - are required.)

Activation and output functions
The nodes in hidden layers use a so called "activation function" to transform the aggregated input from different feeding nodes of the previous layer into one distinct value in a given interval - e.g. between 0 and 1. Again, we should be prepared to have a parameter to choose between different activation functions.

We should be aware of the fact that the nodes of the output layers need special consideration as the "activation function" there produces the final output - which must allow for a distinction of categories. This may lead to a special form - e.g. a kind of probability function. So, output functions should also be regarded as variable.

An ANN class and its interface

I develop my code as a Python module in an PyDEV/Eclipse which uses a virtual Python environment. I described the setup of such an environment in detail in XXXX. In the directory structure there I place a module "myann.py" at the location "...../ml_1/mynotebooks/mycode/myann.py". We need to import some libraries first

Modules and libraries to import

'''
Module to create a simple layered neural network for the MNIST data set
Created on 23.08.2019
@author: ramoe
'''
import numpy as np
import math 
import sys
import time
import tensorflow
from sklearn.datasets import fetch_mldata
from sklearn.datasets import fetch_openml
from keras.datasets import mnist as kmnist
from scipy.special import expit  
from matplotlib import pyplot as plt
#from matplotlib.colors import ListedColormap
#import matplotlib.patches as mpat 
#from keras.activations import relu

 

Why do I import "tensorflow" and "keras"?
Well, only for the purpose to create the input data of MNIST quickly. Sklearn's "fetchml_data" is doomed to end. "fetch_openml" does not use caching in some older versions and is also otherwise terribly slow. But, "keras", which in turn needs tensorflow as a backend, provides its own tool to provide the MNIST data.

We need "scipy" to get an optimized version of the "sigmoid"-function which is one important version of an activation function. "numpy" and "math" are required for fast array- and math-operations. "time" is required to measure the run time of program segments and "mathplotlib" will help to visualize some information during training.

The "__init__"-function of our class

We encapsulate most of the functionality in a class and its methods. Python provides the "__init__"-function, which we can use as a kind of "constructor" - although it technically is not the same as a constructor in other languages. Anyway, we can use it as an interface to provide parameters and initialize variables in a class instance.

We shall build up our "__init__"-function during the next articles step by step in form of well separated segments; in the beginning we only look a attributes and methods required to set up (MNIST) input data, to create the basic network layers and to create arrays with the (MNIST) input data.

Parameters

class MyANN:
    def __init__(self, 
                 my_data_set = "mnist", 
                 n_hidden_layers = 1, 
                 ay_nodes_layers = [0, 100, 0], # array which should have as much elements as n_hidden + 2
                 n_nodes_layer_out = 10,  # number of nodes in output layer 
                 
                 my_activation_function = "sigmoid", 
                 my_out_function        = "sigmoid",   
                 
                 vect_mode = 'cols', 
                 
                 figs_x1=12.0, figs_x2=8.0, 
                 legend_loc='upper right'
                 ):
        '''
        Initialization of MyANN
        Input: 
            data_set: type of dataset; so far only the "mnist", "mnist_784" and the "mnist_keras" datsets are known. 
                      We use this information to prepare the input data and learn about the feature dimension. 
                      This info is used in preparing the size of the input layer.     
            n_hidden_layers = number of hidden layers => between input layer 0 and output layer n 
            ay_nodes_layers = [0, 100, 0 ] : We set the number of nodes in input layer_0 and the output_layer to zero 
                                             Will be set to real number afterwards by infos from the input dataset. 
                                             All other numbers are used for the node numbers of the hidden layers.
            n_nodes_layer_out = expected number of nodes in the output layer (is checked); 
                                this number corresponds to the number of categories to be distinguished   
            
            my_activation_function : name of the activation function to use 
            my_out_function : name of the "activation" function of the last layer whcih produces the output values 
            
            vect_mode: Are 1-dim data arrays (vectors) ordered by columns or rows ?
            
            figs_x1=12.0, figs_x2=8.0 : Standard sizing of plots , 
            legend_loc='upper right': Position of legends in the plots 
            
         '''

 

You see that I defined multiple parameters, which are explained in the doc string. We use a name to choose the dataset to train our ANN on. We assume that special methods to fetch the required data are implemented in our class. This requires that the class knows exactly which data sets it is capable to handle. We provide an list with this information below. The other parameters should be clear from their explanation.

Initialization of class attributes

We first initialize a bunch of class attributes which we shall use to define the network of layers, nodes, weights, to keep our input data and functions.

        # Array (Python list) of known input data sets 
        self.__input_data_sets = ["mnist", "mnist_784", "mnist_keras"]  
        self._my_data_set = my_data_set
        
        # X, y, X_train, y_train, X_test, y_test  
            # will be set by analyze_input_data 
            # X: Input array (2D) - at present status of MNIST image data, only.    
            # y: result (=classification data) [digits represent categories in the case of Mnist]
        self._X       = None 
        self._X_train = None 
        self._X_test  = None   
        self._y       = None 
        self._y_train = None 
        self._y_test  = None
        
        # relevant dimensions 
        # from input data information;  will be set in handle_input_data()
        self._dim_sets     = 0  
        self._dim_features = 0  
        self._n_labels     = 0   # number of unique labels - will be extracted from y-data 
        
        # Img sizes 
        self._dim_img      = 0 # should be sqrt(dim_features) - we assume square like images  
        self._img_h        = 0 
        self._img_w        = 0 
        
        # Layers
        # ------
        # number of hidden layers 
        self._n_hidden_layers = n_hidden_layers
        # Number of total layers 
        self._n_total_layers = 2 + self._n_hidden_layers  
        # Nodes for hidden layers 
        self._ay_nodes_layers = np.array(ay_nodes_layers)
        # Number of nodes in output layer - will be checked against information from target arrays
        self._n_nodes_layer_out = n_nodes_layer_out
        
        
        # Weights 
        # --------
        # empty List for all weight-matrices for all layer-connections
        # Numbering : 
        # w[0] contains the weight matrix which connects layer 0 (input layer ) to hidden layer 1 
        # w[1] contains the weight matrix which connects layer 1 (input layer ) to (hidden?) layer 2 
        self._ay_w = []  
        
        # Known Randomizer methods ( 0: np.random.randint, 1: np.random.uniform )  
        # ------------------
        self.__ay_known_randomizers = [0, 1]

        # Types of activation functions and output functions 
        # ------------------
        self.__ay_activation_functions = ["sigmoid"] # later also relu 
        self.__ay_output_functions     = ["sigmoid"] # later also softmax 
        
        # the following dictionaries will be used for indirect function calls 
        self.__d_activation_funcs = {
            'sigmoid': self._sigmoid, 
            'relu':    self._relu
            
            }
        self.__d_output_funcs = { 
            'sigmoid': self._sigmoid, 
            'softmax': self._softmax
            }    

        # The following variables will later be set by _check_and set_activation_and_out_functions()            
        self._my_act_func = my_activation_function
        self._my_out_func = my_out_function
        self._act_func = None    
        self._out_func = None    

        # Plot handling 
        # --------------
        # Alternatives to resize plots 
        # 1: just resize figure  2: resize plus create subplots() [figure + axes] 
        self._plot_resize_alternative = 1 
        # Plot-sizing
        self._figs_x1 = figs_x1
        self._figs_x2 = figs_x2
        self._fig = None
        self._ax  = None 
                # alternative 2 does resizing and (!) subplots() 
        self.initiate_and_resize_plot(self._plot_resize_alternative)  

        # ***********
        # operations 
        # ***********
        # check and handle input data 
        self._handle_input_data()
        
        print("\nStopping program regularily")
        sys.exit()
       

 
To make things not more complicated as necessary in the beginning I omit the usage of properties and full encapsulation of private attributes. For convenience reasons I also use only one underscore for some attributes and functions/methods to allow for external usage. This is helpful in a testing phase. However, many items can in the end be switched to really private properties or methods.

List of known input datasets
The list of known input data sets is kept in the variable "self.__input_data_sets". The variables

self._X, self._X_train, self._X_test, self._y, self._y_train, self._y_test

will be used to keep all data of the chosen dataset, the training data, the test data for checking of the reliability of the algorithm after training, and the corresponding target data (y_...) for classification. The target data in the MNIST case contain the digit the respective image (X_..) represents.

All of the named attributes will become numpy arrays. A method called "_handle_input_data(self)" will load the (MNIST) input data and fill the numpy arrays.

The input arrays "X_..." will via their dimensions provide the information on the number of data sets (_dim_sets) and the number of features (_dim_features). The target data arrays "Y_..." provide the number of "classes" (MNIST: 10 digits) to distinguish between by the ANN algorithm. We keep this number in the attribute "_n_labels".

It is also useful to keep the pixel dimensions of input image data. At least for MNIST we assume quadratic images (_img_h = img_w = _dim_img).

Layers and weights
The number of total layers ("_n_total_layers") is by 2 bigger than the number of hidden layers (_n_hidden_layers).

We take the number of nodes in the layers from the respective list provided as an input parameter "ay_nodes_layers" to our class. We transform the list into a numpy array "_ay_nodes_layers". The expected number of nodes in the output layer is used for consistency checks and saved in "_n_nodes_layer_out".

The "weights" of a network must be given in form of matrices - as they describe connections between nodes of different adjacent layers. As the number of layers is not fixed but can be set by the user, I use a list "_ay_w" to collect such matrices in the order of layer_0 (input) to layer_n (output).

Weights, i.e. the matrix elements must initially be defined by random numbers. To provide such numbers we have to use randomizer functions. Depending on the kind (floating point numbers, integer numbers) we use at least two randomizers (randint, uniform).

Usable activation and output functions are named in Python dictionaries which point to respective methods. This allows for "indirect addressing" of these functions later on. You may recognize this by the direct refernece of the dictionary elements to defined class methods (no strings are used!).

For the time being we work with the "sigmoid" and the "relu" functions for activation and the "sigmoid" and "softmax" functions for output creation. The attributes "self._act_func" and "self._out_func" are used later on to invoke the functions requested by the respective parameters of the classes interface.

The final part of the code segment given above is used for plot sizing with the help of "matplotlib"; a method "initiate_and_resize_plot()" takes care of this. It can use 2 alternative ways of doing so.

Read and provide the input data

Now let us turn to some methods. We first need to read in and prepare the input data. We use a method "_handle_input_data()" to work on this problem.

For the time being we have only three different ways to load the MNIST dataset from different origins. We need not yet call different methods but deal with each MNIST source within the method:

    # Method to handle different types of input data sets 
    def _handle_input_data(self):    
        '''
        Method to deal with the input data: 
        - check if we have a known data set ("mnist" so far)
        - reshape as required 
        - analyze dimensions and extract the feature dimension(s) 
        '''
        # check for known dataset 
        try: 
            if (self._my_data_set not in self._input_data_sets ): 
                raise ValueError
        except ValueError:
            print("The requested input data" + self._my_data_set + " is not known!" )
            sys.exit()   

        
        # handle the mnist original dataset 
        if ( self._my_data_set == "mnist"): 
            mnist = fetch_mldata('MNIST original')
            self._X, self._y = mnist["data"], mnist["target"]
            print("Input data for dataset " + self._my_data_set + " : \n" + "Original shape of X = " + str(self._X.shape) +
              "\n" + "Original shape of y = " + str(self._y.shape))
            self._X_train, self._X_test, self._y_train, self._y_test = self._X[:60000], self._X[60000:], self._y[:60000], self._y[60000:] 
            
        # handle the mnist_784 dataset 
        if ( self._my_data_set == "mnist_784"): 
            mnist2 = fetch_openml('mnist_784', version=1, cache=True, data_home='~/scikit_learn_data') 
            self._X, self._y = mnist2["data"], mnist2["target"]
            print ("data fetched")
            # the target categories are given as strings not integers 
            self._y = np.array([int(i) for i in self._y])
            print ("data modified")
            print("Input data for dataset " + self._my_data_set + " : \n" + "Original shape of X = " + str(self._X.shape) +
              "\n" + "Original shape of y = " + str(self._y.shape))
            self._X_train, self._X_test, self._y_train, self._y_test = self._X[:60000], self._X[60000:], self._y[:60000], self._y[60000:] 

        # handle the mnist_keras dataset 
        if ( self._my_data_set == "mnist_keras"): 
            (self._X_train, self._y_train), (self._X_test, self._y_test) = kmnist.load_data()
            len_train =  self._X_train.shape[0]
            #print(len_train)
            print("Input data for dataset " + self._my_data_set + " : \n" + "Original shape of X_train = " + str(self._X_train.shape) +
              "\n" + "Original Shape of y_train = " + str(self._y_train.shape))
            len_test =  self._X_test.shape[0]
            #print(len_test)
            print("Original shape of X_test = " + str(self._X_test.shape) +
              "\n" + "Original Shape of y_test = " + str(self._y_test.shape))
            self._X_train = self._X_train.reshape(len_train, 28*28) 
            self._X_test  = self._X_test.reshape(len_test, 28*28) 

        # Common Mnist handling 
        if ( self._my_data_set == "mnist" or self._my_data_set == "mnist_784" or self._my_data_set == "mnist_keras" ): 
            self._common_handling_of_mnist()    
        
        # Other input data sets can not yet be handled 

 
We first check, whether the input parameter fits a known dataset - and raise an error if otherwise. The data come in different forms for the three sources of MNIST. For each set we want to extract arrays

self._X_train, self._X_test, self._y_train, self._y_test.

We have to do this a bit differently for the 3 cases. Note that the "mnist_784" set from "fetch_openml" gives the target
category values in form of strings and not integers. We correct this after loading.

The fastest method is the MNIST dataset based on "keras"; the keras function "kmnist.load_data()" provides already a 60000:10000 ratio for training and test data. However, we get the images in a (60000, 28, 28) array shape; we therefore reshape the "_X_train"-array to (60000, 784) and "_X_test"-array to (10000, 784).

Analyze the input data and encode the provided target labels

The further handling of the MNIST data requires some common analysis.

    # Method for common input data handling of Mnist data sets
    def _common_handling_of_mnist(self):    
        print("\nFinal input data for dataset " + self._my_data_set + 
              " : \n" + "Shape of X_train = " + str(self._X_train.shape) + 
              "\n" + "Shape of y_train = " + str(self._y_train.shape) + 
              "\n" + "Shape of X_test = " + str(self._X_test.shape) + 
              "\n" + "Shape of y_test = " + str(self._y_test.shape) 
              )
        
        # mixing the training indices 
        shuffled_index = np.random.permutation(60000)
        self._X_train, self._y_train = self._X_train[shuffled_index], self._y_train[shuffled_index]
        
        # set dimensions 
        self._dim_sets = self._y_train.shape[0]
        self._dim_features = self._X_train.shape[1] 
        self._dim_img = math.sqrt(self._dim_features)
        # we assume square images 
        self._img_h = int(self._dim_img)
        self._img_w = int(self._dim_img)
         
        # Print dimensions  
        print("\nWe have " + str(self._dim_sets) + " data sets for training") 
        print("Feature dimension is " + str(self._dim_features) + " (= " + str(self._img_w)+ "x" + str(self._img_h) + ")") 

        # we need to encode the digit labels of mnist 
        self._get_num_labels()
        self._encode_all_mnist_labels()

 
As you see we retrieve some of our class attributes which we shall use during training and do some printing. This is trivial. Not so trivial is the handling of the output data.

What shape do we expect for the "_X_train" and "_y_train"? Each element of the input data set is an array with values for all features. Thus the "_X_train.shape" should be (60000, 784). For _y_train we expect a simple integer describing the digit to which the MNIST input image corresponds. Thus we expect a one dimensional array with _y_train.shape = (60000).

However, the output data of our ANN for one input element come as an array of values for different categories - and not as a simple number. This shows that we need to encode the "_y_train"-data, i.e. the target labels, into a usable array form. We use two methods to achieve this:

    # Method to encode mnist labels 
    def _get_num_labels(self):
        self._n_labels = len(np.unique(self._y_train))
        print("The number of labels is " + str(self._n_labels))
        

    # Method to encode all mnist labels 
    def _encode_all_mnist_labels(self, b_print=True):
        '''
        We shall use vectorized input and output - i.e. we process a whole batch of input data sets in parallel
        (see article in the Linux blog) 
        The output array will then have the form OUT(i_out_node, idx) where 
            i_out_node enumerates the node of the last layer (i.e. the category)    
            idx enumerates the data set within a batch,
        After training, if y_train[idx] = 6, we would expect an output value of OUT[6,idx] = 1.0 and OUT[i_node, idx]=0.0 otherwise 
        for a categorization decision in the ideal case. Realistically, we will get a distribution of numbers over the nodes 
        with values between 0.0 and 1.0 - with hopefully the maximum value at the right node OUT[6,idx].  
        
        The following method creates an arrays OneHot[i_out_node, idx] with 
        OneHot[i_node_out, idx] = 1.0, if i_node_out =  int(y[idx])
        OneHot(i_node_out, idx] = 0.0, if i_node_out != int(y[idx])  
        
        This will allow for a vectorized comparison of calculated values and knwon values during training 
        '''
       
        self._ay_onehot = np.zeros((self._n_labels, self._y_train.shape[0]))
        # ay_oneval is just for convenience and printing purposes 
        self._ay_oneval = np.zeros((self._n_labels, self._y_train.shape[0], 2))
        
        if b_print: 
            print("\nShape of y_train = " + str(self._y_train.shape))
            print("Shape of ay_onehot = " + str(self._ay_onehot.shape))
        
        # the next block is just for illustration purposes and a better understanding
        if b_print: 
            values = enumerate(self._y_train[0:12])
            print("\nValues of the enumerate structure for the first 12 elements : = ")

            for iv in values: 
                print(iv)
        
        # here we prepare the array for vectorized comparison 
        print("\nLabels for the first 12 datasets:")
        for idx, val in enumerate(self._y_train):
            self._ay_onehot[val, idx ]   = 1.0
            self._ay_oneval[val, idx, 0] = 1.0 
            self._ay_oneval[val, idx, 1] = val 
       
        if b_print: 
            print("\nShape of ay_onehot = " + str(self._ay_onehot.shape))
            print(self._ay_onehot[:, 0:12])
            #print("Shape of ay_oneval = " + str(self._ay_oneval.shape))
            #print(self._ay_oneval[:, 0:12, :])

 
The first method only determines the number of labels (= number of categories).

We see from the code of the second method that we encode the target labels in the form of two arrays. The relevant one for our optimization algorithm will be "_ay_onehot". This array is 2-dimensional. Why?

Working with mini-batches

A big advantage of the optimization method we shall use later on during training is that we will perform weight adjustment for a whole bunch of training data in one step. Meaning:

We propagate a whole bunch of test data in parallel throughout the grid, get an array with result data (output array) for which we then calculate a value of the total cost function and an array containing the difference of the output array to an array of correct values (error) for all test data of the bunch. The "error" (i.e. the difference) will be back-propagated to determine the derivative of the cost function for corrections of the node weights.

Such a bunch is called a "batch" and if it is significantly smaller than the whole set of training data - a "mini-batch". Working with "mini-batches" during optimization is a compromise between using the full data set for gradient determination during each optimization step ("batch approach") or using just one input data element of the training set for gradient descent correction ("stochastic approach"). See chapter 4 of the book of Geron and chapter 2 in the book of Raschka for some information on this topic.

The advantage of mini-batches is that we can use vectorized linear algebra operations over all elements of the batch. Linear Algebra libraries are optimized to perform such operations on modern CPUs and GPUs.

You really should keep in mind for understanding the code for the propagation and optimization algorithms discussed in forthcoming articles that the cost function is determined over all elements of a batch and the derivative determination is a matrix like operation involving all input elements of each defined batch! This corresponds to the point that the separation interface in the feature space must be adjusted with respect to all given data points in the end. Mini-batches only help in so far as we look at selected samples of data to achieve bigger steps of our gradient guided descent into an optimum in the beginning - with the disadvantage of making some jumpy stochastic turns instead of a smoother approach.

What is the shape of the output array?
A single input element of the bunch is an array of 784 feature values. The corresponding output array is an array with values for 10 categories (here digits). But, what about a whole bunch of test data, i.e. a "batch"?

As I have explained already in my last article
Numpy matrix multiplication for layers of simple feed forward ANNs
that the output array for a bunch of test data will have the form "_ay_a_Out[i_out_node, idx]" with:

  • i_out_node enumerating the node of the last layer, i.e. the various possible category
  • idx enumerating the data set within a batch of training data

We shall construct the output function such that it provides values within the interval [0,1] for each node of the output layer. We define a perfectly working grid after training as one that would produce a "1" for the correct category node (i.e. the expected digit) and "0" otherwise.

For error-back-propagation we need to compare the real results with the correct category values for the test data. To be able to do this we must build up a 2-dim array of the same shape as "_ay_a_Out" with correct output values for all test data of the batch. E.g.: If we expect the digit 7 for an input array of index idx within the set of training data, we need a 2-dim output array with the element [[0,0,0,0,0,0,0,1,0,0], idx].

By using Numpy's zero()-function and Pythons "enumerate()"-function we achieve exactly this for all data elements of the training data set. See the method "_encode_all_mnist_labels()". Thus, _ay_onehot.shape is expected to be (10, 60000). From this 2-dim array we can later cut out bunches of consecutive test data for mini-batches.

We print out the values for the first 12 elements of the input data set.

The array "_ay_oneval" is provided for convenience and print purposes: it provides the expected digit value in addition.

Tests and output via a Jupyter notebook

Let us test the reading of the input data and the label encoding with a Jupyter notebook. In previous articles I have described already how to do so.

I build a Jupyter notebook called "myANN" (in my present working directory "/projekte/GIT/ai/ml1/mynotebooks").

I start it with

myself@mytux:/projekte/GIT/ai/ml1> source bin/activate 
(ml1) myself@mytux:/projekte/GIT/ai/ml1> jupyter notebook
[I 15:07:30.953 NotebookApp] Writing notebook server cookie secret to /run/user/21001/jupyter/notebook_cookie_secret
[I 15:07:38.754 NotebookApp] jupyter_tensorboard extension loaded.
[I 15:07:38.754 NotebookApp] Serving notebooks from local directory: /projekte/GIT/ai/ml1
[I 15:07:38.754 NotebookApp] The Jupyter Notebook is running at:
[I 15:07:38.754 NotebookApp] http://localhost:8888/?token=06c2626c8724f65d1e3c4a50457da0d6db414f88a40c7baf
[I 15:07:38.755 NotebookApp] Use Control-C to stop this server and shut down all kernels (twice to skip confirmation).
[C 15:07:38.771 NotebookApp] 

and add two cells. The first one is for imports of libraries.

By the last line I import my present class code. With the second cell I start my class; the "__init__"-function is automatically executed:

The output is:

Note that the display of "_ay_onehot" shows the categories in vertical direction (rows) and the index for the input data element in horizontal direction (columns)! You see that correspondence of the labels in the enumerate structure correspond to the "1"s in the "_ay_onehot"-array.

Conclusion

Reading the MNIST dataset into Numpy arrays via Keras is simple. We have prepared an array "_ay_onehot"- which we shall use during optimization to calculate a difference between calculated output of the ANN and the expected category for an element of our training data set.

In the next article
A simple program for an ANN to cover the Mnist dataset – II
we shall define initial weights for our ANN.

Literature and links

Referenced Books
"Python machine Learning", Seb. Raschka, 2016, Packt Publishing, Birmingham, UK
"Machine Learning mit Sckit-Learn & TensorFlow", A. Geron, 2018, O'REILLY, dpunkt.verlag GmbH, Heidelberg, Deutschland

Links regarding cost (or loss) functions and logistic regression
https://towardsdatascience.com/introduction-to-logistic-regression-66248243c148
https://cmci.colorado.edu/classes/INFO-4604/files/slides-5_logistic.pdf
Wikipedia article on Loss functions for classification
https://towardsdatascience.com/optimization-loss-function-under-the-hood-part-ii-d20a239cde11
https://stackoverflow.com/questions/32986123/why-the-cost-function-of-logistic-regression-has-a-logarithmic-expression
https://medium.com/technology-nineleaps/logistic-regression-gradient-descent-optimization-part-1-ed320325a67e
https://blog.algorithmia.com/introduction-to-loss-functions/
uni leipzig on logistic regression