A simple Python program for an ANN to cover the MNIST dataset – V – coding the loss function

We proceed with encoding a Python class for a Multilayer Perceptron [MLP] to be able to study at least one simple examples of an artificial neural network [ANN] in detail. During the articles

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 came so far that we could apply the “Feed Forward Propagation” algorithm [FFPA] to multiple data records of a mini-batch of training data in parallel. We spoke of a so called vectorized form of the FFPA; we used special Linear Algebra matrix operations of Numpy to achieve the parallel operations. In the last article

A simple program for an ANN to cover the Mnist dataset – IV – the concept of a cost or loss function

I commented on the necessity of a so called “loss function” for the MLP. Although not required for a proper training algorithm we will nevertheless encode a class method to calculate cost values for mini-batches. The behavior of such cost values with training epochs will give us an impression of how good the training algorithm works and whether it actually converges into a minimum of the loss function. As explained in the last article this minimum should correspond to an overall minimum distance of the FFPA results for all training data records from their known correct target values in the result vector space of the MLP.

Before we do the coding for two specific cost or loss functions – namely the “Log Loss“-function and the “MSE“-function, I will briefly point out the difference between standard “*”-operations between multidimensional Numpy arrays and real “dot”-matrix-operations in the sense of Linear Algebra. The latte one follows special rules in multiplying specific elements of both matrices and summing up over the results.

As in all the other articles of this series: This is for beginners; experts will not learn anything new – especially not of the first section.

Element-wise multiplication between multidimensional Numpy arrays in contrast to the “dot”-operation of linear algebra

I would like to point out some aspects of combining two multidimensional Numpy arrays which may be confusing for Python beginners. At least they were for me 🙂 . As a former physicist I automatically expected a “*”-like operation for two multidimensional arrays to perform a matrix operation in the sense of linear algebra. This lead to problems when I tried to understand Python code of others.

Let us assume we have two 2-dimensional arrays A and B. A and B shall be similar in the sense that their shape is identical, i.e. A.shape = B.shape – e.g (784, 60000):
The two matrices each have the same specific number of elements in their different dimensions.

Whenever we operate with multidimensional Numpy arrays with the same same shape we can use the standard operators “+”, “-“, “*”, “/”. These
operators then are applied between corresponding elements of the matrices. I.e., the mathematical operation is applied between elements with the same position along the different dimensional axes in A and B. We speak of an element-wise operation. See the example below.

This means (A * B) is not equivalent to the C = numpy.dot(A, B) operation – which appears in Linear Algebra; e.g. for vector and operator transformations!

The”dot()”-operation implies a special operation: Let us assume that the shape of A[i,j,v] is

A.shape = (p,q,y)

and the shape of B[k,w,m] is

B.shape = (r,z,s)

with

y = z .

Then in the “dot()”-operation all elements of a dimension “v” of A[i,j,v] are multiplied with corresponding elements of the dimension “w” of B[k,w,m] and then the results summed up.

dot(A, B)[i,j,k,m] = sum(A[i,j,:] * B[k,:,m])

The “*” operation in the formula above is to be interpreted as a standard multiplication of array elements.

In the case of A being a 2-dim array and B being a 1-dimensional vector we just get an operation which could – under certain conditions – be interpreted as a typical vector transformation in a 2-dim vector space.

So, when we define two Numpy arrays there may exist two different methods to deal with array-multiplication: If we have two arrays with the same shape, then the “*”-operation means an element-wise multiplication of the elements of both matrices. In the context of ANNs such an operation may be useful – even if real linear algebra matrix operations dominate the required calculations. The first “*”-operation will, however, not work if the array-shapes deviate.

The “numpy.dot(A, B)“-operation instead requires a correspondence of the last dimension of matrix A with the second to last dimension of matrix B. Ooops – I realize I just used the expression “matrix” for a multidimensional Numpy array without much thinking. As said: “matrix” in linear algebra has a connotation of a transformation operator on vectors of a vector space. Is there a difference in Numpy?

Yes, there is, indeed – which may even lead to more confusion: We can apply the function numpy.matrix()

A = numpy.matrix(A),
B = numpy.matrix(B)

then the “*”-operator will get a different meaning – namely that of numpy.dot(A,B):

A * B = numpy.dot(A, B)

So, better read Python code dealing with multidimensional arrays rather carefully ….

To understand this better let us execute the following operations on some simple examples in a Jupyter cell:

A1 = np.ones((5,3))
A1[:,1] *= 2
A1[:,2] *= 4
print("\nMatrix A1:\n")
print(A1)


A2= np.random.randint(1, 10, 5*3)
A2 = A2.reshape(5,3)
# A2 = A2.reshape(3,5)
print("\n Matrix A2 :\n")
print(A2)


A3 = A1 * A2 
print("\n\nA3:\n")
print(A3)

A4 = np.dot(A1, A2.T)
print("\n\nA4:\n")
print(A4)

A5 = np.matrix(A1)
A6 = np.matrix(A2)

A7 = A5 * A6.T 
print("\n\nA7:\n")
print(A7)

A8 = A5 * A6

We get the following output:


Matrix A1:

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

 Matrix A2 :

[[6 8 9]
 [9 1 6]
 [8 8 9]
 [2 8 3]
 [5 8 8]]


A3:

[[ 6. 16. 36.]
 [ 9.  2. 24.]
 [ 8. 16. 36.]
 [ 2. 16. 12.]
 [
 5. 16. 32.]]


A4:

[[58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]]


A7:

[[58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]
 [58. 35. 60. 30. 53.]]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-10-4ea2dbdf6272> in <module>
     28 print(A7)
     29 
---> 30 A8 = A5 * A6
     31 

/projekte/GIT/ai/ml1/lib/python3.6/site-packages/numpy/matrixlib/defmatrix.py in __mul__(self, other)
    218         if isinstance(other, (N.ndarray, list, tuple)) :
    219             # This promotes 1-D vectors to row vectors
--> 220             return N.dot(self, asmatrix(other))
    221         if isscalar(other) or not hasattr(other, '__rmul__') :
    222             return N.dot(self, other)

<__array_function__ internals> in dot(*args, **kwargs)

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

This example obviously demonstrates the difference of an multiplication operation on multidimensional arrays and a real matrix “dot”-operation. Note especially how the “*” operator changed when we calculated A7.

If we instead execute the following code

A1 = np.ones((5,3))
A1[:,1] *= 2
A1[:,2] *= 4
print("\nMatrix A1:\n")
print(A1)


