Unsupervised Learning: Principal Components and k-Means Clustering

Every method so far in this series had an answer to learn from: a price, a class, a default. Often there is none. A table of crime statistics for 50 states has no “correct output”; the questions are about its structure. Which variables move together? Can four numbers per state be summarised by two? Are there groups of similar states? This is unsupervised learning, and its two classic tools are principal component analysis (PCA), which finds a few directions that capture most of the variation, and k-means clustering, which finds groups.

Without an outcome there is no test error to check against, so mistakes are easy and silent. This post works through both methods on a real data set (USArrests, crime rates in the 50 US states in 1973; McNeil, 1977) and on simulations, and shows what they find, what they assume, and where they fail.

1. Learning Without Labels

Unsupervised learning has only the predictors $X_{1}, \dots, X_{p}$ and looks for structure in them:

Question	Method	Result
Can the data be described with fewer variables?	principal component analysis	a few new variables (components)
Are there groups of similar observations?	clustering (k-means, hierarchical)	a group label per observation

It is often a first, exploratory step, and it is harder to judge than supervised learning: there is no response to check a principal component or a clustering against, so the result has to be judged by whether it makes sense and whether it is stable.

2. Principal Components

PCA (Pearson, 1901; Hotelling, 1933) replaces $p$ correlated variables with new ones, the principal components, that are uncorrelated and ordered by how much of the variation they capture. The first principal component of centred variables $X_{1}, \dots, X_{p}$ is the linear combination

Z_{1} = ϕ_{11} X_{1} + ϕ_{21} X_{2} + \dots + ϕ_{p 1} X_{p}, j = 1 \sum p ϕ_{j 1}^{2} = 1,

with the largest variance. The coefficients $ϕ_{j 1}$ are the loadings, the vector $ϕ_{1}$ the loading vector; without the constraint that it has length 1, the variance could be made arbitrarily large. The values $z_{i 1} = \sum_{j} ϕ_{j 1} x_{ij}$ are the scores, where each observation lands on the component. The second component has the largest variance among all combinations uncorrelated with the first (its loading vector is orthogonal to $ϕ_{1}$ ), and so on, up to $min (n - 1, p)$ components. They are the eigenvectors of the covariance (or correlation) matrix, in practice computed from the singular value decomposition of the data matrix.

Two views of the same direction

The first component has a second, equivalent definition: it is the line closest to the data, minimising the sum of squared perpendicular distances. For every direction, each point’s squared distance from the centre splits into its squared projection on the line plus its squared distance to the line (Pythagoras); the total is fixed, so maximising one minimises the other. For the murder and rape arrest rates of the 50 states (both per 100,000 residents, centred), scanning all directions in steps of 0.05°:

Criterion	Best direction
Largest variance of the projections	73.10°
Smallest sum of squared distances to the line	73.10°
First principal component (exact)	73.11°

3. USArrests

The data: for each of the 50 states, arrests per 100,000 residents for murder, assault and rape, and the share of the population living in urban areas. After standardising each variable:

Variable	PC1 loading	PC2 loading
Murder	0.536	−0.418
Assault	0.583	−0.188
UrbanPop	0.278	0.873
Rape	0.543	0.167
Variance explained	62.0%	24.7%

PC1 weights the three crime rates almost equally and urbanisation less: it is an overall crime axis (Florida 2.98 and Nevada 2.85 at one end, North Dakota −2.96 and Vermont −2.77 at the other). PC2 is mostly urbanisation (Hawaii 1.56). Two components keep 86.8% of the variance of four variables, and a scatter plot of the scores, with the loadings drawn as arrows (a biplot; Gabriel, 1971), shows both the states and the variables:

The signs are arbitrary. Multiplying a loading vector and its scores by −1 describes the same component; software may report either, and the published USArrests loadings appear with both signs. Only the relative signs within a component carry meaning.

Scaling matters. Without standardising, PCA is dominated by whichever variable has the largest variance in its own units. Assault arrests run into the hundreds (variance 6,945), murders into the tens (variance 19). Unscaled, PC1 is essentially Assault (loading 0.995) and explains 96.6% of the variance, a statement about units, not about crime. Variables in different units should be standardised before PCA; variables in the same units (as murder and rape above) may be left as they are.

How many components?

The proportion of variance explained (PVE) by component $m$ is its variance divided by the total, and the cumulative PVE shows how much is kept by the first $m$ components (USArrests: 62.0%, 86.8%, 95.7%, 100%). A scree plot of the PVE often has an “elbow” after which components add little. There is no objective rule: PCA is usually exploratory, and the number of components is a judgement. When the components are used as inputs to a supervised model (principal components regression), that model’s cross-validated error can choose the number.

4. Clustering

