Year: 2021

Abstract:
Unsupervised feature learning often finds low-dimensional embeddings that capture the structure of complex data. For tasks for which expert prior topological knowledge is available, incorporating this into the learned representation may lead to higher quality embeddings. For example, this may help one to embed the data into a given number of clusters, or to accommodate for noise that prevents one from deriving the distribution of the data over the model directly, which can then be learned more effectively. However, a general tool for integrating different prior topological knowledge into embeddings is lacking. Although differentiable topology layers have been recently developed that can (re)shape embeddings into prespecified topological models, they have two important limitations for representation learning, which we address in this paper. First, the currently suggested topological losses fail to represent simple models such as clusters and flares in a natural manner. Second, these losses neglect all original structural (such as neighborhood) information in the data that is useful for learning. We overcome these limitations by introducing a new set of topological losses, and proposing their usage as a way for topologically regularizing data embeddings to naturally represent a prespecified model. We include thorough experiments on synthetic and real data that highlight the usefulness and versatility of this approach, with applications ranging from modeling high-dimensional single-cell data, to graph embedding.

heading: 1 INTRODUCTION
text: Motivation Modern data often arrives in complex forms that complicate their analysis. For example, high-dimensional data cannot be visualized directly, whereas relational data such as graphs lack the natural vectorized structure required by various machine learning models (Bhagat et al., 2011; Kazemi & Poole, 2018; Goyal & Ferrara, 2018). Representation learning aims to derive mathematically and computationally convenient representations to process and learn from such data. However, obtaining an effective representation is often challenging, for example, due to the accumulation of noise in high-dimensional biological expression data (Vandaele et al., 2021b). In other examples such as community detection in social networks, graph embeddings struggle to clearly separate communities with interconnections between them. In such cases, expert prior knowledge of the topological model may improve learning from, visualizing, and interpreting the data. Unfortunately, a general tool for incorporating prior topological knowledge in representation learning is lacking.
In this paper, we introduce such tool under the name of topological regularization. Here, we build on the recently developed differentiation frameworks for optimizing data to capture topological properties of interest (Brüel-Gabrielsson et al., 2020; Solomon et al., 2021; Carriere et al., 2021). Unfortunately, such topological optimization has been poorly studied within the context of representation learning. For example, the used topological losses are indifferent to any structure other than topological, such as neighborhoods, which may be useful for learning. Therefore, topological optimization often destructs natural and informative properties of the data in favor of the topological loss.
Our proposed method of topological regularization effectively resolves this by learning an embedding representation that incorporates the topological prior. As we will see in this paper, these priors can be directly postulated through topological loss functions. For example, if the prior is that the data lies on a circular model, we design a loss function that is lower whenever a more prominent cycle is present in the embedding. By extending the previously suggested topological losses to fit a wider set of models, we show that topological regularization effectively embeds data according to a variety of topological priors, ranging from clusters, cycles, and flares, to any combination of these.
Related Work Certain methods that incorporate topological information into representation learning have already been developed. For example, Deep Embedded Clustering (Xie et al., 2016) simultaneously learns feature representations and cluster assignments using deep neural networks. Constrained embeddings of Euclidean data on spheres have also been studied by Bai et al. (2015). However, such methods often require extensive development for one particular kind of input data and topological model. Contrary to this, incorporating topological optimization into representation learning provides a simple yet versatile approach towards combining any data embedding method with topological priors, that generalizes well to any input data as long as the output is a point cloud.
Topological autoencoders (Moor et al., 2020) already combine topological optimization with a data embedding procedure. The main difference here is that the topological information used for optimization is obtained from the original high-dimensional data, and not passed as a prior. While this may sound as a major advantage—and certainly can be as shown by Moor et al. (2020)— obtaining such topological information heavily relies on distances between observations, which are often meaningless and unstable in high dimensions (Aggarwal et al., 2001). Furthermore, certain constructions such as the α-filtration obtained from the Delaunay triangulation—which we will use extensively and is further discussed in Appendix A—are expensive to obtain from high-dimensional data (Cignoni et al., 1998), and are therefore best computed from the low-dimensional embedding.
Our work builds on a series of recent papers (Brüel-Gabrielsson et al., 2020; Solomon et al., 2021; Carriere et al., 2021), which showed that topological optimization is possible in various settings and developed the mathematical foundation thereof. However, studying the use of topological optimization for data embedding and visualization applications, as well as the new losses we develop therefore and the insights we derive from them in this paper, are to best of our knowledge novel.
Contributions We include a sufficient background on persistent homology—the method behind topological optimization—in Appendix A (note that all of its concepts important for this paper are summarized in Figure 1). We summarize the previous idea behind topological optimization of point clouds (Section 2.1). We introduce a new set of losses to model a wider variety of shapes in a natural manner (Section 2.2). We show how these can be used to topologically regularize embedding methods for which the output is a point cloud (Equation (1)). We include experiments on synthetic and real data that show the usefulness and versatility of topological regularization, and provide additional insights into the performance of data embedding methods (Section 3 & Appendix B). We discuss open problems in topological representation learning and conclude on our work (Section 4).

