Deep Dreams of a CNN trained on MNIST data – I – a first approach based on one selected map of a convolutional layer

Posted on 26. October 2020 by eremo

It is fun to play around with Convolutional Neural Networks [CNNs] on the level of an dedicated amateur. One of the reasons is the possibility to visualize the output of elementary building blocks of this class of AI networks. The resulting images help to understand CNN algorithms in an entertaining way – at least in my opinion. The required effort is in addition relatively limited: You must be willing to invest a bit of time into programming, but on a quite modest level of difficulty. And you can often find many basic experiments which are in within the reach of limited PC capabilities.

A special area where the visualization of CNN guided processes is the main objective is the field of “Deep Dreams“. Anyone studying AI methods sooner or later stumbles across the somewhat psychedelic, but none the less spectacular images which Google presented in 2016 as a side branch of their CNN research. Today, you can download DeepDream generators from GitHub.

When I read a bit more about “DeepDream” experiments, I quickly learned that people use quite advanced CNN architectures, like Google’s Inception CNNs, and apply them to high resolution images (see e.g. the Book of F. Chollet on “Deep Learning with Keras and Python” and ai.googleblog.com, 2015, inceptionism-going-deeper-into-neural). Even if you pick up an already trained version of an Inception CNN, you need some decent GPU power to do your own experiments. Another questionable point for an interested amateur is: What does one actually learn from applying “generators”, which others have programmed, and what from just following a “user guide” without understanding what a DeepDream SW actually does? Probably not much, even if you produce stunning images after some time…

So, I asked myself: Can one study basic methods of the DeepDream technology with self programmed tools and a simple dataset? Could one create a “DeepDream” visualization with a rather simply structured CNN trained on MNIST data?
The big advantage of the MNIST data set is that the individual samples are small; and the amount of numerical operations, which a related simple CNN must perform on input images, fits well to the capabilities of PC technology – even if the latter is some years old.

After a first look into DeepDream algorithms, I think: Yes, it should be possible. In a way DeepDream experiments are a natural extension of the visualization of CNN filters and maps which I have already discussed in depth in another article series. Therefore, DeepDream visualizations might even help us to better understand how the internal filter of CNNs work and what “features” are. However, regarding the creation of spectacular images we need to reduce our expectations to a reasonably low level:

A CNN trained on MNIST data works with gray images, low resolution and only simple feature patterns. Therefore, we will never produce such impressive images as published by DeepDream artists or by Google. But, we do have a solid chance to grasp some basic principles and ideas of this side-branch of AI with very simplistic tools.

As always in this blog, I explore a new field step-wise and let you as a reader follow me through the learning process. Throughout most of this new series of articles we will use a CNN created with the help of Keras and filter visualization tools which were developed in another article series of this blog. The CNN has been trained on the MNIST data set already.

In this first post we are going to pick just a single selected feature or response map of a deep CNN layer and let it “dream” upon a down-scaled image of roses. Well, “dream“, as a matter of fact, is a misleading expression; but this is true for the whole DeepDream
business – as we shall see. A CNN does not dream; “DeepDream” creation is more to be seen as an artistic discipline using algorithmic image enhancement.

The input image which we shall feed into our CNN today is shown below:

As our CNN works on a resolution level of 28×28 pixels, only, the “dreaming” will occur in a coarse way, very comparable to hallucinations on the blurred vision level of a short-sighted, myopic man. More precisely: Of a disturbed myopic man who works the whole day with images of digits and lets this poor experience enter and manipulate his dreamy visions of nicer things :-).

Actually, the setup for this article’s experiment was a bit funny: I got the input picture of roses from my wife, who is very much interested in art and likes flowers. I am myopic and in my soul still a theoretical physicist, who is much more attracted by numbers and patterns than by roses – if we disregard the interesting fractal nature of rose blossoms for a second :-).

What do DeepDreams based on single maps of trained MNIST CNNs produce?

To rouse your interest a bit or to disappoint you from the start, I show you a typical result of today’s exercise: “Dreams” or “hallucinations” based on MNIST and a selected single map of a deep convolutional CNN layer produce gray scale images with ghost-like “apparitions”.

When these images appeared on my computer screen, I thought: This is fun, indeed! But my wife just laughed – and said “physicists” with a known undertone and something about “boys and toys” …. I hope this will not stop you from reading further. Later articles will, hopefully, produce more “advanced” hallucinations. But as I said: It all depends on your expectations.

But, lets focus: How did I create the simple “dream” displayed above?

Requirements – a CNN and analysis and visualization tools described in another article series of this blog

I shall use results and methods, which I have already explained in another article series. You need a basic understanding of how a CNN works, what convolutional layers, kernel based filters and cost functions are, how we can build simple CNNs with the help of Keras, … – otherwise you will be lost from the beginning.
A simple CNN for the MNIST datasets – I – CNN basics
We also need a CNN, already trained on the MNIST data. I have shown how to build and train a very simple, yet suitable CNN with the help of Keras and Python; see e.g.:
A simple CNN for the MNIST datasets – II – building the CNN with Keras and a first test
A simple CNN for the MNIST dataset – III – inclusion of a learning-rate scheduler, momentum and a L2-regularizer
In addition we need some
code to create input image patterns which trigger response maps or full layers of a CNN optimally. I called such pixel patterns “OIPs”; others call them “features”. I have offered a Python class in the other article series which offers an optimization loop and other methods to work on OIPs and filter visualization.
A simple CNN for the MNIST dataset – XI – Python code for filter visualization and OIP detection

