FLASC: a flare-sensitive clustering algorithm

Clustering algorithms

Clustering algorithms have been reviewed on many occasions. For example, Ezugwu et al. (2022) list twelve other surveys to start their comprehensive review. A thought-provoking 1971 publication critically reviews clustering theory and practice, much of which remains relevant today (Cormack, 1971). For example, Cormack (1971) explains that clustering techniques are often based on conflicting ideas of what clusters are due to the absence of an agreed-upon formal definition. Generally, clusters are described as groups of similar observations dissimilar to others (Cormack, 1971; Xu & Wunsch, 2005; Ezugwu et al., 2022). Whether clusters should be spherical or if “multidimensional amoebae” shapes are acceptable was similarly contested at the time (Cormack, 1971).

Cormack’s categorisation of clustering algorithms is also still relevant. Generally, clustering algorithms are divided into those that rely on a cluster hierarchy and those that operate by optimising partitions (Xu & Wunsch, 2005; Ezugwu et al., 2022). Below, we provide a brief overview of these categories. We refer to Xu & Wunsch (2005) and Ezugwu et al. (2022) for more elaborate reviews of clustering research.

Structure learning

Many data types lie not just on a manifold but on a smooth, one-dimensional manifold. Extracting such manifolds can be essential in unsupervised learning applications. For example, road networks can be extracted from GPS measurements (Bonnaire, Decelle & Aghanim, 2022), and cell developmental trajectories can be extracted from gene expression data (Qiu et al., 2017; Vandaele, Saeys & De Bie, 2020). Algorithms for extracting such structures are related to our work because the branch-based subgroups we are interested in can be extracted from them by partitioning the data between their intersections (Chervov, Bac & Zinovyev, 2020).

Most work on extracting smooth, one-dimensional manifolds is based on principal curves: a smooth, self-consistent curve that passes through the middle of the data (Hastie & Stuetzle, 1989). Techniques estimating principal curves, trees, or graphs are often based on expectation maximisation (Dempster, Laird & Rubin, 1977) and optimise the one-dimensional manifold directly (e.g., Bonnaire, Decelle & Aghanim, 2022; Mao et al., 2017). Alternative approaches are more closely related to non-linear dimensionality reduction (DR) algorithms that model the data’s structure as an undirected graph (e.g., Roweis & Saul, 2000; Tenenbaum, de Silva & Langford, 2000; Belkin & Niyogi, 2003; van der Maaten & Hinton, 2008; McInnes, Healy & Melville, 2020). For example, Vandaele, Saeys & De Bie (2020) use (manually) pruned minimum spanning trees over edges weighted by their boundary coefficient to extract a graph’s backbone. Alternatively, Ge et al. (2011) extract graph skeletons using a Reeb Graph. Reeb Graphs track connected components’ existence, merges, and splits in level sets of a continuous function defined on a manifold. Using geodesic distances to an arbitrary eccentric point as the continuous function makes the Reeb Graph capture the manifold’s skeleton. Interestingly, Mapper—an algorithm approximating Reeb Graphs (Singh, Memoli & Carlsson, 2007)—has also been used for detecting branch-based subpopulations (Kamruzzaman, Kalyanaraman & Krishnamoorthy, 2018).

There are several similarities between these methods and our work. Like Vandaele, Saeys & De Bie (2020), our approach detects tree-based branching hierarchies. Like Ge et al. (2011), our approach is topologically inspired. Where they create a Reeb Graph, we compute a join tree. The main difference between these methods and our work is their goal. We aim to identify relevant branch-based subpopulations. Ge et al. (2011), Vandaele, Saeys & De Bie (2020), and the expectation maximisation-based algorithms explicitly model the data’s structure, necessitating a higher computational cost.

HDBSCAN*

HDBSCAN* is a state-of-the-art density-based clustering algorithm (Campello, Moulavi & Sander, 2013; Campello et al., 2015). Informally, density-based clustering specifies clusters as regions of high density separated by regions of lower density. This formulation does not limit clusters to convex shapes and provides a natural way to separate noise points from clusters. The algorithm is well suited for exploring unfamiliar data because HDBSCAN* does not require the number of clusters or the distance between clusters to be specified in advance.

