Topology-Preserving Dimensionality Reduction via Interleaving Optimization

Optimization of persistent homology-based objective functions has attracted much recent attention [29, 17, 19, 22, 6, 5, 26] with a particular focus on applications in computational geometry and deep learning. Of particular relevance is the work of [27] which uses a persistence-based objective to preserve critical edges in the persistent homology of the Vietoris-Rips filtration in a learned latent space of an autoencoder.

The idea of using persistent homology to compare different dimension reduction schemes was initiated by Rieck and Leitte [30, 31], which uses the $2$-Wasserstein distance to compare the persistence diagrams of a data set before and after dimension reduction. Several non-differentiable methods incorporating persistent homology have been developed [12, 36, 13]. With the development of optimization techniques for persistent homology, several differentiable methods have been proposed based on optimization of the $2$-Wasserstein metric on persistence diagrams [20, 34] and an approach based on simulated annealing [37]. While not employed on Vietoris-Rips filtrations, the work of [29] uses the bottleneck distance for optimization of functional maps on shapes.

The interleaving/bottleneck distance has long been a key tool developed for the study of persistent homology under perturbation of the input [9]. The interleaving distance was first introduced by [7], and applied to the study of Vietoris-Rips filtrations in [8] to bound the distance of persistent homology by the Gromov-Hausdorff distance between the input point clouds. The idea of developing confidence regions for persistence pairs was developed by [15] in the context of sampling.

This work presents a novel approach to dimension reduction using optimization of the interleaving/bottleneck distance between the persistent homology of Vietoris-Rips Filtrations of an original data set $X$ and the data set $Y$ after dimensionality reduction.

1. 1.

   We show how the interleaving distance can be used to quantify a scale at which topological features in $X$ and features in $Y$ are in correspondence, and be used to select homological features of $Y$ in correspondence with features in $X$.

2. 2.

   We show how to incorporate the interleaving distance explicitly into the optimization of the embedding $Y$ and prove the existence of descent directions under mild conditions.

3. 3.

   We demonstrate this technique in finding optimal linear projections of the data set $X$ to preserve the bottleneck distance on several examples¹¹1Our implementations are made publicly available at https://github.com/CompTop/Interleaving-DR. with interesting topology.

We are interested in discovering and preserving topological features of a point cloud $X$ together with a notion of dissimilarity $d$, and refer to the combination of these two data as a dissimilarity space $(X,d)$. A dissimilarity is a function $d:X\times X\to\mathbb{R}_{\geq 0}$, with $d(x,x)=0$ for any $x\in X$. We will typically consider dissimilarities that are metrics (in particular which satisfy triangle inequality), but many of the bounds here hold more generality. Examples of topological features of the space $(X,d)$ include clusters formed through single-linkage clustering or “holes” forming loops in the $r$-nearest neighbors graph of $X$.

Vietoris-Rips filtrations are commonly used in conjunction with persistent homology to create features for finite dimensional metric spaces (point clouds) [2]. Given a dissimilarity space $(X,d)$, a Vietoris-Rips complex consists of simplices with a maximum pairwise dissimilarity between vertices is less than some threshold $r$ :

A Vietoris-Rips filtration is a nested sequence of Vietoris-Rips complexes $X_{r}\subseteq X_{s}$ if $r\leq s$.

Homology is a functor from the category of topological spaces and continuous maps to the category of vector spaces and linear maps (for a general introduction see [18]). The dimension of the $k$-dimensional homology vector space $H_{k}(X)$ of a topological space $X$ counts $k$-dimensional topological features of $X$: $\dim H_{0}(X)$ is the number of connected components, $\dim H_{1}(X)$ counts loops, and $\dim H_{k}$ generally counts $k$-dimensional voids. A continuous map $f:X\to Y$ has an induced map $H_{k}(f):H_{k}(X)\to H_{k}(Y)$ which maps vectors associated with topological features in $X$ to vectors associated with topological features in $Y$. The computation of $H_{k}(X)$ begins with the construction of a chain complex $C_{\ast}(X)=\{C_{k}(X),\partial_{k}:C_{k}(X)\to C_{k-1}(X)\}_{k\geq 0}$ where $C_{k}(X)$ is a vector space with a basis element for each $k$-simplex in $X$, and the boundary map $\partial_{k}$ sends each basis element to a linear combination of basis elements of faces in the boundary of the associated simplex. The boundary maps satisfy $\partial_{k}\circ\partial_{k+1}=0$, and $H_{k}(X)$ is the quotient vector space $\ker\partial_{k}/\operatorname{img}\partial_{k+1}$.

