Faster DB-scan and HDB-scan in Low-Dimensional Euclidean Spaces

TL;DR

This paper introduces a new $O(n\,log n)$ algorithm for DBscan in 2D that is faster and less sensitive to parameters, along with an efficient approach for HDBscan and its approximation in fixed dimensions.

Contribution

It presents the first $O(n\,log n)$ algorithms for DBscan and HDBscan in low-dimensional Euclidean spaces, improving efficiency and robustness.

Findings

The new DBscan algorithm is theoretically optimal in 2D.

The simplified algorithm performs well in practice.

The HDBscan algorithm is also optimized to $O(n\,log n)$ in the plane.

Abstract

We present a new algorithm for the widely used density-based clustering method DBscan. Our algorithm computes the DBscan-clustering in $O(n\log n)$ time in $\mathbb{R}^2$, irrespective of the scale parameter $\varepsilon$ (and assuming the second parameter MinPts is set to a fixed constant, as is the case in practice). Experiments show that the new algorithm is not only fast in theory, but that a slightly simplified version is competitive in practice and much less sensitive to the choice of $\varepsilon$ than the original DBscan algorithm. We also present an $O(n\log n)$ randomized algorithm for HDBscan in the plane---HDBscan is a hierarchical version of DBscan introduced recently---and we show how to compute an approximate version of HDBscan in near-linear time in any fixed dimension.

Equations

$$E=\{(p,q)\in D\times D:q\in N_{\epsilon }(p,D)\setminus \{p\}\}.$$

$$O\left(\sum _{j=1}^m\sum _{j^{\prime }=\max (j-2,1)}^{\min (j+2,m)}\left(|{\mathcal{B}}_j|+|{\mathcal{B}}_{j^{\prime }}|\right)\right)=O\left(\sum _{j=1}^m|{\mathcal{B}}_j|\right)=O(n).$$

$$O\left(n_b\cdot \sum _{b^{\prime }\in N_{\epsilon }(b,\mathcal{B})}{n}_{b^{\prime }}\right)=O\left(\text{mbox}\text{sc}MinPts\cdot \sum _{b^{\prime }\in N_{\epsilon }(b,\mathcal{B})}{n}_{b^{\prime }}\right).$$

$$\sum _{b\in {\mathcal{B}}^{\ast }}\sum _{b^{\prime }\in N_{\epsilon }(b,{\mathcal{B}}^{\ast })}O((n_b+n_{b^{\prime }})\log (n_b+n_{b^{\prime }})),$$

$$\sum _{b\in {\mathcal{B}}^{\ast }}\sum _{b^{\prime }\in N_{\epsilon }(b,{\mathcal{B}}^{\ast })}O((1/\delta {)}^{d/3}\cdot (n_b+n_{b^{\prime }}))=O((1/\delta {)}^{d/3}n).$$

$$\begin{array}[]{lll}|rr_{i}|&\leqslant&|rr^{*}|+|r^{*}r_{i}^{*}|+|r_{i}r_{i}^{*}|\\[5.69054pt] &\leqslant&\sin\theta\cdot|pr|+|r^{*}q^{*}|+\sin\theta\cdot|pr_{i}|\\[5.69054pt] &\leqslant&2\sin\theta\cdot|pq|/\cos\theta+|r^{*}q^{*}|\\[5.69054pt] &=&2\sin\theta\cdot|pq|/\cos\theta+|pq^{*}|-|pr^{*}|\\[5.69054pt] &\leqslant&2\sin\theta\cdot|pq|/\cos\theta+|pq|-|pr|\cos\theta\\[5.69054pt] &\leqslant&\left(2\frac{\sin\theta}{\cos\theta}+1-\cos\theta\right)\cdot(1-\delta)\varepsilon\end{array}$$

Full text

1 Introduction

Clustering is one of the most fundamental tasks in data mining. Due to the wide variety of applications where clustering is important, the clustering problem comes in many variants. These variants differ for example in the dimensionality of the data set $D$ and in the underlying metric, but also in the objective of the clustering. Thus a multitude of clustering algorithms has been developed [20], each with their own strengths and weaknesses. We are interested in density-based clustering, where clusters are defined by areas in which the density of the data points is high and clusters are separated from each other by areas of low density.

One of the most popular density-based clustering methods is dbscan; the paper by Ester et al. [11] on dbscan has been cited over 8,800 times, and in 2014 dbscan received the test-of-time award from KDD, a leading data-mining conference. dbscan has two parameters, $\varepsilon$ and MinPts, that together determine when the density around a point $p\in D$ is high enough for $p$ to be part of a cluster (as apposed to being noise); see Section 2 for a precise definition of the dbscan clustering. Typically MinPts is a constant—in the original article [11] it is concluded that $\mbox{{\sc MinPts}}=4$ works well—but finding the right value for $\varepsilon$ is more difficult. The worst-case running time of the original dbscan algorithm is $\Theta(n^{2})$. It is often stated that the running time is $O(n\log n)$ for Euclidean spaces when a suitable indexing structure such as an R-tree is used to support the dbscan algorithm. While this may be true in certain practical cases, it is not true from a theoretical point of view.

Several variants of dbscan algorithm have been proposed, often with the goal to speed up the computation. Some (idbscan [5] and fdbscan [15]) do so at the expense of computing a slightly different, and not clearly defined, clustering. Others (gridbscan [16]) compute the same clustering as dbscan, but without speeding up the worst-case running time.

A fundamental bottleneck of the original dbscan algorithm is that it performs a query with each point $p\in D$ to find $N_{\varepsilon}(p,D)$, the set of points within distance $\varepsilon$ of $p$. Thus $\sum_{p\in D}|N_{\varepsilon}(p,D)|$ is a lower bound on the running time of the dbscan algorithm. In the worst case $\sum_{p\in D}|N_{\varepsilon}(p,D)|=\Theta(n^{2})$, so even with a fast indexing structure the worst-case running time of the original dbscan algorithm is $\Omega(n^{2})$. (Apart from this, the worst-case query time of R-trees and other standard indexing structures is not logarithmic even if we disregard to time to report points.) In most practical instances the dbscan algorithm is much faster than quadratic. The reason is that $\varepsilon$ is typically small so that the sets $N_{\varepsilon}(p,D)$ do not contain many points and the range queries can be answered quickly. However, the fact that the algorithm always explicitly reports the sets $N_{\varepsilon}(p,D)$ makes the running time sensitive to the choice of $\varepsilon$ and the density of the point set $D$. For example, suppose we have a disk-shaped cluster with a Gaussian distribution around the disk center. Then a suitable value of $\varepsilon$ will lead to large sets $N_{\varepsilon}(p,D)$ for points $p$ near the center of the cluster.

Chen et al. [8] overcame the quadratic bottleneck of the standard approach, and designed an algorithm111As described, the algorithm actually computes a variation of the dbscan clustering, but it is easily adapted to compute the true dbscan clustering. with $O(n^{2-\frac{2}{d+2}}\;\mbox{polylog}\;n)$ worst-case running time. Note that for $d=2$ the running time of the exact algorithm is $O(n^{1.5}\;\mbox{polylog}\;n)$. They also present an approximate algorithm that is more practical. Chen *et al.*remark that their exact algorithm is mainly of theoretical interest. The natural question is then whether or not it is possible to to compute the dbscan clustering in subquadratic time in the worst case, irrespective of the value of $\varepsilon$, with a simple and practical algorithm?

Although dbscan is used extensively and performs well in many situations, it has its drawbacks. One is that it produces a flat (non-hierarchial) clustering which heavily depends on the choice of the scale parameter $\varepsilon$. Ankerst et al. [3] therefore introduced optics, which can be seen as a hierarchical version of dbscan. Recently Campello et al. [7] proposed an improved density-based hierarchical clustering method—similar to optics but cleaner—together with a cluster-stability measure that can be used to automatically extract relevant clusters. The new method, called hdbscan, only needs the parameter MinPts, which is much easier to choose than $\varepsilon$. (Campello *et al.*used MinPts=4 in all their experiments.) While hdbscan is very powerful, the algorithm to compute the hdbscan hierarchy runs in quadratic time; not only in the worst-case, but actually also in the best-case. There have been only few papers dealing with speeding up hdbscan or its predecessor optics. A notable recent exception is Poptics [19], a parallel algorithm that computes a similar (though not the same) hierarchy as optics. We do not know of any algorithm that computes the hdbscan or optics hierarchy in subquadratic time. Thus the second question we study is: is it possible to compute the hdbscan hierarchy in subquadratic time?

