DBSCAN, in depth: core, border and noise points, a verified worked example, choosing eps and minPts, and the real complexity

Density as the definition of a cluster

DBSCAN has two parameters. eps (ε) is a radius, and minPts is a count. The ε-neighbourhood of a point p is every point within distance ε of p, including p itself. Conventions differ on that last detail: the original paper and scikit-learn both count the point itself, so with minPts = 4 a point needs three other points within ε.

Every point then gets exactly one of three roles:

• Core point: its ε-neighbourhood holds at least minPts points. Core points are the interior of a cluster.
• Border point: not core, but inside the ε-neighbourhood of at least one core point. Border points sit on the edge of a cluster and belong to it, but they do not extend it.
• Noise point: neither. Noise is a first-class output, labelled -1 by most libraries.

Two relations turn roles into clusters. Point q is directly density-reachable from p if p is core and q lies in p's ε-neighbourhood. q is density-reachable from p if a chain of core points links them, each directly reachable from the previous one. Reachability is not symmetric, because a border point can be reached but cannot reach anything. So DBSCAN uses density-connectedness: p and q are connected if some point o can reach both. A cluster is a maximal set of density-connected points. Read carefully, this says that the core points form a graph (edges between core points within ε), clusters are the connected components of that graph, and each border point is attached to a component it touches.

The algorithm

The algorithm is a breadth-first search over that implicit graph. Visit each point once. If it is not core, provisionally mark it noise. If it is core and unlabelled, start a new cluster and expand outward: every neighbour joins the cluster, and every neighbour that is itself core contributes its own neighbours to the queue. A point first marked noise can be upgraded to a border point when an expansion reaches it.

from collections import deque
import math

NOISE = -1

def region_query(points, i, eps):
    """Indices within eps of points[i], including i. O(n) here; use a spatial index for real data."""
    return [j for j, q in enumerate(points) if math.dist(points[i], q) <= eps]

def dbscan(points, eps, min_pts):
    labels = [None] * len(points)          # None = unvisited
    cluster = -1
    for i in range(len(points)):
        if labels[i] is not None:
            continue
        neighbours = region_query(points, i, eps)
        if len(neighbours) < min_pts:
            labels[i] = NOISE              # may be upgraded to border later
            continue
        cluster += 1
        labels[i] = cluster
        queue = deque(neighbours)
        while queue:
            j = queue.popleft()
            if labels[j] == NOISE:
                labels[j] = cluster        # noise becomes a border point
            if labels[j] is not None:
                continue
            labels[j] = cluster
            nj = region_query(points, j, eps)
            if len(nj) >= min_pts:         # only core points extend the cluster
                queue.extend(nj)
    return labels

Two details carry the correctness. First, a point is expanded only if it is core, which is what stops clusters bleeding through a thin line of border points. Second, the noise label is provisional; a point visited early in the outer loop may later turn out to be on the edge of a cluster. Each point is labelled once and queried at most once.

Worked example: eleven points

Take the eleven points A(1,1), B(2,1), C(1,2), D(2,2), E(3,2), J(4.4,2), F(8,8), G(9,8), H(8,9), I(9,9) and K(5,6), with ε = 1.5 and minPts = 4. Counting neighbourhoods with the point itself included:

The outer loop starts at A, which is core, so cluster 0 begins. The queue pulls in B, C and D, all core; B and D add E; E is core and adds J. J has only two points in its neighbourhood, so it joins cluster 0 as a border point and the expansion stops there. The loop skips to F, which starts cluster 1 and gathers G, H and I. K is visited last, has no neighbours at all, and stays noise. k-means with k = 2 would force K into a cluster and drag that centroid toward it.

Now shrink ε to 1.39, just below the square diagonals of 1.414. Only D keeps four neighbours, so it is the sole core point; cluster 0 becomes D plus border points B, C and E, while A, J, K and all of F, G, H and I become noise (scikit-learn agrees). A 7% change in ε dissolved a whole cluster, because it crossed the plateau of typical neighbour distances. That is why choosing ε deserves a method.

Choosing ε and minPts

minPts first. It controls how much noise you tolerate and how smooth the clusters are. The usual starting rules are minPts ≥ d + 1 for d dimensions, and 2·d as a default for low-dimensional data (Sander et al. 1998); raise it for noisy or very large data sets. minPts = 1 or 2 degenerates to single-linkage clustering cut at height ε, which chains clusters together through any thin bridge of points.

Then ε from a k-distance plot. For every point compute the distance to its k-th nearest neighbour, with k = minPts under the self-included convention, sort the values and plot them. Points inside clusters have small, similar k-distances, which form a flat plateau; noise points have large ones, which form a steep tail. Put ε at the knee where the curve turns upward.

import numpy as np
from sklearn.neighbors import NearestNeighbors

def k_distances(X, min_pts):
    # kneighbors on the training data returns each point as its own first neighbour,
    # which matches the self-included minPts convention used by scikit-learn's DBSCAN.
    nn = NearestNeighbors(n_neighbors=min_pts).fit(X)
    dist, _ = nn.kneighbors(X)
    return np.sort(dist[:, -1])