PCA and clustering both simplify data, in different directions: PCA looks for a few variables that summarise the observations, clustering for a few groups of observations that are similar to each other and different from other groups. A typical use is market segmentation: grouping customers by income, occupation, distance from a city and so on, then treating each segment differently. “Similar” needs a definition, usually a distance, and the choice of distance (and of the variables’ scales) decides the result.

k-means

k-means partitions the observations into a chosen number $K$ of clusters $C_{1}, \dots, C_{K}$ so that the within-cluster variation, summed over clusters, is as small as possible. With squared Euclidean distance, the within-cluster variation of cluster $k$ is

W (C_{k}) = \frac{1}{∣ C _{k} ∣} i, i^{'} \in C_{k} \sum j = 1 \sum p (x_{ij} - x_{i^{'} j})^{2} = 2 i \in C_{k} \sum j = 1 \sum p (x_{ij} - \overset{x}{ˉ}_{kj})^{2},

the sum of all pairwise squared distances divided by the cluster size, which equals twice the sum of squared distances to the cluster’s centroid $\overset{x}{ˉ}_{k}$ . (On the simulated data below, both sides came to 374.313.) So k-means minimises the within-cluster sum of squares (WCSS) around the centroids.

Checking all partitions is impossible. The number of ways to split $n$ observations into $K$ non-empty groups is a Stirling number of the second kind, close to $K^{n} / K!$ : 9,330 for 10 points in 3 clusters, and about $6.2 \times 1 0^{70}$ for 150 points in 3 clusters. Instead, Lloyd’s algorithm (Lloyd, 1982; MacQueen, 1967) alternates two steps:

Update: move each centroid to the mean of the points in its cluster (the point that minimises their squared distances).
Assign: move each point to the cluster of its nearest centroid.

Each step can only lower the WCSS or leave it unchanged, so the algorithm always stops, typically after a few iterations (3 to 12 on the data below), and the WCSS never increased in any of the runs. But it stops at a local minimum, which depends on the random start:

Data (simulated)	Random starts	Result
3 groups, K = 3	50	all within 0.1% of the best (187.06), 24% exactly at it; 94.7% of points in their true group
6 groups, K = 6	100	best 155.1, worst 492.6; 76% of starts within 0.1% of the best

There is no guarantee that most starts reach the global minimum: with more clusters, a quarter of the starts here ended in clearly worse solutions. The standard remedy is to run the algorithm from many random starts and keep the solution with the lowest WCSS, or to start from well-spread centroids (k-means++; Arthur & Vassilvitskii, 2007). Results also differ between runs because the cluster labels are arbitrary: the same partition can be numbered in $K!$ ways.

Choosing K

The WCSS always falls as $K$ grows (with $K = n$ it is zero), so it cannot choose $K$ by itself. The silhouette (Rousseeuw, 1987) compares, for each point, its average distance to its own cluster ( $a$ ) with that to the nearest other cluster ( $b$ ): $s = (b - a) / max (a, b)$ , close to 1 when the point sits well inside its cluster. On the three-group data:

K	Best WCSS	Average silhouette
2	405.6	0.376
3	187.1	0.497
4	147.4	0.420
6	107.5	0.375

The silhouette picked the true $K = 3$ . On real data the best $K$ is rarely this clear, and whether the clusters make sense in the application (or merely cut a continuum) matters at least as much as any index.

5. Where k-Means Struggles

k-means looks for compact, roughly spherical clusters of similar spread, in the units the data are given in. Each part of that sentence is an assumption (Jain, 2010):

Situation (simulated, K = 2, best of 20 starts)	Share of points in the right group	Why
A second feature in much larger units	0.500 (chance)	distances are dominated by it
The same, standardised first	0.890	both features count
Two interleaved half-moons	0.750	clusters must be convex
A wide cluster (170 points) and a tight one (30)	0.735	prefers clusters of similar spread

Outliers pull centroids toward them, and every point, however far, is assigned to some cluster: k-means makes hard, all-or-nothing assignments. Gaussian mixture models give each point a probability of belonging to each cluster and allow different shapes and sizes; hierarchical clustering builds a tree of nested clusters (a dendrogram) for all $K$ at once; density-based methods find clusters of arbitrary shape and leave outliers unassigned.

6. Using Clusters

Clusters are often used as inputs to a later model: a separate model per cluster, or the cluster as a predictor. Two details matter.

A cluster label is a nominal variable. The numbers 1, 2, 3 are names; any numbering of the same clusters is equally valid. Used as a single numeric predictor in a regression, the cluster label gave an R² anywhere from 0.136 to 0.813, depending only on which of the 6 possible numberings was used; as two dummy variables it gave 0.824, whatever the numbering. (Cluster centroids or distances to them are numeric and can be used as they are.)

