# The Hierarchical Structure of Galactic Haloes:
Generalised N-Dimensional Clustering with CluSTAR-ND

William H. Oliver1, Pascal J. Elahi2, and Geraint F. Lewis1
1Sydney Institute for Astronomy, School of Physics A28, The University of
Sydney, NSW, 2006, Australia
2Pawsey Supercomputing Research Centre, 1 Bryce Avenue, Kensington, WA, 6151,
Australia E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present CluSTAR-ND, a fast hierarchical galaxy/(sub)halo finder that
produces Clustering Structure via Transformative Aggregation and Rejection in
N-Dimensions. It is designed to improve upon Halo-OPTICS – an algorithm that
automatically detects and extracts significant astrophysical clusters from the
3D spatial positions of simulation particles – by decreasing run-times,
possessing the capability for metric adaptivity, and being readily applicable
to data with any number of features. We directly compare these algorithms and
find that not only does CluSTAR-ND produce a similarly robust clustering
structure, it does so in a run-time that is at least $3$ orders of magnitude
faster. In optimising CluSTAR-ND’s clustering performance, we have also
carefully calibrated $4$ of the $7$ CluSTAR-ND parameters which – unless
specified by the user – will be automatically and optimally chosen based on
the input data. We conclude that CluSTAR-ND is a robust astrophysical
clustering algorithm that can be leveraged to find stellar satellite groups on
large synthetic or observational data sets.

###### keywords:

galaxies: structure – galaxies: star clusters: general – methods: data
analysis – methods: statistical

††pubyear: 2022††pagerange: The Hierarchical Structure of Galactic Haloes:
Generalised N-Dimensional Clustering with CluSTAR-ND–The Hierarchical
Structure of Galactic Haloes:  Generalised N-Dimensional Clustering with
CluSTAR-ND

##  1 Introduction

The process of identifying clusters – often referred to as clustering111Note
that we use this term throughout this paper to mean both the general process
of finding clusters and the resultant set of classifications from this process
– depending on the context. As such, our use of this term is much broader than
its typical use within astro- and cosmo-related fields, e.g. referring
specifically to large-scale-structure. – from a data set has been an ongoing
data mining problem within the field of machine learning. The function of any
algorithm tasked with undertaking this problem is to determine a set of
statistically coherent groups within the intrinsic feature space of the data.
Many such algorithms exist for this purpose, however due to the innate
subjectivity of the definition of a cluster it is not always obvious as to
which of these may be useful for clustering of a given type.

The groupings found by a clustering algorithm can be categorised into one or
more of the typical models, such as; centroid-based (e.g. K-Means; MacQueen et
al., 1967; Lloyd, 1982), distribution-based (e.g. EM; Dempster et al., 1977),
density-based (e.g. DBSCAN; Ester et al., 1996), and others. The way in which
a clustering algorithm partitions the data is also of importance and provides
a similar categorisation. These algorithms may return a flat or hierarchical
clustering such that clusters are mutually exclusive or can be proper
sub/supersets of one another – this is one of the differences between DBSCAN
and HDBSCAN (Campello et al., 2015) for example. In addition to these, an
algorithm may also return an overlapping set of clusters. A clustering may
place all points within clusters (e.g. K-Means) or more commonly it may leave
some points out of clusters classifying them as noise or outliers to the
clusters that have been predicted. Moreover, these algorithms may return a
hard or soft – also referred to as fuzzy – clustering. This distinguishment
clarifies whether the nature of a point being within a cluster is binary-based
or probability-based (e.g Fuzzy C-Means; Dunn, 1973). Since each of these
separate class systems are codependent – in that an algorithm may be
classified by multiple at a time – it suffices to say that choosing an
appropriate clustering algorithm for the problem at hand is non-trivial.

In cosmology, a typical description of a galaxy and its halo is any spatial
overdensity that is denser than the critical – or mean – density of the
universe by some factor $\Delta$. Two common ways to identify such
overdensities from both cosmological simulations and observed data sets are
via the Spherical-Overdensity method (SO; Press & Schechter, 1974) and the
Friends-Of-Friends algorithm (FOF; Davis et al., 1985). By themselves, these
algorithms perform a density-based flat and hard clustering with noise which
is well suited to the intended definition of galactic haloes – and as an
additional constraint, the SO method can also be considered part of the
distribution-based family due to ensuring that clusters adhere to a particular
volumetric shape.

It is a primary prediction of the $\Lambda$CDM cosmological model that
galaxies are formed hierarchically via continual accretion and merger events
(White & Rees, 1978; Kauffmann et al., 1993; Ghigna et al., 1998). Depending
upon the conditions of the satellite at the time of the merger as well as the
ongoing conditions of the host halo (Bullock & Johnston, 2005; Johnston et
al., 2008), the particles within these infalling groups may become mixed in
with, and indistinguishable from, the surrounding halo222The phase-space
volume of infalling groups is conserved in a fully collisionless Newtonian
gravity simulation.. This prediction is mostly confirmed in nature barring a
few exceptions such as the missing satellite problem (Klypin et al., 1999;
Moore et al., 1999; Reed et al., 2005; Tollerud et al., 2008; Springel et al.,
2008; Ishiyama et al., 2013) and the core-cusp problem (Flores & Primack,
1994; Moore, 1994; Van Den Bosch et al., 2000). As such, an appropriate
flavour of clustering algorithm for application to within galactic haloes is
that of the density-based hierarchical type with noise. Many such algorithms
exist for specifically this purpose, including but not limited to; SUBFIND
(Springel et al., 2001), AHF (Knollmann & Knebe, 2009), ROCKSTAR (Behroozi et
al., 2012), VELOCIraptor (Elahi et al., 2019), and more. Comparisons of these
and others tend to show the outputs of halo-finders to be similar (Knebe et
al., 2011, 2013).

The above halo-finders mostly rely on the SO method and/or the FOF algorithm
which both restrict the clustering to certain types of clusters. The varied
shapes and density profiles of subhaloes often mean that these base algorithms
need to be applied iteratively in order to find a range of subhalo types. In
doing so, the run-times become larger although this does also allow some of
these halo-finders to return hierarchical clusterings and even utilise
adaptive metrics. Despite these measures of fine-tuning, SO-based halo-finders
struggle to find curved structures such as streams and FOF-based halo-finders
can still fall victim to point-point noise and fail to detect sparsely
populated structures. The galaxy/(sub)halo finder Halo-OPTICS (Oliver et al.,
2020) is based on OPTICS (Ankerst et al., 1999), a generalisation of DBSCAN,
and is more robust to point-point noise. Coupled with an adaptive metric,
Halo-OPTICS could be a powerful astrophysical clustering algorithm if not for
the fact that the run-time of OPTICS can become unmanageable when operating in
this way. Alternatives typically only use a set of nearest neighbours to
perform locally adaptive metric computation. One such algorithm, the subhalo-
finder EnLink (Sharma & Johnston, 2009), uses an entropy-based locally
adaptive metric in combination with a group finder that is also robust to
point-point noise in order to find a range of subhalo types. VELOCIraptor and
HSF (Maciejewski et al., 2009) also can attempt to generate locally adaptive
metrics to identify clusters.

Having a fast clustering algorithm is not simply a matter of convenience. It
is critical when clustering over extensive data sets with large feature
spaces. It is also of particular importance when used to find a fuzzy
clustering of an astrophysical data set. Unlike the output from the Fuzzy
C-Means algorithm, which is a fuzzy clustering of a hard data set, the fuzzy
clustering found in this context is itself from a fuzzy data set and needs to
have propagated the uncertainties of a point’s features into the probability
of that point’s membership within a given cluster as well as the probability
of that cluster’s existence. Solving this problem typically requires taking
many samplings of the data and comparing each of the clusterings therefrom
(e.g. Fuentes et al., 2017; Malhan et al., 2022). In order to be able to
perform such a task whilst also maintaining a high-calibre clustering power,
it is becoming necessary to use an algorithm that is: density-based;
hierarchical; robust to point-point noise; equipped with adaptive metric
capabilities; applicable to data sets with any number of features; and
demonstrates exceptionally modest run-times.

We present CluSTAR-ND, a derivative of Halo-OPTICS that possesses the above
qualities, and apply it to synthetic survey data produced by Galaxia (Sharma
et al., 2011) in order to formalise the process of clustering n-dimensional
data sets of galactic haloes. We first provide a relevant outline of the Halo-
OPTICS algorithm (Sec. 2). We then describe the CluSTAR-ND algorithm in Sec.
3. In doing so we give an overview of the concept and motivation behind the
algorithm (Sec. 3.1), describe the root-haloes produced (Sec. 3.2), and the
means by which the substructure therein is found (Secs. 3.3 – 3.6). In Sec. 4
we summarise the details of the synthetic galaxies from Galaxia. Following
this we make a direct comparison between CluSTAR-ND and the Halo-OPTICS
algorithm in Sec. 5. In Sec. 6 we discuss the influence of the algorithm’s
parameters and optimise them by training them on the galaxies from Galaxia. We
then discuss the contrast in the information content available from the
clusterings produced from different clustering scenarios in Sec. 7. Finally,
we present our conclusions and directions for future work in Sec. 8.

##  2 An Overview of the Halo-OPTICS Algorithm

In order to understand CluSTAR-ND it is first necessary to grasp the algorithm
Halo-OPTICS (Oliver et al., 2020). Halo-OPTICS is a hierarchical
galaxy/(sub)halo finder that can be used for the density-based determination
of astrophysical clusters in a 3D spatial data set via a global distance
metric, i.e. the Euclidean distance. It is an extension of the well-known
hierarchical clustering algorithm, OPTICS (Ankerst et al., 1999) – which is
itself an extension of DBSCAN (Ester et al., 1996). DBSCAN is non-hierarchical
and will produce a flat clustering with noise of a given data set such that
the points in each cluster are densely connected with a density greater than
some threshold. All points with a density lower than this threshold are
classified as noise and are not clustered. OPTICS not only connects all points
in this manner, but also keeps track of how the points are connected together
– i.e. in which order and with what measure of local density. This adjustment
allows OPTICS to build a reachability plot which contains information of the
clustering structure.

Both the OPTICS and DBSCAN algorithms require $2$ parameters; $\epsilon$, a
search radius, and $N_{\mathrm{min}}$ (often denoted as MinPts), the minimum
number of points a cluster can have. The algorithms also need a distance
metric to be defined over the feature space of the input data. Conditioned
upon this metric, OPTICS uses the concepts of core distance and reachability
distance to produce its output – the reachability plot. The core distance of
any given point is the distance between that point and it’s
${N_{\mathrm{min}}}^{\mathrm{th}}$ nearest neighbour. It then follows that a
core-point is any point whose core distance is less than or equal to
$\epsilon$. The reachability distance of any point, $q$, with respect to
another, $o$, is the maximum of the core distance of $o$ and the distance
between $q$ and $o$.

Initially, OPTICS computes the core distances of each point and sets the
reachability distances of all points to infinity. The algorithm then
iteratively orders each point in the data set. In each iteration of this
ordering process, the next-to-be-ordered point is chosen as the point with the
smallest reachability distance from the set of all unordered points – for this
reason the first point is chosen at random. The chosen point is then removed
from the list of unordered points and appended to the ordered list. If this
point is a core-point, then all remaining unordered points within a radius of
$\epsilon$ of it are found. The reachability distance of each of these
unordered neighbours is then set to be the minimum of their currently assigned
reachability distance and their reachability distance with respect to the
recently ordered core-point. Following this step, the ordering process
continues on to the next iteration – repeating this cycle until all points are
in the ordered list. It is in this manner that OPTICS creates an ordering of
the data points, whilst concurrently seeking out regions of minimal
reachability distance, i.e. highest density. After the ordering process, the
reachability plot can then be constructed by plotting the final reachability
distances as a function of the corresponding ordered indices, for all points.
Within this plot, mutually clustered points appear as valleys since these
points are both denser than their surrounds (smaller distances between points)
and local to each other (ordered consecutively). More details on the OPTICS
algorithm can be found in Ankerst et al. (1999) and Oliver et al. (2020).

Halo-OPTICS takes the OPTICS algorithm and not only formalises the way in
which the OPTICS parameters are chosen but also provides a robust cluster
extraction technique that operates on the reachability plot in order to
extract a tree-structured hierarchy of suitable astrophysical clusters333The
tree-structure of this hierarchy allows for clusters to be referred to by the
typical terminology, e.g. root-level cluster (largest cluster in the tree),
parent cluster (superset of the child), child cluster (subset of the parent),
and leaf cluster (has no child clusters of its own).. In Halo-OPTICS, the
OPTICS parameter $\epsilon$ is effectively converted to the more physical
parameter $\Delta$, thereby classifying the root-level galaxy haloes using the
factor by which they are denser than the critical (or mean) density of the
Universe. Halo-OPTICS does this by using an approximate mapping between the
FOF linking length $l_{x}$, $N_{\mathrm{min}}$, and $\epsilon$. The technique
of extracting clusters is based on the designs of Sander et al. (2003), Zhang
et al. (2013), and McConnachie et al. (2018). The cluster extraction technique
of Halo-OPTICS creates the hierarchy of clusters by finding the valleys within
valleys that satisfy a range of conditions. These conditions assert that a
cluster must; contain at least $N_{\mathrm{min}}$ points; have a median
density that is a factor of $\rho_{\mathrm{threshold}}$ denser than the
cluster’s surrounds; not be a single child cluster of it’s parent cluster; not
share more than $f_{\mathrm{reject}}$ of it’s points with it’s parent; and it
rejects local outliers that have local-outlier-factors (defined in Eq. 8 and
Breunig et al. (1999)) greater than or equal to $S_{\mathrm{outlier}}$ from
each cluster. The Halo-OPTICS algorithm therefore transforms one OPTICS
parameter and adds $3$ unitless parameters for which near optimal values are
$\rho_{\mathrm{threshold}}=2$, $f_{\mathrm{reject}}=0.9$, and
$S_{\mathrm{outlier}}=2$. More details on the Halo-OPTICS algorithm can be
found in Oliver et al. (2020).

Input the data and parameters.Find all field haloes using 3D FOF or treatinput
data as if it were a field halo.For each field halo, search for
substructurewithin it using the process in Fig. 2.Output the hierarchy of
clusters. Figure 1: An activity of the outer methods of the CluSTAR-ND
process. Before finding the substructure CluSTAR-ND first decides on the root-
level clusters - where the option is given to produce 3D FOF field haloes or
to treat the input data as if it were a field halo.

##  3 CluSTAR-ND: A Hierarchical Galaxy/(Sub)Halo Finder

###  3.1 Concept and Motivation

