# Topological Autoencoders++: Fast and Accurate Cycle-Aware Dimensionality Reduction

Mattéo Clémot, Julie Digne, Julien Tierny

► To cite this version:

Mattéo Clémot, Julie Digne, Julien Tierny. Topological Autoencoders++: Fast and Accurate Cycle-Aware Dimensionality Reduction. 2025. hal-05004858v1

HAL Id: hal-05004858

https://hal.science/hal-05004858v1

Preprint submitted on 25 Mar 2025 (v1), last revised 16 Dec 2025 (v2)

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Topological Autoencoders++: Fast and Accurate Cycle-Aware Dimensionality Reduction

Mattéo Clémot, Julie Digne, Julien Tierny

**Abstract**—This paper presents a novel topology-aware dimensionality reduction approach aiming at accurately visualizing the cyclic patterns present in high dimensional data. To that end, we build on the *Topological Autoencoders (TopoAE)* [1] formulation. First, we provide a novel theoretical analysis of its associated loss and show that a zero loss indeed induces identical persistence pairs (in high and low dimensions) for the 0-dimensional persistent homology (PH<sup>0</sup>) of the Rips filtration. We also provide a counter example showing that this property no longer holds for a naive extension of TopoAE to PH<sup>d</sup> for d ≥ 1. Based on this observation, we introduce a novel generalization of TopoAE to 1-dimensional persistent homology (PH<sup>1</sup>), called TopoAE++, for the accurate generation of cycle-aware planar embeddings, addressing the above failure case. This generalization is based on the notion of *cascade distortion*, a new penalty term favoring an isometric embedding of the 2-chains filling persistent 1-cycles, hence resulting in more faithful geometrical reconstructions of the 1-cycles in the plane. We further introduce a novel, fast algorithm for the exact computation of PH for Rips filtrations in the plane, yielding improved runtimes over previously documented topology-aware methods. Our method also achieves a better balance between the topological accuracy, as measured by the Wasserstein distance, and the visual preservation of the cycles in low dimensions. Our C++ implementation is available at https://github.com/MClemot/TopologicalAutoencodersPlusPlus.

**Index Terms**—Topological data analysis, persistent homology, dimensionality reduction.

## I. INTRODUCTION

DIMENSIONALITY reduction (DR) [9], [10] is a fundamental tool in data science for apprehending high-dimensional data. Such datasets are typically represented as a set of $n$ points $X = \{x_0, x_1, \dots, x_{n-1}\}$ living in a high-dimensional Euclidean space $\mathbb{R}^h$. Then, the goal of DR is to compute an embedding $Z = \phi(X)$ of this input point cloud in a lower-dimensional Euclidean space $\mathbb{R}^l$, with $l \ll h$. This process can be motivated in practice by the need to either (i) simplify the input data (to facilitate further processing and analysis) or to (ii) generate a visualization of the data that can be interpreted by a human observer (e.g., $l = 2$ yields planar views of the input point cloud).

Unless $X$ is already *flat* (e.g., living in an $l$-dimensional hyperplane), its projection to $\mathbb{R}^l$ will necessarily involve a form of *distortion*. Depending on the application needs, several distortion criteria have been introduced in the literature. Then, a vast corpus of DR algorithms [11]–[13] have been proposed,

either (i) to minimize a specific distortion criterion, or (ii) to preserve certain geometrical features (thereby indirectly minimizing the associated distortion criterion). For instance, the seminal multidimensional scaling algorithm [10] aims at minimizing metric distortion, which is achieved in practice by trying to preserve pairwise distances.

In many use-cases, it may be desirable to preserve topological characteristics when projecting down to a visual space. For instance, the preservation of the *clusters* present in high-dimensions enables a quick visual interpretation of the main trends in the dataset. *TopoMap* [7], for example, guarantees that, after its projection, the hierarchical clustering of the dataset is identical in high and low dimensions. This provides users with strong confidence regarding the structural interpretation of the projected data. Technically, the preservation of hierarchical clusters relates to the preservation of the 0-dimensional persistent homology (PH<sup>0</sup>) of the Rips filtration of $X$ (see Sec. III-B for a technical description). In this context, several topology-aware DR methods [1], [7], [14] have been proposed, but these mostly focused in practice on preserving PH<sup>0</sup>, i.e., the data point clusters (see Sec. II-B).

In this work, we investigate the preservation of more advanced topological features than clusters (captured by PH<sup>0</sup>). In particular, we focus on the preservation of the cyclic structures, the *topological handles* captured by PH<sup>1</sup>, when projecting to a planar visual space (i.e., $l = 2$). The presence of such cycles in the input high-dimensional point cloud can imply the existence of periodic phenomena in the data, or the possibility of multiple paths between clusters (e.g., multiple sequences of transition states). As showcased later in Sec. VII-D, the preservation of such cyclic structures can be of paramount importance for the visual interpretation of real-life high-dimensional datasets.

However, preserving PH<sup>1</sup> while projecting to the plane is a significantly more challenging problem than preserving PH<sup>0</sup>. First, from a computational point of view, the 0-dimensional case is a very specific instance of PH computation, which can be efficiently computed with a union-find data structure [15] (with a runtime complexity that is virtually linear with the number of considered edges in the Rips complex). As of dimension 1, generic PH computation algorithms need to be considered, with a worst case runtime complexity that is cubic with the number of considered simplices. Second, as discussed by several authors [1], [7], PH<sup>0</sup> preservation can be obtained by finding a planar layout of $X$ which preserves its minimum spanning tree (MST). Since trees can always be embedded in the plane, the existence of (at least) one solution (i.e., a projection preserving PH<sup>0</sup>) is guaranteed.

Mattéo Clémot is with University Claude Bernard Lyon 1. E-mail: matteo.clemot@univ-lyon1.fr

Julie Digne is with the CNRS and University Claude Bernard Lyon 1. E-mail: julie.digne@cnrs.fr

Julien Tierny is with the CNRS and Sorbonne University. E-mail: julien.tierny@sorbonne-universite.fr
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison of DR methods on three synthetic 3D point clouds (see Sec. VII-A for a description), along with the metric distortion $\mathcal{D}(X, Z)$ between the input $X$ and its planar embedding $Z$, the Wasserstein distance $PDW^1(X, Z) = W_2(D_{\text{Rips}}^1(X), D_{\text{Rips}}^1(Z))$ between their respective 1-dimensional persistence diagrams, and the running time $t$ in seconds. The best value for an indicator is written in bold, the second best is underlined. In the second line, a generator of the most persistent $\text{PH}^1$ pair in the high-dimensional input is projected in the 2D embedding in transparency and slightly smoothed for visualization purposes (to distinguish it from the point cloud). The generator color map depicts its arc-length parameterization. In the third line, the input – that is sampled around the edges of a tetrahedron in 3D – features 3 significantly persistent $\text{PH}^1$ pairs, and a generator for each is represented both in the input and in the 2D embeddings with a specific color. Quantitatively, our approach (TopoAE++) generates planar embeddings with competitive Wasserstein distances. Qualitatively, it produces planar embeddings yielding less crossings of the projected high-dimensional persistent generators (colored curves), thereby producing visualizations that depict more faithfully the topological handles present in high dimensions.</figcaption></figure>

Unfortunately, this observation no longer holds for higher dimensional PH. As discussed later in Sec. VII-D, it is easy to design examples of high-dimensional point clouds whose $\text{PH}^1$ cannot be preserved through planar projections (e.g., a dense sampling of a non-planar graph, see Fig. 13). The existence of such cases, for which an exact preservation of $\text{PH}^1$ cannot be obtained, illustrates the difficulty of the problem, and the need for efficient optimization methods for preserving $\text{PH}^1$ as accurately as possible, which is the topic of this paper.

Specifically, this paper revisits the *Topological Autoencoders* (TopoAE), a prominent technique in topology-aware dimensionality reduction, which we extend into TopoAE++, for the accurate visualization of the cyclic patterns present in the data. We make the following new contributions:

1) A novel theoretical analysis of TopoAE's loss ($\mathcal{L}_{\text{TAE}}^0$):

* We show that $\mathcal{L}_{\text{TAE}}^0$ is an upper bound of the Wasserstein distance between the 0-dimensional persistence diagrams in high and low dimensions, with regard to the Rips filtration. In particular, when $\mathcal{L}_{\text{TAE}}^0 = 0$, the persistence pairs are identical in high and low dimensions for the 0-dimensional persistent homology ($\text{PH}^0$) of the Rips filtration.

* We provide a counter example showing that the above property does not hold for a naive extension of TopoAE to 1-dimensional persistent homology ($\text{PH}^1$), with a loss noted $\mathcal{L}_{\text{TAE}}^1$.

2) A generalization of TopoAE's loss for 1-dimensional persistent homology ($\text{PH}^1$): We introduce a new loss called *cascade distortion* (CD), noted $\mathcal{L}_{\text{CD}}^1$, which addresses the above counter example. This term favors an isometric embedding of the 2-chains filling persistent 1-cycles.

We show that for $\text{PH}^0$ this loss generalizes TopoAE's loss (i.e., $\mathcal{L}_{\text{CD}}^0 = \mathcal{L}_{\text{TAE}}^0$). Extensive experiments show the practical performance of this loss for the accurate preservation of generators through the projection.

3) An efficient computation algorithm for TopoAE++: To accelerate runtime, we provide a new, fast algorithm for the exact computation of PH for Rips filtrations in the plane. This geometric algorithm relies on the fast identification of cycle-killing triangles via local minmax triangulations of the relative neighborhood graph, leveraging duality for the efficient computation of $\text{PH}^1$. We believe this contribution to be of independent interest.

4) Implementation: We provide a C++ implementation of our algorithms that can be used for reproducibility.

## II. RELATED WORK

The literature related to our work can be classified into two categories: general purpose DR techniques (Sec. II-A) and topology-aware techniques (Sec. II-B).

### A. General purpose dimensionality reduction

Numerous DR techniques have been proposed and the related literature has been summarized in several books [9], [10] and surveys [11]–[13]. Principal Component Analysis (PCA) [2] is by far the most popular linear DR technique. Although it is an indispensable tool for data analysis, its linear nature does not always allow it to apprehend complex non-linear phenomena. One of the first non linear DR methods is the multidimensional scaling (MDS) [3]. It aims at preserving as far as possible the pairwise distances in the high- and low-dimensional point clouds. Another approach, particularly related to our work, consists in optimizing an autoencoder neural
network [16]. The *encoder* is used to represent the explicit
projection map from the high-dimensional input space to the
low-dimensional representation space, while the *decoder* tries
to reconstruct the input data from its encoded representation.
We will refer to these methods as *global* methods.

Global methods have been used successfully in many ap-
plications, but they do not take into account the possible
distribution of the input points over some implicit, unknown
manifold. This may lead to the unwanted preservation of
distances between points that are close in the ambient space
but far apart on this manifold. *Locally topology-aware* methods
have therefore been introduced to address this issue. For in-
stance, Isomap [4] preserves geodesic distances on a captured
manifold structure of the input data. Other methods also fea-
ture neighborhood preservation objectives. For example, Local
Linear Embedding (LLE) [17] relies on linear reconstructions
of local neighborhoods. Other methods leverage additional
landmarks [18] or user-provided control points [19].

All these methods try to preserve local Euclidean distances when projecting to a lower dimension. However, this can sometimes lack relevance in the applications, especially with high-dimensional datasets for which the distribution of pairwise Euclidean distances tend to be uniform. For such cases, local distance preservation fails at characterizing well relevant local relations. To alleviate this issue, SNE [20] and later t-SNE [5] use a conditional probability formulation to represent similarities between points and try to have similar distributions both in high- and low-dimension thanks to a Kullback-Leibler divergence minimization. More recently UMAP has been introduced [6] along a theoretical foundation on category theory. It provides results that are similar visually to t-SNE, but in a more scalable way. Variants were later introduced to better preserve the global structure in the embedding, such as TriMAP [21] that constrains the proximity order within triplets of points, or PaCMAP [22] that adds constraints on more distant point pairs. Although these methods succeed in preserving the local topology, they are not explicitly aware of the global structure of the input, which can lead to the loss of noteworthy global, topological features.

B. *Globally topology-aware dimensionality reduction*

Topology-based methods have become popular over the last
two decades in data analysis and visualization [23] and have
been applied to various areas: astrophysics [24], [25], bio-
logical imaging [26]–[28], quantum chemistry [29]–[31], fluid
dynamics [32], [33], material sciences [34]–[36], turbulent
combustion [37], [38]. They leverage tools that define concise
signatures of the data based on its topological properties
and that have been summarized in topological data analysis
reference books [39], [40] and surveys [41].

Several DR methods have been proposed by the visual-
ization community to preserve specific topological signatures
of the input data. For instance, terrain metaphors have been
investigated for the visualization of an input high-dimensional
scalar field, in the form of a three-dimensional terrain, whose
elevation yields an identical contour tree [42] or an identi-
cal set of separatrices [43], [44]. This framework has been

extended to density estimators [45]–[48] for the support of
high-dimensional point clouds. However, such metaphors com-
pletely discard the metric information of the input space [45],
possibly placing next to each other points which are arbitrarily
far in the input space (and reciprocally). Yan et al. [49]
introduced a DR approach driven by the Mapper structure [50],
an approximation of the Reeb graph [51], which can capture in
practice large handles in the data, however without guarantees,
since the number of handles in the considered manifold is only
an upper bound on the number of loops in the Reeb graph [39].

To incorporate the metric information from the input data
while preserving at the same time some of their topological
characteristics, several approaches have focused on driving the
projection by the *persistence diagram* of the Rips filtration of
the point cloud (see Sec. III-B for a technical description).
Carriere et al. [8] presented a generic persistence optimization
framework with an application to dimensionality reduction.
Their approach explicitly minimizes the Wasserstein distance
(Sec. III-B) between the 1-dimensional persistence diagrams
in high and low dimensions. However, this approach solely
focuses on this penalization term. As a result, although the
number and persistence of cycles may be well-preserved,
the solver tends to produce cycles in low dimensions which
involve arbitrary points (e.g., which were not necessarily
located along the cycles in high dimensions), which challenges
visual interpretation, as later detailed in Sec. VII-D.

To enforce a correspondence between the topological fea-
tures at the data point level, additional structures need to be
preserved. For the specific case of 0-dimensional persistent ho-
mology (PH<sup>0</sup>), Doraiswamy et al. introduced *TopoMap* [7], an
algorithm which constructively preserves the *persistence pairs*
(Sec. III-B) through the preservation of the minimum spanning
tree of the data. An accelerated version, with improved layouts,
has recently been proposed [52]. Alternative approaches have
considered the usage of an optimization framework (typically
based on an autoencoder neural network [16]), with the
integration of specific topology-aware losses [1], [14], [53]–
[55]. Among them, a prominent approach is the *Topological
Autoencoders* (TopoAE) [1]. Its loss aims at preserving the
diameter of the simplices involved in persistence pairs, when
going from high to low dimensions and reciprocally. However,
the above techniques focused in practice on the preservation
of PH<sup>0</sup> and did not, to our knowledge, report experiments
regarding the preservation of higher dimensional PH. Specif-
ically, we show in Sec. IV that, while a zero TopoAE loss
indeed implies a preservation of the persistence pairs for PH<sup>0</sup>,
it is not the case for higher dimensional PH. We provide a
counter example for PH<sup>1</sup>, which is addressed by our novel,
generalized loss.

III. BACKGROUND

This section presents the technical background to our work.
We refer the reader to textbooks [39], [40] for comprehensive
introductions to computational topology.

A. *Geometric structures*

We introduce below geometric concepts related to persistent
homology (Sec. III-D) that are relevant to our algorithm.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Left: geometric graphs for a point cloud X in the plane. Its MST is depicted in red. The edges from the RNG which are not in the MST are shown in continuous yellow, while those which are in the UG but not in the RNG are shown in dashed yellow. The remaining Delaunay edges are shown in gray. Three lenses are shown: the purple one is devoid of points of X, therefore the associated yellow edge belongs to the RNG; the blue one – associated with the dashed yellow edge – is devoid of points of the link of that edge but contains another point of X (the blue one), therefore this edge is in the UG but not in the RNG; the green one contains the green point and the associated gray edge is not in RNG nor in UG. Right: within the highlighted bottom RNG-polygon, a MML triangulation has replaced the Delaunay triangulation, with its longest edge highlighted in black. Note that for the two other polygons, the Delaunay triangulation is already a MML triangulation within these polygons.</figcaption></figure>

