3.3 Clustering: Unsupervised Classification

🖥️ Lecture slides — Session 13 (Wed Oct 28)

Clustering is a form of unsupervised classification: the algorithm discovers structure in the data from the features alone, without labels. The goal is to group observations into coherent subgroups.

Whether a group is coherent depends on the distance between data points. The distance metric quantifies how similar or dissimilar two data points are, and every clustering algorithm relies on one.

This tutorial does not cover all possible clustering methods. No single clustering method works best for all scenarios; the right choice depends strongly on the inherent data structure. A simple summary, with toy examples of 2D data structure, is available in the sklearn package.

This lecture focuses on fundamental concepts that are relevant to the geosciences: 1) the definition of distance, illustrated with the most popular clustering algorithm, 2) K-means clustering, and 3) agglomerative clustering.

1. Distance¶

Distance is the basic measure of dissimilarity between data points, and clustering results depend directly on the metric you choose. Distance also returns later in the course as the building block of loss and cost functions for training deep learning models.

There are various ways to estimate and quantify distance between data points. Some of the most commonly used distance metrics include:

• Euclidean Distance: This is the straight-line distance between two data points in a multidimensional space. It is often used when the data features have similar units or scales. $d(\mathbf{x},\mathbf{y}) = \sqrt{ \sum_{i=1}^N (x_i-y_i)^2 }$

• Manhattan Distance: Also known as the ‘L1’ distance, it measures the sum of the absolute differences between corresponding elements of two data points. It’s suitable when movement along axes is restricted, such as in grid-based data. $d(\mathbf{x},\mathbf{y}) = \sum_{i=1}^N |x_i-y_i|$

• Geodesic Distance: Geodesic distance measures the shortest path between two points on the surface of a sphere. It matters for Earth’s geography, including GPS navigation and geodetic measurements.

• Correlation-Based Distances: In geophysical and geospatial data analysis, correlation-based distance metrics like Pearson correlation or Spearman rank correlation are often used to assess relationships between variables.

• Covariance-Based Distances: These distances, which take into account the spatial covariance or variogram models, are prevalent in geostatistics and spatial analysis.

• Cosine Similarity: This metric calculates the cosine of the angle between two data vectors, providing a measure of their similarity, particularly in high-dimensional spaces. It’s frequently used for text or image data.

Scikit-learn contains the most commonly used metrics in the context of Classic ML under the package metrics.DistanceMetric. More details in this scikit-learn documentation.

Relation to PCA Clustering and PCA both simplify the data via a small number of summaries. But the differences are:

• PCA seeks to reduce the dimensionality of the data, to find a low-dimensional representation of the data that explains a good fraction of the data variance,
• Clustering seeks to find homogeneous groups within the observations.

In fact, it is common to combine both for complex and high dimensional data: 1) PCA, 2) clustering on the PCs.

There are two main methods using clustering, k-means clustering and hierarchical clustering.

The toolbox scikit-learn has a collection of clustering algorithms and detailed documentation with tutorials.

2. Tutorial set up¶

Import useful Python packages.

Data¶

Earlier editions of this chapter used a 146 MB SeaFlow flow-cytometry table; we retired it in 2026 in favor of lighter datasets.

Our first working dataset is the Old Faithful geyser record from Yellowstone. Each row pairs the duration of one eruption (current, in minutes) with the duration of the following eruption (next, in minutes). We download it with pooch, which caches the file locally.

K-means relies on Euclidean distance, so features with larger numeric ranges dominate the result. Here both columns happen to be in minutes with similar ranges, but that is a coincidence of this dataset. Standardizing each feature (zero mean, unit variance) is a habit to keep for any distance-based method — in the seismic exercise at the end of this notebook, the features span orders of magnitude and scaling is not optional. We standardize now and store the result in data_faithful, a 2-column numpy array.

3. K-means¶

