# Geometric reconstructions of density based clusterings

A. L. Garcia-Pulido<EMAIL_ADDRESS>Department Of Computer
Science, Ashton Building, University of Liverpool, Liverpool L69 3BX, UK  and
K. P. Samardzhiev<EMAIL_ADDRESS>Geographic Data Science Lab,
University of Liverpool, Roxby Building, Liverpool, L69 7ZT United Kingdom

###### Abstract.

$\operatorname{DBSCAN}^{*}$ and $\operatorname{HDBSCAN}^{*}$ are well
established density based clustering algorithms. However, obtaining the
clusters of very large datasets is infeasible, limiting their use in real
world applications.

By exploiting the geometry of Euclidean space, we prove that it is possible to
systematically construct the $\operatorname{DBSCAN}^{*}$ and
$\operatorname{HDBSCAN}^{*}$ clusters of a finite $X\subset\mathbb{R}^{n}$
from specific subsets of $X$. We are able to control the size of these subsets
and therefore our results make it possible to cluster very large datasets.

To illustrate our theory, we cluster the Microsoft Building Footprint Database
of the US, which is not possible using the standard implementations.

ALGP is a member of the Centre for Topological Data Analysis funded by the
EPSRC grant “New Approaches to Data Science: Application Driven Topological
Data Analysis” EP/R018472/1.

##  1\. Introduction

Let $(X,d)$ be a finite metric space. A _clustering_ of $X$ often refers to
finding disjoint subsets (_clusters_) of $X$ such that points within a subset
are closer to each other than to points in the others. A variant of this is
density based clustering, where in addition, the points in a cluster must
satisfy a given density condition. This clusters $X$ into “high density”
regions. Density based clustering has received interest from many fields
because of its applicability to a vast range of real world problems [1, 16, 5,
18, 2, 15].

A widely used density based clustering algorithm is $\operatorname{DBSCAN}$
[6], we work with $\operatorname{DBSCAN}^{*}$ which is a modification of the
original algorithm. Given parameters $\varepsilon>0$ and $k\in\mathbb{N}$,
$\operatorname{DBSCAN}^{*}$ first considers the set
$\mathcal{C}(\varepsilon)\subset X$ of points $p\in X$ for which the
cardinality of $B_{\varepsilon}(p)$ at least $k+1$. It then clusters
$\mathcal{C}(\varepsilon)$ by adding $p,q\in\mathcal{C}(\varepsilon)$ to the
same cluster if $d(p,q)\leq\varepsilon$ (see Section 2). Formally, we consider
the $\operatorname{DBSCAN}^{*}$ clusters as the connected components of a
graph with vertex set $\mathcal{C}(\varepsilon)$ and $(p,q)$ an edge if
$d(p,q)\leq\varepsilon$ (see 2.3).

Since it is often unclear clear how to choose $\varepsilon$,
$\operatorname{HDBSCAN}^{*}$ evaluates $\operatorname{DBSCAN}^{*}$ clusters
over the full range of $\varepsilon>0$ and, by evaluating how this hierarchy
of clusters evolves, selects the most prevalent clusters of the dataset [3]. A
core part of $\operatorname{HDBSCAN}^{*}$ is the construction of this
hierarchy using a complete weighted graph $(X,\omega)$, where $\omega(p,q)$ is
the maximum of $d(p,q)$ and the least $\varepsilon>0$ such that
$p,q\in\mathcal{C}(\varepsilon)$ (see 2.4). The $\operatorname{DBSCAN}^{*}$
clusters of parameter $\varepsilon$ are then the non-trivial connected
components of a graph $(X,\omega)_{\varepsilon}$ obtained by removing all
edges of weight strictly greater than $\varepsilon$.

An essential component of $\operatorname{DBSCAN}^{*}$ and
$\operatorname{HDBSCAN}^{*}$ are $k$ nearest neighbour searches, but their
complexity grows quadratically and the memory requirements increase with the
size of the dataset. This makes it impractical, or even impossible, to use
these algorithms for cluster analysis of very large datasets. Although there
has been significant effort to improve their performance [7, 14, 8, 10, 9, 17,
13, 4], the tested data and parameters are incomparable to those used in real
world applications. A thorough benchmark [11] shows that, even when running
overnight, it is only feasible to cluster datasets with 5 million points and
with $k\leq 5$. For comparison, the Microsoft Building Footprint Database of
the US [12] has 125 million points and one requires $k\geq 1,000$ to delineate
geographical areas of interest.

In this article we introduce a new framework to obtain the
$\operatorname{DBSCAN}^{*}$ and $\operatorname{HDBSCAN}^{*}$ clusters of very
large datasets. In both cases, we offer an alternative to the explored
approaches of optimising and approximating the steps of the algorithms in
order to do this. Instead, we provide constructions that avoid the lengthy $k$
nearest neighbour searches and prove that they produce the same clusters.

Our starting point is to demonstrate how it is possible to use the natural
partition of Euclidean space into cubes to produce $\operatorname{DBSCAN}^{*}$
clusters subject to minimising the number of distance calculations (see
Section 3). In particular, our construction allows us to skip $k$ nearest
neighbour calculations in highly dense areas of $X$, which are a common
challenge for implementations such as sklearn.DBSCAN. In fact, since our proof
is constructive we use it to derive an algorithm $\operatorname{S-DBSCAN}^{*}$
to construct the $\operatorname{DBSCAN}^{*}$ clusters, which is parallel and
scales to very large datasets.

For $\operatorname{HDBSCAN}^{*}$, we construct a weighted graph
$(X,\omega^{\prime})$ and prove that, for every $\varepsilon>0$, the connected
components of $(X,\omega^{\prime})_{\varepsilon}$ precisely equal those of
$(X,\omega)_{\varepsilon}$, see 5.10. To mitigate the asymptotic complexity of
calculating $k$ nearest neighbours the main idea is to construct this graph
from small independent subsets $S_{i}\subset X$, whilst ensuring that the
connected components are preserved (see Definitions 5.3 and 5.6).

As with $\operatorname{S-DBSCAN}^{*}$, we use our proof to give an algorithm,
$\operatorname{S-HDBSCAN}^{*}$, that constructs the
$\operatorname{HDBSCAN}^{*}$ clusters, see Section 5. We show that it is
possible to iterate through progressively less dense areas of the dataset;
partition the data of each iteration into independent, manageable pieces;
remove parts of the data that are redundant for the processing of subsequent
iterations; and recover the $\operatorname{HDBSCAN}^{*}$ clusters of the full
dataset by processing the small pieces of each iteration.

Reconstructing a global property of a graph relying only on local properties
is a very hard task. In our case, a central factor that makes it possible to
achieve this is, is our careful selection of the $S_{i}$ and the
identification and explicit description of reduced subsets $\partial S_{i}$
that contain the relevant points that interact with other $S_{j}$.

In fact, we develop theory that allows us to study finite subsets of
$\mathbb{R}^{n}$ in a similar way to how one treats open sets in classical
topology, see Section 4. This machinery enables us to incorporate Euclidean
geometry onto finite sets, for example, define notions such as boundary and
neighbourhood of a set; calculate distances between sets and prove that it is
attained near the boundary (see 4.1 and 4.2). This point of view is
instrumental in our construction of $(X,\omega^{\prime})$ and in the proof
that it has the desired property. It proves a great computational asset by
controlling the number of distance calculations when clustering a real world
dataset.

To put our framework into practice, we include examples where we cluster two
real world datasets, using $\operatorname{S-DBSCAN}^{*}$ in Section 6, and
$\operatorname{S-HDBSCAN}^{*}$ in Section 7. Our largest dataset is the entire
Microsoft Building Footprint Database of the United States [12] which we
cluster with $k=1,900$.

##  2\. $\operatorname{DBSCAN}^{*}$ and $\operatorname{HDBSCAN}^{*}$

Given an edge-weighted graph $G$, we write $V(G)$ for the set of vertices of
$G$, $E(G)$ for the set of edges of $G$ and $\omega(G)$ for the weight
function $\omega(G)\colon E(G)\to[0,\infty)$. We will often omit the argument
of $G$ whenever it is clear from the context. We will write $\pi_{0}(G)$ for
the set of connected components of $G$.

Let $G,H$ be two weighted graphs. The _union_ $G\cup H$ is the weighted graph
with vertex set $V(G)\cup V(H)$, edge set $E(G)\cup E(H)$ and weights

$\omega(p,q)=\begin{cases}\min\\{\omega(G)(p,q),\omega(H)(p,q)\\}&(p,q)\in
E(G)\cap E(H)\\\ \omega(G)(p,q)&(p,q)\in E(G)\setminus E(H)\\\
\omega(H)(p,q)&(p,q)\in E(H)\setminus E(G).\end{cases}$

If $\varepsilon\in[0,\infty)$ and $G$ is a weighted graph, the graph
$G_{\varepsilon}$ will denote the subgraph of $G$ with vertex set $V(G)$, set
of edges $E(G_{\varepsilon})=\\{(p,q)\in
E(G):\omega(G)(p,q)\leq\varepsilon\\}$ and with weight
$\omega(G_{\varepsilon}):=\omega(G)|_{E(G_{\varepsilon})}$. For
$k\in\mathbb{N}$, let $G_{\varepsilon,k}$ denote the subgraph of
$G_{\varepsilon}$ induced by vertices of degree at least $k$.

Given a finite metric space $(X,d)$, let $G(X,d)$ be the complete graph on $X$
with weight $\omega(p,q)=d(p,q)$.

######  Remark 2.1.

In the context of topological data analysis, $G(X,d)_{\varepsilon}$ is
sometimes referred to as the $\varepsilon$-neighbourhood graph of $X$ and it
is the $1$-skeleton of the well-known Čech complex of $(X,d)$.

Let $(X,d)$ be a metric space. Given a point $p\in X$ and $\varepsilon\geq 0$,
we denote by $B_{\varepsilon}(p)$ the closed ball of radius $\varepsilon$ with
centre $p$. The cardinality of a set $A$ will be denoted by $|A|$. Given
$A,B\subset X$ and $p\in X$ define the _distance to $A$ from $p$_ as

$d(p,A)=\inf_{a\in A}d(x,a),$

and the _distance between $A$ and $B$_ as

$d(A,B)=\inf_{a\in A}d(a,B).$

For $\varepsilon\geq 0$ define the _$\varepsilon$ -neighbourhood_ of $A\subset
X$ as

$B_{\varepsilon}(A)=\\{p\in X:d(p,A)\leq\varepsilon\\}.$

Throughout this article we work with a fixed $k\in\mathbb{N}$.

In this section we describe two density based clustering algorithms:
$\operatorname{DBSCAN}^{*}$ and $\operatorname{HDBSCAN}^{*}$.

###  2.1. $\operatorname{DBSCAN}^{*}$

In this section, we fix $(X,d)$ a finite metric space and $\varepsilon\geq 0$.
$\operatorname{DBSCAN}^{*}$ is an algorithm that clusters regions of $X$ with
minimum local density $k+1$. The local density of $X$ is measured using the
cardinality of balls with radius $\varepsilon$.

######  Definition 2.2.

A point $p\in X$ is called a _core point_ if $B_{\varepsilon}(p)$ contains at
least $k$ points excluding $p$, that is, if

$|B_{\varepsilon}(p)|>k.$

We denote by $\mathcal{C}(\varepsilon)$ the set of core points and define the
set of _noise points_
$\mathcal{N}(\varepsilon)=X\setminus\mathcal{C}(\varepsilon)$.

Notice that $\mathcal{C}(\varepsilon)=V(G(X,d)_{\varepsilon,k})$.

######  Definition 2.3.

The set of $\operatorname{DBSCAN}^{*}(\varepsilon)$ _clusters_ is

