Title: No title found

Abstract: We propose a novel approach for preserving topological structures of the input space in latent representations of autoencoders. Using persistent homology, a technique from topological data analysis, we calculate topological signatures of both the input and latent space to derive a topological loss term. Under weak theoretical assumptions, we can construct this loss in a differentiable manner, such that the encoding learns to retain multi-scale connectivity information. We show that our approach is theoretically well-founded and that it exhibits favourable latent representations on a synthetic manifold as well as on real-world image data sets, while preserving low reconstruction errors.

1. INTRODUCTION:

While topological features, in particular multi-scale features derived from persistent homology, have seen increasing usage in the machine learning community (Carrière et al., 2019; Guss & Salakhutdinov, 2018; Hofer et al., 2017; 2019a; b; Ramamurthy et al., 2019; Reininghaus et al., 2015; Rieck et al., 2019a; b) , using topology directly as a constraint for current deep learning methods remains a challenge. This is due to the inherently discrete nature of these computations, making backpropagation through the computation of topological signatures immensely difficult or only possible in certain special circumstances (Chen et al., 2019; Hofer et al., 2019b; Poulenard et al., 2018) .

In this work, we present a novel approach that permits obtaining gradients during the computation of topological signatures. This permits employing topological constraints while training deep neural networks as well as building topology-preserving autoencoders based on the following contributions:

1. We develop a novel topological loss term for autoencoders that helps harmonise the topology of the data space and the topology of the latent space. 2. We prove that our approach is stable on the level of mini-batches, resulting in suitable approximations of the persistent homology of a data set. 3. We empirically demonstrate that our novel loss term aids in dimensionality reduction by preserving topological structures in data sets; in particular, the learned latent representations are useful in that the preservation of topological structures can aid interpretability.

2 BACKGROUND: PERSISTENT HOMOLOGY Persistent homology (Barannikov, 1994; Edelsbrunner & Harer, 2008 ) is a method from the field of computational topology, which develops tools for analysing topological features (connectivitybased features such as connected components) of data sets. We first need to introduce the underlying concept of simplicial homology. For a simplicial complex K, i.e. a generalised graph with higher-order connectivity information such as cliques, simplicial homology employs matrix reduction algorithms to assign K a family of groups, the homology groups. The d th homology group H d (K) of K contains d-dimensional topological features, such as connected components (d = 0), cycles/tunnels (d = 1), and voids (d = 2). Homology groups are typically summarised by their ranks, thereby obtaining a simple invariant "signature" of a manifold. For example, a circle in R 2 has one feature with d = 1 (a cycle), and one feature with d = 0 (a connected component). In practice, the underlying manifold M is unknown and we are working with a point cloud X := {x 1 , . . . , x n } ⊆ R d and a metric dist : X × X → R such as the Euclidean distance. Persistent homology extends simplicial homology to this setting: instead of approximating M by means of a single simplicial complex, which would be an unstable procedure due to the discrete nature of X, persistent homology tracks changes in the homology groups over multiple scales of the metric. This is achieved by constructing a special simplicial complex, the Vietoris-Rips complex (Vietoris, 1927) . For 0 ≤ < ∞, the Vietoris-Rips complex of X at scale , denoted by R (X), contains all simplices (i.e. subsets) of X whose elements {x 0 , x 1 , . . . } satisfy dist(x i , x j ) ≤ for all i, j. Given a matrix Aof pairwise distances of a point cloud X, we will use R (A) and R (X) interchangeably because constructing R only requires distances. Vietoris-Rips complexes satisfy a nesting relation, i.e. R i (X) ⊆ R j (X) for i ≤ j , making it possible to track changes in the homology groups as increases (Edelsbrunner et al., 2002) . Figure 1 illustrates this process. Since X contains a finite number of points, a maximum ˜ value exists for which the connectivity stabilises; therefore, calculating R ˜ is sufficient to obtain topological features at all scales.

We write PH(R (X)) for the persistent homology calculation of the Vietoris-Rips complex. It results in a tuple

({D 1 , D 2 , . . .}, {π 1 , π 2 , . . .}) of persistence diagrams (1 st component) and persistence pairings (2 nd component). The d-dimensional persistence diagram D d (Figure 1d) of R (X) contains coordinates of the form (a, b)

, where a refers to a threshold at which a d-dimensional topological feature is created in the Vietoris-Rips complex, and b refers to a threshold at which it is destroyed. When d = 0, for example, the threshold indicates at which distance two connected components in X are merged into one; this is related to spanning trees (Kurlin, 2015) . The d-dimensional persistence pairing contains indices (i, j) corresponding to simplices s i , s j ∈ R (X) that create and destroy the topological feature identified by (a, b) ∈ D d , respectively. Persistence diagrams are known to be stable with respect to small perturbations in the data set (Cohen-Steiner et al., 2007) . Two diagrams D and D can be compared using the bottleneck distance

d b (D, D ) := inf η : D→D sup x∈ D x -η(x) ∞

, where η : D → D denotes a bijection between the points of the two diagrams, and • ∞ refers to the L ∞ norm. We use D X to refer to the set of persistence diagrams of a point cloud X arising from PH(R (X)).

3. PROPOSED METHOD:

We propose a generic framework for constraining autoencoders to preserve topological structures (measured via persistent homology) of the data space in their latent encodings.

3.1. VIETORIS-RIPS COMPLEX CALCULATION:

Given a finite metric space S, we first calculate the persistent homology of the Vietoris-Rips complex of its distance matrix A S . Typically, the usual Euclidean distance is used to calculate A S , but both the persistent homology calculation and our method are not restricted to any particular distance; previous research (Wagner & Dłotko, 2014) shows that even similarity measures that do not satisfy the properties of a metric can be used successfully with PH(•). Subsequently, let := max A S so that R A S is the corresponding Vietoris-Rips complex as described in Section 2. Given a maximum dimension 1 of d ∈ N >0 , we obtain a set of persistence diagrams D S , and a set of persistence pairings π S . The d th persistence pairing π S d contains indices of simplices that are relevant for creating and destroying d-dimensional topological features. We can consider each pairing to represent edge indices, namely the edges that are deemed to be "topologically relevant" by the computation of persistent homology (see below for more details). This works because the Vietoris-Rips complex is a clique complex, i.e. it is fully determined by its edges (Zomorodian, 2010).

For 0-dimensional topological features, it is sufficient to look at the edge indices, i.e. the indices of "destroyer" simplices, contained in π S 0 . Each of these indices corresponds to an edge in the minimum spanning tree of the data set. This calculation is computationally efficient, having a worst-case complexity of O m2 • α m 2 , where m is the batch size and α(•) denotes the extremely slow-growing inverse Ackermann function (Cormen et al., 2009, Chapter 22) . For 1-dimensional features, where edges are paired with triangles, we obtain edge indices by selecting the edge with the maximum weight of the triangle. While this procedure (and thus our method) generalises to higher dimensions, our current implementation supports no higher-dimensional features. Since preliminary experiments showed that including 1-dimensional topological features merely increases runtime, we will focus only on 0-dimensional persistence diagrams in the subsequent experiments. We thus use D S , π S to denote the 0-dimensional persistence diagram and pairing of S, respectively.

3.2. TOPOLOGICAL AUTOENCODER:

In the following, we consider a mini-batch X of size m from the data space X as a point cloud. Furthermore, we define an autoencoder as the composition of two functions h • g, where g : X → Z represents the encoder and h : Z → X represents the decoder. We denote latent codes with Z := g(X). During a forward pass of the autoencoder, we compute the persistent homology of the mini-batch in both the data as well as the latent space, yielding two sets of tuples, i.e. D X , π X := PH(R (X)) and D Z , π Z := PH(R (Z)). The values of the persistence diagram can be retrieved by subsetting the distance matrix with the indices provided by the persistence pairings; we write D X A X π X to indicate that the diagram, which is a set, contains the same information as the distances we retrieve with the pairing. We treat A X π X as a vector in R |π X | . Informally speaking, the persistent homology calculation can thus be seen as a selection of topologically relevant edges of the Vietoris-Rips complex, followed by the selection of corresponding entries in the distance matrix. By comparing both diagrams, we can construct a topological regularisation term L t := L t A X , A Z , π X , π Z , which we add to the reconstruction loss of an autoencoder, i.e.

L := L r (X, h(g(X))) + λ L t (1)

where λ ∈ R is a parameter to control the strength of the regularisation. Next, we discuss how to specify L t : in our case, the PH calculation represents a selection of topologically relevant distances from the distance matrix. Each persistence diagram entry corresponds to a distance between two data points. We assume that the distances are unique so that each entry in the diagram has an infinitesimal neighbourhood that only contains a single point. Given this fixed pairing and a differentiable distance function, the persistence diagram entries are thus also differentiable with respect to the encoder parameter. Hence, the persistence pairing does not change upon a small perturbation of the underlying distances, thereby guaranteeing the existence of the derivative of our loss function. This, in turn, permits the calculation of gradients for backpropagation.

A straightforward approach to impose the data space topology on the latent space would be to directly calculate a loss based on the selected distances in both spaces. Such an approach will not result in informative gradients for the autoencoder, as it merely compares topological features without matching 2 the edges between R (X) and R (Z). A cleaner approach would be to enforce similarity on the intersection of the selected edges in both complexes. However, this would initially include very few edges, preventing efficient training and leading to very biased estimates of the topological alignments between the spaces 3 . To overcome this, we account for the union of all selected edges in X and Z. Our topological loss term decomposes into two components, each handling the "directed" loss occurring as topological features in one of the two spaces remain fixed. Hence, L t = L X →Z + L Z→X , where

EQUATION

respectively. The key idea for both terms is to align and preserve topologically relevant distances from both spaces. By taking the union of all selected edges (and the corresponding distances), we obtain an informative loss term that is determined by at least |X| distances. This loss can be seen as a more generic version of the loss introduced by Hofer et al. (2019b) , whose formulation does not take the two directed components into account and optimises the destruction values of all persistence tuples with respect to a uniform parameter (see also Section 4 for a brief discussion). By contrast, our formulation aims to to align the distances between X and Z (which in turn will lead to an alignment of distances between X and Z). If they are aligned perfectly, L X →Z = L Z→X = 0 because both pairings and their corresponding distances coincide. The converse implication is not true: if L t = 0, the persistence pairings (and their corresponding persistence diagrams) are not necessarily identical.

Letting θ refer to the parameters of the encoder, we have

EQUATION

where π X denotes the cardinality of a persistence pairing and A Z π X i refers to the i th entry of the vector of paired distances. This derivation works analogously for L Z→X (with π X being replaced by π Z ). Furthermore, any derivative of A X with respect to θ must vanish because the distances of the input samples do not depend on the encoding by definition.

These equations presume infinitesimal perturbations. In fact, the persistence diagrams change in a non-differentiable manner during the training phase but for a given update step, a diagram is robust to infinitesimal changes of its entries (Cohen-Steiner et al., 2007) . As a consequence, our topological loss is differentiable for each update step during training. We make our code publicly available4 .

3.3. STABILITY:

With persistence diagrams being stable under small perturbations of the underlying space (Cohen-Steiner et al., 2007) , we still have to analyse our topological approximation on the level of minibatches. The following theorem guarantees that subsampled persistence diagrams are close to the persistence diagrams of the original point cloud.