A2= np.random.randint(1, 10, 5*3)
#A2 = A2.reshape(5,3)
A2 = A2.reshape(3,5)
print("\n Matrix A2 :\n")
print(A2)

A3 = A1 * A2 
print("\n\nA3:\n")
print(A3)

we directly get an error:


Matrix A1:

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

 Matrix A2 :

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

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-12-c4d3ffb1e683> in <module>
     13 
     14 
---> 15 A3 = A1 * A2
     16 print("\n\nA3:\n")
     17 print(A3)

ValueError: operands could not be broadcast together with shapes (5,3) (3,5) 

As expected!

Cost calculation for our ANN

As we want to be able to use different types of cost/loss functions we have to introduce new corresponding parameters in the class’s interface. So we update the “__init__()”-function:


    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', 
                 
                 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 ?
            
            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._ay_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
            }  
        
        
          
        # 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    

        # list for cost values of mini-batches during training 
        # The list will later be split into sections for epochs 
        self._ay_cost_vals = []

        # 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()

        # number of epochs 
        self._n_epochs = n_epochs
        # maximum number of batches to handle (<0 => all!) 
        self._n_max_batches = n_max_batches

        # regularization parameters
        self._lambda2_reg = lambda2_reg
        self._lambda1_reg = lambda1_reg

        # paramter 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 mini-batch index array, too 
        self._prepare_epochs_and_batches()
        
        # perform training 
      
  start_c = time.perf_counter()
        self._fit(b_print=False, 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 way of accessing a method/function by a parameterized “name”-string should already be familiar from other methods. The method with the given name must of course exist in the Python module; otherwise already Eclipse#s PyDev we display errors.

    '''-- Method to set the loss function--'''
    def _check_and_set_loss_function(self):
        # check for known loss functions 
        try: 
            if (self._my_loss_func not in self.__d_loss_funcs ): 
                raise ValueError
        except ValueError:
            print("\nThe requested loss function " + self._my_loss_func + " is not known!" )
            sys.exit()   
             
        # set the function to variables for indirect addressing 
        self._loss_func = self.__d_loss_funcs[self._my_loss_func]
        
        if self._b_print_test_data:
            z = 2.0
            print("\nThe loss function of the ANN/MLP was defined as \""  + self._my_loss_func + '"') 
        '''
        '''
        return None
#

The “Log Loss” function

The “LogLoss”-function has a special form. If “a_i” characterizes the FFPA result for a special training record and “y_i” the real known value for this record then we calculate its contribution to the costs as:

Loss = SUM_i [- y_i * log(a_i) – (1 – y_i)*log(1 – a_i)]

This loss function has its justification in statistical considerations – for which we assume that our output function produces a kind of probability distribution. Please see the literature for more information.

Now, due to the encoded result representation over 10 different output dimensions in the MNIST case, corresponding to 10 nodes in the output layer; see the second article of this series, we know that a_i and y_i would be 1-dimensional arrays for each training data record. However, if we vectorize this by treating all records of a mini-batch in parallel we get 2-dim arrays. Actually, we have already calculated the respective arrays in the second to last article.

The rows (1st dim) of a represent the output nodes (training data records, the columns (2nd dim) of a represent the results of the FFPA-result values, which due to our output function have values in the interval ]0.0, 1.0].

The same holds for y – with the difference, that 9 of the values in the rows are 0 and exactly one is 1 for a training record.

The “*” multiplication thus can be done via a normal element-wise array “multiplication” on the given 2-dim arrays of our code.

a = ay_ANN_out
y = ay_y_enc

Numpy offers a function “numpy.sum(M)” for a multidimensional array M, which just sums up all element values. The result is of course a simple scalar.

This information should be enough to understand the following new method:

    ''' 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

The Mean Square Error [MSE] cost function

Although not often used for classification tasks (but more for regression problems) this loss function is so simple that we encode it on the fly. Here we just calculate something like a mean quadratic error:

Loss = 9.5 * SUM_i [ (y_ia_i)**2 ]

This loss function is convex by definition and leads to the following method code:

    ''' 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
#

Regularization terms

Regularization is a means against overfitting during training. The trick is that the cost function is enhanced by terms which include sums of linear or quadratic terms of all weights of all layers. This enforces that the weights themselves get minimized, too, in the search for a minimum of the loss function. The less degrees of freedom there are the less the chance of overfitting …

In the literature (see the book hints in the last article) you find 2 methods for regularization – one with quadratic terms of the weights – the so called “Ridge-Regression” – and one based on a sum of absolute values of the weights – the so called “Lasso regression”. See the books of Geron and Rashka for more information.

Loss = SUM_i [- y_i * log(a_i) – (1 – y_i)*log(1 – a_i)]
+  lambda_2 * SUM_layer [ SUM_nodes [ (w_layer_nodes)**2 ] ]
+  lambda_1 * SUM_layer [ SUM_nodes [ |w_layer_nodes| ] ]

Note that we included already two factors “lambda_2” and “lamda_1” by which the regularization terms are multiplied and added to the cost/loss function in the “__init__”-method.