$\mathcal{D}^{*}(\varepsilon)=\\{V(C):C\in\pi_{0}(G_{\varepsilon,k})\\}.$

When $\varepsilon$ is clear from the context, we write $\mathcal{C}$,
$\mathcal{N}$ and $\mathcal{D}^{*}$ instead of $\mathcal{C}(\varepsilon)$,
$\mathcal{N}(\varepsilon)$ and $\mathcal{D}^{*}(\varepsilon)$, respectively.

###  2.2. $\operatorname{HDBSCAN}^{*}$

$\operatorname{HDBSCAN}^{*}$ evaluates
$\operatorname{DBSCAN}^{*}(\varepsilon)$ clusters over the full range
$\varepsilon\in[0,\infty)$. By considering how these clusters evolve and
persist as $\varepsilon$ varies, it is possible to quantify the relative
density of a cluster compared to surrounding regions. The final
$\operatorname{HDBSCAN}^{*}$ clusters are $\operatorname{DBSCAN}^{*}$ clusters
which prevail the most through all scales. These clusters may come from
different $\operatorname{DBSCAN}^{*}$ scales.

We recall the construction of $\operatorname{HDBSCAN}^{*}$ and its persistence
score, described as in [10]. Again, for this section, we fix $(X,d)$ a finite
metric space.

######  Definition 2.4.

A _nearest neighbour_ of $p\in X$ is an element $p_{1}\in X$ such that
$d(p,p_{1})=d(p,X\setminus\\{p\\})$. For $j>1$, a _$j$ -th nearest neighbour_
of $p$ is an element $p_{j}\in X$ which satisfies
$d(p,p_{j})=d(p,X\setminus\\{p,p_{1},\cdots,p_{j-1}\\})$, where $p_{i}$ is an
$i$-th nearest neighbour for $1\leq i\leq j-1$. For $p\in X$, we will write
$\operatorname{core}_{k}(p)$ for the distance from $p$ to a $k$-th nearest
neighbour.

We define the _reachability distance_ $\rho\colon X\times X\to\mathbb{R}$ as

$\rho(p,q)=\operatorname{max}\\{\operatorname{core}_{k}(p),\operatorname{core}_{k}(q),d(p,q)\\}$

if $p\neq q$ and $\rho(p,p)=0$.

To describe $\operatorname{HDBSCAN}^{*}$ we first note that for any finite
metric space $(X,d)$ and $\varepsilon>0$,

$E(G(X,\rho)_{\varepsilon})=E(G(X,d)_{\varepsilon,k}).$

Consequently, for $p\neq q\in X$,

(2.1) $p,q\in\Gamma\text{ for some
}\Gamma\in\pi_{0}(G(X,\rho)_{\varepsilon})\Leftrightarrow
p,q\in\tilde{\Gamma}\text{ for some
}\tilde{\Gamma}\in\pi_{0}(G(X,d)_{\varepsilon,k}).$

Thus, by varying $\varepsilon$ from $0$ to $\operatorname{diam}(X)$,
$G(X,\rho)_{\varepsilon}$ produces a hierarchy of the
$\operatorname{DBSCAN}^{*}$ clusters. $\operatorname{HDBSCAN}^{*}$ first
derives a summary of this hierarchy by defining an equivalence relation on
these clusters.

Let $\varepsilon\geq 0$. For $0\leq\gamma\leq\varepsilon$, the inclusion

$i_{\gamma}\colon G(X,d)_{\gamma,k}\hookrightarrow G(X,d)_{\varepsilon,k},$

gives the induced map

$(i_{\gamma})_{*}\colon\pi_{0}(G(X,d)_{\gamma,k})\hookrightarrow\pi_{0}(G(X,d)_{\varepsilon,k}).$

For $m\in\mathbb{N}$, set

$\mathcal{D}^{*}(\varepsilon,m)=\\{C\in\mathcal{D}^{*}(\varepsilon):|C|\geq
m\\}.$

Define an equivalence relation on

$\bigcup_{\varepsilon\in[0,\infty)}\mathcal{D}^{*}(\varepsilon,m)$

by $C\sim D$, with $C\in\mathcal{D}^{*}(\delta,m)$ and
$D\in\mathcal{D}^{*}(\varepsilon,m)$ and without loss of generality
$\delta\leq\varepsilon$, if and only if

$(i_{\delta})_{*}(C)=D\text{ and }|(i_{\gamma})_{*}^{-1}|=1$

for every $\gamma\in[\delta,\varepsilon]$.

Notice that for an equivalence class $[C]$ and $\gamma\in[0,\infty)$

$|[C]\cap\mathcal{D}^{*}(\gamma,m)|\in\\{0,1\\}.$

Set $C_{\gamma}$ to be the unique element of
$[C]\cap\mathcal{D}^{*}(\gamma,m)$ if the latter is non empty, and $\emptyset$
otherwise. Define the _persistence score_ of an equivalence class $[C]$ by

$\Sigma([C])=\int_{0}^{\infty}\frac{|C_{\gamma}|}{\gamma^{2}}.$

This is a weighted average of the density of the clusters $C_{\gamma}$ within
an equivalence class.

The final step in $\operatorname{HDBSCAN}^{*}$ is to choose final clusters
amongst the equivalence classes. For an equivalence class $[C]$, define

$p([C])=\bigcup_{\varepsilon\geq 0}C_{\varepsilon}.$

Note that since $X$ is a finite metric space there are only a finite number of
equivalence classes, say $\\{[C^{1}],\ldots,[C^{\alpha}]\\}$. We now select
$I\subset\\{1,\ldots,\alpha\\}$ which maximises

$\sum_{i\in J}\Sigma[C^{i}]$

amongst all $J\subset\\{1,\ldots,\alpha\\}$ with

$p([C^{i}])\cap p([C^{j}])=\emptyset$

whenever $i\neq j\in J$. The set of $\operatorname{HDBSCAN}^{*}$ _clusters_ is

$\mathcal{H}^{*}=\\{p([C^{i}]):i\in I\\}.$

##  3\. Reconstructing $\operatorname{DBSCAN}^{*}$ using covering cubes

In this section we describe $\operatorname{S-DBSCAN}^{*}$, our parallel
algorithm that produces $\operatorname{DBSCAN}^{*}$ clusters. This algorithm
relies on the natural partition of Euclidean space into cubes.

From now on, $X$ will denote a finite subset of $\mathbb{R}^{n}$ with $n\geq
2$. We will consider subsets of $\mathbb{R}^{n}$ equipped with either the
$l_{2}$ metric, denoted by $d$, or the reachability distance derived from $d$,
denoted by $\rho$. Whilst our ideas are applicable to any norm on
$\mathbb{R}^{n}$ we restrict to this setting for simplicity.

######  Definition 3.1.

For $\varepsilon>0$, define $\mathcal{Q}(\varepsilon)$ to be the collection of
cubes of the form

$\left\\{(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}:\quad
j_{i}\frac{\varepsilon}{2\sqrt{n}}\leq
x_{i}\leq(j_{i}+1)\frac{\varepsilon}{2\sqrt{n}}\right\\},$

with $j\in\mathbb{Z}^{n}$.

We say that $S,T\in\mathcal{Q}(\varepsilon)$ are _adjacent_ if $S\cap
T\neq\emptyset$. For $S\in\mathcal{Q}(\varepsilon)$ and
$m\in\mathbb{N}\cup\\{0\\}$, define the $m$-extension of $S$ by

$S^{m}=\left\\{x\in\mathbb{R}^{n}:\operatorname{max}{|x_{i}-s_{i}|}\leq
m\frac{\varepsilon}{2\sqrt{n}}\text{ for some }s\in S\right\\}.$

Let $\mathfrak{m}$ be the smallest integer such that $\mathfrak{m}\geq
2\sqrt{n}$. Note that $S^{\mathfrak{m}}$ is the smallest extension of $S$
satisfying $B_{\varepsilon}(S)\subset S^{\mathfrak{m}}$.

For $A\subset\mathbb{R}^{n}$, define

$\mathcal{I}(A)=\\{S\in\mathcal{Q}(\varepsilon):S\cap A\neq\emptyset\\}.$

For $Y\subset\mathbb{R}^{n}$ we write $Y_{A}$ for $Y\cap A$.

For the rest of the section we fix $\varepsilon>0$. Recall that $\mathcal{C}$
are the core points of $\operatorname{DBSCAN}^{*}$ from 2.2.

Note that, because of the choice of side length of cubes
$S\in\mathcal{Q}(\varepsilon)$,

$p\in S,\ q\in S^{1}\Longrightarrow d(p,q)\leq\varepsilon.$

In particular,

(3.1) $p,q\in S_{\mathcal{C}}\Longrightarrow p,q\in C\text{ for some
}C\in\mathcal{D}^{*}.$

$\operatorname{S-DBSCAN}^{*}$ constructs a graph
$\mathcal{G}=(\mathcal{I}(\mathcal{C}),\mathcal{E})$ whose connected
components are in bijective correspondence with the
$\operatorname{DBSCAN}^{*}$ clusters: $S,T\in\mathcal{I}(\mathcal{C})$ are in
the same connected component of $\mathcal{G}$ if and only if all of the
elements of $S_{\mathcal{C}}$ and $T_{\mathcal{C}}$ belong to the same
$\operatorname{DBSCAN}^{*}$ cluster. This relies on the observation in eq.
3.1.

The idea behind the construction of this graph arises from algebraic topology.
In terms of covering spaces, a $\operatorname{DBSCAN}^{*}$ cluster is
equivalent to a connected component of $\cup\mathfrak{B}$, where one regards

$\mathfrak{B}=\\{B_{\varepsilon/2}(x)\\}_{x\in\mathcal{C}}$

as a covering of $\mathcal{C}$ in $\mathbb{R}^{n}$. We consider two other
coverings of $\mathcal{C}$, $\mathfrak{U}=\mathcal{I}(\mathcal{C})$ (a
refinement of $\mathfrak{B}$) and
$\mathfrak{W}=\\{S^{\mathfrak{l}}\\}_{S\in\mathcal{I}(\mathcal{C})}$ (refined
by $\mathfrak{B}$), with $\mathfrak{l}$ the smallest integer
$\mathfrak{l}\geq\sqrt{n}$. The connected components of $\cup\mathfrak{U}$ and
$\cup\mathfrak{W}$ are an under and over approximation of the connected
components of $\cup\mathfrak{B}$. In general, both of $\mathfrak{U}$ and
$\mathfrak{W}$ are significantly smaller than $\mathcal{C}$ and
$\mathfrak{B}$, especially around highly dense areas of $\mathcal{C}$.
Moreover, working with a regular tiling by cubes is more efficient than
working with arbitrary balls.

The first step is to identify the core points $\mathcal{C}$ and the cubes
containing them $\mathcal{I}(\mathcal{C})$. For this we:

  * •

Categorise each cube $S\in\mathcal{Q}(\varepsilon)$ as “dense”, “sparse” or
“locally dense” depending on the cardinality of a neighbourhood of $S$. This
is done so that “dense” cubes only contain core points and “sparse” cubes only
contain noise points.

  * •

For “locally dense” $S\in\mathcal{Q}(\varepsilon)$, identify the core points
$p\in S$, obtain their $\varepsilon$-neighbours $B_{\varepsilon}(p)\cap X$,
and identify the cubes containing any $\varepsilon$-neighbours of points in
$S_{\mathcal{C}}$.

In the second step we build the graph $\mathcal{G}$:

  * •

