Barcodes of Towers and a Streaming Algorithm for Persistent Homology

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. Experimental evaluations show that our approach can efficiently handle towers with billions of complexes.


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 (K i ) i=0,...,m and simplicial maps φ i : K i → 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 "unionby-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(∆ • n • log(n 0 )), where ∆ 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 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 ω keeping only up to O(ω) columns of the boundary matrix in memory.This yields a space complexity of O(ω 2 ) and a time complexity of O((∆ • n • log(n 0 ))ω 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 available1 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 co-homology, 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 φ i : ).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 library2 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.

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 (0-simplices) which we denote by u, v, and w.For a k-simplex σ, a simplex τ is a face of σ if τ ⊆ σ.If τ is of dimension , we call it a -face.If < k, we call τ a proper face of σ, and if = k − 1, we call it a facet.For a simplex σ and a vertex v / ∈ σ, we define the join v * σ as the simplex {v} ∪ σ.These definition are inspired by visualizing a k-simplex as the convex hull of k + 1 affinely independent points in R k , but we will not need this geometric interpretation in our arguments.An (abstract) simplicial complex K over V is a set of simplices that is closed under taking faces.We call V the vertex set of K and write V(K) := V .The dimension of K is the maximal dimension of its simplices.For σ, τ ∈ K, we call σ a coface of τ in K if τ is a face of σ.In this case, σ is a cofacet of τ if their dimensions differ by exactly one.A simplicial complex L is a subcomplex of K if L ⊆ K. Given W ⊆ V, the induced subcomplex by W is the set of all simplices σ in K with σ ⊆ W. For a subcomplex L ⊆ K and a vertex v ∈ V(K) \ V(L), we define the join v * L := {v * σ | σ ∈ L}.For a vertex v ∈ K, the star of v in K, denoted by St(v, K), is the set of all cofaces of v in 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 St(v, K).Equivalently, the closed star is the smallest subcomplex of K containing the star of v.The link of v, Lk(v, K), is defined as St(v, K) \ St(v, K).It can be checked that the link is a subcomplex of K.When the complex is clear from context, we will sometimes omit the K in the notation of stars and links.
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 ∈ L, we call φ an elementary inclusion.The simplest example of a non-inclusion simplicial map is and φ is the identity on all remaining vertices of K. We call φ an elementary contraction and write (u, v) u as a shortcut.These notions were introduced by Dey, Fan and Wang in [15] and they also showed that any simplicial map K φ → L can be written as the composition of elementary contractions3 and inclusions.
A tower of length m is a collection of simplicial complexes K 0 , . . ., K m and simplicial maps φ i : K i → K i+1 for i = 0, . . ., m−1.From this initial data, we obtain simplicial maps φ i,j : K i → K j for i ≤ j by composition, where φ i,i is simply the identity map on 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 φ i are inclusion maps.The dimension of a tower is the maximal dimension among the K i , and the width of a tower is the maximal size among the K i .For filtrations, dimension and width are determined by the dimension and size of K m , but this is not necessarily true for general towers.
Homology and Collapses.For a fixed base field F, let H p (K) := H p (K, F) the p-dimensional homology group of K.It is well-known that H p (K) is a F-vector space.Moreover, a simplicial map K φ → L induces a linear map H p (K) φ * → H p (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 K, is a simplex with exactly one proper coface in K.An elementary collapse in K is the operation of removing a free face and its unique coface from K, yielding a subcomplex of K. We say that K collapses to L, if there is a sequence of elementary collapses transforming K into L.The following result is then well-known: Lemma 1.Let K be a complex that collapses into the complex L.Then, the inclusion map L φ → K induces an isomorphism φ * between H p (L) and H p (K).
Barcodes.A persistence module is a sequence vector spaces V 0 , . . ., V m and linear maps f i,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 ≤ d for the persistence module we can write every persistence module as the direct sum I b1,d1 ⊕ . . .⊕ I bs,ds , 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 ), . . ., (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 (K i ) i=0,...,m is a filtration such that K 0 = ∅ and K i+1 has exactly one more simplex than K i .We label the simplices of K m accordingly as σ 1 , . . ., σ m , with The filtration can be encoded as a boundary matrix ∂ of dimension m×m, where the (ij)-entry is 1 if σ i is a facet of σ j , and 0 otherwise.In other words, the j-th column of ∂ encodes the facets of σ j , and the i-th row of ∂ encodes the cofacets of σ i .Moreover, ∂ is upper-triangular because every K i is a simplicial complex.We will sometimes identify rows and columns in ∂ with the corresponding simplex in K m .Adding the k-simplex σ i to 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 σ i and the i-th column of ∂ 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 Z 2 , and interpret the coefficients of ∂ accordingly.In an arbitrary matrix A, a left-to-right column addition is an operation of the form A k ← A k + A with < k, where A k and A 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 ∂ does not have a unique reduction, the pivots of all its reductions are the same.Moreover, the pivots (b 1 , d 1 ), . . ., (b s , d s ) of R are precisely the barcode of the filtration.A direct consequence is that a simplex σ 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]).

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: a tower where, w.l.o.g., K 0 = ∅ and each φ i is either an elementary inclusion or an elementary contraction.Let ∆ denote the dimension and ω the width of the tower, and let n ≤ m denote the total number of elementary inclusions, and n 0 the number of vertex inclusions.Then, there exists a filtration , where the inclusions are not necessarily elementary, such that T and F have the same barcode and the width of the filtration , where C ω is the cost of an operation in a dictionary with ω elements.
The remainder of the section is organized as follows.We define F in Section 3.1 and prove that it yields the same barcode as 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 F and analyze its time and space complexity.