Theorem 1. Let X be a point cloud of cardinality n and X (m) be one subsample of X of cardinality m, i.e. X (m) ⊆ X, sampled without replacement. We can bound the probability of the persistence diagrams of X (m) exceeding a threshold in terms of the bottleneck distance as

EQUATION

)

where d H refers to the Hausdorff distance between the point cloud and its subsample.

Proof. See Section A.1 in the supplementary materials.

For m → n, we have lim m→n d H X, X (m) = 0. Please refer to Section A.2 for an analysis of empirical convergence rates as well as a discussion of a worst-case bound. Given certain independence assumptions, the next theorem approximates the expected value of the Hausdorff distance. The calculation of an exact representation is beyond the scope of this work, though.

Theorem 2. Let A ∈ R n×m be the distance matrix between samples of X and X (m) , where the rows are sorted such that the first m rows correspond to the columns of the m subsampled points with diagonal elements a ii = 0. Assume that the entries a ij with i > m are random samples following a distance distribution F D with supp(f D ) ∈ R ≥0 . The minimal distances δ i for rows with i > m follow a distribution F ∆ . Letting Z := max 1≤i≤n δ i with a corresponding distribution F Z , the expected Hausdorff distance between X and X (m) for m < n is bounded by:

E d H (X, X (m) ) = E Z∼F Z [Z] ≤ +∞ 0 1 -F D (z) (n-1) dz ≤ +∞ 0 1 -F D (z) m(n-m) dz (6)

Proof. See Section A.3 in the supplementary materials.

From Eq. 6, we obtain E[d H (X, X m )] = 0 as m → n, so the expected value converges as the subsample size approaches the total sample size5 . We conclude that our subsampling approach results in point clouds that are suitable proxies for the large-scale topological structures of the point cloud X.

4. RELATED WORK:

Computational topology and persistent homology (PH) have started gaining traction in several areas of machine learning research. PH is often used as as post hoc method for analysing topological characteristics of data sets. Thus, there are several methods that compare topological features of high-dimensional spaces with different embeddings to assess the fidelity and quality of a specific embedding scheme (Khrulkov & Oseledets, 2018; Paul & Chalup, 2017; Rieck & Leitte, 2015; 2017; Yan et al., 2018) . PH can also be used to characterise the training of neural networks (Guss & Salakhutdinov, 2018; Rieck et al., 2019b) , as well as their decision boundaries (Ramamurthy et al., 2019) . Our method differs from all these publications in that we are able to obtain gradient information to update a model while training. Alternatively, topological features can be integrated into classifiers to improve classification performance. Hofer et al. (2017) propose a neural network layer that learns projections of persistence diagrams, which can subsequently be used as feature descriptors to classify structured data. Moreover, several vectorisation strategies for persistence diagrams exist (Adams et al., 2017; Carrière et al., 2015) , making it possible to use them in kernelbased classifiers. These strategies have been subsumed (Carrière et al., 2019) in a novel architecture based on deep sets. The commonality of these approaches is that they treat persistence diagrams as being fixed; while they are capable of learning suitable parameters for classifying them, they cannot adjust input data to better approximate a certain topology. Such topology-based adjustments have only recently become feasible. Poulenard et al. (2018) demonstrated how to optimise real-valued functions based on their topology. This constitutes the first approach for aligning persistence diagrams by modifying input data; it requires the connectivity of the data to be known, and the optimised functions have to be node-based and scalar-valued. By contrast, our method works directly on distances and sidesteps connectivity calculations via the Vietoris-Rips complex. Chen et al. ( 2019) use a similar optimisation technique to regularise the decision boundary of a classifier. However, this requires discretising the space, which can be computationally expensive. Hofer et al. (2019b) , the closest work to ours, also presents a differentiable loss term. They propose regularising the connectivity of the latent space by itself. Specifically, their formulation enforces a single scale on the latent space. The learned encoding is then applied to a classification task. By contrast, our directed loss term aims to preserve the data space topology in the latent space for dimensionality reduction.

5. EXPERIMENTS:

Our main task is to learn a latent space in an unsupervised manner such that topological features of the data space-measured using persistent homology approximations on every batch-are preserved as much as possible. 

5.1. EXPERIMENTAL SETUP:

Data sets We use a SPHERES data set that consists of ten high-dimensional 100-spheres living in 101-dimensional space that are enclosed by one larger sphere that consists of the same number of points as the total of inner spheres (please refer to Section A.4 for more details). We also use three image data sets (MNIST, FASHION-MNIST, and CIFAR-10), which are particularly amenable to our topology-based analysis because real-world images are known to lie on or near low-dimensional manifolds (Lee et al., 2003; Peyré, 2009) .

Baselines & Training:

We compare our approach with several dimensionality reduction techniques, including UMAP (McInnes et al., 2018) , t-SNE (van der Maaten & Hinton, 2008), Isomap (Tenenbaum et al., 2000) , PCA, as well as standard autoencoders (AE). We apply our proposed topological constraint to this standard autoencoder architecture (TopoAE). For comparability and interpretability, each method is restricted to two latent dimensions. We split each data set into training and testing (using the predefined split if available; 90% versus 10% otherwise). Additionally, we remove 15% of the training split as a validation data set for tuning the hyperparameters. We normalised our topological loss term by the batch size m in order to disentangle λ from it. All autoencoders employ batch-norm and are optimized using ADAM (Kingma & Ba, 2014). Since t-SNE is not intended to be applied applying to previously unseen test samples, we evaluate this method only on the train split.