The two related methods are easy to understand:

    ''' 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._ay_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( self._ay_w[idx][:, 1:])) 
        L1 *= 0.5 * self._lambda1_reg
        if b_print:
            print("\nL1: total L1 = " + str(L1))
        return L1 
#

Why do we not start with index “0” in the weight arrays – self._ay_w[idx][:, 1:]?
The reason is that we do not include the Bias-node in these terms. The weight at the bias nodes of the layers is not varied there during optimization!

Note: Normally we would expect a factor of 1/m, with “m” being the number of records in a mini-batch, for all the terms discussed above. Such a constant factor does not hamper the principal procedure – if we omit it consistently also for for the regularization terms discussed below. It can be taken care of by choosing smaller “lambda”s and a smaller step size during optimization.

Inclusion of the loss function calculations within the handling of mini-batches

For our approach with mini-batches (i.e. an approach between pure stochastic and full batch handling) we have to include the cost calculation in our method “_handle_mini_batch()” to handle mini-batches. Method “_handle_mini_batch()” is modified accordingly:

    ''' -- Method to deal with a batch -- '''
    def _handle_mini_batch(self, num_batch = 0, b_print_y_vals = False, b_print = False):
        '''
        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 ay_Z_in_layer and ay_A_out_layer 
        and fill in the first input elements for layer L0  
        '''
        ay_Z_in_layer  = [] # Input vector in layer L0;  result of a matrix operation in L1,...
        ay_A_out_layer = [] # Result of activation function 
    
        #print("num_batch = " + str(num_batch))
        #print("len of ay_mini_batches = " + str(len(self._ay_mini_batches))) 
        #print("_ay_mini_batches[0] = ")
        #print(self._ay_mini_batches[num_batch])
    
        # Step 1: Special treatment of the ANN's input Layer L0
        # Layer L0: Fill in the input vector for the ANN'
s input layer L0 
        ay_idx_batch = self._ay_mini_batches[num_batch]
        ay_Z_in_layer.append( self._X_train[ay_idx_batch] ) # numpy arrays can be indexed by an array of integers
        #print("\nPropagation : Shape of X_in = ay_Z_in_layer = " + str(ay_Z_in_layer[0].shape))           
        if b_print_y_vals:
            print("\n idx, expected y_value of Layer L0-input :")           
            for idx in self._ay_mini_batches[num_batch]:
                print(str(idx) + ', ' + str(self._y_train[idx]) )
        
        # Step 2: Layer L0: We need to transpose the data of the input layer 
        ay_Z_in_0T       = ay_Z_in_layer[0].T
        ay_Z_in_layer[0] = ay_Z_in_0T

        # Step 3: Call the forward propagation method for the mini-batch data samples 
        self._fw_propagation(ay_Z_in = ay_Z_in_layer, ay_A_out = ay_A_out_layer, b_print = b_print) 
        
        if b_print:
            # index range of layers 
            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(ay_Z_in_layer[il].shape))
                print("Shape of A_out of layer L" + str(il) + " = " + str(ay_A_out_layer[il].shape))

        
        # Step 4: To be done: cost calculation for the batch 
        ay_y_enc = self._ay_onehot[:, ay_idx_batch]
        ay_ANN_out = ay_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)
        self._ay_cost_vals.append(total_costs_batch)
        
        # Step 5: To be done: gradient calculation via back propagation of errors 
        # Step 6: Adjustment of weights  
        
        # try to accelerate garbage handling
        if len(ay_Z_in_layer) > 0:
            del ay_Z_in_layer
        if len(ay_A_out_layer) > 0:
            del ay_A_out_layer
        
        return None
#

 

Note that we save the cost values of every batch in the 1-dim array “self._ay_cost_vals”. This array can later on easily be split into arrays for epochs.

The whole process must be supplemented by a method which does the real cost value calculation:

    ''' -- Main Method to calculate costs -- '''
    def _calculate_loss_for_batch(self, ay_y_enc, ay_ANN_out, b_print = False, b_print_details = False ):
        '''
        Method which calculates the costs including regularization terms  
        The cost function is called according to an input parameter of the class 
        '''
        pure_costs_batch = self._loss_func(ay_y_enc, ay_ANN_out, b_print = False)
        
        if ( b_print and b_print_details ): 
            print("Calc_Costs: Shape of ay_ANN_out = " + str(ay_ANN_out.shape))
            print("Calc_Costs: Shape of ay_y_enc = " + str(ay_y_enc.shape))
        if b_print: 
            print("From Calc_Costs: pure costs of a batch =  " + str(pure_costs_batch))
        
        # Add regularitzation terms - L1: linear reg. term, L2: quadratic reg. term 
        # the sums over the weights (squared) have to be performed for each batch again due to intermediate corrections 
        L1_cost_contrib = 0.0
        L2_cost_contrib = 0.0
        if self._lambda1_reg > 0:
            L1_cost_contrib = self._regularize_by_L1( b_print=False )
        if self._lambda2_reg > 0:
            L2_cost_contrib = self._regularize_by_L2( b_print=False )
        
        total_costs_batch = pure_costs_batch + L1_cost_contrib + L2_cost_contrib
        return total_costs_batch
#

Conclusion

By the steps
discussed above we completed the inclusion of a cost value calculation in our class for every step dealing with a mini-batch during training. All cost values are saved in a Python list for later evaluation. The list can later be split with respect to epochs.

In contrast to the FFP-algorithm all array-operations required in this step were simple element-wise operations and summations over all array-elements.

Cost value calculation obviously is simple and pretty fast regarding CPU-consumption! Just test it yourself!

In the next article we shall analyze the mathematics behind the calculation of the partial derivatives of our cost-function with respect to the many weights at all nodes of the different layers. We shall see that the gradient calculation reduces to remarkable simple formulas describing a kind of back-propagation of the error terms [y_ia_i] through the network.

We will not be surprised that we need to involve some real matrix operations again as in the FFPA !

 

A simple Python program for an ANN to cover the MNIST dataset – IV – the concept of a cost or loss function

We continue with our effort to write a Python class for a Multilayer Perceptron [MLP] – a simple form of an artificial neural network [ANN]. In the previous articles of this series

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

I have already explained

  • what parts of an MLP setup we need to parameterize; e.g. the number of layers, the number of nodes per layer, the activation and output functions;
  • how we create node layers and the corresponding weight arrays,
  • how (and also a bit of why) we work with “mini-batches” of test data during training,
  • how we can realize a “vectorized” form of the required “Feed Forward Propagation” algorithm [FFP]. A vectorized form enables us to process all training data records of a mini-batch in parallel. We used Linear Algebra functions provided by Numpy for this purpose; these functions are supported by the the OpenBlas library on a Linux system.

We also set up a basic loop over a number of epochs during training. (Remember: An epoch corresponds to a training step over all training data records). The number of epochs is handled as a parameter to the class’s interface. By artificially repeating the FFP algorithm up to a thousand times, we already got an impression of the code’s performance and its dependence on the number of CPU cores and the size of a mini-batch.

A special method of our class MyANN controls the handling of a mini-batch of multiple input data records via two major steps so far:

  • Step 1: Extract the data records for the mini-batch from the input data.
  • Step 2: Apply FW-propagation to all data records of the mini-batch.

The next natural step would be to encode a training algorithm which optimizes the weight parameters of our MLP. However, in this article we shall not code anything. Instead, I shall discuss some aspects of the so called “cost function” of a MLP. I think this to be useful to get a basic understanding of what training of an ANN actually means and what the differences are in comparison to other ML-algorithms as e.g. the SVM approach. Understanding the cost function’s role for the training of a MLP will also help to better understand the origin and the mathematical form of the back-propagation-algorithm used for training and discussed in a later article.

I simplify a lot below; more details can be found in the literature on machine Learning [ML]; see the section “Links” for some references. Note that if you know all about the theoretical concepts behind ANN training you will not learn anything new here. This is for beginners (and for later reference in this article series).

The concept of a cost function: Turning a classification problem into an optimization problem

What do we mean by training an ANN? Training means to optimize the weights of the ANN such that the “Feed Forward Propagation” in the end delivers correct predictions for new datasets. A cost function is a central concept of the so called “gradient descent method” used for this optimization. By the way: A synonym for cost function is “loss function”. We use both terms
alike below.

The relation between ANN-training based on a loss function and the classification task, which we want to solve with an ANN, is a subtle one. Let us first discuss what we understand by “classification”:

Classification means to separate the input data into categories; i.e.: finding categorical separation surfaces in the multidimensional vector space of input data. In case of the MNIST dataset such separation interfaces should discriminate between 10 different clusters of data points.

I have discussed the problem of finding a separation surface for the case of the moons dataset example in previous articles in this blog. We then used SVM-algorithms to solve this particular problem. Actually, we determined parameters of (non-linear) polynomials to define a separation surface with a (soft) maximum distance from category related clusters of data points in an extended feature space (=input vector space). The extended feature space covered not only basic features of the input data but also powers of it.

All in all we worked directly in an multidimensional extension of the input vector space and optimized parameters describing linear separation interfaces there. If we had several categories instead of 2 we could use a so called “one versus all”-strategy to calculate 10 linear separation interfaces and determine the distance of any new data point towards the separation surfaces as a confidence measure (score) for a prediction. The separation with the highest score would be used to discriminate between the 10 possible solutions and choose the optimal one. Yes, working in an extended input vector space and with parameters of multiple linear separation surfaces was a bit difficult.

Actually, working with ANNs and cost functions corresponds to a more elegant way of optimizing; it starts with measuring distances in the output vector space of the ANN/MLP:

In the context of classification tasks (with known results for training data) a loss function provides a fictitious cost value which weighs the deviations (or distances) of calculated result values (of the ML-algorithm under training) from the already known correct result for training records. I.e. it measures the errors for the training data records in the output space. The optimization task then means to minimize the cost function and thereby minimize a kind of mean error for all input data records.
The hope is that the collection of resulting weight values allows for predictions of other unknown input data, too.

The result of an ANN/MLP for a training data record is the outcome of a complex transformation performed by the ANN. In case of an MLP the transformation of input into output data is done by the “Feed Forward Propagation” algorithm [FFP]. Thus a reasonably designed cost function becomes dependent on the parameters of the FFP-algorithm – predominantly on the weights given at the nodes of the MLP’s layers. We concentrate on this type of parameter below; but note that in special ANN cases there may be additional other parameters to be varied for training and ANN optimization.

The MLP’s weights can in principle be varied continuously during training. The parameter (vector) space thus can be described by multiple real value axes – one for each of the weights. The parameter space of a MLP is a multi-dimensional one with a dimension equal to or bigger than the space of input data – and of course also the result space. (That the dimension is bigger follows from the required node number in the input layer.)

With the help of a suitable cost function we can pose a mathematical optimization problem for the weight parameters:

Find a point in the weight vector space for which the cost
function gives us a minimum, which in turn corresponds to an overall minimum of the deviation distances.

A simple example for a cost function would be a sum of square values for the length of the difference vectors in the output space for all training data.

There are several things to mention:

  1. The result space is a multidimensional vector space (in case of MNIST a 10 dimensional one); so the distance between points there has to be defined via a mathematical norm over components.
  2. The result space in classification problems typically has a much smaller dimension “m” than the dimension “n” of the space of the input data (m < n).
  3. It makes almost no sense to display the cost function over the multidimensional space of input data – as a working ML-algorithm should deliver small cost values for all input data. However, it makes a lot of sense to display the costs over the multi-dimensional vector space of continuous weight values.
  4. We deal with batches of many training data records; it follows that a reasonable cost function in this case must combine deviations of individual records from optimal values. This is very often done via some kind of sum over individual cost contributions from each training record.

A continuous differentiable cost function defines a hyperplane for gradient-descent

In many MLP cases the cost function will be a function of the weight parameters only; this requires a reasonable node independent form of the activation functions. A loss function with a continuous dependency on all ANN parameters (as the weights) provides a multidimensional hyperplane in an (n+1)-dimensional space – with “n” being the number of FFP variables. The (n+1)-th dimension is for the cost values. As the the FFP-algorithm depends on a multitude of linear and non-linear operations we expect that the hyperplane-surface will have a rather complex form – with maxima and minima as well as so called saddle points.

However, if we construct the cost function cleverly the optimum values for the ANN’s weights will lead to a global minimum of this hyperplane – which then in turn corresponds to a minimum of distances between the propagation results and the known values for the training data:

The task to find categorical separation surfaces in the vector space of input data is reformulated as an optimization task in the cost-weight vector space: There it means finding a (global) minimum of the cost hyperplane.

Let us assume we sit at some point on a yet unexplored hyperplane. A quite general way to find the (global) minimum of this hyperplane is to follow a path indicated by the (tangential) gradient vector at the local point: The gradient is vertically oriented with respect to contour lines of constant cost values on the hyperplane. It thus gives us the direction along which a maximum cost change occurs per unit change of some weights. Calculating corrections of the weights translates into following the gradient with small steps. Geometrically speaking:

We follow the direction the overall gradient points to – and translate the movement into to small components along each weight axis – which gives us the individual weight corrections. Our hope is that the overall gradient points into the direction of the global minimum. (In case of local minima or large planes of the hyperplane we would have to adopt the step size somehow.)

This is called the “gradient descent method“. In one of the next contributions to this article series we shall see how this in turn efficiently translates into the backward propagation of errors through the network via matrix operations. Our optimization task is thus reduced to a
systematic variation of the weights during gradient descent with a series of mathematical operations determining gradient components and resulting weight corrections.

Smooth or stochastic gradient descent?

The cost function absorbs complexity stemming from the large amount of all training data rather smoothly by summing up the individual contributions of training data records. Let us look a bit at the gradient: Normally we would have to calculate partial derivatives of all cost contributions off all data sets with respect to all individual weights. For big training data sets this corresponds to a lot of mathematical operations – both matrix operations (linear algebra) and value calculations of nonlinear (activation and output) functions.

What happens if we took not all data records but concentrated on the contributions of selected input data, only? And corrected afterwards again for another disjunctive set of selected data points? I.e. what if we calculated the full required correction only piece-wise for different collections (mini-batches) of input data records?

Then the reduced gradient components would guide us into a direction on the hyperplane which deviates from the overall gradients direction. Taking the next data record would correct this movement a bit into another direction again. If we perform gradient correction for batches of different data records or in the extreme case for individual records we would move somewhat erratically around the overall gradient’s direction; we speak of a “stochastic gradient descent” [SGD].

The erratic movement of SGD helps to overcome local overall minima. But all in all it may take more steps to come to a global minimum or at least close to it – as the a stochastic movement may never converge into the overall minimum’s point in the weight space – but hop instead around it.

The question of how many input data we include in the cost function determining one single weight correction step during an epoch leads to the choice between the following cases:

  • stochastic gradient descent (sequence of weight corrections during an epoch – each based on just one training data record at a time and for all weights),
  • full batch gradient descent (one weight correction per epoch – based on all training data records and for all weights),
  • mini-batch gradient descent (sequence of weight corrections during an epoch – each based on a batch of multiple training data records and for all weights).

A stochastic or mini-bath based gradient descent may mean much faster corrections in terms of a reduced number of (vectorized) mathematical operations and CPU consumption – at least at the beginning of the descent. The CPU time of the training process for large amounts of input data may actually be reduced by factors!

In the case of mini-batches we can, therefore, optimize the performance by varying the mini-batch size. The required matrix operations can be performed vectorized over all data records of the batch; i.e. the operations can be performed “in parallel”. Fortunately, we do not need to care about the necessary CPU register handling whilst coding – optimized libraries will take care of this. As we have seen already in this blog, also threading for a reasonable amount of CPU cores may influence the performance on a specific system a lot.

For our Python class we will therefore provide parameters for the size of a mini-batch – and adapt both the calculation of cost-contributions and respective weight corrections accordingly.

Note that we do not only hope for that the weights determined by gradient descent provide reasonable result values for the training data but also for any other data later on provided to the ANN/MLP. Solving the optimization problem in the end must provide reliable and complex separation surfaces in
the multidimensional input vector space (for MNIST with a dimension of n=784). The mathematical equivalence of the problem of finding separation surfaces in the input vector space to the optimization problem in the result space can be proven for regression problems. (Actually, I do not know whether a mathematical equivalence has been proven for general problems. So, for some ML classification tasks gradient descent may not work sufficiently well.)

Choosing a cost function

Cost functions should be designed carefully. A “cost function” must have certain properties for the so called “gradient descent method” to work successfully:

  • For convenience the global extremum should be a minimum.
  • The cost function must be continuous and differentiable with respect to the ANN’s weights.
  • The requirement of differentiability translates back to the requirement of differentiable activation and output functions – as we shall see in detail in a later article.
  • It should expose a basic convex form in the surroundings of the global minimum (second partial derivatives > 0).
  • The “cost function” must have certain properties for making use of an efficient way to calculate gradients, i.e. partial derivatives. We shall see that some reasonable cost functions turn this task into a back propagation of errors. The efficiency comes via similar matrix operations as those used in the forward propagation algorithm.

Besides choosing a cost function carefully also the choice of the activation function is important for the success of gradient descent. The path to global minimum on a hyperplane may also depend on the starting point (defined by the statistically chosen initial weight values) as well as on an adaptive step size (called learning rate).

Most Machine Learning algorithms can incorporate a variety of reasonable “cost functions. For classification tasks often the following cost functions are used:

  • Categorial Cross-Entropy
  • Log Loss ( = Logistic Regression Loss )
  • Relative Entropy,
  • Exponential Loss
  • MSE (Mean Square Error)

Each of these functions is more or less appropriate for a specific type of classification problem. See the literature for more information on each of these cost functions.

In our code for MNIST-problem we will only include two of these functions as a starting point – Log Loss and MSE. MSE is e.g. used by T. Rashid in his book (see section Links) on building an MLP with Python for the MNIST case. Information on the Log Loss function are provided by the book of Rashka and the book of Geron; see the references in the section “Links” below.

Do we need cost function values at all?

The training of an ANN – i.e. the optimization of weights – does not require the explicit calculation of cost values. The reason for this is of course that gradient descent first of all works with partial derivatives with respect to weights. To calculate them we must use the chain rule with respect to the activation function, the output of lower layers and so on. But the cost values themselves are nowhere required. As a consequence in all of the book of T. Rashid on “Make your Own Neural network” the calculation of costs is never encoded.

Nevertheless, in the next article of this series we shall discuss the code for cost calculations of mini-batches. The reason for this is that we can use the cost values to study the progress of training and the convergence into a minimum: The change of total “costs” provides a way to control and watch the
success of training through its epochs.

Summary and conclusion

The concept of a cost function is central to MLPs and classification tasks: Classification means to separate the input data into categories. The task to find categorical separation surfaces in the vector space of input data is reformulated as an optimization task. This in turn requires us to find a minimum of the cost/loss hyperplane over the multidimensional space of potential weight-parameters. Calculating corrections of the weights during following a gradient guided path to a minimum in turn efficiently translates into the backward propagation of errors through the network via matrix operations.

Links and Literature

https://www.python-course.eu/matrix_arithmetic.php

Gradient descent and cost functions
towardsdatascience.com understanding-the-mathematics-behind-gradient-descent-dde5dc9be06e
ml-cheatsheet readthedocs – gradient-descent.html
page.mi.fu-berlin.de
neural chapter K7.pdf

Regularization
chunml.github.io tutorial on Regularization/

Books
“Neuronale Netze selbst programmieren”, Tariq Rashid, 2017, O’Reilly Media Inc. + dpunkt.verlag GmbH
“Machine Learning mit SciKit-Learn & TensorFlow, Aurelien Geron, 2018, O’Reilly Media Inc. + dpunkt.verlag GmbH
“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

A simple Python program for an ANN to cover the MNIST dataset – III – forward propagation

I continue with my efforts of writing a small Python class by which I can setup and test a Multilayer Perceptron [MLP] as a simple example for an artificial neural network [ANN]. In the last two articles of this series

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

I defined some code elements, which controlled the layers, their node numbers and built weight matrices. We succeeded in setting random initial values for the weights. This enables us to work on the forward propagation algorithm in this article.

Methods to cover training and mini-batches

As we later on need to define methods which cover “training epochs” and the handling of “mini-batches” comprising a defined number of training records we extend our set of methods already now by

An “epoch” characterizes a full training step comprising

  • propagation, cost and derivative analysis and weight correction of all data records or samples in the set of training data, i.e. a loop over all mini-batches.

Handling of a mini-batch comprises

  • (vectorized) propagation of all training records of a mini-batch,
  • cumulative cost analysis for all training records of a batch,
  • cumulative, averaged gradient evaluation of the cost function by back-propagation of errors and summation over all records of a training batch,
  • weight corrections for nodes in all layers based on averaged gradients over all records of the batch data.

Vectorized propagation means that we propagate all training records of a batch in parallel. This will be handled by Numpy matrix multiplications (see below).

We shall see in a forthcoming post that we can also cover the cumulative gradient calculation over all batch samples by matrix-multiplications where we shift the central multiplication and summation operations to appropriate rows and columns.

However, we do not care for details of training epochs and complete batch-operations at the moment. We use the two methods “_fit()” and “_handle_mini_batch()” in this article only as envelopes to trigger the epoch loop and the matrix operations for propagation of a batch, respectively.

Modified “__init__”-function

We change and extend our “__init_”-function of class MyANN a bit:

    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",   
                 
                 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 
                 
                 vect_mode = 'cols', 
                 
                 figs_x1=12.0, figs_x2=8.0, 
                 legend_loc='upper right',
                 
n                 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 
            
            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  
            
            vect_mode: Are 1-dim data arrays (vctors) 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
            
            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._ay_w = []  
        
        # --- New -----
        # Two lists for output of propagation
        # __ay_x_in  : input data of mini-batches on the different layers; the contents is calculated by the propagation algorithm    
        # __ay_a_out : output data of the activation function; the contents is calculated by the propagation algorithm
        # Note that the elements of these lists are numpy arrays     
        self.__ay_X_in  = []  
        self.__ay_a_out = [] 
        
        
        # 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    

        # 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()

        # number of epochs 
        self._n_epochs = n_epochs
        # maximum number of batches to handle (<0 => all!) 
        self._n_max_batches = n_max_batches


        # print 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 mini-batch index array, too 
        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()

 
Readers who have followed me so far will recognize that I renamed the parameter “n_mini_batch” to “n_size_mini_batch” to indicate its purpose a bit more clearly. We shall derive the number of required mini-batches form the value of this parameter.
I have added two new parameters:

  • n_epochs = 1
  • n_max_batches = -1

“n_epochs” will later receive the user’s setting for the number of epochs to follow during training. “n_max_Batches” allows us to limit the number of mini-batches to analyze during tests.

The kind reader will also have noticed that I encapsulated the series of operations for preparing the weight-matrices for the ANN in a new method “_set_ANN_structure()

    
    '''-- Main 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() 

        # check and set activation functions 
        self._check_and_set_activation_and_out_functions()
        
        return None

 
