Diffeomorphic interpolation for efficient persistence-based topological optimization

Authors: mathieu  carrière
marc  theveneau
théo  lacombe
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 11:54:39

Abstract: Topological Data Analysis (TDA) provides a pipeline to extract quantitative and powerful topological descriptors from structured objects. 
This enables the definition of topological loss functions, which assert to which extent a given object exhibits some topological properties. 
One can then use these losses to perform topological optimization via gradient descent routines. 
While theoretically sounded, topological optimization faces an important challenge: gradients tend to be extremely sparse, in the sense that the loss function typically depends (locally) on only very few coordinates of the input object, yielding dramatically slow optimization schemes in practice. 

In this work, focusing on the central case of topological optimization for point clouds, we propose to overcome this limitation using diffeomorphic interpolation, turning sparse gradients into smooth vector fields defined on the whole space. 
In particular, this approach combines efficiently with subsampling techniques routinely used in TDA, as the diffeomorphism derived from the gradient computed on the subsample can be used to update the coordinates of the full and possibly large input object. We then illustrate the usefulness of our approach on black-box autoencoder (AE) regularization, where we aim at applying some topological priors on the latent spaces associated to fixed, black-box AE models without modifying their (unknown) architectures and parameters. We empirically show that, while vanilla topological optimization has to be re-run every time that new data comes out of the black-box models, learning a diffeomorphic flow can be done once and then re-applied to new data in linear time. Moreover, reverting the flow allows us to generate data by sampling the topologically-optimized latent space directly, allowing for better interpretability of the model.

Main paper: 
 introduction: 
 Persistent homology (PH) is a central tool of Topological Data Analysis (TDA) that enables the extraction of quantitative topological information (such as, e.g., the number and sizes of loops, connected components, branches, cavities, etc) about structured objects (such as graphs, times series or point clouds sampled from, e.g., submanifolds), summarized in compact descriptors called persistence diagrams (PDs). PDs were initially used as features in Machine Learning (ML) pipelines; due to their strong invariance and stability properties, they have been proved to be powerful descriptors in the context of classification of time series [42, 19] , graphs [5, 23, 24] , images [2, 25, 16] , shape registration [7, 6, 36] , or analysis of neural networks [22, 3, 26] , to name a few.Another active line or research at the crossroad of TDA and ML is (persistence-based) topological optimization, where one wants to modify an object X so that it satisfies some topological constraints as reflected in its persistence diagram Dgm(X). The first occurrence of this idea appears in [21] , where one wants to deform a point cloud X ∈ R n×d so that Dgm(X) becomes as close as possible (w.r.t. an appropriate metric denoted by W ) to some target diagram D target , hence yielding to the problem of minimizing X → W (Dgm(X), D target ). This idea has then been revisited with different flavors, for instance by adding topology-based terms in standard losses in order to regularize ML models [14, 33, 29] , improving ML model reconstructions by explicitly accounting for topological features [16] , or improving correspondences between 3D shapes by forcing matched regions to have similar topology [36] . Formally, this goes through the minimization of an objective functionL : X → ℓ(Dgm(X)) ∈ R,where ℓ is a user-chosen loss function that quantifies to what extent Dgm(X) reflects some prescribed topological properties inferred from X. Under mild assumptions (see Section 2), the map L is differentiable generically and its gradients are obtained as a byproduct of the computation of Dgm(X). However, these approaches are limited in practice by two major issues: (i) the computation of X → Dgm(X) scales poorly with the size of X (e.g., number of points n in a point cloud, number of nodes in a graph, etc), and (ii) the gradient ∇L(X) tends to be very sparse: if X = (x 1 , . . . , x n ) ∈ R n×d is a point cloud, ∇L(X) i ̸ = 0 for only very few indices i ∈ {1, . . . , n} (the corresponding points are called the critical points of the topological gradient, see Section 2.1). Related works. Several articles have studied topological optimization in the TDA literature. The standard, or vanilla, framework to define and study gradients obtained from topological losses was described in [4, 28] , where the high sparsity and long computation times were first identified. To mitigate this issue, the authors of [34] introduced the notion of critical set that extends the usually sparse set of critical points in order to get a gradient-like object that would update more points in X. In [39] , the authors used an average of the vanilla topological gradients of several subsamples to get a denser and faster gradient. On the theoretical side, the authors of [27] demonstrated that adapting the stratified structure induced by PDs to the gradient definition enables faster convergence.Limitations. Despite proposing interesting ways to accelerate gradient descent, the approaches mentioned above are still limited in the sense that their proposed gradients are not defined on the whole space, but only on a sparse subset of the current observation X, which still prevents their use in different contexts, that we investigate in our experiments (Section 4). First, when the data has more than tens of thousands of points, the number of subsamples needed to capture relevant topological structures (when using [39] ), as well as the critical set computations (when using [34] ), both become practically infeasible. Second, when optimizing the topology of datasets obtained as latent spaces of a black-box autoencoder model (i.e., an autoencoder with forbidden access to its architecture, parameters, and training), then (a) the topological gradients of [39, 34] cannot be reused to process new such datasets, and topological optimization has to be performed from scratch every time that new data comes in, (b) this also impedes their transferability, as re-running gradient descent every time makes it very difficult to guarantee some stability for the final output, and finally (c) one cannot generate new data by sampling the optimized latent spaces directly, as it would require to apply the sequence of reverted gradients (which are not well-defined everywhere).In this article, we propose to replace the standard gradient ∇L(X) of (??) derived formally by a diffeomorphism v : R d → R d that interpolates ∇L(X) on its non-zero entries, that is v(x i ) = ∇L(X) i for all i ∈ I := {j | ∇L(X) j ̸ = 0} and is, in some sense, as smooth as possible. More precisely, our contribution is three-fold:• We introduce a diffeomorphic interpolation of the vanilla topological gradient, which extends this gradient to a smooth and denser vector field defined on the whole space R d , and which is able to move a lot more points in X at each iteration, • We prove that its updates indeed decrease topological losses, and we quantify its smoothness by upper bounding its Lipschitz constant (again, in the context of topological losses), • We showcase its practical efficiency: we show that it compares favorably to the main baseline [34] in terms of convergence speed, that its combination with subsampling [39] allows to process datasets whose sizes are currently out of reach in TDA, and that it can successfully be used for the tasks mentioned above concerning black-box autoencoder models.t = 0 t = 1 t = 2 t = 3 t = 5Section 2 provides necessary background in Topological Data Analysis and diffeomorphic interpolations. Section 3 presents our approach and its corresponding guarantees, and Section 4 showcases our experiments. Limitations and further research directions are discussed in Section 5. 
