Cluster Analysis

You have measured two hundred olive oils, or a week of tablet spectra from a production line, or a drawer full of unknown mineral samples. Nobody has told you which samples belong together. Cluster analysis asks the data themselves: are there natural groups here, and which samples belong to each?

Grouping things by their similarities is as old as science, but the algorithmic discipline of clustering grew up in mid-twentieth-century biology. In 1963, the entomologist Robert R. Sokal and the microbiologist Peter H. A. Sneath published Principles of Numerical Taxonomy [1]. They proposed classifying organisms not by expert judgement or assumed evolutionary history, but by measuring many observable traits and letting an algorithm group organisms by their overall similarity. Many taxonomists objected, but the toolbox they built (distance matrices, linkage rules and dendrograms) proved useful far beyond biology, and it is still the backbone of hierarchical clustering.

The other thread came from engineering and statistics. In 1957, Stuart Lloyd at Bell Labs devised a method for quantising signals in pulse-code modulation: representing a continuous signal by a few discrete levels. His procedure alternated between assigning each value to its nearest level and moving each level to the mean of the values assigned to it. The report circulated inside Bell Labs and was only published in 1982 [2]. Closely related ideas were described by Hugo Steinhaus in 1956 [13], and in 1967 James MacQueen gave the method its name, k-means, and studied its properties [3]. The batch version used today, and in the demo below, is usually called Lloyd's algorithm; MacQueen's original moved a centroid after each sample was assigned. A few years earlier, in 1963, Joe H. Ward Jr. had published a hierarchical method that merges groups so as to keep them as compact as possible [4].

In analytical chemistry, clustering was taken up as soon as computers became part of the laboratory; by 1983 Massart and Kaufman had written a whole book on interpreting analytical data with cluster analysis [14]. Today it is a standard exploratory tool, usually applied after spectral preprocessing and often on PCA scores.

In many chemical problems there are no labels to start with, only measurements and a question: are there natural groups in these data? A few typical situations:

• Origin. Two hundred olive oils from five regions. Before building any classifier, you want to know whether the spectra separate naturally, or whether some regions overlap so much that classification will be hard.
• Process monitoring. A tablet spectrum every 30 seconds for a week, and something drifted. Clustering can show when a new group appeared, and whether it was a sudden jump or a gradual shift.
• Library matching. Five hundred unknown Raman spectra and a library of ten thousand references. Clustering the unknowns first means you only need to match one representative per group, and then check the members of a group individually only when the match is ambiguous.
• Unusual samples. If most samples fall into a few tight groups and three sit far from all of them, those three deserve a closer look.

In all these cases clustering is exploratory: it looks for structure, it does not test a hypothesis.

Before grouping samples, we need to say how similar two samples are. The choice of distance shapes the clusters you find [5],[6].

The ordinary straight-line distance. For two samples $\mathbf{x}$ and $\mathbf{y}$ measured on $p$ variables:

It treats all variables equally and works well when the variables are on comparable scales and the groups are roughly round.

Also called city-block or $L_1$ distance, it adds absolute differences, like walking along the streets of a grid:

Because the differences are not squared, a single very large difference counts less than in the Euclidean distance.

where $\mathbf{S}$ is the covariance matrix of the variables. It accounts for correlations and different variances, so a difference along a direction in which the data naturally vary a lot counts less. It needs a reliable estimate of $\mathbf{S}$ and its inverse, which is impossible with more variables than samples, so in spectroscopy it is used on PCA scores rather than on raw spectra.

In clustering there is a catch. If $\mathbf{S}$ is the covariance of all samples together, the separation between the groups itself inflates the variance along the direction that separates them, and the Mahalanobis distance then shrinks exactly the differences you are looking for. Ideally $\mathbf{S}$ would be the pooled within-group covariance (the covariance inside each group, averaged over the groups), but that is unknown until the groups are found. This is why the Mahalanobis distance is used less in clustering than in classification, where the groups are known in advance.

