A simple Python program for an ANN to cover the MNIST dataset – VII – EBP related topics and obstacles

I continue with my series about a Python program to build simple MLPs:

A simple Python 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

On our tour we have already learned a lot about multiple aspects of MLP usage. I name forward propagation, matrix operations, loss or cost functions. In the last article of this series
A simple program for an ANN to cover the Mnist dataset – VI – the math behind the „error back-propagation“
I tried to explain some of the math which governs “Error Back Propagation” [EBP]. See the PDF attached to the last article.

EBP is an algorithm which applies the “Gradient Descent” method for the optimization of the weights of a Multilayer Perceptron [MLP]. “Gradient Descent” itself is a method where we step-wise follow short tracks perpendicular to contour lines of a hyperplane in a multidimensional parameter space to hopefully approach a global minimum. A step means a change of of parameter values – in our context of weights. In our case the hyperplane is the surface formed by the cost function over the weights. If we have m weights we get a hyperplane in an (m+1) dimensional space.

To apply gradient descent we have to calculate partial derivatives of the cost function with respect to the weights. We have discussed this in detail in the last article. If you read the PDF you certainly have noted: Most of the time we shall execute matrix operations to provide the components of the weight gradient. Of course, we must guarantee that the matrices’ dimensions fit each other such that the required operations – as an element-wise multiplication and the numpy.dot(X,Y)-operation – become executable.

Unfortunately, there are some challenges regarding this point which we have not covered, yet. One objective of this article is to get prepared for these potential problems before we start coding EBP.

Another point worth discussing is: Is there really just one cost function when we use mini-batches in combination with gradient descent? Regarding the descriptions and the formulas in the PDF of the last article this was and is not fully clear. We only built sums there over cost contributions of all the records in a mini-batch. We did NOT use a loss function which assigned to costs to deviations of the predicted result (after forward propagation) from known values for all training data records.

This triggers the question what our
code in the end really does if and when it works with mini-batches during weight optimization … We start with this point.

In the following I try to keep the writing close to the quantity notations in the PDF. Sorry for a bad display of the δs in HTML.

Gradient descent and mini-batches – one or multiple cost functions?

Regarding the formulas given so far, we obviously handle costs and gradient descent batch-wise. I.e. each mini-batch has its own cost function – with fewer contributions than a cost function for all records would have. Each cost function has (hopefully) a defined position of a global minimum in the weights’ parameter space. Taking this into consideration the whole mini-batch approach is obviously based on some conceptually important assumptions:

  • The basic idea is that the positions of the global minima of all the cost-functions for the different batches do not deviate too much from each other in the basic parameter space.
  • If we additionally defined a cost function for all test data records (over all batches) then this cost function should display a global minimum positioned in between the ones of the batches’ cost functions.
  • This also means that there should be enough records in each batch with a really statistical distribution and no specialties associated with them.
  • Contour lines and gradients on the hyperplanes defined by the loss functions will differ from each other. On average over all mini-batches this should not hinder convergence into a common optimum.

To understand the last point let us assume that we have a batch for MNIST dataset where all records of handwritten digits show a tendency to be shifted to the left border of the basic 28×28 pixel frames. Then this batch would probably give us other weights than other batches.

To get a deeper understanding, let us take only two batches. By chance their cost functions may deviate a bit. In the plots below I have just simulated this by two assumed “cost” functions – each forming a hyperplane in 3 dimensions over only two parameter (=weight) dimensions x and y. You see that the “global” minima of the blue and the red curve deviate a bit in their position.

The next graph shows the sum, i.e. the full “cost function”, in green in comparison to the (vertically shifted and scaled) original functions.

Also here you clearly see the differences in the minimas’ positions. What does this mean for gradient descent?

Firstly, the contour lines on the total cost function would deviate from the ones on the cost function hyperplanes of our 2 batches. So would the directions of the different gradients at the point presently reached in the parameter space during optimization! Working with batches therefore means jumping around on the surface of the total cost function a bit erratically and not precisely along the direction of steepest descent there. By the way: This behavior can be quite helpful to overcome local minima.

Secondly, in our simplified example we would in the end not converge completely, but jump or circle around the minimum of the total cost function. Reason: Each batch forces the weight corrections for x,y into different directions namely those of
its own minimum. So, a weight correction induced by one bath would be countered by corrections imposed by the optimization for the other batch. (Regarding MNIST it would e.g. be interesting to run a batch with handwritten digits of Europeans against a batch with digits written by Americans and see how the weights differ after gradient descent has converged for each batch.)