We consider in this section a point cloud $X$ containing $n$ points, in a $d$-dimensional space $\mathbb{R}^d$ endowed with the ambient Euclidean metric. The lens of an edge $(a, b) \in X^2$ is defined as the points of $X$ whose distance to either $a$ or $b$ is smaller than $\|a-b\|_2$. The relative neighborhood graph RNG($X$) [56] is the graph on $X$ that features every edge with an empty lens, i.e., $\text{RNG}(X) = \{e \in X^2 | \text{Lens}(e) = \emptyset\}$. The RNG is known to be a subset of the Delaunay triangulation [56]. In the plane, it can be computed in time $O(n \log n)$ [57]–[59]. The Urquhart graph UG($X$) [60] is another subgraph of the Delaunay triangulation that is obtained by removing the longest edge of each triangle. It can also be thought as the subgraph of $\text{DEL}(X)$ that only keeps edges whose lens does not contain any points of its link<fnref n="1" />. Finally, we write $\text{MST}(X)$ the Euclidean minimum spanning tree of $X$, i.e., the tree spanning $X$ that minimizes the sum of the Euclidean distances of its edges. Under general position hypothesis<fnref n="2" />, we have the following inclusions [56] (Fig. 2, left, shows a case where the three inclusions are strict):

$$ \text{MST}(X) \subseteq \text{RNG}(X) \subseteq \text{UG}(X) \subseteq \text{DEL}(X). \quad (1) $$

Our algorithm also relies on the minmax length (MML) triangulations of a planar point cloud $X \subset \mathbb{R}^2$. It is defined as a triangulation that minimizes the length of its longest edge (Fig. 2, right). Such a triangulation can be computed in quadratic time $O(n^2)$ with an algorithm [61] that breaks down

the construction of such a triangulation over each polygon formed by the RNG and the convex hull of $X$.

We conclude this section with a geometric construction that is often used in the topological data analysis of point clouds: the Rips complex of $X$ with diameter threshold $r \ge 0$ is the simplicial complex featuring all the simplices $\sigma \subset X$ whose diameter $\delta(\sigma) = \max_{x,y \in \sigma} \|x - y\|_2$ is smaller than $r$ (i.e., every edge included in $\sigma$ must be shorter than $r$) (see Fig. 3):

$$ \text{Rips}_r(X) = \{\sigma \subset X \mid \delta(\sigma) \le r\}. \quad (2) $$

We write $\mathcal{K}^k = \{\sigma \in \mathcal{K} | \dim \sigma = k\}$ the set of $k$-simplices of a simplicial complex $\mathcal{K}$, and $\mathcal{K}^{(k)} = \{\sigma \in \mathcal{K} | \dim \sigma \le k\}$ its $k$-skeleton. Then, $\text{Rips}_r^{(k)}(X)$ is the Rips complex of $X$ with diameter $r \ge 0$ and maximum simplex dimension $k$.

## B. Persistent homology

Persistent homology, introduced independently by several authors [62]–[65], focuses on the evolution of the topological properties of an increasing sequence of $m$ nested simplicial complexes $\mathcal{K} = (\mathcal{K}_0, \mathcal{K}_1, \dots, \mathcal{K}_{m-1})$, which is called a filtration. Persistent homology can be applied on various filtrations, like the sublevel sets of a scalar field on a mesh. However a filtration of particular interest for point clouds is the Rips filtration (Fig. 3), which consists in the sequence of Rips complexes with an increasing threshold $r$. More precisely, the Rips filtration of $X$ can be written as $(\text{Rips}_0(X), \text{Rips}_{r_1}(X), \dots, \text{Rips}_{r_{m-1}}(X), \text{Rips}_{\infty}(X))$ with increasing $r_i$'s such that the Rips complex changes between $r_i$ and $r_{i+1}$ (i.e. $\text{Rips}_{r_i}(X) \subsetneq \text{Rips}_{r_{i+1}}(X)$).

Persistent homology groups for a filtration $\mathcal{K}$ are algebraic constructions that contain persistent homology classes and that are defined from homology groups (see Sec. A-A for a brief summary or [39] for a more complete introduction). Informally, these classes represent topological features (that can be geometrically interpreted as clusters for $\text{PH}^0$, cycles for $\text{PH}^1$, voids for the $\text{PH}^2\dots$) which appear at some point in the filtration and may disappear later. Formally, for $0 \le i \le j \le m$, the $k$-th persistent homology group $\mathcal{H}_k^{i,j}(\mathcal{K})$ is the image of the $k$-th homology group $\mathcal{H}_k(\mathcal{K}_i)$ by the inclusion morphism $f_k^{i,j}: \mathcal{H}_k(\mathcal{K}_i) \to \mathcal{H}_k(\mathcal{K}_j)$. Its classes can be understood as those which already exist in $\mathcal{H}_k(\mathcal{K}_i)$ and that still exist in $\mathcal{H}_k(\mathcal{K}_j)$. Therefore, we say that a class $\gamma \in \mathcal{H}_k(\mathcal{K}_i)$ is born entering $\mathcal{K}_i$ when $\gamma \notin \mathcal{H}_k^{i-1,i}(\mathcal{K})$, i.e., when it is not the image of a pre-existing class by the inclusion morphism. We also say that such a class born at $\mathcal{K}_i$ dies entering $\mathcal{K}_j$ when $f_k^{i,j-1}(\gamma) \notin \mathcal{H}_k^{i-1,j-1}(\mathcal{K})$ but $f_k^{i,j}(\gamma) \in \mathcal{H}_k^{i-1,j}(\mathcal{K})$, i.e., when $\gamma$ is merged into a pre-existing class (possibly the $0$ class) when going from $\mathcal{K}_{j-1}$ to $\mathcal{K}_j$. In this case, the birth of $\gamma$ is associated with a $k$-simplex $\sigma_i \in \mathcal{K}_i$ and its death with $(k+1)$-simplex $\sigma_j \in \mathcal{K}_j$. Such a pair $(\sigma_i, \sigma_j)$ is called a persistence pair. If the filtration is defined by scalar values $f$ on the simplices (e.g., the diameter $\delta$ for the Rips filtration), the persistence of such a class is the positive value $f(\sigma_j) - f(\sigma_i)$.

In the specific case of a Rips filtration, we can say more simply that for any $k$, a $\text{PH}^k$ pair of positive persistence $(\sigma_i, \sigma_j)$ is created (resp. killed) by the longest edge of $\sigma_i$ (resp. $\sigma_j$). Indeed, a Rips complex is entirely determined by

<footnote n="1">¹In an abstract simplicial complex $\mathcal{K}$ (i.e., a family of sets that is closed under taking subsets, where sets represent simplices and their subsets represent their faces), the link of a simplex $\sigma \in \mathcal{K}$ is the set $\{\tau \in \mathcal{K} | \sigma \cap \tau = \emptyset \text{ and } \sigma \cup \tau \in \mathcal{K}\}$. It provides a notion of neighborhood boundary for $\sigma$.</footnote>

<footnote n="2">²i.e., we suppose that no $d+2$ points lie on a common $(d-1)$-sphere (for the Delaunay triangulation [39]), and that pairwise distances are unique for the uniqueness of the MST, which can be enforced via small perturbations.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Rips complexes of the same planar point cloud X, for increasing diameter thresholds r (only the 2-skeletons, i.e., the vertices, edges and triangles, of the simplicial complexes are shown). This threshold increase induces a sequence of nested simplicial complexes, whose topology varies along the process. From left to right, the number of connected components is successively n, then 2, then 1, while the number of handles is successively 0, 2, 1, 2, 1 and 0.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Two point clouds (red and blue) in the plane represented with their RNG, along with their respective 1-dimensional persistence diagrams (right). The edges in RNG \ MST are highlighted in bold and yellow. The number of non-diagonal points in each diagram is exactly the number of RNG-polygons in the associated point cloud (3 for the red one, 2 for the blue one), and also the number of RNG \ MST edges. The optimal assignment inducing the Wasserstein distance $W_2$ between them is shown in black.</figcaption></figure>

its 1-skeleton since the diameter of any simplex is the length of the longest edge it features (it is a flag complex). In the following we refer to these edges as k-critical edges.

Persistence diagrams (Fig. 4) are concise topological representations that summarize the persistent homology groups. More precisely, the k-th persistence diagram is a 2D multiset where each k-dimensional persistence pair $(\sigma_i, \sigma_j)$ is represented by a point $(f(\sigma_i), f(\sigma_j)) \in \mathbb{R}^2$ above the diagonal. We write $\mathcal{D}^k(f)$ for the k-th persistence diagram of the filtration given by a scalar function $f: \mathcal{K} \to \mathbb{R}$, or directly $\mathcal{D}_{\text{Rips}}^k(X)$ for the k-th persistence diagram of the Rips filtration of $X$. Many points clouds have the same persistence diagrams: in particular, $\mathcal{D}_{\text{Rips}}^k(X)$ is invariant to point permutations in $X$.

Persistence diagrams (of equal homology dimension) are usually compared with metrics born of optimal transport that come with stability properties [66], [67]. Let $\mathcal{D}, \mathcal{D}' \subset \mathbb{R}^2$ be two persistence diagrams to be compared. We define the augmented diagrams as $\bar{\mathcal{D}} = \mathcal{D} \cup \Delta\mathcal{D}'$ and $\bar{\mathcal{D}}' = \mathcal{D}' \cup \Delta\mathcal{D}$ where $\Delta$ is the map that projects to the diagonal, i.e., $\Delta: (b, d) \mapsto (\frac{b+d}{2}, \frac{b+d}{2})$. This ensures that $|\bar{\mathcal{D}}| = |\bar{\mathcal{D}}'|$. The p-cost $c_p(x, y)$ between two points $x \in \bar{\mathcal{D}}$ and $y \in \bar{\mathcal{D}}'$ is defined as 0 if both $x$ and $y$ are on the diagonal and as $\mathbb{R}^2$'s p-norm $\|x-y\|_p$ otherwise. Then, the $L_q$-Wasserstein distance $W_q^p$ is:

$$W_q^p(\mathcal{D}, \mathcal{D}') = \min_{\psi \in \Psi} \left( \sum_{x \in \bar{\mathcal{D}}} c_p(x, \psi(x))^q \right)^{1/q}, \quad (3)$$

where $\Psi = \text{Bij}(\bar{\mathcal{D}}, \bar{\mathcal{D}}')$ (Fig. 4). Because for any fixed $q$, all $W_q^p$ are bi-Lipschitz equivalent [67], we focus on the distances

$W_p = W_p^p$ with $p < \infty$. In practice, the Wasserstein distances $W_p$ can be approximated thanks to the auction algorithm, an asymmetric and generic approach introduced by Bertsekas [68] that can be specialized to persistence diagrams with dedicated data structures to leverage their geometric structure [69].

### C. Persistence homology computation

Since the introduction of computational PH [65], numerous algorithms and implementations have been proposed, based either on the boundary matrix reduction approach or on its geometric interpretation known as the PairCells algorithm [40], which we detail in the following as our novel loss (introduced in Sec. V) is defined relatively to the cascade, an object computed by this algorithm. Indeed, PairCells computes not only persistent pairs but also representatives of each PH class for any simplicial complex $\mathcal{K}$ with its simplices endowed with distinct scalar values $f(\sigma)$ that are increasing by inclusion (which induces the filtration of sublevel sets). For each k-simplex $\sigma \in \mathcal{K}$, it stores a k-chain cascade($\sigma$) that is initially equal to $\sigma$, and whose $(k-1)$-boundary $\partial(\text{cascade}(\sigma))$ is, after the execution of the algorithm, a representative of the PH class killed by $\sigma$. In the following we will call this boundary a persistent generator of the PH class [70], [71]. In addition, these cascades will enable us in Sec. V to find 2-chains that fill the 1-cycles.

PairCells (Alg. 1) takes one k-simplex $\sigma \in \mathcal{K}$ at a time, by increasing value $f(\sigma)$, and applies the EliminateBoundaries procedure (Alg. 2) on it. This procedure runs a homologous propagation [71] from $\sigma$: the $(k-1)$-boundary $\partial(\text{cascade}(\sigma))$ is expanded step by step, each time selecting its youngest $(k-1)$-simplex $\tau$ (i.e., the one with the highest value $f(\tau)$) and adding modulo 2 to cascade($\sigma$) the cascade associated with the partner of $\tau$, if it exists (see Fig. 5). The procedure stops either when an unpaired simplex $\tau$ is selected, in which case $\sigma$ and $\tau$ are paired together since $\sigma$ kills the $(k-1)$-dimensional homology class introduced by $\tau$; or when the boundary $\partial(\text{cascade}(\sigma))$ becomes empty, in which case we know that $\sigma$ creates a k-dimensional homology class. This algorithm can be understood as a geometric interpretation of the persistence algorithms based on boundary matrix reduction.

In practice, if we are not interested in the cascades themselves, it can be faster to only store the boundary of each simplex's cascade $\partial(\text{cascade}(\sigma))$ (instead of the whole cascade
<figure><img src="image_6.png" /><figcaption>Fig. 5. Some steps of two executions of the *EliminateBoundaries* Alg. 2 procedure run on the 2-simplices $\sigma_0$ (top line) and $\sigma_1$ (bottom line), corresponding to the creation of two PH¹ pairs in the previous planar point cloud (with its smallest cycle removed for clarity). Unpaired edges – which are in RNG($X$) \ MST($X$) – are depicted in yellow. At each step, *cascade*[$\sigma$] is shown in red, and its boundary $\partial$(*cascade*[$\sigma$]) is highlighted in bold, with its longest edge, i.e., $\tau$ = Youngest($\partial$(*cascade*[$\sigma$]))), dotted. When it exists, *cascade*[*partner*[$\tau$]] is shown in blue. A persistent pair ($\tau$, $\sigma$) is created when *partner*[$\tau$] = $\emptyset$, i.e., when $\tau$ coincides with a yellow unpaired edge (rightmost images). In both executions, Rips<sub>$\delta(\sigma)$</sub><sup>1</sup> is shown in light gray.</figcaption></figure>

itself) and to perform a modulo 2 addition operation on those boundaries when expanding a cascade (line 6 of Alg. 2) [71]. Another acceleration is the preliminary pairing of the *apparent pairs*, pairs of simplices with the same filtration value [71], [72] (some of them appear in Fig. 5 as these dotted red edges whose partner's cascade is a single triangle).

<table><caption>**Algorithm 1: PairCells**</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>for</td><td>σ ∈ 𝒦 by increasing f(σ) do</td></tr><tr><td>2</td><td></td><td>partner[σ] ← ∅</td></tr><tr><td>3</td><td></td><td>cascade[σ] ← σ</td></tr><tr><td>4</td><td></td><td>EliminateBoundaries(σ)</td></tr><tr><td>5</td><td></td><td>if ∂(cascade[σ]) ≠ ∅ then</td></tr><tr><td>6</td><td></td><td>τ ← Youngest(∂(cascade[σ]))</td></tr><tr><td>7</td><td></td><td>partner[σ] ← τ</td></tr><tr><td>8</td><td></td><td>partner[τ] ← σ</td></tr></tbody></table>

<table><caption>**Algorithm 2: EliminateBoundaries**</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>while</td><td>∂(cascade[σ]) ≠ ∅ do</td></tr><tr><td>2</td><td></td><td>τ ← Youngest(∂(cascade[σ]))</td></tr><tr><td>3</td><td></td><td>if partner[τ] = ∅ then</td></tr><tr><td>4</td><td></td><td>return</td></tr><tr><td>5</td><td></td><td>else</td></tr><tr><td>6</td><td></td><td>cascade[σ] ←<br/>cascade[σ] + cascade[partner[τ]]</td></tr></tbody></table>

## D. Geometric interpretation of persistent homology

Recent boundary matrix reduction-based open-source implementations include *Gudhi* [73], *PHAT* [74] or *Ripser* [72]. In particular, *Ripser* is specialized in the PH computation of Rips filtrations (given by point clouds or distance matrices) yielding much faster and more memory-efficient execution

than previous implementations. However, it is sometimes possible to leverage the Euclidean embedding of point clouds for which we want to compute the PH of the Rips filtration to speed up the computation. In this section, we highlight some known results that link PH⁰ and PH¹ of the Rips filtration to some geometric structures defined in Sec. III-A.

The first result is that the set of edges that destroy a PH⁰ class of the Rips filtration (i.e., that reduce the number of connected components by 1) is exactly the set of edges of the Euclidean minimum spanning tree (see [7] for example for a proof). It allows to compute the PH⁰ of the Rips filtration by constructing the minimum spanning tree of the input point cloud, e.g., with Kruskal's algorithm based on the union-find data structure [15], in time complexity $O(n^2\alpha(n^2))$ where $\alpha$ is the inverse of the Ackermann function. This approach is largely used in practice, e.g. in *Ripser* [72].

On the contrary, the edges in the complementary of MST($X$) are the edges that create a PH¹ class [75], possibly with zero persistence – in which case this class dies immediately with a triangle of same diameter. More recently, Koyama et al. [76] showed that the edges of RNG($X$) \ MST($X$) are exactly the edges that create a PH¹ class of positive persistence. A corollary is that the quantity $|RNG(X) \setminus MST(X)|$ is exactly the number of PH¹ classes with positive persistence (see Fig. 4). Koyama et al. also introduced *Euclidean PH¹*, a method that specializes the PH¹ computation for Rips filtrations of point clouds in $\mathbb{R}^2$ and $\mathbb{R}^3$, using the above geometric result. Compared to other Rips PH algorithm like *Ripser*, they improve both running time and memory consumption in this specific case thanks to the construction from the RNG of a reduced Vietoris-Rips complex that features far fewer triangles than the original one, but that has the same PH⁰ and PH¹. This permits to get a much smaller boundary matrix to be reduced.

## E. Topological Autoencoders

Our method builds on *Topological Autoencoders* (TopoAE) [1], a dimension reduction technique that combines an autoencoder with a topological feature preservation loss. We therefore briefly summarize this approach below.
<figure><img src="image_7.png" /><figcaption>Fig. 6. An autoencoder with reconstruction term $\mathcal{L}_r$ and topological regularization term $\mathcal{L}_t$. Here both the encoder and the decoder are fully-connected networks with 1 hidden layer.</figcaption></figure>

Let $\mathcal{X}$ be the input, high-dimensional space and $\mathcal{Z}$ the latent, low-dimensional space. For visualization purposes, in the following we take $\mathcal{Z} = \mathbb{R}^2$. An autoencoder [16] is a composition $\text{dec} \circ \text{enc}$ of an encoder $\text{enc} : \mathcal{X} \to \mathcal{Z}$ and a decoder $\text{dec} : \mathcal{Z} \to \mathcal{X}$. We write $Z = \text{enc} X \in \mathcal{Z}$ the representation of an input vector $X \in \mathcal{X}$ and $\tilde{X} = \text{dec}(\text{enc} X)$ its reconstruction. In practice, enc and dec are neural networks, e.g., with one or several fully connected layer(s), possibly with one or several additional convolutional layer(s) if $\mathcal{X}$ represents an image space. The reconstruction loss of an autoencoder is simply the MSE between an input and its reconstruction:

$$ \mathcal{L}_r = \| \text{dec}(\text{enc} X) - X \|^2. \quad (4) $$

For constraining autoencoders to preserve topological structures, one can add to this reconstruction term a topological regularization term $\mathcal{L}_t$ which compares some topological abstraction of the input and its low-dimensional representation (Fig. 6). Moore et al. [1] discussed a regularization term considering directly the lengths of carefully selected edges:

$$ \mathcal{L}_{\text{TAE}}^{d} = \|A^{X}[\pi^{d}(X)] - A^{Z}[\pi^{d}(X)]\|_{2}^{2} + \|A^{Z}[\pi^{d}(Z)] - A^{X}[\pi^{d}(Z)]\|_{2}^{2}, \quad (5) $$

where $A^X[\cdot]$ and $A^Z[\cdot]$ denote evaluating the distance matrix of $X$ and $Z$ respectively, and where $\pi^d(X)$ (resp. $\pi^d(Z)$) denotes all $k$-critical edges, for any $0 \le k \le d$, in $X$ (resp. $Z$), i.e., the edges that either create or kill a $\text{PH}^k$ class of positive persistence (see Sec. III-B).

In practice however, the method is only documented for $\text{PH}^0$, so that $\pi^0(X) = \text{MST}(X)$ and $\pi^0(Z) = \text{MST}(Z)$ (see Sec. III-D). In that case, the regularization term rewrites as:

$$ \mathcal{L}_{\text{TAE}}^{0} = \|A^{X}[\text{MST}(X)] - A^{Z}[\text{MST}(X)]\|_{2}^{2} + \|A^{Z}[\text{MST}(Z)] - A^{X}[\text{MST}(Z)]\|_{2}^{2}. \quad (6) $$

The idea of this loss function is to force the length of topologically critical edges in the input to be preserved in the representation (first line of Eq. 6), and reciprocally (second line). More precisely, the second term tries to push away in $\mathcal{Z}$ the ends of the edges in $\text{MST}(Z) \setminus \text{MST}(X)$ that are, intuitively, too short in $\mathcal{Z}$. The hope is that at the end of the optimization, the length of both MST's edges match in high and low dimension, inducing the same 0-th persistence diagrams (which we show in Sec. IV-A).

<figure><img src="image_8.png" /><figcaption>Fig. 7. Notations for the proof of Lem. 1. In this example $\text{MST}(X) = \{\{v_1, v_2\}, \{v_2, v_3\}, \{v_3, v_4\}\}$ (in red) and $\text{MST}(Z) = \{\{v_1, v_2\}, \{v_2, v_3\}, \{v_2, v_4\}\}$ (in blue). Therefore $\mathcal{K} = \{\{v_1, v_2\}, \{v_2, v_3\}, \{v_2, v_4\}, \{v_3, v_4\}\}$ (in yellow). Here $\mathcal{L}_{\text{TAE}}^{0}(X, Z)$ evaluates as $(\delta_X(\{v_2, v_4\}) - \delta_Z(\{v_2, v_4\}))^2 + (\delta_X(\{v_3, v_4\}) - \delta_Z(\{v_3, v_4\}))^2 = 2(\sqrt{2} - 1)^2$. Note that this provides an example where $\mathcal{L}_{\text{TAE}}^{0}(X, Z) > 0$ but $\mathcal{D}_{\text{Rips}}^{0}(X) = \mathcal{D}_{\text{Rips}}^{0}(Z)$, i.e. the converse of the stated implication (Lem. 1) is not true.</figcaption></figure>

Finally, the weights of the autoencoder are optimized iteratively to minimize the value of the loss function. Note that this involves computing the required critical edges of the current representation $Z$ at each iteration.

## IV. TOPOLOGICAL AUTOENCODER ANALYSIS

In this section we provide a novel theoretical analysis of the TopoAE loss function [1] and its link with the Wasserstein distance between the persistence diagrams in high and low dimensions. Our aim is to better understand its limitations and the configurations that may cause them.

### A. Relation to the Wasserstein distance

We first show below that the topological regularization term $\mathcal{L}_{\text{TAE}}^{0}$ (Eq. 6) upper bounds the $L_2$-Wasserstein distance $W_2$ between both 0-dimensional persistence diagrams.

**Lemma 1.** For any point clouds $X$ and $Z$ of equal size, we have the following inequality:

$$ W_2(\mathcal{D}_{\text{Rips}}^{0}(X), \mathcal{D}_{\text{Rips}}^{0}(Z))^{2} \le \mathcal{L}_{\text{TAE}}^{0}(X, Z). $$

Besides, under general position hypothesis (unique pairwise distances), if $\mathcal{L}_{\text{TAE}}^{0}(X, Z) = 0$, then $\mathcal{D}_{\text{Rips}}^{0}(X) = \mathcal{D}_{\text{Rips}}^{0}(Z)$ and the $\text{PH}^0$ pairs are the same, i.e. $\text{MST}(X) = \text{MST}(Z)$.

*Proof.* Let $\mathcal{K}$ be the (abstract) 1-dimensional simplicial complex $\mathcal{K} = \text{MST}(X) \cup \text{MST}(Z)$ (see Fig. 7). Let $\delta_X : \mathcal{K} \to \mathbb{R}$ (resp. $\delta_Z : \mathcal{K} \to \mathbb{R}$) be the simplex diameter function in $X$ (resp. in $Z$). Because the minimum spanning tree contains exactly all the edges that kill a 0-dimensional persistent homology class, and because $\text{MST}(X) \subset \mathcal{K}$ and $\text{MST}(Z) \subset \mathcal{K}$, we have $\mathcal{D}_{\text{Rips}}^{0}(X) = \mathcal{D}^{0}(\delta_X)$ and $\mathcal{D}_{\text{Rips}}^{0}(Z) = \mathcal{D}^{0}(\delta_Z)$. Hence,
the following inequalities:

$$
\begin{align*}
W_2(\mathcal{D}_{\text{Rips}}^0(X), \mathcal{D}_{\text{Rips}}^0(Z))^2 &= W_2(\mathcal{D}^0(\delta_X), \mathcal{D}^0(\delta_Z))^2 \\
&\le \sum_{\substack{\sigma \in \mathcal{K} \\ \dim(\sigma) \in \{0,1\}}} |\delta_X(\sigma) - \delta_Z(\sigma)|^2 \quad \text{(by Cellular Wasserstein stability theorem [67])} \\
&= \sum_{\substack{\sigma \in \mathcal{K} \\ \dim(\sigma)=1}} |\delta_X(\sigma) - \delta_Z(\sigma)|^2 \quad \text{(Rips filtration)} \\
&\le \sum_{\substack{\sigma \in \text{MST}(X) \\ \dim(\sigma)=1}} |\delta_X(\sigma) - \delta_Z(\sigma)|^2 + \sum_{\substack{\sigma \in \text{MST}(Z) \\ \dim(\sigma)=1}} |\delta_X(\sigma) - \delta_Z(\sigma)|^2 \\
&= \|A^X[\pi^0(X)] - A^Z[\pi^0(X)]\|_2^2 + \|A^Z[\pi^0(Z)] - A^X[\pi^0(Z)]\|_2^2 \\
&= \mathcal{L}_{\text{TAE}}^0(X, Z).
\end{align*}
$$

In particular, when $\mathcal{L}_{\text{TAE}}^0(X, Z) = 0$, then $\mathcal{D}_{\text{Rips}}^0(X) = \mathcal{D}_{\text{Rips}}^0(Z)$. In that case, general position hypothesis (unique pairwise distances) on $X$ and $Z$ – that guarantees the uniqueness of their respective MST – implies $\text{MST}(X) = \text{MST}(Z)$, hence the same $\text{PH}^0$ pairs. $\square$

This result is, to our knowledge, the first theoretical guar-
antee that justifies using the original, 0-dimensional version of
the TopoAE loss function (Eq. 6) as a regularization term for
preserving PH<sup>0</sup>. In particular, if this term is zero, we know
that the representation Z has the same PH<sup>0</sup> – and therefore
the same 0-th persistence diagram – than the input X.

B. *Counter-example for higher dimensional PH*

Unfortunately, Lem. 1 does not generalize to higher-
dimensional PH. Indeed, as soon as we add the edges that
create or destroy PH¹ classes (Sec. III-E), a zero L¹TAE loss
(Eq. 5) no longer implies equality of the persistence diagrams:

$$
\mathcal{L}_{\text{TAE}}^{1}(X, Z) = 0 \iff \left\{ \begin{array}{l} \mathcal{D}_{\text{Rips}}^{0}(X) = \mathcal{D}_{\text{Rips}}^{0}(Z) \\ \mathcal{D}_{\text{Rips}}^{1}(X) = \mathcal{D}_{\text{Rips}}^{1}(Z) \end{array} \right. . \quad (7)
$$

Fig. 8 provides a detailed counter-example for two point clouds $X \subset \mathbb{R}^3$ and $Z \subset \mathbb{R}^2$ which illustrates this non-implication. This example shows that preserving the length of topologically critical edges is not enough to guarantee the preservation of $\text{PH}^k$ when $k \ge 1$.

V. CASCADE DISTORTION

This section presents a new penalty term called *cascade distortion* (CD) which addresses the above counter-example. This term favors an isometric embedding of the 2-chains filling persistent 1-cycles (i.e., that are computed thanks to the *PairCells* algorithm), leading to a more accurate geometrical representation of the persistent 1-cycles in low dimensions.

A. Formulation

Let $\mathcal{P}^k$ be the set of $\text{PH}^k$ pairs with a positive persistence, i.e., $\mathcal{P}^k = \{(\tau, \sigma) \in \mathcal{K}^k \times \mathcal{K}^{k+1} | \text{partner}[\tau] = \sigma \text{ and } \delta(\tau) < \delta(\sigma)\}$. We now define the $d$-dimensional global *cascade*

<figure><img src="image_9.png" /><figcaption>Fig. 8. An example for the non-implication of Eq. 7. The left point cloud Z is in $\mathbb{R}^2$ and is a deformed regular hexagon; the right one X is in $\mathbb{R}^3$ and its leftmost point has been rotated with respect to the dotted yellow axis with angle $\pi/3$ (see Sec. B for raw coordinates). Both RNGs are pictured in black and the length of their edges are 1. The triangle that kills the cycle in X (resp. Z) is shown in blue (resp. red) with its longest edge highlighted in bold. For Z, the remaining cascade edge (see Eq. 8) is also shown (its length is 1.623). The 0-th persistence diagrams are the same: $\mathcal{D}_{\text{Rips}}^0(X) = \mathcal{D}_{\text{Rips}}^0(Z) = \{(0, 1) \text{ with multiplicity } 5, (0, +\infty)\}$; but not the 1-th persistence diagrams: $\mathcal{D}_{\text{Rips}}^1(Z) = \{(1, 1.819)\}$ and $\mathcal{D}_{\text{Rips}}^1(X) = \{(1, 1.732)\}$. However, the critical edges considered in $\mathcal{L}_{\text{TAE}}^1$ (i.e., RNG edges and 1-cycle killing edges, highlighted in yellow) have the same length in both X and Z, therefore $\mathcal{L}_{\text{TAE}}^1(X, Z) = 0$. On the contrary, the two additional edges considered in the 1-dimensional *cascade distortion* (highlighted in purple, see Eq. 9) are longer in Z than in X, hence $\mathcal{L}_{\text{CD}}^1(X, Z) > 0$.</figcaption></figure>

skeleton (GCS<sup>d</sup>) as the set of edges that appear in the cascade
of any PH<sup>k</sup> pair with 0 ≤ k ≤ d:

$$
\text{GCS}^d(X) = \bigcup_{k=0}^{d} \bigcup_{(\tau, \sigma) \in \mathcal{P}^k} \text{cascade}[\sigma]^{(1)}. \quad (8)
$$

In Fig. 5, GCS¹(X) appears in the bottom-right image as the
set of red edges. This leads us to propose the following *d*-dimensional *Cascade Distortion* (CD) loss function, which
looks similar to the original TopoAE loss (Eq. 5), but which
constrains the length of all the edges appearing within the
global cascade of X, instead of its critical edges only:

$$
\mathcal{L}_{\text{CD}}^{d}(X, Z) = \|A^{X}[\text{GCS}^{d}(X)] - A^{Z}[\text{GCS}^{d}(X)]\|_{2}^{2} + \\
\qquad \|A^{Z}[\pi^{d}(Z)] - A^{X}[\pi^{d}(Z)]\|_{2}^{2} \tag{9}
$$

In practice, we focus our study on $\mathcal{L}_{\text{CD}}^1$ for visualization purposes. Intuitively, the role on the first term (top line of Eq. 9) is to favor an isometric embedding of the input cycles and the 2-chains that fill them, so that their shapes are preserved at best. The second term (bottom line of Eq. 9), like in TopoAE, involves pushing away vertices that are too close in $Z$, hence penalizing the cycles that exist in $Z$ but not in $X$. This way, the first term favors a faithful projection of the high-dimensional cycles, while the second term removes spurious low-dimensional cycles from the planar projection. In addition, this loss enables the usage of a dedicated algorithm for the fast computation of $\pi^1(Z)$ at each optimization iteration (Sec. VI).

In practice, GCS¹(X) is computed using a modified version of PairCells, which, similarly to [70], uses a two-step approach. First, the persistence pairs $\mathcal{P}^1$ with a positive persistence are computed with an efficient, dedicated algorithm,
e.g., Ripser [72]. Then, for each pair $(\tau, \sigma) \in \mathcal{P}^1$ in ascending order, we run the *EliminateBoundary* procedure (Alg. 2) from $\sigma$ until the youngest edge of the boundary $\partial(\text{cascade}[\sigma])$ is $\tau$, storing the cascade edges encountered in the meantime. At each step, Youngest($\partial(\text{cascade}[\sigma])$) forms either an apparent pair (with a triangle that is searched at this moment), or a positive persistence pair that was handled previously. Note that it is mandatory to choose the same order as *Ripser* on the triangles with identical diameter (i.e., the reverse colexicographic order). This is more efficient than a naive execution of *PairCells*, since it benefits from the accelerations of *Ripser*.

## B. Relation to TopoAE

The above *Cascade Distortion* term (Eq. 9) can be understood as a generalization of the original TopoAE topological regularization term in that when restricting to $\text{PH}^0$, both formulations are equivalent.

**Lemma 2.** For any point cloud $X$, $\text{GCS}^0(X) = \text{MST}(X)$. Therefore, the 0-dimensional CD and TopoAE regularization terms are equal: $\mathcal{L}_{\text{TAE}}^0(X, Z) = \mathcal{L}_{\text{CD}}^0(X, Z)$.

*Proof.* Let $X$ be a point cloud. First note that for any $\sigma$, we have $\sigma \in \text{cascade}[\sigma]$, hence:

$$ \text{MST}(X) = \bigcup_{(\tau, \sigma) \in \mathcal{P}^0} \{\sigma\} \subset \bigcup_{(\tau, \sigma) \in \mathcal{P}^0} \text{cascade}[\sigma] = \text{GCS}^0(X). $$

Reciprocally, note that for any vertex $\tau$, partner[$\tau$] is an edge that kills a $\text{PH}^0$ class, hence is an MST edge. Then, with line 6 (Alg. 2), cascade[$\sigma$] writes as the sum of $\sigma$ and cascades of edges of the form partner[$\tau$]. Therefore, by recurrence on the edges by increasing length, for any $\sigma \in \text{MST}(X)$, i.e., for any $(\tau, \sigma) \in \mathcal{P}^0$, we have cascade[$\sigma$] $\subset$ MST($X$). Hence $\text{GCS}^0(X) \subset$ MST($X$). In the end, $\text{GCS}^0(X) = \text{MST}(X)$. $\square$

For higher-dimensional PH however, we can only say that $\mathcal{L}_{\text{TAE}}^d(X, Z) \le \mathcal{L}_{\text{CD}}^d(X, Z)$, i.e., the latter is more restrictive than the former. For instance, for $d = 1$, in the counterexample of Fig. 8 where we had $\mathcal{D}_{\text{Rips}}^1(X) \ne \mathcal{D}_{\text{Rips}}^1(Z)$ but $\mathcal{L}_{\text{TAE}}^1(X, Z) = 0$, the 1-dimensional cascade distortion is here indeed positive: $\mathcal{L}_{\text{CD}}^1(X, Z) > 0$. More generally, $\mathcal{L}_{\text{CD}}^1$ additionally constrains the length of the edges of the 2-chains that fill the persistent 1-cycles, hence favoring an isometric embedding of those chains. Intuitively, this can be understood as trying to preserve the whole shape of the 1-cycles, instead of only their birth and death.

In the following, we call TopoAE++ the dimensionality reduction technique obtained by replacing in TopoAE the topological regularization term $\mathcal{L}_{\text{TAE}}^0$ with $\mathcal{L}_{\text{CD}}^1$, and by including the acceleration presented in the next section.

# VI. FAST ALGORITHM FOR RIPS PH IN THE PLANE

Because the computation of the persistence diagram of the low-dimensional representation is executed at each step of the optimization process, it is a critical stage in terms of computation time. Therefore, we propose in this section a

<figure><img src="image_10.png" /><figcaption>Fig. 9. RNG-polygons (green, blue and red polygons) with their minmax length triangulations (gray edges). The arrows show to which critical edges (birth: dashed; death: continuous) is associated each point of the 1-dimensional persistence diagram. Each birth edge is in the RNG and each death edge is the longest edge of the minmax length triangulations of its RNG-polygon. Missing edges from the convex hull are shown in yellow, as well as the remaining polygons that the convex hull forms with the RNG.</figcaption></figure>

new fast algorithm that computes PH for the Rips filtration of a two-dimensional point cloud, which is useful when our embedding space is $\mathbb{R}^2$, e.g., for visualization purpose. This algorithm is purely geometric and involves no matrix reduction step, contrary to other algorithms (e.g. [72]–[74], [76]), which allows it to be significantly faster and less memory-intensive (see Sec. VII-E1 for timing results).

## A. Geometrical results

We rely on the fact that the edges that introduce a cycle are exactly the edges in $\text{RNG}(X) \setminus \text{MST}(X)$ (as previously exploited in the related work [76], see Sec. III-D). Consequently, the number of points in the 1-th persistence diagram is the number of polygons formed by $\text{RNG}(X)$, which is also the cardinal of $\text{RNG}(X) \setminus \text{MST}(X)$. In contrast to previous work [76] (which also exploited boundary matrix reduction), we additionally claim that the critical edges that kill these $\text{PH}^1$ classes can be found efficiently through geometric – rather than algebraic – considerations, namely with a minmax length triangulation of each RNG-polygons (see Fig. 9).

**Lemma 3.** In the plane, an edge that intersects an RNG edge kills no $\text{PH}^1$ class of positive persistence.

*Proof.* Let $ab$ be an edge that intersects an RNG edge $e = p_1p_2$. Suppose that $ab$ is the longest edge of the triangles $abc$ and $abd$, one of which kills a $\text{PH}^1$ class $\gamma$ of which the quadrilateral $acbd$ is a representative 1-cycle (see Fig. 10). Because $\text{Lens}(e) = \emptyset$ (since $e \in \text{RNG}(X)$), $|e| < |ab|$. In addition, $c$ and $d$ are in $\text{Lens}(a, b)$. Therefore, the 2-chain $acp_1 + cbp_1 + bdp_2 + dap_2 + ap_1p_2 + bp_1p_2$, whose boundary is $acbd$, contains only triangles of diameter $< |ab|$. Hence, the class $\gamma$ is killed before the value $|ab|$, therefore not by $ab$. $\square$

Therefore, the edges that kill a $\text{PH}^1$ class of positive persistence can only lie inside either an RNG-polygon, or a polygon formed by the RNG and the convex hull of $X$ (Fig. 9). Then, we have the following lemma for each RNG-polygon.

**Lemma 4.** Let $\Pi$ be a RNG-polygon. Under general position hypothesis (unique pairwise distances), the longest edge $e_{\text{MML}}$,
<figure><img src="image_11.png" /><figcaption>Fig. 10. Notations for proof of Lem. 3. $e = p_1p_2$ (in red) is supposed to be an RNG edge, hence an empty lens. An edge $ab$ (light green) that crosses $e$ is shown to not be able to kill any PH¹ class of positive persistence.</figcaption></figure>

of length noted $\delta_{MML}$, of any minmax length triangulation of $\Pi$ kills a PH¹ class $\gamma$ with positive persistence for which $\Pi$ is a representative.

*Proof.* Let $\Pi$ be a RNG-polygon with point set $P$. Let $e_{MML}$ be the longest edge of the minmax length triangulations of this polygon, which is well-defined with the general position hypothesis. We have the following:

* for diameter threshold $r = \delta_{MML}$, $\text{Rips}_r(P)$ contains a MML triangulation of $\Pi$, which can be seen as a 2-chain whose boundary is $\Pi$;

* for any diameter threshold $r < \delta_{MML}$, $\text{Rips}_r(P)$ contains no 2-chain whose boundary is $\Pi$: otherwise, one could construct a triangulation of $\Pi$ whose longest edge is of length at most $r$, hence smaller than $\delta_{MML}$.

Hence, $e_{MML}$ kills a PH¹ class $\gamma$ with positive persistence for which $\Pi$ is a representative. □

Therefore, in the plane, as the number of RNG-polygons is exactly the number of PH¹ classes with a positive persistence, the edges that destroy the latter are exactly the edges found in the above lemma, i.e., that are the longest of the MML triangulations of an RNG-polygon (see Fig. 9).

Finally, for any RNG-polygon, we state a relation between $\delta_{MML}$ and the length $\delta_{Del}$ of the longest Delaunay edge $e_{Del}$ inside this polygon, enabling the search of $e_{MML}$ among a smaller subset of the polygon's diagonals (see Sec. VI-B).

**Lemma 5.** In the plane, the values $\delta_{MML}$ and $\delta_{Del}$ of an RNG-polygon verify:

$$ \frac{\sqrt{3}}{2} \delta_{Del} \le \delta_{MML} \le \delta_{Del}. $$

*Proof.* Let $\Pi$ be an RNG-polygon and $P$ its point set. Let $e_{Del} \in \text{DEL}(P)$ be the longest Delaunay edge inside $\Pi$, of length $\delta_{Del} = |e_{Del}|$. The right inequality results from the Delaunay triangulation inside $\Pi$ being a 2-chain whose boundary is $\Pi$ and whose maximum diameter is $\delta_{Del}$: hence the associated Rips PH¹ class $\gamma$ is killed at value $\delta_{MML} \le \delta_{Del}$. For the left inequality, because $e_{Del}$ is a Delaunay edge, there exists a circumcircle $C$ to $e_{Del}$ of radius $r(C) \ge |e_{Del}|/2$ that contains no other point of $P$ (see Fig. 11). Now any 2-chain with $\Pi$ as boundary has to contain the center $c(C)$ of that circle. However, the smallest triangles with no vertex within $C$ and that contain $c(C)$ are exactly the equilateral triangles circumscribed by $C$, and their side length is $\sqrt{3}r(C) \ge \frac{\sqrt{3}}{2}|e_{Del}|$.

<figure><img src="image_12.png" /><figcaption>Fig. 11. Notations for the proof of Lem. 5. The RNG-polygon $\Pi$ (bold black), is such that its Delaunay triangulation (in gray) is not a MML triangulation. Hence, its cycle killing edge $e_{MML}$ (bold red) is shorter than its longest Delaunay edge $e_{Del}$ (blue), i.e., $\delta_{MML} < \delta_{Del}$. The circumcircle $C$ of the longest Delaunay edge $e_{Del}$ contains no points of $X$ in its interior.</figcaption></figure>

Therefore, the associated Rips PH¹ class $\gamma$ is killed by a triangle $\tau$ of diameter at least $\frac{\sqrt{3}}{2}|e_{Del}|$, hence $\frac{\sqrt{3}}{2}\delta_{Del} \le \delta_{MML}$. □

## B. Algorithm description

Leveraging the above results, we introduce a fast algorithm that computes PH for the Rips filtration of a 2-dimensional point cloud. It is summarized in Alg. 3 and we describe below more precisely the different stages.

<table><caption>**Algorithm 3:** 2D Rips persistence method overview</caption><thead><tr><th>Input:</th><td>Point cloud $X \subset \mathbb{R}^2$</td></tr><tr><th>Output:</th><td>Persistence diagrams $\mathcal{D}_{\text{Rips}}^0(X)$ and $\mathcal{D}_{\text{Rips}}^1(X)$ with associated simplices</td></tr></thead><tbody><tr><td>1</td><td>Compute the Delaunay triangulation $\text{DEL}(X)$</td></tr><tr><td>2</td><td>Compute the Urquhart graph $\text{UG}(X)$</td></tr><tr><td>3</td><td>Apply Kruskal's algorithm on $\text{UG}(X)$ to compute $\text{MST}(X)$ and obtain $\mathcal{D}_{\text{Rips}}^0(X)$</td></tr><tr><td>4</td><td>Check the emptiness of the lens of the edges in $\text{UG}(X) \setminus \text{MST}(X)$ to compute $\text{RNG}(X)$</td></tr><tr><td>5</td><td>For each RNG-polygon, find its local minmax length triangulation longest edge (Alg. 4)</td></tr><tr><td>6</td><td>Apply Kruskal's algorithm on the dual $\text{RNG}(X)$ graph to obtain $\mathcal{D}_{\text{Rips}}^1(X)$</td></tr></tbody></table>

a) **Delaunay triangulation:** We use CGAL [77] to compute the Delaunay triangulation, with identifiers on vertices and triangles. This implementation deals with the boundary (simplices on the convex hull) with a virtual point at infinity.

b) **Relative neighborhood graph:** To compute $\text{RNG}(X)$, we first compute $\text{UG}(X)$ by keeping Delaunay edges whose link (made of two points) is outside their lens. We then apply Kruskal's algorithm (which uses a union-find data structure [15]) on $\text{UG}(X)$ to compute $\text{MST}(X)$, of which
**Algorithm 4:** Find the longest edge of a minmax length triangulation of an RNG-polygon of *X*