The geometry of the centroids follows the data. If two variables are perfectly correlated, all points lie on a line, and so does every centroid, which is a mean of points on it: with correlation 1, the three centroids were collinear to machine precision (determinant $1 0^{- 15}$ ). With correlation 0 nothing forces them onto a line (determinant 3.1). Of the two statements “with correlation 1 the centroids lie on a line” and “with correlation 0 they lie on a line”, only the first is true.

Finally, clustering that uses the test data, or that is tuned while looking at the outcome, leaks information into a later supervised model, as in the screening example of Cross-Validation and the Bootstrap.

7. Common Claims, Checked

Claim	Verdict
”PC1 is the direction of largest variance.”	Yes, and equivalently the line closest to the data: both criteria gave 73.10°, PC1 73.11°.
”The loadings of a component are unique.”	Up to sign: multiplying a loading vector and its scores by −1 gives the same component, and software differs.
”PCA can be run on the raw data.”	Only if the variables share units. Unscaled, USArrests’ PC1 was Assault alone (loading 0.995, 96.6% of the variance); standardised, a crime axis with 62.0%.
”Two variables projected onto two components is no reduction.”	Not by itself; the reduction comes from dropping components that explain little (USArrests: two of four keep 86.8%).
”There are almost Kⁿ ways to partition n points into K clusters.”	About Kⁿ/K! (a Stirling number): 9,330 for n = 10, K = 3, not 59,049. Labels are interchangeable.
”Lloyd’s algorithm finds the best clustering.”	It never increases the WCSS and always stops, but at a local minimum: with K = 6, starts ended between 155 and 493.
”Most random starts reach the global minimum.”	Not guaranteed: 76% of starts got within 0.1% with K = 6. Use many starts or k-means++.
”Choose K where the WCSS is smallest.”	The WCSS always falls with K (405.6 → 107.5); use the silhouette, the elbow, and judgement. The silhouette picked K = 3 (0.497).
“k-means finds the groups in the data.”	Only compact, similarly spread groups in the given units: 0.500 accuracy with an unscaled feature (0.890 standardised), 0.750 on half-moons.
”A cluster ID can be used as an ordinal predictor.”	It is nominal: as one number, R² ranged from 0.136 to 0.813 with the numbering; as dummies, 0.824.
”If V1 and V2 have correlation 1, the centroids lie on a line; with correlation 0, too.”	Only the first: with correlation 1 every point and centroid lies on a line (det 10⁻¹⁵); with correlation 0 not (det 3.1).

Summary

Idea	In one line
Unsupervised learning	Structure in X without a response; hard to evaluate
PCA	Uncorrelated directions of maximal variance = lines closest to the data
Loadings and scores	Unit-length directions (sign arbitrary) and coordinates along them
Scaling	Standardise variables in different units, or the units decide
How many components	PVE and scree plot; a judgement, or CV for a downstream model
k-means	Minimise within-cluster sum of squares; Lloyd’s algorithm reaches a local minimum
Choosing K	WCSS always falls; silhouette, elbow, meaning
Limits	Scale, convex shapes, similar spreads, outliers; nominal cluster labels

Next: Neural Networks and Backpropagation: How a Network Learns: from one neuron to a trained network, backpropagation derived and checked by hand, and what decides whether training works.

The experiments: unsupervised.py (numpy 2.1, seeded), with the USArrests data in usarrests.py; its PCA matches the published results (standard deviations 1.575, 0.995, 0.597, 0.416; PC1 loadings 0.536, 0.583, 0.278, 0.543). The lab runs k-means in the browser; from the same start it reproduces the numpy result exactly (WCSS 187.155366, cluster sizes 45, 54, 51). The biplot and PVE values are the numpy results.

References

Principal components

Pearson, K. (1901). On lines and planes of closest fit to systems of points in space. Philosophical Magazine, 2(11).
Hotelling, H. (1933). Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24(6).
Gabriel, K. R. (1971). The biplot graphic display of matrices with application to principal component analysis. Biometrika, 58(3).
Jolliffe, I. T. (2002). Principal Component Analysis (2nd ed.). Springer.
Jolliffe, I. T., & Cadima, J. (2016). Principal component analysis: A review and recent developments. Philosophical Transactions of the Royal Society A, 374(2065).

Clustering

MacQueen, J. (1967). Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 1.
Lloyd, S. P. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2).
Rousseeuw, P. J. (1987). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20.
Arthur, D., & Vassilvitskii, S. (2007). k-means++: The advantages of careful seeding. Proceedings of the 18th ACM-SIAM Symposium on Discrete Algorithms.
Jain, A. K. (2010). Data clustering: 50 years beyond K-means. Pattern Recognition Letters, 31(8).

Data

McNeil, D. R. (1977). Interactive Data Analysis. Wiley. (The USArrests data, distributed with R.)

Textbooks

James, G., Witten, D., Hastie, T., & Tibshirani, R. (2021). An Introduction to Statistical Learning (2nd ed.). Springer.
Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.). Springer.

Notes

Explorer