This makes us understand multiple things:

  • Mini-batches should be built with a statistical distribution of records and their composition should be changed statistically from epoch to epoch.
  • We need a criterion to stop iterating over too many epochs senselessly.
  • We should investigate whether the number and thus the size of mini-batches influences the results of EBP.
  • At the end of an optimization run we could invest in some more iterations not for the batches, but for the full cost function of all training records and see if we can get a little deeper into the minimum of this total cost function.
  • We should analyze our batches – if we keep them up and do not create them statistically anew at the beginning of each epoch – for special data records whose properties are off of the normal – and maybe eliminate those data records.

Repetition: Why Back-propagation of 2 dimensional matrices and not vectors?

The step wise matrix operations of EBP are to be performed according to a scheme with the following structure:

  • On a given layer N apply a layer specific matrix “NW.T” (depending on the weights there) by some operational rule on some matrix “(N+1)δS“, which contains some data already calculated for layer (N+1).
  • Take the results and modify it properly by multiplying it element-wise with some other matrix ND (containing derivative expressions for the activation function) until you get a new NδS.
  • Get partial derivatives of the cost function with respect to the weights on layer (N-1) by a further matrix operation of NδS on a matrix with output values (N-1)A.TS on layer (N-1).
  • Proceed to the next layer in backward direction.

The input into this process is a matrix of error-dependent quantities, which are defined at the output layer. These values are then back-propagated in parallel to the inner layers of our MLP.

Now, why do we propagate data matrices and not just data vectors? Why are we allowed to combine so many different multiplications and summations described in the last article when we deal with partial derivatives with respect to variables deep inside the network?

The answer to the first question is numerical efficiency. We operate on all data records of a mini-batch in parallel; see the PDF. The answer to the second question is 2-fold:

  • We are allowed to perform so many independent operations because of the linear structure of our cost-functions with respect to contributions coming from the records of a mini-batch and the fact that we just apply linear operations between layers during forward propagation. All contributions – however non-linear each may be in itself – are just summed up. And propagation itself between layers is defined to be linear.
  • The only non-linearity occurring – namely in the form of non-linear activation functions – is to be applied just on layers. And there it works only node-wise! We do not
    couple values for nodes on one and the same layer.

In this sense MLPs are very simple by definition – although they may look complex! (By the way and if you wonder why MLPs are nevertheless so powerful: One reason has to do with the “Universal Approximation Theorem”; see the literature hint at the end.)

Consequence of the simplicity: We can deal with δ-values (see the PDF) for both all nodes of a layer and all records of a mini-batch in parallel.

Results derived in the last article would change dramatically if we had rules that coupled the Z- or A-values of different nodes! E.g. if the squared value at node 7 in layer X must always be the sum of squared values at nodes 5 an 6. Believe me: There are real networks in this world where such a type of node coupling occurs – not only in physics.

Note: As we have explained in the PDF, the nodes of a layer define one dimension of the NδS“-matrices,
the number of mini-batch records the other. The latter remains constant. So, during the process the δ-matrices change only one of their 2 dimensions.

Some possible pitfalls to tackle before EBP-coding

Now, my friends, we can happily start coding … Nope, there are actually some minor pitfalls, which we have to explain first.

Special cost-, activation- and output-functions

I refer to the PDF mentioned above and its formulas. The example explained there referred to the “Log Loss” function, which we took as an example cost function. In this case the outδS and the 3δS-terms at the nodes of the outermost layer turned out to be quite simple. See formula (21), (22), (26) and (27) in the PDF.

However, there may be other cost functions for which the derivative with respect to the output vector “a” at the outermost nodes is more complicated.

In addition we may have other output or activation functions than the sigmoid function discussed in the PDF’s example. Further, the output function may differ from the activation function at inner layers. Thus, we find that the partial derivatives of these functions with respect to their variables “z” must be calculated explicitly and as needed for each layer during back propagation; i.e., we have to provide separate and specific functions for the provision of the required derivatives.

At the outermost layer we apply the general formulas (84) to (88) with matrix ED containing derivatives of the output-function Eφ(z) with respect to the input z to find EδS with E marking the outermost layer. Afterwards, however, we apply formula (92) – but this time with D-elements referring to derivatives of the standard activation-function φ used at nodes of inner layers.

The special case of the Log Loss function and other loss functions with critical denominators in their derivative