Active and small coning
Coning.We briefly revisit the coning strategy introduced by Dey, Fan and Wang [15].Let φ : K → L be an elementary contraction (u, v) u and define An example is shown in Figure 1.show that L ⊆ L * and that the map induced by inclusion is an isomorphism between H(L) and H(L * ).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 T , we can also obtain an non-zigzag filtration using coning, if we continue the operation on L * instead of going back to L.More precisely, we set K0 := K 0 and if φ i is an inclusion of simplex σ, we set Ki+1 := Ki ∪ {σ}.If φ i is a contraction (u, v) u, we set Ki+1 = Ki ∪ u * St(v, Ki ) .Indeed, it can be proved that ( Ki ) i=0,...,m has the same barcode as T .However, the filtration will not be small, and we will define a smaller variant now.
Our new construction yields a sequence of complexes K0 , . . ., Km with Ki ⊆ Ki+1 .During the construction, we maintain a flag for each vertex in Ki , 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 Ki , let ActSt(u, Ki ) denote its active closed star, which is the set of active simplices in Ki in the closed star of u.
and mark u as inactive.Otherwise, we set and mark v as inactive.This ends the description of the construction.We write F for the filtration ( Ki ) i=0,...,m .
There are two major changes compared to the construction of ( Ki ) i=0,...,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 Ki induced by the active vertices is isomorphic to K i .As a consequence, we add the same number of simplices by passing from Ki to Ki+1 as in the approach by Dey et al. does when passing from K to L * .
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 is the elementary contraction of u and v.It is obvious that both choices yield identical complexes K i+1 up to renaming of vertices.However, the choices make a difference in terms of the size of Ki+1 , because the active closed star of u to v in Ki might differ in size.Our construction simply choose the representative which causes the smaller Ki+1 .