Our results. We present an $O(n\log n)$ algorithm to compute the dbscan clustering for a set $D$ of $n$ points in the plane, irrespective of the setting of the parameter $\varepsilon$ used to define the dbscan clustering. (Here, and in our other results, we assume that the parameter MinPts is a fixed constant. As mentioned this is the case in practice, where one typically uses $\mbox{{\sc MinPts}}=4$.) We remark that our algorithm is not only fast in theory, but a slightly simplified version is also competitive in practice and much less sensitive to the choice of $\varepsilon$ than the original dbscan algorithm. Some basic experimental results are provided in Section 6.

We also present a new algorithm for planar hdbscan: we show how to compute the hdbscan hierarchy in ${\mathbb{R}}^{2}$ in $O(n\log n)$ expected time, thus obtaining the first subquadratic algorithm for the problem.

Finally, we provide a slightly improved version of the approximate dbscan clustering algorithm by Chen et al. [8] and by Gan and Tao [12] (their results are discussed in more detail below). Specifically we improve the dependency on the approximation parameter $\delta$. We then extend the concept of an approximate dbscan clustering as defined to the hierarchical version. We thus obtain $\delta$-approximate hdbscan, an approximate version of the hdbscan hierarchy of Campello et al. [7], where the parameter $\delta$ specifies the accuracy of the approximation. (Intuitively, a $\delta$-approximate hdbscan hierarchy has the same clusters as the standard hdbscan hierarchy at any level $\varepsilon$, except that clusters at distance $(1-\delta)\cdot\varepsilon$ from each other may be merged. See Section 5 for precise definition.) We show that a $\delta$-approximate hdbscan hierarchy in ${\mathbb{R}}^{d}$ can be computed in $O((n/\delta^{(d-1)/2})\log n)$ time.

Further related work. This work should be viewed as the journal publication of the (so far unpublished) masters thesis of the second author [14], which contained the results on dbscan, extended with results on hdbscan. In the meantime, Gan and Tao [12] published a paper in which they extend the work from the masters thesis to ${\mathbb{R}}^{d}$, resulting in an algorithm for dbscan with a running time of $O(n^{2-\frac{2}{\lceil d/2\rceil+1}+\gamma})$; we briefly comment on how this is done at the end of Section 3. Gan and Tao also prove that computing the dbscan clustering in ${\mathbb{R}}^{d}$ for $d\geqslant 3$ is at least as hard as the so-called unit-spherical emptiness problem, which is believed to require $\Omega(n^{4/3})$ time [10]. Finally, Gan and Tao show that a $\delta$-approximate dbscan clustering can be computed in $O(n/\delta^{d-1})$ expected time, using a modified version of the exact algorithm. Their approximate clustering is the same as the approximate clustering defined by Chen et al. [8], who already showed how to compute it in $O(n\log n+n/\delta^{d-1})$ time deterministically. (Gan and Tao were unaware of the paper by Chen et al..) As we show in Section 5 our algorithm can also be used to obtain a deterministic algorithm with $O(n\log n+n/\delta^{d/3})$ running time.

2 Preliminaries on dbscan and dbscan∗

Let $D$ be a set of points in ${\mathbb{R}}^{d}$. dbscan distinguishes three types of points: core points (points in the “interior” of a cluster), border points (points on the boundary of a cluster), and noise (points not in any cluster). The distinction is based on two global parameters, $\varepsilon$ and MinPts. Define $N_{\varepsilon}(p,D):=\{q\in D:|pq|\leqslant\varepsilon\}$ to be the neighborhood of a point $p$, where $|pq|$ denotes the (Euclidean) distance between $p$ and $q$; the point $p$ itself is included in $N_{\varepsilon}(p,D)$. A point $p\in D$ is a core point if $|N_{\varepsilon}(p,D)|\geqslant\mbox{{\sc MinPts}}$, and a non-core point $q$ in the neighborhood of a core point is a border point. We denote the set of core points by $D_{\mathrm{core}}$, and the set of border points by $D_{\mathrm{border}}$. The remaining points are noise. In dbscan∗ [7] border points are not part of a cluster but are considered noise.

Ester et al. [11] define the dbscan clusters based on the concept of density-reachability (a detailed description is given below). Equivalently, we can define the clusters as the connected components of a certain graph.

To this end, define the neighborhood graph $\mathcal{G}(D,E)$ as the (undirected) graph with node set $D$ and edges connecting pairs of points within distance $\varepsilon$; see Fig. 1. In other words,

[TABLE]

Note that a point $p\in D$ is a core point if and only if its degree in $\mathcal{G}$ is at least $\mbox{{\sc MinPts}}-1$, since then its neighborhood contains at least MinPts points (including $p$ itself). Now consider the subgraph $\mathcal{G}_{\mathrm{core}}(D_{\mathrm{core}},E_{\mathrm{core}})$ induced by the core points, that is, $\mathcal{G}_{\mathrm{core}}$ is the graph whose nodes are the core points and whose edges connect two core points when they are within distance $\varepsilon$ from each other. We call $\mathcal{G}_{\mathrm{core}}$ the core graph. The connected components of $\mathcal{G}_{\mathrm{core}}$ are the clusters in dbscan∗. The clusters in dbscan are the same, except that they also contain border points. Formally, a border point $q$ belongs to a cluster $C$ if $q$ has an edge (in $\mathcal{G}$) to a core point $p\in C$. Thus a border point can belong to multiple clusters. The original dbscan algorithm assigns a border point $p$ to the first cluster that finds $p$ (clusters are constructed one by one); we assign border points to the cluster of their nearest core point.

2.1 The original definition of the dbscan clustering.

For comparison purposes only, we restate the original definition of the dbscan clustering. Ester et al. [11] define when two points are in the same cluster based on the concept of density-reachability, as explained next. A point $q\in D$ is directly density-reachable from a point $p\in D$ if $q\in N_{\varepsilon}(p,D)$ and $p$ is a core point. We denote this by $p\rightarrow q$. A core point is always directly density-reachable from itself, since a point $p\in D$ is always in its own neighborhood. A point $q$ is density-reachable from a core point $p$, denoted by $p\leadsto q$, if there is a sequence $p=r_{0},\ldots,r_{k}=q$ (for some $k\geqslant 0$) such that $r_{0}\rightarrow r_{1}\rightarrow\cdots\rightarrow r_{k}$. Observe that if two points $p$ and $q$ are both core points, then $p\leadsto q$ if and only if $q\leadsto p$. Two points $p$ and $q$ are density-connected if there is a point $r$ such that $r\leadsto p$ and $r\leadsto q$.

A cluster is now defined as a subset $C\subseteq D$ such that (i) if a core point $p$ is in $C$ then all points $q$ that are reachable from $p$ are in $C$, and (ii) any two points in $C$ are density-connected to each other. As observed by Ester et al., a cluster must contain at least one core point, and a cluster is uniquely defined by any of its core points. More precisely, if $p\in C$ is a core point then $C=\{q\in D:p\leadsto q\}$. Each core point belongs to exactly one dbscan cluster. Under the above definition border points can belong to multiple clusters, however. This is typically undesirable, so the original dbscan algorithm assigns each border point $q$ to only one cluster, namely the cluster from which $q$ is discovered first by their algorithm. This implies that the computed clustering depends on the order in which their algorithm happens to handle the points.

3 A fast algorithm for dbscan

The original dbscan algorithm reports, while generating and exploring the clusters, for each point $p\in D$ all its neighbors. In other words, it spends time on every edge in the neighborhood graph. Our new algorithm avoids this by working with a smaller graph, the box graph $\mathcal{G}_{\mathrm{box}}$. Its nodes are disjoint rectangular boxes with a diameter of at most $\varepsilon$ that together contain all the points in $D$, and its edges connect pairs of boxes within distance $\varepsilon$; see Fig. 2.

The boxes are generated such that (i) any two points in the same box are in each other’s neighborhood, and (ii) the degree of any node in the box graph is $O(1)$. Property (i) allows us to immediately classify all points in a box as core points when it contains at least MinPts points, and property (ii) allows us to quickly retrieve the neighbors of any given point in a box. Next we describe the algorithm, which consists of four easy steps, in detail.