We shall extend this class by further methods throughout our forthcoming work. To develop and run the codes you should have a working Jupyter environment, a virtual Python environment, an IDE like Eclipse with PyDev for building larger code segments and a working Cuda installation for a NVidia graphics card. My 960GTX proved to be fully sufficient for what we are going to do.

Deep “Dream” – or some funny image manipulation?

As it unfortunately happens so often with AI topics: Also in case of the term “DeepDream” the vocabulary is exaggerated and thoroughly misleading. A simple CNN neither thinks nor “dreams” – it is a software manifestation of the results of an optimization algorithm applied to and trained on selected input data. If applied to new input, it will only detect patterns for which it was optimized before. You could also say:

A CNN is a manifestation of learned prejudices.

CNNs and other types of AI networks filter input according to encoded rules which serve a specific purpose and which reflect the properties of the selected training data set. If you ever used the CNN of my other series on your own hand-written images after a training only on the (US-) MNIST images you will quickly see what I mean. The MNIST dataset reflects an American style of writing digits – a CNN trained on MNIST will fail relatively often when confronted with image samples of digits written by Europeans.

Why do I stress this point at all? Because DeepDreams reveal such kinds of “prejudices” in a visible manner. DeepDream technology extracts and amplifies patterns within images, which fit the trained filters of the involved CNN. F. Chollet correctly describes “DeepDream” as an image manipulation technique which makes use of algorithms for the visualization of CNN filters.

The original algorithmic concept for DeepDreams consists of the following steps:

Extend your algorithm for CNN filter visualization (= OIP creation) from a specific map to the optimization of the response of complete layers. Meaning: Use the total response of all maps of a layer to define contributions to your cost function. Then mix these contributions in a defined weighted way.
Take some image of whatever motive you like and prepare 4 or more down-scaled versions of this image, i.e. versions with different levels of size and resolution below the original size and resolution.
Offer the image with the lowest resolution to the CNN as an input image.
Loop over all prepared image sizes :
- Apply your algorithm for filter visualization of all maps and layers to the input image – but only for a very limited amount of epochs.
- Upscale the resulting output image (OIP-image) to the next level of higher resolution.
- Add details of the original image with the same resolution to the upscaled OIP-image.
- Offer the resulting image as a new input image to your CNN.

Readers who followed me through my last series on “a simple CNN for MNIST” should already raise their eyebrows: What if the CNN expects a certain fixed size of of the input image? Well, a good question. I’ll come back to it in a second. For the time being, let us say that we will concentrate more on resolution than on an
actual image size.

The above steps make it clear that we manipulate an image multiple times. In a way we transform the image slowly to improve a layer’s response and repeat the process with growing resolution. I.e., we apply pattern detection and amplification on more and more details – in the end using all available large and small scale filters of the CNN in a controlled way without fully eliminating the original contents.

What to do about the low resolution of MNIST images and the limited capability of a CNN trained on them?

MNIST images have a very low resolution, real images instead a significantly higher one. With our CNN specialized on MNIST input the OIP-creation algorithm only works on (28×28)-images (and with some warnings, maybe, on smaller ones). What to do about it when we work with input images of a size of e.g. 560×560 pixels?

Well, we just work on the given level of resolution! We have three options:

We can downsize the input image itself or parts of it to the MNIST dimensions – with the help of a bicubic interpolation. Then our OIP-algorithm has the chance to detect OIPs on the coarse scale and to change the downsized image accordingly. Then we can upscale the result again to the original image size – and add details again.
We can split the input image into tiles of size (28×28) and offer these tiles as input to the CNN.
We can combine both of the above options.

Its like what a shortsighted human would do: Work with a blurred impression of the full scale image or look at parts of it from a close distance and then reassemble his/her impressions to larger scales.

A first step – apply only one specific map of a convolutional layer on a down-scaled image version

In this article we have a very limited goal for which we do not have to change our tools, yet:

Preparation:
- We choose a map.
- We downscale the original image to (28×28) by interpolation, upscale the result again by interpolating again (with loss) and calculate the difference to the original full resolution image (all interpolations done in a bicubic way).
Loop (4 times or so):
- We apply the OIP-algorithm on the downscaled input image for a fixed amount of epochs
- We upscale the result by bicubic interpolation to the original size.
- We re-add the difference in details.
- We downscale the result again.

With this approach I try to apply some of the elements of the original algorithm – but just on one scale of coarse resolution. I shall discuss the code for realizing the recipe given above with Python and Jupyter in the next article. For today let us look at some of the ghost like apparitions in the dreams for selected maps of the 3rd convolutional layer; see:
A simple CNN for the MNIST dataset – IX – filter visualization at a convolutional layer

DeepDreams based on selected maps of the 3rd convolutional layer of a CNN trained on MNIST data