The called functions have remained unchanged in comparison to the last article.

Preparing epochs and batches

We can safely assume that some steps must be performed to prepare epoch- and batch handling. We, therefore, introduced a new function “_prepare_epochs_and_batches()”. For the time being this method only calculates the number of mini-batches from the input parameter “n_size_mini_batch”.

We use the Numpy-function “array_split()” to split the full range of input data into batches.

 
    ''' -- Main Method to prepare epochs -- '''
    def _prepare_epochs_and_batches(self):
        # set number of mini-batches and array with indices of input data sets belonging to a batch 
        self._set_mini_batches()
        return None
##    
    ''' -- Method to set the number of batches based on given batch size -- '''
    def _set_mini_batches(self, variant=0): 
        # number of mini-batches? 
        self._n_mini_batches = math.ceil( self._y_train.shape[0] / self._n_size_mini_batch )
        print("num of mini_batches = " + str(self._n_mini_batches))
        
        # create list of arrays with indices of batch elements 
        self._ay_mini_batches = np.array_split( range(self._y_train.shape[0]), self._n_mini_batches )
        print("\nnumber of batches : " + str(len(self._ay_mini_batches)))
        print("length of first batch : " + str(len(self._ay_mini_batches[0])))
        print("length of last batch : "  + str(len(self._ay_mini_batches[self._n_mini_batches - 1]) ))
        return None

 
