Unsupervised Learning: Clustering and Dimensionality Reduction

Most of the world's data has no labels. Customer logs, sensor streams, documents and images arrive without anyone saying which group they belong to or which features matter. Unsupervised learning finds structure in that data on its own: clustering groups similar points, and dimensionality reduction compresses many features into a few that keep the important variation. Practitioners use these tools for customer segmentation, anomaly screening, compressing features before a supervised model, deduplicating data, and, very often, simply looking at a high-dimensional dataset before deciding what to do with it.

Intuition: sorting a pile without instructions

Imagine someone dumps a thousand photos on a table and asks you to "organize them." Nobody tells you the categories. You might group them by who is in them, by location, or by colour. Each choice is a different notion of similarity, and each produces a different, equally valid grouping. That is the central fact of unsupervised learning: the algorithm finds the structure that its similarity measure and assumptions can see. There is no ground truth to score against, so judgment about the distance metric, the scaling of features, and the shape of clusters you expect matters more than in supervised learning.

Dimensionality reduction is the same pile viewed differently. If every photo is described by a million pixel values, but most of the variation comes from a few factors (lighting, pose, who is in the frame), you can describe each photo with a handful of numbers and lose little. PCA finds the directions of greatest variation; t-SNE and UMAP try to keep neighbours close in a 2D picture.

k-means clustering

The objective

Given $n$ points $x_1, \dots, x_n \in \mathbb{R}^d$ and a number of clusters $k$, k-means looks for cluster assignments $c_i \in \{1,\dots,k\}$ and centroids $\mu_1,\dots,\mu_k \in \mathbb{R}^d$ that minimize the within-cluster sum of squares (WCSS, called inertia in scikit-learn):

$J = \sum_{i=1}^{n} \lVert x_i - \mu_{c_i} \rVert^2 .$

Finding the global minimum is NP-hard in general, so we use Lloyd's algorithm, which alternates two steps:

1. Assignment step: give each point to its nearest centroid, $c_i = \arg\min_j \lVert x_i - \mu_j \rVert^2$.
2. Update step: move each centroid to the mean of its assigned points, $\mu_j = \frac{1}{|S_j|}\sum_{i \in S_j} x_i$, where $S_j$ is the set of points assigned to cluster $j$.

Each step can only lower $J$ or leave it unchanged: the assignment step picks the best centroid for each point given fixed centroids, and the mean is the point that minimizes squared distance to a set. Since there are finitely many assignments, the algorithm stops after a finite number of iterations, but only at a local minimum. That is why implementations run several random restarts (n_init) and keep the best, and why k-means++ initialization (pick each new initial centroid with probability proportional to its squared distance from the nearest centroid already chosen) is the standard default.

Worked iteration

Six points in 2D: $A=(1,1)$, $B=(2,1)$, $C=(1,2)$, $D=(5,4)$, $E=(6,5)$, $F=(5,6)$. Use $k=2$ with a deliberately poor start: $\mu_1=(1,1)$, $\mu_2=(2,1)$.

Iteration 1, assignment. Squared distances to $(\mu_1, \mu_2)$:

$J = 0 + 0 + 1 + 18 + 32 + 34 = 85$.

Iteration 1, update. $\mu_1 = \text{mean}(A, C) = (1, 1.5)$. $\mu_2 = \text{mean}(B, D, E, F) = \left(\tfrac{2+5+6+5}{4}, \tfrac{1+4+5+6}{4}\right) = (4.5, 4)$.

Iteration 2, assignment. For B: $\lVert B-\mu_1\rVert^2 = (2-1)^2 + (1-1.5)^2 = 1.25$, $\lVert B-\mu_2\rVert^2 = 2.5^2 + 3^2 = 15.25$, so B moves to cluster 1. A, C stay in cluster 1; D, E, F stay in cluster 2. With these assignments, evaluated against the current centroids, $J = 0.25+1.25+0.25+0.25+3.25+4.25 = 9.5$.

Iteration 2, update. $\mu_1 = \text{mean}(A,B,C) = (4/3, 4/3) \approx (1.33, 1.33)$, $\mu_2 = \text{mean}(D,E,F) = (16/3, 5) \approx (5.33, 5)$.

Iteration 3. No point changes cluster, so the algorithm has converged. Final $J = 4.0$. The objective fell $85 \to 9.5 \to 4.0$, exactly as the monotonicity argument promises.

Choosing k

There is no universally correct $k$; use several signals together.