With the image sections displayed below I have tried to collect results for different maps which focus on certain areas of the input image (with the exception of the first image section).

The first two images of each row display the detected OIP-patterns on the (28×28) resolution level with pixel values encoded in a (viridis) color-map; the third image in gray scale. The fourth image reveals the dream on the blurry vision level – up-scaled and interpolated to the original image size. You may still detect traces of the original rose blossoms i these images. The last two images of each row display the results
after re-adding details of the original image and an adjustment of the width of the value distribution. The detected and enhanced pattern then turns into a whitey, ghostly shadow.

I have given each section a fancy headline.

I never promised you a rose garden …

“Getting out …”

“Donut …”

“Curls to form a 3 …”

“Two of them …”

“The creepy roots of it all …”

“Look at me …”

“A hidden opening …”

“Soft is something different …”

“Central separation …”

Conclusion: A CNN
detects patterns or parts of patterns it was trained for in any kind of offered input …

You can compare the results to some input patterns (OIPs) which strongly trigger individual maps on the 3rd convolutional layer; you will detect similarities. E.g. four OIP- or feature patterns to which map 56 reacts strongly, look like:

This explains the basic shape of the “apparition” in the first “dream”:

This proves that the filters of a trained CNN actually detect patterns, which proved to be useful for a certain training purpose, in any kind of input which shows some traces of such patterns. A CNN simply does not “know” better: If you only have a hammer to interact with the world, everything becomes a nail to you in the end – this is the level of stupidity on which a CNN algorithm works. And it actually is a fundamental ingredient of DeepDream image manipulation – a transportation of learned patterns or prejudices to an environment outside the original training context.

In the next article
Deep Dreams of a CNN trained on MNIST data – II – some code for pattern carving
I provide the code for creating the above images.

Further articles in this series

Deep Dreams of a CNN trained on MNIST data – II – some code for pattern carving
Deep Dreams of a CNN trained on MNIST data – III – catching dream patterns at smaller length scales

A simple CNN for the MNIST dataset – X – filling some gaps in filter visualization

Posted on 4. October 2020 by eremo

I continue my series of articles on a CNN for the MNIST dataset.

In the last article we harvested the fruits of previous efforts. We produced a variety of input image patterns which triggered certain maps of the innermost convolutional layer of our CNN and the filters behind it maximally. I called such pixel-patterns OIPs. (I am still careful to avoid the expression “feature”, which is used by many authors as a term describing a physical entity identified by the human brain. A connotation which I do not like …)

Although my code for creating OIPs allows for a variation of fluctuations on 4 different length scales it proved to be hard to find OIPs for quite a lot of maps at the highest convolutional layer and for their filters. I could not always find a combination of initial random pixel fluctuations which led to loss values > 0. More precisely: I did not find a pattern by trial and error on a reasonable short timescale of some minutes.

We already know that OIP pixel patterns for the innermost convolutional layer are a bit artificial:

They are idealizations; the displayed pattern may not be present in the same form in the real input images which were used during training.
They may contain repetitions of sub-patterns at different positions.

The images in the last articles showed us in addition that some maps at the innermost Conv-layer are sensitive to rather complicated and specific patterns on relatively large length scales.

This is to be expected as at least filters on the higher convolutional levels accumulate information on a coarse level of image coverage: Filtered information coming from original large areas of the image are in the end convoluted down to grids of (3×3) neurons, i.e. onto low resolution maps. Filters at this convolutional depth can therefore require relatively large scale patterns to be passed. So, it is no real wonder that some maps do not react at all to random statistical fluctuations on very short length scales, e.g. on a two pixel scale. The activation of the respective neurons may stay at zero.

In this article I shall describe a simple method which allowed me to create OIP patterns for around half of the maps (at the 3d Conv-level) for which I did not have any luck before. This is done by a “precursor run” which tests the reaction of a map to a large sample of input images with pattern variations on relatively long length scales.

Systematic analysis? My simple approach …

What is the basic idea? Let us assume that we cover the original image surface by a a grid of squares, e.g. by 16 (7×7)-squares – where (7×7) means a square with a side length of 7 pixels. We then assign a constant grey-value to all the pixels in such a square. Instead of picking random values in the range [0, 255] we use a limited number of N distinct normalized values in an interval [0, 1].

We then investigate all combinations of distinct value distributions across the 16 squares. For each combination we construct an input image with bicubic interpolation on the original (28×28) scale and standardize the resulting distribution of pixel values. (The standardization is required for my CNN, which was trained on standardized images). Thus we produce a (huge) number of input images with smooth large scale fluctuations, which we then can present to our OIP-optimization algorithm for a test of a map’s reaction.

Do we have a chance to find OIPs by systematic trials?

Now, do we have a real chance to work a bit more systematically with the described approach? Unfortunately, the answer is: Not really, at least not with a PC equipment and for fluctuations at middle length scales with more than 5 distinct pixel values. The reason in terms of Fourier series’ is that the number of amplitude combinations for multiple sinus waves grows exponentially with a defined number of discrete amplitudes values out of an interval.

