DBSCAN, hierarchical clustering, and Gaussian mixtures

1. Intuition

DBSCAN (density-based). A point is core if at least min_samples points (including itself) lie within distance ε. Clusters grow by chaining core points that are within ε of each other, then adding border points that are near a core point. Everything else is noise. You never choose K, and clusters can be any shape.

Hierarchical (agglomerative). Start with every point as its own cluster and repeatedly merge the two closest clusters. The full merge history is a tree, the dendrogram. Cut it at a height to get a clustering; different heights give coarser or finer groups.

Gaussian mixture model (GMM). Assume the data comes from K overlapping bell-shaped (Gaussian) clouds, each with its own center, spread, and orientation. Each point gets a probability of belonging to each cloud: soft membership. A point between two clusters might be 60/40.

2. Visualization

Synthetic data. On two moons plus 15 random points, DBSCAN finds both moons and labels 10 points as noise; K-means can only draw a straight boundary.

Interactive: the K-means lab's "Rings" and "Stretched" presets show the shapes that motivate these methods. Open it.

3. The math

Symbols

Gaussian mixture

Expectation-maximization (EM) alternates: compute responsibilities $\gamma_{ik}$ (E-step), then re-estimate $\pi_k, \mu_k, \Sigma_k$ as responsibility-weighted averages (M-step). K-means is the limiting case with hard assignments and identical round clusters.

Linkage for hierarchical clustering

The distance between clusters $A$ and $B$ can be the closest pair (single linkage, finds chains), the farthest pair (complete), the average pair (average), or the increase in within-cluster variance from merging (Ward, which produces compact clusters similar to K-means).

Worked example: DBSCAN by hand

Points on a line: $1.0, 1.5, 2.0, 2.4, 6.0, 6.3, 6.8, 12.0$, with ε = 0.6 and min_samples = 3.

Result: cluster 0 = {1.0, 1.5, 2.0, 2.4}, cluster 1 = {6.0, 6.3, 6.8}, and 12.0 is noise. scikit-learn returns exactly these labels.

Worked example: soft membership

A two-component GMM fit to synthetic data (a tight cluster at 0 and a wider one at 4) gives the point $(2.2, 0)$ responsibilities $[0.240, 0.760]$. Although it is closer to the left center, the right cluster's larger spread makes the point more plausible under it. K-means would assign it to the left cluster with full confidence.

4. Implementation

Build it from scratch

DBSCAN and the EM loop for a Gaussian mixture each fit in about 20 lines of NumPy. The full script, which checks both against scikit-learn, is code/07-unsupervised/clustering.py.

DBSCAN, part 1: neighbors and core points. Compute all pairwise distances once, list each point's ε-neighbors (including itself, as in the worked example), and mark a point core when it has at least min_samples of them.

def dbscan(X, eps, min_samples):
    """Label each point with a cluster id, or -1 for noise. Also return the core-point mask."""
    D = np.sqrt(((X[:, None, :] - X[None, :, :]) ** 2).sum(axis=2))   # [n, n] pairwise distances
    nbrs = [np.flatnonzero(row <= eps) for row in D]                 # neighbors, including itself
    core = np.array([len(nb) >= min_samples for nb in nbrs])
    labels = np.full(len(X), -1)
    cluster = 0

Predict: on the line example with ε = 0.6, how many entries does nbrs hold for the point 12.0, and is it core?

DBSCAN, part 2: grow clusters. Start a new cluster at each unlabeled core point and flood-fill: every neighbor joins, but only core neighbors are pushed on the stack to expand further. Border points join a cluster and stop there; whatever is never reached keeps the label −1 (noise).

    for i in np.flatnonzero(core):
        if labels[i] != -1:                                          # already reached from another core point
            continue
        labels[i], stack = cluster, [i]
        while stack:                                                 # flood-fill through core points
            j = stack.pop()
            for q in nbrs[j]:
                if labels[q] == -1:
                    labels[q] = cluster                              # core or border point joins
                    if core[q]:
                        stack.append(q)                              # only core points keep expanding
        cluster += 1
    return labels, core

Predict: if you removed the if core[q] check and pushed every neighbor, what would happen to 1.0 and 2.4 in the worked example, and to the two moons when a few noise points sit between them?

Gaussian density. EM needs $\mathcal{N}(x \mid \mu, \Sigma)$ for every point. The einsum computes $d^\top \Sigma^{-1} d$ row by row (the squared Mahalanobis distance), so an elongated covariance makes distance along the long axis cheap.

def gauss_pdf(X, mean, cov):
    """Multivariate normal density N(x | mean, cov) for every row of X."""
    d = X - mean                                                     # [n, p]
    q = np.einsum("ij,jk,ik->i", d, np.linalg.inv(cov), d)           # squared Mahalanobis distance
    return np.exp(-0.5 * q) / np.sqrt((2 * np.pi) ** X.shape[1] * np.linalg.det(cov))

Predict: with cov the identity matrix, what does q reduce to?

EM loop. The E-step turns weighted densities into responsibilities $\gamma_{ik}$ (rows of R, each summing to 1). The M-step is K-means' update step with soft weights: $N_k$ effective points, new weights, weighted means, weighted covariances. The tiny ridge 1e-6 * I keeps a component from collapsing onto a single point.