• Elbow method. Plot $J$ against $k$. The best achievable $J$ never increases as $k$ grows (with $k=n$ distinct points it is zero), so look for the "elbow" where extra clusters stop buying much reduction. Elbows are often ambiguous.
• Silhouette score. For point $i$, let $a(i)$ be its mean distance to other points in its own cluster and $b(i)$ the mean distance to points in the nearest other cluster. Then $s(i) = \frac{b(i) - a(i)}{\max\{a(i), b(i)\}} \in [-1, 1].$ Values near 1 mean well-separated, near 0 mean on a boundary, negative means probably misassigned. The mean silhouette over all points is a standard model-selection score. In the example above, for A: $a(A) = \text{mean}(1, 1) = 1$, $b(A) = \text{mean}(\sqrt{25}, \sqrt{41}, \sqrt{41}) \approx 5.94$, so $s(A) \approx (5.94-1)/5.94 \approx 0.83$.
• Gap statistic, BIC (for mixtures), stability across resamples, and most importantly the use case: a marketing team that can run four campaigns needs four segments, whatever the elbow says.

import numpy as np
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

X, _ = make_blobs(n_samples=600, centers=4, cluster_std=0.8, random_state=0)
X = StandardScaler().fit_transform(X)        # k-means is distance-based: scale first

for k in range(2, 8):
    km = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    print(k, round(km.inertia_, 1), round(silhouette_score(X, km.labels_), 3))

The loop prints inertia and mean silhouette for each $k$; on this synthetic data the silhouette peaks at the true value of 4. Note StandardScaler: k-means uses Euclidean distance, so a feature measured in dollars will swamp one measured in years unless you scale.

What k-means assumes

k-means implicitly assumes clusters are roughly spherical, similar in size and similar in spread, because it assigns by distance to a single centre. It fails on elongated clusters, nested rings, and clusters of very different densities, and it must assign every point, including outliers, which drag centroids around.

Hierarchical clustering

Agglomerative clustering starts with every point as its own cluster and repeatedly merges the two closest clusters until one remains. The full merge history is a tree called a dendrogram; cutting it at a height gives a flat clustering, so you can inspect many values of $k$ from one run. The key choice is the linkage, the definition of distance between clusters:

Agglomerative clustering needs pairwise distances, so naive implementations cost $O(n^2)$ memory and at least $O(n^2)$ time; it is practical for thousands to tens of thousands of points, not millions. Its advantages are the dendrogram, the freedom to use any distance (for example edit distance between strings), and determinism.

DBSCAN: density-based clustering

DBSCAN (Ester, Kriegel, Sander and Xu, 1996) defines clusters as dense regions separated by sparse ones. It has two parameters: a radius $\varepsilon$ and a count min_samples. A core point has at least min_samples points within distance $\varepsilon$ (counting itself in scikit-learn). Clusters are formed by connecting core points that lie within $\varepsilon$ of one another, plus the non-core border points within reach of them. Everything else is labelled noise (label -1 in scikit-learn).

Strengths: no need to choose $k$, arbitrary cluster shapes (rings, crescents), and explicit outlier detection. Weaknesses: one global $\varepsilon$ cannot fit clusters of very different densities, and choosing $\varepsilon$ is fiddly. A common heuristic is the k-distance plot: sort every point's distance to its min_samples-th nearest neighbour and pick $\varepsilon$ near the knee. HDBSCAN, available as sklearn.cluster.HDBSCAN in recent scikit-learn versions, removes the single-$\varepsilon$ limitation by building a hierarchy over density levels and keeping the most stable clusters; it is usually the better first try for density-based clustering.

Gaussian mixture models

A Gaussian mixture model (GMM) is a probabilistic model that says each point was generated by first picking a component $j$ with probability $\pi_j$, then drawing from a Gaussian $\mathcal{N}(\mu_j, \Sigma_j)$:

$p(x) = \sum_{j=1}^{k} \pi_j \, \mathcal{N}(x \mid \mu_j, \Sigma_j), \qquad \sum_j \pi_j = 1 .$

It is fitted with the expectation-maximization (EM) algorithm. The E-step computes each point's responsibility $r_{ij} = \pi_j \mathcal{N}(x_i\mid\mu_j,\Sigma_j) / p(x_i)$, a soft assignment. The M-step re-estimates $\pi_j$, $\mu_j$, $\Sigma_j$ as responsibility-weighted averages. Like k-means it converges to a local optimum of its objective (here the log-likelihood). In fact k-means is the limiting case of a GMM with equal, spherical covariances shrinking to zero, so soft assignments become hard.

