A simple Python program for an ANN to cover the MNIST dataset – VIII – coding Error Backward Propagation

I continue with my series on a Python program for coding small “Multilayer Perceptrons” [MLPs].

A simple program for an ANN to cover the Mnist dataset – VII – EBP related topics and obstacles
A simple program for an ANN to cover the Mnist dataset – VI – the math behind the „error back-propagation“
A simple program for an ANN to cover the Mnist dataset – V – coding the loss function
A simple program for an ANN to cover the Mnist dataset – IV – the concept of a cost or loss function
A simple program for an ANN to cover the Mnist dataset – III – forward propagation
A simple program for an ANN to cover the Mnist dataset – II – initial random weight values
A simple program for an ANN to cover the Mnist dataset – I – a starting point

After all the theoretical considerations of the last two articles we now start coding again. Our objective is to extend our methods for training the MLP on the MNIST dataset by methods which perform the “error back propagation” and the correction of the weights. The mathematical prescriptions were derived in the following PDF:
My PDF on “The math behind EBP”

When you study the new code fragments below remember a few things:
We are prepared to use mini-batches. Therefore, the cost functions will be calculated over the data records of each batch and all the matrix operations for back propagation will cover all batch-records in parallel. Training means to loop over epochs and mini-batches – in pseudo-code:

  • Loop over epochs
    1. adjust learning rate,
    2. check for convergence criteria,
    3. Shuffle all data records in the test data set and build new mini-batches
  • Loop over mini-batches
    1. Perform forward propagation for all records of the mini-batch
    2. calculate and save the total cost value for each mini-batch
    3. calculate and save an averaged error on the output layer for each mini-batch
    4. perform error backward propagation on all records of the mini-batch to get the gradient of the cost function with respect to all weights
    5. adjust all weights on all layers

As discussed in the last article: The cost hyperplane changes a bit with each mini-batch. If there is a good mixture of records in a batch then the form of its specific cost hyperplane will (hopefully) resemble the form of an overall cost hyperplane, but it will not be the same. By the second step in the outer loop we want to avoid that the same data records always get an influence on the gradients at the same position in the correction procedure. Both statistical elements help a bit to overcome dominant records and a non-equal distribution of test records. If we had only pictures for number 3 at the end of our MNIST data set we
may start learning “3” representations very well, but not other numbers. Statistical variation also helps to avoid side minima on the overall cost hyperplane for all data records of the test set.

We shall implement the second step and third step in the epoch loop in the next article – when we are sure that the training algorithm works as expected. So, at the moment we will stop our training only after a given number of epochs.

More input parameters

In the first articles we had build an __init__()-method to parameterize a training run. We have to include three more parameters to control the backward propagation.

learn_rate = 0.001, # the learning rate (often called epsilon in textbooks)
decrease_const = 0.00001, # a factor for decreasing the learning rate with epochs
mom_rate = 0.0005, # a factor for momentum learning

The first parameter controls by how much we change weights with the help of gradient values. See formula (93) in the PDF of article VI (you find the Link to the latest version in the last section of this article). The second parameter will give us an option to decrease the learning rate with the number of training epochs. Note that a constant decrease rate only makes sense, if we can be relatively sure that we do not end up in a side minimum of the cost function.

The third parameter is interesting: It will allow us to mix the presently calculated weight correction with the correction term from the last step. So to say: We extend the “momentum” of the last correction into the next correction. This helps us not to follow indicated direction changes on the cost hyperplanes too fast.

Some hygienic measures regarding variables

In the already written parts of the code we have used a prefix “ay_” for all variables which represent some vector or array like structure – including Python lists and Numpy arrays. For back propagation coding it will be more important to distinguish between lists and arrays. So, I changed the variable prefix for some important Python lists from “ay_” to “li_”. (I shall do it for all lists used in a further version). In addition I have changed the prefix for Python ranges to “rg_”. These changes will affect the contents and the interface of some methods. You will notice when we come to these methods.

The changed __input__()-method

Our modified __init__() function now looks like this:

    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,  # expected number of nodes in output layer 
                                                  
                 my_activation_function = "sigmoid", 
                 my_out_function        = "sigmoid",   
                 my_loss_function       = "LogLoss",   
                 
                 n_size_mini_batch = 50,  # number of data elements in a mini-batch 
                 
                 n_epochs      = 1,
                 n_max_batches = -1,  # number of mini-batches to use during epochs - > 0 only for testing 
                                      # a negative value uses all mini-batches 
                 
                 lambda2_reg = 0.1,     # factor for quadratic regularization term 
                 lambda1_reg = 0.0,     # factor for linear regularization term 
                 
                 vect_mode = 'cols', 
                 
                 learn_rate = 0.001,        # the learning rate (often called epsilon in textbooks) 
                 decrease_const = 0.00001,  # a factor for decreasing the learning rate with epochs
                 mom_rate  
 = 0.0005,       # a factor for momentum learning
                 
                 figs_x1=12.0, figs_x2=8.0, 
                 legend_loc='upper right',
                 
                 b_print_test_data = True
                 
                 ):
        '''
        Initialization of MyANN
        Input: 
            data_set: type of dataset; so far only the "mnist", "mnist_784" 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_out_layer = expected number of nodes in the output layer (is checked); 
                                this number corresponds to the number of categories NC = number of labels 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 which produces the output values 
            my_loss_function : name of the "cost" or "loss" function used for optimization 
            
            n_size_mini_batch : Number of elements/samples in a mini-batch of training data 
                                The number of mini-batches will be calculated from this
            
            n_epochs : number of epochs to calculate during training
            n_max_batches : > 0: maximum of mini-batches to use during training 
                            < 0: use all mini-batches  
            
            lambda_reg2:    The factor for the quadartic regularization term 
            lambda_reg1:    The factor for the linear regularization term 
            
            vect_mode: Are 1-dim data arrays (vctors) ordered by columns or rows ?
            
            learn rate :     Learning rate - definies by how much we correct weights in the indicated direction of the gradient on the cost hyperplane.
            decrease_const:  Controls a systematic decrease of the learning rate with epoch number 
            mom_const:       Momentum rate. Controls a mixture of the last with the present weight corrections (momentum learning)
            
            figs_x1=12.0, figs_x2=8.0 : Standard sizing of plots , 
            legend_loc='upper right': Position of legends in the plots
            
            b_print_test_data: Boolean variable to control the print out of some tests data 
             
         '''
        
        # 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._li_w = []
        
        # Arrays for encoded output labels - will be set in _encode_all_mnist_labels()
        # -------------------------------
        self._ay_onehot = None
        self._ay_oneval = None
        
        # 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 
        
        # Types of cost functions 
        # ------------------
        self.__ay_loss_functions = ["LogLoss", "MSE" ] # later also othr types of cost/loss functions  


        # 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
            }  
        self.__d_loss_funcs = { 
            'LogLoss': self._loss_LogLoss, 
            'MSE': self._loss_MSE
            }  
        # Derivative functions 
        self.__d_D_activation_funcs = {
            'sigmoid': self._D_sigmoid, 
            'relu':    self._D_relu
            }
        self.__d_D_output_funcs = { 
            'sigmoid': self._D_sigmoid, 
            'softmax': self._D_softmax
            }  
        self.__d_D_loss_funcs = { 
            'LogLoss': self._D_loss_LogLoss, 
            'MSE': self._D_loss_MSE
            }  
        
        
        # 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._my_loss_func = my_loss_function
        self._act_func = None    
        self._out_func = None    
        self._loss_func = None    
        
        # number of data samples in a mini-batch 
        self._n_size_mini_batch = n_size_mini_batch
        self._n_mini_batches = None  # will be determined by _get_number_of_mini_batches()

        # maximum number of epochs - we set this number to an assumed maximum 
        # - as we shall build a backup and reload functionality for training, this should not be a major problem 
        self._n_epochs = n_epochs
        
        # maximum number of batches to handle ( if < 0 => all!) 
        self._n_max_batches = n_max_batches
        # actual number of batches 
        self._n_batches = None

        # regularization parameters
        self._lambda2_reg = lambda2_reg
        self._
lambda1_reg = lambda1_reg
        
        # parameter for momentum learning 
        self._learn_rate = learn_rate
        self._decrease_const = decrease_const
        self._mom_rate   = mom_rate
        self._li_mom = [None] *  self._n_total_layers
        
        # book-keeping for epochs and mini-batches 
        # -------------------------------
        # range for epochs - will be set by _prepare-epochs_and_batches() 
        self._rg_idx_epochs = None
        # range for mini-batches 
        self._rg_idx_batches = None
        # dimension of the numpy arrays for book-keeping - will be set in _prepare_epochs_and_batches() 
        self._shape_epochs_batches = None    # (n_epochs, n_batches, 1) 

        # list for error values at outermost layer for minibatches and epochs during training
        # we use a numpy array here because we can redimension it
        self._ay_theta = None
        # list for cost values of mini-batches during training 
        # The list will later be split into sections for epochs 
        self._ay_costs = None
        
        
        # Data elements for back propagation
        # ----------------------------------
        
        # 2-dim array of partial derivatives of the elements of an additive cost function 
        # The derivative is taken with respect to the output results a_j = ay_ANN_out[j]
        # The array dimensions account for nodes and sampls of a mini_batch. The array will be set in function 
        # self._initiate_bw_propagation()
        self._ay_delta_out_batch = None
        

        # parameter to allow printing of some test data 
        self._b_print_test_data = b_print_test_data

        # 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()
        # set the ANN structure 
        self._set_ANN_structure()
        
        # Prepare epoch and batch-handling - sets ranges, limits num of mini-batches and initializes book-keeping arrays
        self._rg_idx_epochs, self._rg_idx_batches = self._prepare_epochs_and_batches()
        
        # perform training 
        start_c = time.perf_counter()
        self._fit(b_print=True, b_measure_batch_time=False)
        end_c = time.perf_counter()
        print('\n\n ------') 
        print('Total training Time_CPU: ', end_c - start_c) 
        print("\nStopping program regularily")
        sys.exit()