topological data analysis: 
 In this section, we recall the basic materials of Topological Data Analysis (TDA), and refer the interested reader to [20, 35] for a thorough overview. We restrict the presentation to our case of interest: extracting topological information from a point cloud using the standard Vietoris-Rips (VR) filtration. A more extensive presentation of the TDA machinery is provided in Appendix A.Let X = (x 1 , . . . , x n ) ∈ R n×d . The Vietoris-Rips filtration consists of building an increasing sequence of simplicial complexes (K t ) t≥0 over X by inserting a simplex σ = (x i1 , . . . , x ip ) whenever ∀j, j ′ ∈ {1, . . . , p}, ∥x ijx i j ′ ∥ ≤ t. Each time a simplex σ is inserted, it either creates a topological feature (e.g., inserting a face creates a cavity, that is a 2-dimensional topological feature) or destroy a pre-existing feature (e.g., the face insertion fills a loop, that is a 1-dimensional feature, making it topologically trivial). The persistent homology machinary tracks the apparition and destruction of such features in the so-called persistence diagram (PD) of X, denoted by Dgm(X). Therefore, Dgm(X) is a set of points in R 2 of the form (t b , t d ) with t d ≥ t b , where each such point accounts for the presence of a topological feature inferred from X that appeared at time t b following the insertion of an edge (x i1 , x i2 ) with ∥x i1x i2 ∥ = t b and disappeared at time t d following the insertion of an edge (x i3 , x i4 ) with ∥x i3x i4 ∥ = t d . Figure 1 illustrates this construction. From a computational standpoint, computing the VR diagram of X ∈ R n×d that would reflect topological features of dimensiond ′ ≤ d runs in O(n d ′ +2), making the computation quickly unpractical when n increases, even when restricting to low-dimensional features such as connected components(d ′ = 0), loops (d ′ = 1) or cavities (d ′ = 2).Topological optimization. PDs are made to be used in downstream pipelines, either as static features (e.g., for classification purposes) or as intermediate representations of X in optimization schemes. In this work, we focus on the second problem. We formally consider the minimization of objective functions of the formEQUATIONHere, ℓ represents a loss function taking value from the space of PDs, denoted by D in the following.The space D can be equipped with a canonical metric, denoted by W and whose formal definition is not required in this work, for which a central result is that the map X → Dgm(X) is stable (Lipschitz continuous) [17, 18, 38] . Therefore, if ℓ is Lipschitz continuous, so is L, hence it admits a gradient almost everywhere by Rademacher theorem. Building on these theoretical statements, one can consider the "vanilla" gradient descent update X k+1 := X k -λ∇L(X k ) for a given learningrate λ > 0 and iterate it in order to minimize (1) . Theoretical properties of this seminal scheme (and natural extensions, e.g., stochastic gradient descent) have been studied in [4, 27] , where convergence (to a local minimum of L) is proved.From a computational perspective, deriving ∇L(X) comes in two steps. Let µ := Dgm(X), written as µ = {(b i , d i ) | i ∈ I} for some finite set of indices I. To each i ∈ I correspond four (possibly coinciding) points x i1 , x i2 , x i3 , x i4 in the input point cloud X. Intuitively, minimizing µ → ℓ(µ) boils down to prescribe a descent direction (δb i , δd i ) ∈ R2 to each (b i , d i ) for i ∈ I, where δb i = ∂ℓ ∂bi (µ) and δd i = ∂ℓ ∂di (µ). Backpropagating this perturbation to X will move the corresponding points x i1 , x i2 , x i3 , x i4 in order to increase or decrease the distances ∥x i1x i2 ∥ = b i and ∥x i3x i4 ∥ = d i accordingly. This yields to the following formula:EQUATIONwhere the notation x → (b i , d i ) means that x ∈ X appears in (at least) one of the four points yielding the presence of (b i , d i ) in the diagram µ = Dgm(X). A fundamental contribution of [28, §3.3 ] is to prove that the chain rule formula (2) is indeed valid 2 . Most of the time, a point x ∈ X will not belong to any critical pair (σ b , σ d ) and the above gradient coordinate is 0. Therefore, the gradient of L depends on very few points of X, yielding the sparsity phenomenon discussed in Section 1.Examples of common topological losses. Let X = (x 1 , . . . , x n ) ∈ R n×d be a point cloud and Dgm(X) = {(b i , d i ) | i ∈ I} be its PD. There are several natural losses that have been introduced in the TDA literature:• Topological simplification losses: typically of the form i∈ Ĩ (b id i ) 2 , where Ĩ ⊆ I. Such losses push (some of the) points in Dgm(X) toward the diagonal ∆ = {b = d}, hence destroying the corresponding topological features appearing in X. • Topological augmentation losses [4] : similar to simplification losses, but typically attempting to push points in Dgm(X) away from ∆, i.e., minimizing -i∈ Ĩ (b id i ) 2 , to make topological features of X more salient. As such losses are not coercive, they are usually coupled with regularization terms to prevent points in X going to infinity. • Topological registration losses [21] : given a target diagram D target , one minimizes W (Dgm(X), D target ) where W denotes a standard metric between persistence diagrams. This loss attempts to modify X so that it exhibits a prescribed topological structure.In order to overcome the sparsity of gradients appearing in TDA, we rely on diffeomorphic interpolations (see, e.g., [43, Chapter 8] ). Let X = (x 1 , . . . , x n ) ∈ R n×d , let I ⊆ {1, . . . , n} denote the set of indices on which ∇L(X) is non-zero and let a i := (∇L(X)) i ∈ R d for i ∈ I. Our goal is to find a smooth vector field ṽ : R d → R d such that, for all i ∈ I, ṽ(x i ) = a i . To formalize this, we consider a Hilbert spaceH ⊂ (R d ) R d for which the map δ α x : f → ⟨α, f (x)⟩ R d = α T f (x) is continuous for any (α, x) ∈ R d × R d .Such a space is called a Reproducing Kernel Hilbert Space (RKHS) 3 . A crucial property is that there exists a matrix-valued kernel operator K : R d ×R d → R d×d whose outputs are symmetric and positive definite, and related to H through the relation ⟨k αx , k β y ⟩ H = α T K(x, y)β for all x, y, α, β ∈ R d , where k αx ∈ H is the unique vector provided by the Riesz representation theorem such that ⟨k α x , f ⟩ = ⟨α, f (x)⟩. Conversely, any such kernel K induces a (unique) RKHS H (of which K is the reproducing kernel). Now, we can consider the following problem:EQUATIONthat is, we are seeking for the smoothest (lowest norm) element of H that solves our interpolation problem. The solution ṽ of this problem is the projection of 0 onto the affine set{v ∈ H | v(x i ) = a i , ∀i ∈ I}. Observe that ṽ belongs to the orthogonal of {v ∈ H | v(x i ) = 0, ∀i ∈ I}, and thus of {v ∈ H | ⟨k αi xi , v⟩ H = 0, ∀i ∈ I, α i ∈ R d },and therefore ṽ ∈ span({k αi xi | i ∈ I}). This justifies to search for ṽ in the form of ṽ(x) = i∈I K(x, x i )α i , and the interpolation that it must satisfy yields ṽ(x) = i∈I K(x, x i )(K -1 a) i , where K is the block matrix (K(x i , x j )) i,j∈I and a = (a i ) i∈I . See also [43, Theorem 8.8] . In particular, it is important to note that ṽ inherits from the regularity of K and will typically be a diffeomorphism in this work. If K is the Gaussian kernel defined by K(x, y) := exp -∥x-y∥ 2 2σ 2 I d for some bandwidth σ > 0, a choice to which we stick to in the rest of this work, the expression of ṽ reduces toEQUATIONwhere ρ σ (u) := e -u 2 2σ 2 , andα i := (K -1 a) i with K = (ρ σ (∥x i -x j ∥)I d ) i,j∈I .Note that ṽ can be understood as the convolution of a with a Gaussian kernel, but involving a correction K -1 guaranteeing that after the convolution, the interpolation constraint is satisfied. We will call ṽ the diffeomorphic interpolation associated to the vectors and indices {a i | i ∈ I}.3 Diffeomorphic interpolation of the vanilla topological gradient 1.0 Point cloud L(X) (Vanilla) v(X) (Diffeo)Figure 2: (blue) A point cloud X, and (black) the negative gradient -∇L(X) of a simplification loss which aims at destroying the loop by collapsing the circle (reduce the loop's death time) and tearing it (increase the birth time). While ∇L(X) only affects four points in X, the diffeomorphic interpolation ṽ(X) (orange, σ = 0.1) is defined on R d , hence extends smoothly to other points in X.We aim at minimizing a loss function L : X → ℓ(Dgm(X)) as in ( 1), starting from some initialization X 0 , and assuming that L is lower bounded (typically by 0) and locally semi-convex. This assumption is typically satisfied by the topological losses ℓ introduced in Section 2.1. Gradient descents implemented in practice are (explicit) discretization of the gradient flowEQUATIONwhere ∂L(X):= {v | L(Y ) ≥ L(X) + v • (Y -X) + o(Y - X)for all X, Y } denotes the subdifferential of L at X. Note that a topological loss L is typically not differentiable everywhere, since the map X → Dgm(X) is differentiable almost everywhere but not in C In this work, we propose to use the dynamic described by the diffeomorphism ṽt introduced in (4) interpolating the current vanilla topological gradient ∇L(X t ) at each time t, formally considering solutions X ofEQUATIONHere, slightly overloading notation, ṽt ( X(t)) denotes the n×d matrix where the i-th line is given by ṽt ( X(t) i ). The flow at time T associated to ( 6) is the mapEQUATIONwhich can inverted by simply following the flow backward (i.e., by following ṽt instead of -ṽ t ). We now guarantee that at each time t, following ṽt instead of the vanilla topological gradient ∇L(X t ) still provides a descent direction for the topological loss L.Proposition 3.1. For each t ≥ 0, it holds that dL( X(t)) dt = -∥∇L( X(t))∥ 2 ≤ 0.Proof. One has dL( X(t)) dt = -⟨∇L( X(t)), ṽt ( X(t))⟩ = - n i=1 (∇L( X(t))) i • (ṽ t ( X(t))) i .Since ∇L( X(t)) i = 0 for i ̸ ∈ I, and ṽt ( X(t)) i = -∇L( X(t)) i for i ∈ I, the result follows. Moreover, it is also possible to upper bound the smoothness, i.e., the Lipschitz constant, of ṽ: Proposition 3.2. Let L be the simplification or augmentation loss computed with k = | Ĩ| PD points, as defined at the end of Section 2.1. Let ṽ = ṽt be the diffeomorphic interpolation associated to the vanilla topological gradient at time t ≥ 0. Then, one has, ∀x, y ∈ R d and t ≥ 0:∥ṽ(x) -ṽ(y)∥ 2 ≤ ∥ṽ(x) -ṽ(y)∥ 1 ≤ C d • σ d-1 • κ(K) • Pers k (Dgm( X(t))) • ∥x -y∥ 2 ,whereC d = √ d • 2 3+ d+1 2 • π d-1 2 , κ(K)is the condition number of K, and Pers k (Dgm( X(t))) is the sum of the k largest distances to the diagonal in Dgm( X(t)).The proof is deferred to Appendix B. This upper bound can be used to quantify how smooth the diffeomorphic interpolation is (as characterized with its Lipschitz constant) based on the parameters it is computed from. In the case of the Gaussian kernel, we found that our upper bound on the Lipschitz constant depends on the kernel bandwidth σ: indeed, the more spread the Gaussian function is, the more critical points can influence other data points potentially far from them, introducing unwanted distortions and larger Lipschitz constant. Similarly, if the condition number κ(K) is large, inverting the kernel matrix might introduce instabilities in Equation ( 4), and thus a larger Lipschitz constant as well. Finally, the total persistence also appears, as the more PD points one has to optimize, the more critical pairs will appear, and thus the more constrained Equation ( 3) is, leading to more complex diffeomorphic interpolation solutions with larger Lipschitz constants. 
subsampling techniques to scale topological optimization: 
 As a consequence of the limited scaling of the Vietoris-Rips filtration with respect to the number of points n of the input point cloud X, it often happens in practical applications that computing the VR diagram Dgm(X) of a large point set X (a fortiori its gradient) turns out to be intractable. A natural workaround is to randomly sample s-points from X, with s ≪ n, yielding a smaller point cloud X ′ ⊂ X. Provided that the Hausdorff distance between X ′ and X is small, the stability theorem [18, 17, 8] ensures that Dgm(X ′ ) is close to Dgm(X). See [11, 12, 9] for an overview of subsampling methods in TDA.However, the sparsity of vanilla topological gradients computed from topological losses strikes further when relying on subsampling: only a tiny fraction of the seminal point cloud X is likely to be updated at each gradient step. In contrast, using the diffeomorphic interpolation ṽ (of the vanilla topological gradient) computed on the subsample X ′ still provides a vector field defined on the whole input space R d , in particular on each point of X and the update can then be performed in linear time with respect to n. This yields Algorithm 1. Figure 3 illustrates the qualitative benefits offered by the joint use of subsampling and diffeomorphic interpolations when compared to vanilla topological gradients. A larger-scale experiment is provided in Section 4.Input: Initial X 0 ∈ R n×d , loss function ℓ, learning rate λ > 0, subsampling size s ∈ {1, . . . , n}, max. epoch T ≥ 1, stopping criterion. Set L : X → ℓ(Dgm(X)) (+ possibly a regularization term in X).for k = 1, . . . , T do Subsample X ′ k-1 = {x ′ 1 , . . . , x ′ s } uniformly from X k-1 . Compute ∇L(X ′ k-1 ) (vanilla topological gradient) Compute the diffeomorphic interpolation ṽ(X ′ k-1 ) from ∇L(X ′ k-1 ) using (4). Set X k := X k-1 -λṽ(X k-1 ). if stopping criterion is reached then Return X k end if end for Return X TStopping criterion. A natural stopping criterion for Algorithm 1 is to assess whether the loss L(X t ) = ℓ(Dgm(X t )) is smaller than some ε > 0. However, computing Dgm(X t ) can be intractable if X t is large. Therefore, a tractable loss to consider is L(X t ) := E[ℓ(Dgm(X ′ t ))], where X ′ t is a uniform s-sample from X t . Under that perspective, Algorithm 1 can be re-interpreted as a kind of stochastic gradient descent on L, for which two standard stopping criteria can be used: (a) compute an exponential moving average of the loss on individual samples X ′ t over iterations, or (b) compute a validation loss, i.e., sample X ′ t,(1) , . . . , X ′ t,(K) and estimate L byK -1 K k=1 ℓ(Dgm(X ′ t,(k) )). Empirically, we observe that the latter approach with K = n/s (more repetitions for smaller sample sizes to mitigate variance) yields the most satisfactory results (faster convergence toward a better objective X t ) overall, and thus stick to this choice in our experiments. 
numerical experiments: 
 We provide numerical evidence for the strength of our diffeomorphic interpolations. PHrelated computations relies on the library Gudhi [40] and automatic differentiation relies on tensorflow [1] . The "big-step gradient" baseline [34] implementation is based on oineus4 . The first two experiments were run on a 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, the last one on a 2x Xeon SP Gold 5115 @ 2.40GHz. Our code is publicly available at https: //github.com/tlacombe/topt.Convergence speed and running times. We sample uniformly n = 200 points on a unit circle in R 2 with some additional Gaussian noise, and then minimize the simplification loss L : X → (b,d)∈Dgm(X) |d| 2 , which attempts to destroy the underlying topology in X by reducing the death times of the loops by collapsing the points. The respective gradient descents are iterated over a maximum of 250 epochs, possibly interrupted before if a loss of 0 is reached (Dgm(X) is empty), with a same learning rate λ = 0.1. The bandwidth of the Gaussian kernel in (4) is set to σ = 0.1. We include the competitor oineus [34] , as-even though relying on a fairly different construction-this method shares a key idea with ours: extending the vanilla gradient to move more points in X. We stress that both approaches can be used in complementarity: compute first the "big-step gradient" of [34] using oineus, and then extend it by diffeomorphic interpolation. Results are displayed in Figure 4 5 . In terms of loss decrease over iterations, both "big-step gradients" and our diffeomorphic interpolations significantly outperform vanilla topological gradients, and their combined use yields the fastest convergence (by a slight margin over our diffeomorphic interpolations alone). However, in terms of raw running times, the use of oineus involves a significant computational overhead, making our approach the fastest to reach convergence by a significant margin.Subsampling. We now showcase how using our diffeomorphic interpolation jointly with subsampling routines (Algorithm 1) allows to perform topological optimization on point clouds with thousands of points, a new scale in the field. For this, we consider the vertices of the Stanford Bunny [41] , yielding a point cloud X 0 ∈ R n×d with n = 35, 947 and d = 3. We consider a topolog- ical augmentation loss (see Section 2.1) for two-dimensional topological features, i.e., we aim at increasing the persistence of the bunny's cavity. The size of n makes the computation of Dgm(X 0 ) untractable (recall that it scales in O(n 4 )); we thus rely on subsampling with sample size s = 100 and compare the vanilla gradient descent scheme and our scheme described in Algorithm 1. Results are displayed in Figure 5 . Because it only updates a tiny fraction of the initial point cloud at each iteration, the vanilla topological gradient with subsampling barely changes the point cloud (nor decreases the loss) in 1,000 epochs. In sharp contrast, as our diffeomorphic interpolation computed on subsamples is defined on R 3 , it updates the whole point cloud at each iteration, making possible to decrease the objective function where the vanilla gradient descent is completely stuck. Note that a step of diffeomorphic interpolation, in that case, takes about 10 times longer than a vanilla step. An additional subsampling experiment can be found in Appendix C.Black-box autoencoder models. Next, we apply our diffeomorphic interpolations to black-box autoencoder models. In their simplest formulation, autoencoders (AE) can be summarized as two maps E : R d → R d ′ and D : R d ′ → R d called encoder and decoder respectively. The intermediate space R d ′ in which the encoder E is valued is referred to as a latent space (LS), with typically d ′ ≪ d. In general, without further care, there is no reason to expect that the LS of a point cloud X, E(X) = {E(x 1 ), . . . , E(x n )}, reflects any geometric or topological properties of X. While this can be mitigated by adding a topological regularization term to the loss function during the training of the autoencoder [29, 4] , this cannot work in the setting where one is given a black-box, pretrained AE. However, replacing (E, D) by (φ • E, D • φ -1 ) for any invertible map φ : R d ′ → R d ′ yields an AE producing the same outputs yet changing the LS E(X), without explicit access to the AE's model. Hence, we propose to learn such a φ with diffeomorphic interpolations: given some latent space E(X), we apply T steps of our diffeomorphic gradient descent algorithm to X → ℓ(Dgm(X)) initialized at E(X). We thus get a sequence of smooth displacements -ṽ 1 , . . . , -ṽ T of R d ′ that discretizes the flow (7) via φ :x 0 → x 0 - T k=1 ṽk (x k-1 ) where x k -x k-1 = -ṽ k (x k-1), and such that Dgm(φ(E(X))) is more topologically satisfying. This fixed diffeomorphism φ can then be re-applied to any new data coming out of the encoder in a deterministic way. Moreover, any random sample from the topologically-optimized LS can be inverted without further computations by following ṽT , ṽT -1 , . . . , ṽ1 , which allows to push the new sample back to the initial LS, and then apply the decoder on it. Again, this cannot be achieved with baselines [39, 34] .In order to illustrate these properties, we trained a variational autoencoder (VAE) to project a family of datasets of images representing rotating objects, named COIL [32] , to two-dimensional latent spaces. Given that every dataset in this family is comprised of 288 pictures of the same object taken with different angles, one can impose a prior on the topology of the corresponding LSs, namely that they are sampled from circles. However, the VAE architecture is shallow (the encoder has one fullyconnected layer (100 neurons), and the decoder has two (50 and 100 neurons), all layers use ReLu activations), and thus the learned latent spaces, although still looking like curves thanks to continuity, do not necessarily display circular patterns. This makes generating new data more difficult, as latent spaces are harder to interpret. To improve on this, we learn a flow φ as described above with an augmentation loss associated to the 1-dimensional PD point which is the most far away from the diagonal, in order to force latent spaces to have a significant one-dimensional topological feature, i.e., loop. As the datasets are small, we do not use subsampling, and we use learning rate λ = 0.1, Gaussian kernels with bandwidth σ = 0.3 and an increase of at least 3. in the topological loss (from an iteration to the next) to stop the algorithm 6 . We provide some qualitative results in Figure 7 (see also Appendix C, Figure 10 ). In order to quantify the improvement, we also computed the Pearson correlation scores between the ground-truth angles θ i and the angles θi computed from the topologically-optimized LSs with θi := ∠(φ•E(x i )-Ê[φ•E(X)], φ•E(x 1 )-Ê[φ•E(X)]),where Ê[φ • E(X)] := n -1 n i=1 φ • E(x i ), and φ denotes our flow. In Table 1 , we provide an average of these scores over 100 test sets obtained by randomly perturbing the training set with uniform noise of amplitude 0.05. As expected, correlation becomes better after forcing the latent spaces to have the topology of a circle. This better interpretability is also illustrated in Figure 6 , in which four angles are specified, which are mapped to the topologically-optimized LS, then pushed to the initial LS of the black-box VAE by following the reverted flow of our learned diffeomorphism φ, and finally decoded back into realistic, COIL-like images.Single-cell data. We also deployed our proposed topological correction of LSs from AEs on a realworld dataset of single cells. Specifically, we designed an experiment on single cell HiC (scHiC) data inspired from [27] . The dataset is comprised of single cells characterized by chromatin folding, that is, each cell is encoded by the spatial distance matrix of its DNA fragments. The dataset we focus on is taken from [31] , in which it was shown that cells are sampled at different phases of the cell cycle. Thus, similar to COIL images, we expect latent embeddings of this dataset to exhibit a circular shape, that we can constrain with diffeomorphic topological optimization. Specifically, we processed this single cell dataset of 1, 171 cells with the stratum-adjusted correlation coefficient (SCC) with 500kb and convolution parameter h = 1 on chromosome 10. Then, we run kernel PCA on the SCC matrix to obtain a preprocessed dataset in R 100 , on which we applied the same VAE architecture than the one described above for COIL images. Finally, we optimized two losses, the first was the same augmentation loss than for the COIL images, the second was the following registration loss:L : X ∈ R n×d → W 2 (Dgm 1 (X), D target ),where Dgm(X) contains the points of the PD of X with distance-to-diagonal at least τ = 1, D target is a target PD with only one point p * = [-3.5, 3.5], and W is the 2-Wasserstein distance between PDs. We used σ = 0.2 for the augmentation loss and σ = 0.025 for the registration loss (σ is set to a smaller value for the registration loss in order to mitigate the effects of matching instability), λ = 0.1 on 500 epochs, with subsampling of size s = 300 for computational efficiency (as computing VR diagrams without radius thresholding on 1, 171 points already takes few minutes on a laptop CPU, which becomes hardly tractable if done repetitively as in gradient descent), and loss increase of 3. as a stopping criterion. Qualitative results are displayed in Appendix C, Figure 11 , and we also measured quantitative improvement with the correlation scores between latent space angles and repli scores7 in Table 1 , in which improvements can be observed. An additional experiment on the influence of the bandwidth parameter σ over these correlation scores, as well as over convergence, can also be found in Appendix C. 
conclusion: 
 In this article, we have presented a way to turn sparse topological gradients into dense diffeomorphisms with quantifiable Lipschitz constants, and showcased practical benefits of this approach in terms of convergence speed, scaling, and applications to black-box AE models on several datasets. Several questions are still open for future work.In terms of theoretical results, we plan on working on the stability between the diffeomorphic interpolations computed on a dataset and its subsamples. This requires some control over the locations of the critical points, which we expect to be possible in statistical estimation; indeed sublevel sets of density functions are know to have stable critical points [10, Lemma 17] . We also plan to look at adaptive kernels whose parameters (like the bandwidth σ) depend on the input point cloud σ = σ(X) (instead of using a fixed kernel and parameters at every iteration), and understand the convergence properties of our proposed diffeomorphic gradient descent. Finally, applying our diffeomorphic interpolation to sparse gradients computed with multiparameter persistent homology is another natural research direction, provided that differentiability properties have recently been proved in that setting [30, 37] .Concerning the AE experiment, we plan to investigate the limitations presented in the figure below: as the initial LSs (colored with the ground-truth angles) have zero (left) or two (right) loops, it is impossible to unfold them with diffeomorphisms; instead the optimized latent spaces either also exhibit no topology, or mixes different angles. In future work, we plan on investigating other losses or gradient descent schemes for diffeomorphic topological optimization, including stratified procedures similar to [27] that allow for local topological changes during training. 
 
 Appendix: The starting point of TDA is to extract quantitative topological information from structured objects-for example graphs, points sampled on a manifold, time series, etc. Doing so relies on a piece of algebraic machinery called persistent homology (PH), which informally detects the presence of underlying topological properties in a multiscale way. Here, estimating topological properties should be understood as inferring the number of connected components (topology of dimension 0), the presence of loops (topology of dimension 1), of cavities (dimension 2), and so on in higher dimensional settings.Simplicial filtrations. Given a finite simplicial complex8 K, a filtration over K is a map t → K t ⊆ K that is non-decreasing (for the inclusion). For each σ ∈ K, one can record the value t(σ) := inf{t | σ ∈ K t } at which the simplex σ is inserted in the filtration. From a topological perspective, the insertion of σ has exactly one of two effects: either it creates a new topological feature in K t (e.g., the insertion of an edge can create a loop, that is, a one-dimensional topological feature) or it destroys an existing feature of lower dimension (e.g., two independent connected components are now connected by the insertion of an edge). Relying on a matrix reduction algorithm [20, §IV.2], PH identifies, for each topological feature appearing in the filtration, the critical pair of simplices (σ b , σ d ) that created and destroyed this feature, and as a byproduct the corresponding birth and death times t(σ b ), t(σ d ). 9 The collection of intervals(t(σ b ), t(σ d )) is a (finite) subset of the open half-plane {(b, d) ∈ R 2 | b < d}, called the persistence diagram (PD) of the filtration (K t ) t . The distance of such a point (t b , t d ) to the diagonal {b = d}, namely 2 -1 2 |t(σ b ) -t(σ d )|, is called the persistence of the corresponding topological feature, as an indicator of "for how long" could this feature be detected in the filtration (K t ) t .Note that if (σ b , σ d ) is a critical pair for our filtration (K t ) t , it holds that |σ b | = |σ d |+1. The quantity |σ b | -1 is the dimension of the corresponding topological feature (e.g., loops, which are created by the insertion of an edge and killed by the insertion of a triangle, are topological features of dimension one). From a computational perspective, deriving the PD of a filtration (K t ) t is empirically 10 quasilinear with respect to the number of simplices in K (which can still be extremely high in the case of Vietoris-Rips filtration-see below-where K = 2 V with |V | = n being typically quite large).The Vietoris-Rips filtration. A particular instance of simplicial filtration that will be extensively used in this work is the Vietoris-Rips (VR) one. Given X = (x 1 , . . . , x n ) ∈ R n×d a point cloud of n points in dimension d, one consider the simplicial complex K = 2 X and then the filtration (K t ) t defined by σ = {x i1 . . . x ip } ∈ K t ⇐⇒ ∀j, j ′ ∈ {1, . . . , p}, ∥x ij -EQUATIONThe corresponding persistence diagram will be denoted, for the sake of simplicity, by Dgm(X).Note that for t < 0, K t = ∅, when t ≥ diam(X), K t = K, and there is always a point with coordinates (0, +∞) in Dgm(X) accounting for the remaining connected component when t → ∞. This is the unique point in Dgm(X) for which the second coordinate is +∞ (and is often discarded in practice, as it does not play any significant role). Intuitively, Dgm(X) reflects the topological properties that can be inferred from the geometry of X; this can be formalized by various results which state, roughly, that if the x i are i.i.d. samples from a regular measure µ supported on a submanifold M ⊂ R d , then with high probability the topological properties of M are reflected in Dgm(X) (see [13, 9] ). From a computational perspective, note that the VR filtration only depends on X through the pairwise distance matrix (∥x i -x j ∥) 1≤i,j≤n , and thus the complexity of computing Dgm(X) depends only linearly in d11 . On the other hand, since Dgm(X) scales (at least) linearly ). This loss favors the apparition of topological features (loops) while the regularization term penalizes points that would go to infinity otherwise.-Experiment reproduced following the setting of [4] , using code available at https://github.com/GUDHI/TDA-tutorial/blob/ master/Tuto-GUDHI-optimization.ipynb.with respect to the number of simplices in K, computing the whole VR diagram of a point cloud X ∈ R n×d can take up to O(2 n ) operations. Even if one restricts topological features of dimension d ′ ≤ d (e.g. d ′ = 1 if one only considers loops)-as commonly done-the complexity is of order O(n d ′ +2 ), which becomes quickly intractable if n is large, even if d ′ = 1 or 2.Proof of Proposition 3.2. One has: ∥ṽ(x) -ṽ(y)∥ 1 = ∥ i∈I (K(x, x i ) - K(y, x i ))(-K -1 ∇L( X(t))) i ∥ 1 ≤ i∈I |ρ σ (∥x -x i ∥) -ρ σ (∥y -x i ∥)| • ∥(-K -1 ∇L( X(t))) i ∥ 1 , since we are using Gaussian kernels. As |ρ σ (∥x -x i ∥) -ρ σ (∥y -x i ∥)| ≤ C d,σ ∥x -y∥ 2 , with C d,σ = 2 d+1 2 π d-1 2 σ d-1 (see [2, Theorem 8]), it follows that ∥ṽ(x) -ṽ(y)∥ 1 ≤ C d,σ ∥x -y∥ 2 • ∥K -1 ∥ 1 • ∥∇L( X(t))∥ 1 .Let us upper bound the term ∥∇L( X(t))∥ 1 . Let us start with the simplification loss, one has ∂ℓ ∂bi = 2(b id i ) and, writing b i = ∥x i1x i2 ∥ 2 (for some critical points x i1 , x i2 ∈ X(t)), one has∂bi ∂xi 1 = xi 1 -xi 2 ∥xi 1 -xi 2 ∥2 and ∂bi ∂xi 2 = - xi 1 -xi 2 ∥xi 1 -xi 2 ∥2. Similarly, one has ∂ℓ ∂di = -2(b id i ), and writing d i = ∥x i3x i4 ∥ 2 provides the corresponding partial derivatives.Applying (2) , this gives:EQUATION→) means that x appears as left point (resp. right point) in the computation of the birth filtration value b i of one of the k PD points (b i , d i ) ∈ Dgm( X(t)) associated to the loss, and similarly for death filtration values. A brutal majoration finally gives ∥∇L(X(t))∥ 1 ≤ 2 √ d x∈ X(t) x→(bi,di) |b i -d i | ≤ 8 √ d • Pers k (Dgm( X(t))), as there are at most four points associated to every (b i , d i ) ∈ Dgm( X(t)). One can easily see that the same bound applies to the augmentation loss.Let us finally bound ∥K -1 ∥ 1 , one has ∥K -1 ∥ 1 = κ(K)/∥K∥ 1 ≤ κ(K). Indeed, as we are using Gaussian kernels,∥K∥ 1 = max 1≤i≤n n j=1 ρ σ (∥x i -x j ∥ 2 ) ≥ 1.Subsampling and improving over [4] . We reproduce the experiment of [4, §5] , see also Figure 8 , but starting from an initial point cloud X 0 of size n = 2, 000 instead of n = 300. This makes the phic interpolation. The increase in running time with respect to the previous experiment mostly lies on instantiating and applying the n × d (d = 3) vector field ṽ (requires to compute ρ i (xx i ) for each new x (n of them) and sampled x i (|I| of them, which is typically very small), hence a ∼ O(n) complexity).• Training the VAE for the COIL and scHiC datasets are the most computationally expensive parts of this work: it takes about 3 hours per image (20 of them), and 3 hours also for the scHiC dataset. In contrast, performing the topological optimization take few dozen of minutes (less than one hour) for each image. Applying it is done in few seconds at most. Recall that our method is designed to handle pre-trained models (which may be way more sophisticated than the one we used!); and its running time does not depend on the complexity of the model.Influence of the bandwidth σ. We reproduce the same experimental setting as in Figure 3 , i.e., sample points uniformly on a circle of radius 1 plus additional noise ∼ N (0, 0.05I 2 ), and consider minimizing the total persistence of the point cloud. We take σ ∈ {0, 0.1, 0.2, 0.3, 0.5, 0.7, 1, 2, 3, 4, 5} (with the convention that σ = 0 corresponds to the vanilla topological gradient) and learning rate λ = 0.1. We also rely on a subsampling with batch size s = 50. To quantify the variability of the scheme with respect to the randomness induced by the subsampling step, we run each gradient descent 50 times with a fixed initialization X 0 , up to a maximum of 200 iterations, stopped earlier if a loss of 0 (no topology left, global minimum has been reached) is measured.Figure 12 displays the results of this illustrative experiment. We report the median of both running time and number of iterations to reach convergence (or reach the 200 iterations limit), along with the 10 and 90 percentiles. The conclusions are:• For σ = 0 (vanilla) and σ ≥ 3, the gradient descent never converges in less than 200 steps. Since the radius of the diameter of the circle is 2, it is not surprising that taking a bandwidth larger than that hinders convergence.• For σ ∈ (0.1, 1], the convergence occurs within the same order of magnitude (between 0.49 and 1.74s), the best performance being reached at σ = 0.3 (recall that we used in the paper, testifying that we did not rely on hyperparameter tuning). It suggests that, on regular structure, the approach is smooth with respect to σ. Empirically, we observe that a good proxy is to take σ < median({|x ix j |} i,j ). Note that even though in theory, σ → 0 should recover the vanilla topological gradients, one is limited by numerical accuracy when evaluating the Gaussian kernel.• The variation around the median over 50 runs is very small: the randomness of the samples at each iteration (hence of the trajectory) barely impacts the decrease of the loss and thus the convergence time. We also studied how the bandwidth σ influences the correlation scores of Table 1 in Figure 13 . We observe oscillations for values that are roughly on the sides (very small or very large), and more stable scores for middle range values. Note that these oscillations could also come from how the correlation score itself is computed. This is not trivial, because the intermediate space D is only a metric space.In many applications, RKHS are restricted to spaces of functions valued in R or C, but the theory adapts faithfully to the more general setting of vector-valued maps.https://github.com/anigmetov/oineusNote that the loss computed with oineus has not been normalized in the figure, which is why its values are larger than the others. This has no influence over its minimum number of iterations needed to reach 0 though.Indeed, we noticed that 3. was a consistent threshold for detecting whether the representative cycle of the most persistent PD point changed between iterations.The repli-score is a real-valued proxy for the cell cycle, which was introduced in[31], and which can be computed out of the copy numbers of genome regions associated to the early phases of the cell cycle. Thus, a large correlation with this score indicates that the circular shape in the optimized latent space is indeed representative of the cell cycle itself.A simplicial complex is a combinatorial object generalizing graphs and triangulations. Given a finite set of vertices V = {v1, . . . , vn}, a finite simplicial complex K is a subset of 2 V (whose elements are called simplices) such that σ ∈ K ⇒ τ ∈ K, ∀τ ⊆ σ (if a simplex is in the complex, its faces must be in it as well).It may happen that a topological feature appears at some time t b and is never destroyed, in which case the death time is set to +∞. However, in the context of the Vietoris-Rips filtration, extensively studied in this work, this (almost) never happens. See the next paragraph.The theoretical worst case yields a cubic complexity, but the matrix that has to be reduced is typically very sparse, enabling this practical speed up.However, the statistical efficiency of Dgm(X) when it is used as an estimator for the topology of an underlying manifold M deteriorate when the intrinsic dimension of M increases. 


