“Thank God for giving you a glimpse of heaven, but do not imagine yourself a bird because you can flap your wings.”
— Alfred de Musset, in “The confession of a child of the century”
In this tutorial, we will try to classify hand-written digits using the tools seen in previous chapters.
Precompiling packages...
Images Being precompiled by another process (pid: 60149, pidfile: /home/vituri/.julia/compiled/v1.12/Images/H8Vxc_fPNkl.ji.pidfile)
2537.5 ms ✓ Images
1 dependency successfully precompiled in 3 seconds. 224 already precompiled.
Precompiling packages...
Plots Being precompiled by another process (pid: 60149, pidfile: /home/vituri/.julia/compiled/v1.12/Plots/ld3vC_fPNkl.ji.pidfile)
1327.8 ms ✓ Plots
1 dependency successfully precompiled in 2 seconds. 176 already precompiled.
Precompiling packages...
1662.0 ms ✓ QuartoNotebookWorkerPlotsExt (serial)
1 dependency successfully precompiled in 2 seconds
14.2 The dataset
MNIST is a dataset consisting of 70.000 hand-written digits. Each digit is a 28x28 grayscale image, that is: a 28x28 matrix of values from 0 to 1. To get this dataset, run
If the console asks you to download some data, just press y.
Notice that we only get the first n_train images so this notebook doesn’t take too much time to run. You can increase n_train to 60000 if you like to live dangerously and have enough RAM memory.
Next, we transpose the digits and save them in a vector
figs = [mnist_digits[:, :, i]' |> Matrix for i ∈ 1:size(mnist_digits)[3]];
The first digit, for example, is the following matrix:
n =10figs_plot = [fig .|> Gray for fig in figs[1:n^2]]mosaicview(figs_plot, nrow = n, rowmajor =true)
14.3 Preparing for war
What topological tools can be useful to distinguish between different digits?
Persistence homology with Vietoris-Rips filtration won’t be of much help: all digits are connected, so the 0-persistence is useless; for the 1-dimensional persistence,
1, 3, 5, 7 do not contain holes;
2 and 4 sometimes contain one hole (depending on the way you write it);
0, 6, 9 contain one hole each;
8 contains two holes.
What if we starting chopping the digits with sublevels of some functions? The excentricity function is able to highlight edges. Doing a sublevel filtration with the excentricity function will permit us to separate digits by the amount of edges they have. So 1 and 3 and 7, for example, will have different persistence diagrams.
14.3.1 From square matrices to points in the plane
In order to calculate the excentricity, we need to convert the “image digits” (28x28 matrices) to points in \(\mathbb{R}^2\) (matrices with 2 columns, one for each dimension, which we will call pointclouds). A simple function can do that:
Getting into details: the excentricity of a metric space \((X, d)\) is a measure of how far a point is from the “center”. It is defined as follows for each \(x \in X\):
\[
e(x) = \sum_{y \in X} \frac{d(x, y)}{N}
\]
where \(N\) is the amount of points of \(X\).
Define a function that takes a digit in \(\mathbb{R}^2\) and return the excentricity as an 28x28 image
We store all the excentricities in the excs vector
excs =excentricity.(figs);
and plot a digit with it’s corresponding excentricity
i =5fig = figs[i]exc = excs[i]heatmap(exc |> rotr90)
Looks good! Time to chop it.
14.3.3 Persistence images
Now we calculate all the persistence diagrams using sublevel filtration. This can take some seconds. Julia is incredibly fast, but does not perform miracles (yet!).
pds =map(excs) do ex m =maximum(ex) ex = m .- exripserer(Cubical(ex), cutoff =0.5)end;
We check the first one
pd = pds[i]pd |> barcode
Compare it with the corresponding heatmap above. There are 3 main edges (and one really small one). It seems ok!
We can see the “step-by-step” creation of these connected components in the following mosaic.
r =range(minimum(exc), maximum(exc), length =25) |> reversefigs_filtration =map(r) do vreplace(x -> x ≤ v ? 0:1, exc) .|> Grayendmosaicview(figs_filtration..., rowmajor =true, nrow =5, npad =20)
Now we create the persistence images of all these barcodes in dimension 0 and 1. We pass the entire collection of barcodes to the PersistenceImage function, and it will ensure that all of them are comparable (ie. are on the same grid).
Top left: the barcode of a digit with respect to sublevels using the excentricity function. Top right: the corresponding persistence diagram. Bottom: 0 and 1 dimensional persistence images. They create a pixelated view of the persistence diagram, using a gaussian blur.
14.4 Fitting a model
In order to use these persistence images in a machine learning model, we first need to vectorize them, ie, transform them into a vector. Machine learning models love vectors! The easist way is to just concatenate the persistence images as follows:
functionconcatenate_pds(imgs_0, pds_0, imgs_1, pds_1) persims = [ [vec(imgs_0(pds_0[i])); vec(imgs_1(pds_1[i])) ] for i in1:length(pds) ] X =reduce(hcat, persims)' XendX =concatenate_pds(imgs_0, pds_0, imgs_1, pds_1)y = mnist_labels .|> string;
We can see that X is a matrix with 10000 rows (the amount of digits) and 128 columns (the persistence images concatenated).
It was also important to convert the mnist_labels to strings, because we want to classify the digits (and not do a regression on them).
We now have a vector for each image. What can we do? We need a model that takes a large vector of numbers and try to predict the digit. Neural networks are excellent in finding non-linear relations on vectors. Let’s try one!
Create the layers
functionnn_model(X) model =Chain(Dense(size(X)[2] =>64) ,Dense(64=>10) )endmodel =nn_model(X)
How to separate “6” and “9”? They are isometric! For some people, “2” and “5” are also isometric (just mirror on the x-axis). Functions that only “see” the metric (like the excentricity) will never be able to separate these digits. In digits, the position of the features is important, so let’s add more slicing filtrations to our arsenal.
To avoid writing all the above code-blocks again, we encapsulate the whole process into a function
functionwhole_process( mnist_digits, mnist_labels, f ; imgs_0 =nothing, imgs_1 =nothing , dim_max =1, sigma =1, size_persistence_image =8 ) figs = [mnist_digits[:, :, i]' |> Matrix for i ∈ 1:size(mnist_digits)[3]] excs =f.(figs); pds =map(excs) do ex m =maximum(ex) ex = m .- exripserer(Cubical(ex), cutoff =0.5, dim_max = dim_max)end; pds_0 = pds .|> first pds_1 = pds .|> last if isnothing(imgs_0) imgs_0 =PersistenceImage(pds_0; sigma = sigma, size = size_persistence_image) end if isnothing(imgs_1) imgs_1 =PersistenceImage(pds_1; sigma = sigma, size = size_persistence_image) end persims = [ [vec(imgs_0(pds_0[i])); vec(imgs_1(pds_1[i])) ] for i ineachindex(pds) ] X =reduce(hcat, persims)' y = mnist_labels .|> string return X, y, pds_0, pds_1, imgs_0, imgs_1end;
We now create the sideways filtrations: from the side and from above.
set_value(x, threshold =0.5, value =0) = x ≥ threshold ? value :0functionfiltration_sideways(fig; axis =1, invert =false) fig2 =copy(fig)if axis ==2 fig2 = fig2'|>Matrixendfor i ∈1:28if invert k =29- i else k = i end fig2[i, :] .=set_value.(fig2[i, :], 0.5, k)end fig2end;
and calculate all 4 persistence diagrams. Warning: this can take a few seconds if you are using 60000 digits!
Let’s explore a bit where the model is making mistakes. Collect all the errors
errors =findall(pred_y .!= y);
and plot the first 3
i = errors[1]println("The model predicted a $(pred_y[i]) but it was a $(y[i])")plot_digit(figs[i])
The model predicted a 2 but it was a 1
i = errors[2]println("The model predicted a $(pred_y[i]) but it was a $(y[i])")plot_digit(figs[i])
The model predicted a 5 but it was a 2
i = errors[3]println("The model predicted a $(pred_y[i]) but it was a $(y[i])")plot_digit(figs[i])
The model predicted a 5 but it was a 2
We can make a mosaic with the first 100 errors
n =10figs_plot = [figs[i] .|> Gray for i in errors[1:n^2]]mosaicview(figs_plot, nrow = n, rowmajor =true)
Many of these digits are really ugly! This makes them hard to classify with our sublevel filtrations. Some other functions could be explored.
14.6 Getting new data
Now we want to see if our model really learned something, or if it just repeated what he saw in the training data. To check data, we need to get new data and calculate the accuracy of the same model on this new data.
Even though we used heavy machinery from topology, at the end our persistence images were vectors that indicated the birth and death of edges. Apart from that, the only machine learning algorithm we used was a simple dense neural network to fit these vectors to the correct labels in a non-linear way. State-of-art machine learning models on the MNIST dataset usually can get more than 99% of accuracy, but they use some complicated neural networks with many layers, and the output prediction are hard to explain. These methods, however, are not excludent of each other: we can use the persistence images (and any other vectorized output from TDA) together with other algorithms.
A curious exercise to the reader is to check if a neural network with two parallel inputs - one for the digits images, followed by convolutional layers - other for the vector of persistence images, followed by dense layers can achieve a better result than the convolutional alone.