heading: 2 METHODS
text: The main purpose of this paper is to present a method to incorporate prior topological knowledge in a point cloud embedding E (dimensionality reduction, graph embedding, . . . ) of a data set X. As will become clear below, these topological priors can be directly postulated through topological loss functions Ltop. Then, the goal is to find an embedding that that minimizes a total loss
Ltot(E,X) := Lemb(E,X) + λtopLtop(E), (1)
where Lemb is a loss that aims to preserve structural attributes of the original data, and λtop > 0 controls the strength of topological regularization. Note that, X itself is not required to be a point cloud, or reside in the same space as E, which is especially useful for representation learning.
In this section, we mainly focus on topological optimization of point clouds, that is, the loss Ltop. The basic idea behind this recently introduced method—as presented by Brüel-Gabrielsson et al. (2020)—is illustrated in Section 2.1. We show that direct topological optimization may neglect important structural information such as neighborhoods, which can effectively be resolved through (1). Hence, as we will also see in Section 3, while representation learning may benefit from topological losses for incorporating prior topological knowledge, topological optimization itself may also benefit from other structural losses to represent the topological prior in a more truthful manner. Nevertheless, some topological models remain difficult to naturally represent through topological optimization. Therefore, we introduce a new set of topological losses, and provide an overview of how different models can be postulated through them in Section 2.2. Experiments with and comparisons to topological regularization of embeddings through (1) will be presented in Section 3.

heading: 2.1 BACKGROUND ON TOPOLOGICAL OPTIMIZATION OF POINT CLOUDS
text: Topological optimization is performed through a topological loss function evaluated on one or more persistence diagrams (Barannikov, 1994; Carlsson, 2009). These diagrams—obtained through persistent homology as formally discussed in Appendix A—summarize all from the finest to coarsest topological holes (connected components, cycles, voids, . . . ) in the data, as illustrated in Figure 1.
While methods that learn from persistent homology are now both well developed and diverse (Pun et al., 2018), optimizing the data representation for the persistent homology thereof only gained recent attention (Brüel-Gabrielsson et al., 2020; Solomon et al., 2021; Carriere et al., 2021). Persistent homology has a rather abstract mathematical foundation within algebraic topology (Hatcher, 2002), and its computation is inherently combinatorial (Barannikov, 1994; Zomorodian & Carlsson, 2005). This complicates working with usual derivatives for optimization. To accommodate for this, topological optimization makes use of Clarke subderivatives (Clarke, 1990), whose applicability to persistence builds on arguments from o-minimal geometry (van den Dries, 1998; Carriere et al., 2021). Fortunately, thanks to the recent work of Brüel-Gabrielsson et al. (2020) and Carriere et al. (2021), powerful tools for topological optimization have been developed for software libraries such as PyTorch and TensorFlow, allowing their application without deeper knowledge of these subjects.
Topological optimization optimizes the data representation with respect to topological information summarized by its persistence diagram(s) D. We will use the approach by Brüel-Gabrielsson et al. (2020), where (birth, death) tuples (b1, d1), (b2, d2), . . . , (b|D|, d|D|) in D are first ordered by decreasing persistence dk − bk. The points (b,∞), usually plotted on top of the diagram such as in Figure 1b, form the essential part of D. The points with finite coordinates form the regular part of D. In case of point clouds, one and only one topological hole, i.e., a connected component born at time α − 0, will always persist indefinitely. Other gaps and holes will eventually be filled (Figure 1). Thus, we only optimize for the regular part in this paper. This is done through a topological loss function, which for a choice of i ≤ j (which, along with the dimension of topological hole, will specify our topological prior as we will see below) and a function g : R2 → R, is defined as
Ltop(D) := j∑
k=i,dk<∞
g(bk, dk), where d1 − b1 ≥ d2 − b2 ≥ . . . . (2)
It turns out that for many useful definitions of g, Ltop(D) has a well-defined Clarke subdifferential with respect to the parameters defining the filtration from which the persistence diagram D is ob-
tained. In this paper, we will consistently use the α-filtration as shown in Figure 1a (see Appendix A for its formal definition), and these parameters are entire point clouds (in this paper embeddings) E ∈ (Rd)n of size n in the d-dimensional Euclidean space. Ltop(D) can then be easily optimized with respect to these parameters through stochastic subgradient algorithms (Carriere et al., 2021).
As it directly measures the prominence of topological holes, we let g : R2 → R : (b, d) 7→ µ(d− b) be proportional to the persistence function. By ordering the points by persistence, Ltop is a function of persistence, i.e., it is invariant to permutations of the points in D (Carriere et al., 2021). The factor of proportionality µ ∈ {1,−1} indicates whether we want to minimize (µ = 1) or maximize (µ = −1) persistence, i.e, the prominence of topological holes. Thus, µ determines whether more clusters, cycles, . . . , should be present (µ = −1) or missing (µ = 1). The loss (2) then reduces to
Ltop(E) := Ltop(D) = µ j∑
k=i,dk<∞
(dk − bk) , where d1 − b1 ≥ d2 − b2 ≥ . . . . (3)
Here, the data matrix E (in this paper the embedding) defines the diagram D through persistent homology of the α-filtration of E, and a persistence (topological hole) dimension to optimize for.
For example, consider (3) with i = 2, j = ∞, µ = 1, restricted to 0-dimensional persistence (measuring the prominence of connected components) of the α-filtration. Figure 2 shows the data from Figure 1 optimized for this loss function for various epochs. The optimized point cloud quickly resembles a single connected component for smaller numbers of epochs. This is the single goal of the current loss, which neglects all other structural structural properties of the data such as its underlying cycles (e.g., the circular hole in the ‘R’) or local neighborhoods. Larger numbers of epochs mainly affect the scale of the data. While this scale has an absolute effect on the total persistence, thus, the loss, the point cloud visually represents a single connected topological component equally well. We also observe that while local neighborhoods are preserved well during the first epochs simply by nature of topological optimization, they are increasingly distorted for a larger number of epochs.

