Each row \(i\) is a point \(x^{(i)} \in \mathbb{R}^d\), here with \(d = 2\) columns. In supervised learning every row also had a label \(y^{(i)}\). Here there is none: the task is to find structure in the \(x^{(i)}\) alone.

Photograph from Color quantization for image reduction using k-means (Analytics Vidhya). The colours are computed here, in your browser.

A Gaussian mixture says how likely each point is to come from each blob. It is the subject of the next lecture; k-means, today, is its simplest relative.

Source: Towards Data Science, Hierarchical clustering explained

The height of each join is the distance between the two groups it merges. Cutting the tree at the red line gives the coloured groups.

Source: Towards Data Science, Hierarchical clustering explained

A cut at \(y_1\), \(y_2\) or \(y_3\) gives a different number of clusters from one tree, so the method does not need the number of clusters in advance.

Source: Towards Data Science, Hierarchical clustering explained

Cats and dogs at the top; big ears, calm behaviour, small or big size below. K-means gives one level only, and you choose which.

K-means describes each group by one point, its centre \(\mu_j\), and puts every \(x^{(i)}\) in the group of the nearest centre. The number of groups \(k\) is chosen by us.

• Assign. With the centres fixed, \(J\) is a sum of one term per point, and choosing the nearest centre makes each term as small as it can be. \( \lVert x^{(i)} - \mu_{c^{(i)}} \rVert^2 = \min_j \lVert x^{(i)} - \mu_j \rVert^2 \)
• Update. With the clusters fixed, the part of \(J\) for cluster \(j\) is a quadratic in \(\mu_j\). Its gradient is zero at the mean of the cluster. \( \nabla_{\mu_j} \sum_{i:\,c^{(i)} = j} \lVert x^{(i)} - \mu_j \rVert^2 = -2 \sum_{i:\,c^{(i)} = j} (x^{(i)} - \mu_j) = 0 \)
• So it stops. \(J\) never rises, and there are finitely many ways to assign \(n\) points to \(k\) clusters, so the assignments stop changing after a finite number of rounds.

See the CS229 notes, chapter 10 (cs229.stanford.edu), for the same argument.

Do it on paper with your neighbour before the next slide.

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)

# define the dataset
n = 50
c1 = np.array([-3, -3])
c2 = np.array([2, 2])
sig = 1.2
data = np.concatenate([c1 + sig*np.random.randn(n, 2), c2 + sig*np.random.randn(n, 2)])

# initialize the centroids (as floats: an integer array would round every mean)
centroids = np.array([[-2, 3], [2, -3]], dtype=float)

max_iterations = 100           # define the maximum number of iterations
convergence_threshold = 1e-4   # define the convergence threshold
iteration = 0                  # initialize the iteration counter

# previous centroids: a random value different from the current centroids
prev_centroids = np.random.rand(2, 2)

Two clouds of 50 points around \((-3, -3)\) and \((2, 2)\). Run this block, then the loop on the next slide.

while np.sum(np.abs(prev_centroids - centroids)) > convergence_threshold and iteration < max_iterations:
    prev_centroids = centroids.copy()          # store the previous centroids
    clusters = np.zeros(len(data), dtype=int)  # initialize the cluster assignments

    # assign each data point to the closest centroid
    for i, point in enumerate(data):
        distances = np.sum((point - centroids)**2, axis=1)
        clusters[i] = np.argmin(distances)

    # update the centroids using the mean of the points in each cluster
    for j in range(2):
        centroids[j] = np.mean(data[clusters == j], axis=0)

    iteration += 1
    J = np.sum((data - centroids[clusters])**2)
    print(f'iteration {iteration}: J = {J:.2f}')

fig, ax = plt.subplots()
ax.scatter(data[:, 0], data[:, 1], c=['C0' if c == 0 else 'C1' for c in clusters], s=12)
ax.plot(centroids[:, 0], centroids[:, 1], 'x', color='C2', markersize=10)
ax.set_title('Last result')
plt.show()

It prints iteration 1: J = 1177.29; iteration 2: J = 298.71; iteration 3: J = 298.71.

Share of the 342 birds in the cluster matched to their species, best of 50 starts.

import numpy as np
import pandas as pd

d = pd.read_csv('penguins.csv').dropna(subset=['bill_length_mm', 'bill_depth_mm'])
X = d[['bill_length_mm', 'bill_depth_mm']].values
Z = (X - X.mean(axis=0)) / X.std(axis=0)          # standardize: mean 0, spread 1

def kmeans(Z, k, rng):
    mu = Z[rng.choice(len(Z), k, replace=False)]   # k random birds as the first centres
    while True:
        D = ((Z[:, None, :] - mu[None, :, :])**2).sum(axis=2)   # n by k squared distances
        c = D.argmin(axis=1)                                     # assignment step
        new = np.array([Z[c == j].mean(axis=0) for j in range(k)])  # update step
        if np.allclose(new, mu):
            return c, mu, D[np.arange(len(Z)), c].sum()          # J, the distortion
        mu = new

rng = np.random.default_rng(0)
runs = [kmeans(Z, 3, rng) for _ in range(20)]    # 20 restarts, keep the lowest J
c, mu, J = min(runs, key=lambda r: r[2])
print(f'J = {J:.1f}')
print(pd.crosstab(c, d.species))                 # the species were never used above

It prints J = 186.9 and the table of clusters against species (317 of 342 on the diagonal of the best matching).

Replace each round centre by a Gaussian with its own covariance, a tilted ellipse of any size, and assign each bird by probability instead of distance. On the same standardized bills, a Gaussian mixture puts 330 of 342 birds (96%) with their species, against 317 for k-means, and the specimen goes back to the Chinstraps.

How a mixture is fitted (the EM algorithm, k-means with soft assignments) is the next lecture. The notation stays: \(x^{(i)}\), \(\mu_j\), \(k\); the assignment \(c^{(i)}\) becomes a probability.

• K-means alternates two steps, assign each point to its nearest centre and move each centre to its mean, and each step lowers the distortion J.
• Distances add up the columns, so standardize first: on the penguins that alone moves the match with the species from 75% to 93%.
• The elbow and the silhouette measure separation, not meaning: the penguins have a good answer at 2, 3 and 6.
• K-means finds round clusters of similar size, from wherever it starts: restart it, and expect it to miss tilted or small groups.