Topological equivalence
We make the following simplifying assumption for T .Let φ i be an elementary contraction of u and v.If our construction of Ki+1 turns v inactive, we assume that φ i (u) = φ i (v) = u.Otherwise, we assume φ i (u) = φ i (v) = v.Indeed, this is without loss of generality because it corresponds to a renaming of the simplices in each 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 K i is equal to the set of active vertices in Ki .This allows us to interpret K i and Ki as simplicial complexes defined over a common vertex set.In fact, K i is the subcomplex of Ki induced by the active vertices: Proof.We use induction on i.The statement is true for i = 0, because K 0 = ∅ = K0 .So assume first φ i : K i → K i+1 is an elementary inclusion that adds a d-simplex σ = (v 0 , . . ., v d ) to K i+1 .If σ is a vertex, it is set active in Ki+1 by construction.Otherwise, v 0 , . . ., v d are active in Ki by induction and stay active in Ki+1 .In any case, σ is active in Ki+1 .The equivalence for the remaining simplices is straight-forward.
If φ i is an elementary contraction (u, v) u, we prove both directions of the equivalence separately.For "⇒", fix a d-simplex σ ∈ K i+1 .It suffices to show that σ ∈ Ki+1 , as in this case, it is also active by Lemma 3. If σ ∈ K i , this follows at once by induction because ∈ K i , u must be a vertex of σ.Moreover, writing σ = {u, v 1 , ..., v d } and σ = {v, v 1 , . . ., v d } we have that σ ∈ K i and φ i (σ ) = σ.In particular, the vertices v 1 , . . ., v d are active in Ki by induction, hence {v 1 , . . ., v d } is in the active closed star of v in Ki .By construction, {u, v 1 , ..., v d } = σ is in Ki+1 .
For "⇐", let σ ∈ Ki+1 \ K i+1 .We show that σ is an inactive simplex in Ki+1 .By Lemma 3, this is equivalent to show that σ contains a vertex not in K i+1 .
Case 1: σ ∈ Ki .If σ is inactive in Ki , it stays inactive in Ki+1 .So, assume that σ is active in Ki and thus σ ∈ K i by induction.But σ / ∈ K i+1 , so σ must have v as a vertex and is therefore inactive in Ki+1 .
Case 2: σ ∈ Ki+1 \ Ki .By construction, σ is of the form {u, v 1 , . . ., v d } such that {v 1 , . . ., v d } is in the active closed star of v in Ki .Assume for a contradiction that v = v j for all j = 1, . . ., d.Then, σ = {v, v 1 , . . ., v d } is active in Ki and thus, by induction, a simplex in K i .But then, φ i (σ ) = σ ∈ K i+1 which is a contradiction to our choice of σ.It follows that v is a vertex of σ which proves our claim.
Lemma 5.For every 0 ≤ i ≤ m, the complex Ki collapses to K i .
Proof.We use induction on i.For i = 0, K 0 = K0 , and the statement is trivial.Suppose that the statement holds for Ki and K i using the sequence s i of elementary collapses.Note that these collapses only concern inactive simplices in Ki .For an inactive vertex v ∈ Ki , the construction of Ki+1 ensures that v does not gain any additional coface.This implies that s i is still a sequence of elementary collapses for Ki+1 , yielding a complex K * i+1 with K i+1 ⊆ K * i+1 ⊆ Ki+1 .In particular, K * i+1 only contains vertices that are still active in Ki .If φ i is an elementary inclusion, K * i+1 = K i+1 , because any all vertices in Ki remain active in Ki+1 .For φ i being an elementary contraction (u, v) u, set S := K * i+1 \ K i+1 as the remaining set of simplices that still need to be collapsed to obtain 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 , . . ., v d } with v 1 , . . ., v d active in Ki+1 .We split S = S u ∪ S ¬u where S u ⊂ S are the simplices in S that contain u as vertex, and S ¬u = S \ S u .
We now define a sequence of elementary collapses from K * i+1 to K i+1 .Choose a simplex σ = {v, v 1 , . . ., v d } ∈ S ¬u of maximal dimension, and let τ = {u, v, v 1 , . . ., v d } denote the corresponding simplex in S u .Then σ is indeed a free face in K * i+1 , because if there was another coface τ = τ , it takes the form {w, v, v 1 , . . ., v d } with w = u active.So, τ ∈ S ¬u , and τ has larger dimension than σ, a contradiction.Therefore, the pair (σ, τ ) defines an elementary collapse in K * i+1 .We proceed with this construction, always collapsing a remaining pair in S ¬u × S u of maximal dimension, until all elements of S have been collapsed.Proposition 6. T and F have the same barcode.
Proof.Let φi : Ki → Ki+1 denote the inclusion map from Ki to Ki+1 .By combining Lemma 5 with Lemma 1, we have an isomorphism inc * i : H(K i ) → H( Ki ), for all 0 ≤ i ≤ m, induced by the inclusion maps inc i : K i → Ki , and therefore the following diagram connecting the persistence modules induced by T and F: The Persistence Equivalence Theorem [17, p.159] asserts that (K j ) j and ( Kj ) j , with j = 0, ..., m, have the same barcode if (1) Two contiguous maps are known to be homotopic [28,Theorem 12.5.]and thus equal at homology level, that is, φ * = ψ * .We show that inc i+1 • φ i and φi • inc i are contiguous.This implies that (1) commutes, because, by functoriality, To show contiguity, fix σ ∈ K i and observe that ( φi • inc i )(σ) = σ because φi and inc i are inclusions.If φ i (σ) = σ, (inc i+1 • φ i )(σ) = σ as well, and the contiguity condition is clearly satisfied.So, let φ i (σ) = σ.Then φ i is an elementary contraction (u, v) u, and σ is of the form {v, v 1 , . . ., v d }, where one of the v j might be equal to u.Then, is in the active closed star of v in Ki , and by construction {u, v, v 1 , ..., v d } ∈ Ki+1 , which proves contiguity of the maps.

