Persistent homology for high-dimensional data based on spectral methods

Dear reviewer BdHd,

thank you for you detailed review. We are happy that you found the problem setting we tackle "significant", our approach "novel", and our analysis "comprehensive". We will address your concerns in the following.

Outlier sensitivity:
This is an interesting question, but we disagree about outliers being a problem for kNN-graph-related methods. First, all methods but Euclidean distance and Fermat distances, in particular DTM-based filtrations which are designed to handle outliers, rely on the kNN graph in some way. However, they do not use the raw distances as edge weights, but aggregate these in a way that avoids sensitivity to outliers. As an illustrative example, consider effective resistance and the 2-NN graph. In case of single outlier $o$, it would have its nearest neighbors $p_1$ and $p_2$ on the main data structure, say the noisy circle. These neighbors might be far away in the Euclidean distance (which is why Euclidean distance can suffer from outliers). But the 2-NN graph on the circle is densely connected, such that the effective resistance distance between $p_1$ and $p_2$ is small. Thus, the triangle $op_1p_2$ is very slim and does not lead to a highly persistent 1D hole.

To validate this, we now conducted additional experiments with uniformly distributed outliers in the bounding box around the noised circle (Appendix G). In low ambient dimension, outliers hurt the Euclidean distance at low noise levels a lot. Effective resistance remained much more robust against outliers.

Second, as we describe in Section 4, outliers cease to be a big problem in high-dimensional spaces, as they are necessarily extremely sparse. It is highly unlikely that an outlier point falls, say, into the middle of circle in high ambient dimensionality. We backed this up by additional experiments, showing that even the Euclidean distance is robust to adding up to 10% uniformly distributed outliers in $\mathbb R^{50}$. Effective resistance continues to outperform the Euclidean distance even in the presence of outliers in this setting.

Usefulness of our detection score:
We completely agree that our detection score is only useful if ground truth on the number of holes in the data is available. This is akin to most metrics in machine learning, which only work in reference to a ground truth. If ground truth is available, our metric works for any number of holes in the data as it targets the gap between the $m$-th and $(m+1)$-st highest persistence for $m$ ground truth holes. For the synthetic datasets ground truth is known by construction. For the single-cell datasets, we made an effort to only benchmark datasets for which the number of loops is clearly determined by the biology of the sample. For the malaria dataset there are two cycles: the parasite reproduction cycle in red blood cells and the parasite transmission cycle between humans and mosquitos as described in Section 7.2. The other five datasets are created so that the gene expression variation comes from the cell division cycle (mytosis), that is, the process of a cell dividing into two copies of itself. Therefore, we expect one loop in these datasets.

Downstream single-cell experiment:
We are happy that your appreciate our diverse set of experiments! An experiment like the one you suggested is indeed a natural next step. Thanks for making this suggestion. However, this type of experiments is future work and could not be conducted during the rebuttal period (and would probably be of more interest for biologists, as opposed to ML researchers). In addition to using an effective form of persistent homology in a classification pipeline, we also believe that looking for holes in complex single-cell datasets would be an interesting and important research direction. But in order to trust persistent homology as a tool for exploring the structure of single-cell data, we first had to validate which approach to PH is trustworthy. This is the aim of our paper.

Restructuring suggestions:
Thank you for these constructive suggestions! The section on the curse of dimensionality for persistent homology is a new contribution and not background. Since the flow from persistent homology (Section 3) to its curse of dimensionatlity (Section 4) and finally the resolution with spectral methods (Section 5) seems fitting, we kept these chapters as they were. We did, however, move the more theoretical section on the relation between spectral distances up and combined the two experimental sections, as you suggested.

If you have are any further concerns, please do let us know and we will be happy to answer them!

Dear reviewer qoEF,

thank you for writing your review! We appreciate that you found our description of the curse of dimensionality for persistent homology "very clear and easy to understand". We address your concerns below:

Novelty:
Please confer the general comment on novelty. Briefly put, we believe that solving a new, relevant problem with a novel combination of existing methods can be a more valuable contribution than inventing a new method specific to the problem. Please note that the corrected effective resistance had NEVER been combined with persistent homology before (as we now explicitly state in the related work section). Also, all these methods have never been applied to high-dimensional data with high-dimensional Gaussian noise, which is a setting very important in many practical applications. Finally, our theoretical analysis of corrected effective resistance and its relationship to diffusion distances is also novel.

