Hierarchical and Density-Based Clustering

Hierarchical clustering builds a sequence of nested groups, while density-based methods connect observations according to local neighborhood rules and can leave some unassigned. They offer alternatives to the compact nearest-center groups favored by k-means. They still require choices about distance, scale, and which groups to report. A returned cluster count reflects those choices as well as the data.

Linkage and the shape of a hierarchy

Agglomerative clustering starts with one cluster per observation and repeatedly merges two clusters. Linkage defines how to compare a pair of clusters; the pointwise distance, feature scaling, optional connectivity restrictions, and the stopping or cutting rule also affect the result. With two groups \(A=\{0,1\}\) and \(B=\{3,4\}\), the cross-group distances are 3, 4, 2, and 3. Single linkage uses 2, complete uses 4, and average uses 3.

from sklearn.cluster import AgglomerativeClustering, KMeans
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

X, y = make_moons(n_samples=1000, noise=0.06, random_state=0)
for link in ("ward", "complete", "average", "single"):
    lab = AgglomerativeClustering(n_clusters=2, linkage=link).fit_predict(X)
    print(f"{link:9s} ARI {adjusted_rand_score(y, lab):.4f}")
print(f"{'k-means':9s} ARI "
      f"{adjusted_rand_score(y, KMeans(2, n_init=10, random_state=0).fit_predict(X)):.4f}")
# ward      ARI 0.3452
# complete  ARI 0.3666
# average   ARI 0.4457
# single    ARI 1.0000
# k-means   ARI 0.2533

For these two crescents, single linkage agrees with the generating labels at ARI 1.0000. ARI compares pairwise cluster membership while accounting for chance agreement and does not require matching numeric cluster IDs. Short within-crescent links connect each crescent before the final cross-crescent merge. Ward gives 0.3452 and k-means 0.2533. Their preference for compact groups is poorly matched to this sample, but these two ARIs do not measure how closely the two fitted partitions agree with each other.

Ward chooses the merge with the smallest increase in within-cluster sum of squared distances: \(\Delta(A,B)=\frac{n_A n_B}{n_A+n_B}\|\bar x_A-\bar x_B\|^2\). The counts \(n_A,n_B\) and means \(\bar x_A,\bar x_B\) describe the two groups. This is the same sum-of-squares criterion used by k-means, but Ward makes irreversible greedy merges instead of reassigning points among centers. It need not reproduce a k-means partition. Ward requires Euclidean geometry in these implementations.

What single linkage costs

Single linkage’s short-hop connectivity can join groups through a thin bridge. The next experiment adds twelve observations between two Gaussian clouds and compares the result with single linkage before the bridge is added.

import numpy as np
from sklearn.cluster import AgglomerativeClustering
from sklearn.metrics import adjusted_rand_score

rng = np.random.default_rng(0)
a = rng.normal([0, 0], 0.5, size=(300, 2))
b = rng.normal([6, 0], 0.5, size=(300, 2))
bridge = np.column_stack([np.linspace(1.2, 4.8, 12),      # twelve points
                          rng.normal(0, 0.05, 12)])
X = np.vstack([a, b, bridge])
y = np.r_[np.zeros(300), np.ones(300), np.full(12, -1)]
real = y >= 0

base_labels = AgglomerativeClustering(n_clusters=2, linkage="single").fit_predict(X[real])
print(f"single without bridge ARI {adjusted_rand_score(y[real], base_labels):.4f}")

for link in ("single", "average", "complete", "ward"):
    lab = AgglomerativeClustering(n_clusters=2, linkage=link).fit_predict(X)
    print(f"{link:9s} ARI on the real clusters {adjusted_rand_score(y[real], lab[real]):.4f}"
          f"   sizes {np.bincount(lab)}")
# single without bridge ARI 1.0000
# single    ARI on the real clusters 0.0000   sizes [611   1]
# average   ARI on the real clusters 1.0000   sizes [306 306]
# complete  ARI on the real clusters 1.0000   sizes [311 301]
# ward      ARI on the real clusters 1.0000   sizes [310 302]