**Input:** Points $P \subset X$ of an RNG-polygon
**Input:** Length $\delta_{\text{Del}}$ of the longest Delaunay edge
**Output:** Edge that destroy the associated 1-cycle
1 candidates $\leftarrow \{\}$
2 **for** $e = v_1v_2 \in \binom{P}{2}$ **s.t.** $|e| \in [\sqrt{3} \delta_{\text{Del}}, \delta_{\text{Del}}]$ **do**
3   // check if $e$ is a 2-edge
4   **if** LeftLens($e$) $\neq \emptyset$ **and** RightLens($e$) $\neq \emptyset$ **then**
5     // check if $e$ is expandable
6     **for** $x \in$ LeftLens($e$) **do**
7       **if** RightLens($v_1x$) $= \emptyset$ **and**
8         RightLens($xv_2$) $= \emptyset$ **then**
9         LeftExpandable $\leftarrow$ true
10        break
11   **for** $y \in$ RightLens($e$) **do**
12     **if** RightLens($v_2y$) $= \emptyset$ **and**
13       RightLens($yv_1$) $= \emptyset$ **then**
14         RightExpandable $\leftarrow$ true
15        break
16   **if** LeftExpandable **and** RightExpandable **then**
17     candidates $\leftarrow$ candidates $\cup \{e\}$
18 **return** $\min(\text{candidates}, |\cdot|)$

$D_{\text{Rips}}^0(X)$ can be deduced (see Sec. III-A and Sec. III-D). Then, checking the emptiness of the lens of the edges in $UG(X) \setminus \text{MST}(X)$ using disk queries on a precomputed k-d tree, we determine whether they belong to $\text{RNG}(X)$, and construct a data structure that represents $\text{RNG}(X)$ and its polygons (see next paragraph).

c) *Data structure for RNG-polygons:* During the computation of $UG(X)$ and $\text{RNG}(X)$, we maintain a union-find data structure on Delaunay triangles that permits to deduce efficiently the RNG-polygons. More precisely, as soon as we encounter a non-RNG edge (either during $UG(X)$ computation or during $\text{RNG}(X)$ refinement), we merge (if not already merged) the two classes associated with its two adjacent triangles. We also keep for each current class the longest deleted Delaunay edge, which needs to be maintained at each merge. Besides, we store $UG(X)$ and $\text{RNG}(X)$ edges as quad-edges (i.e., an edge is stored as the identifiers of its two extremities and the identifiers of its two adjacent triangles). This allows an implicit representation of the RNG-polygons, each of which is associated with a root of the union-find data structure over triangles. For the next stage, we also recover the vertices associated with each polygon.

d) *Minmax length triangulations:* Now all we need is to find the edges that kill persistent cycles: using Lem. 4, we search the longest edge $e_{\text{MML}}$ of a MML triangulation of each polygon. For that we partially follow the algorithm described in [61] that computes MML triangulations in quadratic time, even though there is no need here to compute a full triangulation since we only need its longest edge. Therefore