For spectra, the angle between two spectra seen as vectors is often useful:

It ignores overall intensity: two spectra with the same shape but different intensities, for example because of different path lengths, have an angle of zero. Use it when shape matters more than size. Its close relative, the correlation coefficient between two spectra, is the same idea applied after centring each spectrum.

Hierarchical clustering builds a tree of nested groups. The usual variant is agglomerative (bottom-up): every sample starts as its own cluster, and the two closest clusters are merged again and again until only one remains.

The result is drawn as a dendrogram, a tree whose branches join at the height (distance) at which the merge happened. Read it from the bottom: samples that join low down are very similar; samples that only join near the top are very different.

The key question in step 4 is how to measure the distance between two groups of samples. Each answer gives a different method.

Single linkage (nearest neighbour) uses the closest pair of members:

It can follow long, irregular groups, but it suffers from chaining: a few intermediate samples act as stepping stones and merge two distinct groups early.

Complete linkage (farthest neighbour) uses the most distant pair:

It produces compact groups and resists chaining, but it can split a large, elongated group.

Average linkage (UPGMA, for "unweighted pair group method with arithmetic mean", introduced by Sokal and Michener in 1958 [18]) uses the mean of all pairwise distances:

where $|A|$ is the number of samples in $A$. It is a compromise between single and complete linkage.

Ward's method [4] does not measure a distance directly. It merges the pair of clusters whose union causes the smallest increase in the total within-cluster sum of squares:

where ${\bar{\mathbf{x}}}_A$ and ${\bar{\mathbf{x}}}_B$ are the cluster means. It tends to produce compact groups of similar size, and with the Euclidean distance it is one of the most popular choices in chemometrics. Software usually draws Ward's dendrogram at the height $\sqrt{2\Delta }$ rather than at $\Delta$ itself, so that two single samples merge at their Euclidean distance; this is what SciPy's ward and R's ward.D2 report, and what the demo below shows.

A dendrogram does not say how many clusters there are; it shows every possible grouping at once. To get a single partition you cut the tree at a chosen height: every branch below the cut becomes a cluster. The panel below lets you do it, and compare linkages.

Where to cut can be guided by:

• Gaps in the tree. A long vertical stretch without merges means that groups only joined at a large distance; cutting there separates them.
• Domain knowledge. If you know there should be three regions of origin, look at the three-cluster cut.
• Validation measures, such as the silhouette, described below.

Hierarchical clustering builds the whole tree. k-means directly divides the samples into a chosen number $k$ of groups. It is fast, simple and one of the most widely used clustering algorithms.

Try it step by step:

k-means minimises the within-cluster sum of squares (WCSS), the total squared distance from every sample to its centroid:

The two steps take turns at lowering $W$. With the centroids fixed, assigning each sample to its nearest centroid is the best possible assignment. With the assignment fixed, the mean is the point that minimises the sum of squared distances to the members (a small least squares problem). Since $W$ never increases and there are only finitely many ways to assign the samples, the algorithm always stops, usually after a few dozen iterations or fewer.

But it stops at a local minimum, which depends on where it started. A poor start can leave two centroids sharing one natural group while two other groups are forced together, and the algorithm cannot escape. The standard remedy is to run k-means many times from different starts and keep the solution with the lowest $W$; k-means++ starts make bad solutions less likely.

Because it works with distances to centroids, k-means implicitly prefers groups that are:

• Roughly round (convex) rather than elongated or curved.
• Of similar spread. A large, diffuse group can be split while small, tight groups are merged.
• Not too different in size.

Centring, scaling and PCA often bring data closer to these conditions, which is one reason preprocessing matters so much for clustering.

Neither k-means nor an uncut dendrogram tells you the "right" number of groups. Several criteria help.