Single linkage has ARI 1.0 before the bridge. After adding it, the two-cluster cut has sizes 611 and 1, with displayed ARI 0.0000 on the 600 original observations. Average, complete, and Ward retain ARI 1.0 on those observations, but their full cluster sizes differ because they assign bridge points differently. This score excludes the bridge; it does not say that every assignment is unchanged.

This effect is called chaining. At a distance threshold, single-linkage clusters correspond to connected components of the graph joining pairs within that distance. A path can join distant regions even if most cross-region pairs are far apart. Whether one added observation is sufficient depends on the gap and threshold; this experiment uses twelve.

Complete linkage bases a merge on the farthest cross pair, so one short bridging pair alone cannot trigger a low-height merge of two extended groups. It favors bounded-diameter groups but does not mathematically prohibit non-convex arrangements of observations. Average linkage uses every cross pair. These criteria produce different compromises; no linkage is guaranteed to recover every connected shape or to resist every kind of noise.

The dendrogram

A fully built merge tree records all successive merges. Cutting at different heights gives nested partitions, so groups at a lower cut are contained in groups at a higher cut. Equal merge heights can make a horizontal cut skip some cluster counts. A library fit stopped early may not contain the full tree. The next block uses the bridge dataset from the preceding block; run the body examples in order.

import numpy as np
from scipy.cluster.hierarchy import linkage, fcluster

Z = linkage(X, method="ward")                  # the merge tree, computed once
print(f"{'cut height':>11} {'clusters':>9}")
for h in (2, 5, 10, 20, 40):
    print(f"{h:11.1f} {len(np.unique(fcluster(Z, h, criterion='distance'))):9d}")
print("last four merge heights:", np.round(Z[-4:, 2], 3))
#  cut height  clusters
#         2.0        25
#         5.0        11
#        10.0         3
#        20.0         2
#        40.0         2
# last four merge heights: [  8.807   9.155  12.027 103.541]

The last four Ward merge heights are about 8.807, 9.155, 12.027, and 103.541. After the merge near 12.027 and before the final merge near 103.541, there are two groups. This explains why cuts at both 20 and 40 give two. The five cut counts alone would not identify the location or width of that gap, which is why the code also prints the actual merge heights.

SciPy’s Ward height is \(\sqrt{2\Delta}\), not the raw sum-of-squares increase. For two singleton points at 0 and 2, the merge increases squared error by 2 and has height 2. Heights depend on scaling, linkage, and, for Ward, group sizes. A large gap can guide inspection of a cut, but is not an independent test that the resulting groups are natural populations.

Density-based clustering

For unweighted DBSCAN, an observation is a core point when its epsilon-radius neighborhood contains at least min_samples observations, including itself. Connected core points form groups. A non-core point within epsilon of a core point can join as a border point; other points receive label −1, or noise. Border points do not propagate connectivity the way core points do, and a border point adjacent to two groups can be assigned differently when input order changes.

For points 0, 0.1, 0.2, and 0.4 on a line, epsilon 0.15 and min_samples=3 make 0.1 a core point: its neighborhood contains the first three observations. Points 0 and 0.2 are border points, and 0.4 is noise. Thus low local neighbor count does not automatically mean noise. The returned number of groups, including possibly zero, follows these rules and parameter settings.

import numpy as np
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

X, y = make_moons(n_samples=1000, noise=0.06, random_state=0)
for eps in (0.05, 0.10, 0.15, 0.20, 0.30, 0.50):
    lab = DBSCAN(eps=eps, min_samples=5).fit_predict(X)
    n_clusters = len(set(lab)) - (1 if -1 in lab else 0)
    print(f"eps {eps:4.2f}  clusters {n_clusters:3d}  noise {np.mean(lab == -1):.3f}"
          f"  ARI {adjusted_rand_score(y, lab):.4f}")
# eps 0.05  clusters  13  noise 0.116  ARI 0.1683
# eps 0.10  clusters   2  noise 0.002  ARI 0.9960
# eps 0.15  clusters   2  noise 0.000  ARI 1.0000
# eps 0.20  clusters   2  noise 0.000  ARI 1.0000
# eps 0.30  clusters   1  noise 0.000  ARI 0.0000
# eps 0.50  clusters   1  noise 0.000  ARI 0.0000

