Unsupervised Learning

Supervised learning asks the model to imitate a labelled target; unsupervised learning asks the model to find structure in the inputs alone, with no teacher signal to define what "right" looks like. This lesson covers the three workhorse algorithms of that regime — k-means clustering, DBSCAN, and principal component analysis — together with the most important diagnostic question the field still has not solved: when is the structure even real? Most of what follows reduces, in one way or another, to a geometric claim: points in high-dimensional space are assumed to live near some lower-dimensional manifold, and the algorithm is asked to recover that manifold. When the assumption is wrong, the algorithm will happily fabricate structure that has no business being there.

Learning Objectives

1. Formulate the k-means objective as a sum of squared distances to centroids, run the Lloyd assignment–update loop to a fixed point, and explain why random initialisation can land in a bad local minimum and why k-means++ initialisation provably mitigates that.
2. Select the number of clusters $k$ using the elbow method (within-cluster sum of squares plotted against $k$) and the silhouette score $s(i) = (b(i) - a(i)) / \max\{a(i), b(i)\}$ as a complement that does not require a clear elbow.
3. Apply DBSCAN by defining a point's $\varepsilon$-neighbourhood and the min_samples threshold, classify points as core, border, or noise, and describe why DBSCAN recovers non-convex cluster shapes and needs no preset $k$.
4. Derive principal component analysis as variance maximisation: centre the data, form the covariance matrix, decompose it into eigenvectors and eigenvalues, and project onto the top components using the explained variance ratio $\lambda_i / \sum_j \lambda_j$.
5. Distinguish dimensionality reduction (lossy projection that preserves geometry for downstream models) from compression (storage-saving reconstruction with explicit reconstruction error), and contrast PCA with t-SNE and UMAP as non-linear neighbours.
6. Diagnose when clustering is meaningless — no natural cluster structure, an imposed $k$ that the data does not support, the curse of dimensionality collapsing pairwise distances, and cluster artefacts driven by preprocessing rather than signal.

1. The Unsupervised Setting

In supervised learning the dataset is $\{(x_i, y_i)\}_{i=1}^n$ and the loss function involves both the inputs and the labels. In unsupervised learning only the inputs are observed: $\mathcal{D} = \{x_i\}_{i=1}^n$ with $x_i \in \mathbb{R}^d$, and the model must produce some representation useful enough to act on without ever being told what the right answer was.

Three task families dominate:

• Clustering partitions the data into groups whose points are more similar to each other than to points in other groups. k-means and DBSCAN are the workhorses.
• Dimensionality reduction projects the data into a lower-dimensional space that retains as much of the original structure as possible. PCA is the linear workhorse; t-SNE and UMAP are the non-linear neighbours.
• Density estimation and representation learning model the distribution $p(x)$ directly (Gaussian mixtures, autoencoders, contrastive methods).

All three reduce to the same question: how are the points arranged in $\mathbb{R}^d$? The answer depends entirely on whether the data actually has the geometric structure the algorithm is looking for. A clustering run on a perfectly uniform cloud will return clusters — every algorithm must return clusters when asked for $k$ clusters — but those clusters are artefacts of the algorithm, not discoveries about the data.

2. k-Means Clustering

k-means takes a dataset $X = \{x_1, \ldots, x_n\}$ and a target $k$, and returns $k$ centroids $\{\mu_1, \ldots, \mu_k\}$ plus an assignment $c(i) \in \{1, \ldots, k\}$ that puts every point in the cluster of its nearest centroid. The objective is the within-cluster sum of squares (WCSS), also called the inertia:

Minimising $J$ over both the centroids and the assignments is jointly non-convex, but it admits a clean coordinate-descent algorithm when you fix one side and optimise the other.

2.1 The assignment step

Given fixed centroids $\mu_1, \ldots, \mu_k$, the optimal assignment for each point is the nearest centroid:

This step is a nearest-centroid lookup; for each point the cost is $O(kd)$.

2.2 The update step