K-means is an unsupervised clustering method. The main idea is to separate the data into K distinct clusters. We then have two problems to solve. First, we need to find the k centroids of the k clusters. Then, we need to affect each data point to the cluster which centroid is the closest to the data point.

The goal is to partition $n$ data points into $k$ clusters. Each observation is labeled to a cluster with the nearest mean.

K-means is iterative:

1. assume initial values for the mean of each of the $k$ clusters
2. Compute the distance for each observation to each of the $k$ means
3. Label each observation as belonging to the nearest means
4. Find the center of mass (mean) of each group of labeled points. These are new means to step 1.

In the following, we denote $n$ the number of data points, and $p$ the number of features for each data point.

K-means only works with the Euclidean distance metric. In fact, its core concept is to use the Euclidean distance to measure and minimize inertia or within cluster variation.

Let us define a function to initialize the centroid of the clusters. We choose random points within the range of values taken by the data.

To be able to affect each data point to the closest centroid, we need to define the distance between two data points. The most common distance is the Euclidean distance:

$d(x,y) = \sqrt{\sum_{i = 1}^p (x_i - y_i)^2}$

where $x$ and $y$ are two data observation points with $p$ variables.

We then define a function to compute the distance between each data point and each centroid.

We now define a function to affect each data point to the cluster which centroid is the closest to the point. We also define an objective function that will be minimized until we reach convergence.

Our objective is to minimize the sum of the square of the distance between each point and the closest centroid:

$obj = \sum_{j = 1}^k \sum_{i = 1}^{N_j} d(x^{(i)} , x^{(j)}) ^2$

where $x^{(i)}$ is the $i^{th}$ point in the cluster $j$, $x^{(j)}$ is the centroid of the cluster $j$, and $N_j$ is the number of points in the cluster $j$.

After all points are assigned to a cluster, compute the new location of the centroid. It is just the value of the mean of all the points affected to that cluster:

For $1 \leq j \leq k$, $x_p^{(j)} = \frac{1}{N_j} \sum_{i = 1}^{N_j} x_p^{(i)}$

We can now code the K-means algorithm by assembling all these functions. We stop the computation when the objective function no longer decreases.

The Old Faithful data shows two regimes — short eruptions followed by short eruptions, and long eruptions followed by long eruptions — so we set k = 2 and run our from-scratch k-means.

K-means with Scikit learn¶

We now use the scikit-learn implementation on the same scaled data.

Two defaults worth knowing: init='k-means++' is the default initialization scheme, and since recent scikit-learn versions n_init defaults to 'auto', which runs several k-means++ restarts and keeps the best. You no longer need to set either parameter by hand.

Practical tips for k-means¶

1. How to evaluate the success of the clustering? How well separated are the clusters? There exist many ways to evaluate the quality of the clusters. Scikit-learn has summarized and packaged these tools into the module sklearn.metrics. See documentation here. Usually, high scores are better. Here is a summary:
   • If the data has ground-truth labels (e.g., the seismic source types in the exercise below), we can use several metrics:
     • homogeneity (each cluster contains only members of a given class, metrics.homogeneity_score(clusterID,true_label)), completeness (all members of a given class are assigned the same cluster, metrics.completeness_score), V-measure (2 x homogeneity x completeness / (homogeneity+completeness), metrics.v_measure_score) and all three metrics.homogeneity_completeness_v_measure(clusterID,true_label).
     • Fowlkes-Mallows index, FMI, that uses TP (True Positive), FP (False Positive), FN (False Negative). Is 0.0 for random cluster assignment and 1.0 for perfect label assignments.
   • If the data does not have ground truth label, you can use:
     • silhouette coefficient using the module metrics.silhouette_score, a high score or coefficient is better. It quantifies how tight the data within the clusters are and how well separated the clusters are. More details on implementation and visualization of the silhouette score here.