Decision: Topological optimization suffers from very sparse gradients which leads to slow convergence. This paper suggests a method to alleviate this difficulty by a smoothening procedure. The reviewers were all convinced by the merit and practical benefits  of the method and recommended acceptance.

Reviews: review: paper_summary: paper_summary: This paper proposes a new topological optimization scheme mainly due to the sparsity of topological gradients. The authors introduce the notion of diffeomorphic interpolation and use this to create a smooth vector field over the whole space, which gives a gradient descent algorithm. The authors show some theoretical results which guarantee that the smoothness of the interpolation can be controlled. Further, the authors show numerical experimental evidence for the quick convergence of their method. The authors show, using experiments, that this scheme can be used to interpolate gradient on the whole space by considering a gradient on a small sample yielding expected results. The authors show that this framework can be used to infuse topological information in a pre-trained autoencoder model.
strength_weakness: strength_weakness: 1. Section 3 is notation intensive. For example, Proposition 3.2, concepts like condition number are not defined earlier and hence, the expression can be moved to the appendix with the main text containing a line which says that the Lipschitz constant can be bounded which guarantees the smoothness. 

2. Proof of Proposition 3.1 can be moved to the appendix as it is not hard to see, given the way \tilde v is defined/constructed to improve the readability. 

3. For the subsampling experiment, the comparison with vanilla topological gradient seems slightly unfair because there is no way for the network to update gradient over the whole space when a small percentage of points are being sampled for gradient update. Updating gradient over the whole space is not feasible, but is there a threshold on the number of points until which the vanilla gradient update is feasible and if the final output, in that case, is comparable to the output from diffeomorphic gradient descent? 