The extended method _set_ANN_structure()

I do not change method “_handle_input_data()”. However, I extend method “def _set_ANN_structure()” by a statement to initialize a list with momentum matrices for all layers.

    '''-- Method to set ANN structure --''' 
    def _set_ANN_structure(self):
        # check consistency of the node-number list with the number of hidden layers (n_hidden)
        self._check_layer_and_node_numbers()
        # set node numbers for the input layer and the output layer
        self._set_nodes_for_input_output_layers() 
        self._show_node_numbers() 
        
        # create the weight matrix between input and first hidden layer 
        self._create_WM_Input() 
        # create weight matrices between the 
hidden layers and between tha last hidden and the output layer 
        self._create_WM_Hidden() 
        
        # initialize momentum differences
        self._create_momentum_matrices()
        #print("\nLength of li_mom = ", str(len(self._li_mom)))
        
        # check and set activation functions 
        self._check_and_set_activation_and_out_functions()
        self._check_and_set_loss_function()
        
        return None

 
The following box shows the changed functions _create_WM_Input(), _create_WM_Hidden() and the new function _create_momentum_matrices():

    '''-- Method to create the weight matrix between L0/L1 --'''
    def _create_WM_Input(self):
        '''
        Method to create the input layer 
        The dimension will be taken from the structure of the input data 
        We need to fill self._w[0] with a matrix for conections of all nodes in L0 with all nodes in L1
        We fill the matrix with random numbers between [-1, 1] 
        '''
        # the num_nodes of layer 0 should already include the bias node 
        num_nodes_layer_0 = self._ay_nodes_layers[0]
        num_nodes_with_bias_layer_0 = num_nodes_layer_0 + 1 
        num_nodes_layer_1 = self._ay_nodes_layers[1] 
        
        # fill the matrix with random values 
        #rand_low  = -1.0
        #rand_high = 1.0
        rand_low  = -0.5
        rand_high = 0.5
        rand_size = num_nodes_layer_1 * (num_nodes_with_bias_layer_0) 
        
        randomizer = 1 # method np.random.uniform   
        
        w0 = self._create_vector_with_random_values(rand_low, rand_high, rand_size, randomizer)
        w0 = w0.reshape(num_nodes_layer_1, num_nodes_with_bias_layer_0)
        
        # put the weight matrix into array of matrices 
        self._li_w.append(w0)
        print("\nShape of weight matrix between layers 0 and 1 " + str(self._li_w[0].shape))
        
#
    '''-- Method to create the weight-matrices for hidden layers--''' 
    def _create_WM_Hidden(self):
        '''
        Method to create the weights of the hidden layers, i.e. between [L1, L2] and so on ... [L_n, L_out] 
        We fill the matrix with random numbers between [-1, 1] 
        '''
        
        # The "+1" is required due to range properties ! 
        rg_hidden_layers = range(1, self._n_hidden_layers + 1, 1)

        # for random operation 
        rand_low  = -1.0
        rand_high = 1.0
        
        for i in rg_hidden_layers: 
            print ("Creating weight matrix for layer " + str(i) + " to layer " + str(i+1) )
            
            num_nodes_layer = self._ay_nodes_layers[i] 
            num_nodes_with_bias_layer = num_nodes_layer + 1 
            
            # the number of the next layer is taken without the bias node!
            num_nodes_layer_next = self._ay_nodes_layers[i+1]
            
            # assign random values  
            rand_size = num_nodes_layer_next * num_nodes_with_bias_layer   
            
            randomizer = 1 # np.random.uniform
            
            w_i_next = self._create_vector_with_random_values(rand_low, rand_high, rand_size, randomizer)   
            w_i_next = w_i_next.reshape(num_nodes_layer_next, num_nodes_with_bias_layer)
            
            # put the weight matrix into our array of matrices 
            self._li_w.append(w_i_next)
            print("Shape of weight matrix between layers " + str(i) + " and " + str(i+1) + " = " + str(self._li_w[i].shape))
#
    '''-- Method to create and initialize matrices for momentum learning (differences) '''
    def _create_momentum_matrices(self):
        rg_layers = range(0, self._n_total_layers - 1)
        for i in rg_layers: 
r
            self._li_mom[i] = np.zeros(self._li_w[i].shape)
            #print("shape of li_mom[" + str(i) + "] = ", self._li_mom[i].shape)

 

The modified functions _fit() and _handle_mini_batch()

The _fit()-function is modified to include a systematic decrease of the learning rate.

    ''' -- Method to set the number of batches based on given batch size -- '''
    def _fit(self, b_print = False, b_measure_batch_time = False):
        
        rg_idx_epochs  = self._rg_idx_epochs 
        rg_idx_batches = self._rg_idx_batches
        if (b_print):    
            print("\nnumber of epochs = " + str(len(rg_idx_epochs)))
            print("max number of batches = " + str(len(rg_idx_batches)))
       
        # loop over epochs
        for idxe in rg_idx_epochs:
            if (b_print):
                print("\n ---------")
                print("\nStarting epoch " + str(idxe+1))
                
                self._learn_rate /= (1.0 + self._decrease_const * idxe)
            
            # loop over mini-batches
            for idxb in rg_idx_batches:
                if (b_print):
                    print("\n ---------")
                    print("\nDealing with mini-batch " + str(idxb+1))
                if b_measure_batch_time: 
                    start_0 = time.perf_counter()
                # deal with a mini-batch
                self._handle_mini_batch(num_batch = idxb, num_epoch=idxe, b_print_y_vals = False, b_print = False)
                if b_measure_batch_time: 
                    end_0 = time.perf_counter()
                    print('Time_CPU for batch ' + str(idxb+1), end_0 - start_0) 
                
                #if idxb == 100: 
                #    sys.exit() 
        
        return None

 
Note that the number of epochs is determined by an external parameter as an upper limit of the range “rg_idx_epochs”.

Method “_handle_mini_batch()” requires several changes: First we define lists which are required to save matrix data of the backward propagation. And, of course, we call a method to perform the BW propagation (see step 6 in the code). Some statements print shapes, if required. At step 7 of the code we correct the weights by using the learning rate and the calculated gradient of the loss function.

Note, that we mix the correction evaluated at the last batch-record with the correction evaluated for the present record! This corresponds to a simple form of momentum learning. We then have to save the present correction values, of course. Note that the list for momentum correction “li_mom” is, therefore, not deleted at the end of a mini-batch treatment !

In addition to saving the total costs per mini-batch we now also save a mean error at the output level. The average is done by the help of Numpy’s function numpy.average() for matrices. Remember, we build the average over errors at all output nodes and all records of the mini-batch.


    ''' -- Method to deal with a batch -- '''
    def _handle_mini_batch(self, num_batch = 0, num_epoch = 0, b_print_y_vals = False, b_print = False, b_keep_bw_matrices = True):
        '''
        For each batch we keep the input data array Z and the output data A (output of activation function!) 
        for all layers in Python lists
        We can use this as input variables in function calls - mutable variables are handled by reference values !
        We receive the A and Z data from propagation functions and proceed them to cost and gradient calculation functions
        
        As an initial step we define the Python lists li_Z_
in_layer and li_A_out_layer 
        and fill in the first input elements for layer L0  
        
        Forward propagation:
        --------------------
        li_Z_in_layer : List of layer-related 2-dim matrices for input values z at each node (rows) and all batch-samples (cols).
        li_A_out_layer: List of layer-related 2-dim matrices for output alues z at each node (rows) and all batch-samples (cols).
                        The output is created by Phi(z), where Phi represents an activation or output function 
        
        Note that the matrices in ay_A_out will be permanently extended by a row (over all samples) 
        to account for a bias node of each inner layer. This happens during FW propagation. 
        
        Note that the matrices ay_Z_in will be temporarily extended by a row (over all samples) 
        to account for a bias node of each inner layer. This happens during BW propagation.
        
        Backward propagation:
        -------------------- 
        li_delta_out:  Startup matrix for _out_delta-values at the outermost layer 
        li_grad_layer: List of layer-related matrices with gradient values for the correction of the weights 
        
        Depending on parameter "b_keep_bw_matrices" we keep 
            - a list of layer-related matrices D with values for the derivatives of the act./output functions
            - a list of layer-related matrices for the back propagated delta-values 
        in lists during back propagation. This can support error analysis. 
        
        All matrices in the lists are 2 dimensional with dimensions for nodes (rows) and training samples (cols) 
        All these lists be deleted at the end of the function to accelerate garbadge handling
        
        Input parameters: 
        ----------------
        num_epoch:     Number of present epoch
        num_batch:    Number of present mini-batch 
        '''
        # Layer-related lists to be filled with 2-dim Numpy matrices during FW propagation
        # ********************************************************************************
        li_Z_in_layer  = [None] * self._n_total_layers # List of matrices with z-input values for each layer; filled during FW-propagation
        li_A_out_layer = li_Z_in_layer.copy()          # List of matrices with results of activation/output-functions for each layer; filled during FW-propagation
        li_delta_out   = li_Z_in_layer.copy()          # Matrix with out_delta-values at the outermost layer 
        li_delta_layer = li_Z_in_layer.copy()          # List of the matrices for the BW propagated delta values 
        li_D_layer     = li_Z_in_layer.copy()          # List of the derivative matrices D containing partial derivatives of the activation/ouput functions 
        li_grad_layer  = li_Z_in_layer.copy()          # List of the matrices with gradient values for weight corrections
        
        if b_print: 
            len_lists = len(li_A_out_layer)
            print("\nnum_epoch = ", num_epoch, "  num_batch = ", num_batch )
            print("\nhandle_mini_batch(): length of lists = ", len_lists)
            self._info_point_print("handle_mini_batch: point 1")
        
        # Print some infos
        # ****************
        if b_print:
            self._print_batch_infos()
            self._info_point_print("handle_mini_batch: point 2")
        
        # Major steps for the mini-batch during one epoch iteration 
        # **********************************************************
        
        # Step 0: List of indices for data records in the present mini-batch
        # ******
        ay_idx_batch = self._ay_mini_batches[num_batch]
        
        # Step 1: Special preparation of the Z-input to the MLP's input Layer L0
        # ******
        # Layer L0: Fill in the input vector for the ANN's input layer L0 
        li_
