Choosing the Cluster Count

You ran k-Means Clustering for k = 1 through 10 and got ten inertia values, each lower than the last. That is not a puzzle, it is arithmetic: adding a cluster can never make the within-cluster sum of squares larger, because the optimizer for k = 4 includes every partition k = 3 can produce, plus better ones. At the extreme, k = n gives every point its own centroid and inertia 0. Inertia is monotonically non-increasing in k, so "the lowest inertia wins" selects n clusters and no clustering at all.

The elbow heuristic works around this by looking for the k where extra clusters stop buying much. The intuition: while k is below the true number of groups, each new centroid lands in a genuinely separate region and inertia plunges; once k passes that point, extra centroids subdivide groups that were already coherent, and the curve flattens.

import numpy as np
from sklearn.cluster import KMeans

ks = range(1, 11)
inertias = [KMeans(n_clusters=k, n_init=10, random_state=42)
            .fit(X_scaled).inertia_ for k in ks]

# elbow by the 'farthest from the chord' rule
p0, p1 = np.array([1, inertias[0]]), np.array([10, inertias[-1]])
line = lambda t: p0 + t * (p1 - p0)
ts = np.linspace(0, 1, len(inertias))
dists = [np.linalg.norm(np.array([k, e]) - line(t))
         for k, e, t in zip(ks, inertias, ts)]
elbow_k = list(ks)[int(np.argmax(dists))]

The max-distance-to-chord rule formalizes "the bend": fit the straight line from the (k=1, k=10) endpoints and take the k farthest below it. It is a defensible programmatic answer, but be honest about its limits in an interview:

Two practical details. Run the curve on scaled features — inertia on unscaled data is dominated by the largest-magnitude column, so the bend reflects that column, not the data's structure. And because each k is a fresh set of local minima, use the same n_init and random_state protocol for every k, or the curve's wobble is initialization noise you are about to misread as structure.

The elbow gives you a shortlist, not a decision. The next two methods put numbers on the candidates.

Before you read an elbow, normalize the objective to the same dataset and preprocessing, and record the random seed and initialization strategy. A single dramatic bend can be a lucky local optimum. Fit several seeds for each candidate k and plot the median inertia with a spread band; a bend that appears across runs is more credible. Also compare the gain from k to k+1 as a percentage of the previous reduction. If the improvement is small relative to the cost of maintaining another segment, the extra cluster needs a business reason. The elbow is a diagnostic for diminishing returns, not a statistical test that reveals a uniquely correct answer.

Worked numbers make the read concrete. Suppose inertias for k = 1..5 come back as [1420, 980, 705, 688, 679] on the same scaled matrix. The drop from k=1 to k=2 removes 440 units, k=2 to k=3 removes 275, then the gains collapse: 17, then 9. The bend sits at k = 3, and the normalized gains tell the same story — 31% and 28% reductions before the elbow, then 2.4% and 1.3% after it. Reporting the curve as percent reductions per step, rather than raw inertia units, is also what makes the plot readable to stakeholders who have never seen the objective.

Inertia only looks inward — how tight each cluster is. The silhouette coefficient adds the missing question: how far is each point from the next cluster over?

For a point i in cluster C, define two quantities:

• a(i) — the average distance from i to every other point in C (its intra-cluster cohesion). For a point alone in its cluster, a(i) = 0.
• b(i) — the smallest average distance from i to the points of any other cluster (its separation from the nearest neighbor cluster).

Then the coefficient for that point is:

s(i) = (b(i) - a(i)) / max(a(i), b(i))

Work the arithmetic once. Suppose i has a(i) = 0.80 (it sits comfortably among its own points) and b(i) = 1.60 (the nearest other cluster is twice as far). Then s(i) = (1.60 - 0.80) / 1.60 = 0.50 — solidly assigned. If instead a(i) = 1.50 and b(i) = 1.60, s(i) = 0.06: i is almost equidistant between clusters and its membership is a coin flip the algorithm resolved arbitrarily. Negative s(i) means a(i) > b(i) — the point is closer to another cluster than to its own, a direct misassignment.

The dataset-level score is the mean of all s(i). sklearn splits it into two calls for good reasons:

from sklearn.metrics import silhouette_score, silhouette_samples

score = silhouette_score(X_scaled, labels)          # one number for the whole run
per_pt = silhouette_samples(X_scaled, labels)       # per point, for profiling

for k in range(2, 8):
    lab = KMeans(n_clusters=k, n_init=10, random_state=42).fit_predict(X_scaled)
    print(k, round(silhouette_score(X_scaled, lab), 3))