Step 1: Compute the box graph $\mathcal{G}_{\mathrm{box}}$. To compute $\mathcal{G}_{\mathrm{box}}$, we first construct a collection of vertical strips that together cover all the points. Let $p_{1},\ldots,p_{n}$ be the points in $D$ sorted by $x$-coordinate, with ties broken arbitrarily. The first strip has $p_{1}$ on its left boundary. We go through the remaining points from left to right, adding them to the first strip as we go, until we encounter a point $p_{i}$ whose distance to the left strip boundary is more than $\varepsilon/\sqrt{2}$. We then start a new strip with $p_{i}$ on its left boundary, and we add points to that strip until we encounter a point whose distance to the left strip boundary is more than $\varepsilon/\sqrt{2}$, and so on, until we have handled all the points. Constructing the strips takes $O(n)$ time, after sorting the points by $x$-coordinate.

Within each strip we perform a similar procedure, going over the points within the strip in order of increasing $y$-coordinate and creating boxes instead of strips. Thus the first box in the strip has the lowest point on its bottom edge, and we keep adding points to this box (enlarging it so that the new point fits, ensuring a tight bounding box) until we encounter a point whose vertical distance to the bottom edge is more than $\varepsilon/\sqrt{2}$. We then start a new box, and so on, until we handled all points in the strip. If the number of points in the $j$-th strip is $n_{j}$, then the time needed to handle all the strips is $\sum_{j}O(n_{j}\log n_{j})=O(n\log n)$.

Let $m$ be the number of strips and $\mathcal{B}_{j}$ the set of boxes in the $j$-th strip. We sometimes refer to a set $\mathcal{B}_{j}$ as a strip, even though formally $\mathcal{B}_{j}$ is a set of boxes. Let $\mathcal{B}:=\mathcal{B}_{1}\cup\cdots\cup\mathcal{B}_{m}$. The nodes of the box graph $\mathcal{G}_{\mathrm{box}}$ are the boxes in $\mathcal{B}$ and there is an edge $(b,b^{\prime})$ when $\mathrm{dist}(b,b^{\prime})\leqslant\varepsilon$, where $\mathrm{dist}(b,b^{\prime})$ denote the minimum distance between $b$ and $b^{\prime}$. Two boxes $b,b^{\prime}$ are neighbors when they are connected by an edge in $\mathcal{G}_{\mathrm{box}}$. Let $N_{\varepsilon}(b,\mathcal{B})$ be the set of neighbors $b$.

Lemma 3.1

$\mathcal{G}_{\mathrm{box}}$* has at most $n$ nodes, each having $O(1)$ neighbors.*

Proof 3.2

The box graph obviously has at most $n$ nodes. Next we give a precise analysis of the number of neighbors a box $b$ in some strip $B_{j}$ can have. Consider the node corresponding to the box $b$. If there are two or more strips in between $b$ and some other box $b^{\prime}$ then $\mathrm{dist}(b,b^{\prime})>2(\varepsilon/\sqrt{2})>\varepsilon$, so $b$ can only have neighbors in $\mathcal{B}_{j-2}$, $\mathcal{B}_{j-1}$, $\mathcal{B}_{j}$, $\mathcal{B}_{j+1}$, or $\mathcal{B}_{j+2}$. We bound the number of neighbors in each of these strips separately.

• •

The number of neighbors in $\mathcal{B}_{j}$ is at most four. Indeed, the boxes in $\mathcal{B}_{j}$ can be ordered vertically, and if there are more than two boxes in between $b$ and some other box $b^{\prime}\in\mathcal{B}_{j}$, then the vertical distance between $b$ and $b^{\prime}$ is more than $\varepsilon$.

• •

*The number of neighbors in $\mathcal{B}_{j-1}$ (and similarly in $\mathcal{B}_{j+1}$) is at most five. Suppose for a contradiction that $b$ has six or more neighbors in $B_{j-1}$. Let $b^{\prime}$ and $b^{\prime\prime}$ be the lowest and highest of these neighbors, respectively. Then the vertical distance between the top boundary of $b^{\prime}$ and the bottom boundary of *$b^{\prime\prime}$is more than $4\varepsilon/\sqrt{2}$. Since the height of $b$ is at most $\varepsilon/\sqrt{2}$, this implies that the vertical distance between $b$ and either $b^{\prime}$ or $b^{\prime\prime}$ is more than $(4\varepsilon/\sqrt{2}-\varepsilon/\sqrt{2})/2>\varepsilon$, contradicting that both $b^{\prime}$ and $b^{\prime\prime}$ are neighbors of $b$.

• •

The number of neighbors in $\mathcal{B}_{j-2}$ (and similarly in $\mathcal{B}_{j+2}$) is at most four. Suppose for a contradiction that $b$ has five or more neighbors in $B_{j-2}$. Let $b^{\prime}$ and $b^{\prime\prime}$ be the lowest and highest of these neighbors, respectively. Then the vertical distance between the top boundary of $b^{\prime}$ and the bottom boundary of $b^{\prime\prime}$ is more than $3\varepsilon/\sqrt{2}$. Since the height of $b$ is at most $\varepsilon/\sqrt{2}$, this implies the vertical distance between $b$ and either $b^{\prime}$ or $b^{\prime\prime}$ is more than $(3\varepsilon/\sqrt{2}-\varepsilon/\sqrt{2})/2=\varepsilon/\sqrt{2}$. The horizontal distance is at least $\varepsilon/\sqrt{2}$ as well, because there is a strip in between. Hence, the total distance between $b$ and either $b^{\prime}$ or $b^{\prime\prime}$ is more than $\varepsilon$, contradicting that both $b^{\prime}$ and $b^{\prime\prime}$ are neighbors of $b$.

Adding up the maximum number of neighbors in each of the strips gives us a maximum of 22 neighbors in total.

This also gives us an easy way to compute the edge set $E_{\mathrm{box}}$ of the box graph, because the edges between boxes in strips $\mathcal{B}_{j}$ and $\mathcal{B}_{j^{\prime}}$ with $|j-j^{\prime}|\leqslant 2$ can be computed in $O(|\mathcal{B}_{j}|+|\mathcal{B}_{j^{\prime}}|)$ time in total by scanning the boxes in $\mathcal{B}_{j}$ and $\mathcal{B}_{j^{\prime}}$ in a coordinated manner. The total time to compute all edges of the box graph is thus

[TABLE]

Adding the time to construct the strips and boxes, we see that Step 1 takes $O(n\log n)$ time and we obtain the following lemma.

Lemma 3.3

The box graph $\mathcal{G}_{\mathrm{box}}(\mathcal{B},E_{\mathrm{box}})$ can be computed in $O(n\log n)$ time.

An alternative for Step 1. An alternative approach is to define the boxes as the non-empty cells in a grid whose cells have height and width $\varepsilon/\sqrt{2}$. If we store the boxes in a hash-table based on the coordinates of their lower left corners, then finding the neighbors of a box $b$ can be done by checking each potential neighbor cell for existence in the hash-table—we do not need to store the box graph explicitly. Creating the boxes (with their corresponding point sets) can be done in $O(n)$ time if the floor function can be computed in $O(1)$ time.

Step 2: Find the core points. The graph $\mathcal{G}_{\mathrm{box}}$ allows us to determine the core points in a simple and efficient manner. The key observation is that the maximum distance between any two points in the same box is at most $\varepsilon$. Hence, if a box contains more than MinPts points, then all of them are core points. Thus the following simple algorithm suffices to determine the core points.

For a box $b\in\mathcal{B}$, let $D(b):=D\cap b$ be the set of point inside $b$, and let $n_{b}:=|D(b)|$. If $n_{b}\geqslant\mbox{{\sc MinPts}}$ then label all points in $b$ as core points. Otherwise, for each point $p\in D(b)$, count the number of points $q$ in neighboring boxes of $b$ for which $|pq|\leqslant\varepsilon$. If this number is at least $\mbox{{\sc MinPts}}-n_{b}$, then label $p$ as core point. The counting is done brute-force, by checking all points in neighboring boxes. Hence, this takes $O(\sum_{b^{\prime}\in N_{\varepsilon}(b,\mathcal{B})}n_{b^{\prime}})$ time for each point $p\in b$.

