Clustering is the task of grouping data points such that points within the same group are more similar to each other than to points in other groups — without using any predefined labels. It is one of the most widely used forms of unsupervised learning, applied to customer segmentation, anomaly detection, document organisation, image compression, and as a preprocessing step for supervised learning. The challenge: clustering is inherently subjective (the “best” number of clusters depends on the use case) and evaluation without ground truth labels requires indirect metrics. This guide covers the four main clustering families with their mathematical foundations, practical code, and when to use each.
Clustering connects to the dimensionality reduction techniques in our Dimensionality Reduction guide (PCA and UMAP are commonly applied before clustering to improve quality and speed). It is examined in our Machine Learning Interview Q&A. The evaluation metrics described here are a subset of the broader model evaluation framework in our Model Evaluation guide. For customer segmentation case studies that apply clustering, see our Data Science Case Study Interview guide. The probability distributions underlying Gaussian Mixture Models connect to our Probability Distributions guide.
K-Means — Fast, Scalable, Assumes Spherical Clusters
K-Means partitions n data points into k clusters by minimising the within-cluster sum of squared distances (inertia): minimize sum over clusters C_j of sum over x_i in C_j of ||x_i – mu_j||^2. The algorithm alternates between two steps until convergence: (1) Assign each point to the nearest centroid; (2) Update each centroid to the mean of its assigned points. K-Means is guaranteed to converge (the objective decreases each step) but not to a global optimum — multiple random initialisations (n_init=10 by default) mitigate local minima.
K-Means++ initialisation: Rather than random initialisation, K-Means++ selects the first centroid randomly, then selects each subsequent centroid with probability proportional to its squared distance from the nearest existing centroid. This produces better initialisations and converges faster to better solutions. Always prefer KMeans++ over random initialisation (it is sklearn’s default).
from sklearn.cluster import KMeans, MiniBatchKMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score, davies_bouldin_score
import numpy as np
import matplotlib.pyplot as plt
# Always standardise before clustering — distance-based methods
# are scale-sensitive
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# --- Choosing k: Elbow method + Silhouette score ---
inertias, silhouettes = [], []
K_range = range(2, 12)
for k in K_range:
km = KMeans(n_clusters=k, init='k-means++', n_init=10,
random_state=42, max_iter=300)
labels = km.fit_predict(X_scaled)
inertias.append(km.inertia_)
silhouettes.append(silhouette_score(X_scaled, labels, sample_size=5000))
# Plot elbow curve
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(K_range, inertias, 'bo-')
ax1.set_xlabel('Number of Clusters k'); ax1.set_ylabel('Inertia')
ax1.set_title('Elbow Method')
ax2.plot(K_range, silhouettes, 'rs-')
ax2.set_xlabel('Number of Clusters k')
ax2.set_ylabel('Silhouette Score (higher = better)')
ax2.set_title('Silhouette Analysis')
plt.tight_layout(); plt.show()
# Optimal k: elbow in inertia + peak in silhouette
optimal_k = K_range[np.argmax(silhouettes)]
km_final = KMeans(n_clusters=optimal_k, init='k-means++',
n_init=20, random_state=42)
labels = km_final.fit_predict(X_scaled)
# For very large datasets: MiniBatchKMeans (10-50x faster)
mbkm = MiniBatchKMeans(n_clusters=optimal_k, batch_size=1024,
n_init=10, random_state=42)
K-Means limitations: Assumes spherical, equal-size clusters (the Voronoi partitioning constraint). Sensitive to outliers — a single outlier can pull a centroid far from its natural position. Requires specifying k in advance. Fails on non-convex cluster shapes (rings, crescents). For these cases, use DBSCAN or Spectral Clustering.
DBSCAN — Density-Based, Discovers Arbitrary Shapes
DBSCAN (Density-Based Spatial Clustering of Applications with Noise) defines clusters as dense regions separated by low-density regions. Two parameters control the density definition: epsilon (the neighbourhood radius) and min_samples (minimum points within epsilon for a core point). Points are classified as: core points (have at least min_samples neighbours within epsilon), border points (within epsilon of a core point but fewer than min_samples neighbours), and noise points (neither — labelled -1). DBSCAN requires no k specification, automatically discovers the number of clusters, and natively handles noise/outliers.
| Property | K-Means | DBSCAN | Hierarchical | Gaussian Mixture |
|---|---|---|---|---|
| Number of clusters | Must specify k | Automatic | Choose by cutting dendrogram | Must specify k |
| Cluster shape | Spherical only | Arbitrary | Arbitrary | Elliptical |
| Outlier handling | Forced into nearest cluster | Explicit noise label (-1) | Forced into cluster | Low probability |
| Scales to large n | Very well (MiniBatch) | Moderate (O(n log n) with index) | Poorly (O(n^2)) | Moderate |
| Soft assignments | No (hard) | No (hard) | No (hard) | Yes — probabilities |
| Best for | Large data, known k, spherical | Arbitrary shapes, outlier detection | Hierarchical structure, small n | Overlapping clusters, density estimation |
from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors
# --- Choosing epsilon with k-distance graph ---
# Plot sorted distances to the k-th nearest neighbour
# The "elbow" in this plot is a good epsilon estimate
k = 5 # typically min_samples
nbrs = NearestNeighbors(n_neighbors=k).fit(X_scaled)
distances, _ = nbrs.kneighbors(X_scaled)
k_distances = np.sort(distances[:, -1])[::-1]
plt.plot(k_distances)
plt.xlabel('Points sorted by distance')
plt.ylabel('Distance to 5th nearest neighbour')
plt.title('k-distance graph — elbow = good epsilon')
plt.show()
# Fit DBSCAN
db = DBSCAN(eps=0.5, min_samples=5, metric='euclidean', n_jobs=-1)
labels = db.fit_predict(X_scaled)
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
n_noise = (labels == -1).sum()
print('Clusters found:', n_clusters)
print('Noise points:', n_noise, '(' + str(round(n_noise/len(labels)*100, 1)) + '%)')
# If too many noise points: decrease min_samples or increase eps
# If clusters merge incorrectly: increase min_samples or decrease eps
Gaussian Mixture Models — Soft, Probabilistic Clustering
A Gaussian Mixture Model (GMM) models the data as a weighted mixture of k multivariate Gaussian distributions. Each component has a mean (mu_k), covariance matrix (Sigma_k), and mixing weight (pi_k). Fitting uses the Expectation-Maximisation (EM) algorithm: the E-step computes the posterior probability that each point belongs to each Gaussian; the M-step updates the parameters to maximise expected log-likelihood. GMMs produce soft cluster assignments — every point gets a probability of belonging to each cluster — and can model elliptical clusters of varying shapes, sizes, and orientations (by choosing the covariance type).
from sklearn.mixture import GaussianMixture
from sklearn.model_selection import cross_val_score
# Select number of components using BIC (Bayesian Information Criterion)
# Lower BIC = better balance of fit and complexity
bic_scores = []
K_range = range(2, 12)
for k in K_range:
gmm = GaussianMixture(n_components=k, covariance_type='full',
n_init=5, random_state=42)
gmm.fit(X_scaled)
bic_scores.append(gmm.bic(X_scaled))
optimal_k = K_range[np.argmin(bic_scores)]
print('Optimal components by BIC:', optimal_k)
# Fit final model
gmm = GaussianMixture(n_components=optimal_k, covariance_type='full',
n_init=10, random_state=42)
gmm.fit(X_scaled)
# Soft assignments (probabilities)
probs = gmm.predict_proba(X_scaled) # shape (n_samples, k)
labels = gmm.predict(X_scaled) # hard assignment = argmax
# Anomaly detection: low-probability points are anomalies
log_probs = gmm.score_samples(X_scaled)
threshold = np.percentile(log_probs, 2) # bottom 2% = anomalies
anomalies = X_scaled[log_probs < threshold]
print('Anomalies detected:', len(anomalies))
For dimensionality reduction before clustering (especially for high-dimensional data), our Dimensionality Reduction guide covers PCA, UMAP, and t-SNE — UMAP in particular preserves cluster structure and is the standard preprocessing step before DBSCAN or K-Means on high-dimensional data. The statistical foundations of GMMs and the EM algorithm connect to the Bayesian thinking in our Bayesian Statistics guide. For customer segmentation case studies using these techniques, our Data Science Case Study Interview guide covers how to present clustering results to business stakeholders.



