Large Scale computation of Means and Clusters for Persistence Diagrams using Optimal Transport
Persistence diagrams (PDs) are now routinely used to summarize the underlying topology of complex data. Despite several appealing properties, incorporating PDs in learning pipelines can be challenging because their natural geometry is not Hilbertian. Indeed, this was recently exemplified in a string of papers which show that the simple task of averaging a few PDs can be computationally prohibitive. We propose in this article a tractable framework to carry out standard tasks on PDs at scale, notably evaluating distances, estimating barycenters and performing clustering. This framework builds upon a reformulation of PD metrics as optimal transport (OT) problems. Doing so, we can exploit recent computational advances: the OT problem on a planar grid, when regularized with entropy, is convex can be solved in linear time using the Sinkhorn algorithm and convolutions. This results in scalable computations that can stream on GPUs. We demonstrate the efficiency of our approach by carrying out clustering with diagrams metrics on several thousands of PDs, a scale never seen before in the literature.
# Update after rebuttal  I thank the authors for their comments in the rebuttal. In conjunction with discussions with other reviewers, they helped me understand the salient points of the paper better. I would still add the following suggestions:  1. It would strengthen the paper if advantages over existing methods that are not Wasserstein-based could be showed (persistence landscapes, persistence images, ...), in particular since those methods have the disadvantage of suffering from descriptors that are not easily interpretable and do not yield 'proper' means (in the sense that the mean landscape, for example, cannot be easily converted back into a persistence diagram.  2. The experiment could be strengthened by discussing a more complex scenario that demonstrates the benefits of TDA _and_ the benefits of the proposed method. This could be a more complex data set, for example, that is not easily amenable to non-topological approaches, such as natural image data. I realize that this a balancing act, however.  Having discussed at length how to weigh the originality of the paper, I would like to see the method published in the hopes that it will further promulgate TDA. I am thus raising my score.  # Summary  This paper presents an algorithm for calculating approximate distances between persistence diagrams, i.e. feature descriptors from topological data analysis. The approximation is based on a reformulation of the distance calculation as an optimal transport problem. The paper describes how to apply previous research, most notably entropic regularization, to make calculations efficient and parallelizable. In addition to showing improved performance for distance calculations, the paper also shows how to calculate barycentres in the space of persistence diagrams. Previously, this involved algorithm that are computationally infeasible at larger scales, but the new formulations permits faster calculations. This is demonstrated by deriving a $k$-means algorithm.  # Review  This is an interesting paper and I enjoyed reading it. I particularly appreciated the step-by-step 'instructions' for how to build the algorithm, which increases accessibility even though the topic is very complex with respect to terminology. However, I lean towards rejecting the paper for now: while the contribution is solid and reproducible, the paper appears to be a (clever!) application of previous works, in particular the ones by Altschuler et al., Benamou et al., Cuturi, as well as Cuturi and Doucet. In contrast to these publications, which have an extremely general scope (viz. the whole domain of optimal transport), this publication is only relevant within the realm of topological data analysis.  Moreover, the experimental section is currently rather weak and does not give credence to the claim in the abstract that 'clustering is now demonstrated at scales never seen before in the literature'. While I agree that this method is capable of handling larger persistence diagrams, I think the comparison with other algorithms needs to have more details: the 'Hera' package, as far as I know, also permits the calculation of approximations to the Wasserstein distance and to the bottleneck distance between persistence diagrams. Were these used in the experiments? And if so, what were the parameters such as the 'relative error'? And for that matter, how many iterations of the new algorithm were required and what are the error bounds?  In addition, it might be fairer to use a single-core implementation for the comparison here---this would also show the benefits of parallelization! Such plots could be shown in the supplementary materials of the paper. At present, the observations with respect to speed are somewhat 'obvious': it seems clear that a multi-core algorithm (on a P100!) can always beat a single-core implementation. The same goes for the barycentre calculations, which are not even written in the same programming language I assume. Instead of showing that the algorithm is faster here (which I will definitely believe) can the authors say more about the qualitative behaviour? Are the barycentres essentially the same? How long does convergence take?  Statements of this sort would also strengthen reproducibility; the authors could show all required parameters in addition to their experiments.  Minor comments:  - L. 32: If different kernels for persistence diagrams are cited, the first one, i.e. the multi-scale kernel by Reininghaus et al. in "A Stable Multi-Scale Kernel for Topological Machine Learning" should also be cited. - In the introduction, works about approximations of persistence diagram distances (such as the one by Kerber et al.) could be briefly discussed as well. Currently, they only appear in the experimental section, which is a little bit late, I think. - L. 57: "verifying" --> "satisfying" (this happens at other places as well) - L. 114: $\mathcal{M}_{+}$ should be briefly defined/introduced before usage - How are parameters tuned? While I understand how the 'on-the-fly' error control is supposed to work, I am not sure how reasonable defaults can be obtained here. Or is the algorithm always iterated until the approximation error is sufficiently small? - L. 141: Please add a citation for the C-transform, such as Santambrogio (who is already cited in the bibliography) - L. 155: "criterion" --> "threshold" - L. 168: What does $|D_1|$ refer to? Some norm of the diagram points? This is important in order to understand the discretization error. - L. 178: From what I can tell, the algorithm presented in the paper generalizes to arbitrary powers; does the statement about the convolutions only hold for the squared Euclidean distance? The authors should clarify this. - L. 231: "Turner et al." is missing a citation - Timing information for the first experiment in Section 5 should be included; it is important to know the dependencies between the 'fidelity' of the approximation and the computational speed; in particular since the paper speaks about 'smearing' in L. 251, which prohibits the recovery of 'spiked' diagrams