Lemma 3.4

Given $\mathcal{G}_{\mathrm{box}}$, we can find all core points in $D$ in $O(n)$ time.

Proof 3.5

The total time spent to handle boxes $b$ with $n_{b}\geqslant\mbox{{\sc MinPts}}$ is clearly $O(n)$. The time needed to handle a box $b$ with $n_{b}<\mbox{{\sc MinPts}}$ is

[TABLE]

Now charge $O(\mbox{{\sc MinPts}})=O(1)$ time to each point in every $b^{\prime}\in N_{\varepsilon}(b,\mathcal{B})$. Because any box $b^{\prime}$ is the neighbor of $O(1)$ other boxes by Lemma 3.1, each point is charged $O(1)$ times, proving the lemma.

Step 3: Compute the cluster cores. In the previous step, we determined the core points. Next we wish to determine the clusters or, more precisely the cluster cores. The core of a cluster is the set of core points in that cluster. In Step 3 we assign to each core point a cluster-id so that core points in the same cluster have the same cluster-id. Again, this can be done in an efficient manner using $\mathcal{G}_{\mathrm{box}}$. To this end, we first remove certain boxes and edges from $\mathcal{G}_{\mathrm{box}}$ to obtain a reduced box graph $\mathcal{G}_{\mathrm{box}}^{*}$. More precisely, we keep only the boxes with at least one core point, and we keep only the edges $(b,b^{\prime})$ for which there are core points $p\in b$, $p^{\prime}\in b^{\prime}$ with $|pp^{\prime}|\leqslant\varepsilon$. Because any two core points in a given box $b$ are connected in $\mathcal{G}_{\mathrm{core}}$, we have the following lemma.

Lemma 3.6

The connected components in $\mathcal{G}_{\mathrm{box}}^{*}$ correspond one-to-one to the connected components in the core graph $\mathcal{G}_{\mathrm{core}}$ and, hence, to the dbscan∗ clusters.

Thus the cluster cores can be computed by computing the connected components in $\mathcal{G}_{\mathrm{box}}^{*}$. The latter can be done in $O(n)$ time using DFS [9]. After computing the connected components we give every core point $p$ a cluster-id corresponding to the connected component of the box $b$ that contains $p$.

To construct $\mathcal{G}_{\mathrm{box}}^{*}$, we need to decide for two given boxes $b,b^{\prime}$ whether there are core points $p\in D(b)$, $p^{\prime}\in D(b^{\prime})$ with $|pp^{\prime}|\leqslant\varepsilon$. If $n_{b}<\mbox{{\sc MinPts}}$ or $n_{b^{\prime}}<\mbox{{\sc MinPts}}$ then we simply check all pairs of core points. We can argue as in the proof of Lemma 3.4 that this takes $O(\mbox{{\sc MinPts}}\cdot n)=O(n)$ time in total. If $n_{b}\geqslant\mbox{{\sc MinPts}}$ and $n_{b^{\prime}}\geqslant\mbox{{\sc MinPts}}$ we have to be more careful,since checking all pairs of points can lead to a quadratic running time. For example, if both cells contain $n/2$ points, then checking all pairs would take $\Omega(n^{2})$ time. When $n_{b}\geqslant\mbox{{\sc MinPts}}$ and $n_{b^{\prime}}\geqslant\mbox{{\sc MinPts}}$ we therefore compute the Delaunay triangulation of $D(b)\cup D(b^{\prime})$ in $O((n_{b}+n_{b^{\prime}})\log(n_{b}+n_{b^{\prime}}))$ time [6], and check whether it has an edge $(p,p^{\prime})$ of length at most $\varepsilon$ between points $p\in D(b),p^{\prime}\in D(b^{\prime})$. This is sufficient because the pair of points defining the closest distance between $D(b)$ and $D(b^{\prime})$ must define an edge in the Delaunay triangulation. This leads to the following lemma.

Lemma 3.7

Computing the cluster cores can be done in $O(n\log n)$ time.

Proof 3.8

The most time consuming part of the construction of $\mathcal{G}_{\mathrm{box}}^{*}$ is to determine for each pair of neighboring boxes in $\mathcal{B}$ whether there are core points $p\in b$, $p^{\prime}\in b^{\prime}$ with $|pp^{\prime}|\leqslant\varepsilon$. As mentioned, the total time spent on pairs with $n_{b}<\mbox{{\sc MinPts}}$ or $n_{b^{\prime}}<\mbox{{\sc MinPts}}$ is $O(n)$. Let $\mathcal{B}^{*}$ be the set of boxes containing at least MinPts points. Then the total time spent on the pairs of boxes from $\mathcal{B}^{*}$ is

[TABLE]

which is $O(n\log n)$ because $|N_{\varepsilon}(b,\mathcal{B}^{*})|=O(1)$ for any box $b$ and $\sum_{b\in\mathcal{B}^{*}}n_{b}\leqslant n$.

Remark. In practice we can also use a brute-force algorithm when $n_{b}\geqslant\mbox{{\sc MinPts}}$ and $n_{b^{\prime}}\geqslant\mbox{{\sc MinPts}}$, because the number of points in boxes with more than MinPts points is typically still not very large. Moreover, if both $b$ and $b^{\prime}$ contain many points, then there are often many pairs of points within distance $\varepsilon$ from each other, and we can stop when we find such a pair.

Step 4: Assigning border points to clusters. It remains to decide for non-core points $p$ whether $p$ is a border point or noise. (For dbscan∗ we can skip Step 4, since in dbscan∗ border points are considered noise.) If $p$ is a border point, it has to be assigned to the nearest cluster. Again, a brute-force method suffices: for each box $b\in B$ and each non-core point $p\in b$, we check all points in $b$ and its neighboring boxes to find $p$’s nearest core point, $p^{\prime}$. If $|pp^{\prime}|\leqslant\varepsilon$, then $p$ is a border point in the same cluster as $p^{\prime}$, otherwise $p$ is noise. We only need to consider boxes $b$ with $n_{b}<\mbox{{\sc MinPts}}$—otherwise all points in $b$ are core points—so the argument from the proof of Lemma 3.4 shows that this takes $O(n)$ time.

Putting it all together. Steps 1 and 3 take $O(n\log n)$ time and Steps 2 and 4 take $O(n)$ time. We thus obtain the following theorem.

Theorem 3.9

Let $D$ be a set of $n$ points in ${\mathbb{R}}^{2}$, and $\varepsilon$ and MinPts be given constants. Then we can compute a dbscan clustering on $D$ according to $\varepsilon$ and MinPts for the Euclidean metric in $O(n\log n)$ time.

Remark: extension to higher dimensions. The algorithm just described can easily be extended to ${\mathbb{R}}^{d}$ for $d>2$, as already observed by Gan and Tao [12]. For completeness we describe the extension and the resulting bound. One trivial modification is that in ${\mathbb{R}}^{d}$ we need to use boxes in $\mathcal{G}_{\mathrm{box}}$ of width at most $\varepsilon/\sqrt{d}$ along each axis to ensure their diameter is at most $\varepsilon$. The only other difference is in Step 3, where we have to decide for all pairs $b,b^{\prime}$ of neighboring boxes whether there are core points $p\in D(b)$ and $p^{\prime}\in D(b^{\prime})$ with $|pp^{\prime}|\leqslant\varepsilon$. When $n_{b}\geqslant\mbox{{\sc MinPts}}$ and $n_{b^{\prime}}\geqslant\mbox{{\sc MinPts}}$ we can no longer use the Delaunay triangulation for this, as it may have quadratic size. Instead we can use a known algorithm for bichromatic closest pair [2], which gives a running time of $O(n^{2-\frac{2}{\lceil d/2\rceil+1}+\gamma})$, where $\gamma>0$ is an arbitrarily small constant. (For $d=3$ the $n^{\gamma}$ factor can be replaced by a polylogarithmic factor.)

4 A fast algorithm for hdbscan in the plane

Campello et al. [7] introduced hdbscan, a hierarchical version of dbscan∗ (similar to optics [3]). The algorithm described by Campello *et al.*to compute the hdbscan hierarchy runs in quadratic time. We show that in ${\mathbb{R}}^{2}$ and under the Euclidean metric, the hdbscan hierarchy can be computed in $O(n\log n)$ time.