Size analysis
The contracting forest.We associate a rooted labeled forest W j to a prefix ∅ = K 0 φ0 −→ . . .φj−1 − −− → K j of T inductively as follows: For j = 0, W 0 is the empty forest.Let W j−1 be the forest of K 0 → . . .→ K j−1 .If φ j−1 is an elementary inclusion of a d-simplex, we have two cases: if d > 0, set W j := W j−1 .If a vertex v is included, W j := W j−1 ∪ {x}, with x a single node tree labeled with v.If φ j−1 is an elementary contraction contracting two vertices u and v in K j−1 , there are two trees in W j−1 , whose roots are labeled u and v.In 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 W i are labeled with the vertices of the complex K i , for every i = 0, . . ., m.Moreover, each leaf corresponds to the inclusion of a vertex in K 0 → . . .→ K i , and each internal node corresponds to a contraction of two vertices.In particular, W i is a full forest, that is, every node has 0 or 2 children.
Let W := W m denote the forest of the tower T .Let Σ denote the set of all simplices that are added at elementary inclusions in T , and recall that n = |Σ|.A d-simplex σ ∈ Σ added by φ i is formed by d + 1 vertices, which correspond to d + 1 roots of W i+1 , and equivalently, to d + 1 nodes in W. For a node x in W, we denote by E(x) ⊆ Σ the subset of simplices with at least one vertex that appears as label in the subtree of 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: We say that the set N of nodes in W is independent, if there are no two nodes x 1 = x 2 in N , such that x 1 is an ancestor of x 2 in W. A vertex in K i appears as label in at most one W-subtree rooted at a vertex in the independent set N .Thus, a d-simplex σ 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 W.Then, The cost of contracting.Recall that a contraction K i φi → K i+1 yields an inclusion Ki → Ki+1 that potentially adds more than one simplex.Therefore, in order to bound the total size of Km , we need to bound the number of simplices added in all these contractions.
We define the cost of a contraction φ i as | Ki+1 \ Ki |, that is, the number of simplices added in this step.Since each contraction corresponds to a node x in W, we can associate these costs to the internal nodes in the forest, denoted by c(x).The leaves get cost 0. Proof.Let φ i : K i → 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 K i that are contracted by φ By Lemma 4, St(w 1 , K i ) = ActSt(w 1 , Ki ), 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, For every d-simplex σ ∈ K i , there is some d-simplex τ ∈ Σ that has been added in an elementary inclusion φ j with j < i, such that We call τ an origin of σ.There might be more than one origin of a simplex, but two distinct simplices in K i cannot have a common origin.Moreover, for every vertex v of σ, the tree in W i whose root is labeled with v contains exactly one vertex v of τ 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 ) \ E(y 2 )| such that no simplex is charged more than twice.Note that then fix an origin τ of σ.Then τ has a vertex that is a label in the subtree rooted at y 1 , so τ ∈ E(y 1 ).At the same time, since w 2 is not a vertex of σ, τ has no vertex in the subtree rooted at y 2 , so τ / ∈ E(y 2 ).We charge τ for the existence of σ.Because different simplices have different origins, every element in E(y 1 ) \ E(y 2 ) is charged at most once for St(w 1 , K i ).If σ ∈ Lk(w 1 , K i ), σ := w 1 * σ ∈ St(w 1 , K i ), and we can choose an origin τ of σ .As before, τ ∈ E(y 1 ) \ E(y 2 ) and we charge τ for the existence of σ.Again, each element in E(y 1 ) \ 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 ≥ 1, is a path in a forest such that x i+1 is the parent of x i , for 1 ≤ i < L. We call L the length of the path and x L its endpoint.For ascending paths in 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 .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 , . . ., x i ) denote the sub-path ending at x i , for i = 1, . . ., 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 ) ≤ 2 • |E(x i )|.For i = 1, this follows because c(p 1 ) = 0.For i = 2, . . ., L, x i is an internal node, and its two children are x i−1 and some other node x i−1 .Using induction and Lemma 8, we have that 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 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 W. Set W (0) ← W. In iteration i, we obtain the forest W (i) from W (i−1) by deleting the only-child-paths of W (i−1) .We stop when W (i) is empty; this happens eventually because at least the leaves of 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 W (1) arises from W (0) by deleting the leaves of W (0) .Note that the intermediate forests W (1) , W (2) , . . .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 W: Proof.The space complexity follows at once from the invariant, because the size of D is at most O(∆|K i |) during the i-th iteration.
For the time bound, we set S := | Km | for convenience and show that the algorithm finishes in O(∆•S) steps, assuming dictionary operations to be of constant cost.We have one simplex insertion per elementary inclusion which requires O(∆) operations.Thus, all inclusions are bounded by O(n∆), which is subsumed by our bound as n ≤ S. For the contraction case, we need O(c i ) to identify the smaller star, O(c i + ∆) to get a sorted list of St(u, ¬v) (or vice versa), and O(∆ • c i ) to add new vertices to D. Moreover, we delete the star of u from D; the cost for that is O(∆ • d i ), where d i is the number of deleted simplices.Thus, the complexity of a contraction is O(∆•(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 Km 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(∆ • S) as required.
Note that the dictionary D has lists of identifiers of length up to ∆+1 as keys, so that comparing two keys has cost O(∆).Therefore, using balanced binary trees as dictionary, we get a complexity of C ω = O(∆ log ω).Using hash tables, we get an expected worst-case complexity of C ω = O(∆).
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. 4o test its performance, we compared it to the software Simpers5 (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•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 • 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.

Tightness and lower bounds
The conversion theorem (Theorem 2) yields an upper bound of O(∆ • 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 ∆ 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 Ω(n log n 0 ).
Proof.Let p = 2 k for some k ∈ N. Consider a graph with p edges (a i , b i ), with a 1 , . . ., a p , b 1 , . . ., 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 0 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, . . ., k.This means that the sum of the costs of all level i contractions is exactly p 2 .Summing up over all i yields a cost of at least k 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 Ω(b), because adding a simplex causes either the birth or the death of exactly one bar in the 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 0-dimensional class, so add a vertex v 0 to the complex.From now, every time a 0-dimensional homology class is born add a new vertex to the complex.When a 0-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.

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: 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.

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 σ.In the addition case, this is followed by a list of simplex identifiers specifying the facets of σ.In other words, the element encodes the next column of the boundary matrix.For the inactive case, it means that σ has become inactive in the complex.In particular, no subsequent simplex in the stream will contain σ 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 as input and clears out all entries in row from the matrix.For that, let j be the column with as pivot.Traverse all non-zero columns of the matrix except column j.If a column i = j has a non-zero entry at row , 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 , we report ( , j) as a persistence pair and check whether is active.If not, we call remove_row( ).If the input stream specifies that simplex becomes inactive, we check whether j appears as pivot in the matrix and call remove_row( ) 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 of the column is inactive, no subsequent column in the stream will have an entry in row .Moreover, the reduction process cannot introduce an entry in row 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 , . . ., 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.

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 L with Ki ⊆ L ⊆ Ki+1 for some i = 0, . . ., m.Moreover, the active simplices define a subcomplex L ⊆ L and there is a moment during the algorithm where L = Ki and L = K i , for every i = 0, . . ., m.We call this the i-th checkpoint.We will make frequent use of the following simple observation.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 L = K i , there are not more than |K i | ≤ ω active simplices.Since each column has a different active pivot, their number is also bounded by ω.If we are between checkpoint i and i + 1, there have been not more than ω 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ω.
Proof.Consider a row index and a time in the algorithm where M represents L. By the same argument as in the previous lemma, we can argue that there are at most 2ω active row indices at any time.Therefore, we restrict our attention to the case that is inactive and distinguish three cases.If represents a negative simplex, we observe that its row should have been was removed due to the compression optimization, and after became inactive, no simplex can have it as a facet either.It follows that row is empty in this case.If is positive and was paired with another index j during the algorithm, then remove_row was called on , either at the moment the pair was formed, or when became inactive.Since the procedure removes the row, we can conclude that row is empty also in this case.The final case is that is positive, but has not been paired so far.It is well-known that in this case, is the generator of an homology class of L. Let For the 4.6 • 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.

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.

Figure 1 :
Figure 1: Construction example of L * , were u and v in K are contracted to w in L Dey et al.show that L ⊆ L * and that the map induced by inclusion is an isomorphism between H(L) and H(L * ).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 T , we can also obtain an non-zigzag filtration using coning, if we continue the operation on L * instead of going back to L.More precisely, we set K0 := K 0 and if φ i is an inclusion of simplex σ, we set Ki+1 := Ki ∪ {σ}.If φ i is a contraction (u, v) u, we set Ki+1 = Ki ∪ u * St(v, Ki ) .Indeed, it can be proved that ( Ki ) i=0,...,m has the same barcode as

Lemma 9 .
An ascending path with endpoint x has cost at most 2 • |E(x)|.An independent set of ascending paths in W has cost at most 2 • (∆ + 1) • n.

Figure 2 :
Figure 2: Iterations of the pruning procedure; the only-child-paths are marked in color y 1 , y 2 and r(y 1 ) = r(y 2 ) max{r(y 1 ), r(y 2 )}, if x has children y 1 , y 2 and r(y 1 ) = r(y 2 ) Proposition 13.The algorithm requires O(∆ • ω) space and O(∆ • | Km | • C ω ) time, where ω = max i=0,...,m |K i | and C ω is the cost of an operation in a dictionary with at most ω elements.

Figure 3 :
Figure 3: [In proof of Proposition 14] Sequence of contractions in the described construction for p = 8 (right) and the corresponding contracting forest (left), whose nodes contains the cost of the contractions For constant dimension, O(b) is also an upper bound: Lemma 15.For a barcode with b bars and maximal dimension ∆, there exists a filtration of size ≤ 2 ∆+2 b realizing this barcode.

Figure
Figure Examples of the influence of contractions on the homology classes: (a) destruction of four 1-homology generators, (b) creation of two 2-homology generators, and (c) no influence at all

Figure 5 :
Figure 5: Evolution of processing time (left Y-axis in sec) and process memory peak (right Y-axis in kB) depending on the chunk size (logarithmic X-axis)

Table 1 :
Experimental results.The symbol ∞ means that the calculation time exceeded 12 hours.