Note that the approach may lead to smaller batch sizes than requested by the user.
array_split() cuts out a series of sub-arrays of indices of the training data. I.e., “_ay_mini_batches” becomes a 1-dim array, whose elements are 1-dim arrays, too. Each of the latter contains a collection of indices for selected samples of the training data – namely the indices for those samples which shall be used in the related mini-batch.

Preliminary elements of the method for training – “_fit()”

For the time being method “_fit()” is used for looping over the number of epochs and the number of batches:

 
    ''' -- Method to set the number of batches based on given batch size -- '''
    def _fit(self, b_print = False, b_measure_batch_time = False):
        # range of epochs
        ay_idx_epochs  = range(0, self._n_epochs)
        
        # limit the number of mini-batches
        n_max_batches = min(self._n_max_
batches, self._n_mini_batches)
        ay_idx_batches = range(0, n_max_batches)
        if (b_print):
            print("\nnumber of epochs = " + str(len(ay_idx_epochs)))
            print("max number of batches = " + str(len(ay_idx_batches)))
        
        # looping over epochs
        for idxe in ay_idx_epochs:
            if (b_print):
                print("\n ---------")
                print("\nStarting epoch " + str(idxe+1))
            
            # loop over mini-batches
            for idxb in ay_idx_batches:
                if (b_print):
                    print("\n ---------")
                    print("\n Dealing 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, b_print_y_vals = False, b_print = b_print)
                if b_measure_batch_time: 
                    end_0 = time.perf_counter()
                    print('Time_CPU for batch ' + str(idxb+1), end_0 - start_0) 
        
        return None
#

 
We limit the number of mini_batches. The double-loop-structure is typical. We tell function “_handle_mini_batch(num_batch = idxb,…)” which batch it should handle.

Preliminary steps for the treatment of a mini-batch

We shall build up the operations for batch handling over several articles. In this article we clarify the operations for feed forward propagation, only. Nevertheless, we have to think a step ahead: Gradient calculation will require that we keep the results of propagation layer-wise somewhere.

As the number of layers can be set by the user of the class we save the propagation results in two Python lists:

  • ay_Z_in_layer = []
  • ay_A_out_layer = []

The Z-values define a collection of input vectors which we normally get by a matrix multiplication from output data of the last layer and a suitable weight-matrix. The “collection” is our mini-batch. So, “ay_Z_in_layer” actually is a 2-dimensional array.

For the ANN’s input layer “L0”, however, we just fill in an excerpt of the “_X”-array-data corresponding to the present mini-batch.

Array “ay_A_out_layer[n]” contains the results of activation function applied onto the elements of “ay_Z_in_layer[n]” of Layer “Ln”. (In addition we shall add a value for a bias neutron; see below).

Our method looks like:

 
    ''' -- Method to deal with a batch -- '''
    def _handle_mini_batch(self, num_batch = 0, b_print_y_vals = False, b_print = False):
        '''
        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 ay_Z_in_layer and ay_A_out_layer 
        and fill in the first input elements for layer L0  
        '''
        ay_Z_in_layer  = [] # Input vector in layer L0;  result of a matrix operation in L1,...
        ay_A_out_layer = [] # Result of activation function 
    
        #print("num_batch = " + str(num_batch))
        #print("len of ay_mini_batches = " + str(len(self._ay_mini_batches))) 
        #print("_ay_mini_batches[0] = ")
        #print(self._ay_mini_batches[num_batch])
    
        # Step 1: Special treatment of the ANN's input Layer L0
        # Layer L0: Fill in the input vector for the ANN's input layer L0 
       
 ay_Z_in_layer.append( self._X_train[(self._ay_mini_batches[num_batch])] ) # numpy arrays can be indexed by an array of integers
        #print("\nPropagation : Shape of X_in = ay_Z_in_layer = " + str(ay_Z_in_layer[0].shape))           
        if b_print_y_vals:
            print("\n idx, expected y_value of Layer L0-input :")           
            for idx in self._ay_mini_batches[num_batch]:
                print(str(idx) + ', ' + str(self._y_train[idx]) )
        
        # Step 2: Layer L0: We need to transpose the data of the input layer 
        ay_Z_in_0T       = ay_Z_in_layer[0].T
        ay_Z_in_layer[0] = ay_Z_in_0T

        # Step 3: Call the forward propagation method for the mini-batch data samples 
        self._fw_propagation(ay_Z_in = ay_Z_in_layer, ay_A_out = ay_A_out_layer, b_print = b_print) 
        
        if b_print:
            # index range of layers 
            ilayer = range(0, self._n_total_layers)
            print("\n ---- ")
            print("\nAfter propagation through all layers: ")
            for il in ilayer:
                print("Shape of Z_in of layer L" + str(il) + " = " + str(ay_Z_in_layer[il].shape))
                print("Shape of A_out of layer L" + str(il) + " = " + str(ay_A_out_layer[il].shape))

        
        # Step 4: To be done: cost calculation for the batch 
        # Step 5: To be done: gradient calculation via back propagation of errors 
        # Step 6: Adjustment of weights  
        
        # try to accelerate garbage handling
        if len(ay_Z_in_layer) > 0:
            del ay_Z_in_layer
        if len(ay_A_out_layer) > 0:
            del ay_A_out_layer
        
        return None

 
Why do we need to transpose the Z-matrix for layer L0?
This has to do with the required matrix multiplication of the forward propagation (see below).

The function “_fw_propagation()” performs the forward propagation of a mini-batch through all of the ANN’s layers – and saves the results in the lists defined above.

Important note:
We transfer our lists (mutable Python objects) to “_fw_propagation()”! This has the effect that the array of the corresponding values is referenced from within “_fw_propagation()”; therefore will any elements added to the lists also be available outside the called function! Therefore we can use the calculated results also in further functions for e.g. gradient calculations which will later be called from within “_handle_mini_batch()”.

Note also that this function leaves room for optimization: It is e.g. unnecessary to prepare ay_Z_in_0T again and again for each epoch. We will transfer the related steps to “_prepare_epochs_and_batches()” later on.

Forward Propagation

In one of my last articles in this blog I already showed how one can use Numpy’s Linear Algebra features to cover propagation calculations required for information transport between two adjacent layers of a feed forward “Artificial Neural Network” [ANN]:
Numpy matrix multiplication for layers of simple feed forward ANNs

The result was that we can cover propagation between neighboring layers by a vectorized multiplication of two 2-dim matrices – one containing the weights and the other vectors of feature data for all mini-batch samples. In the named article I discussed in detail which rows and columns are used for the central multiplication with weights and summations – and that the last dimension of the input array should account for the mini-batch samples. This requires the transpose operation on the input array of Layer L0. All other intermediate layer results (arrays) do already get the right form for vectorizing.

“_fw_propagation()” takes the following form:

 
    ''' -- Method to handle FW propagation for a mini-batch --'''
    def _fw_propagation(self, ay_Z_in, ay_A_out, b_print= False):
        
        b_internal_timing = False
        
        # index range of layers 
        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(ay_Z_in[il].shape))
            
            # Step 1: Take input of last layer and apply activation function 
            if il == 0: 
                A_out_il = ay_Z_in[il] # L0: activation function is identity 
            else: 
                A_out_il = self._act_func( ay_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     
            ay_A_out.append(A_out_il)
            if b_print: 
                print("Shape of A_out of layer L" + str(il) + " (with bias) = " + str(ay_A_out[il].shape))
            
            # Step 3: Propagate by matrix operation 
            Z_in_ilp1 = np.dot(self._ay_w[il], A_out_il) 
            ay_Z_in.append(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(ay_Z_in[il].shape))
        A_out_il = self._out_func( ay_Z_in[il] ) # use the output function 
        ay_A_out.append(A_out_il)
        if b_print:
            print("Shape of A_out of last layer L" + str(il) + " = " + str(ay_A_out[il].shape))
        
        return None
#

 
First we set a range for a loop over the layers. Then we apply the activation function. In “step 2” we add a bias-node to the layer – compare this to the number of weights, which we used during the initialization of the weight matrices in the last article. In step 3 we apply the vectorized Numpy-matrix multiplication (np.dot-operation). Note that this is working for layer L0, too, because we already transposed the input array for this layer in “_handle_mini_batch()”!

Note that we need some special treatment for the last layer: here we call the out-function to get result values. And, of course, we do not add a bias neuron!

It remains to have a look at the function “_add_bias_neuron_to_layer(A_out_il, ‘row’)”, which extends the A-data by a constant value of “1” for a bias neuron. The function is pretty simple:

    ''' Method to add values for a bias neuron to A_out '''
    def _add_bias_neuron_to_layer(self, A, how='column'):
        if how == 'column':
            A_new = np.ones((A.shape[0], A.shape[1]+1))
            A_new[:, 1:] = A
        elif how == 'row':
            A_new = np.ones((A.shape[0]+1, A.shape[1]))
            A_new[1:, :] = A
        return A_new    

A first test

We let the program run in a Jupyter cell with the following parameters:

This produces the following output ( I omitted the output for initialization):

 
Input data for dataset mnist_keras : 
Original shape of X_train = (60000, 28, 28)
Original Shape of y_train = (60000,)
Original shape of X_test = (10000, 28, 28)
Original Shape of y_test = (10000,)

Final input data for dataset mnist_keras : 
Shape of X_train = (60000, 784)
Shape of y_train = (60000,)
Shape of X_test = (10000, 784)
Shape of y_test = (10000,)

We have 60000 data sets for training
Feature dimension is 784 (= 28x28)
The number of labels is 10

Shape of y_train = (60000,)
Shape of ay_onehot = (10, 60000)

Values of the enumerate structure for the first 12 elements : 
(0, 6)
(1, 8)
(2, 4)
(3, 8)
(4, 6)
(5, 5)
(6, 9)
(7, 1)
(8, 3)
(9, 8)
(10, 9)
(11, 0)

Labels for the first 12 datasets:

Shape of ay_onehot = (10, 60000)
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]
 [0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
 [1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 1. 0. 1. 0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 1. 0.]]