Three traps the interview panel may probe. First, silhouette is undefined at k = 1 (there is no other cluster to compare against), which is one reason k = 1 never appears in these curves. Second, it is defined on the same distance metric you clustered with — compute it on unscaled features after scaling for clustering and you are measuring a different geometry. Third, a huge cluster swallowing many small ones can still score decently if the big cluster is tight; the mean hides the tail, which is why silhouette_samples matters — plot the distribution per cluster, not just the average, and look for clusters whose median s(i) sits near zero.

Cost note: silhouette needs all pairwise distances, O(n^2) time and memory, so on large datasets score a stratified sample (silhouette_score(..., sample_size=10000)) rather than the full table.

Silhouette values also depend on the chosen distance and the representation. Euclidean distance on unscaled features can reward units rather than structure, while cosine distance may be more appropriate for normalized embeddings. Compute the score on a representative sample when the full pairwise calculation is too expensive, and keep the sampling procedure fixed across candidate values of k. Inspect both the global mean and per-cluster distributions: one large, well-separated group can conceal a tiny cluster with many negative silhouettes. A negative score means the observation is, on average, closer to another cluster; it does not prove the data point is erroneous.

The per-cluster view changes decisions. With silhouette_samples you can compute a mean per cluster instead of one global number: a global 0.51 might decompose into cluster means of 0.71, 0.68, 0.55, and 0.09. That 0.09 cluster is the finding — its points sit on a boundary, so either k is one too high, that cluster needs its own features, or it is genuinely transitional and should be reported as such. When you compare candidate values of k, track the minimum per-cluster mean alongside the global mean; two candidates with identical global scores can differ sharply in their weakest cluster, and shipping a segmentation with one boundary-soup cluster is how trust in the whole analysis erodes.

The elbow is a picture and the silhouette is an average — both inherit k-Means' own definition of "good". The gap statistic adds a control group: it asks how much structure your data has compared with featureless data.

The construction is a hypothesis test in disguise. Fit k-Means for each k and record log(W_k), where W_k is inertia at that k. Then generate B reference datasets (typically 10–50) by sampling uniformly within each feature's range (or, better, uniform along each PCA axis so uncorrelated features are respected). Cluster each reference set the same way, and compute:

gap(k) = mean over b of log(W_k on reference b) - log(W_k on your data)

A large gap means your data is far more clustered than unstructured noise at that k. The selection rule, from Tibshirani et al., handles the error bars: compute sd_k = standard deviation of log(W_k) across references, scaled by sqrt(1 + 1/B), and pick the smallest k such that gap(k) >= gap(k+1) - sd_{k+1}. The one-standard-error logic stops you from chasing noise at larger k.

B trades variance against cost: with B = 5 the sd estimate is noisy and the rule fires inconsistently across seeds, while B = 50 stabilizes the band but multiplies the k-Means fits (B datasets times every candidate k). Ten to twenty references is the practical default; go higher only when two candidate k values sit within the error bands of each other and the choice has real money riding on it.

from sklearn.cluster import KMeans
import numpy as np

def log_inertia(X, k):
    return np.log(KMeans(n_clusters=k, n_init=10).fit(X).inertia_)

B, ks = 20, range(2, 11)
logW = np.array([log_inertia(X, k) for k in ks])
refs = np.array([[log_inertia(np.random.uniform(X.min(0), X.max(0), X.shape), k)
                  for k in ks] for _ in range(B)])
gap = refs.mean(0) - logW
sd = refs.std(0) * np.sqrt(1 + 1 / B)
pick = next(k for i, k in enumerate(ks)
            if i == len(ks) - 1 or gap[i] >= gap[i + 1] - sd[i + 1])

Now the part no formula can do. Three methods, three opinions — elbow says 4, silhouette peaks at 3, gap's rule returns 5. This is normal, and here is the tiebreaker hierarchy interviewers want to hear:

Statistical evidence narrows k to a shortlist. Business interpretability decides. A marketing team can run four campaigns; a segment of 0.4% of users is not worth a campaign regardless of its silhouette; two clusters that differ only on a feature nobody can act on are not segments. Concretely: check that every candidate cluster has enough members to matter (say >5% of the data), that each is describable in one sentence from its feature profile, and that the differences between clusters map to actions. If k = 3 wins the stats and the stakeholders can name exactly three distinct treatments, that is your answer — and the profiling work that turns a k into stakeholder-ready segments is the subject of Cluster Evaluation and Profiling. If the stats are split, present the shortlist with evidence and let the decision owner break the tie; a k nobody acts on is worth less than a slightly worse k that ships.