2. The number of clusters $k$ is a tunable parameter. To find the optimal number of clusters, we discuss below a few strategies.
3. The optimization of k-means may have local minima. The result therefore may differ due to the initialization of the clusters. It is recommended to use:
   • repeat random initialization, repeat k-means, and use the best set of cluster (one that has the lowest final error)
   • choose K-means++ initialization scheme using the sklearn parameter init='k-means++'. The algorithm selects initial cluster centroids using sampling based on an empirical probability distribution of the points’ contribution to the overall inertia.
4. The data variance of some of the features may affect the results. It may be difficult to find clusters if some of the data features (axis) are much larger than other. The data may need to be pre-processed (centered and scaled) or pre-conditioned (e.g., PCA).

Silhouette analysis

The silhouette coefficient and silhouette score are metrics used to assess the quality of clustering in unsupervised learning.

1. Silhouette Coefficient: The silhouette coefficient for a data sample measures how similar it is to its own cluster (cohesion) compared to other clusters (separation). It is calculated for each data sample and ranges from -1 to 1. A high silhouette coefficient indicates that the object is well matched to its own cluster and poorly matched to neighboring clusters. Conversely, a low silhouette coefficient suggests that the object may be in the wrong cluster.

   The formula for the silhouette coefficient (s) for a single data point is given by:

   $s(i) = \frac{b(i) - a(i)}{\max\{a(i), b(i)\}}$ where:

   • $a(i)$ is the average distance from the $i$-th data point to the other data points in the same cluster.
   • $b(i)$ is the smallest average distance from the $i$-th data point to data points in a different cluster, minimized over clusters.

2. Silhouette Score: The silhouette score is the average silhouette coefficient across all data samples in the dataset. It provides a global measure of how well-separated clusters are. The silhouette score ranges from -1 to 1 as well, where a high score indicates good clustering and a low score suggests overlapping or misclassified clusters.

   The formula for the silhouette score is given by:

$S = \frac{\sum_{i=1}^{N} s(i)}{N}$ where:

• $N$ is the number of data points in the dataset.

Interpretation:

• A silhouette coefficient close to 1 indicates that the data point is well-matched to its own cluster and poorly matched to neighboring clusters, signifying a robust and distinct cluster.
• A silhouette coefficient close to -1 suggests that the data point is possibly misclassified, as it is better matched to a neighboring cluster.
• A silhouette coefficient around 0 indicates overlapping clusters.

For the silhouette score:

• A score close to 1 implies well-defined clusters.
• A score around 0 suggests overlapping clusters.
• A negative score indicates that the majority of data points may be assigned to the wrong clusters.

In short, silhouette analysis helps select the number of clusters and assess the overall quality of a clustering.

Choice of the number of clusters: the Elbow Method¶

The elbow method is designed to find the optimal number of clusters. It runs the clustering algorithm with an increasing number of clusters $k$, and measures the average distance between data points and the cluster centroids. There are two typical metrics in the Elbow method:

• Distortion: It is calculated as the average of the squared distances from the cluster centers of the respective clusters. Typically, the Euclidean distance metric is used.
• Inertia: It is the sum of squared distances of samples to their closest cluster center.

For each value of $k$, we compute the mean of the square of the distance between the data points and the centroid of the cluster to which they belong. We then plot this value as a function of $k$. Hopefully, it decreases and then reaches a plateau. The optimal number of clusters is the value at the kink — the “elbow” of the curve.

We illustrate the method on the standardized Old Faithful data.

Compute the value of E for k between 1 and 8 and plot it. Look for the value of $k$ where the curve bends.

The elbow method does not always work very well. For example, see what happens when the points get closer to each other. Because data_faithful is standardized and centered at zero, we shrink the data toward the origin.

Let us see what happens when we decrease the number of data.

Both stress tests show the limits of the elbow criterion: when the clusters get closer together, or when the data get sparse, the curve loses its clear kink and the choice of $k$ becomes ambiguous. Use the elbow method as a guide, not a rule, and cross-check with the silhouette score.

Repeat K-means¶

The result is very sensitive to the location of the initial centroids. Repeat the clustering N times and choose the clustering with the best objective function.