Originally, our intention was to create an optimised version of the Halo-
OPTICS algorithm (refer to Sec. 2 for details) that implemented a locally
adaptive metric and could be applied to a data set with an arbitrary number of
features with minimal additional input from the user. Ultimately, since Halo-
OPTICS iteratively performs a radial search about each of the $n$
$d$-dimensional points in the data set, the cost of having it compute the
distances via a locally adaptive metric is too great. For $m$ points within
some radius of the query point, this algorithm would effectively need to have
computed the inverse and determinant of $m$-$d\times d$ covariance matrices
which makes the time complexity of the radial search and distance calculations
$\mathcal{O}(m\mathrm{log}(n))$ and $\mathcal{O}(md^{3})$ respectively (each
with different constant factors). In practice this increases the run-times
dramatically such that even applying this to small galactic haloes becomes
challenging.

We have designed CluSTAR-ND with the intention of mimicking Halo-OPTICS along
with the additional benefits of having; faster run-times, a variety of
adaptive metric settings, and the capability of performing on a data set with
any number of features. In addition to the input data, CluSTAR-ND has $7$
parameters that it uses to construct clusters:

  1. 1.

$l_{x}$ ($\in\mathbb{R}_{>0}$), the spatial linking length that is used to
find root-level haloes as discussed in Sec. 3.2. May be given by the user
although is set to $\infty$ by default – which makes CluSTAR-ND treat the
input data as a root-level halo.

  2. 2.

$adaptive$ ($\in\\{0,1,2\\}$), a flag that defines the behaviour of the metric
used to calculate the Mahalanobis distances between points as discussed in
Sec. 3.3. May be given by the user although is set to $1$ by default – which
provides consistently robust results as shown in Sec. 7.

  3. 3.

$k_{\mathrm{den}}$ ($\in\mathbb{N}_{\geq 7}$), the number of nearest
neighbours that are used to calculate the measure of local density for each
point as discussed in Sec. 3.4. May be given by the user although is set to
$20$ by default – which provides consistently robust results as shown in Secs.
6 and 7.

  4. 4.

$k_{\mathrm{link}}$ ($\in\mathbb{N}_{\geq 7\ \land\ \leq k_{\mathrm{den}}}$),
the number of nearest neighbours that are used to densely connect the points
in the data as discussed in Sec. 3.5. May be given by the user although is
automatically calculated to be optimal by default as in Sec. 6.1 – these
optimal values are based on the input data and the value of $k_{\mathrm{den}}$
and ensure maximal cluster completeness and inter-neighbourhood connectivity.

  5. 5.

$\rho_{\mathrm{threshold}}$ ($\in\mathbb{R}_{\geq 1}$), the factor by which
the median density of a cluster must be denser than that cluster’s surrounds
as discussed in Sec. 3.6. May be given by the user although is automatically
calculated to be optimal by default as in Sec. 6.2 – these optimal values are
based on the input data and the value of $k_{\mathrm{den}}$ and ensure that
the leaf level of the returned clusters are satellite-like overdensities.

  6. 6.

$f_{\mathrm{reject}}$ ($\in\mathbb{R}_{\geq 0\ \land\ \leq 1}$), the maximum
fraction of points that can be shared by parent-child clusters in the
hierarchy as discussed in Sec. 3.6. May be given by the user although is set
to $0.9$ by default – which provides a simple to interpret hierarchy that
agrees with analyses in Oliver et al. (2020) as discussed in Sec. 6.4.

  7. 7.

$S_{\mathrm{outlier}}$ ($\in\mathbb{R}_{\geq 1}$), the local-outlier-factor
that is used to reject outlier points from a candidate cluster as discussed in
Sec. 3.6. May be given by the user although is set to $2.5$ by default – which
provides a moderate level of outlier removal whilst maintaining good
clustering results as is shown in Sec. 6.3.

Since optimised values for $4$ of these parameters are found, the user only
needs to consider choosing values for $l_{x}$, $adaptive$ ,and
$k_{\mathrm{den}}$. The details of the algorithm CluSTAR-ND and how to choose
these three parameters is outlined in the subsections of Secs. 3 and 6.

###  3.2 Defining Root-Level Clusters

One major distinction between the group-finding component of the CluSTAR-ND
algorithm and that of the Halo-OPTICS algorithm is that the implementation of
CluSTAR-ND does away with having to perform a radial search about each point.
This is the largest factor in why CluSTAR-ND is so much faster than Halo-
OPTICS and is also why CluSTAR-ND does not have an equivalent parameter to
that of the OPTICS parameter $\epsilon$. In OPTICS this parameter not only
aids in the ordering of points but also defines the maximum reachability
distance that a clustered point can have – thereby prescribing an approximate
minimum density that a cluster can have. The functionality of $\epsilon$ that
leads to the detection and extraction of clusters can be approximated without
its presence (refer to Sec. 3.5 for details on this), however the effect of
defining a cluster’s minimum density requires another separate step.

In Halo-OPTICS, the choice of $\epsilon$ is redirected to the choice of
$\Delta$ – the overdensity factor. This is the factor by which a field halo
(root-level cluster) is denser than the critical density of the Universe,
$\rho_{\mathrm{crit}}$. This is accomplished via a mapping between the FOF
linking length, $l_{x}$, and $\epsilon$ and is used to approximate the FOF
field haloes. In CluSTAR-ND’s implementation, we offer the functionality of
directly computing the 3DFOF field haloes from the positions444If this
functionality is used, then the first $3$ features of the input data must be
the Cartesian spatial coordinates of the points. of the data set as a whole.
This splits the input data into some number of field haloes and opens the
opportunity for trivial parallelisation over each of these haloes for the
remaining steps.

These 3DFOF field haloes will be found by CluSTAR-ND whenever
$l_{x}\neq\infty$. A common choice for $l_{x}$ when finding field haloes from
cosmological simulations is $l_{x}=0.2L_{\mathrm{box}}/N$ – where
$L_{\mathrm{box}}$ is the side length of the simulation box and $N$ is the
number of particles within it – which corresponds to field halo overdensities
of $\gtrsim 100\bar{\rho}$ (Elahi et al., 2019). Alternatively when
$l_{x}=\infty$ (default), CluSTAR-ND will neglect finding field haloes in this
way which results in the root-level cluster becoming the entire data set. The
choice of whether to find field haloes using the 3DFOF approach or not must be
made under consideration of the input data and the intended output. This
choice is outlined in the first few nodes of Fig. 1.

Other hierarchical astrophysical clustering algorithms such as VELOCIraptor
(Elahi et al., 2019) also perform this step in order to define the root-level
clusters before then finding the substructure within them. However due to the
possibility of not having to perform this step in CluSTAR-ND, it is also
possible for the user to apply some other algorithm to the data – or to not do
so at all – in order to perform this step prior to using CluSTAR-ND. By
default $l_{x}$ is set to $\infty$ which forces CluSTAR-ND to treat the input
data as a single root-level halo from which the substructure therein is
subsequently found.

Once the root-level cluster(s) have been found, we add them to a list of
unsearched clusters in order to find the substructure within them. Each of the
unsearched clusters are independent of each other, meaning that we can search
for substructure in parallel. This remains true even when we set $adaptive=2$,
which ensures that CluSTAR-ND iteratively searches within the top level of
substructure until there are no more clusters found – as shown in Figs. 1, 2,
and 3. In this way, setting $adaptive=2$ creates a locally adaptive metric
since the data is transformed before each instance of finding substructure.
The following subsection provides details of how these transformations are
performed.

Input parent cluster and parameters.Take a PCA transformation of the data
andthen rescale each component to unit variance.Find the $k_{\mathrm{den}}$
nearest neighboursfor each point from the transformed data.Find the density of
each point using Eqs. 2 & 3.Produce a clustering, $\mathcal{C}$, using the
process in Fig. 3. if $\mathcal{C}\neq\varnothing\ \land\ adaptive=True$ Keep
only first levelof substructure in $\mathcal{C}$For each $c\in\mathcal{C}$,
search for substructurewithin it using the process in Fig. 2.Output the
hierarchy of substructure. Figure 2: An activity chart of the methods
concerned with the set up for finding substructure. A transformation of the
data is taken, the nearest neighbour lists are found, and then from them the
density of each point is calculated. This generalises the approach to
clustering the substructure. If such a clustering is found and $adaptive$
parameter is set to $2$, then the process in this activity chart is invoked
again on first level of substructure. This occurs iteratively until there are
no more significant clusters found, at which point the process produces a
cascade of returns of the adaptive hierarchy of substructure. If $adaptive$ is
set to $0$ or $1$, the process returns the hierarchy of substructure that has
been found using a parent cluster (globally) defined metric.

###  3.3 Data Transformation

When searching for substructure within a previously unsearched cluster, we
must implement a strategy that scales the various dimensions so that the
clusters found are not overly dependent on a subset of the features. Such a
strategy needs to be robust against the unit choice and coordinate
orientation. So that we produce the desired effects, we choose to use a
Principle Component Analysis (PCA) to first transform the data within the
cluster. We then scale each PCA component to unit variance. If we calculate
distances via the Euclidean distance metric on the transformed space, we are
then effectively using a Mahalanobis distance metric (Mahalanobis, 1936) such
that

$s^{2}(x_{i},x_{j})=(x_{i}-x_{j})^{T}\Sigma^{-1}(x_{i}-x_{j}),$ (1)

where $\Sigma=E[({\bf X}-\mu_{\bf X})({\bf X}-\mu_{\bf X})^{T}]$ is the
covariance matrix of the cluster before the transformation. This guarantees
that substructure found within a cluster will be dense with regards to the
cluster. When $adaptive=2$, CluSTAR-ND approximates a locally adaptive
Mahalanobis metric since a new covariance matrix is redefined at every level
of the hierarchy555It should be noted that this is not equivalent to the
locally adaptive Mahalanobis metric of EnLink (Sharma & Johnston, 2009), which
is constructed via a entropy-based binary space partitioning algorithm and
defines a unique covariance matrix for each point a priori to the clustering
of the data. While this method is very effective for defining a locally
adaptive metric based on the local distribution of points, it has a
substantially larger run-time than the top-down method we implement in
CluSTAR-ND.. Hence the choice of whether to set $adaptive$ to be $0$, $1$ or
$2$ must be made in consideration of the trade-off between run-time
constraints and the additional clustering power that an adaptive metric gives
over that of a global metric.

###  3.4 Density Estimation

Following the transformation of the data, we now seek to find a measure of the
local density surrounding each point. We do this by first finding the
$k_{\mathrm{den}}$ nearest neighbours for each point within the transformed
space. The density for each point is then approximated using the neighbour
list and a multivariate kernel within a balloon estimator such that

$\rho_{i}\propto\frac{1}{h_{i}^{d}}\sum_{j=1}^{k_{\mathrm{den}}}K\left(\frac{s(x_{i},x_{j})}{h_{i}}\right).$
(2)

Here $s(x_{i},x_{j})$ is defined in Eq. 1, $d$ is the dimensionality of the
feature space, $K$ is a kernel function, and $h_{i}$ is a smoothing length
corresponding to point $p_{i}$ – which we choose to be the distance from
$p_{i}$ to it’s ${k_{\mathrm{den}}}^{\mathrm{th}}$-most nearest neighbour.
Choosing the smoothing length in this way is not atypical (e.g. Sain, 2002),
however this also has the added bonus of being equal to the core distance from
OPTICS. Together with an appropriate kernel function, $K$, this choice allows
CluSTAR-ND to compute a determinable and smoothed analogue of the rigid
density estimator that is the reachability distance. For this we implement an
Epanechnikov kernel (Epanechnikov, 1969), which is defined as

$K(u)\propto(1-u^{2}),$ (3)

and is $0$ for $u>1$. Typically there is a normalisation constant that is
defined such that the integral of the kernel over the space is $1$. However,
as we clarify in the following subsection, the process of determining
significant substructure with CluSTAR-ND only relies on density by processing
points in order of descending density and by comparing the relative densities
of points. For this reason, we do not need to compute the constant factors of
Eqs. 2 and 3.

###  3.5 Densely Connecting Points via Neighbourly Aggregation

Once the data has been transformed, the nearest neighbour lists found, and the
local densities computed for any given unsearched cluster, the algorithm now
seeks to find the substructure of that cluster by aggregating the points that
are densely connected. The key to creating a similar output to that of Halo-
OPTICS without performing any radial search is to densely connect points via
their nearest neighbour lists. The OPTICS process allows for the algorithm to
gain knowledge about the locality of points nearby to (and from the
perspective of) the points that have already been appended to the ordered
list. It uses this knowledge to constantly be seeking out regions of higher
denser since the next-to-be-ordered point is always the point with the
smallest reachability distance i.e. largest density.

CluSTAR-ND mimics this process by further restricting the nearest neighbour
lists of each point to the nearest $k_{\mathrm{link}}$ points666The reason for
doing this is that the number of neighbours needed to densely connect points
is not necessarily equal to number of neighbours that should be used to
calculate the local density of a point. In fact, by restricting the nearest
neighbour lists, CluSTAR-ND is able to return a greater resolution on the
clustering structure – as is described in Sec. 3.5., finding all points that
are local density maxima with respect to their surrounds, and then using these
points as seeds from which to densely connect other points to. Seeding the
list of dense connections with each of the local maxima ensures that the
aggregation process within CluSTAR-ND does not need to seek out regions with a
higher density than the points that have already been connected together.
Instead it simply connects each point (that isn’t a local density maximum) in
order of decreasing density to a set of already densely connected points –
this process is similar to those of the core-search in VELOCIraptor (Elahi et
al., 2019) as well as the group-finders of SUBFIND (Springel et al., 2001) and
EnLink (Sharma & Johnston, 2009).