The node numbers for the 4 layers are : 
[784 100  50  10]

Shape of weight matrix between layers 0 and 1 (100, 785)
Creating weight matrix for layer 1 to layer 2
Shape of weight matrix between layers 1 and 2 = (50, 101)
Creating weight matrix for layer 2 to layer 3
Shape of weight matrix between layers 2 and 3 = (10, 51)

The activation function of the standard neurons was defined as "sigmoid"
The activation function gives for z=2.0:  0.8807970779778823

The output function of the neurons in the output layer was defined as "sigmoid"
The output function gives for z=2.0:  0.8807970779778823
num of mini_batches = 300

number of batches : 300
length of first batch : 200
length of last batch : 200

number of epochs = 1
max number of batches = 2

 ---------

Starting epoch 1

 ---------

 Dealing with mini-batch 1

Starting propagation between L0 and L1
Shape of Z_in of layer L0 (without bias) = (784, 200)
Shape of A_out of layer L0 (with bias) = (785, 200)

Starting propagation between L1 and L2
Shape of Z_in of layer L1 (without bias) = (100, 200)
Shape of A_out of layer L1 (with bias) = (101, 200)

Starting propagation between L2 and L3
Shape of Z_in of layer L2 (without bias) = (50, 200)
Shape of A_out of layer L2 (with bias) = (51, 200)

