Barcodes of Towers and a Streaming Algorithm for Persistent Homology

TL;DR

This paper introduces a streaming algorithm to efficiently compute persistent homology from towers of simplicial complexes, producing a filtration with the same barcode while optimizing space and computational resources.

Contribution

It presents a novel streaming algorithm that constructs a filtration equivalent to a tower's barcode with minimal space complexity, combining coning strategies and matrix reduction techniques.

Findings

Algorithm is asymptotically space-efficient

Can compute barcode directly from towers in practice

Space complexity independent of tower length

Abstract

A tower is a sequence of simplicial complexes connected by simplicial maps. We show how to compute a filtration, a sequence of nested simplicial complexes, with the same persistent barcode as the tower. Our approach is based on the coning strategy by Dey et al. (SoCG 2014). We show that a variant of this approach yields a filtration that is asymptotically only marginally larger than the tower and can be efficiently computed by a streaming algorithm, both in theory and in practice. Furthermore, we show that our approach can be combined with a streaming algorithm to compute the barcode of the tower via matrix reduction. The space complexity of the algorithm does not depend on the length of the tower, but the maximal size of any subcomplex within the tower.

Equations

$$\underset{\text{b−1 times}}{\underbrace{0\stackrel{0}{\to }\dots \stackrel{0}{\to }0}}\stackrel{0}{\to }\underset{\text{d−b+1 times}}{\underbrace{\mathbb{F}\stackrel{\mathrm{id}}{\to }\dots \stackrel{\mathrm{id}}{\to }\mathbb{F}}}\stackrel{0}{\to }\underset{\text{m−d times}}{\underbrace{0\stackrel{0}{\to }\dots \stackrel{0}{\to }0}},$$

$${\mathbb{L}}^{\ast }=\mathbb{K}\cup \left(u\ast \overline{\mathrm{St}}(v,\mathbb{K})\right).$$

$${\hat{\mathbb{K}}}_{i+1}={\hat{\mathbb{K}}}_i\cup \left(v\ast \mathrm{Act}\overline{\mathrm{St}}(u,{\hat{\mathbb{K}}}_i)\right)$$

$${\hat{\mathbb{K}}}_{i+1}={\hat{\mathbb{K}}}_i\cup \left(u\ast \mathrm{Act}\overline{\mathrm{St}}(v,{\hat{\mathbb{K}}}_i)\right)$$

$$\begin{array}{ccccccc}H({\mathbb{K}}_0)& \stackrel{{\varphi }_0^{\ast }}{\to }& H({\mathbb{K}}_1)& \stackrel{{\varphi }_1^{\ast }}{\to }& ...& \stackrel{{\varphi }_{m-1}^{\ast }}{\to }& H({\mathbb{K}}_m)\\ \downarrow {\mathrm{inc}}_0^{\ast }& & \downarrow {\mathrm{inc}}_1^{\ast }& & \downarrow {\mathrm{inc}}_m^{\ast }& & \\ H({\hat{\mathbb{K}}}_0)& \stackrel{{\hat{\varphi }}_0^{\ast }}{\to }& H({\hat{\mathbb{K}}}_1)& \stackrel{{\hat{\varphi }}_1^{\ast }}{\to }& ...& \stackrel{{\hat{\varphi }}_{m-1}^{\ast }}{\to }& H({\hat{\mathbb{K}}}_m)\end{array}$$

$$|E(x)|$$

$$\sum _{x\in N}|E(x)|\le (\Delta +1)\cdot n,$$

$$c(x)\le \min \{|C_1|,|C_2|\}.$$

$$c(p_i)=c(p_{i-1})+c(x_i)\le 2\cdot (|E(x_{i-1})|+|E(x_{i-1}^{\prime })\setminus E(x_{i-1})|)\le 2\cdot |E(x_i)|,$$