We see the limitations already in our simplified approach. Let us assume that we pick 5 distinct and normalized pixel values in the interval ]0, 1[. Let us further assume that we cover the whole image area (of a size 28×28) with a grid of (9×9) squares, each with a constant pixel value. (We ignore the fact that 28/9 = 3.11 and trust in bicubic interpolation to take care of the 0.11 🙂 )

Even under these conditions we get already around 5**9 = 1.953 million pattern variations for which we have to test a map’s response. For each of these patterns we would have to follow at least 4 to 10 iterations (= epochs) of the optimization loop to identify the input image as a promising candidate for a full optimization or not. This means that we need 8 to 20 million full gradient calculations on a parameter space with 784 (=28×28) dimensions. Hard work even for graphic cards (at the consumer level). At smaller length scales, as induced e.g. by a (7×7) coverage grid, we go far beyond standard calculation capacities on PC hardware.

What was/is within my computational reach? If we reduced the number of distinct pixel values to 3 (e.g.: 0.2, 0.5, 0.8) and concentrated on large scale fluctuations – e.g. based on (9×9) squares – we would have around 19700 combinations to investigate.

I had to work a bit on the respective Python code to get it run under Tensorflow 2.2.1/Cuda 10.2 with a reasonable percentage of around 30% on a graphics card. With 10 epochs for an optimization trial run for each image a complete run over 20000 images takes less than 5 minutes on my old GTX960 graphics card. With a modern card, a recent CPU and some more optimization one would probably arrive at values significantly below a minute. So, systematically working with large scale fluctuations is possible on PC-equipment!

Note that a coverage of the image surface with (9×9)-squares corresponds to the resolution which the (3×3)-maps of the innermost Conv-layers represent! (Including padding in the filter definitions). I assume that it is reasonable to work and test on this scale.

Selection of 8 promising candidates out of 19683

So, I added some “precursor” methods to my class My_OIP to prepare and perform systematic runs for all of the maps for which I had not found a pattern, yet. I chose the following 3 discrete pixel values: [0.2, 0.5, 0.8].

The resulting (9×9) images where rescaled with a bicubic interpolation to the MNIST size of (28×28) and afterwards standardized.

During the test of a map’s response I selected 8 candidates out of 19683 for subsequent thorough optimization runs with 1200 epochs. The selection was done by looking at the largest loss values reached after 10 optimization epochs. I then applied the procedure to all 67 maps for which I had not got an OIP, yet.

Note that we cannot be sure whether 10 optimization iterations really give us a perfect indication of the highest loss values which would be reached after a full optimization run with 600 to 1200 epochs. So, it is worthwhile to play around with the number of test epochs; the selection of the input fluctuation patterns may change.

I got results then for 33 of the maps which had no OIP. I present the OIP-images below. The results cover almost 50% of the missing OIPs. So, my computational efforts were well invested.

OIPs for additional 33 maps

Below I present MNIST-related OIP-images for the maps with the following indices on the 3rd Conv-layer of my (!) CNN ):
1, 3, 9, 15, 16, 29, 35, 37, 38, 44, 46, 50, 51, 55, 59, 60, 70, 78, 81, 91, 93, 94, 96, 97, 101, 104, 108, 109, 111, 112, 116, 122, 125.
This time without contrast enhancement – as I was too lazy to integrate it into the precursor code.

CNNs are fun, aren’t they?

Unique OIPs?

An interesting question is: Are the OIPs really unique per map in the sense that we produce the same OIP-images for different fluctuations on the input images? Well, there is some unique overall form of the patterns, but there may occur translations in position and there may be differences in details. Below, I show you different derived images for some maps:

Map 46:

Map 81:

Map 116:

Map 125:

Map 127:

Conclusion

By a systematic investigation of the reaction of maps to large scale fluctuations we can find OIPs for maps where a trial and error approach is not sufficient to trigger a response of the filters. I will discuss the Python code for a related “precursor run” in my next post

A simple CNN for the MNIST dataset – XI – Python code for filter visualization and OIP detection.

A simple Python program for an ANN to cover the MNIST dataset – XIV – cluster detection in feature space

Posted on 14. April 2020 by eremo

We extend our studies of a program for a Multilayer perceptron and gradient descent in combination with the MNIST dataset:

In this article we shall work a bit on the following topic: How can we reduce the computational time required for gradient descent runs of our MLP?

Readers who followed my last articles will have noticed that I sometimes used 1800 epochs in a gradient descent run. The computational time including

costly intermediate print outs into Jupyter cells,
a full determination of the reached accuracy both on the full training and the test dataset at every epoch

lay in a region of 40 to 45 minutes for our MLP with two hidden layers and roughly 58000 weights. Using an Intel I7 standard CPU with OpenBlas
support. And I plan to work with bigger MLPs – not on MNIST but other data sets. Believe me: Everything beyond 10 minutes is a burden. So, I have a natural interest in accelerating things on a very basic level already before turning to GPUs or arrays of them.

Factors for CPU-time

This introductory question leads to another one: What basic factors beyond technical capabilities of our Linux system and badly written parts of my Python code influence the consumption of computational time? Four points come to my mind; you probably find even more:

One factor is certainly the extra forward propagation run which we apply to all samples of both the test and training data seat the end of each epoch. We perform this propagation to make predictions and to get data on the evolution of the accuracy, the total loss and the ratio of the regularization term to the real costs. We could do this in the future at every 2nd or 5th epoch to save some time. But this will reduce CPU-time only by less than 22%. 76% of the CPU-time of an epoch is spent in batch-handling with a dominant part in error backward propagation and weight corrections.
The learning rate has a direct impact on the number of required epochs. We could enlarge the learning rate in combination with input data normalization; see the last article. This could reduce the number of required epochs significantly. Depending on the parameter choices before by up to 40% or 50%. But it requires a bit of experimenting ….
Two other, more important factors are the frequent number of matrix operations during error back-propagation and the size of the involved matrices. These operations depend directly on the number of nodes involved. We could therefore reduce the number of nodes of our MLP to a minimum compatible with the required accuracy and precision. This leads directly to the next point.
The dominant weight matrix is of course the one which couples layer L0 and layer L1. In our case its shape is 784 x 70; it has almost 55000 elements. The matrix for the next pair of layers has only 70×30 = 2100 elements – it is much, much smaller. To reduce CPU time for forward propagation we should try to make this matrix smaller. During error back propagation we must perform multiple matrix multiplications; the matrix dimensions depend on the number of samples in a mini-batch AND on the number of nodes in the involved layers. The dimensions of the the result matrix correspond to the those of the weight matrix. So once again: A reduction of the nodes in the first 2 layers would be extremely helpful for the expensive backward propagation. See: The math behind EBP.

We shall mainly concentrate on the last point in this article.

Reduction of the dimensions of the dominant matrix”requires a reduction of input features

The following numbers show typical CPU times spend for matrix operations during error back propagation [EBP] between different layers of our MLP and for two different batches at the beginning of gradient descent:

Time_CPU for BW layer operations (to L2) 0.00029015699965384556
Time_CPU for BW layer operations (to L1) 0.0008645610000712622
Time_CPU for BW layer operations (to L0) 0.006551215999934357

Time_CPU for BW layer operations (to L2) 0.00029157400012991275
Time_CPU for BW layer operations (to L1) 0.0009575330000188842
Time_CPU for BW layer operations (to L0) 0.007488838999961445

The operations involving layer L0 cost a factor of 7 more CPU time than the other operations! Therefore, a key to the reduction of the number of mathematical operations is obviously the reduction of the number of nodes in the input layer! We cannot reduce the numbers in the hidden layers much, if we
do not want to hamper the accuracy properties of our MLP too much. So the basic question is

Can we reduce the number of input nodes somehow?

Yes, maybe we can! Input nodes correspond to “features“. In case of the MNIST dataset the relevant features are given by the gray-values for the 784 pixels of each image. A first idea is that there are many pixels within each MNIST image which are probably not used at all for classification – especially pixels at the outer image borders. So, it would be helpful to chop them off or to ignore them by some appropriate method. In addition, special significant pixel areas may exist to which the MLP, i.e. its weight optimization, reacts during training. For example: The digits 3, 5, 6, 8, 9 all have a bow within the lower 30% of an image, but in other regions, e.g. to the left and the right, they are rather different.

If we could identify suitable image areas in which dark pixels have a higher probability for certain digits then, maybe, we could use this information to discriminate the represented digits? But a “higher density of dark pixels in an image area” is nothing else than a description of a “cluster” of (dark) pixels in certain image areas. Can we use pixel clusters at numerous areas of an image to learn about the represented digits? Is the combination of (averaged) feature values in certain clusters of pixels representative for a handwritten digit in the MNIST dataset?

If the number of such pixel clusters could be reduced below lets say 100 then we could indeed reduce the number of input features significantly!

Cluster detection

To be able to use relevant “clusters” of pixels – if they exist in a usable form in MNIST images at all – we must first identify them. Cluster identification and discrimination is a major discipline of Machine Learning. This discipline works in general with unlabeled data. In the MNIST case we would not use the labels in the “y”-data at all to identify clusters; we would only use the “X”-data. A nice introduction to the mechanisms of cluster identification is given in the book of Paul Wilcott (see Machine Learning – book recommendations for the reference). The most fundamental method – called “kmeans” – iterates over 3 major steps [I simplify a bit :-)]:

We assume that K clusters exist and start with random initial positions of their centers (called “centroids”) in the multidimensional feature space
We measure the distance of all data points to he centroids and associate a point with that centroid to which the distance is smallest
We determine the “center of mass” (according to some distance metric) of the identified data point groups and assume it as a new position of the centroids and move the old positions (a bit) in this direction.

We iterate over these steps until the centroids’ positions hopefully get stable. Pretty simple. But there is a major drawback: You must make an assumption on the number “K” of clusters. To make such an assumption can become difficult in the complex case of a feature space with hundreds of dimensions.

You can compensate this by executing multiple cluster runs and comparing the results. By what? Regarding the closure or separation of clusters in terms of an appropriate norm. One such norm is called “cluster inertia“; it measures the mean squared distance to the center for all points of a cluster. The theory is that the sum of the inertias for all clusters drops significantly with the number of clusters until an optimal number is reached and the inertia curve flattens out. The point where this happens in a plot of inertia vs. number of clusters is called “elbow“.
Identifying this “elbow” is one of the means to find an optimal number of clusters. However, this recipe does not work under all circumstances. As the number of clusters get big we may be confronted with a smooth decline of the inertia sum.