Shape of Z_in of layer L3 = (10, 200)
Shape of A_out of last layer L3 = (10, 200)

 ---- 

After propagation through all layers: 
Shape of Z_in of layer L0 = (784, 200)
Shape of A_out of layer L0 = (785, 200)
Shape of Z_in of layer L1 = (100, 200)
Shape of A_out of layer L1 = (101, 200)
Shape of Z_in of layer L2 = (50, 200)
Shape of A_out of layer L2 = (51, 200)
Shape of Z_in of layer L3 = (10, 200)
Shape of A_out of layer L3 = (10, 200)

 ---------

 Dealing with mini-batch 2

Starting propagation between L0 and L1
Shape of Z_in of layer L0 (without bias) = (784, 200)
Shape of A_out of layer L0 (with bias) = (785, 200)

Starting propagation between L1 and L2
Shape of Z_in of layer L1 (without bias) = (100, 200)
Shape of A_out of layer L1 (with bias) = (101, 200)

Starting propagation between L2 and L3
Shape of Z_in of layer L2 (without bias) = (50, 200)
Shape of A_out of layer L2 (with bias) = (51, 200)

Shape of Z_in of layer L3 = (10, 200)
Shape of A_out of last layer L3 = (10, 200)

 ---- 

After propagation through all layers: 
Shape of Z_in of layer L0 = (784, 200)
Shape of A_out of layer L0 = (785, 200)
Shape of Z_in of layer L1 = (100, 200)
Shape of A_out of layer L1 = (101, 200)
Shape of Z_in of layer L2 = (50, 200)
Shape of A_
out of layer L2 = (51, 200)
Shape of Z_in of layer L3 = (10, 200)
Shape of A_out of layer L3 = (10, 200)


 ------