Preliminaries on hdbscan. Recall that dbscan∗ is the version of dbscan in which border points are considered noise. The hdbscan hierarchy is a tree structure encoding the clusterings of dbscan∗ (for a fixed MinPts) that arise as $\varepsilon$ increases from $\varepsilon=0$ to $\varepsilon=\infty$. Initially, when $\varepsilon=0$, all points are noise. As $\varepsilon$ increases, three types of events can happen to the dbscan∗ clustering:

• •

Type (i): the status of a point changes. In this event, a point changes from being noise to being a core point. The value of $\varepsilon$ at which this happens for a point $p$ is called the core distance of $p$; we denote it by $d_{\mathrm{core}}(p)$.

• •

Type (ii): a new cluster starts. This event is triggered by a type (i) event, when a point becoming a core point forms a new (singleton) cluster.

• •

Type (iii): two clusters merge. This event can be triggered by a type (i) event or it can happen when $\varepsilon=|pq|$ for core points $p,q$ from different clusters.

Note that all events happen at values of $\varepsilon$ such that $\varepsilon=|pq|$ for some pair of points $p,q\in D$. Also note that several events may happen simultaneously, not only when two pairwise distances are equal, but also when an event triggers other events. This process can be modeled as a dendrogram: a tree whose leaves correspond to the points in $D$ and whose nodes correspond to clusters arising during the process. This dendrogram, where each node stores the value of $\varepsilon$ at which the corresponding cluster was created, is the hdbscan hierarchy. Notice that we can easily extract the dbscan∗ clustering for any desired value of $\varepsilon$ (with respect to the given MinPts) from the dendrogram in linear time. Campello *et al.*compute the hdbscan hierarchy as follows.

For two points $p,q\in D$, define $d_{\mathrm{mr}}(p,q):=\max\left(d_{\mathrm{core}}(p),d_{\mathrm{core}}(q),|pq|\right)$ to be the mutual reachability distance of $p$ and $q$. The mutual reachability graph $\mathcal{G}_{\mathrm{mr}}$ is defined as the complete graph with node set $D$ in which each edge $(p,q)$ has weight $d_{\mathrm{mr}}(p,q)$. Campello *et al.*observe that hdbscan hierarchy can easily be computed from a minimum spanning tree (mst) on $\mathcal{G}_{\mathrm{mr}}$. (Indeed, the cluster-growing process corresponds to the process of computing an mst on $\mathcal{G}_{\mathrm{mr}}$ using Kruskal’s algorithm [9].) Hence, they compute the hdbscan hierarchy as follows.

Compute the core distances $d_{\mathrm{core}}(p)$ for all points $p\in D$. 2. 2.

Compute an mst $\mathcal{T}$ of the mutual reachability graph $\mathcal{G}_{\mathrm{mr}}$. 3. 3.

Convert $\mathcal{T}$ into a dendrogram where each internal node stores the value of $\varepsilon$ at which the corresponding cluster is formed.

Our planar algorithm. The most time-consuming parts in the algorithm above are Steps 1 and 2; Step 3 takes $O(n)$ time after sorting the edges of $\mathcal{T}$ by weight.

For Step 1 we observe that $d_{\mathrm{core}}(p)$ is the distance of point $p$ to its $\ell$-th nearest neighbor for $\ell=\mbox{{\sc MinPts}}-1$. Hence, to compute all core distances it suffices to compute for each point its $k$ nearest neighbors. This can be done (in any fixed dimension) in $O(n\ell\log n)$ time [21]. Since $\ell=\mbox{{\sc MinPts}}-1=O(1)$ this implies that Step 1 takes $O(n\log n)$ time.

Step 2 is more difficult to do in subquadratic time. The main problem is that we cannot afford to look at all edges of $\mathcal{G}_{\mathrm{mr}}$ when computing $\mathcal{T}$. To overcome this problem we need the following generalization of Delaunay triangulations, introduced by Gudmundsson et al. [13]. Recall that a pair of points $p,q\in D$ forms an edge in the Delaunay triangulation of $D$ if and only if there is a circle with $p$ and $q$ on its boundary and no points from $D$ in its interior [6]. We say that the pair $p,q\in D$ forms a $k$-th order Delaunay edge, or $k$-OD edge for short, if and only if there exists a circle with $p$ and $q$ on its boundary and at most $k$ points from $D$ in its interior [13]. (Thus the 0-OD edges are precisely the edges of the Delaunay triangulation.) The $k$-OD edges are useful for us because of the following lemma.

Lemma 4.1

Let $\overline{\mathcal{G}}_{\mathrm{mr}}$ be the subgraph of $\mathcal{G}_{\mathrm{mr}}$ that contains only the $k$-OD edges, where $k=\max(\mbox{{\sc MinPts}}-3,0)$. Then an mst of $\overline{\mathcal{G}}_{\mathrm{mr}}$ is also an mst of $\mathcal{G}_{\mathrm{mr}}$.

Proof 4.2

Imagine computing an mst $\mathcal{T}$ on $\mathcal{G}_{\mathrm{mr}}$ using Kruskal’s algorithm [9]. This algorithm treats the edges $(p,q)$ of $\mathcal{G}_{\mathrm{mr}}$ in order of increasing weight, that is, increasing values of $d_{\mathrm{mr}}(p,q)$. When it processes $(p,q)$ it checks if $p$ and $q$ are already in the same connected component—in our application this component corresponds to a cluster at the current value of $\varepsilon$—and, if not, merges these components. We will argue that whenever we process an edge $(p,q)$ that is not in $\overline{\mathcal{G}}_{\mathrm{mr}}$, that is, an edge that is not a $k$-OD edge, then $p$ and $q$ are already in the same connected component. Hence, there is no need to process $(p,q)$, which proves that an mst of $\overline{\mathcal{G}}_{\mathrm{mr}}$ is also an mst of $\mathcal{G}_{\mathrm{mr}}$.

*Let $C_{pq}$ be the circle such that $p$ and $q$ form a diametrical pair of $C$, and let $D(C_{pq})\subset D$ be the set of points lying in the interior of $C_{pq}$. If $|D(C_{pq})|\leqslant k$, then $(p,q)$ is a $k$-OD edge, so assume $|D(C_{pq})|\geqslant k+1$. Note that $d_{\mathrm{core}}(r)<|pq|$ for all $r\in D(C_{pq})$. Indeed, since $p,q$ is a diametrical pair of $C_{pq}$, the distance from $r$ to any other point in $C_{pq}$ (including $p$ and $q$) is smaller than $|pq|$. Hence, for $\varepsilon=|pq|$ we have $|N_{\varepsilon}(r,D)|\geqslant|D(C_{pq})|+2=k+3=\mbox{{\sc MinPts}}$. Thus all points $r\in C_{pq}$ are core points when we process $(p,q)$. Moreover, for all edges $(s,t)$ with $s,t\in D(C_{pq})\cup\{p,q\}$ we have $d_{\mathrm{mr}}(s,t)\leqslant|pq|$. Hence, it suffices to prove the following.

Claim: Let $C$ be a circle with two points $p,q$ on its boundary and let $D(C)\subset D$ be the set of points from $D$ in the interior of $C$. Then there is a path from $p$ to $q$ in $\overline{\mathcal{G}}_{\mathrm{mr}}$ that uses only points in $D(C)\cup\{p,q\}$.

We prove this claim by induction on $|D(C)|$. If $|D(C)|\leqslant k$ then $(p,q)$ is a $k$-OD edge itself and we are done. Otherwise, pick any point $r\in D(C)$. Now shrink $C$, while keeping $p$ in its boundary, until we obtain a circle $C_{1}$ that also has $r$ on its boundary, as shown in Figure 3. By induction, there is a path from $p$ to $r$ in $\overline{\mathcal{G}}_{\mathrm{mr}}$ that uses only points in $D(C_{1})\cup\{p,r\}\subset D(C)\cup\{p,q\}$. A similar argument shows that there is a path from $r$ to $q$ that uses only points in $D(C)\cup\{p,q\}$. This proves the claim and, hence, the lemma.*