Begin with $\mathcal{G}_{0}=(\mathcal{I}(\mathcal{C}),\mathcal{E}_{0})$, where
$(S,T)\in\mathcal{E}_{0}$ if $S\cap T\neq\emptyset$; that is, the connected
components of $\mathcal{G}_{0}$ correspond to those of $\cup\mathfrak{U}$. We
use the results from the first step and the cubes from $\mathfrak{W}$ to
remove and add edges to $\mathcal{G}_{0}$ and obtain $\mathcal{G}$.

This approach, and our specific choice of side length of the cubes in
$\mathcal{Q}(\varepsilon)$, allows us to identify large portions of $X$ for
which the $\operatorname{DBSCAN}^{*}$ clusters can be partially constructed
without the need for any pairwise distance calculations.

The first step only requires us to consider a neighbourhood of a cube and can
be computed in parallel. Most of the second step only requires either the
neighbourhood of a cube, or a connected component of a graph, and each can be
processed in parallel.

Note for computational reasons it is more convenient to use the covers
$\\{S^{1}\\}_{S\in\mathcal{I}(\mathcal{C})}$ and
$\\{S^{\mathfrak{\mathfrak{m}}}\\}_{S\in\mathcal{I}(\mathcal{C})}$, instead of
$\mathfrak{U}$ and $\mathfrak{W}$, which we do below.

###  3.1. Categorising cubes

######  Observation 3.2.

Let $S\in\mathcal{Q}(\varepsilon)$. If $|S_{X}^{1}|>k$, then
$S_{X}\subset\mathcal{C}$. If $|S^{\mathfrak{m}}_{X}|\leq k$ or
$S_{X}=\emptyset$, then $S_{X}\subset\mathcal{N}$.

###### Proof.

The first assertion follows from the fact that if $S$ and $T$ are adjacent
cubes in $\mathcal{Q}(\varepsilon)$, then $d(p,q)\leq\varepsilon$ for any
$p\in S$ and $q\in T$. The second is satisfied since, for any $p\in S$,
$B_{\varepsilon}(p)\subset S^{\mathfrak{m}}$. ∎

This observation motivates the following definition:

######  Definition 3.3.

We say that $S\in\mathcal{Q}(\varepsilon)$ is

  1. (1)

_dense_ if $S_{X}\neq\emptyset$ and $|S^{1}_{X}|>k$,

  2. (2)

_sparse_ if $|S^{\mathfrak{m}}_{X}|\leq k$ or $S_{X}=\emptyset$,

  3. (3)

_locally dense_ if $S_{X}\neq\emptyset$, $|S^{1}_{X}|<k$ and
$|S^{\mathfrak{m}}_{X}|>k$.

By 3.2, we only have to determine the core points of locally dense cubes,
which significantly reduces the number of pairwise distance calculations.

###  3.2. Identifying core points in locally dense cubes and their
$\varepsilon$-neighbours

We use the following observation to identify the core points in a locally
dense cube and their $\varepsilon$-neighbours.

######  Observation 3.4.

Let $S\in\mathcal{Q}(\varepsilon)$ be locally dense. A $p\in S$ is a core
point if and only if

$|B_{\varepsilon}(p)\cap(S_{X}^{\mathfrak{m}}\setminus S^{1})|>k-|S^{1}_{X}|.$

This implies that for every locally dense cube $S\in\mathcal{Q}(\varepsilon)$
and $p\in S$ we only need to calculate $d(p,q)$ for $q\in
S^{\mathfrak{m}}\setminus S^{1}$. If we deduce that $p\in S_{\mathcal{C}}$, we
store these $\varepsilon$-neighbours of $p$ in a list $\mathcal{B}_{p}$ and
add the elements of $S^{1}_{X}\setminus\\{p\\}$. We also store a list
$\mathcal{B}_{S}$ of the cubes $T$ such that $T_{X}\cap
B_{\varepsilon}(S_{\mathcal{C}})\neq\emptyset$.

###  3.3. Building the graph

Let us first consider the graph
$\mathcal{G}_{1}=(\mathcal{I}(\mathcal{C}),\mathcal{E}_{1})$ with
$(S,T)\in\mathcal{E}_{1}$ if

  1. (1)

$S,T$ are adjacent or,

  2. (2)

$S,T$ are locally dense, $T\in\mathcal{B}_{S}$ and $S\in\mathcal{B}_{T}$ or

  3. (3)

$S$ is locally dense, $T$ is dense and $T\in\mathcal{B}_{S}$.

This graph is easy to construct using the previous results of this section.
However, $\pi_{0}(\mathcal{G}_{1})$ only approximates the
$\operatorname{DBSCAN}^{*}$ clusters.

The first problem is if $S,T\in\mathcal{I}(\mathcal{C})$ are as in Item 2, we
cannot guarantee the existence of $p\in S_{\mathcal{C}}$ and $q\in
T_{\mathcal{C}}$ with $d(p,q)\leq\varepsilon$. Compare this to Items 1 and 3,
where finding such $p,q$ is guaranteed. We prune $\mathcal{G}_{1}$ to correct
this.

Precisely, we calculate $\pi_{0}(\mathcal{G}_{1})$ and for each
$\Gamma\in\pi_{0}(\mathcal{G}_{1})$ we remove edges $(S,T)$ with $S,T$ non
adjacent and such that for every $p\in S$,

$B_{\varepsilon}(p)\cap T_{\mathcal{C}}=\emptyset.$

We do this in parallel across the connected components of $\mathcal{G}_{1}$ to
construct $\mathcal{G}_{2}=(\mathcal{I}(\mathcal{C}),\mathcal{E}_{2})$. By
construction, if $(S,T)\in\mathcal{E}_{2}$, then there are $p\in
S_{\mathcal{C}}$ and $q\in T_{\mathcal{C}}$ with $d(p,q)\leq\varepsilon$. In
particular, $p,q\in C$ for some $C\in\mathcal{D}^{*}$. By eq. 3.1,
$S_{\mathcal{C}}\cup T_{\mathcal{C}}\subset C$.

Conversely, suppose that $p\in S_{C}$ and $q\in T_{\mathcal{C}}$ with
$d(p,q)\leq\varepsilon$. If $S,T$ are adjacent or one of them is locally
dense, then by definition $(S,T)\in\mathcal{E}_{1}$ and this edge is not
pruned when constructing $\mathcal{E}_{2}$.

The remaining problem is if $S,T\in\mathcal{I}(\mathcal{C})$ are both dense
and belong to different connected components of $\mathcal{G}_{2}$. For every
pair of such cubes with $T\subset S^{\mathfrak{m}}\setminus S^{1}$, we
calculate pairwise distances to determine if there are $p\in
S_{X}\subset\mathcal{C}$ and $q\in T_{X}\subset\mathcal{C}$ with
$d(p,q)\leq\varepsilon$. If so, we add $(S,T)$ to form a final graph
$\mathcal{G}$ and stop processing the connected components containing $S$ and
$T$.

We handle this case last to avoid as many pairwise distance calculations
involving dense cubes as possible; this is crucial as their cardinality can be
very large. By doing this step last we increase the likeliness that two dense
cubes will be connected by a path constructed earlier in the process.

##  4\. Geometric constructions relative to covering cubes

In this section, we fix $\varepsilon>0$. Recall from 3.1 the set of cubes
$\mathcal{Q}=\mathcal{Q}(\varepsilon)$ and its properties. We also fix
$A\subset X$ and recall that $\mathcal{I}(A)=\\{S\text{
cube}:S_{A}\neq\emptyset\\}$. We use the elements of $\mathcal{Q}$ to define
standard notions from topology (interior, boundary and closure) of $A$
relative to $\mathcal{Q}$. These constructions are a fundamental component of
our alternate construction of $\operatorname{HDBSCAN}^{*}$ clusters.

At the end of the section we demonstrate how these constructions allow us to
easily extend $\operatorname{S-DBSCAN}^{*}$ to $\operatorname{S-DBSCAN}$, a
parallel algorithm that obtains the $\operatorname{DBSCAN}$ clusters, whilst
minimising pairwise distance calculations.

######  Definition 4.1.

A cube $S\in\mathcal{I}(A)$ is an _interior cube_ of $A$ if $S^{1}_{X}\subset
A$ and $T\in\mathcal{I}(A)$ for every cube $T\subset S^{1}$. The _interior of
$A$_ is the union of interior cubes,

$\operatorname{int}(A)=\bigcup_{S\text{ interior}}S.$

A cube $S\in\mathcal{I}(A)$ is a _boundary cube_ of $A$ if $S$ is not an
interior cube. The _boundary of $A$_ is defined as

$\partial A=\bigcup\limits_{S\text{ boundary}}S,$

and the _closure of $A$_ is

$\bar{A}=\operatorname{int}(A)\cup\partial
A=\bigcup\limits_{S\in\mathcal{I}(A)}S.$

Let $Y\subset\mathbb{R}^{n}$ and $N\geq 0$. We define the _$N$ -extension_ of
$Y$ to be

$Y^{N}=\bigcup\limits_{S\in\mathcal{I}(Y)}S^{N}$

From now on we will write $\mathfrak{n}$ for the smallest integer satisfying
$\mathfrak{n}\geq\sqrt{n}-1$. Notice that $\mathfrak{n}\geq 1$.

We now show that the distance to a point outside $A$ is attained near the
boundary of $A$.

######  Proposition 4.2.

Let $x\in X\setminus A$ and $r\in\mathbb{N}\cup\\{0\\}$. For any $p\in
A\setminus(\partial A)^{\mathfrak{n}+r}_{A}$, there exists $a\in(\partial
A)_{A}$ such that

$d(p,x)>d(a,x)+r\frac{\varepsilon}{2\sqrt{n}}.$

In particular, for any $x\in X\setminus A$ and $p\in A$, there exists
$a\in(\partial A)^{\mathfrak{n}}_{A}$ such that

$d(a,x)\leq d(p,x).$

###### Proof.

Let $x,p$ be as in the hypothesis and $S\in\mathcal{Q}$ with $p\in S$.

First note that

(4.1) $S^{\mathfrak{n}+r}\subset\operatorname{int}(A).$

Indeed $S\subset\operatorname{int}(A)$ because $S\in\mathcal{I}(A)$ and
$p\not\in\partial A$. If $i\geq 0$ is maximal such that
$S^{i}\subset\operatorname{int}(A)$ and $i<\mathfrak{n}+r$, there exists
$T\in\mathcal{Q}$ with $T\subset S^{i+1}\cap\partial A$ and this would imply

$S\subset T^{i+1}\subset(\partial A)^{i+1}\subset(\partial
A)^{\mathfrak{n}+r},$

contradicting the hypothesis. Note that eq. 4.1 implies

(4.2) $T\in\mathcal{Q},\ T\subset S^{\mathfrak{n}+r+1}\Longrightarrow
T\in\mathcal{I}(A)\text{ and }T\not\in\mathcal{I}(X\setminus A).$

Let $[p,x]$ denote the line segment from $p$ to $x$ and let
$0=t_{0}<t_{1}<t_{2}<\ldots<t_{j}=1$ and $p_{i}=(1-t_{i})p+t_{i}x$ be such
that

$[p_{i-1},p_{i}]\in\\{S\cap[p,x]:S\in\mathcal{I}([p,x])\\}.$

We write $Q_{i}$ for the cube satisfying $[p_{i-1},p_{i}]=Q_{i}\cap[p,x]$.

Let $i$ be maximal such that $Q_{i}\in\mathcal{I}(A)$ and
$Q_{i}\not\in\mathcal{I}(X\setminus A)$. In particular, $Q_{i}\subset\partial
A$ and by eq. 4.2 $Q_{i+1}\not\subset S^{\mathfrak{n}+r+1}$. Since
$p\not\in(\partial A)^{\mathfrak{n}+r}\supset(Q_{i})^{\mathfrak{n}+r}$ we have
the strict inequality