Input the points ($P$), $k_{\mathrm{den}}$ nearestneighbours, densities, and
parameters.Construct the $k_{\mathrm{link}}\ (\leq k_{\mathrm{den}})$ nearest
neighboursfor each point from their $k_{\mathrm{den}}$ nearest neighbours.Find
the set of local density maxima, $M$, and seed the listsof densely connected
points with them, $D=\\{\\{p\\}\mid\forall p\in M\\}$.Initialise the sets
unprocessed and clustered points,$U=P\setminus M$ and
$\mathcal{C}=\varnothing$ respectively.Retrieve the next point, $p_{i}$,
where$i=\operatorname*{arg\,max}_{i}\\{\rho_{i}\in\rho\mid p_{i}\in U\\}$.
Remove $p_{i}$ from $U$.Find the set of neighbourly connections,$T=\\{d\in
D\mid\exists p_{j}\in N_{i}^{\prime}\setminus U,\ p_{j}\in d\\}$. if $|T|>1$
Append $p_{i}$ to $d\in D$,where $d\in T$.Remove all $d\in T$from $D$.Append
the new set$\bigcup T\cup\\{p_{i}\\}$ to $D$.Filter $T$ usingEqs. 3.6 and
5.Append the elementsof $T$ to $\mathcal{C}$.While $U\neq\varnothing$Filter
$D$ usingEqs. 3.6 and 5. if $|D|>1$ Append the elementsof $D$ to
$\mathcal{C}$.Filter $\mathcal{C}$ using Eq. 3.6.Remove any $p$ from any $c$
in $\mathcal{C}$ not satisfying Eq. 9.Output $\mathcal{C}$, the subclusters of
$P$. Figure 3: The activity chart for the aggregation and rejection process
that constructs a hierarchy of subclusters similar to those from Halo-OPTICS.
The nearest neighbour lists are reduced and from them the local density maxima
are found. The lists of densely connected points are then seeded with these
local maxima. In order of decreasing local density, the points are either
appended to an existing list of densely connected points or used to merge
multiple into a new one. In the event of the latter and subject to the
conditions of Eqs. 3.6 and 5, these lists are considered to be potential
clusters. Following this inner loop, if not all points were densely connected,
the remaining connected lists are subject to these same conditions in order to
become potential clusters. The hierarchy of clusters is then cleaned using Eq.
3.6, and each remaining cluster is also cleaned subject to Eq. 9.

The key condition that CluSTAR-ND exploits here is that even though the OPTICS
algorithm gains knowledge of all unordered points within a radius of
$\epsilon$, the reachability distance of those points will typically be large
if they are not directly connectable via the $N_{\mathrm{min}}$ nearest
neighbours of either itself or the already ordered points. Hence at any given
time during OPTICS ordering process, the next-to-be-ordered point is almost
always a neighbour of a point (or vice versa) that has already been ordered.
The only time that it isn’t as such, is if the next-to-be-ordered point cannot
be directly connected through the $N_{\mathrm{min}}$-sized neighbourhoods of
those points that have already been ordered – or of such a neighbourhood of
it’s own. This does not mean that such a point can never be connected through
$N_{\mathrm{min}}$-sized neighbourhoods, however doing so will require
traversing points that are less dense than itself – densely connecting points
concurrently. Such a scenario where the point with the smallest reachability
distance cannot be densely connected through $N_{\mathrm{min}}$-sized
neighbourhoods is possible in OPTICS, although whether two points belonging to
entirely separate sets of densely connected points are at all associated is
somewhat uncertain in an astrophysical context. This distinction between
CluSTAR-ND and Halo-OPTICS is investigated further in Sec. 5.

The process of aggregating points is iteratively performed; by retrieving the
unprocessed point with the highest density ($p_{i}$), finding the set(s) of
densely connected points ($T$) that have members within the
$k_{\mathrm{link}}$ nearest neighbours of $p_{i}$, then either appending
$p_{i}$ to an existing list of densely connected points or constructing
clusters from $T$ and merging the lists of $T$ to create a new list of densely
connected points (with $p_{i}$ included). The condition for whether clusters
are to be constructed and a new list of densely connected points created, is
if $|T|>1$ i.e. if the point $p_{i}$ can be densely connected to multiple sets
of densely connected points. In this scenario, the previous sets of $T$ are
only retained as potential clusters if they satisfy Eq. 3.6 – a relative
density condition dependent upon the trainable parameter
$\rho_{\mathrm{threshold}}$. More details on this condition and others that
are responsible for cluster rejection are given in Sec. 3.6.

The circumstances under which $|T|>1$ for any given point depend strongly upon
the value of $k_{\mathrm{link}}$ since this determines how many neighbours are
considered for connecting the points. Decreasing $k_{\mathrm{link}}$ ensures
that a greater number of points can be densely connected to the potential
clusters before eventually finding a connecting point whose neighbourhood
satisfies $|T|>1$. In effect, this gives a greater resolution to the clusters
and, since $k_{\mathrm{link}}$ can be much less than $k_{\mathrm{den}}$, is
typically greater even than those from Halo-OPTICS. Of course, there is a
reasonable lower limit to the value of $k_{\mathrm{link}}$ that can be found
by ensuring that all points be densely connected under a uniform distribution
of points. This lower limit is investigated further in Sec. 6.1. Finally, if
all points within the previously unsearched parent cluster do not become
densely connected through this aggregation process and if more than $1$ of the
densely connected subsets have at least $k_{\mathrm{link}}$ points, then we
append those that do to the list of substructures as well.

###  3.6 Cluster and Cluster Member Rejection

The cluster extraction algorithm of Halo-OPTICS is implemented in CluSTAR-ND,
however it is not entirely its own separate process that simply takes place
after the aggregation/ordering of points as it is in Halo-OPTICS. To detect
and extract clusters from the reachability plot, Halo-OPTICS performs a set of
routines – each of which are described in Sec. 3.2 of Oliver et al. (2020).
First, the algorithm determines all local maxima of the reachability plot from
OPTICS and then using this, builds the hierarchy of clusters by taking
contiguous subsets of the ordered list that are situated on either side of
each local maxima – ensuring that each subset only contains points with
reachability distances less than that at the corresponding local maxima.
CluSTAR-ND performs these steps by ensuring that points are aggregated in
order of decreasing local density and then by considering lists of densely
connected points that are connected to one another via the neighbourhood of a
less dense point as potential clusters.

Halo-OPTICS then rejects all potential clusters that contain less than
$N_{\mathrm{min}}$ points or that have median local densities less than
$\rho_{\mathrm{threshold}}$ times the surrounding density of those potential
clusters. Following this, the algorithm rejects all single leaf potential
clusters. These criteria are fulfilled by CluSTAR-ND in a single step during
the aggregation process when multiple lists of densely connected points are
considered to be potential clusters. As such, CluSTAR-ND conducts the
following:

$\displaystyle\mathrm{reject\ any}\ d\ \mathrm{from}\ T\ \mathrm{if}$
$\displaystyle\ \mathrm{it\ has\ no\ children\ and\ either;}$
$\displaystyle|d|<k_{\mathrm{link}}\ \mathrm{or}$
$\displaystyle\mathrm{median}\\{\rho_{j}\mid p_{j}\in
d\\}/\rho_{i}<\rho_{\mathrm{threshold}}.$ (4)

Here, $\rho_{i}$ is the local density of $p_{i}$ – the point responsible for
densely connecting the otherwise disconnected sets in $T$. Following the
rejection of all $d\in T$ that do not meet these criteria, the set $T$ is then
subject to the condition that

$\mathrm{if}\ |T|=1,\ \mathrm{set}\ T=\varnothing.$ (5)

The conditions in Eqs. 3.6 and 5 are related to steps $3$ and $4$ from the
cluster extraction process of Halo-OPTICS. Next, for all parent-child cluster
pairs sharing at least $f_{\mathrm{reject}}$ of the parent’s points; CluSTAR-
ND rejects the child if it has child clusters of its own, otherwise CluSTAR-ND
rejects the parent. To re-iterate, given a set of substructure $S$ and any
parent-child pair $s_{\mathrm{parent}}$ and $s_{\mathrm{child}}$,

$\displaystyle\mathrm{if}\ |s_{\mathrm{child}}|/|s_{\mathrm{parent}}|$
$\displaystyle\geq f_{\mathrm{reject}}\ \mathrm{then},$
$\displaystyle\mathrm{reject}\ s_{\mathrm{child}}\ \mathrm{if;}\ \exists s\in
S\ \mathrm{such\ that}\ s\subset s_{\mathrm{child}},$
$\displaystyle\mathrm{else;\ reject}\ s_{\mathrm{parent}}.$ (6)

This is equivalent to step $5$ of the extraction process in Halo-OPTICS. Steps
$6$ and $7$ of the Halo-OPTICS extraction process are concerned with removing
outliers based on how their local density compares with the rest of their
neighbours. While this approach works in the context of using a reachability
distance as the measure of inverse density, it can be improved upon with a
more robust measure of density such as that that CluSTAR-ND computes. The
measure used by Halo-OPTICS is the local-outlier-factor which relies upon the
local-reachability-density (Breunig et al., 1999), both of which are defined
in Eqs $3$ and $4$ of Oliver et al. (2020). In effect, the local-outlier-
factor compares how reachable a point is from its neighbours to how reachable
its neighbours from their own. We provide CluSTAR-ND with a comparable measure
with which to define outliers that incorporates the estimate of local density
it uses from Eq. 2.

We define the kernel density estimate analogue of the local-reachability-
density of a point $p_{i}$ with respect to transformed space as

$\mathrm{lrd}(p_{i})\propto\frac{1}{\sum\limits_{p_{j}\in
N_{k_{\mathrm{den}}}}\mathrm{min}\left(\rho_{i},\ \rho_{j}\right)^{-1/d}},$
(7)

where $N_{k_{\mathrm{den}}}$ is the $k_{\mathrm{den}}$ nearest neighbours of
$p_{i}$ within the previously unsearched cluster. CluSTAR-ND then uses this to
similarly compute an analogue of the local-outlier-factor,

$\mathrm{lof}(p_{i})=\frac{\sum\limits_{p_{j}\in
N_{k_{\mathrm{den}}}}\frac{\mathrm{lrd}(p_{j})}{\mathrm{lrd}(p_{i})}}{k_{\mathrm{den}}}.$
(8)

In Halo-OPTICS, a point is considered an outlier from a cluster if the local-
outlier-factor (calculated with respect to that cluster) is larger than
$S_{\mathrm{outlier}}$. In CluSTAR-ND we adjust this rule by first defining a
cut-off density,

$\rho_{\mathrm{cut}}=\mathrm{min}\\{\rho_{i}\mid\mathrm{lof}(p_{i})<S_{\mathrm{outlier}}\\},$
(9)

which we then use to reject all points from a cluster that have a density less
than this. This ensures that only local-outliers at the outskirts of a cluster
are rejected and not those that are embedded deeper within it – a more
appropriate regime for outlier detection and rejection within a hierarchical
set of clusters.

##  4 Synthetic Data

Within the remainder of this paper we compare (Sec. 5), train (Sec. 6), and
scrutinise (Sec. 7) our algorithm against synthetic survey data of Milky Way
(MW) type galaxies produced by Galaxia (Sharma et al., 2011). Given certain
survey restrictions, such as one or more colour-magnitude bounds, a survey
size, and geometry, Galaxia is able to return a synthetic catalogue of stars
in accordance with a given model of the MW. The code is generalised enough so
that it can also accept star formation rates, age–metallicity relations,
age–velocity-dispersion relations, and analytic density distribution functions
from which it can also use to produce stellar catalogues.

Among the synthetic data that can be produced with Galaxia are resamplings of
the $11$ $\Lambda$CDM stellar haloes from Bullock & Johnston (2005) and the
complementary $6$ artificial stellar haloes from Johnston et al. (2008). The
original $\Lambda$CDM haloes are simulated using a hybrid semi-analytical plus
hydrodynamic $N$-body approach that replicates a density profile and satellite
distribution that is similar to the MW. The simulation model assumes a
$\Lambda$CDM cosmology with parameters $\Omega_{m}=0.3$,
$\Omega_{\Lambda}=0.7$, $\Omega_{b}h^{2}=0.024$, $h=0.7$, and
$\sigma_{8}=0.9$. The simulations generate, track, and evolve a number of
individual satellites that are first modelled as $N$-body dark matter systems
within a parent galaxy whose disk, bulge, and halo are represented by time
dependent semi-analytic functions. Semi-analytical prescriptions are then used
to assign star formation histories and leaky accreting boxes to each. A
chemical enrichment model is also used to calculate the metallicities as a
function of age for the stellar populations (Robertson et al., 2005; Font et
al., 2006). Ultimately, the dark matter distributions of the satellites follow
NFW profiles (Navarro et al., 1996) while the stellar distributions follow
King profiles (King, 1962) – the latter of which is constructed in order to
reproduce an agreement with the structural properties of the Local Group’s
dwarf galaxies.

Each satellite created in this way has three main model parameters; the time
since accretion ($t_{\mathrm{acc}}$), the luminosity ($L_{\mathrm{sat}}$), and
the orbital circularity ($\epsilon=J/J_{\mathrm{circ}}$). The distribution of
these parameters specify the accretion history of a halo and as such a further
$6$ artificial haloes were created in Johnston et al. (2008) for the purpose
of studying the effects of different accretion histories on the properties of
haloes. These $6$ artificial haloes have accretion events that are
predominantly; radial ($\epsilon<0.2$); circular ($\epsilon>0.7$); old
($t_{\mathrm{acc}}>11$ Gyr); young ($t_{\mathrm{acc}}<8$ Gyr); high luminosity
($L_{\mathrm{sat}}>10^{7}L_{\mathrm{\odot}}$); and low luminosity
($L_{\mathrm{sat}}<10^{7}L_{\mathrm{\odot}}$). Importantly, the output from
Galaxia also bestows each star with a label that corresponds to which
satellite group it belonged to at the time that that satellite was created
within the simulation – this allows us to test how well CluSTAR-ND performs in
Sec. 6. In addition to this, Galaxia maintains a list of satellite properties
that includes information on whether the satellite is self-bound or not. More
details on these simulations can be found in Sec. 3.4 of Sharma et al. (2011)
and references therein.

We use Galaxia to produce a random sample from each synthetic galaxy in both
sets of the stellar haloes. From these we use various combinations of the
spatial (denoted by ${\bf x}=(x,y,z)$ in the later sections), kinematic
(denoted by ${\bf v}=(v_{x},v_{y},v_{z})$ in the later sections), and chemical
(denoted by ${\bf m}=($[Fe/H],[$\alpha$/Fe]$)$ in the later sections)
information in order to compare, optimise, and apply our algorithm. For each
galaxy, the points are contained inside the $282$ kpc virial radius of the
corresponding host dark matter halo – these have an overdensity factor of
$\Delta\approx 337$ times the mean dark matter density of the universe. As
such we treat these galaxies as field haloes and hence do not use CluSTAR-ND
to find field haloes before finding substructure. Since these field haloes and
those that are found via 3D FOF are well defined, the analysis in the
following sections is an assessment of the substructure finding component of
CluSTAR-ND.