Scale features before any of this. ε is one radius in every direction, so a feature measured in thousands swamps one measured in fractions. Standardise, or better, choose units in which one unit of distance means the same thing on every axis. For geographic data skip scaling and use a real distance: scikit-learn accepts metric='haversine' with algorithm='ball_tree', latitude and longitude in radians, and ε expressed in radians (kilometres divided by the Earth's mean radius, about 6371 km).

Complexity: what the index buys you

The cost is n region queries plus the total size of all neighbourhoods. With a linear scan that is Θ(n²) distance computations, painful beyond a few tens of thousands of points. With a spatial index (k-d tree, ball tree, R-tree or a grid) each query costs roughly the logarithm of n plus the size of the answer, which is where the often-quoted O(n log n) comes from.

That figure is an average for well-chosen parameters in low dimensions, not a guarantee. Three things break it. If ε is large relative to the data, every neighbourhood holds a large fraction of the points and the output size alone is quadratic. In high dimensions tree indexes degrade toward a linear scan, and distances concentrate so that the k-distance plot loses its knee. And in theory, Gan and Tao (SIGMOD 2015) showed that exact DBSCAN in three or more dimensions cannot beat roughly n to the power 4/3 time unless a long-standing hardness conjecture fails; Schubert, Sander, Ester, Kriegel and Xu revisited the question in ACM TODS 2017 and concluded that with sensible parameters the original algorithm with an index remains competitive in practice.

Memory is the trap that bites first in scikit-learn. Its implementation computes all neighbourhoods up front and stores them, so a large ε on a few million points can exhaust RAM long before CPU time is a problem. Mitigations: a smaller ε, deduplicating points and passing sample_weight for the counts, precomputing a sparse radius-neighbours graph in chunks and passing metric='precomputed', or moving to a grid-based or GPU implementation such as the one in RAPIDS cuML.

Border points and determinism

DBSCAN is deterministic for core points and noise: whether a point is core depends only on its neighbourhood, and the clusters of core points are connected components, which do not depend on visit order. Border points are different. A border point within ε of core points from two different clusters joins whichever cluster's expansion reaches it first, so shuffling the input can move it. That breaks naive tests that compare label arrays.

There are three responses. Compare clusterings with a permutation-invariant score such as the adjusted Rand index rather than raw labels. Sort the input deterministically before clustering if you need reproducible output. Or use the variant Campello and colleagues call DBSCAN*, which treats border points as noise so that the result is fully order-independent; this is also the foundation HDBSCAN builds on.

Using scikit-learn in practice

In practice you will call a library. The scikit-learn version exposes the same two parameters, with min_samples playing minPts, and reports which samples were core:

import numpy as np
from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler

X = StandardScaler().fit_transform(raw_features)          # make eps mean the same on every axis
model = DBSCAN(eps=0.3, min_samples=10, n_jobs=-1).fit(X)

labels = model.labels_                                    # -1 = noise
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
noise_ratio = float(np.mean(labels == -1))
core_mask = np.zeros(len(X), dtype=bool)
core_mask[model.core_sample_indices_] = True
border_mask = (labels != -1) & ~core_mask

print(f"{n_clusters} clusters, {noise_ratio:.1%} noise, {border_mask.sum()} border points")

Log those three numbers on every run; a jump in the noise ratio signals that ε no longer fits the data. DBSCAN has no predict method, because assigning a new point needs the core points: to label streaming data, store the core samples and assign a new point to the cluster of any core point within ε, or noise otherwise.

DBSCAN, OPTICS, HDBSCAN or k-means

DBSCAN's real limitation is the single ε. A dense city-centre cluster and a sparse suburban one cannot both be found with one radius: a small ε fragments the suburb into noise, a large one merges the city into one blob. OPTICS orders points by reachability so that every ε below a maximum can be read off one plot. HDBSCAN builds the full hierarchy of DBSCAN* clusterings over all ε and keeps the most stable clusters, which is usually the better default when you do not know the density in advance. DBSCAN remains right when the density scale is known physically, such as GPS fixes within 50 metres, because then ε is a requirement rather than a guess.

Failure modes

• Everything in one cluster: ε too large, or features unscaled so one axis dominates. Check the k-distance plot on the scaled data.
• Everything noise: ε too small or minPts too high for the sample size; often happens after subsampling, which thins density.
• Bridges between clusters: a thin chain of points connects two groups. Raise minPts so the chain points stop being core.
• High-dimensional embeddings: distances concentrate and no ε works. Reduce dimensionality first or use cosine distance with HDBSCAN.
• Out-of-memory in fit: neighbourhoods are being materialised; shrink ε, deduplicate with sample weights, or precompute a sparse graph.
• Flaky tests: border points moved with input order; compare with the adjusted Rand index.

What to do next

1. Type in the thirty-line implementation and reproduce the eleven-point example, then check it against scikit-learn's labels.
2. Generate two interleaved half-moons with noise and compare DBSCAN with k-means on the same data.
3. Plot the k-distance curve for one of your own data sets and pick ε at the knee; record the noise ratio.
4. Replace the linear region query with a spatial index, using k-d trees as a guide, and time both at 10,000 and 100,000 points.
5. Rewrite the cluster step as connected components of core points with union-find, and compare it with the breadth-first version.
6. Run HDBSCAN on the same data and decide which result you would ship, and why.