heading: 2.2 NEWLY PROPOSED TOPOLOGICAL LOSS FUNCTIONS
text: In this paper, the prior topological knowledge incorporated into the point cloud data matrix embedding E is directly postulated through a topological loss function. For example, letting D be the 0-dimensional persistence diagram of E (H0 in Figure 1b), and i = 2, j = ∞, and µ = 1 in (3), corresponds to the criterion that E should represents one closely connected component, as illustrated in Figure 2. Therefore, we often regard a topological loss as a topological prior, and vice versa.
Although persistent homology effectively measures the prominence of topological holes, topological optimization is often ineffective for representing such holes in a natural manner, for example, through sufficiently many data points. An extreme example of this are clusters, despite the fact that they are captured through the simplest form of persistence, i.e., 0-dimensional. This is shown in Figure 3. Optimizing the point cloud in Figure 3a for (at least) two clusters can be done by defining Ltop(E) as in (3), letting D be the 0-dimensional persistence diagram of E, i = j = 2, and µ = −1. However, we observe that topological optimization simply displaces one point away from all other points (Figure 3b). Purely topologically, this is indeed a correct representation of two clusters.
To represent topological holes through more points, we propose topological optimization for the loss
L̃top(E) := E (Ltop ({x ∈ S : S is a random sample of E with sampling fraction fS})) , (4) where Ltop is defined as in (3). In practice, during each optimization iteration, L̃top is approximated by the mean of Ltop evaluated over nS random samples of E. Figure 3 shows the result
for a sampling fraction fS = 0.1 and repeats nS = 1. The new data representation visualizes the clusters already well and far more naturally. An added benefit of the new loss (4) is that topological optimization can be conducted significantly faster, as the α-filtration and persistent homology are evaluated on smaller samples (a computational analysis can be found at the end of Appendix A).
In summary, various topological priors can now be formulated through topological losses as follows.
k-dimensional holes Optimizing for k-dimensional holes (k = 0 for clusters), can generally be done through (3) or (4), by letting D be the corresponding k-dimensional persistence diagram. The terms i and j in the summation are used to express how many holes one exactly, at least, or at most wants. Finally, µ can be chosen to either decrease (µ = 1) or increase (µ = −1) persistence.
Flares Persistent homology is invariant to certain topological changes. For example, both a linear ‘I’-structured model and a bifurcating ‘Y’-structured model consist of one connected component, and no higher-dimensional holes. These models are indistinguishable based on the (persistent) homology thereof, even though they are topologically different in terms of their singular points.
Capturing such additional topological phenomena is possible through a refinement of persistent homology under the name of functional persistence, also well discussed and illustrated by Carlsson (2014). The idea is that instead of evaluating persistent homology on a data matrix E, we evaluate it on a subset {x ∈ E : f(x) ≤ τ} for a well chosen function f : E → R and hyperparameter τ . Inspired by this approach, for a diagram D of a point cloud E, we propose the topological loss
L̃top(E) := Ltop ({x ∈ E : fE(x) ≤ τ}) , informally denoted [Ltop(D)]f−1E ]−∞,τ ] , (5)
where f is a real-valued function on E, possibly dependent on E—which changes during optimization—itself, τ a hyperparameter, and Ltop is an ordinary topological loss as defined by (3). In particular, we will focus on the case where f equals a scaled centrality measure on E:
fE ≡ EE :≡ 1− gE
max gE , where gE(x) := ∥∥∥∥∥∥x− 1|E| ∑ y∈E y ∥∥∥∥∥∥ . (6) For τ ≥ 1, L̃top(E) = Ltop(E). For sufficiently small τ > 0, L̃top evaluates Ltop on the points ‘far away’ from the center of E. As we will see in the experiments, this is especially useful in conjunction with 0-dimensional persistence to optimize for flares, i.e., star-structured shapes.
Combinations Naturally, through linear combination of loss functions, different topological priors can be combined, e.g., if we want the represented model to both be connected and include a cycle.

heading: 3 EXPERIMENTS
text: In this section, we show how our proposed topological regularization of data embeddings (1) leads to a powerful and versatile approach for representation learning. In particular, we show that
• embeddings benefit from prior topological knowledge through topological regularization; • conversely, topological optimization may also benefit from incorporating structural infor-
mation as captured through embedding losses, leading to more qualitative representations; • subsequent learning tasks may benefit from expert prior topological knowledge.
In Section 3.1, we show how topological regularization improves standard PCA dimensionality reduction and allows better understanding of its performance when noise is accumulated over many dimensions. In Section 3.2, we present applications to high-dimensional single-cell trajectory data and graph embedding. Quantitative results are discussed in Section 3.3.
Topological optimization was performed in Pytorch, using code adapted from Brüel-Gabrielsson et al. (2020). Data sizes, hyperparameters, losses, and optimization times are summarized in Tables 1 & 2. Code for this project is provided with our submission and will be made available on github.
Appendix B presents supplementary experiments that show topological regularization is consistent with earlier topological data analysis of the Harry Potter network (Section B.1), a domain application on real biological cells in the ‘cell cycle’ differentiation model (Section B.2), and how topological regularization reacts to different topological losses and hyperparameters (Section B.3).