<figure><img src="image_13.png" /><figcaption>Fig. 12. Illustration of the expandability – tested by Alg. 4 – of the 2-edge $e = v_1v_2$. Indeed there exists $x, y \in \text{Lens}(e)$ (red area) on either side of $e$ such that the “inner” half-lenses (green areas) of the edges $v_1x, xv_2, v_2y, yv_1$ are empty. Thus, there exists a planar triangulation of $\Pi$ whose longest edge is $e$. Besides, $e$ happens to be the smallest of the expandable 2-edges in $\Pi$: $e$ is therefore the longest edge of the MML triangulations of $\Pi$.</figcaption></figure>

for each RNG-polygon, we enumerate the polygon diagonals with length between $\frac{\sqrt{3}}{2} \approx 0.866$ and 1 times the length of the longest deleted Delaunay edge within the polygon (see Lem. 5). Among them, we search the smallest that is an *expandable 2-edge* (a property introduced in [61] that allows an edge to be the longest of a triangulation of the polygon, see Sec. A-B for more details) using disk queries on a k-d tree on the vertices of the polygon (see Fig. 12 and Alg. 4).

e) *Persistent pairs:* The last step consists in pairing edges that create cycles with those that kill them (a cycle-creating edge is not necessarily associated with a killer edge of an adjacent polygon, in case of nested cycles). We solve this relying on duality, an overall strategy also employed in the seminal persistence algorithm [65]. Specifically, for our setup, we compute $\text{PH}^0$ of the dual of $\text{RNG}(X)$ with opposite filtration values, which is the $\text{PH}^1$ of the input. In other words, we apply Kruskal’s algorithm on a dual graph where dual nodes are RNG-polygons $\Pi$ with filtration value $-\delta_{\text{MML}}(\Pi)$, and where dual edges are cycle-creating edges $e \in \text{RNG}(X) \setminus \text{MST}(X)$ with filtration value $-|e|$.

### C. Parallelization

Several stages can be parallelized for accelerating the running time. The most critical stage is the computation of minmax length triangulations (it represents around 50% of the computation time on average in sequential for uniform point clouds). It can be naively parallelized on the different RNG-polygons. However, these polygons can be of varying sizes, leading to a poor speed-up in some cases. Instead, we develop a slightly more efficient parallelization, which consists in storing in a single array every candidate edge (of length in $[\frac{\sqrt{3}}{2}\delta_{\text{Del}}, \delta_{\text{Del}}]$) of every polygon, and to parallelize the expandability checks on this array.

Finally, other stages can be easily parallelized with a less noticeable effect, such as using a parallel sorting implementa-
tion for sorting the UG($X$) edges by length (before applying Kruskal's algorithm on them to get MST($X$)).

## VII. RESULTS

This section presents both qualitative and quantitative experimental results. All results were obtained on a laptop computer with 64 GB of RAM, a Core i7-13850HX CPU (8 cores at 5.3 GHz and 12 cores at 3.8 GHz), and a RTX 2000 GPU. We implemented our method in C++ (with the libraries OpenMP, LibTorch [78] and CGAL [77]) as modules for TTK [79], [80]. We used as encoder (resp. decoder) a fully connected network with two hidden layers of size 128 and 32 (resp. 32 and 128), with ReLU activation functions, and batch normalization. For each experiment, we performed 1000 iterations of the Adam optimizer [81] with a learning rate fixed at 0.01.

### A. Test data

We present results on synthetic, three-dimensional datasets in Fig. 1, which have been specifically designed to include clearly salient topological features. The first one consists in 3 Gaussian clusters, and features 3 significantly persistent $PH^0$ pairs. The second one consists in points sampled around an ellipse twisted along its major axis, and features one significantly persistent $PH^1$ pair. The third one consists in points sampled around the edges of a tetrahedron (i.e., the complete $K_4$ graph which is planar) and features 3 significantly persistent $PH^1$ pairs. We also consider a stress case in Fig. 13, with a dense point cloud sampling the complete $K_5$ graph embedded in $\mathbb{R}^3$ (i.e., 5 vertices in 3D, all pairwise connected by an edge). Since $K_5$ is a non-planar graph, the resulting 3D persistent generators (which follow the edges of $K_5$) cannot be projected to the plane without intersections. Then, it is not possible to project this stress case to the plane faithfully with respect to $PH^1$, as the intersections of the projected 3D generators create additional, short cycles in 2D.

We present results on real-life, acquired or simulated, high-dimensional, cycle-featuring datasets in Fig. 14. COIL-20 [82] is an image dataset with 20 classes, each one containing 72 grayscale 64 × 64 images. Each class consists in a single object (e.g., a rubber duck) viewed from 72 evenly spaced angles, so that the images within each class live on a manifold that is homeomorphic to a circle. This leads to significantly persistent pairs when computing the $PH^1$ of the input endowed with the Euclidean metric between the images. Cyclic features may also appear in motion capture data [83]. Hence, we also evaluate our approach on a dataset extracted from the CMU motion capture database [84], which records a human subject doing some periodic movement (e.g., walking or running), each frame consisting in 62 skeleton joint angles measurements. Finally, data from biological processes can also feature meaningful cycles like in single-cell omics [85], a field that aims at sequencing single cells. This enables the observation of gene expression within individual cells at different stage of the cell cycle (resting, growth, DNA replication, mitosis). Specifically, we used the simulated dataset "cyclic_2" [86], which features a significantly persistent $PH^1$ pair. See Tab. I for a summary of the presented datasets.

<table><caption>TABLE I
SUMMARY OF DATASETS.</caption><thead><tr><th>Dataset</th><th>Nature</th><th>Dimension (h)</th><th>Size (n)</th><th>Fig.</th></tr></thead><tbody><tr><td>3Clusters</td><td>Synthetic</td><td>3</td><td>800</td><td rowspan="3">Fig. 1</td></tr><tr><td>Twist</td><td>Synthetic</td><td>3</td><td>100</td></tr><tr><td>K_4</td><td>Synthetic</td><td>3</td><td>300</td></tr><tr><td>K_5</td><td>Synthetic</td><td>3</td><td>500</td><td>Fig. 13</td></tr><tr><td>COIL-20-1</td><td>Acquired</td><td>4096</td><td>72</td><td rowspan="3">Fig. 14</td></tr><tr><td>MoCap</td><td>Acquired</td><td>62</td><td>138</td></tr><tr><td>SingleCell</td><td>Simulated</td><td>1170</td><td>243</td></tr></tbody></table>

### B. Quantitative criteria

In order to evaluate the topological accuracy of the representations generated by our approach, we compute the $L_2$-Wasserstein distance (see Sec. III-B) between the 1-dimensional persistence diagrams of the input $X$ and the low-dimensional representation $Z$:

$$PDW^1(X, Z) = W_2(D_{\text{Rips}}^1(X), D_{\text{Rips}}^1(Z)).$$

This quantity, which we want to be as low as possible, summarizes in some measure the extent to which similar cycles exist in the two point clouds. However, as persistence diagrams lose information – they are in particular invariant to permutations in the point clouds – it does not convey whether the set of vertices involved in these cycles are the same, in contrast to TopoAE-like loss functions.

We also compute the metric distortion $\mathcal{D}$, which should also be as low as possible, expressed as the root mean squared error (RMSE) between the pairwise distances in $X$ and those in $Z$ – which is what MDS-like methods typically minimize:

$$\mathcal{D}(X, Z) = \sqrt{\frac{1}{n} \sum_{i<j} (\|X_i - X_j\|_2 - \|Z_i - Z_j\|_2)^2}.$$

We refer the reader to Sec. C-A for a description of other quality scores often considered in the DR literature (see Sec. C-B for additional experiments with these indicators).

### C. Competing approaches

We compare our approach to global methods (PCA [2], MDS [3]) that do not take topology into account. We also compare to several locally topology-aware methods (Isomap [4], t-SNE [5], UMAP [6]) that are able to consider the local manifold structure of the data but that might fail to project faithfully the overall structure. Regarding previously documented globally topology-aware methods, although TopoMap [7] preserves exactly $PH^0$ and TopoAE [1] preserves $PH^0$ in the sense of Lem. 1, both methods have no guarantee about $PH^1$. In addition, we show the results of the method suggested in [8], which we denote TopoAE+$W^1$ and which consists in adding $PDW^1(X, Z)$ to the TopoAE loss $L_{\text{TAE}}^0$. Both terms in the loss are given the same weights. It is then supposed to give the best results in terms of $PDW^1(X, Z)$ values and therefore make the cycles in $X$ and $Z$ as similar (i.e., in number, in size) as possible. However, as explained in Sec. VII-B, $PDW^1(X, Z)$ is not sensible to which simplices are involved in the cycles, and the cycles in $Z$ might be created between simplices that do not belong to a cycle in $X$.
<table><caption>Fig. 13. Comparison of DR methods on a stress case, i.e., a point cloud sampled around the edges of the complete $K_5$ graph embedded in 3D. It features 6 significantly persistent $PH^1$ pairs, and a generator for each is represented with a curve of distinct color. As $K_5$ is a non-planar graph, this point cloud cannot be projected faithfully in 2D in a topological sense (i.e., the high-dimensional persistent generators cannot be projected to the plane without intersections). Despite this, our method (TopoAE++) achieves a competitive Wasserstein distance.</caption><thead><tr><th rowspan="2">Input</th><th colspan="2">Global methods</th><th colspan="4">Locally topology-aware methods</th><th colspan="4">Globally topology-aware methods</th></tr><tr><th>PCA [2]</th><th>MDS [3]</th><th>Isomap [4]</th><th>t-SNE [5]</th><th>UMAP [6]</th><th>TopoMap [7]</th><th>TopoAE [1]</th><th>TopoAE+W<sup>1</sup> [8]</th><th>TopoAE++</th></tr></thead><tbody><tr><td>K_5 (n = 500)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PDW<sup>1</sup>(X, Z)</td><td>2.8e-01</td><td>3.5e-01</td><td>1.3e-01</td><td>4.7e+01</td><td>9.7e-01</td><td>4.8e-01</td><td>1.6e-01</td><td><b>2.2e-03</b></td><td>8.0e-02</td></tr><tr><td>D(X, Z)</td><td>3.1e-01</td><td><b>2.4e-01</b></td><td>4.7e-01</td><td>2.1e+01</td><td>9.8e+00</td><td>1.2e+00</td><td>9.2e-01</td><td>8.6e-01</td><td>3.4e-01</td></tr><tr><td>t(s)</td><td>0.00043</td><td>0.59</td><td>0.088</td><td>0.63</td><td>0.14</td><td>0.063</td><td>29</td><td>170</td><td>6.0</td></tr></tbody></table>

<table><caption>Fig. 14. Comparison of DR methods on three real-life, acquired or simulated, high-dimensional datasets (see Sec. VII-A for a description) that all feature a single significantly persistent $PH^1$ pair. A generator for this pair is projected down to 2D, in transparency and slightly smoothed (with a color map depicting its arc-length parameterization). Quantitatively, our approach (TopoAE++) achieves a competitive Wasserstein distance. Qualitatively, it generates fewer crossings of the projected high-dimensional persistent generator.</caption><thead><tr><th rowspan="2">Input</th><th colspan="2">Global methods</th><th colspan="4">Locally topology-aware methods</th><th colspan="4">Globally topology-aware methods</th></tr><tr><th>PCA [2]</th><th>MDS [3]</th><th>Isomap [4]</th><th>t-SNE [5]</th><th>UMAP [6]</th><th>TopoMap [7]</th><th>TopoAE [1]</th><th>TopoAE+W<sup>1</sup> [8]</th><th>TopoAE++</th></tr></thead><tbody><tr><td>COIL-20-1 (n = 72)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PDW<sup>1</sup>(X, Z)</td><td>1.2e+01</td><td>1.0e+01</td><td>6.0e+02</td><td>1.3e+01</td><td>1.2e+01</td><td>2.7e+01</td><td>7.5e+00</td><td>9.1e-02</td><td><b>1.6e-02</b></td></tr><tr><td>D(X, Z)</td><td>2.9e+00</td><td><b>2.0e+00</b></td><td>2.1e+01</td><td>4.9e+00</td><td>6.2e+00</td><td>7.2e+00</td><td>8.2e+00</td><td>6.7e+00</td><td>1.9e+01</td></tr><tr><td>t(s)</td><td>0.0071</td><td>0.086</td><td>0.083</td><td>0.30</td><td>2.5</td><td>0.14</td><td>2.0</td><td>4.0</td><td>2.4</td></tr><tr><td>MoCap (n = 138)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PDW<sup>1</sup>(X, Z)</td><td>2.1e+03</td><td>1.9e+03</td><td>1.0e+04</td><td>1.6e+03</td><td>1.5e+03</td><td>2.1e+03</td><td>1.2e+03</td><td><b>9.4e+00</b></td><td>1.2e+01</td></tr><tr><td>D(X, Z)</td><td>2.4e+01</td><td><b>1.7e+01</b></td><td>1.3e+02</td><td>9.9e+01</td><td>1.0e+02</td><td>8.5e+01</td><td>2.6e+01</td><td>5.6e+01</td><td>2.9e+01</td></tr><tr><td>t(s)</td><td>0.041</td><td>0.25</td><td>0.11</td><td>0.33</td><td>2.6</td><td>0.056</td><td>3.1</td><td>12</td><td>3.4</td></tr><tr><td>SingleCell (n = 243)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PDW<sup>1</sup>(X, Z)</td><td>3.0e+03</td><td>3.1e+03</td><td>1.8e+05</td><td>1.8e+03</td><td>1.8e+03</td><td>2.8e+04</td><td>1.8e+03</td><td><b>6.4e+02</b></td><td>8.6e+02</td></tr><tr><td>D(X, Z)</td><td>3.1e+01</td><td><b>2.3e+01</b></td><td>4.8e+02</td><td>9.6e+01</td><td>1.0e+02</td><td>1.2e+03</td><td>9.8e+01</td><td>3.2e+02</td><td>1.5e+02</td></tr><tr><td>t(s)</td><td>0.0089</td><td>0.77</td><td>0.019</td><td>0.33</td><td>2.5</td><td>0.023</td><td>7.6</td><td>27</td><td>6.0</td></tr></tbody></table>

For PCA, MDS, Isomap and t-SNE, we used the implementations in scikit-learn [87] (with default parameters except Isomap for which the number of neighbors is set to 8); for UMAP, the Python package umap-learn provided by the authors; and for TopoMap, the existing implementation in TTK. For TopoAE-like approaches, we used our own implementation, with Ripser [72] as latent space $PH^1$ algorithm for TopoAE and TopoAE+$W^1$, and Alg. 3 for TopoAE++. For non-deterministic methods (i.e., TopoAE, TopoAE+$W^1$ and TopoAE++ which rely on a stochastic optimization of a neural network with a random initialization), we performed 10 runs and kept the result with the best $PDW^1(X, Z)$.

## D. Result analysis

1) **Qualitative analysis:** Our visual validation approach involves computing a generator for each significantly persistent pair in the input and observing their low-dimensional projections. These generators should ideally be projected without self-intersection and coincide with the cycles in the low-dimensional representation. In such a case, the output would be considered to faithfully represent the high-dimensional cycles.

In the depicted examples, TopoAE++ succeeds in projecting the cycles from the input in this visual sense, whenever

possible (with the exception of $K_5$, for which a faithful planar projection is not possible, see Sec. VII-A). It also tends to preserve the *shape* of the cycles. For example, *Twist* (Fig. 1) is projected with the original ellipse shape preserved, but without torsion. Similarly, in *COIL-20-1* (Fig. 14), the high-dimensional cycle has a stretched shape (since the rubber duck has a similar silhouette when viewed from the front and back, but not when viewed from the left and right), which TopoAE++ takes into account (views from the left and right corresponds to the "ends" of the stretched planar cycle). *MoCap* (Fig. 14) comes from the capture of a subject performing a translation in addition to its periodic movement (running). It is therefore projected to a spiral shape, where the longest edge of the projection of the input cycle (appearing in lavender) corresponds to the same body configuration during the running cycle, but at two different locations. Finally, *SingleCell* (Fig. 14) is difficult to project in the plane as it presents points that lie on a cycle but are distant from each other due to the high dimensionality: TopoAE++ takes that into account and its result can be interpreted as a compromise between preserving $PH^0$ and $PH^1$.