Persistent homology [14, 38] is an algebraic invariant of filtrations which captures how the topology of a filtration changes using homology. The output of persistent homology is a persistence vector space $V_{\ast}$, consisting of vector spaces $\{V_{r}=H_{k}(X_{r})\}_{r\in\mathbb{R}}$ and linear maps induced by inclusion $\{\iota^{V}_{r,s}:V_{r}\to V_{s}\}_{r\leq s\in\mathbb{R}}$ which satisfy a consistency condition $\iota^{V}_{r,t}=\iota^{V}_{r,s}\iota^{V}_{s,t}$ for all $r\leq s\leq t$.

Persistence vector spaces are classified up to isomorphism by birth-death pairs $\{(b_{i},d_{i})\}_{i\in I}$, or equivalently their persistence barcode [38] or interval indecomposables [3]. Each pair $(b,d)$ is associated to the appearance of a new homology vector at filtration parameter $b$ (meaning it is not in the image of an induced map), which maps through the persistence vector space until it enters the kernel of an induced map at filtration parameter $d$. The length of the pair $(b,d)$ is the difference $|d-b|$.

Every birth and death in persistent homology is associated with the addition of a particular simplex in the filtration. This follows from the definition of homology of the quotient vector space $H_{k}(X)=\ker\partial_{k}/\operatorname{img}\partial_{k+1}$. The addition of a new $k$-simplex increases the dimension of $C_{k}(X)$ by one, and either increases the dimension of $\ker\partial_{k}$ by one, causing a birth in $H_{k}(X)$, or increases the dimension of $\operatorname{img}\partial_{k}$ by one, causing a death in $H_{k-1}(X)$. Because the Vietoris-Rips filtration is determined by its edges, the filtration value of every simplex can be mapped to the largest pairwise distance. This provides a way to map the gradient of a function with respect to births and deaths to a gradient with respect to each pairwise distance – see [17] for additional details.

Interleavings allow for the comparison of two persistence vector spaces [7], as well as other objects filtered by some partially ordered set [1]. Let $V_{\ast}$ and $W_{\ast}$ be 1-parameter persistence vector spaces. An $\epsilon$-shift map $f_{\ast}:V_{\ast}\to W_{\ast}$ is a collection of maps $f_{r}:V_{r}\to W_{r+\epsilon}$ so that the following diagram commutes for all parameters $r$

An $\epsilon$-interleaving between $V_{\ast}$ and $W_{\ast}$ is a pair of $\epsilon$-shift maps $f_{\ast}:V_{\ast}\to W_{\ast}$ and $g_{\ast}:W_{\ast}\to V_{\ast}$ so that $g_{r+\epsilon}f_{r}=\iota^{V}_{r,r+2\epsilon}$ and $f_{r+\epsilon}g_{r}=\iota^{W}_{r,r+2\epsilon}$ for all parameters $r$.

The interleaving distance [7] on persistence modules $V_{\ast}$ and $W_{\ast}$ is

This notion of distance satisfies triangle inequality through the composition of interleavings. Note that the addition or removal an arbitrary number of zero-length pairs to a persistence vector space $V_{\ast}$ to obtain $V^{\prime}_{\ast}$ results in $d_{I}(V_{\ast},V^{\prime}_{\ast})=0$.

The construction of general interleavings, let alone those that would realize the interleaving distance, can be a daunting task. Fortunately, for 1-parameter persistent homology the interleaving distance $d_{I}$ is equivalent to the geometric (and easily computable) bottleneck distance $d_{B}$ on persistence diagrams [25].

The bottleneck distance considers the birth-death pairs $\{(b_{i},d_{i})\}_{i\in I}$ as points in the 2-dimensional plane. The persistence diagram $\operatorname{dgm}(V_{\ast})$ is the union of this discrete multi-set of the points $\{(b_{i},d_{i})\}$ with the diagonal $\Delta=\{(x,x)\mid x\in\mathbb{R}\}$ where points in $\Delta$ are counted with infinite multiplicity.