def gmm_em(X, k, iters=200, seed=0, tol=1e-8):
    """Fit a full-covariance Gaussian mixture with EM. Returns weights, means, covariances, responsibilities."""
    n, p = X.shape
    r = np.random.default_rng(seed)
    pi, mu = np.full(k, 1 / k), X[r.choice(n, k, replace=False)]    # equal weights, random points as means
    S = np.array([np.cov(X.T) for _ in range(k)])                    # start every component at the data covariance
    prev = -np.inf
    for _ in range(iters):
        dens = np.stack([pi[j] * gauss_pdf(X, mu[j], S[j]) for j in range(k)], axis=1)   # [n, k]
        ll = np.log(dens.sum(axis=1)).sum()                          # log-likelihood never decreases
        R = dens / dens.sum(axis=1, keepdims=True)                   # E-step: responsibilities, rows sum to 1
        Nk = R.sum(axis=0)                                           # M-step: effective points per component
        pi, mu = Nk / n, (R.T @ X) / Nk[:, None]
        S = np.array([(R[:, j, None] * (X - mu[j])).T @ (X - mu[j]) / Nk[j] + 1e-6 * np.eye(p)
                      for j in range(k)])                            # small ridge stops collapse onto one point
        if ll - prev < tol:
            break
        prev = ll
    return pi, mu, S, R, ll

Predict: if you replaced each row of R with a one-hot vector at its largest entry fixed every S[j] to the identity, and kept every weight at 1/k, which algorithm would this loop become?

Hierarchical clustering has no from-scratch version here: the naive algorithm (recompute all cluster distances after every merge) is $O(n^3)$, and SciPy's linkage is the standard tool.

Use the library

In production, use scikit-learn for DBSCAN and Gaussian mixtures and SciPy for the dendrogram. Scale first, because every method here uses distances.

from scipy.cluster.hierarchy import dendrogram, fcluster, linkage
from sklearn.cluster import DBSCAN
from sklearn.mixture import GaussianMixture
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

X = StandardScaler().fit_transform(X_raw)

# DBSCAN: pick eps from the "knee" of sorted k-th neighbor distances
dists, _ = NearestNeighbors(n_neighbors=5).fit(X).kneighbors(X)
labels = DBSCAN(eps=0.3, min_samples=5).fit_predict(X)      # -1 marks noise

Z = linkage(X, method="ward")
labels_h = fcluster(Z, t=3, criterion="maxclust")            # cut into 3 clusters

gmm = GaussianMixture(n_components=3, covariance_type="full", random_state=0).fit(X)
probs = gmm.predict_proba(X)                                  # soft memberships
anomaly_score = -gmm.score_samples(X)                         # low likelihood = unusual

GaussianMixture stops when the average per-point log-likelihood improves by less than tol, which defaults to 1e-3. On the two-blob example that stops after 3 iterations, with $P(\text{left} \mid (2.2, 0)) = 0.402$ instead of the converged 0.240. Check gmm.converged_ and gmm.n_iter_, and lower tol when you read individual probabilities.

Run it

From guide/:

uv run code/07-unsupervised/clustering.py

Key lines of output:

DBSCAN eps=0.6 min_samples=3
  from scratch labels: [0, 0, 0, 0, 1, 1, 1, -1]  core: [1, 2, 5]
  scikit-learn labels: [0, 0, 0, 0, 1, 1, 1, -1]  core: [1, 2, 5]
GMM on two synthetic blobs (tight at x=0, wide at x=4)
  from scratch: weights [0.499 0.501], means x [0.003 3.879], log-likelihood -1524.10
  scikit-learn: weights [0.499 0.501], means x [0.004 3.88 ], log-likelihood -1524.10
  P(cluster | x=(2.2, 0)): from scratch [0.24 0.76], scikit-learn [0.24 0.76] (left, right)
moons: DBSCAN found 2 clusters and 10 noise points; from scratch agrees on every point: True

The DBSCAN labels match the table in the worked example: indices 1, 2, and 5 (the points 1.5, 2.0, 6.3) are core, and 12.0 is noise. EM from scratch reaches the same likelihood and the same soft membership as scikit-learn.

5. Engineering

Anomaly detection. Points DBSCAN labels as noise, points with low GMM likelihood, or points far from their K-means centroid are candidate anomalies. Isolation Forest is another common, scalable choice. Anomaly scores need a threshold, and that threshold should be set with whatever labeled incidents you have.

Choosing K for GMMs. Use BIC or AIC (they penalize extra components) together with domain review.

Watch out

Failure modes. Unscaled features changing every distance; one ε that can't fit clusters of different densities; reading dendrogram structure in random data; GMM components collapsing onto single points (use regularization or a different covariance type).

Common mistakes

• Running DBSCAN with default ε on unscaled data.
• Treating DBSCAN noise as errors to delete rather than signals to inspect.
• Using hierarchical clustering on 1 million points.
• Reading GMM probabilities from a fit that stopped early: check converged_ and lower tol.

6. Knowledge check

Practice exercise. With the same points and ε = 0.6, what changes if min_samples = 2?

Implementation challenge. Generate three blobs with very different densities. Show that one DBSCAN ε can't separate all three well, then try HDBSCAN (sklearn.cluster.HDBSCAN) and compare.

Summary

• DBSCAN grows clusters from dense core points, finds arbitrary shapes, and labels noise; tune ε and min_samples, and scale first.
• Hierarchical clustering builds a merge tree; cut the dendrogram to choose the number of clusters.
• Gaussian mixtures give soft memberships and likelihoods via EM; K-means is a hard, round special case.
• Low density or low likelihood points are natural anomaly candidates.

Next: PCA, t-SNE, and UMAP

Related: K-means · Model card: DBSCAN · Model card: Gaussian mixture

Interview angle