At epsilon 0.15 and 0.20, DBSCAN matches both crescents with no noise labels, although no parameter directly requests two groups. At 0.05 it returns thirteen clusters and rejects 11.6% of observations; at 0.30 it returns one cluster. The ARI is computed on every row and treats all −1 labels as one predicted group. That is an evaluation convention, not a claim that rejected points form a coherent cluster.

Of the radii tested, 0.10 through 0.20 perform well on these labels. The sweep does not locate every successful radius. Epsilon is measured in the units of the chosen distance, and min_samples changes which points can connect groups. A sorted neighbor-distance curve is one way to propose candidate radii; the second exercise shows why it still needs investigation.

DBSCAN provides a rejection category as part of its clustering rule. Its clusters partition the assigned subset, while noise rows remain outside those clusters. Treating all noise as one additional label yields a partition for bookkeeping, but does not make it a meaningful group. Downstream use should specify how rejected observations are handled and report their fraction.

Different density scales

import numpy as np
from sklearn.cluster import DBSCAN, HDBSCAN
from sklearn.metrics import adjusted_rand_score

rng = np.random.default_rng(0)
X = np.vstack([rng.normal([0, 0], 0.12, size=(700, 2)),       # tight
               rng.normal([2.0, 0], 0.45, size=(700, 2)),     # medium
               rng.normal([6.0, 0], 1.30, size=(700, 2))])    # diffuse
y = np.r_[np.zeros(700), np.ones(700), 2 * np.ones(700)]

print(f"{'eps':>6} {'clusters':>9} {'noise':>7} {'ARI':>8}")
for eps in (0.15, 0.25, 0.40, 0.70, 1.00, 1.50):
    lab = DBSCAN(eps=eps, min_samples=10).fit_predict(X)
    n = len(set(lab)) - (1 if -1 in lab else 0)
    print(f"{eps:6.2f} {n:9d} {np.mean(lab == -1):7.3f} {adjusted_rand_score(y, lab):8.4f}")
for mcs in (20, 50):
    lab = HDBSCAN(min_cluster_size=mcs, min_samples=10, copy=True).fit_predict(X)
    n = len(set(lab)) - (1 if -1 in lab else 0)
    print(f"HDBSCAN mcs={mcs:3d}  clusters {n:2d}  noise {np.mean(lab == -1):.3f}"
          f"  ARI {adjusted_rand_score(y, lab):.4f}")
    print("  retained by generating group:", np.round([np.mean(lab[y == j] >= 0) for j in range(3)], 3))
#    eps  clusters   noise      ARI
#   0.15         4   0.375   0.8399
#   0.25         6   0.177   0.8161
#   0.40         4   0.050   0.9092
#   0.70         1   0.007   0.0002
#   1.00         1   0.001   0.0000
#   1.50         1   0.000   0.0000
# HDBSCAN mcs= 20  clusters  3  noise 0.071  ARI 0.9036
#   retained by generating group: [1.    0.993 0.794]
# HDBSCAN mcs= 50  clusters  3  noise 0.071  ARI 0.9036
#   retained by generating group: [1.    0.993 0.794]

These generating groups have standard deviations 0.12, 0.45, and 1.30. The six tested DBSCAN radii yield 4, 6, 4, 1, 1, and 1 groups. Small radii reject or fragment sparse regions, while larger radii connect groups through intervening observations. This sweep does not prove that no radius can return three groups. Different densities can make a global radius difficult to choose, but sufficiently separated unequal-density groups can still be recovered by DBSCAN.

HDBSCAN forms a hierarchy using mutual-reachability distance: for the default scaling, it takes the maximum of the two points’ core distances and their ordinary distance. A core distance is the neighbor radius associated with min_samples. The hierarchy is condensed using min_cluster_size, and the default excess-of-mass selection favors stable groups using membership and persistence across density levels. It is not simply choosing the widest interval on the epsilon axis. Both parameters matter. The code fixes min_samples=10 while changing the minimum cluster size; leaving it unset in scikit-learn would change it with that size.