Total training Time_CPU:  0.010270356000546599

Stopping program regularily
stopped

 
We see that the dimensions of the Numpy arrays fit our expectations!

If you raise the number for batches and the number for epochs you will pretty soon realize that writing continuous output to a Jupyter cell costs CPU-time. You will also notice strange things regarding performance, multithreading and the use of the Linalg library OpenBlas on Linux system. I have discussed this extensively in a previous article in this blog:
Linux, OpenBlas and Numpy matrix multiplications – avoid using all processor cores

So, for another tests we set the following environment variable for the shell in which we start our Jupyter notebook:

export OPENBLAS_NUM_THREADS=4

This is appropriate for my Quad-core CPU with hyperthreading. You may choose a different parameter on your system!

We furthermore stop printing in the epoch loop by editing the call to function “_fit()”:

self._fit(b_print=False, b_measure_batch_time=False)

We change our parameter setting to:

Then the last output lines become:

The node numbers for the 4 layers are : 
[784 100  50  10]

Shape of weight matrix between layers 0 and 1 (100, 785)
Creating weight matrix for layer 1 to layer 2
Shape of weight matrix between layers 1 and 2 = (50, 101)
Creating weight matrix for layer 2 to layer 3
Shape of weight matrix between layers 2 and 3 = (10, 51)

The activation function of the standard neurons was defined as "sigmoid"
The activation function gives for z=2.0:  0.8807970779778823

The output function of the neurons in the output layer was defined as "sigmoid"
The output function gives for z=2.0:  0.8807970779778823
num of mini_batches = 150

number of batches : 150
length of first batch : 400
length of last batch : 400


 ------
Total training Time_CPU:  146.44446582399905

Stopping program regularily
stopped

Good !
The time required to repeat this kind of forward propagation for a network with only one hidden layer with 50 neurons and 1000 epochs is around 160 secs. As backward propagation is not much more complex than forward propagation this already indicates that we should be able to train such a most simple MLP with 60000 28×28 images in less than 10 minutes on a standard CPU.

Conclusion

In this article we saw that coding forward propagation is a pretty straight-forward exercise with Numpy! The tricky thing is to understand the way numpy.dot() handles vectorizing of a matrix product and which structure of the matrices is required to get the expected numbers!

In the next article

A simple program for an ANN to cover the Mnist dataset – IV – the concept of a cost or loss function

we shall start working on cost and gradient calculation.