Figure 4: A clustering comparison of the outputs of Halo-OPTICS and CluSTAR-ND
after having been applied to the 3D spatial dimensions of the synthetic
galaxies from Galaxia. Each panel belongs to one of three pairs and each pair
corresponds to a different comparison wherein CluSTAR-ND has been used with a
different metric adaptivity setting – which is annotated on the relevant
panels. Plotted within each panel are the recovery and purity fractions of
each best-fitting cluster to either; the Halo-OPTICS predicted clusters (left
panel in each pair) or the CluSTAR-ND predicted clusters (right panel in each
pair). Best-fitting clusters are those that produce the maximum Jaccard index
in Eq. 5.1 (i.e. $T^{*}=\operatorname*{arg\,max}_{T\in S_{2}}J(C,T)$), the
recovery and purity are then found according to Eq. 5.1. To clarify, the label
at the top of each panel indicates which clustering is used as $S_{1}$ within
these equations. Every cluster, except for the root-level clusters (L = $0$),
appear in these panels as a marker coloured by its level within the hierarchy
and shaped by whether or not its best-fitting cluster is a mutually best-
fitting – note that BFNM here stands for best-fitting not mutual. The axes are
extended beyond the possible recovery and purity values for easier
readability, this extended area is coloured grey. Within this area on each
panel and for both the recovery and purity distributions, there is a box and
whisker plot that denotes the $Q_{0}$ – $Q_{4}$ quartiles with mean values
also indicated with diamonds markers.

##  5 A Comparison with Halo-OPTICS

Our first assessment of CluSTAR-ND’s performance is done by comparing it to
Halo-OPTICS. We compare both the clustering and the run-times of these
algorithms when applied to the synthetic haloes outlined in Sec. 4. Since
Halo-OPTICS is designed for clustering on 3D spatial data, we only apply these
algorithms to the 3D positions of the points in these haloes.

We apply CluSTAR-ND in a way that treats each synthetic galaxy as a field halo
($l_{x}=\infty$) – due to each galaxy only containing points within its virial
radius. Similarly, for Halo-OPTICS we find the smallest value of $\epsilon$
that will guarantee that every point in these reduced galaxies is a part of at
least one core-point’s $N_{\mathrm{min}}$ nearest neighbours – thereby
ensuring that every point will be given a reachability distance less than or
equal to $\epsilon$. By applying these algorithms to ready-made galactic data
sets, we can ensure that we are only comparing the substructure found and not
the root-level clusters – which in both cases are equal to the entire data set
due to the above.

In order to assess how adequately CluSTAR-ND reproduces Halo-OPTICS, we set
$k_{\mathrm{den}}=20$ and similarly for Halo-OPTICS we set
$N_{\mathrm{min}}=20$. For CluSTAR-ND we set $k_{\mathrm{link}}=7$, as is
described in Sec. 6.1. For both codes we use the near-optimal parameters found
in Oliver et al. (2020) such that $\rho_{\mathrm{threshold}}=2$,
$f_{\mathrm{reject}}=0.9$. However since the outlier detection is different
between the codes, we set $S_{\mathrm{outlier}}=\infty$ which ensures that
there are no outliers removed from clusters. We review the optimal choice of
these parameters in Secs. 6.2 – 6.4.

###  5.1 Output Similarity

To appropriately compare the output of the two codes, we find the best-fitting
match from the Halo-OPTICS catalogue of clusters for each of the CluSTAR-ND
clusters (ignoring the root-level clusters which are fixed to be equal). As a
measure of which pair is best-fitting we use the maximum Jaccard index
(Jaccard, 1912), which – given two clusterings of the same data set ($S_{1}$
and $S_{2}$) – allows us to not only find which clusters from $S_{1}$ are
best-fitting to the clusters from $S_{2}$ but to also compare how well the
clusters in $S_{2}$ are matched by those in $S_{1}$. Hence for two such
clusterings the maximum Jaccard index is given by

$\displaystyle J_{\mathrm{max}}(C)=\mathrm{max}\\{J(C,T)\mid\forall T\in
S_{2}\\},\ C\in S_{1},\ \mathrm{where}$ $\displaystyle J(C,T)=\frac{|C\cap
T|}{|C\cup T|}.$ (10)

Once each cluster from $S_{1}$ has been assigned a best-fitting cluster from
$S_{2}$, we then compute the recovery and purity fractions of the pairs. We
define these such that

$\displaystyle R(C)=\frac{|C\cap T^{*}|}{|T^{*}|},\ \mathrm{and}$
$\displaystyle P(C)=\frac{|C\cap T^{*}|}{|C|}.$ (11)

where $T^{*}\in S_{2}$ is the best-fitting match to $C\in S_{1}$ – i.e.
$T^{*}=\operatorname*{arg\,max}_{T\in S_{2}}J(C,T)$. Neatly, this implies that

$\displaystyle\frac{1}{J_{\mathrm{max}}(C)}=\frac{1}{R(C)}+\frac{1}{P(C)}-1,\
\mathrm{and\ hence}$ $\displaystyle
J_{\mathrm{max}}(C)\leq\mathrm{min}\\{R(C),\ P(C)\\}.$ (12)

Fig. 4 depicts the recovery and purity fractions determined in this way for
all of the best-fitting cluster pairs found within the synthetic haloes from
Galaxia for various values of the $adaptive$ parameter. Here we see that the
clusterings from CluSTAR-ND are best-fitted by those from Halo-OPTICS in a way
that most commonly leads to clusters from CluSTAR-ND being encompassed by
those from Halo-OPTICS – indicated in each comparison by high purity and
varied recovery in the left panel of each pair and the opposite in the right
panel of each pair. This matches the trends seen in both Figures 6 and 10 of
Oliver et al. (2020) where Halo-OPTICS was typically shown over-encompassing
both mock cluster sets and predicted clusters found by VELOCIraptor (Elahi et
al., 2019). Regardless, this property of the best-fitting clusters from
CluSTAR-ND and Halo-OPTICS suggests that CluSTAR-ND may benefit from an
additional method to allocate points to the clusters that it already finds.
Such a task could be achieved via an expectation-maximisation-like technique
that allocates previously noisy points to already existing clusters from
CluSTAR-ND depending on how well they would assimilate within them.

A small number of clusters are also shown in Fig. 4 to have high purity and
very low recovery ($<0.2$). These clusters are those that have been found by
one algorithm but not the other – hence they are not mutually best-fitting –
and are consequences of: 1) CluSTAR-ND’s density estimate being more reliable
than that of Halo-OPTICS; 2) CluSTAR-ND being able to robustly detect clusters
with a smaller number of points; and 3) the algorithmic differences between
connecting points. We do also see a few high recovery and high purity clusters
whose best-fitting cluster is not mutually best-fitting returned by both
codes. These are clusters from one code that effectively sit between two
levels of the hierarchy in the other code – this can produce a good fit to a
cluster that is not necessarily a mutually best-fitting cluster. These
clusters arise for the same reasons as above.

Overall, it is clear that CluSTAR-ND is able to reproduce the clusters and
hierarchy of Halo-OPTICS with a similar clustering power. We see from Fig. 4
and the box and whisker plots therein that typically the clusters from
CluSTAR-ND (Halo-OPTICS) are mutually best-fitting to those from Halo-OPTICS
(CluSTAR-ND) with both high recovery and high purity with medians of $\sim
0.860$ ($\sim 0.990$) and $\sim 0.998$ ($\sim 0.926$) respectively.

For completeness, we also discuss the clustering power of these outputs from
Halo-OPTICS in comparison with the optimised outputs of CluSTAR-ND in Sec. 7.
In Sec. 7 we see that overall CluSTAR-ND is expected to provide an equal
clustering power to Halo-OPTICS, however the clustering power does tend to
differ slightly on a per galaxy basis. This difference is due to the competing
effects of CluSTAR-ND having a more robust density estimation, while the
ordering process of Halo-OPTICS can more easily gather like-satellite points
before joining child clusters into their shared parent cluster. We will now
define this measure of clustering power as it is used throughout both Secs. 6
and 7.

####  5.1.1 Measuring Clustering Power

In order to measure the clustering power of CluSTAR-ND we must construct a
scalar function that reduces the complexity of comparing the quality of fit
between an entire CluSTAR-ND clustering output when applied to a synthetic
galaxy and the ground truth labels of that synthetic galaxy. Since each galaxy
has been constructed through the accretion of satellites (refer to Sec. 4 for
details), the ground truth labels form a flat clustering without noise, i.e.
each data point in the galaxy is a member of exactly one satellite. This is a
fundamentally distinct clustering type from the hierarchical clustering with
noise that CluSTAR-ND produces – so we only find the goodness of fit between
the leaf clusters produced by CluSTAR-ND and the satellite labels within each
galaxy.

In order to compare a leaf-clustering produced by CluSTAR-ND ($C_{g}$) to the
ground truth clustering ($T_{g}$) defined over the same galaxy ($g$), we
construct a mutual-information-based objective function similar to that built
by Vinh et al. (2009). To do this we must first append an additional cluster
to those in $C_{g}$ to create $C_{g}^{*}$ so that the sets $C_{g}^{*}$ and
$T_{g}$ are defined over the same set of data points. The additional cluster
is made of the remaining points in the data set that have not been clustered
into the leaf clusters returned by CluSTAR-ND. This means that $C_{g}^{*}$ is
an artificial construction of a flat clustering without noise using $C_{g}$
which we can then easily match to the flat clustering without noise that is
$T_{g}$.

The mutual-information between $C_{g}^{*}$ and $T_{g}$, $I(C_{g}^{*};T_{g})$,
is then the amount of information obtained about the true clusters by having
observed the predicted clusters (Shannon, 1948).

$\displaystyle I(X;Y)=H(X)-H(Y|X),\ \mathrm{where}$ $\displaystyle
H(X)\equiv-\sum_{x\in X}P(x)\ \mathrm{log}(P(x)),\ \mathrm{and}$
$\displaystyle H(Y|X)\equiv-\sum_{x\in X,\ y\in Y}P(x,y)\
\mathrm{log}\bigg{(}\frac{P(x,y)}{P(x)}\bigg{)}.$ (13)

This measure is non-negative and will have a different theoretical maximum
depending on the synthetic galaxy in question. For our purposes, we normalise
$I(C_{g}^{*};T_{g})$ between $0$ and $1$ such that these values represent the
absolute worst and best values that a clustering algorithm can be expected to
have respectively (see below). We then also take the average of this
normalised value over each of the synthetic galaxies so that we may compare
sets of CluSTAR-ND clusterings to sets of ground truth clusterings. Let
$\mathcal{S}$ represent some feature space combination on which the predicted
clusters are dependent, then the resultant objective function we use is
defined as

$\displaystyle
F(\mathcal{S})=\frac{1}{N_{\mathrm{galaxies}}}\sum_{g\in\mathrm{galaxies}}F_{g}(\mathcal{S}),\
\mathrm{where}$ $\displaystyle
F_{g}(\mathcal{S})=\frac{I(C_{g}^{*}(\mathcal{S});T_{g})-E[I(R_{g}^{*}(\mathcal{S});T_{g})]}{H(T_{g})-E[I(R_{g}^{*}(\mathcal{S});T_{g})]}.$
(14)

Here $H(T_{g})$ is the entropy of $T_{g}$ (this normalises the maximum of
$F_{g}(\mathcal{S})$ to $1$ since $I(T_{g};T_{g})=H(T_{g})$) and
$E[I(R_{g}^{*}(\mathcal{S});T_{g})]$ is the expected mutual-information that
arises when a clustering, $R_{g}^{*}(\mathcal{S})$, with equal cluster sizes
as $C_{g}^{*}(\mathcal{S})$ has been created via random assignment (this
normalises the minimum of $F_{g}(\mathcal{S})$ to $0$ since no clustering
algorithm can ever do worse than random assignment). Unlike Vinh et al. (2009)
(who relies on a hypergeometric model), we calculate
$E[I(R_{g}^{*}(\mathcal{S});T_{g})]$ empirically by taking the average of the
mutual-information that is found when points have been randomly assigned to a
mock clustering, $R_{g}^{*}(\mathcal{S})$, with the number of clusters and
number of points within each cluster being equal to those in
$C_{g}^{*}(\mathcal{S})$. This average is calculated using $10$ random
realisations of $R_{g}^{*}(\mathcal{S})$ which we find to be sufficient since
the variance of $I(R_{g}^{*}(\mathcal{S});T_{g})$ is small compared to its
expectation.

The normalisation adjustments we have made mean that $F_{g}(\mathcal{S})$ can
now be interpreted as the proportion of relevant information obtained about
the true clusters of a galaxy by having observed the predicted leaf clusters
produced from CluSTAR-ND when applied to that galaxy. Similarly,
$F(\mathcal{S})$ is the average of this proportion across the galaxies. We use
these mutual-information-based measures of clustering power throughout both
Secs. 6 and 7 as they provide the means to compare entire clusterings
($F_{g}(\mathcal{S})$) and sets of entire clusterings ($F(\mathcal{S})$).

While it is possible to construct an objective function with a similar purpose
from the easy-to-interpret measures of recovery/purity/Jaccard index, it is
non-trivial to reduce these to an appropriate scalar measure of clustering
power as these are designed to be used to compare one cluster to another
(rather than sets of clusters). Most critically, these measures do not account
for true-negative predictions and so optimising the parameters of CluSTAR-ND
(as is done throughout Sec. 6 with $F(\mathcal{S})$) by using an objective
function built from such measures will likely create an artificial incentive
for CluSTAR-ND to make fewer predictions on clusters. These measures also do
not appropriately reward CluSTAR-ND for having matched a true cluster with two
predicted clusters – clearly this is not as nice as a perfect match, but it is
better than a $J_{\mathrm{max}}(C)=0.5$ match for example. Nevertheless,
maximising $F(\mathcal{S})$ does also maximise the recovery, purity, and
Jaccard index that the predicted clusters will have – it just does so without
promoting fewer predictions. In fact, a value of $F_{g}(\mathcal{S})=f$ is
approximately equivalent to having $\sim f\times|T_{g}|$ true clusters matched
perfectly ($J_{\mathrm{max}}(C)=1$) by CluSTAR-ND, of course, this relation is
not one-to-one as $F_{g}(\mathcal{S})=f$ could also occur via a larger number
of predicted clusters with smaller $J_{\mathrm{max}}(C)$ values – although, by
adjusting $F_{g}(\mathcal{S})$ for random chance allocation we have limited
this trade-off to meaningful predictions.

###  5.2 Run-time Disparity

We now compare run-times of CluSTAR-ND and Halo-OPTICS using synthetic
galaxies from Galaxia in Fig. 5. All runs were performed using a single core
on an Intel i5 vPro processor and both codes are written using Python3 but do
make use of optimised numerical packages such as NumPy (Harris et al., 2020),
SciPy (Virtanen et al., 2020), Scikit-learn (Pedregosa et al., 2011), and
Numba (CluSTAR-ND only; Lam et al., 2015). We see in Fig. 5 that for these
galaxies the run-times of CluSTAR-ND are at least $3$ orders of magnitude
faster than those of Halo-OPTICS and this disparity grows for larger data
sets. CluSTAR-ND follows an $\mathcal{O}(n\mathrm{log}(n))$ time complexity
while Halo-OPTICS appears to follow a $\mathcal{O}(n^{2.17})$ time complexity.
This is not the true time complexity of Halo-OPTICS and is partially an
artefact of the changing clustering structuring within the galaxies – a
feature that affects both algorithms albeit affecting CluSTAR-ND to a much
lesser degree.