Repeat k-means 20 times with k = 2 and keep the best run.

K-means can be a slow process to calculate and there are avenues to speed this up. One solution is to do mini-batch K-means, which takes a subset of data at each iteration to construct the centroids.

4. Hierarchical Clustering¶

In K-means, we use the Euclidean distance and prescribe the number of clusters K.

In hierarchical clustering, we choose difference distance metrics, visualize the data structure, and then decide on the number of clusters. There are two approaches to building the hierarchy of clusters:

• Agglomerative: each point starts in each unique cluster. data is merged in pairs as on creates a hierarchy of clusters.
• Divisive: initially, all data is into 1 cluster. The data is recursively split into smaller and smaller clusters.

There are several types of linkages. sklearn has detailed documentation, mostly for agglomerative: The different linkages methods are:

• Ward minimizes the sum of squared differences within all clusters. It is a variance-minimizing approach and in this sense is similar to the k-means objective function but tackled with an agglomerative hierarchical approach.
• Maximum or complete linkage minimizes the maximum distance between observations of pairs of clusters.
• Average linkage minimizes the average of the distances between all observations of pairs of clusters.
• Single linkage minimizes the distance between the closest observations of pairs of clusters.

We first import relevant packages.

Here, we create a fake data set that has 2 clusters that intermingle in a few data points.

Combine these two data sets as one.

First we explore the dendrograms.

Now that we have explored the structure of the data and built some intuition on how many clusters and the distribution of the data samples in the cluster.

Next, we choose a distance threshold and assign each data point to a cluster ID.

In Sci-kit learn, the entire algorithm is incorporated in the function AgglomerativeClustering sklearn doc here. The function still needs either a threshold distance or a number of clusters to perform the clustering and assign a cluster ID to each data sample.

5. PCA before clustering¶

Let us generate synthetic data.

Let us now do k-means clustering with 3 clusters.

What happens if we apply PCA + normalization before the clustering?

Exercise: clustering volcanic and tectonic seismic sources¶

Volcanoes like Mt Hood, an ice-capped volcano in the Cascades, produce a mix of seismic sources: volcanic and tectonic earthquakes, surface events such as rockfalls and avalanches, and background noise. Analysts label these events by hand. Can clustering find the source types from waveform features alone?

We use a curated dataset of Pacific Northwest seismic events — earthquake, explosion, surface event, and noise, 1000 of each — described by physical waveform features (spectral shape, envelope statistics, kurtosis, band energies). The dataset is archived on Zenodo: DOI 10.5281/zenodo.14025693. It returns in notebook 3.5 for supervised classification.

Step 1: load the data. The loader below downloads and caches the four class files, concatenates them, and drops the one feature column with missing values.

Step 2: standardize the 61 features. The features span orders of magnitude, so scaling is mandatory before any distance-based method.

Step 3: PCA. Plot the cumulative explained variance and keep enough principal components to explain about 80% of the variance.

Step 4: k-means with 4 clusters on the principal components. We know there are four source types, so we set k = 4.

Step 5: compare the clusters to the true labels. This dataset has ground-truth labels, so we can score the clustering with three metrics:

• Homogeneity is 1 when each cluster contains members of a single class only.
• Completeness is 1 when all members of a class land in the same cluster.
• V-measure is the harmonic mean of the two. All three range from 0 (random assignment) to 1 (perfect match).

Step 6: cross-tabulate clusters against source types.

Without ever seeing a label, k-means recovers a good part of the structure: one cluster is dominated by earthquakes and another by noise. Look at your crosstab for the pair of sources that mixes the most — here, explosions and surface events largely share a cluster. That is physically plausible: both are shallow, near-surface sources (quarry blasts, rockfalls, avalanches), so their waveforms carry little of the deep high-frequency energy that separates earthquakes, and their envelope and spectral features overlap. The supervised classifiers in notebook 3.5 do much better on the same features, which is the point of labels.