Z_in_layer[0] = self._X_train[ay_idx_batch] # numpy arrays can be indexed by an array of integers
        if b_print:
            print("\nPropagation : Shape of X_in = li_Z_in_layer = ", li_Z_in_layer[0].shape)           
            #print("\nidx, expected y_value of Layer L0-input :")           
            #for idx in self._ay_mini_batches[num_batch]:
            #    print(str(idx) + ', ' + str(self._y_train[idx]) )
            self._info_point_print("handle_mini_batch: point 3")
        
        # Step 2: Layer L0: We need to transpose the data of the input layer 
        # *******
        ay_Z_in_0T       = li_Z_in_layer[0].T
        li_Z_in_layer[0] = ay_Z_in_0T
        if b_print:
            print("\nPropagation : Shape of transposed X_in = li_Z_in_layer = ", li_Z_in_layer[0].shape)           
            self._info_point_print("handle_mini_batch: point 4")
        
        # Step 3: Call forward propagation method for the present mini-batch of training records
        # *******
        # this function will fill the ay_Z_in- and ay_A_out-lists with matrices per layer
        self._fw_propagation(li_Z_in = li_Z_in_layer, li_A_out = li_A_out_layer, b_print = b_print) 
        
        if b_print:
            ilayer = range(0, self._n_total_layers)
            print("\n ---- ")
            print("\nAfter propagation through all " + str(self._n_total_layers) + " layers: ")
            for il in ilayer:
                print("Shape of Z_in of layer L" + str(il) + " = " + str(li_Z_in_layer[il].shape))
                print("Shape of A_out of layer L" + str(il) + " = " + str(li_A_out_layer[il].shape))
                if il < self._n_total_layers-1:
                    print("Shape of W of layer L" + str(il) + " = " + str(self._li_w[il].shape))
                    print("Shape of Mom of layer L" + str(il) + " = " + str(self._li_mom[il].shape))
            self._info_point_print("handle_mini_batch: point 5")
        
        
        # Step 4: Cost calculation for the mini-batch 
        # ********
        ay_y_enc = self._ay_onehot[:, ay_idx_batch]
        ay_ANN_out = li_A_out_layer[self._n_total_layers-1]
        # print("Shape of ay_ANN_out = " + str(ay_ANN_out.shape))
        
        total_costs_batch = self._calculate_loss_for_batch(ay_y_enc, ay_ANN_out, b_print = False)
        # we add the present cost value to the numpy array 
        self._ay_costs[num_epoch, num_batch] = total_costs_batch
        if b_print:
            print("\n total costs of mini_batch = ", self._ay_costs[num_epoch, num_batch])
            self._info_point_print("handle_mini_batch: point 6")
        print("\n total costs of mini_batch = ", self._ay_costs[num_epoch, num_batch])
        
        # Step 5: Avg-error for later plotting 
        # ********
        # mean "error" values - averaged over all nodes at outermost layer and all data sets of a mini-batch 
        ay_theta_out = ay_y_enc - ay_ANN_out
        if (b_print): 
            print("Shape of ay_theta_out = " + str(ay_theta_out.shape))
        ay_theta_avg = np.average(np.abs(ay_theta_out)) 
        self._ay_theta[num_epoch, num_batch] = ay_theta_avg 
        
        if b_print:
            print("\navg total error of mini_batch = ", self._ay_theta[num_epoch, num_batch])
            self._info_point_print("handle_mini_batch: point 7")
        print("avg total error of mini_batch = ", self._ay_theta[num_epoch, num_batch])
        
        
        # Step 6: Perform gradient calculation via back propagation of errors
        # ******* 
        self._bw_propagation( ay_y_enc = ay_y_enc, 
                              li_Z_in = li_Z_in_layer, 
                              li_A_out = li_A_out_layer, 
                              li_delta_out = li_delta_out, 
                              li_delta = li_delta_
layer,
                              li_D = li_D_layer, 
                              li_grad = li_grad_layer, 
                              b_print = b_print,
                              b_internal_timing = False 
                              ) 
        
        
        # Step 7: Adjustment of weights  
        # *******
        rg_layer=range(0, self._n_total_layers -1)
        for N in rg_layer:
            delta_w_N = self._learn_rate * li_grad_layer[N]
            self._li_w[N] -= ( delta_w_N + (self._mom_rate * self._li_mom[N]) )
            # save momentum
            self._li_mom[N] = delta_w_N
        
        # try to accelerate garbage handling
        # **************
        if len(li_Z_in_layer) > 0:
            del li_Z_in_layer
        if len(li_A_out_layer) > 0:
            del li_A_out_layer
        if len(li_delta_out) > 0:
            del li_delta_out
        if len(li_delta_layer) > 0:
            del li_delta_layer
        if len(li_D_layer) > 0:
            del li_D_layer
        if len(li_grad_layer) > 0:
            del li_grad_layer
            
        return None

 

Forward Propagation

The method for forward propagation remains unchanged in its structure. We only changed the prefix for the Python lists.

    ''' -- Method to handle FW propagation for a mini-batch --'''
    def _fw_propagation(self, li_Z_in, li_A_out, b_print= False):
        
        b_internal_timing = False
        
        # index range of layers 
        #    Note that we count from 0 (0=>L0) to E L(=>E) / 
        #    Careful: during BW-propgation we may need a correct indexing of lists filled during FW-propagation
        ilayer = range(0, self._n_total_layers-1)

        # propagation loop
        # ***************
        for il in ilayer:
            if b_internal_timing: start_0 = time.perf_counter()
            
            if b_print: 
                print("\nStarting propagation between L" + str(il) + " and L" + str(il+1))
                print("Shape of Z_in of layer L" + str(il) + " (without bias) = " + str(li_Z_in[il].shape))
            
            # Step 1: Take input of last layer and apply activation function 
            # ******
            if il == 0: 
                A_out_il = li_Z_in[il] # L0: activation function is identity 
            else: 
                A_out_il = self._act_func( li_Z_in[il] ) # use real activation function 
            
            # Step 2: Add bias node
            # ****** 
            A_out_il = self._add_bias_neuron_to_layer(A_out_il, 'row')
            # save in array     
            li_A_out[il] = A_out_il
            if b_print: 
                print("Shape of A_out of layer L" + str(il) + " (with bias) = " + str(li_A_out[il].shape))
            
            # Step 3: Propagate by matrix operation
            # ****** 
            Z_in_ilp1 = np.dot(self._li_w[il], A_out_il) 
            li_Z_in[il+1] = Z_in_ilp1
            
            if b_internal_timing: 
                end_0 = time.perf_counter()
                print('Time_CPU for layer propagation L' + str(il) + ' to L' + str(il+1), end_0 - start_0) 
        
        # treatment of the last layer 
        # ***************************
        il = il + 1
        if b_print:
            print("\nShape of Z_in of layer L" + str(il) + " = " + str(li_Z_in[il].shape))
        A_out_il = self._out_func( li_Z_in[il] ) # use the output function 
        li_A_out[il] = A_out_il
        if b_print:
            print("Shape of A_out of last layer L" + str(il) + " = " + str(li_A_out[il].shape))
        
        return None

 
nAddendum, 15.05.2020:
We shall later learn that the treatment of bias neurons can be done more efficiently. The present way of coding it reduces performance – especially at the input layer. See the article series starting with
MLP, Numpy, TF2 – performance issues – Step I – float32, reduction of back propagation
for more information. At the present stage of our discussion we are, however, more interested in getting a working code first – and not so much in performance optimization.