GMMs add three things k-means lacks: soft membership (a customer can be 70% one segment), elliptical clusters of different sizes and orientations via full covariance matrices, and a likelihood, which allows model selection with the Bayesian information criterion (BIC) and lets you flag low-density points as anomalies. The price is more parameters (a full covariance has $d(d+1)/2$ entries per component), so in high dimensions use covariance_type="diag" or reduce dimensionality first.

Dimensionality reduction with PCA

Derivation

Centre the data so each feature has mean zero, and stack it in a matrix $X \in \mathbb{R}^{n \times d}$. The sample covariance matrix is $S = \frac{1}{n} X^\top X$ (many texts use $\frac{1}{n-1}$; the directions found are identical). We want the unit vector $w$ such that projecting the data onto it, $z_i = w^\top x_i$, has maximum variance:

$\max_{w} \; \operatorname{Var}(z) = w^\top S w \quad \text{subject to} \quad w^\top w = 1 .$

Add a Lagrange multiplier $\lambda$: $\mathcal{L}(w,\lambda) = w^\top S w - \lambda (w^\top w - 1)$. Setting the gradient with respect to $w$ to zero gives $2Sw - 2\lambda w = 0$, that is

$S w = \lambda w .$

So $w$ must be an eigenvector of $S$, and the variance it captures is $w^\top S w = \lambda w^\top w = \lambda$. The best direction is the eigenvector with the largest eigenvalue; this is the first principal component. Repeating the argument with the extra constraint that each new direction is orthogonal to the earlier ones yields the remaining eigenvectors in decreasing order of eigenvalue. Because $S$ is symmetric and positive semidefinite, its eigenvectors can be chosen orthogonal and its eigenvalues are real and non-negative.

The SVD route. Write the centred data as $X = U \Sigma V^\top$ (singular value decomposition). Then $X^\top X = V \Sigma^2 V^\top$, so the columns of $V$ are the principal directions and the eigenvalues are $\lambda_j = \sigma_j^2 / n$. Libraries compute PCA through the SVD of $X$ rather than forming $X^\top X$, because it is more numerically stable; for very large data they use randomized SVD.

Explained variance

The fraction of total variance captured by component $j$ is

$\text{EVR}_j = \frac{\lambda_j}{\sum_{l=1}^{d} \lambda_l}.$

A common rule is to keep the smallest number of components whose cumulative ratio passes a target such as 90% or 95%, or to look for an elbow in the scree plot of $\lambda_j$.

Worked example

Four centred points: $(2,2), (-2,-2), (1,-1), (-1,1)$. Using $\frac{1}{n}$:

• $\operatorname{Var}(x_1) = (4+4+1+1)/4 = 2.5$, $\operatorname{Var}(x_2) = 2.5$,
• $\operatorname{Cov}(x_1,x_2) = (4 + 4 - 1 - 1)/4 = 1.5$.

So $S = \begin{pmatrix} 2.5 & 1.5 \\ 1.5 & 2.5 \end{pmatrix}$. The characteristic equation $(2.5-\lambda)^2 - 1.5^2 = 0$ gives $\lambda = 2.5 \pm 1.5$, so $\lambda_1 = 4$, $\lambda_2 = 1$. The eigenvectors are $w_1 = \frac{1}{\sqrt 2}(1,1)$ and $w_2 = \frac{1}{\sqrt 2}(1,-1)$. The first component explains $4/(4+1) = 80\%$ of the variance. Projecting $(2,2)$ onto $w_1$ gives $z = (2+2)/\sqrt 2 \approx 2.83$; projecting $(1,-1)$ gives $0$, because that point lies entirely along $w_2$.

import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA

X, y = load_digits(return_X_y=True)          # 1,797 images, 64 pixel features
pca = PCA().fit(X)
cum = np.cumsum(pca.explained_variance_ratio_)
k95 = int(np.argmax(cum >= 0.95)) + 1
print(f"components for 95% variance: {k95} of {X.shape[1]}")

Z = PCA(n_components=k95).fit_transform(X)   # compressed features
print(Z.shape)

On the digits dataset roughly 30 of 64 components keep 95% of the variance (the exact count prints when you run it). Two cautions: PCA is scale-sensitive, so standardize features measured in different units, and PCA is linear, so it cannot unroll curved structure.

t-SNE and UMAP: visualization, not measurement