Given two data sets with an equal number of points, the run-time of Halo-
OPTICS will be larger for the data set that has a larger range of densities
within the virial radius i.e. is more cusp-like rather than core-like. This is
because there will be a larger fraction of points whose mutual distance is
less than or equal to $\epsilon$. Ordinarily, the best-case run-time of Halo-
OPTICS is $\mathcal{O}(n\mathrm{log}(n))$, however as the larger synthetic
galaxies we use in this paper are more cusp-like, we find that the time
complexity of Halo-OPTICS approaches $\mathcal{O}(n^{2})$ with increasing $n$.
So while it is likely that the run-time of Halo-OPTICS applied to each galaxy
individually is sub-quadratic, the overall trend of increasing cupsyness gives
rise to a super-quadratic time complexity of $\mathcal{O}(n^{2.17})$.

Figure 5: The run-times of both CluSTAR-ND and Halo-OPTICS when applied to the
3D spatial dimensions of the synthetic galaxies from Galaxia. The lines of
best fit correspond to an $O(n\mathrm{log}(n))$ time complexity for CluSTAR-ND
and an $O(n^{2.17})$ time complexity for Halo-OPTICS. We see that even when
applied to the smallest size galaxy ($\sim 2\times 10^{5}$ points) in the
suite CluSTAR-ND is $3$ orders of magnitude faster than Halo-OPTICS, and due
to the difference in time complexities, this speed increase balloons outwards
to nearly $5$ orders of magnitude when applied to the largest of these
galaxies ($\sim 2.5\times 10^{6}$ points). While some of this run-time
disparity can be attributed to changing clustering structure (see Sec. 5.2),
the modest run-time of CluSTAR-ND makes it ideally suited for application to
large data sets.

As CluSTAR-ND does not perform a radial search it does not suffer from this
same drawback and moreover, the clustering structure does not strongly affect
the run-time of the $k_{\mathrm{den}}$ nearest neighbour search and density
computation either. In reality the run-time of CluSTAR-ND in Fig. 2 is
dominated by the $k_{\mathrm{den}}$ nearest neighbour search – taking up a
constant $\sim 89$% of the run-time for this number of dimensions and this
value of $k_{\mathrm{den}}$. In this implementation, CluSTAR-ND uses SciPy’s
cKDTree (Bentley, 1975; Maneewongvatana & Mount, 1999; Virtanen et al., 2020)
which has an expected build and search run-time of
$\mathcal{O}(n\mathrm{log}(n))$. Such a run-time complexity is theoretically
the fastest that can be achieved for this problem and hence CluSTAR-ND’s
overall run-time only stands to improve by some constant factor. Faster run-
times can be achieved this way by running the $k_{\mathrm{den}}$ nearest
neighbour search in parallel, by using a faster implementation of the kdtree,
or by using an approximate nearest neighbour algorithm rather than an exact
one.

In the event that these run-time improvements are used and the dimensionality
of the data set is sufficiently low, it is possible that the run-time of
$k_{\mathrm{den}}$ nearest neighbour search becomes small enough that it is no
longer the most time consuming component of the algorithm. The next most time
consuming component of the CluSTAR-ND algorithm is the aggregation and
rejection process in Fig. 3 which takes up a constant $\sim 10$% of the total
run-time for each of the runs in Fig. 5. This run-time is predominantly
affected by the merging of potential clusters. In the best-case run-time
scenario, there will be no mergers and the time complexity for this part of
the algorithm would be $\mathcal{O}(n)$. While technically it is possible for
a data set to be structured in such a way that gives rise to no mergers during
the aggregation process, for astrophysical data, it would be highly atypical.
Most if not all points in the data set will be densely connected together
during this process. As the mergers are discovered they are combined in a
tree-like structure. By merging the densely connected lists in a compact
binary-tree structure the aggregation process will be performed in an
$\mathcal{O}(n\mathrm{log}(n))$ time. This is both the expected and worst case
time complexity. The time complexity of the aggregation process is not only
affected by the way in which the densely connected lists are merged but also
by how many times this occurs. The values $k_{\mathrm{den}}$ and
$k_{\mathrm{link}}$ both affect this number and will increase the run-time by
some constant factor the smaller that they are.

The remaining $\lesssim 1$% of CluSTAR-ND’s run-time can be attributed to the
transformation of the data and general data manipulation – the former of which
has an $\mathcal{O}(nd^{2}+d^{3})$ time complexity. Overall, the these
constituents make the run-time complexity of CluSTAR-ND
$\mathcal{O}(n\mathrm{log}(n))$. In Fig. 5 we do not show the run-time of
CluSTAR-ND when $adaptive=0$ since this is only slightly different ($\lesssim
1$% different as above) from when $adaptive=1$. We see that when $adaptive=2$
the run-time increases by some factor per data set. In practice, this factor
that is less than the number of levels found within the clustering hierarchy.
When applied to the galaxies we analyse in this paper, the depth of the
hierarchy is typically capped at $3$ or $4$ levels. CluSTAR-ND’s fast run-time
makes it ideal for application to large data sets such as the Gaia DR3
catalogue.

##  6 Parameter Influence and Optimisation

The CluSTAR-ND algorithm has $7$ parameters. Parameters $l_{x}$, $adaptive$,
and $k_{\mathrm{den}}$ may be chosen by the user and are responsible for the
spatial linking length that finds field haloes, the adaptivity of the distance
metric, and the number of neighbours used to find the density of the points
respectively. A typical choice for $l_{x}$ is $0.2$ times the inter-particle
spacing and then whether $adaptive$ is set to $0$, $1$, or $2$ is decided in
consideration of the run-time vs clustering power trade-off as well as whether
the user wishes to apply a PCA transform to the input data.

The parameter $k_{\mathrm{den}}$ has a more loosely constrained range – if it
is too small then the density fluctuations between neighbouring points can be
large and noisy; if it too large then significant differences in local density
may be smoothed out enough that some structures may not be detected by
CluSTAR-ND. Typically, similar parameters to $k_{\mathrm{den}}$ in other
clustering algorithms are chosen with values ranging from $20$ to $500$. It is
commonly accepted that $k_{\mathrm{den}}$ should increase when using data with
a larger a feature space – however this depends on the intended resolution of
the density estimate and, in a clustering context, the expected number of
points within the true clusters present in the data. A more thorough
investigation of the effects of the choice of $k_{\mathrm{den}}$ is carried
out in Secs. 6.2 – 6.4.

Although the remaining $4$ parameters in CluSTAR-ND may be chosen by the user,
these do have properties that allow for their optimal values to be determined.
It should be mentioned that these parameters do have some co-dependence with
regards to the state of the output, though this influence is small. Due to the
fact that none of these parameters act within the same step of the algorithm,
we are able to isolate their effects in order to find their optimal values.
Furthermore and since the labels within the synthetic data provide a flat
clustering without noise, we opt to assess the performance of CluSTAR-ND – and
hence the corresponding optimal values of some of the cluster extraction
parameters ($\rho_{\mathrm{threshold}}$ and $S_{\mathrm{outlier}}$) – by
determining how well the leaf clusters of CluSTAR-ND are able to match a
single label within the synthetic data sets.

While the affects of the cluster extraction parameters in CluSTAR-ND are
expected to be similar to those in Halo-OPTICS, the lists of points that they
act upon are generally slightly different. The near-optimal values of Halo-
OPTICS parameters were determined by maximising the recovery and purity
fractions of a set of hand-crafted 3D Gaussian distributions in both Cartesian
and spherical coordinate systems. While this technique can be used to mimic
returning the optimal set of astrophysical clusters, it of course is not
equivalent to doing so. Furthermore, this investigation did not analyse the
effects of higher dimensions and/or other subspaces on these parameters. We
now investigate the effects of these parameters and deduce their preferred
values.

###  6.1 Values for $\boldsymbol{k_{\mathrm{link}}}$

The parameter $k_{\mathrm{link}}$ is the number of neighbours used to
aggregate the points together in order of decreasing density so that the
substructure may be found. The consequence of decreasing $k_{\mathrm{link}}$
is that the clusters will be more complete up until their true boundary
density. Irrespective of this however, a practical lower limit exists. During
the aggregation process, multiple potential clusters will be connected
together into a new potential parent cluster whenever the current point being
processed has multiple of its $k_{\mathrm{link}}$ nearest neighbours in
different potential clusters. This has the consequence that at most
$k_{\mathrm{link}}-1$ potential clusters can be connected together during this
step (since the point being processed is one of its own neighbours and does
not yet belong to any potential cluster). The number of potential clusters
that can be joined depends on the specific arrangement of points – although
with a sufficiently large value of $k_{\mathrm{link}}$, the number of these
clusters that are observed to be joined together in these steps is most
commonly $2$, occasionally $3$, and rarely $4$. The is also a subtle
dependency on $k_{\mathrm{den}}$ as well, however the number of merging
clusters only starts to increase for small values of this parameter i.e.
$<20$. Given that typically $k_{\mathrm{den}}\geq 20$, we suggest that for the
choice of $k_{\mathrm{link}}$ not to affect this behaviour, its value needs to
be at least as large as $5$ or $6$.

Furthermore, we must also consider that certain arrangements of points may
give arise to the issue whereby the $k_{\mathrm{link}}$ nearest neighbours of
each point is not an extensive enough neighbourhood for the process to be able
to seamlessly aggregate most, if not all, points together. Densely connecting
all points via their neighbourhood in a data set is not always feasible. In
particular, if all points in the data set are contained within dense clusters
that are each sparsely distributed with respect to each other then connecting
all points may require a value of $k_{\mathrm{link}}$ that approaches the size
of the data set itself.

To manage this issue, we choose practical values for $k_{\mathrm{link}}$ by
ensuring that CluSTAR-ND is able to densely connect all points given that they
are from a randomly sampled uniform distribution in $d$-dimensions when it is
using a particular value for $k_{\mathrm{den}}$. We vary both $d$
consecutively from $1$ to $9$ and $k_{\mathrm{den}}$ within the set $\\{20,\
30,\ 40,\ 60,\ 80,\ 120,\ 160,\ 240,\ 360,\ 480\\}$. For each
$d$-$k_{\mathrm{den}}$ combination we create $300$ random samples of $10^{5}$
points from a $d$-dimensional uniform distribution within the unit hypercube.
For each of these we find the smallest value of $k_{\mathrm{link}}$ that will
densely connect all $10^{5}$ points together by the end of the aggregation
process. The histograms of these smallest values of $k_{\mathrm{link}}$ can be
seen in the top panel of Fig. 6 for various combinations of $d$ and
$k_{\mathrm{den}}$. Within these histograms a value of $k_{\mathrm{link}}=x$
will contribute to the height of the bar in the interval $(x-1,x]$
777Constructing the histogram in this way helps to provide a continuous
analogue of the probability distribution since in such an analogue the value
of $k_{\mathrm{link}}=x$ will have been drawn from this interval..

We then fit a skew-normal such that the sum of squared differences between the
probability within each of the histogram’s intervals and the probability of
the skew-Gaussian within those same intervals is minimised. The skew-normal
distribution we fit has the standard form of

$\displaystyle
f(k_{\mathrm{link}})=\frac{2}{\omega}\phi\bigg{(}\frac{k_{\mathrm{link}}-\xi}{\omega}\bigg{)}\Phi\bigg{(}\alpha\frac{k_{\mathrm{link}}-\xi}{\omega}\bigg{)},\
\mathrm{where}$
$\displaystyle\phi(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}},\
\mathrm{and}$
$\displaystyle\Phi(x)=\int_{-\infty}^{x}\phi(y)dy=\frac{1}{2}\bigg{[}1+\mathrm{erf}\bigg{(}\frac{x}{\sqrt{2}}\bigg{)}\bigg{]}.$
(15)

Here, $\xi$, $\omega$, and $\alpha$ are the location, shape, and scale
parameters respectively. The continuous distributions shown in the bottom
panel of Fig. 6 depicts the fitted skew-normal distributions for various
combinations of $d$ and $k_{\mathrm{den}}$.

We now use the fitted skew-normal distributions to find the continuous
analogue value for $k_{\mathrm{link}}$ that predicts that all $10^{5}$ points
within each uniform distribution will be densely connected together $95$% of
the time. Once found, we then fit a function of the form $k_{\mathrm{link},\
.95}=\alpha d^{\beta}+\gamma{k_{\mathrm{den}}}^{\delta}+\epsilon$ by
minimising the squared differences between these distributional
$95^{\mathrm{th}}$ percentile values and the function’s output – finding that
$\alpha\approx 12.0$, $\beta\approx-2.2$, $\gamma\approx-23.0$,
$\delta\approx-0.6$, and $\epsilon\approx 9.0$. The predicted continuous
analogue value of $k_{\mathrm{link},\ .95}$ are shown as dashed vertical lines
in the bottom panel of Fig. 6 for various combinations of $d$ and
$k_{\mathrm{den}}$. We then convert these back to the original discrete regime
by rounding up to the nearest integer. Likewise, the discrete values of
$k_{\mathrm{link},\ .95}$ are shown as dashed vertical lines in the top panel
of Fig. 6 for various combinations of $d$ and $k_{\mathrm{den}}$.

The fitted values of $k_{\mathrm{link},\ .95}$ in the continuous analogue do
not always correspond well to the true $95^{\mathrm{th}}$ percentile of the
fitted skew-normal distributions, however we find that they do correspond well
to the empirical $95^{\mathrm{th}}$ percentile of histograms produced over
each combination of $d$ and $k_{\mathrm{den}}$. To be more certain that the
automation of $k_{\mathrm{link}}$ will not artificially break densely
connected regions we add $1$ to these discrete values of the
$95^{\mathrm{th}}$ percentile and – to keep the potential cluster merging
behaviour unaffected as discussed above – always ensure that it is larger than
$7$. To be clear, unless the user specifies a value for $k_{\mathrm{link}}$
CluSTAR-ND calculates it using automatically using

$k_{\mathrm{link}}=\mathrm{max}\\{\mathrm{ceil}(12.0d^{-2.2}-23.0{k_{\mathrm{den}}}^{-0.6}+10.0),\
7\\}.$ (16)