heading: 3.1 SYNTHETIC DATA
text: We sampled 50 points uniformly from the unit circle in R2. We then added 500-dimensional noise to the resulting data matrix X , where the noise in each dimension is sampled uniformly from [−0.45, 0.45]. Since the additional noisy features are irrelevant to the topological (circular) model, an ideal projection embedding X is its restriction to its first two data coordinates (Figure 4a).
However, it is probabilistically unlikely that that the irrelevant features will have a zero contribution to a PCA embedding of the data (Figure 4b). Measuring the feature importance of each feature as the sum of its two absolute contributions (the loadings) to the projection, we observe that most of the 498 irrelevant features have a small nonzero effect on the PCA embedding (Figure 5). Intuitively, each added feature slightly shifts the projection plane away from the plane spanned by the first two coordinates. As a result, the circular hole is less prominent in the PCA embedding of the data. Note that we observed this to be a notable problem for ‘small n large p’ data sets, as similar to other machine learning models and as also recently studied by Vandaele et al. (2021a), more data can significantly accommodate for the effect of noise and result in a better embedding model on its own.
We can regularize this embedding using a topological loss function Ltop measuring the persistence of the most prominent 1-dimensional hole in the embedding (i = j = 1 in (3)). For a simple Pytorch compatible implementation, we used Lemb(W ,X) := MSE ( XWW T ,X ) , as to minimize the reconstruction error between X and its linear projection obtained through W . To this, we added the loss 104L⊥(W ), where L⊥(W ) := ∥W TW − I2∥2 is used to encourage orthonormality of the matrix W to be optimized, initialized with the PCA-loadings. The resulting embedding is shown in Figure 4d (here L⊥(W ) ∼ 0.03), which better captures the circular hole. Furthermore, we see that irrelevant features now more often contribute less to the embedding according to W (Figure 5).
For comparison, Figure 4c shows the optimization of W without accounting for the reconstruction loss (L⊥(W ) still included, here L⊥(W ) ∼ 0.2). From this and also from Figure 5, we observe that W struggles more to converge to the correct projection, resulting in a less prominent hole.

heading: 3.2 REAL DATA
text: Circular Cell Trajectory Data We considered a single-cell trajectory data set of 264 cells in a 6812-dimensional gene expression space (Cannoodt et al., 2018; Saelens et al., 2019). This data can be considered a snapshot of the cells at a fixed time. The ground truth is a circular model connecting three cell groups through cell differentiation (see also Section B.2). It has been shown by Vandaele (2020) that real single-cell data with such model are difficult to embed in a circular manner.
To explore this, we repeated the experiment with the same losses as in Section 3.1 on this data, where the (expected) topological loss is now modified through (4) with fS = 0.25, and nS = 1. From Figure 6a, we see that while the ordinary PCA embedding does somehow respect the positioning of the cell groups (marked by their color), it indeed struggles to embed the data in a manner that visualizes the present cycle. However, as shown in Figure 6c, by topologically regularizing the embedding (here L⊥(W ) ∼ 4e−3) we are able to embed the data much better in a circular manner. Figure 6b shows the optimization of W without the reconstruction loss (L⊥(W ) ∼ 4e−3). The embedding is similar to the one in Figure 6c, with the pink colored cell group slightly more dispersed.
Bifurcating Cell Trajectory Data We considered a second cell trajectory data set of 154 cells in a 1770-dimensional expression space (Cannoodt et al., 2018). The ground truth here is a bifurcating model connecting four different cell groups through cell differentiation. However, this time we used the UMAP loss for the embeddings. We used a topological loss Ltop ≡ Lconn − Lflare, where Lconn measures the total (sum of) finite 0-dimensional persistence in the embedding to encourage connectedness of the representation, and Lflare is as in (5), measuring the persistence of the third most prominent 0-dimensional hole in {y ∈ E : EE(y) ≤ 0.75}, where EE is as in (6). Thus, Lflare is used to optimize for a ‘flare’ with (at least) three clusters away from the embedding mean. We observe that while the ordinary UMAP embedding is more ‘blobby’ (Figure 7a), the topologically regularized embedding is more constrained towards a connected bifurcating shape (Figure 7c).
For comparison, we conducted topological optimization for the loss Lflare of the initialized UMAP embedding without the UMAP embedding loss. The resulting embedding is now more fragmented (Figure 7b). We thus see that topological optimization may also benefit from the embedding loss.
Graph Embedding The topological loss in (1) can be evaluated on any embedding, and does not require a point cloud as original input. We can thus use topological regularization for embedding a graph G, to learn a representation of the nodes of G in Rd that well respects properties of G.
To explore this, we considered the Karate network (Zachary, 1977), a well known and studied network within graph mining that consists of two different communities. The communities are represented by two key figures (John A. and Mr. Hi), as shown in Figure 8a. To embed the graph, we used a DeepWalk variant adapted from Dagar et al. (2020). While the ordinary DeepWalk em-
bedding (Figure 8b) well respects the ordering of points according to their communities, the two communities remained close to each other. We thus regularized this embedding for (at least) two clusters using the topological loss as defined by (4), where Ltop measures the persistence of the second most prominent 0-dimensional hole, and fS = 0.25, nS = 10. The resulting embedding (Figure 8d) now nearly perfectly separates the two ground truth communities present in the graph.
Topological optimization of the initialized DeepWalk embedding with the same topological loss but without the DeepWalk loss creates some natural community structure, but also results in a few outliers (Figure 8c). Thus, although our introduced loss (4) enables more natural topological modeling to some extent, we again observe that using this in conjunction with embedding losses, i.e., our proposed method of topological regularization, leads to the best qualitative results.

