Diffeomorphic interpolation for efficient persistence-based topological optimization
Authors: Mathieu Carrière
Marc Theveneau
Théo Lacombe
Venue: NeurIPS 2024 poster (2024)
Keywords: Persistent Homology
Persistence Diagrams
Optimization
Topological Data Analysis
Primary area: optimization
TL;DR: We propose a diffeomorphic interpolation of (typically sparse) gradients appearing in Topological Data Analysis, yielding substantially faster and smoother optimization schemes.

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.

Diffeomorphic interpolation for efficient
persistence-based topological optimization
Mathieu Carrière
Mare Theveneau"
DataShape
Shape Analysis Group
Centre Inria d' Université Côte d' Azur Sophia Antipolis, France mathieu. www.iereeiniiii..rr.....
Computer Science department, McGill Montréal, Quebec, Canada marc. twwveheaummaillll Lateverneau@maa..mgglll......
Théo Lacombe
Laboratoire d' d'Informatique Gaspard Monge, Univ. Gustave Eiffel, CNRS, LIGM, F-77454
Marne-la- France theo. laaaoolacccmbeemmyyeleffe.....
Abstract
Topological Data Analysis (TDA) provides pipeline extract quantitative topo- logical descriptors from structured objects. This enables the definition of topo- logical loss functions, which assert what extent a given object exhibits some topological properties. These losses can then be used to perform topological op- timization 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 on only very few coordinates of the input object, yielding dramatically slow optimization schemes in practice.
Focusing on the central case of topological optimization for point clouds, pro- pose in this work to overcome this limitation using diffeomorphic interpolation. turning sparse gradients into smooth vector fields defined on the whole space, with quantifiable Lipschitz constants. In particular, we show that our approach combines efficiently with subsampling techniques routinely used in TDA. the diffeomorphism derived from the gradient computed on a subsample can used update the coordinates of the full input object, allowing to perform topo- logical optimization on point clouds at an unprecedented scale. Finally, we also showease the relevance of our approach for black-box autoencoder (AE) regular- ization, where we aim at enforcing topological priors on the latent spaces associ- ated fixed, pre-trained, black-box AE models, and where we show that learning a diffcomorphic flow can be done once and then re-applied to new data in linear time (while vanilla topological optimization has be re-run from scratch). More- over, reverting the flow allows to generate data by sampling the topologically- optimized latent space directly, yielding better interpretability of the model.
Introduction
Persistent homology (PH) is central tool of Topological Data Analysis (TDA) that enables the extraction quantitative topological information (such as, e.g., the number and sizes of loops, con- nected components, branches, cavities, etc) about structured objects (such graphs, times series or
Part this work was done when MT was doing an intermship at the Laboratoire d'Informatique Gaspard Monge and student at Université Paris-Saclay.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

point clouds sampled from, e.g., submanifolds), summarized in compact descriptors called persis- tence 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 few.
Another active line or research the crossroad of TDA and ML is (persistence-basse -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 E IRxxd that Dgm(X (X) becomes as close as possible (w.r.t an appropriate metric denoted by W) to some target diagram Diarget. hence yielding to the problem of minimizing X -> W (Dgmm(X (x). Ptarget): ). This idea has then been revisited with different flavors, for instance by adding topology-base 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 function
L ((Dgm(XX) R,
where l is user-chosen loss function that quantifies to what extent Dgm(X) reflects some pre- scribed topological properties inferred from X. Under mild assumptions (see Section 2), the map L is differentiable generically and its gradients are obtained as byproduct of the computation of Dgm(X). ()). However, these approaches are limited in practice by two major issues: (i) the com- putation of X + Dgm () scales poorly with the size of X (e.g., number of points in a point cloud, number of nodes in a graph, etc), and (ii) the gradient (T(X) tends to be very sparse: if X = (Ii,..., In) E "x is a point cloud, VX()); + 0 for only very few indices i E (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 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 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 used 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 cannol 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).
Contributions and Outline. In this article, we propose to replace the standard gradient VL(X) of (??) derived formally by diffeomorphism v : Rd + Rd that interpolates V VL(X) on its non-zero entries, that is v(z;) = VL()); for all i E I = (j LX))), # 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 ex- tends this gradient to smooth and denser vector field defined on the whole space Rdd, and which is able to move a lot more points in X at each iteration,
2

---

Figure I: Illustration of the Vietoris- filtration point cloud in B focusing on one-dimensional topo- logical features (loops). When the filtration parameter increases, loops appear and disappcar the filtration. These values are accounted in the resulting persistence diagram (right).
We prove that its updates indeed decrease topological losses, and we quantify its smooth- ness by upper bounding its Lipschitz constant (again, the context of topological losses),
We showcase its practical efficiency: we show that it compares favorably the main base- line [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 can success- fully be used for the tasks mentioned above concerning black-box autoencoder models.
Section 2 provides necessary background in Topological Data Analysis and diffeomorphic interpo- lations. Section 3 presents our approach and its corresponding guarantees, and Section 4 showcases our experiments. Limitations and further research directions are discussed in Section
Background
2.1 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 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 (#1...., xn) E R The Vietoris- Rips filtration consists of building an increasing sequence of simplicial complexes (K)) over X by inserting a simplex o = (Ii, Tsp) when- ever Yj,j' E (1....,p). ||:;; Ii,y < I. Each time a simplex a 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 loop, that is a 1-dimensional feature, making it topologically trivial). The persistent homology machinary tracks the apparition and de- struction of such features in the so-called persistence diagram (PD) of X, denoted by Dgm(X), Therefore, Dgm(X] (X) is a set of points in R² of the form (4,, with La 2 l.. where each such point accounts for the presence of topological feature inferred from X that appeared time I, following the insertion of an edge (Ma,, Tt2) with || Figll = 14 and disappeared at time la following the insertion of an edge (Iig, 24,)) with ||:r Iz, = td. Figure illustrates this construction. From computational standpoint, computing the VR diagram of X E R that would reflect topolog- ical features of dimension d' < d runs in (( +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) 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 form
L: E [a ((Dgm(X)) E R.
(I)
Here, É 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 - Dgin (X) is stable (Lipschitz continuous) [17, 18, 38). Therefore, if is Lipschitz continuous, so is L., hence it admits gradient almost everywhere by Rademacher theorem. Building on these theoretical statements, one
3

---

can consider the "vanilla" gradient descent update Xx+1 == Xx - XVL(X)) for given learning- rate À > 0 and iterate it 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 local minimum of I.) is proved.
From a computational perspective, deriving (L(X) comes in two steps. Let == Dgm(X) (X), written as μ {b,,, d i E } for some finite set of indices 1. To each E I correspond four (possibly coinciding) points Tig Zizz -ig. Ti, in the input point cloud X. Intuitively, minimizing μ + {(u) boils down to prescribe descent direction (ōb, E IR² each (b,, (b. for i I, where ōb, at (p) and ôd, ad: (n)) Backpropagating this perturbation X will move the corresponding 06,
86
points Tiy Fig Fig Fi, in order to increase or decrease the distances |,,, - Fial ll = b, and || :ial di accordingly. This yields to the following formula:
aT. dr
ae
bb. ae
0d,
(X)
Σ
+
(X),
(2)
bb.
Or
0d, Ox
x d.)
where the notation x -> (b;, d,) means that x E X appears in (at least) one of the four points yielding the presence of (b;, d,) in the diagram f (xm(().. A fundamental contribution of [28, $3.3] is to prove that the chain rule formula (2) is indeed valid?. Most of the time, a point x E X will not belong to any critical pair (0%, σa) and the above gradient coordinate is O. Therefore, the gradient of L depends on very few points of X, yielding the sparsity phenomenon discussed in Section 1.
(..,.... ..)) E Raxd be point cloud and Dgm(X) X) = [(h;, d) i E /} be its PD. There are several natural losses that have been introduced in the TDA literature:
Examples of common topological losses. Let X
Topological simplification losses: typically of the form Σiej(b; - d,)², where is I. Such losses push (some of the) points in Dgm(X) toward the diagonal A = {b = d),
hence destroying the corresponding topological features appearing in X.
Topological augmentation losses [4]: similar to simplification losses, but typically attempt- ing to push points in Dgm( (X) away from A, i.e., minimizing - ib; - -d,²² 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 Dearget: one minimizes W (Dgmm(X) Drarget) Dearge) where W denotes standard metric between persistence diagrams.
This loss attempts to modify X so that exhibits a prescribed topological structure.
2.2 Diffeomorphic interpolations
In order to overcome the sparsity of gradients appearing in TDA, we rely on diffeomorphic interpo- lations (see, e.g., [43, Chapter 8J). Let X (.,,,..,...)) of indices on which (XX) is non-zero and let a, = (VL(X)); Rd for i E I. Our goal is to find smooth vector field 6 IRd + Rd such that, for all i E T, ((x;) = a, To formalize this, we consider a Hilbert space H C (R=))** for which the map 8 f + (a. fz)) g. a f(x) is continuous for any (a, .2) E Rd Rd. Such space is called a Reproducing Kernel Hilbert Space (RKHS)². A cru- cial property is that there exists a matrix-valued kernel operator K Rdd xR + Idd whose outputs are symmetric and positive definite, and related to H through the relation (kg k.g) = a" K (z, y)B
Ra let I S (I,...,,) nn denote the set
for all x. y, a, E Rd, where R. E H is the unique vector provided by the Riesz representation theorem such that (ke, f) = (a, f()). Conversely, any such kernel induces a (unique) RKHS Hl (of which K is the reproducing kernel). Now, we can consider the following problem:
minimize [|e]] s.t. v(z,) ag, Vi E I,
(3)
that is, we are secking for the smoothest (lowest norm) element of И that solves our interpolation problem. The solution ì of this problem is the projection of 0 onto the affine set {u E H v(x;)) =
Thii not trivial. because the intermediate space only a metric space. 'In many applications. RKHS are restricted to spaces of functions valued in IR or C. but the theory adapts faithfully to the more general setting of vector-valued maps.

---

a, Vi E I]. Observe that ë belongs to the orthogonal of {or E H v(zz)) = 0, Vi E I), and thus of {u E H| (k2, v) = 0, Vi E 1, 04 E Rd), and therefore E span( (1x2 i E T)). This justifies to search for ē in the form of ((x) = Eier K (a, z;)as, and the interpolation that it must satisfy yields 0(x) = Σ (x, [K a)), where IK is the block matrix See also [43, Theorem 8.8]. In particular,
2y)),ye and a = (a:)vez- important to note that inherits from the regularity of K and will typically be diffeomorphism in this work. If K is the Gaussian kernel defined by
K (x. == exp
| Ia for some bandwidth o > 0, a choice to which we stick to in the rest 2
of this work, the expression of ō reduces
5(x) = Σ Po (11= xz(l)as: iE7
(4)
2 20 and a; (K- with K (po( - 23 2zll) Note that
where Pa
can be understood the convolution of with Gaussian kernel, but involving correction K- guaranteeing that after the convolution, the interpolation constraint is satisfied. We will call the diffeomorphic interpolation associated to the vectors and indices {a; I}.
3
Diffeomorphic interpolation of the vanilla topological gradient
3.1 Methodology
We aim at minimizing a loss function L X (DDg(XX) as (I), starting from some initialization Xo. and assuming that L is lower bounded (lypically by 0) and locally semi-convex. This assumption typically satisfied by the topological losses l introduced in Section 2.1. Gradient descents implemented in practice are (explicit) discretization of the gradient flow
dX dt E (t)( X(0) = Xo.
(5)
where (T.(X)= = [wl ((Y) 2 T(X)+ X) for all X, denotes the subdifferential of at X. Note that a topological loss I. is typically not differentiable every- where, since the map X + Dgm(X) is differentiable almost everywhere but not in CL. However, uniqueness of the gra- is In this work, we propose to use the dynamic described by the diffeomorphism in introduced in (4) interpolating the current vanilla topological gradient VI((X,) each time i, formally considering solutions X of
+ v. ()-
dient flow on a maximal interval [0, + is guaranteed lower bounded and locally semi-convex [15, $B.1].
Figure (blue) A point cloud X, and (black) the negative gradi- ent -VL(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 (T.(X) only affects four points in X, the diffcomorphic in- terpolation ō(X) (orange, = 0.1) is defined on Re hence extends smoothly to other points in X.
dX dt
i. (X((), X())= Xn-
(6)
Here, slightly overloading notation, 8e (X)) denotes the n nxd matrix where the i-th line is given by Tz (X()))) The flow at time T associated to (6) is the map
ΦT : Do > 3o L. 8 (=()))dt
*(t) - (z(()), x00) = To
which can inverted by simply following the flow backward (i.e., by following i instead of -G)) We now guarantee that at each time following ün instead of the vanilla topological gradient VL(X,) still provides descent direction for the topological loss L.
Proposition 3.1. For each 2 0, it holds that d
(VL(( So.
de
Proof. One has (L((())) = - ( E (o, Since de
VL(((t)), 0 for i & I. and (x)())) - -VL( for i I, the result follows.
D
5

---

Moreover, it is also possible to upper bound the smoothness, i.c., the Lipschitz constant, of i:
Proposition 3.2. Let L be the simplification or augmentation loss computed with k [7] PD points,
as defined the end of Section 2. Let û = Oy be the diffeomorphic interpolation associated to the vanilla topological gradient at time 1. 2 0. Then, one has, Vx. y E TR. and 1. 2 0:
[0(x)) - 81y)l2 10(x)) +(v) la I Ca ad-1_ k()) Porss Perssx(Denn (Dgm (X())) || -yl2,
where Ca = add23+ T d-1 the sum of the k largest distances to the diagonal in Dgm(X(t)) (t)).
k(K is the condition number of K, and Perss (Dgm( (f)) is
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 Lip- schitz constant depends on the kernel bandwidth a: indeed, the more spread the Gaussian function is, the more critical points can influence other data points potentially far from them, introducing un- wanted 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 larger Lipschitz constant as well. Finally, the total persistence also appears, as the more PD points one has to opti- mize, 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.
3.2 Subsampling techniques to scale topological optimization
As consequence of the limited scaling of the Victoris-Rip 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 < n, yielding a smaller point cloud X' X. Provided that the Hausdorff distance between X' and X 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 fur- ther 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 7 (of the vanilla topological gradient) computed on the subsample X' still provides vector field defined on the whole input space IRd, in particular on each point of X and the update can then be performed in linear time with respect n. This yields Algorithm 1. Figure illustrates the qualitative benefits offered by the joint use of subsampling and diffeomorphic interpolations when compared to vanilla topological gradients. larger-scale experiment provided in Section 4.
Algorithm 1 Diffeomorphic gradient descent for topological loss functions with subsampling
Input: Initial X, E Randd loss function e, learning rate À > 0, subsampling size s E nn......).
max. epoch T 2 1, stopping criterion.
Set L : X + {(Dmm(X) (+ possibly a regularization term in X).
for k 1.....7 do
Subsample X (x', (x',....,=') uniformly from XR-1-
Compute VL(X (vanilla topological gradient)
Compute the diffeomorphic interpolation o(X)_,) from VL(X;_)) using (4).
Set X Xa-1 A6(Xx-1)
if stopping criterion is reached then Return X
end if end for Return XT
Stopping criterion. A natural stopping criterion for Algorithm is to assess whether the loss L(X;) {(ggm(X() is smaller than some > O. However, computing Dgm(X,) can be in- tractable if X, is large. Therefore, tractable loss to consider is i(x.) = E[{(Dgm(X[)), where

---

(a) L =
(b) = 100 (Vanilla)
(c) = 500 (Vanilla)
(d) 100 (Diffeo)
Figure 3: Showcase of the usefulness of subsampling combined with diffeomorphic interpolations to minimize topological simplification loss, with parameters À = 0.1, s 50, n 500. (a) Initial point cloud X (blue), subsample X' (red). vanilla topological gradient on the subsample (black) and corresponding diffeomorphic interpolation (orange). (b) and (c), the point cloud X, after running
100 and 500 steps vanilla gradient descent. (d) the point cloud X₂ after running = 100 steps of diffcomorphic gradient descent.
X a uniform s-sample from X₁. Under that perspective, Algorithm can be re-interpreted as a kind of stochastic gradient descent on i. for which two standard stopping criteria can be used: (a) compute an exponential moving average of the loss individual samples X; over iterations, or (b) compute validation loss, i.e., sample X ".. XX,,,, and estimate Ż by E e(Dgm(Xi,.e)
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 better objective X,) overall, and thus stick to this choice in our experiments.
Numerical experiments
We provide numerical evidence for the strength of diffeomorphic interpolations. PH- related computations relies on the library Gudhi [40] and automatic differentiation relies on tensorflow [1]. The "big-step gradient" baseline [34] implementation based on oineus". The first two experiments were run on a 11th Gen Intel (R) Core (TM) i5-1135G7 @ 2. 20GGHZ
the last one on a 2x Xeon SP Gold 5115 @ 2. 2.40GHz. Our code is publicly available https: //gi github.co. com com/tlacombe/ton
Convergence speed and running times. We sample uniformly n = 200 points on unit circle in I with some additional Gaussian noise, and then minimize the simplification loss L X + D(bd))Dgm(X) |d]?, 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 loss of is reached (Dgm(X) X) is empty), with a same learning rate À = 0.1. The bandwidth of the Gaussian kernel in (4) is set to a = 0.1. We include the competitor oineus [34], 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 4s, 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 diffcomorphic 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 subsam- pling routines (Algorithm 1) allows to perform topological optimization on point clouds with thou- sands of points, new scale in the field. For this, consider the vertices of the Stanford Bunny [41], yielding a point cloud X, E Raxd with n = 35, 947 and d = 3. We consider a topolog-
'nttps://gii http:/gg.tt..ooooooo.oooooocoo.t wtt..//aaaantoo//oiinnoeee
'Note 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 influence over its minimum number of iterations needed reach 0 though.
7

---

Init.
Vanilla output
Diffeo output
Oineus output
Oineus+ Diffeo outout
35
33
3.8
-30
00
16
00
21
71-
13
11-
15
15-
Figure 4: (Top) From left to right initial point cloud, and final point cloud for the different flows. (Bottom) Evolution of the losses with respect to the number of iterations and with respect to running time.
Evalution over iterations
Dilleo
Epoch
Figure 5: From left to right: initial Stanford bunny Xo. the point cloud after 1, 000 epachs of vanilla topological gradient descent (barely any changes), the point cloud after 200 epochs of diffeomorphic gradient descent, after 1,000 epochs, and eventually the evolution of losses for both methods over iterations.
ical augmentation loss (see Section 2.1) for two-dimensional topological features, i.c., we aim at increasing the persistence of the bunny's 's cavity. The size of n makes the computation of Dgm (Xo) untractable (recall that scales in O(n")) )); we thus rely on subsampling with sample size 8 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 de- creases the loss) in 1,000 epochs. In sharp contrast, as our diffeomorphic interpolation computed on subsamples is defined on R$, it updates the whole point cloud each iteration, making possible to decrease the objective function where the vanilla gradient descent is completely stuck. Note that step of diffcomorphic interpolation. in that case, takes about 10 times longer than vanilla step. An additional subsampling experiment can be found in Appendix C.
Black-box autoencoder models. Next, we apply our diffcomorphic interpolations to black-box autoencoder models. In their simplest formulation, autoencoders (AE) can be summarized as two maps E Rd + Rd and D Rd + Rd called encoder and decoder respectively. The intermediate space TRd" in which the encoder E is valued is referred to as latent space (LSS), with typically d' & d. In general, without further care, there no reason to expect that the L.S of a point cloud X, E(X) = (E(x1),.... E(xx)), 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 black-box, pre- trained AE. However, replacing (E, D) by (fp o E, D o ip ") for any invertible map (p 12 + Rd yields an AE producing the same outputs yet changing the LS E(X). without explicit access to the AE's model. Hence, propose to learn such & with diffeomorphic interpolations: given some latent space EXX), we apply T steps of our diffeomorphic gradient descent algorithm to X +> {(ggm(X) initialized at E(X). We thus get a sequence of smooth displacements -U,...., -Ū of Rdd that discretizes the flow (7) via yp : 2o + 2o EE- 0R (x where Ty --t-1 - (m)- and such that (g is more topologically satisfying. This fixed diffeomorphism &o can
8

---

Figure 7: COIL images, their corresponding initial LSs in blue and final LSs obtained with diffeomorphic gradient descent in orange, and the corresponding topological losses, for both vase (left) and duck (right).
then be re-applied to any new data coming out of the encoder in a deterministic way. Moreover, any random sample from the aoplooicalll-ootimmzze by following Ur, UT-1, 01, 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].
LS can be inverted without further computations
In order to illustrate these properties, we trained variational autoencoder (VAE) to project family of datasets of images representing rotating objects, named COIL [32], to two-dimensional latent spaces. Given that every dataset 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 shallow (the encoder has one fully- connected layer (100 neurons), and the decoder has two (50 and 100 neurons), all layers use ReL.u activations), and thus the learned latent spaces, although still looking like curves thanks to continuity, do not necessarily display circular patterss. This makes generating new data more difficult, as latent spaces are harder to interpret. To improve on this, we learn a flow yp as described above with an augmentation loss associated to the I-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, Gaussian kernels with bandwidth o = 0.3 and an increase of at least 3. in the topological loss (from an iteration to the next) stop the algorithm
do not subsampling, and we use learning rate À 0.1,
Ostreand HelLL
We provide some qualitative results in Figure 7 (see also Appendix C, Figure 10). In order to quantify the improve- ment, we also computed the Pearson correlation scores between the ground truth angles 0g and the angles 6, com- puted from the wwwpllooicalllooppmmmmieo
LSs with
0. (x= Blipo E(X)), so E E(x1) F((o))()))
where Blpo E(X) == n Σm= 8o E(x;), and y de- notes our flow. In Table 1, we provide an average of these scores over 100 test sets obtained by randomly per- turbing the training set with uniform noise of amplitude 0.05. As expected, correlation becomes better after forc- ing 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 popologcalll-oopimmiee
Figure 6: As the four samples from the dapologically-optimzzz LS (blue) are far from the initial (orange), the decoded im- ages are fuzzy. However, reverting yo and following the corresponding green trajecto- ries allows render good-looking images.
LS, then pushed to the initial LS of the black-box box VAE by following the reverted flow of our learned diffcomorphism so, and finally decoded back into realistic, COIL-like images.
Single-cell data. We also deployed our proposed topological correction of I.s from AEs on a real- world dataset of single cells. Specifically, 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.
'Indeed, we noticed that 3. was consistent threshold for detecting whether the representative cycle of the most persistent PD point changed between iterations.

---

Dataset
Duck
Cat
Pig
Vase
Тeapot
No optim.
0.56 1 7.5e-04
0.78 1 1.4e-03
0.17 1 5.7e-04
0.86 7.2e-03
0.32 1 1.6e-03
Diffeo
0.61 + 3.1e-03
0.83 1.2e-03
0.76 + 2.1e-04
0.93 + 9.8e-04
0.39 + 3.4e-03
Dataset
scHiC (augmentation)
scHiC (registration)
No optim. Diffco
0.79 8.
0.792 8.
0.84 + 4.3e-03
0.794 + 8.4e-03
Table 1: Means and variances of correlation scores computed over 100 test sets, for both COIL and scHiC.
Specifically, we processed this single cell dataset of 1, 171 cells with the sratum-adjjusted correlation coefficient (SCC) with 500kb and convolution parameter h = on chromosome 10. Then, we run kernel PCA on the SCC matrix to obtaín a preprocessed dataset in IR 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:
: Rnxd 1 W2 (Dgm' (X), Dtarget),
where Dgm(X) contains the points of the PD of X with distance-to-de to-daaoonl at least 1, Dtrrget is target PD with only one point p" -3.5, 3.5). and W is the 2- 2-Wasserstein distance between PDs. We used a = 0.2 for the augmentation loss and smaller value for the registration loss in order to mitigate the effects of matching instability), À = 0.1 on 500 epochs, with subsampling of size = 300 for computational efficiency (as computing VR diagrams without radius thresholding on 171 points already takes few minutes 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 scores' in Table 1, in which improvements can be observed. An additional experiment the influence of the bandwidth parameter a over these correlation scores, as well as over convergence, also be found in Appendix C.
0.025 for the registration loss (a is set to a
5 Conclusion
In this article, we have presented a way to turn sparse topological gradients into dense diffeomor- phisms 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 in- terpolations computed on dataset and its subsamples. This requires some control over the loca- tions 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 fixed kernel and parameters at every iteration), and understand the convergence properties of our proposed diffeomorphic gradient descent. Finally, applying our dif- feomorphic interpolation to sparse gradients computed with multiparameter persistent homology 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 the figure below: as the initial LSS (colored with the ground-truth angles) have zero (left) or two (right) loops, is impossible to unfold them with diffeomor- phisms; instead the optimized latent spaces either also exhibit no topol- ogy, or mixes different angles. future work, we plan on investigating other losses gradient descent schemes for diffeomorphic topological optimization. including strat- ified procedures similar to [27] that allow for local topological changes during training.
"The repli- spore a real-valued proxy for cell cycle, which was introduced [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.
10

---

Acknowledgements
M.C. supported by ANR grant "TopModel", ANR-23-CE23-0014, The authors are grateful to the OPAL infrastructure from Université Côte Azur for providing resources and support.
References
[1] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Good- fellow, Andrew Harp, Geoffrey Irving, Michael Isard. Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg. Dandelion Mané, Rajat Monga. Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, llya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Lcrge-scale machine learning on heterogeneous systems, 2015. Software avail- able from tensorflow.org.
[2] Henry Adams, Tegan Emerson, Michael Kirby, Rachel Neville, Chris Peterson, Patrick Ship- man, Sofya Chepushtanova, Eric Hanson, Francis Motta, and Lori Ziegelmeier. Persistence images: stable vector representation of persistent homology. Journal of Machine Learning Research, 18(8):1335, 2017.
[3] Tolga Birdal, Aaron Lou, I Leonidas Guibas, and Umut Simsekli. Intrinsic dimension, persistent homology and generalization in neural networks. In Advances Neural Information Process- ing Systems 34 (NeurIPS 2021), volume 34, pages 6779. Curran Associates, Inc., 2021.
[4] Mathieu Carrière, Frédéric Chazal, Mare Glisse, Yuichi Ike, Hariprasad Kannan, and Yuhei Umeda. Optimizing persistent homology based functions. In 38th International Conference on Machine Learning (ICML. 202/), volume 139, pages 1294-1303. PMLR, 2021.
[5] Mathieu Carrière, Frédéric Chazal, Yuichi Ike, Théo Lacombe, Martin Royer, and Yuhei Umeda. PersLay: a neural network layer for persistence diagrams and new graph topological signatures. 23rd International Conference on Artificial Intelligence and Statistics (AISTATS 2020), pages 2786-2796. PMLR, 2020.
[6] Mathieu Carrière, Steve Oudot, and Maks Ovsjanikov. Stable topological signatures for points on 3d shapes. Computer Graphics Forum, 34(5):1-12, 2015.
[7] Frédéric Chazal, David Cohen-Steiner, Leonidas Guibas, Facundo Mémoli, and Steve Oudot. Hrusdo-Ha stable signatures for shapes using persistence. Computer Graphies Forum, 28(5):133311403 2009.
[8] Frédéric Chazal, Vin De Silva, and Steve Oudot. Persistence stability for geometric complexes. Geometriae Dedicata, 173(1):193-21 191-214, 2014.
[9] Frédéric Chazal, Brittany Fasy, Fabrizio Lecci, Bertrand Michel, Alessandro Rinaldo, and Larry Wasserman. Subsampling methods for persistent homology. In 32nd International Con- ference on Machine Learning (ICML 2015), volume 37, pages 2143-2151. PMLR, 2015.
[10] Frédéric Chazal, Brittany Fasy, Fabrizio Lecci, Bertrand Michel, Alessandro Rinaldo, and Larry Wasserman. Robust topological inference: distance to a measure and kernel distance. Journal of Machine Learning Research, 18(159):1-40) 2018.
[11] Frédéric Chazal, Brittany Fasy, Fabrizio Lecci, Alessandro Rinaldo, Aarti Singh, and Larry Wasserman. On the bootstrap for persistence diagrams and landscapes. Modelirovanie Analiz Informatsionnykk Sistem, (0006))1111120 2013.
[12] Frédéric Chazal, Brittany Fasy, Fabrizio Lecci, Alessandro Rinaldo, and Larry Wasserman. Stochastic convergence of persistence landscapes and silhouettes. Journal of Computational Geometry, 662):140-161, 2015.
[13] Frédéric Chazal and Steve Yann Oudot. Towards persistence-based reconstruction in euclidean spaces. In Proceedings of the twenty-fourth annual symposiun on Computational geometry, pages 232-241. ACM, 2008
[14] Chao Chen, Xiuyan Ni, Qinxun Bai, and Yusu Wang. A topological regularizer for classifiers via persistent homology. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2573-2582, PMLR, 2019.

---

[15] Lenaic Chizat and Francis Bach. On the global convergence of gradient descent for over- parameterized models using optimal transport. Advances in neural information processing systems, 31, 2018.
[16] James R Clough, Nicholas Byrne IIkay Oksuz, Veronika A Zimmer, Julia A Schnabel, and Andrew P King. A topological loss function for deep-learning based image segmentation using persistent homology. IEEE transactions pattern analysis and machine intelligence, 44( 44111)88766-877 2020.
[17] David Cohen- Herbert Edelsbrunner, and John Harer. Stability of persistence diaggams. Discrete & Compatational Geometry, 37(1):103-120, 2007.
[18] David Cohen-Steiner, Herbert Edelsbrunner, John Harer, and Yuriy Mileyko. Lipschitz func- tions have p-stable persistence. Foundations of computational mathematics, 10(2):12--139 2010.
[19] Meryll Dindin, Yuhei Umeda, and Frederic Chazal. Topological data analysis for arrbyth- mia detection through modular neural networks. In Advances in Artificial Intelligence: 33rd Canadian Conference on Artificial Intelligence, Canadian Al 2020. Ottawa, ON, Canada, May
13-15, 2020, Proceedings 33, pages 177-188. Springer, 2020.
[20] Herbert Edelsbrunner and John Harer. Computational topology: an introduction. American Mathematical Soc., 2010.
[21] Marcio Gameiro, Yasuaki Hiraoka, and Ippei Obayashi. Continuation of point clouds via persistence diagrams. Physica D: Nonlinear Phenomena, 334: 18-132, 2016.
[22] Thomas Gebhart and Paul Schrater. Adversary detection in neural networks via persistent homology. arXiv preprint 17Xiv:1711.100555 2017.
[23] Christoph Hofer, Florian Graf, Bastian Rieck, Marc Niethammer, and Roland Kwitt. Graph fil- tration learning. In 37th International Conference on Machine Learning (C20), volume 119, pages 4314-4323. PMLR, 2020.
[24] Max Horn, Edward de Brouwer, Michael Moor, Bastian Rieck, and Karsten Borgwardt. Topo- logical Graph Neural Networks. In 10th International Conference on Learning Representations (ICLR 2022). OpenReviews.net, 2022.
[25] Xiaoling Hu, Li Fuxin, Dimitris Samaras, and Chao Chen. Topology-preserving deep image segmentation. In Advances in Neural Information Processing Systems 32 (NeurIPS 2019). pages 5657-5668. Curran Associates, Inc., 2019.
[26] Théo Lacombe, Yuichi Ike, Mathieu Carrière, Frédéric Chazal, Mare Glisse, and Yuhei Umeda Topological Uncertainty: monitoring trained neural networks through persistence of activation graphs. In 30th International Joint Conference on Artificial Intelligence (IJCAI 2021), pages 2666-2672. International Joint Conferences on Artificial Intelligence Organization, 2021.
[27] Jacob Leygonie, Mathicu Carrière, Théo Lacombe, and Steve Oudot. A gradient sampling algorithm for stratified maps with applications to topological data analysis. Mathematical Programming, pages 1-41, 2023.
[28] Jacob Leygonie, Steve Oudot, and Ulrike Tillmann. A framework for differential calculus on persistence barcodes. Foundations of Computational Mathematics, pages 1-63, 2021.
[29] Michael Moor, Max Horn, Bastian Rieck, and Karsten Borgwardt. Topological autoencoders. In 37th International Conference on Machine Learning (ICML 2020), volume 119, pages 7045-7054. PMLR, 2020.
[30] Soham Mukherjee, Shreyas Samaga, Cheng Xin, Steve Oudot, and Tamal Dey. D-gril: End- to-end topological learning with 2-parameter persistence. In CoRR. arXiv:2406.07 100, 2024.
[31] Takashi Nagano, Yaniv Lubling, Csilla Várnai, Carmel Dudley, Wing Leung, Yael Baran, Netta Mendelson-Cohen, Steven Wingett, Peter Fraser, and Amos Tanay. Cell-cycle dynamics of chromosomal organization at single-cell resolution. Nature, 547:61-67, 2017.
[32] S. A. Nene, S. K. Nayar, and H. Murase. Columbia Object Image Library (COII.-100). In Technical Report CUCS-005-96, 1996.
[33] Arnur Nigmetov, Aditi S Krishnapriyan, Nicole Sanderson, and Dmitriy Morozov. Topological regularization via persistence-sensitii optimization. arXiv preprint 27Xiv:2011.05290 2020.
12

---

[34] Arnur Nigmetov and Dmitrly Morozov. Topological optimization with big steps. Discrete & Computational Geometry, pages 1-35, 2024.
[35] Steve Y Oudot. Persistence theory: from quiver represeniations data analysis, volume 209. American Mathematical Society, 2015.
[36] Adrien Poulenard, Primoz Skraba, and Maks Ovsjanikov. Topological function optimization for continuous shape matching, In Computer Graphics Forum, volume 37, pages 13-25, Wiley Online Library, 2018
[37] Luis Scoccola, Siddharth Setlur, David Loiseaux, Mathieu Carrière, and Steve Oudot. Dif- ferentiability and convergence of filtration learning with multiparameter persistence. In 41st International Conference on Machine Learning (ICMI. 2024). PMLR, 2024.
[38] Primoz Skraba and Katharine Turner. Wasserstein stability for persistence diagrams. arXiv
preprint arXiv:2006.16824 2020.
[39] Yitzchak Solomon, Alexander Wagner, and Paul Bendich. A fast and robust method for global topological functional optimization. In International Conference on Artificial Intelligence and Statistics, pages 109-117. PMLR, 2021.
[40] The GUDHI Project. GUDHI User and Reference Manual. GUDHI Editorial Board, 3.6.0 edition, 2022.
[41] Greg Turk and Mare Levoy. Zippered polygon meshes from range images. In Proceedings of the 21st annual conference on Computer graphics and interactive techniques, pages 311-318, 1994.
[42] Yuhei Umeda. Time series classification via topological data analysis. Information and Media Technologies,
2017.
[43] Laurent Younes. Shapes and diffeomorphisms. Springer-Verlags 2010.
13

---

A A more extensive presentation of TDA
The starting point of TDA is to extract quantitative topologiaal information from structured objects-foo 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 pres- ence of underlying topological properties in a multiscale way. Here, estimating topological proper- ties 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 finite simplicial complex² K, a filtration over K is a map t + K; C K that is non-decreasing (for the inclusion). For each o E K, one can record the value t(o) == inf{t] lo E Kr} at which the simplex is inserted in the filtration. From a topological perspective, the insertion of a has exactly one of two effects: either it creates a new topological feature in K2 (e.g., the insertion of an edge destroys an existing feature of lower dimension (e.g., two independent connected components are now connected by the insertion of edge). Relying on matrix reduction algorithm [20, $IV.2), PH identifies, for each topological feature appearing in the filtration, the critical pair of simplices (ob,oa) that created and destroyed this feature, and as byproduct the corresponding birth and ? The collection of intervals (t(os), t(ox)) is a (finite) subset of the open half-plane {(b, d) R²| 6 < d}. called the persistence diagram (PD) of the filtration (Ke).. The distance of such a point (to; (a) to the diagonal (b- = d], namely 2 (((o)) - ((ra)), , 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,).-
create loop, that is, a one-dimensiona. topological feature) it
death times t(r,), aaa).)
Note that if (ab, ..)) is a critical pair for our filtration (Kt)., it holds that l [0%] 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 filtration (Ke)s is empirically" quasi- linear 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 with V n being typically quite large).
= Io The quantity
The Vietoris-Rips Rips filtration. A particular instance of simplicial filtration that will be extensively used in this work the Vietoris-Rip (VR) one. Given X = (Z1..... ...))) E Raxd a point cloud of points in dimension d, one consider the simplicial complex K 22 and then the filtration (Kt). defined by
[ ... .., E K2 S Vj.j' E (l,....p). PP, I| - -. و S L.
(8)
The corresponding persistence diagram will be denoted, for the sake of simplicity, by Dgm(X).
Note that for t < 0. K₁ = 0. when & > dXam(X). K₁ = K, and there is always a point with coordinates (00) in Dgm(X) accounting for the remaining connected component when & + 00. This is the unique point in Dgm(X) for which the second coordinate is +oo (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; are i.i.d. samples from regular measure μr supported on a submanifold M C IRd, 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 (||æ;- (it-2) (2yl)) gi..gn. and thus the complexity of computing Dgm(X) X) depends only linearly in d'l. On the other hand, since Dgm X) scales (at least) linearly
"A simplicial complex is a combinatorial object generalizing graphs and triangulations. Given finite set of vertices V = (FI,.... V,.), a finite simplicial complex is subset 2 (whose elements are called simplices) such that a E K + E K, VT C a (if a simplex is in the complex, its faces must be in it as well). 'It may happen that topological feature appears some time ts and is never destroyed, in which case the death time set to +oc. However, 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 cubic complexity, but the matrix that has to be reduced is typically very sparse, enabling this practical speed up.
Howeeerr 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.
14

---

Secal ext Sad pol del
Sanae and buil prokoees dagras
Eodaate of Co cas beneow
1.00
07
4.29
0.0
Figure 8: Topological optimization an initial point cloud X (in red) by minimizing X H Dggx))))a(X) - +d [-1, This loss favors the apparition of topological features (loops) while the regularization term penalizes points that would go infinity otherwise. Experimmet reproduced following the setting of [4], using code available at https:// http:///gtttbt..coo CAa//UDHII/DDA tutoouttttoiia/blb.. master/Tuto- GUDHI - opt.iizatiiiiin.ooommmii ipynb.
with respect the number of simplices in K, computing the whole VR diagram of point cloud X E RRddd can take up to O(2") operations. Even if one restricts topological features of dimension d' d (e.g. d' = if one only considers loops) as commonly done-thh the complexity is of order O2) which becomes quickly intractable if n is large, even d' oor 2.
B Delayed proofs
Proof of Proposition 3.2. One has:
l(ū(z) - 7(y))1 = (KEE,
K K(y, (y,,x:))(-KK --))((KK-VL(S) (1)):ll I s Eier Ipo (|x (1<x0xa x;ll)
po (ly xill))
since using Gaussian kernels. As ( 2lll) poll =11i| I Co || vll..
2+1 d-1 gd-l (see [2, Theorem 8)), follows that [[((z) - 8(1)112 I 2 π
with Cd.o = 2
Co.o || - vlle |K 11a (L((()))l
Let us upper bound the term
Let us start with the simplification loss, one has at
;b; 24y Tig 2|2 (for some critical points Fig Tia E (t)), one has at
2(b; - d4) and, writing b; = |r,,
Fay Bea | -.g 12 and ab._ || 32g Similarly, one has dd, 22bb d.), and writing 8m.g
ab_ =
21, -ey
dj = | rig 2 provides the corresponding partial derivatives.
Applying (2), this gives:
VV(((())
Σ
Σ
2(b, - di)
T. Fie
Σ
2(b: - di)
Di, 36
1
Nig 112
| x:
xEX(1)
(b,,d))
2 (b.
Σ
2(b; di)
Lig
2(b: di)
Fig
x2|2 11
[x
Fig 112
x (b,,d))
(br,d)
where } (resp. (3) means that x appears as left point (resp. right point) in the computation of the birth filtration value b, of one of the PD points (b,, (b, E (t)(X(t)) associated to the loss, and similarly for death filtration values. A brutal majoration finally gives (V)(((t))l s 2v Ex-(bid,)) 1b, d,l S 8yd Persy (erss(Dgmm (Dgm(((1)). as there are at most four points associated to every (b,, dc) E Dgm (gm(((()) (t)). One can easily see that the same bound applies to the augmentation loss.
Let us finally bound ||K [|,, one has | = <(K)/ S <(K). Indeed, as we are using Gaussian kernels, ||||] = maxigign Σf=1 Pa (1|2 - 25122) 1.
C Complementary experimental results and details
Subsampling and improving over [4]. We reproduce the experiment of [4, $5), see also Figure 8, but starting from initial point cloud Xo of size = 2, 000 instead of n 300. This makes the
15

---

UNilA rpors 750
D Oif ezech 752
Evoution of
Pho
075
333
C22
32
-222
-02-
-222
-03
-1.93
3100
-132
200th
Figure 9: Topological optimization with subsampling. From left right, the initial point cloud Xu. the point cloud after 750 steps of vanilla gradient descent (+subsampling), the point cloud after 750 steps of diffcomor- phic interpolation gradient descent (+subsampling), loss evolution over epochs. Parameters: À = 0.1, o = 0.1.
I I NY
LL
Figure 10: Topologically-optinni LSs and losses for duck, cat, pig, vase and teapot.
computation of Dgm(X+) each gradient step unpractical. Following Section 3.2 rely subsampling with sample size = 100 and apply Algorithm 1. Results are summarized in Figure 9. While relying on vanilla gradients and subsampling barely changes the point cloud even after 750 epochs, the diffeomorphic interpolation gradient with subsampling manages to decrease the loss.
More extensive reports of running time and comments. On the hardware used in our exper- iments (the first two experiments were run on 11th Gen Intel (R) Core (TM) 15-1135G7 2. 40GH2, the last one 2x Xeon SP Gold 5115 2. 2.40GH..), report the approximate fol- lowing running times:
Small point cloud optimization without subsampling (see Figure 4, n 200 points): one gradient descent iteration takes about 1s for the vanilla topological gradient and our diffeo- morphic interpolation. The of oineus integrated in pipeline raises the running time (per iteration) to 10 to 20 seconds. Note that the diffeomorphic interpolation and oineus may converge in less steps than the vanilla topological gradient, preserving competitive advantage. We also believe that oineus has significant room for improvement in terms of running times and may be a promising method in the future be used jointly with approach.
Iterating over the stanford bunny with subsampling (n 35, 947, & = 100) takes about seconds per iteration for the vanilla topological gradient and 20 second for diffcomor-
16

---

Figure 11: Decrease of augmentation (up left) and registration (up right) losses for the scHiC dataset, well as the corresponding initial (orange) and optimized (blue) LSs displayed below.
phic interpolation. The increase in running time with respect to the previous experiment mostly lies instantiating and applying the (d = 3) vector field (requires com- pute p=( - x:)) for each new x (n of them) and sampled o; (77) of them, which is typically very small), hence ~ (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 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 a. We reproduce the same experimental setting as in Fig- ure 3, i.c., sample points uniformly on a circle of radius 1 plus additional noise ~ N 0.05/2). and consider minimizing the total persistence of the point cloud. We take (0. 00, 0.2, 0.3, 0.5, 0.7, 2, 4,5) 5} (with the convention that a = 0 corresponds to the vanilla topological gradient) and learning rate 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 Xo- up max- imum of 200 iterations, stopped earlier if loss of 0 (no topology left, global minimum has been reached) measured.
0.1. We also rely on subsampling with batch size
Figure 12 displays the results of this illustrative experiment. We report the median both running time and number of iterations reach convergence (or reach the 200 iterations limit), along with the 10 and 90 percentiles. The conclusions are:
For 0 (vanilla) and 2 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 bandwidth larger than that hinders convergence.
For E (0.1, 1], the convergence occurs within the same order of magnitude (between 0.49 and 1.74s), the best performance being reached at a 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 a. Empirically, we observe that a good proxy is to take < median( [[=;-1, =y)): Note that even though in theory, - 0 should recover the vanilla topological gradients, one limited by numerical accuracy when evaluating the Gaussian kemel.
17

---

The variation around the median over 50 runs is very small: the randomness of the samples each iteration (hence of the trajectory) barely impacts the decrease of the loss and thus the convergence time.
Impact of and sampling convergence
#lterations to converge stop
200
15.0
12.5
Total running time
150
10.0
7.5
100
5.0
2.5
2 bandwidth
bandwidth
Figure Topological simplification, point cloud diameter with median pairwise distance 10-90 percentiles over 50 runs. (Leftt Time to converge for different values of a E [0, (a 0 corre- sponds to Vanilla), (Right) #iterations to converge (or stop after 200 iterations, indicated by the dashed red line).
Median
We also studied how the bandwidth a influences the correlation scores of Table in Figure 13. We observe oscillations for values that roughly the sides (very small 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.
0.7
apea.
some
Figure 13: Influence of the kemel bandwidth correlation scores for few datasets. The values of are evenly spaced between 0.05 and for the CDIL datascts, and between 0.025 and 0.5 for the scHiC dataset.
18

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
Our abstract is quite extensive, and there is a clear "Contribution and Outline" paragraph at the end of our introduction along with "Related works" and "LLimitatinss" (of previous works) paragraphs.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: This is discussed in the Conclusion section.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Our two theoretical results followed by proof. The first pretty short and included in the main material, the second one is longer and technical and thus deferred the appendix.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
We mentioned important hyper-parameters for our PoC experiments (learning rate X, bandwidth a, number of epochs, stopping criterion if any.). The code is available in the supplementary material as well.
Justification:
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Code provided in the supplementary material with (hopefully) clear Readme.md file which provides some ready-to-use quick start examples, is also available as public Github repository, see https: /tttps:/githhtt conn/tlacombeetop....
6. Experimental Setting/eetalls
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
We hopefully provided the important experimental details to understand our experiments.
Justification:
7. Experiment Statistical Significance
Question:
Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments?
Answer: [Ys]]
19

---

Justification: We did not report error bars in our PoC experiments as they were meant for illustrative purposes. We did report error bars in the form of standard deviations for our main experiments (topological correction of latent spaces of VAEs for images and single cells datasets). These were obtained by generating multiple test sets by adding noise the training set.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory. time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Hardware mentioned in the first paragraph of Section 4, and we reported running times either in the plots in the appendix.
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https://neurip:: https//ccccppcccccccccceeeeeeeee
Answer: [Yes]
Justification: We not use sensitive nor closed datasets or other aspect that may enter in conffict with NeurlPS Code of Ethics.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: We do expect that-in its current state-our article may have any neg- ative societal impact. It is a work that introduced novel method to perform topological optimization. This recent research field does not have direct applications that would have negative societal impact so far. Of course, it falls in the same pitfall as any machine learn- ing related research problem, but we do not see something specific to our work that would be worth discussing.
Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: We used single generative model (VAE) to reproduce images from the COIL dataset and single cells from the Hi-C dataset. We trained it by ourselves.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We used three non-synthetic datasets: the Stanford Bunny, the COIL dataset and the Hi-C dataset, for which we give the proper credits.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documenta- tion provided alongside the assets?
Answer: [NA]
Justification: No new asset.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsoureing experiments and research with human subjects, does the pa- per include the full text instructions given participants and screenshots, applicable, well as details about compensation (if any)?
20

---

Answer: [NA]
Justification: the article does not involve crowdsourcing nor research with human subjects.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: the article does not involve crowdsourcing nor research with human subjects.
21


Reviews:
[{"review_id": "z34CJndwoM", "reviewer": "Reviewer_fiZD", "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.", "questions": "Can this method work for differentiable multiparameter vectorization methods?", "limitations": "Yes, the authors have discussed limitations.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Overall, the paper is well-written and organized. \n\n2. The proposed method converges faster than its competitors, as shown in Figure 4.\n \n3. The idea of diffeomorphically interpolating the gradient on a subsample to the entire sample, to the best of my knowledge, is novel and quite useful in practice. \n\n4. The concept of inducing topological information into the latent space of a pre-trained autoencoder is also new to the best of my knowledge.", "weaknesses": "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. \n\n2. 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. \n\n3. 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? \n\nMinor: \n\nIn Algorithm 1, should it be “Set $X_k \\coloneqq X_{k-1} - \\lambda \\tilde v(X_{k-1})$\"?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> Weaknesses 1. and 2. (Section 3 is notation intensive. // Proof of Proposition 3.1 can be moved to the appendix)\n\nThank you for the suggestion. We will alleviate Section 3, stating Proposition 3.2. more informally and deferring the complete, technical statement to the appendix. As far as the proof of Proposition 3.1 is concerned, as it is very short, we consider keeping it in the main paper as it shows in a simple way how $\\tilde v$ interplays with the loss, which is helpful for developing intuition. \n\n> Weakness 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)\n\nWe agree. One could say that it is precisely the purpose of this experiment: since the vanilla topological gradients are not defined on the whole space (not even the whole input point cloud), and because of the inherent computational limitation of (Vietoris-Rips) persistent homology, it is expected that vanilla gradient descent will perform poorly on this task.\n\nSo yes, to some extent, the comparison is indeed unfair, but this is exactly the point we wanted to emphasize: extending gradients via diffeomorphic interpolations has a crucial impact in terms of computational efficiency, opening the way for new experiments involving topological optimization that are not currently accessible with vanilla gradient descent. \n\n> Minor (Should $X_t$ be $X_k$ ?) \n\nIndeed, that’s a typo. Thank you for catching it!\n\n> Question : Can this method work for differentiable multiparameter vectorization methods?\n\nOur approach may be generalized easily: as long as one is given a (sparse) gradient on a point cloud $X \\in \\mathbb{R}^{n \\times d}$, one can build its diffeomorphic interpolation and use the resulting $\\tilde{v}$ to update $X$ instead. \n\nTherefore, any methods that returns such gradients can benefit from our approach. This includes recent pipelines involving differentiable multi-parameter PH such as D-GRIL [_D-GRIL: End-to-end topological learning with 2-parameter persistence._ Mukherjee et al. arXiv:2406.07100, 2024] and differentiable signed measures [_Differentiability and convergence of filtration learning with multiparameter persistence._ Scoccola et al. **ICML**, 2024]. In both methods, gradients can be obtained by computing free resolutions of multi-parameter persistence modules, and are thus also very sparse (similar to the 1-parameter case studied in our work) as they depend only on the multigraded Betti numbers of the module.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "zhgCpWtCXM", "reviewer": "Reviewer_GJy4", "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.", "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?\n\n(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)\n\n(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?\n\n(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?\n\n(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\n\n(Q6) Do you see application in the real-world of your method, and how big datasets can your approach handle?\n\n(Q7) Can you please elaborate on how exactly your method offers better interpretability for the black-box autoencoder application?\n\n(Q8) How robust is your method to taking different samples and to different sampling strategies?", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "(S1) The paper addresses a relevant problem in the TDA field\n\n(S2) The paper is comprehensive and well-written\n\n(S3) The paper proposes interesting applications of the proposed method", "weaknesses": "(W1) The experimental setup could be strengthen by providing a real-world application", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> (Q1) color swap in Figures 6 and 7\n\nThis was accidental. Thank you for catching this error!\n\n> (Q2) elaborate on (…) black-box AE models, (…) compare to other topological optimisation techniques\n\nSorry: we made a misleading mistake by accidentally naming “Vanilla” the competitor of “Diffeo” (ours) in the Table in page 9. Here, “Vanilla” is supposed to mean _no topological optimization on the latent space_ (like “Vanilla VAE”, not “Vanilla topological gradient descent”). We will change this in the revised version and are deeply sorry for the confusion it yielded. \n\nWe do not compare our method with existing topological optimization methods (Oineus, Vanilla, etc.) in this experiment. That is for a good reason: the experiment we propose can only be made using our method. Indeed, given a pre-trained VAE $D \\circ E$, we design $\\varphi : \\mathbb{R}^{d’} \\to \\mathbb{R}^{d’}$ such that $\\varphi \\circ E$ induces a “better” latent representation (in terms of topological properties) by concatenating our diffeomorphic interpolations $\\tilde{v}_k, k=1,\\dots,T$.\n\nIn contrast, **existing methods do not provide maps defined on the latent space**. Their gradients are only defined on training observations. Given (say) a new latent representation $z$, there is no way to “correct it” using vanilla or oineus gradients.\n\nAside from computational efficiency, this is the biggest edge of our method over pre-existing topological optimization techniques: having access to maps defined on the whole space opens the way for new applications.\n\n> (Q3) ran experiment multiple times?\n\nSee Global Rebuttal and companion pdf, Table 1. As for the Pig result, the large margin is because the initial space (from the VAE) was particularly bad (topologically), see Figure 2(b).\n\n> (Q4) Figure 3 (…) expected end result of topological optimisation? The final shape (…) is very distorted, (…) how is that good?\n\nThe loss that we minimize in this PoC experiment is proportional to the topological information contained in the input point cloud. That is, we want to destroy the topological structure (here, the presence of a loop). A good output (low loss) is a point cloud with no underlying loop. The take-home messages of this experiment are:\n- The vanilla gradient only updates a few points (sparse gradient) of the (small) subsample. This yields a low decrease in the loss: after $t=100$ steps, the point cloud has barely changed. \n- In sharp contrast, our diffeomorphic interpolation updates _the whole point cloud_ at each step. It yields a much faster optimization (showcased also in Figure 4), and the whole point cloud is deformed smoothly from the input one. \n- In terms of raw topological loss, there is no “better” final output: both methods manage to eventually provide a point cloud with trivial topology. Nonetheless, having a smooth, invertible, deformation between the input and final objects can be useful and more interpretable; this is for instance leveraged in the black-box AE experiment.\n\n> (W1, real-world experiment) + (Q5) real-world uses of topological optimisation?\n\nSee Global Rebuttal and the companion pdf (Figure 2(b) and Table 1) for a real-world application on single-cell data. In a nutshell, topological optimization can be useful if you have a topological prior (e.g., \"my model should exhibit a periodic / cyclic structure\"). One can also consider the references given in lines 40--42 in the main paper for other possible applications of topological optimization. \n\n> (Q6) how big datasets can your approach handle?\n\nDiffeomorphic interpolation updates computed from a subsample of the dataset modify the whole input point cloud, a sharp difference with the vanilla gradients that only update a fraction of the original point cloud at each step. Our PoC experiments illustrate this (in particular Figure 5 shows that our approach can seamlessly perform topological optimization on a point cloud with $\\sim 36,000$ points, a scale completely out of reach of previous methods), and we confirmed this on our new real-world application, which involves Vietoris-Rips PH on $1,171$ single cells. While this is already hardly tractable for vanilla gradients, subsampling combined with our diffeomorphic interpolations allows to perform the task fastly and efficiently.\n\nIn terms of computational complexity, the overhead induced by our method compared to vanilla gradients are: \n- inverting the kernel matrix $\\mathbf{K}\\in \\mathbb{R}^{|I| \\times |I|}$ where $I$ is the set of indices for which the vanilla gradient is non-zero. Since the vanilla gradient is usually sparse, $I$ is small and this matrix inversion is basically negligible.\n- Computing $\\tilde{v}$ on the whole input point cloud $X$, which scales _linearly in the number of points_.\n\n> (Q7) elaborate on (…) interpretability for black-box autoencoder?\n\nThe idea showcased in this experiment is the following: each dataset is a collection of images representing a rotating object, hence it is natural to assume that there must be a 2-dimensional underlying circle (topological prior).\nWhen training a VAE on such a dataset, we obtain an encoder $E$ and a decoder $D$ that enables the generation of new data. However, the latent space does not exhibit the expected circle structure (it has no reason to do so).\nWe change the latent space representation using the diffeomorphism $\\varphi$ optimized so that $(\\varphi \\circ E(x_i), i=1,\\dots,n )$ satisfies our topological prior ($(x_i)_{i}$ is our training set). Since $\\varphi$ is invertible, we can generate new observations using $D \\circ \\varphi^{-1}$.\n\nThis increases interpretability: after topological correction, the latent representation is organized on a circle, and angles in the latent space reflect angles between rotated images. You can now sample \"an image rotated by 30°” (wrt some reference image) by simply taking an angle of 30° in the latent space.\n\n> (Q8) How robust (…) to (…) different samples\n\nSee Global Rebuttal and companion pdf, Figure 1.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0HSlYMYteb", "reviewer": "Reviewer_v7wD", "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.", "questions": "1. What is the intuition behind the proof of Theorem 2? (Cf. Weakness 1)\n2. 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?\n3. 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?\n4. What is the impact of different choices of sigma on the results of the experiments? How sensitive are the experiments to varying sigma?\n5. 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? \n6. 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.\n7. 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?\n\nI would like to thank the authors in advance for their answers!\n\nComments:\n* Line 20: \"can be done once and then re-applied to new data in linear time\" -> What does that even mean?\n* ll 33: \"active line OR research\"\n* ll 37: \"yielding to the problem\"\n* 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:\n* ll 100: \"machinary\"\n*  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\n* 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?\n* ll 193: \"i-th line\"\n* ll 214: \"s-points\"\n* ll 476: The grammar of this sentence is confusing.\n* page 13, footnote 5: only non-empty simplices. footnote 8: \"deteriorate\"", "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...", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper presents a novel combination of topological gradients and diffeomorphic interpolation. The authors do a good job of embedding this approach into the existing literature and highlighting its differences and relevancy.\n2. The submission is technically sound and provides theoretical guarantees which are proven. The application to the latent space of variational auto encoders provides a nice and simple application adding to the theoretical contributions.\n3. The proposed method significantly enhances the applicability of optimisation of point clouds w.r.t. to topological loss functions. The idea to solve some of the problems of topological gradients using diffeomorphic interpolation is new and very good.\n(4. The gif is very nice :) )", "weaknesses": "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!\n2. 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.\n3. 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.)\n4. 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.)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. What is the intuition behind the proof of Theorem 2?\n\nOur goal with this result was to quantify the smoothness of our proposed diffeomorphic interpolation (as characterized with its Lipschitz constant) based on the kernel it is computed from. In the case of the Gaussian kernel, we found that the Lipschitz constant depends on the kernel bandwidth $\\sigma$: indeed, the more spread the Gaussian is, the more critical points (identified during the computation of PH) can influence other data points potentially far from them, introducing unwanted distortions and larger Lipschitz constant. Similarly, if the condition number is large, inverting the kernel matrix might introduce instabilities in the $\\alpha_i$ 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 the minimization problem in section 2.2 is, leading to more complex diffeomorphic interpolation solutions with larger Lipschitz constants.\n\n> 2. What do the theoretical guarantees tell us in practice? (…) evolution of loss with discrete step-size? (…) what makes it much harder?\n\nWe agree that the theoretical guarantees remain of limited use in practice: they only tell us that our diffeomorphic interpolations also provide meaningful descent directions for the loss; and thus that for a sufficiently small step-size, the loss must decrease. Empirically, we showcase that overall, the loss decreases substantially faster using our method. Proving this formally in some reasonable settings (continuous-time, etc.) remains quite hard though, for various reasons:\n- the losses considered are typically not convex (with respect to the input point cloud $X$), and thus it is hard to quantify how large the learning rate can be taken and how much the loss decreases after a discrete-time step,\n- Our $\\tilde{v}$ that dictates the dynamic of the flow is not a gradient (at least for Euclidean geometry), so we cannot apply faithfully the gradient flow/descent literature to get theoretical guarantees,\n- The dependence of the Vietoris--Rips filtration on the point cloud $X$ is cumbersome. We can derive theoretical guarantees as long as the ordering of the pairwise distances $\\|x_i - x_j\\|$ remains unchanged (i.e., as long as the points are moved by less than $\\min_{ij,kl} |\\|x_i - x_j\\| - \\|x_k - x_l\\| |$), but this is fairly restrictive. \n\nIt is clear that studying this model is of interest, but we believe that it may be postponed for future work. \n\n> 3. and 4. How did you choose the bandwidth in practice? (…) heuristic ? (…) vary sigma locally ? (…) different choices of sigma ?\n\nWe chose $\\sigma$ with no parameter tuning (see Global Rebuttal and companion pdf, which actually show that the results of our method in Figure 3 can actually be improved if we pick $\\sigma = 0.3$ instead of $0.1$); the heuristic being that $\\sigma$ should be a fraction of the median distance between input points (this heuristic is common in calibration of Gaussian kernel in other statistical scenarios). \n\nIt may be possible to consider a space-dependent $\\sigma$ or other variant: as long as $(x,y) \\mapsto K(x,y)$ remains valued in positive semidefinite matrices, the theory presented in Section 2.2. adapts faithfully (see also line 154 in the main paper). \n\nWe included in Global Rebuttal and companion pdf (Figures 1 and 2(a)) some experiments for varying $\\sigma$, showcasing the robustness of our approach w.r.t. $\\sigma$. \n\n> 5. Oineus loss larger\n\nGood catch, thanks! We forgot a normalization factor in the way the oineus loss is displayed. Fortunately, as our stopping criterion is that a loss of exactly $0$ should be reached (see line 248 in the main paper), this does not change the number of iterations (e.g., 49 iterations for oineus alone) nor the running times, and thus the conclusions of this proof-of-concept experiment remain unchanged. We will correct this in the next version. \n\n> 6. In Figure 3, how did you pick the hyperparameters? (…) It might be the case that different hyperparameter regimes are optimal.\n\nWe actually pick $\\lambda = 0.1$ for every experiments, with no specific tuning. Given that our approach is directly extrapolating the vanilla topological gradient, we believe that picking the same learning rate for both models is the natural choice. \n\n> Figure 9, (a) 3rd from left (…) purple points? (b) 4th (…) loss spikes for the diffeo algorithm be explained?\n\n(a) For 3rd Figure, transparent blue points simply represent the initial point cloud as a reference. Orange points represent the final state of the point cloud. \n\n(b) We report the loss on the subsample batch (because evaluating the loss on the whole point cloud is computationally prohibitive, see line 227). Because of the stochastic nature of our approach, it may happen that this estimated loss is, on a few occasions, quite high because the subsample does not reflect the expected topological structure. \n\nThis does not happen with vanilla gradient descent for a simple reason: the initial point cloud is broadly uniform, so are the subsamples (with high probability) yielding a loss around $-0.03$, and because _vanilla gradient descent fails to perform any significant update of the point cloud_, this situation remains through iterations. \nIn contrast, because our diffeomorphic interpolations do change the point cloud, they reach configurations where some subsamples may yield a high loss, though the loss decreases overall. \n\nIn a nutshell, these spikes are not related to the optimization technique involved (if the vanilla gradient descent was initialized with, say, the epoch 600 reached by our approach based on diffeomorphisms, one would observe similar spikes), but to the “topological variance” in the subsample for a given configuration, the point being that vanilla gradient descent does not reach such configurations.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vDH2eYNmaL", "reviewer": "Reviewer_uH7P", "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.", "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": "See above.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. In practical applications, optimizing topological loss is often slow. The proposed method significantly accelerates this process.\n2. Diffeomorphic interpolation is proven to provide a descent direction for the topological loss, with a provided smoothness bound.\n3. The proposed method is especially beneficial for subsampling techniques, leading to both faster convergence and reduced computation time.", "weaknesses": "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> topological registration losses is not demonstrated\n\nSee Global Rebuttal and companion pdf (Figure 2(b) and Table 1), where we showcased our methodology on a topological registration loss on a real-world dataset of single cells. While we do observe oscillations, we still note a global decrease of the loss. As for the correlation score, the improvement is very marginal, which we believe is due to the fact that correlation score only measures the correlation between latent space angles and ground truth. Thus, while the optimized latent space does exhibit a better delineated circular shape with the registration loss (see Figure 2(b)), the size of that circle is prescribed by the coordinates of the target PD point $[-3.5, 3.5]$ (whereas this size is made as large as possible with the augmentation loss), which in turn does not change the latent space angles very much (as well as the correlation score overall) compared to the augmentation loss. \n\n> 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?\n\nThank you for the suggestion. Proving a convergence theorem is indeed one of the future directions we aim at pursuing. We believe that such convergence results should be accessible provided that existing results in this vein (such as, e.g, [_Optimizing persistent homology based functions_, Carrière et al., **ICML**, 2021]) can be thought of as particular instances of our setting when $\\sigma \\to 0$. Hence, convergence should also occur for $\\sigma$ and learning rates that are sufficiently small, in order to ensure that they do not interfere too much with the different terms involved by the high topological complexity. \n\nGuaranteeing the convergence for fixed bandwidth $\\sigma$ and learning rate $\\lambda$ in the greatest generality seems too demanding: we designed in Global Rebuttal and the companion pdf some cases where convergence seems to fail numerically because $\\sigma$ is much too large with respect to the diameter of the point cloud. As commented in Global Rebuttal, a good heuristic in our experiments is to take $\\sigma$ as a fraction of the median distance in the point cloud.  \nTherefore, an appealing generalization suggested by reviewer v7wD would be to use in practice an adaptive / locally varying $\\sigma$ that would in some sense, adapt to small vs large topological features. Note that one should do it in a way such that $K(x,y)$ remains positive semi-definite for any $x,y$ though, so that the theory of diffeomorphic interpolations presented in Section 2.2 still applies.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
# Global rebuttal

We thank the reviewers for their constructive reviews and overall positive comments on our work. We will take into account all comments on clarity, typo, etc., that have been reported by the reviewers in the revised version of our work. 
In addition to individual responses, we provide complementary experimental results hopefully addressing remaining concerns in our attached companion pdf:

- 1. A proof-of-concept study showcasing the robustness of our method to (i) the choice of $\sigma$ as suggested by Reviewer v7wD and (ii) the randomness induced by the subsampling as suggested by Reviewer GJy4 (Figure 1). The influence of $\sigma$ over the correlation scores is also provided (Figure 2(a)).  
- 2. A real-world experiment (suggested by Reviewers GJy4, v7wD) on single-cell genomics that involves minimizing a _topological registration loss_ as suggested by Reviewer uH7P (Figure 2(b)).
- 3. The results of the black-box VAE experiments averaged over 100 runs as suggested by Reviewer GJy4 (Table 1), in which we compute the mean scores and standard deviations after adding for each run random uniform noise (with amplitude $0.05$ for COIL and $0.01$ for the single-cell dataset).

We describe below the experiments we ran (for 1. and 2.; 3. is exactly the one in the paper but averaged over 100 runs).

**1. Influence of $\sigma$ and subsampling (PoC).** We reproduce the same experimental setting as in Figure 3 of the main paper, i.e. sample points uniformly on a circle of radius $1$ plus additional noise $\sim \mathcal{N}(0,0.05 I_2)$, and consider minimizing the total persistence of the point cloud. 
We take $\sigma \in \{ 0, 0.1, 0.2, 0.3, 0.5, 0.7, 1, 2, 3, 4, 5 \}$ (with the convention that $\sigma = 0$ corresponds to the vanilla topological gradients and learning rate $\lambda = 0.1$ (as in the paper). 
We rely on a subsampling with batch size $s=50$ (as in the paper). 
To quantify the variability of the scheme (suggestion by reviewer GJy4) 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 1 in the companion pdf 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 $\sigma = 0$ (vanilla) and $\sigma \geq 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 $\sigma \in (0.1,1]$, the convergence occurs within the same order of magnitude (between $0.49$ and $1.74$s), the best performance being reached at $\sigma = 0.3$ (recall that we used $\sigma = 0.1$ 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 $\sigma$. Empirically, we observe that a good proxy is to take $\sigma < \mathrm{med}((\|x_i - x_j\|)_{ij})$. Note that even though it theory, $\sigma \to 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.

As for correlation scores in Figure 2(a), we observe oscillations for $\sigma$ 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. 

**2. Real-world application on single-cell data.** We also designed an experiment on single cell HiC data inspired from [_A Gradient Sampling Algorithm for Stratified Maps with Applications to Topological Data Analysis_. Leygonie et al. **Math. Prog.**, 2023]. The dataset is comprised of single cells characterized by chromatin folding, that is, each cell is encoded by the distance matrix of its DNA fragments. The dataset we focus on is taken from [_Cell-cycle dynamics of chromosomal organization at single-cell resolution_. Nagano et al. **Nature**, 2017], in which it was shown that cells are sampled at different phases of the cell cycle. Thus, 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 $500$kb and convolution parameter $h=1$ on chromosome $10$. Then, we run kernel PCA on the SCC matrix to obtain a preprocessed dataset of dimension $100$, on which we applied the same VAE architecture than the one used for COIL data in the main paper. Finally, we optimized the same augmentation loss than for the COIL data, as well as the following registration loss:

$$L :  X\mapsto W_{2}^2({\rm Dgm}^1(X), {\rm Dgm}^{1, t}),$$

where ${\rm Dgm}^1(X)$ contains the PD points of $X$ with distance-to-diagonal at least $\tau=1$, ${\rm Dgm}^{1,t}$ is a target PD with only one point $p^*=[-3.5,3.5]$, and $W_{2}$ is the 2-Wasserstein distance. We used $\sigma=0.25$ (aug.) and $\sigma=0.025$ (reg., sigma is smaller to mitigate the effects of matching instability), $\lambda=0.1$ on $500$ epochs, with subsampling of size $s=300$ for computational efficiency (as computing Vietoris-Rips PH 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 we measured correlation between latent space angles and repli scores in Table 1.

Meta-review: 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.
Decision: Accept (poster)