t-SNE (van der Maaten and Hinton, 2008) converts distances in the original space into neighbour probabilities, using a Gaussian around each point whose width is set by the perplexity (roughly the effective number of neighbours, commonly 5 to 50). It then places points in 2D so that a heavy-tailed Student-t distribution over 2D distances matches those probabilities, minimizing the KL divergence between the two. The heavy tail lets moderately distant points spread apart, which produces the characteristic well-separated islands.

UMAP (McInnes, Healy and Melville, 2018) builds a weighted k-nearest-neighbour graph, interprets it as a fuzzy topological structure, and optimizes a low-dimensional layout with a cross-entropy objective. Its main knobs are n_neighbors (local versus global emphasis) and min_dist (how tightly points pack). It is usually faster than exact t-SNE, scales to millions of points, can transform new data, and tends to keep somewhat more global arrangement, although neither method guarantees global geometry.

Both are excellent for looking: spotting that labelled classes separate, finding mislabelled examples, seeing whether embeddings from a model group by topic. Both are dangerous for measuring:

• Cluster sizes mean little. t-SNE expands dense clusters and contracts sparse ones.
• Distances between clusters mean little. Two islands far apart on the plot may not be far apart in the data.
• Hyperparameters change the picture. Different perplexities or n_neighbors values can produce visibly different layouts from the same data, and random noise can form apparent clusters at low perplexity.
• Runs differ. Different seeds give different layouts. Always try several.
• Do not cluster on the 2D output and present the clusters as findings; cluster in the original or PCA space and use the plot to illustrate.

A standard recipe is PCA to about 50 dimensions first (removes noise and speeds things up), then t-SNE or UMAP to 2D.

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE

X, y = load_digits(return_X_y=True)
X50 = PCA(n_components=50, random_state=0).fit_transform(X)
for perp in (5, 30, 50):                       # compare several perplexities
    emb = TSNE(n_components=2, perplexity=perp, init="pca",
               random_state=0).fit_transform(X50)
    print(perp, emb.shape)                     # plot emb coloured by y

UMAP lives in the separate umap-learn package (import umap; umap.UMAP(n_neighbors=15, min_dist=0.1).fit_transform(X50)).

Evaluating clusters

Without labels, use internal metrics: silhouette, Davies-Bouldin (lower is better), Calinski-Harabasz (higher is better). They reward compact, separated clusters and therefore favour k-means-like shapes, so they can mislead for DBSCAN on non-convex data. When some labels exist, external metrics such as the adjusted Rand index and normalized mutual information compare a clustering with them. The most important check is qualitative: inspect samples from each cluster and ask whether a domain expert would find the grouping useful, and test stability by re-running on bootstrap samples.

Common pitfalls

• Not scaling features before k-means, hierarchical clustering, DBSCAN or PCA, so the feature with the biggest units dominates.
• Trusting one k-means run. Use multiple initializations (n_init) and a fixed random_state for reproducibility.
• Forcing k-means onto non-spherical data. Rings and crescents need DBSCAN, HDBSCAN, spectral clustering or a GMM.
• Treating the elbow as objective. It is a heuristic; combine it with silhouette, stability and business constraints.
• Reading t-SNE or UMAP plots as maps. Cluster sizes and inter-cluster gaps are not reliable quantities.
• Fitting PCA on all data before a train/test split. PCA is a learned transformation; fit it on training data only, ideally inside a pipeline.
• High-dimensional distances. In hundreds of dimensions, distances concentrate and nearest neighbours become less meaningful. Reduce dimensions or use a representation (such as learned embeddings with cosine similarity) built for the task.

In practice

• Start by looking. PCA to 2 components and a UMAP or t-SNE plot are cheap first diagnostics on any new dataset or set of embeddings.
• Segmentation: standardize, try k-means for a few $k$ values with silhouette scores, then check with a GMM if soft membership or elliptical segments make sense. Profile each segment with feature means so stakeholders can name it.
• Anomaly screening: DBSCAN/HDBSCAN noise labels and low GMM likelihoods are useful first-pass flags, covered further in the anomaly detection lesson.
• Embeddings: modern pipelines cluster vectors from neural encoders (text or image embeddings). Normalize them and use cosine distance (for k-means, L2-normalizing first makes Euclidean distance equivalent in ranking to cosine).
• Compression: PCA before a supervised model can remove noise and collinearity, but tree ensembles rarely benefit; test it with cross-validation rather than assuming.
• Scale: for millions of points use MiniBatchKMeans, randomized PCA, and approximate nearest-neighbour libraries.