In addition, significant computational scaling issues forced us to forgo running a hyperparameter search for Isomap on real-world data sets, so we only compare this algorithm on the synthetic data set. Please refer to Section A.5 for more details on architectures and hyperparameters.

Evaluation We evaluate the quality of latent representations in terms of (1) low-dimensional visualisations, (2) dimensionality reduction quality metrics (evaluated between input data and latent codes), and (3) reconstruction errors (Data MSE; evaluated between input and reconstructed data), provided that invertible transformations are available6 . For (2), we consider several measures (please refer to Section A.6 for more details). First, we calculate KL σ , the Kullback-Leibler divergence between the density estimates of the input and latent space, based on density estimates (Chazal et al., 2011; 2014b) , where σ ∈ R >0 denotes the length scale of the Gaussian kernel, which is varied to account for multiple data scales. We chose minimising KL 0.1 as our hyperparameter search objective. Furthermore, we calculate common non-linear dimensionality reduction (NLDR) quality metrics, which use the pairwise distance matrices of the input and the latent space (as indicated by the " " in the abbreviations), namely the root mean square error ( -RMSE), which-despite its name-is not The reported measures are computed on the testing splits (except for t-SNE where no transformation between splits is available, so we report the measures on a random subsample of the train split).

5.2. RESULTS:



Quantitative Results:

Table 1 reports the quantitative results. Overall, we observe that our method shines in preserving the data density over multiple length scales (as measured by KL). Furthermore, we find that TopoAE displays competitive continuity values ( -Cont) and reconstruction errors (Data MSE). The remaining classical measures favour the baselines (foremost the train (!) performance of t-SNE). However, we will subsequently see that those classic measures fail at detecting the relevant structural information, as exemplified with known ground truth manifolds, such as the SPHERES data set.

Visualisation of latent spaces For the SPHERES data set (Figure 2 ), we observe that only our method is capable of assessing the nesting relationship of the high-dimensional spheres correctly. By contrast, t-SNE "cuts open" the enclosing sphere, distributing most of its points around the remaining spheres. We see that the KL-divergence confirms the visual assessment that only our proposed method preserves the relevant structure of this dataset. Several classical evaluation measures, however, favour t-SNE, even though this method fails to capture the global structure and nesting relationship of the enclosing sphere manifold accounting for half of the dataset. On FASHION-MNIST (Figure 3a -e), we see that, as opposed to AE, which is purely driven by the reconstruction error, our method has the additional objective of preserving structure. Here, this constraint helps the autoencoder to "organise" the latent space, resulting in a comparable pattern as in UMAP, which is also topologically motivated. Furthermore, we observe that t-SNE tends to fragment certain classes (dark orange, red) into multiple distinct subgroups. This is very likely not to reflect the underlying manifold structure, but constitutes an artefact frequently observed with this method. For MNIST, the latent embeddings (Figure 3f-j ) demonstrate that the non-linear competitors-mostly by pulling apart distinct classes-lose some of the relationship information between clusters when comparing against our proposed method or PCA. Finally, we observe that CIFAR-10 (Figure 3k -o), is challenging to embed in two latent dimensions in a purely unsupervised manner. Interestingly, our method (repeatedly) identified a linear substructure separating the latent space in two additional groups of classes.

6. DISCUSSION AND CONCLUSIONS:

We presented topological autoencoders, a novel method for preserving topological information, measured in terms of persistent homology, of the input space when learning latent representations with deep neural networks. Under weak theoretical assumptions, we showed how our persistent homology (PH) calculations can be combined with backpropagation; moreover, we proved that approximating PH on the level of mini-batches is theoretically justified.

In our experiments, we observed that our method is uniquely able to capture spatial relationships between nested high-dimensional spheres. This is relevant, as the ability to cope with several manifolds in the domain of manifold learning still remains a challenging task. On real-world data sets, we observed that our topological loss leads to competitive performance in terms of numerous quality metrics (such as a density preservation metric), while not adversely affecting the reconstruction error. In both synthetic and real-world data sets, we obtain interesting representations, as our method does not merely pull apart different classes, but tries to spatially arrange them meaningfully. Thus, we do not observe mere distinct "clouds", but rather entangled structures, which we consider to constitute a more meaningful representation of the underlying manifolds (an auxiliary analysis in Supplementary Section A.9 confirms that our method influences topological features, measured using PH, in a beneficial manner).

Our topological loss formulation is highly generic; it only requires the existence of a distance matrix between individual samples (either globally, or on the level of batches). As a consequence, our topological loss term can be directly integrated into a variety of different architectures and is not limited to standard autoencoders. For instance, we can also apply our constraint to variational setups (see Figure A.7 for a sketch).

Employing our generic constraint to more involved architectures will be an exciting route for future work. One issue with the calculation is that, given the computational complexity of calculating R (•), for higher-dimensional features we would scale progressively worse with increasing batch size. However, in our low-dimensional setup, we observed that runtime tends to grow with decreasing batch-size, i.e. the mini-batch speed-up still dominates runtime (for more details concerning the effect of batch sizes, see Supplementary Section A.7). In future work, scaling to higher dimensions could be mitigated by approximating the calculation of persistent homology (Choudhary et al., 2018; Kerber & Sharathkumar, 2013; Sheehy, 2013) or by exploiting recent advances in parallelising it (Bauer et al., 2014; Lewis & Morozov, 2015) . 

A APPENDIX:

A.1 PROOF OF THEOREM 1

Theorem 1. Let X be a point cloud of cardinality n and X (m) be one subsample of X of cardinality m, i.e. X (m) ⊆ X, sampled without replacement. We can bound the probability of X (m) exceeding a threshold in terms of the bottleneck distance as

EQUATION