Methods for Error Backward Propagation

In contrast to the recipe given in my PDF on the EBP-math we cannot calculate the matrices with the derivatives of the activation functions “ay_D” in advance for all layers. The reason was discussed in the last article VII: Some matrices have to be intermediately adjusted for a bias-neuron, which is ignored in the analysis of the PDF.

The resulting code of our method for EBP looks like given below:

 
    ''' -- Method to handle error BW propagation for a mini-batch --'''
    def _bw_propagation(self, 
                        ay_y_enc, li_Z_in, li_A_out, 
                        li_delta_out, li_delta, li_D, li_grad, 
                        b_print = True, b_internal_timing = False):
        
        # List initialization: All parameter lists or arrays are filled or to be filled by layers 
        # Note: the lists li_Z_in, li_A_out were already filled by _fw_propagation() for the present batch 
        
        # Initiate BW propagation - provide delta-matrices for outermost layer
        # *********************** 
        # Input Z at outermost layer E  (4 layers -> layer 3)
        ay_Z_E = li_Z_in[self._n_total_layers-1]
        # Output A at outermost layer E (was calculated by output function)
        ay_A_E = li_A_out[self._n_total_layers-1]
        
        # Calculate D-matrix (derivative of output function) at outmost the layer - presently only D_sigmoid 
        ay_D_E = self._calculate_D_E(ay_Z_E=ay_Z_E, b_print=b_print )
        
        # Get the 2 delta matrices for the outermost layer (only layer E has 2 delta-matrices)
        ay_delta_E, ay_delta_out_E = self._calculate_delta_E(ay_y_enc=ay_y_enc, ay_A_E=ay_A_E, ay_D_E=ay_D_E, b_print=b_print) 
        
        # We check the shapes
        shape_theory = (self._n_nodes_layer_out, self._n_size_mini_batch)
        if (b_print and ay_delta_E.shape != shape_theory):
            print("\nError: Shape of ay_delta_E is wrong:")
            print("Shape = ", ay_delta_E.shape, "  ::  should be = ", shape_theory )
        if (b_print and ay_D_E.shape != shape_theory):
            print("\nError: Shape of ay_D_E is wrong:")
            print("Shape = ", ay_D_E.shape, "  ::  should be = ", shape_theory )
        
        # add the matrices to their lists ; li_delta_out gets only one element 
        idxE = self._n_total_layers - 1
        li_delta_out[idxE] = ay_delta_out_E # this happens only once
        li_delta[idxE]     = ay_delta_E
        li_D[idxE]         = ay_D_E
        li_grad[idxE]      = None    # On the outermost layer there is no gradient ! 
        
        if b_print:
            print("bw: Shape delta_E = ", li_delta[idxE].shape)
            print("bw: Shape D_E = ", ay_D_E.shape)
            self._info_point_print("bw_propagation: point bw_1")
        
        
        # Loop over all layers in reverse direction 
        # ******************************************
        # index range of target layers N in BW direction (starting with E-1 => 4 layers -> layer 2))
        if b_print:
            range_N_bw_layer_test = reversed(range(0,
 self._n_total_layers-1))   # must be -1 as the last element is not taken 
            rg_list = list(range_N_bw_layer_test) # Note this exhausts the range-object
            print("range_N_bw_layer = ", rg_list)
        
        range_N_bw_layer = reversed(range(0, self._n_total_layers-1))   # must be -1 as the last element is not taken 
        
        # loop over layers 
        for N in range_N_bw_layer:
            if b_print:
                print("\n N (layer) = " + str(N) +"\n")
            # start timer 
            if b_internal_timing: start_0 = time.perf_counter()
            
            # Back Propagation operations between layers N+1 and N 
            # *******************************************************
            # this method handles the special treatment of bias nodes in Z_in, too
            ay_delta_N, ay_D_N, ay_grad_N = self._bw_prop_Np1_to_N( N=N, li_Z_in=li_Z_in, li_A_out=li_A_out, li_delta=li_delta, b_print=False )
            
            if b_internal_timing: 
                end_0 = time.perf_counter()
                print('Time_CPU for BW layer operations ', end_0 - start_0) 
            
            # add matrices to their lists 
            li_delta[N] = ay_delta_N
            li_D[N]     = ay_D_N
            li_grad[N]= ay_grad_N
            #sys.exit()
        
        return

 

We first handle the necessary matrix evaluations for the outermost layer. We use two helper functions there to calculate the derivative of the output function with respect to the a-term [ _calculate_D_E() ] and to calculate the values for the “delta“-terms at all nodes and for all records [ _calculate_delta_E() ] according to the prescription in the PDF:

    
    ''' -- Method to calculate the matrix with the derivative values of the output function at outermost layer '''
    def _calculate_D_E(self, ay_Z_E, b_print= True):
        '''
        This method calculates and returns the D-matrix for the outermost layer
        The D matrix contains derivatives of the output function with respect to local input "z_j" at outermost nodes. 
        
        Returns
        ------
        ay_D_E:    Matrix with derivative values of the output function 
                   with respect to local z_j valus at the nodes of the outermost layer E
        Note: This is a 2-dim matrix over layer nodes and training samples of the mini-batch
        '''
        if self._my_out_func == 'sigmoid':
            ay_D_E = self._D_sigmoid(ay_Z=ay_Z_E)
        
        else:
            print("The derivative for output function " + self._my_out_func + " is not known yet!" )
            sys.exit()
        
        return ay_D_E

    ''' -- Method to calculate the delta_E matrix as a starting point of the backward propagation '''
    def _calculate_delta_E(self, ay_y_enc, ay_A_E, ay_D_E, b_print= False):
        '''
        This method calculates and returns the 2 delta-matrices for the outermost layer 
        
        Returns
        ------
        delta_E:     delta_matrix of the outermost layer (indicated by E)
        delta_out:   delta_out matrix => elements are local derivative values of the cost function 
                     with respect to the output "a_j" at an outermost node  
                     !!! delta_out will only be returned if calculable !!!
        
        Note: these are 2-dim matrices over layer nodes and training samples of the mini-batch
        '''
        
        if self._my_loss_func == 'LogLoss':
            # Calculate delta_S_E directly to avoid problems with zero denominators
            ay_delta_E = ay_A_E - ay_y_enc
            # delta_out is fetched but may be None 
            ay_delta_out, ay_D_
numerator, ay_D_denominator = self._D_loss_LogLoss(ay_y_enc, ay_A_E, b_print = False)
            
            # To be done: Analyze critical values in D_denominator 
            
            # Release variables explicitly 
            del ay_D_numerator
            del ay_D_denominator
            
        
        if self._my_loss_func == 'MSE':
            # First calculate delta_out and then the delta_E
            delta_out = self._D_loss_MSE(ay_y_enc, ay_A_E, b_print=False)
            # calculate delta_E via matrix multiplication 
            ay_delta_E = delta_out * ay_D_E
                    
        return ay_delta_E, ay_delta_out

 

Further required helper methods to calculate the cost functions and related derivatives are :

    ''' method to calculate the logistic regression loss function '''
    def _loss_LogLoss(self, ay_y_enc, ay_ANN_out, b_print = False):
        '''
        Method which calculates LogReg loss function in a vectorized form on multidimensional Numpy arrays 
        '''
        b_test = False

        if b_print:
            print("From LogLoss: shape of ay_y_enc =  " + str(ay_y_enc.shape))
            print("From LogLoss: shape of ay_ANN_out =  " + str(ay_ANN_out.shape))
            print("LogLoss: ay_y_enc = ", ay_y_enc) 
            print("LogLoss: ANN_out = \n", ay_ANN_out) 
            print("LogLoss: log(ay_ANN_out) =  \n", np.log(ay_ANN_out) )

        # The following means an element-wise (!) operation between matrices of the same shape!
        Log1 = -ay_y_enc * (np.log(ay_ANN_out))
        # The following means an element-wise (!) operation between matrices of the same shape!
        Log2 = (1 - ay_y_enc) * np.log(1 - ay_ANN_out)
        
        # the next operation calculates the sum over all matrix elements 
        # - thus getting the total costs for all mini-batch elements 
        cost = np.sum(Log1 - Log2)
        
        #if b_print and b_test:
            # Log1_x = -ay_y_enc.dot((np.log(ay_ANN_out)).T)
            # print("From LogLoss: L1 =   " + str(L1))
            # print("From LogLoss: L1X =  " + str(L1X))
        
        if b_print: 
            print("From LogLoss: cost =  " + str(cost))
        
        # The total costs is just a number (scalar)
        return cost
#
    ''' method to calculate the derivative of the logistic regression loss function 
        with respect to the output values '''
    def _D_loss_LogLoss(self, ay_y_enc, ay_ANN_out, b_print = False):
        '''
        This function returns the out_delta_S-matrix which is required to initialize the 
        BW propagation (EBP) 
        Note ANN_out is the A_out-list element ( a 2-dim matrix) for the outermost layer 
        In this case we have to take care of denominators = 0 
        '''
        D_numerator = ay_ANN_out - ay_y_enc
        D_denominator = -(ay_ANN_out - 1.0) * ay_ANN_out
        n_critical = np.count_nonzero(D_denominator < 1.0e-8)
        if n_critical > 0:
            delta_s_out = None
        else:
            delta_s_out = np.divide(D_numerator, D_denominator)
        return delta_s_out, D_numerator, D_denominator