Constructing $k_{\mathrm{link}}$ this way standardises the degree with which
it is possible to connect the points of a data set whilst also ensuring
maximal cluster completeness and as such these calculated values serve as a
practical lower limit for $k_{\mathrm{link}}$. The user may choose different
values for $k_{\mathrm{link}}$, although there are some considerations. The
value of $k_{\mathrm{link}}$ must always be smaller than or equal to that of
$k_{\mathrm{den}}$ – which applies to the user’s choice of $k_{\mathrm{den}}$
just as it does to that of $k_{\mathrm{link}}$. In the event that
$k_{\mathrm{link}}$ is chosen automatically by CluSTAR-ND the above optimal
values are effectively the lower limits of $k_{\mathrm{den}}$ given the number
of dimensions $d$. The value of $k_{\mathrm{link}}$ could be chosen to be
larger than the automated values for the purpose of decreasing the run-time,
however this comes at the detriment of cluster completeness. Other than these
considerations, the optimal values of $k_{\mathrm{link}}$ do not depend on the
choice of $k_{\mathrm{den}}$ nor any other parameter from CluSTAR-ND as they
have been determined via the simple criterion that all points in a
$d$-dimensional uniform distribution be densely connected.

Figure 6: The empirical distributional relationships between the minimum
$k_{\mathrm{link}}$ required to connect all $10^{5}$ points of a
$d$-dimensional uniform distribution contained to the unit hypercube given
that $k_{\mathrm{den}}$ nearest neighbours are used to estimate the local
density. The value of $d$ is varied from $1$ to $9$ and the value of
$k_{\mathrm{den}}$ is varied within the set $\\{20,\ 30,\ 40,\ 60,\ 80,\ 120,\
160,\ 240,\ 360,\ 480\\}$ – however only specific values of the latter are
shown. The histograms in the top panel are the binned minimum
$k_{\mathrm{link}}$ values and are created from $300$ re-samplings of the
aforementioned data distribution for each $d$ and $k_{\mathrm{den}}$
combination. The skew-Gaussian distributions in the bottom panel are the
fitted continuous analogues of the histograms in the top panel. The dashed
lines in the lower panel are the $95^{\mathrm{th}}$ percentile values of the
fitted skew-Gaussian distributions and the dashed lines in the top panel are
the discrete analogue of the latter. Figure 7: The clustering power of
CluSTAR-ND as it depends on the various feature subspaces, $k_{\mathrm{den}}$,
and $\rho_{\mathrm{threshold}}$. The objective function used here is defined
in Eq. 5.1.1 however the plots show a Gaussian-smoothed version of this for
simplicity. The dashed lines indicate the fitted values of
$\rho_{\mathrm{threshold}}$ as a function the number of features, $d$, and the
value of $k_{\mathrm{den}}$. These fitted values are calculated using Eq. 17.

###  6.2 Values for $\boldsymbol{\rho_{\mathrm{threshold}}}$

The parameter $\rho_{\mathrm{threshold}}$ is used to provide the functionality
of determining which groupings of points are significantly dense with respect
to their surrounds. More specifically, it is responsible for ensuring that the
median density of any cluster is at least $\rho_{\mathrm{threshold}}$ times
the density of that cluster’s surrounds – refer to Eq. 3.6.

For its use in CluSTAR-ND, we take a closer look at the
$\rho_{\mathrm{threshold}}$ parameter by optimising it on the synthetic
galaxies described in Sec. 4. Since the labels within these data sets only
provide a flat clustering of the synthetic galaxies without noise – instead of
a hierarchical clustering with noise as is the case with CluSTAR-ND’s output –
we only use the leaf clusters found by CluSTAR-ND in order to determine its
performance with varying $\rho_{\mathrm{threshold}}$. This approach then
requires $f_{\mathrm{reject}}$ to be sufficiently large – to ensure that the
largest leaf cluster shares less than $f_{\mathrm{reject}}$ of it’s points
with the root cluster – so that the leaf clusters will be unaffected by the
choice of $f_{\mathrm{reject}}$. In fact, this poses a practical lower limit
to $f_{\mathrm{reject}}$, however we revisit this with greater detail in Sec.
6.4. Here, we set $f_{\mathrm{reject}}=1$. Since the effect of
$S_{\mathrm{outlier}}$ on the final clusters is only minimal, we simply set
$S_{\mathrm{outlier}}=\infty$ to remove its influence – we take a closer look
at $S_{\mathrm{outlier}}$ in Sec. 6.3. The $k_{\mathrm{link}}$ parameter is
automatically chosen and accordingly to the scheme outlined in Sec. 6.1. For
simplicity, we also optimise these parameters using $adaptive=1$.

By setting the above parameters in this way, we are able to find the optimal
values of $\rho_{\mathrm{threshold}}$ as they depend on $k_{\mathrm{den}}$ and
the number of features, $d$. We vary $\rho_{\mathrm{threshold}}$ from $1$ to
$4.5$ in intervals of $0.1$ and $k_{\mathrm{den}}$ again from the set of
$\\{20,\ 30,\ 40,\ 60,\ 80,\ 120,\ 160,\ 240,\ 360,\ 480\\}$. For each
combination of $\rho_{\mathrm{threshold}}$ and $k_{\mathrm{den}}$, we run
CluSTAR-ND over each of the various feature subspace combinations out of the
spatial ($\bf{x}$), kinematic ($\bf{v}$), and chemical ($\bf{m}$)
subspaces888Note that we do not analyse the clusterings over the
$\bf{m}$-subspace alone as typically each data point shares the same values of
both [Fe/H] and [$\alpha$/Fe] with a number of other data points. Producing
clusterings over the $\bf{m}$-subspace thereby artificially predicts many
thousands of spurious clusters that only exist in this context because each
point they contain is effectively a duplicate – this dramatically increases
the estimation of their local density and the density contrast between nearby
neighbourhoods. It is likely that this will also affect the clusterings in
higher dimensional feature spaces that contain the $\bf{m}$-subspace to some
extent. However, since the number of points per value of the ([Fe/H],
[$\alpha$/Fe]) tuple is small (compared to the total number of points) in each
synthetic galaxy, the global Mahalanobis metric used in CluSTAR-ND here is
effectively blind to the fact the data sets in higher dimensions are made up
many hyper-planes separated along the [Fe/H] and [$\alpha$/Fe] directions.
Specifically, the blindness that CluSTAR-ND has to this bias in higher
dimensional settings is a result of the PCA transformation where the resultant
inter-point spacing within each local hyper-plane is on the order of the
inter-plane spacing between each hyper-plane.. Then for every run, we compare
the leaf clusters found by CluSTAR-ND to the ground truth labels from the
relevant synthetic galaxy. We use the mutual-information-based objective
function, $F(\mathcal{S})$, introduced in Sec. 5.1.1 and Eq. 5.1.1 to assess
the clustering power of CluSTAR-ND for each $\rho_{\mathrm{threshold}}$ –
$k_{\mathrm{den}}$ – subspace combination.

Fig. 7 contains a panel for each feature subspace combination whereby
$F(\mathcal{S})$ is plotted against the values of $\rho_{\mathrm{threshold}}$
for each value of $k_{\mathrm{den}}$ (with specific values shown in colour).
As is mentioned in Sec. 4, the spatial and kinematic features (${\bf x}$ and
${\bf v}$) are each $d=3$ subspaces and the chemical features (${\bf m}$) is a
$d=2$ subspace. From these plots we can clearly see a inter-dependency between
$\rho_{\mathrm{threshold}}$, $k_{\mathrm{den}}$, and the number of features,
$d$. As such, we wish to find a function that describes the optimal values for
$\rho_{\mathrm{threshold}}$ such that this function is monotonically
increasing with $d$ and monotonically decreasing with $k_{\mathrm{den}}$. We
also need for $\rho_{\mathrm{threshold}}$ to always be at least $1$ as there
is no reason for choosing otherwise.

In light of these factors and $\rho_{\mathrm{threshold}}$’s coupled dependency
on $k_{\mathrm{den}}$ and $d$ (unlike the inter-dependency between these
variables and $k_{\mathrm{link}}$ investigated in Sec. 6.1), we choose to fit
a function of the form $\rho_{\mathrm{threshold}}(k_{\mathrm{den}},d)=1+\alpha
f(d)/g(k_{\mathrm{den}})$ where $f$ and $g$ are monotonically increasing and
strictly non-zero for positive inputs. We consider each power-law/logarithmic
form combination for $f$ and $g$, but ultimately find that setting
$f(d)=d^{\beta}$ and $g(k_{\mathrm{den}})=\mathrm{ln}(k_{\mathrm{den}})$ gives
the best overall fit. In fitting, we maximise the function’s output on all
Gaussian-smoothed objective function curves (for all values of
$k_{\mathrm{den}}$ not just those that are coloured) shown in Fig. 7
simultaneously with an equal weighting placed on each. Following this, we find
that the optimal values of $\rho_{\mathrm{threshold}}$ are well-described by
the function

$\rho_{\mathrm{threshold}}(k_{\mathrm{den}},d)=1+0.81\frac{d^{0.81}}{\mathrm{ln}(k_{\mathrm{den}})}.$
(17)

Some fitted values produced by Eq. 17 are shown in Fig. 7 as vertical dashed
lines with colours that correspond to the appropriate value of
$k_{\mathrm{den}}$. While these fitted values are not always at the global
maxima of the objective function curves, the function does not sacrifice much
relevant information for the added convenience of having this parameter be
chosen automatically from the user’s choice of $k_{\mathrm{den}}$ and the
number of features, $d$, in the data set that the user wishes to produce a
clustering from.

In addition to being able to justify a set of optimal choices for
$\rho_{\mathrm{threshold}}$, we can also see from Fig. 7 the difference in
information content as it depends on the specific feature subspace combination
present within the input data. Here we can easily see that including extra
informative dimensions will increase clustering power. We can also see that
smaller $k_{\mathrm{den}}$ values produce better clusterings in an
astrophysical context since they increase the resolution with which CluSTAR-ND
can resolve structure. This is particularly prevalent in the larger feature
spaces where fine-structure can be separated from background noise more easily
– however choosing small values for $k_{\mathrm{den}}$ comes with diminishing
returns as this also increases the effects of Poisson noise on the density
estimation which can artificially introduce spurious structures. We take a
closer look at the robustness of the CluSTAR-ND output by comparing recovery
and purity statistics of specific clusters across these same feature subspace
combinations in Sec. 7.

Figure 8: The clustering power of CluSTAR-ND as it depends on the various
feature subspaces, $k_{\mathrm{den}}$, and $S_{\mathrm{outlier}}$. The
objective function used here is defined in Eq. 5.1.1 and is re-normalised
between each curve’s respective minimum and maximum over this range for
$S_{\mathrm{outlier}}$ to better depict the global influence of this
parameter. As with Fig. 7, these plots show a Gaussian-smoothed version of the
objective function. We see that for small $S_{\mathrm{outlier}}$ values the
clustering power decreases dramatically and for large values the clustering
power asymptotically approaches a maximum (values of which are those shown in
Fig. 7 at the relevant $\rho_{\mathrm{threshold}}$ values for each
$k_{\mathrm{den}}$ and number of features, $d$). The dashed lines indicate the
value of $S_{\mathrm{outlier}}=2.5$ which we set as the standard choice for
CluSTAR-ND. This value is both small enough to remove distinct outliers from
clusters while being large enough to not have the resultant clustering
breakdown. Ultimately however, any value above about
$S_{\mathrm{outlier}}=1.5\sim 2$ will produce good clustering results –
adjusting this further is dependent on how strict the user wishes to be on the
removal of potential outliers from clusters.

###  6.3 Values for $\boldsymbol{S_{\mathrm{outlier}}}$

The parameter $S_{\mathrm{outlier}}$ is used to categorise the points of the
input data as either local inliers or local outliers depending on whether
their local-outlier-factor (or rather its density kernel analogue in Eq. 8) is
less than or greater than $S_{\mathrm{outlier}}$. Following the aggregation
process described in Sec. 3.5, a cut-off density is then defined for each
cluster according to Eq. 9. This classification is then used to remove all
points from each cluster whose local density is less than the cut-off density,
$\rho_{\mathrm{cut}}$, of that cluster – i.e. only local outliers whose
density is less than that of all local inliers are removed.

To better understand the effect $S_{\mathrm{outlier}}$ has on the final
clustering produced and gauge an optimal value for it from the synthetic
galaxies discussed in Sec. 4, we compute the objective function,
$F(\mathcal{S})$, in Eq. 5.1.1 from the CluSTAR-ND output for various
combinations of $S_{\mathrm{outlier}}$, $k_{\mathrm{den}}$, and feature
subspace in largely the same way as in Sec. 6.2. We vary
$S_{\mathrm{outlier}}$ from $1$ to $4.5$ in intervals of $0.1$. We similarly
vary $k_{\mathrm{den}}$ from the set of $\\{20,\ 30,\ 40,\ 60,\ 80,\ 120,\
160,\ 240,\ 360,\ 480\\}$ and the feature subspace combinations of the spatial
(x), kinematic (v), and chemical (m) subspaces. Given these combinations of
$k_{\mathrm{den}}$ and the number of features, $d$, we allow for
$k_{\mathrm{link}}$ and $\rho_{\mathrm{threshold}}$ to be automatically chosen
according to Eqs. 16 and 17 respectively. In our investigation of
$S_{\mathrm{outlier}}$ we again set $f_{\mathrm{reject}}=1$ and $adaptive=1$.

The relationship between the clustering power of CluSTAR-ND and
$S_{\mathrm{outlier}}$ is shown in Fig. 8 for each combination of
$k_{\mathrm{den}}$ and $d$. More simplistic than that between the clustering
power and $\rho_{\mathrm{threshold}}$, the relationship here is very similar
between these combinations and suggests that so long as $S_{\mathrm{outlier}}$
is large enough the CluSTAR-ND output will be robust999The $F({\bf v})$ panel
of Fig. 8 appears to include a oscillation-like pattern for
$S_{\mathrm{outlier}}\gtrsim 1.5$ and $k_{\mathrm{den}}=20$. This is an
artefact of the objective function minimum being closer to its maximum for
this combination. As such, the scale of the noise appears much larger and
noticeable compared to the other $k_{\mathrm{den}}$ and $d$ combinations. The
exact reasoning for this noise is that for a given change in
$S_{\mathrm{outlier}}$ the quality of some clusters will increase while others
may decrease. These competing effects take place for every curve in Fig. 8 but
are simply more noticeable on this curve due to the imposed scale between the
minimum and maximum.. For values larger than about $1.5$ or $2$, the choice of
$S_{\mathrm{outlier}}$ has a minimal effect and will only serve to remove a
select few or perhaps none of the points from each cluster. Of course changing
$S_{\mathrm{outlier}}$ from $3$ and $4$ will still have an effect on which
points are apart of the final clusters, however this effect is minimal and
cannot be robustly defined beyond the resultant clustering include more
potential outliers per cluster than beforehand. In light of this ambiguity, we
choose to use $S_{\mathrm{outlier}}=2.5$ as the standard value within CluSTAR-
ND as this will provide strong clustering results with a moderate level of
outlier detection.

