Density-based clustering (DBSCAN)

Deep dive

How it works

Parts:

• Metric — distance function: Euclidean, cosine, haversine for lat/lon.
• Region index — answers query(i) = all j with d(i,j) <= eps: uniform grid, kd-tree, r-tree, approximate graph index (HNSW).
• Labels — per point: unseen, noise (-1) or cluster id.
• Seed list — frontier of the cluster being grown.

Operations:

1. Take the next unseen point p; run query(p).
2. If fewer than minPts results: mark p noise (it may later become a border point) and continue.
3. Otherwise open a new cluster: label p, put its neighbours in the seed list.
4. Pop a seed q. If q was noise, relabel it as this cluster's border point; do not expand.
5. If q is unseen: label it, run query(q); if q is core, append its neighbours to the seeds.
6. Repeat 4-5 until the seed list is empty; the cluster is the set of labelled points.
7. Go to 1. Each point is labelled once and queried once.

 .  n        o o o o      C core: >= minPts within eps
      b C C C C C b       b border: within eps of a core, not core itself
 n     b C C C C b   n    n noise: neither
 cores chain to one cluster when each step is <= eps

Choosing eps (k-distance plot). For each point take the distance to its k-th nearest neighbour (k tied to minPts), sort descending, plot. The knee height is a candidate eps; points above it become noise. Ester et al. 1996 propose this sorted k-dist graph.

Choosing minPts. Larger values tolerate more noise; scikit-learn advises raising it on noisy and large data. With minPts = 2 clusters are the connected components of the eps-graph (single-linkage cut); with 1 every point is core. Schubert et al. 2017 discuss heuristics and indicators of bad settings; start from 5 and inspect the result.

Metrics. For lat/lon use haversine (scikit-learn takes [lat, lon] in radians, so eps = km / 6371, mean Earth radius); a degree grid is wrong because longitude degrees shrink with latitude, so project to local metres or use a tree. For text embeddings L2-normalise and use 1 - dot.

Stories from embeddings. Cluster article embeddings in a sliding time window (bounded n) after collapsing near-duplicates (minhash); region queries come from an HNSW index (approximate: some edges are missed) or a blocked matrix product. Pick eps from a sampled k-distance plot; score on labelled stories (evaluation-metrics).

Varying density. OPTICS orders points by reachability, from which clusterings for every eps up to a bound can be read; HDBSCAN considers all eps and keeps the stable clusters, leaving a minimum cluster size to set.

Code

Rust, std only: DBSCAN over any region-query function, a grid index, and brute force for comparison, on two blobs, a ring and 1,000 uniform noise points.

use std::collections::HashMap;

type P = [f64; 2];
fn d2(a: &P, b: &P) -> f64 { (a[0] - b[0]).powi(2) + (a[1] - b[1]).powi(2) }

/// Uniform grid, cell side = eps: every eps-neighbour is in the 3x3 cells around a point.
struct Grid { eps: f64, cells: HashMap<(i64, i64), Vec<usize>> }

impl Grid {
    fn cell(&self, q: &P) -> (i64, i64) { ((q[0] / self.eps).floor() as i64, (q[1] / self.eps).floor() as i64) }
    fn new(p: &[P], eps: f64) -> Self {
        let mut g = Grid { eps, cells: HashMap::new() };
        for (i, q) in p.iter().enumerate() { g.cells.entry(g.cell(q)).or_default().push(i); }
        g
    }
    /// Region query: all j with |p_i - p_j| <= eps, including i.
    fn query(&self, p: &[P], i: usize, out: &mut Vec<usize>) {
        out.clear();
        let (cx, cy) = self.cell(&p[i]);
        for (dx, dy) in (-1..=1).flat_map(|x| (-1..=1).map(move |y| (x, y))) {
            if let Some(v) = self.cells.get(&(cx + dx, cy + dy)) {
                out.extend(v.iter().copied().filter(|&j| d2(&p[i], &p[j]) <= self.eps * self.eps));
            }
        }
    }
}

/// DBSCAN over any region query. Labels: -1 noise, 0.. cluster id. minPts counts the point itself.
fn dbscan(n: usize, min_pts: usize, mut query: impl FnMut(usize, &mut Vec<usize>)) -> Vec<i32> {
    const UNSEEN: i32 = -2;
    let mut label = vec![UNSEEN; n];
    let (mut nb, mut nb2, mut c) = (Vec::new(), Vec::new(), 0);
    for i in 0..n {
        if label[i] != UNSEEN { continue; }
        query(i, &mut nb);
        if nb.len() < min_pts { label[i] = -1; continue; } // may still become a border point
        label[i] = c;
        let mut seeds = nb.clone();
        let mut k = 0;
        while k < seeds.len() {
            let j = seeds[k];
            k += 1;
            if label[j] == -1 { label[j] = c; continue; } // border: claimed, not expanded
            if label[j] != UNSEEN { continue; }
            label[j] = c;
            query(j, &mut nb2);
            if nb2.len() >= min_pts { seeds.extend_from_slice(&nb2); } // j is core: grow
        }
        c += 1;
    }
    label
}