heading: 3.3 QUANTITATIVE EVALUATION
text: Table 3 summarizes the embedding and topological losses we obtained for the ordinary embeddings, the topologically optimized embeddings (initialized with the ordinary embeddings, but not using the embedding loss), as well as for the topologically regularized embeddings. As one would expect, topological regularization balances the embedding losses between the embedding losses of the ordinary and topologically optimized embeddings. More interestingly, topological regularization may actually result in a more optimal, i.e., lower topological loss than topological optimization only, here in particular for the synthetic cycle data and Harry Potter graph. This suggest that combining topological information with other structural information may facilitate convergence to the correct embedding model, as we also qualitatively confirmed for these data sets (see also Section B.1). We also observe that there are more significant differences in the obtained topological losses than in the embedding losses with and without regularization. This suggests that the optimum region for the embedding loss may be somewhat flat with respect to the corresponding region for the topological loss. Thus, slight shifts in the local embedding optimum, e.g., as caused by noise, may result in much worse topological embedding models, which can be resolved through topological regularization.
We evaluated the quality of the embedding visualizations presented in this section, by assessing how informative they are for predicting ground data truth labels. For the Synthetic Cycle data, these labels are the 2D coordinates of the noise-free data on the unit circle in R2, and we used a multi-ouput support vector regressor model. For the cell trajectory data and Karate network, we used the ground truth cell groupings and community assignments, respectively, and a support vector machine model. All points in the 2D embeddings were then split into 90% points for training and 10% for testing. Consecutively, we used 5-fold CV on the training data to tune the regularization hyperparameter C ∈ {1e − 2, 1e − 1, 1, 1e1, 1e2}. Other settings were the default from SCIKIT-LEARN. The performance of the final tuned and trained model was then evaluated on the test data, through the r2 coefficient of determination for the regression problem, and the accuracy for all classification problems. The averaged test performance metrics and their standard deviations, obtained over 100 random train-test splits, are summarized in Table 4. Although topological regularization consistently leads to more informative visualization embeddings, quantitative differences may be minor, as we can see from the figures in this section that most local similarities are preserved with regularization.

heading: 4 DISCUSSION AND CONCLUSION
text: We proposed a new approach for representation learning under the name of topological regularization, which builds on the recently developed differentiation frameworks for topological optimization. This led to a versatile and effective way for embedding data according to expert prior topological knowledge, directly postulated through (some newly introduced) topological loss functions.
A clear limitation of topological regularization is that prior topological knowledge is not always available. How to select the best from a list of priors is thus open to further research (see also Section B.3). Furthermore, designing topological loss functions currently requires some understanding of persistent homology. It may be useful to study how to facilitate that design process for lay users. From a foundational perspective, our work provides new research opportunities into extending the developed theory for topological optimization (Carriere et al., 2021) to our newly introduced losses and their integration into data embeddings. Finally, topological optimization based on combinatorial structures other than the α-complex may be of theoretical and practical interest. For example, point cloud optimization based on graph-approximations such as the minimum spanning tree (Vandaele et al., 2021b), or varying the functional threshold τ in the loss (5) alongside the filtration time (Chazal et al., 2009), may lead to new topological loss functions with fewer hyperparameters.
Nevertheless, through our approach, we already provided new and important insights into the performance of embedding methods, such as their potential inability to converge to the correct topological model due to the flatness of the embedding loss near its (local) optimum, with respect to the topological loss. Furthermore, we showed that including prior topological knowledge provides a promising way to improve consecutive—even non-topological—learning tasks (see also Section B.2). In conclusion, topological regularization enables both improving and better understanding representation learning methods, for which we provided and thoroughly illustrated the first directions in this paper.