Both HDBSCAN settings return three groups, reject 7.1% of rows, and give ARI 0.9036. Retention is about 100%, 99.3%, and 79.4% in the tight, medium, and diffuse generating groups. The count is right while many diffuse-group observations are excluded. DBSCAN at epsilon 0.40 has four groups and slightly higher all-row ARI, 0.9092. Count, pairwise agreement, and retention answer different questions. Report them together; reporting ARI only on retained rows would also need the retained fraction and an explanation of which rows were excluded.

Exercises

1. Cost and memory. Time repeated agglomerative and DBSCAN fits as sample size grows. Separately calculate square and condensed pairwise-distance storage, without treating either calculation as measured peak memory.

Which algorithm and representation are actually being timed? Which large arrays are only hypothetical calculations?

import time
import numpy as np
from sklearn.cluster import AgglomerativeClustering, DBSCAN
from threadpoolctl import threadpool_limits

rng = np.random.default_rng(0)
with threadpool_limits(limits=1):
    for n in (1000, 2000, 4000, 8000):
        X = rng.normal(size=(n, 10))
        times = {}
        for name, model in (
            ("agglomerative", AgglomerativeClustering(n_clusters=5)),
            ("DBSCAN", DBSCAN(eps=2.0, min_samples=5))):
            model.fit(X)
            t0 = time.perf_counter()
            for _ in range(3):
                model.fit(X)
            times[name] = (time.perf_counter() - t0) * 1000 / 3
        print(f"n={n:5d} agglomerative {times['agglomerative']:8.1f} ms"
              f" DBSCAN {times['DBSCAN']:8.1f} ms")

for n in (1000, 2000, 4000, 8000):
    print(f"n={n:5d} square {n*n*8/2**20:7.1f} MiB"
          f" condensed {n*(n-1)//2*8/2**20:7.1f} MiB")
for n in (100_000, 1_000_000):
    print(f"n={n:9,d} square {n*n*8/2**30:10,.1f} GiB"
          f" condensed {n*(n-1)//2*8/2**30:10,.1f} GiB")
# n= 1000 agglomerative     12.8 ms DBSCAN     25.1 ms # varies by machine
# n= 2000 agglomerative     64.3 ms DBSCAN     94.9 ms # varies by machine
# n= 4000 agglomerative    371.1 ms DBSCAN    359.5 ms # varies by machine
# n= 8000 agglomerative   1957.3 ms DBSCAN   1363.7 ms # varies by machine
# n= 1000 square     7.6 MiB condensed     3.8 MiB
# n= 2000 square    30.5 MiB condensed    15.3 MiB
# n= 4000 square   122.1 MiB condensed    61.0 MiB
# n= 8000 square   488.3 MiB condensed   244.1 MiB
# n=  100,000 square       74.5 GiB condensed       37.3 GiB
# n=1,000,000 square    7,450.6 GiB condensed    3,725.3 GiB

2. A candidate radius from neighbor distances. Compute a chord-based knee on the moons example and compare DBSCAN results on both sides of the suggested radius.

Does the candidate preserve the two crescents? Check rejection rates as well as the number of groups.

import numpy as np
from sklearn.neighbors import NearestNeighbors
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score

X, y = make_moons(n_samples=1000, noise=0.06, random_state=0)
k = 5

d, _ = NearestNeighbors(n_neighbors=k).fit(X).kneighbors(X)
kdist = np.sort(d[:, -1])                      # query itself is included among these k neighbors
print("k-distance percentiles:",
      {p: round(float(np.percentile(kdist, p)), 4) for p in (50, 75, 90, 95, 99)})

i = np.arange(len(kdist))                      # largest vertical gap below the endpoint chord
line = kdist[0] + (kdist[-1] - kdist[0]) * i / (len(kdist) - 1)
knee = kdist[np.argmax(line - kdist)]
print(f"knee suggests eps = {knee:.4f}")

for eps in (knee * 0.5, knee * 0.75, knee, knee * 1.5, knee * 2.0):
    lab = DBSCAN(eps=eps, min_samples=k).fit_predict(X)
    n = len(set(lab)) - (1 if -1 in lab else 0)
    print(f"  eps {eps:.4f}  clusters {n:2d}  noise {np.mean(lab == -1):.3f}"
          f"  ARI {adjusted_rand_score(y, lab):.4f}")