Given fixed assignments, the optimal centroid of cluster $j$ is the mean of its assigned points:

This is the step that gives the algorithm its name: every centroid becomes the k-means (mean) of its cluster. The closed-form optimality follows from setting the gradient of $J$ with respect to $\mu_j$ to zero.

The two steps together strictly decrease (or leave unchanged) $J$ at every iteration. $J$ is bounded below by zero, so the sequence must converge, and because there are only finitely many possible assignments the algorithm must reach a fixed point in a finite number of steps.

3. The Lloyd Algorithm and Convergence

The Lloyd algorithm alternates the assignment and update steps until neither changes the partition:

import numpy as np

def kmeans(X, k, max_iter=300, tol=1e-6, seed=0):
    """Plain Lloyd k-means with random init (for contrast with k-means++)."""
    rng = np.random.default_rng(seed)
    #Random init: pick k distinct rows of X as initial centroids.
    init_idx = rng.choice(X.shape[0], size=k, replace=False)
    centroids = X[init_idx].astype(float).copy()

    for _ in range(max_iter):
        #Assignment step: each point joins its nearest centroid.
        dists = np.linalg.norm(X[:, None, :] - centroids[None, :, :], axis=2)
        #`dists` has shape (n, k); argmin over axis=1 gives the cluster id.
        labels = np.argmin(dists, axis=1)

        #Update step: each centroid becomes the mean of its cluster.
        new_centroids = np.array([
            X[labels == j].mean(axis=0) if np.any(labels == j) else centroids[j]
            for j in range(k)
        ])

        #Check for convergence: centroids moved by less than `tol`.
        shift = np.linalg.norm(new_centroids - centroids, axis=1).max()
        centroids = new_centroids
        if shift < tol:
            break

    return centroids, labels

The empty-cluster guard (if np.any(labels == j)) handles the rare case where a centroid loses all of its points; in that case the previous centroid is kept, but a more aggressive implementation relaunches that centroid at the point furthest from any current centroid.

3.1 Local minima and seed sensitivity

The Lloyd algorithm is not guaranteed to find the global minimum of $J$. Different initialisations routinely produce different final inertia, often with very different cluster geometries. A common diagnostic is to run k-means ten or twenty times with different seeds and keep the run with the lowest inertia. k-means++ is the principled way to make the initialisation itself less of a lottery.

4. k-Means++ and Why Random Init Fails

Random initialisation samples $k$ points uniformly from the dataset and uses them as the starting centroids. Two failure modes are common:

1. Collapsed seeds. Several seeds land close together. The algorithm splits that region into multiple clusters and leaves another dense region without a centroid at all.
2. Outlier seeds. A seed lands in a sparse tail. The cluster it grows is a thin sliver collecting outliers; the bulk of the data is left to fight over the remaining $k-1$ centroids.

Both failure modes are unrecoverable: Lloyd's monotonic-decrease property prevents the algorithm from ever moving a centroid by more than a few units in a single step, so a bad seed tends to stay bad.

k-means++ replaces the random sample with a distance-proportional sample:

def kmeans_plus_plus_init(X, k, rng):
    """k-means++ seeding: each new centroid is sampled with probability
    proportional to its squared distance from the nearest already-chosen
    centroid."""
    n = X.shape[0]
    centroids = np.empty((k, X.shape[1]), dtype=float)

    #Pick the first centroid uniformly at random.
    centroids[0] = X[rng.integers(n)]

    for j in range(1, k):
        #Squared distance from each point to its nearest existing centroid.
        dists = np.min(
            np.linalg.norm(X[:, None, :] - centroids[None, :j, :], axis=2),
            axis=1,
        )
        probs = dists ** 2
        probs /= probs.sum()      #Normalise to a probability distribution.
        centroids[j] = X[rng.choice(n, p=probs)]

    return centroids