Several studies have implemented and adapted the HDBSCAN* algorithm. McInnes & Healy (2017) improved the algorithm’s computational performance by using space trees to find the data points’ nearest neighbours and provided an efficient Python implementation (McInnes, Healy & Astels, 2017). Stewart & Al-Khassaweneh (2022) created a Java implementation with a novel prediction technique for unseen data points. Jackson, Qiao & Xing (2018) presented an approximate HDBSCAN* algorithm that uses NN-descent (Dong, Moses & Li, 2011) to find the nearest neighbours with fast distributed performance. Malzer & Baum (2020) introduced a cluster selection distance threshold that effectively creates a hybrid between DBSCAN’s (Ester et al., 1996) and HDBSCAN*’s cluster selection, improving the algorithm’s performance on data sets with small clusters and a large density variability. Neto et al. (2021) showed how relative neighbourhood graphs (RNGs) (Toussaint, 1980) can be used to efficiently compute HDBSCAN* cluster hierarchies for multiple min. cluster size values. Their follow-up work presented MustaCHE, a visualisation tool for the resulting meta-cluster hierarchy (Neto et al., 2018). To our knowledge, no previous study has adapted HDBSCAN* for detecting flares.

Because our work builds on HDBSCAN*, it is relevant to explain how the algorithm works in more detail. The remainder of this section describes HDBSCAN* following Campello et al. (2015)’s explanation. We refer the reader to McInnes & Healy (2017) for a more formal, statistically motivated description.

The HDBSCAN* algorithm

HDBSCAN* is based on density-based clustering concepts pioneered by Wishart (1969) and formalised by Hartigan (1975). We demonstrate these ideas using a 2D point cloud adapted from McInnes, Healy & Astels (2022), shown in Fig. 1A. In general, let $\mathbf{X}=\{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\}$ be a data set consisting of N feature vectors ${\mathbf{x}}_{(\cdot )}$ and a distance metric $d({\mathbf{x}}_i,{\mathbf{x}}_j)$. Then, HDBSCAN* estimates the density at point ${\mathbf{x}}_i$ as (Campello et al., 2015):

(1) $${\lambda }_k({\mathbf{x}}_i)=1/\kappa ({\mathbf{x}}_i),$$where the point’s core distance $\kappa ({\mathbf{x}}_i)$ is the distance to its $k$-nearest neighbour.

Figure 1B illustrates the example’s density profile as contours in a height map. Density contour clusters intuitively correspond to peaks in the density profile, for example, the clusters indicated in colour in Fig. 1C. More formally, the density contour clusters at some threshold ${\lambda }_t$ are a collection of maximal, connected subsets in a level set $\{\mathbf{x}|\lambda (\mathbf{x})\ge {\lambda }_t\}$ (Hartigan, 1975). In other words, density contour clusters are the connected components of points with a density higher than some threshold. Density contour trees capture the hierarchy in which density contour clusters merge as the density threshold decreases. From a topological perspective, density contour trees are a join tree of the data’s density profile.

Data sets generally do not have an inherent notion of connectivity between their data points. Such connectivity is needed to determine whether two points are part of the same density contour cluster at a threshold ${\lambda }_t$. HDBSCAN* solves this problem by considering points to be connected if the distance between them is smaller than or equal to $1/{\lambda }_t$. This solution is possible because density is defined in terms of distance in Eq. (1). HDBSCAN* uses a mutual reachability distance between points for this purpose, which is defined as (Campello et al., 2015):