Explicit formula for corrected effective resistance:
Our explicit formula enabled us to relate corrected effective resistance to diffusion distances. Using it, we found, for instance, that the corrected effective resistance decays eigenvectors with high eigenvalues more aggressively than the uncorrected version, but less aggressively than diffusion distances. We consider it an important theoretical result of our paper.

Theoretical treatment of curse of dimensionality:
This is a very interesting, but difficult question, and is beyond the scope of this rebuttal. Our Section 4 provides a qualitative treatment of this curse of dimensionality (which we have not seen explicitly discussed in the literature before). A quantitative treatment remains for future work. In fact, motivated by your question, we started working on it and have some promising ideas for a theoretical treatment, but would need more time to work on it than the rebuttal period allows.

If you have are any further concerns, we will be happy to answer them!

Dear reviewer k35W,
thank you for reading our paper and writing the review! Please note many of the details you asked about are all present in the manuscript. Our overall approach is conceptionally standard and close to previous works such as Bendrich et al. 2011, Anai et al. 2020, Fermandez et al. 2023: We compute a pairwise distance matrix, build the filtered Vietoris-Rips complex based on that distance matrix, and compute persistent homology of this filtered simplicial complex. Computationally, this was done using the ripser library. We only needed to compute various distances matrices and pass them to ripser. We have described this approach in detail in Section 3.

That said, we will make an effort to resolve any confusion in our replies below!

But it is not clear what does persistences p_m of the m-th most persistent features mean?

We describe in Section 3 that the persistence of a feature is the difference between its death and birth time: "Each hole has associated birth and death times $(\tau_b, \tau_d)$, i.e., the first and last filtration value $\tau$ at which that hole exists. Their difference $p = \tau_d - \tau_b$ is called the persistence or life time of the hole and quantifies its prominence."

The m-th most persistent feature is that with the m-th highest persistence.

How are these persistent features related to the loops/holes?

We describe persistent homology in Section 3. A $\delta$-dimensional feature that persists from $\tau_b$ to $\tau_d$ corresponds to $\delta$-dimensional hole in the simplicial complex between filtration values $\tau_b$ and $\tau_d$ which in turn can be thought of $\delta$-dimensional holes that persist in a union of balls around the data points of radius $\tau_b$ to $\tau_d$. We stress that we explain this in Sec.3.

Given a distance metric, how are these features and the detection score computed?

We use the package ripser to compute the persistence diagram, as stated in Section 3 and Appendix F. This is a standard package for the computation of persistent homology. Given a pairwise distance matrix, ripser computes the persistent homology of the associated Vietoris-Rips complex. The output, a persistence diagram, contains the information about detected holes / loops. When we talk of the persistence of a feature, say, for effective resistance, we mean a feature that was computed by the ripser package when providing it with the effective resistance distance as input.
We describe the computation of our detection score in Section 7 (paragraph "Performance metric"): If the ground truth number of holes of dimension $\delta$ is $m$, we take the persistences $p_m$ and $p_{m+1}$ of the $m$-th and the $(m+1)$-st most persistent feature in dimension $\delta$ and compute the ratio $s_m=(p_m -p_{m+1}) / p_m$ as our detection score. In particular, all these details of our approach are present in the manuscript.

How are the holes/loops detected and what is the persistence (birth-death of holes) with respect to. Is the resolution scale with respect to the different distance metrics considered?

Holes are detected with the ripser package (Section 3). Indeed, the resolution scale is with respect to the different distance metrics considered. This is explained in Section 3, where we write that the ball-growing procedure is made computationally tractable by the use of a filtered simplicial complex, that we only work with the Vietoris-Rips complex, and that for this complex it suffices to specify a pairwise distance metric. Later in Section 7, paragraph "Distance measures" we write that we consider twelve distances as input to persistent homology.

If so, what is the role of the K-NN graph?

All distances other than Euclidean and Fermat depend on the kNN graph in various ways. We described the all distances briefly in Section 7 and in detail in Appendix C. Note that the input is always a point cloud, never a predefined graph. If a distance relies on the kNN graph, we compute it from the point cloud (usually with a parallelized brute-force approach implemented with the package pykeops). The exact role of the kNN graph depends on the distance. For instance, the distance we call Geodesics is the shortest path distance on the kNN graph.

The graph connection is predefined to find these distances in this case.

We do not quite understand what you mean by this sentence. Could you please elaborate?

How does the proposed method differ from them is not clearly described.