#
    ''' method to calculate the MSE loss function '''
    def _loss_MSE(self, ay_y_enc, ay_ANN_out, b_print = False):
        '''
        Method which calculates LogReg loss function in a vectorized form on multidimensional Numpy arrays 
        '''
        if b_print:
            print("From loss_MSE: shape of ay_y_enc =  " + str(ay_y_enc.shape))
            print("From loss_MSE: shape of ay_ANN_out =  " + str(ay_ANN_out.shape))
            #print("LogReg: ay_y_enc = ", ay_y_enc) 
            #print("LogReg: ANN_out = \n", ay_
ANN_out) 
            #print("LogReg: log(ay_ANN_out) =  \n", np.log(ay_ANN_out) )
        
        cost = 0.5 * np.sum( np.square( ay_y_enc - ay_ANN_out ) )

        if b_print: 
            print("From loss_MSE: cost =  " + str(cost))
        
        return cost
#
    ''' method to calculate the derivative of the MSE loss function 
        with respect to the output values '''
    def _D_loss_MSE(self, ay_y_enc, ay_ANN_out, b_print = False):
        '''
        This function returns the out_delta_S - matrix which is required to initialize the 
        BW propagation (EBP) 
        Note ANN_out is the A_out-list element ( a 2-dim matrix) for the outermost layer
        In this case the output is harmless (no critical denominator) 
        '''
        delta_s_out = ay_ANN_out - ay_y_enc
        return delta_s_out

 
You see that we are a bit careful to avoid zero denominators for the Logarithmic loss function in all of our helper functions.

The check statements for shapes can be eliminated in a future version when we are sure that everything works correctly. Keeping the layer specific matrices during the handling of a mini-batch will be also good for potentially required error analysis in the beginning. In the end we only may keep the gradient-matrices and the layer specific matrices required to process the local calculations during back propagation.