heading: B SUPPLEMENTARY EXPERIMENTS
text: In this part of the appendix, we include additional experiments that illustrate the usefulness and effectiveness of topologically regularization. Section B.1 includes an additional qualitative evaluation on the Harry Potter network, which shows that its topologically regularized graph embedding, optimized for a circular prior, is consistent with earlier and independent topological data analysis of this network (Vandaele et al., 2020). In Section B.2, we show how including expert prior topological information in the embedding procedure may facilitate automated pseudotime inference in our real cyclic single-cell expression data presented in the main paper. Finally, in Section B.3, we show how topological regularization reacts to different topological loss functions, either designed to model the same prior topological information, or different and potentially wrong prior information.
B.1 ADDITIONAL QUALITATIVE EVALUATION ON THE HARRY POTTER NETWORK
We considered an additional experiment on the Harry Potter graph obtained from https:// github.com/hzjken/character-network. This graph is composed of characters from the Harry Potter novel (these constitute the nodes in the graph), and edges marking friendly relationships between them (Figure 11). Only the largest connected component is used. This graph has previously been analyzed by Vandaele et al. (2020), who identified a circular model therein that transitions between the ‘good’ and ‘evil’ characters from the novel.
We furthermore regularized this embedding using a topological loss function Ltop that measures the persistence of the most prominent 1-dimensional hole in the embedding (see also Table 2 in the main paper), the result of which is shown in Figure 12c. Interestingly, the topologically regularized embedding now better captures the circularity of the model identified by Vandaele et al. (2020), and focuses more on distributing the characters along it. Note that although this model is included in the visualizations, it is not used to derive the embeddings, nor is it derived from them.
For comparison, Figure 12b shows the result of optimizing the initialized ordinary graph embedding for the same topological loss, but without the graph embedding loss. We observe that this results in a sparse enlarged cycle. Most characters are positioned poorly along the circular model, and concentrate near a small region. Interestingly, even though we only optimized for the topological loss here, it is actually less optimal, i.e., higher, than when we applied topological regularization (see also Table 3 in the main paper). This is also a result from the sparsity of the circle, which constitutes to a larger birth time, and thus a lower persistence, of the corresponding hole.
To embed the Harry Potter graph, we used a simple graph embedding model where the sigmoid of the inner product between embedded nodes captures the (Bernoulli) probability of an edge occurrence (Rendle et al., 2020). Thus, this probability will be high for nodes close to each other in the embedding, and low for more distant nodes. These probabilities are then optimized to match the binary edge indicator vector. Figure 12a shows the result of this embedding, along with the circular model presented by Vandaele et al. (2020). For clarity, character labels are only annotated for a subset of the nodes (the same as by Vandaele et al. (2020)).
B.2 APPLICATION: PSEUDOTIME INFERENCE IN CELL TRAJECTORY DATA
Single-cell omics data include various types of data collected on cellular level, such as transcriptomics, proteomics and epigenomics data. Studying the topological model underlying the data may lead to better understanding of the dynamic processes of biological cells and the regulatory interactions involved therein. Such dynamic processes can be modeled through trajectory inference (TI) methods, also called pseudotime analysis, which order cells along a trajectory based on the similarities in their expression patterns (Saelens et al., 2019).
For example, the cell cycle is a well known biological differentiation model that takes place in a cell as it grows and divide. The cell cycle consists of different stages, namely growth (G1), DNA synthesis (S), growth and preparation for mitosis (G2), and mitosis (M), where the latter two are often grouped together in a G2M stage. Hence, by studying expression data of cells that follow such differentiation model, one may identify the genes involved in and between particular stages of the cell cycle (Liu et al., 2017). Pseudotime analysis allows such study by assigning to each cell a time during the differentiation process in which it occurs, and thus, the relative positioning of all cells within the cell cycle model.
Hence, the analysis of single cell cycle data constitutes a problem where prior topological information is available. As the signal-to-noise ratio is commonly low in high-dimensional expression data (Libralon et al., 2009; Zhang et al., 2021), this data is usually preprocessed through a dimensionality reduction method prior to automated pseudotime inference (Cannoodt et al., 2016; Saelens et al., 2019). Topological regularization provides a tool to enhance the desired topological signal during the embedding procedure, and as such, facilitate automated inference that depends on this signal.
To explore this, we used persistent homology for an automated (cell) cycle and pseudotime inference method, with and without topological regularization during the PCA embedding of the data. The data we use for this experiment is our real cell cycle data presented in the main paper, also previously analyzed by Buettner et al. (2015). The topological regularization procedure we followed here is the same as in the main paper, i.e., we regularized for a circular prior (see also Tables 1 & 2 in the main paper). Our automated pseudotime inference method now consists of the following steps.
(a) The representation of the most prominent cycle obtained through persistent homology. (b) Orthogonal projection of the embedded data onto the cycle representation. (c) The pseudotimes inferred from the orthogonal projection, quantified on a continuous color scale.
(a) The representation of the most prominent cycle obtained through persistent homology. (b) Orthogonal projection of the embedded data onto the cycle representation. (c) The pseudotimes inferred from the orthogonal projection, quantified on a continuous color scale.
Figure 14: Automated pseudotime inference of real cell cycle data through persistent homology, from the topologically regularized PCA embedding of the data.
1. First, a representation of the most prominent cycle in the embedding is obtained through persistent homology from the α-filtration, using the Python software library Dionysus (https://pypi.org/project/dionysus/). It can be seen as a circular representation—discretized in edges between data points—of the point in the 1stdimensional persistence diagram that corresponds to the most persisting cycle (Figures 13a & 14a).
2. An orthogonal projection of the embedded data onto the representative cycle is obtained. This is an intermediate step to derive continuous pseudotimes from a discretized topological representation, as earlier described by Saelens et al. (2019) (Figures 13b & 14b).
3. The lengths between points on the orthogonal projection are used to obtain circular pseudotimes between 0 and 2π (Figures 13c & 14c).
Note that within the scope of the current paper, we do not advocate this to be a new and generally applicable automated single cell cycle and pseudotime inference method. However, we chose this particular method because it illustrates well what persistent homology identifies in the data, and what we aim to optimize through topological regularization.
For example, we observe that the (representation of) the most prominent cycle in the ordinary PCA embedding is rather spurious, and mostly linearly separates G1 and G2M cells. Projecting cells onto the edges of this cycle, places the majority of the cells on a single edge (Figure 13b). The resulting pseudotimes are mostly continuous only for cells projected onto this edge, whereas they are very discretized for all other cells (Figure 13c). However, by incorporating our prior topological knowledge into the PCA embedding, the (representation of) the most prominent cycle in the embedding now much better characterizes the transmission between the G1, G2M, and S stages in the cell cycle model (Figure 14a). The automated procedure for pseudotime inference now also reflects a more continuous transmission between the cell stages (Figures 14b & 14c).
B.3 TOPOLOGICAL REGULARIZATION FOR DIFFERENT LOSS FUNCTIONS
In its current stage, postulating topological prior information is directly done by designing a topological loss function. This requires some understanding of persistent homology, and how to facilitate that design process for lay users is open to further research. Therefore, it may be useful to explore how topological regularization reacts to different (and potentially even wrong) prior topological information and loss functions, which is the topic of this section.
Different Loss Functions for the same Topological Prior (varying fS ,nS , or τ ) As we have introduced new topological loss functions (4) and (5) in the main paper, it may be useful to investigate how topological regularization reacts to different choices of the hyperparameters that are important for these losses. These are the sampling fraction fS and number of repeats nS in (4), and the functional threshold τ in (5). Note that although changing these hyperparameters may result in changing the geometrical characteristics that are specified, the topological characteristics may be considered to roughly remain the same with respect to the chosen topological loss function.
Varying fS and nS Figure 15 shows the topologically regularized PCA embedding of our real single-cell data following the cell cycle differentiation model from the main paper and also discussed in Section 13c, for varying choices of the values fS and nS . The topological loss function and all other hyperparameters are identical to those in the main paper for this data (see also Tables 1 & 2 in the main paper).
We observe that both by increasing fS and nS , the cycle tends to be more well defined, i.e., we can better visually deduce a subset of the data that captures a nearly perfect circle in the plane. Nevertheless, we also conclude that using more data during topological regularization, i.e., higher values of fS , does not necessarily result in better embeddings. Indeed, this gives a higher chance that the topological loss function will focus more on spurious holes during optimization, for example as illustrated by Figure 13a. From Figure 15—although more difficult to deduce due to overlapping points—we see that higher values of fS generally result in slightly worse clustering, most notably distributing the pink colored cells (or as in Section B.2 discussed, G1 cells) more across the entire circular model.
An added advantage of our introduced topological loss function (4), which is approximated through a sampling strategy, is that we can use lower sampling fractions fS not only to avoid optimizing for spurious holes, but also conduct this optimization significantly faster, as also discussed in our computational analysis in Section A. This can also be seen from Figure 16, which compares the computation times in seconds of the topologically regularized embeddings for different sampling fractions fS and number of repeats nS .
Varying τ Figure 17 shows the topologically regularized UMAP embedding of our real single-cell data following a bifurcating cell differentiation model from the main paper, for varying choices of the functional threshold τ . Recall that 0 ≤ τ ≤ 1 was a threshold on the normalized distance of the embedded points to their mean as defined in (6) in the main paper. Intuitively, τ specified how close we can look to our embedding mean, with higher values of τ allowing more points (thus closer to the embedding mean) to be included when optimizing for three clusters away from the center.
This is consistent with what we observe from Figure 17. While little differences can be observed between τ = 0.25 and τ = 0.5, we observe that for τ = 0.75 (which was the value chosen in our main paper), points near the center of bifurcation are more stretched towards the leaves, i.e., endpoints, in the embedded model. This agrees with the fact that for higher values of τ , more points that are close to the center are included when optimizing for three clusters away from the center.
Different Loss Functions for the same Topological Prior (varying g) For one particular type of prior topological information, we may also design different topological loss functions by varying the real-valued function g evaluated on the points in the persistence diagram in the loss (2) from the main paper. In particular, our topological loss function (3) in the main paper, where we let g :≡ bk − dk, is inspired by the topological loss function
Ltop(E) := Ltop(D) = µ |Dreg|∑ i=k
(bk,dk)∈Dreg
(dk − bk)p ( dk + bk
2
)q , (7)
introduced by Brüel-Gabrielsson et al. (2020), and more formally investigated within an algebraic setting by Adcock et al. (2013). Here, p and q are hyperparameters that control the strength of penalization, whereas i and the choice of persistence diagram (or homology dimension) is used to postulate the topological information of interest. The choice of µ ∈ {1,−1} determines whether we want to increase (µ = −1) or decrease (µ = 1) the topological information for which we designed the topological loss function in the data. The term (dk−bk)p can be used to control the prominence, i.e., the persistence, of the most prominent topological holes, whereas the term ( dk+bk
2
)q can be
used to control whether prominent topological features persist either early or later in the filtration.
Our topological loss function (3) in the main paper is identical to (7) with p = 1 and q = 0, and the additional change that we sum to the hyperparameter j instead of |Dreg|, which allows for even more flexible prior topological information. Note that this one of the most simple and intuitive, yet flexible, manners to postulate topological loss functions. Indeed, (dk − bk)p=1 directly measures the the persistence of the k-th most prominent topological hole for a chosen dimension of interest. The role of q ̸= 0 is to be further investigated within the context of topological regularization and formulating prior topological information.
It may already be clear that the same topological information could be postulated through different nonzero choices for p > 0. To explore this, we considered two other topologically regularized UMAP embeddings of our real bifurcating single cell data presented in the main paper, for the same topological loss function where only the parameter p is varied, i.e.,
Ltop(E) = ∞∑ k=2 (dk − bk)p − [(d3 − b3)p]E−1E ]−∞,0.75] . (8)
Figure 18 shows the resulting embeddings for p = 12 and p = 2, as well as the original embedding for p = 1 from the main paper for easy comparison.
In the first term of our topological loss (8), higher values of p will more strongly penalize larger gaps between neighboring points. This can be seen in Figure 18 for the majority of the embedded points, which are more pulled towards the leaves (endpoints) of the represented topological model. However, this results in a few more displaced points near the center of bifurcation, which try to remain connected as we simultaneously fit the UMAP loss function. Indeed, in the main paper we have seen that without the UMAP loss function, the model represented in the topologically optimized embedding will be more fragmented. Hence, in the current experiment, larger values of p in the term ∑∞ k=2(dk − bk)p seem to better motivate many smaller gaps (at the cost of a few larger gaps), whereas lower values appear to motivate more evenly spaced points. The fact that more points are pulled from the center towards the leaves for larger values of p and not from the leaves towards the center, is consistent with the fact that in the original high-dimensional space, more points are far away from each other and the center than they are close, as previously shown by Vandaele (2020) for this particular data set.
Nevertheless, although different values of p appear to affext the overall spacing between points, we see the embedded topological shape remains recognizable and consistent (Figure 18).
Different Loss Functions for different Topological Priors (example for PCA) As a user may not always have strong expert prior topological knowledge available, it may be useful to study how topological regularization reacts to different topological loss functions that are designed to model weaker or even wrong topological information. Indeed, if the impact of topological regularization on embeddings may be too strong, the user might wrongfully conclude that a particular topological model is present in the data when it can be easily fitted according to any specified shape.
To explore this, we considered various additional topological functions to topologically regulariz the PCA embedding of our synthetic data set of points lying on a circle with added noise in high dimension, presented in the main paper. Note that here we know that there is exactly one topological model, i.e., a circle, truthfully present in the data. We performed the topologically regularized PCA embedding of this data for the exact same hyperparameters as summarized in Table 1 in the main paper, but now for the three other topological loss functions as summarized in Table 5. The intuition of these topological loss functions, in order of appearance in Table 5, is as follows.
• The first topological loss function can be considered to specify a weaker form of prior topological information than in our original embedding, where we restricted to the term −(d1− b1). Now, the topological loss function states that the sum of persistence of the two
Figure 19 shows the topologically regularized PCA embeddings for these topological loss functions in order. While we observe small changes between the embeddings, which somewhat do agree with their respective designed topological loss function, in all cases, the embeddings do not strongly deviate from the ground truth topological model, which is one (and only one) prominent circle.
Different Loss Functions for different Topological Priors (example for UMAP) Note that in case of an orthogonal projection, which we encouraged in Figure 19 as we did in the main paper, it may be difficult to provide an embedding that well captures nontrivial topological models such as in Figure 18, when they are not truthfully present. However, non-linear dimensionality reduction methods such as UMAP may be more prone to model wrong prior topological information, as they directly optimize for the data point coordinates in the embedding space. To explore this, we considered an experiment where we embed our real bifurcating single cell data for a circular prior through the topological loss function −(d1−b1) evaluated on the 1st-dimensional persistence diagram. Thus, we topologically regularized the UMAP embedding of this data set for the exact same topological loss function as we originally used for our synthetic data lying on a circle (see also Figure 4 in the main paper). The result of topological regularization, using the exact same hyperparameters that we originally used for this data set as summarized in Table 1 in the main paper, is shown in Figure 20a.
We observe that when the embedding optimization is ran for 100 epochs (as was the embedding in the main paper), a cycle indeed starts to appear in the embedding. To further explore this, we increased the number of epochs to 1000, which showed to not increase the prominence of the cycle
much further (Figure 20a). However, by omitting the UMAP loss during the same topological optimization procedure of the embedding, we observe that a cycle becomes much more prominently present in the embedding (Figure 20b). Hence, inclusion of the structural UMAP embedding loss function through topological regularization prevents the cycle to be prominently and falsely represented in the entire data, and the cycle remains spurious in the embedding.
Naturally, higher topological regularization strengths λtop will more significantly impact the topological representation in the data embedding, so caution will need to be maintained when such information is absent. Figure 21 shows the topologically regularized UMAP embeddings for the circular prior for different topological regularization strengths λtop. Here, we again ran for 1000 epochs to explore long-term effects of topological regularization during the optimization.
When multiplying our original regularization strength λtop = 10 with a factor of 10, we see that the circular prior becomes much more prominently present in the topologically regularized UMAP embedding (λtop = 100). For lower topological regularization strengths, the topological prior struggles to become represented in the embedded (for λtop = 10 as in Figure 20a), or even appears to have no effect on the embedding at all (for λtop = 1). This can also be observed from Figure 22, where the embedding and topological loss functions are plotted according to the number of epochs and for the different topological regularization strengths λtop. We observe that for λtop ∈ {1, 10}, the topological loss functions eventually increase again and stagnate. For λtop = 100, the topological loss function appears to decrease indefinitely, at the cost of the embedding loss function which eventually increases again.
For all three cases, we also see that the incorrectly specified topological prior is not represented naturally in the data (Figure 21). For example, for λtop = 100, the cycle is prominently present, but the majority of points remain clustered together and neglect the circular shape that we aim to model through our topological loss functions. This artefact of topological optimization as well as the evolution of the loss functions (which are more generally applicable for higher dimensional embeddings) may provide some visual feedback to the user that the bias to the topological prior may be too weak or too strong. How these observations can be quantified, as to validate the passed topological prior and tune hyperparameters such as topological the regularization strength in an automated manner, is open to further research.

