Choosing k

Week 6, Problem 4 (live)

This page belongs to Problem 4 of Week 6. It runs Python in your browser, so there is nothing to install. The first visit downloads Python, numpy, scikit-learn, and matplotlib which can take a minute; wait until the loading message disappears.

This is new, and I’m trialing it. I’m not entirely sure this will work on all devices, heaven knows what this is going to look like on mobile. I heartily welcome feedback about this assignment, feel free to send me an email at . I have confirmed that this is properly running on your machine, so if you do feel funny and write an infinite loop, it’s your own browser on the hook ¬‿¬ .

Below, each grey box is a small editor. Press Run Code to run it, then change it and run it again as often as you like. Boxes marked Exercise contain blanks ______ for you to fill in, and running them also checks your answer. Theoretically (?) your edits are saved in this browser, and Start Over restores the original code.

Every box can use the helpers in kmeans_tools.py, which build the datasets and draw the plots. You do not need to read that file, that’s mainly to keep the boilerplate down, but you may. The numbered Questions are things to think about as you go.

NoteWhat is already imported

Every box on this page starts with the following already imported, so you do not need to import anything yourself.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

For reference, the documentation pages you are most likely to want are:

The cost curve

Recall the k-means cost: if point \(i\) is in cluster \(c(i)\) and cluster \(\ell\) has center \(\boldsymbol{\mu}_\ell\), the cost is \[ \sum_{i=1}^n\|\mathbf{x}_i-\boldsymbol{\mu}_{c(i)}\|_2^2 . \] In code, the data points are the rows of X, labels[i] is \(c(i)\), and row l of centers is \(\boldsymbol{\mu}_\ell\). Python counts from \(0\), so the clusters are numbered \(0,\ldots,k-1\) rather than \(1,\ldots,k\). NumPy lets us index an array by an array, so centers[labels] is the matrix whose \(i\)-th row is \(\boldsymbol{\mu}_{c(i)}\). After fitting, scikit-learn stores the labels and centers it found as model.labels_ and model.cluster_centers_.

X - centers[labels] is the matrix whose \(i\)-th row is \(\mathbf{x}_i-\boldsymbol{\mu}_{c(i)}\). Square every entry, then add them all up.

The elbow method computes the cost of the best clustering it can find for each \(k\) and plots it against \(k\). The cost always goes down as \(k\) grows (more clusters never hurt the cost), so we look for the \(k\) after which it stops going down by much. If that sounds like a heuristic, it absolutely is! The helper report(X, title) does this for \(k=1,\ldots,8\), running k-means++ with ten restarts for each \(k\). It also plots the silhouette score from the next section, and colors the data by the clustering with the best silhouette score.

The four groups of four_blobs sit at the corners of a square of side \(7\), and spread is the standard deviation of each group. For now, don’t worry about the Silhouette score, we’ll talk about that shortly.

Rhetorical Question 1. Where is the elbow when spread=0.5? Change spread to 1.5 and then to 2.5, and run the box again each time. What happens to the elbow? Explain what you see in terms of how tight the groups are and how far apart they are.

With spread=0.5 there is a sharp elbow at \(k=4\): the cost drops from about \(1300\) to about \(100\), then barely moves. At spread=1.5 the elbow is still at \(k=4\) but softer, and at spread=2.5 the curve bends so gradually that there is no clear elbow at all. What matters is the spread relative to the separation \(7\). When the groups are tight, each one is far from the others, so splitting off a true group removes a lot of cost and splitting within a group removes very little. As the spread grows the groups overlap, the gap between those two kinds of splits shrinks, and the elbow smooths out.

Silhouette scores

