Stable Volumes for Persistent Homology

Figure 1 can be used to intuitively explain how a PD is determined from data. The input data in the example constitute a pointcloud with five points, as shown in Fig. 1(a). The points appear as a regular triangle and a square. Let $a$ be the length of edges of the regular triangle and the square. Since the five points themselves do not carry any interesting topological information, we construct a topological structure by putting circles with values of radii $r$ as indicated in Fig. 1(b)(c)(d)(e).

We now face the problem of how to determine the proper size of the circles. If the circles are too small, as in Fig. 1(b), the topology is simply equivalent to the five points. On the other hand, if the circles are too large, as in Fig. 1(e), the shape becomes acyclic, which makes it impossible to uncover any interesting topological information. PH solves this problem by considering the appearance and disappearance of homology generators associated with radii changes; that is, PH considers the increasing sequence from (b) to (e).

In the Fig. 1 example, two holes (homology generators in the 1st homology) appear at (c), one of which disappears at (d). The other hole disappears at (e). The pair of radii²²2The squares of radii are also often used in the literature of PH. of the appearance and disappearance of each homology generator is called a birth-death pair, and the multiset of all birth-death pairs is called a persistence diagram. In this example, the PD for the 1st homology is $\{(a/2,a/\sqrt{3}),(a/2,a/\sqrt{2})\}$. The persistence diagram is visualized by a scatter plot or a 2D histogram (Fig. 1(f)). The two pairs correspond to a regular triangle and a square. We can extract the geometric information of the pointcloud in Fig. 1(a) using the PD.

It is beneficial in data analysis using PH to detect the original homological structures corresponding to each birth-death pair. This is sometimes called “inverse analysis on PH”. Some applications of PH[18, 19] already use inverse analysis on PH. However, detection is not an easy task, as the representative cycle corresponding to a homology generator is not unique.

To solve this problem, various methods involving the solution of an optimization problem in homology algebra has been proposed[9, 17, 12, 9, 16, 11, 24, 32, 29] in various settings. For PH, optimal cycles[17], volume-optimal cycles[24], and persistent 1-cycles[12] have been proposed. These studies give formalizations suitable to PH. The “tightest” or “minimal” cycle is considered the best to interpret, and solving the optimization problem in homology algebra produces the tightest cycle. As an example, consider the simplicial complex in Fig. 2. The 1st homology group of the simplicial complex $H_{1}$ is isomorphic to $\mathbb{Z}$. $z_{1}$ and $z_{2}$ give the same generator of $H_{1}$, but $z_{1}$ is superior since the $z_{1}$ surrounds the hole more tightly than $z_{2}$. The above methods find $z_{1}$ by using the prescribed optimization technique.

It should be noted that there are two choices of functions to be optimized. One is the size of the cycle; the other is the internal volume of the cycle. In Fig. 2, $z_{1}$ is tighter than $z_{2}$ in both senses; however, it is possible for the solutions to the two optimization problems to differ. The better choice depends on the problem. This paper mainly uses the internal volume as an optimization function since it is suitable to define stable volumes in Section 3.

The existing methods have the following two problems:

• •

  The methods do not always give the minimal building blocks of the data

• •

  The result is unstable to noise

The example in Fig. 1 can be used to demonstrate these problems.

The pointcloud in Fig. 1 has two minimal building blocks, a triangle and a square. However, the optimization technique sometimes fails to find these minimal building blocks when a small noise is added. Figure 3(a) shows the pointcloud in Fig. 1(a) when a small noise is added. This figure also shows circles, and, as can be seen, a pentagon appears before either a triangle or a square appears. By applying the homology optimization technique to the data in Fig. 3(a), we produce a triangle and a pentagon, as shown in Fig. 3(b). On the other hand, consider the pointcloud in Fig. 3(c). The pointcloud here is also very close to that in Fig. 1(a); however, while a square appears first, the homology optimization technique gives a square and a triangle, as shown in Fig. 3(d).

The example demonstrates the following problems: (1) A small noise can change the result, and (2) The optimization technique sometimes fails to give minimal building blocks of the data. The first problem is related to the reliability of the analysis; the second problem is related to the interpretability of the result.

Previous studies [10, 8, 2, 22] have proved the stability theorem for persistence diagrams. The literature indicates that a PD is continuously changed by a small noise in the input data if we consider reasonable metrics. The stability of other PH outputs has also been studied in the context of machine learning and PH[4, 21, 1] or pointcloud summary method[20, 30]. Although these types of stability play an important role in the study of PH, such a stability theorem does not hold for the solution of homology optimization problems.