On contrary, global methods (PCA, MDS) do not guarantee the preservation of these cycles (see, e.g., their results on
Twist, Fig. 1, and COIL-20-1, Fig. 14), while locally topology-aware methods sometimes fail to project them correctly (e.g., t-SNE and UMAP break apart the cycles in K_4, Fig. 1, and SingleCell, Fig. 14, when projecting). In addition, even if the cycle is visually correctly projected, these methods tend to forget the shape of the cycle (e.g., Isomap projects Twist, Fig. 1, and COIL-20-1, Fig. 14, as almost perfect circles). Finally, as expected, previously documented globally topology-aware methods generally fail to faithfully project the cycles, as they only incorporate constraints on $PH^0$ (TopoMap, TopoAE) or on the Wasserstein distance $PDW^1$ (TopoAE+$W^1$), but not on $PH^1$ (see the next paragraph).

2) Quantitative analysis: From a more quantitative point of view, in the examples shown, TopoAE++ generally produces the second best $PDW^1$ after $TopoAE+W^1$, both significantly better than all other methods for this metric. However, as $TopoAE+W^1$ directly optimizes $PDW^1$ and forgets the simplices involved in $PH^1$ pairs, the vertices involved in planar cycles (i.e., corresponding to the points of $D_{\text{Rips}}^1(Z)$) do not match those involved in input cycles (i.e., corresponding to the points of $D_{\text{Rips}}^1(X)$). This leads to a good $PDW^1$ metric but a poor projection. For instance, in Twist (Fig. 1), $TopoAE+W^1$ creates a cycle (the top one) that has the right birth and death – while the bottom loop has a persistence close to zero – to cancel out $PDW^1$; yet the original cycle is poorly projected, with a self-intersection. Similarly, on K_5 (Fig. 13), $TopoAE+W^1$ places the points in order to have 6 persistent cycles with the right birth and death to cancel out $PDW^1$, while $TopoAE++$ tries to preserve $PH^1$, which is impossible here due to the non-planarity of $K_5$.

The distortion metric is higher with $TopoAE++$, than in global methods (PCA, MDS) which directly minimize it (by design for MDS). It can be seen as a necessary compromise to faithfully project cycles (e.g., Twist, Fig. 1, requires distortion to be untwisted in 2D; the cycle in COIL-20-1, Fig. 14, which is folded up in high dimension, requires distortion to be unfolded). See Fig. 19 of Sec. C for additional metrics commonly used in DR, along with a companion discussion.

Fig. 15 provides further evaluations of the internal aspects of our approach. Specifically, for the counter-example of

<figure><img src="image_14.png" /><figcaption>Fig. 15. Topological accuracy comparison between projections obtained with the loss $L_{\text{TAE}}^1$ (a naive extension of [1] to $PH^1$, Eq. 5) and our cascade distortion ($L_{\text{CD}}^1$, Eq. 9) on the counter-example $X \in \mathbb{R}^3$ of Fig. 8, producing the 2D embeddings $Z_{\text{TAE}}$ and $Z_{\text{CD}}$ respectively (left). In this example, the optimization achieved a virtually zero value for both losses. However, in the space of persistence diagrams, $D_{\text{Rips}}^1(Z_{\text{TAE}})$ does not match $D_{\text{Rips}}^1(X)$, while $D_{\text{Rips}}^1(Z_{\text{CD}})$ coincides with it (right).</figcaption></figure>

Fig. 8, it compares projections obtained by minimizing a naive extension of [1] to $PH^1$ ($L_{\text{TAE}}^1$, Eq. 5) and our novel loss ($L_{\text{CD}}^1$, Eq. 9). In this example, while the optimization achieves a virtually zero value for both losses, a clear gap occurs between the persistence diagrams $D_{\text{Rips}}^1(Z_{\text{TAE}})$ and $D_{\text{Rips}}^1(X)$, while $D_{\text{Rips}}^1(Z_{\text{CD}})$ coincides with $D_{\text{Rips}}^1(X)$. This indicates that $L_{\text{CD}}^1$ enables a better preservation of the persistent homology of the input $Z$. This gain in topological accuracy is further evaluated in Fig. 16, which compares the resulting Wasserstein distances ($PDW^1(X, Z)$) for our real-life datasets. In particular, we observe in this figure that, while both losses do not produce intersections in the projection of the high-dimensional generator (colored curve), projections based on our novel loss improve topological accuracy by 48% on average. This illustrates that our novel loss produces projections which are more faithful to the input data, in terms of the number, size and shape of its cyclic patterns.

## E. Time performance

1) Fast 2D persistence computation: We compared our persistence algorithm (Alg. 3) for point clouds in $\mathbb{R}^2$ to the following more generic Rips $PH^0$ and $PH^1$ implementations: Gudhi [73] (a generic implementation that handles various filtrations), Ripser [72] (a specialization to Rips filtrations that can be applied to any point cloud in $\mathbb{R}^d$ or distance matrix), and Euclidean $PH1$ [76] (that works theoretically for any point cloud in $\mathbb{R}^d$, but implemented only for $d=2$ or $d=3$). We observe that on uniformly sampled point clouds, the running time asymptotic behavior is roughly cubic (with regard to the number of input points) for Gudhi and quadratic for Ripser, while it remains close to being linear for Euclidean $PH1$ and our algorithm. However, our approach offers a speedup of about 2 orders of magnitudes compared to Euclidean $PH1$ (see Fig. 17). The speedup seems even more noticeable with points clouds sampled next to a circle which creates one high-persistence pair in the 1-dimensional persistence diagram and particularly slows down Euclidean $PH1$ (see Fig. 18). Finally,

<figure><img src="image_15.png" /><figcaption>Fig. 16. Topological accuracy comparison ($PDW^1(X, Z)$) between projections obtained with the loss $L_{\text{TAE}}^1$ and our cascade distortion $L_{\text{CD}}^1$, for our real-life datasets. For $L_{\text{TAE}}^1$, we report the runtime with the gold standard for Rips persistence computation [72] as well as, in parentheses, with our novel planar Rips persistence computation algorithm (Sec. VI). On average, projections using our novel loss ($L_{\text{CD}}^1$), combined with our fast planar Rips persistence algorithm, improve topological accuracy by 48%, with a runtime improvement of 71%.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 17. Comparison of Rips persistence computation times between *Gudhi*, *Ripser*, *Euclidean PH1* and our algorithm (on a single thread, 1T, and on 20 threads, 20T), on point clouds randomly, uniformly distributed over $[0, 1]^2$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 18. Comparison of Rips persistence computation times between *Gudhi*, *Ripser*, *Euclidean PH1* and our algorithm (on a single thread, 1T, and on 20 threads, 20T), on a stress case (for *Ripser* and *Euclidean PH1*): a point cloud randomly distributed over the unit circle, with Gaussian noise ($\sigma = 0.1$).</figcaption></figure>

we have observed empirically that the use of Lem. 5 in Alg. 4 (i.e., checking expandability only for the edges of length in $[\sqrt{3}2\delta_{\text{Del}}, \delta_{\text{Del}}]$) enables discarding 84% of the possible polygon diagonals on average, which significantly contributes to the reported speedup.

2) Overall DR approach: TopoAE++ can hardly compete with conventional DR methods (PCA, MDS, Isomap, t-SNE) in terms of computation time. However, this must be put into perspective. First, we construct a complete projection from the input space $\mathcal{X}$ to the latent space $\mathbb{R}^2$, and not just a projection of the point cloud $\mathcal{X}$ (i.e., we can project a new point, at the low price of one evaluation of the encoder network). Second, we incorporate constraints from persistent homology, which is computationally expensive by nature.

However, our fast 2D Rips PH algorithm (Sec. VI) manages to overcome the runtime overhead due to the computation of $PH^1$ of $\mathcal{Z}$ at each iteration of the optimization. Indeed, comparing TopoAE++ with TopoAE+$\mathcal{W}^1$ – which also requires this computation, implemented with *Ripser* – we observe a speedup by a factor of 2 to 20 in the presented examples. The experiments from Fig. 16 (comparing TopoAE++ to a naive extension of $[1]$ to $PH^1$) report compatible speedups (3.7 on average). In addition, TopoAE++ is sometimes even faster than TopoAE although the latter only computes $PH^0$, both initially for the input $\mathcal{X}$ and at each optimization step for $\mathcal{Z}$.

## F. Limitations

The first obvious limitation is of a fundamental nature: there are points clouds for which there is no planar embedding that is faithful with regard to the $PH^1$ (see Fig. 13 and its description, Sec. VII-A). Besides, like many other topological methods, we are limited in the size of the input due to the initial PH computation, if only in terms of memory. Specifically, *Ripser*, which is commonly considered as the gold standard implementation for Rips filtrations in high dimensions, runs out of memory on a 64 GB machine when computing the $PH^1$ of around 50,000 points. This limitation could be addressed by batch processing or subsampling approaches (see, e.g., [1]). In addition, as with the original TopoAE, unsatisfactory results may occur when the minimization of the topological regularization term is not successful (i.e., stuck in a “bad” local minimum). In our experiments, as for TopoAE and TopoAE+$\mathcal{W}^1$, this was mitigated by running 10 times the projection with a random initialization, and considering as an output the projection minimizing $PD\mathcal{W}^1(X, \mathcal{Z})$. More sophisticated strategies could be considered in the future, in particular for an early detection of potentially bad local minima. Finally, we have no theoretical guarantee that our cascade distortion loss function upper bounds the Wasserstein distance between the persistence diagrams, despite its experimental ability to accurately preserve the cycles. However, we have not found a counter-example yet, similar to the example presented in Fig. 8, and this question remains open for the moment.

## VIII. CONCLUSION

We have introduced TopoAE++, a dimensionality reduction method aiming at accurately visualizing the cyclic patterns present in high diemnsional data. For this, we revisited TopoAE [1], provided a novel theoretical analysis of its original formulation for $PH^0$, and introduced a new generalization to $PH^1$. We have shown experimentally that our novel projection provides an improved balance between the topological accuracy and the visual preservation of the input 1-cycles. As a side benefit of our work, to overcome the computational overhead due to the $PH^1$ computation, we have presented a novel, fast, geometric algorithm that computes Rips PH for planar point clouds, which may be of independent interest.

In future work, the possibility to handle automatically bad local minima during the minimization, e.g., by penalizing self-intersections of input generators, could be investigated. Besides, an interesting extension would be to generalize the method to either higher-dimensional latent spaces (e.g., $\mathcal{Z} = \mathbb{R}^3$, still for visualization purpose), or constraints on higher-dimensional homology (e.g., constraining $PH^2$ to preserve the *cavities*). However, this might prove challenging in terms of computational cost, since the generalization of our fast planar Rips $PH^1$ algorithm does not seem trivial to such a context. Finally, another interesting extension would be to consider other filtrations of specific interest for the applications.

## ACKNOWLEDGMENTS