Formula (21) shows something interesting for the quantity outδS, which is a starting point for backward propagation: a denominator depending on critical factors, which directly involve output “a” at the outer nodes or “a” in a difference term. But in our one-hot-approach “a” may become zero or come close to it – during training by accident or by convergence! This is a dangerous thing; numerically we absolutely want to avoid any division by zero or by small numbers close to the numerical accuracy of a programming language.

What mathematically saves us in the special case of Log Loss are formulas (26) and (27), where due to some “magic” the dangerous denominator is cancelled by a corresponding factor in the numerator when we evaluate EδS.

In the general case, however,
we must investigate what numerical dangers the functional form of the derivative of the loss function may bring with it. In the end there are two things we should do:

  • Build a function to directly calculate EδS and put as much mathematical knowledge about the involved functions and operations into it as possible, before employing an explicit calculation of values of the cost function’s derivative.
  • Check the involved matrices, whose elements may appear in denominators, for elements which are either zero or close to it in the sense of the achievable accuracy.

For our program this means: Whether we calculate the derivative of a cost function to get values for “outδS” will depend on the mathematical nature of the cost function. In case of Log Loss we shall avoid it. In case of MSE we shall perform the numerical operation.

Handling of bias nodes

A further complication of our aspired coding has its origin in the existence of bias nodes on every inner layer of the MLP. A bias node of a layer adds an additional degree of freedom whilst adjusting the layer’s weights; a bias node has no input, it produces only a constant output – but is connected with weights to all normal nodes of the next layer.

Some readers who are not so familiar with “artificial neural networks” may ask: Why do we need bias nodes at all?

Well, think about a simple matrix operation on a 2 dim-vector; it changes its direction and length. But if we want to approximate a function for regression or a separation hyperplanes for classification by a linear operation then we need another element which corresponds to a constant translation part in a linear transformation: z = w1*x1 + w2*x2 + const.. Take a simple function y=w*x + c. The “c” controls where the line crosses the y axis. We need such a parameter if our line should separate clusters of points separably distributed somewhere in the (x,y)-plane; the w is not sufficient to orientate and position the hyperplane in the (x,y)-plane.

This very basically what bias neurons are good for regarding the basically linear operation between two MLP-layers. They add a constant to an otherwise linear transformation.

Do we need a bias node on all layers? Definitely on the input layer. However, on the hidden layers a trained network could by learning evolve weights in such a way that a bias neuron comes about – with almost zero weights on one side. At least in principle; however, we make it easier for the MLP to converge by providing explicit “bias” neurons.

What did we do to account for bias nodes in our Python code so far? We extended the matrices describing the output arrays ay_A_out of the activation function (for input ay_Z_in) on the input and all hidden layers by elements of an additional row. This was done by the method “add_bias_neuron_to_layer()” – see the codes given in article III.

The important point is that our weight matrices already got a corresponding dimension when we built them; i.e. we defined weights for the bias nodes, too. Of course, during optimization we must calculate partial derivatives of the cost function with respect to these weights.

The problem is:

We need to back propagate a delta-matrix Nδ for layer N via
( (NW.T).dot(Nδ) ). But then we can not apply a simple element-wise matrix multiplication with the (N-1)D(z)-matrix at layer N-1. Reason: The dimensions do not fit, if we calculate the elements of D only for the existing Z-Values at layer N-1.

There are two solutions for coding:

  • We can add a row artificially and intermediately to the Z-matrix to calculate the D-matrix, then calculate NδS as
    ( (NW.T).dot(Nδ) ) * (N_1)D
    and eliminate the first artificial row appearing in NδS afterwards.
  • The other option is to reduce the weight-matrix (NW) by a row intermediately and restore it again afterwards.

What we do is a matter of efficiency; in our coding we shall follow the first way and test the difference to the second way afterwards.

Check the matrix dimensions

As all steps to back-propagate and to circumvent the pitfalls require a bit of matrix wizardry we should at least check at every step during EBP backward-propagation that the dimensions of the involved matrices fit each other.

Outlook

Guys, after having explained some of the matrix math in the previous article of this series and the problems we have to tackle whilst programming the EBP-algorithm we are eventually well prepared to add EBP-methods to our Python class for MLP simulation. We are going to to this in the next article:
A simple program for an ANN to cover the Mnist dataset – VIII – coding Error Backward Propagation

Literature

“Machine Learning – An Applied Mathematics Introduction”, Paul Wilmott, 2019, Panda Ohana Publishing

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.

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 !