Both in Fig. 3(a) and Fig. 3(c), the minimal building blocks are a triangle and a square, and it is more natural that an inverse analysis technique gives a triangle and a square for both data. Therefore, we would hope to detect such structures even if a small noise is added.

These problems are practically significant, especially when we analyze crystalline structures using PH. We can demonstrate the problems using the synthetic crystalline data in Fig. 4. Here, the pointcloud consists of 27 points, arranged in a 3x3x3 cubical lattice (Fig. 4(a)). The distance between two vertices on the cube is $1$. A small noise is added to the pointcloud. Figure 4(b) shows the 1st PD of the pointcloud. As can be seen here, some birth-death pairs are concentrated around $(1/4,1/\sqrt{2})\approx(0.5,0.7)$. These pairs correspond to 1x1 squares in the lattice. By applying volume-optimal cycles[24], we found some loops corresponding to the pairs in Fig. 4(c) and (d). Some of the cycles shown in (c) are squares, as expected; however, others in (d) are larger structures resembling chairs. These larger structures were detected with the same mechanism in Fig. 3.

The purpose of this paper is to present a new method to solve these problems. That is, we propose a method to find minimal building blocks that is robust to noise.

Bendich, Bubenik, and Wagner [3] proposed a statistical approach for the problem, which can be outlined as follows:

1. 1.

   Computing a PD from the input data and choosing a birth-death pair to analyze

2. 2.

   Adding a small noise to the input data, computing a PD, and applying inverse analysis

3. 3.

   Repeating (2) multiple times

4. 4.

   Computing an average of the results of the inverse analysis

We can flexibly change the method of inverse analysis and the way in which the average of the results is computed. We can combine the statistical approach with optimal or volume-optimal cycles are as follows:

By applying this method to the birth-death pair $(1/2,1/\sqrt{2})$ in Fig. 1, we produce the result shown in Fig. 5. The result indicates that the four points on the square are robust to noise and that the leftmost point is less robust. This result is consistent with the fact that the pair $(1/2,1/\sqrt{2})$ corresponds to a square.

The stability of the method in a probabilistic sense was also proved by Bendich, Bubenik, and Wagner. Notably, the method provides a more reliable inverse analysis. The idea is quite clever, simple, and easy to implement. However, the method has a high computation cost since it requires the user to compute PDs and optimal cycles a large number of times. In [3], the authors repeat the computation of the generators 1,000 or 10,000 times. We suspect that fewer repetitions may be sufficient, but ten or a hundred trials will be necessary. Multiple trials and errors are typically needed to tune the noise bandwidth in order to apply the method, and thus the cost is not ignorable.

The “reconstructed shortest cycles” functionality of Ripserer.jl[34]³³3https://github.com/mtsch/Ripserer.jl by Čufar gives another solution of the problem. The functionality reconstructs the tightest 1-cycle using the shortest path algorithm and the representative of persistent cohomology. The functionality accepts a noise bandwidth parameter and computes a tighter loop.

The advantage of reconstructed shortest cycles is its efficiency. Since it uses the shortest path algorithm, the computational complexity is small. However, reconstructed shortest cycles can be applied only to 1st PH since it uses the shortest path algorithm. Another disadvantage is the lack of mathematical justification of the functionality. Now⁴⁴4Dec. 4, 2021 the functionality is declared as experimental and gives no mathematical documentation. We explain the algorithm in Appendix A.

In this paper, we propose stable volumes and a variant of stable volumes, called stable sub-volumes. The proposed method is based on the volume-optimal cycles and optimal volumes shown in [24]. Stable volumes produce minimal building blocks with a lower computation cost than the statistical approach. Below, we provide an outline of stable volumes. The exact formalization will be introduced in Section 2 and further developed in subsequent sections.

Let $X$ be a simplicial complex and $X^{(k)}$ the set of all $k$-simplices of $X$. We define a level function as follows.

Here, for simplicity, we assume that a level function $\hat{r}$ can distinguish every simplex; that is, we assume the following:

From the definition of the level function, $X_{\hat{r},t}=\{\sigma\mid\hat{r}(\sigma)<t\}$ is a subcomplex of $X$ and $t<t^{\prime}$ implies $X_{\hat{r},t}\subseteq X_{\hat{r},t^{\prime}}$. This means that $\{X_{\hat{r},t}\}_{t\in\mathbb{R}}$ is a filtration of simplicial complexes. From the filtration, the $k$th PD, $\mathrm{PD}_{k}(\{X_{\hat{r},t}\})$, is defined. Any birth-death pair in the diagram can be written as $(b,d)=(\hat{r}(\tau_{b}),\hat{r}(\omega_{d}))$, where $\tau_{b}$ is a $k$-simplex and $\omega_{d}$ is a $(k+1)$-simplex due to the theory of PH.