where d H refers to the Hausdorff distance between the point cloud and its subsample, i.e. 

for a baseline distance dist(x, y) such as the Euclidean distance.

Proof. The stability of persistent homology calculations was proved by Chazal et al. (2014a) for finite metric spaces. More precisely, given two metric spaces X and Y , we have

EQUATION

where d GH (X, Y ) refers to the Gromov-Hausdorff distance (Burago et al., 2001, p. 254) of the two spaces. It is defined as the infimum Hausdorff distance over all isometric embeddings of X and Y . This distance can be employed for shape comparison (Chazal et al., 2009; Mémoli & Sapiro, 2004) , but is hard to compute. In our case, with X = X and Y = X (m) , we consider both spaces to have the same metric (for Y , we take the canonical restriction of the metric from X to the subspace Y ). By definition of the Gromov-Hausdorff distance, we thus have We can see that the convergence rate is roughly similar, but shown on different absolute levels that depend on the ambient dimension. While bounding the convergence rate of this expression is feasible (Chazal et al., 2015a; b) , it requires more involved assumptions on the measures from which X and X (m) are sampled. Additionally, we can give a simple bound using the diameter diam(X) := sup{dist(x, y) | x, y ∈ X}. We have d H X, X (m) ≤ diam(X) because the supremum is guaranteed to be an upper bound for the Hausdorff distance. This worst-case bound does not account for the sample size (or mini-batch size) m, though (see Theorem 2 for an expression that takes m into account).

EQUATION

A.3 PROOF OF THEOREM 2:

Prior to the proof we state two observations that arise from our special setting of dealing with finite point clouds. Observation 1. Since X (m) ⊆ X, we have sup x ∈X (m) inf x∈X dist(x, x ) = 0. Hence, the Hausdorff distance simplifies to:

EQUATION

In other words, we only have to consider a "one-sided" expression of the distance because the distance from the subsample to the original point cloud is always zero. Observation 2. Since our point clouds of interest are finite sets, the suprema and infima of the Hausdorff distance coincide with the maxima and minima, which we will subsequently use for easier readability.

Hence, the computation of d H (X, X (m ) can be divided into three steps.

1. Using the baseline distance dist(•, •), we compute a distance matrix A ∈ R n×m between all points in X and X (m) .

2. For each of the n points in X, we compute the minimal distance to the m samples of X (m) by extracting the minimal distance per row of A and gather all minimal distances in δ ∈ R n .

3. Finally, we return the maximal entry of δ as d H X, X (m) .

In the subsequent proof, we require an independence assumption of the samples.

Proof. Using Observations 1 and 2 we obtain a simplified expression for the Hausdorff distance, i.e.

EQUATION

The minimal distances of the first m rows of A are trivially 0. Hence, the outer maximum is determined by the remaining nm row minima {δ i | m < i ≤ n } with δ i = min 1≤j≤m (a ij ). Those minima follow the distribution F ∆ (y) with

EQUATION

)

Next, we consider Z := max 1≤i≤n δ i . To evaluate the density of Z, we first need to derive its distribution F Z :

F Z (z) = P(Z ≤ z) = P max m<i≤n δ i ≤ z = P   m<i≤n δ i ≤ z   (15)

Next, we approximate Z by Z by imposing i.i.d sampling of the minimal distances δ i from F ∆ . This is an approximation because in practice, the rows m + 1 to n are not stochastically independent because of the triangular inequality that holds for metrics. However, assuming i.i.d., we arrive at

EQUATION

. Since Z has positive support its expectation can then be evaluated as: The independence assumption leading to Z results in overestimating the variance of the drawn minima δ i . Thus, the expected maximum of those minima, E[Z ], is overestimating the actual expectation of the maximum E[Z], which is why Eq. 17 to Eq. 18 constitute an upper bound of E[Z], and equivalently, an upper bound of E d H (X, X (m) ) . When increasing m, E[d H (X, X m )] decreases monotonically since for a particular m, we draw nm samples from the minimal distance distribution F ∆ , and their maximum determines the Hausdorff distance. In contrast, our preliminary upper bound on the left-hand side of Eq. 18 forms a downwards-facing parabola due to the quadratic form in the exponent. This indicates that a tighter bound is achieved for m = n by using the minimal subsample size of m = 1.

EQUATION

A.4 SYNTHETIC DATASET SPHERES consists of eleven high-dimensional 100-spheres living in 101-dimensional space. Ten spheres of radius r = 5 are each shifted in a random direction (by adding the same Gaussian noise vector per sphere). To this end, we draw ten d-dimensional Gaussian vectors following N (0, I( 10 / √ d))

for d = 101. Crucially, to add interesting topological information to the data set, the ten spheres are enclosed by an additional larger sphere of radius 5r. The spheres were generated using the library scikit-tda.

A.5 ARCHITECTURES AND HYPERPARAMETER TUNING:

Architectures for synthetic data set For the synthetically generated dataset we use a simple Multilayer perceptron architecture consisting of two hidden layer with 32 neurons each both encoder and decoder and a bottleneck of two neurons such that the sequence of hidden-layer neurons is 32 -32 -2 -32 -32. ReLU non-linearities and batch normalization were applied between the layers excluding the output layer and the bottleneck layer. The networks were fit using mean squared error loss. Architectures for real world data sets For the MNIST, FASHION-MNIST, and CIFAR-10 datasets, we use an architecture inspired by DeepAE (Hinton & Salakhutdinov, 2006) . This architecture is composed of 3 layers of hidden neurons of decreasing size (1000 -500 -250) for the encoder part, a bottleneck of two neurons, and a sequence of three layers of hidden neurons in decreasing size (250 -500 -1000) for the decoder. In contrast to the originally proposed architecture, we applied ReLU non-linearities and batch normalization between the layers as we observed faster and more stable training. For the non-linearities of the final layer, we applied the tanh non-linearity, such that the image of the activation matches the range of input images scaled between -1 and 1. Also here, we applied mean squared error loss.