Review summaries:
This paper proposes a topological regularization method for incorporating topological prior knowledge for shaping data embeddings. This is achieved by introducing a new family of loss functions based on the characteristics of persistence diagrams. The results of the empirical evaluation suggest that the embeddings produced with the proposed method better capture the topological aspects of input data, by taking into account the prior topological knowledge.

The paper argues for including prior topological information about the structure between data points in some lower dimensional embedding space with the intention of improving the quality of the embeddings produced by tasks such as dimensionality reduction. The paper accomplishes this by adding a novel topological regularisation term to the standard embedding loss. The topological regularisation term incorporates the persistent homology of filtrations of the embedded data points. How this information is manipulated within the term encodes the prior knowledge the investigator wishes to use, such as the embedding should comprise a single cycle. The paper supports its method an array of qualitative and quantitative analyses on synthetic and real datasets. 

The paper develops a way to construct low-dimensional embeddings that combine an embedding cost (e.g. PCA, t-SNE, UMAP, etc.) with the topological cost, where the latter forces the embedding to have a particular topological structure as represented by the persistent homology. As an example, one can construct a data embedding and force it to be a circle, or a figure-8, or have three clusters.  The paper is based on a series of several 2020-2021 papers, with which I am unfamiliar, so it was difficult for me to judge on the novelty. I found the paper interesting, and the method seems to work well. My biggest criticism is that the goal of the method appears somewhat artificial: low-dimensional embeddings are typically done for data exploration, but one does not want to enforce any particular structure for the purposes of data exploration, and the topological structure is typically not known a priori. That's why I am giving a borderline score.