Run k-means for $k=1,2,\dots$ and plot the best $W$ against $k$. The best $W$ never increases as $k$ grows (with as many clusters as samples it is zero), but after the true number of groups the improvement becomes small. The "elbow" where the curve flattens suggests a reasonable $k$. Real data do not always have a sharp elbow, so treat it as a hint.

The silhouette of sample $i$ compares how close it is to its own cluster with how close it is to the nearest other cluster [9]:

where $a(i)$ is the average distance from sample $i$ to the other members of its cluster, and $b(i)$ is the average distance to the members of the nearest other cluster. For a cluster with a single sample, $s(i)$ is set to 0. The value lies between $-1$ and $+1$:

• Close to $+1$: the sample sits well inside its cluster, far from the others.
• Close to $0$: the sample is on the border between two clusters.
• Negative: the sample is, on average, closer to another cluster, and may be misassigned.

The average silhouette width over all samples summarises the whole partition. Comparing it for different $k$, or different cuts of a dendrogram, and choosing the largest is a common rule.

The gap statistic [8] compares the observed $\log W$ with what you would expect from data with no cluster structure, typically samples spread uniformly over the same range. The chosen $k$ is the smallest one for which the gap at $k$ is at least the gap at $k+1$ minus one standard error. It is more principled than the elbow, but needs a set of simulated reference data sets (10 to 50 are typical).

Finding clusters is easy; any algorithm will return them, even for random data. Finding meaningful clusters is harder.

Internal measures use only the data:

• Silhouette plots. Plot the silhouette of every sample, sorted by cluster. Good clusters show wide bars and few negative values.

• Within-cluster sum of squares, compared with what the same algorithm gives on structureless data (for example, data with each variable shuffled independently).

• Davies-Bouldin index [11]. For each cluster it finds the most similar other cluster, measuring similarity as the sum of their spreads divided by the distance between their centroids, and averages over clusters:

  $$DB=\frac{1}{k}\sum _{c=1}^k\underset{l\ne c}{\max }\frac{{\sigma }_c+{\sigma }_l}{d({\mathbf{\mu }}_c,{\mathbf{\mu }}_l)}$$

  where ${\sigma }_c$ is the average distance of the members of cluster $c$ to its centroid. Lower is better.

When true labels are known, even for some samples, you can compare the clusters with them:

• Adjusted Rand index (ARI) [12]: the agreement between two partitions, corrected so that random assignments score about 0 and identical partitions score 1.
• Normalised mutual information (NMI): how much knowing the cluster tells you about the true label, from 0 (nothing) to 1 (everything). Unlike the ARI, it is not corrected for chance, and it tends to grow with the number of clusters [19].

External validation is the strongest check, but it needs labels, which exploratory work often lacks.

Clustering is very sensitive to how the data are prepared, because everything depends on distances.

Apply the usual corrections first:

• Baseline correction (for example asymmetric least squares or polynomial fitting) to remove slowly varying backgrounds.
• Scatter correction (SNV or MSC) to remove multiplicative effects of particle size or path length.
• Smoothing (for example Savitzky-Golay) if the noise is strong.

Otherwise, baseline and intensity differences can dominate the distances and hide the chemical groups.

Mean centring focuses on differences between samples rather than on the common signal (it does not change Euclidean distances, but it matters for PCA and for the spectral angle). Autoscaling gives every variable the same weight. Without it, variables with large numerical ranges dominate. With spectra, where all variables share the same unit, autoscaling is usually a bad idea: it gives the noise of empty spectral regions the same weight as the bands. Mean-centre spectra; autoscale tables that mix different units (see centring and scaling).

With hundreds or thousands of variables, clustering in the full space has two problems:

1. The curse of dimensionality. In very high dimensions, distances between points become more and more alike, which makes groups harder to distinguish.
2. Noise. Many spectral variables carry mostly noise, which adds to every distance and dilutes the informative differences.

The usual solution is to compute a PCA first and cluster the scores on the first few components. This removes most of the noise and keeps the structure in a handful of dimensions. Typically, 2 to 10 components are enough. The order of operations is: correct the spectra, centre (and scale if needed), compress with PCA, then cluster.