(2) $$d_{mreach}({\mathbf{x}}_i,{\mathbf{x}}_j)=\{\begin{array}{cc}max\left\{\kappa ({\mathbf{x}}_i),\kappa ({\mathbf{x}}_j),d({\mathbf{x}}_i,{\mathbf{x}}_j)\right\}\hfill & \mathrm{i}\mathrm{f}{\mathbf{x}}_i\ne {\mathbf{x}}_j,\hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\hfill \end{array}$$where the value of $k$, as used in $\kappa$, is specified manually and acts as a smoothing factor for the density estimation.

We can now recover the density contour tree using the mutual reachability distance to provide connectivity. The edges that change the connectivity between density contour clusters are exactly those edges in the data’s Minimum Spanning Tree (MST) (as cited in Cormack, 1971). HDBSCAN* uses an MST to efficiently compute a single linkage clustering hierarchy (Sibson, 1973). The resulting dendrogram is simplified using a manually specified minimum cluster size $m_c$ to recover a condensed cluster hierarchy that resembles the data’s density profile, as shown in Fig. 1D. From the root down, only the sides of a split containing more than $m_c$ points are considered to represent clusters. Sides with fewer points are interpreted as “falling out of the parent cluster” (McInnes & Healy, 2017) or the cluster disappearing completely.

HDBSCAN* provides two strategies for selecting clusters from the condensed hierarchy: the excess of mass (EOM) strategy and the leaf strategy (Campello et al., 2015). The leaf strategy selects all leaf segments in the condensed hierarchy, typically resulting in multiple small clusters. The EOM strategy maximises relative cluster stability while preventing any data point from being a member of more than one selected cluster. A cluster $C_j$’s relative stability ${\sigma }_k(C_j)$ is defined as Campello et al. (2015):

(3) $${\sigma }_k(C_j)={\sum }_{{\mathbf{x}}_i\in C_j}{\lambda }_{k,max}^{C_j}({\mathbf{x}}_i)-{\lambda }_{k,min}^{C_j},$$where ${\lambda }_{k,max}^{C_j}({\mathbf{x}}_i)$ is the density at which ${\mathbf{x}}_i$ falls out of $C_j$ or $C_j$ separates into two clusters, and ${\lambda }_{k,min}^{C_i}$ is the minimum density at which $C_j$ exists. In other words, the stability of a cluster is the sum of density ranges in which points are part of the cluster, corresponding to the area of the cluster’s icicle in Fig. 1D. HDBSCAN*’s cluster selection epsilon parameter can be used to specify a minimum persistence for EOM clusters (Malzer & Baum, 2020).

Stability

Stability is an important property of algorithms, indicating that their output differs only slightly when the input changes slightly. Two notions of stability are relevant for FLASC: (1) the algorithm has to provide similar results when run repeatedly on (different) samples of an underlying distribution, and (2) the detected branch hierarchies have to represent the clusters’ underlying topology accurately. The deterministic density-based design of FLASC provides stability in the first sense.

Vandaele et al. (2021) analysed the second notion of stability for graph-based branch detection, explaining that the graph approximation should accurately represent the underlying shape and the graph-based centrality function should accurately describe the points’ centrality in a cluster’s metric space $(C_j,d_{mreach})$. For the normalised centrality used by Carlsson (2014), Vandaele et al. (2021) show that the bound on the bottleneck distances between true and empirical persistence diagrams is tight if the metric distortion induced by the graph and its maximum edge weight are small.

Both the full and core cluster approximation graphs used by FLASC satisfy the low maximum edge weight requirement, as their largest edge weight is the minimum mutual reachability distance required for all points in the cluster to be connected in the graph. Additionally, the metric distortion should be small as only edges in the local neighbourhood of data points are included because the clusters do not contain noise points.

Our eccentricity function (Eq. (4)) also meets Vandaele et al.’s requirement to be a $\mathrm{c}$-Lipschitz-continuous function when considered over the cluster centrality graph’s edges:

(6) $$|max\{e({\mathbf{x}}_i),e({\mathbf{x}}_j)\}-max\{e({\mathbf{x}}_k),e({\mathbf{x}}_l)\}|\le \mathrm{c}{d}_{mreach}({\mathbf{x}}_i,{\mathbf{x}}_l),$$where $\mathrm{c}$ is a constant describing the continuity, $({\mathbf{x}}_i,{\mathbf{x}}_j)\in {\mathbf{G}}_k^{C_j}$, $({\mathbf{x}}_k,{\mathbf{x}}_l)\in {\mathbf{G}}_k^{C_j}$, and the mutual reachability between ${\mathbf{x}}_i$ and ${\mathbf{x}}_l$ is the largest of the four points. This property, however, does not guarantee that the detected hierarchy accurately represents the cluster’s topology because Eq. (4) is sensitive to an interplay between the cluster’s shape and the position of its centroid. Consider, for example, a U-shaped cluster. Topologically, this cluster is equivalent to an I-shaped cluster. Equation (4), however, will contain two local eccentricity maxima and three local eccentricity minima because the centroid is located between the U-shape’s arms. As a result, the detected branching hierarchy is indistinguishable from a Y-shaped cluster. We aim to show that the current approach strikes a good balance between computational cost and stability in the experiments presented in the next section.

Computational complexity

The algorithm’s most computationally expensive steps are constructing the full cluster approximation graphs and computing their single linkage hierarchies. Naively, the worst-case complexity for creating a cluster approximation graph is $\mathcal{O}(n_c^2)$, where $n_c$ is the number of points in the cluster. Usually, the average case is much better because the approximation graphs are rarely fully connected. After all, HDBSCAN*’s noise classification limits the density range within the clusters. Furthermore, the space tree re-used from the HDBSCAN* clustering step provides fast asymptotic performance for finding the graph’s edges. The exact run-time bounds depend on the data properties. They are challenging to describe (as explained in McInnes & Healy, 2017). However, an average complexity proportional to $\mathcal{O}(n_e\log N)$ is expected, where $n_e$ is the number of edges in the approximation graph. Computing single linkage hierarchies from the cluster centrality graphs is possible in $\mathcal{O}(n_e\alpha (n_e))$ using the Union-Find implementation from McInnes, Healy & Astels (2017) ( $\alpha$ is the inverse Ackermann function). Like McInnes & Healy (2017), we feel confident that FLASC achieves sub-quadratic complexity on average, which we demonstrate in the ‘Experiments’ section.

Branch detection ability

This first synthetic benchmark compares the branch detection ability of several clustering algorithms. The benchmark is designed to answer the following research question: how well can FLASC detect branches compared to other clustering algorithms?

Six algorithms with Python implementation were selected from different clustering algorithm categories: HDBSCAN* (Campello, Moulavi & Sander, 2013; McInnes & Healy, 2017) represents density-based clustering algorithms and serves as a baseline; Single Linkage Clustering (SLINK) (Sibson, 1973) represents hierarchical clustering algorithms; kMeans (MacQueen, 1967; Lloyd, 1982; Arthur & Vassilvitskii, 2006) represents partitional clustering algorithms; Spectral clustering (Damle, Minden & Ying, 2019) represents graph-based clustering algorithms; and CDC (Peng et al., 2022) represents boundary-seeking clustering algorithms⁴ .

We expect that kMeans, Spectral clustering, CDC, and FLASC will be able to detect branching structures to varying degrees. In contrast, we expect SLINK and HDBSCAN* will struggle due to the lack of low-density regions separating the branches from the clusters. The following two subsections explain how the data for this benchmark was generated and how we measured the algorithms’ branch detection accuracies, respectively.

Datasets

The datasets for this benchmark contain four star-shaped clusters and uniformly distributed noise points (Fig. 4). Each star is parameterised by its number of branches ( $3$, $4$, $5$, $10$), branch lengths ( $1.8$, $2.3$, $2.0$, $3.5$), noise levels ( $0.2$, $0.2$, $0.02$, $0.1$), and points per branch ( $40$, $150$, $20$, $100$). The stars were generated and positioned in 2D. Branches were evenly spread and contain points exponentially spaced from the centre outwards, ensuring density is highest at the stars’ centres and lowest at the branch ends. This configuration prevents the branches from reliably containing local density maxima. Normally distributed noise ( $\mu =0$) was added to the points in $2$, $8$, or $16$ dimensions. The noise level was converted to a standard deviation $\sigma =n_r\sqrt{1/d}$ to correct for the number of dimensions $d$. Uniformly distributed noise points were added to simulate outliers. The number of uniform noise points was $5\mathrm{\%}$ of the clean data set size. All structure is present in the first two dimensions. The additional dimensions contain only noise.

Table 1:

Optimal parameter values per algorithm as found by the grid search described in the ‘Evaluation and settings’ sub-section.

Algorithm	Parameter	Dimensions
		2	8	16
FLASC	Approximation graph	Core	Core	Core
	Cluster selection method	Leaf	Leaf	Leaf
	Min. samples ( $k$)	6	6	2
	Min. cluster size ( $m_c$)	70	100	84
	Min. branch size ( $m_b$)	10	8	14
HDBSCAN*	Cluster selection method	EOM	EOM	Leaf
	Min. samples ( $k$)	2	20	2
	Min. cluster size ( $m_c$)	20	20	58
CDC	Num. neighbours ( $k$)	2	N/A	N/A
	Ratio	0.61	N/A	N/A
kMeans	Num. clusters	27	27	27
SLINK	Num. clusters	27	27	27
Spectral	Num. clusters	27	23	18

DOI: 10.7717/peerj-cs.2792/table-1

Note:

These parameter values were used to compare the algorithms in this benchmark. The CDC implementation only supports 2D data.

Results

Figure 4 shows which clusters the algorithms detected with their optimal parameters (Table 1). The figure also lists the average ARI scores over the five generated data samples. Only FLASC achieved ARI scores greater than $0.61$ and assigned each branch to a distinct group. The other algorithms achieved ARI scores lower than $0.36$. SLINK and HDBSCAN* assigned whole stars to single groups, thereby not detecting most branches. CDC, kMeans, and Spectral clustering did assign multiple groups to the stars, but these groups did not correspond to the branches. Some groups spanned multiple branches, and some branches were given multiple groups. Generally, these findings agree with our expectation that clustering algorithms struggle to detect branches as clusters and highlight the benefit FLASC brings for detecting this type of pattern.

Next, we explore FLASC’s parameter sensitivity using results from the 2D grid search. Figure 5 shows the average ARI curves for the min. samples ( $k$), min. cluster size ( $m_c$), and min. branch size ( $m_b$) parameters. The curves do not show the optimal performance that can be reached by combining the best parameter values—as indicated by the black crosses—because the averages are computed over five data samples and the evaluated (non-optimal) values for the other parameters. Instead, the curves summarise how performance changes by varying the parameters, indicating the algorithm’s sensitivity. Shaded areas around the curves indicate a 95% confidence interval around means.

The average ARI curve for min. samples is shown in Fig. 5A. The curve contains a peak around $k=10$, slightly higher than the optimal value of $k=6$. This pattern indicates that too large values of $k$ reduce the algorithm’s performance. Generally, $k$ should be set so that the $k$-closest neighbours of points within clusters are noticeably closer than the $k$-closest neighbours of noisy points between clusters. Higher values reduce variation and peaks in the modelled density profile, leading to larger clusters. In addition, higher values increase connectivity in the cluster approximation graphs, which reduces the algorithm’s sensitivity to branches. As a result, we tend to set $k$ between $2$ and $25$.

The curves for min. cluster size and min. branch size are shown in Figs. 5B and 5C, respectively. For these parameters, the averages were computed with $k=6$. Both parameters show stable performance once they reach a large enough value to exclude noisy clusters ( $m_c>50$) and branches ( $m_b>3$). Similarly, both parameters should not be set too large, as that would exclude true clusters and branches ( $m_b<37$). Generally, these parameters should be set to the smallest sizes one is interested in finding. The exact sizes vary on a case-by-case basis and may require some domain knowledge to deduce. When larger thresholds are used, the cluster and branching hierarchies contain fewer and larger leaf clusters. So, these parameters can fine-tune clustering results when using the leaf selection strategy.

The different cluster approximation graphs and cluster selection methods have similar curves in Fig. 5. We provide the following guidelines for tuning these parameters. The core cluster approximation graph is appropriate when clusters span large density ranges or contain branching structures at multiple densities. The full approximation graph is more suitable for finding branches at the lowest density within a cluster. The leaf selection strategy works well for finding all density peaks in the data. It tends to find many smaller clusters and can indicate many points as noise. The EOM strategy is more suitable for finding larger clusters containing multiple density peaks. The additional selection parameters can filter out low-persistent clusters and branches when needed.

Finally, we explore FLASC’s practical stability in the sense of its output similarity on the five sampled datasets. The output similarity is explored through unweighted average positions of the detected groups. We refer to these positions as the groups’ centroids. For each ground-truth centroid, we selected the closest centroid FLASC detected with its optimal parameters (Table 1). The centroids are shown in Fig. 6, with the black crosses indicating the ground-truth centroids and the coloured dots indicating their closest detected centroid. The 95 percentile distance to the closest detected centroid is shown as a line to indicate the centroids’ spread and describe the variation in FLASC’s output. Generally, the detected centroids are located close to the ground-truth centroids, and the 95 percentile distances are considerably smaller than the longest branches.

Computational performance

Now, we turn to the computational performance. This second synthetic benchmark is designed to answer the following research question: how does FLASC’s compute cost compare to other clustering algorithms?

The algorithms selected in the previous benchmark are compared on their run time scaling over the data set size and number of dimensions. Given the challenges in accurately benchmarking computational performance (Kriegel, Schubert & Zimek, 2017), we limit this comparison to the trends in run time scaling over data set size and number of dimensions for specific implementations. The following two subsections explain how the data for this benchmark was generated and how we measured the algorithms’ compute cost.

Results

Figure 7 shows the measured run times in seconds over the data set size and number of dimensions. Quadratic regression lines are drawn to visualise the algorithms’ scaling behaviour. The shaded areas around each line indicate the regression’s $95\mathrm{\%}$ confidence interval.

There are several patterns of note in Fig. 7. CDC and Spectral clustering show the steepest quadratic scaling trends in all evaluated dimensions. SLINK exhibits a shallower quadratic scaling trend. kMeans achieved the shallowest scaling trend in all three conditions. HDBSCAN* and both FLASC variants regression lines fall between SLINK and kMeans. In the two dimensional case, their trends are most similar to kMeans. As more dimensions are added, their trends approach SLINK’s trend. This pattern is expected because this dataset is designed to challenge the algorithms’ use of space trees to accelerate the clustering process.

The difference in compute time between the two FLASC variants and HDBSCAN* is small compared to the total compute time, indicating that the branch detection post-processing step is relatively cheap compared to the initial HDBSCAN* clustering step, making it a viable option in practical data analyses workflows.

Use case: diabetes types

Next, we present a data exploration case in which identifying branch-based subgroups is essential to understand the data. Reaven & Miller (1979) attempted to clarify a “horseshoe”-shaped relation between glucose levels and insulin responses in diabetes patients. Three of the metabolic variables they measured were very informative in a 3D scatterplot, showing a dense core with two less-dense branches, which they considered unlikely to be a single population. Seeing that plot was instrumental in their understanding of the data (Miller, 1985).

More recently, Singh, Memoli & Carlsson (2007) used the same data set to demonstrate how Mapper visualises these flares without manually specifying which dimensions to plot. Their analysis leveraged the flares’ lower density, allowing them to be detected with a density-based lens function instead of a centrality measure. In general, though, local density minima do not always relate to branches, especially for data sets with multiple branching clusters.

In this use case, we show how FLASC can detect the branching pattern in this data set and classify the observations by their branch without manually extracting the flares from a visualisation.

Results

Figure 8 shows the detected subgroups encoded using colour on the 3D scatterplot. FLASC’s classification distinguishes the branches from the central core (Fig. 8A). The algorithm also finds a low-persistent flare representing the central core’s bottom. This flare could be ignored by specifying a persistence threshold. In contrast, HDBSCAN*’s classification does not find the branches (Fig. 8B). Instead, it finds part of the left branch as a small low-persistent cluster and merges most of the right branch with the central core.

This case study demonstrated how FLASC detects branch-based subgroups that do not contain local density maxima without having to specify the relevant features in advance or extract the subgroups visually. Practically, FLASC would have made it easier for researchers to detect the three groups in this data set, which was relevant for understanding diabetes and its causes.

Use case: cell development

Finally, we demonstrate a use case where detecting a branch hierarchy is important for understanding the data set. Specifically, we analyse a cell development atlas for the C. Elegans—a small roundworm often used in biological studies—by Packer et al. (2019). They analysed gene expressions in C. Elegans embryos to uncover the trajectories along which cells develop. Broadly speaking, this data set describes what happens in cells as they develop from a single egg cell into all the different tissues within fully grown C. Elegans worms.

After pre-processing, the data set appears to contain clusters and branching structures when viewed in a 3D projection (Packer et al., 2019). In this use case, we demonstrate that FLASC’s branch hierarchy provides interesting information about the data set’s structure even when the main subgroups can be detected as clusters.

Results

Figures 9A and 9B show 2D UMAP projections of the data. Data points are coloured to indicate the detected clusters and branches, respectively. There are two main differences between the branches and the clusters. Firstly, two regions where branches merge are detected as clusters, namely clusters 10 and 14. These regions do not have a distinct branch label but are identifiable in the branch hierarchy (Fig. 9D). Secondly, only branch 13 was not detected as a cluster in this dataset, indicating that the branches represent regions with high local density. Considering that the branches correspond to developmental end-states, it is unsurprising that local density maxima occur within them. One could imagine that the variation in gene expression is higher during development and that fully developed cells are observed more frequently. Both scenarios could cause these local density maxima.

More interesting are the differences between the cluster and branch hierarchies. Figures 9C and 9D visualise these hierarchies as an icicle plot. These designs are adapted from McInnes, Healy & Astels (2017) to indicate the selected branches and clusters using colour and a label. Segment widths encode the number of points in the tree below the segment. The hierarchies highlight that while HDBSCAN* detects the branches as clusters, it does not capture the trajectories. For example, the hierarchy does not reflect that clusters 0 and 1 connect to the whole data set through cluster 10. FLASC’s branch hierarchy, on the other hand, more closely resembles the data’s shape. For example, the hierarchy describes the embedding’s left and right sides with five and six branches, respectively. Generally, branches that merge into the cluster near each other are also close in the branch-condensed tree. The branch-condensed tree captures only the most eccentric connection for branches connected to multiple other branches in the cluster approximation graph.

This use case demonstrated that FLASC’s branch hierarchy provides information about the data’s shape that may not be obvious from cluster hierarchies. In addition, we found it beneficial to suppress noisy connectivity using dimensionality reduction techniques when detecting branches.

FLASC’s practical value

As demonstrated by the cell development use case, the argument that branches are not detectable as clusters only applies when they do not contain local density maxima. Subpopulations tend to have some location in feature space where observations are more likely. These locations are detectable as local density maxima, allowing data points surrounding them to be classified as a particular cluster. If one is only interested in the existence of subgroups, then FLASC only provides a benefit on datasets where relevant subgroups are sparse (e.g., the diabetes types use case). If one is also interested in the clusters’ shapes, then FLASC’s branch hierarchy provides information that cannot be extracted from a cluster hierarchy. We envision FLASC as a valuable tool for exploring unfamiliar data, providing guidance into which subpopulations exist and informing follow-up questions. Knowing that a cluster may represent multiple subgroups can be very relevant.

Alternative eccentricity metrics

The presented FLASC algorithm uses a geometric distance-to-centroid metric to describe how eccentric data points are within a cluster (Fig. 2A, Eq. (4)). An interesting alternative is an unweighted geodesic eccentricity, which measures path lengths between each data point and the cluster’s root point in the cluster approximation graph. Here, the root point can be chosen as the data point closest to the cluster’s centroid, as we did for the branch-membership vectors (Fig. 3). Such a geodesic eccentricity would agree with the notion that distances in high dimensional data may not accurately reflect distances along the intrinsic structure of a data set, which was one of the motivations for Reversed Graph Embeddings (Mao et al., 2017). It would also be closer to the maximum shortest-path centrality metric used by Vandaele et al. (2021).

Several trade-offs between the geometric and geodesic eccentricity metrics made us choose the geometric one:

• Computing the geodesic eccentricity is more expensive because it requires an additional traversal over the entire cluster approximation graph. The extra cost, however, should be low compared to other parts of the algorithm.

• The resolution of the unweighted geodesic eccentricity is lower, as it expresses the number of edges to the root point. As a result, zero-persistent branches are more likely to occur. In addition, it reduces the detectability of small branches that are well connected. On the other hand, that can be seen as beneficial noise suppression. In addition, the branch selection persistence parameter becomes more interpretable and would represent the traversal depth of a branch in the approximation graph.

• The cluster’s centroid may lie outside the cluster, resulting in a root point and eccentricity values that do not accurately describe its centre. For example, imagine a U-shaped cluster. The centroid would lie in between the two arms, and the root would lie in one of the arms. As a result, the geodesic metric would find one smaller and one larger branch rather than two equal branches. On the other hand, the geometric eccentricity finds the two branches and the connecting bend as three separate groups. Confusingly, it also contains two regions with a local eccentricity maximum, which FLASC gives a single label. Placing the root at an arbitrary eccentric location, as in Ge et al. (2011), avoids this issue but necessitates a different interpretation of the branch hierarchy and membership probability.

FLASC’s general process can also be used with measures that capture other aspects than eccentricity. At its core, FLASC consists of two filtrations, one to determine the connectivity between data points and one to analyse a signal on the resulting graph. The process would then describe how many distinct local minima (or maxima) of the measure exist within the clusters. The resulting interpretation does not have to relate to the cluster’s shape. For example, one could use the boundary coefficient (Vandaele, Saeys & De Bie, 2020) or direction centrality metric (Peng et al., 2022) to create an efficient and principled boundary-seeking clustering algorithm.

One could even interpret FLASC as two applications of HDBSCAN*: one over the density and one over the eccentricity. This perspective raises a possible improvement to the algorithm by translating the mutual reachability concept to the centrality metric. The idea of ‘pushing away points in low-density regions’ can also be applied to the centrality, emphasising the centrality difference between the centre and branch ends. Additionally, smoothing the centrality profile by incorporating neighbouring values could improve the algorithm’s robustness to noise. The additional computational cost should be low, as points’ neighbours are already known when the centrality is computed. Another way to improve noise robustness could be to implement the mutual $k$-nearest neighbour approach Dalmia & Sia (2021) used to improve cluster separation in UMAP projections. It would provide a subgraph of the core approximation graph that better reflects the cluster’s connectivity in high dimensional data sets. We leave evaluating these ideas for future work.

Visually summarising data’s shape

A strength of Mapper (Singh, Memoli & Carlsson, 2007) and Reversed Graph Embeddings (Mao et al., 2017) is that they can summarise the data’s shape using intuitive visualisations. While FLASC’s branch-condensed tree provides some information about the clusters’ shapes, interpreting the shape is not trivial. Studying how well two-dimensional layouts of FLASC’s cluster approximation graphs work as shape summarising visualisations would be an interesting future research direction. These graphs directly encode the connectivity used by the algorithm. Another benefit is that—unlike in Mapper—all (non-noise) data points are represented in the graphs once. Directly visualising the graphs, however, probably does not scale to larger sizes in terms of computational cost for the layout algorithm and visual interpretability. Ways to summarise the networks would have to be found, which could be based on kMeans centroids like in Reversed Graph Embeddings (Mao et al., 2017), local density maxima in the cluster, or a Reeb-Graph approach similar to Ge et al. (2011).