Overall, our paper is about evaluating which (possibly existing) methods perform well for detecting holes in the presence of high-dimensional noise, a problem setting that has not been addressed as we state in the first paragraph of the related work section. We added a sentence to the related work section clarifying that persistent homology had so far never been combined with the corrected effective resistance of von Luxburg et al. 2010a, one of the methods that we find performed best in our problem setting.

t-SNE and UMAP seem to perform better at detection holes and these method should also have lower computational cost.

Please see the general comment on the performance of t-SNE and UMAP. t-SNE and UMAP are typically slower that effective resistance or diffusion distance. Both first compute the kNN graph of the data, but t-SNE and UMAP then perform an iterative optimization scheme, while effective resistance only needs to compute the pseudoinverse of a matrix and diffusion distances only need to power a matrix.

Furthermore, t-SNE and UMAP are typically used for 2D embeddings, and it is not possible to detect higher-order homologies (e.g. H2 homologies such as voids) using a 2D embedding.

Here cycles/holes only seem to consider edges, and not higher order simplices. Is this correct?

No, this is not correct. We always consider simplices of the dimension needed to compute the type of holes (loops, voids) that we are interested in. For instance, if we compute 2D holes (voids), we consider simplices of dimensions 1, 2, and 3. This is handled by the ripser package.

How is the Vetoris-Rips complex constructed?

The Vietoris-Rips complex is a standard construction in topology. We describe it in Section 3: "... the Vietoris–Rips complex, which includes an n-simplex (v0, v1, ... , vn) at filtration time $\tau$ if the distance between all pairs vi, vj is at most $\tau$. Therefore, to build a Vietoris–Rips complex, it suffices to provide pairwise distances between all pairs of points."

In particular, the Vietoris-Rips complex on $N$ points can contain simplices of dimension up to $N-1$. Whether a particular simplex, even a high-dimensional one, is contained at time $\tau$ only depends on pairwise distances. The construction of the Vietoris-Rips complex is handled by the risper package. One of its arguments is the maximal feature dimension that should be computed. Based on this, it only creates those simplices necessary, so that in practice very high-dimensional simplices are not constructed.

Given n points, the complex can have large number of higher order simplices, and detecting high dimensional holes is very expensive.

As mentioned above, this part is abstracted away by the ripser package. However, since we always set a maximal dimension for the topological features that should be computed (2 for void detection and 1 for loop detection), no simplices of dimension 4 or higher are needed. Indeed, we observe that then we compute persistent homology of 2D holes (voids) the run time increases a lot. We discussed this in Section 8: "One limitation of persistent homology is its computational complexity. ...." We also reported exemplary run times in Table S3, again showing that void detection is much more costly than loop detection.

I thank the authors for their thorough responses to all reviewers' comments. These are very helpful, and make many of the details related to the methodology now more clear.

However, I think the current version of the paper falls short due to the following reasons, and I encourage the authors to consider these points when revising the draft:

(a) It appears the readers of the paper are expected to know what the ripser library package does, and should be familiar with the "standard techniques" presented in the existing literature listed in this response to understand and follow the paper. Many of the details in the responses to the reviewers are clearly missing.

(b) In the response above, authors say "We have described this approach in detail in Section 3." However, I do not see many of the details that are discussed here in the response, in section 3 (which is half a page) of the paper. In the section, there is just a single sentence, "We compute persistent homology via the ripser package", and the readers are expected to know what this package does.

(c) Also, in section 3 and the response above, it is suggested, "we use the Vietoris–Rips complex, which includes an n-simplex ..." As noted by the authors above, given $N$ points, a Vietoris–Rips complex can include simplices of order up to $N-1$. However, it appears the full Vietoris–Rips complex (with simplices of all orders) is never constructed, and only a simplicial complex with simplices of order up to $k=2$ or $3$ is considered. This makes sense (due to the cost), but was never clear in the paper. Similarly, the paper just says that persistent homology involves computing the homology groups of the union of all balls in order to find the holes. The dimension of these holes are not specified (mention of $\delta$-dimensional feature and holes only appear here in the response). Note that the general definition of persistent homology is constructing the full simplicial complex and finding/counting holes of all dimensions that persist across different scales $t$. The authors' response above suggests in the paper $\delta =1$ or $2$ are considered. Again, this is perfectly fine (since higher order holes are expensive to estimate), but is never discussed in the paper.

Since all these details are currently missing in the paper, a reader would just assume the full complex and holes of all dimensions are computed. Since constructing the full Vietoris–Rips complex and computing all Betti numbers will be very expensive (could be of exponential cost), a natural conclusion is to assume that the overall cost of the proposed method will be extremely high.

