DBSCAN

In 2015–2025 the Earth had 19,073 earthquakes of magnitude 5 and above. The map shows they are not scattered at random: they stretch in chains along the boundaries of lithospheric plates, thicken into swarms and now and then happen alone in the middle of a plate. How many groups are there? k-means has to be told the number in advance, and it will cut the points into round pieces. DBSCAN does not ask for the number of groups. It looks for dense places, where a point has many neighbours, joins them in chains into clusters of any shape, and leaves lone shocks aside as grey noise. With a radius of 100 km and a threshold of ten points it gives 162 clusters, and 14 % of the earthquakes remain noise.

Move the radius. At 25 km more than half the earthquakes become noise, and only the densest swarms remain. At 400 km the clusters merge into a few giant belts. The switch “and within 30 days” makes neighbours only of shocks that were close both in place and in time. Then instead of plate boundaries you see swarms: series of aftershocks around strong earthquakes. The k-means button colours the same points into as many groups by nearest centre, and the difference from DBSCAN is visible at once: the belts are cut into pieces of roughly the same span, and lone shocks in mid-ocean get a group too.

The idea

A cluster is a place where points are dense, and between clusters it is empty. Density is measured in the simplest way: how many points lie in a circle of radius ε around a point. If many, the point is inside a cluster, and its neighbours are in the cluster too. If few, it is at the edge of a cluster or entirely alone.

There are two parameters: the radius ε and the threshold minPts. The neighbours of a point p are all points at a distance of at most ε from it, including itself. Denote them N(p). Each point gets one of three roles.

• A core is a point with |N(p)| ≥ minPts. Cores are the dense interior of clusters.
• A border point is not a core, but lies within ε of some core. A border point belongs to the cluster, but the cluster does not grow further through it.
• Noise is everything else: the point is not a core, and there is no core nearby.

A cluster is cores linked by chains in which every two neighbouring cores are neighbours, together with the border points of these cores. The algorithm gathers them like this. Take a point not yet looked at and count its neighbours. If there are fewer than minPts, the point is counted as noise for now: later a core may pick it up, and then it becomes a border point. If there are enough neighbours, a new cluster begins: all the neighbours join a queue. Each point from the queue is sent to this cluster, and if it is a core too, its neighbours join the queue as well. When the queue is empty, the cluster is ready, and we move on to the next point not looked at. In essence this is breadth-first search in the graph whose edges go from cores to their neighbours.

The problem about clusters on a map of epicentres rests on this theorem: the number of clusters is the number of connected components of the graph “core – neighbouring core”, and it can be gathered any way you like, even with a disjoint set union. If minPts = 2, any point with at least one neighbour becomes a core. Then the clusters of DBSCAN are the same pieces into which the network of edges shorter than ε falls apart, that is, single linkage from the article on the minimum spanning tree, and points without neighbours go to noise. A threshold minPts above two protects against chains: a thin bridge of sparse points does not join two clusters if the bridge's points are not cores.

Step by step

The steps cannot be made out on the whole world, so take one arc of islands in the north-west Pacific: 185 earthquakes in 2018–2021 on a flat map of this region. DBSCAN goes through them in time order. The circle around the current point is its radius ε. Cores are large circles, border points are circles with an empty middle, noise is grey dots.

With ε = 100 km and minPts = 4 there are 6 clusters: two long ones, of 67 and 59 points, stretched along the arc, another of 35 points, and three small ones of 4 points each. Of the 185 points 167 are cores, 6 border points and 12 noise. Increase ε to 120 km, and the northern cluster of 67 points merges with the middle one of 35 and one of the small groups into a cluster of 107 points: with the larger radius a chain of cores was found between them. Raise minPts to 5, and the small groups fall apart into noise, leaving three clusters. The “k-means” switch colours the same points into as many groups by nearest centre. It cuts the long arc into pieces of roughly the same span, it has no noise at all, and it gathers one of its six groups from scattered loners that DBSCAN counted as noise.

Code

Earthquakes lie on a sphere, and the distance between them is an arc along the Earth's surface. Latitude and longitude are awkward for distances: at the sixtieth parallel a degree of longitude is half as long as at the equator, and at the 180th meridian longitude jumps from 180 to −180. So each point is turned into a unit vector (x, y, z), as in the article on the k-d tree. An arc d on an Earth of radius R corresponds to the chord 2·sin(d/2R), and the chord grows with the arc. So “the arc is at most ε” is the same as “the chord is at most the chord of ε”, and ordinary distances in three-dimensional space can be compared.

A grid finds the neighbours: space is cut into cubes with a side equal to the chord of ε, and the points are put into cubes in a dictionary. A neighbour is at most one cube side away from the point in each coordinate, so it lies in the point's cube or in one of the 26 adjacent ones. The other cubes need not be looked at. This is the same idea as in a hash table: the cube's key tells at once where to look.