All neural network architectures were fit using Adam and weight-decay of 10 -5 .

Hyperparameter tuning For hyperparameter tuning we apply random sampling of hyperparameters using the scikit-optimize library (scikit-optimize contributers, 2018) with 20 calls per method on all datasets. We select the best model parameters in terms of KL 0.1 on the validation split and evaluate and report it on the test split. To estimate performance means and standard deviations, we repeated the evaluation on an independent test split 5 times by using the best parameters (as identified in the hyperparameter search on the validation split) and refitting the models by resampling the train-validation split.

Neural networks For the neural networks we sample the learning rate log-uniformly in the range 10 -4 -10 -2 , the batch-size uniformly between 16 and 128, and for the TopoAE method we sample the regularization strength log-uniformly in the range 10 -1 -3. Each model was allowed to train for at most 100 epochs, whereas we applied early stopping with patience = 10 based on the validation loss.

Competitor methods For t-SNE, we sample the perplexity uniformly in the range 5 -50 and the learning rate log-uniformly in the range 10 -1000. For Isomap and UMAP, the number of neighbors included in the computation was varied between 15 -500. For UMAP, we additionally vary the min_dist parameter uniformly between 0 and 1.

A.6 MEASURING THE QUALITY OF LATENT REPRESENTATIONS:

Next to the reconstruction error (if available; please see the paper for a discussion on this), we use a variety of NLDR metrics to assess the quality of our method. Our primary interest concerns the quality of the latent space because, among others, it can be used to visualise the data set. We initially considered classical quality metrics from non-linear dimensionality reduction (NLDR) algorithms (see Bibal & Frénay (2019) ; Gracia et al. (2014); van der Maaten et al. (2009) for more in-depth descriptions), namely

(1) the root mean square error ( -RMSE) between the distance matrix of the original space and the latent space (as mentioned in the main text, this is not related to the reconstruction error), (2) the mean relative rank error ( -MRRE), which measures the changes in ranks of distances in the original space and the latent space (Lee & Verleysen, 2009) , (3) the trustworthiness ( -Trust) measure (Venna & Kaski, 2006) , which checks to what extent the k nearest neighbours of a point are preserved when going from the original space to the latent space, and (4) the continuity ( -Cont) measure (Venna & Kaski, 2006) , which is defined analogously to -Trust, but checks to what extent neighbours are preserved when going from the latent space to the original space.

All of these measures are defined based on comparisons of the original space and the latent space; the reconstructed space is not used here. As an additional measure, we calculate the Kullback-Leibler divergence between density distributions of the input space and the latent space. Specifically, for a point cloud X with an associated distance dist, we first use the distance to a measure density estimator (Chazal et al., 2011; 2014b) , defined as f σ X (x) := y∈X exp -σ -1 dist(x, y) 2 , where σ ∈ R >0 represents a length scale parameter. For dist, we use the Euclidean distance and normalise it between 0 and 1. Given σ, we evaluate KL σ := KL f σ X f σ Z , which measures the similarity between the two density distributions. Ideally, we want the two distributions to be similar because this implies that density estimates in a low-dimensional representation are similar to the ones in the original space.

A.7 ASSESSING THE BATCH SIZE:

As we used fixed architectures for the hyperparameter search, the batch-size remains the main determinant for the runtime of TopoAE. In Figure A .6, we display trends (linear fits) on how loss measures vary with batch size. Addtionally, we draw runtime estimates. As we applied early stopping, for better comparability, we approximated the epoch-wise runtime by dividing the execution time of a run by its number of completed epochs. Interestingly, these plots suggest that the runtime grows with decreasing batch-size (even though the topological computation is more costly for larger batch-sizes!).

In these experiments, sticking to 0-dimensional topological features we conclude that the benefit of using mini-batches for neural network training still dominate the topological computations. The few steep peaks most likely represent outliers (the corresponding runs stopped after few epochs, which is why the effective runtime could be overestimated).

For the loss measures, we see that reconstruction loss tends to decrease with increasing batch-size, while our topological loss tends to increase with increasing batch-size (despite normalization). The second observation might be due to larger batch-size enabling more complex data point arrangements and corresponding topologies.

A.8 EXTENDING TO VARIATIONAL AUTOENCODERS:

In Figure A.7 we sketch a preliminary experiment, where we apply our topological constraint to variational autoencoders for the SPHERES dataset. Also here, we observe that our constraint helps identifying the nesting structure of the enclosing sphere.

A.9 TOPOLOGICAL DISTANCE CALCULATIONS:

To assess the topological fidelity of the resulting latent spaces, we calculate several topological distances between the test data set (full dimensionality) and the latent spaces obtained from each method (two dimensions). More precisely, we calculate (i) the 1 st Wasserstein distance (W 1 ), (ii) the 2 nd Wasserstein distance (W 2 ), and (iii) the bottleneck distance (W ∞ ) between the persistence diagrams obtained from the test data set of the SPHERES data and their resulting 2D latent representations. Even though our loss function is not optimising this distance, we observe in Table A .1 that the topological distance of our method ("TopoAE") is always the lowest among all the methods. In particular, it is always smaller than the topological distance of the latent space of the autoencoder architecture; this is true for all distance measures, even though W ∞ , for example, is known to be susceptible to outliers. Said experiment serves as a simple "sanity check" as it demonstrates that the changes induced by our method are beneficial in that they reduce the topological distance of the latent space to the original data set. For a proper comparison of topological features between the two sets of spaces, a more involved approach would be required, though.