###  6.4 Values for $\boldsymbol{f_{\mathrm{reject}}}$

The parameter $f_{\mathrm{reject}}$ defines the maximum fraction of points
(relative to the parent cluster) that can be shared between any of the parent-
child cluster pairs that exist immediately following the aggregation process.
A cluster is rejected from a parent-child pair if the child shares more than
$f_{\mathrm{reject}}$ of the parent’s points – according to Eq. 3.6 – i.e.
CluSTAR-ND removes the smallest cluster of the pair that has child clusters of
its own, unless it would be a root-level cluster in which case the child from
the pair is removed. These rules mean that $f_{\mathrm{reject}}$ is implicitly
responsible for determining the size and shape of the hierarchy between the
terminating clusters of the hierarchy. Refer to Sec. 3.6 for more details on
the rules surrounding the $f_{\mathrm{reject}}$ parameter.

Choosing extreme $f_{\mathrm{reject}}$ values of $0$ and $1$ would force
CluSTAR-ND to return $(\leq$101010The $\leq$ sign here indicates that the
number of leaf clusters found in the $f_{\mathrm{reject}}=0$ case may be less
than or equal to the number of leaf clusters found in the
$f_{\mathrm{reject}}=1$. The inequality will exist if there is one or more
leaf clusters whose parent is the root-level cluster.$)l$ and $2(l-1)$
clusters per root-level cluster respectively – where $l$ is the number of leaf
clusters within that root-level cluster. In the same scenarios, the hierarchy
returned by CluSTAR-ND would have depths of $1$ and somewhere in the range of
[$\mathrm{log}_{2}(l-1)$, $l-1$] per root-level cluster respectively. To use
the standard terminology, these values of $f_{\mathrm{reject}}$ create
hierarchies that form perfect $(\leq)l$-ary trees and full binary trees
respectively – with the exact shape of the latter dependent upon $l$’s
representation in base-$2$ and the clustering structure within the data set.
Varying $f_{\mathrm{reject}}$ between these values transforms the shape of the
hierarchy in a way is difficult to predict.

We should expect that a reasonable lower limit for $f_{\mathrm{reject}}$ is
$0.5$, since anything lower would suggest that equally sized groups can not be
merged into a parent cluster. However, further constraining this within this
paper would be a difficult task since the satellite labels within the
synthetic data from Galaxia (Sec. 4) provide a flat clustering without noise
of each galaxy and this information alone does not indicate the optimal
hierarchy shape. To properly determine the optimal $f_{\mathrm{reject}}$ from
this data, we would need to assess the physics of the galaxies and their
constituents – to determine how bound each satellite is to each other, and
prepare an objective function that can be used to infer the most physical
sense of what a satellite within a satellite should look like in an
astrophysical context. Specifically, this would require analysis of the mass,
orbits of the points, backwards integration, and analysis of the star
formation histories within satellites and how this relates to each other
satellite within a galaxy. Additionally, we can not even be certain from the
outset that we would see a single most-optimal value for $f_{\mathrm{reject}}$
that is able to consistently provide such a hierarchy because the rules that
govern this may not even be suitable to do so.

In light of these issues – and since we expect the hierarchy produced by
CluSTAR-ND to be similar to that of Halo-OPTICS– we choose to follow the
verdicts with respect to the parameter $f_{\mathrm{reject}}$ in Oliver et al.
(2020). Considering the results from figure 6 therein, a value of
$f_{\mathrm{reject}}=0.9$ is near-optimal and ensures good quality clustering
results (high Jaccard index as well as high recovery and purity fractions) of
the hierarchical combinations of hand-crafted leaf clusters. From a hierarchy
interpret-ability stand point, running CluSTAR-ND with this value of
$f_{\mathrm{reject}}$ over the synthetic galaxies discussed in Sec. 4 creates
consistently shaped hierarchies with levels $L=1$, $2$, and $3+$ containing
$\sim 70$%, $\sim 20$%, and $\sim 10$% of all galactic substructure clusters
($L>0$). Furthermore, these percentages do not seem to change more than $\sim
15$% regardless of the value of $k_{\mathrm{den}}$, the feature space of the
input data, nor the choice synthetic galaxy – implying that we should expect
the resultant hierarchy shape to remain simple to interpret and predictable
when CluSTAR-ND is applied to astrophysical data sets.

##  7 Information Content of Different Clustering Scenarios

In the preceding sections we have presented the CluSTAR-ND algorithm (Sec. 3),
showing that it produces comparable clusters to Halo-OPTICS in a far more
computationally efficient manner (Sec. 5), and optimised its parameters for
generalised application to n-dimensional data sets (Sec. 6). In this section
we will further demonstrate the clustering power of CluSTAR-ND and quantify
the difference between the clusters that result when CluSTAR-ND is applied to
different galactic environments (Sec. 7.1), to different feature spaces (Sec.
7.2), and using metric adaptivity settings (Sec. 7.3).

Table 1: A summary of the minima, medians, means, and maxima of the normalised mutual information values ($F_{g}(\mathcal{S})$) for the various feature spaces ($\mathcal{S}$) and $adaptive$ settings shown within Fig. 9. The mean values are equivalent to the objective function defined $F(\mathcal{S})$ in Eq. 5.1.1. For completeness, we also include these statistics for the clusterings produced with Halo-OPTICS using the $({\bf x})$ feature space in Sec. 5. $\mathcal{S}$ | $adaptive$ | Min | Median | Mean | Max
---|---|---|---|---|---
$({\bf x})$ | Halo-OPTICS | 0 | 0.027 | 0.043 | 0.198
$0$ | 0.001 | 0.028 | 0.043 | 0.197
$1$ | 0.006 | 0.027 | 0.043 | 0.193
$2$ | 0 | 0.023 | 0.039 | 0.193
$({\bf v})$ | $0$ | 0 | 0.009 | 0.029 | 0.261
$1$ | 0 | 0.008 | 0.029 | 0.266
$2$ | 0 | 0.006 | 0.028 | 0.265
$({\bf x},{\bf m})$ | $1$ | 0.024 | 0.062 | 0.079 | 0.364
$2$ | 0.017 | 0.066 | 0.075 | 0.341
$({\bf v},{\bf m})$ | $1$ | 0 | 0.063 | 0.093 | 0.473
$2$ | 0 | 0.070 | 0.102 | 0.469
$({\bf x},{\bf v})$ | $1$ | 0 | 0.065 | 0.089 | 0.445
$2$ | 0 | 0.067 | 0.092 | 0.443
$({\bf x},{\bf v},{\bf m})$ | $1$ | 0 | 0.077 | 0.110 | 0.485
$2$ | 0 | 0.083 | 0.109 | 0.486
Figure 9: The normalised mutual information $F_{g}(\mathcal{S})$ from Eq.
5.1.1 of various CluSTAR-ND clustering scenarios involving different synthetic
galaxies, feature spaces, and metric adaptivity settings. The
$F_{g}(\mathcal{S})$ values of both the Halo-OPTICS clusterings from Sec. 5
and a hypothetical set of clusterings that perfectly classify
($J_{\mathrm{max}}=1$) only the self-bound satellites within each MW type
galaxy are also shown for reference. The vertical axes indicate the names of
the synthetic galaxies as well as whether they are from the original MW type
galaxies of Bullock & Johnston (2005) or from the additional artificial
galaxies of Johnston et al. (2008). It is easy to see that the galactic
environ, the number and type of informative features, and the metric
adaptivity all affect the power of the resultant clustering.

The clustering power of the now-optimised CluSTAR-ND is summarised in Tab. 1 &
Fig. 9: Fig. 9 illustrates the breakdown of the objective function from Eq.
5.1.1 (i.e. $F_{g}(\mathcal{S})$) over each synthetic galaxy from Sec. 4 – the
measure of normalised mutual information is each galaxy’s contribution to the
objective function prior to having averaged over them; Tab. 1 displays the
minima, medians, means, and maxima of the distribution of $F_{g}(\mathcal{S})$
as it depends on the metric adaptivity setting and the feature space. With
Fig. 9 we can also see how the clustering power differs depending on the
galactic environment as well. While these values may seem low, both Tab. 1 and
Fig. 9 show that CluSTAR-ND has performed well given the circumstances of this
data. As a reference point, Fig. 9 also shows the clustering power of a
hypothetical clustering that exclusively and perfectly classifies
($J_{\mathrm{max}}=1$) the self-bound satellites from within each MW type
galaxy. We see that in comparison, CluSTAR-ND performs well and that the
galaxies are composed of mostly disrupted tidal debris from previous mergers –
which is mostly impossible to aggregate together without joining multiple
unbound groups and effectively reducing our understanding of the catalogue of
true satellites. We see that even with CluSTAR-ND only using the positions
($\mathcal{S}=({\bf x})$) it will consistently return the self-bound
satellites and some tidal debris to the effect of having an
$F_{g}(\mathcal{S})$ value that is at least a factor of $\sim 1$–$2$ times
larger than for just classifying the self-bound satellites. This is
unsurprising since Halo-OPTICS was shown to be able to do the same thing with
only positional information (Oliver et al., 2020). As the feature space size
increases CluSTAR-ND’s clustering power increases too such that when
$\mathcal{S}=({\bf x},{\bf v},{\bf m})$, $F_{g}(\mathcal{S})$ is typically at
least a factor of $\sim 3$–$6$ times larger than that for self-bound
satellites only.

All clusters analysed in this section have been produced using
$k_{\mathrm{den}}=20$ – as this consistently produced the best results in Sec.
6. All other parameters are chosen automatically according to the analysis in
Sec. 6 unless specified otherwise. We also highlight the impact of feature
space and adaptivity setting using two satellites, shown in Fig. 10 and 11.

Figure 10: A satellite and the corresponding best-fitting clusters predicted
by CluSTAR-ND when using various feature spaces and metric adaptivity
settings. The projection shown of the satellite in each panel is of the two
spatial PCA components with the largest variances. The feature space used for
the satellite’s prediction is shown in the top left corner of each panel. The
metric adaptivity setting is indicated by colour; where blue and red
corresponds to an $adaptive$ value of $1$ and $2$ respectively – the purple
indicates their intersection. Large points belong to the satellite, small
points are incorrectly associated to the satellite. In the top right corner of
each panel are the maximum Jaccard index (Eq. 5.1), recovery, and purity (Eq.
5.1) of the best-fitting clusters from the two $adaptive$ values. These panels
show the influence of the different feature spaces as well as that of the
different metric adaptivity settings.

###  7.1 The Effects of Galactic Environment

Fig. 9 shows that there are distinct differences in CluSTAR-ND’s clustering
power that are dependent upon the type of galaxy that the algorithm is applied
to. This is not unexpected, the constituent satellites of any particular
galaxy will be more/less difficult to classify for any clustering algorithm if
a large number of these satellites are disrupted/intact. In essence, if a
large number of a satellite’s points’ $k_{\mathrm{link}}$ neighbourhoods
contain points from multiple satellites then CluSTAR-ND’s ability to provide a
well-matched cluster to that satellite diminishes as any matching predicted
cluster will also contain contamination points from other satellites. Such
galactic environments will occur if the constituent satellites have
experienced strong tidal forces throughout their history or if they are
sparsely distributed at the time of their infall such that they encompass a
large volume of the feature space. It is for these reasons that we see the
clustering power of CluSTAR-ND to be less when applied to galaxies created by
old satellites, satellites on radial trajectories, and low luminosity
satellites.

Contrarily, if most $k_{\mathrm{link}}$ neighbourhoods contain only points
from a single satellite then the aggregation process can likely link together
many points from a single satellite in sequence before the group becomes
overly contaminated by members of other satellites. Of course, CluSTAR-ND’s
ability to create meaningful clusters from these neighbourhoods also depends
on the density of the points and whether these neighbourhoods are isolated or
connected to other like-neighbourhoods. Galactic environments that harbour
such conditions will contain a larger proportion of intact satellites compared
to those that do not. Accordingly, Fig. 9 shows that galaxies made of young
and (for some feature spaces) high luminosity satellites are more easily
categorised into their respective satellites.

As the MW type galaxies are essentially a mix of disrupted and intact
satellites we befittingly see a clustering power that lies somewhere in the
middle of these extremes. Nevertheless, from these observations we can gain
some understanding about what to expect if CluSTAR-ND were to be applied to a
data set of stars in the galactic bulge compared to if it were applied to the
stars in the stellar halo.

###  7.2 The Effects of Feature Space

The effects of group-mixing in a data set are not only prevalent between
different galactic environments, but are noticeable between different feature
spaces as well. Generally speaking, Tab. 1 portrays that the CluSTAR-ND
clustering power is greater for larger feature spaces. One exception to this
is that we find the $({\bf v},{\bf m})$ feature space more informative than
the $({\bf x},{\bf v})$ feature space. Amongst equal sized feature spaces, we
see that spatial coordinates are generally more informative than kinematic
coordinates although the opposite is true when chemical abundances are added
in amongst these feature spaces.

From Fig. 9 however, we see that the above observations from Tab. 1 are not
consistent across galaxies – i.e. the information content of any particular
galaxy clustering is coupled to both the galactic environment and the
underlying feature space of the data. Out of the synthetic galaxies we
investigate in this paper, the clustering scenario that gives the most
information content for each galaxy is most commonly derived from the largest
feature space, $({\bf x},{\bf v},{\bf m})$, however the title of most
informative feature space is also sometimes held by the $({\bf x},{\bf m})$,
$({\bf v},{\bf m})$, or $({\bf x},{\bf v})$ feature spaces.

Adding in extra features will often increase clustering power by creating some
separation between multiple satellites within the extra volume that is
provided by the additional dimensions. However, this will not always be the
case as it is possible that by adding in extra uninformative features the
points belonging to multiple satellites will be brought closer together
(relative to other points) – making it more difficult for CluSTAR-ND to
disentangle them from each other and their surrounds. These effects are why
the clustering of the satellite depicted in Fig. 10 is worse for the $({\bf
x},{\bf v},{\bf m})$ feature space than it is for the $({\bf x},{\bf v})$
feature space, even though the former feature space is more voluminous. The
same effect can be also be seen in Fig. 11, not only between these feature
spaces but between the $({\bf x})$ and $({\bf x},{\bf m})$ feature spaces as
well.