Gudmundsson *et al.*showed that the number of $k$-OD edges is $O(n(k+1))$ and that the set of all $k$-OD edges can be computed in $O(n(k+1)\log n)$ time with a randomized incremental algorithm. Lemma 4.1 implies that after computing the core distances and the $k$-OD edges in $O(n\log n)$ time—recall that $k=\max(\mbox{{\sc MinPts}}-3,0)=O(1)$—we can compute the mst for $\mathcal{G}_{\mathrm{mr}}$ by considering only $O(n)$ edges. Thus computing the mst can be done in $O(n\log n)$ time [9]. Since the rest of the algorithm takes linear time, we obtain the following theorem.

Theorem 4.3

Let $D$ be a set of $n$ points in ${\mathbb{R}}^{2}$ and MinPts be a given constant. We can compute the hdbscan hierarchy on $D$ for the Euclidean metric with a randomized algorithm in $O(n\log n)$ expected time.

5 Approximate dbscan∗ and hdbscan

The approach from the previous section for computing the hdbscan hierarchy will not give a subquadratic bound in higher dimensions, since the number of Delaunay edges can be quadratic. The running time of the algorithm from Section 3 to compute a single dbscan clustering is subquadratic in any fixed dimension, but the running time quickly deteriorates as the dimension increases. In this section we introduce an approximate version of dbscan∗ and hdbscan, both of which can be computed in linear time. An approximation for dbscan∗ that runs in expected linear time was already provided by Chen et al. [8] Gan and Tao [12], however our approach runs has a slightly better dependency on the approximation factor $\delta$. To the best of our knowledge this is the first linear time approximation algorithm for ${\sc hdbscan}$.

5.1 Approximate dbscan∗

First we define what exactly is an approximate dbscan∗ clustering. Our definition of approximate dbscan∗ is essentially the same as the definitions by Chen et al. [8] and Gan and Tao [12]. The main difference is that we base our definition on dbscan∗ instead of dbscan, which avoids some technical difficulties in the definition.

Let MinPts be a fixed constant. Let $\mathcal{C}_{\varepsilon}(D)$ denote the set of clusters in the dbscan∗ clustering for a given value of $\varepsilon$. We call a clustering $\mathcal{C}_{1}$ a refinement of a clustering $\mathcal{C}_{2}$, denoted by $\mathcal{C}_{1}\prec\mathcal{C}_{2}$, when for every cluster $C_{1}\in\mathcal{C}_{1}$ there is a cluster $C_{2}\in\mathcal{C}_{2}$ with $C_{1}\subseteq C_{2}$. Recall that, as $\varepsilon$ increases, the dbscan∗ clusters merge or expand and new singleton clusters may appear, but clusters do not shrink or disappear. Hence, if $\varepsilon<\varepsilon^{\prime}$ then222Here it is important that we consider dbscan∗ and not dbscan. Indeed, in dbscan border points can “flip” between clusters as $\varepsilon$ increases, and so we do not necessarily have $\mathcal{C}_{\varepsilon}(D)\prec\mathcal{C}_{\varepsilon^{\prime}}(D)$. $\mathcal{C}_{\varepsilon}(D)\prec\mathcal{C}_{\varepsilon^{\prime}}(D)$. An approximate dbscan∗ clustering is now defined as follows.

Definition 5.1

A $\delta$-approximate dbscan∗ clustering of a data set $D$, for given parameters $\varepsilon$ and MinPts, and a given error $\delta>0$, is now defined as a clustering $\mathcal{C}^{*}$ of $D$ into clusters and noise such that $\mathcal{C}_{(1-\delta)\varepsilon}(D)\prec\mathcal{C}^{*}\prec\mathcal{C}_{\varepsilon}(D)$.

Thus if we choose $\delta$ sufficiently small, then a $\delta$-approximate dbscan∗ clustering is very similar to the exact dbscan∗ clustering for the given parameter values.

The algorithm. As mentioned in the introduction, both Chen et al. [8] and Gan and Tao [12] already presented algorithms for this. We obtain a slightly better dependency on $\delta$ than Gan and Tao [12] by plugging in a better algorithm for approximate bichromatic closest pair.

Recall that the bottleneck in computing a dbscan∗ clustering lies in checking, for pairs $b,b^{\prime}$ of neighboring boxes, whether there is a pair $(p,p^{\prime})\in D(b)\times D(b^{\prime})$ with $|pp^{\prime}|\leqslant\varepsilon$. We can perform this check approximately by computing an approximate bichromatic closest pair $(p,p^{\prime})\in D(b)\times D(b^{\prime})$ such that $|pp^{\prime}|\leqslant(1+\alpha)\cdot\mathrm{dist}(D(b),D(b^{\prime}))$ for $\alpha=\delta/(1-\delta)$, where $\mathrm{dist}(D(b),D(b^{\prime}))$ denotes the distance between the points of the true closest pair. This can be done in $O((1/\alpha)^{d/3}(n_{b}+n_{b^{\prime}}))=O((1/\delta)^{d/3}(n_{b}+n_{b^{\prime}}))$ time [4]. We add the edge $(b,b^{\prime})$ to the box graph $\mathcal{G}_{\mathrm{box}}$ when $|pp^{\prime}|\leqslant\varepsilon$. This way we can obtain the following result.

Theorem 5.2

Let $D$ be a set of $n$ points in ${\mathbb{R}}^{d}$, and let $\varepsilon$ and MinPts be given constants. Then, for any given $\delta>0$, we can compute a $\delta$-approximate dbscan∗ clustering on $D$ with respect to $\varepsilon$ and MinPts for the Euclidean metric in $O(n\log n+(1/\delta)^{d/3}n)$ time.

Proof 5.3

Let $\mathcal{C}$ be the clustering computed using the approximate bichromatic closest pairs. When $\mathrm{dist}(D(b),D(b^{\prime})\leqslant\varepsilon$, then the reported approximate bichromatic closest pair $(p,p^{\prime})$ has $|pp^{\prime}|\leqslant\varepsilon$, so the set of edges added to the box graph is a subset of the edge set of the actual box graph at the given value of $\varepsilon$. Hence, $\mathcal{C}\prec\mathcal{C}_{\varepsilon}(D)$. On the other hand, when $\mathrm{dist}(D(b),D(b^{\prime}))\leqslant(1-\delta)\varepsilon$ then the approximate closest pair $(p,p^{\prime})$ has $|pp^{\prime}|\leqslant(1+\alpha)(1-\delta)\varepsilon=\varepsilon$. Hence, we are guaranteed to add all edges of the box graph for $(1-\delta)\varepsilon$ and so $\mathcal{C}_{(1-\delta)\varepsilon}(D)\prec\mathcal{C}$. Thus $\mathcal{C}$ is a $\delta$-approximate dbscan∗ clustering.

Recall that the running time of our algorithm is $O(n\log n)$, plus the time needed to compute the edges of the box graph. In the approximate version the latter step takes time

[TABLE]

Approximate hdbscan. Our definition of an approximate hdbscan hierarchy is based on the definition of $\delta$-approximate dbscan∗ clusterings: we say that a hierarchy is a $\delta$-approximate hdbscan hierarchy if, for any value of $\varepsilon$, the clustering extracted from the hierarchy is a $\delta$-approximate dbscan∗ clustering for that value of $\varepsilon$. Next we show how to compute a $\delta$-approximate hdbscan hierarchy in $O(n\log n)$ time, in any fixed dimension.

As in Section 4 we follow the algorithm by Campello et al. [7], and we speed up Step 2 of the algorithm by computing an mst on a subgraph of the mutual reachability graph $\mathcal{G}_{\mathrm{mr}}$ rather than on the whole graph. (Steps 1 and 3 can still be done in $O(n\log n)$ and $O(n)$ time, respectively.) The difference with the exact algorithm of Section 4 is that we will select the edges of the subgraph in a different manner, using ideas from so-called $\theta$-graphs [18].

Let $p\in D$ be a point. We partition ${\mathbb{R}}^{d}$ into simplicial cones with apex $p$ and whose angular diameter is $\theta$, where $\theta$ will be specified later. (The angular diameter of a cone $c$ with apex $p$ is the maximum angle between any two vectors emanating from $p$ and inside $c$.) Let $\Gamma_{p}$ be the resulting collection of cones and consider a cone $c\in\Gamma_{p}$. Let $D(c)\subseteq D$ denote the set of points inside $c$. (If a point lies on the boundaries of several cones we can assign it to one of these cones arbitrarily.) Pick a half-line $\ell_{c}$ with endpoint $p$ that lies inside $c$. A $\theta$-graph would now be obtained by projecting all points from $D(c)$ orthogonally onto $\ell_{c}$, and adding an edge from $p$ to the point closest to $p$ in this projection, with ties broken arbitrarily. We do the same, except that we add edges to the $k$ closest points for $k:=2\cdot\mbox{{\sc MinPts}}-3$. If $c$ contains fewer than $k$ points, we simply connect $p$ to all points in $D(c)$. Doing this for all the cones $c\in\Gamma_{p}$ gives us a set $E_{p}$ of $O(k/\theta)=O(1/\theta)$ edges for point $p$. Let $E(\theta):=\bigcup_{p\in D}E_{p}$. The set $E(\theta)$ can be computed by making a straightforward adaptation to the algorithm to compute a $\theta$-graph in ${\mathbb{R}}^{d}$ [18, Chapter 5], leading to the following result.

Lemma 5.4

$E(\theta)$* has $O(n/\theta^{d-1})$ edges and can be computed in $O((n/\theta^{d-1})\log^{d-1}n)$ time.*

The set $E(\theta)$, where $\theta$ is chosen such that $\cos\theta\geqslant 1-\delta$, defines the subgraph $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)$ on which we compute the mst in Step 2. Since $\cos\theta>1-\theta^{2}/2$, we have $\cos\theta\geqslant 1-\delta$ when $\theta:=\sqrt{2\delta}$. Next we show that an mst on $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)$ defines a $\delta$-approximate hdbscan clustering.