The optimal volume of the pair $(b,d)$ is defined by the solution to the following minimization problem:

where $C_{k+1}(X;\Bbbk)$ is a $(k+1)$th chain complex whose coefficient field is $\Bbbk$, $\mathcal{F}_{k}=\{\sigma\in X^{(k)}\mid\hat{r}(\tau_{d})<\hat{r}(\sigma)<\hat{r}(\sigma_{d})\}$, $\tau^{*}$ is an element of the dual basis of the cochain complex $C^{k}(X;\Bbbk)$, and $\|z\|_{0}=\{\omega\mid\alpha_{\omega}\not=0\}$ is the $\ell^{0}$ norm of $z$⁵⁵5The $\ell^{0}$ norm is, in fact, a norm since the triangle inequality does not hold. However, this is often called the $\ell^{0}$ norm in the context of machine learning and sparse modeling.. The optimal volume minimizes the internal volume surrounded by the loop. The optimal volume is defined on $\Bbbk=\mathbb{Z}/2\mathbb{Z}$; however, in practice, optimal volumes are determined by using linear programming with a change of coefficient field from $\mathbb{Z}/2\mathbb{Z}$ to $\mathbb{R}$ and an approximation of the optimization of the $\ell^{0}$ norm by the $\ell^{1}$ norm. For optimal volume $z$, $\partial z$ is called a volume-optimal cycle. Reference [24] shows that optimal volumes and volume-optimal cycles have good mathematical properties and provide good representations of a birth-death pair on a PD.

The stable volume for the pair $(b,d)$ with a noise bandwidth parameter $\epsilon$ is defined by the solution to the following minimization problem:

where $\mathcal{F}_{\epsilon,k}=\{\sigma\in X^{(k)}\mid\hat{r}(\tau_{d})+\epsilon\leq\hat{r}(\sigma)<\hat{r}(\omega_{d})\}$. This is quite similar to the definition of optimal volumes. The following theorem states that the stable volume is considered to be the “robust part” against noise.

We can regard $\mathcal{R}_{\epsilon}$ as the set of all possible level functions with small noise $\epsilon/2$. Then the theorem states that the stable volume is the common area of all possible optimal volumes under small noise. This theorem does not hold for the PH of another dimension, but it suggests that stable volumes give a better result than optimal volumes.

For practical applications, assumption (1) is overly strict. Vietoris-Rips, Čech, and alpha filtrations do not satisfy the condition. Therefore, we introduce an order with levels in Section 2.2 and formalize all of them using the order with levels.

As the stable volume method has already been implemented in HomCloud [25], stable volumes can now be used to analyze the data.

We can demonstrate stable volumes by using the same input data as in Fig. 4. Figure 6 shows the optimal volumes (Fig. 6(a)) and stable volumes (Fig. 6(b)) of the same two birth-death pairs. The birth-death pairs correspond to 1x1 squares in the lattice points, but Fig. 6(a) shows larger loops than expected. In contrast Fig. 6(b) shows the expected 1x1 squares. This example suggests that stable volumes give better results and help promote a better intuitive understanding of PDs.

Figure 7 shows another example of stable volume. Here, we consider 2D lattice points with defects (Fig. 7 (a)). The configuration of the points is somewhat perturbed from the complete lattice, and some points are removed randomly. Figure 7(b) shows the 1st PD of the points, and Fig. 7(c) shows the optimal volume of (0.496, 4.371). Figure 7(d) shows the stable volume of the same birth-death pair. The optimal volume in Fig. 7(d) seems to surround the holes in the pointcloud more naturally than is the case in Fig. 7(c).

Additional demonstrations of stable volumes using synthetic and real data are given in Section 6.

The remainder of the paper is organized as follows: Section 2 introduces the concept of PH, optimal volumes, and volume-optimal cycles. Section 3 defines stable volumes for the $(n-1)$th PH embedded in $\mathbb{R}^{n}$. In this section we also discuss the limitation of stable volumes. Section 4 shows an alternative formalization of stable volumes using mathematical optimization. Section 4.1 provides a proof that the two formalizations are consistent for the $(n-1)$th PH when the simplicial complex is embedded in $\mathbb{R}^{n}$. Section 4.2 extends the formalization to other cases. Section 5 discusses the software implementation of stable volumes. Section 6 gives examples using synthetic and real data. This section also provies the comparison with previous methods, the statistical method and the reconstructed shortest cycles. Section 7 discusses the way in which the noise bandwidth parameter can be tuned. Section 8 compares stable volumes and sub-volumes. Section 9 summarizes the paper and offers concluding remarks. In this section we also discuss the differences between the proposed method and the methods previous works.