• Classification of spectra. NIR, mid-IR or Raman spectra of unknown samples, clustered with Ward's method on PCA scores, often reveal groups that correspond to different compositions, origins or processing.
• Origin and variety. Wines, coffees, teas, oils and pharmaceutical raw materials are routinely grouped by origin or variety from their spectra or chemical profiles. A classic example is Forina's set of 572 Italian olive oils from nine growing areas, described by the proportions of eight fatty acids [17]. Clustering shows how naturally the groups separate before any supervised model is built.
• Authenticity and adulteration. Genuine and adulterated products can form separate groups; samples in between may indicate partial adulteration, or simply natural variability. Only reference samples can tell.
• Process monitoring. Spectra collected over time form one tight group while the process is stable. A new group signals a change. Clustering complements PCA-based control charts by characterising the different states.

The examples simulate spectra of three sample types, mean-centre them (all variables are absorbances, so autoscaling would only inflate the noise of the empty regions), compress them with PCA, and then apply Ward's hierarchical clustering and k-means.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.metrics import silhouette_score, adjusted_rand_score
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
import matplotlib.pyplot as plt

# --- Synthetic spectra: three groups with distinct profiles ---
rng = np.random.default_rng(42)
n_per_group, n_wavelengths = 40, 200
wavelengths = np.linspace(400, 2500, n_wavelengths)

def band(center, width2):
    return np.exp(-((wavelengths - center) ** 2) / width2)

profiles = [0.8 * band(800, 5000),
            0.7 * band(1400, 6000),
            0.5 * band(800, 5000) + 0.6 * band(1700, 4000)]
X = np.vstack([p + rng.normal(0, 0.05, (n_per_group, n_wavelengths)) for p in profiles])
true_labels = np.repeat([0, 1, 2], n_per_group)

# --- Preprocessing: mean centring and PCA (no autoscaling: all variables are absorbances) ---
X_c = X - X.mean(axis=0)
pca = PCA(n_components=5)
scores = pca.fit_transform(X_c)
print(f"Variance explained by 5 PCs: {pca.explained_variance_ratio_.sum():.1%}")   # about 85%

# --- Hierarchical clustering with Ward's method ---
Z = linkage(scores, method='ward')
hc_labels = fcluster(Z, t=3, criterion='maxclust')   # cut the tree into 3 clusters

plt.figure(figsize=(12, 4))
dendrogram(Z, no_labels=True)
plt.title("Dendrogram (Ward's method on PCA scores)")
plt.ylabel('Merge height')
plt.show()

# --- k-means: elbow and silhouette for k = 2..8 ---
ks = range(2, 9)
wcss, sil = [], []
for k in ks:
    km = KMeans(n_clusters=k, n_init=20, random_state=0).fit(scores)
    wcss.append(km.inertia_)
    sil.append(silhouette_score(scores, km.labels_))

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(list(ks), wcss, 'o-'); axes[0].set_xlabel('k'); axes[0].set_ylabel('WCSS')
axes[0].set_title('Elbow plot')
axes[1].plot(list(ks), sil, 'o-'); axes[1].set_xlabel('k'); axes[1].set_ylabel('Average silhouette')
axes[1].set_title('Silhouette')
plt.tight_layout(); plt.show()

# --- Final k-means with k = 3, compared with the true groups ---
km = KMeans(n_clusters=3, n_init=20, random_state=0).fit(scores)
print(f"Silhouette (k=3): {silhouette_score(scores, km.labels_):.3f}")
# Well-separated groups: both methods should give ARI = 1
print(f"ARI k-means vs truth: {adjusted_rand_score(true_labels, km.labels_):.3f}")
print(f"ARI Ward vs truth:    {adjusted_rand_score(true_labels, hc_labels):.3f}")