Minor: 

In Algorithm 1, should it be “Set $X_k \coloneqq X_{k-1} - \lambda \tilde v(X_{k-1})$"?
questions: questions: Can this method work for differentiable multiparameter vectorization methods?
limitations: limitations: Yes, the authors have discussed limitations.
paper_summary: paper_summary: The paper proposes a novel way based on diffeomorphic interpolation to deal with the problem of sparse gradients when performing topological optimisation, which is relevant in the intersection between Topological Data Analysis (TDA) and ML.
strength_weakness: strength_weakness: (W1) The experimental setup could be strengthen by providing a real-world application
questions: questions: (Q1) In Figure 7, the initial LS is in blue and the final is in orange, while in Figure 6, the initial LS is in orange and the final is in orange, I find this confusing, is there a reason for it?

(Q2) Can you please elaborate more on the importance of your suggested application to black-box autoencoder models, and how does your method compare to other topological optimisation techniques such as Oineus (as in the table in this section you have only compared your approach to the Vanilla, and you reported Oineus for the rest of the experiments)

(Q3) The performance between Vanilla and Diffeo is very similar in the table in section Black-box autoencoder models (apart from Pig), have you ran your experiments multiple times and what is the mean and standard deviation? Also, do you have an explanation as to why your method outperforms Vanilla by a large margin for Pig?