The first centroid is uniform; each subsequent centroid is sampled with probability proportional to $D(x)^2$, where $D(x)$ is the distance from $x$ to the nearest already-chosen centroid. Points that are far from any existing centroid are over-represented in the draw, so the seeds are spread out across the data instead of clumping.

The result is a theorem: Arthur and Vassilvitskii (2007) showed that k-means++ produces an initialisation whose expected $J$ is within a factor of $O(\log k)$ of the optimum. The constant in the big-O is small in practice, so k-means++ almost always matches or beats the best of many random restarts. scikit-learn's KMeans uses k-means++ by default; pass init="random" to force the bad old behaviour when teaching the lesson.

5. Choosing k: Elbow and Silhouette

Choosing $k$ is the part of k-means with no closed-form answer. Two diagnostics are the standard tools.

5.1 The elbow method

The within-cluster sum of squares $J$ is a monotonically non-increasing function of $k$. Adding more centroids can only lower $J$, and the decrease is steep while the new centroids are still capturing real structure and shallow afterwards. The "elbow" is the $k$ at which the marginal decrease sharply bends.

A worked sketch with synthetic numbers:

The big drop is at $k = 3$; the curve flattens past it. The elbow suggests three clusters, which in this case is the truth. The weakness is that real curves often do not have a clean elbow — the marginal decreases fade smoothly, and the elbow is in the eye of the beholder.

5.2 The silhouette score

The silhouette score gives a per-point answer to "how well does this point fit its cluster?" and an aggregate answer for the whole clustering. For point $i$, let

• $a(i)$ = mean distance from $i$ to the other points in $i$'s cluster.
• $b(i)$ = mean distance from $i$ to the points in the nearest other cluster (the cluster that minimises the mean distance to $i$).

Then the silhouette coefficient is

A score near $+1$ means the point is far from neighbouring clusters and tightly inside its own; a score near $0$ means the point sits on the boundary; a score near $-1$ means the point may belong in a different cluster. The aggregate silhouette score is the mean of $s(i)$ over all points, and the rule of thumb is to pick $k$ that maximises it.

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

scores = []
for k in range(2, 11):
    km = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    scores.append(silhouette_score(X, km.labels_))

best_k = 2 + int(np.argmax(scores))

The silhouette score is more honest than the elbow on data without a clean break in the inertia curve, because it directly measures cluster separation rather than curvature. It is also expensive: $O(n^2 d)$ per $k$ from the pairwise distance matrix.

6. DBSCAN: Density-Based Clustering

DBSCAN (Density-Based Spatial Clustering of Applications with Noise, Ester et al., 1996) does not ask for $k$. Instead it asks for two parameters: a neighbourhood radius $\varepsilon$ and a minimum sample count min_samples. The algorithm classifies every point as one of three kinds:

• Core point. At least min_samples other points (including the point itself in some conventions) lie within distance $\varepsilon$. A core point is the interior of a cluster.
• Border point. Fewer than min_samples points lie within $\varepsilon$, but it lies within $\varepsilon$ of some core point. A border point is the fringe.
• Noise point (also called outlier). Fewer than min_samples neighbours and not within $\varepsilon$ of any core point. Noise points are not assigned to any cluster.

Two core points in the same $\varepsilon$-neighbourhood belong to the same cluster; connectivity is propagated transitively. The result is a clustering in which the algorithm decides how many clusters exist, if any, based on density alone.

6.1 Worked example

Take $\varepsilon = 1.0$ and min_samples = 4. The figure: a dense crescent-shaped region has 80 points, a second dense round blob has 60, and 7 outliers sit far from either. DBSCAN will:

1. Identify the 80 + 60 = 140 points whose $\varepsilon$-neighbourhoods contain at least 4 points as core points.
2. Stitch the crescents' cores together via $\varepsilon$-connectivity into one cluster, and the blob's cores into a second.
3. Classify any point that touches a core point by $\varepsilon$ but is not itself a core as a border point of its cluster.
4. Leave the 7 far-away outliers as noise.