(4.3) $d(p,p_{i})>(\mathfrak{n}+r+1)\frac{\varepsilon}{2\sqrt{n}}.$

Let $a\in(\partial A)_{A}\cap Q_{i}$, so that

(4.4) $d(p_{i},a)\leq\frac{\varepsilon}{2}.$

Since $p,p_{i},x$ all lie in the same line segment, we have

$d(p,x)=d(p,p_{i})+d(p_{i},x).$

Combining this and eqs. 4.4 and 4.3 and using the reverse triangle inequality,

$\displaystyle d(p,x)$
$\displaystyle>(\mathfrak{n}+r+1)\frac{\varepsilon}{2\sqrt{n}}+d(x,a)-d(a,p_{i})$
$\displaystyle\geq(\mathfrak{n}+r+1)\frac{\varepsilon}{2\sqrt{n}}+d(x,a)-\frac{\varepsilon}{2}$
$\displaystyle\geq r\frac{\varepsilon}{2\sqrt{n}}+d(x,a),$

using the fact that $\mathfrak{n}+1\geq\sqrt{n}$. ∎

We next show that in order to find the closure of the extension of a set, we
only need to extend its boundary.

######  Lemma 4.3.

If $r\in\mathbb{N}\cup\\{0\\}$, then
$(\bar{A})^{r}=\operatorname{int}(A)\cup(\partial A)^{r}.$

###### Proof.

Certainly $(\bar{A})^{r}\supset\operatorname{int}(A)\cup(\partial A)^{r}.$ To
prove the other containment, let $S\in\mathcal{I}(A)$. If $S\subset\partial
A$, then $S^{r}\subset(\partial A)^{r}$. Otherwise, suppose
$S\subset\operatorname{int}(A)$ and let $T\in\mathcal{Q}$ be such that
$T\subset S^{r}$. Notice that if $T\in\mathcal{I}(A)$, then
$T\subset\bar{A}\subset\operatorname{int}(A)\cup(\partial A)^{r}$.

Now suppose $T\not\in\mathcal{I}(A)$. Then there exists $t>0$ such that
$T^{t-1}_{A}=\emptyset$ and $T^{t}_{A}\neq\emptyset$. Since $S\subset T^{r}$
and $S\in\mathcal{I}(A)$, $t<r$. Let $Q\subset T^{t}$ be such that
$Q\in\mathcal{I}(A)$. Then $Q\subset\partial A$ and so $T\subset(\partial
A)^{t}\subset(\partial A)^{r}$. ∎

Finally, we note that to find the $\varepsilon$-neighbours of $A$ it is
sufficient to extend the boundary $\mathfrak{m}$ times.

######  Corollary 4.4.

$B_{\varepsilon}(A)\subset\bar{A}^{\mathfrak{m}}=\operatorname{int}(A)\cup(\partial
A)^{\mathfrak{m}}.$

###### Proof.

For any $S\in\mathcal{I}(A)$, $B_{\varepsilon}(S)\subset S^{\mathfrak{m}}$ and
so $S^{\mathfrak{m}}\subset\bar{A}^{\mathfrak{m}}$. In consequence,

$B_{\varepsilon}(A)\subset\bigcup\limits_{S\in\mathcal{I}(A)}S^{\mathfrak{m}}\subset\bar{A}^{\mathfrak{m}}.$

∎

###  4.1. $\operatorname{S-DBSCAN}$

The original density-based clustering algorithm $\operatorname{DBSCAN}$ was
first introduced in [6]. In contrast to $\operatorname{DBSCAN}^{*}$,
$\operatorname{DBSCAN}$ makes a distinction between non-core points which are
$\varepsilon$-neighbours of core points and those that are not.

######  Definition 4.5.

An $x\in\mathcal{N}(\varepsilon)$ is a _border point_ if $x\in
B_{\varepsilon}(C)$ for some $C\in\mathcal{D}^{*}(\varepsilon)$.
$\operatorname{DBSCAN}$ assigns any border point $x$ to the first cluster
$C\in\mathcal{D}^{*}(\varepsilon)$ that the algorithm finds such that $x\in
B_{\varepsilon}(C)$. We denote the set of $\operatorname{DBSCAN}(\varepsilon)$
clusters by $\mathcal{D}(\varepsilon)$.

Using $\operatorname{S-DBSCAN}^{*}$ and the results from Section 4, it is
straightforward to obtain $\operatorname{DBSCAN}$ clusters.

######  Remark 4.6.

Fix $C\in\mathcal{D}^{*}$. By 4.4, the border points of $C$ are contained in
$(\partial C)^{\mathfrak{m}}_{\mathcal{N}}$. By 4.2, $x\in(\partial
C)^{\mathfrak{m}}_{\mathcal{N}}$ is a border point of $C$ if there exists
$p\in(\partial C)_{C}^{\mathfrak{n}}$ such that $d(p,x)\leq\varepsilon$. For
locally dense $S\in\mathcal{Q}$, we previously calculated the
$\varepsilon$-neighbours of $p\in S$ and so we can immediately add those that
belong to $\mathcal{N}$. For dense cubes $S\in\mathcal{Q}$, we can immediately
add $S^{1}_{\mathcal{N}}$ to $C$. Thus we only need to restrict to pairwise
distance calculations between core points in a dense cube $S\subset(\partial
C)^{\mathfrak{n}}$ and noise points in $T\subset S^{\mathfrak{m}}\setminus
S^{1}$.

##  5\. Reconstructing $\operatorname{HDBSCAN}^{*}$ using cluster boundaries

In this section we use the theory from Section 4 to give an alternate
construction to obtain the $\operatorname{HDBSCAN}^{*}$ clusters of $X$. In
fact, we construct a weighted graph $\mathcal{F}$ satisfying

(5.1)
$\\{V(C):C\in\pi_{0}(\mathcal{F}_{\alpha})\\}=\\{V(C):C\in\pi_{0}(G(X,\rho)_{\alpha})\\},$

for every $\alpha\geq 0$. The definition of the persistence score and the
process used to choose the final clusters from Section 2.2 can be applied to
$\mathcal{F}$ to produce the $\operatorname{HDBSCAN}^{*}$ clusters. Our proofs
are constructive and we use them in Section 7 to describe our algorithm
$\operatorname{S-HDBSCAN}^{*}$.

We construct $\mathcal{F}$ by constructing smaller weighted graphs restricted
to subsets $C_{i}\subset X$ and combine them to form $\mathcal{F}$. Since our
central requirement is that $\mathcal{F}$ satisfies eq. 5.1 and that its
simple construction is computationally feasible for very large $X$, we must
ensure that

  1. (1)

we can prescribe a set $N_{k}(C_{i})$ (the search space for $k$ nearest
neighbours) that is simple, comparable in size to $C_{i}$ and contains the
$k$-nearest neighbours of elements in $C_{i}$;

  2. (2)

we can independently construct graphs for each $C_{i}$ and have minimal
interaction between different $C_{i}$;

  3. (3)

we can prescribe reduced subsets of $C_{i}$ that contain the relevant points
that interact with other $C_{j}$.

In order to meet the above constraints we

  * •

select sets $C_{i}$ with a uniform bound
$\operatorname{core}_{k}\leq\varepsilon$ so that taking $N_{k}(C_{i})$ to be
an extension of $C_{i}$ (as in Section 4); that is,
$C_{i}\subset\mathcal{C}(\varepsilon)$ item 1;

  * •

iteratively construct our graph using an increasing sequence of $\varepsilon$
and choose the $C_{i}$ at each iteration so that the $N_{k}(C_{i})$ do not
intersect other $C_{j}$ item 2.

  * •

ensure $C_{i}$ is a union of $\mathcal{D}^{*}(\varepsilon)$ clusters.
Otherwise, since all the elements in a $\mathcal{D}^{*}(\varepsilon)$ cluster
are in the same connected component at scale $\varepsilon$, we would need to
search other $C_{j}$ for edges of weight in $[0,\varepsilon]$ item 2.

  * •

ensure $C_{i}$ _is_ a $\mathcal{D}^{*}(\varepsilon)$ cluster. Otherwise,
$C_{i}$ could be dispersed throughout $X$ and its extension could be the whole
of $X$ item 2.

  * •

prove that an extension of $\partial C_{i}$, for each $C_{i}$ in an iteration,
contains all the points that interact with previous iterations. By
construction, the $C_{i}$ of one iteration do not interact with each other
item 3.

Having established that $\mathcal{F}$ must be constructed iteratively using
$\operatorname{DBSCAN}^{*}$ clusters to partition $X$, we begin by considering
an initial clustering of $X$, $\mathcal{D}^{*}(\varepsilon_{1})$, and obtain
the graph $G(C,\rho)$ for every cluster
$C\in\mathcal{D}^{*}(\varepsilon_{1})$. In practice, one picks
$\varepsilon_{1}$ so that calculating $\rho$ in each cluster is feasible.

We then wish to repeat this for a choice of $\varepsilon_{2}>\varepsilon_{1}$.
However, increasing the scale can result in clusters that are too large for
local calculations to be efficient or feasible. This would also not take into
account the graphs obtained in the initial step. Instead, we consider a subset
of $X$ consisting of $\mathcal{N}(\varepsilon_{1})$ and an extension of the
boundary of each cluster in $\mathcal{D}^{*}(\varepsilon_{1})$.

The key result is to show that the core points that play a role in forming
$\operatorname{DBSCAN}^{*}$ clusters at any subsequent scale
$\alpha\geq\varepsilon_{1}$ (i.e. points that merge two clusters or that
connect a cluster to a noise point) lie near the boundary of clusters in
$\mathcal{D}^{*}(\varepsilon_{1})$. Indeed, in 5.2 we determine an explicit
size of extension of the boundary (depending only upon $n$) which always
contains these points.

We define the subset $X_{2}$ containing this extension of each cluster in
$\mathcal{D}^{*}(\varepsilon_{1})$ and $\mathcal{N}(\varepsilon_{1})$ of $X$,
see 5.3. We make sure that the reachability distance _with respect to_ $X_{2}$
coincides with that of $X$ in all points that play a role in forming future
$\operatorname{DBSCAN}^{*}$ clusters, see 5.4. We then cluster $X_{2}$ using
$\varepsilon_{2}$ and obtain a graph $G(C,\rho)$ for each cluster
$C\in\mathcal{D}^{*}(\varepsilon_{2})$. These graphs are then combined with
the graphs of the initial step.

This process can be repeated for any increasing sequence of $\varepsilon_{i}$,
see 5.6. The main inductive argument is 5.7. To complete the hierarchical
clustering, we define $\mathcal{F}$ by combining the constructed graph with
the subgraph of $G(X_{i+1},\rho)$ with edge set

$\\{(p,q)\in E(G(X_{i+1},\rho)):\omega(p,q)\geq\varepsilon_{i}\\},$

see 5.9. In 5.10, we prove that eq. 5.1 is satisfied.

######  Remark 5.1.

One could also construct approximations of the $\operatorname{HDBSCAN}^{*}$
clusters by choosing to terminate the construction of $\mathcal{F}$ at
iteration $i$ based on a given condition and discard the remaining noise
points. For example, we could terminate the core points exceed a given
cardinality threshold, or simply skip the last step for a large
$\varepsilon_{i}$. The discarded points will have little effect on the overall
clustering because they are the least dense and hence contribute least to the
score.

Since we work with subsets of $X$ clustered using different values of
$\varepsilon$, we first introduce some notation to accommodate this.