19073 earthquakes of magnitude 5+, 2015-2025
100 km, 10 points: 162 clusters, 15224 cores, 1177 border points, 2672 noise (14.0%); largest 2917, 2473, 2325; 2373291 distances instead of 181880128
100 km and 30 days, 10 points: 126 clusters, 3499 cores, 362 border points, 15212 noise (79.8%); largest 465, 313, 215; 2373291 distances instead of 181880128

Instead of 181.9 million pairs the grid computed 2.4 million distances, 77 times fewer. The three largest clusters are belts along deep-sea trenches, where one plate goes under another. The largest, 2917 earthquakes, stretches for 4779 km over the islands of the south-west Pacific. The second, 2473 earthquakes, is an arc of 7469 km along the north-west and north of the Pacific: several island arcs in a row, linked by cores into one chain. The third, 2325, is a narrow trench about 2500 km long south of the equator, on both sides of the 180th meridian.

Neighbours in time

The second line of the output is the same algorithm with a different notion of neighbour: two points are neighbours if they are at most 100 km and at most 30 days apart. DBSCAN does not care what distance is. It needs only an answer to the question “are these two points neighbours”, and any rule that answers it gives its own clustering. Here the clusters are no longer plate boundaries but swarms: a strong earthquake and a series of aftershocks after it, dense both in space and in time. 79.8 % of the points become noise: most strong earthquakes do not have even nine neighbours closer than 100 km and 30 days. The largest swarm is 465 shocks of magnitude 5 and above over five months of 2025 in the north-west Pacific; the strongest of them was of magnitude 8.8.

The grid does not change: neighbours in time are a subset of neighbours in place, so the cubes still cut off distant points, and time is checked as a second condition. The number of comparisons is the same, 2,373,291.

How to choose ε and minPts

The threshold minPts is easier to choose: it smooths the density estimate, and in noisy data it is raised. The authors of DBSCAN suggested 4 for points on a plane. In 1998 Jörg Sander and co-authors proposed twice the dimension of the data, 2·d, and for large and noisy data sets even more is advised. ε is harder. If there is knowledge of the subject, that is what is used. Erich Schubert and co-authors in a 2017 paper give this example: for GPS coordinates a specialist may decide that neighbours are points closer than a kilometre. Our 100 km is just a round number, and the graph below shows what it means for this data.

For each point find the distance to its (minPts − 1)-th nearest neighbour: this is its core distance. Minus one because the point itself is already the first of the minPts. The point is a core for exactly those ε that are at least this distance. Line the points up by decreasing core distance: this gives the sorted graph of distances to the neighbours, proposed already in the original DBSCAN paper. A horizontal line at height ε divides it into two parts. Everything below the line is cores, everything above is border points or noise.

The authors of DBSCAN advised looking for the first “valley” on this graph, where the curve bends sharply, and, if the share of noise in the data is known, marking it off on the left and taking ε at that place. On real data there is often no bend. There is none here either: on a logarithmic scale the curve descends smoothly. With minPts = 10, half the earthquakes have their ninth neighbour no further than 40 km, and the thousand loneliest further than 225 km. With ε = 100 km, 15,224 points become cores, 79.8 %, and the graph says the same at this place as the listing. Another 1177 points turn out to be border points, so 14 % go to noise.

In the same paper Schubert and co-authors name two signs of bad parameters. The first is the share of noise: a reasonable one is usually between 1 and 30 %. The second is the size of the largest cluster: if it holds more than 20–50 % of all points that fell into clusters, ε is most likely too large. On the earthquakes at 100 km the noise is 14 %, and the largest cluster is 17.8 % of the clustered points. At 150 km the largest cluster is already 35 %, at 300 km 64 %: the belts merge. At 50 km the noise is 31.6 %, slightly above the upper limit. Both signs narrow the choice to about 50–150 km. The same paper points out a detail that people often get wrong. Neighbours in nearest-neighbour search do not include the point itself, while in DBSCAN's density count they do. So the threshold minPts needs the (minPts − 1)-th neighbour. This off-by-one shift is the catch in the problem about the density threshold.

DBSCAN and k-means

Both algorithms divide points into groups, but they understand a group differently. k-means looks for k centres and gives each point to the nearest, so its groups are convex regions around centres, of roughly the same span. DBSCAN looks for connected dense regions, and their shape can be anything. On earthquakes this is the difference between a chain along a trench and a circle around a centre.

In the map scene k-means gets the same number of groups, 162, and starts with k-means++. The diameter of a typical k-means group (the median) is 1114 km, of a typical DBSCAN cluster 211 km. But the longest DBSCAN cluster stretches for 7469 km, and the longest k-means group for 2542 km: a belt that DBSCAN sees whole, k-means cuts into pieces. It splits 19 DBSCAN clusters into several sizeable parts, and 52 of its groups mix points of different DBSCAN clusters. And k-means has no noise: it handed all 2672 lone shocks to the nearest groups. Their median distance to the centre of their group is 375 km, and for 64 of them it is over a thousand.

This does not mean DBSCAN is better. k-means answers a different question: how to replace the points with k representatives with the least error. For the palette of a photograph that is just what is needed. DBSCAN answers the question of where points lie densely. It does not need to know the number of groups in advance, and it can say “this point belongs to nothing”. But it has two parameters to which the result is sensitive, while k-means has one.