$$r(x)=\{\begin{array}{ll}1,& \text{ if x is a leaf}\\ r(y_1)+1,& \text{ if x has children y1,y2 and r(y1)=r(y2)}\\ \max \{r(y_1),r(y_2)\},& \text{ if x has children y1,y2 and r(y1)=r(y2)}\end{array}$$

$$s(x)=1+s(y_1)+s(y_2)\ge 1+\left(2^{r(y_1)}-1\right)+\left(2^{r(y_2)}-1\right)=2^{r(x)}-1.$$

$$\beta (\hat{\mathbb{L}}):=\sum _{i=0}^{\Delta }{\beta }_i(\hat{\mathbb{L}})$$

Full text

1 Introduction

Motivation and problem statement.

Persistent homology [18, 7, 17] is a paradigm to analyze how topological properties of general data sets evolve across multiple scales. Thanks to the success of the theory in finding applications (see, e.g., [30, 22] for recent enumerations), there is a growing demand for efficient computations of the involved topological invariants.

In this paper, we consider a sequence of simplicial complexes $(\mathbb{K}_{i})_{i=0,\ldots,m}$ and simplicial maps $\phi_{i}:\mathbb{K}_{i}\rightarrow\mathbb{K}_{i+1}$ connecting them, calling this data a (simplicial) tower of length $m$. Applying the homology functor with an arbitrary field, we obtain a persistence module, a sequence of vector spaces connected by linear maps. Such a module decomposes into a barcode, a collection of intervals, each representing a homological feature in the tower that spans over the specified range of scales.

Our computational problem is to compute the barcode of a given tower efficiently. The most prominent case of a tower is when all maps $f_{i}$ are inclusion maps. In this case one obtains a filtration, a sequence of nested simplicial complexes. A considerable amount of work went into the study of fast algorithms for the filtration case, which culminated in fast software libraries for this task. The more general case of towers recently received growing interest in the context of sparsification technique for the Vietoris-Rips and Čech complexes; see the related work section below for a detailed discussion.

Results.

As our first result, we show that any tower can be efficiently converted into a small filtration with the same barcode. Using the well-known concept of mapping cylinders from algebraic topology [21], it is easy to see that such a conversion is possible in principle. Dey, Fan, and Wang [15] give an explicit construction, called “coning”, for the generalized case of zigzag towers. Using a simple variant of their coning strategy, we obtain a filtration whose size is only marginally larger than the length of the tower in the worst case. Furthermore, we experimentally show that the size is even smaller on realistic instances.

To describe our improved coning strategy, we discuss the case that a simplicial map in the tower contracts two vertices $u$ and $v$. The coning strategy by Dey et al. proposes to join $u$ with the closed star of $v$, making all incident simplices of $v$ incident to $u$ without changing the homotopy type. The vertex $u$ is then taken as the representative of the contracted pair in the further processing of the tower. We refer to the number of simplices that the join operation adds to the complex as the cost of the contraction. Quite obviously, the method is symmetric in $u$ and $v$, and we have two choices to pick the representative, leading to potentially quite different costs. We employ the self-evident strategy to pick the representative that leads to smaller costs (somewhat reminiscent of the “union-by-weight” rule in the union-find data structure [13, §21]). Perhaps surprisingly, this idea leads to an asymptotically improved size bound on the filtration. We prove this by an abstraction to path decompositions on weighted forest which might be of some independent interest. Altogether, the worst-case size of the filtration is $O(\Delta\cdot n\cdot\log(n_{0}))$, where $\Delta$ is the maximal dimension of any complex in the tower, and $n$/$n_{0}$ is the number of simplices/vertices added to the tower.

We also provide a conversion algorithm whose time complexity is roughly proportional to the total number of simplices in the resulting filtration. One immediate benefit is a generic solution to compute barcodes of towers: just convert the tower to a filtration and apply one of the efficient implementations for barcodes of filtrations. Indeed, we experimentally show that on not-too-large towers, our approach is competitive with, and sometimes outperforms Simpers, an alternative approach that computes the barcode of towers with annotations, a variant of the persistent cohomology algorithm.

Our second contribution is a space-efficient version of the just mentioned algorithmic pipeline that is applicable to very large towers. To motivate the result, let the width of a tower denote the maximal size of any simplicial complex among the $\mathbb{K}_{i}$. Consider a tower with a very large length (say, larger than the number of bytes in main memory) whose width remains relatively small. In this case, our conversion algorithm yields a filtration that is very large as well. Most existing implementations for barcode computation read the entire filtration on initialization and must be converted to streaming algorithm to handle such instances. Moreover, algorithms based on matrix reduction are required to keep previously reduced columns because they might be needed in subsequent reduction steps. This leads to a high memory consumption for the barcode computation.

We show that with minor modifications, the standard persistent algorithm can be turned into a streaming algorithm with smaller space complexity in the case of towers. The idea is that upon contractions, simplices become inactive and cannot get additional cofaces. Our approach makes use of this observation by modifying the boundary matrix such that columns associated to inactive simplices can be removed. Combined with our conversion algorithm, we can compute the barcode of a tower of width $\omega$ keeping only up to $O(\omega)$ columns of the boundary matrix in memory. This yields a space complexity of $O(\omega^{2})$ and a time complexity of $O((\Delta\cdot n\cdot\log(n_{0}))\omega^{2})$ in the worst case. We implemented a practically improved variant that makes use of additional heuristics to speed up the barcode computation in practice and resembles the chunk algorithm presented in [1].

We tested our implementation on various challenging data sets. The source code of the implementation is available111at https://bitbucket.org/schreiberh/sophia/ and the software was named Sophia.

Related work.

Already the first works on persistent homology pointed out the existence of efficient algorithm to compute the barcode invariant (or equivalently, the persistent diagram) for filtrations [18, 33]. As a variant of Gaussian elimination, the worst-case complexity is cubic. Remarkable theoretical follow-up results are a persistence algorithm in matrix multiplication time [27], an output-sensitive algorithm to compute only high-persistent features with linear space complexity [11], and a conditional lower bound on the complexity relating the problem to rank computations of sparse matrices [19].

On realistic instances, the standard algorithm has shown a quasi-linear behavior in practice despite its pessimistic worst-case complexity. Nevertheless, many improvements of the standard algorithm have been presented in the last years which improve the runtime by several orders of magnitude. One line of research exploits the special structure of the boundary matrix to speed up the reduction process [10]. This idea has led to efficient parallel algorithms for persistence in shared [1] and distributed memory [2]. Moreover, of same importance as the reduction strategy is an appropriate choice of data structures in the reduction process as demonstrated by the PHAT library [3]. A parallel development was the development of dual algorithms using persistent cohomology, based on the observation that the resulting barcode is identical [14]. The annotation algorithm [15, 4] is an optimized variant of this idea realized in the Gudhi library [25]. It is commonly considered as an advantage of annotations that only a cohomology basis must be kept during the reduction process, making it more space efficient than reduction-based approaches. We refer to the comparative study [29] for further approaches and software for persistence on filtrations.

Moreover, generalizations of the persistence paradigm are an active field of study. Zigzag persistence is a variant of persistence where the maps in the filtration are allowed to map in either direction (that is, either $\phi_{i}:\mathbb{K}_{i}\hookrightarrow\mathbb{K}_{i+1}$ or $\phi_{i}:\mathbb{K}_{i}\hookleftarrow\mathbb{K}_{i+1}$). The barcode of zigzag filtrations is well-defined [8] as a consequence of Gabriel’s theorem [20] on decomposable quivers – see [30] for a comprehensive introduction. The initial algorithms to compute this barcode [9] has been improved recently [26]. Our case of towers of complexes and simplicial maps can be modeled as a zigzag filtration and therefore sits in-between the standard and the zigzag filtration case.

Dey et al. [15] described the first efficient algorithm to compute the barcode of towers. Instead of the aforementioned coning approach explained in their paper, their implementation handles contractions with an empirically smaller number of insertions, based on the link condition. Recently, the authors have released the SimPers library222http://web.cse.ohio-state.edu/~tamaldey/SimPers/Simpers.html that implements their annotation algorithm from the paper.

The case of towers has received recent attention in the context of approximate Vietoris-Rips and Čech filtrations. The motivation for approximation is that the (exact) topological analysis of a set of $n$ points in $d$-dimensions requires a filtration of size $O(n^{d+1})$ which is prohibitive for most interesting input sizes. Instead, one aims for a filtration or tower of much smaller size, with the guarantee that the approximate barcode will be close to the exact barcode (“close” usually means that the bottleneck distance between the barcodes on the logarithmic scale is upper bounded; we refer to the cited works for details). The first such type of result by Sheehy [31] resulted in a approximate filtration; however, it has been observed that passing to towers allows more freedom in defining the approximation complexes and somewhat simplifies the approximation schemes conceptually. See [15, 6, 24, 12] for examples. Very recently, the SimBa library [16] brings these theoretical approximation techniques for Vietoris-Rips complexes into practice. The approach consists of a geometric layer to compute a tower, and an algebraic layer to compute its barcode, for which they use SimPers. Our approach can be seen as an alternative realization of this algebraic layer.

This paper is a more complete version of the conference paper [23]. It provides missing proof details from [23] and a conclusion, both omitted in the former version for space restrictions. In addition, the experimental results were redone with the most recent versions of the corresponding software libraries. Furthermore, the new subsection 3.5 discusses the tightness of the complexity bound of our first main result.

Outline.

We introduce the necessary basic concepts in Section 2. We describe our conversion algorithm from general towers to barcodes in Section 3. The streaming algorithm for persistence is discussed in Section 4.

2 Background

Simplicial Complexes.

Given a finite vertex set $V$, a simplex is merely a non-empty subset of $V$; more precisely, a $k$-simplex is a subset consisting of $k+1$ vertices, and $k$ is called the dimension of the simplex. Throughout the paper, we will denote simplices by small Greek letters, except for vertices ([math]-simplices) which we denote by $u$, $v$, and $w$. For a $k$-simplex $\sigma$, a simplex $\tau$ is a face of $\sigma$ if $\tau\subseteq\sigma$. If $\tau$ is of dimension $\ell$, we call it a $\ell$-face. If $\ell<k$, we call $\tau$ a proper face of $\sigma$, and if $\ell=k-1$, we call it a facet. For a simplex $\sigma$ and a vertex $v\notin\sigma$, we define the join $v\ast\sigma$ as the simplex $\{v\}\cup\sigma$. These definition are inspired by visualizing a $k$-simplex as the convex hull of $k+1$ affinely independent points in $\mathbb{R}^{k}$, but we will not need this geometric interpretation in our arguments.

An (abstract) simplicial complex $\mathbb{K}$ over $V$ is a set of simplices that is closed under taking faces. We call $V$ the vertex set of $\mathbb{K}$ and write $\mathcal{V}(\mathbb{K}):=V$. The dimension of $\mathbb{K}$ is the maximal dimension of its simplices. For $\sigma,\tau\in\mathbb{K}$, we call $\sigma$ a coface of $\tau$ in $\mathbb{K}$ if $\tau$ is a face of $\sigma$. In this case, $\sigma$ is a cofacet of $\tau$ if their dimensions differ by exactly one. A simplicial complex $\mathbb{L}$ is a subcomplex of $\mathbb{K}$ if $\mathbb{L}\subseteq\mathbb{K}$. Given $\mathcal{W}\subseteq\mathcal{V}$, the induced subcomplex by $\mathcal{W}$ is the set of all simplices $\sigma$ in $\mathbb{K}$ with $\sigma\subseteq\mathcal{W}$. For a subcomplex $\mathbb{L}\subseteq\mathbb{K}$ and a vertex $v\in\mathcal{V}(\mathbb{K})\setminus\mathcal{V}(\mathbb{L})$, we define the join $v\ast\mathbb{L}:=\{v\ast\sigma\mid\sigma\in\mathbb{L}\}$. For a vertex $v\in\mathbb{K}$, the star of $v$ in $\mathbb{K}$, denoted by $\mathrm{St}(v,\mathbb{K})$, is the set of all cofaces of $v$ in $\mathbb{K}$. In general, the star is not a subcomplex, but we can make it a subcomplex by adding all faces of star simplices, which is denoted by the closed star $\overline{\mathrm{St}}(v,\mathbb{K})$. Equivalently, the closed star is the smallest subcomplex of $\mathbb{K}$ containing the star of $v$. The link of $v$, $\mathrm{Lk}(v,\mathbb{K})$, is defined as $\overline{\mathrm{St}}(v,\mathbb{K})\setminus\mathrm{St}(v,\mathbb{K})$. It can be checked that the link is a subcomplex of $\mathbb{K}$. When the complex is clear from context, we will sometimes omit the $\mathbb{K}$ in the notation of stars and links.

Simplicial maps.

A map $\mathbb{K}\stackrel{{\scriptstyle\phi}}{{\rightarrow}}\mathbb{L}$ between simplicial complexes is called simplicial if with $\sigma=\{v_{0},\ldots,v_{k}\}\in\mathbb{K}$, $\phi(\sigma)$ is equal to $\{\phi(v_{0}),\ldots,\phi(v_{k})\}$ and $\phi(\sigma)$ is a simplex in $\mathbb{L}$. By definition, a simplicial map maps vertices to vertices and is completely determined by its action on the vertices. Moreover, the composition of simplicial maps is again simplicial.

A simple example of a simplicial map is the inclusion map $\mathbb{L}\stackrel{{\scriptstyle\phi}}{{\hookrightarrow}}\mathbb{K}$ where $\mathbb{L}$ is a subcomplex of $\mathbb{K}$. If $\mathbb{K}=\mathbb{L}\cup\{\sigma\}$ with $\sigma\notin\mathbb{L}$, we call $\phi$ an elementary inclusion. The simplest example of a non-inclusion simplicial map is $\mathbb{K}\stackrel{{\scriptstyle\phi}}{{\rightarrow}}\mathbb{L}$ such that there exist two vertices $u,v\in\mathbb{K}$ with $\mathcal{V}(\mathbb{L})=\mathcal{V}(\mathbb{K})\setminus\{v\}$, $\phi(u)=\phi(v)=u$, and $\phi$ is the identity on all remaining vertices of $\mathbb{K}$. We call $\phi$ an elementary contraction and write $(u,v)\leadsto u$ as a shortcut. These notions were introduced by Dey, Fan and Wang in [15] and they also showed that any simplicial map $\mathbb{K}\stackrel{{\scriptstyle\phi}}{{\rightarrow}}\mathbb{L}$ can be written as the composition of elementary contractions333 They talk about ”collapses” instead of ”contractions”, but this notion clashes with the standard notion of simplicial collapses of free faces that we use later. Therefore, we decided to use ”contraction”, even though the edge between the contracted vertices might not be present in the complex.

and inclusions.

A tower of length $m$ is a collection of simplicial complexes $\mathbb{K}_{0},\ldots,\mathbb{K}_{m}$ and simplicial maps $\phi_{i}:\mathbb{K}_{i}\rightarrow\mathbb{K}_{i+1}$ for $i=0,\ldots,m-1$. From this initial data, we obtain simplicial maps $\phi_{i,j}:\mathbb{K}_{i}\to\mathbb{K}_{j}$ for $i\leq j$ by composition, where $\phi_{i,i}$ is simply the identity map on $\mathbb{K}_{i}$. The term “tower” is taken from category theory, where it denotes a (directed) path in a category with morphisms from objects with smaller indices to objects with larger indices. Indeed, since simplicial complexes form a category with simplicial maps as their morphisms, the specified data defines a tower in this category. A tower is called a filtration if all $\phi_{i}$ are inclusion maps. The dimension of a tower is the maximal dimension among the $\mathbb{K}_{i}$, and the width of a tower is the maximal size among the $\mathbb{K}_{i}$. For filtrations, dimension and width are determined by the dimension and size of $\mathbb{K}_{m}$, but this is not necessarily true for general towers.

Homology and Collapses.

For a fixed base field $\mathbb{F}$, let $H_{p}(\mathbb{K}):=H_{p}(\mathbb{K},\mathbb{F})$ the $p$-dimensional homology group of $\mathbb{K}$. It is well-known that $H_{p}(\mathbb{K})$ is a $\mathbb{F}$-vector space. Moreover, a simplicial map $\mathbb{K}\stackrel{{\scriptstyle\phi}}{{\rightarrow}}\mathbb{L}$ induces a linear map $H_{p}(\mathbb{K})\stackrel{{\scriptstyle\phi^{\ast}}}{{\rightarrow}}H_{p}(\mathbb{L})$. In categorical terms, the equivalent statement is that homology is a functor from the category of simplicial complexes and simplicial maps to the category of vector spaces and linear maps.

We will make use of the following homology-preserving operation: a free face in $\mathbb{K}$, is a simplex with exactly one proper coface in $\mathbb{K}$. An elementary collapse in $\mathbb{K}$ is the operation of removing a free face and its unique coface from $\mathbb{K}$, yielding a subcomplex of $\mathbb{K}$. We say that $\mathbb{K}$ collapses to $\mathbb{L}$, if there is a sequence of elementary collapses transforming $\mathbb{K}$ into $\mathbb{L}$. The following result is then well-known:

Lemma 1**.**

Let $\mathbb{K}$ be a complex that collapses into the complex $\mathbb{L}$. Then, the inclusion map $\mathbb{L}\stackrel{{\scriptstyle\phi}}{{\hookrightarrow}}\mathbb{K}$ induces an isomorphism $\phi_{\ast}$ between $H_{p}(\mathbb{L})$ and $H_{p}(\mathbb{K})$.

Barcodes.

A persistence module is a sequence vector spaces $\mathbb{V}_{0},\ldots,\mathbb{V}_{m}$ and linear maps $f_{i,j}:\mathbb{V}_{i}\to\mathbb{V}_{j}$ for $i<j$ such that $f_{i,i}=\mathrm{id}_{\mathbb{V}_{i}}$ and $f_{i,k}=f_{j,k}\circ f_{i,j}$ for $i\leq k\leq j$. As primary example, we obtain a persistence module by applying the homology functor on any simplicial tower. Persistence modules admit a decomposition into indecomposable summands in the following sense. Writing $I_{b,d}$ with $b\leq d$ for the persistence module

[TABLE]

we can write every persistence module as the direct sum $I_{b_{1},d_{1}}\oplus\ldots\oplus I_{b_{s},d_{s}}$, where the direct sum of persistence modules is defined component-wise for vector spaces and linear maps in the obvious way. Moreover, this decomposition is uniquely defined up to isomorphisms and re-ordering, thus the pairs $(b_{1},d_{1}),\ldots,(b_{s},d_{s})$ are an invariant of the persistence module, called its barcode. When the persistence module was generated by a tower, we also talk about the barcode of the tower.

Matrix reduction.

In this paragraph, we assume that $(\mathbb{K}_{i})_{i=0,\ldots,m}$ is a filtration such that $\mathbb{K}_{0}=\emptyset$ and $\mathbb{K}_{i+1}$ has exactly one more simplex than $\mathbb{K}_{i}$. We label the simplices of $\mathbb{K}_{m}$ accordingly as $\sigma_{1},\ldots,\sigma_{m}$, with $\mathbb{K}_{i}\setminus\mathbb{K}_{i-1}=\{\sigma_{i}\}$. The filtration can be encoded as a boundary matrix $\partial$ of dimension $m\times m$, where the $(ij)$-entry is $1$ if $\sigma_{i}$ is a facet of $\sigma_{j}$, and 0 otherwise. In other words, the $j$-th column of $\partial$ encodes the facets of $\sigma_{j}$, and the $i$-th row of $\partial$ encodes the cofacets of $\sigma_{i}$. Moreover, $\partial$ is upper-triangular because every $\mathbb{K}_{i}$ is a simplicial complex. We will sometimes identify rows and columns in $\partial$ with the corresponding simplex in $\mathbb{K}_{m}$. Adding the $k$-simplex $\sigma_{i}$ to $\mathbb{K}_{i-1}$ either introduces one new homology class (of dimension $k$), or turns a non-trivial homology class (of dimension $k-1$) trivial. We call $\sigma_{i}$ and the $i$-th column of $\partial$ positive or negative, respectively (with respect to the given filtration).

For the computation of the barcode, we assume for simplicity homology over the base field $\mathbb{Z}_{2}$, and interpret the coefficients of $\partial$ accordingly. In an arbitrary matrix $A$, a left-to-right column addition is an operation of the form $A_{k}\leftarrow A_{k}+A_{\ell}$ with $\ell<k$, where $A_{k}$ and $A_{\ell}$ are columns of the matrix. The pivot of a non-zero column is the largest non-zero index of the corresponding column. A non-zero entry is called a pivot if its row is the pivot of the column. A matrix $R$ is called a reduction of $A$ if $R$ is obtained by a sequence of left-to-right column additions from $A$ and no two columns in $R$ have the same pivot. It is well-known that, although $\partial$ does not have a unique reduction, the pivots of all its reductions are the same. Moreover, the pivots $(b_{1},d_{1}),\ldots,(b_{s},d_{s})$ of $R$ are precisely the barcode of the filtration. A direct consequence is that a simplex $\sigma_{i}$ is positive if and only if the $i$-th column in $R$ is zero.

The standard persistence algorithm processes the columns from left to right. In the $j$-th iteration, as long as the $j$-th column is not empty and has a pivot that appears in a previous column, it performs a left-to-right column addition. In this work, we use a simple improvement of this algorithm that is called compression: before reducing the $j$-th column, it first scans through the non-zero entries of the column. If a row index $i$ corresponds to a negative simplex (i.e., if the $i$-th column is not zero at this point in the algorithm), the row index can be deleted without changing the pivots of the matrix. After this initial scan, the column is reduced in the same way as in the standard algorithm. See [1, §. 3] for a discussion (we remark that this optimization was also used in [33]).

3 From towers to filtrations

We phrase now our first result which says that any tower can be converted into a filtration of only marginally larger size with a space-efficient streaming algorithm:

Theorem 2** (Conversion Theorem).**

Let $\mathcal{T}:\,\mathbb{K}_{0}\xrightarrow{\phi_{0}}\mathbb{K}_{1}\xrightarrow{\phi_{1}}\ldots\xrightarrow{\phi_{m-1}}\mathbb{K}_{m}$ be a tower where, w.l.o.g., $\mathbb{K}_{0}=\emptyset$ and each $\phi_{i}$ is either an elementary inclusion or an elementary contraction. Let $\Delta$ denote the dimension and $\omega$ the width of the tower, and let $n\leq m$ denote the total number of elementary inclusions, and $n_{0}$ the number of vertex inclusions. Then, there exists a filtration $\mathcal{F}:\,\hat{\mathbb{K}}_{0}\xhookrightarrow{}\hat{\mathbb{K}}_{1}\xhookrightarrow{}\ldots\xhookrightarrow{}\hat{\mathbb{K}}_{m}$, where the inclusions are not necessarily elementary, such that $\mathcal{T}$ and $\mathcal{F}$ have the same barcode and the width of the filtration $|\hat{\mathbb{K}}_{m}|$ is at most $O(\Delta\cdot n\log n_{0})$. Moreover, $\mathcal{F}$ can be computed from $\mathcal{T}$ with a streaming algorithm in $O(\Delta\cdot|\hat{\mathbb{K}}_{m}|\cdot C_{\omega})$ time and space complexity $O(\Delta\cdot\omega)$, where $C_{\omega}$ is the cost of an operation in a dictionary with $\omega$ elements.

The remainder of the section is organized as follows. We define $\mathcal{F}$ in Section 3.1 and prove that it yields the same barcode as $\mathcal{T}$ in Section 3.2. In Section 3.3, we prove the upper bound on the width of the filtration. In Section 3.4, we explain the algorithm to compute $\mathcal{F}$ and analyze its time and space complexity.

3.1 Active and small coning

Coning.

We briefly revisit the coning strategy introduced by Dey, Fan and Wang [15]. Let $\phi:\mathbb{K}\to\mathbb{L}$ be an elementary contraction $(u,v)\leadsto u$ and define

[TABLE]

An example is shown in Figure 1.

Dey et al. show that $\mathbb{L}\subseteq\mathbb{L}^{\ast}$ and that the map induced by inclusion is an isomorphism between $H(\mathbb{L})$ and $H(\mathbb{L}^{\ast})$. By applying this result at any elementary contraction, this implies that every zigzag tower can be transformed into a zigzag filtration with identical barcode.

Given a tower $\mathcal{T}$, we can also obtain an non-zigzag filtration using coning, if we continue the operation on $\mathbb{L}^{\ast}$ instead of going back to $\mathbb{L}$. More precisely, we set $\tilde{\mathbb{K}}_{0}:=\mathbb{K}_{0}$ and if $\phi_{i}$ is an inclusion of simplex $\sigma$, we set $\tilde{\mathbb{K}}_{i+1}:=\tilde{\mathbb{K}}_{i}\cup\{\sigma\}$. If $\phi_{i}$ is a contraction $(u,v)\leadsto u$, we set $\tilde{\mathbb{K}}_{i+1}=\tilde{\mathbb{K}}_{i}\cup\left(u\ast\overline{\mathrm{St}}(v,\tilde{\mathbb{K}}_{i})\right)$. Indeed, it can be proved that $(\tilde{\mathbb{K}}_{i})_{i=0,\ldots,m}$ has the same barcode as $\mathcal{T}$. However, the filtration will not be small, and we will define a smaller variant now.

Our new construction yields a sequence of complexes $\hat{\mathbb{K}}_{0},\ldots,\hat{\mathbb{K}}_{m}$ with $\hat{\mathbb{K}}_{i}\subseteq\hat{\mathbb{K}}_{i+1}$. During the construction, we maintain a flag for each vertex in $\hat{\mathbb{K}}_{i}$, which marks the vertex as active or inactive. A simplex is called active if all its vertices are active, and inactive otherwise. For a vertex $u$ and a complex $\hat{\mathbb{K}}_{i}$, let $\mathrm{Act\overline{\mathrm{St}}}(u,\hat{\mathbb{K}}_{i})$ denote its active closed star, which is the set of active simplices in $\hat{\mathbb{K}}_{i}$ in the closed star of $u$.

The construction is inductive, starting with $\hat{\mathbb{K}}_{0}:=\emptyset$. If $\mathbb{K}_{i}\stackrel{{\scriptstyle\phi_{i}}}{{\to}}\mathbb{K}_{i+1}$ is an elementary inclusion with $\mathbb{K}_{i+1}=\mathbb{K}_{i}\cup\{\sigma\}$, set $\hat{\mathbb{K}}_{i+1}:=\hat{\mathbb{K}}_{i}\cup\{\sigma\}$. If $\sigma$ is a vertex, we mark it as active. It remains the case that $\mathbb{K}_{i}\stackrel{{\scriptstyle\phi_{i}}}{{\to}}\mathbb{K}_{i+1}$ is an elementary contraction of the vertices $u$ and $v$. If $|\mathrm{Act\overline{\mathrm{St}}}(u,\hat{\mathbb{K}}_{i})|\leq|\mathrm{Act\overline{\mathrm{St}}}(v,\hat{\mathbb{K}}_{i})|$, we set

[TABLE]

and mark $u$ as inactive. Otherwise, we set

[TABLE]

and mark $v$ as inactive. This ends the description of the construction. We write $\mathcal{F}$ for the filtration $(\hat{\mathbb{K}}_{i})_{i=0,\ldots,m}$.

There are two major changes compared to the construction of $(\tilde{\mathbb{K}}_{i})_{i=0,\ldots,m}$. First, to counteract the potentially large growth of the involved cones, we restrict coning to active simplices. We will show below that the subcomplex of $\hat{\mathbb{K}}_{i}$ induced by the active vertices is isomorphic to $\mathbb{K}_{i}$. As a consequence, we add the same number of simplices by passing from $\hat{\mathbb{K}}_{i}$ to $\hat{\mathbb{K}}_{i+1}$ as in the approach by Dey et al. does when passing from $\mathbb{K}$ to $\mathbb{L}^{\ast}$.

A second difference is that our strategy exploits that an elementary contraction of two vertices $u$ and $v$ leaves us with a choice: we can either take $u$ or $v$ as the representative of the contracted vertex. In terms of simplicial maps, these two choices correspond to setting $\phi_{i}(u)=\phi_{i}(v)=u$ or $\phi_{i}(u)=\phi_{i}(v)=v$, if $\phi_{i}$ is the elementary contraction of $u$ and $v$. It is obvious that both choices yield identical complexes $\mathbb{K}_{i+1}$ up to renaming of vertices. However, the choices make a difference in terms of the size of $\hat{\mathbb{K}}_{i+1}$, because the active closed star of $u$ to $v$ in $\hat{\mathbb{K}}_{i}$ might differ in size. Our construction simply choose the representative which causes the smaller $\hat{\mathbb{K}}_{i+1}$.

3.2 Topological equivalence

We make the following simplifying assumption for $\mathcal{T}$. Let $\phi_{i}$ be an elementary contraction of $u$ and $v$. If our construction of $\hat{\mathbb{K}}_{i+1}$ turns $v$ inactive, we assume that $\phi_{i}(u)=\phi_{i}(v)=u$. Otherwise, we assume $\phi_{i}(u)=\phi_{i}(v)=v$. Indeed, this is without loss of generality because it corresponds to a renaming of the simplices in each $\mathbb{K}_{i}$ and yields equivalent persistence modules. The advantage of this convention is the following property, which follows from a straight-forward inductive argument.

Lemma 3**.**

For every $i$ in $\{0,...,m\}$, the set of vertices of $\mathbb{K}_{i}$ is equal to the set of active vertices in $\hat{\mathbb{K}}_{i}$.

This allows us to interpret $\mathbb{K}_{i}$ and $\hat{\mathbb{K}}_{i}$ as simplicial complexes defined over a common vertex set. In fact, $\mathbb{K}_{i}$ is the subcomplex of $\hat{\mathbb{K}}_{i}$ induced by the active vertices:

Lemma 4**.**

A simplex $\sigma$ is in $\mathbb{K}_{i}$ if and only if $\sigma$ is an active simplex in $\hat{\mathbb{K}}_{i}$.

Proof.

We use induction on $i$. The statement is true for $i=0$, because $\mathbb{K}_{0}=\emptyset=\hat{\mathbb{K}}_{0}$. So assume first $\phi_{i}:\mathbb{K}_{i}\to\mathbb{K}_{i+1}$ is an elementary inclusion that adds a $d$-simplex $\sigma=(v_{0},\ldots,v_{d})$ to $\mathbb{K}_{i+1}$. If $\sigma$ is a vertex, it is set active in $\hat{\mathbb{K}}_{i+1}$ by construction. Otherwise, $v_{0},\ldots,v_{d}$ are active in $\hat{\mathbb{K}}_{i}$ by induction and stay active in $\hat{\mathbb{K}}_{i+1}$. In any case, $\sigma$ is active in $\hat{\mathbb{K}}_{i+1}$. The equivalence for the remaining simplices is straight-forward.

If $\phi_{i}$ is an elementary contraction $(u,v)\leadsto u$, we prove both directions of the equivalence separately. For “$\Rightarrow$”, fix a $d$-simplex $\sigma\in\mathbb{K}_{i+1}$. It suffices to show that $\sigma\in\hat{\mathbb{K}}_{i+1}$, as in this case, it is also active by Lemma 3. If $\sigma\in\mathbb{K}_{i}$, this follows at once by induction because $\mathbb{K}_{i}\subseteq\hat{\mathbb{K}}_{i}\subseteq\hat{\mathbb{K}}_{i+1}$. If $\sigma\notin\mathbb{K}_{i}$, $u$ must be a vertex of $\sigma$. Moreover, writing $\sigma=\{u,v_{1},...,v_{d}\}$ and $\sigma^{\prime}=\{v,v_{1},\ldots,v_{d}\}$ we have that $\sigma^{\prime}\in\mathbb{K}_{i}$ and $\phi_{i}(\sigma^{\prime})=\sigma$. In particular, the vertices $v_{1},\ldots,v_{d}$ are active in $\hat{\mathbb{K}}_{i}$ by induction, hence $\{v_{1},\ldots,v_{d}\}$ is in the active closed star of $v$ in $\hat{\mathbb{K}}_{i}$. By construction, $\{u,v_{1},...,v_{d}\}=\sigma$ is in $\hat{\mathbb{K}}_{i+1}$.

For “$\Leftarrow$”, let $\sigma\in\hat{\mathbb{K}}_{i+1}\setminus\mathbb{K}_{i+1}$. We show that $\sigma$ is an inactive simplex in $\hat{\mathbb{K}}_{i+1}$. By Lemma 3, this is equivalent to show that $\sigma$ contains a vertex not in $\mathbb{K}_{i+1}$.

Case 1: $\sigma\in\hat{\mathbb{K}}_{i}$. If $\sigma$ is inactive in $\hat{\mathbb{K}}_{i}$, it stays inactive in $\hat{\mathbb{K}}_{i+1}$. So, assume that $\sigma$ is active in $\hat{\mathbb{K}}_{i}$ and thus $\sigma\in\mathbb{K}_{i}$ by induction. But $\sigma\notin\mathbb{K}_{i+1}$, so $\sigma$ must have $v$ as a vertex and is therefore inactive in $\hat{\mathbb{K}}_{i+1}$.

Case 2: $\sigma\in\hat{\mathbb{K}}_{i+1}\setminus\hat{\mathbb{K}}_{i}$. By construction, $\sigma$ is of the form $\{u,v_{1},\ldots,v_{d}\}$ such that $\{v_{1},\ldots,v_{d}\}$ is in the active closed star of $v$ in $\hat{\mathbb{K}}_{i}$. Assume for a contradiction that $v\neq v_{j}$ for all $j=1,\ldots,d$. Then, $\sigma^{\prime}=\{v,v_{1},\ldots,v_{d}\}$ is active in $\hat{\mathbb{K}}_{i}$ and thus, by induction, a simplex in $\mathbb{K}_{i}$. But then, $\phi_{i}(\sigma^{\prime})=\sigma\in\mathbb{K}_{i+1}$ which is a contradiction to our choice of $\sigma$. It follows that $v$ is a vertex of $\sigma$ which proves our claim. ∎

Lemma 5**.**

For every $0\leq i\leq m$, the complex $\hat{\mathbb{K}}_{i}$ collapses to $\mathbb{K}_{i}$.

Proof.

We use induction on $i$. For $i=0$, $\mathbb{K}_{0}=\hat{\mathbb{K}}_{0}$, and the statement is trivial. Suppose that the statement holds for $\hat{\mathbb{K}}_{i}$ and $\mathbb{K}_{i}$ using the sequence $s_{i}$ of elementary collapses. Note that these collapses only concern inactive simplices in $\hat{\mathbb{K}}_{i}$. For an inactive vertex $v\in\hat{\mathbb{K}}_{i}$, the construction of $\hat{\mathbb{K}}_{i+1}$ ensures that $v$ does not gain any additional coface. This implies that $s_{i}$ is still a sequence of elementary collapses for $\hat{\mathbb{K}}_{i+1}$, yielding a complex $\hat{\mathbb{K}}^{*}_{i+1}$ with $\mathbb{K}_{i+1}\subseteq\hat{\mathbb{K}}^{*}_{i+1}\subseteq\hat{\mathbb{K}}_{i+1}$. In particular, $\hat{\mathbb{K}}^{*}_{i+1}$ only contains vertices that are still active in $\hat{\mathbb{K}}_{i}$. If $\phi_{i}$ is an elementary inclusion, $\hat{\mathbb{K}}^{*}_{i+1}=\mathbb{K}_{i+1}$, because any all vertices in $\hat{\mathbb{K}}_{i}$ remain active in $\hat{\mathbb{K}}_{i+1}$. For $\phi_{i}$ being an elementary contraction $(u,v)\leadsto u$, set $S:=\hat{\mathbb{K}}^{*}_{i+1}\setminus\mathbb{K}_{i+1}$ as the remaining set of simplices that still need to be collapsed to obtain $\mathbb{K}_{i+1}$. All simplices of $S$ have $v$ as vertex. More precisely, $S$ is the set of all simplices of the form $\{v,v_{1},\ldots,v_{d}\}$ with $v_{1},\ldots,v_{d}$ active in $\hat{\mathbb{K}}_{i+1}$. We split $S=S_{u}\cup S_{\neg u}$ where $S_{u}\subset S$ are the simplices in $S$ that contain $u$ as vertex, and $S_{\neg u}=S\setminus S_{u}$.

We claim that the mapping that sends $\{u,v,v_{1},\ldots,v_{d}\}\in S_{u}$ to $\{v,v_{1},\ldots,v_{d}\}\in S_{\neg u}$ is a bijection. This map is clearly injective. If $\sigma=\{v,v_{1},\ldots,v_{d}\}\in S_{\neg u}$, then $\sigma\in\hat{\mathbb{K}}_{i}$ (because every newly added simplex in $\hat{\mathbb{K}}_{i+1}$ contains $u$). Also, $\sigma\in\hat{\mathbb{K}}^{*}_{i+1}$, and is therefore active in $\hat{\mathbb{K}}_{i}$. By construction, $\{u,v,v_{1},\ldots,v_{d}\}\in\hat{\mathbb{K}}_{i+1}$, proving surjectivity.

We now define a sequence of elementary collapses from $\hat{\mathbb{K}}^{*}_{i+1}$ to $\mathbb{K}_{i+1}$. Choose a simplex $\sigma=\{v,v_{1},\ldots,v_{d}\}\in S_{\neg u}$ of maximal dimension, and let $\tau=\{u,v,v_{1},\ldots,v_{d}\}$ denote the corresponding simplex in $S_{u}$. Then $\sigma$ is indeed a free face in $\hat{\mathbb{K}}^{*}_{i+1}$, because if there was another coface $\tau^{\prime}\neq\tau$, it takes the form $\{w,v,v_{1},\ldots,v_{d}\}$ with $w\neq u$ active. So, $\tau^{\prime}\in S_{\neg u}$, and $\tau^{\prime}$ has larger dimension than $\sigma$, a contradiction. Therefore, the pair $(\sigma,\tau)$ defines an elementary collapse in $\hat{\mathbb{K}}^{*}_{i+1}$. We proceed with this construction, always collapsing a remaining pair in $S_{\neg u}\times S_{u}$ of maximal dimension, until all elements of $S$ have been collapsed. ∎

Proposition 6**.**

$\mathcal{T}$* and $\mathcal{F}$ have the same barcode.*

Proof.

Let $\hat{\phi}_{i}:\hat{\mathbb{K}}_{i}\to\hat{\mathbb{K}}_{i+1}$ denote the inclusion map from $\hat{\mathbb{K}}_{i}$ to $\hat{\mathbb{K}}_{i+1}$. By combining Lemma 5 with Lemma 1, we have an isomorphism $\mathrm{inc}_{i}^{*}:H(\mathbb{K}_{i})\to H(\hat{\mathbb{K}}_{i})$, for all $0\leq i\leq m$, induced by the inclusion maps $\mathrm{inc}_{i}:\mathbb{K}_{i}\to\hat{\mathbb{K}}_{i}$, and therefore the following diagram connecting the persistence modules induced by $\mathcal{T}$ and $\mathcal{F}$:

[TABLE]

The Persistence Equivalence Theorem [17, p.159] asserts that $(\mathbb{K}_{j})_{j}$ and $(\hat{\mathbb{K}}_{j})_{j}$, with $j=0,...,m$, have the same barcode if (1) commutes, that is, if $\mathrm{inc}_{i+1}^{*}\circ\phi_{i}^{*}=\hat{\phi}_{i}^{*}\circ\mathrm{inc}_{i}^{*}$, for all $0\leq i<m$.

Two simplicial maps $\phi:\mathbb{K}\to\mathbb{L}$ and $\psi:\mathbb{K}\to\mathbb{L}$ are contiguous if, for all $\sigma\in\mathbb{K}$, $\phi(\sigma)\cup\psi(\sigma)\in\mathbb{L}$. Two contiguous maps are known to be homotopic [28, Theorem 12.5.] and thus equal at homology level, that is, $\phi^{\ast}=\psi^{\ast}$. We show that $\mathrm{inc}_{i+1}\circ\phi_{i}$ and $\hat{\phi}_{i}\circ\mathrm{inc}_{i}$ are contiguous. This implies that (1) commutes, because, by functoriality, $\mathrm{inc}_{i+1}^{*}\circ\phi_{i}^{*}=(\mathrm{inc}_{i+1}\circ\phi_{i})^{*}=(\hat{\phi}_{i}\circ\mathrm{inc}_{i})^{*}=\hat{\phi}_{i}^{*}\circ\mathrm{inc}_{i}^{*}$.

To show contiguity, fix $\sigma\in\mathbb{K}_{i}$ and observe that $(\hat{\phi}_{i}\circ\mathrm{inc}_{i})(\sigma)=\sigma$ because $\hat{\phi}_{i}$ and $\mathrm{inc}_{i}$ are inclusions. If $\phi_{i}(\sigma)=\sigma$, $(\mathrm{inc}_{i+1}\circ\phi_{i})(\sigma)=\sigma$ as well, and the contiguity condition is clearly satisfied. So, let $\phi_{i}(\sigma)\neq\sigma$. Then $\phi_{i}$ is an elementary contraction $(u,v)\leadsto u$, and $\sigma$ is of the form $\{v,v_{1},\ldots,v_{d}\}$, where one of the $v_{j}$ might be equal to $u$. Then, $(\mathrm{inc}_{i+1}\circ\phi_{i})(\sigma)=\{u,v_{1},...,v_{d}\}$. Consequently, $(\mathrm{inc}_{i+1}\circ\phi_{i})(\sigma)\cup(\hat{\phi}_{i}\circ\mathrm{inc}_{i})(\sigma)=\{u,v,v_{1},...,v_{d}\}$. By Lemma 4, $\sigma=\{v,v_{1},\ldots,v_{d}\}$ is in the active closed star of $v$ in $\hat{\mathbb{K}}_{i}$, and by construction $\{u,v,v_{1},...,v_{d}\}\in\hat{\mathbb{K}}_{i+1}$, which proves contiguity of the maps. ∎

3.3 Size analysis

The contracting forest.

We associate a rooted labeled forest $\mathcal{W}_{j}$ to a prefix $\emptyset=\mathbb{K}_{0}\xrightarrow{\phi_{0}}\ldots\xrightarrow{\phi_{j-1}}\mathbb{K}_{j}$ of $\mathcal{T}$ inductively as follows: For $j=0$, $\mathcal{W}_{0}$ is the empty forest. Let $\mathcal{W}_{j-1}$ be the forest of $\mathbb{K}_{0}\rightarrow\ldots\rightarrow\mathbb{K}_{j-1}$. If $\phi_{j-1}$ is an elementary inclusion of a $d$-simplex, we have two cases: if $d>0$, set $\mathcal{W}_{j}:=\mathcal{W}_{j-1}$. If a vertex $v$ is included, $\mathcal{W}_{j}:=\mathcal{W}_{j-1}\cup\{x\}$, with $x$ a single node tree labeled with $v$. If $\phi_{j-1}$ is an elementary contraction contracting two vertices $u$ and $v$ in $\mathbb{K}_{j-1}$, there are two trees in $\mathcal{W}_{j-1}$, whose roots are labeled $u$ and $v$. In $\mathcal{W}_{j}$, these two trees are merged by making their roots children of a new root, which is labeled with the vertex that $u$ and $v$ are mapped to.

We can read off from the construction immediately that the roots of $\mathcal{W}_{i}$ are labeled with the vertices of the complex $\mathbb{K}_{i}$, for every $i=0,\ldots,m$. Moreover, each leaf corresponds to the inclusion of a vertex in $\mathbb{K}_{0}\rightarrow\ldots\rightarrow\mathbb{K}_{i}$, and each internal node corresponds to a contraction of two vertices. In particular, $\mathcal{W}_{i}$ is a full forest, that is, every node has [math] or $2$ children.

Let $\mathcal{W}:=\mathcal{W}_{m}$ denote the forest of the tower $\mathcal{T}$. Let $\Sigma$ denote the set of all simplices that are added at elementary inclusions in $\mathcal{T}$, and recall that $n=|\Sigma|$. A $d$-simplex $\sigma\in\Sigma$ added by $\phi_{i}$ is formed by $d+1$ vertices, which correspond to $d+1$ roots of $\mathcal{W}_{i+1}$, and equivalently, to $d+1$ nodes in $\mathcal{W}$. For a node $x$ in $\mathcal{W}$, we denote by $E(x)\subseteq\Sigma$ the subset of simplices with at least one vertex that appears as label in the subtree of $\mathcal{W}$ rooted at $x$. If $y_{1}$ and $y_{2}$ are the children of $x$, $E(y_{1})$ and $E(y_{2})$ are both subsets of $E(x)$, but not disjoint in general. However, the following relation follows at once:

[TABLE]

We say that the set $N$ of nodes in $\mathcal{W}$ is independent, if there are no two nodes $x_{1}\neq x_{2}$ in $N$, such that $x_{1}$ is an ancestor of $x_{2}$ in $\mathcal{W}$. A vertex in $\mathbb{K}_{i}$ appears as label in at most one $\mathcal{W}$-subtree rooted at a vertex in the independent set $N$. Thus, a $d$-simplex $\sigma$ can only appear in up to $d+1$ $E$-sets of vertices in $N$. That implies:

Lemma 7**.**

Let $N$ be an independent set of vertices in $\mathcal{W}$. Then,

[TABLE]

The cost of contracting.

Recall that a contraction $\mathbb{K}_{i}\stackrel{{\scriptstyle\phi_{i}}}{{\to}}\mathbb{K}_{i+1}$ yields an inclusion $\hat{\mathbb{K}}_{i}\hookrightarrow\hat{\mathbb{K}}_{i+1}$ that potentially adds more than one simplex. Therefore, in order to bound the total size of $\hat{\mathbb{K}}_{m}$, we need to bound the number of simplices added in all these contractions.

We define the cost of a contraction $\phi_{i}$ as $|\hat{\mathbb{K}}_{i+1}\setminus\hat{\mathbb{K}}_{i}|$, that is, the number of simplices added in this step. Since each contraction corresponds to a node $x$ in $\mathcal{W}$, we can associate these costs to the internal nodes in the forest, denoted by $c(x)$. The leaves get cost [math].

Lemma 8**.**

Let $x$ be an internal node of $\mathcal{W}$ with children $y_{1}$, $y_{2}$. Then, $c(x)\leq 2\cdot|E(y_{1})\setminus E(y_{2})|$.

Proof.

Let $\phi_{i}:\mathbb{K}_{i}\to\mathbb{K}_{i+1}$ denote the contraction that is represented by the node $x$, and let $w_{1}$ and $w_{2}$ the labels of its children $y_{1}$ and $y_{2}$, respectively. By construction, $w_{1}$ and $w_{2}$ are vertices in $\mathbb{K}_{i}$ that are contracted by $\phi_{i}$. Let $C_{1}=\overline{\mathrm{St}}(w_{1},\mathbb{K}_{i})\setminus\overline{\mathrm{St}}(w_{2},\mathbb{K}_{i})$ and $C_{2}=\overline{\mathrm{St}}(w_{2},\mathbb{K}_{i})\setminus\overline{\mathrm{St}}(w_{1},\mathbb{K}_{i})$. By Lemma 4, $\overline{\mathrm{St}}(w_{1},\mathbb{K}_{i})=\mathrm{Act\overline{\mathrm{St}}}(w_{1},\hat{\mathbb{K}}_{i})$, and the same for $w_{2}$. So, because the simplices the two active closed stars have in common will not influence the cost of the contraction, we have,

[TABLE]

In particular, $c(x)\leq|C_{1}|$. We will show that $|C_{1}|\leq 2\cdot|E(y_{1})\setminus E(y_{2})|$.

For every $d$-simplex $\sigma\in\mathbb{K}_{i}$, there is some $d$-simplex $\tau\in\Sigma$ that has been added in an elementary inclusion $\phi_{j}$ with $j<i$, such that $\phi_{i-1}\circ\phi_{i-2}\circ\ldots\circ\phi_{j+1}(\tau)=\sigma$. We call $\tau$ an origin of $\sigma$. There might be more than one origin of a simplex, but two distinct simplices in $\mathbb{K}_{i}$ cannot have a common origin. Moreover, for every vertex $v$ of $\sigma$, the tree in $\mathcal{W}_{i}$ whose root is labeled with $v$ contains exactly one vertex $v^{\prime}$ of $\tau$ as label. We omit the proof which works by simple induction.

We prove the inequality by a simple charging argument. For each vertex in $C_{1}$, we charge a simplex in $|E(y_{1})\setminus E(y_{2})|$ such that no simplex is charged more than twice. Note that $C_{1}=(\mathrm{St}(w_{1},\mathbb{K}_{i})\cup\mathrm{Lk}(w_{1},\mathbb{K}_{i}))\setminus\overline{\mathrm{St}}(w_{2},\mathbb{K}_{i})$. If $\sigma\in\mathrm{St}(w_{1},\mathbb{K}_{i})$, then fix an origin $\tau$ of $\sigma$. Then $\tau$ has a vertex that is a label in the subtree rooted at $y_{1}$, so $\tau\in E(y_{1})$. At the same time, since $w_{2}$ is not a vertex of $\sigma$, $\tau$ has no vertex in the subtree rooted at $y_{2}$, so $\tau\notin E(y_{2})$. We charge $\tau$ for the existence of $\sigma$. Because different simplices have different origins, every element in $E(y_{1})\setminus E(y_{2})$ is charged at most once for $\mathrm{St}(w_{1},\mathbb{K}_{i})$. If $\sigma\in\mathrm{Lk}(w_{1},\mathbb{K}_{i})$, $\sigma^{\prime}:=w_{1}\ast\sigma\in\mathrm{St}(w_{1},\mathbb{K}_{i})$, and we can choose an origin $\tau^{\prime}$ of $\sigma^{\prime}$. As before, $\tau^{\prime}\in E(y_{1})\setminus E(y_{2})$ and we charge $\tau^{\prime}$ for the existence of $\sigma$. Again, each element in $E(y_{1})\setminus E(y_{2})$ is charged at most once among all elements in the link. This proves the claim. ∎

An ascending path $(x_{1},...,x_{L})$, with $L\geq 1$, is a path in a forest such that $x_{i+1}$ is the parent of $x_{i}$, for $1\leq i<L$. We call $L$ the length of the path and $x_{L}$ its endpoint. For ascending paths in $\mathcal{W}$, the cost of the path is the sum of the costs of the nodes. We say that the set $P$ of ascending paths is independent, if the endpoints in $P$ are pairwise different and form an independent set of nodes. We define the cost of $P$ as the sum of the costs of the paths in $P$.

Lemma 9**.**

An ascending path with endpoint $x$ has cost at most $2\cdot|E(x)|$. An independent set of ascending paths in $\mathcal{W}$ has cost at most $2\cdot(\Delta+1)\cdot n$.

Proof.

For the first statement, let $p=(x_{1},...,x_{L})$ be an ascending path with $v_{L}=v$. Without loss of generality, we can assume the the path starts with a leaf $x_{1}$, because otherwise, we can always extend the path to a longer path with at least the same cost. We let $p_{i}=(x_{1},\ldots,x_{i})$ denote the sub-path ending at $x_{i}$, for $i=1,\ldots,L$, so that $p_{L}=p$. We let $c(p_{i})$ denote the cost of the path $p_{i}$ and show by induction that $c(p_{i})\leq 2\cdot|E(x_{i})|$. For $i=1$, this follows because $c(p_{1})=0$. For $i=2,\ldots,L$, $x_{i}$ is an internal node, and its two children are $x_{i-1}$ and some other node $x^{\prime}_{i-1}$. Using induction and Lemma 8, we have that

[TABLE]

where the last inequality follows from (2). The second statement follows from Lemma 7 because the endpoints of the paths form an independent set in $\mathcal{W}$. ∎

Ascending path decomposition.

An only-child-path in a binary tree is an ascending path starting in a leaf and ending at the first encountered node that has a sibling, or at the root of the tree. An only-child-path can have length $1$, if the starting leaf has a sibling already. Examples of only-child-paths are shown in Figure 2. We observe that no node with two children lies on any only-child-path, which implies that the set of only-child-paths forms an independent set of ascending paths.

Consider the following pruning procedure for a full binary forest $\mathcal{W}$. Set $\mathcal{W}_{(0)}\leftarrow\mathcal{W}$. In iteration $i$, we obtain the forest $\mathcal{W}_{(i)}$ from $\mathcal{W}_{(i-1)}$ by deleting the only-child-paths of $\mathcal{W}_{(i-1)}$. We stop when $\mathcal{W}_{(i)}$ is empty; this happens eventually because at least the leaves of $\mathcal{W}_{(i-1)}$ are deleted in the $i$-th iteration. Because we start with a full forest, the only-child-paths in the first iteration are all of length $1$, and consequently $\mathcal{W}_{(1)}$ arises from $\mathcal{W}_{(0)}$ by deleting the leaves of $\mathcal{W}_{(0)}$. Note that the intermediate forests $\mathcal{W}_{(1)},\mathcal{W}_{(2)},\ldots$ are not full forests in general. Figure 2 shows the pruning procedure on an example.

To analyze this pruning procedure in detail, we define the following integer valued function for nodes in $\mathcal{W}$:

[TABLE]

Lemma 10**.**

A node $x$ of a full binary forest $\mathcal{W}$ is deleted in the pruning procedure during the $r(x)$-th iteration.

Proof.

We prove the claim by induction on the tree structure. If $x$ is a leaf, it is removed in the first iteration, and $r(x)=1$. If $x$ is an internal node with children $y_{1}$ and $y_{2}$, let $r_{1}:=r(y_{1})$ and $r_{2}:=r(y_{2})$. By induction, $y_{1}$ is deleted in the $r_{1}$-th iteration and $y_{2}$ is deleted in the $r_{2}$-th iteration. There are two cases: if $r_{1}=r_{2}$, $x$ still has two children after $r_{1}-1$ iterations. This implies that $x$ does not lie on an only-child- path in the forest $\mathcal{W}_{(r_{1}-1)}$ and is therefore not deleted in the $r_{1}$-th iteration. But because its children are deleted, $x$ is a leaf in $\mathcal{W}_{(r_{1})}$ and therefore deleted in the $(r_{1}+1)$-th iteration. It remains the second case that $r_{1}\neq r_{2}$. Assume without loss of generality that $r_{1}>r_{2}$, so that $r(x)=r_{1}$. In iteration $r_{1}$, $y_{1}$ lies on an only-child-path in $\mathcal{W}_{(r_{1}-1)}$. Because $y_{2}\notin\mathcal{W}_{(r_{1}-1)}$, $y_{1}$ has no sibling in $\mathcal{W}_{(r_{1}-1)}$, so the only-child-path extends to $x$. Consequently, $x$ is deleted in the $r_{1}$-th iteration. ∎

Lemma 11**.**

For a node $x$ in a full binary forest, let $s(x)$ denote the number of nodes in the subtree rooted at $x$. Then $s(x)\geq 2^{r(x)}-1$. In particular, $r(x)\leq\log_{2}(s(x)+1)$.

Proof.

We prove the claim by induction on $r(x)$. Note that $r(x)=1$ if and only if $v$ is a leaf, which implies the statement for $r(x)=1$. For $r(x)>1$, it is sufficient to prove the statement assuming that $v$ has a minimal $s$-value among all nodes with the same $r$-value. Since $x$ is not a leaf, it has two children $y_{1}$ and $y_{2}$. They satisfy $r(y_{1})=r(y_{2})=r(x)-1$, because otherwise $r(x)=r(y_{1})$ or $r(x)=r(y_{2})$, contradicting the minimality of $x$. By induction hypothesis, we obtain that

[TABLE]

∎

Proposition 12**.**

$|\hat{\mathbb{K}}_{m}|\leq n+2\cdot(\Delta+1)\cdot n\cdot(1+\log_{2}(n_{0}))=O(n\cdot\Delta\cdot\log_{2}(n_{0}))$, where $n_{0}$ is the number of vertices included in $\mathcal{T}$.

Proof.

Applying the pruning procedure to the contraction forest $\mathcal{W}$ of $\mathcal{T}$, we obtain in every iteration a set of independent ascending paths of $\mathcal{W}$, and the cost of this set is bounded by $2\cdot(\Delta+1)\cdot n$ with Lemma 9. Because $\mathcal{W}$ has at most $2\cdot n_{0}-1$ nodes, any node $x$ satisfies $r(x)\leq\log_{2}(2\cdot n_{0})$ by Lemma 11. It follows that the pruning procedure ends after $1+\log_{2}(n_{0})$ iterations, so the total cost of the contraction forest is at most $2\cdot(\Delta+1)\cdot n\cdot(1+\log_{2}(n_{0}))$. By definition, this cost is equal to the number of simplices added in all contraction steps. Together with the $n$ simplices added in inclusion steps, the bound follows. ∎

3.4 Algorithm

We will make frequent use of the following concept: a dictionary is a data structure that stores a set of items of the form (k,v), where k is called the key and v is called the value of the item. We assume that all keys stored in the dictionary are pairwise different. The dictionary support three operations: insert(k,v) adds a new item to the dictionary, delete(k) removes the item with key k from the dictionary (if it exists), and search(k) returns the item with key k, or returns that no such item exists. Common realizations are balanced binary search trees [13, §12] and hash tables [13, §11].

Simplicial complexes by dictionaries.

The main data structure of our algorithm is a dictionary $D$ that represents a simplicial complex. Every item stored in the dictionary represents a simplex, whose key is the list of its vertices. Every simplex $\sigma$ itself stores an dictionary $CoF_{\sigma}$. Every item in $CoF_{\sigma}$ is a pointer to another item in $D$, representing a cofacet $\tau$ of $\sigma$. The key of the item is a vertex identifier (e.g., an integer) for $v$ such that $\tau=v\ast\sigma$.

How large is $D$ for a simplicial complex $\mathbb{K}$ with $n$ simplices of dimension $\Delta$? Observe that $D$ stores $n$ items, and each key is of size $\leq\Delta+1$. That yields a size of $O(n\Delta)$ plus the size of all $CoF_{\sigma}$. Since every simplex is the cofacet of at most $\Delta$ simplices, the size of all these inner dictionaries is also bounded by $O(n\Delta)$ (assuming that the size of a dictionary is linear in the number of stored elements).

We can insert and delete simplices efficiently in $D$ using dictionary operations. For instance, to insert a simplex $\sigma$ given as a list of vertices, we insert a new item in $D$ with the key. Then we search for the $\Delta$ facets of $\sigma$ (using dictionary search in $D$), and notify each facet $\tau$ of $\sigma$ about the insertion by adding a pointer to $\sigma$ to $CoF_{\tau}$, using the vertex $\sigma\setminus\tau$ as the key. The deletion procedure works similarly. Each simplex insertion and deletion requires $O(\Delta)$ dictionary operations. In what follows, it is convenient to assume that dictionary operations have unit costs; we multiply the time complexity with the cost of a dictionary operation at the end to compensate for this simplification.

The conversion algorithm.

We assume that the tower $\mathcal{T}$ is given to us as a stream where each element represents a simplicial map $\phi_{i}$ in the tower. Specifically, an element starts with a token {INCLUSION, CONTRACTION} that specifies the type of the map. In the first case, the token is followed by a non-empty list of vertex identifiers specifying the vertices of the simplex to be added. In the second case, the token is followed by two vertex identifiers $u$ and $v$, specifying a contraction of type $(u,v)\leadsto u$.

The algorithm outputs a stream of simplices specifying the filtration $\mathcal{F}$. Specifically, while handling the $i$-th input element, it outputs the simplices of $\hat{\mathbb{K}}_{i+1}\setminus\hat{\mathbb{K}}_{i}$ in increasing order of dimension (to ensure that every prefix is a simplicial complex). For simplicity, we assume that output simplices are also specified by a list of vertices – the algorithm can easily be adapted to return the boundary matrix of the filtration in sparse list representation with the same complexity bounds.

We use an initially empty dictionary $D$ as above, and maintain the invariant that after the $i$-th iteration, $D$ represents the active subcomplex of $\hat{\mathbb{K}}_{i}$, which is equal to $\mathbb{K}_{i}$ by Lemma 4.

If the algorithm reads an inclusion of a simplex $\sigma$ from the stream, it simply adds $\sigma$ to $D$ (maintaining the invariant) and writes $\sigma$ to the output stream.

If the algorithm reads a contraction of two vertices $u$ and $v$, from $\mathbb{K}_{i}$ to $\mathbb{K}_{i+1}$, we let $c_{i}$ denote the cost of the contraction, that is, $c_{i}=|\hat{\mathbb{K}}_{i+1}\setminus\hat{\mathbb{K}}_{i}|$. The first step is to determine which of the vertices has the smaller active closed star in $\hat{\mathbb{K}}_{i}$, or equivalently, which vertex has the smaller closed star in $\mathbb{K}_{i}$. The size of the closed star of a vertex $v$ could be computed by a simple graph traversal in $D$, starting at a vertex $v$ and following the cofacet pointers recursively, counting the number of simplices encountered. However, we want to identify the smaller star with only $O(c_{i})$ operations, and the closed star can be much larger. Therefore, we change the traversal in several ways:

First of all, observe that $|\overline{\mathrm{St}}(u)|\leq|\overline{\mathrm{St}}(v)|$ if and only if $|\mathrm{St}(u)|\leq|\mathrm{St}(v)|$ (in $\mathbb{K}_{i}$). Now define $\mathrm{St}(u,\neg v):=\mathrm{St}(u)\setminus\mathrm{St}(v)$. Then, $|\mathrm{St}(u)|\leq|\mathrm{St}(v)|$ if and only if $|\mathrm{St}(u,\neg v)|\leq|\mathrm{St}(v,\neg u)|$, because we subtracted the intersection of the stars on both sides. Finally, note that $\min\{|\mathrm{St}(u,\neg v)|,|\mathrm{St}(v,\neg u)|\}\leq c_{i}$, as one can easily verify. Moreover, we can count the size of $\mathrm{St}(u,\neg v)$ by a cofacet traversal from $u$, ignoring cofacets that contain $v$ (using the key of $CoF_{\ast}$), in $O(|\mathrm{St}(u,\neg v)|)$ time. However, this is still not enough, because counting both sets independently gives a running time of $\max\{|\mathrm{St}(u,\neg v)|,|\mathrm{St}(v,\neg u)|\}$. The last trick is that we count the sizes of $\mathrm{St}(u,\neg v)$ and $\mathrm{St}(v,\neg u)$ at the same time by a simultaneous graph traversal of both, terminating as soon as one of the traversal stops. The running time is then proportional to $2\cdot\min\{|\mathrm{St}(u,\neg v)|,|\mathrm{St}(v,\neg u)|\}=O(c_{i})$, as required.

Assume w.l.o.g. that $|\overline{\mathrm{St}}(u)|\leq|\overline{\mathrm{St}}(v)|$. Also in time $O(c_{i})$, we can obtain $\mathrm{St}(u,\neg v)$. We sort its elements by increasing dimension, which can be done in $O(c_{i}+\Delta)$ using integer sort. For each simplex $\sigma=\{u,v_{1},\ldots,v_{k}\}\in\mathrm{St}(u,\neg v)$ in order, we check whether $\{v,v_{1},\ldots,v_{k}\}$ is in $D$. If not, we add it to $D$ and also write it to the output stream. Then, we write $\{u,v,v_{1},\ldots,v_{k}\}$ to the output stream (note that we do not have to check it existence in $\hat{\mathbb{K}}_{i}$, because it does not by construction, and there is no need to add it to $D$ because of the next step). At the end of the loop, we wrote exactly the simplices in $\mathbb{K}_{i+1}\setminus\mathbb{K}_{i}$ to the output stream, which proves correctness.

It remains to maintain the invariant on $D$. Assuming still that $|\overline{\mathrm{St}}(u)|\leq|\overline{\mathrm{St}}(v)|$, $u$ turns inactive in $\hat{\mathbb{K}}_{i+1}$. We simply traverse over all cofaces of $u$ and remove all encountered simplices from $D$. After this operations, the invariant holds. This ends the description of the algorithm.

Complexity analysis.

By applying the operation costs on the above described algorithm, we obtain the following statement. Combined with Propositions 6 and 12, it completes the proof of Theorem 2.

Proposition 13**.**

The algorithm requires $O(\Delta\cdot\omega)$ space and $O(\Delta\cdot|\hat{\mathbb{K}}_{m}|\cdot C_{\omega})$ time, where $\omega=\max_{i=0,\ldots,m}|\mathbb{K}_{i}|$ and $C_{\omega}$ is the cost of an operation in a dictionary with at most $\omega$ elements.

Proof.

The space complexity follows at once from the invariant, because the size of $D$ is at most $O(\Delta|\mathbb{K}_{i}|)$ during the $i$-th iteration.

For the time bound, we set $S:=|\hat{\mathbb{K}}_{m}|$ for convenience and show that the algorithm finishes in $O(\Delta\cdot S)$ steps, assuming dictionary operations to be of constant cost. We have one simplex insertion per elementary inclusion which requires $O(\Delta)$ operations. Thus, all inclusions are bounded by $O(n\Delta)$, which is subsumed by our bound as $n\leq S$. For the contraction case, we need $O(c_{i})$ to identify the smaller star, $O(c_{i}+\Delta)$ to get a sorted list of $\mathrm{St}(u,\neg v)$ (or vice versa), and $O(\Delta\cdot c_{i})$ to add new vertices to $D$. Moreover, we delete the star of $u$ from $D$; the cost for that is $O(\Delta\cdot d_{i})$, where $d_{i}$ is the number of deleted simplices. Thus, the complexity of a contraction is $O(\Delta\cdot(c_{i}+d_{i}))$.

Since $c_{i}$ is the number of simplices added to the filtration at step $i$, the sum of all $c_{i}$ is bounded by $O(S)$. Moreover, because every simplex that ever appears in $D$ belongs to $\hat{\mathbb{K}}_{m}$ and every simplex is inserted only once, the sum of all $d_{i}$ is bounded by $O(S)$ as well. Therefore, the combined cost over all contractions is $O(\Delta\cdot S)$ as required. ∎

Note that the dictionary $D$ has lists of identifiers of length up to $\Delta+1$ as keys, so that comparing two keys has cost $O(\Delta)$. Therefore, using balanced binary trees as dictionary, we get a complexity of $C_{\omega}=O(\Delta\log\omega)$. Using hash tables, we get an expected worst-case complexity of $C_{\omega}=O(\Delta)$.

Experimental results.

The following tests where made on a 64-bit Linux (Ubuntu) HP machine with a 3.50 GHz Intel processor and 63 GB RAM. The programs were all implemented in C++ and compiled with optimization level -O2. Our algorithm was implemented in the software Sophia.444https://bitbucket.org/schreiberh/sophia/

To test its performance, we compared it to the software Simpers555http://web.cse.ohio-state.edu/~tamaldey/SimPers/Simpers.html (downloaded in August 2017), which is the implementation of the Annotation Algorithm from Dey, Fan and Wang described in [15]. Simpers computes the persistence of the given filtration, so we add to our time the time the library PHAT (version 1.5) needs to compute the persistence from the generated filtration. PHAT was used with its default parameters and its ’–ascii –verbose’ options activated. Simpers also used its default parameters except for the dimension parameter which was set to 5.

The results of the tests are in Table 1. The timings for File IO are not included in any process time except the input reading of Sophia. The memory peak was obtained via the ’/usr/bin/time -v’ Linux command. Each command was repeated 10 times and the average was token. The first three towers in the table, data1-3, were generated incrementally on a set of $n_{0}$ vertices: In each iteration, with $90\%$ probability, a new simplex is included, that is picked uniformly at random among the simplices whose facets are all present in the complex, and with $10\%$ probability, two randomly chosen vertices of the complex are contracted. This is repeated until the complex on the remaining $k$ vertices forms a $k-1$-simplex, in which case no further simplex can be added. The remaining data was generated from the SimBa (downloaded in June 2016) library with default parameters using the point clouds from [16]. To obtain the towers that SimBa computes internally, we included a print command at a suitable location in the SimBa code.

To verify that the space consumption of our algorithm does not dependent on the length of the tower, we constructed an additional example whose size exceeds our RAM capacity, but whose width is small: we took 10 random points moving on a flat torus in a randomly chosen fixed direction. When two points get in distance less than $t_{1}$ to each other, we add the edge between them (the edge remains also if the points increase their distance later on). When two points get in distance less than $t_{2}$ from each other with $t_{2}<t_{1}$, we contract the two vertices and let a new moving point appear somewhere on the torus. This process yields a sequence of graphs, and we take its flag complex as our simplicial tower. In this way, we obtain a tower of length about $3.5\cdot 10^{9}$ which has a file size of about $73$ GB, but only has a width of $367$. Our algorithm took about $2$ hours to convert this tower into a filtration of size roughly $4.6\cdot 10^{9}$. During the conversion, the virtual memory used was constantly around $22$ MB and the resident set size about $3.8$ MB only, confirming the theoretical prediction that the space consumption is independent of the length of the tower. The information about the memory use was taken from the ’/proc/<pid>/stat’ system file every 100 000 insertions/contractions during the process.

3.5 Tightness and lower bounds

The conversion theorem (Theorem 2) yields an upper bound of $O(\Delta\cdot n\log n_{0})$ for the size of a filtration equivalent to a given tower. It is natural to ask whether this bound can be improved. In this section, we assume for simplicity that the maximal dimension $\Delta$ is a constant. In that case, it is not difficult to show that our size analysis cannot be improved:

Proposition 14**.**

There exist an example of a tower with $n$ simplices and $n_{0}$ vertices such that our construction yields a filtration of size $\Omega(n\log n_{0})$.

Proof.

Let $p=2^{k}$ for some $k\in\mathbb{N}$. Consider a graph with $p$ edges $(a_{i},b_{i})$, with $a_{1},\ldots,a_{p}$, $b_{1},\ldots,b_{p}$ $2p$ distinct vertices. Our tower first constructs this graph with inclusions (in an arbitrary order). Then, the $a$-vertices are contracted in a way such that the contracting forest is a fully balanced binary tree – see Figure 3 for an illustration.

To bound the costs, it suffices to count the number of edges added between $a$-vertices and $b$-vertices in each step. We call them $ab$-edges from now on. Define the level of a contraction to be its level in the contracting tree, where [math] is the level of the leaves, and $k$ is the level of the root. On level $1$, a contraction of $a_{i}$ and $a_{j}$ yields exactly one new $ab$-edge, either $(a_{i},b_{j})$ or $(a_{j},b_{i})$. The resulting contracted vertex has two incident $ab$-edges. A contraction on level $2$ yields two novel $ab$-edges, and a vertex with $4$ incident $ab$-edges. By a simple induction, we observe that a contraction on level $i$ introduces $2^{i-1}$ new $ab$-edges, and hence has a cost of at least $2^{i-1}$, for $i=1,\ldots,k$. This means that the sum of the costs of all level $i$ contractions is exactly $\frac{p}{2}$. Summing up over all $i$ yields a cost of at least $k\frac{p}{2}$. The result follows because $k=\log(n_{0}/2)$ and $p=n/3$. ∎

Another question is whether a different approach could convert a tower into a filtration with an (asymptotically) smaller number of simplices. For this question, consider the inverse persistence computation problem: given a barcode, find a filtration of minimal size which realizes this barcode. Note that solving this problem results in a simple solution for the conversion problem: compute the barcode of the tower first using an arbitrary algorithm; then compute a filtration realizing this barcode. While useful for lower bound constructions, we emphasize the impracticality of this solution, as the main purpose of the conversion is a faster computation of the barcode.

Let $b$ be the number of bars of a barcode. An obvious lower bound for the filtration size is $\Omega(b)$, because adding a simplex causes either the birth or the death of exactly one bar in the barcode. For constant dimension, $O(b)$ is also an upper bound:

Lemma 15**.**

For a barcode with $b$ bars and maximal dimension $\Delta$, there exists a filtration of size $\leq 2^{\Delta+2}b$ realizing this barcode.

Proof.

Begin with an empty complex. Now consider the birth and death times represented by the barcode one by one. The first birth will be the one of a [math]-dimensional class, so add a vertex $v_{0}$ to the complex. From now, every time a [math]-dimensional homology class is born add a new vertex to the complex. When a [math]-dimensional homology class dies, link the corresponding vertex to $v_{0}$ with an edge.

When a $k$-dimensional homology class is born, with $k>0$, add the boundary of a $(k+1)$-simplex to the complex that is incident to $v_{0}$ and to $k$ novel vertices. When this homology class dies, add the corresponding $(k+1)$-simplex. This way, the resulting filtration realizes the barcode. For a bar of the barcode in dimension $k$, we have to add all proper faces of a $(k+1)$-simplex (except for one vertex). Since that number is at most $2^{k+2}-3$, the result follows. ∎

If a tower has length $m$, what is the maximal size of its barcode? If the size of the barcode is $O(m)$, the preceding lemma implies that a conversion to a filtration of linear size is possible (for constant dimension). On the other hand, any example of a tower yielding a super-linear lower bound would imply a lower bound for any conversion algorithm, because a filtration has to contain at least one simplex per bar in its barcode.

To approach the question, observe that a single contraction might destroy many homology classes at once: consider a “fan” of $t$ empty triangles, all glued together along a common edge $ab$ (see Figure 4 (a)). Clearly, the complex has $t$ generators in $1$-homology. When $a$ and $b$ are contracted, the complex transforms to a star-shaped graph which is contractible. Moreover, a contraction can also create many homology classes at once: consider a collection of $t$ disjoint triangulated spheres, all glued together along a common edge $ab$. For every sphere, remove one of the triangles incident to $ab$ (see Figure 4 (b)). The resulting complex is acyclic. The contraction of the edge $ab$ “closes” each of the spheres, and the resulting complex has $t$ generators in $2$-homology. Finally, a contraction might not not affect the homology at all – see Figure 4 (c).

The above examples show that a single contraction can create and destroy many bars. For a super-linear bound on the barcode size, however, we would have to construct an example where sufficiently many contractions create a large number of bars. So far, we did neither succeed in constructing such an example, nor are we able to show that such an example does not exist.

4 Persistence by Streaming

Even if the complex we need to maintain in memory during the conversion is relatively small, at the end, the algorithm to compute the final persistence still needs to memorize the whole filtration we give it as input. If the complex in the original tower has an interesting maximum size during the whole process, we should be able to compute its persistence even if the tower is extremely long. So we design here a streaming variation of the reduction algorithm that computes the barcode of filtrations, such that it has an efficient memory use. More precisely, we will prove the following theorem:

Theorem 16**.**

With the same notation as in Theorem 16, we can compute the barcode of a tower $\mathcal{T}$ in worst-case time $O(\omega^{2}\cdot\Delta\cdot n\cdot\log n_{0})$ and space complexity $O(\omega^{2})$

We describe the algorithm in Section 4.1 and prove the complexity bounds in Section 4.2. The described algorithm requires various adaptations to become efficient in practice, and we describe an improved variant in Section 4.3. We present some experiments in Section 4.4.

4.1 Algorithmic description

On a high level, our algorithm converts the tower into an equivalent filtration and computes the barcode of that filtration. We focus on the second part, which we describe as a streaming algorithm. The input to the algorithm is a stream of elements, each starting with a token {ADDITION, INACTIVE} followed by a simplex identifier which represents a simplex $\sigma$. In the addition case, this is followed by a list of simplex identifiers specifying the facets of $\sigma$. In other words, the element encodes the next column of the boundary matrix. For the inactive case, it means that $\sigma$ has become inactive in the complex. In particular, no subsequent simplex in the stream will contain $\sigma$ as a facet. It is not difficult to modify the algorithm from Section 3.4 to return a stream as required, within the same complexity bounds.

The algorithm uses a matrix data type $M$ as its main data structure. We realize $M$ as a dictionary of columns, indexed by a simplex identifier. Each column is a sorted linked list of identifiers corresponding to the non-zero row indices of the column. In particular, we can access the pivot of the column in constant time and we can add two columns in time proportional to the maximal size of the involved lists. Note that most algorithms store the boundary matrix as an array of columns, but we use dictionaries for space efficiency.

There are two secondary data structures that we mention briefly: given a row index $r$, we have to identify the column index $c$ of the column that has $r$ as pivot in the matrix (or to find out that no such column exists). This can be done using a dictionary with key and value type both simplex identifiers. Finally, we maintain a dictionary representing the set of simplex identifiers that represent active simplices of the complex, plus a flag denoting whether the corresponding simplex is positive or negative. It is straight-forward to maintain these structures during the algorithm, and we will omit the description of the required operations.

The algorithm uses two subroutines. The first one, called reduce_column, takes a column identifier $j$ as input is defined as follows: iterate through the non-zero row indices of $j$. If an index $i$ is the index of an inactive and negative column in $M$, remove the entry from the column $j$ (this is the “compression” described at the end of Section 2). After this pre-processing, reduce the column: while the column is non-empty, and its pivot $i$ is the pivot of another column $k<j$ in the matrix, add column $k$ to column $j$.

The second subroutine, remove_row, takes a index $\ell$ as input and clears out all entries in row $\ell$ from the matrix. For that, let $j$ be the column with $\ell$ as pivot. Traverse all non-zero columns of the matrix except column $j$. If a column $i\neq j$ has a non-zero entry at row $\ell$, add column $j$ to column $i$. After traversing all columns, remove column $j$ from $M$.

The main algorithm can be described easily now: if the input stream contains an addition of a simplex, we add the column to $M$ and call reduce_column on it. If at the end of that routine, the column is empty, it is removed from $M$. If the column is not empty and has pivot $\ell$, we report $(\ell,j)$ as a persistence pair and check whether $\ell$ is active. If not, we call remove_row($\ell$). If the input stream specifies that simplex $\ell$ becomes inactive, we check whether $j$ appears as pivot in the matrix and call remove_row($\ell$) in this case.

Proposition 17**.**

The algorithm computes the correct barcode.

Proof.

First, note that removing a column from $M$ within the procedure remove_row does not affect further reduction steps: Since the pivot $\ell$ of the column is inactive, no subsequent column in the stream will have an entry in row $\ell$. Moreover, the reduction process cannot introduce an entry in row $\ell$ because the routine has removed all such entries.

Note that remove_row might also include right-to-left column additions, and we also have to argue that they do not change the pivots. Let $R$ denote the matrix obtained by the standard compression algorithm that does not call remove_row (as described in Section 2). Just before our algorithm calls reduce_column on the $j$-th column, let $S$ denote the matrix with $(j-1)$ columns that is represented by $M$. It is straight-forward to verify by an inductive argument that every column of $S$ is a linear combination of $R_{1},\ldots,R_{j-1}$. reduce_column adds a subset of the columns of $S$ to the $j$-th column. Thus, the reduced column can be expressed by a sequence of left-to-right column additions in $R$, and thus yields the same pivot as the standard compression algorithm. ∎

4.2 Complexity analysis

We analyze how large the structure $M$ can become during the algorithm. After every iteration, the matrix represents the reduced boundary matrix of some intermediate complex $\hat{\mathbb{L}}$ with $\hat{\mathbb{K}}_{i}\subseteq\hat{\mathbb{L}}\subseteq\hat{\mathbb{K}}_{i+1}$ for some $i=0,\ldots,m$. Moreover, the active simplices define a subcomplex $\mathbb{L}\subseteq\hat{\mathbb{L}}$ and there is a moment during the algorithm where $\hat{\mathbb{L}}=\hat{\mathbb{K}}_{i}$ and $\mathbb{L}=\mathbb{K}_{i}$, for every $i=0,\ldots,m$. We call this the $i$-th checkpoint. We will make frequent use of the following simple observation.

Lemma 18**.**

$|\hat{\mathbb{K}}_{i+1}\setminus\hat{\mathbb{K}}_{i}|\leq|\mathbb{K}_{i}|\leq\omega$**

Lemma 19**.**

At every moment, the number of columns stored in $M$ is at most $2\omega$.

Proof.

It can be verified easily that throughout the algorithm, a column is stored in $M$ only if not zero, and its pivot is active. So, assume first that we are at the $i$-th checkpoint for some $i$. Since $\mathbb{L}=\mathbb{K}_{i}$, there are not more than $|\mathbb{K}_{i}|\leq\omega$ active simplices. Since each column has a different active pivot, their number is also bounded by $\omega$. If we are between checkpoint $i$ and $i+1$, there have been not more than $\omega$ columns added to $M$ since the $i$-th checkpoint from Lemma 18. The bound follows. ∎

The number of rows is more difficult to bound because we cannot guarantee that each column in $M$ corresponds to an active simplex. Still, the number of rows is asymptotically the same as for columns:

Lemma 20**.**

At every moment, the number of rows stored in $M$ is at most $4\omega$.

Proof.

Consider a row index $\ell$ and a time in the algorithm where $M$ represents $\hat{\mathbb{L}}$. By the same argument as in the previous lemma, we can argue that there are at most $2\omega$ active row indices at any time. Therefore, we restrict our attention to the case that $\ell$ is inactive and distinguish three cases. If $\ell$ represents a negative simplex, we observe that its row should have been was removed due to the compression optimization, and after $\ell$ became inactive, no simplex can have it as a facet either. It follows that row $\ell$ is empty in this case. If $\ell$ is positive and was paired with another index $j$ during the algorithm, then remove_row was called on $\ell$, either at the moment the pair was formed, or when $\ell$ became inactive. Since the procedure removes the row, we can conclude that row $\ell$ is empty also in this case. The final case is that $\ell$ is positive, but has not been paired so far. It is well-known that in this case, $\ell$ is the generator of an homology class of $\hat{\mathbb{L}}$. Let

[TABLE]

denote the sum of the Betti numbers of the complex. Then, it follows that the number of such row indices is at most $\beta(\hat{\mathbb{L}})$.

We argue that $\beta(\hat{\mathbb{L}})\leq 2\omega$ which proves our claim. Assume that $\hat{\mathbb{K}}_{i}\subseteq\hat{\mathbb{L}}\subset\hat{\mathbb{K}}_{i+1}$. We have that $\beta(\hat{\mathbb{K}}_{i})=\beta(\mathbb{K}_{i})$ by Lemma 5, and since $\mathbb{K}_{i}$ has at most $\omega$ simplices, $\beta(\mathbb{K}_{i})\leq\omega$. Since we add at most $\omega$ simplices to get from $\hat{\mathbb{K}}_{i}$ to $\hat{\mathbb{L}}$, and each addition can increase $\beta$ by at most one, we have indeed that $\beta(\hat{\mathbb{L}})\leq 2\omega$. ∎

Proposition 21**.**

The algorithm has time complexity $O(\omega^{2}\cdot\Delta\cdot n\cdot\log n_{0})$ and space complexity $O(\omega^{2})$.

Proof.

The space complexity is immediately clear from the preceding two lemmas, as $M$ is the dominant data structure in terms of space consumption. For the time complexity, we observe that both subroutines reduce_column and remove_row need $O(\omega)$ column additions and $O(\omega)$ dictionary operations in the worst case. A column addition costs $O(\omega)$, and a dictionary operation is not more expensive (since the dictionaries contain at most $O(\omega)$ elements and their keys are integers). So, the complexity of both methods is $O(\omega^{2})$. Since each routine is called at most once per input element, and there are $O(\Delta\cdot n\cdot\log n_{0})$ elements by Theorem 2, the bound follows. ∎

4.3 Implementation

The algorithm described in Section 4.1 is not efficient in practice for three reasons: First, it has been observed as a general rule that the clearing optimization [3] significantly improves the standard algorithm. However, in the above form, that optimization is not usable because it requires to process the columns in a non-incremental way. Second, the remove_row routine scans the entire matrix $M$; while not affecting the worst-case complexity, frequently scanning the matrix should be avoided in practice. Finally, the above algorithm uses lists to represent columns, but it has been observed that this is a rather inefficient way to perform matrix operations [3].

We outline a variant of the above algorithm that partially overcomes these drawbacks and behaves better in practice. In particular, our variant can be implemented with all column representations available in the Phat library. The idea is to perform a “batch” variant of the previous algorithm: We define a chunk size $C$ and read in $C$ elements from the stream; we insert added columns in the matrix, removing row entries of already known inactive negative columns as before, but not reducing the columns yet. After having read $C$ elements, we start the reduction of the newly inserted columns using the clearing optimization. That is, we go in decreasing dimension and remove a column as soon as its index becomes the pivot of another column; see [10] for details. After the reduction ends, except for the last chunk, we go over the columns of the matrix and check for each pivot whether it is active. If it is, we traverse its row entries in decreasing order, but skipping the pivot. Let $\ell$ be the current entry. If $\ell$ is the inactive pivot of some column $j$, we add $j$ to the current column. If $\ell$ is inactive and represents a negative column, we delete $\ell$ from the current column. After performing these steps for all remaining columns of the matrix, we go over all columns again, deleting every column with inactive pivot. (By the way, also cleaning up the secondary data structures described in 4.1.)

It remains the question of how to choose the parameter $C$. The chunk provides a trade-off between time and space efficiency. Roughly speaking, the matrix can have up to $O(\omega+C)$ columns during this reduction, but the larger the chunks are, the more benefit one can draw from the clearing optimization (the clearing optimization fails for pairs where the simplices are in different chunks). We therefore recommend to choose $C$ rather large, but making sure that the matrix will still fit into memory.

4.4 Experimental evaluation

The tests were made with the same setup as in Section 3.4. Figure 5 shows the effect of the chunk size parameter $C$ on the runtime and memory consumption of the algorithm. The data used is S3 (see Section 3.4); we also performed the tests on the other examples from Table 1, with similar outcome. The File IO operations are included in the measurements. Confirming the theory, as the chunk size decreases, our implementation needs less space but more computation time (while the running time seems to increase slightly again for larger chunk sizes).

For the $4.6\cdot 10^{9}$ inclusions tower from Section 3.4, with $C=200\,000$, the algorithm took around 4.5 hours, the virtual memory used was constantly around $68$ MB and the resident set size constantly around $49$ MB, confirming the theoretical statement that the memory size does not depend on the length of the filtration.

5 Conclusion

In the first part of the paper, we have presented an efficient algorithm to reduce the computation of the barcode of a simplicial tower to the computation of a barcode of a filtration with slightly larger size. With our approach, every algorithmic improvement for persistence computation on the filtration case becomes immediately applicable to the case of towers as well. In the second part of the paper, we present a streaming variant of the classical persistence algorithm for the case of towers. In here, the extra information provided by the towers allows a more space efficient storage of the boundary matrix.

There are various theoretical and practical questions remaining for further work: as already exposed in Section 3.5, the question of how large can the barcode of a tower become has immediate consequences on the conversion from towers to filtrations. Moreover, while we focused on constant dimension in Section 3.5, we cannot exclude the possibility that our algorithm achieves a better asymptotic bound for non-constant dimensions.

We made our software publicly available in the Sophia library. There are several open questions regarding practical performance. For instance, our complex representation, based on hash table, could be replaced with other variants, such as the Simplex tree [5]. Moreover, it would be interesting to compare our streaming approach for the barcode computation with a version that converts to a filtration and subsequently computes the barcode with the annotation algorithm. The reason is that the latter algorithm only maintains a cohomology basis of the currently active complex in memory and therefore avoids the storage of the entire boundary matrix.

Since both the simplex tree and annotation algorithm are part of the Gudhi library [32], we plan to integrate our conversion algorithm in an upcoming version of the Gudhi library, both to increase the usability of our software and to facilitate the aforementioned comparisons.

Acknowledgements

The authors are supported by the Austrian Science Fund (FWF) grant number P 29984-N35.