A B C d 1 d 3 d 2 (a) X A B C d 2 d 3 d 1 (b) Z Figure A

.8: An undesirable configuration of the latent space of three non-collinear points. Pairwise distances are shown as dotted lines. This configuration results in equal persistence diagrams for X and Z. We prevent this by not explicitly minimising the distances between persistence diagrams but by instead including the persistence pairing information.

A.10 ALTERNATIVE LOSS FORMULATIONS:

Our choice of loss function was motivated by the observation that only aligning the persistence diagrams between mini-batches of X and Z can lead to degenerate or 'meaningless' latent spaces. As a simple example (see Figure A .8 for a visualisation), imagine three non-collinear points in the input space and the triangle they are forming. Now assume that the latent space consists of the same triangle (in terms of its side lengths) but with permuted labels. A loss term of the form

EQUATION

only measures the distance between persistence diagrams (which would be zero in this situation) and would not be able to penalise such a configuration. 

:

This means that we do not have to consider higher-dimensional topological features, making the calculation more efficient.

We use the term "matching" only to build intuition. Our approach does not calculate a matching in the sense of a bottleneck or Wasserstein distance between persistence diagrams.

When initialising a random latent space Z, the persistence pairing in the latent space will select random edges, resulting in only 1 expected matched edge (independent of mini-batch size) between the two pairings. Thus, only one edge (referring to one pairwise distance between two latent codes) could be used to update the encoding of these two data points.

https://osf.io/abuce/?view_only=f16d65d3f73e4918ad07cdd08a1a0d4b

For m = n, the two integrals switch their order as m(n -m) = 0 < n -1 (for n > 1).

Invertible transformations are available for PCA and all autoencoder-based methods.

Figures and Tables:

Figure 1: The Vietoris-Rips complex R (X) of a point cloud X at different scales 0 , 1 , and 2 . As the distance threshold increases, the connectivity changes. The creation and destruction of d-dimensional topological features is recorded in the d th persistence diagram D d .

Figure 2: Latent representations of the SPHERES data set. Only our method is capable of representing the complicated nesting relationship inherent to the data. TopoAE used a batch-size of 28. See Figure A.2 in the supplementary materials for an enlarged version.

Figure 3: Latent representations of the FASHION-MNIST (top row), MNIST (middle row), CIFAR-10 (bottom row) data sets. TopoAE batch-sizes (same data set order): (95, 126, 82). Please refer to Figures A.3, A.4, and A.5 in the supplementary materials for enlarged versions.

Figure A.1: Empirical convergence rate (mean) of the Hausdorff distance for a subsample of size m of 100 points in a d-dimensional space, following a standard normal distribution.

10) from which the original claim from Eq. 7 follows by taking probabilities on both sides.A.2 EMPIRICAL CONVERGENCE RATES OF d H X, X(m)

Figure A.1 depicts the mean of the convergence rate (mean) of the Hausdorff distance for a subsample of size m of 100 points in a d-dimensional space, following a standard normal distribution. We can see that the convergence rate is roughly similar, but shown on different absolute levels that depend on the ambient dimension. While bounding the convergence rate of this expression is feasible(Chazal et al., 2015a;b), it requires more involved assumptions on the measures from which X and X (m) are sampled. Additionally, we can give a simple bound using the diameter diam(X) := sup{dist(x, y) | x, y ∈ X}. We have d H X, X (m) ≤ diam(X) because the supremum is guaranteed to be an upper bound for the Hausdorff distance. This worst-case bound does not account for the sample size (or mini-batch size) m, though (see Theorem 2 for an expression that takes m into account).

Figure A.2: A depiction of all latent spaces obtained for the SPHERES data set. TopoAE used batch-size 28. This is an enlarged version of the figure shown in Section 5.2.

Figure A.3: Latent representations of the FASHION-MNIST data set. TopoAE used batch-size 95. This is a larger extension of the figure shown in Section 5.2.

Figure A.4: Latent representations of the MNIST data set. TopoAE used batch-size 126. This is a larger extension of the figure shown in Section 5.2.

Figure A.5: Latent representations of the CIFAR-10 data set. TopoAE used batch-size 82. This is a larger extension of the figure shown in Section 5.2.

Figure A.6: From our hyperparameter searches we scatter batch-sizes against three measures of interest: Topological Loss, Reconstruction Loss, and KL 0.1 , our objective for hyperparameter search. Additionally, we draw epoch-wise runtime estimates.

Figure A.7: A depiction of latent spaces obtained for the SPHERES data set with variational autoencoders (VAEs). Here, VAE represents a standard MLP-based VAE, whereas TopoVAE represents the same architecture plus our topological constraint.

Multiple evaluation metrics (Section 5.1). The hyperparameters of all tunable methods were selected to minimise the objective KL 0.1 . The winner is shown in bold and underlined, the runner-up in bold. Please refer to Supplementary TableA.2 for more σ scales and variance estimates.