PH is defined on a filtration of topological spaces. Let $X$ be a finite simplicial complex⁶⁶6To simplify the explanation in this paper, we use a simplicial complex. However, the content of the paper can be easily extended to cubical or cell complexes. and

be a filtration of subcomplexes of $X$. The $k$th persistent homology $H_{k}(\mathbb{X};\Bbbk)$ with a coefficient field $\Bbbk$ is defined as

where $H_{k}(X_{i};\Bbbk)$ is the $k$th homology $\Bbbk$-vector space of $X_{i}$ and $\to$ is the linear map induced by the inclusion maps. The fundamental theorem of PH is as follows.

On an intuitive level, the map $0\to\Bbbk$ indicates the appearance of a homology generator, $\Bbbk\to\Bbbk$ indicates the persistence of the generator, and $\Bbbk\to 0$ indicates the disappearance of the generator.

We note that we should use a field as a coefficient ring of the homology module since the theorem holds only when $\Bbbk$ is a field.

We earlier introduced the idea of a level function. We now introduce the concept of total order consistent with the level function. Total order is used to rigorously describe the algorithms; the level function is needed to rigorously describe Theorem 2. Simplices with the same level often appear in an alpha filtration or a Vietoris-Rips filtration; thus, the concept of total order is required to deal with such filtrations. The pair consisting of the total order and the level function is called an order with level.

Let $\prec$ be a total order on the set of all simplices of $X$. We assume that

We can number all simplices in $X$ by this ordering as follows:

From the condition (14), $X_{i}=\{\sigma_{1},\ldots,\sigma_{i}\}$ is always a subcomplex of $X$ and

is a filtration of a simplicial complex. In this filtration, the number of simplices is increased one by one, which simplifies the mathematical description.

To include additional information in the filtration, we consider an order with level.

From the total order $\preceq_{r}$, a filtration $\mathbb{X}_{r}$ is defined by (15) and (16). Theorem A gives the unique decomposition $H_{k}(\mathbb{X}_{r};\Bbbk)\simeq\bigoplus_{m=1}^{L}I(b_{m},d_{m})$. For each $(b_{m},d_{m})$, the simplex $\sigma_{b_{m}}$ is called a birth simplex, $\sigma_{d_{m}}$ is called a death simplex, and the pair $(\sigma_{b_{m}},\sigma_{d_{m}})$ is called a birth-death simplices pair. When $d_{m}=\infty$, we write the pair as $(\sigma_{b_{m}},\star)$ using the special symbol $\star$. $D_{k}(\mathbb{X}_{r})$ denotes the set of all birth-death simplices pairs. Moreover, $\hat{r}(\sigma_{b_{m}})$ is called a birth time, $\hat{r}(\sigma_{d_{m}})$ is called a death time, and the pair of the birth time and death time, $(\hat{r}(\sigma_{b_{m}}),\hat{r}(\sigma_{d_{m}}))$, is called a birth-death pair. When $\sigma_{d_{m}}=\star$, $\hat{r}(\star)$ is defined as $+\infty$. The multiset of birth-death pairs is called a persistence diagram:

We write the persistence diagram as $\mathrm{PD}_{k}(X,r)$.

We can readily show the following fact.

By using the level information, we can describe the stability theorem of persistence diagrams [10, 8, 2, 22]. Let $\mathrm{PD}_{1},\mathrm{PD}_{2}$ be two persistence diagrams. The bottleneck distance between two diagrams $d_{B}(\mathrm{PD}_{1},\mathrm{PD}_{2})$ is defined as

where $\gamma$ ranges over all bijections from $\mathrm{PD}_{1}\cup\Delta$ to $\mathrm{PD}_{2}\cup\Delta$, $\Delta$ is the diagonal, and $\|\cdot\|_{\infty}$ is the $\ell^{\infty}$ norm. In our setting, the stability theorem of PH is described as below.

The theorem ensures the robustness of PDs to noise.

Optimal volume is proposed in [24] as a way to extract homological structures corresponding to a birth-death pair. For an order with levels $r=(\hat{r},\preceq_{r})$ and a birth-death simplices pair $(\tau_{0},\omega_{0})\in D_{k}(\mathbb{X}_{r})$ with $\omega_{0}\not=\star$, the optimal volume for the pair is formalized as the solution to the following minimization problem:

where $C_{k+1}(X;\Bbbk)$ is a chain complex whose coefficient field is $\Bbbk$, $\mathcal{F}_{k}=\{\sigma\in X^{(k)}\mid\tau_{0}\prec_{r}\sigma\prec_{r}\omega_{0}\}$, $\tau^{*}$ is an element of the dual basis of cochain complex $C^{k}(X;\Bbbk)$, and $\|z\|_{0}=\{\omega\mid\alpha_{\omega}\not=0\}$ is the $\ell^{0}$ norm of $z$. For solution $z$, $\partial z$ is called a volume-optimal cycle.

Reference [24] shows the following fact, which indicates that the volume-optimal cycle is suitable for representing a birth-death pair.

In this section, we consider a triangulation of a convex set in $\mathbb{R}^{n}$ and its $(n-1)$th PH with $n\geq 2$. More precisely, we assume the following:

Schweinhart [29] pointed out that the $(n-1)$th PH was isomorphic to the 0th persistent cohomology of the dual filtration by Alexander duality. Zeroth cohomology is deeply related to the connected components in the dual filtration. This gives rise to another formalization of $(n-1)$th PH. Reference [24] generalizes the idea shown in [29].

We now examine an order with levels $r$ on $X$. Consider the filtration $\mathbb{X}_{r}$ in (15) and (16) given by $\preceq_{r}$. For simplicity, we will use $\mathbb{Z}_{2}$ as the coefficient field. Then the following proposition and theorems hold [29, 24].

From Theorem E, we know that $D_{n-1}(\mathbb{X}_{r})$ can be regarded as a forest (i.e., the union of distinct trees) by the inclusion relation. In [29], the forest is called persistence trees. Moreover, we can compute the persistence trees by the merge-tree algorithm (Algorithm 2).

To describe the algorithm, consider the one-point compactification space $\mathbb{R}^{n}\cup\{\infty\}\simeq S^{n}$. The following facts are well known.

We can extend the order with levels $r$ onto $X\cup\{\omega_{\infty}\}$ by regarding $\omega_{\infty}$ as the maximum element and $\hat{r}(\omega_{\infty})=+\infty$. To describe the algorithm, we consider a directed graph $\bar{G}$ whose nodes are $n$-cells in $X\cup\{\omega_{\infty}\}$. An edge has extra data in $X^{(n-1)}$, and we can write the edge from $\omega$ to $\omega^{\prime}$ with extra data $\tau$ as $(\omega\xrightarrow{\tau}\omega^{\prime})$. The directed graph $\bar{G}$ is increasingly updated in Algorithm 2.

Since the graph is always a forest throughout the algorithm[24], we can find a root of a tree that contains an $n$-cell $\omega$ in the graph $\bar{G}$ by recursively following the edges from $\omega$. We call this procedure Root($\omega,\bar{G}$).

The following theorem gives the interpretation of the result of the algorithm as persistence information.

We can prove the above theorems using Alexander duality. Appendix A in [23]⁸⁸8 This paper is the preprint version of [24]. The contents of Appendix A in [23] are omitted from [24]; thus, we sometimes refer to the preprint version. provides a detailed discussion.

To show the theorem, we introduce a dual graph of $X$ and its subgraphs.

We call the graph $G_{X}$ the dual graph of $X$. We define the subgraphs of $G_{X}$, $G_{X}(\succeq_{r}\sigma)$ and $G_{X}(\hat{r}\geq s)$, as

where $\sigma$ is a simplex in $X$ and $s$ is a real number. The condition (14) ensures that $G_{X}(\succeq_{r}\sigma)$ and $G_{X}(\hat{r}\geq s)$ are subgraphs of $G_{X}$. We also define $G_{X}(\succ_{r}\sigma)$ by replacing $\succeq_{r}$ with $\succ_{r}$.

We introduce the notation $\bar{G}_{\sigma}$ as $\bar{G}$ in Algorithm 2 when the inside of the (LOOP) is finished at $\sigma$. The following propositions are essential for the proof of Theorem 11. The facts are shown as Fact 17 and Fact 18 in [23].

The next proposition describes the role of the birth and death simplices. $K(G,\omega)$ denotes the connected component of a graph $G$ which contains a vertex $\omega$. $K_{\mathrm{V}}(G,\omega)$ and $K_{\mathrm{E}}(G,\omega)$ denote all the vertices and all the edges of $K(G,\omega)$.

We can easily prove the proposition by considering Algorithm 2 and Facts 15 and 16.

The following equality holds from the definition of stable volume (25), Theorem F, Facts 15 and 16, and Proposition 17.

The following equality also holds from Proposition 17.

Therefore, we can show the following relationship to prove Claim 12.

The following lemma is essential to prove Claim 12.