Then we turn to loop over all other layers down to layer L0. The matrix operation to be done for all these layers are handled in a further method:

    
    ''' -- Method to calculate the BW-propagated delta-matrix and the gradient matrix to/for layer N '''
    def _bw_prop_Np1_to_N(self, N, li_Z_in, li_A_out, li_delta, b_print=False):
        '''
        BW-error-propagation bewtween layer N+1 and N 
        Inputs: 
            li_Z_in:  List of input Z-matrices on all layers - values were calculated during FW-propagation
            li_A_out: List of output A-matrices - values were calculated during FW-propagation
            li_delta: List of delta-matrices - values for outermost ölayer E to layer N+1 should exist 
        
        Returns: 
            ay_delta_N - delta-matrix of layer N (required in subsequent steps)
            ay_D_N     - derivative matrix for the activation function on layer N 
            ay_grad_N  - matrix with gradient elements of the cost fnction with respect to the weights on layer N 
        '''
        
        if b_print:
            print("Inside _bw_prop_Np1_to_N: N = " + str(N) )
        
        # Prepare required quantities - and add bias neuron to ay_Z_in 
        # ****************************
        
        # Weight matrix meddling betwen layer N and N+1 
        ay_W_N = self._li_w[N]
        shape_W_N   = ay_W_N.shape # due to bias node first dim is 1 bigger than Z-matrix 
        if b_print:
            print("shape of W_N = ", shape_W_N )

        # delta-matrix of layer N+1
        ay_delta_Np1 = li_delta[N+1]
        shape_delta_Np1 = ay_delta_Np1.shape 

        # !!! Add intermediate row (for bias) to Z_N !!!
        ay_Z_N = li_Z_in[N]
        shape_Z_N_orig = ay_Z_N.shape
        ay_Z_N = self._add_bias_neuron_to_layer(ay_Z_N, 'row')
        shape_Z_N = ay_Z_N.shape # dimensions should fit now with W- and A-matrix 
        
        # Derivative matrix for the activation function (with extra bias node row)
        #    can only be calculated now as we need the z-values
        ay_D_N = self._calculate_D_N(ay_Z_N)
        shape_D_N = ay_D_N.shape 
        
        ay_A_N = li_A_out[N]
        shape_A_N = ay_A_N.shape
        
        # print shapes 
        if b_print:
            print("shape of W_N = ", shape_W_N)
            print("
shape of delta_(N+1) = ", shape_delta_Np1)
            print("shape of Z_N_orig = ", shape_Z_N_orig)
            print("shape of Z_N = ", shape_Z_N)
            print("shape of D_N = ", shape_D_N)
            print("shape of A_N = ", shape_A_N)
        
        
        # Propagate delta
        # **************
        if li_delta[N+1] is None:
            print("BW-Prop-error:\n No delta-matrix found for layer " + str(N+1) ) 
            sys.exit()
            
        # Check shapes for np.dot()-operation - here for element [0] of both shapes - as we operate with W.T !
        if ( shape_W_N[0] != shape_delta_Np1[0]): 
            print("BW-Prop-error:\n shape of W_N [", shape_W_N, "]) does not fit shape of delta_N+1 [", shape_delta_Np1, "]" )
            sys.exit() 
        
        # intermediate delta 
        # ~~~~~~~~~~~~~~~~~~
        ay_delta_w_N = ay_W_N.T.dot(ay_delta_Np1)
        shape_delta_w_N = ay_delta_w_N.shape
        
        # Check shapes for element wise *-operation !
        if ( shape_delta_w_N != shape_D_N ): 
            print("BW-Prop-error:\n shape of delta_w_N [", shape_delta_w_N, "]) does not fit shape of D_N [", shape_D_N, "]" )
            sys.exit() 
        
        # final delta 
        # ~~~~~~~~~~~
        ay_delta_N = ay_delta_w_N * ay_D_N
        # reduce dimension again 
        ay_delta_N = ay_delta_N[1:, :]
        shape_delta_N = ay_delta_N.shape
        
        # Check dimensions again - ay_delta_N.shape should fit shape_Z_in_orig
        if shape_delta_N != shape_Z_N_orig: 
            print("BW-Prop-error:\n shape of delta_N [", shape_delta_N, "]) does not fit original shape Z_in_N [", shape_Z_N_orig, "]" )
            sys.exit() 
        
        if N > 0:
            shape_W_Nm1 = self._li_w[N-1].shape 
            if shape_delta_N[0] != shape_W_Nm1[0] : 
                print("BW-Prop-error:\n shape of delta_N [", shape_delta_N, "]) does not fit shape of W_Nm1 [", shape_W_Nm1, "]" )
                sysexit() 
        
        
        # Calculate gradient
        # ********************
        #     required for all layers down to 0 
        # check shapes 
        if shape_delta_Np1[1] != shape_A_N[1]:
            print("BW-Prop-error:\n shape of delta_Np1 [", shape_delta_Np1, "]) does not fit shape of A_N [", shape_A_N, "] for matrix multiplication" )
            sys.exit() 
        
        # calculate gradient             
        ay_grad_N = np.dot(ay_delta_Np1, ay_A_N.T)
        
        # regularize gradient (!!!! without adding bias nodes in the L1, L2 sums) 
        ay_grad_N[:, 1:] += (self._li_w[N][:, 1:] * self._lambda2_reg + np.sign(self._li_w[N][:, 1:]) * self._lambda1_reg) 
        
        #
        # Check shape 
        shape_grad_N = ay_grad_N.shape
        if shape_grad_N != shape_W_N:
            print("BW-Prop-error:\n shape of grad_N [", shape_grad_N, "]) does not fit shape of W_N [", shape_W_N, "]" )
            sys.exit() 
        
        # print shapes 
        if b_print:
            print("shape of delta_N = ", shape_delta_N)
            print("shape of grad_N = ", shape_grad_N)
            
            print(ay_grad_N)
        
        return ay_delta_N, ay_D_N, ay_grad_N

 
This function does more or less exactly what we have requested by our theoretical analysis in the last two articles. Note the intermediate handling of bias nodes! Note also that bias nodes are NOT regarded in regularization terms L1 and L2! The function to calculate the derivative of the activation function is:

   
#
    ''' -- Method to calculate the matrix with the derivative values of the output function at outermost layer '''

    def _calculate_D_N(self, ay_Z_N, b_print= False):
        '''
        This method calculates and returns the D-matrix for the outermost layer
        The D matrix contains derivatives of the output function with respect to local input "z_j" at outermost nodes. 
        
        Returns
        ------
        ay_D_E:    Matrix with derivative values of the output function 
                   with respect to local z_j valus at the nodes of the outermost layer E
        Note: This is a 2-dim matrix over layer nodes and training samples of the mini-batch
        '''
        if self._my_out_func == 'sigmoid':
            ay_D_E = self._D_sigmoid(ay_Z = ay_Z_N)
        
        else:
            print("The derivative for output function " + self._my_out_func + " is not known yet!" )
            sys.exit()
        
        return ay_D_E
 

 

The methods to calculate regularization terms for the loss function are:

   
#
    ''' method do calculate the quadratic regularization term for the loss function '''
    def _regularize_by_L2(self, b_print=False): 
        '''
        The L2 regularization term sums up all quadratic weights (without the weight for the bias) 
        over the input and all hidden layers (but not the output layer)
        The weight for the bias is in the first column (index 0) of the weight matrix - 
        as the bias node's output is in the first row of the output vector of the layer 
        '''
        ilayer = range(0, self._n_total_layers-1) # this excludes the last layer 
        L2 = 0.0
        for idx in ilayer:
            L2 += (np.sum( np.square(self._li_w[idx][:, 1:])) ) 
        L2 *= 0.5 * self._lambda2_reg
        if b_print: 
            print("\nL2: total L2 = " + str(L2) )
        return L2 
#
    ''' method do calculate the linear regularization term for the loss function '''
    def _regularize_by_L1(self, b_print=False): 
        '''
        The L1 regularization term sums up all weights (without the weight for the bias) 
        over the input and all hidden layers (but not the output layer
        The weight for the bias is in the first column (index 0) of the weight matrix - 
        as the bias node's output is in the first row of the output vector of the layer 
        '''
        ilayer = range(0, self._n_total_layers-1) # this excludes the last layer 
        L1 = 0.0
        for idx in ilayer:
            L1 += np.sum(np.abs( self._li_w[idx][:, 1:]))
        L1 *= 0.5 * self._lambda1_reg
        if b_print:
            print("\nL1: total L1 = " + str(L1))
        return L1 

 
Addendum, 15.05.2020:
Also the BW-propagation code presented here will later be the target of optimization steps. We shall see that it – despite working correctly – can be criticized regarding efficiency at several points. See again the article series starting with
MLP, Numpy, TF2 – performance issues – Step I – float32, reduction of back propagation.

Conclusion

We have extended our set of methods quite a bit. At the core of the operations we perform matrix operations which are supported by the Openblas library on a Linux system with multiple CPU cores. In the next article

A simple program for an ANN to cover the Mnist dataset – IX – First Tests

we shall test the convergence of our training for the MNIST data set. We shall see that a MLP with two hidden layers with 70 and 30 nodes can give us a convergence
of the averaged relative error down to 0.006 after 1000 epochs on the test data. However, we have to analyze such results for overfitting. Stay tuned …

Links

My PDF on “The math behind EBP”

A simple Python program for an ANN to cover the MNIST dataset – VI – the math behind the “error back-propagation”

I continue with my article series on how to program a training algorithm for a multi-layer perceptron [MLP]. In the course of my last articles

A simple program for an ANN to cover the Mnist dataset – V – coding the loss function
A simple program for an ANN to cover the Mnist dataset – IV – the concept of a cost or loss function
A simple program for an ANN to cover the Mnist dataset – III – forward propagation
A simple program for an ANN to cover the Mnist dataset – II – initial random weight values
A simple program for an ANN to cover the Mnist dataset – I – a starting point

we have already created code for the “Feed Forward Propagation” algorithm [FFPA] and two different cost functions – “Log Loss” and “MSE”. In both cases we took care of a vectorized handling of multiple data records in mini-batches of training data.

Before we turn to the coding of the so called “error back-propagation” [EBP], I found it usefull to clarify the math behind behind this method for ANN/MLP-training. Understanding the basic principles of the gradient descent method for the optimization of MLP-weights is easy. But comprehending

  • why and how gradient descent method leads to the back propagation of error terms
  • and how we cover multiple training data records at the same time

is not – at least not in my opinion. So, I have discussed the required analysis and resulting algorithmic steps in detail in a PDF which you find attached to this article. I used a four layer MLP as an example for which I derived the partial derivatives of the “Log Loss” cost function for weights of the hidden layers in detail. I afterwards generalized the formalism. I hope the contents of the PDF will help beginners in the field of ML to understand what kind of matrix operations gradient descent leads to.

PDF on the math behind Error Back_Propagation

In the next article we shall encode the surprisingly compact algorithm for EBP. In the meantime I wish all readers Merry Christmas …

Addendum 01.01.2020 / 23.02.2020 : Corrected a missing “-” for the cost function and resulting terms in the above PDF.

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 784×100. 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

(784×100) matrix on (784)-vector, (100×50) matrix on (100)-vector, (50×10) 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

(784×100) matrix on an (784, N)-array, (100×50) matrix on an (100, N)-array, (50×10) 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 Python program for an ANN to cover the MNIST dataset – I – a starting point

For beginners both in Python and Machine Learning [ML] the threshold to do some real programming and create your own Artificial Neural Network [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 in the future, I have some ideas along this line …

So I thought, just let me set up a small Python3 and Numpy based program to create a simple kind of ANN – a “Multilayer Perceptron” [MLP] – and train it for the MNIST dataset. I expected that my readers and I myself would earn something on various methods used in ML during our numerical experiments. Well, we shall see ..-

Regarding ANN-theory I take a brutal shortcut and assume that my readers are already acquainted with the following topics:

  1. what simple neural networks with “hidden” layers look like and what its weights are,
  2. what a “cost function” is and how the gradient descent method works,
  3. what logistic regression is and what the cost function J for it looks like,
  4. why we need the back propagation of deviations from known values and why it gives you the required partial derivatives of a typical cost function with respect to an ANN’s “weights”,
  5. what a mini-batch approach to weight optimization is.

I cannot spare you the effort of studying most of these topics in advance. Otherwise I would have to write an introductory book on ML myself. [I would do, but you need to give me a sponsor 🙂 ] But even if you are not fully acquainted with all the named topics: I shall briefly comment on each and every of the points in the forthcoming articles. The real basics are, however, much better and more precisely documented in the literature; e.g. in the books of Geron and Rashka (see the references at the end of this article). I recommend to read one of these books whilst we move on with a sequence of steps to build the basic code for our MLP.

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 images of handwritten numbers with 28×28 pixels [px]. It is a standard data set used in many elementary courses on ML. The challenge for our 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 a fully fletched multi-layer MLP. “Stochastic Gradient Descent”-approaches for a pure binary classificator 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 S. Raschka in his book (see the second to last section for a reference). However, at multiple points I take my 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 and synonyms:

Category: Each input data element is associated with a category to which it belongs. In our case a category
corresponds to a “digit”. The classification algorithm (here: the MLP) may achieve an ability to predict the association of an unknown MNIST like input data sample with its correct category. It should – after some training – detect the (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“. Another term which we I shall frequently use for such an element is a “sample“. The MNIST input set of training data consists of 60000 records or samples – which we provide via an array _X_train. The array is two dimensional as each sample consists of values more multiple properties.

Feature: A sample of the input data set may be equivalent to a mathematical vector, whose elements specify (numerical) values for multiple properties – so called “features” – of a sample. Thus, input samples correspond to points in a multidimensional feature space.

Output data set: A complete set of output data after a so called “propagation” through the ANN for the input data set. “Propagation”means a series of defined mathematical transformations of the original features data of the input sample. The number of samples or records in the output data sets is equal to the number of records in the input data set. The output set will be represented by a Numpy array “_ay_ANN_out“.

A data record or sample 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. E.g., we may be given a “probability” for each category – which allows us to decide with which of the categories we should associate the output element.

A simple MLP network – layers, nodes, weights

An AN network is composed of a series of horizontally and/or vertically arranged layers with nodes. The nodes represent the artificial neurons. A MLP is an ANN which has a rather simple structure: It consists of an input layer, multiple sequential intermediate “hidden” layers, and an output layer. All nodes of a specific layer are connected with all nodes of neighboring (!) layers, only. We speak of a “dense” or fully “connected layer” structure.

The simplifying sketch below displays an ANN with just three sequentially arranged layers – an input layer, a “hidden” middle layer and an output layer.

Note that in general there can be (many) more hidden layers than just one. Note also that modern ANNs (e.g. Convoluted Networks) may have a much more complicated topological structure with hundreds of layers. There may also be
cascaded networks where a specific layer has connections to many more than just the neighbor layers.

Input layer and its number of nodes
To feed input data into the MLP we need an “input layer” with sufficient input nodes. How many? Well, this depends on the number of features your data set represents. In the MNIST case a sample image contains 28×28 pixels. For each pixel with a gray value (integer number between 0 and 256) given. So a typical image represents 28×28 = 768 different “features” – i.e. 786 numbers for “gray”-values between 0 and 255. We need as many input nodes in our MLP to represent the full image information by the input layer.

For other input data the number of features may be different; in addition features my follow a multidimensional order or organization which first must be “flattened” out in one dimension. The number of input nodes must then be adjusted accordingly. The number of input nodes should, therefore, be a parameter or be derived from information on the type of input data. The way of how you map complicated and structured features to input layers and whether you map all data to a one dimensional input vector is a question one should think about carefully. (Most people today treat e.g. a time dimension of input data as just a special form of a feature – I regard this as questionable in some cases, but this is beyond this article series …)

For our MLP we always assume a mapping of features to a flat one dimensional vector like structure.

Output layer and its number of nodes
We shall use our MLP for classification tasks in the beginning. We, therefore, assume that the output of the ANN should allow for the distinction 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 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 of our program, too.

How we indicate the association of a (transformed) sample at the output layer 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 in later articles.

The numbers of “hidden layers” and their nodes
We want the numbers of nodes on “hidden layers” to be parameters for our program. For simple data as MNIST images we do not need big networks, but we want to be able to play around a bit with 1 up to 3 layers. (For an ANN to recognize hand written MNIST digits an input layer “L0” and only one hidden layer “L1” before an output layer “L2” are fully sufficient. Nevertheless in most of our experiments we will actually use 2 hidden layers. There are three reasons: You can approximate any continuous function with two hidden layers (with a special non-linear activation function; see below) and an output layer (with just a linear output function). The other reason is that the full mathematical complexity of “learning” of a MLP appears with two hidden layers (see a later article).

Activation and output functions
The nodes in hidden layers use a so called “activation function” to transform aggregated input from different feeding nodes of the previous layer into one distinct value within a defined interval – e.g. between -1 and 1. Again, we should be prepared to have a program 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 in turn must allow for a distinction of categories. This may lead to a special form – e.g. a kind of probability function. So, the type of the “output function” should also be regarded as variable parameter.

A Python class for our ANN and its interface

I develop my code as a Python module in an Eclipse/PyDev IDE, which itself uses a virtual Python3 environment. I described the setup of such a development environment in detail in another previous article of this blog. In the resulting directory structure of the PyDev project I place a module “myann.py” at the location “…../ml_1/mynotebooks/mycode/myann.py”. This file shall contain the code of class “MyANN” for our ANN.

Modules and libraries to import

We need to import some libraries at the head of our Python program first:

'''
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. The alternative “fetch_openml” does not use caching in some older versions and is also in general 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 so called “sigmoid“-function – which is an important version of an activation function. We shall use it most of the time. “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 us to visualize some information gathered during and after training.

The “__init__”-function of our class MyANN

We encapsulate most of the required 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 feed in parameters and to initialize variables of a class instance.

We shall build up our “__init__()”-function during the next articles step by step. In the beginning we shall only focus on attributes and methods of our class required to import the MNIST data and put them into Numpy arrays and to create the basic network layers.

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 Python “doc”-string. We use a “string” to choose the dataset to train our ANN on. To be able to work on other data sets later on we assume that specific methods for importing a variety of special input data sets are implemented in our class. This requires that the class knows exactly which kinds of data sets it is capable to handle. We provide an list with this information below. The other parameters should be clear from their inline documentation.

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 I omit the usage of “properties” and a full encapsulation of private attributes. For convenience reasons I 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 sample data of the chosen dataset – i.e. the training data, the test data for checking of the reliability of the algorithm after training and the corresponding target data (y_…) for classification – in distinct array variables during code execution. The target data in the MNIST case contain the digit a specific sample image ( of _X_train or _X_test..) 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 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). Numpy provides the various dimensions of multidimensional arrays in form of a tuple.

The target data arrays “_y_…” provides the number of “categories” (MNIST: 10 digits) the ANN must distinguish after training. 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 the variable “_n_nodes_layer_out”.

The “weights” of an ANN must be given in form of matrices: A weight describes a connection between two nodes of different adjacent layers. So we have as many connections as there are node combinations (nodex_(N+1), nodey_N), with “nodex_N” meaning a node on layer L_N and nodey_(N+1) a node on layer L_(N+1).

As the number of layers is not fixed, but can be set by the user, I use a Python 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 set as random numbers. To provide such numbers we have to use randomizer functions. Depending on the kind (floating point numbers, integer numbers) of random numbers to produce we use at least two randomizers (randint, uniform). For the weights we use the uniform-randomizer.

Allowed activation and output function names are listed in Python dictionaries which point to respective methods. This allows for an “indirect addressing” of these functions later on. You may recognize this by the direct reference 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.

    # 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 the 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 directly after loading.

The fastest method for importing the MNIST dataset is 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).

Analysis of the input data and the one-hot-encoding of the target labels

A 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 as trivial is, however, 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). So far, so good …

But: The output data of our ANN for one input element will be provided as an array of values for our 10 different categories – and not as a simple number. To account for this we need to encode the “_y_train”-data, i.e. the target labels, into an usable array form. We use two methods to achieve this:

    # Method to get the number of target 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 weight optimization method we shall use later on during training of our MLP is that we will perform weight adjustment for a whole bunch of training samples in one step. Meaning:

We propagate a whole bunch of test data samples in parallel through the grid to get an array with result data (output array) for all samples.

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 the training and learning phase of an ANN is a compromise between

  • using the full data set (for gradient evaluation) during each training step (“batch approach“)
  • and using just one sample of the training data set for gradient descent and weight corrections (“stochastic approach“).

See chapter 4 of the book of Geron and chapter 2 in the book of Raschka for some thorough 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 the resulting vector and matrix operations on modern CPUs and GPUs.

You really should keep the following point in mind to understand the code for the propagation and optimization algorithms discussed in forthcoming articles:
The so called “cost function” will be determined as some peculiar sum over all elements of a batch and the usual evaluation of partial derivatives during gradient descent will be based on matrix operations involving all input elements of a defined batch!

Mini-batches also will help during training in so far as we look at a bunch of multiple selected samples in parallel to achieve bigger steps of our gradient guided descent into an minimum of the cost hyperplane in the beginning – with the disadvantage of making some jumpy stochastic turns on the cost hyperplane instead of a smoother approach.

I probably lost you now 🙂 . The simpler version is: Keep in mind that we later on will work with batches of training data samples in parallel!

However, the separation interface for our categories in the feature space must in the end be adjusted with respect to all given data points of the training set. This means we must perform the training successively for a whole sequence of mini-batches which together cover all available training samples.

What is the shape of the output array?
A single element of the batch 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
the output array for a batch 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. a possible category
  • idx enumerating the data sample within a batch of training data

We shall construct the output function such that it provides something like “probability” values within the interval [0,1] for each node of the output layer. We define a perfectly working MLP as one which – after training – produces a “1.0” at the correct category node (i.e. the expected digit) and “0.0” at all other output nodes.

One-hot-encoding of labels
Later on we must compare the real results for the training samples with the expected correct values. 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 samples of the batch. E.g.: If we expect the digit 7 for the input array of a sample with 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]. The derivation of such an array from a given category label is called “one-hot-encoding“.

By using Numpy’s zero()-function and Pythons “enumerate()”-function we can achieve such an encoding for all data elements of the training data set. See the method “_encode_all_mnist_labels()”. Thus, the array “_ay_onehot” will have a shape of (10, 60000). From this 2-dim array we can later slice out bunches of consecutive test data for mini-batches. The array “_ay_oneval” is provided for convenience and print purposes, only: it provides the expected digit value in addition.

First tests via a Jupyter notebook

Let us test the import of the input data and the label encoding with a Jupyter notebook. In previous articles I have described already how to use such a notebook. I set up 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 the import of libraries.

By the last line I import my present class code. With the second cell I create an instance of my class; the “__init__()”-function is automatically executed ad calls the other methods defined so far:

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 the labels in the enumerate structure correspond to the “1”s in the “_ay_onehot”-array.

Conclusion

Importing the MNIST dataset into Numpy arrays via Keras is simple – and has a good performance. We have learned a bit about “one-hot-encoding” and prepared an array “_ay_onehot”, which we shall use during ANN training and weight optimization. It will allow us to calculate a difference between the actual output values of the ANN at the nodes of the output layer and a “1.0” value at the node for the expected sample category and “0.0” otherwise.

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

Further articles in this series

A simple Python program for an ANN to cover the MNIST dataset – XIV – cluster detection in feature space
A simple Python program for an ANN to cover the MNIST dataset – XIII – the impact of regularization
A simple Python program for an ANN to cover the MNIST dataset – XII – accuracy evolution, learning
rate, normalization

A simple Python program for an ANN to cover the MNIST dataset – XI – confusion matrix
A simple Python program for an ANN to cover the MNIST dataset – X – mini-batch-shuffling and some more tests
A simple Python program for an ANN to cover the MNIST dataset – IX – First Tests
A simple Python program for an ANN to cover the MNIST dataset – VIII – coding Error Backward Propagation
A simple Python program for an ANN to cover the MNIST dataset – VII – EBP related topics and obstacles
A simple Python program for an ANN to cover the MNIST dataset – VI – the math behind the „error back-propagation“
A simple Python program for an ANN to cover the MNIST dataset – V – coding the loss function
A simple Python program for an ANN to cover the MNIST dataset – IV – the concept of a cost or loss function
A simple Python program for an ANN to cover the MNIST dataset – III – forward propagation
A simple Python program for an ANN to cover the MNIST dataset – II – initial random weight values

Numpy matrix multiplication for layers of simple feed forward ANNs

I move a small step forward on my trials as a beginner in machine learning. This is once again going to be an article from which experts on “Artificial Neural networks” [ANNs] will learn nothing new. But some of my readers may also be newbies both to Python, to numpy and to ANNs. For them this article may be helpful.

I shall have a brief look at some matrix (or array) operations which constitute important parts of the information propagation on ANNs. How would such operations look like in Python with Numpy? One should have a clear understanding of the so called “dot”-operation between multidimensional arrays in this context. Such operations can be performed highly optimized on GPUs – and whole sets of data samples can be handled in form of vectorized code instructions. I will try to explain this in form of matrix terms required to simulate information transport between some layers of a very simple network.

Propagation operations on ANNs

Basic Artificial Neural Networks consist of layers with neurons. Each layer has a defined number of neurons. Each neuron is connected with all neurons of the next layer (receiving neurons). The output “O(k, L)” of a neuron N(k,L) of Layer “L” to a neuron N(m,L+1) is multiplied by some “weight” W[(m, L+1),(k, L),].

The contributions of all neurons on layer “L” to N(m, L+1) are summed up. The result is then transformed by some local “activation” function and afterwards presented by N(m, L+1) as output which shall be transported to the neurons of layer L+2. The succession of such operations is called “propagation“.

At the core of propagation we, therefore, find operations which we can express in the following form :

Input vector at N(m, L+1) : I(m, L+1) = SUM(over_all_k) [ W[(m, L+1), (k, L)] * O(k, L) ]
Output vector at N(m, L+1) : O(m, L+1) = f_act( I(m, L+1) )

To make things simpler: Let us set f_act = 1. Then we get

O(m, L+1) = SUM(over all k) [ W[(m, L+1), (k,N)] * O(k, L) ]

The SUM is basically a matrix operations on vector O(k, L). If output vectors were always given by singular values in “k” rows (!) then we could write

O(L+1) = W(L, L+1) * O(L) )

Batches and numerical optimization of the matrix operations for propagation

In the past I have programmed gas flow simulations in exploding stars in Fortran/C and simulated cascaded production networks for Supply Chains in PHP. There, you work with deep and complicated layer structures, too. Therefore I felt well prepared for writing a small program to set up a simple artificial neural network with some hidden layers.

Something, you will find in literature (e.g. in the book of S. Rashka, “Python machine Learning”, 2016, PACKT Publishing] are remarks on batch based training. Rashka e.g. rightfully claims in his discussion of an example code that the use of mini-batches of input data sets is advantageous (see chapter 12, page 379 in the book mentioned above):

“Mini-batch learning has the advantage over online learning that we can make use of our vectorized implementations to improve computational efficiency.”

I swallowed this remark first without much thinking. As a physicist I had used highly optimized libraries for Linear Algebra operation on supercomputers with vector registers already 35 years ago. Ok, numpy can perform a matrix operation on some vector very fast, especially on GPUs – so what?

But when I tried to follow the basic line of programming of
Rashka for an ANN in Python, I suddenly found that I did not fully understand what the code was doing whilst handling batches of input data sets. The reason was that I had not at all grasped the full implications of the quoted remark. It indicated not just the acceleration of a matrix operation on one input vector; no, instead, it indicated that we could and should perform matrix operations in parallel on all input vectors given by a full series of many different data sets in a mini-batch collection.

So, we then want to perform an optimized matrix operation on another matrix of at least 2 dimensions. When you think about it, you will quickly understand that any kind of operation of one n-dim matrix on another matrix with more and different dimensions must be defined very precisely to avoid misunderstandings and miscalculations.

To make things simple:
Let us take a (3,5)-matrix A and a (10,3)-matrix B. What should a “dot“-like operation A * B mean in this case? Would it be suitable at all, regarding the matrix element arrangements in rows and columns?

I studied the numpy documentation on “numpy.dot(A, B)”. See e.g. here. Ooops, not fully trivial: “If A is an N-D array and B is an M-D array (where M>=2), it is a sum product over the last axis of A and the second-to-last axis of B”. The matrices must obviously be arranged in a suitable way ….

So, I decided to performing some simple matrix experiments to get a clear understanding. I assumed a fictitious ANN of three layers and some direct propagation between the layers. What does propagation mean in such a case, how is it expressed in matrix operation terms in Python and what does batch-handling mean in this context?

Three simple layers and a batch of input data sets

The following graphics shows my simple network:

Our batch has 10 data sets of test data for our ANN. Each data set describes 3 properties; thus the input layer “L1” has 3 nodes. Our “hidden” layer “L2” shall have 5 nodes and our output layer “L3” only 2 nodes.

Weight matrices and propagation operations

Now, let us define some simple matrices with input vectors for transportation and weight-matrices for the layer-connections with Python. All activation functions shall in our simple example just perform a multiplication by 1.

The input to layer “L1” – i.e. a batch of 10 data sets – shall be given by a matrix of 10 lines (dimension 1) and 3 columns (dimension 2)

import numpy as np
import scipy
import time 

A1 = np.ones((10,3))
A1[:,1] *= 2
A1[:,2] *= 4
print("\nMatrix of input vectors to layer L1:\n")
print(A1)


print("\nMatrix of input vectors to layer 1:\n")
print(A1)

Matrix of input vectors to layer 1:

Matrix of input vectors to layer L1:

[[1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]
 [1. 2. 4.]]

We now define the weight matrixW1” for transport between layer L1 and layer L2 as:

W1= np.random.randint(1, 10, 5*3)
W1 = w1.reshape(5,3)
print("\nMatrix of weights W1(L1, L2) :\n")
print(W1)

Matrix of weights W1(L1, L2) :

[[3 8 5]
 [4 9 7]
 [4 1 3]
 [8 8 4]
 [3 1 9]]


We set the first dimension to the number of nodes in the target layer (L2) and the second dimension to
the number of nodes in the lower preceding layer (here: L1). Expressed in the shape notation for an array “(dimension 1, dimension 2, …) “:

Shape of W : (dim 1, dim 2) = (node number of target layer, node number of preceding layer)

If we now tried a numpy “dot”-operation for our matrices

A2 = np.dot(W1, A1)

we would get an error:

ValueError: shapes (5,3) and (10,3) not aligned: 3 (dim 1) != 10 (dim 0)

Well, this is actually consistent with the numpy documentation: The last dimension of “W1” should be consistent with the second to last (=first) dimension of “strong>A1“.

What we need to do here is to transpose the matrix of our input vectors for our multiple data sets:

A1 = np.transpose(A1)
print("\nMatrix of output vectors of layer 1 to layer 2:\n")
print(A1)

Matrix of output vectors of layer 1 to layer 2:

[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
 [4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]]

Then

A2 = np.dot(W1, A1)
print("\n\nA2:\n")
print(A2)

A2:

[[39. 39. 39. 39. 39. 39. 39. 39. 39. 39.]
 [50. 50. 50. 50. 50. 50. 50. 50. 50. 50.]
 [18. 18. 18. 18. 18. 18. 18. 18. 18. 18.]
 [40. 40. 40. 40. 40. 40. 40. 40. 40. 40.]
 [41. 41. 41. 41. 41. 41. 41. 41. 41. 41.]]

Now, we need a weight matrix W2 between layer L1 and layer L2:

W2= np.random.randint(1, 3, 2*5)
W2 = W2.reshape(2,5)
print("\n\nW2:\n")
print(W2)

W2:
[[2 2 1 2 2]
 [2 2 2 2 1]]

We followed the same policy – and used the number of nodes in target layer L2 as the first dimension of the weight array.

Now, do we have to to anything with array A2, if we want to use it as an input for the next matrix operation ?

A3 = np.dot(W2, A2)

No, we do not need to change anything! “A2” has the required form, already. The second to last dimension in “A2” is 5 in our example – as is the last dimension in “W2“!

So, the following will work directly :

A3 = np.dot(W2, A2)
print("\n\nA3:\n")
print(A3)

A3:

[[358. 358. 358. 358. 358. 358. 358. 358. 358. 358.]
 [335. 335. 335. 335. 335. 335. 335. 335. 335. 335.]]

You can try out our complete “propagation” code in a Jupyter notebook cell.

import numpy as np
import scipy
import time 

A1 = np.ones((10,3))
A1[:,1] *= 2
A1[:,2] *= 4
print("\nMatrix of input vectors to layer L1:\n")
print(A1)


W1= np.random.randint(1, 10, 5*3)
W1 = W1.reshape(5,3)
print("\nMatrix of weights W1(L1, L2) :\n")
print(W1)

A1 = np.transpose(A1)
print("\nMatrix of output vectors of layer L1 to layer L2:\n")
print(A1)

A2 = np.dot(W1, A1)
print("\n\nA2:\n")
print(A2)

W2= np.random.randint(1, 3, 2*5)
W2 = W2.reshape(2,5)
print("\n\nW2:\n")
print(W2)

A3 = np.dot(W2, A2)
print("\n\nA3:\n")
print(A3)

Conclusion

The first thing we learned is that matrix operations on simple vectors can be extended to a set of vectors, i.e. a matrix. In the case of ANNs one dimension of such an (numpy) array would cover the number of data sets of a typical input mini-batch. The other dimension would cover the “properties” of the input data.

The second thing we saw is that the input matrix to an ANN often must be transposed to work together with weight matrices, if such matrices follow the policy that the first dimension is given by the number of cells in the target layer and the second dimension by the number of nodes of the
preceding layer.

The weighted input to the first target layer is then given by a matrix “dot”-operation on the transposed input matrix. The outcome is a matrix where the first dimension is defined by the number of nodes of the target layer and the second dimension the number of data sets. This has already the correct form for a further propagation to the next layer. Note that the application of an activation function to each of the matrix elements would not change the required arrangement of matrix data!

After the input layer we, therefore, can just continue to directly apply weight matrices to output matrices by a numpy-“dot”-operation. The weight matrix structure just should reflect our shape policy of

(dim1, dim2) = (node number of target layer, node number of preceding layer)

.

In the next articles we shall use these insights to build a Python class for the setup and training of a simple ANN with one or two hidden layers for the analysis of the famous MNIST dataset.