lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added material

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2021-07-28 23:27:13 +02:00
parent 82dae187cc
commit 76009c48f8
4 changed files with 1591 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/src/Clustering/Clustering.do.txt
+974
View File
@@ -0,0 +1,974 @@
======= Clustering Analysis =======
In this chapter we will concern ourselves with the study of _cluster analysis_.
In general terms cluster analysis, or clustering, is the task of grouping a
data-set into different distinct categories based on some measure of equality of
the data. This measure is often referred to as a _metric_ or _similarity
measure_ in the literature (note: sometimes we deal with a _dissimilarity
measure_ instead). Usually, these metrics are formulated as some kind of
distance function between points in a high-dimensional space.
There exists a lot of such distance measures. The simplest, and also the most
common is the _Euclidean distance_ (i.e. Pythagoras). A good source for those of
you wanting a thorough overview is the article (DOI:10.5120/ijca2016907841
Irani, Pise, Phatak). A few other metrics mentioned there are: *cosine
similarity*, *Manhattan distance*, *Chebychev distance* and the *Minkowski
distance*. The Minkowski distance is a general formulation which encapsulates a
range of metrics. All of these, and many more, can be used in clustering. There
exists different categories of clustering algorithms. A few of the most
common are: *centroid-*, *distribution-*, *density-* and *hierarchical-
clustering*. We will concern ourselves primarily with the first one.
===== Basic Idea of the K-means Clustering Algorithm =====
The simplest of all clustering algorithms is the aptly named _k-means algorithm_
, sometimes also referred to as *Lloyds algorithm*. It is the simplest and also
the most common. From its simplicity it obtains both strengths and weaknesses.
These will be discussed in more detail later. The k-means algorithm is a
_centroid based_ clustering algorithm.
Assume, we are given $n$ data points and we wish to split the data into $K < n$
different categories, or clusters. We label each cluster by an integer $k\in\{
1, \cdots, K \}$. In the basic k-means algorithm each point is assigned to only
one cluster $k$, and these assignments are *non-injective* i.e. many-to-one. We
can think of these mappings as an encoder $k = C(i)$, which assigns the $i$-th
data-point $\bf x_i$ to the $k$-th cluster. Before we jump into the mathematics
let us describe the k-means algorithm in words:
o We start with guesses / random initializations of our $k$ cluster centers /
centroids
o For each centroid the points that are most similar are identified
o Then we move / replace each centroid with a coordinate average of all the
points that were assigned to that centroid.
o Iterate this points 2, 3) until the centroids no longer move (to some
tolerance)
Now we consider the method formally. Again, we assume we have $n$ data-points
(vectors)
!bt
\begin{equation}\label{eq:kmeanspoints}
\bm{x_i} = \{x_{i, 1}, \cdots, x_{i, p}\}\in\mathbb{R}^p.
\end{equation}
!et
which we wish to group into $K < n$ clusters. For our dissimilarity measure we
will use the *squared Euclidean distance*
!bt
\begin{equation}\label{eq:squaredeuclidean}
d(\bm{x_i}, \bm{x_i'}) = \sum_{j=1}^p(x_{ij} - x_{i'j})^2
= ||\bm{x_i} - \bm{x_{i'}}||^2
\end{equation}
!et
Next we define the so called *within-cluster point scatter* which gives us a
measure of how close each data point assigned to the same cluster tends to be to
the all the others.
!bt
\begin{equation}\label{eq:withincluster}
W(C) = \frac{1}{2}\sum_{k=1}^K\sum_{C(i)=k}
\sum_{C(i')=k}d(\bm{x_i}, \bm{x_{i'}}) =
\sum_{k=1}^KN_k\sum_{C(i)=k}||\bm{x_i} - \bm{\overline{x_k}}||^2
\end{equation}
!et
where $\bm{\overline{x_k}}$ is the mean vector associated with the $k$-th
cluster, and $N_k = \sum_{i=1}^nI(C(i) = k)$, where the $I()$ notation is
similar to the Kronecker delta (*Commonly used in statistics, it just means that
when $i = k$ we have the encoder $C(i)$*). In other words, the within-cluster
scatter measures the compactness of each cluster with respect to the data points
assigned to each cluster. This is the quantity that the $k$-means algorithm aims
to minimize. We refer to this quantity $W(C)$ as the within cluster scatter
because of its relation to the *total scatter*.
!bt
\begin{equation}\label{eq:totalscatter}
T = W(C) + B(C) = \frac{1}{2}\sum_{i=1}^n
\sum_{i'=1}^nd(\bm{x_i}, \bm{x_{i'}})
= \frac{1}{2}\sum_{k=1}^K\sum_{C(i)=k}
\Big(\sum_{C(i') = k}d(\bm{x_i}, \bm{x_{i'}})
+ \sum_{C(i')\neq k}d(\bm{x_i}, \bm{x_{i'}})\Big)
\end{equation}
!et
Which is a quantity that is conserved throughout the $k$-means algorithm. It can
be thought of as the total amount of information in the data, and it is composed
of the aforementioned within-cluster scatter and the *between-cluster scatter*
$B(C)$. In methods such as principle component analysis the total scatter is not
conserved.
Given a cluster mean $\bm{m_k}$ we define the _total cluster variance_
!bt
\begin{equation}\label{eq:totalclustervariance}
\min_{C, \{\bm{m_k}\}_1^K}\sum_{k=1}^KN_k\sum||\bm{x_i} - \bm{m_k}||^2
\end{equation}
!et
Now we have all the pieces necessary to formally revisit the k-means algorithm.
If you at this point feel like some of the above definitions came a bit out of
no-where, don't fret, the method does get a whole lot simpler once we start
programming.
===== The K-means Clustering Algorithm =====
The k-means clustering algorithm goes as follows (note in my opinion this
description is a bit complicated and is lifted directly out of ESL HASTIE for
deeper understanding purposes)
o For a given cluster assignment $C$, and $k$ cluster means
$\{m_1, \cdots, m_k\}$. We minimize the total cluster variance with respect to
the cluster means $\{m_k\}$ yielding the means of the currently assigned
clusters.
o Given a current set of $k$ means $\{m_k\}$ the total cluster variance is
minimized by assigning each observation to the closest (current) cluster mean.
That is $$C(i) = \underset{1\leq k\leq K}{\mathrm{argmin}}
||\bm{x_i} - \bm{m_k}||^2$$
o Steps 1 and 2 are repeated until the assignments do not change.
As previously stated the above formulation can be a bit difficult to understand,
*at least the first time*, due to the dense notation used. But all in all the
concept is fairly simple when explained in words. The math needs to be
understood but to help you along the way we summarize the algorithm as follows
(try to look at the terms above to match with the summary).
o Before we start we specify a number $k$ which is the number of clusters we
want to try to separate our data into.
o We initially choose $k$ random data points in our data as our initial
centroids, *or means* (this is where the name comes from).
o Assign each data point to their closest centroid, based on the squared
Euclidean distance.
o For each of the $k$ cluster we update the centroid by calculating new mean
values for all the data points in the cluster.
o Iteratively minimize the within cluster scatter by performing steps (3, 4)
until the new assignments stop changing (can be to some tolerance) or
until a maximum number of iterations have passed.
That's it, nothing magical happening.
===== Writing Our Own Code =====
In the following section we will work to develop a deeper understanding of the
previously discussed mathematics through developing codes to do k-means cluster
analysis.
=== Basic Python ===
Let us now program the most basic version of the algorithm using nothing but
Python with numpy arrays. This code is kept intentionally simple to gradually
progress our understanding. There is no vectorization of any kind, and even most
helper functions are not utilized. Throughout our implementation process it will
be helpful to keep in mind both the mathematical description of the algorithm
*and* our summary from above. In addition, try to think of ways to optimize this
while reading the next section. We will get to it, take it as a challenge to see
if your optimizations are better.
First of all we need a dataset to do our cluster analysis on, for clarity (and
lack of googling beforehand) we generate it ourselves using Gaussians. First we
import
!bc pycod
import time
import numpy as np
import tensorflow as tf
from matplotlib import image
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from IPython.display import display
np.random.seed(2021)
!ec
Next we define functions, for ease of use later, to generate Gaussians and to
set up our toy data set.
!bc pycod
def gaussian_points(dim=2, n_points=1000, mean_vector=np.array([0, 0]),
sample_variance=1):
"""
Very simple custom function to generate gaussian distributed point clusters
with variable dimension, number of points, means in each direction
(must match dim) and sample variance.
Inputs:
dim (int)
n_points (int)
mean_vector (np.array) (where index 0 is x, index 1 is y etc.)
sample_variance (float)
Returns:
data (np.array): with dimensions (dim x n_points)
"""
mean_matrix = np.zeros(dim) + mean_vector
covariance_matrix = np.eye(dim) * sample_variance
data = np.random.multivariate_normal(mean_matrix, covariance_matrix,
n_points)
return data
def generate_simple_clustering_dataset(dim=2, n_points=1000, plotting=True,
return_data=True):
"""
Toy model to illustrate k-means clustering
"""
data1 = gaussian_points(mean_vector=np.array([5, 5]))
data2 = gaussian_points()
data3 = gaussian_points(mean_vector=np.array([1, 4.5]))
data4 = gaussian_points(mean_vector=np.array([5, 1]))
data = np.concatenate((data1, data2, data3, data4), axis=0)
if plotting:
fig, ax = plt.subplots()
ax.scatter(data[:, 0], data[:, 1], alpha=0.2)
ax.set_title('Toy Model Dataset')
plt.show()
if return_data:
return data
data = generate_simple_clustering_dataset()
!ec
Now that we are our, albeit very simple, dataset we are ready to start
implementing the k-means algorithm.
!bc pycod
n_samples, dimensions = data.shape
n_clusters = 4
# we randomly initialize our centroids
np.random.seed(2021)
centroids = data[np.random.choice(n_samples, n_clusters, replace=False), :]
distances = np.zeros((n_samples, n_clusters))
# first we need to calculate the distance to each centroid from our data
for k in range(n_clusters):
for n in range(n_samples):
dist = 0
for d in range(dimensions):
dist += np.abs(data[n, d] - centroids[k, d])**2
distances[n, k] = dist
# we initialize an array to keep track of to which cluster each point belongs
# the way we set it up here the index tracks which point and the value which
# cluster the point belongs to
cluster_labels = np.zeros(n_samples, dtype='int')
# next we loop through our samples and for every point assign it to the cluster
# to which it has the smallest distance to
for n in range(n_samples):
# tracking variables (all of this is basically just an argmin)
smallest = 1e10
smallest_row_index = 1e10
for k in range(n_clusters):
if distances[n, k] < smallest:
smallest = distances[n, k]
smallest_row_index = k
cluster_labels[n] = smallest_row_index
!ec
Let's plot and see
!bc pycod
fig = plt.figure()
ax = fig.add_subplot()
unique_cluster_labels = np.unique(cluster_labels)
for i in unique_cluster_labels:
ax.scatter(data[cluster_labels == i, 0],
data[cluster_labels == i, 1],
label = i,
alpha = 0.2)
ax.scatter(centroids[:, 0], centroids[:, 1], c='black')
ax.set_title("First Grouping of Points to Centroids")
plt.show()
!ec
So what do we have so far? We have 'picked' $k$ centroids at random from our
data points. There are other ways of more intelligently choosing their
initializations, however for our purposes randomly is fine. Then we have
initialized an array 'distances' which holds the information of the distance,
*or dissimilarity*, of every point to of our centroids. Finally, we have
initialized an array 'cluster_labels' which according to our distances array
holds the information of to which centroid every point is assigned. This was the
first pass of our algorithm. Essentially, all we need to do now is repeat the
distance and assignment steps above until we have reached a desired convergence
or a maximum amount of iterations.
!bc pycod
max_iterations = 100
tolerance = 1e-8
for iteration in range(max_iterations):
prev_centroids = centroids.copy()
for k in range(n_clusters):
# this array will be used to update our centroid positions
vector_mean = np.zeros(dimensions)
mean_divisor = 0
for n in range(n_samples):
if cluster_labels[n] == k:
vector_mean += data[n, :]
mean_divisor += 1
# update according to the k means
centroids[k, :] = vector_mean / mean_divisor
# we find the dissimilarity
for k in range(n_clusters):
for n in range(n_samples):
dist = 0
for d in range(dimensions):
dist += np.abs(data[n, d] - centroids[k, d])**2
distances[n, k] = dist
# assign each point
for n in range(n_samples):
smallest = 1e10
smallest_row_index = 1e10
for k in range(n_clusters):
if distances[n, k] < smallest:
smallest = distances[n, k]
smallest_row_index = k
cluster_labels[n] = smallest_row_index
# convergence criteria
centroid_difference = np.sum(np.abs(centroids - prev_centroids))
if centroid_difference < tolerance:
print(f'Converged at iteration {iteration}')
break
elif iteration == max_iterations:
print(f'Did not converge in {max_iterations} iterations')
!ec
And thats it! We now have an extremely barebones, un-optimized k-means
clustering implementation. Lets plot the final result
!bc pycod
fig = plt.figure()
ax = fig.add_subplot()
unique_cluster_labels = np.unique(cluster_labels)
for i in unique_cluster_labels:
ax.scatter(data[cluster_labels == i, 0],
data[cluster_labels == i, 1],
label = i,
alpha = 0.2)
ax.scatter(centroids[:, 0], centroids[:, 1], c='black')
ax.set_title("Final Result of K-means Clustering")
plt.show()
!ec
And here is an animation of the progression of the algorithm.
MOVIE: [clustering_example_images/simple_clustering.gif]
The completed code, up to this point, is wrapped in a function for convenience
later.
!bc pycod
def naive_kmeans(data, n_clusters=4, max_iterations=100, tolerance=1e-8):
start_time = time.time()
n_samples, dimensions = data.shape
n_clusters = 4
#np.random.seed(2021)
centroids = data[np.random.choice(n_samples, n_clusters, replace=False), :]
distances = np.zeros((n_samples, n_clusters))
for k in range(n_clusters):
for n in range(n_samples):
dist = 0
for d in range(dimensions):
dist += np.abs(data[n, d] - centroids[k, d])**2
distances[n, k] = dist
cluster_labels = np.zeros(n_samples, dtype='int')
for n in range(n_samples):
smallest = 1e10
smallest_row_index = 1e10
for k in range(n_clusters):
if distances[n, k] < smallest:
smallest = distances[n, k]
smallest_row_index = k
cluster_labels[n] = smallest_row_index
for iteration in range(max_iterations):
prev_centroids = centroids.copy()
for k in range(n_clusters):
vector_mean = np.zeros(dimensions)
mean_divisor = 0
for n in range(n_samples):
if cluster_labels[n] == k:
vector_mean += data[n, :]
mean_divisor += 1
centroids[k, :] = vector_mean / mean_divisor
for k in range(n_clusters):
for n in range(n_samples):
dist = 0
for d in range(dimensions):
dist += np.abs(data[n, d] - centroids[k, d])**2
distances[n, k] = dist
for n in range(n_samples):
smallest = 1e10
smallest_row_index = 1e10
for k in range(n_clusters):
if distances[n, k] < smallest:
smallest = distances[n, k]
smallest_row_index = k
cluster_labels[n] = smallest_row_index
centroid_difference = np.sum(np.abs(centroids - prev_centroids))
if centroid_difference < tolerance:
print(f'Converged at iteration {iteration}')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
print(f'Did not converge in {max_iterations} iterations')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
!ec
Now there are a few glaring improvements to be done here. First of all is
organizing things into functions for better readability. Second is getting rid
of the small inefficiencies like manually calculating distances and argmin. And
finally, we need to optimize for better run-time. It's like we always say: the
best way of looping in Python is to not loop in Python. Let us tackle the first
two improvements.
=== Towards a More Numpythonic Code ===
!bc pycod
def get_distances_to_clusters(data, centroids):
"""
Function that for each cluster finds the squared Euclidean distance
from every data point to the cluster centroid and returns a numpy array
containing the distances such that distance[i, j] means the distance between
the i-th point and the j-th centroid.
Inputs:
data (np.array): with dimensions (n_samples x dim)
centroids (np.array): with dimensions (n_clusters x dim)
Returns:
distances (np.array): with dimensions (n_samples x n_clusters)
"""
n_samples, dimensions = data.shape
n_clusters = centroids.shape[0]
distances = np.zeros((n_samples, n_clusters))
for k in range(n_clusters):
for i in range(n_samples):
dist = 0
for j in range(dimensions):
dist += np.abs(data[i, j] - centroids[k, j])**2
distances[i, k] = dist
return distances
def assign_points_to_clusters(distances):
"""
Function to assign each data point to the cluster to which it is the closest
based on the squared Euclidean distance from the get_distances_to_clusters
method.
Inputs:
distances (np.array): with dimensions (n_samples x n_clusters)
Returns:
cluster_labels (np.array): with dimensions (n_samples)
"""
cluster_labels = np.argmin(distances, axis=1)
return cluster_labels
def k_means(data, n_clusters=4, max_iterations=100, tolerance=1e-8):
"""
Naive implementation of the k-means clustering algorithm. A short summary of
the algorithm is as follows: we randomly initialize k centroids / means.
Then we assign, using the squared Euclidean distance, every data-point to a
cluster. We then update the position of the k centroids / means, and repeat
until convergence or we reach our desired maximum iterations. The method
returns the cluster assignments of our data-points and a sequence of
centroids.
Inputs:
data (np.array): with dimesions (n_samples x dim)
n_clusters (int): hyperparameter which depends on dataset
max_iterations (int): hyperparameter which depends on dataset
tolerance (float): convergence measure
Returns:
cluster_labels (np.array): with dimension (n_samples)
centroid_list (list): list of centroids (np.array)
with dimensions (n_clusters x dim)
"""
start_time = time.time()
n_samples, dimensions = data.shape
#np.random.seed(2021)
centroids = data[np.random.choice(len(data), n_clusters, replace=False), :]
distances = get_distances_to_clusters(data, centroids)
cluster_labels = assign_points_to_clusters(distances)
for iteration in range(max_iterations):
prev_centroids = centroids.copy()
for k in range(n_clusters):
vector_mean = np.zeros(dimensions)
mean_divisor = 0
for n in range(n_samples):
if cluster_labels[n] == k:
vector_mean += data[n, :]
mean_divisor += 1
# And update according to the new means
centroids[k, :] = vector_mean / mean_divisor
distances = get_distances_to_clusters(data, centroids)
cluster_labels = assign_points_to_clusters(distances)
centroid_difference = np.sum(np.abs(centroids - prev_centroids))
if centroid_difference < tolerance:
print(f'Converged at iteration: {iteration}')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
print(f'Did not converge in {max_iterations} iterations')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
cluster_labels, centroids = k_means(data)
!ec
Notice we added timing to our code. Although, the result of timing once will not
be optimal, it can still give us hints about the order of magnitude of our
subsequent improvements
So we see an improvement from just switching to numpy's argmin function. There
is a very nice tool (or category of tools) called profilers. These can be
utilized to make clearer which improvements to our code we should care most
about here is an "excellent source": "https://ipython-books.github.io/42-profiling-your-code-easily-with-cprofile-and-ipython/"
on the topic. Even before optimizing we can understand which parts of our code
will be taking the most of the run-time. It will be the longest Python loop,
i.e. the loop over all the samples. Nonetheless, let us do some profiling!
!bc pycod
test_data = generate_simple_clustering_dataset(n_points=10000, plotting=False)
%prun -l 10 cluster_labels, centroids = k_means(test_data)
!ec
Here we can see the reason for profiling. We now know for certain a lot can be
gained just by vectorizing our distance function. Ideally we wish to perform
most of our loops in numpy, i.e. C. To do this we need our array shapes to match
and clever reshaping will let us do so.
!bc pycod
def np_get_distances_to_clusters(data, centroids):
"""
Squared Euclidean distance between all data-points and every centroid. For
the function to work properly it needs data and centroids to be numpy
broadcastable. We sum along the dimension axis.
Inputs:
data (np.array): with dimensions (samples x 1 x dim)
centroids (np.array): with dimensions (1 x n_clusters x dim)
Returns:
distances (np.array): with dimensions (samples x n_clusters)
"""
distances = np.sum(np.abs((data - centroids))**2, axis=2)
return distances
def np_assign_points_to_clusters(distances):
"""
Assigning each data-point to a cluster given an array distances containing
the squared Euclidean distance from every point to each centroid. We do
np.argmin along the cluster axis to find the closest cluster. Returns a
numpy array with corresponding labels.
Inputs:
distances (np.array): with dimensions (samples x n_clusters)
Returns:
cluster_labels (np.array): with dimensions (samples x None)
"""
cluster_labels = np.argmin(distances, axis=1)
return cluster_labels
def np_k_means(data, n_clusters=4, max_iterations=100, tolerance=1e-8):
"""
Numpythonic implementation of the k-means clusting algorithm.
Inputs:
data (np.array): with dimesions (n_samples x dim)
n_clusters (int): hyperparameter which depends on dataset
max_iterations (int): hyperparameter which depends on dataset
tolerance (float): convergence measure
Returns:
cluster_labels (np.array): with dimension (n_samples)
centroid_list (list): list of centroids (np.array)
with dimensions (n_clusters x dim)
"""
start_time = time.time()
n_samples, dimensions = data.shape
#np.random.seed(2021)
centroids = data[np.random.choice(len(data), n_clusters, replace=False), :]
distances = np_get_distances_to_clusters(np.reshape(data,
(n_samples, 1, dimensions)),
np.reshape(centroids,
(1, n_clusters, dimensions)))
cluster_labels = np_assign_points_to_clusters(distances)
for iteration in range(max_iterations):
prev_centroids = centroids.copy()
for k in range(n_clusters):
points_in_cluster = data[cluster_labels == k]
mean_vector = np.mean(points_in_cluster, axis=0)
centroids[k] = mean_vector
distances = np_get_distances_to_clusters(np.reshape(data,
(n_samples, 1, dimensions)),
np.reshape(centroids,
(1, n_clusters, dimensions)))
cluster_labels = np_assign_points_to_clusters(distances)
centroid_difference = np.sum(np.abs(centroids - prev_centroids))
if centroid_difference < tolerance:
print(f'Converged at iteration: {iteration}')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
print(f'Did not converge in {max_iterations} iterations')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
!ec
When working towards becoming a data scientist using Python this last step is
arguably one of the most important. Thinking of ways to avoid explicitly looping
by adding dimensions to our arrays in such a way that they become broadcastable
using numpy (also tensorflow and many others). Let us take a look at our the
fruits of our labor.
!bc pycod
cluster_labels, centroids = np_k_means(data)
!ec
So our new code is around _two orders of magnitude_ faster on this limited test
example!
So let us recap what we have learned. The $k$-means algorithm works by
iteratively updating $k$ means according to the squared Euclidean distance to
our data points. One way of thinking about it is in terms of compressing the
original high dimensional data to a lower dimension. The process of which,
*hopefully*, captures something about the structure of the original data.
Phrased according to our example, we start of with 4000 '*groups*' which we
'compress' to just four '*groups*'. We have first hand seen the unreasonable
inefficiency of pure Python loops and learned to think about loops as just
another way of looking at the dimensionality of our data.
=== Onto Bigger and Better Things ===
As is often the case, someone has already done everything we have just done
(often but not always better). Let us take a look at "scikit-learn's
implementation": "https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html".
We define a wrapper function for timing and benchmarking purposes later.
!bc pycod
def skl_kmeans(data, n_clusters=4, max_iterations=100, tolerance=1e-8):
start_time = time.time()
#np.random.seed(2021)
kmeans = KMeans(n_clusters=n_clusters, max_iter=max_iterations,
tol=tolerance).fit(data)
print(f'Converged at iteration {kmeans.n_iter_}')
print(f'Runtime: {time.time() - start_time} seconds')
return kmeans.labels_, kmeans.cluster_centers_
!ec
Now let us get an idea how we hold up against the professionals.
!bc pycod
cluster_labels, centroids = skl_kmeans(data)
!ec
Here we get another important lesson, there is a cost for convenience. In
addition, there might be various tests and other things enforcing robustness in
the sklearn implementation which are not present in our numpy code. However,
the result still speaks for themselves.
The final implementation we will look at is one in tensorflow. Now while
tensorflow is not strictly meant for things like this. It is natively highly
parallelized with functionality out of the box to run on GPU's. The following
code is based on, and adapted from this "blog-post": "https://www.altoros.com/blog/using-k-means-clustering-in-tensorflow/"
!bc pycod
def tf_kmeans(data, n_clusters=4, max_iterations=100, tolerance=1e-8):
start_time = time.time()
data = tf.constant(data)
centroids = tf.constant(tf.slice(tf.random.shuffle(data),
begin=[0, 0],
size=[n_clusters, -1]))
@tf.function
def update_centroids(data, centroids):
data_expanded = tf.expand_dims(data, 0)
centroids_expanded = tf.expand_dims(centroids, 1)
distances = tf.reduce_sum(
tf.square(
tf.subtract(data_expanded, centroids_expanded)), 2)
cluster_labels = tf.argmin(distances, 0)
means = []
for k in range(n_clusters):
temp = tf.reshape(tf.where(tf.equal(cluster_labels, k)),
shape=[1, -1])
temp = tf.gather(data, temp, validate_indices=None)
temp = tf.reduce_mean(temp, axis=[1])
means.append(temp)
updated_centroids = tf.concat(means, 0)
return cluster_labels, updated_centroids
for iteration in range(max_iterations):
prev_centroids = tf.identity(centroids)
cluster_labels, centroids = update_centroids(data, centroids)
if tf.reduce_sum(
tf.abs(
tf.subtract(prev_centroids, centroids))) < tolerance:
print(f'Converged at iteration {iteration}')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
print(f'Did not converge in {max_iterations} iterations')
print(f'Runtime: {time.time() - start_time} seconds')
return cluster_labels, centroids
!ec
!bc pycod
cluster_labels, centroids = tf_kmeans(data)
!ec
_Note_ tensorflow has a different seed than numpy, so direct comparison of just
one run of the code (bad practice) is not doable. Also notice the @tf.function.
This was introduced in tensorflow 2 and signifies that the function is
'compiled'. We will look at benchmarking all our methods next.
But first we plot or result to see we get something which looks correct. We also
take this opportunity to define a function for our plotting.
!bc pycod
def make_plot(data, cluster_labels, centroids, method_string):
"""
Simple plot function
Inputs:
data (np.array like) with dimensions: (n_samples x dim)
cluster_labels (np.array like) with dimensions: (n_samples)
centroids (np.array like) with dimensions: (n_clusters x dim)
"""
fig = plt.figure()
ax = fig.add_subplot()
unique_cluster_labels = np.unique(cluster_labels)
for i in unique_cluster_labels:
ax.scatter(data[cluster_labels == i, 0],
data[cluster_labels == i, 1],
label = i,
alpha = 0.2)
ax.scatter(centroids[:, 0], centroids[:, 1], c='black')
ax.set_title(f"K Means Clustering using {method_string} Method")
plt.show()
!ec
!bc pycod
make_plot(data, cluster_labels, centroids, "Tensorflow")
!ec
Looks good! To summarize, we have now seen how to code the most basic Python
version, how to optimize and use a profiler, we ended up with a really quick
numpythonic version of our code. Then we looked at two implementations of the
$k$ means algorithm using higher level libraries scikit-learn and tensorflow.
===== Benchmarking and Testing Our Code =====
Now that we have developed our code, the next logical step is testing how well
our different implementations fare against each other in various tests. We start
with doing a Monte Carlo run over all our methods to produce numerical estimates
for the performance of our implementations.
!bc pycod
monte_carlo_cycles = 10
n_data_points = np.logspace(10, 100, num=2)
method_dict = {
'Python': naive_kmeans,
'Numpy': np_k_means,
'Scikit-Learn': skl_kmeans,
'Tensorflow': tf_kmeans
}
times_array = np.zeros((4, len(n_data_points), monte_carlo_cycles))
for i, (method_string, method) in enumerate(method_dict.items()):
for j, n_points in enumerate(n_data_points):
data = generate_simple_clustering_dataset(n_points=n_points,
plotting=False)
for k in range(monte_carlo_cycles):
toc = time.time()
cluster_labels, centroids = method(data)
tic = time.time() - toc
times_array[i, j, k] = tic
!ec
!bc pycod
means = np.mean(times_array, axis=2)
std = np.std(times_array, axis=2)
fig = plt.figure()
ax = fig.add_subplot()
for i, method_string in enumerate(method_dict):
print(i, method_string)
ax.errorbar(n_data_points, means[i, :], yerr=std[i, :],
label=method_string,
markersize=10,
capsize=5)
ax.set_xlabel('$n$ data points')
ax.set_ylabel('$t$ seconds')
fig.legend()
plt.show()
!ec
The naive Python implementation is too slow to see the nuance in the others, we
plot again without it
!bc pycod
fig = plt.figure()
ax = fig.add_subplot()
# this turned out a bit hacky
for i, method_string in enumerate(list(method_dict.keys())[1:]):
ax.errorbar(n_data_points, means[i+1, :], yerr=std[i+1, :],
label=method_string,
markersize=10,
capsize=5)
ax.set_xlabel('$n$ data points')
ax.set_ylabel('$t$ seconds')
fig.legend()
plt.show()
!ec
_THERE NEEDS TO BE A DISCUSSION OF RESULTS HERE WHEN I PRODUCED THEM_
_POTENTIALLY NEED TO DO SOMETHING LIKE THE SKL CLUSTERING EXAMPLE_
_THAT COULD LEAD TO DISCUSSION OF STRENGTHS AND WEAKNESSES OF METHOD_
===== Vector Quantization: Actual Example Use Case =====
Before we wrap up this topic we will consider the topic of _vector quantization_
. It is a technique which comes from signal processing theory. And is a basis
for some kinds of lossy compression. This is obvious because the $k$-means
algorithm is non-injective (many-to-one). Take notice of the fact that while the
$k$-means algorithm is an example of a vector quantization, there also exists
other vector quantization algorithms. We will be doing compression on an image,
but in principle this method can be applied to any form of data. We leave it to
the interested reader to try doing compression on something else.
We load in our image. And do some array manipulation to make the image size more
manageable. The way this will work is our centroids will represent colors, and
clustering means that we 'force' pixels that originally were not that color to
'turn' that color. Some cool art can probably be made using this. Do not
hesitate to try for yourself on another image!!
!bc pycod
example_image = image.imread('clustering_example_images/some_image.jpg') / 255
print(f'Initial shape: {example_image.shape}')
x, y, rgb = example_image.shape
example_image = example_image[0:500, 200:800, :]
reshaped_image = np.reshape(example_image, (500 * 600, rgb))
print(f'New shape: {reshaped_image.shape}')
!ec
FIGURE: [clustering_example_images/some_image.jpg]
Let's first try clustering with different numbers of centroids and look at the
effect. We will use the numpy implementation because it was the fastest.
!bc pycod
cluster_labels_4, centroids_4 = np_k_means(reshaped_image, n_clusters=4,
tolerance=1e-4)
cluster_labels_6, centroids_6 = np_k_means(reshaped_image, n_clusters=6,
tolerance=1e-4)
cluster_labels_8, centroids_8 = np_k_means(reshaped_image, n_clusters=8,
tolerance=1e-4)
compressed_image_4 = reshaped_image.copy()
compressed_image_6 = reshaped_image.copy()
compressed_image_8 = reshaped_image.copy()
for i in range(len(reshaped_image)):
compressed_image_4[i] = centroids_4[cluster_labels_4[i]]
compressed_image_6[i] = centroids_6[cluster_labels_6[i]]
compressed_image_8[i] = centroids_8[cluster_labels_8[i]]
!ec
!bc pycod
fig, ax = plt.subplots(2, 2, figsize=(8, 5))
ax[0, 0].imshow(example_image)
ax[0, 0].set_title('Original')
ax[0, 0].axis('off')
ax[0, 1].imshow(np.reshape(compressed_image_4, (500, 600, 3)))
ax[0, 1].set_title('$n_c = 4$')
ax[0, 1].axis('off')
ax[1, 0].imshow(np.reshape(compressed_image_6, (500, 600, 3)))
ax[1, 0].set_title('$n_c = 6$')
ax[1, 0].axis('off')
ax[1, 1].imshow(np.reshape(compressed_image_8, (500, 600, 3)))
ax[1, 1].set_title('$n_c = 8$')
ax[1, 1].axis('off')
plt.show()
!ec
And we have randomly stumbled upon how stencils for spray painting can be made.
Definitely a lot of potential for cool art here!
doc/src/Clustering/clean.sh
Executable
+3
View File
@@ -0,0 +1,3 @@
#!/bin/sh
doconce clean
rm -rf *.pdf *.tex ipynb*.tar.gz *.html ._*.html *~ reveal.js Trash README.txt
doc/src/Clustering/make.sh
Executable
+79
View File
@@ -0,0 +1,79 @@
#!/bin/sh
set -x
function system {
"$@"
if [ $? -ne 0 ]; then
echo "make.sh: unsuccessful command $@"
echo "abort!"
exit 1
fi
}
if [ $# -eq 0 ]; then
echo 'bash make.sh slides1|slides2'
exit 1
fi
name=$1
rm -f *.tar.gz
opt="--encoding=utf-8"
# Note: Makefile examples contain constructions like ${PROG} which
# looks like Mako constructions, but they are not. Use --no_mako
# to turn off Mako processing.
opt="--no_mako"
rm -f *.aux
html=${name}-reveal
system doconce format html $name --pygments_html_style=perldoc --keep_pygments_html_bg --html_links_in_new_window --html_output=$html $opt
system doconce slides_html $html reveal --html_slide_theme=beige
# Plain HTML documents
html=${name}-solarized
system doconce format html $name --pygments_html_style=perldoc --html_style=solarized3 --html_links_in_new_window --html_output=$html $opt
system doconce split_html $html.html --method=space10
html=${name}
system doconce format html $name --pygments_html_style=default --html_style=bloodish --html_links_in_new_window --html_output=$html $opt
system doconce split_html $html.html --method=space10
# Bootstrap style
html=${name}-bs
system doconce format html $name --html_style=bootstrap --pygments_html_style=default --html_admon=bootstrap_panel --html_output=$html $opt
system doconce split_html $html.html --method=split --pagination --nav_button=bottom
# IPython notebook
system doconce format ipynb $name $opt
# Publish
dest=../../pub
if [ ! -d $dest/$name ]; then
mkdir $dest/$name
mkdir $dest/$name/html
mkdir $dest/$name/ipynb
fi
cp -r ${name}*.html ._${name}*.html reveal.js $dest/$name/html
# Figures: cannot just copy link, need to physically copy the files
if [ -d fig-${name} ]; then
if [ ! -d $dest/$name/html/fig-$name ]; then
mkdir $dest/$name/html/fig-$name
fi
cp -r fig-${name}/* $dest/$name/html/fig-$name
fi
cp ${name}.ipynb $dest/$name/ipynb
ipynb_tarfile=ipynb-${name}-src.tar.gz
if [ ! -f ${ipynb_tarfile} ]; then
cat > README.txt <<EOF
This IPython notebook ${name}.ipynb does not require any additional
programs.
EOF
tar czf ${ipynb_tarfile} README.txt
fi
cp ${ipynb_tarfile} $dest/$name/ipynb
doc/src/GANS/GenerativeAdversarialNetworks.do.txt
+535
View File
@@ -0,0 +1,535 @@
======= Generative Models =======
_Generative models_ describe a class of statistical models that are a contrast
to _discriminative models_. Informally we say that generative models can
generate new data instances while discriminative models discriminate between
different kinds of data instances. A generative model could generate new photos
of animals that look like 'real' animals while a discriminative model could tell
a dog from a cat. More formally, given a data set $x$ and a set of labels /
targets $y$. Generative models capture the joint probability $p(x, y)$, or
just $p(x)$ if there are no labels, while discriminative models capture the
conditional probability $p(y | x)$. Discriminative models generally try to draw
boundaries in the data space (often high dimensional), while generative models
try to model how data is placed throughout the space.
===== Generative Adversarial Networks =====
_Generative Adversarial Networks_ are a type of unsupervised machine learning
algorithm proposed by "Goodfellow et. al": "https://arxiv.org/pdf/1406.2661.pdf"
in 2014 (Read the paper first it's only 6 pages). The simplest formulation of
the model is based on a game theoretic approach, *zero sum game*, where we pit
two neural networks against one another. We define two rival networks, one
generator $g$, and one discriminator $d$. The generator directly produces
samples
!bt
\begin{equation}
x = g(z; \theta^{(g)})
\end{equation}
!et
The discriminator attempts to distinguish between samples drawn from the
training data and samples drawn from the generator. In other words, it tries to
tell the difference between the fake data produced by $g$ and the actual data
samples we want to do prediction on. The discriminator outputs a probability
value given by
!bt
\begin{equation}
d(x; \theta^{(d)})
\end{equation}
!et
indicating the probability that $x$ is a real training example rather than a
fake sample the generator has generated. The simplest way to formulate the
learning process in a generative adversarial network is a zero-sum game, in
which a function
!bt
\begin{equation}
v(\theta^{(g)}, \theta^{(d)})
\end{equation}
!et
determines the reward for the discriminator, while the generator gets the
conjugate reward
!bt
\begin{equation}
-v(\theta^{(g)}, \theta^{(d)})
\end{equation}
!et
During learning both of the networks maximize their own reward function, so that
the generator gets better and better at tricking the discriminator, while the
discriminator gets better and better at telling the difference between the fake
and real data. The generator and discriminator alternate on which one trains at
one time (i.e. for one epoch). In other words, we keep the generator constant
and train the discriminator, then we keep the discriminator constant to train
the generator and repeat. It is this back and forth dynamic which lets GANs
tackle otherwise intractable generative problems. As the generator improves with
training, the discriminator's performance gets worse because it cannot easily
tell the difference between real and fake. If the generator ends up succeeding
perfectly, the the discriminator will do no better than random guessing i.e.
50\%. This progression in the training poses a problem for the convergence
criteria for GANs. The discriminator feedback gets less meaningful over time,
if we continue training after this point then the generator is effectively
training on junk data which can undo the learning up to that point. Therefore,
we stop training when the discriminator starts outputting $1/2$ everywhere.
At convergence we have
!bt
\begin{equation}
g^* = \underset{g}{\mathrm{argmin}}\hspace{2pt}
\underset{d}{\mathrm{max}}v(\theta^{(g)}, \theta^{(d)})
\end{equation}
!et
The default choice for $v$ is
!bt
\begin{equation}
v(\theta^{(g)}, \theta^{(d)}) = \mathbb{E}_{x\sim p_\mathrm{data}}\log d(x)
+ \mathbb{E}_{x\sim p_\mathrm{model}}
\log (1 - d(x))
\end{equation}
!et
The main motivation for the design of GANs is that the learning process requires
neither approximate inference (variational autoencoders for example) nor
approximation of a partition function. In the case where
!bt
\begin{equation}
\underset{d}{\mathrm{max}}v(\theta^{(g)}, \theta^{(d)})
\end{equation}
!et
is convex in $\theta^{(g)} then the procedure is guaranteed to converge and is
asymptotically consistent
( "Seth Lloyd on QuGANs": "https://arxiv.org/pdf/1804.09139.pdf" ). This is in
general not the case and it is possible to get situations where the training
process never converges because the generator and discriminator chase one
another around in the parameter space indefinitely. A much deeper discussion on
the currently open research problem of GAN convergence is available
"here": "https://www.deeplearningbook.org/contents/generative_models.html". To
anyone interested in learning more about GANs it is a highly recommended read.
Direct quote: "In this best-performing formulation, the generator aims to
increase the log probability that the discriminator makes a mistake, rather than
aiming to decrease the log probability that the discriminator makes the correct
prediction." "Another interesting read": "https://arxiv.org/abs/1701.00160"
===== Writing Our First Generative Adversarial Network =====
Let us now move on to actually implementing a GAN in tensorflow. We will study
the performance of our GAN on the MNIST dataset. This code is based on and
adapted from the
"google tutorial": "https://www.tensorflow.org/tutorials/generative/dcgan"
First we import our libraries
!bc pycod
import os
import time
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
from tensorflow.keras import layers
from tensorflow.keras.utils import plot_model
!ec
Next we define our hyperparameters and import our data the usual way
!bc pycod
BUFFER_SIZE = 60000
BATCH_SIZE = 256
EPOCHS = 30
data = tf.keras.datasets.mnist.load_data()
(train_images, train_labels), (test_images, test_labels) = data
train_images = np.reshape(train_images, (train_images.shape[0],
28,
28,
1)).astype('float32')
# we normalize between -1 and 1
train_images = (train_images - 127.5) / 127.5
training_dataset = tf.data.Dataset.from_tensor_slices(
train_images).shuffle(BUFFER_SIZE).batch(BATCH_SIZE)
!ec
Let's have a quick look
!bc pycod
plt.imshow(train_images[0], cmap='Greys')
plt.show()
!ec
Now we define our two models. This is where the 'magic' happens. There are a
huge amount of possible formulations for both models. A lot of engineering and
trial and error can be done here to try to produce better performing models. For
more advanced GANs this is by far the step where you can 'make or break' a
model.
We start with the generator. As stated in the introductory text the generator
$g$ upsamples from a random sample to the shape of what we want to predict. In
our case we are trying to predict MNIST images ($28\times 28$ pixels).
!bc pycod
def generator_model():
"""
The generator uses upsampling layers tf.keras.layers.Conv2DTranspose() to
produce an image from a random seed. We start with a Dense layer taking this
random sample as an input and subsequently upsample through multiple
convolutional layers.
"""
# we define our model
model = tf.keras.Sequential()
# adding our input layer. Dense means that every neuron is connected and
# the input shape is the shape of our random noise. The units need to match
# in some sense the upsampling strides to reach our desired output shape.
# we are using 100 random numbers as our seed
model.add(layers.Dense(units=7*7*BATCH_SIZE,
use_bias=False,
input_shape=(100, )))
# we normalize the output form the Dense layer
model.add(layers.BatchNormalization())
# and add an activation function to our 'layer'. LeakyReLU avoids vanishing
# gradient problem
model.add(layers.LeakyReLU())
model.add(layers.Reshape((7, 7, BATCH_SIZE)))
assert model.output_shape == (None, 7, 7, BATCH_SIZE)
# even though we just added four keras layers we think of everything above
# as 'one' layer
# next we add our upscaling convolutional layers
model.add(layers.Conv2DTranspose(filters=128,
kernel_size=(5, 5),
strides=(1, 1),
padding='same',
use_bias=False))
model.add(layers.BatchNormalization())
model.add(layers.LeakyReLU())
assert model.output_shape == (None, 7, 7, 128)
model.add(layers.Conv2DTranspose(filters=64,
kernel_size=(5, 5),
strides=(2, 2),
padding='same',
use_bias=False))
model.add(layers.BatchNormalization())
model.add(layers.LeakyReLU())
assert model.output_shape == (None, 14, 14, 64)
model.add(layers.Conv2DTranspose(filters=1,
kernel_size=(5, 5),
strides=(2, 2),
padding='same',
use_bias=False,
activation='tanh'))
assert model.output_shape == (None, 28, 28, 1)
return model
!ec
And there we have our 'simple' generator model. Now we move on to defining our
discriminator model $d$, which is a convolutional neural network based image
classifier.
!bc pycod
def discriminator_model():
"""
The discriminator is a convolutional neural network based image classifier
"""
# we define our model
model = tf.keras.Sequential()
model.add(layers.Conv2D(filters=64,
kernel_size=(5, 5),
strides=(2, 2),
padding='same',
input_shape=[28, 28, 1]))
model.add(layers.LeakyReLU())
# adding a dropout layer as you do in conv-nets
model.add(layers.Dropout(0.3))
model.add(layers.Conv2D(filters=128,
kernel_size=(5, 5),
strides=(2, 2),
padding='same'))
model.add(layers.LeakyReLU())
# adding a dropout layer as you do in conv-nets
model.add(layers.Dropout(0.3))
model.add(layers.Flatten())
model.add(layers.Dense(1))
return model
!ec
Let us take a look at our models. _Note_: double click images for bigger view.
!bc pycod
generator = generator_model()
plot_model(generator, show_shapes=True, rankdir='LR')
!ec
!bc pycod
discriminator = discriminator_model()
plot_model(discriminator, show_shapes=True, rankdir='LR')
!ec
Next we need a few helper objects we will use in training
!bc pycod
cross_entropy = tf.keras.losses.BinaryCrossentropy(from_logits=True)
generator_optimizer = tf.keras.optimizers.Adam(1e-4)
discriminator_optimizer = tf.keras.optimizers.Adam(1e-4)
!ec
The first object, *cross_entropy* is our loss function and the two others are
our optimizers. Notice we use the same learning rate for both $g$ and $d$. This
is because they need to improve their accuracy at approximately equal speeds to
get convergence (not necessarily exactly equal). Now we define our loss
functions
!bc pycod
def generator_loss(fake_output):
loss = cross_entropy(tf.ones_like(fake_output), fake_output)
return loss
!ec
!bc pycod
def discriminator_loss(real_output, fake_output):
real_loss = cross_entropy(tf.ones_like(real_output), real_output)
fake_loss = cross_entropy(tf.zeros_liks(fake_output), fake_output)
total_loss = real_loss + fake_loss
return total_loss
!ec
Next we define a kind of seed to help us compare the learning process over
multiple training epochs.
!bc pycod
noise_dimension = 100
n_examples_to_generate = 16
seed_images = tf.random.normal([n_examples_to_generate, noise_dimension])
!ec
Now we have everything we need to define our training step, which we will apply
for every step in our training loop. Notice the @tf.function flag signifying
that the function is tensorflow 'compiled'. Removing this flag doubles the
computation time.
!bc pycod
@tf.function
def train_step(images):
noise = tf.random.normal([BATCH_SIZE, noise_dimension])
with tf.GradientTape() as gen_tape, tf.GradientTape() as disc_tape:
generated_images = generator(noise, training=True)
real_output = discriminator(images, training=True)
fake_output = discriminator(generated_images, training=True)
gen_loss = generator_loss(fake_output)
disc_loss = discriminator_loss(real_output, fake_output)
gradients_of_generator = gen_tape.gradient(gen_loss,
generator.trainable_variables)
gradients_of_discriminator = disc_tape.gradient(disc_loss,
discriminator.trainable_variables)
generator_optimizer.apply_gradients(zip(gradients_of_generator,
generator.trainable_variables))
discriminator_optimizer.apply_gradients(zip(gradients_of_discriminator,
discriminator.trainable_variables))
return gen_loss, disc_loss
!ec
Next we define a helper function to produce an output over our training epochs
to see the predictive progression of our generator model. _Note_: I am including
this code here, but comment it out in the training loop.
!bc pycod
def generate_and_save_images(model, epoch, test_input):
# we're making inferences here
predictions = model(test_input, training=False)
fig = plt.figure(figsize=(4, 4))
for i in range(predictions.shape[0]):
plt.subplot(4, 4, i+1)
plt.imshow(predictions[i, :, :, 0] * 127.5 + 127.5, cmap='gray')
plt.axis('off')
plt.savefig(f'./images_from_seed_images/image_at_epoch_{str(epoch).zfill(3)}.png')
plt.close()
#plt.show()
!ec
Setting up checkpoints to periodically save our model during training so that
everything is not lost even if the program were to somehow terminate while
training.
!bc pycod
# Setting up checkpoints to save model during training
checkpoint_dir = './training_checkpoints'
checkpoint_prefix = os.path.join(checkpoint_dir, 'ckpt')
checkpoint = tf.train.Checkpoint(generator_optimizer=generator_optimizer,
discriminator_optimizer=discriminator_optimizer,
generator=generator,
discriminator=discriminator)
!ec
Now we define our training loop
!bc pycod
def train(dataset, epochs):
generator_loss_list = []
discriminator_loss_list = []
for epoch in range(epochs):
start = time.time()
for image_batch in dataset:
gen_loss, disc_loss = train_step(image_batch)
generator_loss_list.append(gen_loss.numpy())
discriminator_loss_list.append(disc_loss.numpy())
#generate_and_save_images(generator, epoch + 1, seed_images)
if (epoch + 1) % 15 == 0:
checkpoint.save(file_prefix=checkpoint_prefix)
print(f'Time for epoch {epoch} is {time.time() - start}')
#generate_and_save_images(generator, epochs, seed_images)
loss_file = './data/lossfile.txt'
with open(loss_file, 'w') as outfile:
outfile.write(str(generator_loss_list))
outfile.write('\n')
outfile.write('\n')
outfile.write(str(discriminator_loss_list))
outfile.write('\n')
outfile.write('\n')
!ec
To train simply call this function. _Warning_: this might take a long time so
there is a folder of a pretrained network already included in the repository.
!bc pycod
train(train_dataset, EPOCHS)
!ec
And here is the result of training our model for 100 epochs
MOVIE: [images_from_seed_images/generation.gif]
Now to avoid having to train and everything, which will take a while depending
on your computer setup we now load in the model which produced the above gif.
!bc pycod
checkpoint.restore(tf.train.latest_checkpoint(checkpoint_dir))
restored_generator = checkpoint.generator
restored_discriminator = checkpoint.discriminator
print(restored_generator)
print(restored_discriminator)
!ec
===== Exploring the Latent Space =====
So we have successfully loaded in our latest model. Let us now play around a bit
and see what kind of things we can learn about this model. Our generator takes
an array of 100 numbers. One idea can be to try to systematically change our
input. Let us try and see what we get
!bc pycod
def generate_latent_points(number=100, scale_means=1, scale_stds=1):
latent_dim = 100
means = scale_means * tf.linspace(-1, 1, num=latent_dim)
stds = scale_stds * tf.linspace(-1, 1, num=latent_dim)
latent_space_value_range = tf.random.normal([number, number],
means,
stds,
dtype=tf.float64)
return latent_space_value_range
def generate_images(latent_points):
# notice we set training to false because we are making inferences
generated_images = restored_generator(latent_space_value_range,
training=False)
return generated_images
!ec
!bc pycod
def plot_result(generated_images, number):
# obviously this assumes sqrt number is an int
fig, axs = plt.subplots(int(np.sqrt(number)), int(np.sqrt(number)),
figsize=(10, 10))
for i in range(int(np.sqrt(number))):
for j in range(int(np.sqrt(number))):
axs[i, j].imshow(generated_images[i*j], cmap='Greys')
axs[i, j].axis('off')
plt.show()
!ec
!bc pycod
generated_images = generate_images(generate_latent_points())
plot_result(generated_images, number)
!ec
Interesting! We see that the generator generates images that look like MNIST
numbers: $1, 4, 7, 9$. Let's try to tweak it a bit more to see if we are able
to generate a similar plot where we generate every MNIST number. Let us now try
to 'move' a bit around in the latent space. _Note_: decrease the plot number if
these following cells take too long to run on your computer.
!bc pycod
plot_number = 225
generated_images = generate_images(generate_latent_points(number=plot_number,
scale_means=5,
scale_stds=1))
plot_result(generated_images, plot_number)
generated_images = generate_images(generate_latent_points(number=plot_number,
scale_means=-5,
scale_stds=1))
plot_result(generated_images, plot_number)
generated_images = generate_images(generate_latent_points(number=plot_number,
scale_means=1,
scale_stds=5))
plot_result(generated_images, plot_number)
!ec
Again, we have found something interesting. *Moving* around using our means
takes us from digit to digit, while *moving* around using our standard
deviations seem to increase the number of different digits! In the last image
above, we can barely make out every MNIST digit. Let us make on last plot using
this information by upping the standard deviation of our Gaussian noises.
!bc pycod
plot_number = 400
generated_images = generate_images(generate_latent_points(number=plot_number,
scale_means=1,
scale_stds=10))
!ec
A pretty cool result! We see that our generator indeed has learned a
distribution which qualitatively looks a whole lot like the MNIST dataset.
===== Interpolating Between MNIST Digits =====
Another interesting way to explore the latent space of our generator model is by
interpolating between the MNIST digits. This section is largely based on
"this excellent blogpost": https://machinelearningmastery.com/how-to-interpolate-and-perform-vector-arithmetic-with-faces-using-a-generative-adversarial-network/
by Jason Brownlee.
So let us start