What it costs

DBSCAN asks for the neighbours of each point exactly once: a point leaves the “not looked at” ones as soon as its neighbours have been counted. So the running time is n neighbour queries plus the time to go through all the neighbours found: O(n·Q + Σ|N(p)|), where Q is the cost of one query. Without an index a query is a scan of all points, Q = O(n), and everything costs O(n²). With a grid a query looks at only 27 cubes, and on the earthquakes that is 77 times fewer distances. Instead of a grid a k-d tree with a search for all points within a radius will do too. In scikit-learn DBSCAN offers a choice: a k-d tree, a ball tree or brute force.

The scan of pairs grows as n², as it should. The grid in Python is dozens of times faster than the scan, but its line also rises more steeply than direct proportion. This is not a mistake. The earthquakes are taken in time order, and over twice the span twice as many neighbours accumulate within the same 100 km. The term Σ|N(p)| grows as a square, and no index removes it: the neighbours found still have to be gone through. The grid removes only comparisons with distant points. In Go the gap is smaller, 5–16 times: there the scan of pairs is a tight loop, while the grid pays for looking up cubes in a dictionary.

The worst case stays quadratic with any index: if all points lie within ε of each other, each query returns almost all of them. The original 1996 paper estimated the time of DBSCAN with an R*-tree as O(n log n) on average, counting on small neighbourhoods. In 2013 Ade Gunawan noticed in a master's thesis that in the worst case the algorithm is quadratic, and for the plane built an exact algorithm in O(n log n) on a grid with cells of ε/√2. Any two points of such a cell are neighbours, and if a cell has gathered minPts points, all of them are cores without a single comparison. In 2015 Junhao Gan and Yufei Tao showed that in three or more dimensions exact DBSCAN apparently needs time of the order of n^4/3. This bound can be got around only by solving one of the long-standing hard problems of computational geometry faster, and few believe that possible. Instead they proposed an approximate DBSCAN in O(n) on average. Two years later Schubert, Sander, Ester, Kriegel and Xu replied that the fault lay not with DBSCAN but with the hope placed in spatial indexes. On one of the data sets from that work, at the smallest ε used in Gan and Tao's experiments, more than 99 % of the points already lie in one cluster. With sensibly chosen ε, ordinary DBSCAN with a good index works no worse than the new approximate method.

Where the algorithm goes wrong

• One density for all. The radius ε is one for all the data. If there is a dense swarm next to a sparse belt, an ε chosen for the belt will merge the swarm with its neighbours, and one chosen for the swarm will turn the belt into noise. This shows on the map: at 50 km the densest trenches still hold together as whole clusters, while the mid-ocean ridges have gone almost entirely to noise. This weakness is treated by OPTICS (1999), which in one pass orders the points so that clusters for all ε up to a given limit can be read from the order, and HDBSCAN (2013), which chooses clusters from a whole hierarchy of densities.
• Chains of cores. If two dense places are joined by a string of cores, they become one cluster, although to the eye they are two. This is how the second-largest earthquake cluster, the arc of 7469 km, was gathered. minPts protects against bridges of sparse points, but not against a bridge of cores.
• Border points depend on order. A border point with cores of two clusters nearby goes to whichever reached it first. So a cluster can even turn out smaller than minPts: some neighbours of its first core were already taken by a cluster gathered earlier. In PostGIS the documentation of the function ST_ClusterDBSCAN warns outright that such border points are assigned arbitrarily, and advises setting the order of traversal.
• Many dimensions. In a space with dozens of features the distances between any points become similar, and a “dense place” loses its meaning. Before DBSCAN the dimension is often reduced, for example by principal component analysis, and a grid is useless in high dimensions: there are 3^d adjacent cells.
• The scale of features. Distance mixes all coordinates, and a feature in large units drowns out the others. In neighbours “in place and time” this problem is solved by two separate thresholds, 100 km and 30 days, rather than one formula like “kilometres plus days”.

History

DBSCAN was proposed by Martin Ester, Hans-Peter Kriegel, Jörg Sander and Xiaowei Xu at the KDD conference in 1996. The paper was called “A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise”, and the algorithm was tested on points of the plane, synthetic ones and real geographic ones from the SEQUOIA 2000 data set. In 1998 Sander, Ester, Kriegel and Xu generalised it to GDBSCAN: there neighbourhood is given by any symmetric rule “neighbours or not”, and not only points can be clustered but also extended objects, such as polygons on a map. In 1999 Ankerst, Breunig, Kriegel and Sander described OPTICS, and in 2013 Campello, Moulavi and Sander HDBSCAN. In 2014 the original paper received the SIGKDD Test of Time Award for its lasting influence. DBSCAN is in scikit-learn, ELKI, R and Weka, and since version 2.3 also in PostGIS, right in database queries.

Problems

The first problem is DBSCAN itself on the plane with a grid: how many clusters, cores, border points and noise. The second is about choosing ε: how to say at once how many cores there will be at any radius.