From 76009c48f89eb3cb4a3901d584430c925ea812fd Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Wed, 28 Jul 2021 23:27:13 +0200 Subject: [PATCH] added material --- doc/src/Clustering/Clustering.do.txt | 974 ++++++++++++++++++ doc/src/Clustering/clean.sh | 3 + doc/src/Clustering/make.sh | 79 ++ .../GANS/GenerativeAdversarialNetworks.do.txt | 535 ++++++++++ 4 files changed, 1591 insertions(+) create mode 100644 doc/src/Clustering/Clustering.do.txt create mode 100755 doc/src/Clustering/clean.sh create mode 100755 doc/src/Clustering/make.sh create mode 100644 doc/src/GANS/GenerativeAdversarialNetworks.do.txt diff --git a/doc/src/Clustering/Clustering.do.txt b/doc/src/Clustering/Clustering.do.txt new file mode 100644 index 000000000..91bb2cd3e --- /dev/null +++ b/doc/src/Clustering/Clustering.do.txt @@ -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! diff --git a/doc/src/Clustering/clean.sh b/doc/src/Clustering/clean.sh new file mode 100755 index 000000000..2e5da2c72 --- /dev/null +++ b/doc/src/Clustering/clean.sh @@ -0,0 +1,3 @@ +#!/bin/sh +doconce clean +rm -rf *.pdf *.tex ipynb*.tar.gz *.html ._*.html *~ reveal.js Trash README.txt diff --git a/doc/src/Clustering/make.sh b/doc/src/Clustering/make.sh new file mode 100755 index 000000000..4722a0918 --- /dev/null +++ b/doc/src/Clustering/make.sh @@ -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 <