What data do we use for gradient descent after cluster detection?

How could we measure whether an image shows certain clusters? We could e.g. measure distances (with some appropriate metric) of all image points to the clusters. The “fit_transform()”-method of KMeans and MiniBatchKMeans provide us with with some distance measure of each image to the identified clusters. This means our images are transformed into a new feature space – namely into a “cluster-distance space”. This is a quite complex space, too. But it has less dimensions than the original feature space!

Note: We would of course normalize the resulting distance data in the new feature space before applying gradient descent.

Application of “KMeansBatch” to MNIST

There are multiple variants of “KMeans”. We shall use one which is provided by SciKit-Learn and which is optimized for large datasets: “MiniBatchKMeans“. It operates batch-wise without loosing too much of accuracy and convergence properties in comparison to KMeans (or a comparison see here). “MiniBatchKMeans”has some parameters you can play with.

We could be tempted to use 10 clusters as there are 10 digits to discriminate between. But remember: A digit can be written in very many ways. So, it is much more probable that we need a significant larger number of clusters. But again: How to determine on which K-values we should invest a bit more time? “Kmeans” and methods alike offer another quantity called “silhouette” coefficient. It measures how well the data points are within, at or outside the borders of a cluster. See the book of Geron referenced at the link given above on more information.

Variation of CPU time, inertia and average silhouette coefficients with the number of clusters “K”

Let us first have a look at the evolution of CPU time, total inertia and averaged silhouette with the number of clusters “K” for two different runs. The following code for a Jupyter cell gives us the data:

    
# *********************************************************
# Pre-Clustering => Searching for the elbow 
# *********************************************************
from sklearn.cluster import KMeans
from sklearn.cluster import MiniBatchKMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score
X = np.concatenate((ANN._X_train, ANN._X_test), axis=0)
y = np.concatenate((ANN._y_train, ANN._y_test), axis=0)
print("X-shape = ", X.shape, "y-shape = ", y.shape)
num = X.shape[0]

li_n = []
li_inertia = []
li_CPU = []
li_sil1 = []

# Loop over the number "n" of assumed clusters 
rg_n = range(10,171,10)
for n in rg_n:
    print("\nNumber of clusters: ", n)
    start = time.perf_counter()
    kmeans = MiniBatchKMeans(n_clusters=n, n_init=500, max_iter=1000, batch_size=500 )  
    X_clustered = kmeans.fit_transform(X)
    sil1 = silhouette_score(X, kmeans.labels_)
    #sil2 = silhouette_score(X_clustered, kmeans.labels_)
    end = time.perf_counter()
    dtime = end - start
    print('Inertia = ', kmeans.inertia_)
    print('Time_CPU = ', dtime)
    print('sil1 score = ', sil1)
    li_n.append(n)    
    li_inertia.append(kmeans.inertia_)    
    li_CPU.append(dtime)    
    li_sil1.append(sil1)    

    
# Plots         
# ******
fig_size = plt.rcParams["figure.figsize"]
fig_size[
0] = 14
fig_size[1] = 5
fig1 = plt.figure(1)
fig2 = plt.figure(2)

ax1_1 = fig1.add_subplot(121)
ax1_2 = fig1.add_subplot(122)

ax1_1.plot(li_n, li_CPU)
ax1_1.set_xlabel("num clusters K")
ax1_1.set_ylabel("CPU time")

ax1_2.plot(li_n, li_inertia)
ax1_2.set_xlabel("num clusters K")
ax1_2.set_ylabel("inertia")

ax2_1 = fig2.add_subplot(121)
ax2_2 = fig2.add_subplot(122)

ax2_1.plot(li_n, li_sil1)
ax2_1.set_xlabel("num clusters K")
ax2_1.set_ylabel("silhoutte 1")

You see that I allowed for large numbers of initial centroid positions and iterations to be on the safe side. Before you try it yourself: Such runs for a broad variation of K-values are relatively costly. The CPU time rises from around 32 seconds for 30 clusters to a little less than 1 minute for 180 clusters. These times add up to a significant sum after a while …

Here are some plots:

The second run was executed with a higher resolution of K_(n+1) – K_n 5 = 5.

We see that the CPU time to determine the centroids’ positions varies fairly linear with “K”. And even for 170 clusters it does not take more than a minute! So, CPU-time for cluster identification is not a major limitation.

Unfortunately, we do not see a clear elbow in the inertia curve! What you regard as a reasonable choice for the number K depends a lot on where you say the curve starts to flatten. You could say that this happens around K = 60 to 90. But the results for the silhouette-quantity indicate for our parameter setting that K=40, K=70, K=90 are interesting points. We shall look at these points a bit closer with higher resolution later on.

Reduction of the regularization factor (for Ridge regularization)

