The ubiquitous deployment of electronic and computer based data collection systems has created massive data sets that defy human analysis. As a result, a broad collection of mechanized data analysis/data mining techniques have emerged. One approach for analyzing data is based on the mathematical field of topology and is called Topological Data Analysis (TDA) [1]- [5]. Topology is a branch of mathematics that characterizes the properties of a space that are preserved under continuous deformation [6]. TDA leverages this aspect of topology to extract knowledge based on the shape and form of the data. There are two main techniques that apply TDA techniques for data analysis, namely: Persistent Homology (PH) and mapper [1]. TDA, and specifically PH, has been Support for this work was provided in part by the National Science Foundation under grant IIS-1909096. demonstrated as an effective data mining technique for a number of fields. For example, PH has been applied to analyze networks [7]- [9], brain artery trees [10], digital images [11]- [16], protein structures [17]- [19] gene sequences [20]- [22], and cardiovascular diseases [23], [24] to name a few.
Computational persistent homology explores the shape of the data as the data in the point cloud is interconnected at different spatial distances (often called ǫ distances). The computation records the ǫ distances when topological features (holes, loops, and voids) appear (called the birth) and disappear (called the death). The < birth, death > pair defines the persistence interval of each topological feature and they can be displayed in a variety of different forms, including: barcodes, persistence diagrams, persistence landscapes, and persistence images [5]. In most cases, persistence intervals are paired with the dimension that they exist and so the persistence interval becomes < dimension >, < birth >, < death >.
Unfortunately, the computation of PH is exponential in both time and space [5]. This prevents the direct application of PH on big data. For example, the most efficient libraries for computing PH (currently Ripser [25], GUDHI [26], and Eirene [27]) can only process a few thousand data points in R 3 (unless strict limits are otherwise placed on the PH processing parameters. This limitation has motivated studies to explore the application of data reduction techniques to permit the application of PH on larger data sets [28]- [31]. These techniques provide good approximations of the PH of the large features in the point cloud and achieve the PH computation at 3-4 orders of magnitude faster than using the entire point cloud (when the entire point cloud can be analyzed). The data reduction technique of Chazal et al [28] uses repeated trials of random selections of data from the original point cloud; the technique of de Silva and Carlsson [30] uses either random or Maxmin (which finds dispersed points across the point cloud) to select "landmark" points that have other nearby members in the point cloud for use; and Moitra et al [29], [31] use the centroids of k-means++ clusters. While experiments show that k-means++ centroids provide good results, more accurate results might be possible with other topologically preserving data reduction methods.
This paper compares different methods to reduce large point cloud data sets for computing PH. The study explores multiple scalings of data reduction using several different data reduction methods and evaluates how well each method preserves the topological features in the point cloud. The key objective is to discover how different methods of data reduction sample the data such that the resulting sample P ′ is a topologically faithful representation of the original point cloud P. By definition, the sample is not able to maintain all of the topological features of the original. However, many uses of TDA are interested only in the presence and shape of the larger (longer lived) topological features while considering the small (short lived) topological features insignificant. Thus, the specific goal of data reduction is to maintain the structure of the large topological features while removing the less concerning small topological features. This paper characterizes the different types of losses that can occur due to data reduction.
In addition, an experimental assessment that expands the study initiated by Moitra et al [29] with additional clustering algorithms and the random and Maxmin methods of [28], [30] is presented. The samples are composed of the cluster centroids from several different clustering algorithms, namely: k-means++ [32], Agglomerative single-linkage [33], and Agglomerative (Ward's) [34]. In addition comparison to random sampling (as performed in [28]) and the Maxmin algorithm developed to select landmark points by de Silva and Carlsson [30] are included in this study. The PH result from each sampling method is compared to the PH result from the full data set where the size of the data permits. While the theory shows many possible losses from sampling, the experimental results show that in practice the large topological features are well-preserved by many of the clustering and random sampling data reductions. However, as the data reductions increase in scale, the PH results are better preserved by only a few of the data clustering algorithms.
The remainder of this paper is organized as follows. Section II presents some of the background on PH and topology preservation. Section III briefly describes related work. Section IV outlines the general strategy for data reduction with cluster centroids. Section V provides the motivation for data reduction using different clustering algorithms. Section VI presents the experimental results on several different data sets. Finally, we conclude the paper with some remarks in Section VII.
This section contains a brief overview of persistent homology and highlights some of the formal theories on the preservation of topological structure by various data clustering algorithms [35]. A detailed overview is available at [1], [5].
Homology is a means to characterize the features of a topological space. Homology groups at different dimensions d represent the topological features in that dimension. The Betti number at each dimension d is the rank of the d th homology group, which corresponds to the number of holes in d [1].
Computing homology on a finite metric space is achieved by approximating the space with a representative complex formed from the points in that space. While there are several types of complexes (cubical, CW, or simplicial [5]) simplicial complexes are the most widely used and they will be used in this paper. Simplicial complexes are composed of simplicies; generalizations of a triangle to any number of dimensions. For example, a 0-simplex is a point, an edge is a 1-simplex, a triangle is a 2-simplex, tetrahedron a 3-simplex, and so on. However, examining a point cloud statically does not yield any meaningful topological information [5]. Instead, the point cloud must be observed at multiple scales via geometric filtrations in order to discern topological features. Features that exist over multiple scales are "persistent" and therefore meaningful from a topology standpoint. Persistent Homology (PH) is then the notion of homology applied to multiple geometric scales on a finite metric space.
One of the most common and computationally feasible complexes to construct filtrations is the Vietoris-Rips (VR) complex [36]. The filtration is a set of subcomplexes with each subcomplex based on a distance ǫ i . The ǫ distances begin at 0 and increases until all points in the point cloud are connected (although in practice a maximum ǫ distance is often defined). As ǫ i increases, topological features (holes, voids, and loops) appear and disappear in the space. Topologically speaking, the ǫ distance that a topological feature first appears in the filtration is called the birth. The death of a topological feature occurs at the ǫ distance where that topological feature no longer exists in the filtration. Each < birth, death > tuple for a topological feature is called a persistence interval. Typically longer persistence intervals represent topologically meaningful features of the data and shorter ones represent noise [1]. However, in some cases the shorter persistence intervals are also of interest [10], [37].
A theoretical framework for the usage of clustering algorithms to preserve topology of a space is developed by Carlsson and Memoli [35]. Preliminary work by Niyogi et al. [38] shows that homological information of a manifold can be inferred from a random sampling of points distributed around it. In this case the manifold M is characterized as a low-dimensional underlying geometric space of the points in question. Assumptions are made that all probability distributions of points around the manifold are supported by it and any noise points are distributed via Gaussians. Given these strict assumptions Niyogi et al. [38] were able to show it is possible to infer higher order homological information from a sampling and reconstruction of connected components. In other words, low dimensional topological features can be used to infer higher dimensional features. In practice it is nearly impossible to satisfy the assumptions of Niyogi et al; especially for experimental data [39]. In order to obtain topological information from noisy data that does not lie around a manifold persistence must be incorporated with any homology analysis. Persistence allows one to determine which homology groups of a space are not created by noise. Persistence also identifies which features continue existing as the dimension of the homology groups increases. This solves the issue of not being able to rely on an underlying manifold for a data set. Carlsson extends the work of [38] by exploring what occurs when a more robust sampling method beyond random sampling is used [39]. Carlsson and Memoli later found that clustering algorithms in particular enable persistent homology to be computed on reduced data sets with provable preservation of persistent features [35].
As outlined in [35], [39] a clustering algorithm C will take an input set of points and create a mapping which separates the points into a set of output points in partitions P. Carlsson notes that clustering is a statistical method of sampling and mapping the connected components of a topological space to a partitioned space [39]. This mapping is unique because it acts as a function between two disparate mathematical areas -topology and set theory. This allows Carlsson to use category theory to characterize how clustering algorithms change topological spaces. A brief description of category theory is given below; see [40], [41] for a detailed background.
Category theory allows for any entity that fulfills certain conditions to be meaningfully compared to another entity that fulfills the same conditions. Generally these entities are mathematical in nature, however, they do not have to be. These entities are known as categories. For something to be referred to as a category, it must be a set of objects with morphisms (mappings/transformations) between pairs of objects and possess composition and identity properties. In this case, the original data set and the reduced data set are both categories and the clustering algorithm C is a functor between them. By definition, functors preserve the composition and identity properties of the categories they map from. This mapping is continuous. This usage of clustering algorithms as functors allowed Carlsson and Memoli [35] to develop theories about how they preserve topological structures of a space and persistence of topological features in the mapping they create.
The computation of PH via the VR complex becomes intractable as the size of a small point cloud extends beyond a few thousand points in R 3 . A number of efforts have been made to simplify the computation wherever possible. The primary bottlenecks of computing PH are the size of the simplicial complex and the size and reduction of the boundary matrix [5]. Early work to alleviate these bottlenecks was achieved by modifications to the simplicial complex, primarily concerned with sparsifying the complex. In particular, Sheehy developed two theoretical methods to sparsify the Vietoris Rips complex using net trees [42]. Sheehy shows that the use of net trees to remove points and their incident simplices from the VR complex does not change the topology. However, this sparsification did not scale well to larger point clouds [43]. To attack this problem, Dey et al. [43] implements a method to approximate VR filtrations on much larger point clouds than Sheehy through batch collapse of simplices. Dey's usage of batch and cluster set distances as opposed to vertex distances when merging complexes resulted in increased scalability from the initial approach described in [44]. More recently, Brehm and Hardering [45] were able to implement a more scalable version of Sheehy's sparsification method in the Julia library Sparips. Sparips first builds a contraction tree (similar to a cover tree [46]) over the raw data before the construction of the boundary matrix. The boundary matrix is then sparsified using information from the contraction tree. Sparips provides a tighter bound on approximation of PH than [43] and obtains performance comparable to GUDHI. All of these approaches reduce the number of complexes constructed from the original point cloud in order to identify significant topological features on larger point clouds.
Data reduction through sampling has also been explored to expand the processing capabilities of PH libraries. Chazal et al. [28] utilize repeated random sampling of the original point cloud and compute PH from the average landscape of those samples. However, as noted by Sheehy [42], the combination of different persistent diagrams into a single diagram did not yield the same accuracy; strict assumptions had to be made on the data for the best results. Moitra et al. [29] used a similar approach that sampled the point cloud using k-means++. Because the k-means++ algorithm samples the data through multiple iterations to reduce the WCSS error, PH only needs to be calculated on one sample, as opposed to the multiple samples in Chazal et al. Moitra et al. shows that sampling the data using k-means++ preserves significant topological features and has similar results to Chazal et al [28]. In addition, k-means++ also allows for the use of upscaling of the reduced data to increase the accuracy of the persistent intervals computed from the sampled point cloud [31]. However, these studies did not explore the impact of data reduction at increasingly large reduction percentages or for the broader impact that various other data reduction methods would have on the results of computing PH on reduced data.
Finally, some studies have attempted to use random projection to enable the computation of PH on high dimensional data sets. Random projection allows higher dimensional data to be mapped to lower dimensions, while preserving distances between points with bounded error [47]. Sheehy uses this idea to prove that random projection preserves the persistent homology of the point cloud to a comparable bound [48]. From the theoretical results, Ramamurthy [49] conducted experiments on PH of randomly projected point clouds and showed that the persistence diagrams were similar for a variety of random projections.
The focus of this paper is the transformation of a point cloud P that is too large for computing PH to another point cloud P ′ with fewer total points such that the PH can be computed on a representative point cloud. Ideally, the transformation should be such that P and P ′ maintain the structure of large topological structures, homeomorphic to some degree. An example of a suitable mapping is illustrated in Figure 1. The leftmost image has 2, 000 points and represents the original point cloud P; the center and rightmost images represent two topologically preserving reductions of 500 and 250 points and are possible representations of the reduced point cloud P ′ . The general approach is to partition the original point
• r i , the distance from the partition centroid, P ′ i ∈ P ′ , to the most distant point in that partition, and • r max = max(r i ), the maximum r i of all the partitions. Using the topologically similar (but smaller) point could P ′ to approximate the PH of P will identify the large topological features (as defined by the bounds of Section IV-B). In particular, let B be the boundary of points in the complex defining a d ≥ 2-dimensional topological feature and let s B = max(distance(b i , b j )) ∀b i , b j ∈ B. Then define the term "large topological feature" to be any feature with diameter s B > 2r max . That is, a large topological feature has a diameter that is not contained within the largest partition of P. Depending on its location in the partitions, any topological feature with a diameter smaller than 2r max may or may not be identified during this step. In particular, any topological feature that falls within the boundary of a partition will be lost; any topological feature that extends beyond the boundaries of the partitions are likely to be retained. The degree to which features are lost may be significantly impacted by the mechanism/algorithm used to define the partitions.
The partitioning step is used to create a mapping from the original, large input point cloud P to a topologically similar but smaller (in terms of total points) point cloud P ′ . Ideally this mapping will be performed in a manner that preserves the larger topological features of the original point cloud. Formally ||P ′ || < M < ||P|| so that PH can be computed on P ′ in the allotted time and space. This step is performed as follows. The PH of the original point cloud is estimated by a computation of PH on the partition centroids P ′ . When computing PH on P ′ instead of P, the smaller topological features that lie within a radius of any partition will be hidden from this PH computation. For the large topological features in P, approximating the PH using P ′ can result in the following deviations from the PH results of P:
1) The persistence interval < birth, death > may occur at different (but bounded) ǫ distances. 2) A false topological feature not actually in P may be identified. 3) A topological feature of P in dimension n > 2 may present itself in P ′ at a different dimension m. 4) A topological feature might be lost by the data reduction step.
The frequency and significance of these deviations is influenced by the partitioning methods used to define P and consequently P ′ . For the remainder of this section, a worst case characterizations of these deviations will be presented.
As developed by Moitra et al [29].
Theorem 1. The shift in the < birth, death > persistence interval values arising when computing PH from P ′ instead of P is bounded by 2r max .
Fortunately, this error can be reduced by an upscaling step [31]. More precisely, the approximate < birth, death > interval for any topological feature identified from the estimated PH can be refined by recomputing the PH using all of (and only) the points from the partitions containing the centroids that form the boundary of the feature. This process is called upscaling. If the upscaled point cloud contains too many points for computing PH, an iterative repartitioning and upscaling can be performed to refine the approximation of the feature boundary. Of course there is a limit; if the convex hull of points on the boundary exceeds M , then the improvements by upscaling may be further limited.
False topological features can arise due to false voids that lie between the centroids of the reduced point cloud P ′ . That is, the gaps between the centroids due to the removal of the partition points can be such that the centroids of P ′ define a complex around a false topological feature. Formally, these false features can occur only when: Theorem 2 (False voids from centroid gaps). False voids can appear when ǫ min < r max < ǫ max .
Proof. Consider a 2-dimensional space of 4 square partitions with radius r uniformly filled with points in P, where ǫ min <r< ǫ max and where the minimum pairwise distance between any two points in P is less than ǫ min . Then P ′ would consist of the centroids of these squares. Computing PH on P ′ would result in the discovery of a topological feature not present in the original point cloud P.
False topological features can be pruned by upscaling [31].
Estimating the PH of P using P ′ can also cause a topological feature to shift dimensions. Shifts into higher dimensions occur when the points in P ′ stretch a feature with a void space in new dimension. Shifts into lower dimensions occur when a topological feature that has a convex hull in k-dimensions in P loses a cover in one (or more) of the dimensions in P ′ so that the convex hull only occurs in j < k-dimensions. That is: Theorem 3 (Feature shift to higher dimensions). The PH computation in P ′ may shift an identified topological feature into a higher dimension.
Proof. Consider a point cloud in R 3 composed of N points that contains a single 2-dimensional circle of radius r > 2ǫ max in the xy-plane at z = j. Consider a partitioning of the space such that one partition is a square with sides of length r and located immediately above the 2-dimensional circle and with all points outside of the square defining a partition. The data reduction step will introduce a void space in the z dimension above and including the circle.
Feature shifts to higher dimensions are not a significant issue. They will be pushed back into the proper dimension with the upscaling computation.
Theorem 4 (Feature shift to lower dimensions). The PH computation in P ′ may shift an identified topological feature into a lower dimension.
Proof. Consider a 3-dimensional space with a sphere at the origin with radius r s (ǫ min < r s < ǫ max ) and with a uniform distribution of points extending some finite distance beyond the surface of the sphere. Consider a partitioning such that (a) a collection of partitions that lie in the xy-plane at z = 0 with the maximum radius r max for all of these partitions is such that ǫ min < r max < ǫ max the remaining space lies in two partitions: one covering all points in P at (x, y, z > 0) and with a centroid at (x, y, z > ǫ max ); and covering all points in P at (x, y, z < 0) and with with a centroid at (x, y, z < ǫ max ). Then the resulting centroid points would be such that the sphere would only appear to the PH algorithm as a 2-dimensional circle and no connections would be made to the centroids above and below the z = 0 axis as they would lie outside the range of the ǫ values examined by the PH algorithm.
Feature shifts into lower dimensions are more problematic. They cannot easily be restored to their proper dimension and the feature will be lost.
In rare cases, a topological feature can be lost when the points defining the boundary for the convex hull in P surrounding the topological feature are insufficient in P ′ for a corresponding convex hull to be defined. This issue may motivate the use of multiple partitioning steps to estimate PH with strategic partitioning methods.
The motivation of this work is to be able to compute PH on big data. Preliminary data [29], [31] suggests that spherical clustering methods such as k-means++ can present a suitable partitioning of the data. However, it is unclear if this is the best method for data reduction or if other algorithms or heuristics provide better topologically preserving abstractions.
This experimental study examines the use of other clustering methods for data reduction. In particular hierarchical clustering algorithms are of interest due to their theoretical guarantees of preserving persistent features [35]. More precisely, a comparative analysis of k-means++, single-linkage agglomerative clustering, ward-linkage agglomerative clustering, and Silva's Witness maxmin sampling [30] is performed.
Density-based cluster algorithms such as DBSCAN [50], HDBSCAN [51], and mean-shift were evaluated but performed poorly due to the lack of ability to set k, the number of centroids to generate. This led to classifications at different bandwidths and parameters to attempt to approximate the point clouds of suitable size. That said, several suitably sized point clouds were derived using these methods. However, ultimately the reduced point clouds P ′ extracted using these densitybased clustering algorithms failed to preserve the topological features on par with other methods. As a result these algorithms are excluded from the results in this section.
The experimental study measures the accuracy of the persistence intervals computed with the various approaches against a reference computation of the persistence intervals on the original point cloud. Accuracy at this stage is determined by the Heat Kernel Distance (HKD) [52], [53] between the original persistence intervals computed from P and the persistence intervals computed from P ′ . The HKD gives a stable heatkernel metric for classification applications and is a way to obtain topological inferences about an object using Gaussian kernel density estimates. Additionally heat kernel distance is robust to noise and outliers in data, which makes it a valuable analysis metric for comparing persistence intervals.
Each algorithm under comparison was implemented as a preprocessing step to reduce the input point cloud. In the case of clustering algorithms such as k-means++, the cluster centroids were output; for hierarchical cases the dendrogram was cut at a specific threshold. After each reduced point cloud was obtained, the PH library Ripser [25] was used to compute PH on the original and reduced point clouds. Ripser is currently a state-of-the-art library for computing PH with Vietoris-Rips complexes, efficient in both speed and memory performance. However, similar results can be obtained with GUDHI [26] or Eirene [27]. All experiments were conducted using an Ryzen Threadripper 1950X with 128GB of RAM.
The HKD comparisons are computed on persistence intervals separated by dimension. While an aggregate comparison of the HKD to the complete set of persistence intervals (independent of dimension) was performed the results did not provide any significant insights and, due to space considerations, they are not presented here. Accuracy measured by persistence interval dimension gives a notion of the preservation of all persistence intervals in the original data set. However, many of the persistence intervals contributing to the HKD results in this instance can be attributed to shorter persistence intervals or noise that are by definition going to be missing in the reduced point cloud. Since the focus on this paper is to how well these reduction methods preserve the significant, larger topological features, a deeper analysis is necessary. As a result, several methods to filter the persistence intervals by length are explored and compared. That is, the persistence intervals from the original point cloud are examined and cutoff lengths were established using several filtering methods in an attempt to isolate the significant persistence intervals for comparison. For example, one filter finds, by dimension, the shortest persistence interval from the longest 10% found in the original point cloud. This length is then used to filter persistence intervals computed from the original and reduced data sets. The HKD for these filtered results are then computed and reported. This will provide insight on how well the reduced data sets preserve the significant, larger topological features. In cases where the persistent homology of the original data set could not be computed due to resource constraints, the maximum number and mechanized approach to filtering the top-most intervals.
Unfortunately with filtering the top 10% of persistence intervals, the resultant intervals may include shorter intervals based on the number of persistence intervals generated. All three filters refine the compared persistence intervals provides a measurement of preservation of the salient topological features in the point cloud. An analysis of all data sets chosen for VI-A are included in Figure II to present a further comparison of k-means++ and Agglomerative Ward and their preservation of these longer persistence intervals at various levels of reduction.
There are several notable findings in the filtered comparison as shorter barcodes contribute less noise to the heat kernel distance metric. In the water treatment dataset, the H 1 features in all three filterings displays a significant loss of identified features when reducing from 200 to 100 points. This dropoff indicates a large feature has been lost during the reduction and may indicate the reduction's limits for that specific dataset. Utilizing the persistence interval filtering can provide insight to the accuracy of each algorithm in extreme cases of reduction to quantify only the large persistence intervals for comparison.
Overall partitioning algorithms perform with similar accuracy up to a significant percentage of data reduction. However, several of the algorithms better preserve the salient topological features at more significant reduction percentages, namely k-means++ and aggWard. These algorithms provide similar accuracy results and can give the most accurate persistence intervals under significant reduction. Between k-means++ and aggWard, there are slight differences in their ability to preserve long persistence intervals depending on the structure of the point cloud. Both are nearly identical in the twoCircles and twoMoons dataset; aggWard seems to perform better with triangulated mesh point clouds, while k-means++ has higher accuracy with the categorical point clouds seeds and water. VII. CONCLUSIONS Persistent homology presents a novel approach to analyzing data. Unfortunately the computation of persistent homology on big data is not currently possible due to its exponential complexity. Data reduction is a classic approximation technique used by the data mining/machine learning communities to attack computational complexity issues and various explorations to use data reduction for the computation of persistent homology have been performed. This paper has explored how well different data reduction strategies functioned to preserve the large topological features present in a point cloud.
While the experiments performed in this analysis focused on several notable clustering algorithms, only slight improvements in some algorithms were identified to provide better partitioning results. Algorithm complexity can play a large role in determining what partitioning algorithm is suitable to an application. In general, k-means++ and Agglomerative Ward tend to consistently provide better overall data reduction results. While Agglomerative Single Link provides comparable results, at larger reductions it tended to lose more persistence intervals in the higher dimensions than k-means++ and Agglomerative Ward. Random sampling and the Maxmin [30] sampling sometimes provide good results, but they tended to be unpredictable and sometimes produce wildly inaccurate results. An additional benefit from the use of clustering is that data around a sampled point can be restored to support the concept of upscaling [31] to restore more accurate persistence intervals from the reduced point cloud.
Big data continues to be on the horizon for persistent homology as techniques for data reduction, simplicial complex collapses and optimizations, and boundary matrix reduction continue to increase performance. Approximations in the hundreds of thousands of points should be possible with a well designed partitioning and upscaling library, even those in higher dimensions. Bringing TDA to big data analysis will provide automated tools for analyzing the connectivity of point clouds beyond current applications and should be continued to be explored in all domains.