The matching in the bottleneck distance actually gives an interleaving which maps a vector associated to a persistence pair in $\operatorname{dgm}(V_{\ast})$ to the vector associated with the matched pair in $\operatorname{dgm}(W_{\ast})$ which realizes the interleaving distance.

In practice, persistent homology of the Vietoris-Rips filtration can be quickly approximated using sub-sampling. Bounds on this approximation come from the Hausdorff or Gromov-Hausdorff distance on between point clouds [8].

The theorem above can also provide us an interleaving/bottleneck distance bound for sub-sampled data sets only if we can find their Gromov-Hausdorff distance.

The proof only needs the triangle inequality of metrics.

Let $H_{k}(X;r)=H_{k}(\mathcal{R}(X,d_{X};r))$ denote the $k$-dimensional homology of the Vietoris-Rips complex at parameter $r$, and $H_{k}(X)$ denote the $k$-dimensional persistent homology of the Vietoris-Rips filtration. Similarly, we have $H_{k}(Y;r)$ and $H_{k}(Y)$. We refer to each persistence pair in $H_{k}(X)$ or $H_{k}(Y)$ as a homological feature of $(X,d_{X})$ or $(Y,d_{Y})$ respectively. We would like to select features of $(Y,d_{Y})$ which are in correspondence with features of $(X,d_{X})$. A simple selection procedure is to compute the interleaving distance $\epsilon=d_{\mathrm{I}}(H_{k}(\mathcal{R}(Y,d_{Y}),H_{k}(\mathcal{R}(X,d_{X}))$, and to select any homology class $(b,d)\in H_{k}(\mathcal{R}(Y,d_{Y}))$ with $|d-b|>2\epsilon$.

As a result, not only does this selection procedure guarantee that we will make no topological type-I errors, but we also will not make any type-II errors involving a homological feature of $(X,d_{X})$ of sufficient size. Note that if we lower our threshold for selecting pairs of $H_{k}(Y)$ then we could introduce the possibility of homological type-I errors, and that there is always the possibility of homological type-II errors when neither pair in the correspondence has sufficient length.

This selection procedure can be used to assess how well any transformation of point cloud data $X$ preserves topology. In particular, we can compute the bottleneck distance $\epsilon=d_{\mathrm{I}}(H_{k}(X),H_{k}(Y))$, the number of features of $H_{k}(Y)$ with $|d-b|>2\epsilon$ and the number of features of $H_{k}(X)$ with $|d-b|>4\epsilon$. In the context of dimension reduction, it would be desirable to minimize the interleaving distance in order to maximize the number of features we can identify which are in correspondence with features in the original data set. This is the approach we pursue in Section 4.

In table Table 1 we compare several algorithms for dimension reduction on a set of images from the Columbia Object Image Library (COIL-100) [28] taken of a tomato at various angles in a circle which has a single large $H_{1}$ homological feature in the original data. Methods compared include PCA, MDS [23], and ISOMAP [32] with a method based on minimizing the bottleneck distance developed in Section 4, PH, and a hybrid PH + PCA. Because every Vietoris-Rips filtration has a single $H_{0}$ pair with death at $\infty$, at least one $H_{0}$ feature will always be selected. Only PH + PCA allows for the selection of an $H_{1}$ feature. Visualization of each embedding can be found in the appendix, and a more detailed visualization of the PH + PCA embedding can be found in Figure 1.

Some care must be taken in interpreting the interleaving as presented here. Importantly, the correspondence between persistence pairs of $H_{k}(X)$ and $H_{k}(Y)$ is entirely algebraic; there are not necessarily topological shift maps $f:\mathcal{R}(X;r)\to\mathcal{R}(Y;r+\epsilon)$ and $g:\mathcal{R}(Y;r)\to\mathcal{R}(X;r+\epsilon)$ which induce the interleaving on homology. Even if the interleaving distance between $X$ and $Y$ is zero, there is no guarantee that there is a natural topological map between the two spaces. There are many possible spaces which have identical persistent homology [10], and to maintain some level of geometric interpretability of persistence pairs of $Y$ in terms of the persistence pairs of $X$, it is desirable to incorporate additional constraints onto the embedding $Y$, as is often the case in dimension reduction algorithms.

Gradient-based optimization techniques can be applied to persistent homology by backpropagating the gradient of a function of the persistence pairs back to the input values of a filtration. This is often done by considering a featurization of the persistence pairs such as algebraic functions of the pairs [17] or persistence landscapes [22, 6], but in our situation, we will use the bottleneck distance $\mathrm{d}_{\mathrm{B}}(\operatorname{dgm}(X),\operatorname{dgm}(Y))$.

Optimization with Vietoris-Rips filtrations is described in detail in [17], which we summarize here. The key is that every simplex addition in the filtration either creates or destroys homology and the corresponding birth or death takes that filtration value. If $f$ is a function of persistence pairs, this allows for the mapping of $\partial f/\partial b$ or $\partial f/\partial d$ to $\partial f/\partial w_{\sigma}$ where $w_{\sigma}$ is the filtration value of the unique simplex $\sigma$. In the case of Vietoris-Rips filtrations, the filtration value of a simplex $(x_{0},\dots,x_{k})$ is the maximum pairwise distance $d_{Y}(x_{i},x_{j})$ where $x_{i},x_{j}$ are vertices in the simplex, so this can then be backpropagated to a gradient $\partial f/\partial d_{Y}(x_{i},x_{j})$. There is a potential issue here which is that a single edge may map to multiple higher-order simplices, but Theorem 4.2 indicates that this is not generally a problem. Finally, if we choose a differentiable metric on $Y$ such as the Euclidean metric, we can backpropagate the gradient to point locations in the embedding.

We are interested in optimizing the embedding $Y$ to minimize the interleaving distance via the bottleneck distance $\mathrm{d}_{\mathrm{B}}(\operatorname{dgm}(Y),\operatorname{dgm}(X))$. Because the original data set $X$ is fixed, we have a function $f(\operatorname{dgm}(Y))=\mathrm{d}_{\mathrm{B}}(\operatorname{dgm}(Y),\operatorname{dgm}(X))$ which we can fit into the optimization framework for persistent homology.

First, we recall that the bottleneck distance uses a matching between persistence pairs of $Y$ and $X$ and the diagonal $\Delta$ representing potential zero-length pairs, and that the distance is computed from the maximum-weight matching. This leads to three possibilities.

1. 1.

   The maximum weight matching occurs between two non-diagonal pairs.

2. 2.

   The maximum weight matching occurs between a pair in $\operatorname{dgm}(Y)$ and the diagonal in $\operatorname{dgm}(X)$.

3. 3.

   The maximum weight matching occurs between the diagonal in $\operatorname{dgm}(Y)$ and and a pair in $\operatorname{dgm}(X)$.

The bottleneck distance does not consider matchings between two diagonal points. In the first two cases, it is possible to find a descent direction. In the third case, $Y$ is at a local saddle point of the bottleneck distance.

Proposition 4.3 implies that we have great freedom to perturb the embedding $Y$ while minimizing the bottleneck distance. In the third case above, when the diagonal of $\operatorname{dgm}(Y)$ is involved in the maximum weight matching, we have the freedom to perturb $Y$ in any direction. This allows for easy optimization of the bottleneck distance in conjunction with a secondary objective.

We generate a data set, 5-Tendril consisting of 500 points in 5 dimensions sampled along 5 tendrils, each of which consists of 50 randomly generated points on a canonical basis vector $e_{i}$ of Euclidean space. Here, our PH optimization method will first sample 100 points and then optimize the bottleneck distance on H0 persistent diagram. We perform 3 different dimension reduction algorithms: PCA, MDS and PH optimization. In Figure 3, the results show that PCA and MDS can only see 3 branches, while PH optimization can see 5.

We generate a data set of 500 samples in 5 dimensions with $\binom{5}{2}=10$ cycles, each of which consists of 50 points and lies in a 2-dimensional plane spanned by two canonical bases $e_{i}$ and $e_{j}$ with center $e_{i}$ + $e_{j}$ and radius one. For speedup, our PH optimization method samples 200 points and then optimizes bottleneck distance on $H_{1}$.

Since the dataset consists of cycles in orthogonal planes and our PH optimization method will also purse an orthonormal projection, there is no way to find a projection that can show and divide all cycles apart within the orthogonality constraint. In Figure 4, we show the dimension reduction results with 4 different methods: PCA will lead to a tangle where cycles cross with each other, and if without labels, one cannot tell how many cycles exist in this plot; ISOMAP will provide us a five-star with half of orthogonal cycles missed, but without labels, it is hard to determine if there are cycles exists; MDS can provide a 5-ring, but some of are not comprised of a single orthogonal cycle, for example, the yellow orthogonal cycle lies on two rings; our PH method can provide 3 cycles and the three straight lines between can indicate the collapse of orthogonal structure in the projection.

The Columbia Object Image Library (COIL-100) [28] dataset contains 7200 colorful images of 100 objects, where each object has 72 $128\times 128$ images with 3 color channels taken at pose intervals of 5 degrees. We pick the tomato dataset (see Figure 2) with shape $72\times 49152(=128*128*3)$ from COIL-100 and perform our PH optimization algorithm on both $H_{1}$ and $H_{0}$ persistent diagrams. Although the number of points is only 72, the high feature dimension makes the direct running of our procedure on Pytorch prohibitive due to the memory issue of Cayley transform. To solve the problem, we devise two procedures: a) Optimizing the projection directly using the Cayley transform to handle orthogonality; b) first we perform PCA on the tomato dataset to reduce dimension to 10 and then keep reducing the dimension to 2 by our Pytorch PH optimization algorithm, denoted by PH+PCA.