Now, I want to discuss an important point which I did not find in the literature:
In my last article we saw that regularization plays a significant but also delicate role in reaching top accuracy values for the test dataset. We saw that Lambda2 = 0.2 was a good choice for a normalized input of the MNIST data. It corresponded to a certain ratio of the regularization term to average batch costs.
But when we reduce the number of input nodes we also reduce the number of total weights. So the weight values themselves will automatically become bigger if we want to get to similar good values at the second layer. But as the regularization term depends in a quadratic way on the weights we may assume that we roughly need a linear reduction of Lambda2. So, for K=100 clusters we may shrink Lambda2 to (0.2/784*100) = 0.025 instead of 0.2. In general:

Lambda2_cluster = Lambda2_std * K / (number of input nodes)

I applied this rule of a thumb successfully throughout experiments with clustering befor gradient descent.

Reference run without clustering

We saw at the end of article XII that we could reach an accuracy of around 0.975 after 500 epochs under optimal circumstances. But in the case I presented ten I was extremely lucky with the statistical initial weight distribution and the batch composition. In other runs with the same parameter setup I got smaller accuracy values. So, let us take an ad hoc run with the following parameters and results:
Parameters: learn_rate = 0.001, decrease_rate = 0.00001, mom_rate = 0.00005, n_size_mini_batch = 500, n_epochs = 600, Lambda2 = 0.2, weights at
all layers in [-2*1.0/sqrt(num_nodes_layer), 2*1.0/sqrt(num_nodes_layer)]
Results: acc_train: 0.9949 , acc_test: 0.9735, convergence after ca. 550-600 epochs

The next plot shows (from left to right and the down) the evolution of the costs per batch, the averaged error of the last mini-batch during an epoch, the ratio of regularization to batch costs and the total costs of the training set, respectively .

The following plot summarizes the evolution of the total costs of the traaining set (including the regularization contribution) and the evolution of the accuracy on the training and the test data sets (in orange and blue, respectively).

The required computational time for the 600 epochs was roughly 18,2 minutes.

Results of gradient descent based on a prior cluster identification

Before we go into a more detailed discussion of code adaption and test runs with things like clusters in unnormalized and normalized feature spaces, I want to show what we – without too much effort – can get out of using cluster detection ahead of gradient descent. The next plot shows the evolution of a run for K=70 clusters in combination with a special normalization:

and the total cost and accuracy evolution

The dotted line marks an accuracy of 97.8%! This is 0.5% bigger then our reference value of 97.3%. The total gain of %gt; 0.5% means however 18.5% of the remaining difference of 2.7% to 100% and we past a value of 97.8% already at epoch 600 of the run.

What were the required computational times?

If we just wanted 97.4% as accuracy we need around 150 epochs. And a total CPU time of 1.3 minutes to get to the same accuracy as our reference run. This is a factor of roughly 14 in required CPU time. For a stable 97.73% after epoch 350 we were still a factor of 5.6 better. For a stable accuracy beyond 97.8% we needed around 600 epochs – and still were by a factor of 3.3 faster than our reference run! So, clustering really brings some big advantages with it.

Conclusion

In
this article I discussed the idea of introducing cluster identification in the (unnormalized or normalized) feature space ahead of gradient descent as a possible means to save computational time. A preliminary trial run showed that we indeed can become significantly faster by at least a factor of 3 up to 5 and even more. This is just due to the point that we reduced the number of input nodes and thus the number of mathematical calculations during matrix operations.

In the next article we shall have a more detailed look at clustering techniques in combination with normalization.

A simple Python program for an ANN to cover the MNIST dataset – I – a starting point

Posted on 29. September 2019 by eremo

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

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

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

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

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

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

Note that this task does NOT require the use of a fully fletched multi-layer MLP. “Stochastic Gradient Descent”-approaches for a pure binary classificator to determine (linear) separation surfaces in combination with a “One-versus-All” strategy for multi-category-classification may be sufficient. See chapters 3 to 5 in the book of Geron for more information.

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

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

Wording

To avoid confusion, I use the following wording and synonyms:

Category: Each input data element is associated with a category to which it belongs. In our case a category
corresponds to a “digit”. The classification algorithm (here: the MLP) may achieve an ability to predict the association of an unknown MNIST like input data sample with its correct category. It should – after some training – detect the (non-linear) separation interfaces for categories in a multidimensional feature space. In the case of MNIST we speak about ten categories corresponding to 10 digits, including zero.

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

Input data set: A complete “set” of input data. Such a set consists of individual “elements” or “records“. Another term which we I shall frequently use for such an element is a “sample“. The MNIST input set of training data consists of 60000 records or samples – which we provide via an array _X_train. The array is two dimensional as each sample consists of values more multiple properties.

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

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

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

A data record of the output data set: One distinct element of the output data set. Note that such an element itself may be an array covering all possible categories in a distinct form. E.g., we may be given a “probability” for each category – which allows us to decide with which of the categories we should associate the output element.

A simple MLP network – layers, nodes, weights

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

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

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

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

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

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

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

How we indicate the association of a (transformed) sample at the output layer to a category numerically – by a probability number between “0” and “1” or just a “1” at the right category and zeros otherwise – can be a matter of discussion. It is also a question of the cost function we wish to use. We will come back to this point in later articles.

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

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

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

A Python class for our ANN and its interface

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