Table content: <table><tr><td>Data set</td><td colspan="2">Method KL 0.01 KL 0.1</td><td colspan="4">KL 1 -MRRE -Cont -Trust -RMSE Data MSE</td></tr><tr><td/><td colspan="3">Isomap 0.181 0.420 0.00881</td><td>0.246 0.790 0.676</td><td>10.4</td><td>-</td></tr><tr><td/><td>PCA</td><td colspan="2">0.332 0.651 0.01530</td><td>0.294 0.747 0.626</td><td>11.8</td><td>0.9610</td></tr><tr><td>SPHERES</td><td>TSNE UMAP</td><td colspan="2">0.152 0.527 0.01271 0.157 0.613 0.01658</td><td>0.217 0.773 0.679 0.250 0.752 0.635</td><td>8.1 9.3</td><td>--</td></tr><tr><td/><td>AE</td><td colspan="2">0.566 0.746 0.01664</td><td>0.349 0.607 0.588</td><td>13.3</td><td>0.8155</td></tr><tr><td/><td colspan="3">TopoAE 0.085 0.326 0.00694</td><td>0.272 0.822 0.658</td><td>13.5</td><td>0.8681</td></tr><tr><td/><td>PCA</td><td colspan="2">0.356 0.052 0.00069</td><td>0.057 0.968 0.917</td><td>9.1</td><td>0.1844</td></tr><tr><td/><td>TSNE</td><td colspan="2">0.405 0.071 0.00198</td><td>0.020 0.967 0.974</td><td>41.3</td><td>-</td></tr><tr><td>F-MNIST</td><td>UMAP</td><td colspan="2">0.424 0.065 0.00163</td><td>0.029 0.981 0.959</td><td>13.7</td><td>-</td></tr><tr><td/><td>AE</td><td colspan="2">0.478 0.068 0.00125</td><td>0.026 0.968 0.974</td><td>20.7</td><td>0.1020</td></tr><tr><td/><td colspan="3">TopoAE 0.392 0.054 0.00100</td><td>0.032 0.980 0.956</td><td>20.5</td><td>0.1207</td></tr><tr><td/><td>PCA</td><td colspan="2">0.389 0.163 0.00160</td><td>0.166 0.901 0.745</td><td>13.2</td><td>0.2227</td></tr><tr><td/><td>TSNE</td><td colspan="2">0.277 0.133 0.00214</td><td>0.040 0.921 0.946</td><td>22.9</td><td>-</td></tr><tr><td>MNIST</td><td>UMAP</td><td colspan="2">0.321 0.146 0.00234</td><td>0.051 0.940 0.938</td><td>14.6</td><td>-</td></tr><tr><td/><td>AE</td><td colspan="2">0.620 0.155 0.00156</td><td>0.058 0.913 0.937</td><td>18.2</td><td>0.1373</td></tr><tr><td/><td colspan="3">TopoAE 0.341 0.110 0.00114</td><td>0.056 0.932 0.928</td><td>19.6</td><td>0.1388</td></tr><tr><td/><td>PCA</td><td colspan="2">0.591 0.020 0.00023</td><td>0.119 0.931 0.821</td><td>17.7</td><td>0.1482</td></tr><tr><td/><td>TSNE</td><td colspan="2">0.627 0.030 0.00073</td><td>0.103 0.903 0.863</td><td>25.6</td><td>-</td></tr><tr><td>CIFAR</td><td>UMAP</td><td colspan="2">0.617 0.026 0.00050</td><td>0.127 0.920 0.817</td><td>33.6</td><td>-</td></tr><tr><td/><td>AE</td><td colspan="2">0.668 0.035 0.00062</td><td>0.132 0.851 0.864</td><td>36.3</td><td>0.1403</td></tr><tr><td/><td colspan="3">TopoAE 0.556 0.019 0.00031</td><td>0.108 0.927 0.845</td><td>37.9</td><td>0.1398</td></tr></table>

1: Various topological distances calculated between the test data set and the corresponding latent space. For computational efficiency reasons, we used subsamples of size m = 500 and repeated the process 10 times. Individual cells thus contain a mean and a standard deviation.

Table content: <table><tr><td>Method</td><td>W 1</td><td>W 2</td><td>W ∞</td></tr><tr><td>Isomap</td><td colspan="3">4.32±0.037 0.477±0.0045 0.165±0.00096</td></tr><tr><td>PCA</td><td colspan="3">4.42±0.053 0.476±0.0046 0.158±0.00108</td></tr><tr><td>t-SNE</td><td colspan="3">4.38±0.038 0.478±0.0045 0.164±0.00094</td></tr><tr><td>UMAP</td><td colspan="3">4.47±0.042 0.478±0.0045 0.160±0.00092</td></tr><tr><td>AE</td><td colspan="3">3.99±0.037 0.469±0.0053 0.154±0.00128</td></tr><tr><td colspan="4">TopoAE 3.73±0.076 0.459±0.0055 0.152±0.00268</td></tr></table>

2: Extended version of the table from the main paper, showing more length scales and variance estimates.

Table content: <table><tr><td>0.01</td></tr><tr><td>KL</td></tr><tr><td>0.001</td></tr><tr><td>KL</td></tr><tr><td>Method</td></tr><tr><td>Data set</td></tr></table>

Peer review (iclr2020):
This paper shows how to train an autoencoder that preserves topologically relevant distance across varying length scales. It presents a new differentiable loss term for topological distance between input and latent space. The topological loss has a number of nice features: compatibility with minibatch training, and extension to higher input and latent space dimensionality fairly easily. Since this trainable topological loss can be applied in more general scenarios than simply auto-encoders or visualization, I think this well-written article is of general interest and worth publishing. The main ideas and related work are presented clearly in sections 1--4, and the experiments compare TopoAE with a variety of low-dimensional visualization methods. Embedding quality is evaluated with a wide variety of metrics on real and synthetic datasets highlighting the preservation of global and local topology. My only real issue, easily fixable, was that after the minibatch stability theory (Sec. 3.3), I really wanted to know what minibatch sizes were used in the experimental figures, and spent twenty minutes before admitting defeat. Please include the actual TopoAE minibatch size when presenting the visualization figures in the article and appendices. The appendices contain a wealth of useful reading and experiments, and the provided source code was clearly organized and useful to browse. I did not go line-by-line through the appendix proofs.