plt.figure(figsize=(6, 5))
plt.scatter(scores[:, 0], scores[:, 1], c=km.labels_, cmap='viridis', edgecolors='k', linewidths=0.3)
plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%})')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%})')
plt.title('k-means clusters on PCA scores')
plt.show()

% --- Synthetic spectra: three groups with distinct profiles ---
rng(42);
n_per_group = 40;
n_wl = 200;
wavelengths = linspace(400, 2500, n_wl);
band = @(c, w2) exp(-((wavelengths - c).^2) / w2);

base1 = 0.8 * band(800, 5000);
base2 = 0.7 * band(1400, 6000);
base3 = 0.5 * band(800, 5000) + 0.6 * band(1700, 4000);
X = [repmat(base1, n_per_group, 1); repmat(base2, n_per_group, 1); ...
     repmat(base3, n_per_group, 1)] + 0.05 * randn(3 * n_per_group, n_wl);
true_labels = repelem((1:3)', n_per_group);

% --- Preprocessing: mean centring and PCA (pca centres by default) ---
[~, scores, latent] = pca(X);
scores = scores(:, 1:5);
fprintf('Variance explained by 5 PCs: %.1f%%\n', 100 * sum(latent(1:5)) / sum(latent));

% --- Hierarchical clustering (Ward) and a 3-cluster cut ---
Z = linkage(scores, 'ward');
hc_labels = cluster(Z, 'maxclust', 3);
figure; dendrogram(Z, 0); title("Dendrogram (Ward's method on PCA scores)");

% --- k-means: elbow and silhouette ---
ks = 2:8; wcss = zeros(size(ks)); sil = zeros(size(ks));
for i = 1:numel(ks)
    [idx, ~, sumd] = kmeans(scores, ks(i), 'Replicates', 20);
    wcss(i) = sum(sumd);
    sil(i) = mean(silhouette(scores, idx));
end
figure;
subplot(1, 2, 1); plot(ks, wcss, 'o-'); xlabel('k'); ylabel('WCSS'); title('Elbow plot');
subplot(1, 2, 2); plot(ks, sil, 'o-'); xlabel('k'); ylabel('Average silhouette'); title('Silhouette');

% --- Final k-means with k = 3 ---
labels = kmeans(scores, 3, 'Replicates', 20);
figure;
gscatter(scores(:, 1), scores(:, 2), labels);
xlabel(sprintf('PC1 (%.1f%%)', 100 * latent(1) / sum(latent)));
ylabel(sprintf('PC2 (%.1f%%)', 100 * latent(2) / sum(latent)));
title('k-means clusters on PCA scores'); grid on;

library(cluster)

# --- Synthetic spectra: three groups with distinct profiles ---
set.seed(42)
n_per_group <- 40
n_wl <- 200
wavelengths <- seq(400, 2500, length.out = n_wl)
band <- function(center, width2) exp(-((wavelengths - center)^2) / width2)

profiles <- rbind(0.8 * band(800, 5000),
                  0.7 * band(1400, 6000),
                  0.5 * band(800, 5000) + 0.6 * band(1700, 4000))
X <- profiles[rep(1:3, each = n_per_group), ] +
     matrix(rnorm(3 * n_per_group * n_wl, sd = 0.05), nrow = 3 * n_per_group)
true_labels <- rep(1:3, each = n_per_group)

# --- Preprocessing: mean centring and PCA (no autoscaling: all variables are absorbances) ---
pca_result <- prcomp(X, center = TRUE, scale. = FALSE)
scores <- pca_result$x[, 1:5]
cat(sprintf("Variance explained by 5 PCs: %.1f%%\n",
            100 * sum(pca_result$sdev[1:5]^2) / sum(pca_result$sdev^2)))

# --- Hierarchical clustering (Ward) and a 3-cluster cut ---
hc <- hclust(dist(scores), method = "ward.D2")
plot(hc, labels = FALSE, main = "Dendrogram (Ward's method on PCA scores)",
     xlab = "", sub = "", ylab = "Merge height")
rect.hclust(hc, k = 3, border = c("#2563eb", "#d97706", "#059669"))
hc_labels <- cutree(hc, k = 3)

# --- k-means: elbow and silhouette ---
ks <- 2:8
fits <- lapply(ks, function(k) kmeans(scores, centers = k, nstart = 20))
wcss <- sapply(fits, function(f) f$tot.withinss)
sil <- sapply(fits, function(f) mean(silhouette(f$cluster, dist(scores))[, "sil_width"]))
par(mfrow = c(1, 2))
plot(ks, wcss, type = "b", pch = 19, xlab = "k", ylab = "WCSS", main = "Elbow plot")
plot(ks, sil, type = "b", pch = 19, xlab = "k", ylab = "Average silhouette", main = "Silhouette")

# --- Final k-means with k = 3 ---
km3 <- kmeans(scores, centers = 3, nstart = 20)
par(mfrow = c(1, 1))
plot(scores[, 1], scores[, 2], col = c("#2563eb", "#d97706", "#059669")[km3$cluster],
     pch = 19, xlab = "PC1", ylab = "PC2", main = "k-means clusters on PCA scores")
print(table(cluster = km3$cluster, truth = true_labels))

Hierarchical clustering and k-means are the workhorses, but a few other methods are worth knowing:

• k-medoids (PAM, "partitioning around medoids") [5] works like k-means, but each centre is a real sample (the medoid) rather than a mean. It accepts any distance, including the spectral angle or the correlation, and it is less sensitive to outliers.
• Fuzzy c-means [15] gives each sample a degree of membership in every cluster instead of a single label. It is widely used in chemometrics when groups overlap.
• Gaussian mixture models are a probabilistic version of k-means: each group is a normal distribution with its own covariance, so groups can be elliptical, and each sample gets a probability of belonging to each group.
• DBSCAN [16] (density-based spatial clustering of applications with noise) grows clusters from dense regions of the data. It finds groups of any shape, does not need $k$, and labels isolated samples as noise.

Clustering is a good tool for:

• Discovering natural groups in unlabelled data.
• Grouping samples by spectral similarity (origin, type, quality).
• Spotting off-specification batches or process changes.
• Reducing a large set of samples to a few representative groups.
• Exploring the data before building a supervised classifier.

It is not the right tool for:

• Predicting a continuous property. Use regression.
• Heavily overlapping classes. Hard clustering will cut them arbitrarily; consider soft or probabilistic methods such as Gaussian mixture models, or supervised methods such as PLS-DA if labels exist.
• Raw, very high-dimensional data. Preprocess and compress with PCA first.

Keep in mind that clustering always returns clusters, even for structureless data, that results depend strongly on preprocessing, distance and linkage, and that small changes in the data or the starting point can move samples between clusters. Validation is not optional. Computationally, the naive agglomerative algorithm needs time up to $n^3$, although modern implementations bring the common linkages (single, complete, average and Ward) down to about $n^2$ [7],[20]. The memory needed for the distance matrix, proportional to $n^2$, remains the limit for very large data sets; k-means scales much better.

Clustering groups samples using nothing but their similarities. Hierarchical methods build a whole tree of groupings and let you choose where to cut it; the linkage rule decides what "close" means for groups. k-means alternates between assigning samples to the nearest centre and moving the centres to the mean, which lowers the within-cluster sum of squares at every step but can get stuck, so it should be run from many starts. The number of clusters is a choice, guided by elbows, silhouettes and chemistry. And because any algorithm will find groups, the real work is checking that they mean something.

From here:

• Principal component analysis is the usual step before clustering, and its scores plot gives a first visual impression of the groups.
• SIMCA (soft independent modelling of class analogy) builds a PCA model for each known class and is used when labelled classes exist.
• Discriminant analysis (LDA, PLS-DA) uses known labels to find the directions that best separate classes.