fn main() {
    let mut s = 0x9E3779B97F4A7C15u64; // xorshift64
    let mut u = move || { s ^= s << 13; s ^= s >> 7; s ^= s << 17; (s >> 11) as f64 / (1u64 << 53) as f64 };
    let mut p: Vec<P> = Vec::new();
    for cx in [-5.0, 5.0] { // two Gaussian blobs, sd 1 (Box-Muller)
        for _ in 0..12_000 {
            let (r, t) = ((-2.0 * (1.0 - u()).ln()).sqrt(), 6.283185307 * u());
            p.push([cx + r * t.cos(), r * t.sin()]);
        }
    }
    for _ in 0..5_000 { // ring: radius 8, width 0.1, centre (0, 12)
        let (a, r) = (6.283185307 * u(), 8.0 + 0.1 * (u() - 0.5));
        p.push([r * a.cos(), 12.0 + r * a.sin()]);
    }
    for _ in 0..1_000 { p.push([-12.0 + 24.0 * u(), -6.0 + 28.0 * u()]); } // uniform noise

    let (eps, min_pts) = (0.3, 10);
    let t = std::time::Instant::now();
    let g = Grid::new(&p, eps);
    let a = dbscan(p.len(), min_pts, |i, o| g.query(&p, i, o));
    let t_grid = t.elapsed();
    let t = std::time::Instant::now();
    let b = dbscan(p.len(), min_pts, |i, o| {
        o.clear();
        o.extend((0..p.len()).filter(|&j| d2(&p[i], &p[j]) <= eps * eps));
    });
    let t_brute = t.elapsed();

    let (ka, kb) = (*a.iter().max().unwrap() + 1, *b.iter().max().unwrap() + 1);
    let noise = a.iter().filter(|&&x| x == -1).count();
    println!("n={} clusters {} noise {} grid {:?} brute {:?}", p.len(), ka, noise, t_grid, t_brute);
    assert!(ka == 3 && kb == 3);
    assert!(noise > 800 && noise <= 1_000);
    // noise and core partition are order-independent; border labels may differ
    assert!(a.iter().zip(&b).all(|(x, y)| (*x == -1) == (*y == -1)));
}

Measured on this Mac (Apple M4 Pro, single thread, cargo build --release, 13 runs, machine not idle): n = 30,000, eps 0.3, minPts 10 gives 3 clusters and 982 noise points; the grid took 54-120 ms (typically about 60), brute force 400-800 ms. Noise and cluster count match; border labels may differ (Pitfalls).

Complexity and constants

• n region queries, each costing the index lookup plus its output size; in dense data the average neighbour count k̄ dominates: about n·(q + k̄).
• Brute force: n² distance evaluations (9·10⁸ at n = 30,000). Grid: 9 cell visits per query in 2-D, 3^d in d dimensions: it stops paying beyond a few dimensions.
• Memory: labels 4 B/point; the seed list can hold duplicates (up to the sum of neighbourhood sizes). scikit-learn stores all neighbourhoods in bulk: O(n·k̄), worst case O(n²) with large eps and low min_samples (its documentation).
• Worst case is O(n²) whatever the index; what counts is whether the index works for your data and eps (Schubert et al. 2017).

Efficiency

• Cost: measured here (M4 Pro, one core, n = 30,000, 2-D): about 2-4 µs per point with the grid, 13-27 µs brute force. Grid: bound by cell memory access and output size; brute force: compute, growing with n.
• Levers:
  • Spatial index, cell side eps — about 7x typical (4-14x range) measured at n = 30,000 (this page's code) — low dimension, small eps.
  • Collapse duplicates and pass weights — scikit-learn accepts sample_weight (docs) — repeated readings, snapped coordinates.
  • Window by time or region before clustering — cost falls with n — event streams, news, GPS days.
  • Reduce dimension (PCA) or use an approximate index — fewer distance operations — embeddings of hundreds of dimensions.
  • Parallel queries (scikit-learn n_jobs) or GPU (cuML) — not measured here — large n.
• Naive → standard → state of the art:
  • Full pairwise distance matrix · O(n²) time and memory.
  • Region queries through a kd-tree, r-tree or grid · near O(n log n) when queries are cheap, O(n²) worst case.
  • Parallel or GPU DBSCAN with an index, or HDBSCAN when density varies · same asymptotics, better constants or robustness; the approximate DBSCAN of SIGMOD 2015 offered no practical benefit with well-chosen parameters (Schubert et al.).

Implementation notes

• Count the point itself in minPts (original paper, scikit-learn); other libraries differ by one.
• Hash grid cells by integer coordinates; for large n sort points by cell for cache locality.
• Scale features first; eps is one number for all axes.

Pitfalls

• Border-point order: a border point within eps of two clusters goes to whichever reaches it first; cores and noise are order-independent, cluster ids are not (scikit-learn documents this). Treat border labels as soft.
• Chaining: a thin bridge of cores merges two clusters; lower eps or raise minPts.
• High dimension: pairwise distances concentrate, so one eps links everything or nothing; reduce dimension, use cosine on normalised embeddings, or prefer HDBSCAN.
• Degrees as metres: Euclidean on lat/lon distorts; use haversine or project.
• Loose eps on dense data: neighbour lists explode (O(n²) memory worst case).

Variants and alternatives

In production

• scikit-learn (DBSCAN, HDBSCAN, OPTICS), ELKI, and cuML on GPUs; typical uses are listed under Use when.

Implementations

Sources

• Ester, Kriegel, Sander, Xu, A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise, KDD 1996. paper
• Schubert, Sander, Ester, Kriegel, Xu, DBSCAN Revisited, Revisited: Why and How You Should (Still) Use DBSCAN, ACM TODS 42, 2017, 1-21. paper
• Ankerst, Breunig, Kriegel, Sander, OPTICS: Ordering Points To Identify the Clustering Structure, SIGMOD 1999, 49-60. paper
• Campello, Moulavi, Sander, Density-Based Clustering Based on Hierarchical Density Estimates, PAKDD 2013, 160-172. paper
• McInnes, Healy, Astels, hdbscan: Hierarchical density based clustering, JOSS 2, 205, 2017. paper · code
• scikit-learn documentation: DBSCAN · clustering guide