For $Y\subset\mathbb{R}^{n}$ and $\varepsilon\geq 0$, we will write
$\mathcal{D}^{*}(Y,\varepsilon)$ for the set of $\operatorname{DBSCAN}^{*}$
clusters of the metric space $(Y,d)$ and write $\mathcal{C}(Y,\varepsilon)$
and $\mathcal{N}(Y,\varepsilon)$ for the corresponding sets of core and noise
points, respectively. For $p,q\in Y$, we write
$\operatorname{core}^{Y}_{k}(p)$ and $\rho^{Y}(p,q)$ for
$\operatorname{core}_{k}(p)$ and $\rho(p,q)$ relative to the metric space
$(Y,d)$. Notice
$\operatorname{core}^{X}_{k}(p)\leq\operatorname{core}^{Y}_{k}(p)$ and
$\rho^{X}(p,q)\leq\rho^{Y}(p,q)$.

From now on, let $\mathfrak{N}=\mathfrak{n}+\mathfrak{m}$. Recall that
$\mathfrak{n}$ and $\mathfrak{m}$ are the least integers such that
$\mathfrak{n}\geq\sqrt{n}-1$ and $\mathfrak{m}\geq 2\sqrt{n}$ as defined in
Section 4 and 3.1, respectively.

The following is an immediate corollary of 4.2.

######  Corollary 5.2.

Let $\varepsilon>0$, $Y\subset X$ and $C\in\mathcal{D}^{*}(Y,\varepsilon)$.
For any $p\in C$ and $x\in Y\setminus C$, there exists $a\in(\partial
C)^{\mathfrak{n}}_{C}$ such that

$d(a,x)\leq d(p,x).$

Next we identify a set for which the relative reachability distance coincides
with the global reachability distance for all points that are found using the
previous corollary and all noise points.

######  Definition 5.3.

For $\varepsilon\geq 0$ and $Y\subset X$, define

$J(Y,\varepsilon)=\bigcup\limits_{C\in\mathcal{D}^{*}(Y,\varepsilon)}(\partial
C)_{C}^{\mathfrak{n}}\bigcup\mathcal{N}(Y,\varepsilon)$

and

$F(Y,\varepsilon)=\bigcup\limits_{C\in\mathcal{D}^{*}(Y,\varepsilon)}(\partial
C)_{C}^{\mathfrak{N}}\bigcup\mathcal{N}(Y,\varepsilon).$

######  Lemma 5.4.

Let $Y\subset X$ and $\varepsilon>0$. If $p,q\in J(Y,\varepsilon)$, then

$\rho^{F(Y,\varepsilon)}(p,q)=\rho^{Y}(p,q).$

###### Proof.

We show that for every $p\in J(Y,\varepsilon)$ and $1\leq i\leq k$, any $i$-th
nearest neighbour of $p$ in $Y$ is contained in $F(Y,\varepsilon)$, from which
the result follows.

First suppose that $p\in(\partial C)^{\mathfrak{n}}_{C}$ for some
$C\in\mathcal{D}^{*}(Y,\varepsilon)$. We will show that

(5.2) $B_{\varepsilon}(p)_{Y}\subset F(Y,\varepsilon).$

Indeed let $y\in B_{\varepsilon}(p)_{Y}$, then

$\operatorname{max}_{1\leq j\leq n}{|y_{j}-p_{j}|}\leq
d(y,p)\leq\varepsilon\leq\frac{\mathfrak{m}}{2\sqrt{n}}\varepsilon.$

Since $p\in(\partial C)^{\mathfrak{n}}$, there exists $x\in\partial C$ with

$\operatorname{max}_{1\leq j\leq
n}{|x_{j}-p_{j}|}\leq\frac{\mathfrak{n}}{2\sqrt{n}}\varepsilon.$

By the triangle inequality,

$\operatorname{max}_{1\leq j\leq
n}{|y_{j}-x_{j}|}\leq\frac{\mathfrak{N}}{2\sqrt{n}}\varepsilon,$

and therefore $y\in(\partial C)^{\mathfrak{N}}\subset F(Y,\varepsilon)$,
proving eq. 5.2. Since $p\in\mathcal{C}(Y,\varepsilon)$, any $i$-th nearest
neighbour of $p$ lies in $B_{\varepsilon}(p)$.

Now suppose that $p\in\mathcal{N}(Y,\varepsilon)$. Let $y\in Y\setminus
F(Y,\varepsilon)$, say $y\in C\setminus(\partial C)^{\mathfrak{N}}_{C}$ for
some $C\in\mathcal{D}^{*}(Y,\varepsilon)$. Let $a\in(\partial C)_{C}$ be given
by 4.2 with

$r=\mathfrak{N}-\mathfrak{n}=\mathfrak{m}\geq 2\sqrt{n}$

such that

$d(p,y)>d(p,a)+\mathfrak{m}\frac{\varepsilon}{2\sqrt{n}}\geq
d(p,a)+\varepsilon.$

If $a_{i}$ is an $i$-th nearest neighbour of $a$, for some $1\leq i\leq k$,
then

$d(p,y)>d(p,a)+d(a,a_{i})\geq d(p,a_{i}).$

By eq. 5.2 applied to $a$, we have $a_{i}\in F(Y,\varepsilon)$. That is, for
any $y\in Y\setminus F(Y,\varepsilon)$, we have found $k$ closer points to $p$
that lie in $F(Y,\varepsilon)$, as required. ∎

We next reformulate 5.2 in terms of the reachability distance.

######  Corollary 5.5.

Let $Y\subset X$ and $\varepsilon>0$. Suppose $p\in C$ for some
$C\in\mathcal{D}^{*}(Y,\varepsilon)$ and $x\in Y\setminus C$. There exists
$a\in(\partial C)_{C}^{\mathfrak{n}}\subset J(Y,\varepsilon)$ such that

$\rho^{Y}(a,x)\leq\rho^{Y}(p,x).$

###### Proof.

Let $a\in(\partial C)^{\mathfrak{n}}_{C}$ be given by 5.2 so that $d(a,x)\leq
d(p,x)$. Since $a\in C$ and $x\not\in C$,

$\operatorname{core}^{Y}_{k}(a)\leq\varepsilon\leq d(a,x).$

Therefore,

$\rho^{Y}(a,x)=\operatorname{max}\\{\operatorname{core}^{Y}_{k}(x),d(a,x)\\}\leq\operatorname{max}\\{\operatorname{core}^{Y}_{k}(x),d(p,x)\\}\leq\rho^{Y}(p,x).$

∎

We now consider a graph constructed using $F$. It depends on an increasing
sequence $0=\varepsilon_{0}<\varepsilon_{1}<\ldots$, which we now fix.

######  Definition 5.6.

Set $X_{1}=X$ and $X_{i+1}=F(X_{i},\varepsilon_{i})$ for $i\geq 1$. Let $H(0)$
be the weighted graph with vertex set $X_{1}$ and no edges. For $i\geq 0$
define

$H(i+1)=\bigcup\limits_{C\in\mathcal{D}^{*}(X_{i+1},\varepsilon_{i+1})}G(B_{\varepsilon_{i+1}}(C),\rho^{B_{\varepsilon_{i+1}}(C)})_{\varepsilon_{i+1}}\bigcup
H(i).$

######  Proposition 5.7.

Let $p_{0},q_{0}\in X$. For every $i\in\mathbb{N}$, there exist
$p_{i},q_{i}\in J(X_{i},\varepsilon_{i})$ such that

  * •

$p_{i-1}=p_{i}$ or $p_{i-1},p_{i}\in A_{i}$ for
$A_{i}\in\mathcal{D}^{*}(X_{i},\varepsilon_{i})$,

  * •

$q_{i-1}=q_{i}$ or $q_{i-1},q_{i}\in B_{i}$ for
$B_{i}\in\mathcal{D}^{*}(X_{i},\varepsilon_{i})$,

and

(5.3) $\rho^{X_{i}}(p_{i},q_{i})\leq\rho^{X_{i-1}}(p_{i-1},q_{i-1}).$

In particular, there exist $p_{i},q_{i}\in J(X_{i},\varepsilon_{i})$ with

$\rho^{X_{i+1}}(p_{i},q_{i})\leq\rho^{X}(p_{0},q_{0})$

and paths in $H(i)$ that join $p_{0}$ to $p_{i}$ and $q_{0}$ to $q_{i}$.

###### Proof.

Set $p_{-1}=p_{0}$, $q_{-1}=q_{0}$ and $X_{-1}=X_{0}=X$. Note that with such
definitions eq. 5.3 is satisfied for $i=0$.

Let $i\in\mathbb{N}$ and suppose that the conclusion holds for $i-1$ and let

$p_{i-1},q_{i-1}\in J(X_{i-1},\varepsilon_{i-1})\subset
F(X_{i-1},\varepsilon_{i-1})=X_{i}.$

By 5.4 applied to $p=p_{i-1}$, $q=q_{i-1}$, $\varepsilon=\varepsilon_{i-1}$
and $Y=X_{i-1}$, we have

(5.4) $\rho^{X_{i}}(p_{i-1},q_{i-1})=\rho^{X_{i-1}}(p_{i-1},q_{i-1}).$

There are four cases to consider.

First suppose that
$p_{i-1},q_{i-1}\in\mathcal{N}(X_{i},\varepsilon_{i})\subset
J(X_{i},\varepsilon_{i})$. Set $p_{i}=p_{i-1}$ and $q_{i}=q_{i-1}$, so that
eq. 5.4 gives eq. 5.3 in this case.

Now suppose that $p_{i-1}\not\in\mathcal{N}(X_{i},\varepsilon_{i})$ and
$q_{i-1}\in\mathcal{N}(X_{i},\varepsilon_{i})$. Let
$A_{i}\in\mathcal{D}^{*}(X_{i},\varepsilon_{i})$ with $p_{i-1}\in A_{i}$ and
set

$q_{i}=q_{i-1}\in\mathcal{N}(X_{i},\varepsilon_{i})\subset
J(X_{i},\varepsilon_{i}).$

Let $p_{i}\in(\partial A_{i})^{\mathfrak{n}}_{A_{i}}\subset
J(X_{i},\varepsilon_{i})$ be as in 5.5 applied to $Y=X_{i}$, $p=p_{i-1}$ and
$x=q_{i-1}=q_{i}$ so that

$\rho^{X_{i}}(p_{i},q_{i})\leq\rho^{X_{i}}(p_{i-1},q_{i})=\rho^{X_{i}}(p_{i-1},q_{i-1}).$

Applying eq. 5.4 gives eq. 5.3 as required.

If $q_{i-1}\not\in\mathcal{N}(X_{i},\varepsilon_{i})$ and
$p_{i-1}\in\mathcal{N}(X_{i},\varepsilon_{i})$, exchange $p_{i-1}$ and
$q_{i-1}$ and apply the previous case.

Finally, if $p_{i-1},q_{i-1}\not\in\mathcal{N}(X_{i},\varepsilon_{i})$, there
exist $A_{i},B_{i}\in\mathcal{D}^{*}(X_{i},\varepsilon_{i})$ with $p_{i-1}\in
A_{i}$ and $q_{i-1}\in B_{i}$. If $A_{i}=B_{i}$, setting
$p_{i}=q_{i}\in(\partial A_{i})_{A_{i}}$ suffices. Otherwise, (similarly to
the previous two cases) by two applications of 5.5, there exist
$p_{i}\in(\partial A_{i})^{\mathfrak{n}}_{A_{i}}\subset
J(X_{i},\varepsilon_{i})$ such that

$\rho^{X_{i}}(p_{i},q_{i-1})\leq\rho^{X_{i}}(p_{i-1},q_{i-1}),$

and $q_{i}\in(\partial B_{i})^{\mathfrak{n}}_{B_{i}}\subset
J(X_{i},\varepsilon_{i})$ such that