Irrespective of the above, Tab. 1 does confirm that this effect is outweighed
by the opposite effect when including additional features to cluster over –
i.e. that including additional features brings together $k_{\mathrm{link}}$
neighbourhoods made of points from the same satellite while separating those
that contain points from other satellites. Specifically, we can observe this
occur in Fig. 10 where the addition of the chemical features in the $({\bf
x},{\bf m})$ feature space has dramatically improved the best-fitting
predicted cluster to the satellite when compared to that found using the
spatial features alone. The same can be said about the $({\bf v},{\bf m})$ and
$({\bf v})$ feature spaces, and not only in Fig. 10 but also in Fig. 11.
Moreover, we see that the amalgamation of the spatial and kinematic features
allows CluSTAR-ND to produce far better clusterings of the respective
satellites in both Figs. 10 and 11.

One method of finding the most informative feature space is by comparing the
Shannon entropy (Shannon, 1948) of the available feature spaces. The smaller
the Shannon entropy of a particular feature space the less uniform the
distribution of points within it and hence more clustered – this does not
necessary guarantee that the feature space is appropriately clustered but it
does hint to this possibility. Such a concept is already utilised within
clustering algorithms such as EnLink (Sharma & Johnston, 2009) in order to
create a meaningful locally adaptive metric and could be utilised when using
CluSTAR-ND to narrow the list of possible feature spaces before applying the
algorithm.

Figure 11: A satellite and the corresponding best-fitting clusters predicted
by CluSTAR-ND when using various feature spaces and metric adaptivity
settings. This satellite is different to that shown in Fig. 10, however the
description of each panel’s information remains the same here. As with Fig.
10, these panels show the influence of the different feature spaces as well as
that of the different metric adaptivity settings – however the $({\bf v},{\bf
m})$ panel in the lower middle is a pronounced case of the clustering
difference that can occur between the globally adaptive metric setting
($adaptive=1$) and the locally adaptive metric setting ($adaptive=2$).

###  7.3 The Effects of Metric Adaptivity Settings

Another way to bring like-satellite points together while separating dislike-
satellite points from each other is via the means of a adaptive metric.
CluSTAR-ND offers three metric adaptivity settings however only the $adaptive$
values of $1$ and $2$ should be considered as adaptive metrics since there is
no transformation performed by CluSTAR-ND when $adaptive=0$. It is for this
reason that the latter should only be used if the feature space consists
solely of similar coordinates such as the $({\bf x})$ or $({\bf v})$ feature
spaces. The $adaptive=0$ setting can also be used on compound feature spaces
so long as some meaningful transformation of the data has been performed. The
settings of $adaptive=1$ and $2$ are implementations of a globally adaptive
metric and a locally adaptive metric respectively.

In Tab. 1 we can see that running CluSTAR-ND with the $adaptive=1$ setting
produces similar results to the $adaptive=0$ in the feature spaces where this
setting is applicable. This gives us confidence that the $adaptive=1$ setting
is appropriate for producing a globally adaptive metric for the larger
compound feature spaces. The locally adaptive metric provides a similar
clustering power overall to the globally adaptive metric. The medians and
means of this measure appear to be higher for the locally adaptive metric when
CluSTAR-ND is applied to compound feature spaces, but it is difficult to say
whether using it provides a better clustering in general. However, in Fig. 9
we see that blanket statements about the use of the locally adaptive metric
over the globally adaptive metric are not always true. With this it is clear
that the optimal value of $adaptive$ has a subtle dependency on the galactic
environment as well as the feature space.

Fig. 11 depicts one such case where the use of the $adaptive=2$ setting has
caused CluSTAR-ND to produce a distinctly better match to a satellite. We see
here that by using the locally adaptive metric setting on the $({\bf v},{\bf
m})$ feature space, the best-fitting match to the satellite has both an
improved recovery and purity. However it is possible that in other feature
spaces, $({\bf x},{\bf v})$ and $({\bf x},{\bf v},{\bf m})$ for example, the
lower arm (on the satellite projection) has been found as a sibling cluster(s)
to the best-fitting cluster(s) shown in those panels. In this sense, it is
likely that the locally adaptive metric has simply created a means to condense
the lower density particle bridge between these regions.

Intuitively, the $adaptive=2$ setting is only preferable to the $adaptive=1$
if the most meaningful child cluster (in any parent-child cluster pair) is
dense with respect to the shape of it’s parent cluster. It is not obvious as
to whether this will be the case before using CluSTAR-ND on a particular data
set. Moreover, it is not obvious after the clustering has been found whether
this technique has yielded the a more meaningful clustering overall. Without a
recognisable way of appropriately using the $adaptive=2$ setting111111Other
than perhaps when using the $({\bf v},{\bf m})$ and $({\bf x},{\bf v})$
feature spaces since it is evident from Fig. 9 that the clustering power of
CluSTAR-ND is either similar or increased when using the $adaptive=2$ setting
compared to the $adaptive=1$ on these feature spaces., we can only expect to
see a different clustering by using it over the $adaptive=1$ setting.

Given that using CluSTAR-ND with the locally adaptive metric is more
computationally demanding and requires a longer run-time (by a factor of $\sim
2$) with varied clustering results, it is left to the user’s discretion as to
whether or not they wish to user it and by default the $adaptive$ parameter
will be set to $1$ ensuring a globally adaptive metric. This aspect of the
CluSTAR-ND algorithm can be certainly be improved and we will endeavour to
include a more appropriately determined locally adaptive metric in a future
work. It is also possible to use CluSTAR-ND on the $adaptive=0$ setting in
conjunction with other external manifold learning codes such as UMAP (McInnes
et al., 2018) – however clustering over these results can be unpredictable and
may also require a different $\rho_{\mathrm{threshold}}$ value than the
optimal one we report in this paper.

##  8 Conclusions

We have presented the CluSTAR-ND algorithm and have shown it to be well-suited
to the hierarchical classification of stellar satellite groups within
galaxies. Compared to its predecessor Halo-OPTICS, CluSTAR-ND is not only
capable of producing a similarly robust clustering output but it does so in an
exceptionally efficient manner – outperforming the Halo-OPTICS run-times by
$3$ orders of magnitude at a minimum. Importantly, CluSTAR-ND is readily
applicable to any point-based astrophysical data set that is defined over an
arbitrary number of features. In optimising its clustering performance, we
have removed the need for user defined parameter values such that the most
optimal parameter settings will be chosen automatically based on the input
data (unless the user specifies otherwise).

In this paper we also investigate the capacity of a computationally cheap
design for producing a locally adaptive metric to improve CluSTAR-ND’s
clustering power. The design iteratively re-applies CluSTAR-ND’s substructure
searching component to each consecutive level of the hierarchy in a top-down
fashion in order to adaptively modify the underlying distance metric and
ensure that each level of clusters found are dense with respect to their
parent cluster in the hierarchy. We find that while this approach does produce
better results in a small number of clustering scenarios (predominantly on
data defined with a combination of kinematic and chemical abundance features),
it is mostly inconsistent with regards to the resultant clustering power and
tends to simply provide a different clustering when compared to the globally
adaptive metric scheme. Since the locally adaptive metric scheme requires more
run-time than the globally adaptive metric scheme, we leave using it as a
situational choice for the user to make.

In a future work we will develop a more clustering-appropriate locally
adaptive metric scheme for CluSTAR-ND to make use of. We will also establish a
method of producing fuzzy clusterings from fuzzy data such that the point-
based uncertainties of the data may be propagated into cluster-based
uncertainties. As of this paper however, CluSTAR-ND is ideally suited to large
astrophysical data sets both synthetic and observational, and hence our
further research will also include the application of CluSTAR-ND to
observational data sets of the MW.

## Acknowledgements

WHO gratefully acknowledges financial support through the Paulette Isabel
Jones PhD Completion Scholarship at the University of Sydney. GFL received no
funding to support this research.

## Data Availability

The data underlying this article may be made available on reasonable request
to the corresponding author.

## References

  * Ankerst et al. (1999) Ankerst M., Breunig M. M., Kriegel H.-P., Sander J., 1999, in ACM Sigmod record. ACM, pp 49–60, doi:10.1145/304182.304187
  * Behroozi et al. (2012) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2012, The Astrophysical Journal, 762, 109
  * Bentley (1975) Bentley J. L., 1975, Commun. ACM, 18, 509–517
  * Breunig et al. (1999) Breunig M. M., Kriegel H.-P., Ng R. T., Sander J., 1999. Principles of Data Mining and Knowledge Discovery. Springer Berlin Heidelberg, pp 262–270, doi:10.1145/342009.335388
  * Bullock & Johnston (2005) Bullock J. S., Johnston K. V., 2005, The Astrophysical Journal, 635, 931
  * Campello et al. (2015) Campello R. J. G. B., Moulavi D., Zimek A., Sander J., 2015, ACM Trans. Knowl. Discov. Data, 10
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D., 1985, The Astrophysical Journal, 292, 371
  * Dempster et al. (1977) Dempster A. P., Laird N. M., Rubin D. B., 1977, Journal of the Royal Statistical Society: Series B (Methodological), 39, 1
  * Dunn (1973) Dunn J. C., 1973, Journal of Cybernetics, 3, 32
  * Elahi et al. (2019) Elahi P. J., Canas R., Poulton R. J. J., Tobar R. J., Willis J. S., Lagos C. d. P., Power C., Robotham A. S. G., 2019, Publications of the Astronomical Society of Australia, 36, e021
  * Epanechnikov (1969) Epanechnikov V. A., 1969, Theory of Probability & Its Applications, 14, 153
  * Ester et al. (1996) Ester M., Kriegel H.-P., Sander J., Xu X., 1996, in Kdd. pp 226–231
  * Flores & Primack (1994) Flores R. A., Primack J. R., 1994, ApJ, 427, L1
  * Font et al. (2006) Font A. S., Johnston K. V., Bullock J. S., Robertson B. E., 2006, The Astrophysical Journal, 638, 585
  * Fuentes et al. (2017) Fuentes S. S., De Ridder J., Debosscher J., 2017, Astronomy & Astrophysics, 599, A143
  * Ghigna et al. (1998) Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 1998, Monthly Notices of the Royal Astronomical Society, 300, 146
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Ishiyama et al. (2013) Ishiyama T., et al., 2013, The Astrophysical Journal, 767, 146
  * Jaccard (1912) Jaccard P., 1912, New phytologist, 11, 37
  * Johnston et al. (2008) Johnston K. V., Bullock J. S., Sharma S., Font A., Robertson B. E., Leitner S. N., 2008, The Astrophysical Journal, 689, 936
  * Kauffmann et al. (1993) Kauffmann G., White S. D. M., Guiderdoni B., 1993, Monthly Notices of the Royal Astronomical Society, 264, 201
  * King (1962) King I., 1962, The Astronomical Journal, 67, 471
  * Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, The Astrophysical Journal, 522, 82
  * Knebe et al. (2011) Knebe A., et al., 2011, Monthly Notices of the Royal Astronomical Society, 415, 2293
  * Knebe et al. (2013) Knebe A., et al., 2013, Monthly Notices of the Royal Astronomical Society, 428, 2039
  * Knollmann & Knebe (2009) Knollmann S. R., Knebe A., 2009, The Astrophysical Journal Supplement Series, 182, 608
  * Lam et al. (2015) Lam S. K., Pitrou A., Seibert S., 2015, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC. pp 1–6
  * Lloyd (1982) Lloyd S., 1982, IEEE transactions on information theory, 28, 129
  * MacQueen et al. (1967) MacQueen J., et al., 1967, in Proceedings of the fifth Berkeley symposium on mathematical statistics and probability. pp 281–297
  * Maciejewski et al. (2009) Maciejewski M., Colombi S., Springel V., Alard C., Bouchet F. R., 2009, MNRAS, 396, 1329
  * Mahalanobis (1936) Mahalanobis P. C., 1936.
  * Malhan et al. (2022) Malhan K., et al., 2022, The Astrophysical Journal, 926, 107
  * Maneewongvatana & Mount (1999) Maneewongvatana S., Mount D. M., 1999, in Center for geometric computing 4th annual workshop on computational geometry. pp 1–8
  * McConnachie et al. (2018) McConnachie A. W., et al., 2018, The Astrophysical Journal, 868, 55
  * McInnes et al. (2018) McInnes L., Healy J., Melville J., 2018, arXiv preprint arXiv:1802.03426
  * Moore (1994) Moore B., 1994, Nature, 370, 629
  * Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, The Astrophysical Journal, 524, L19
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Oliver et al. (2020) Oliver W. H., Elahi P. J., Lewis G. F., Power C., 2020, Monthly Notices of the Royal Astronomical Society, 501, 4420
  * Pedregosa et al. (2011) Pedregosa F., et al., 2011, Journal of Machine Learning Research, 12, 2825
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, The Astrophysical Journal, 187, 425
  * Reed et al. (2005) Reed D., Governato F., Quinn T., Gardner J., Stadel J., Lake G., 2005, Monthly Notices of the Royal Astronomical Society, 359, 1537
  * Robertson et al. (2005) Robertson B., Bullock J. S., Font A. S., Johnston K. V., Hernquist L., 2005, The Astrophysical Journal, 632, 872
  * Sain (2002) Sain S. R., 2002, Computational Statistics & Data Analysis, 39, 165
  * Sander et al. (2003) Sander J., Qin X., Lu Z., Niu N., Kovarsky A., 2003, in Whang K.-Y., Jeon J., Shim K., Srivastava J., eds, Advances in Knowledge Discovery and Data Mining. Springer Berlin Heidelberg, pp 75–87
  * Shannon (1948) Shannon C. E., 1948, The Bell System Technical Journal, 27, 379
  * Sharma & Johnston (2009) Sharma S., Johnston K. V., 2009, The Astrophysical Journal, 703, 1061
  * Sharma et al. (2011) Sharma S., Bland-Hawthorn J., Johnston K. V., Binney J., 2011, The Astrophysical Journal, 730, 3
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, Monthly Notices of the Royal Astronomical Society, 328, 726
  * Springel et al. (2008) Springel V., et al., 2008, Monthly Notices of the Royal Astronomical Society, 391, 1685
  * Tollerud et al. (2008) Tollerud E. J., Bullock J. S., Strigari L. E., Willman B., 2008, The Astrophysical Journal, 688, 277
  * Van Den Bosch et al. (2000) Van Den Bosch F. C., Robertson B. E., Dalcanton J. J., De Blok W., 2000, The Astronomical Journal, 119, 1579
  * Vinh et al. (2009) Vinh N. X., Epps J., Bailey J., 2009, in Proceedings of the 26th Annual International Conference on Machine Learning. ICML ’09. Association for Computing Machinery, New York, NY, USA, p. 1073–1080, doi:10.1145/1553374.1553511, https://doi.org/10.1145/1553374.1553511
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, Monthly Notices of the Royal Astronomical Society, 183, 341
  * Zhang et al. (2013) Zhang A. X., Noulas A., Scellato S., Mascolo C., 2013, in 2013 International Conference on Social Computing. IEEE, pp 69–74