The crescent is not convex, which is the central capability: k-means would cut the crescent in half with a Voronoi boundary; DBSCAN follows the density. The catch is that DBSCAN needs the user to pick $\varepsilon$ and min_samples, and the two parameters are coupled — doubling $\varepsilon$ roughly compensates for halving min_samples — so a sensible heuristic is to set min_samples to $d+1$ (where $d$ is the data dimension) and use a $k$-distance plot to pick $\varepsilon$ at the knee.

6.2 k-means vs DBSCAN

DBSCAN's main failure mode is that it cannot cope with clusters of very different densities — the same $\varepsilon$ that finds a sparse cluster will dissolve a dense one into noise. HDBSCAN extends DBSCAN to varying densities by building a hierarchy of clusters and extracting the most stable ones.

7. Principal Component Analysis

PCA takes a centred data matrix $X \in \mathbb{R}^{n \times d}$ (rows are points, columns are features) and finds orthogonal directions $w_1, \ldots, w_d \in \mathbb{R}^d$ such that projecting onto the first $q < d$ directions captures as much of the data's variance as possible. The directions $w_j$ are called the principal components; they are the eigenvectors of the data's covariance matrix.

7.1 Centring

PCA is sensitive to translation. The first step is to subtract the column mean so each feature has zero empirical mean:

The centred matrix $\tilde{X}$ has column sums of zero. From here on "the data" means $\tilde{X}$.

7.2 Covariance matrix

The (empirical) covariance matrix is

$\Sigma$ is symmetric and positive semi-definite, so its eigenvalues are real and non-negative, and its eigenvectors are mutually orthogonal.

7.3 Eigen-decomposition

The spectral theorem gives $\Sigma = W \Lambda W^{\top}$ where $W = [w_1, \ldots, w_d]$ is orthogonal and $\Lambda = \text{diag}(\lambda_1, \ldots, \lambda_d)$ with $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_d \geq 0$. The eigenvectors $w_j$ are the principal components; the eigenvalues $\lambda_j$ are the variance of the data along $w_j$.

The explained variance ratio of component $j$ is

and the cumulative explained variance of the first $q$ components is $\sum_{j=1}^{q} \text{EVR}_j$. A common target is to keep enough components to explain 90–95% of the variance.

7.4 Projection

Projecting the data onto the first $q$ components gives the reduced representation

where $W_q = [w_1, \ldots, w_q]$. Each row $z_i$ is the low-dimensional code for point $x_i$. The reconstruction back to $\mathbb{R}^d$ is

and the reconstruction error is the sum of squared residuals over the discarded components:

7.5 Implementation

scikit-learn's PCA performs the eigen-decomposition under the hood and exposes the components and the explained variance ratio:

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

#Standardise so each feature has unit variance; PCA's scale-invariance
#is only over rotation, not over per-axis scaling.
X_std = StandardScaler().fit_transform(X)

pca = PCA(n_components=0.95)        #keep enough PCs for 95% variance
Z = pca.fit_transform(X_std)        #shape (n, q), q chosen automatically
print("components retained:", pca.n_components_)
print("explained variance ratio:", pca.explained_variance_ratio_)

n_components=0.95 is the most common production setting; n_components can also be an integer (keep the first $q$ components) or None (keep them all and inspect the spectrum).

8. PCA: Derivation from Variance Maximisation

The PCA objective can be derived two equivalent ways: variance maximisation (sequential) or reconstruction minimisation (joint). The variance-maximisation derivation is the most common.

8.1 First component

Find the unit vector $w_1$ that maximises the variance of the projection $\tilde{X} w_1$:

By the Rayleigh quotient, the maximum of $w^{\top} \Sigma w$ over $\|w\|=1$ is the largest eigenvalue $\lambda_1$, attained at the corresponding eigenvector $w_1$. This is the first principal component.

8.2 Subsequent components

For the second component, impose orthogonality to $w_1$ and maximise again:

The Lagrangian has solution $w = w_2$, the eigenvector of $\Sigma$ associated with $\lambda_2$, the second-largest eigenvalue. Continuing, $w_j$ is the eigenvector of $\Sigma$ associated with $\lambda_j$, the $j$-th largest eigenvalue.

8.3 Reconstruction view

The other derivation starts from "minimise the squared reconstruction error of $X$ from its rank-$q$ approximation". The optimal low-rank approximation in Frobenius norm is the projection onto the subspace spanned by the first $q$ eigenvectors of $\tilde{X}^{\top} \tilde{X}$, which is the Eckart–Young–Mirsky theorem. The two derivations arrive at the same $w_j$'s.

8.4 Singular value decomposition

In practice, the eigen-decomposition of $\Sigma$ is computed via the SVD of $\tilde{X}$:

where $V$'s columns are the principal components $w_j$ and $S^2 / n = \Lambda$. SVD is numerically more stable than forming $\tilde{X}^{\top} \tilde{X}$ explicitly when $d$ is large.

9. Dimensionality Reduction vs Compression

PCA is called both a dimensionality-reduction method and a compression method, but the two goals are different and the success criteria differ.

9.1 Dimensionality reduction

The goal is to produce a representation $Z \in \mathbb{R}^{n \times q}$ that is useful for downstream tasks — classification, regression, clustering, visualisation. The low-dimensional coordinates need not be recoverable into the original space; what matters is that distances and neighbourhood structure are preserved well enough for the downstream model to perform. A 2D PCA projection used to scatter-plot high-dimensional data is dimensionality reduction; no one will try to reconstruct the original points from the plot.

9.2 Compression

The goal is to store the data in fewer bits while permitting recovery that is approximately faithful. The success metric is the reconstruction error $\| X - \hat{X} \|_F^2$ or a perceptual equivalent. Compression is PCA plus a code book: project onto $W_q$, store the projection coefficients (or a quantised version of them), and reconstruct on demand.

A pipeline that trains a classifier on the PCA scores is using PCA for dimensionality reduction. A pipeline that stores the PCA scores, transmits them, and reconstructs at the receiver is using PCA for compression. The mathematics is the same; the engineering question is different.

10. PCA vs t-SNE vs UMAP

PCA is the linear baseline. Two non-linear methods dominate modern visualisation: t-SNE and UMAP. Both work by constructing a probability distribution over pairs of points in the original space and another in the embedded space, and minimising a divergence between them.

The PCA → t-SNE/UMAP pipeline is standard: PCA first to 50 dimensions to denoise and speed up the neighbour search, then t-SNE or UMAP to 2 or 3 dimensions for plotting. Both t-SNE and UMAP distort global distances: the apparent distance between clusters in a t-SNE plot is not a faithful report of the original distance. They are tools for showing which points are near each other, not for measuring how far.

from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt

Z50 = PCA(n_components=50).fit_transform(X_std)   #denoise + speed up
Z2  = TSNE(n_components=2, perplexity=30, random_state=0).fit_transform(Z50)

plt.scatter(Z2[:, 0], Z2[:, 1], c=y, s=4, cmap="tab10")
plt.title("t-SNE projection (preceded by PCA(50))")
plt.show()

A common mistake is to read structure in a t-SNE plot that is actually an artefact of the perplexity parameter: changing perplexity can split, merge, or move clusters without changing the underlying data. UMAP is less sensitive to its hyperparameters but not immune.

11. When Clustering Is Meaningless

The deepest pitfall in unsupervised learning is not algorithmic; it is epistemological. Most datasets do not have cluster structure, and most clustering algorithms will return clusters anyway because that is what they were asked to do. The four common failure modes:

11.1 No natural structure

If the data is sampled from a uniform or smoothly varying distribution, no partition into discrete clusters is correct. k-means will still return $k$ groups; the clusters are artefacts of the algorithm's geometry, not discoveries about the data. A silhouette score near zero across the board is a strong hint that the data does not cluster. Hierarchical clustering and density-based methods are not immune either; they will still return something.

11.2 Forced $k$