$\rho^{X_{i}}(p_{i},q_{i})\leq\rho^{X_{i}}(p_{i},q_{i-1}).$

By combining these two inequalities with eq. 5.4, we obtain eq. 5.3.

For the in particular statement, since $p_{i},q_{i}\in
J(X_{i},\varepsilon_{i})$, eq. 5.4 implies

$\rho^{X_{i+1}}(p_{i},q_{i})=\rho^{X_{i}}(p_{i},q_{i})\leq\rho^{X_{i-1}}(p_{i-1},q_{i-1})\leq\ldots\leq\rho^{X}(p_{0},q_{0}).$

Also, for any $1\leq j\leq i$, $A\in\mathcal{D}^{*}(X_{j},\varepsilon_{j})$,
and $a,b\in A$, there is a path in $H(j)$ joining $a$ and $b$ (its vertices
are contained in $A$). Since $E(H(j))\subset E(H(i))$, this path is contained
in $H(i)$. Thus, there exists a path from $p_{0}$ to $p_{i}$ in $H(i)$.
Similarly, there exists a path from $q_{0}$ to $q_{i}$ in $H(i)$. ∎

Using the previous Proposition we are able to construct a graph whose
connected components are identical to those produced by
$\operatorname{HDBSCAN}^{*}$ up to a given scale.

######  Lemma 5.8.

Let $i\in\mathbb{N}$. For every $\varepsilon_{i-1}<\alpha\leq\varepsilon_{i}$,

  1. (1)

$E(H(i)_{\alpha})\subset E(G(X,\rho^{X})_{\alpha})$

  2. (2)

if $(p,q)\in E(G(X,\rho^{X})_{\alpha})$ then there is a path in
$H(i)_{\alpha}$ joining $p$ and $q$.

###### Proof.

First notice that, for any $m\in\mathbb{N}$ and
$\varepsilon_{m-1}<\alpha\leq\varepsilon_{m}$,

(5.5)
$H(m)_{\alpha}=\bigcup\limits_{C\in\mathcal{D}^{*}(X_{m},\varepsilon_{m})}G(B_{\varepsilon_{m}}(C),\rho^{B_{\varepsilon_{m}}(C)})_{\alpha}\bigcup
H(m-1).$

We prove Item 1 by induction. For any $C\in\mathcal{D}^{*}(Y,\varepsilon)$ and
$p,q\in C$, we have
$\operatorname{core}_{k}^{B_{\varepsilon}(C)}(p)=core_{k}^{Y}(p)$ and
$\rho^{B_{\varepsilon}(C)}(p,q)=\rho^{Y}(p,q)$. This implies

$E(G(X,\rho^{X})_{\alpha})=\bigcup\limits_{C\in\mathcal{D}^{*}(X,\varepsilon_{1})}E(G(B_{\varepsilon_{1}}(C),\rho^{B_{\varepsilon_{1}}(C)})_{\alpha})=:E(H(1)_{\alpha})$

proving the lemma when $i=1$.

Now assume that item 1 is satisfied for $i=m$ and let
$\varepsilon_{m}<\alpha\leq\varepsilon_{m+1}$. Given
$C\in\mathcal{D}^{*}(X_{m+1},\varepsilon_{m+1})$, the containment

$E(G(B_{\varepsilon_{m+1}}(C),\rho^{B_{\varepsilon_{m+1}}(C)})_{\alpha})\subset
E(G(X,\rho^{X})_{\alpha})$

follows from the fact that $\rho^{X}\leq\rho^{A}$ for any subset $A\subset X$.
Also by the induction hypothesis,

$E(H(m)_{\varepsilon_{m}})\subset E(G(X,\rho^{X})_{\varepsilon_{m}})\subset
E(G(X,\rho^{X})_{\alpha}).$

Thus Equation 5.5 concludes the proof of item 1.

To prove item 2 let $m\in\mathbb{N}$,
$\varepsilon_{m-1}<\alpha\leq\varepsilon_{m}$ and $(p,q)\in
E(G(X,\rho^{X})_{\alpha})$. Applying 5.7 with $p_{0}=p$ and $q_{0}=q$, we
obtain points $p_{m-1},q_{m-1}\in J(X_{m-1},\varepsilon_{m-1})\subset X_{m}$
such that

$\rho^{X_{m}}(p_{m-1},q_{m-1})\leq\rho^{X}(p,q)\leq\alpha\leq\varepsilon_{m}$

and two paths in $E(H(m-1))\subset E(H(m)_{\alpha})$ joining $p_{0}$ and
$p_{m-1}$ and $q_{0}$ and $q_{m-1}$. In particular, $(p_{m-1},q_{m-1})\in
E(G(B_{\varepsilon_{m}}(C),\rho^{B_{\varepsilon_{m}}(C)})_{\alpha})$ for some
$C\in\mathcal{D}^{*}(X_{m},\varepsilon_{m})$ and so, by Equation 5.5,
$(p_{m-1},q_{m-1})\in E(H(m)_{\alpha})$. Combining this edge and the two
previous paths gives the required path from $p$ to $q$ in $H(m)_{\alpha}$.

∎

Finally we define a graph whose connected components coincide with those of
$\operatorname{HDBSCAN}^{*}$ at all scales.

######  Definition 5.9.

For fixed $i\in\mathbb{N}$ define

$\mathcal{F}(i)=G(X_{i+1},\rho^{X_{i+1}})\cup H(i).$

######  Theorem 5.10.

For any $i\in\mathbb{N}$, $\varepsilon_{i-1}<\nu\leq\varepsilon_{i}$,
$\alpha\geq 0$ and any $p\neq q\in X$,

(5.6) $p,q\in V(\Gamma)\text{ for
}\Gamma\in\pi_{0}(H(i)_{\nu})\Longleftrightarrow p,q\in
V(\tilde{\Gamma})\text{ for }\tilde{\Gamma}\in\pi_{0}(G(X,\rho)_{\nu}).$

and

(5.7) $p,q\in V(\Gamma)\text{ for
}\Gamma\in\pi_{0}(\mathcal{F}(i)_{\alpha})\Longleftrightarrow p,q\in
V(\tilde{\Gamma})\text{ for }\tilde{\Gamma}\in\pi_{0}(G(X,\rho)_{\alpha}).$

###### Proof.

Fix $i\in\mathbb{N}$. Equation 5.6 follows immediately from 5.8.

To prove Equation 5.7 when $\alpha\leq\varepsilon_{i}$, first notice that

$E(G(X_{i+1},\rho^{X_{i+1}})_{\alpha})\subset E(H(i)_{\alpha}).$

Consequently, the result follows from Equation 5.6.

Now suppose that $\alpha>\varepsilon_{i}$ and note that

$\mathcal{F}(i)_{\alpha}=G(X_{i+1},\rho^{X_{i+1}})_{\alpha}\cup H(i).$

Since $\rho^{X_{i+1}}\leq\rho^{X}$ wherever they are both defined,

$E(G(X_{i+1},\rho^{X_{i+1}})_{\alpha})\subset E(G(X,\rho^{X})_{\alpha}).$

By 5.8 item 1,

$E(H(i)_{\alpha})\subset E(G(X,\rho^{X})_{\alpha}).$

Therefore

$E(\mathcal{F}_{\alpha})\subset E(G(X,\rho^{X})_{\alpha})$

and so the first implication in Equation 5.7 holds.

Finally, assume $(p_{0},q_{0})\in E(G(X,\rho^{X})_{\alpha})$. Let
$p_{i},q_{i}\in J(X_{i},\varepsilon_{i})\subset X_{i+1}$ be as in 5.7 so that

$\rho^{X_{i+1}}(p_{i},q_{i})\leq\rho^{X}(p_{0},q_{0})\leq\alpha.$

In particular, $(p_{i},q_{i})\in E(G(X_{i+1},\rho^{X_{i+1}})_{\alpha})$. Since
there exist paths in $H(i)$ connecting $p_{0}$ to $p_{i}$ and $q_{0}$ to
$q_{i}$, there exists a path from $p_{0}$ to $q_{0}$ in
$\mathcal{F}(i)_{\alpha}$. Consequently, the second implication in Equation
5.7 holds.

∎

##  6\. Example: $\operatorname{S-DBSCAN}^{*}$ to cluster building data of the
US

To illustrate our methods, we cluster datasets of building footprints in the
United States. We use the Microsoft Open Buildings dataset [12] which is
freely available.

We cluster two different datasets using $\operatorname{S-DBSCAN}^{*}$ and
sklearn.$\operatorname{DBSCAN}^{*}$ and compare the runtime in each instance.
Our first dataset corresponds to the state of Utah which consists of 1,004,734
data points. The second dataset corresponds to the entire United States and
has size 124,828,547. For each of these datasets, we produce the
$\operatorname{DBSCAN}^{*}$ clusters for $\varepsilon=$ 3,000 and
$\varepsilon=$12,000, each with $k=$1,900. In each instance, the run time is
recorded in seconds.

These experiments were performed using an Intel Xeon CPU running at 2.30GHz
with 64GB of memory on 64-bit Ubuntu 20.04.3 LTS.

Tables 1 and 2 show the runtime of $\operatorname{S-DBSCAN}^{*}$ and
sklearn.$\operatorname{DBSCAN}^{*}$, the resulting number of core points
$|\mathcal{C}|$, number of clusters $|\mathcal{D}^{*}|$ and maximum and mean
cluster sizes. It should be noted that the implementation of
$\operatorname{S-DBSCAN}^{*}$ is not optimized and relies purely on python
code which negatively affects performance. On the other hand,
sklearn.$\operatorname{DBSCAN}^{*}$ is optimized for speed and not memory
efficiency. Five workers are used for all runs of the
$\operatorname{S-DBSCAN}^{*}$ algorithm.

| $\varepsilon$ | Time | $|\mathcal{C}|$ | $|\mathcal{D}^{*}|$ | $\operatorname{max}|C|$ | mean $|C|$
---|---|---|---|---|---|---
$\operatorname{S-DBSCAN}^{*}$ | 2,000 | 34 | 670,611 | 30 | 305,024 | 22,354
$\operatorname{S-DBSCAN}^{*}$ | 8,000 | 6 | 877,112 | 19 | 741,710 | 46,164
$\operatorname{DBSCAN}^{*}$ | 2,000 | 91 | 670,611 | 30 | 305,024 | 22,354
$\operatorname{DBSCAN}^{*}$ | 8,000 | - | - | - | - | -
Table 1. $\operatorname{DBSCAN}^{*}$ for Utah | $\varepsilon$ | Time | $|\mathcal{C}|$ | $|\mathcal{D}^{*}|$ | $\operatorname{max}|C|$ | mean $|C|$
---|---|---|---|---|---|---
$\operatorname{S-DBSCAN}^{*}$ | 3,000 | 1,752 | 79,919,086 | 3,244 | 6,237,660 | 24,636
$\operatorname{S-DBSCAN}^{*}$ | 12,000 | 961 | 120,976,900 | 429 | 94,517,493 | 281,997
$\operatorname{DBSCAN}^{*}$ | 3,000 | - | - | - | - | -
$\operatorname{DBSCAN}^{*}$ | 12,000 | - | - | - | - | -
Table 2. $\operatorname{DBSCAN}^{*}$ for USA

In Table 1 it can be seen that, for the state of Utah,
$\operatorname{S-DBSCAN}^{*}$ has a lower runtime than
sklearn.$\operatorname{DBSCAN}^{*}$ for $\varepsilon=2,000$ . Furthermore, for
$\varepsilon=$ 8,000 sklearn.$\operatorname{DBSCAN}^{*}$ ran out of memory. In
contrast, the runtime of $\operatorname{S-DBSCAN}^{*}$ _decreases_ for the
larger value of $\varepsilon$. This is because more cubes of the partition
become classified as dense, and consequently, the required number of distance
calculations significantly decreases.