We include visualization of the results using methods PH, PH+PCA, PCA, ISOMAP, MDS into the Appendix and show the result of PH+PCA in Figure 1. In Figure 1a), we can see that the 72 points in $\mathbb{R}^{2}$ constitute a great circle, which demonstrate that the dataset consists of pictures taken at 360 degrees at pose interval 5 around the tomato. In Figure 1b), we draw the persistence diagram of the projected data set and a unconfident band with width 9.065, which is the twice of bottleneck distance between the original and projected tomato dataset. By the discussion in Section 3, there is only one H1 class in the projected data that we can ensure also exists in the original one. In Figure 1c), we then draw the certain H1 representative in red and the longest uncertain H1 representative in blue.

Natural image patches are a well-studied data set with interesting topological structures at various densities [4]. We follow the data processing procedure of [24] to sample $3\times 3$ patches from the van Hateren natural images database [33]. We further refine a sub-sample of 50,000 patches using the co-density estimator of [11] with $k=5,p=40\%$ to obtain a data set of 20,000 patches which resembles the “three-circle” model of [11].

In Figure 5 we apply our procedure to two initial projections. In both cases, we use a greedy subsampling of 100 points in the data which contains 5 robust $H_{1}$ pairs, agreeing with the three-circle model. We first initialize the projection with the first two principal components of the data and then refine by optimizing the bottleneck distance on $H_{1}$ pairs. The initial projection onto principal components displays a clear “primary circle,” with the two secondary circles projecting onto two chords. Our procedure decreases the $H_{1}$ bottleneck distance on the 100 sampled points from $5.6$ to $4.7$ but there is minimal visual difference between the two projections, indicating that the initialization was near a local optimum.

The second row of Figure 5 starts with a random projection into two dimensions. We again optimize to minimize the bottleneck distance on the $H_{1}$ pairs, and decrease the bottleneck distance on the sampled points from $7.9$ to $3.7$. In this case, there is a noticeable visual difference between the two projections. In the first projection, a noisy projection of the primary circle is visible, and in the second projection, we see a clear visualization of one of the two secondary circles, with the primary circle and other secondary circle collapsed to a chord.

In both these experiments, our $H_{1}$ bottleneck distance bounds do not allow us to confidently select any features in the visualization. As with the orthogonal cycle data, the reason is fundamental: the three circle model embeds into a Klein bottle [4], and any projection to fewer than 4-dimensions (necessary to embed the Klein bottle) will result in spurious intersections of at least two of the three circles. Despite these limitations, our method is still able to present a subset of the important $H_{1}$ features of this data.