(Q4) Figure 3 is difficult to interpret. What is considered to be the expected end result of topological optimisation? The final shape after 100 iterations of your algorithm is very distorted, can you please add more details on how is that good?

(Q5) What are some real-world uses of topological optimisation? Can you please elaborate on when one should care about this, and what possible applications are there

(Q6) Do you see application in the real-world of your method, and how big datasets can your approach handle?

(Q7) Can you please elaborate on how exactly your method offers better interpretability for the black-box autoencoder application?

(Q8) How robust is your method to taking different samples and to different sampling strategies?
limitations: limitations: Yes
paper_summary: paper_summary: In "Diffeomorphic interpolation for efficient persistence-based topological optimization", the authors provide a novel method to compute interpolations of gradients of persistence diagrams. This approach makes optimisation of point clouds w.r.t. loss functions defined on the associated persistence diagrams feasible, as it provides generalised gradients which smoothly vary over the underlying space. In summary, the authors use some form of Gaussian kernels and ideas from stochastic gradient descent to extend the standard gradients only defined on a few points to the entire  base space. They compare the running time of their algorithm favourably to other existing topological gradients, and perform qualitative experiments on synthetic data and the latent space of an autoencoder.
strength_weakness: strength_weakness: 1. There is no intuition provided behind theoretical result number 2. One can either just believe it, or read the very brief and technical derivation in the appendix; a brief summary and intuition in the main text would be great!
2. While there are theoretical guarantees for the derivative of the loss at a single time point (and thus in the step-size -> 0 case), any guarantees which also apply in the case for discrete steps are missing. To the reviewer's understanding, the diffeomorphic interpolation of the loss function does not need to be continuous in t, even in the case of infinitesimal small step sizes. Thus the relevancy of the theoretical guarantees in practice is not clear.
3. The choice of the bandwidth of the Gaussian kernel used for the diffeomorphic interpolation, a seemingly very important hyperparameter, is not discussed in the paper. The impact of different choices of sigma is not looked at in the experiments. (See the question section.)
4. The experiments presented are either toy examples or very artificial. An application on a down-stream task relevant in some field of practice or real-world application would significantly enhance the strength of the paper. (This is a caveat applicable to much of the TDA literature in general though.)
questions: questions: 1. What is the intuition behind the proof of Theorem 2? (Cf. Weakness 1)
2. What do the theoretical guarantees tell us in practice? Do they tell us something? Can you provide any guarantees on the evolution of the loss over a number of steps with a discrete step size? Or even in the continuous case? If not, what are the theoretical reasons these cases are so much harder?
3. How did you choose the bandwidth in practice? If someone knows nothing about their dataset, what would be a good heuristic to choose sigma? Do we need to fix a global sigma or can we vary sigma locally?
4. What is the impact of different choices of sigma on the results of the experiments? How sensitive are the experiments to varying sigma?
5. In the bottom part of Figure 4, why is the loss of Oineus significantly larger than the loss of the other methods for small x values? 
6. In Figure 3, how did you pick the hyperparameters? Did you tune them for the diffeomorphic interpolation, or for the vanilla case? It might be the case that different hyperparameter regimes are optimal.
7. In Figure 9, 3rd from the left: what are the purple points? 4th from the left: How can the large loss spikes for the diffeo algorithm be explained? How do the hyperparameter choices influence the loss stability?