Table 2 shows the running time of $\operatorname{S-DBSCAN}^{*}$ and
sklearn.$\operatorname{DBSCAN}^{*}$ when clustering building data of the US.
For all parameter choices sklearn.$\operatorname{DBSCAN}^{*}$ runs out of
memory. Similarly to the Utah dataset, for the larger value of $\varepsilon$,
the $\operatorname{S-DBSCAN}^{*}$ runtime decreases.

##  7\. Example: $\operatorname{S-HDBSCAN}^{*}$ to cluster building data of
the US

Recall that, for a fixed value of $\varepsilon$, the connected components of
$\operatorname{DBSCAN}^{*}$ (with more than one element) coincide with the
connected components of $G(X,\rho)_{\varepsilon}$, used to construct the
$\operatorname{HDBSCAN}^{*}$ clusters, see eq. 2.1. This fact allows us to
construct $\operatorname{HDBSCAN}^{*}$ in steps by first separating the data
into $\operatorname{DBSCAN}^{*}$ clusters and then obtaining the connected
components of $G(X,\rho)_{\varepsilon}$ restricted to each
$\operatorname{DBSCAN}^{*}$ cluster.

More explicitly the steps of $\operatorname{S-HDBSCAN}^{*}$ are:

  1. (1)

Cluster $X_{1}=X$ using $\operatorname{S-DBSCAN}^{*}$ with an initial choice
of $\varepsilon_{1}$.

  2. (2)

For each cluster $C\in\mathcal{D}^{*}(X_{1},\varepsilon_{1})$, we calculate
the reachability distance $\tilde{\rho}$ using _only_ points in
$C^{\mathfrak{m}}$ and obtain the weighted graph
$G(C,\tilde{\rho})_{\varepsilon_{1}}$.

  3. (3)

Combine an extended boundary of each cluster and
$\mathcal{N}(X_{1},\varepsilon_{1})$ to form $X_{2}$. Repeat steps 1 and 2 on
$X_{2}$ for a choice of $\varepsilon_{2}>\varepsilon_{1}$.

  4. (4)

Combine an extended boundary of each cluster in
$\mathcal{D}^{*}(X_{2},\varepsilon_{2})$ and
$\mathcal{N}(X_{2},\varepsilon_{2})$ to form $X_{3}$.

  5. (5)

Calculate $\tilde{\rho}$ in $X_{3}$ and the weighted graph
$G(X_{3},\tilde{\rho})$.

  6. (6)

Combine all weighted graphs produced at each step to produce our final graph
$G$. 5.10 guarantees that, for any $\varepsilon>0$, the connected components
of $G_{\varepsilon}$ agree with those of $G(X,\rho)_{\varepsilon}$.

This procedure can be iterated any number of times, we use $3$ for
illustration only. At each step, each $\operatorname{S-DBSCAN}^{*}$ cluster is
processed independently and in parallel.

One of the fastest implementations of $\operatorname{HDBSCAN}^{*}$ is
sklearn.$\operatorname{HDBSCAN}^{*}$, [10]. The speed of this algorithm comes
from combining several high performance algorithms used to optimise the steps
of $\operatorname{HDBSCAN}^{*}$ with high time complexity. A central tool used
in these optimisations are _$k$ -d trees_, which are state-of-the-art in
nearest neighbour type searches.

A $k$-d tree partitions $\mathbb{R}^{n}$ into non-overlapping regions, indexed
by a given _reference tree_ $\tau_{r}\subset X\subset\mathbb{R}^{n}$. Given a
_query point_ $p\in\tau_{q}\subset X$, one can use the partition to
efficiently find the nearest neighbours in $\tau_{r}$ to $p$. To further
optimise performance one can build a second partition based on the _query
tree_ $\tau_{q}$.

In [10], query and reference trees are used to calculate
$\operatorname{core}_{k}$ and a _minimum spanning tree_. In practice, given a
weighted graph $H$, to obtain the connected components of $H_{\varepsilon}$,
for every $\varepsilon$, a standard and efficient technique is to store a
minimum spanning tree of $H$. This is done using the dual tree Borůvka
algorithm.

The performance and memory usage of calculating $\operatorname{core}_{k}$ and
the Borůvka algorithm worsens as the size of the reference and query trees
increase. In fact, in [11] it can be seen that clustering very large datasets
using sklearn.$\operatorname{HDBSCAN}^{*}$ becomes infeasible. An underlying
reason for this is that these trees are constructed using the entire dataset.
A common solution is to process the points in the query tree in batches.
However, given the size of our dataset and our choice of $k$, the reference
trees are simply too large and batch processing the query trees alone does not
make the task tangible.

To solve this issue $\operatorname{S-HDBSCAN}^{*}$ partitions $X$ to create
many smaller trees. The difficulty lies in partitioning the reference trees
whilst ensuring that the correct $\operatorname{core}_{k}$ and the correct
weights of a spanning graph are produced. Since $\operatorname{DBSCAN}^{*}$
clusters naturally partition $X$ they can be used to achieve this. We take
this further and iteratively use $\operatorname{DBSCAN}^{*}$ clusters to
incrementally construct a spanning graph. After each iteration we can remove
points from the dataset that are now redundant, namely the interiors of
clusters of previous iterations. Moreover, the calculations for each cluster
are independent and are processed in parallel. A crucial point is that, the
minimum spanning trees of each cluster are combined into a graph with a
relatively small amount of edges that, by 5.10, allows us to recover the
connected components required for $\operatorname{HDBSCAN}^{*}$.

In practice, from one iteration to the next we only add edges $(p,q)$ with
weight $\varepsilon_{i-1}<\omega(p,q)\leq\varepsilon_{i}$ because this is
sufficient to complete a small spanning graph. Further, we only add an edge
$(p,q)$ if $p$ and $q$ were not in the same connected component in a previous
iteration.

We also note that $\operatorname{S-HDBSCAN}^{*}$ can be used recursively: if a
particular cluster cannot be processed directly, the same approach can be
applied to the cluster.

To demonstrate our theory, we built an implementation of
$\operatorname{S-HDBSCAN}^{*}$ and ran it on the same example datasets than
Section 6 (using the same computer). Since this was for demonstration purposes
only it consists of unoptimised Python scripts. As a point of comparison, we
also ran sklearn.$\operatorname{HDBSCAN}^{*}$.

The results are shown in Table 3, with time measured in seconds. The three
final columns refer to the final $\operatorname{HDBSCAN}^{*}$ clusters, after
the hierarchical clusters have been scored across all scales. In the case of
Utah, where sklearn.$\operatorname{HDBSCAN}^{*}$ completes, we see that
$\operatorname{S-HDBSCAN}^{*}$ has a comparable running time. However for the
entire dataset, sklearn.$\operatorname{HDBSCAN}^{*}$ does not terminate,
whereas $\operatorname{S-HDBSCAN}^{*}$ does.

|  | Time | $|\mathcal{H}^{*}|$ | $\operatorname{max}|C|$ | mean $|C|$
---|---|---|---|---|---
$\operatorname{S-HDBSCAN}^{*}$ | Utah | 649 | 40 | 307,149 | 20,386
$\operatorname{HDBSCAN}^{*}$ | Utah | 415 | 40 | 307,149 | 20,386
$\operatorname{S-HDBSCAN}^{*}$ | USA | 89,670 | 5,222 | 44,515,947 | 23,904
$\operatorname{HDBSCAN}^{*}$ | USA | - | - | - | -
Table 3. Final $\operatorname{HDBSCAN}^{*}$ clusters, $k=m=1900$

###  7.1. Tree sizes at each iteration

We now illustrate the reduction of the query and reference trees, at each
iteration, granted by our theory. For the first iteration, we use
$\varepsilon_{1}=3,000$ and $\operatorname{S-DBSCAN}^{*}$ produces 3,244
clusters of the entire dataset $|X_{1}|=124,828,547$. For the second iteration
we use $\varepsilon_{2}=12,000$ and $\operatorname{S-DBSCAN}^{*}$ produces 429
clusters of $|X_{2}|=110,943,965$. The last step of the process uses
$|X_{3}|=66,145,965$. We refer the reader to table 2 for details of the
$\operatorname{DBSCAN}^{*}$ clusters of $X_{1}$, relevant to iterations 1 and
2.

Since our datasets are contained in $\mathbb{R}^{2}$, we have
$\mathfrak{n}=1$, $\mathfrak{m}=3$ and $\mathfrak{N}=4$. Table 4 contains the
size distribution of the query tree $\tau_{q}(C)=C$ and reference tree
$\tau_{r}(C)$ of clusters $C\in\mathcal{D}^{*}(X_{1},\varepsilon_{1})$ used to
calculate $\operatorname{core}_{k}$. By 4.4, it suffices to take
$\tau_{r}(C)=C^{3}$ as a reference tree to calculate
$\operatorname{core}_{k}(p)$ for points $p\in C$. The query and reference tree
required to build the minimum spanning tree of $C$ is $C$. Table 5 contains
the sizes of these trees for the 10 largest clusters.

| $\tau_{q}(C)$ | $\tau_{r}(C)$
---|---|---
mean | 24,636 | 29,059
std | 179,487 | 190,162
min | 1 | 1,956
25% | 1,908 | 3,702
50% | 3,204 | 5,692
75% | 7,080 | 10,742
max | 6,237,660 | 6,745,531
count | 3,244
Table 4. First iteration: query and reference tree sizes. $\tau_{q}(C)$ | $\tau_{r}(C)$ | $\tau_{q}(C)$ | $\tau_{r}(C)$
---|---|---|---
6,237,660 | 6,745,531 | 5.00% | 5.40%
4,055,724 | 4,090,661 | 3.25% | 3.28%
3,089,031 | 3,204,504 | 2.47% | 2.57%
2,162,629 | 2,285,680 | 1.73% | 1.83%
1,975,669 | 2,202,149 | 1.58% | 1.76%
1,721,895 | 1,814,330 | 1.38% | 1.45%
1,680,441 | 1,760,533 | 1.35% | 1.41%
1,661,300 | 1,740,421 | 1.33% | 1.39%
1,453,234 | 1,574,444 | 1.16% | 1.26%
1,450,173 | 1,556,962 | 1.16% | 1.25%
25,487,756 | 26,975,215 | 20.42% | 21.61%
Table 5. First iteration: tree sizes of largest clusters and their percentage
of the whole dataset.

This iteration of $\operatorname{S-HDBSCAN}^{*}$ processes
$\operatorname{core}_{k}$ and the minimum spanning tree of approximately 80
million points, see Table 2. Tables 4 and 5 show that this is achieved using
reference trees with less than 7 million points, 5% of the entire dataset. In
fact, about 55 million of these points (44% of the dataset) are processed with
query and reference trees smaller than 1.5 million points. Table 4 shows that
75% of the clusters have reference trees with less than 11,000 points, four
orders of magnitude smaller than the entire dataset.

This massive reduction of the sizes of these reference trees is what makes it
possible for us to cluster a dataset of this size, which we do in 24hrs. In
data of this magnitude, it is appropriate to model the complexity of
sklearn.$\operatorname{HDBSCAN}^{*}$ algorithm as $O(n^{2})$, [11]. Under this
assumption, dividing the dataset into such small pieces is a significant
improvement.