Lemma 5.5

Let $\mathcal{T}$ be an mst of $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)$ and let $\varepsilon>0$. Let $\mathcal{C}(\mathcal{T},\varepsilon)$ be the clustering induced by $\mathcal{T}$. Then $\mathcal{C}$ is a $\delta$-approximate dbscan∗ clustering for the given $\varepsilon$.

Proof 5.6

For a weighted graph $\mathcal{G}$ and threshold weight $\tau$, let $\mathcal{G}[\tau]$ denote the subgraph obtained by removing all edges of weight greater than $\tau$. In order to show that $\mathcal{C}(\mathcal{T},\varepsilon)\prec\mathcal{C}_{\varepsilon}(D)$ we must show that any connected component of $\mathcal{T}[\varepsilon]$ is contained in a connected component of $\mathcal{G}_{\mathrm{mr}}[\varepsilon]$. Since $\mathcal{T}$ is a subgraph of $\mathcal{G}_{\mathrm{mr}}$ this is obviously the case.

Next we prove that $\mathcal{C}_{(1-\delta)\varepsilon}(D)\prec\mathcal{C}(\mathcal{T},\varepsilon)$. For this we must prove that any connected component of $\mathcal{G}_{\mathrm{mr}}[(1-\delta)\varepsilon]$ is contained in a connected component of $\mathcal{T}[\varepsilon]$. Since $\mathcal{T}$ is an mst of $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)$, the connected components of $\mathcal{T}[\varepsilon]$ are the same as the connected components of $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)[\varepsilon]$. It thus suffices to show the following: for any edge $(p,q)\in\mathcal{G}_{\mathrm{mr}}[(1-\delta)\varepsilon]$, there is a path from $p$ to $q$ in $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)[\varepsilon]$. We show this by induction on $|pq|$, similarly to the way in which it is shown that a $\theta$-graph has a small dilation.

*Let $(p,q)$ be an edge in $\mathcal{G}_{\mathrm{mr}}[(1-\delta)\varepsilon]$. Consider the set $\Gamma_{p}$ of cones with apex $p$ that was used to define the edge set $E_{p}$, and let $c\in\Gamma_{p}$ be the cone containing $q$. Recall that we added an edge from $p$ to the $k$ points in $c$ that are closest to $p$ when projected onto the half-line $\ell_{c}$, where $k:=2\cdot\mbox{{\sc MinPts}}-3$. Hence, when $q$ is one of these $k$ closest points we are done. Otherwise, let $r\in D(c)$ be the $(\mbox{{\sc MinPts}}-1)$-th closest point.

Claim: (i) $d_{\mathrm{core}}(r)\leqslant(1-\delta)\varepsilon$, (ii) $|pr|\leqslant\varepsilon$, and (iii) $|rq|<|pq|$.

Before we prove this claim, we first we argue that the claim allows us to finish our inductive proof. Since $(p,q)$ is an edge in $\mathcal{G}_{\mathrm{mr}}[(1-\delta)\varepsilon]$ we have $d_{\mathrm{mr}}(p,q)\leqslant(1-\delta)\varepsilon$. Thus $|pq|\leqslant(1-\delta)\varepsilon$ and $d_{\mathrm{core}}(q)\leqslant(1-\delta)\varepsilon$. Together with parts (i) and (iii) of the claim this implies that $(r,q)$ is an edge in $\mathcal{G}_{\mathrm{mr}}[(1-\delta)\varepsilon]$ with $|rq|<|pq|$.*

In the base case of our inductive proof, where $(p,q)$ is the shortest edge in $\mathcal{G}_{\mathrm{mr}}[(1-\delta)\varepsilon]$, this cannot occur. Thus $q$ must be one of the $k$ closest points in the cone $c$, and we have an edge between $p$ and $q$ in $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)[\varepsilon]$ by construction.

If we are not in the base case, then we have a path from $r$ to $q$ in $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)[\varepsilon]$ by the induction hypothesis. Moreover, $(p,r)$ is an edge in $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)$ by construction. Since $|pr|\leqslant\varepsilon$ by part (ii) of the claim, we have a path from $p$ to $q$ in $\overline{\mathcal{G}}_{\mathrm{mr}}(\delta)[\varepsilon]$.

*It remains to prove the claim. For this we use the following fact [18, Lemma 4.1.4], which is also used to prove that a $\theta$-graph has small dilation.

Fact: Let $s,t$ be any two points in a cone $c\in\Gamma_{p}$ such that, when projected onto the half-line $\ell_{c}$, the distance from $p$ to $s$ is smaller than the distance from $p$ to $t$. Then $|ps|\leqslant|pt|/\cos\theta$ and $|st|<|pt|$.

Part (iii) of the claim immediately follows from this fact by taking $s:=r$ and $t:=q$. Part (ii) follows again by taking $s:=r$ and $t:=q$, using that $|pq|\leqslant(1-\delta)\varepsilon$ and that we have chosen $\delta$ such that $\cos\theta=1-\delta$. For part (i) we must prove that there are at least $\mbox{{\sc MinPts}}-1$ points within distance $(1-\delta)\varepsilon$ from $r$. Recall that $r$ is the $(\mbox{{\sc MinPts}}-1)$-th closest point to $p$ in the cone $c$, measured in the projection onto the half-line $\ell_{c}$. Let $r_{1},\ldots,r_{k}$ be the $k$ closest points; thus $r=r_{i}$ for $i=\mbox{{\sc MinPts}}-1$. We distinguish two cases: $|pr|\leqslant(1-\delta)\varepsilon$ and $|pr|>(1-\delta)\varepsilon$. See also Fig. 4.*

In the former case we can conclude that $|r_{i}r|\leqslant(1-\delta)\varepsilon$ for all $1\leqslant i\leqslant\mbox{{\sc MinPts}}-2$ by setting $s:=r_{i}$ and $t:=r$ and using $|pr|\leqslant(1-\delta)\varepsilon$. Thus, including the point $p$, we know that $r$ has at least $\mbox{{\sc MinPts}}-1$ points within distance $(1-\delta)\varepsilon$.

In the latter case we will argue that $|r_{i}r|\leqslant(1-\delta)\varepsilon$ for all $\mbox{{\sc MinPts}}\leqslant i\leqslant 2\cdot\mbox{{\sc MinPts}}-3$. Since by part (iii) of the claim we have $|rq|\leqslant(1-\delta)\varepsilon$, we conclude that also in the latter case $r$ has at least $\mbox{{\sc MinPts}}-1$ points within distance $(1-\delta)\varepsilon$. To argue that $|r_{i}r|\leqslant(1-\delta)\varepsilon$ we first note that for any point $s\in c$ we have $|ss^{*}|\leqslant\sin\theta\cdot|ps|$, where $s^{*}$ denotes the orthogonal projection of $s$ onto $\ell_{c}$. Thus