Modules and libraries to import

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

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

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

We need “scipy” to get an optimized version of the so called “sigmoid“-function – which is an important version of an activation function. We shall use it most of the time. “numpy” and “math” are required for fast array- and math-operations. “time” is required to measure the run time of program segments and “mathplotlib” will help us to visualize some information gathered during and after training.

The “init”-function of our class MyANN

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

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

Parameters

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

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

Initialization of class attributes

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

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

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

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

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

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

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

List of known input datasets

The list of known input data sets is kept in the variable “self.__input_data_sets”. The variables

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

will be used to keep all sample data of the chosen dataset – i.e. the training data, the test data for checking of the reliability of the algorithm after training and the corresponding target data (y_…) for classification – in distinct array variables during code execution. The target data in the MNIST case contain the digit a specific sample image ( of _X_train or _X_test..) represents.

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

The input arrays “X_…” will via their dimensions provide the information on the number of data sets (_dim_sets) and the number of features (_dim_features). Numpy provides the various dimensions of multidimensional arrays in form of a tuple.

The target data arrays “_y_…” provides the number of “categories” (MNIST: 10 digits) the ANN must distinguish after training. We keep this number in the
attribute “_n_labels”. It is also useful to keep the pixel dimensions of input image data. At least for MNIST we assume quadratic images (_img_h = img_w = _dim_img).

Layers and weights

The number of total layers (“_n_total_layers”) is by 2 bigger than the number of hidden layers (_n_hidden_layers).

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

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

As the number of layers is not fixed, but can be set by the user, I use a Python list “_ay_w” to collect such matrices in the order of layer_0 (input) to layer_n (output).

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

Allowed activation and output function names are listed in Python dictionaries which point to respective methods. This allows for an “indirect addressing” of these functions later on. You may recognize this by the direct reference of the dictionary elements to defined class methods (no strings are used!).

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

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

Read and provide the input data

Now let us turn to some methods. We first need to read in and prepare the input data. We use a method “_handle_input_data()” to work on this problem. For the time being we have only three different ways to load the MNIST dataset from different origins.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

The first method only determines the number of labels (= number of categories). We see from the code of the second method that we encode the target labels in the form of two arrays. The relevant one for our optimization algorithm will be “_ay_onehot”. This array is 2-dimensional. Why?

Working with mini-batches

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

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

Such a bunch is called a “batch” and if it is significantly smaller than the whole set of training data – a “mini-batch“.

Working with “mini-batches” during the training and learning phase of an ANN is a compromise between

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

See chapter 4 of the book of Geron and chapter 2 in the book of Raschka for some thorough information on this topic.

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

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

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

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

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

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

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

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

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

One-hot-encoding of labels
Later on we must compare the real results for the training samples with the expected correct values. To be able to do this we must build up a 2-dim array of the same shape as “_ay_a_out” with correct output values for all test samples of the batch. E.g.: If we expect the digit 7 for the input array of a sample with index idx within the set of training data, we need a 2-dim output array with the element [[0,0,0,0,0,0,0,1,0,0], idx]. The derivation of such an array from a given category label is called “one-hot-encoding“.

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

First tests via a Jupyter notebook

Let us test the import of the input data and the label encoding with a Jupyter notebook. In previous articles I have described already how to use such a notebook. I set up a Jupyter notebook called “myANN” (in my present working directory “/projekte/GIT/ai/ml1/mynotebooks”).

I start it with

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

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

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

The output is:

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

Conclusion

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

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

Literature and links

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

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

Further articles in this series

What do DeepDreams based on single maps of trained MNIST CNNs produce?

Requirements – a CNN and analysis and visualization tools described in another article series of this blog

Deep “Dream” – or some funny image manipulation?

What to do about the low resolution of MNIST images and the limited capability of a CNN trained on them?

A first step – apply only one specific map of a convolutional layer on a down-scaled image version

DeepDreams based on selected maps of the 3rd convolutional layer of a CNN trained on MNIST data

Conclusion: A CNN detects patterns or parts of patterns it was trained for in any kind of offered input …

Further articles in this series

Systematic analysis? My simple approach …

Do we have a chance to find OIPs by systematic trials?

Selection of 8 promising candidates out of 19683

OIPs for additional 33 maps

Unique OIPs?

Conclusion

Factors for CPU-time

Reduction of the dimensions of the dominant matrix”requires a reduction of input features

Cluster detection

What data do we use for gradient descent after cluster detection?

Application of “KMeansBatch” to MNIST

Variation of CPU time, inertia and average silhouette coefficients with the number of clusters “K”

Reduction of the regularization factor (for Ridge regularization)

Reference run without clustering

Results of gradient descent based on a prior cluster identification

Conclusion

Wording

A simple MLP network – layers, nodes, weights

A Python class for our ANN and its interface

Modules and libraries to import

The “__init__”-function of our class MyANN

Initialization of class attributes

List of known input datasets

Layers and weights

Read and provide the input data

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

Working with mini-batches

First tests via a Jupyter notebook

Conclusion

Literature and links

Further articles in this series

Conclusion: A CNN
detects patterns or parts of patterns it was trained for in any kind of offered input …

The “init”-function of our class MyANN