This work is partially supported by the European Commission grant ERC-2019-COG “TORI” (ref. 863464, https:
//erc-tori.github.io/); and by ANR-23-PEIA-0004 (PDE-AI
https://pde-ai.math.cnrs.fr/).

REFERENCES

[1] M. Moor, M. Horn, B. Rieck, and K. Borgwardt, "Topological Autoencoders," in International conference on machine learning. PMLR, 2020.

[2] K. Pearson, "On lines and planes of closest fit to systems of points in space," *The London, Edinburgh, and Dublin philosophical magazine and journal of science*, 1901.

[3] W. S. Torgerson, "Multidimensional Scaling: I. theory and method," *Psychometrika*, 1952.

[4] J. B. Tenenbaum, V. d. Silva, and J. C. Langford, "A global geometric framework for nonlinear dimensionality reduction," *Science*, 2000.

[5] L. Van der Maaten and G. Hinton, "Visualizing data using t-SNE," *JMLR*, 2008.

[6] L. McInnes, J. Healy, and J. Melville, "UMAP: Uniform Manifold Approximation and Projection for dimension reduction," *arXiv preprint arXiv:1802.03426*, 2018.

[7] H. Doraiswamy, J. Tierny, P. J. Silva, L. G. Nonato, and C. Silva, "TopoMap: A 0-dimensional homology preserving projection of high-dimensional data," *IEEE TVCG*, 2020.

[8] M. Carriere, F. Chazal, M. Glisse, Y. Ike, H. Kannan, and Y. Umeda, "Optimizing persistent homology based functions," in *ICML*, 2021.

[9] I. Borg and G. P., *Modern Multidimensional Scaling - Theory and Applications*. Springer Series in Statistics, 1997.

[10] J. A. Lee and M. Verleyse, *Nonlinear Dimensionality Reduction*. Springer, 2007.

[11] L. van der Maaten, E. Postma, and J. van den Herik, "Dimensionality reduction: A comparative review," Tilburg University, Tech. Rep., 2007.

[12] Z. Cunningham, J. P. Ghahramani, "Linear dimensionality reduction: Survey, insights, and generalizations." *JMLR*, 2015.

[13] L. G. Nonato and M. Aupetit, "Multidimensional projection for visual analytics: Linking techniques with distortions, tasks, and layout enrichment," *IEEE TVCG*, 2019.

[14] I. Trofimov, D. Cherniavskii, E. Tulchinskii, N. Balabin, E. Burnaev, and S. Barannikov, "Learning topology-preserving data representations," in *ICLR*, 2023.

[15] T. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, *Introduction to Algorithms*. MIT Press, 2009.

[16] G. E. Hinton and R. R. Salakhutdinov, "Reducing the dimensionality of data with neural networks," *Science*, 2006.

[17] S. T. Roweis and L. K. Saul, "Nonlinear dimensionality reduction by locally linear embedding," *Science*, 2000.

[18] V. D. Silva and J. B. Tenenbaum, "Global versus local methods in nonlinear dimensionality reduction," in *NeurIPS*, 2003.

[19] P. Joia, D. Coimbra, J. Cuminato, F. Paulovich, and L. Nonato, "Local affine multidimensional projection," *IEEE TVCG*, 2011.

[20] G. E. Hinton and S. Roweis, "Stochastic neighbor embedding," *Advances in neural information processing systems*, 2002.

[21] E. Amid and M. K. Warmuth, "TriMap: Large-scale dimensionality reduction using triplets," 2022.

[22] Y. Wang, H. Huang, C. Rudin, and Y. Shaposhnik, "Understanding how dimension reduction tools work: an empirical approach to deciphering t-SNE, UMAP, TriMap, and PaCMAP for data visualization," *JMLR*, 2021.

[23] C. Heine, H. Leitte, M. Hlawitschka, F. Iuricich, L. De Floriani, G. Scheuermann, H. Hagen, and C. Garth, "A survey of topology-based methods in visualization," *CGF*, 2016.

[24] T. Sousbie, "The persistent cosmic web and its filamentary structure: Theory and implementations," *Royal Astronomical Society*, 2011.

[25] N. Shivashankar, P. Pranav, V. Natarajan, R. van de Weygaert, E. P. Bos, and S. Rieder, "Felix: A topology based framework for visual exploration of cosmic filaments," *IEEE TVCG*, 2016.

[26] K. Anderson, J. Anderson, S. Palande, and B. Wang, "Topological data analysis of functional MRI connectivity in time and space domains," in *MICCAI Workshop on Connectomics in NeuroImaging*, 2018.

[27] H. A. Carr, J. Snoeyink, and M. van de Panne, "Simplifying Flexible Isosurfaces Using Local Geometric Measures," in *IEEE VIS*, 2004.

[28] A. Bock, H. Doraiswamy, A. Summers, and C. T. Silva, "Topoangler: Interactive topology-based extraction of fishes," *IEEE TVCG*, 2018.

[29] D. Guenther, R. Alvarez-Boto, J. Contreras-Garcia, J.-P. Piquemal, and J. Tierny, "Characterizing molecular interactions in chemical systems," *IEEE TVCG*, 2014.

[30] H. Bhatia, A. G. Gyulassy, V. Lordi, J. E. Pask, V. Pascucci, and P.-T. Bremer, "Topoms: Comprehensive topological exploration for molecular and condensed-matter systems," *J. Comput. Chem.*, 2018.

[31] M. Olejniczak and J. Tierny, "Topological data analysis of vortices in the magnetically-induced current density in lih molecule," *PCCP*, 2023.

[32] J. Kasten, J. Reininghaus, I. Hotz, and H. Hege, "Two-dimensional time-dependent vortex regions based on the acceleration magnitude," *IEEE TVCG*, 2011.

[33] F. Nauleau, F. Vivodtzev, T. Bridel-Bertomeu, H. Beaugendre, and J. Tierny, "Topological analysis of ensembles of hydrodynamic turbulent flows an experimental study," in *IEEE LDAV*, 2022.

[34] A. Gyulassy, M. A. Duchaineau, V. Natarajan, V. Pascucci, E. Bringa, A. Higginbotham, and B. Hamann, "Topologically clean distance fields," *IEEE TVCG*, 2007.

[35] A. Gyulassy, A. Knoll, K. Lau, B. Wang, P. Bremer, M. Papka, L. A. Curtiss, and V. Pascucci, "Interstitial and interlayer ion diffusion geometry extraction in graphitic nanosphere battery materials," *IEEE TVCG*, 2016.

[36] M. Soler, M. Petitfrere, G. Darche, M. Plainchault, B. Conche, and J. Tierny, "Ranking viscous finger simulations to an acquired ground truth with topology-aware matchings," in *IEEE LDAV*, 2019.

[37] A. Gyulassy, P. Bremer, R. Grout, H. Kolla, J. Chen, and V. Pascucci, "Stability of dissipation elements: A case study in combustion," *CGF*, 2014.

[38] D. E. Laney, P. Bremer, A. Mascarenhas, P. Miller, and V. Pascucci, "Understanding the structure of the turbulent mixing layer in hydrodynamic instabilities," *IEEE TVCG*, 2006.

[39] H. Edelsbunner and J. L. Harer, *Computational topology: an introduction*. AMS, 2010.

[40] A. Zomorodian, "Computational topology," in *Algorithms and theory of computation handbook: special topics and techniques*, 2nd ed., 2010.

[41] F. Chazal and B. Michel, "An introduction to topological data analysis: fundamental and practical aspects for data scientists," *Frontiers in artificial intelligence*, 2021.

[42] G. Weber, P.-T. Bremer, and V. Pascucci, "Topological Landscapes: A Terrain Metaphor for Scientific Data," *IEEE TVCG*, 2007.

[43] S. Gerber, P. Bremer, V. Pascucci, and R. Whitaker, "Visual Exploration of High Dimensional Scalar Functions," *IEEE TVCG*, 2010.

[44] S. Gerber, O. Rübel, P.-T. Bremer, V. Pascucci, and R. T. Whitaker, "Morse–Smale Regression," *J. Comput. Graph. Stat.*, 2013.

[45] P. Oesterling, C. Heine, H. Jänicke, and G. Scheuermann, "Visual analysis of high dimensional point clouds using topological landscapes," in *Proc. PacificVis*, 2010.

[46] P. Oesterling, G. Scheuermann, S. Teresniak, G. Heyer, S. Koch, T. Ertl, and G. H. Weber, "Two-stage framework for a topology-based projection and visualization of classified document collections," in *Proc. IEEE VAST*, 2010.

[47] P. Oesterling, C. Heine, H. Jänicke, G. Scheuermann, and G. Heyer, "Visualization of High Dimensional Point Clouds Using their Density Distribution's Topology," *IEEE TVCG*, 2011.

[48] P. Oesterling, C. Heine, G. H. Weber, and G. Scheuermann, "Visualization and point clouds as topological landscape profiles to guide local data analysis," *IEEE TVCG*, 2013.

[49] L. Yan, Y. Zhao, P. Rosen, C. Scheidegger, and B. Wang, "Homology-preserving dimensionality reduction via manifold landmarking and tearing," *CoRR*, 2018.

[50] G. Singh, F. Mémoli, and G. E. Carlsson, "Topological methods for the analysis of high dimensional data sets and 3D object recognition," in *PBG*. Eurographics Association, 2007.

[51] G. Reeb, "Sur les points singuliers d'une forme de Pfaff complètement intégrable ou d'une fonction numérique," *Comptes-rendus de l'Académie des Sciences*, 1946.

[52] V. Guardieiro, F. I. de Oliveira, H. Doraiswamy, L. G. Nonato, and C. Silva, "TopoMap++: A faster and more space efficient technique to compute projections with topological guarantees," *IEEE TVCG*, 2024.

[53] S. Barannikov, I. Trofimov, N. Balabin, and E. Burnaev, "Representation topology divergence: A method for comparing neural network representations," *arXiv preprint arXiv:2201.00058*, 2021.

[54] B. J. Nelson and Y. Luo, "Topology-preserving dimensionality reduction via interleaving optimization," *arXiv preprint arXiv:2201.13012*, 2022.

[55] S. T. Schönenberger, A. Varava, V. Polianskii, J. J. Chung, D. Kragic, and R. Siegwart, "Witness Autoencoder: Shaping the latent space with witness complexes," in *TDA & Beyond*, 2020.

[56] G. T. Toussaint, "The relative neighbourhood graph of a finite planar set," *Pattern recognition*, 1980.

[57] K. J. Supowit, "The relative neighborhood graph, with an application to minimum spanning trees," *JACM*, 1983.
[58] J. W. Jaromczyk and M. Kowaluk, "A note on relative neighborhood graphs," in *Proceedings of the third annual symposium on Computational geometry*, 1987.

[59] A. Lingas, "A linear-time construction of the relative neighborhood graph from the Delaunay triangulation," *Computational Geometry*, 1994.

[60] R. Urquhart, "Algorithms for computation of relative neighbourhood graph," *Electronics Letters*, 1980.

[61] H. Edelsbrunner and T. S. Tan, "A quadratic time algorithm for the minmax length triangulation," *SIAM Journal on Computing*, 1993.

[62] S. Barannikov, "Framed Morse complexes and its invariants." *Adv. Soviet Math.*, 1994.

[63] P. Frosini and C. Landi, "Size theory as a topological tool for computer vision," *Pattern Recognition and Image Analysis*, 1999.

[64] V. Robins, "Toward computing homology from finite approximations," *Topology Proceedings*, 1999.

[65] H. Edelsbrunner, D. Letscher, and A. Zomorodian, "Topological persistence and simplification," *DCG*, 2002.

[66] D. Cohen-Steiner, H. Edelsbrunner, and D. Morozov, "Vines and vineyards by updating persistence in linear time," in *Proceedings of the twenty-second annual symposium on Computational geometry*, 2006.

[67] P. Skraba and K. Turner, "Wasserstein stability for persistence diagrams," *arXiv preprint arXiv:2006.16824*, 2023.

[68] D. P. Bertsekas, "A new algorithm for the assignment problem," *Mathematical Programming*, 1981.

[69] M. Kerber, D. Morozov, and A. Nigmetov, "Geometry helps to compare persistence diagrams," 2017.

[70] F. Iuricich, "Persistence cycles for visual exploration of persistent homology," *IEEE TVCG*, 2022.

[71] P. Guillou, J. Vidal, and J. Tierny, "Discrete morse sandwich: Fast computation of persistence diagrams for scalar data–an algorithm and a benchmark," *IEEE TVCG*, 2023.

[72] U. Bauer, "Rips: efficient computation of Vietoris–Rips persistence barcodes," *Journal of Applied and Computational Topology*, 2021.

[73] C. Maria, J.-D. Boissonnat, M. Glisse, and M. Yvinec, "The Gudhi library: Simplicial complexes and persistent homology," in *ICMS*, 2014.

[74] U. Bauer, M. Kerber, J. Reininghaus, and H. Wagner, "PHAT – Persistent homology algorithms toolbox," *Journal of symbolic computation*, 2017.

[75] P. Skraba, G. Thoppe, and D. Yogeshwaran, "Randomly weighted $d$-complexes: Minimal spanning acycles and persistence diagrams," *arXiv preprint arXiv:1701.00239*, 2017.

[76] M. A. Koyama, V. Robins, K. Turner, and F. Memoli, "Reduced Vietoris–Rips complexes: New methods to compute Vietoris–Rips persistent homology," *arXiv preprint arXiv:2307.16333*, 2023.

[77] A. Fabri and S. Pion, "CGAL: the Computational Geometry Algorithms Library," in *Proceedings of the 17th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems*, 2009.

[78] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Z. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, "PyTorch: An Imperative Style, High-Performance Deep Learning Library," in *NeurIPS*, 2019.

[79] J. Tierny, G. Favelier, J. A. Levine, C. Gueunet, and M. Michaux, "The Topology ToolKit," *IEEE TVCG*, 2017.

[80] T. B. Masood, J. Budin, M. Falk, G. Favelier, C. Garth, C. Gueunet, P. Guillou, L. Hofmann, P. Hristov, A. Kamakshidasan et al., "An overview of the Topology ToolKit," *Topological Methods in Data Analysis and Visualization VI: Theory, Applications, and Software*, 2021.

[81] D. P. Kingma and J. Ba, "Adam: A method for stochastic optimization," *arXiv preprint arXiv:1412.6980*, 2017.

[82] S. A. Nene, S. K. Nayar, and H. Murase, "Columbia Object Image Library (COIL-20)," 1996.

[83] B. Wang, B. Summa, V. Pascucci, and M. Vejdemo-Johansson, "Branching and circular features in high dimensional data," *IEEE TVCG*, 2011.

[84] "Carnegie Mellon University - CMU Graphics Lab - Motion capture library." [Online]. Available: http://mocap.cs.cmu.edu/

[85] W. Saelens, R. Cannoodt, H. Todorov, and Y. Saeys, "A comparison of single-cell trajectory inference methods," *Nature biotechnology*, 2019.

[86] R. Cannoodt, W. Saelens, H. Todorov, and Y. Saeys, "Single-cell omics datasets containing a trajectory," 2018. [Online]. Available: https://zenodo.org/records/1443566

[87] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg et al., "Scikit-learn: Machine learning in python," *JMLR*, 2011.

[88] J. A. Lee and M. Verleysen, "Quality assessment of dimensionality reduction: Rank-based criteria," *Neurocomputing*, 2009.

[89] J. Venna and S. Kaski, "Local multidimensional scaling," *Neural Networks*, 2006.
# APPENDIX A
## ADDITIONAL BACKGROUND

### A. Homology

Homology groups, introduced by Henri Poincaré in the early 20th century in his Analysis Situs, are algebraic constructions that describe "holes" in a topological space. Here, we provide a brief summary of the notions used to define them.

Let $\mathcal{K}$ be a simplicial complex. A $k$-chain $c$ is a formal sum of $k$-simplices with modulo 2 coefficients<fnref n="3" />: $c = \sum \alpha_i \sigma_i$ with $\alpha_i \in \mathbb{Z}/2\mathbb{Z}$. Two $k$-chains can be summed simplex by simplex, which permits to define the group of $k$-chains, noted $C_k(\mathcal{K})$. The boundary of a $k$-simplex is the sum of its $(k-1)$-dimensional faces, i.e. more formally, if $\sigma = (v_0, \dots, v_k)$, then its boundary is

$$ \partial_k \sigma = \sum_{i=0}^{k} (v_0, \dots, \hat{v}_i, \dots, v_k) $$

where the hat means that $v_i$ is omitted. Then the boundary of a $k$-chain $c$ is defined as the modulo 2 sum of the boundaries of its simplices, i.e. $\partial_k c = \sum \alpha_i \partial_k \sigma_i$. This defines the $k$-th boundary map $\partial_k$, which is a morphism $\partial_k : C_k(\mathcal{K}) \to C_{k-1}(\mathcal{K})$.

A $k$-cycle is a $k$-chain whose boundary is 0. The group of $k$-cycles is noted $Z_k(\mathcal{K}) = \ker \partial_k$. A $k$-boundary is a $k$-chain which is the boundary of some $(k+1)$-chain. The group of $k$-boundaries is noted $B_k(\mathcal{K}) = \text{Im } \partial_{k+1}$. The fundamental lemma of homology states that $\partial_{k-1} \circ \partial_k = 0$, i.e., that the boundary of a boundary is always an empty chain. In particular this implies the inclusions $B_k(\mathcal{K}) \subseteq Z_k(\mathcal{K}) \subseteq C_k(\mathcal{K})$.

Two $k$-cycles $a, b \in Z_k(\mathcal{K})$ are homologous whenever $a = b + \partial c$ for some $c \in C_{k+1}(\mathcal{K})$. The set of the equivalence classes for this relation forms the $k$-th homology group of $\mathcal{K}$: $\mathcal{H}_k(\mathcal{K}) = Z_k(\mathcal{K})/B_k(\mathcal{K})$. Its rank, i.e. the maximal number of linearly independent classes (called *generators*), is called the $k$-th *Betti number* of $\mathcal{K}$: $\beta_k(\mathcal{K}) = \text{rank} \mathcal{H}_k(\mathcal{K}) = \log_2(|\mathcal{H}_k(\mathcal{K})|)$. The Betti numbers can be easily interpreted geometrically for small values of $k$: for a simplicial complex $\mathcal{K}$ embedded in $\mathbb{R}^3$, $\beta_0(\mathcal{K})$ is the number of connected components, $\beta_1(\mathcal{K})$ is the number of cycles, and $\beta_2(\mathcal{K})$ is the number of voids.

### B. Minmax length triangulations

In this section, we give more details on the computation of minmax length (MML) triangulations, that are all taken from [61]. General position is supposed, i.e., no two edges are equally long. A 2-edge is an edge whose right and left lenses are both non-empty. Conversely, a 1-edge is an edge such that one of its half-lens is empty but not the other one. By definition, 1- and 2-edges are not RNG edges. Note that every triangulation of an RNG polygon has a 2-edge (its longest edge is always a 2-edge).

An *expandable* 2-edge is a 2-edge $e = pq$ for which there exists $x \in \text{LeftLens}(e)$ and $y \in \text{RightLens}(e)$ such that the edges $px, xq, qy, yp$ are either RNG edges, or 1-edges with

an empty right lens (see Fig. 12). Expandability is a property that somehow allows an edge to be the longest edge of a triangulation. Indeed, let $e$ be an expandable 2-edge inside an RNG-polygon $\Pi$. Then it is possible to construct by recurrence a triangulation of $\Pi$ such that the only 2-edge – and thus longest edge – in that triangulation is $e$ itself.

Besides, the 2-edge lemma (Sec. 5.3 of [61]) states that there exists a MML triangulation of $\Pi$ that contains an expandable 2-edge. Therefore, to find a MML triangulation of an RNG-polygon, it suffices to find the shortest expandable 2-edge $e$, which is the longest edge of any MML-triangulation. This is the only information we use in our application. However, if a complete MML triangulation is wanted, $e = pq$ has to be found together with $x$ and $y$ defined as above. Then, the incomplete polygons defined by $px, xq, qy, yp$ have to be triangulated (Sec. 5.1 of [61]).

# APPENDIX B
## RAW DATA FOR FIG. 8

The point clouds $X \subset \mathbb{R}^3$ and $Z \subset \mathbb{R}^2$ used in the counter-example depicted in Fig. 8 are given by the following coordinates:

$$ X = \begin{bmatrix} 1 & 0 & 0 \\ 0.5 & 0.866 & 0 \\ -0.48 & 0.667 & 0 \\ -0.73 & -0.199 & 0.433 \\ -0.48 & -1.065 & 0 \\ 0.5 & -0.866 & 0 \end{bmatrix} \quad Z = \begin{bmatrix} 1 & 0 \\ 0.5 & 0.866 \\ -0.48 & 0.667 \\ -0.98 & -0.199 \\ -0.48 & -1.065 \\ 0.5 & -0.866 \end{bmatrix}. $$

# APPENDIX C
## ADDITIONAL QUANTITATIVE ANALYSIS

See Fig. 19 for topological metrics and some other commonly used DR indicators computed on the examples presented in the paper.

### A. Description of the indicators

a) Pairwise distances-based indicators: In addition to the metric distortion $\mathcal{D}$, one can measure the linear correlation $\text{LC} \in [-1, 1]$ between the pairwise distances in $X$ and in $Z$. Having $\text{LC} = 1$ means that these pairwise distances are perfectly correlated in high and low dimension.

b) Triplets-based indicator: The triplet accuracy $\text{TA} \in [0, 1]$ is the proportion of triangles whose three edges have the same relative order (in terms of length) both in high and low dimension [22]. It measures to some extent the preservation of the global structure of $X$, which is preserved when $\text{TA} = 1$.