(d) The role of k-NN graph as related to persistent feature computation is still not clear to me, and section 7 does not contain much details related to this as suggested in the authors' response above.

(e) In the response above, authors say "our paper is about evaluating which (possibly existing) methods perform well for detecting holes in the presence of high-dimensional noise". This should have been the main premise of the presentation, but in the abstract and introduction, this fact is never made clear. Therefore, naturally reviewers have asked about novelty and relation to existing methods.

Overall, the current version of the paper is very difficult to follow and the main contributions, methodology, and significance are all not really clear, as also suggested by other reviewers' comments. Hence, I am keeping my current score.

I thank the authors for the response and the changes made to the paper.

But, I think the authors are missing the main point reviewers are trying to make, that the current paper presentation does not clearly convey different aspects (e.g., main contributions/premise, methodology used, and significance) of the presented work.

One of the key duties of reviewers is to provide suggestions to help improve the paper.

Mentioning something in some part of the paper, is not equivalent to a clear explanation or discussion, right? Readers should not be excepted to guess, or go back and forth/search for basic key details. For example, as the authors above say, the fact that only 1D and 2D holes are computed, is just mentioned (and one had to search for it) as a pass-by sentence, but in the main methodology section, this was not discussed. Same with the k-NN graph, indeed it is mentioned at several places, but the nice discussion presented in the responses to reviewers (that a $k$-NN graph is constructed from point cloud using certain package, then the spectral distance methods are applied to this graph, etc) was missing, and its role was not clear to the readers.

Regarding the ripser library, I agree a detailed explanation of the package is certainly not necessary. But, since the study heavily depends on this package, at least a few sentences on what does it take as inputs and arguments (e.g., the maximal feature dimension), how are the resolution scales defined for persistence (like how many levels does it consider), what does it output for persistent homology (is it a list of holes, how do we get m-th persistent features, etc) should be provided. Right now, we just have a single sentence about this. Without these basic details, it is very difficult to understand the methodology used, right?

Same with the main premise of the paper. The two sentences picked by the authors above might allude that the paper does some kind of comparison. However, the paper is not set up as "a comparison study of existing methods for detecting holes in the presence of high-dimensional noise". There are other contradicting sentences which imply differently. For examples, in the intro, it says "we suggest to use persistent homology with distances based on the spectral decomposition of the kNN graph Laplacian, such as the effective resistance ....", which can be interpreted as a new proposal which was not explored before. In the abstract, the term "As a remedy, .." (which is conveniently skipped above) suggests a new solution is proposed. I feel some of the comments by the reviewers are rooted in this fact, that this main premise (which is well explained in the general comment/response by the authors) was not really clear in the paper.

(Q23) Multiparameter persistence:
Multiparameter persistence is an interesting method and might offer benefits in the high-dimensional setting. It has been employed to spatial transcriptomics data in Benjamin, Katherine, et al. "Multiscale topology classifies and quantifies cell types in subcellular spatial transcriptomics." arXiv preprint arXiv:2212.06505 (2022).

However, its output is not as readily interpretable as the persistence diagrams in single parameter persistent homology, which is why we limited out study to this more established method.

Minor 1. Miscellaneous:
We added the feature dimension to all persistence diagram, changed the filtration parameter to $\tau$, changed the notation in Appendix C, replaced "(neg control)" by "(0 loops)", and inserted a line break after “, and then added isotropic Gaussian noise samples from …”. Thank you for these suggestions!

Minor 2. MDS of highly noised circle:
It is true that one can still see a hole in the middle of Figure 3d. But this is a 2D MDS and cannot capture the exact shape of the data in $\mathbb R^{50}$. Indeed, we see in Figure 3c that the feature for the correct loop disappears into the cloud of noise features (green) and in Figure 1a,b we see that the most persistent loop actually does not wrap around the noisy circle.

Minor 3. Definition of hitting times:
The hitting time $H_{ij}$ is the expected number of edges that a random walker starting at node $i$ traverses before reaching node $j$ for the first time, as stated in Definition 2. We changed the ambiguous term "average time" in the main text to "average number of edges".

Minor 4. Notation for identity matrix:
By $\mathbb I_d$ we always denote the $d\times d$ identity matrix and never the matrix of all ones.

Minor 5. Overloaded variable $D$:
We used $D$ both to denote the embedding dimension of Laplacian Eigenmaps and the degree matrix of the symmetric kNN graph. We changed the former to $\tilde{d}$. Thanks for your keen eye!