I would like to thank the authors in advance for their answers!

Comments:
* Line 20: "can be done once and then re-applied to new data in linear time" -> What does that even mean?
* ll 33: "active line OR research"
* ll 37: "yielding to the problem"
* ll 84: "allows to process datasets whose sizes are currently out of reach in TDA," This claim is a bit bold: Subsampling and taking landmarks was a thing even before this paper. Maybe this applies to optimisation of point clouds in TDA:
* ll 100: "machinary"
*  section 2.2: The defintion of RKHS is really confusing. If I didn't know the definition beforehand, I would not have understood it there
* ll 181: where does v live? can X and Y be any arbitrary point clouds? With which definitions? How can X be the argument on the left side of the definition, and still be referenced as all X in the set builder notation on the RHS?
* ll 193: "i-th line"
* ll 214: "s-points"
* ll 476: The grammar of this sentence is confusing.
* page 13, footnote 5: only non-empty simplices. footnote 8: "deteriorate"
limitations: limitations: As detailed above, a discussion on hyperparameter choices, the stability of the algorithm, and the relevance of the guarantees in practice is missing. For possible negative societal impact: In Figure 5 the authors blow up a synthetic bunny and make the forgivable mistake to fail to mention that this experiment should not be performed on real-world bunnies...
paper_summary: paper_summary: The paper introduces a novel approach to topological optimization using diffeomorphic interpolation to tackle the sparsity and inefficiency of gradient descent in topological data analysis (TDA). The authors focus on optimizing point clouds by transforming sparse gradients into smooth vector fields. They achieve this by assigning exponentially decaying weights to each point with a gradient, allowing the influence to span the entire space. This approach enables a larger set of points to be optimized in each iteration, thereby accelerating convergence. The effectiveness of the method is demonstrated through experiments that showcase its fast convergence. Additionally, the method proves to be particularly efficient when combined with subsampling techniques, facilitating large-scale topological optimization.
strength_weakness: strength_weakness: The authors describe three types of topological loss where the proposed method can be applied. However, the application to topological registration losses is not demonstrated in their experiments. My concern is that, in topological registration loss, the match between the output diagram and the target diagram is often unstable. Any mismatch could result in an increased registration loss. With diffeomorphic interpolation, this increment could be even larger, potentially making the overall optimization more unstable.
questions: questions: The experiments on convergence rate are only designed for simple topological situations. When the topological complexity increases, the optimization becomes unstable. Will the proposed method introduce additional instability to the optimization? In other words, could a convergence theorem for the proposed method be derived?
limitations: limitations: See above.