c) Rank-based indicators: The rank of a point $X_i$ relative to another point $X_j$ is the integer $\rho_{ij} \in \mathbb{N}$ such that $X_j$ is the $\rho_{ij}$-th nearest neighbor of $X_i$. There exists indicators based on these ranks [88], [89] that measure the preservation of the nearest neighbors. In particular, the *trustworthiness* ($\text{Trust} \in [0, 1]$) is penalized when a high rank $\rho_{ij}$ in the input becomes low in the representation, i.e., when faraway points in high-dimension become neighbors in low-dimension. On the contrary, the *continuity* ($\text{Cont} \in [0, 1]$) is penalized when a low rank $\rho_{ij}$ in the input becomes high in the representation, i.e.,

<footnote n="3">³The coefficients can be chosen in another field or ring (which yields different homology groups), but $\mathbb{Z}/2\mathbb{Z}$ is the most common choice.</footnote>
<table><caption>Fig. 19. PH-based metrics and other common indicators used in DR (see Sec. C for a description). The best value (before rounding) for an indicator is written in bold, the second best (before rounding) is underlined.</caption><thead><tr><th rowspan="2">Dataset</th><th rowspan="2">Indicator</th><th colspan="2">Global methods</th><th colspan="3">Locally topology-aware methods</th><th colspan="4">Globally topology-aware methods</th></tr><tr><th>PCA [2]</th><th>MDS [3]</th><th>Isomap [4]</th><th>t-SNE [5]</th><th>UMAP [6]</th><th>TopoMap [7]</th><th>TopoAE [1]</th><th>TopoAE+<sup>W</sup> [8]</th><th>TopoAE++</th></tr></thead><tbody><tr><td rowspan="7">3Clusters</td><td>PDW<sup>0</sup></td><td>2.0e+01</td><td>1.5e+01</td><td>2.6e+01</td><td>6.7e+02</td><td>2.1e+02</td><td><strong>4.2e-04</strong></td><td>2.0e+01</td><td>1.2e+01</td><td><u>1.0e+01</u></td></tr><tr><td>PDW<sup>1</sup></td><td>4.5e-01</td><td>5.0e-01</td><td>4.7e-01</td><td>2.3e+01</td><td>7.2e-01</td><td>4.3e+00</td><td>5.8e-01</td><td><strong>1.1e-01</strong></td><td><u>4.4e-01</u></td></tr><tr><td>DD</td><td><strong>2.7e-01</strong></td><td><u>4.1e-01</u></td><td>1.1e+00</td><td>3.5e+01</td><td>8.5e+00</td><td>7.0e+00</td><td>1.3e+00</td><td>1.5e+00</td><td><u>6.4e-01</u></td></tr><tr><td>LC</td><td><strong>1.00</strong></td><td><u>0.99</u></td><td>0.98</td><td>0.85</td><td>0.89</td><td>0.86</td><td>0.95</td><td>0.93</td><td>0.97</td></tr><tr><td>TA</td><td><strong>0.89</strong></td><td><u>0.88</u></td><td>0.72</td><td>0.38</td><td>0.42</td><td>0.48</td><td>0.70</td><td>0.66</td><td>0.78</td></tr><tr><td>Trust</td><td>0.93</td><td>0.94</td><td>0.92</td><td><strong>0.99</strong></td><td>0.98</td><td>0.96</td><td>0.97</td><td>0.97</td><td>0.96</td></tr><tr><td>Cont</td><td><strong>0.99</strong></td><td>0.98</td><td>0.98</td><td>0.98</td><td>0.98</td><td>0.92</td><td>0.98</td><td>0.97</td><td>0.98</td></tr><tr><td rowspan="7">Twist</td><td>PDW<sup>0</sup></td><td>1.1e-01</td><td>1.3e-01</td><td>4.5e-01</td><td>2.9e+00</td><td>6.6e-01</td><td><strong>7.5e-06</strong></td><td>8.6e-02</td><td>2.8e-02</td><td><u>5.2e-03</u></td></tr><tr><td>PDW<sup>1</sup></td><td>1.1e+00</td><td>2.2e-01</td><td>1.9e+01</td><td>3.3e+00</td><td>1.5e+01</td><td>1.8e-01</td><td>1.6e-01</td><td><strong>1.8e-05</strong></td><td><u>1.7e-03</u></td></tr><tr><td>DD</td><td><u>2.2e-01</u></td><td><strong>2.1e-01</strong></td><td>2.4e+00</td><td>3.9e+00</td><td>2.8e+00</td><td>8.2e-01</td><td>3.8e-01</td><td>7.3e-01</td><td><u>7.8e-01</u></td></tr><tr><td>LC</td><td><strong>0.99</strong></td><td><u>0.99</u></td><td>0.86</td><td>0.99</td><td>0.58</td><td>0.86</td><td>0.98</td><td>0.99</td><td>0.99</td></tr><tr><td>TA</td><td><strong>0.91</strong></td><td><u>0.89</u></td><td>0.63</td><td>0.87</td><td>0.42</td><td>0.64</td><td>0.85</td><td>0.87</td><td>0.88</td></tr><tr><td>Trust</td><td>0.98</td><td>0.99</td><td><strong>1.00</strong></td><td>0.99</td><td>1.00</td><td>0.99</td><td>0.99</td><td>0.96</td><td><u>1.00</u></td></tr><tr><td>Cont</td><td>0.99</td><td>0.99</td><td><strong>1.00</strong></td><td>0.99</td><td>1.00</td><td>0.98</td><td>0.99</td><td>0.98</td><td><u>1.00</u></td></tr><tr><td rowspan="7">K_4</td><td>PDW<sup>0</sup></td><td>1.2e-01</td><td>6.9e-02</td><td>1.1e-01</td><td>3.0e+01</td><td>2.0e+01</td><td><strong>2.3e-07</strong></td><td>1.9e-01</td><td>7.8e-02</td><td><u>5.8e-02</u></td></tr><tr><td>PDW<sup>1</sup></td><td>2.7e-01</td><td>2.8e-01</td><td>3.8e-01</td><td>4.7e+01</td><td>3.9e-01</td><td>1.4e-01</td><td>6.5e-02</td><td><strong>7.9e-04</strong></td><td><u>2.4e-02</u></td></tr><tr><td>DD</td><td><u>2.7e-01</u></td><td><strong>1.8e-01</strong></td><td>5.1e-01</td><td>1.6e+01</td><td>9.1e+00</td><td>5.7e-01</td><td>1.0e+00</td><td>6.8e-01</td><td><u>3.4e-01</u></td></tr><tr><td>LC</td><td><u>0.83</u></td><td><strong>0.89</strong></td><td>0.65</td><td>0.66</td><td>0.53</td><td>0.57</td><td>0.75</td><td>0.68</td><td>0.77</td></tr><tr><td>TA</td><td><u>0.61</u></td><td><strong>0.63</strong></td><td>0.39</td><td>0.48</td><td>0.34</td><td>0.38</td><td>0.51</td><td>0.47</td><td>0.56</td></tr><tr><td>Trust</td><td>0.97</td><td>0.98</td><td>0.95</td><td><strong>1.00</strong></td><td>1.00</td><td>0.99</td><td>0.98</td><td>0.98</td><td><u>1.00</u></td></tr><tr><td>Cont</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.99</td><td>0.99</td><td>0.98</td><td>0.99</td><td>0.99</td><td><u>1.00</u></td></tr><tr><td rowspan="7">K_5</td><td>PDW<sup>0</sup></td><td>2.7e-01</td><td>1.5e-01</td><td>1.8e-01</td><td>9.3e+01</td><td>1.8e+01</td><td><strong>7.3e-07</strong></td><td>2.2e-01</td><td>1.3e-01</td><td><u>9.0e-02</u></td></tr><tr><td>PDW<sup>1</sup></td><td>2.8e-01</td><td>3.5e-01</td><td>1.3e-01</td><td>4.7e+01</td><td>9.7e-01</td><td>4.8e-01</td><td>1.6e-01</td><td><strong>2.2e-03</strong></td><td><u>8.0e-02</u></td></tr><tr><td>DD</td><td><u>3.1e-01</u></td><td><strong>2.4e-01</strong></td><td>4.7e-01</td><td>2.1e+01</td><td>9.8e+00</td><td>1.2e+00</td><td>9.2e-01</td><td>8.6e-01</td><td><u>3.4e-01</u></td></tr><tr><td>LC</td><td><u>0.81</u></td><td><strong>0.86</strong></td><td>0.67</td><td>0.76</td><td>0.58</td><td>0.41</td><td>0.73</td><td>0.72</td><td>0.81</td></tr><tr><td>TA</td><td><strong>0.63</strong></td><td><u>0.61</u></td><td>0.43</td><td>0.55</td><td>0.41</td><td>0.33</td><td>0.50</td><td>0.50</td><td>0.57</td></tr><tr><td>Trust</td><td>0.93</td><td>0.97</td><td>0.96</td><td><strong>1.00</strong></td><td>1.00</td><td>0.99</td><td>0.97</td><td>0.97</td><td>1.00</td></tr><tr><td>Cont</td><td>0.99</td><td>0.99</td><td><strong>1.00</strong></td><td>0.99</td><td>0.98</td><td>0.99</td><td>0.99</td><td>0.99</td><td><u>1.00</u></td></tr><tr><td rowspan="7">COIL-20-1</td><td>PDW<sup>0</sup></td><td>1.2e+02</td><td>8.1e+01</td><td>3.0e+00</td><td>1.5e+02</td><td>1.5e+02</td><td><strong>2.3e-06</strong></td><td>1.9e+00</td><td>4.5e+00</td><td><u>1.2e+00</u></td></tr><tr><td>PDW<sup>1</sup></td><td>1.2e+01</td><td>1.0e+01</td><td>6.0e+02</td><td>1.3e+01</td><td>1.2e+01</td><td>2.7e+01</td><td>7.5e+00</td><td><u>9.1e-02</u></td><td><strong>1.6e-02</strong></td></tr><tr><td>DD</td><td><u>2.9e+00</u></td><td><strong>2.0e+00</strong></td><td>2.1e+01</td><td>4.9e+00</td><td>6.2e+00</td><td>7.2e+00</td><td>8.2e+00</td><td>6.7e+00</td><td><u>1.9e+01</u></td></tr><tr><td>LC</td><td><strong>0.94</strong></td><td><u>0.93</u></td><td>0.83</td><td>0.88</td><td>0.85</td><td>0.74</td><td>0.81</td><td>0.62</td><td>0.88</td></tr><tr><td>TA</td><td><strong>0.81</strong></td><td>0.74</td><td>0.47</td><td><strong>0.77</strong></td><td>0.65</td><td>0.44</td><td>0.56</td><td>0.38</td><td>0.73</td></tr><tr><td>Trust</td><td>0.97</td><td>0.97</td><td><u>1.00</u></td><td><strong>1.00</strong></td><td>0.99</td><td>0.92</td><td>0.96</td><td>0.91</td><td>0.99</td></tr><tr><td>Cont</td><td>0.96</td><td>0.97</td><td>0.98</td><td><strong>1.00</strong></td><td>0.99</td><td>0.96</td><td>0.97</td><td>0.94</td><td>0.99</td></tr><tr><td rowspan="7">MoCap</td><td>PDW<sup>0</sup></td><td>4.3e+03</td><td>2.2e+03</td><td>3.9e+02</td><td>6.8e+03</td><td>6.8e+03</td><td><strong>2.4e-03</strong></td><td>1.6e+03</td><td>7.9e+02</td><td><u>1.9e+03</u></td></tr><tr><td>PDW<sup>1</sup></td><td>2.1e+03</td><td>1.9e+03</td><td>1.0e+04</td><td>1.6e+03</td><td>1.5e+03</td><td>2.1e+03</td><td>1.2e+03</td><td><strong>9.4e+00</strong></td><td><u>1.2e+01</u></td></tr><tr><td>DD</td><td><u>2.4e+01</u></td><td><strong>1.7e+01</strong></td><td>1.3e+02</td><td>9.9e+01</td><td>1.0e+02</td><td>8.5e+01</td><td>2.6e+01</td><td>5.6e+01</td><td><u>2.9e+01</u></td></tr><tr><td>LC</td><td><u>0.94</u></td><td><strong>0.95</strong></td><td>0.67</td><td>0.90</td><td>0.47</td><td>0.43</td><td>0.89</td><td>0.61</td><td>0.83</td></tr><tr><td>TA</td><td><u>0.72</u></td><td><strong>0.73</strong></td><td>0.53</td><td>0.64</td><td>0.34</td><td>0.35</td><td>0.68</td><td>0.42</td><td>0.59</td></tr><tr><td>Trust</td><td>0.98</td><td>0.98</td><td><u>1.00</u></td><td>1.00</td><td><strong>1.00</strong></td><td>0.97</td><td>0.98</td><td>0.93</td><td>1.00</td></tr><tr><td>Cont</td><td>0.98</td><td>0.99</td><td><strong>1.00</strong></td><td>1.00</td><td>0.99</td><td>0.98</td><td>0.99</td><td>0.98</td><td>1.00</td></tr><tr><td rowspan="7">SingleCell</td><td>PDW<sup>0</sup></td><td>4.6e+05</td><td>4.6e+05</td><td>3.5e+05</td><td>4.6e+05</td><td>4.6e+05</td><td><strong>1.2e-01</strong></td><td>4.6e+05</td><td><u>1.6e-05</u></td><td><u>3.8e+05</u></td></tr><tr><td>PDW<sup>1</sup></td><td>3.0e+03</td><td>3.1e+03</td><td>1.8e+05</td><td>1.8e+03</td><td>1.8e+03</td><td>2.8e+04</td><td>1.8e+03</td><td><strong>6.4e+02</strong></td><td><u>8.6e+02</u></td></tr><tr><td>DD</td><td><u>3.1e+01</u></td><td><strong>2.3e+01</strong></td><td>4.8e+02</td><td>9.6e+01</td><td>1.0e+02</td><td>1.2e+03</td><td>9.8e+01</td><td>3.2e+02</td><td><u>1.5e+02</u></td></tr><tr><td>LC</td><td><strong>0.99</strong></td><td><u>0.98</u></td><td>0.98</td><td>0.88</td><td>0.81</td><td>0.62</td><td>0.72</td><td>0.74</td><td>0.85</td></tr><tr><td>TA</td><td><strong>0.94</strong></td><td><u>0.86</u></td><td>0.84</td><td>0.68</td><td>0.54</td><td>0.40</td><td>0.55</td><td>0.52</td><td>0.61</td></tr><tr><td>Trust</td><td>0.99</td><td>0.99</td><td>0.99</td><td><strong>0.99</strong></td><td>0.99</td><td>0.96</td><td>0.97</td><td>0.96</td><td>0.97</td></tr><tr><td>Cont</td><td>0.99</td><td>0.99</td><td>0.99</td><td><strong>0.99</strong></td><td>0.99</td><td>0.95</td><td>0.94</td><td>0.96</td><td>0.96</td></tr></tbody></table>

when neighbors in high-dimension are projected to faraway points in low-dimension. These indicators are computed for a number *K* of nearest neighbors (in our experiments, *K* = 10).

## B. Observations

TopoMap has the best results for PDW<sup>0</sup>, which is expected since it preserves by design exactly PH<sup>0</sup>. Although in theory this quantity should be exactly 0, the approximation performed by the auction algorithm when estimating the Wasserstein distance (to compute PDW<sup>0</sup>) makes it slightly positive. Our approach (TopoAE++) has competitive results for this metric as well, which was expected since, as TopoAE, it constrains the preservation of PH<sup>0</sup> in the sense of Lem. 1.

Global indicators, namely the metric distortion $\mathcal{D}$, the linear correlation LC between pairwise distances and the triplet accuracy TA, are best preserved by global methods, i.e., PCA

and MDS. Our approach (TopoAE++) presents competitive results for these indicators when comparing to other locally topology-aware methods, while it clearly outperforms globally topology-aware methods.

Finally, TopoAE++ also presents good results for neighborhood quality indicators (Trust and Cont). Indeed, in our datasets, the correct embedding of the input cycle(s) favors the preservation of the neighborhood. More precisely, a self-intersection in the projection of an input cycle penalizes the trustworthiness (since faraway vertices would be projected as neighbors near this intersection), while a broken input cycle penalizes the continuity (since two neighbors in the input generator would be embedded faraway in low-dimension).