# k-distance percentiles: {50: 0.0413, 75: 0.0513, 90: 0.0674, 95: 0.0812, 99: 0.1063}
# knee suggests eps = 0.0608
#   eps 0.0304  clusters 45  noise 0.657  ARI 0.0073
#   eps 0.0456  clusters 29  noise 0.171  ARI 0.0923
#   eps 0.0608  clusters  4  noise 0.053  ARI 0.5734
#   eps 0.0913  clusters  2  noise 0.006  ARI 0.9881
#   eps 0.1217  clusters  2  noise 0.000  ARI 1.0000

3. What internal scores and stability measure. Compare silhouette, subsampling stability, and a uniform-reference gap calculation on four Gaussian blobs and a uniform square.

Can a repeatable geometric partition arise without separated generating groups? What reference distribution does the gap comparison use?

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import adjusted_rand_score, silhouette_score

rng = np.random.default_rng(0)

def stability(X, k, reps=25, frac=0.8):
    labels, idxs = [], []
    for _ in range(reps):
        idx = rng.choice(len(X), size=int(frac * len(X)), replace=False)
        labels.append(KMeans(k, n_init=10, random_state=0).fit_predict(X[idx]))
        idxs.append(idx)
    scores = []
    for i in range(reps):
        for j in range(i + 1, reps):
            common = np.intersect1d(idxs[i], idxs[j])
            ai, aj = dict(zip(idxs[i], labels[i])), dict(zip(idxs[j], labels[j]))
            scores.append(adjusted_rand_score([ai for c in common],
                                              [aj for c in common]))
    return float(np.mean(scores))

def gap(X, k, B=20):
    lo, hi = X.min(0), X.max(0)
    obs = np.log(KMeans(k, n_init=10, random_state=0).fit(X).inertia_)
    ref = [np.log(KMeans(k, n_init=10, random_state=0)
                  .fit(rng.uniform(lo, hi, X.shape)).inertia_) for _ in range(B)]
    return np.mean(ref) - obs, np.std(ref, ddof=1) * np.sqrt(1 + 1/B)

real, _ = make_blobs(n_samples=1200, centers=4, cluster_std=1.0, random_state=0)
noise = rng.uniform(0, 1, (1200, 2))

for name, X in (("four real clusters", real), ("uniform noise", noise)):
    for k in (2, 3, 4, 6):
        m = KMeans(k, n_init=10, random_state=0).fit(X)
        stab = stability(X, k)
        g, s = gap(X, k)
        print(f"{name:20s} k={k}  silhouette {silhouette_score(X, m.labels_):.4f}"
              f"   stability {stab:.4f}   gap {g:+.4f}   ref spread {s:.4f}")
    print()
# four real clusters   k=2  silhouette 0.4599   stability 0.9613   gap +0.7917   ref spread 0.0163
# four real clusters   k=3  silhouette 0.4777   stability 0.9920   gap +0.9961   ref spread 0.0217
# four real clusters   k=4  silhouette 0.5024   stability 0.9889   gap +1.0563   ref spread 0.0202
# four real clusters   k=6  silhouette 0.4016   stability 0.7190   gap +0.8808   ref spread 0.0134
# 
# uniform noise        k=2  silhouette 0.3640   stability 0.9435   gap +0.0036   ref spread 0.0219
# uniform noise        k=3  silhouette 0.3860   stability 0.8705   gap -0.0099   ref spread 0.0182
# uniform noise        k=4  silhouette 0.4075   stability 0.9689   gap -0.0232   ref spread 0.0251
# uniform noise        k=6  silhouette 0.3809   stability 0.7255   gap -0.0227   ref spread 0.0161

References

• Ester, M., Kriegel, H.-P., Sander, J., & Xu, X. (1996). A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. KDD.
• Campello, R. J. G. B., Moulavi, D., & Sander, J. (2013). Density-Based Clustering Based on Hierarchical Density Estimates. PAKDD.