On ε Approximations of Persistence Diagrams
MMATHEMATICS OF COMPUTATIONVolume 00, Number 0, Pages 000–000S 0025-5718(XX)0000-0 ON ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS
JONATHAN JAQUETTE AND MIROSLAV KRAM ´AR
Abstract.
Biological and physical systems often exhibit distinct structures atdifferent spatial/temporal scales. Persistent homology is an algebraic tool thatprovides a mathematical framework for analyzing the multi-scale structuresfrequently observed in nature. In this paper a theoretical framework for thealgorithmic computation of an arbitrarily good approximation of the persistenthomology is developed. We study the filtrations generated by sub-level setsof a function f : X → R , where X is a CW-complex. In the special case X = [0 , N , N ∈ N we discuss implementation of the proposed algorithms.We also investigate a priori and a posteriori bounds of the approximationerror introduced by our method. Introduction
The formation of complex patterns can be encountered in almost all sciences.Many patterns observed in nature are extremely complicated and it is usually im-possible to completely understand them using analytical methods. Often a coarserbut computationally tractable description is needed. In recent years computationaltopology [8, 14] has become widely recognized as an important tool for quantify-ing complex structures. Computational homology has been used to quantitativelystudy spatio-temporal chaos [9] as well as to analyze complicated flow patternsgenerated by Rayleigh-B´enard convection [19, 20], time evolution of the isotropi-cally compressed granular material [15] and complex microstructures generated inbinary metal alloys through a process called spinodal decomposition [10, 21]. Un-fortunately homology groups can be extremely sensitive to small perturbations [17].Persistent homology [7] overcomes this problem [4] and provides a more powerfultool for studying evolution of complex patterns in the granular media [17, 16] andfluid dynamics [18].In the applications mentioned above the pattern is given by a sub- or super-level set of a real valued function f . Therefore analyzing the patterns usuallyinvolves numerical study based on a suitable discretization. It is important to knowwhen the homology groups of a sub-level set of f can be inferred from a discreteapproximation of this set. Probabilistic results, for a function f : [ a, b ] ⊂ R → R , arepresented in [24]. An alternative approach is to construct a cubical approximationof the sub-level set with the same homology groups [5]. The algorithm presented in[5] can handle functions defined on a unit square and uses a uniform grid on [0 , .If the resolution of a grid is fine enough, then the sub-level set can be approximated The first author’s research was funded in part by AFOSR Grant FA9550-09-1-0148 and NSFGrant DMS-0915019.The second author’s research was funded in part by NSF Grants DMS-1125174 and DMS-0835621. c (cid:13) XXXX American Mathematical Society a r X i v : . [ m a t h . A T ] J a n JONATHAN JAQUETTE AND MIROSLAV KRAM´AR by a collection of the grid elements that intersect the sub-level set. This often leadsto an unnecessary large complex. The size of the complex can be reduced by using arandomized version of the algorithm [3]. Using a regular CW complex to constructthe nodal approximation [6] simplifies the homology computations.In this paper we develop a method for approximating the sub-level sets of acontinuous function f : X → R , where X is a regular CW complex. Note thatthe sub-level set X t = f − ( −∞ , t ] = { x ∈ X : f ( x ) ≤ t } does not have to be a CWcomplex. The following definition provides conditions on a CW structure E of X under which the singular homology groups of the CW complex X t = (cid:91) { e ∈ E : f (cl( e )) ≤ t } , are isomorphic to the homology groups H ∗ ( X t ). Throughout the paper, we use thenotation “ f (cl( e )) ≤ t ” to mean “ f ( x ) ≤ t for all x ∈ cl ( e )”. Definition 1.1.
Let X be a finite regular complex with a CW structure E . Supposethat f : X → R is a continuous function and t ∈ R . We say that the CW structure E is compatible with X t if every cell e ∈ E satisfies one of the following conditions:(1) f (cl( e )) > t ,(2) f (cl( e )) ≤ t ,(3) there exists a deformation retraction h ( x , s ) : X t ∩ cl( e ) × [0 , → X t ∩ cl( e )of the set X t ∩ cl( e ) onto X t ∩ bd( e ).According to the following theorem, proved in Section 3, the singular homologygroups H ∗ ( X t ) can be obtained by computing the cellular homology groups of thecomplex X t . Theorem 1.2.
Let X be a finite regular CW complex. If the CW decomposition E of X is compatible with X t , then the map i ∗ : H ∗ ( X t ) → H ∗ ( X t ) , induced by theinclusion i : X t → X t , is an isomorphism. Due to computational complexity of computing homology, it is extremely impor-tant that the number of cells in the complex X t is as small as possible. To achievethis goal we use a multi-scale CW structure E . The diameter (size) of the cells in E varies with local complexity of the pattern. Larger cells are used at places withsmaller complexity while small cells might be necessary at the places where spacialscales of the pattern are small. To produce a multi-scale CW structure we startfrom a coarse CW structure of X and keep refining the cells that do not satisfyany of the conditions in Definition 1.1. For X = [0 , the number of cells in X t issmaller than the number of cells in the cubical complex used in [5, 6].To facilitate the refinement process we introduce self similar grids. The size ofthe grid elements might vary but all the grid elements of G are homeomorphic to X .The resolution of the self similar grid G can be dynamically adjusted at differentparts of X . To increase the resolution, each grid element can be refined into a unionof smaller complexes homeomorphic to X . We represent a self similar grid by atree structure. The root of this tree corresponds to the whole complex X . The i -thlevel nodes represent the complexes obtained by i refinements of X . Finally, theleaves correspond to the grid elements of G .The tree representation of G provides a basis for a memory efficient data struc-ture for storing the complex X t and its boundary operator. We stress that theboundary operator is not stored as a large matrix. Instead it is computed using the N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 3 boundary operator ∂ : H ∗ ( X ) → H ∗− ( X ) corresponding to the coarse (unrefined)CW structure, the tree structure representing the self similar grid and geometricinformation about intersections of the grid elements.While constructing a compatible CW structure for an arbitrary regular CWcomplex X might be challenging, it is straight forward when X = [0 , N , where N ∈ N . The coarse CW structure of [0 , N , employed in this paper, is given bystandard decomposition of a cube into vertices, edges and higher dimensional faces(cells). A self similar grid G on [0 , N consists of dyadic cubes of different sizes.If every cube (grid element) in G satisfies the following definition, then every cellin the CW structure E generated by the grid G satisfies one of the conditions inDefinition 1.1. Moreover, if a cell e ∈ E s i satisfies Condition (3) in Definition 1.1,then the function f ◦ h ( x , s ) can be taken to be non-increasing in s . Definition 1.3.
Let f ∈ C ([0 , N , R ) and C be an N dimensional dyadic cube.We say that C is ( f, t )- verified if 0 (cid:54)∈ ∂f∂x li ( C ) for some coordinate vectors x l , . . . , x l n ,where 0 ≤ n ≤ N , and for every N − n dimensional cell e of the cube C which isorthogonal to the vectors x l , . . . , x l n , then either f (cl( e )) > t or f (cl( e )) ≤ t .A cell e is considered to be orthogonal to a collection of vectors v , . . . , v n ifthe projection of e onto span { v , . . . , v n } is zero dimensional. The condition inDefinition 1.3 can be checked using interval arithmetics. Therefore the homologygroups H ∗ ( X t ) can be evaluated on a computer. However, if t is close to a criticalvalue, then using interval arithmetics to guarantee that the cubes surrounding thecritical point satisfy Definition 1.3 might require an extremely fine resolution ormight not be possible at all. Another problem with this approach is that thetopology of the set X t depends on the choice of the threshold t , and small changesof t can lead to large changes of the topology as demonstrated in [17]. So theparticular choice of the threshold t might be hard to justify. Both of these problemscan be overcome by using the persistent homology [7, 8].Persistent homology relates the homology groups of X t for different values of t . It reduces the function f to a collection of points in the plane. This collectionof points is called a persistence diagram and denoted by PD ( f ). Each point inthe persistence diagram encodes a well defined geometric feature of f . Hence thepersistent homology provides a finer description of the pattern than the homol-ogy groups of any one sub-level set. The set of all persistence diagrams can beturned into a complete metric space by using a Bottleneck distance d B between thediagrams [4, 22].The main contribution of this paper is a theoretical framework for constructingan ε approximation of the persistence diagram PD ( f ). For ε > f : X → R we define an ε approximate filtration of f as follows. Definition 1.4.
Let X be a finite regular CW complex and f : X → R a continuousfunction. Suppose that a sequence of real numbers { s i } m +1 i = − has the followingproperties:(1) s − = −∞ and s < min x ∈ X f ( x ),(2) s m > max x ∈ X f ( x ) and s m +1 = ∞ ,(3) | s i − s i − | < ε for 1 ≤ i ≤ m .Then the collection of sub-level sets { X s i } m +1 i = − is called an ε approximate filtration of f . JONATHAN JAQUETTE AND MIROSLAV KRAM´AR
It follows from the stability results [4] that the persistence diagram PD (cid:48) of any ε approximate filtration of f is close to PD ( f ). More precisely(1.1) d B ( PD ( f ) , PD (cid:48) ) < ε. Construction of the persistence diagram PD (cid:48) requires determining homology groupsof the sub-level sets in some ε approximate filtration of f . Freedom in choosing an ε approximate filtration of f allows us to stay away from the critical values. In orderto compute the persistence diagram of an ε approximate filtration { X s i } m +1 i = − onemight be tempted to replace this filtration by its cellular counterpart {X s i } m +1 i = − .However there is no guarantee that X s i ⊆ X s i +1 and so {X s i } m +1 i = − may not be afiltration. (See Figure 2 for a counterexample.) In order to construct a cellularfiltration we have to require that the CW structures {E s i } m +1 i = − are commensurable. Definition 1.5.
Let {E s i } m +1 i = − be a collection of CW structures on X . We saythat they are commensurable if for every e, e (cid:48) ∈ E (cid:48) := (cid:83) E s i such that e ∩ e (cid:48) (cid:54) = ∅ then either e ⊆ e (cid:48) or e (cid:48) ⊆ e .The CW structures {E s i } m +1 i = − constructed using a self similar grid are alwayscommensurable. The CW complexes (cid:8) ¯ X s i (cid:9) m +1 i = − define by ¯ X s i = (cid:84) j ≥ i X s j form an ε approximate cellular filtration of f . This filtration has the following property. Theorem 1.6.
Let (cid:8) ¯ X s i (cid:9) m +1 i = − be an ε approximate cellular filtration of f corre-sponding to {E s i } m +1 i = − . Suppose that for every cell e ∈ E s i which satisfies Condi-tion (3) in Definition 1.1, the function f ◦ h ( x , s ) is non-increasing in s . Then thepersistence diagrams of the filtrations { X s i } m +1 i = − and (cid:8) ¯ X s i (cid:9) m +1 i = − are equal. Corollary 1.7.
Let (cid:8) ¯ X s i (cid:9) m +1 i = − be an ε approximate cellular filtration of f . If PD (cid:48) is the persistence diagram of this filtration, then d B ( PD ( f ) , PD (cid:48) ) < ε .The persistence diagram PD (cid:48) can be computed even if we fail to construct CWstructures compatible with X s i for some of the thresholds s i . In this case thedistance d B ( PD ( f ) , PD (cid:48) ) can be evaluated a posteriori . The required resources forcomputing the persistence homology can be further reduced if zig-zag persistence[1] is used.The remainder of this paper is organized as follows. In Section 2 we recallthe definition of a CW-complex, homology and persistent homology. Theorem 1.2is proven in Section 3. In Section 4 we address the Inequality (1.1). A precisedefinition of a self similar grid is given in Section 5. These grids are employed inSection 6 to build a compatible CW structure. A memory efficient data structure forstoring the CW complex X t and its boundary operator is presented in Section 7. Inthe last section we use our method to compute the approximate persistence diagramof a two dimensional Fourier series with 5 modes and discuss the a posteriori errorestimates. 2. Background
In this section we recall definitions of a CW-complex [12], homology [14] andpersistent homology [7]. Notation introduced in this section is used throughout thepaper. We start by introducing a regular CW-complex.Let X be a Hausdorff space. A (finite) regular CW decomposition of X is a(finite) set E of subspaces of X with the following properties: N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 5 (1) E is a partition of X , that is X = (cid:83) e ∈E e and e (cid:54) = e (cid:48) ⇒ e ∩ e (cid:48) = ∅ , (2) Every e ∈ E is homeomorphic to some euclidean space R dim( e ) , (3) For every n -cell e ∈ E there exists a continuous embedding Φ e ( B n , S n − ) → ( X n − ∪ e, X n − ) such that Φ e ( B n ) = cl( e ) where cl( e ) is the closure of e , B n = { x ∈ R n : (cid:107) x (cid:107) ≤ } and S n = (cid:8) x ∈ B n +1 : (cid:107) x (cid:107) = 1 (cid:9) .The number dim( e ) is well determined and is called the dimension of e . The sets e ∈ E which are homeomorphic to R n are the n -cells . The boundary of a cell e is defined by bd( e ) := cl( e ) \ e . The union X n = (cid:83) dim ( e ) ≤ n e is the n -skeleton of the CW decomposition. The map Φ e is called a characteristic map for e and ϕ e = Φ e | S n − : S n − → X n − an attaching map. The space X is called a regularCW complex if there exists a regular CW decomposition of X .Cellular homology is useful for computing the homology groups of CW com-plexes. For a regular CW-complex X the n -chains are defined by C n ( X ) = H n ( X n , X n − ) . There is a one-to-one correspondence between the basis elements of C n ( X ) and the n -cells of X . For every n dimensional cell e the map Φ e generates a monomorphism(Φ e ) ∗ : H n ( B n , S n − ) → C n ( X ). For a fixed generator a ∈ H n ( B n , S n − ) the setΛ n = { (Φ e ) ∗ ( a ) : e ∈ E and dim( e ) = n } is a basis of C n ( X ). For every λ ∈ Λ = (cid:83) n ∈ Z Λ n we denote its correspondinggeometric cell by | λ | .The boundary operator ∂ n : C n ( X ) → C n − ( X ) is defined by the followingcomposition H n ( X n , X n − ) d n → H n − ( X n − ) j n − → H n − ( X n − , X n − )where d n is the connecting homomorphism from the long exact sequence of the pair( X n , X n − ) and j n − is the quotient map. In the above chosen basis of C n ( X ), theboundary operator ∂ n : C n ( X ) → C n − ( X ) is completely determined by the valueson the basis elements Λ n and can be uniquely expressed as ∂ n ( λ ) = (cid:88) µ ∈ Λ n − [ λ : µ ] µ. The coefficient [ λ : µ ] is called an incidence number of the cells. Regularity of thecomplex X implies that the incidence number is non-zero only for the pairs of cells λ ∈ Λ n , τ ∈ Λ n − such that bd( | λ | ) ∩ | τ | (cid:54) = ∅ . Moreover [ λ : τ ] = ± X is defined by C ( X ) = { C n ( X ) , ∂ n } n ∈ Z .The singular homology H ∗ ( X ) of X is isomorphic to the cellular homology of C ( X ),which is given by H ∗ ( C ( X )) = ker( ∂ ∗ ( C ∗ ( X )))im( ∂ ∗ +1 ( C ∗ +1 ( X ))) . We close this section with a short overview of persistent homology, following thepresentation given in [7]. Let f : X → R be a function defined on a CW complex X . We denote its sub-level sets by X t = f − ( −∞ , t ]. The function f is called tame if the homology groups of every sub-level set have finite ranks and there are onlyfinitely many values t across which the homology groups are not isomorphic. Let t < t . . . < t m be the values for which the homology of X t i − (cid:15) and X t i + (cid:15) differ forall 0 < (cid:15) <<
1. We define an interleaved sequence { s i } m +1 i = − such that s − = −∞ , s m +1 = ∞ and s i − < t i < s i , for 1 ≤ i ≤ m . For each − ≤ i ≤ j ≤ m + 1 JONATHAN JAQUETTE AND MIROSLAV KRAM´AR there is a natural inclusion of the set X s i into X s j and we denote the inducedhomomorphism between the corresponding homology groups by f i,jn : H n ( X s i ) → H n ( X s j ) . Persistent homology makes use of the above mentioned maps to compare the ho-mology groups of X t for different values of t . We say that the homology class α ∈ H n ( X s i ) is born at t i if it does not come from a class in H n ( X s i − ), that is α (cid:54)∈ im( f i − ,in ). The class α dies at t j if α is in im( f i,j − n ) but not in im( f i,jn ). Wedenote the birth-death pair for α by ( b α , d α ). The n -th persistence diagram PD n ( f )is defined to be the multi set of the points containing:(1) one point for each pair ( b α , d α ).(2) infinitely many copies of the points ( s, s ) at the diagonal.If X has dimension n , then the persistence diagram for the function f : X → R is the collection PD ( f ) = { PD i ( f ) } ni =0 . Actually, persistence diagrams can bedefined for continuous function f : X → R , where X is the realization of a simplicialcomplex as a topological space[2]. In this case the persistence diagrams can containaccumulation points at the diagonal. In other words there might be a sequence ofoff diagonal points converging to the diagonal. Using the same reasoning as in [2]the same is true if X is a finite regular CW-complex.Differences between the persistence diagrams can be assessed using a bottleneckdistance defined as follows. Definition 2.1.
Let PD = { PD i } ni =0 and P D (cid:48) = (cid:8) PD (cid:48) i (cid:9) ni =0 be persistence dia-grams. The bottleneck distance between PD and PD (cid:48) is defined to be d B ( PD , PD (cid:48) ) = max i inf γ : PD i → PD (cid:48) i sup p ∈ PD i (cid:107) p − γ ( p ) (cid:107) ∞ , where (cid:107) ( a , b ) − ( a , b ) (cid:107) ∞ := max {| a − a | , | b − b |} and γ ranges over all bi-jections. 3. Homology of sub-level sets
In this section we present a method for computing the homology of sub-level sets.Let X be a CW-complex and f : X → R a continuous function; we are interestedin the homology of the sub-level set X t := f − ( −∞ , t ]. Suppose that E is a CWstructure on X , then we define the cellular approximation to be(3.1) X t := (cid:91) { e ∈ E : f (cl( e )) ≤ t } . Note that X t is a closed subset of X t , and since X t is a finite union of cells, thenit is a sub-complex of X . For the function f shown in Figure 1 the set X ishomologous to X . According to Theorem 1.2, the homology groups of X t and X t are isomorphic if the CW-structure E is compatible with X t . We close this sectionby proving Theorem 1.2. Proof. If H ∗ ( X t , X t ) ∼ = 0 (using singular homology), then the fact that i ∗ is anisomorphism follows from the long exact sequence of the pair. We prove that H ∗ ( X t , X t ) ∼ = 0 by induction on the dimension of X .For a zero dimensional complex X , the set X t is equal to X t and so trivially H ∗ ( X t , X t ) ∼ = 0. Now we show that H ∗ ( X t , X t ) ∼ = 0 for an n dimensional complex X under the assumption that the statement is true for any complex with dimensionless than n . Let us first suppose that X contains a single n dimensional cell e . Then N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 7 (a) (b) (c)
Figure 1. (a) A simple function f : [0 , → R . The sub-levelset X = f − (( −∞ , , at whichthe value of f is below the semitransparent plane. (b) The set X is shown in purple. The black vertices, edges and squares boundedby the edges define a CW structure E on [0 , . (c) The set X corresponding to the decomposition E . X = X n − ∪ e . Define X n − t to be the n − X t . While X t is not a CWcomplex in general, we define an analog of its n − X n − t := X t ∩ X n − .By the inductive hypothesis, the homology groups H ∗ ( X n − t , X n − t ) are trivial.The CW structure E is compatible with X t , so the cell e satisfies one of theconditions in Definition 1.1. If f (cl( e )) > t , then X t = X n − t and X t = X n − t .Hence we obtain the set equality H ∗ ( X t , X t ) = H ∗ ( X n − t , X n − t ) and so by theinductive assumption the relative homology groups stay trivial.Now suppose that f (cl( e )) ≤ t . In this case both cl( e ) ⊆ X t and cl( e ) ⊆ X t . Ingeneral, we have the equality ( X t \ e, X t \ e ) = ( X n − t , X n − t ), so it follows from theinductive assumption that H ∗ ( X t \ e, X t \ e ) ∼ = 0. We would like to apply the excisiontheorem to obtain H ∗ ( X t , X t ) ∼ = H ∗ ( X t \ e, X t \ e ) ∼ = 0. However, the assumptionsof the singular excision theorem are not actually satisfied because cl( e ) is not inthe interior of X t . This problem can be overcome by expanding X t and X t inside e using a small ε neighborhood of bd( e ). Note that for ε small enough the homologygroups of X t and X t do not change.If e satisfies Condition (3) in Definition 1.1, then there exists a deformationretraction h ( x , s ) of the set X t ∩ cl( e ) onto X t ∩ bd( e ). Note that h is a continuousmap on the closed set X t ∩ cl( e ) × [0 , X t \ e × [0 ,
1] is also continuous. Since these two maps agree on the intersectionof their domains, we can extend h to a deformation retraction of X t onto X n − t .Existence of h also implies that there exists some t (cid:48) > t such that t (cid:48) ∈ f (cl( e )),hence e ∩ X t = ∅ and X t = X n − t . By induction H ∗ ( X t , X t ) ∼ = 0 again.Finally suppose that X is an arbitrary n dimensional complex. Let e , . . . , e m be its n dimensional cells. We repeat the above process m times. At each step weadd a cell e k to X n ∪ e , ∪ . . . ∪ e k − . (cid:3) JONATHAN JAQUETTE AND MIROSLAV KRAM´AR Persistent homology
In order to exactly compute the persistence diagrams PD ( f ) of a continuousbounded function, one must identify all thresholds t i for which the homology of X t i − (cid:15) and X t i + (cid:15) differ for all 0 < (cid:15) <<
1. If f is a Morse function, then thethresholds t i are the critical values of f , which one may identify after locating everypoint where the Jacobian Df vanishes. Such a root finding problem is frequentlystudied in nonlinear optimization [11]. In applications though, it may not alwaysbe feasible to rigorously find every root of Df . In such a situation, it would bedesirable to still be able to compute the persistence diagrams of f , even at the costof losing some accuracy.In this section we present a framework for computing an approximation of thepersistence diagram PD ( f ). For a given ε > f . The persistence diagram PD (cid:48) of this approxi-mation is ε close to PD ( f ), that is d B ( PD ( f ) , PD (cid:48) ) < ε . We start by choosing anapproximate filtration (cf. Definition 1.4). The following lemma shows that thebottleneck distance between the persistent diagram for the ε approximate filtrationof f and the persistence diagram of f is at most ε . Lemma 4.1.
Let X be a CW complex and f : X → R a continuous boundedfunction. If PD (cid:48) is a persistence diagram of any ε approximate filtration of f then d B ( PD ( f ) , PD (cid:48) ) < ε .Proof. To prove the lemma we have to show that for every n -th persistence diagram PD n ( f ) there exists a bijection γ from PD n ( f ) to PD (cid:48) n such that (cid:107) p − γ ( p ) (cid:107) ∞ < ε for every p ∈ PD ( f ).If p = ( b, d ) and d − b < ε , then we define γ ( p ) = ( b, b ). Clearly (cid:107) p − γ ( p ) (cid:107) ∞ < ε inthis case. The collection of sets { X s i } m +1 i = − is an (cid:15) filtration. Therefore, if d − b ≥ ε ,then there exist indices i and j such that s i − < b ≤ s i ≤ s j − < d ≤ s j . Sothe topological feature corresponding to p = ( b, d ) is not present in the sub-levelset X s i − but it can be found in X s i . Hence for the ε approximate filtration thereexists a unique topological feature q , corresponding to p , born at s i . By the samereasoning the feature q dies at s j and we define γ ( p ) = ( s i , s j ). Again the distance (cid:107) p − γ ( p ) (cid:107) ∞ < ε , and the number of points in PD (cid:48) is smaller or equal to the numberof points in PD ( f ). So γ is a bijection. (cid:3) It is appealing to construct a cellular approximate filtration by replacing the sets X s i with their cellular surrogates X s i introduced in the previous section. Howeverthe collection {X s i } m +1 i = − might not be a filtration, as shown in Figure 2. Namely,the set X is not a subset of X because of different choices of the compatible CWstructures for the sub-level sets. To avoid this problem we construct a commonrefinement of the CW structures. Our refinement process relies on the fact that theCW structures E s i are commensurable (cf. Definition 1.5). Definition 4.2.
Let {E s i } m +1 i = − be a commensurable collection of CW structureson X . Then the set E := { e ∈ E (cid:48) : either e ∩ e (cid:48) = ∅ or e ⊆ e (cid:48) for every e (cid:48) ∈ E (cid:48) } is called the common refinement of the CW structures {E s i } m +1 i = − .The CW structures E and E shown in Figures 2 (b-c) are commensurable andtheir common refinement is E . In this paper we will use the common refinement of N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 9 (a) (b) (c) (d)
Figure 2. (a) A contour plot of a function f : [0 , → R . (b) Thesub-level set X and its cellular approximation X obtained by thecompatible CW structure with four squares (2 dimensional cells.)(c) The sub-level set X and its cellular approximation X obtainedby the compatible CW structure with a single 2 dimensional cell.(d) The sub-level set X with the refined approximation ¯ X , whichis equal to X ∩ X .all the structures {E s i } m +1 i = − to compute the persistent homology. If zig-zag persis-tence [1] is used, then it suffices to consider the common refinements of consecutiveCW structures E s i and E s i +1 . We stress that the common refinement does notneed to be compatible with any of the sub-level sets X s i . Hence we need to de-velop a more sophisticated approach for building a cellular representation of theapproximate filtration. Before introducing a cellular representation of the filtration { X s i } m +1 i = − let us prove that the common refinement E is a regular CW structureon X . Lemma 4.3.
Let E be a common refinement of the commensurable CW structures {E s i } m +1 i = − on X . Then E is a CW structure on X .Proof. First we show that every point x ∈ X is contained in a unique cell e ∈ E .Note that the set C ( x ) = { e ∈ E (cid:48) : x ∈ e } is not empty. If e, e (cid:48) ∈ C ( x ), then byDefinition 1.5 either e ⊆ e (cid:48) or e (cid:48) ⊆ e . This implies that the cell e x := (cid:84) { e ∈ C ( x ) } is an element of E (cid:48) . Moreover for any e ∈ C ( x ) the cell e x ⊆ e and e ∩ e x = ∅ if e (cid:54)∈ C ( x ). Therefore e x ∈ E . Suppose by contradiction that e x is not a unique cellin E containing x . Then there exist e, e (cid:48) ∈ E such that x ∈ e ∩ e (cid:48) . It follows fromDefinition 4.2 that e ⊆ e (cid:48) and e (cid:48) ⊆ e . So e = e (cid:48) .We showed that E is a disjoint cellular decomposition of X . Naturally, everycell e ∈ E ⊆ E (cid:48) is homeomorphic to some euclidean space R dim( e ) . Now we needto define the characteristic maps for the cells. Every cell in e ∈ E belongs to atleast one E s i . We identify the characteristic map for e ∈ E with the characteristicmap given by one (arbitrary but fixed) CW structure E s i . Since every cell in E (cid:48) iscovered by a collection of smaller cells in E , the boundary of a cell e ∈ E is containedin the union of cells in E of dimension less than dim( e ). Hence the set of cells E equipped with these maps is a CW structure on X . (cid:3) Since the common refinement E does not have to be compatible with any of thesub-level sets, we cannot use it to generate a cellular approximation of the filtration { X s i } m +1 i = − . Instead we approximate the sub-level sets X s i by(4.1) ¯ X s i = (cid:92) j ≥ i X s j . Definition 4.4.
Let { X s i } m +1 i = − be an ε approximate filtration of a continuousfunction f and {E s i } m +1 i = − a collection of commensurable CW structures, where E s i is compatible with X s i for every s i . Then the collection of sets (cid:8) ¯ X s i (cid:9) m +1 i = − is calledan ε approximate cellular filtration of f corresponding to {E s i } m +1 i = − . Remark . Every set ¯ X s i is a sub-complex of X , where the CW structure of X isgiven by the common refinement E . Moreover (cid:8) ¯ X s i (cid:9) m +1 i = − is a filtration, that is ∅ = ¯ X s − ⊆ ¯ X s ⊆ · · · ⊆ ¯ X s m +1 = X. Our goal is to show that the persistence diagrams of the filtrations { X s i } m +1 i = − and (cid:8) ¯ X s i (cid:9) m +1 i = − are equal. This is true if the diagram shown below in Figure 3,with the maps induced by inclusion, commutes and the maps j i ∗ are isomorphisms.It follows from Remark 4.5 and the functorality of homology that the diagramcommutes. To prove that j i ∗ are isomorphisms we will have to strengthen thedefinition of a compatible CW structure. . . . −−−−→ H ∗ ( X s i − ) −−−−→ H ∗ ( X s i ) −−−−→ H ∗ ( X s i +1 ) −−−−→ . . . (cid:120) j i − ∗ (cid:120) j i ∗ (cid:120) j i +1 ∗ . . . −−−−→ H ∗ ( ¯ X s i − ) −−−−→ H ∗ ( ¯ X s i ) −−−−→ H ∗ ( ¯ X s i +1 ) −−−−→ . . . Figure 3.
The persistence diagrams of the filtrations { X s i } m +1 i = − and (cid:8) ¯ X s i (cid:9) m +1 i = − are equal if the above diagram commutes and themaps j i ∗ , induced by inclusion, are isomorphisms.First, we illustrate the idea on a simple example shown in Figure 2. The cells e ⊂ X that are not in ¯ X are inside of the single 2 dimensional cell e (cid:48) := (0 , ∈ E and additionally ¯ X = X ∩ bd( e (cid:48) ). Furthermore e (cid:48) ∩ X = ∅ because the image f (cl( e (cid:48) )) contains both values above and below 1. The fact that the CW structure E is compatible with X implies that there is a deformation retraction h ( x , s ) from X ∩ cl( e (cid:48) ) onto X ∩ bd( e (cid:48) ). If f ◦ h ( x , s ) is non-increasing in s , then it can berestricted to a deformation retraction from X ∩ cl( e (cid:48) ) onto X ∩ bd( e (cid:48) ). By the sameargument as in the proof of Theorem 1.2 the inclusion map from ¯ X = X ∩ bd( e (cid:48) )to X induces an isomorphism on the homology level. We now prove Theorem 1.6,which treats the general case. Proof.
The persistence diagrams of the filtrations { X s i } m +1 i = − and (cid:8) ¯ X s i (cid:9) m +1 i = − areequal if the diagram in Figure 3 commutes and the maps j ∗ , induced by inclusion,are isomorphisms. We already showed that the diagram computes. Now we haveto prove that every j i ∗ : H ∗ ( ¯ X s i ) → H ∗ ( X s i ) is an isomorphism. According toTheorem 1.2, the inclusion of each X s i into X s i induces an isomorphism. Hence it N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 11 suffices to show that the inclusion from ¯ X s i into X s i also induces an isomorphismin homology.The set ¯ X s i is a sub-complex of X s i . So¯ X s i = X s i \ { e ∪ e ∪ . . . ∪ e N } , for some cells e a ∈ E . We denote E = { e , e , . . . , e N } . First we divide the cellsin E into disjoint groups. Then we show that by removing all of the cells inthe first group we obtain a CW sub-complex of X s i and the map induced by theinclusion of this complex into X s i is an isomorphism. We prove that the maps¯ j i ∗ : H ∗ ( ¯ X s i ) → H ∗ ( X s i ) are isomorphisms by inductively removing the cells in allthe other groups.Before grouping the cells into disjoint classes let us show that for every cell e a ∈ E there exists e (cid:48) a ∈ E s i and e (cid:48)(cid:48) a ∈ E s j , where j > i , such that:(1) e a ⊆ e (cid:48) a (cid:40) e (cid:48)(cid:48) a ,(2) e (cid:48)(cid:48) a ∩ X s j = ∅ .By the construction of E , each cell e a is a subset of a unique cell e (cid:48) a ∈ E s i and e (cid:48) a ⊆ X s i . The fact that e a ∩ ¯ X s i = ∅ implies e a ∩ X s j = ∅ for some j > i. It followsfrom the construction of the cellular approximation X s j that e (cid:48) a (cid:54)∈ E s j . Otherwisewe would have e a ⊆ e (cid:48) a ⊆ X s j . Therefore there exists a cell e (cid:48)(cid:48) a ∈ E s j which satisfiesthe properties (1) and (2).Let { e (cid:48)(cid:48) , . . . , e (cid:48)(cid:48) N } be the cells with the properties (1) and (2) corresponding to thecells in E and let E (cid:48)(cid:48) = { e (cid:48)(cid:48) , . . . , e (cid:48)(cid:48) M } ⊆ { e (cid:48)(cid:48) , . . . , e (cid:48)(cid:48) N } be a set of the maximal cells.That is, a cell e (cid:48)(cid:48) i is removed from the set { e (cid:48)(cid:48) , . . . , e (cid:48)(cid:48) N } if there is another cell e (cid:48)(cid:48) j suchthat e (cid:48)(cid:48) i ⊆ e (cid:48)(cid:48) j . We order the cells in E (cid:48)(cid:48) by dimension so that dim( e (cid:48)(cid:48) i ) ≥ dim( e (cid:48)(cid:48) j ) for i ≤ j . We say that the cells e i , e j ∈ E are in the same group if there exists some l such that both e i and e j are subsets of e (cid:48)(cid:48) l .We will remove the cells e ∈ E from X s i in M steps. During the l -th step all ofthe cells e ∈ E such that e ⊂ e (cid:48)(cid:48) l are removed. Note that ¯ X s i = X s i \ ∪ Ma =1 e (cid:48)(cid:48) a . Dueto the ordering of the cells in E (cid:48)(cid:48) we get a CW-complex at each step. To show thislet e ∈ E be contained in some e (cid:48)(cid:48) l ∈ E (cid:48)(cid:48) . We show that e is not in the boundary ofany cell in X s i \ ∪ (cid:96)a =1 e (cid:48)(cid:48) a . The cell e (cid:48)(cid:48) l ∈ E s j for some j > i and e (cid:48)(cid:48) l (cid:54)∈ X s j . It followsfrom the construction of X s j that every cell c ∈ E s j such that e (cid:48)(cid:48) l ∩ cl( c ) (cid:54) = ∅ is notin X s j . So every cell in X s i whose closure intersects the cell e is contained in somecell e (cid:48)(cid:48) k ∈ E (cid:48)(cid:48) such that dim( e (cid:48)(cid:48) k ) > dim( e (cid:48)(cid:48) l ) and was removed before e . Thereforeat each step we get a closed subset of the cells from X s i which is a CW-complex.Moreover, at each step we can express the complex as a finite union of cells from E s i , so it also has a CW structure induced by E s i .First we suppose that l = 1 and prove that the map (cid:96) ∗ : H ∗ ( X s i \ e (cid:48)(cid:48) ) → H ∗ ( X s i ),induced by inclusion, is an isomorphism. To do so we apply the five lemma to thefollowing diagram. −−−−→ H ∗ ( X s i ∩ bd( e (cid:48)(cid:48) )) −−−−→ H ∗ ( X s i \ e (cid:48)(cid:48) ) ⊕ H ∗ ( X s i ∩ cl( e (cid:48)(cid:48) )) −−−−→ H ∗ ( X s i ) −−−−→ (cid:120) (cid:120) (cid:120) (cid:96) ∗ −−−−→ H ∗ ( X s i ∩ bd( e (cid:48)(cid:48) )) −−−−→ H ∗ ( X s i \ e (cid:48)(cid:48) ) ⊕ H ∗ ( X s i ∩ bd( e (cid:48)(cid:48) )) −−−−→ H ∗ ( X s i \ e (cid:48)(cid:48) ) −−−−→ Horizontal maps are given by the cellular Mayer-Vietoris sequences for the CWcomplexes given by the middle terms. The vertical maps are induced by inclusion.
The first two vertical maps are induced by identity maps and are thereby isomor-phisms. To satisfy the assumptions of the five lemma we have to show that themap k ∗ : H ∗ ( X s i ∩ bd( e (cid:48)(cid:48) )) → H ∗ ( X s i ∩ cl( e (cid:48)(cid:48) )) is also an isomorphism.We recall that the cell e (cid:48)(cid:48) ∈ E s j and e (cid:48)(cid:48) ∩ X s j = ∅ for some j > i . Thereforethe image f ( e ) contains points both above and below s j and there has to exista deformation retraction h ( x , s ) : X s j ∩ cl( e (cid:48)(cid:48) ) → X s j ∩ bd( e (cid:48)(cid:48) ). We have assumedthat f ◦ h ( x , s ) is non-increasing in s . So h restricts to a deformation retractionfrom X s i ∩ cl( e (cid:48)(cid:48) ) onto X s i ∩ bd( e (cid:48)(cid:48) ). Hence the top horizontal map in the diagrambelow is an isomorphism. H ∗ ( X s i ∩ bd( e (cid:48)(cid:48) )) ∼ = −−−−→ H ∗ ( X s i ∩ cl( e (cid:48)(cid:48) )) (cid:120) ∼ = (cid:120) ∼ = H ∗ ( X s i ∩ bd( e (cid:48)(cid:48) )) k ∗ −−−−→ H ∗ ( X s i ∩ cl( e (cid:48)(cid:48) ))There exist regular CW decompositions of bd( e (cid:48)(cid:48) ) and cl( e (cid:48)(cid:48) ) comprised of cellsfrom E s i . It follows from Theorem 1.2 that the vertical arrows in the diagramabove are isomorphisms. By the commutativity of the diagram, it follows thatthe map k ∗ is an isomorphism. So by the five lemma, it follows that the map (cid:96) ∗ : H ∗ ( X s i \ e (cid:48)(cid:48) ) → H ∗ ( X s i \ e (cid:48)(cid:48) ) is an isomorphism.Using mathematical induction we suppose that the inclusion map from X s i \∪ la =1 e (cid:48)(cid:48) a to X s i induces an isomorphism on the homology level. To show that this isalso true for l + 1 we repeat the above process replacing the set X s i by X s i \ ∪ la =1 e (cid:48)(cid:48) a and e (cid:48)(cid:48) by e (cid:48)(cid:48) l +1 . (cid:3) Cellular complexes on grid
In order to build a cellular approximation X t of X t we need to construct com-mensurable CW structures on the space X . We do this in an adaptive manner,starting with a coarse CW structure E on X . To obtain a cellular approximationwhich is compatible with a given sub-level set, the cells that do not satisfy any ofthe conditions in Definition 1.1 are repeatedly refined. If the critical values of f arebounded away from t , then this process terminates in a finite number of steps.To formalize this process we introduce the notion of a self similar grid. Usingthis type of grid also facilitates a memory efficient method for storing the complexand its boundary operator, see Section 7. At the end of this section we show howto compute the boundary operator for the refined CW structure on X in terms ofthe incidence numbers for the basis Λ given by E .Before defining a self similar grid in full generality we consider the case of a unitcube [0 , N . We choose E to be the standard CW structure on the cube consistingof the vertices, edges, higher dimensional faces and a single N-dimensional cell. Aself similar grid for the unit cube is built out of dyadic cubes defined as follows. Definition 5.1.
Let n be an arbitrary natural number. For every sequence { m i } Ni =1 ,where 0 ≤ m i < n , we define an N dimensional dyadic cube as a closed set C := N (cid:89) i =1 (cid:20) m i n , m i + 12 n (cid:21) . N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 13
Figure 4. (a) A dyadic grid G on a unit square. (b) CW structureson the dyadic cubes. (c) CW structure E G induced by G .Every dyadic cube is homeomorphic to [0 , N . Compositions of this homeomor-phism with the characteristic maps for the cells e ∈ E define a CW structure onthe dyadic cubes. A dyadic grid is a particular self similar grid consisting of dyadiccubes. Definition 5.2.
Let G = { C i } Mi =0 be a finite collection of N dimensional dyadiccubes. Then G is a dyadic grid if the following conditions are satisfied:(1) [0 , N = (cid:83) ni =1 C i ,(2) int ( C i ) ∩ int ( C j ) = ∅ if i (cid:54) = j .Figure 4(a) shows a simple example of dyadic grid on a unit square. The CWstructures of the dyadic cubes are depicted in Figure 4(b). Note that orientationof a cell depends on the precise choice of the homeomorphism h i : [0 , N → C i . Inour example we reversed orientation of the cells in two cubes at the bottom leftcorner.We generalize the idea of dyadic grids to CW complexes which are homeomorphicto a closed ball. In the rest of this paper we assume that there is a regular CWstructure E on X with a single N dimensional cell e such that X = cl( e ), thusimplying that X is homeomorphic to an N dimensional closed ball. The buildingblocks of a self similar grid are homeomorphic to X and their CW structure isagain induced by these homeomorphisms. They also have to cover X in a similarmanner as the dyadic cubes cover the unit cube. A self similar grid on X is a specialcollection of the building blocks. Definition 5.3.
Let X be an N dimensional regular CW complex with a single N dimensional cell e such that X = cl( e ). The set B is called a collection of buildingblocks of X if it satisfies the following conditions:(1) An element C ∈ B is a CW complex homeomorphic to X , and C ⊆ X .(2) Let C , C ∈ B and e i ⊂ C i for i = 1 ,
2. If e ∩ e (cid:54) = ∅ , then either e ⊆ e or e ⊆ e .(3) X ∈ B . Definition 5.4.
Let B be a collection of building blocks of X . Then G ⊆ B is a self similar grid on X if the following conditions are satisfied:(1) X = (cid:83) C ∈G C , (2) For all C , C ∈ G the intersection int ( C ) ∩ int ( C ) = ∅ if C (cid:54) = C .Using the same arguments as in the proof of Lemma 4.3, one can show that byremoving the redundant cells, a self similar grid G on X generates a CW structure E G . For our simple example of a dyadic grid, this structure is shown in Figure 4(c).If we fix a collection of building blocks B , then a collection of CW structures {E G i } generated by the self similar grids {G i } is always commensurable. The followingtwo lemmas formalize the above statements. Lemma 5.5.
Let G = { C i } Mi =0 be a self similar grid on X and E (cid:48)G := { e : e is a cell of C i for some i = 0 , . . . , M } . Then E G := (cid:8) e ∈ E (cid:48)G : either e ∩ e (cid:48) = ∅ or e ⊆ e (cid:48) for every e (cid:48) ∈ E (cid:48)G (cid:9) is a CW structure on X generated by G . Lemma 5.6.
Let B be a collection of building blocks of X and {G i } m +1 i = − be acollection of self similar grids on X such that G i ⊆ B for i ∈ {− , . . . , m + 1 } . Thenthe CW structures {E G i } m +1 i = − , generated by the grids {G i } m +1 i = − , are commensurable. In the second part of this section we analyze the boundary operator associatedwith the CW structure E G . Our goal is to express this operator using the incidencenumbers given by the basis Λ corresponding to E . We start by building a basis forthe chain complex C ∗ ( X ) induced by the cellular structure E G . Every cell e ∈ E G belongs to some C i ∈ G . We associate this cell with one of these (arbitrary butfixed) building blocks C i . Let h i be the homeomorphism from X to C i . Thenthere is a cell e (cid:48) ∈ E such that e = h i ( e (cid:48) ). Now the characteristic map for e isgiven by h i ◦ Φ e (cid:48) where Φ e (cid:48) is the characteristic map for e (cid:48) . The map h i inducesa chain isomorphism between the chain complexes, ( h i ) ∗ : C ∗ ( X ) → C ∗ ( C i ). Fix λ ∈ Λ to be the basis element associated with the cell e (cid:48) ∈ E , and (cid:96) i : C i → X theinclusion map. We choose ( (cid:96) i ◦ h i ) ∗ ( λ ) to be a basis element in C ∗ ( X ) with thecellular structure E G , corresponding to the geometric cell e . We denote ( (cid:96) i ◦ h i ) ∗ ( λ )by ( i, λ ) and the corresponding basis by Λ G . As before | i, λ | is the topological cellassociated with ( i, λ ).To define a boundary operator in the basis Λ G using the incidence numbers forΛ, we need to understand the relationship between the cells belonging to differentbuilding blocks. Namely, we need a method for comparing the orientations of twointersecting cells with the same dimension, and we do so using their local homologygroups.Let us consider cells e i = | i, λ | ⊂ C i and e j = | j, µ | ⊂ C j such that n = dim( e i ) =dim( e j ) and e i ∩ e j (cid:54) = ∅ . Without loss of generality we can suppose that e i ⊆ e j andselect a point q ∈ e i . Since bd( e i ) and bd( e j ) are deformation retracts of cl( e i ) \ q and cl( e j ) \ q respectively, then it follows that the maps induced by inclusion r ∗ : H ∗ (cl( e i ) , bd( e i )) → H ∗ (cl( e i ) , cl( e i ) \ q ) s ∗ : H ∗ (cl( e j ) , bd( e j )) → H ∗ (cl( e j ) , cl( e j ) \ q )are isomorphisms. By excision, the map t ∗ : H ∗ (cl( e i ) , cl( e i ) \ q ) → H ∗ (cl( e j ) , cl( e j ) \ q )induced by inclusion is also an isomorphism. Thereby we obtain the following com-mutative diagram where all of the maps are isomorphisms: N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 15 H ∗ (cl( e i ) , bd( e i )) H ∗ (cl( e i ) , cl( e i ) \ q ) H ∗ ( B n , S n − ) H ∗ (cl( e j ) , bd( e j )) H ∗ (cl( e j ) , cl( e j ) \ q ) r ∗ t ∗ ( h j ◦ Φ | µ | ) ∗ ( h i ◦ Φ | λ | ) ∗ s ∗ Hence for a fixed generator a ∈ H ∗ ( B n , S n − ) then(5.1) r − ∗ ◦ t − ∗ ◦ s ∗ ◦ ( h j ◦ Φ | µ | ) ∗ ( a ) = { ( i, λ ) : ( j, µ ) } ( h i ◦ Φ | λ | ) ∗ ( a )where { ( i, λ ) : ( j, µ ) } is a unit. If we are working with Z coefficients then { ( i, λ ) : ( j, µ ) } = ± e i and e j have the same orientation if { ( i, λ ) : ( j, µ ) } = 1 andthe opposite otherwise.First we evaluate a boundary for some cells in the complex shown in Figure 4.The boundary of the two dimensional cell in the top left square C consists of fiveone dimensional cells. Three of them belong to C . However the bottom edge wasreplaced by two smaller edges. Orientation of the removed edge agrees with thesmaller edges. Hence we can simply replace this edge by those smaller edges. Forthe bottom right square the problem is a little more complicated. The missing leftedge is divided into two smaller edges but orientation of the bottom smaller edgeis opposite. Therefore we need to reverse orientation of this cell by multiplying itby −
1. The following lemma formalizes this procedure.
Lemma 5.7.
Let ( i, λ ) ∈ Λ n G . Then (5.2) ∂ n ( i, λ ) = (cid:88) µ ∈ Λ n − [ λ : µ ] (cid:88) ( j,τ ) ∈N ( i,µ ) { ( i, µ ) : ( j, τ ) } ( j, τ ) , where N ( i, µ ) = (cid:8) ( j, τ ) ∈ Λ n − G : | i, µ | ∩ | j, τ | (cid:54) = ∅ (cid:9) .Proof. First we consider the boundary of ( i, λ ) inside C i equipped with the basis { ( i, λ ) : λ ∈ Λ } . This boundary can be written as (cid:80) µ ∈ Λ n − [ λ : µ ]( i, µ ). Howeverthis is not necessarily the boundary of ( i, λ ) in the basis Λ G . This is becausesome elements ( i, µ ) might be represented by basis elements ( j, τ ) correspondingto a different building block, or they can be broken up to smaller pieces also fromdifferent building blocks. The set N ( i, µ ) ⊂ Λ G is a representation of ( i, µ ) inthe basis Λ G . Geometrically, this means that cl( | i, µ | ) = (cid:83) ( j,τ ) ∈N ( i,µ ) cl( | j, τ | ). Toobtain ∂ n ( i, λ ) we have to replace each ( i, µ ) by the elements in N ( i, µ ). In order toaccommodate possible differences in orientation, the sign of the incidence number[ λ : µ ] has to be changed if ( j, τ ) ∈ N ( i, µ ) and ( i, µ ) have opposite orientation.Therefore each term contained in the second sum is multiplied by { ( i, µ ) , ( j, τ ) } . (cid:3) Construction of compatible CW structures on a unit cube
By now the results presented in this paper have been mostly theoretical. Inthis section we restrict our attention to functions f ∈ C ([0 , N , R ) and presentan algorithm for constructing a CW structure compatible with the sub-level set X t = f − (( −∞ , t ]). It is complicated to check the conditions in Definition 1.1algorithmically, especially Condition (3). Instead we introduce analytical conditionsfor dyadic cubes that force their cells to satisfy the conditions in Definition 1.1. These conditions can be checked using interval arithmetic, and the algorithm canbe carried out by a computer. We divide this section in two parts. In the first partwe present the conditions posed on the dyadic cubes. The second part describesthe algorithm for constructing a compatible CW structures on a unit cube.6.1.
Verified Cubes.
The dyadic grid G is said to be ( f, t )- verified if every dyadiccube C ∈ G is ( f, t )-verified (cf. Definition 1.3). We will show that if G is ( f, t )-verified, then the associated CW complex E G is compatible with X t . Note thatevery cell in E G belongs to at least one dyadic cube C ∈ G . Hence it sufficesto show that every cell in a ( f, t )-verified cube C ∈ G satisfies one of the threeconditions in Definition 1.1. The following lemma allows us to prove this for justthe top dimensional cell. Lemma 6.1.
If a cube C is ( f, t ) -verified, then for every cell e ⊆ C , the cube cl( e ) is ( f, t ) -verified.Proof. Fix vectors x l , . . . , x l n such that 0 / ∈ ∂f∂x li ( C ) for all i ≤ n , with 0 ≤ n ≤ N .Select any cell e ⊆ C . There exists a k ≤ n such that if we reorder the sequence { l i } , then 0 / ∈ ∂f∂x li (cl( e )) for all i ≤ k and cl( e ) is orthogonal to x l i for any i > k .Select any dim( e ) − k dimensional cell e (cid:48) ⊆ cl( e ) such that e (cid:48) is orthogonal to thevectors x l , . . . , x l k . Since C is a cube, then there exists a N − n cell e (cid:48)(cid:48) ⊆ C suchthat e (cid:48) ⊆ e (cid:48)(cid:48) and e (cid:48)(cid:48) is orthogonal to the vectors x l , . . . , e l n . Since C is ( f, t )-verified, then either f (cl( e (cid:48)(cid:48) )) > t or f (cl( e (cid:48)(cid:48) )) ≤ t . Thereby either f (cl( e (cid:48) )) > t or f (cl( e (cid:48) )) ≤ t , thus showing that cl( e ) is ( f, t )-verified. (cid:3) There are three special cases of a cube C being ( f, t )-verified corresponding to thecases when n = 0 and n = N :(1) f ( C ) > t (2) f ( C ) ≤ t (3) 0 / ∈ ∂f∂x i ( C ) for 1 ≤ i ≤ N .If a ( f, t )-verified cube C satisfies one of Conditions (1) or (2) above, then all of itscells satisfy Condition (1) or (2) in Definition 1.1 respectively. On the other hand, ifa cube C satisfies neither of Conditions (1) or (2) above, then we must construct adeformation retraction required by Condition (3) in Definition 1.1. Using rescalingwe can further restrict our attention to a unit cube.Before constructing the deformation retraction in full generality, we demonstratethe basic ideas on two examples shown in Figure 5. The unit cube depicted inFigure 5(a) satisfies Condition (3) above. In this case both partial derivatives arenegative. For ¯ x ∈ int([0 , \ X t ) we define a vector field ξ ( x ) = x − ¯ x . The set X t can be retracted onto X t ∩ ∂ [0 , along the flow lines, indicated by the arrows inFigure 5(a), of the flow generated by ξ ( x ). To ensure that the function f ◦ h ( x , s )is non-increasing in s we have to choose ¯ x sufficiently close to the origin.For the example shown in Figure 5(b) the partial derivative ∂f∂x ([0 , ) < ∂f∂x is both positive and negative on [0 , . In this case the cells e = { }× (0 ,
1) and e = { }× (0 ,
1) are orthogonal to the coordinate vector x , and satisfy f (cl( e )) > t and f (cl( e )) ≤ t . Hence cube [0 , is ( f, t )-verified. In this case we define thevector field ξ ( x , x ) = ( x − ¯ x , ε ( x − )) for some point (¯ x , ) ∈ int( X t \ [0 , )and ε small enough. Again we retract the set X t along the flow lines indicated bythe arrows. The first coordinate of the vector field ξ is similar to the one used in N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 17
Figure 5. (a) A deformation retract of a sub-level set to theboundary when the cube satisfies Condition (3). (b) A deformationretract of a sub-level set to the boundary for a generic ( f, t )-verifiedcube.the previous example. Part of the set X t ∩ ∂ [0 , is contained inside the horizontaledges of the cube. The second coordinate of the vector field ξ ensures that someflow lines reach this part of the set X t ∩ ∂ [0 , . The following two lemmas formalizethe above discussion. Lemma 6.2. If / ∈ ∂f∂x i ([0 , N ) for ≤ i ≤ N and f passes through t , then thereexists a deformation retraction h from X t onto X t ∩ ∂ [0 , N . Moreover f ◦ h ( x , s ) is non-increasing in s .Proof. Without loss of generality we can suppose that ∂f∂x i ( x ) < x ∈ [0 , N and 1 ≤ i ≤ N . Thereby the value of f at the origin has to be greater than t . To prove the existence of a deformation retraction h : X t → X t ∩ ∂ [0 , N weuse the Wa˙zewski principal. We will show that X t is a Wa˙zewski set for a flow ϕ corresponding to an appropriate linear vector field ξ ( x ) on R N .The vector field is given by ξ ( x ) = x − ¯ x where ¯ x ∈ int([0 , N \ X t ) such that( x − ¯ x ) · ∇ f ( x ) < x ∈ X t . To show that there exists such ¯ x we analyzethe continuous function x · ∇ f ( x ). For every x ∈ [0 , N the function x · ∇ f ( x ) ≤ ε > x · ∇ f ( x ) ≤ − ε on X t . Moreover for δ > x · ∇ f ( x ) > − ε for all x ∈ B δ := (cid:8) y ∈ [0 , N : (cid:107) y (cid:107) < δ (cid:9) . Wechoose ¯ x to be an arbitrary but fixed point in B δ ∩ int([0 , N ). Clearly ¯ x (cid:54)∈ X t and( x − ¯ x ) · ∇ f ( x ) < x ∈ X t .Let ϕ be a flow given by the vector field ξ ( x ). We define the exit set and theimmediate exit set for X t as follows: W = { x ∈ X t : ∃ s > ϕ ( x , s ) / ∈ X t } ,W − = { x ∈ X t : ∀ s > ϕ ( x , [0 , s )) (cid:54)⊂ X t } . The unstable fixed point of the linear system x (cid:48) = ξ ( x ) is not in X t . Thereforeall the orbits have to leave X t and W = X t . By the chain rule dds f ◦ ϕ ( x , s ) =( x − ¯ x ) · ∇ f ( x ) which is negative for all x ∈ X t and f ◦ ϕ ( x , s ) is strictly decreasingwith respect to s for all x ∈ X t . Hence for every x ∈ X t ∩ (0 , N there exists s x > such that ϕ ( x , [0 , s x )) ⊆ X t and the immediate exit set satisfies W − ⊆ X t ∩ ∂ [0 , N .The point ¯ x is in int([0 , N ). Thus for every x ∈ ∂ [0 , N the vector field ξ ( x ) pointsout of the unit cube and ϕ ( x , s ) / ∈ [0 , N for s >
0. Thereby W − = X t ∩ ∂ [0 , N .The sets X t and X t ∩ ∂ [0 , N are closed so it follows that X t is a Wa˙zewskiset. By the Wa˙zewski principal there exists a deformation retraction h ( x , s ) from W = X t onto W − = X t ∩ ∂ [0 , N . The set X t is retracted along the flow lines of ϕ . So f ◦ h ( x , s ) is non-increasing in s . (cid:3) Lemma 6.3. If [0 , N is ( f, t ) -verified and f passes through t , then there existsa deformation retraction h from X t onto X t ∩ ∂ [0 , N . Moreover f ◦ h ( x , s ) isnon-increasing in s .Proof. Without loss of generality we can suppose that ∂f∂x i ( x ) < x ∈ [0 , N and 1 ≤ i ≤ k . We denote by π A projection on the first k coordinates and by π B projection on the last N − k coordinates. The cube B = π B ([0 , N ) is the closureof a cell in [0 , N which is orthogonal to the vectors x , . . . , x k , so either f ( B ) > t or f ( B ) ≤ t . Since t ∈ f ([0 , N ) and ∂f∂x i ([0 , N ) < ≤ i ≤ k then it followsthat f ( B ) > t .We apply the Wa˙zewski principal to a linear flow ϕ radiating from some point ¯ x .By the same argument as in the previous proof one can show that for some ε > x = (¯ x , . . . , ¯ x k , / , . . . , / ∈ int([0 , N \ X t ) such that π A (¯ x − x ) ·∇ f ( x ) < − ε for all x ∈ X t . We now define the following vector field ξ ( x ) = π A ( x − ¯ x ) + ε η π B ( x − ¯ x )where η = sup x ∈ X t | π B ( x − ¯ x ) · ∇ f ( x ) | . The value η has to be finite because it issupremum of a continuous function on the compact set X t . Define ϕ to be the flowgenerated by the vector field ξ .The function f is non-increasing along the flow lines of ϕ as corroborated by thefollowing calculation ddt f ◦ ϕ ( x , s ) = (cid:18) π A ( x − ¯ x ) + ε η π B ( x − ¯ x ) (cid:19) · ∇ f ( ϕ ( x , s )) < − (cid:15) + (cid:15) η η < . Applying the Wa˙zewski principal to the same sets W and W − as in the proof ofthe previous lemma guarantees that there exists a deformation retraction h from X t onto X t ∩ ∂ [0 , N and f ◦ h ( x , s ) is non-increasing in s . (cid:3) Finally we are in the position to show that an ( f, t )-verified grid G induces a CWstructure E G compatible with X t . Theorem 6.4.
Let f ∈ C ([0 , N , R ) . If the dyadic grid G on [0 , N is ( f, t ) -verified, then the CW structure E G is compatible with X t .Proof. If e is a cell in E G , then there exists some ( f, t )-verified cube C ∈ G whichcontains e . By Lemma 6.1 it follows that cl( e ) is ( f, t )-verified. If either f (cl( e )) >t or f (cl( e )) ≤ t , then cl( e ) will satisfy Condition (1) or (2) in Definition 1.1,respectively. Otherwise f passes through t , and so by Lemma 6.3 there exists a N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 19 C C C C C C C C C C C C C C Figure 6. (a) A dyadic cubical grid G with seven building blocks.(b) Tree representation of G . A building block C i ∈ G correspondsto the leaf with same label. (c) Basis Λ G chosen in Figure 4(c) isstored in the leaves of the tree corresponding to G .deformation retraction h of the set X t ∩ cl( e ) onto X t ∩ bd( e ). Moreover the function f ◦ h ( x , s ) is non-increasing in s . Thereby E G is compatible with X t . (cid:3) Remark . The fact that the function f ◦ h ( x , s ) is non-increasing in s impliesthat the sets X t obtained using the dyadic grid can be used to produce a cellularfiltration (see Theorem 1.6).6.2. Algorithm for Constructing a Verified Grid.
Recall that a dyadic grid G is ( f, t )-verified if every cube C ∈ G is ( f, t )-verified. To construct a verified grid westart with a grid G containing a single cube [0 , N . Then we iteratively refine thecubes C ∈ G that are not ( f, t )-verified. If C is not ( f, t )-verified, then we replaceit by 2 N dyadic cubes inside C whose edge length equals to half of the edge lengthof C . We call those smaller cubes the offspring of C . If the following algorithmterminates, then the returned grid G is ( f, t )-verified. G = ConstructVerifiedGrid ( f, t ) G = (cid:8) [0 , N (cid:9) while there exists a cube C ∈ G which is not ( f, t )-verified do Remove C from G Add the offspring of C to G end whilereturn G Data Structures and Algorithms
In order to compute the homology of the sub-level set X t ⊆ X we need todevelop data structures for storing the basis C ∗ ( X t ) and the boundary operator. InSection 7.1 we introduce a representation of the basis Λ G associated with the grid G . A special structure of self similar grid is exploited to devise a memory efficientdata structure. We do not store the boundary operator for the basis Λ G which isusually a large (sparse) matrix and requires a lot of memory to be stored. Only theincidence numbers for the coarse CW structure E G need to be stored. In the case ofa unit cube the incidence numbers can be computed and the memory requirementsare further reduced.In Section 7.2, the boundary operator ∂ n : Λ n G → Λ n − G is computed using For-mula (5.2). Evaluating the terms in (5.2) involves manipulating geometric objects.This might be complicated for a general self similar grid. However for a dyadic grid it can be performed using a special representation of the cells, as explained inSection 7.3.7.1. Representation of Λ G . A collection of building blocks B , given by Defini-tion 5.3, can be partially ordered by inclusion. For C i , C j ∈ B we say that C i precedes C j and write C i ≺ C j if C j ⊆ C i . The fact that X ≺ C for all C ∈ B implies that the graph representing the partial ordering is a tree T B with the rootcorresponding to the set X . Therefore a self similar grid G = { C i } Mi =0 can be repre-sented by a subtree T G of T B containing the vertices corresponding to the buildingblocks C i ∈ G and their predecessors.Figure 6 shows an example of a self similar grid with seven building blocks and itstree representation T G . There is a one to one correspondence between the buildingblocks C i ∈ G and the leaves of the tree T G . The nodes of the tree correspond to thebuilding blocks that were refined in the process of creating the grid. In our examplethe root corresponds to the unit square and the unlabeled node in the middle layerto the subdivided cube at the bottom left corner. Definition 7.1.
For every node (leaf) s in T G we denote the corresponding buildingblock by B ( s ). If s is not a leaf, then B ( s ) (cid:54)∈ G . For a leaf s we define Λ n ( s ) := (cid:8) ( i, λ ) ∈ Λ n G stored in the leaf s (cid:9) .The basis Λ G is stored in the leaves of the tree T G . Figure 6(c) shows how thebasis of the CW structure chosen in Figure 4(c) is stored. In general, the leafcorresponding to the building block C i ∈ G stores the basis elements ( i, λ ) ∈ Λ G .7.2. Boundary Operator.
Evaluation of Formula (5.2) requires construction ofthe set N ( i, µ ). In this section we present a tree traversing algorithm for construct-ing N ( i, µ ). The algorithm consists of two parts. In the first part, the algorithmmoves up the tree. It stops at the node s such that | i, µ | ⊆ int( B ( s )). We willprove that the elements of N ( i, µ ) are stored in the leaves that can be reached bydescending down form the node s . In the second part, the algorithm descends tothe leaves that can contain elements of N ( i, µ ) and retrieves them. Before statingthe algorithm we recall some basic terminology. Definition 7.2.
Let T be a tree. Suppose that s , s (cid:48) are two nodes of T connectedby an edge. We say that s is a predecessor (child) of s (cid:48) if the path from s to theroot is shorter (longer) than the path from s (cid:48) to the root. N = ConstructN ( i, µ )Let s be a leaf corresponding to C i ∈ G . while | i, µ | (cid:54)⊆ int( B ( s )) and s is not the root do Set s to the predecessor of send while Put s to the stack H while H (cid:54) = ∅ do Remove one element s from H . if s is a leaf thenfor every ( j, τ ) ∈ Λ dim( µ ) ( s ) such that | i, µ | ∩ | j, τ | (cid:54) = ∅ do Put ( j, τ ) to N ( i, µ ) end forelse N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 21 [0 ,
0] [0 , , , (0) (1)(2) (3) (0 , , ,
2) (0 ,
3) (1 , , , , ,
0) (2 , ,
2) (2 , , , , , Figure 7. (a) Binary (b) decimal indices for the dyadic cubesobtained by one subdivision of a unit square. (c) Binary indicesfor the dyadic cubes obtained by two subdivision of a unit square. for every child s (cid:48) of s such that | i, µ | ∩ B ( s (cid:48) ) (cid:54) = ∅ do Add s (cid:48) to H end forend ifend whilereturn N The following lemma guarantees that the set N retrieved by ConstructN ( i, µ ) isequal to the set N ( i, µ ) given by (5.2). Lemma 7.3.
For every basis element ( i, µ ) ∈ Λ G the algorithm ConstructN ( i, µ ) returns the set N ( i, µ ) .Proof. Let s be the node reached at the end of the ascending part (first whilecycle) of the algorithm. We show that every leaf containing an element of N ( i, µ )can be reached by descending from the node s . If s is the root, then this is trivial.If s is not the root, then | i, µ | ⊆ int( B ( s )). For every ( j, τ ) ∈ N ( i, µ ) the cell | j, τ | ⊆ | i, µ | ⊆ int( B ( s )). So int( C i ) ∩ int( B ( s )) (cid:54) = ∅ for every building block C j ∈ G , such that | j, τ | ⊆ C j . The interior of every building block is formedby a single N dimensional cell. Hence by Condition (2) of Definition 5.3 eitherint( C j ) ⊆ int( B ( s )) or int( B ( s )) ⊆ int( C j ). If int( B ( s )) ⊆ int( C j ), then thereexists a building block C k ∈ G , corresponding to some leaf below s , such thatint( C j ) ∩ int( C k ) (cid:54) = ∅ . This violates Condition (2) of Definition 5.4. Therefore C j ⊆ B ( s ) and the leaf, containing ( j, τ ) can be reached by descending from thenode s .In the descending part, the algorithm reaches every leaf s corresponding to abuilding block C j that intersects | i, µ | . All the elements of Λ dim( µ ) ( s ) which intersect | i, µ | are added to N . Hence N = N ( i, µ ). (cid:3) Dyadic Grid.
In this section we deal with implementation of the algorithmfor a dyadic grid. First we label the dyadic cubes. Our labeling schema is motivatedby the refinement process applied to the unit cube. During a single refinement thecube [0 , N is divided into 2 N cubes with edge lengths equal to 1 /
2. Each of thesecubes contains exactly one vertex of the original cube. Let C be a cube containingthe vertex ( s , s , . . . , s N − ). The coordinates of the vertex can be interpreted as a binary number. The number s = (cid:80) N − n =0 s n n is the index of the cube C . Note thatthe index contains complete information about position of the cube. Coordinatesof the smallest vertex (in the lexicographical order) of C are given by x i = s i , seeFigure 7.Further refinement of the cube C , considered above, produces another collectionof 2 N dyadic cubes. Each of these cubes contains exactly one vertex of C . Again,the position of each cube inside C is encoded by a binary number with N digits.Therefore every cube obtained by two consecutive refinements of [0 , N is indexedby a sequence ( s, s (cid:48) ), see Figure 7(c). In general a dyadic cube obtained by n sub-divisions is uniquely indexed by s = ( s , s , . . . , s n ) where s i ∈ (cid:8) , , . . . , N − (cid:9) .The position of the cube corresponding to s is given by(7.1) C s = N − (cid:89) i =0 n (cid:88) j =1 − j s ij , n (cid:88) j =1 − j (1 + s ij ) , where s ij is the i -th binary digit of s j .We choose a CW structure for the unit cube to be the standard cellular structureon [0 , N . Every cell e is a product of N generalized intervals and can be uniquelyindexed by λ = ( λ , λ , . . . , λ N − ) ∈ Z N . In particular the cell corresponding tothe index λ is | λ | = I × I × . . . × I N − where I i = λ i = 0 , λ i = 1 , (0 ,
1) if λ i = 2.If the cell e (cid:48) is a translation of the cell e , then Φ e (cid:48) = T ◦ Φ e , where T is thetranslation map. The characteristic map of the cell e whose closure contains theorigin can be chosen arbitrarily. The orientation of the cells e (cid:48) that do not containthe origin in their closure is given by Φ e (cid:48) .To define a CW structure on C s we use the homeomorphism h s : [0 , N → C s defined coordinate-wise by(7.2) h i ( x ) = 2 − n x i + n (cid:88) j =1 − j s ij , where h s = ( h ( x ) , h ( x ) , . . . , h N − ( x )). Our choice of orientation and homeomor-phisms h s guarantees that two intersecting cells have the same orientation. SoFormula (5.2) reduces to ∂ n ( i, λ ) = (cid:88) µ ∈ Λ n − [ λ : µ ] (cid:88) ( j,τ ) ∈N ( i,µ ) ( j, τ ) . Implementation of the algorithm
ConstructN ( i, µ ) requires evaluation of intersec-tions of the cells. For the dyadic grid the cells are cubes and their intersection canbe computed from the cell index ( s , λ ) using (7.1) and (7.2). In practice we do notcompute the intersections. Instead we formally manipulate indices of the cells toobtain the index of the cell corresponding to the intersection. N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 23 Application
Error Estimates for Persistence Diagrams.
The algorithm we have de-scribed for computing persistence diagrams largely depends on constructing ( f, s i )-verified grids at a variety of thresholds. Depending on how many ( f, s i )-verifiedgrids we are able to construct, the following corollary provides an a posteriori bound on the bottleneck distance between our calculated persistence diagrams andthe true persistent diagrams for our function. Corollary 8.1.
Fix a sequence of real numbers { s i } m +1 i = − with the following prop-erties:(1) s − = −∞ and s < min x ∈ X f ( x ),(2) s m +1 = ∞ and s m > max x ∈ X f ( x ).Suppose there exists a sequence of grids {G s i } m +1 i = − such that each grid G s i is ( f, s i )-verified. Define {E s i } m +1 i = − to be the associated CW complexes, define { ¯ X s i } m +1 i = − tobe the corresponding approximate filtration, and define P D (cid:48) to be the persistencediagram associated with { ¯ X s i } m +1 i = − . Then d B ( P D ( f ) , P D (cid:48) ) ≤ ε where ε = max ≤ i ≤ m | s i − s i − | . Proof.
By Theorem 6.4 it follows that each CW structure E s i is compatible with X s i . Moreover, for every cell e ∈ E s i which satisfies Condition (3) in Definition 1.1,the function f ◦ h ( x , s ) is non-increasing in s . It follows from Theorem 1.6 that { ¯ X s i } m +1 i = − and { X s i } m +1 i = − have the same persistence diagrams. For every δ > { X s i } m +1 i = − is an ( ε + δ )-approximate filtration. Taking the limit as δ goesto zero, it follows from Lemma 4.1 that d B ( P D ( f ) , P D (cid:48) ) ≤ ε . (cid:3) In practice, we do not know a priori the thresholds for which we can constructan ( f, s i )-verified grid. Consequently, we can only give a posteriori error estimatesfor our persistence diagrams. For example, fix a sequence of real numbers { s i } m +1 i = − as above and suppose that it is uniformly spaced. That is ∆ = | s i − s i − | for1 ≤ i ≤ m . If we attempt to create ( f, s i )-verified grids corresponding to { s i } m +1 i = − then our resulting persistence diagram will satisfy(8.1) d B ( P D ( f ) , P D (cid:48) ) ≤ ∆( F + 1)where F is the largest number of consecutive thresholds s i which we failed toconstruct ( f, s i )-verified grids. While we have a priori control on choosing ∆, thefinal error estimate can only be determined after finishing the calculation.The algorithm we employ for creating an ( f, t )-verified grid will either terminatesuccessfully or halt after reaching one of two failure conditions. The first failurecondition is reached if the algorithm has exceeded a predetermined maximum res-olution of the domain. This condition is imposed to prevent the algorithm fromsubdividing grid elements ad infinitum . The second failure condition occurs if, byusing interval arithmetic, a vertex has an image which is neither strictly greaterthan nor strictly less than the threshold t . If this conditions is satisfied, then weare unable to decide whether or not to include that vertex in the approximation.In general, we only need to worry about the first failure condition, and regionsof the sub-level set which are geometrically complex. In particular, if a cube C - - - - - - - - Figure 8. (a) Plot of a trigonometric polynomial with randomcoefficients. (b) β persistence diagram and (c) β persistence di-agram.contains a critical point of the function f , then it can only be ( f, t )-verified if itsatisfies either t > f ( C ) or t ≤ f ( C ). Consequently, if f has a critical value veryclose to the threshold t , then a dyadic grid will need to have a very fine resolution inorder to be ( f, t )-verified. Additionally, it is sometimes the case that the offspringof a cube which is ( f, t )-verified are themselves not ( f, t )-verified. For these reasons,having very precise interval arithmetic can help prevent runaway subdivisions ofthe domain.8.2. Computational Results.
We implemented our algorithm for studying thesub-level sets of functions defined on the unit square: [0 , . The source code canbe downloaded from [13]. In particular, we studied random Fourier series of theform f ( x, y ) := N (cid:88) i =1 N (cid:88) j =1 a i,j, sin(2 πix ) sin(2 πjy ) + a i,j, sin(2 πix ) cos(2 πjy )+ a i,j, cos(2 πix ) sin(2 πjy ) + a i,j, cos(2 πix ) cos(2 πjy )where the coefficients a i,j,k were taken from a normal distribution of mean 0 andstandard deviation 1. Depicted in Figure 8 is one such computation, where we used N = 5 as the maximum number of modes.For the filtration we tried to construct ( f, s i )-verified grids for thresholds be-tween −
15 and 15 using spacings of 0 .
05. Out of the 601 thresholds, we failed toconstruct ( f, s i )-verified grids for 11 thresholds. The largest number of consecutivefailures occurred when both − . − .
15 failed to be verified. Hence the largestgap between any two verified thresholds in our filtration was 0 .
15. Therefore thepersistent diagrams we calculated are within a bottleneck distance of 0 .
15 of thetrue persistent diagrams for our given function.The computation took 193 hours and used 16 MB of memory, and was performedon a single Intel i7-2600 processor. The algorithm created a filtered CW structureon [0 , which contained 145,721 cells, which was then processed by the Perseus software program to compute its persistent homology [23]. The (cid:15) -approximatefiltration was constructed sequentially by computing ¯ X s i := X s i ∩ ¯ X s i +1 , so atmost three CW complexes were stored in memory at a given time. At the costof using more memory, one could reduce the computation time by computing the N ε APPROXIMATIONS OF PERSISTENCE DIAGRAMS 25 cellular approximations {X s i } in parallel, and subsequently creating the filtration { ¯ X s i } . Additionally, the one point of departure our implementation took from thealgorithm here described is that the cells were represented using a lookup table, asopposed to the product of generalized intervals.The reason our algorithm failed to verify the 11 thresholds is because when at-tempting to verify them, it surpassed 25 subdivisions of the domain. In general,the size of geometric features in these sub-level sets was much larger than 2 − .However, the error bounds produced by interval arithmetic have a risk of creatinga cascade of cubes which fail to become ( f, t )-verified. In order to compute intervalarithmetic more precisely, whenever we computed a function on a cube, we subdi-vided the cube seven times and computed the function using interval arithmetic oneach piece.When we repeated the computation but instead used six subdivisions, the num-ber of thresholds we failed to verify increased 4 . . . . Acknowledgments
The authors would like to thank Konstantin Mischaikow for his insightful com-ments and conversations.
References
1. Gunnar Carlsson, Vin de Silva, and Dmitriy Morozov,
Zigzag persistent homology and real-valued functions , Proceedings of the Twenty-fifth Annual Symposium on Computational Ge-ometry (New York, NY, USA), SCG ’09, ACM, 2009, pp. 247–256.2. Frederic Chazal, Vin de Silva, Marc Glisse, and Steve Oudot,
The structure and stability ofpersistence modules , arXiv (2013), 1207.3674.3. Gregory Scott Cochran, Thomas Wanner, and Pawe(cid:32)l D(cid:32)lotko,
A randomized subdivision al-gorithm for determining the topology of nodal sets , Computational Methods in Science andEngineering (2013), no. 5, B1034–B1054.4. David Cohen-Steiner, Herbert Edelsbrunner, and John Harer, Stability of persistence dia-grams , Discrete Comput. Geom. (2007), no. 1, 103–120. MR 2279866 (2008i:68130)5. Sarah Day, William D. Kalies, and Thomas Wanner, Verified homology computations for nodaldomains , Multiscale Modeling and Simulation (2009), 1695–1726.6. Pawe(cid:32)l D(cid:32)lotko, Tomasz Kaczynski, Marian Mrozek, and Thomas Wanner, Coreduction homol-ogy algorithm for regular CW-complexes , Discrete & Computational Geometry (2010),361–388.7. Herbert Edelsbrunner and John Harer, Persistent homology-a survey , Contemporary mathe-matics (2008), 257–282.8. Herbert Edelsbrunner and John L. Harer,
Computational topology , American MathematicalSociety, Providence, RI, 2010, An introduction. MR 2572029 (2011e:00001)9. Marcio Gameiro, Konstantin Mischaikow, and William Kalies,
Topological characterization ofspatial-temporal chaos , Phys. Rev. E (2004), 035203.10. Marcio Gameiro, Konstantin Mischaikow, and Thomas Wanner, Evolution of pattern com-plexity in the Cahn-Hilliard theory of phase separation , Acta Materialia (2005).11. Eldon Hansen and G William Walster, Global optimization using interval analysis: revisedand expanded , vol. 264, CRC Press, 2003.12. Allen Hatcher,
Algebraic topology , Cambridge University Press, 2001.13. J. Jaquette and M. Kramar, C++ code available at: http://chomp.rutgers.edu/Projects/Rigorous_Computational_Dynamics/Approximating_Persistence_Diagrams.html .
14. Tomasz Kaczynski, Konstantin Mischaikow, and Marian Mrozek,
Computational homology ,Springer, 2004.15. L. Kondic, A. Goullet, C.S. O’Hern, M. Kramar, K. Mischaikow, and R.P. Behringer,
Topologyof force networks in compressed granular media , Europhys. Lett. (2012), 54001.16. M. Kramar, A. Goullet, L. Kondic, and K. Mischaikow, Persistence of force networks incompressed granular media , Phys. Rev. E (2013), 042207.17. , Quantifying force networks in particulate systems , Physica D (2014), 37–55.18. M. Kramar, R. Levanger, M. Shcatz, M. Paul, M. Xu, J. Tithof, B. Suri, and K. Mischaikow,
Detecting recurrent dynamics of Kolmogorov flow and Rayleigh-Benard convection , (2014),In preparation.19. Kapilanjan Krishan, Huseyin Kurtuldu, Michael Schatz, Marcio Gameiro, Konstantin Mis-chaikow, and Santiago Madruga,
Homology and symmetry breaking in Rayleigh-Benard con-vection: Experiments and simulations , Physics of Fluids (2007).20. H. Kurtuldu, K. Mischaikow, and M. F. Schatz, Measuring the departures from the Boussinesqapproximation in Rayleigh-Benard convection experiments , Journal of Fluid Mechanics (2011), 543–557.21. R. Mendoza, K. Thornton, I. Savin, and P.W. Voorhees,
The evolution of interfacial topologyduring coarsening , Acta Materialia (2006), 743–750.22. Yuriy Mileyko, Sayan Mukherjee, and John Harer, Probability measures on the space of per-sistence diagrams , Inverse Problems (2011), no. 12, 124007.23. Konstantin Mischaikow and Vidit Nanda, Morse theory for filtrations and efficient computa-tion of persistent homology , Discrete & Computational Geometry (2013), no. 2, 330–353.24. Konstantin Mischaikow and Thomas Wanner, Topology-guided sampling of nonhomogeneousrandom processes , Annals of Applied Probability (2010), no. 3, 1068–1097. Department of Mathematics, Hill Center-Busch Campus, Rutgers University, 110Frelinghusen Rd, Piscataway, NJ 08854-8019, USA
E-mail address : [email protected] Department of Mathematics, Hill Center-Busch Campus, Rutgers University, 110Frelinghusen Rd, Piscataway, NJ 08854-8019, USA
E-mail address ::