The elbow is judged by eye. The silhouette score is a number and therefore, a little easier to quantify. Fix a clustering and a point \(i\) in a cluster \(C\) with at least two points. Let \[ a(i)=\frac{1}{|C|-1}\sum_{j\in C,\ j\neq i}\|\mathbf{x}_i-\mathbf{x}_j\|_2 \] be the average distance from \(\mathbf{x}_i\) to the other points of its own cluster, and let \[ b(i)=\min_{C'\neq C}\ \frac{1}{|C'|}\sum_{j\in C'}\|\mathbf{x}_i-\mathbf{x}_j\|_2 \] be the average distance from \(\mathbf{x}_i\) to the points of the nearest other cluster. The silhouette of point \(i\) is \[ s(i)=\frac{b(i)-a(i)}{\max\{a(i),b(i)\}}\in[-1,1]. \] It is close to \(1\) when \(\mathbf{x}_i\) is much closer to its own cluster than to any other, close to \(0\) when it sits on a border, and negative when it looks like it belongs to another cluster. (A point alone in its cluster gets \(s(i)=0\).) The silhouette score of the clustering is the average of \(s(i)\) over all points, and the heuristic picks the \(k\) with the largest silhouette score. It needs at least two clusters, so \(k\geq2\).

dist[own] contains the distance from \(\mathbf{x}_i\) to itself, which is \(0\). It adds nothing to the sum, but it should not be counted when you divide. own.sum() is the number of points in the cluster.

Rhetorical Question 2. Rerun the four_blobs box above with spread=0.5 and with spread=2.5, and read the silhouette column of the printed table. Which \(k\) does the silhouette score pick in each case? Does it agree with the elbow, and how confident should you be in its choice when spread=2.5?

It picks \(k=4\) in both cases. With spread=0.5 the score at \(k=4\) is about \(0.87\), far above every other \(k\), so it agrees with the elbow and the choice is convincing. With spread=2.5 it still picks \(k=4\), but with a score of only about \(0.40\), against about \(0.37\) to \(0.38\) for most other \(k\). A difference that small is kind of tight and could flip with a bad random sample.

Four datasets

Now try both heuristics on data that is harder to read. pairs_of_blobs has four groups arranged as two far-apart pairs. uniform_square has 200 points spread uniformly over a square, so it has no clusters at all. two_rings is the pair of concentric circles from the lecture.

You can also color any dataset by a clustering of your choice, for example report(pairs_of_blobs(), "pairs of blobs", k_show=4).

Rhetorical Question 3. For pairs_of_blobs, where are the elbows, and which \(k\) does the silhouette score pick? Which \(k\) is “right”?

There are two elbows: a huge one at \(k=2\), where the cost falls from about \(20000\) to about \(500\), and a smaller one at \(k=4\), where it falls to about \(60\). The silhouette score picks \(k=2\) (about \(0.91\)). Both \(k = 2\) and \(k = 4\) is pretty defensible. Which one is “right” depends on what you want the clusters for, and no heuristic can decide that for you.

Rhetorical Question 4. For uniform_square, is there an elbow? Multiply each cost by its \(k\); what do you notice? Which \(k\) does the silhouette score pick, and are the scores large? k-means still splits the square into \(k\) clusters: does that tell you the data has clusters?

There is no elbow; the cost decreases smoothly. The products \(k\times\text{cost}\) are all roughly the same (between about \(30\) and \(40\)), so the cost behaves like \(C/k\). That is what happens when \(k\) centers tile a region evenly: each cell has area about \(1/k\), and the squared distance to its center is proportional to that area. The silhouette score picks a large \(k\) (about \(8\) here), but all scores sit near \(0.4\) and barely change with \(k\). Finally, no! One million times no! k-means always returns \(k\) clusters, whatever the data, so the existence of its output is no evidence that clusters exist!

Rhetorical Question 5. For two_rings, which \(k\) do the heuristics suggest, and do the clusters mean anything?

The cost has no clear elbow, and the silhouette score keeps creeping up and picks the largest \(k\) we try. The clusters are pie-slice wedges that cut across both rings, at least to me they mean nothing. Morally, the rings are two clusters, so the only candidate ought to be \(k=2\), but k-means isn’t strong enough on its own.

Units matter

In units_example, the first feature is a length measured in decimeters. Measuring it in centimeters instead multiplies the first column by \(10\) and changes nothing about the objects being measured.

The usual fix is to standardize each feature: subtract its mean and divide by its standard deviation, so every column has mean \(0\) and standard deviation \(1\).

X.mean(axis=0) is the vector of column means, and X.std(axis=0) is the vector of column standard deviations. Without axis=0, NumPy averages over the whole matrix.

Rhetorical Question 6. Which \(k\) do the heuristics choose for the data in decimeters, and which for the data in centimeters? Which answer is right?

In decimeters both heuristics choose \(k=4\), the four groups. In centimeters both choose \(k=2\): the first feature now ranges over about \(60\) units against about \(4\) for the second, so distances are dominated by the first feature and k-means splits left from right only. The decimeter answer is right, in the sense that the objects form four groups whatever units we measure in.

Seeding, in practice

Problem 1 showed that one seed in each group is enough for Lloyd’s algorithm to succeed, and k-means++ seeding tries to spread the seeds out so that this happens. Here we run Lloyd’s algorithm once from each kind of seed (n_init=1), \(100\) times each, on ten groups of \(30\) points, and count how often a single run reaches the best cost we know of.

costs <= 1.01 * best is an array of True and False, and the mean of such an array is the fraction of True entries.

Rhetorical Question 7. Report the two success rates. For ten groups of equal size, random seeding puts exactly one seed in each group with probability only \(10!/10^{10}\approx0.0004\) (you do not need to show this). Why is the success rate of random seeding so much higher than that? Using your success rate for random seeding, about how many restarts (n_init) would you need for at least one run to reach the best cost with probability \(0.99\)? Hint: If one run succeeds with probability \(p\), then \(r\) independent runs all fail with probability \((1-p)^r\).

The success rates should come out around \(0.08\) for random seeding and nearly \(1\) (about \(0.99\)) for k-means++. One seed per group is sufficient for success but not necessary. If a group gets two seeds and a neighboring group gets none, Lloyd’s algorithm can often push one of the two centers over into the empty group, so many bad seedings still repair themselves. For the restarts we need \((1-p)^r\leq0.01\), that is \(r\geq\ln(0.01)/\ln(1-p)\). With \(p=0.08\) this is \(r\geq4.61/0.083\approx56\), so about \(55\) to \(60\) restarts, against a single run of k-means++.