By the end of this lecture, you should be able to:
To restate:
While clustering and classification both assign labels to observations, there is a key difference.
| Mapping | Supervised Learning Problem | Unsupervised Learning Problem |
|---|---|---|
| \(\hat f_{\mathcal D}(X): \mathbb R^p \to \mathbb R^k\) | Regression | Dimensionality Reduction |
| \(\hat f_{\mathcal D}(X): \mathbb R^p \to \{1, 2, \ldots, k\}\) | Classification | Clustering |
Consider the following synthetic dataset of \(n=140\) set of \(p=2\)-dimensional covariates: \(\mathcal D \in \{ X_1, \ldots X_{140} \}, X_i \in \mathbb R^2\).
library(Stat406)
library(mvtnorm)
library(ggplot2)
set.seed(406406406)
X1 <- rmvnorm(50, c(-1, 2), sigma = matrix(c(1, .5, .5, 1), 2))
X2 <- rmvnorm(40, c(2, -1), sigma = matrix(c(1.5, .5, .5, 1.5), 2))
X3 <- rmvnorm(40, c(4, 4))
data <- rbind(X1, X2, X3)
tibble(
x1 = data[, 1],
x2 = data[, 2],
) |> ggplot(aes(x = x1, y = x2)) +
geom_point(size = 2) +
theme(legend.position = "none")
# data: 140 x 2 matrix with our covariates
assignments <- kmeans(data, centers = 3)$cluster
tibble(
x1 = data[, 1],
x2 = data[, 2],
cluster = as.factor(assignments)
) |> ggplot(aes(x = x1, y = x2, colour = cluster)) +
geom_point(size = 2) +
scale_colour_manual(values = c("blue", "orange", "green"))
As with (almost) all other learning procedures we’ve discussed, we will formalize this clustering procedure as an optimization problem.
Assume we have fixed \(k\), i.e. the number of clusters we hope to partition our data into. (We’ll discuss how to choose \(k\) momentarily.)
What we are trying to learn: \(\hat f_{\mathcal D}(X) : \mathbb R^p \to \{ 1, \ldots, k \}\), or some function that maps covariates onto clusters.
For any given \(\hat f_{\mathcal D}(X) : \mathbb R^p \to \{ 1, \ldots, k \}\), let \(C_1, \ldots, C_k \subset \mathcal D\) represent the clusters induced by this mapping; i.e.
\[C_j = \{ X_i \in \mathcal D : \hat f_{\mathcal D}(X_i) = j \}\]
In other words, \(C_j\) is the set of all training covariates assigned to cluster \(j\).
The within-distance variation of this clustering is defined as:
\[ W_j := \frac{1}{|C_j|^2} \sum_{i, {i'} \in C_j} \| X_i - X_{i'} \|_2^2 \]
where \(|C_j|\) is the number of points in cluster \(j\). It is the average pairwise-distance between points in cluster \(j\).
We want to minimize the total within-cluster variation, but weight each cluster by its size so larger clusters contribute more.
Mathematically, this gives us:
\[ \mathrm{argmin}_{\hat f_{\mathcal D}(X)} \frac{1}{2} \sum_{j=1}^k \vert C_j \vert W_j = \frac{1}{2} \sum_{j=1}^k \frac{1}{\vert C_j \vert} \sum_{i, i' \in C_j} \Vert X_i - X_{i'} \Vert^2_2. \]
The fraction \(1/2\) is for mathematical convenience (it does not affect the maximization), as will become clear shortly.
This optimization problem is unwieldy, especially because we have to consider all pairwise distances between points in each cluster. (This is, worst case scenario, \(O(n^2)\) pairwise distances to consider!)
Fortunately, we can make use of the following identity to simplify from a pair-wise summation to a per-data-point summation:
\[ \frac{1}{2} \sum_{i, i' \in C_j} \| X_i - X_{i'} \|_2^2 = |C_j| \sum_{i \in C_j} \| X_i - \mu_j \|_2^2 \]
where \(\mu_j\) is the centroid of cluster \(j\), or the empirical average of the covariates within that cluster:
\[ \mu_j = \frac{1}{|C_j|} \sum_{i \in C_j} X_i. \]
If you’re confused about where this identity comes from, consider our favourite trick that we used when deriving the bias-variance tradeoff: adding and subtracting zero:
\[ \sum_{i, i' \in C_j} \| X_i - X_{i'} \|_2^2 = \sum_{i, i' \in C_j} \| (X_i {\color{blue}- \mu_j}) + ({\color{blue}\mu_j -} X_{i'}) \|_2^2. \]
Just like we did with the bias-variance tradeoff, expand the square and you’ll find that the cross terms vanish!
While we could go through this derivation step-by-step, it’s perhaps simpler if we replace the summations with expectations over the empirical distribution of points in cluster \(C_j\).
\[ \sum_{i, i' \in C_j} \| (X_i {\color{blue}- \mu_j}) + ({\color{blue}\mu_j -} X_{i'}) \|_2^2 = \vert C_j \vert^2 \mathbb E_{X, X' \sim \hat P_{C_j}} \left[ \| (X_i {\color{blue}- \mu_j}) + ({\color{blue}\mu_j -} X_{i'}) \|_2^2 \right] \]
We also note that \(\mu_j\) is the expectation of \(X\) under the empirical distribution \(\hat P_{C_j}\):
\[ \begin{aligned} &\vert C_j \vert^2 \mathbb E_{X, X' \sim \hat P_{C_j}} \left[\Vert (X_i {\color{blue}- \mu_j}) + ({\color{blue}\mu_j -} X_{i'}) \Vert_2^2 \right] \\ =& \vert C_j \vert^2 \mathbb E_{X, X' \sim \hat P_{C_j}} \left[\Vert (X_i - \mathbb E[X]) + (\mathbb E[X] - X_{i'}) \Vert_2^2 \right] \end{aligned} \]
and now this derivation should look exactly like what we did in the bias-variance tradeoff! Follow those steps to see the cross terms disappear!Intuitively, these centroids represent the “average” point in each cluster. They are represented on the plots below by black crosses.
centroids <- as.data.frame(kmeans(data, centers = 3)$centers)
ggplot() +
geom_point(data = tibble(
x1 = data[, 1],
x2 = data[, 2],
cluster = as.factor(assignments)
), aes(x = x1, y = x2, colour = cluster), size = 2) +
geom_point(data = centroids, aes(x = V1, y = V2), colour = "black", size = 5, shape = 3) +
scale_colour_manual(values = c("blue", "orange", "green"))
Even after simplifying the optimization problem using centroids, it is still impossible to solve it exactly. Instead, we rely on an iterative algorithm to approximate a solution.
We produce a series of clustering assignments
\[ f_{\mathcal D}^{(0)}(X), f_{\mathcal D}^{(1)}(X), f_{\mathcal D}^{(2)}(X), \ldots f_{\mathcal D}^{(t)}(X), \ldots \]
so that \(\sum_{j=1}^k |C_j| W_j\) decreases as \(t\) increases.
\(f_{\mathcal D}^{(0)}(X)\) will essentially be random:
\[ f_{\mathcal D}^{(0)}(X_i) = \mathrm{argmin}_{j \in \{1, \ldots, k\}} \| X_i - \mu_j \|_2^2. \]
\[ \mu_j = \frac{1}{|C_j|} \sum_{i \in C_j} X_i. \]
To get \(f_{\mathcal D}^{(1)}(X)\) from \(f_{\mathcal D}^{(0)}(X)\), we repeat steps 2 and 3 above.
More generally, we repeat steps 2 and 3 until the cluster centroids are stable (i.e. when the distance that the cluster centroids change is \(< \epsilon\) for some pre-defined constant \(\epsilon\). Most implementations will define this constant for you.)
Because k-means starts with a random initialization of the cluster centroids, it may converge to different clustering assignments on different runs. Most implementations of k-means (including R’s built-in kmeans function) will try many random initializations and return the clustering assignment with the lowest within-cluster variation.
However, note that you may get different clustering assignments on different runs of k-means, especially if the clusters are not well-separated.
In general, there are two key desiderata for our clusters. We’ve already talked about one of them:
\[ W := \frac{1}{2} \sum_{j=1}^k \vert C_j \vert W_j = \frac{1}{2}\sum_{j=1}^k \vert C_j \vert \sum_{i \in C_j} \left\| X_i - \mu_j \right\|_2^2 \]
\[ B := \frac{1}{2} \sum_{j=1}^k \sum_{j'=1}^k \frac{|C_j| |C_{j'}|}{n} \| \mu_j - \mu_{j'} \|_2^2. \]
As with the within-cluster variation, this summation over all pairs can be simplified to a summation over individual clusters:
\[ B = \sum_{j=1}^k |C_j| \| \mu_j - \mu \|_2^2 \]
where \(\mu = \frac{1}{n} \sum_{i=1}^n X_i\) is the overall mean of the covariates.
kmeans function in R automatically computes all of these variables
kmeans(data, centers = 3)K-means clustering with 3 clusters of sizes 49, 39, 42
Cluster means:
[,1] [,2]
1 -0.9528507 2.141750
2 1.8193830 -1.531834
3 4.0791847 3.836696
Clustering vector:
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[38] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2
[75] 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[112] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Within cluster sum of squares by cluster:
[1] 98.81053 111.78974 81.12076
(between_SS / total_SS = 80.2 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" You can guess what each of these variable names mean ;)
We want to minimize within-cluster variation and maximize between-cluster variation. However, these quantities will often be at odds with one another.
As the total number of clusters \(k\) increases…
A simple heuristic for choosing \(k\) is to maximize the (weighted) ratio between these two quantities:
\[ \mathrm{CH} := \frac{B/(k-1)}{W/(n-k)}. \]
For all \(k \in \{2, 3, \ldots, k_{\max}\}\) (where \(k_{\max} < n\) is some maximum number of clusters we want to consider):
Choose the number of clusters \(\hat k\) that maximizes the CH index.
Example clustering problem: You are studying birds, and you have collected measurements of their beak length, beak depth, wing length, and weight. You suspect there might be several distinct species of birds in your dataset. Since species are distinct categories, it makes sense to use clustering to group the birds into categories that may correspond to species or sub-species.
Example dimensionality reduction problem: You are studying mental health data. You have collected survey responses from individuals on various aspects of their mental health, including stress levels, anxiety, depression, and overall well-being. You want to reduce these multiple dimensions into a smaller set of underlying factors that capture the main variations in mental health. Since many mental health factors exist on a spectrum, dimensionality reduction is more appropriate here.
There are many alternative flavours of clustering that can be more appropriate for different applications/data types:
Hierarchical clustering organizes clusters in a tree-like structure, creating hierarchical relationships between the different clusters. This algorithm is especially popular for biological data.
Gaussian mixture models are a generalization of k-means where the cluster distances are generalized and cluster membership becomes “fuzzy” or probabilistic. You may learn about this method in a more advanced machine learning class.
Spectral clustering is similar to k-means clustering, except it is designed for graphical data (i.e. if your data is a social network, a “musical influence” chart, etc.)
You don’t need to know any of these, but be aware that they’re out there. Also, if k-means doesn’t immediately work for your data, with some googling you’ll likely find a variant that does!