[TABLE]

where the last inequality uses $|pq|\leqslant(1-\delta)\varepsilon$ and that we are now considering the case $|pr|>(1-\delta)\varepsilon$. Since we can assume that $\theta$ is small enough to ensure $2\sin\theta<\cos^{2}\theta$, we conclude that, indeed, $|rr_{i}|\leqslant(1-\delta)\varepsilon$. This finishes the proof for part (i) of the claim and hence, of the lemma.

Combining the previous two lemmas we obtain the following theorem.

Theorem 5.7

Let $D$ be a set of $n$ points in ${\mathbb{R}}^{d}$, and let $\varepsilon$ and MinPts be given constants. Then, for any given $\delta>0$, we can compute a $\delta$-approximate hdbscan clustering on $D$ with respect to $\varepsilon$ and MinPts for the Euclidean metric in $O((n/\delta^{(d-1)/2})\log^{d-1}n)$ time.

6 Experimental evaluation

In this section we experimentally investigate the efficiency of our new algorithm and compare it to the original algorithm. The only goal of these experiments is to serve as a proof of concept to illustrate that indeed for very basic point distributions the original algorithm has a bad running time, whereas the new algorithm performs much better. We first describe some implementation details and we discuss our implementation of the original algorithm. We then describe the data sets and parameters for the tests. Finally we present the results and conclusions.

Implementation details. We implemented two versions of our new algorithm: the strip-based approach for Step 1 and the grid-based approach. In Step 3, where we have to decide for all pairs $b,b^{\prime}$ of neighboring boxes whether there are core points $p\in D(b)$ and $p^{\prime}\in D(b^{\prime})$ with $|pp^{\prime}|\leqslant\varepsilon$, we use the following randomized brute-force approach. (This is instead of the theoretically better Delaunay triangulations or spherical emptiness queries.) Without loss of generality assume $n_{b}\leqslant n_{b^{\prime}}$. For each core point $p\in D(b)$ we first test if $\mathrm{dist}(p,b^{\prime})\leqslant\varepsilon$. If not, no point of $D(b^{\prime})$ can be within distance $\varepsilon$ of $p$. If so, we test each point $p^{\prime}\in D(b^{\prime})$ to see if $|pp^{\prime}|\leqslant\varepsilon$. If during this procedure we find two core points within distance $\varepsilon$ from each other, we can stop. The randomization is obtained by considering the points in each box in random order; it ensures that if there are many pairs within distance $\varepsilon$, we expect to find such a pair early.

The original dbscan algorithm performs a spherical range query for each point $p\in D$ to find all $q\in D$ with $|pq|\leqslant\varepsilon$. To this end an indexing structure such as an R-tree is typically used. In our implementation we use the box-graph to answer the queries. Note that an R-tree also groups the points into boxes at the leaf level; the tree structure is then used find the leaf boxes intersecting the query range, after which the points inside these boxes are tested. The boxes in our box-graph can be seen as being optimized for the radius $\varepsilon$ of the query range, and so the box-graph should be at least as efficient as a general-purpose R-tree.

Experimental set-up. We ran the algorithms on several synthetic data sets in 2D and 4D, each consisting of four clusters. (For other numbers of clusters the results are similar.) The clusters either have a uniform or Gaussian distribution, and their centers are placed roughly 700 units apart in a hypercube with edge length 1,000 units. Uniform clusters are generated within a ball of radius 300 around the cluster center, Gaussian clusters are generated with a standard deviation of 100. For several data sets we added 5% noise to the input, uniformly inside a slightly expanded bounding box around the clusters.

Parameters. We analyse the efficiency of the algorithms with respect to two parameters: the input size and the density within the clusters. As a measure of the density we use $n(r/\varepsilon)^{d}$, where $r$ is the cluster radius; thus, for the uniform distribution, the density represents the expected number of points within distance $\varepsilon$ from a core point, and also the expected number of points in the boxes of the box-graph within the clusters.

Measurements. We compare the algorithms in two different ways: we measure the actual execution times, and the number of pairs of points for which a distance computation is done. The latter is the main operation for finding the clusters in both the new and the original algorithm, so it provides a good implementation-independent measure. For the original algorithm we also count the sum of neighborhood sizes of all points. Following earlier work, we call this the number of seeds. This is a lower bound for the number of operations needed in the original algorithm, independent of the indexing structure used to find the neighborhoods.

Result for fixed input size. In these experiments we fix the input size and run the algorithms with different values of $\varepsilon$. In 2D we ran the algorithms on a data set of each type—uniform or Gaussian and with or without noise—in which the clusters contain 500,000 points each. In 4D we use the same types of data sets, but with 200,000 points per cluster. The results are shown in Fig. 5.

Result for fixed density. In our second set of experiments we use data sets of different sizes, where for each data set we pick $\varepsilon$ such that the density (as defined above) remains constant. For each type of data we generated data sets ranging from 20,000 to 500,000 points per cluster in 2D, and from 20,000 to 300,000 points per cluster for 4D. We ran these experiments for two different values of $\varepsilon$: one that is roughly the smallest value needed to find the four clusters, and one that is roughly the highest value for which the four clusters do not start to merge. The results are given in Fig. 6.

Discussion. The running times in Fig. 5 (top row) show a lot of fluctuation, while the results on the number of distance computations are very stable. We suspect the fluctuation in running time is related to memory issues, or to other processes claiming resources. Hence, most of our conclusions will rely on the number of distance calculations and the sum of neighborhood-sizes.

The results from Fig. 5 clearly show the different dependency on density between the two algorithms. For the original algorithm the number of distance calculations (shown in the second row) grows linearly with the density. The sum of the neighborhood sizes (seeds) is slightly over 50% of the number of distance computations, so the linear dependency is inherent in the original algorithm, which reports the complete neighborhood of all points. The new algorithm, on the other hand, is insensitive to the density. Note that the bottom row in Fig. 5 shows that there is a fairly large range of values of $\varepsilon$ at which the correct clusters are found. (The smallest range occurs for Gaussian clusters with noise, and even there the four clusters are still found for $\varepsilon^{2}$ roughly between 4 and 20.) The original algorithm can compete with the new algorithm only when $\varepsilon$ is chosen very near the lower end of the correct range of $\varepsilon$; for a “safe” value of $\varepsilon$ in the middle of the range, our algorithm is significantly faster. This is important, as determining the right value of $\varepsilon$ is hard and may require running the algorithm several times with different values of $\varepsilon$. Moreover, different clusters may have different densities. Choosing a value for $\varepsilon$ that is near the low end of the range for some clusters may then either fail to find some other clusters or it may lead to an explosion in running time for the other clusters.

Fig. 6 shows that the strip-based approach tends to make fewer comparison than the grid-based approach, especially in the low-density setting shown in the upper row in the figure. (In Fig. 5 the low-density setting is in the far left in each graph, and therefore less visible.) This is caused by two things: the boxes are created in a data-driven way in the strip-based approach, and the boxes are smaller (being bounding boxes of the points in it). The former means that points within distance $\varepsilon$ end up in the same box more often (in which case they are not compared), while the latter means that some pairs of bounding boxes can be at distance more than $\varepsilon$, while the “corresponding” grid cells have distance just below $\varepsilon$. Except for the low-density setting, the actual computation time (see the top row in Fig. 5) of the grid-based approach is better in 2D (even though the number of comparisons is not), due to the smaller constant factors in the approach.

7 Concluding remarks

We presented a new algorithm for dbscan in ${\mathbb{R}}^{2}$, which runs in $O(n\log n)$ time in the worst case. We also presented an $O(n\log n)$ algorithm for hdbscan in ${\mathbb{R}}^{2}$—this is the first subquadratic algorithm for this problem—and we presented near-linear algorithms for approximate hdbscan. It will be interesting to do a more thorough experimental evaluation of our new algorithms. Specifically to compare the approach to other implementations that are current available. It would also be interesting to implement the approximation algorithms and compare them to the exact ones. From the theoretical side, a main open problem is to see if we can compute the exact hdbscan hierarchy in subquadratic time in dimensions $d\geqslant 3$.