Choosing $k$ is a modelling decision, and the wrong $k$ yields wrong clusters no matter how good the algorithm. With $k = 2$ on a dataset with three real clusters, the algorithm will merge two of the real clusters into one and call it a cluster. The elbow method and the silhouette score mitigate this but cannot eliminate it; the only honest fix is to ask whether the problem needs discrete clusters at all.

11.3 Curse of dimensionality

In high dimensions, pairwise distances between random points concentrate around their mean. The contrast between "near" and "far" disappears, and cluster boundaries built on Euclidean distance become meaningless. PCA preprocessing, cosine distance, or a domain-specific metric can rescue some problems, but clustering in raw hundreds-of-dimensions space is usually hopeless.

11.4 Preprocessing artefacts

A common failure mode: standardise the features, run k-means, find beautiful clusters — and discover that the clusters are an artefact of the standardisation. Per-feature scaling weights features by their inverse variance, so a feature measured in noisy units (low variance) gets disproportionate weight and dominates the clustering. The same data without standardisation clusters completely differently. The fix is to justify the preprocessing choice independently of the clustering result and to check that the clusters are stable across reasonable preprocessing perturbations.

11.5 A diagnostic checklist

Before reporting a clustering result:

1. Is there an independent reason to expect clusters? If the labels of a downstream task are known to be the ground truth, validate the clustering against them with adjusted Rand index or mutual information; if not, you are clustering on faith.
2. Is the silhouette score meaningfully above zero across most points? A mean silhouette of 0.05 is "no structure found", not "weak structure".
3. Are the clusters stable across random restarts and across reasonable perturbations of the preprocessing? A clustering that rearranges itself when one feature is removed is not a clustering of the data; it is a clustering of the noise.
4. Do the clusters mean something? The strongest signal of real structure is that the clusters correspond to a downstream property that was not used in forming them. Two clusters of customers that differ in churn rate are evidence of structure; two clusters of customers that differ only in cluster id are not.

Key Takeaways

• Unsupervised learning recovers structure from $\{x_i\}$ alone: clustering partitions points into groups, dimensionality reduction finds lower-dimensional coordinates, and density estimation models $p(x)$ directly.
• k-means minimises the within-cluster sum of squares $J = \sum_i \| x_i - \mu_{c(i)} \|^2$ by alternating nearest-centroid assignment and mean-of-cluster updates; the algorithm decreases $J$ monotonically and converges in finite steps to a local minimum.
• k-means++ seeds centroids with probability proportional to squared distance from the nearest already-chosen centroid, provably yielding expected cost within $O(\log k)$ of the optimum and avoiding the collapse-and-outlier failures of pure random initialisation.
• The elbow method inspects inertia's bend as $k$ grows; the silhouette score $s(i) = (b(i) - a(i)) / \max\{a(i), b(i)\}$ measures per-point cluster fit and picks $k$ that maximises the mean.
• DBSCAN defines core, border, and noise points via $\varepsilon$ and min_samples, discovers the number of clusters from density, follows arbitrary shapes, and explicitly labels outliers — at the cost of having to pick the two coupled parameters.
• PCA finds orthogonal directions that maximise variance; equivalently, they are the eigenvectors of the centred covariance matrix $\Sigma = \frac{1}{n} \tilde{X}^{\top} \tilde{X}$, sorted by descending eigenvalue.
• Dimensionality reduction (keep enough structure for downstream tasks) and compression (encode-and-reconstruct with explicit error budget) use the same PCA mathematics but answer different engineering questions.
• t-SNE and UMAP are non-linear neighbours that preserve local structure for visualisation but distort global distances; both should be preceded by a PCA denoising step, and both should be read with their hyperparameters in mind.
• Clustering is meaningless when the data has no natural structure, when $k$ is forced on data that does not support it, when the curse of dimensionality collapses the distance metric, or when the clusters are artefacts of preprocessing rather than signal. Validate against an independent criterion before reporting a clustering.