There are two possible ways to continue constructing a spanning graph from
scales $\varepsilon_{1}$ to $\varepsilon_{2}$. The first is to repeat the
first iteration with $\varepsilon_{2}$ but only processing noise points from
the first iteration. That is, for each
$C\in\mathcal{D}^{*}(X_{1},\varepsilon_{2})$ the trees
$\tau^{c}_{q}(C)=C\cap\mathcal{N}(X_{1},\varepsilon_{1})$ and
$\tau^{c}_{r}(C)=C^{3}$ calculate $\operatorname{core}_{k}$. The tree
$\tau^{m}=C$ is used as a query and reference tree to calculate the minimum
spanning tree. The sizes of these trees for the 10 largest clusters are
contained in the first three columns of Table 6. Table 7 shows the same values
as a percentage of the size of the entire dataset.

We consider an alternative way that reduces the sizes of the trees even
further by first removing an interior from each cluster of the first
iteration. More precisely, set $X_{2}=F(X_{1},\varepsilon_{1})$ as in 5.3
(keeping only the 4th extension of the boundary of each cluster) and, for each
$C\in\mathcal{D}^{*}(X_{2},\varepsilon_{2})$, set
$\tau^{c}_{q}(C)=C\cap\mathcal{N}(X_{1},\varepsilon_{1})$ and
$\tau^{c}_{r}(C)=C^{3}$. To calculate the minimum spanning tree, the query and
reference trees agree and equal $\tau^{m}(C)=C\cap J(X_{1},\varepsilon_{1})$
(as defined in 5.3, keeping only the 1st extension of the boundary of previous
iteration clusters). The sizes of these trees for the 10 largest clusters are
shown in columns 1, 4 and 5 of Tables 6 and 7.

By 5.10 and 5.5, the spanning graphs produced by either of these methods are
equivalent to $G(X,\rho)_{\varepsilon}$ for any
$\varepsilon\leq\varepsilon_{2}$.

Both | With interior | Without
---|---|---
$\tau^{c}_{q}(C)$ | $\tau^{c}_{r}(C)$ | $\tau^{m}(C)$ | $\tau^{c}_{r}(C)$ | $\tau^{m}(C)$
36,039,363 | 95,728,364 | 94,517,493 | 85,297,674 | 63,373,014
984,007 | 10,580,949 | 10,495,579 | 8,697,199 | 4,699,615
806,604 | 3,575,609 | 3,530,441 | 3,157,678 | 2,039,457
242,434 | 1,553,354 | 1,519,012 | 1,241,211 | 740,889
208,895 | 1,530,487 | 1,518,726 | 1,157,362 | 537,396
170,270 | 813,466 | 797,390 | 740,581 | 425,973
105,852 | 578,156 | 576,816 | 530,175 | 362,543
89,095 | 545,200 | 498,704 | 477,756 | 356,340
83,919 | 477,756 | 443,445 | 391,985 | 231,348
67,557 | 454,799 | 441,073 | 369,047 | 222,249
38,797,996 | 115,838,140 | 114,338,679 | 102,060,668 | 72,988,824
Table 6. Second iteration: tree sizes of largest clusters Both | With interior | Without
---|---|---
$\tau^{c}_{q}(C)$ | $\tau^{c}_{r}(C)$ | $\tau^{m}(C)$ | $\tau^{c}_{r}(C)$ | $\tau^{m}(C)$
28.87% | 76.69% | 75.72% | 68.33% | 50.77%
0.79% | 8.48% | 8.41% | 6.97% | 3.76%
0.65% | 2.86% | 2.83% | 2.53% | 1.63%
0.19% | 1.24% | 1.22% | 0.99% | 0.59%
0.17% | 1.23% | 1.22% | 0.93% | 0.43%
0.14% | 0.65% | 0.64% | 0.59% | 0.34%
0.08% | 0.46% | 0.46% | 0.42% | 0.29%
0.07% | 0.44% | 0.40% | 0.38% | 0.29%
0.07% | 0.38% | 0.36% | 0.31% | 0.19%
0.05% | 0.36% | 0.35% | 0.30% | 0.18%
31.08% | 92.80% | 91.60% | 81.76% | 58.47%
Table 7. Second iteration: percentage of total dataset of tree sizes of
largest clusters

In this second iteration there is a very large cluster containing 95 million
points, which has large reference and query trees as a result. However, this
cluster consists only of 37 million new points and, at most, 3,244 clusters
already formed in the iteration 1. Consequently, our calculations of
$\operatorname{core}_{k}$ restrict to the new 37 million points and the
Borůvka algorithm only requires 37 million iterations, instead of 95 million,
to complete the spanning graph of this cluster. We can see in Tables 6 and 7
that the trees corresponding to the remaining 428 $\operatorname{DBSCAN}^{*}$
clusters have a comparable distribution to those of the first iteration.

Table 7 shows that removing the interior of the clusters from iteration 1
reduces the size of the reference trees by about 10% of the total dataset when
calculating $\operatorname{core}_{k}$ and by 30% when calculating the spanning
graphs.

For the final iteration, we calculate core distances with a single query tree
$\tau^{c}_{q}=\mathcal{N}(X_{1},\varepsilon_{2})$. If we keep the interiors of
clusters from the previous iterations, we must use $\tau^{c}_{r}=X_{1}$ as a
reference tree. However, we can reduce the reference tree significantly by
removing the interiors of clusters giving $\tau^{c}_{r}=X_{3}$. For the
minimum spanning tree we have query and reference trees equal to
$\tau^{m}=X_{1}$ in the first case and $\tau^{m}=J(X_{2},\varepsilon_{2})$ in
the second case.

Both | With interior | Without
---|---|---
$\tau^{c}_{q}$ | $\tau^{c}_{r}$ | $\tau^{m}$ | $\tau^{c}_{r}$ | $\tau^{m}$
3,851,647 | 124,828,547 | 124,828,548 | 66,145,965 | 28,931,402
3.09% | 100% | 100% | 52.99% | 23.18%
Table 8. Final iteration: tree sizes and their percentage of the whole
dataset.

This iteration highlights the importance of removing the interiors of previous
iterations. In Table 8 we see that calculating $\operatorname{core}_{k}$ has a
relatively small query tree and that the impact of removing the interiors of
previous iterations is significant, reducing the size of the reference tree by
almost 50%. For calculating the spanning graph, the reduction from removing
the interiors is even greater and reduces the size of both the query and the
reference trees by almost 80% of the total dataset.

###  7.2. Final remarks

There is a lot of potential to implement $\operatorname{S-HDBSCAN}^{*}$ in
parallel and to employ sophisticated algorithms that optimise its performance.
Moreover, the number of points that we removed in this final iteration
motivates further study in how to optimally choose the $\varepsilon_{i}$ to
have a greater impact by removing larger interiors in each iteration.

## References

  * [1] J. Ajayakumar, A. Curtis, and J. Curtis.  A clustering environment for real-time tracking and analysis of COVID-19 case clusters.  In Proceedings of the 2nd ACM SIGSPATIAL International Workshop on Spatial Computing for Epidemiology (SpatialEpi 2021), SpatialEpi’21, page 1–9, New York, NY, USA, 2021. Association for Computing Machinery.  doi:10.1145/3486633.3491091.
  * [2] F. Azimlu, S. Rahnamayan, and M. Makrehchi.  House price prediction using clustering and genetic programming along with conducting a comparative study.  In Proceedings of the Genetic and Evolutionary Computation Conference Companion, page 1809–1816, New York, NY, USA, 2021. Association for Computing Machinery.  doi:10.1145/3449726.3463141.
  * [3] R. J. G. B. Campello, D. Moulavi, and J. Sander.  Density-based clustering based on hierarchical density estimates.  In Advances in Knowledge Discovery and Data Mining, pages 160–172. Springer Berlin Heidelberg, 2013.  doi:10.1007/978-3-642-37456-2_14.
  * [4] M. de Berg, A. Gunawan, and M. Roeloffzen.  Faster DBSCAN and HDBSCAN in low-dimensional Euclidean spaces.  International Journal of Computational Geometry & Applications, 29(01):21–47, 2019.  doi:10.1142/S0218195919400028.
  * [5] G. Dong, L. Cai, D. Datta, S. Kumar, L. E. Barnes, and M. Boukhechba.  Influenza-like symptom recognition using mobile sensing and graph neural networks.  In Proceedings of the Conference on Health, Inference, and Learning, CHIL ’21, page 291–300, New York, NY, USA, 2021. Association for Computing Machinery.  doi:10.1145/3450439.3451880.
  * [6] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu.  A density-based algorithm for discovering clusters in large spatial databases with noise.  In Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, KDD’96, page 226–231. AAAI Press, 1996\.  doi:10.5555/3001460.3001507.
  * [7] S. F. Galán.  Comparative evaluation of region query strategies for DBSCAN clustering.  Information Sciences, 502:76–90, 2019.  doi:10.1016/j.ins.2019.06.036.
  * [8] Y. He, H. Tan, W. Luo, S. Feng, and J. Fan.  MR-DBSCAN: a scalable MapReduce-based DBSCAN algorithm for heavily skewed data.  Frontiers of Computer Science, 8(1):83–99, dec 2013.  doi:10.1007/s11704-013-3158-3.
  * [9] J. Jackson, A. Qiao, and E. P. Xing.  Scaling HDBSCAN clustering with kNN graph approximation.  Proceedings of the SysML Conference 2018, pages 14–16, 2018.  URL: https://mlsys.org/Conferences/doc/2018/105.pdf.
  * [10] L. McInnes and J. Healy.  Accelerated hierarchical density based clustering.  In 2017 IEEE International Conference on Data Mining Workshops (ICDMW), pages 33–42, 2017.  doi:10.1109/ICDMW.2017.12.
  * [11] L. McInnes, J. Healy, and S. Astels.  Benchmarking performance and scaling of python clustering algorithms, 2016\.  URL: https://hdbscan.readthedocs.io/en/latest/performance_and_scalability.html#.
  * [12] Microsoft.  Microsoft maps: US building footprints, 2018.  URL: https://github.com/Microsoft/USBuildingFootprints.
  * [13] A. C. A. Neto, J. Sander, R. J. Campello, and M. A. Nascimento.  Efficient computation of multiple density-based clustering hierarchies.  In 2017 IEEE International Conference on Data Mining (ICDM), pages 991–996. IEEE, 2017.  doi:10.1109/ICDM.2017.127.
  * [14] M. M. A. Patwary, D. Palsetia, A. Agrawal, W.-k. Liao, F. Manne, and A. Choudhary.  A new scalable parallel DBSCAN algorithm using the disjoint-set data structure.  In SC’12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, pages 1–11. IEEE, 2012\.  doi:10.1109/SC.2012.9.
  * [15] T. Statham, L. Wolf, and S. Fox.  Applications of gridded population datasets: delineating urban areas.  In Proceedings of the 28th Geographical Information Science Research UK Conference, 2020.  URL: http://london.gisruk.org/gisruk2020_proceedings/GISRUK2020_paper_63.pdf.
  * [16] N. Wang and Y. Zhang.  Adaptive and fast image superpixel segmentation approach.  Image and Vision Computing, 116:104315, 2021.  doi:10.1016/j.imavis.2021.104315.
  * [17] Y. Wang, S. Yu, Y. Gu, and J. Shun.  Fast parallel algorithms for Euclidean minimum spanning tree and hierarchical spatial clustering.  Proceedings of the ACM SIGMOD International Conference on Management of Data, pages 1982–1995, 2021.  doi:10.1145/3448016.3457296.
  * [18] N. Zhao, J. Chen, X. Peng, H. Wang, X. Wu, Y. Zhang, Z. Chen, X. Zheng, X. Nie, G. Wang, Y. Wu, F. Zhou, W. Zhang, K. Sui, and D. Pei.  Understanding and handling alert storm for online service systems.  In 2020 IEEE/ACM 42nd International Conference on Software Engineering: Software Engineering in Practice (ICSE-SEIP), pages 162–171, 2020\.