Comments: comment: I thank the authors for their efforts and the elaborative response. I have read the response and would like to stick with my score.
_I believe that adding the explanations and intuition you gave me in the rebuttal would be a valuable addition to the (appendix of the) paper (...)_

_I appreciate your study on the effect of the bandwidth on the performance of your algorithm (...) to fully appreciate this I think you should use (a) the $\sigma$ computed by this heuristic in your experiments or (b) at least report the suggested $\sigma$ in comparison to the $\sigma$ used in your experiments._

Thank you for your suggestions. We will include both our theorem's discussion and highlight our heuristic and the corresponding $\sigma$ value in the experiment section. 

_I am not convinced by the additional experiments on the VAE (...) did you check whether this particular circle, and the individual location of the points on the circle correspond to any meaningful real-world features?_

This is a rightful question, and we apologize for not getting into much details about it in the rebuttal (mostly due to lack of space). In order to make sure that the circular shape that we end up with (after topological optimization) is biologically relevant and not artifactual, we measured the correlation between the latent angles (computed with the same formula than the one in the main paper), and a quantity called the **repli-score**, which is computed out of the copy numbers of genome regions associated to the early phases of the cell cycle. It was introduced and provided in [_Cell-cycle dynamics of chromosomal organization at single-cell resolution_. Nagano et al. **Nature**, 2017] as an efficient, computable, and real-valued proxy for the cell cycle. Thus, a large correlation indicates that the circular shape in the optimized latent space is indeed representative of the cell cycle itself.  

