What it actually does
Stand on any point and count how many others are within a fixed radius ε. If at least minPts are (the point itself included), you are standing somewhere dense, and the point is a core point. Two core points within ε of each other are in the same cluster, and so is anything that chains through a sequence of such steps — the cluster is the crowd you can walk across without ever crossing a gap wider than ε in a crowded place. Points within ε of a core point but not dense themselves are border points: they join the cluster, but you cannot walk onward through them. Everything left over is noise. Nothing in this asks for the number of clusters, and nothing assumes a shape — a ring, a crescent and a blob are all just connected regions of density — which is exactly where k-means, with its convex Voronoi cells, fails. What it asks for instead is a single, global notion of "close", which is the whole difficulty: ε is a distance, so it is only as meaningful as the metric that measures it, and one ε cannot describe two clusters of different density.
- A single density threshold describes every cluster, so all clusters are of roughly comparable density — the assumption HDBSCAN relaxes and the one most often violated in real data
- The distance is meaningful and the features are on comparable scales; ε is only as good as d, and an unscaled feature quietly decides every neighbourhood
- The dimension is low enough for distance to discriminate. In high dimension distances concentrate — the ratio between the nearest and farthest neighbour tends to 1 — so ε either captures nothing or everything and the single-radius idea stops working
- Clusters are separated by regions of lower density, rather than by overlapping or touching densities; thin bridges of points chain clusters together
- Noise is meaningful: points outside every dense region are genuinely unlike the rest, not just under-sampled areas of a legitimate sparse mode
Fixed-point condition
There is no loss here, and pretending otherwise would be the wrong reading. DBSCAN does not minimize anything: no two parameter settings can be compared by a number that goes down, which is why ε and minPts cannot be tuned by the usual train-and-validate loop. What it defines is a closure. Start from one core point, repeatedly add the entire ε-neighbourhood of every core point you hold, and stop when nothing new is added. That stopping point is a fixed point of the expansion operator — apply it once more and the set does not change — and because the operator only grows sets, the smallest such fixed point containing the seed is well defined and reached in finitely many steps. That is why the objective is stated as a fixed-point condition rather than as a loss, in the same family as the Bellman equation for Q-learning: the answer is whatever the operator leaves unchanged. Equivalently, the clusters are the connected components of the graph whose vertices are the core points and whose edges join two cores within ε of each other. Border points attach to a component through one core neighbour without extending it, and noise is whatever no component reaches. The consequence worth stating plainly is that correctness is about the definition, not about optimality: there is no sense in which one valid clustering is better than another, only whether (ε, m) described the density you meant.
- the distance between two points — the single place where the data representation enters, and the choice that decides the shape of the neighbourhood (Euclidean ball, Manhattan diamond, cosine cone, great-circle cap)
- the neighbourhood radius, in the units of d — a distance, not a count, so it is only meaningful relative to the metric and the feature scaling
- the ε-neighbourhood of p: every point within ε of it, p itself included
- minPts, the density threshold: how many points must fall inside the ε-ball for its centre to count as dense
- the set of core points — the dense ones. Core points and noise are determined entirely by (ε, m); only border points are left to ordering
- the set reached after t expansion steps from a seed core point p₀; it only ever grows, so the sequence converges
- the cluster containing p₀: the smallest set that is closed under "add the whole neighbourhood of every core point already in it". Noise is every point in no such set
How it is minimized
Breadth-first closure of density-reachability, with a region query per point
- assignment, applied as the queue is drained: pop a point, claim its neighbourhood for the current cluster, and enqueue the new ones
- the id of the cluster currently being grown; it increments each time a fresh, unvisited core point seeds a new one
- the frontier queue. The "only if p is core" condition is what stops expansion at a border point: it is labelled but never queried onward
- the region query — a range search under the chosen metric, and the operation the whole runtime is made of
- each point is expanded at most once, which is what bounds the algorithm at n region queries rather than at the (much larger) number of edges in the neighbourhood graph
The closure is computed the obvious way: pick an unvisited point, query its ε-neighbourhood, and if it is dense start a cluster and grow it breadth-first, querying each newly reached point exactly once. Because every point is expanded at most once the algorithm performs at most n region queries, so its cost is entirely the cost of the range search — which is why the choice of spatial index, not any property of the clustering, decides whether it is fast. The same structure shows what is deterministic and what is not. Core points and noise are fixed by (ε, m) regardless of visiting order, because the closure is order-independent; only border points can differ, since a border point lying within ε of cores from two different clusters is claimed by whichever expansion reaches it first. That is a known, minor non-determinism, and a consequence of the definition rather than a bug. The other consequence is that the whole procedure is equivalent to connected components over the core points, which is why it parallelizes: the region queries are independent, and the merge is a union-find. There is no gradient, no learning rate and no iteration to convergence — which is also why there is nothing to warm-start and no incremental fit in the standard implementations.
eps (ε)
The neighbourhood radius, and the parameter that matters. Too small and everything is noise; too large and distinct clusters merge into one. Choose it from the k-distance plot: sort every point’s distance to its k-th nearest neighbour (k = minPts) and take the knee, the point where distances jump from "inside a cluster" to "isolated". It is a distance, so it is invalid if the features are unscaled or the metric is wrong (read off the k-distance knee; never a default. Rescaling features changes it, so fix the scaler first)
minPts (min_samples)
The density threshold. Larger values demand denser cores, absorb noise, and suppress the thin-bridge chaining that merges clusters; smaller values find small clusters and are noisier. Mind the convention: in scikit-learn it counts the point itself, and in some papers it does not (at least d + 1; 2·d is the common rule of thumb for d features, and 4 to 10 for low-dimensional spatial data. Raise it for noisy data)
metric
The distance d, which defines the shape of the ε-ball. Euclidean by default; Manhattan gives a diamond and is more tolerant of one wildly different coordinate; cosine gives a cone and compares direction only; haversine gives a great-circle cap for latitude and longitude; precomputed accepts any distance matrix, such as dynamic time warping, Jaccard or edit distance (chosen from the data representation — see the similarity notes in the applications below. Changing it invalidates ε)
algorithm / leaf_size
The spatial index behind the region query: kd-tree, ball-tree, or brute force. A cost decision, not a modelling one — except that tree indexes lose their advantage past roughly 10 to 20 dimensions and above all for non-Euclidean or non-vector metrics (kd-tree or ball-tree in low dimension; brute force or approximate neighbours in high dimension)
min_cluster_size (HDBSCAN)
HDBSCAN drops ε entirely and keeps one intuitive parameter: the smallest group worth calling a cluster. It replaces the single global ε with a hierarchy over all ε values and selects the most stable clusters from it, which is what handles clusters of different density (a domain decision — the smallest group you would act on)
Terminates after at most n region queries, deterministically for core points and noise and with a small ordering dependence for border points. There is no local optimum because there is no objective: the output is the unique closure, not the best of several. The failure modes are failures of the parameters and the metric, not of convergence, and three deserve names. One global ε: if clusters differ in density, no single ε is right — a value that separates the dense ones fragments the sparse one into noise, and a value that keeps the sparse one whole merges the dense ones. This is the problem HDBSCAN and OPTICS exist to solve. Chaining: a thin line of points bridging two clusters, each within ε of the next, connects them into one — the same single-linkage effect hierarchical clustering has, and raising minPts is the direct defence, since a bridge of border-dense points stops being core. And the curse of dimensionality, below. DBSCAN cannot fail to stop; it can only answer a question you did not mean to ask.
n region queries, so O(n · q) where q is the cost of one range search. With a kd-tree or ball-tree in low dimension q is O(log n) plus the neighbourhood size, giving O(n log n) overall. The worst case is O(n²): with brute force, with a dimension high enough that the tree prunes nothing, or with an ε so large that every neighbourhood contains a constant fraction of the data — the work is bounded below by the total size of all neighbourhoods, which an index cannot reduce. Memory is O(n) with an index, but implementations that precompute and store every neighbourhood use O(n · average neighbourhood size), which is the quiet cause of out-of-memory failures on dense data; a precomputed distance matrix, needed for non-vector metrics such as dynamic time warping, is O(n²) outright.
How, where, and why it is used
Time-Series Forecasting
AdaptedIt does not forecast, and the honest use is upstream of the forecaster. Cluster the series by shape, then pool: fit one model per cluster, or use the cluster as a feature, and route every noise series to its own model or a robust fallback, since those are exactly the ones that resemble nothing else. What makes DBSCAN better than k-means for this job is a single fact: it only needs pairwise distances, never a mean. Dynamic time warping, which handles phase shifts and local stretching, has no meaningful centroid and so cannot drive k-means, but it feeds DBSCAN directly as a precomputed distance matrix.
- Grouping SKUs or stores by demand profile so a pooled model is fitted per group, with irregular products flagged as noise instead of dragging a group’s model
- Smart-meter and energy-load segmentation by daily-profile shape, where the noise points are atypical households that need their own treatment
- Clustering sensors or machines by their operating-signal embeddings before fitting per-cluster baselines for forecasting or monitoring
The noise label is the useful part here. A partitioning method must assign every series to some group and then fits a pooled model that the odd series contaminate; DBSCAN lets the series that do not resemble anything opt out. The costs are specific. The distance matrix is O(n²) in the number of series and each DTW entry is itself quadratic in series length, so it does not scale past a few thousand series without an approximate lower-bound or a learned embedding. A single ε assumes one similarity scale across the whole catalogue, which is rarely true when fast and slow movers coexist. And it has no predict method, so a new series is placed by nearest-core assignment, which is a rule you must build. Where those bite, cluster learned embeddings with k-means or HDBSCAN instead.
- Z-normalize each series so the distance compares shape, otherwise the clusters recover magnitude and ε means nothing across scales
- Compute DTW (or a shape-based distance) with a Sakoe-Chiba band, and feed the matrix as metric="precomputed"; an unbanded DTW is both slow and prone to pathological alignments
- Pick ε from the k-distance curve of that same matrix, and re-check it whenever the series length or normalization changes
- Fit the clustering on the training window only — clustering over the full history leaks future shape into the pooled-model assignment
Rolling-origin backtest of the pooled forecasts against per-series and global-model baselines, not a clustering metric. A clustering that scores well on silhouette and does not lower forecast error described the data without helping it, and the noise series should be scored separately so that their (usually worse) error is visible rather than averaged away.
- Euclidean distance on unaligned series, where a one-period phase shift looks like a large difference — the reason a warping distance is needed here at all
- An O(n²) DTW matrix on tens of thousands of series, which fails on memory long before time
- A large noise share, which means ε was set for the dense minority and most of the catalogue got no pooled model
- Cluster ids that change between refits, silently changing which pooled model a series is scored by
Anomaly Detection
Primary fitThe output of the algorithm is the detector. Run it on the feature representation of normal-plus-unknown data, and the points labelled noise (-1) are the ones with no dense neighbourhood — the outliers — while the clusters describe what normal looks like, with no assumption about how many normal modes there are or what shape they take. Two readings are available and both are useful: the noise points are the point anomalies, and a cluster that is unusually tight, small and appearing suddenly is often a coordinated group of anomalies that found each other.
- Spatial and sensor outlier screening, such as GPS fixes or readings that belong to no dwell or operating cluster
- Transaction and event screening in a space of amount, time, location and device features, where the benign flows form dense modes and unusual ones are isolated
- Network-traffic baselining, where normal flows cluster by feature profile and scans or exfiltration fall outside every dense region
- Manufacturing and inspection data, separating spatially clustered defects (a process problem) from isolated random ones
Noise is a first-class output rather than an afterthought, which is the structural advantage over k-means (it assigns every point to a cluster and then needs a distance threshold bolted on) and over a Gaussian mixture (it needs a density model with a chosen number of components). The method makes no assumption about how many normal modes there are or what shape they have, and it needs no labelled anomalies. Its limits are as important. The output is binary: there is no score, so ranking alerts needs a second step — LOF, HDBSCAN’s GLOSH outlier scores, or simply the k-distance — and a hard noise threshold is brittle. One global ε again: a normal point in a sparse but legitimate region is flagged, and a tight anomaly cluster is called normal. And it fails in high dimension, where distances concentrate and every point looks equally isolated; there, an isolation forest, which uses no distances at all, is the better detector. Metric choice is part of the model: Euclidean on standardized features, cosine for embeddings, haversine for coordinates, and a wrong choice silently re-draws what counts as an outlier.
- Standardize numeric features so that ε means the same thing on every axis, and make the scaler part of the model artefact
- Reduce dimension first (PCA or a learned embedding), since a single ε is meaningless in dozens of raw dimensions
- Choose the metric from the representation: Euclidean for standardized coordinates, cosine for embeddings (L2-normalize and use Euclidean), haversine for latitude and longitude in radians, a precomputed matrix for sets or sequences
- Set minPts from how large a group of events you are willing to call normal, then read ε from the k-distance plot for that k
- Clean the fit window, or a dense cluster of contamination is learned as a normal mode
Precision@k and PR-AUC against confirmed incidents, never accuracy, with the noise points ranked by a continuous score (k-distance, LOF or GLOSH) so that a threshold can be chosen. Report the noise rate alongside precision: an alert rate you cannot staff is a failure regardless of how many true anomalies it contains. Compare against an isolation forest on the same features — if it wins clearly, the dimension or the varying density is the thing costing you.
- Treating the noise label as a calibrated anomaly score, when it is a yes/no verdict at one (ε, minPts)
- A coordinated, dense group of anomalies forming its own cluster and being reported as perfectly normal
- One global ε over regions of different density, flagging every point of a sparse-but-legitimate segment
- Unscaled or high-dimensional features, where the k-distance plot has no knee and every point looks equally isolated
- No predict method: scoring a new event requires a deliberate nearest-core rule, and doing it differently in training and serving shifts the verdicts
Optimization
Not applicableDBSCAN solves no optimization problem: its output is a closure under a density-reachability operator, not the minimizer of an objective, so it is not a solver for routing, scheduling or allocation. The only optimization attached to it is tuning ε and minPts, which is external to the algorithm — done by reading a k-distance plot or by an application-level metric — and there is nothing to optimize inside it.
Risk & Fraud
ViableTwo readings of one output. The noise points are the isolated, unlike-anyone events worth a look, and the clusters are where the structure is: a dense cluster of accounts that share a device fingerprint, an address pattern and a creation-time burst is a fraud ring, which is a cluster in a feature space where the analyst has chosen what "close" means. The signal is often the cluster, not the noise.
- Account-farm and bot-ring discovery, clustering accounts by device, network and behavioural features
- Transaction clustering to find unusually tight groups of near-identical payments, as in card testing
- Reducing a review queue by labelling the dense, well-understood bulk and surfacing the leftover for investigation
Fraud is adversarial and clustered by construction — attackers reuse infrastructure — so a method that finds arbitrary-shaped dense groups without a preset count is a natural fit, and the noise output gives a second, independent signal. It is an exploratory and triage tool rather than a decision engine: it has no labels and no calibrated score, so it complements a supervised model rather than replacing it. Metric choice does the real work, since similarity here is mixed — Jaccard on sets of shared identifiers, cosine on behaviour embeddings, Euclidean on standardized numeric features — and a precomputed combination is common.
- Build the representation from the investigator’s notion of "the same actor": shared devices, IPs, addresses and timing, encoded so that distance reflects that notion
- Use Jaccard distance on sets of shared identifiers and cosine on embeddings, via a precomputed matrix or L2-normalized vectors
- Down-weight ubiquitous identifiers such as a shared corporate NAT address, or the cluster is the whole company
- Fix minPts at the smallest group you would investigate, and read ε from the k-distance plot at that k
Precision of the surfaced clusters and noise points against confirmed investigations, and the investigator hours per confirmed finding, rather than any internal clustering metric. Hold out a recent time slice to confirm the clusters are not simply last week’s known cases.
- A shared infrastructure identifier merging unrelated accounts into one enormous cluster through a chaining effect
- Treating cluster membership as evidence of fraud, when it is only evidence of similarity
- A fixed ε drifting out of validity as the population or the feature scaling changes
Control & Operations
ViableGeographic, trajectory and spatial-sensor clustering, which is where DBSCAN is most at home. Latitude-longitude fixes are clustered under the haversine distance, so ε is a real-world distance on the sphere; a stay-point detector clusters the fixes of a trajectory to find where it dwelt; and a robot perception stack clusters a LiDAR point cloud, after ground removal, into obstacle candidates under Euclidean distance. Spatio-temporal variants use two radii, one in space and one in time.
- Stay-point and place detection from GPS trajectories, and hotspot discovery for pickups, incidents or demand
- Obstacle segmentation in LiDAR point clouds, where each dense cluster is a candidate object
- Vessel, vehicle and aircraft position data, clustering stops and anchorages and flagging off-pattern movements as noise
- Spatial siting and service-area analysis, where the real shape of a hotspot is not a circle
This is the problem the algorithm was designed for: spatial data of low dimension, a physically meaningful ε (metres or kilometres), no known count of clusters, arbitrary shapes — a road, a coastline, a crescent-shaped crowd — and noise that is real (a stray fix, a one-off) rather than a modelling nuisance. The metric is not a detail: Euclidean distance on raw latitude and longitude is wrong, because a degree of longitude shrinks with latitude, so the data must be projected or the haversine metric used. The trade-offs are the usual ones. A single ε assumes uniform density, so a city centre and its suburbs need either different ε or HDBSCAN, and a ball-tree is the index that supports haversine.
- Use the haversine metric on latitude and longitude in radians, or project to a local planar coordinate system and use Euclidean — never Euclidean on raw degrees
- Convert ε from a physical distance to radians (divide kilometres by the Earth radius) when using haversine, and keep the unit in the parameter name
- For trajectories, thin the fixes by time or distance first, so a slow-moving point does not look dense simply because it was sampled often
- For spatio-temporal clustering, use two radii, one for space and one for time, rather than mixing units in a single distance
- For point clouds, remove ground points first and consider a voxel-grid downsample, since the cost is the number of region queries
Task-level: place-detection precision and recall against labelled visits, obstacle recall at a fixed false-positive rate for perception, and the downstream decision the cluster feeds. Check stability under a small perturbation of ε, since a hotspot that vanishes when ε moves by ten percent is an artefact.
- Euclidean distance on latitude and longitude in degrees, which distorts every neighbourhood away from the equator
- Forgetting that haversine takes radians, and that its ε is an angle on the unit sphere rather than a distance
- A fixed ε across a dense city centre and sparse outskirts, which fragments the outskirts into noise
- Sampling-rate artefacts in trajectories, where the number of fixes rather than the time spent produces the density
What it costs to run
n range queries, so O(n log n) with a spatial index in low dimension and O(n²) in the worst case; on 2-D spatial data, millions of points cluster in seconds to minutes. The index build is O(n log n) and is a one-off. With a precomputed distance matrix the cost and the memory are both O(n²), which is the practical ceiling for DTW or edit-distance clustering at a few tens of thousands of items.
There is no predict method in the standard implementation, because the model is the data. A new point is placed by a nearest-core rule: find the closest core point, and if it is within ε adopt that core’s cluster, otherwise label the point noise. That is one nearest-neighbour query against the stored core points, so the artefact to ship is the core points themselves plus ε, the metric and the scaler — which can be large — and an approximate index when it is. Note that this rule is a deliberate, separate piece of code, and it must match the training-time definition exactly or the labels disagree.
Batch refit; there is no incremental update in the standard implementation, and the refit is cheap enough that it is the reason not to bother with incremental variants. The reason to refit rarely is stability rather than cost: cluster ids are arbitrary and change on every run, and a drift in the data moves ε, so anything keyed to a cluster id must be remapped or re-approved. Re-derive ε from the k-distance plot on each refit instead of carrying a stale value forward.
- Track the noise rate over time — a rising rate is either a real wave of anomalies or an ε that no longer matches the density, and the k-distance plot distinguishes them
- Watch the number and size distribution of clusters; one cluster swallowing the rest points at chaining or a drifted ε, and a collapse into noise points at the opposite
- Re-plot the k-distance curve on recent data and compare the knee to the deployed ε — a moved knee is drift measured directly
- Monitor the share of new points that fall within ε of a core point at serving time, which is the coverage of the current model
- The noise label is -1 and is easily consumed as a valid cluster id by downstream code that indexes with it
- The meaning of min_samples differs between implementations — whether the point itself counts — and a port between libraries silently changes which points are core
- The scaler is part of the model: ε is a distance in scaled units, so refitting the scaler changes every neighbourhood without any code change
- Haversine requires radians and returns angles on the unit sphere, so an ε in kilometres must be divided by the Earth radius or every point is noise
- Dense data with a large ε makes implementations that store every neighbourhood run out of memory long before they run out of time
- Duplicate points inflate neighbourhood counts, so a heavily repeated record can make a point core on its own; collapse duplicates with sample weights
- Border points can differ between runs and between implementations; compare core points and noise, not raw labels
Pros and cons, in context
No number of clusters to choose, and arbitrary cluster shapes
A ring, a crescent and a blob are all connected regions of density, so the method finds them where k-means, with its convex Voronoi cells, splits them down the middle. The saving is real but it is exchanged, not removed: the two parameters ε and minPts replace K, and they are harder to guess.
Noise is a first-class output
Outliers are labelled instead of being forced into the nearest cluster, which is both a better clustering (they do not drag a centroid) and a ready-made anomaly detector. It is the right property when the data contain genuine strays and the wrong one when every point must be assigned, in which case the noise needs a second step.
Works with any distance, including non-vector ones
It needs only pairwise distances, never a mean, so it runs on a precomputed matrix of dynamic-time-warping, Jaccard or edit distances where k-means has no centroid to compute. The price is O(n²) memory in exactly those cases.
Deterministic for core points and noise, and robust to outliers
There is no initialization and no local optimum to restart from; the same data and parameters give the same dense structure. Only border points depend on ordering, which seldom matters.
Fast in low dimension
O(n log n) with a kd-tree or ball-tree, which makes it the standard tool for geographic and spatial data of two or three dimensions. The advantage disappears with dimension, where the index stops pruning.
One global ε cannot describe clusters of different density
A radius that separates the dense clusters fragments the sparse one into noise, and one that keeps the sparse one whole merges the dense ones. This is the characteristic failure on real data, and the reason HDBSCAN and OPTICS exist; if densities clearly differ, start there.
Fails in high dimension
Distances concentrate as the dimension grows, so the contrast between near and far collapses and no ε separates dense from sparse. Reduce dimension first or use a method that does not rely on distance; this is the main place an isolation forest beats it as a detector.
Sensitive to ε and the metric, with no loss to tune against
There is no objective to minimize, so ε and minPts cannot be validated the way a supervised model’s parameters are. The k-distance plot gives a defensible ε, but it is a heuristic with a judgement call at the knee, and a rescaling or a metric change invalidates it.
No predict method, and no native anomaly score
The model is the data: new points need a hand-built nearest-core rule, and the noise label is binary. Both are addressable (nearest-core assignment, LOF or GLOSH) and both are routinely forgotten until deployment.
Chaining through thin bridges
A sparse line of points connecting two clusters merges them, just as in single-linkage hierarchical clustering. Raising minPts is the defence, at the cost of discarding genuinely small clusters as noise.
Intuitive, idiomatic, optimized
The same algorithm in three languages, each progressing through Kent Beck's three stages. Every optimization names the lever it pulls and the trade-off it accepts — cost profiles are illustrative, not measured benchmarks.
The math, transcribed literally. No dependencies.
"""DBSCAN - the definition, transcribed.
A point is CORE if at least min_pts points (itself included) lie within eps of
it. A cluster is grown from a core point by repeatedly claiming the whole
eps-neighbourhood of every core point it reaches. A point that is within eps of
a core point but is not core itself is a BORDER point: it joins the cluster
but the growth does not continue through it. Everything unreached is NOISE.
The only place the data representation enters is distance(): change it and the
shape of the neighbourhood changes with it.
"""
from collections import deque
NOISE = -1
UNVISITED = -2
def distance(a, b):
"""Euclidean distance. This is the choice that makes the eps-ball round."""
total = 0.0
for j in range(len(a)):
difference = a[j] - b[j]
total += difference * difference
return total ** 0.5
def region_query(X, i, eps):
"""Every index within eps of point i, point i itself included."""
neighbours = []
for j in range(len(X)):
if distance(X[i], X[j]) <= eps:
neighbours.append(j)
return neighbours
def fit(X, eps, min_pts):
n = len(X)
labels = [UNVISITED] * n
cluster = 0
for seed in range(n):
if labels[seed] != UNVISITED:
continue
neighbours = region_query(X, seed, eps)
if len(neighbours) < min_pts:
# Not dense. Tentatively noise - a later cluster may still
# claim it as a border point.
labels[seed] = NOISE
continue
# seed is a core point: start a cluster and grow it.
labels[seed] = cluster
queue = deque(neighbours)
while queue:
candidate = queue.popleft()
if labels[candidate] == NOISE:
# Was written off as noise but sits within eps of a core
# point: it is a border point. Label it, do not expand it.
labels[candidate] = cluster
continue
if labels[candidate] != UNVISITED:
continue
labels[candidate] = cluster
reachable = region_query(X, candidate, eps)
if len(reachable) >= min_pts:
# candidate is core too, so its neighbourhood is reachable.
queue.extend(reachable)
cluster += 1
return labels
def k_distances(X, k):
"""Sorted distance from every point to its k-th nearest neighbour.
Plot this and read eps off the knee: distances inside a cluster are small
and flat, and the curve jumps where points become isolated.
"""
result = []
for i in range(len(X)):
distances = sorted(distance(X[i], X[j]) for j in range(len(X)))
result.append(distances[min(k, len(distances) - 1)])
return sorted(result)O(n^2) distance evaluations in interpreter loops: every region query scans all n points.
How this language actually wants it written.
"""DBSCAN - typed, validated, tree-backed region queries, a real predict()."""
from __future__ import annotations
from collections import deque
from dataclasses import dataclass, field
from typing import Literal
import numpy as np
from numpy.typing import NDArray
from scipy.spatial import cKDTree
Matrix = NDArray[np.float64]
Labels = NDArray[np.int64]
NOISE = -1
UNVISITED = -2
Metric = Literal["euclidean", "manhattan", "chebyshev", "cosine"]
# Minkowski p for each metric the kd-tree can serve directly. The metric is not
# a detail: it decides the SHAPE of the eps-ball - a sphere, a diamond, a cube.
_MINKOWSKI_P: dict[str, float] = {
"euclidean": 2.0,
"manhattan": 1.0,
"chebyshev": float("inf"),
"cosine": 2.0,
}
def _normalise(X: Matrix, metric: Metric) -> Matrix:
"""Cosine distance compares direction only, so project onto the unit sphere.
On unit vectors ||a - b||^2 = 2 (1 - cos(a, b)), so a cosine radius is just
a Euclidean radius and the same tree serves it.
"""
if metric != "cosine":
return X
norms = np.linalg.norm(X, axis=1, keepdims=True)
if np.any(norms == 0.0):
raise ValueError("cosine distance is undefined for zero vectors")
return X / norms
def _radius(eps: float, metric: Metric) -> float:
return float(np.sqrt(2.0 * eps)) if metric == "cosine" else eps
@dataclass(frozen=True)
class DbscanModel:
"""A fitted model. The 'model' is the data: core points, eps, the metric."""
labels: Labels
core_indices: NDArray[np.int64]
n_clusters: int
eps: float
metric: Metric
_core_points: Matrix = field(repr=False, compare=False)
_core_tree: cKDTree | None = field(repr=False, compare=False)
def predict(self, X: Matrix) -> Labels:
"""Nearest-core rule: adopt the cluster of the closest core point if it
is within eps, otherwise NOISE. sklearn has no such method, and a port
that implements it differently in serving shifts every verdict."""
if X.ndim != 2 or X.shape[1] != self._core_points.shape[1]:
raise ValueError(
f"expected (n, {self._core_points.shape[1]}) input, got {X.shape}"
)
out = np.full(X.shape[0], NOISE, dtype=np.int64)
if self._core_tree is None:
return out
points = _normalise(X, self.metric)
distance, nearest = self._core_tree.query(
points,
k=1,
p=_MINKOWSKI_P[self.metric],
distance_upper_bound=_radius(self.eps, self.metric),
)
reachable = np.isfinite(distance)
out[reachable] = self.labels[self.core_indices[nearest[reachable]]]
return out
def fit(
X: Matrix,
eps: float,
min_samples: int,
metric: Metric = "euclidean",
) -> DbscanModel:
"""Cluster X. min_samples counts the point itself, as scikit-learn does.
Raises ValueError on malformed input.
"""
if X.ndim != 2 or X.shape[0] == 0:
raise ValueError(f"X must be a non-empty 2-D array, got shape {X.shape}")
if not np.isfinite(X).all():
raise ValueError("X contains NaN or infinity")
if eps <= 0.0:
raise ValueError(f"eps must be positive, got {eps}")
if min_samples < 1:
raise ValueError(f"min_samples must be at least 1, got {min_samples}")
if metric not in _MINKOWSKI_P:
raise ValueError(f"metric must be one of {sorted(_MINKOWSKI_P)}, got {metric}")
points = _normalise(X, metric)
radius = _radius(eps, metric)
p = _MINKOWSKI_P[metric]
tree = cKDTree(points)
n = points.shape[0]
labels = np.full(n, UNVISITED, dtype=np.int64)
is_core = np.zeros(n, dtype=bool)
frontier: deque[int] = deque()
cluster = 0
for seed in range(n):
if labels[seed] != UNVISITED:
continue
neighbours = tree.query_ball_point(points[seed], r=radius, p=p)
if len(neighbours) < min_samples:
labels[seed] = NOISE
continue
is_core[seed] = True
labels[seed] = cluster
frontier.clear()
frontier.extend(neighbours)
while frontier:
candidate = frontier.popleft()
if labels[candidate] == NOISE:
labels[candidate] = cluster # border: label, never expand
continue
if labels[candidate] != UNVISITED:
continue
labels[candidate] = cluster
reachable = tree.query_ball_point(points[candidate], r=radius, p=p)
if len(reachable) >= min_samples:
is_core[candidate] = True
frontier.extend(reachable)
cluster += 1
core_indices = np.flatnonzero(is_core).astype(np.int64)
core_points = points[core_indices]
return DbscanModel(
labels=labels,
core_indices=core_indices,
n_clusters=cluster,
eps=eps,
metric=metric,
_core_points=core_points,
_core_tree=cKDTree(core_points) if core_indices.size else None,
)What changed: The brute-force region query becomes a kd-tree range search, which is the change that moves the cost from O(n^2) to roughly O(n log n) in low dimension — the algorithm is n region queries, so the index is the algorithm. The metric is now an explicit, validated choice rather than a hard-wired Euclidean: Manhattan and Chebyshev are the same tree with a different Minkowski p (a diamond and a cube instead of a sphere), and cosine is handled by projecting onto the unit sphere, where a cosine radius becomes a Euclidean one. The fitted result is a frozen dataclass that keeps the core points, because a DBSCAN model is its data, and it gains the predict() the standard implementation lacks — a nearest-core rule with a hard radius — so that serving and training use one definition. Inputs are validated up front, with the min_samples convention (self-inclusive) stated.
O(n log n) region queries in low dimension, but still one Python call per point; degrades toward O(n^2) as dimension grows and the tree stops pruning.
Faster, and the justification for each change.
"""DBSCAN as connected components - batched queries, no Python loop over points.
DBSCAN clusters are the connected components of the graph whose vertices are
the core points and whose edges join two cores within eps. That restatement
removes the queue-and-expand loop entirely: count neighbours in one batched
call, find the cores, build the core-core edge list, and hand it to a sparse
connected-components routine. Border points then attach to their nearest core.
"""
import numpy as np
from numpy.typing import NDArray
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from scipy.spatial import cKDTree
Matrix = NDArray[np.float64]
Labels = NDArray[np.int64]
NOISE = -1
def k_distance_curve(X: Matrix, k: int, p: float = 2.0) -> NDArray[np.float64]:
"""Sorted k-th neighbour distance for every point, in one batched query."""
design = np.ascontiguousarray(X, dtype=np.float64)
distances, _ = cKDTree(design).query(design, k=k, p=p, workers=-1)
return np.sort(distances[:, -1])
def dbscan(X: Matrix, eps: float, min_samples: int, p: float = 2.0) -> Labels:
points = np.ascontiguousarray(X, dtype=np.float64)
n = points.shape[0]
# 1. Every neighbourhood SIZE at once. return_length avoids materializing
# n Python lists of indices, which is where the time and memory went.
tree = cKDTree(points)
counts = tree.query_ball_point(
points, r=eps, p=p, return_length=True, workers=-1
)
is_core = counts >= min_samples
labels = np.full(n, NOISE, dtype=np.int64)
core_index = np.flatnonzero(is_core)
if core_index.size == 0:
return labels
# 2. Core-core edges, then connected components over the core graph.
core_points = points[core_index]
core_tree = cKDTree(core_points)
pairs = core_tree.query_pairs(r=eps, p=p, output_type="ndarray")
m = core_index.size
graph = coo_matrix(
(np.ones(len(pairs), dtype=np.int8), (pairs[:, 0], pairs[:, 1])),
shape=(m, m),
).tocsr()
_, core_labels = connected_components(graph, directed=False)
labels[core_index] = core_labels
# 3. Border points: nearest core within eps, else stay noise. One batched
# nearest-neighbour query, with the radius as a hard upper bound.
border = np.flatnonzero(~is_core)
if border.size:
distance, nearest = core_tree.query(
points[border], k=1, p=p, distance_upper_bound=eps, workers=-1
)
reachable = np.isfinite(distance)
labels[border[reachable]] = core_labels[nearest[reachable]]
return labelsWhat changed: The per-point queue-and-expand loop is gone. DBSCAN is the connected components of the core graph, so the algorithm becomes four array operations: one batched call that returns every neighbourhood size (with return_length, so no list of indices is ever built), a boolean mask for the cores, one query_pairs for the core-core edges, and a sparse connected-components pass that replaces the breadth-first growth. Border points are attached by a single batched nearest-core query with the radius as a hard bound. A side effect is that border assignment becomes deterministic (nearest core) instead of depending on visiting order. Every step runs in compiled code across worker threads.
Eliminate Python-level loops over samples
The seed loop and the frontier queue become a batched neighbour count, a sparse connected-components call and a batched nearest-core query — three compiled passes instead of one interpreter iteration and one tree call per point.
Trade-off: The edge list is materialized: its size is the number of core-core pairs, which in a dense region approaches the square of that region’s point count, so a large eps on dense data can exhaust memory where the queue version only ever held one frontier.
Batch work to amortize interpreter overhead
query_ball_point with return_length and workers=-1 sizes every neighbourhood in one call across all cores, and the k-distance curve is one batched k-nearest-neighbour query rather than n sorted full scans.
Trade-off: Counting neighbourhoods and then collecting edges queries the tree twice; the single-pass loop would query once. The extra pass is paid for by being parallel and vectorized, and it is a loss on small inputs where overhead dominates.
Ensure contiguous memory layout and a single dtype
The points are forced to one contiguous float64 buffer before the tree is built, so the tree and every query read memory sequentially and never copy or convert per call.
Trade-off: It copies the input when it arrives in another layout or dtype, doubling peak memory for the largest array in the problem; float32 would halve it and cost precision near the eps boundary.
Pre-allocate output arrays and use in-place operations
The label array is allocated once and filled by two masked assignments — cores, then reachable borders — with noise as the default, so no per-cluster list is ever built or concatenated.
Trade-off: Masked fancy-index assignment still allocates temporary index arrays on the right-hand side, so it is cheaper than a loop but not allocation-free.
O(n log n) tree queries run across all cores plus a linear-time connected-components pass; memory O(n + core-core edges). Illustrative, not a measured benchmark.