Minor 6. Complexity of persistent homology:
We added Myers, Audun D., et al. "Persistent homology of coarse-grained state-space networks." Physical Review E 107.3 (2023): 034303 as reference for the computational complexity of persistent homology. The term $d+1$ is correct: The complexity of the reduction algorithm is $\mathcal{O}(N^3)$, where $N$ is the number of simplices in the simplicial complex. A Vietoris-Rips complex can contain at most $\binom{n}{\delta}$ simplices of dimension $\delta$. To compute $\delta$-dimensional homology one needs the simplices of dimension $\delta-1, \delta, \delta+1$. Finally,

Minor 7. Capitalization of references:
We have originally capitalized the titles of references as in the original publications, but we are going to change it to a unified format in the next days.

Minor 8. Dimensions of the eyeglasses dataset:
You are right about the eyeglasses dataset. The length of the straight segments is actually $1.06$. We corrected the text.

Dear reviewer TJmS,
many thanks for your very detailed and helpful review! We really appreciate all the concrete suggestions and questions. They prompted many improvements to the manuscript. We will address your questions below.

(W1) Clarifications:
We clarified numerous aspects of the paper, and maintain that our statements are correct. If any questions remain, we are happy to address them as well.

(W2) Performance of $t$-SNE and UMAP:
Please see our general comment on the performance of $t$-SNE and UMAP (as this point came up in several reviews).

(Q7a) Exclusion of embedding methods from Figure 5:
We excluded the $t$-SNE and UMAP results from Figure 5 partly for the above reasons and partly because this already busy figure would have become messy otherwise.
We ran a large number of experiments on several datasets. Therefore, we deliberately chose different methods and datasets for the figures in the main text to illustrate the key take-aways of our analysis in a concise way. Note that we show the results of all methods on all dataset, inlcuding $t$-SNE and UMAP in the added supplementary figures S10-S21.

(Q21) $t$-SNE performance with low perplexity:
The reason for the dip of the $t$-SNE embedding method for low noise and low perplexity is a curvy embedding of the circle with various bottlenecks that give rise to additional features of high persistence. We included a new Figure S9 explaining this in more detail.

(W3) Novelty:
Please, see statement on novelty to all reviewers.

(Q1) Focus of the paper:
Our work addresses high-dimensional Gaussian noise, however we added an experiment with outliers in the revised version, following reviewer BdHd's concern (Appendix G). As stated at the end of Section 4 and the beginning of Section 9, moving from low to high ambient dimension, outliers cease to be an issue for persistent homology, instead Gaussian noise becomes problematic. For this reason, we focus on Gaussian noise. As can be seen in the new Fig. S6, adding up to 10% outliers impacts all methods only minimally in high ambient dimension. While not the focus of our work, please note that effective resistance also handles outliers better than Euclidean distance in low ambient dimesion (Figure S6).

(Q7b) Adequate presentation of DTM and Fermat:
You are right that DTM and Fermat distances were introduced to handle outliers, as we write in the related work section and Section 7. We have added a sentence in Section 7 stressing that we apply them to a different type of noise in our paper.

(Q2) Using $\ell_\infty$ to deal with noise:
Thank you for this suggestion. We ran experiments on our toy datasets for $\ell_p$ distances with $p=3,5,7, \infty$. The results are close to or worse than for the Euclidean distance (Figure S21). We believe that the reason for this is that for any $p\neq 2$ the $\ell_p$ distance favors directions and is not isotropic anymore. However, we place our toy datasets in a random orientation in ambient $\mathbb R^{50}$, so that the directions in which the dataset varies meaningfully are not aligned with the coordinate dimensions.

(Q3) Term "ring":
Thank you for spotting our imprecision. When speaking of "ring" we always meant a circle (plus Gaussian noise). We have changed all uses of "ring" to a more precise phrasing, as suggested.

Moreover, we have included Figure S4 that shows the circle with different levels of Gaussian noise in 2D as well as 2D MDS embeddings of the circle with different levels of Gaussian noise in 50D.

(Q4) Node weights in DTM:
We deliberately omitted this detail since, as you acknowledge, it does not matter for our experiments and omitting it simiplified the exposition. In the interest of accurarcy, we have added this aspect to the detailed description of DTM in Appendix C. We tweaked the relevant sentence in the related work to "The main idea of most of these suggestions is to replace the Euclidean distance with a different distance matrix, before running persistent homology. "