_What would you need to do to generalise your topological optimisation to work with DTM filtrations? (...)_

Thank you for this suggestion. As far as our method is concerned, one only requires a sparse gradient on a point cloud in order to run. Hence, the main question is how to compute such gradients for DTM-based filtrations. For instance, DTM-Rips filtrations should be possible, as the equation for computing the DTM is easy to differentiate, and the gradient for Rips is already known, but DTM-Alpha filtrations should be more difficult, as it is less clear what the gradient of the Alpha filtration is. See also our answer to the last question of reviewer fiZD in the context of multiparameter persistence (which can also be used for dealing with noisy data).
Thank you very much for your detailed rebuttal! I believe your changes and proposed changes to the manuscript are very valuable. I still have some brief comments:

1. I believe that adding the explanations and intuition you gave me in the rebuttal would be a valuable addition to the (appendix of the) paper. For example, sometimes it is very relevant and interesting information why a theorem takes the form it has, and why it is hard to prove something more difficult.
2. I appreciate your study on the effect of the bandwidth $\sigma$ on the performance of your algorithm. I also value the heuristic for choosing $\sigma$. However, to fully appreciate this I think you should use __(a)__ the $\sigma$ computed by this heuristic in your experiments or __(b)__ at least report the suggested $\sigma$ in comparison to the $\sigma$ used in your experiments.
3. I am not convinced by the additional experiments on the VAE. (Or I might have misunderstood.) I know that you have added a loss-term to optimise for the embedding of the VAE having a circular structure. However, did you check whether this particular circle, and the individual location of the points on the circle correspond to any meaningful real-world features? I know that, because of the setup of the experiments and the studied cell-cycle, it is not unreasonable to expect some circle to form. I just don't understand how to quantitatively evaluate the meaningfulness of your result. (Maybe using circular coordinates, for example from [1])
4. I know that this is a complicated question: What would you need to do to generalise your topological optimisation to work with DTM filtrations? It seems that many topological features in noisy real-world settings are much more likely to be uncovered by DTM filtrations, for example in the scihc embedding. [2]

Despite all the above remarks, I believe that this is a very cool paper even in its current state. I am looking forward to exciting applications of this in the world of ML. I will thus be raising my score and recommending the paper for acceptance

[1] Perea, Jose A. "Sparse circular coordinates via principal ℤ-bundles." Topological Data Analysis: The Abel Symposium 2018. Cham: Springer International Publishing, 2020.
[2] Anai, H., Chazal, F., Glisse, M., Ike, Y., Inakoshi, H., Tinarrage, R., & Umeda, Y. (2020). DTM-based filtrations. In Topological Data Analysis: The Abel Symposium 2018 (pp. 33-66). Springer International Publishing.
Dear reviewer v7wD,

Thank you for your review. The authors have tried to address your concerns in the rebuttal. Please carefully read their rebuttal, and let them know what you think, and whether there is any more clarifications you require. Note that author-reviewer discussion ends on August 13th.
Thanks!
the AC
Thank you for the response. It’s encouraging to hear that global convergence has been observed in experiments. I look forward to seeing a convergence theorem for diffeomorphic interpolation. Overall, this is a strong paper.
I thank the authors for their efforts and the response. I have read the response and would like to stick with my score.