Pawe

Thomas Wanner

Rigorous cubical approximation and persistent homology of continuous functions

The interaction between discrete and continuous mathematics lies at the heart of many fundamental problems in applied mathematics and computational sciences. In this paper we discuss the problem of discretizing vector-valued functions defined on finite-dimensional Euclidean spaces in such a way that the discretization error is bounded by a pre-specified small constant. While the approximation scheme has a number of potential applications, we consider its usefulness in the context of computational homology. More precisely, we demonstrate that our approximation procedure can be used to rigorously compute the persistent homology of the original continuous function on a compact domain, up to small explicitly known and verified errors. In contrast to other work in this area, our approach requires minimal smoothness assumptions on the underlying function.

Introduction

In many applied situations it is necessary to replace a given smooth function by a discretized version which in some sense is close to the given original mapping. In this paper we discuss a number of algorithms for approximating a continuous function f : D → R m , defined on a compact rectangular domain D ⊂ R n , using piecewise constant functions. In order to achieve approximations with mathematically verified approximation bounds, our approach relies on the use of rigorous computer arithmetic. One such example is interval arithmetic [START_REF] Moore | Introduction to Interval Analysis[END_REF], and it will form the foundation of the algorithms described in the present paper. However, our approach can also be used with any other rigorous computer arithmetic which provides rigorous enclosures of function values. For example, our algorithms can readily be adapted to employ alternative tools such as affine arithmetic [START_REF] De Figueiredo | Affine arithmetic: Concepts and applications[END_REF], generalized interval arithmetic [START_REF] Hansen | A generalized interval arithmetic[END_REF], or in fact, any other rigorous arithmetic capable of providing range enclosures.

Upon successful completion, the algorithms presented in this paper will return a decomposition of the rectangular domain D into a collection of compact n-dimensional rectangles, any two of which only intersect in a subset of their topological boundaries. Furthermore, each rectangle R in the collection is assigned a value val(R) in R m . The approximating function f is then defined in the following way. For every rectangle R in the constructed decomposition, the function value of f in the interior of R is given by val(R). On the other hand, if x ∈ D lies in the intersection of k rectangles R 1 , . . . , R k , then we define f (x) = min i=1,...,k val(R i ). In this way, the decomposition of D produced upon successful completion of the algorithm gives rise to a lower semi-continuous piecewise constant approximation f of the given continuous function f . We will see later that the algorithm in fact guarantees the inequality f -f ∞ < ε , where ε > 0 denotes the input parameter of the algorithm, and where • ∞ denotes the usual maximum norm. This procedure is illustrated in the two panels of Figure 1. In the left panel, the approximating function f is shown for the one-dimensional function f (x) = 2 -25x + 108x 2 -162x 3 + 81x 4 on the interval D = [0, 1], and with approximation parameter ε = 1/2. The right panel contains the approximation for the so-called Ackley function, which is one of the standard test functions in nonlinear optimization and which will be discussed in more detail later in this paper.

The approximation algorithms presented in this paper can in principle be applied in a number of situations. For example, in the context of global optimization they allow one to rigorously narrow down the search for global minima or maxima of a real-valued function, as can easily be seen in the examples of Figure 1. Rather than pursuing such applications, we confine ourselves to the computation of the persistent homology of a function. One-dimensional persistent homology was introduced in [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF], and since then has become one of the main topological tools in the applied sciences for studying the shape of functions or spaces. Its popularity can be attributed to a variety of factors:

• Persistent homology provides quantitative information about the object of interest, including information about its shape,

• the obtained information is both metric and coordinate free, and

• the information is preserved under continuous deformations, e.g., it is invariant with respect to isometries. • Persistence captures topological information at different scales,

• it is robust with respect to the influences of random noise, and

• efficient algorithms are available for its computation.

At present, persistent homology can only be determined algorithmically for the case of scalar-valued functions, i.e., only one-dimensional persistence is amenable to a computational treatment. For this, one needs to find a suitable filtration of a finite cell complex, and the results of this paper will enable its efficient construction. Due to deep results from quiver theory [START_REF] Carlsson | The theory of multidimensional persistence[END_REF], there are no finite invariants which completely characterize multi-dimensional persistent homology. Nevertheless, there are many groups working on providing partial invariants for multi-dimensional persistent homology that can be used in practice.

The algorithms presented in this paper for the case of functions f : D ⊂ R n → R m provide discretizations of f with rigorous error bounds. These discretizations can immediately be used in the computation of one-dimensional persistent homology, i.e., in the scalar-valued case m = 1. We would like to point out, however, that the approximations can potentially also be used in the context of recently developed algorithms which compute partial invariants of multi-dimensional persistent homology. In order to specifically accommodate the computation of persistent homology, we consider specific approximations in which, in addition to the rectangular decomposition of the domain D, we also keep track of all the intersections of the rectangles in the decomposition. This is achieved by constructing a weighted cell complex which will later be used to compute persistent homology or some invariants of multi-dimensional persistent homology.

Approximating continuous functions is a part of the classical simplicial approximation theorem [START_REF] Wylie | Homology Theory. An Introduction to Algebraic Topology[END_REF]. However, we are not aware of any attempts to create lower semi-continuous approximations of continuous functions using rigorous computer arithmetic. In the context of differentiable functions, the work [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF] comes close in spirit to the present paper. The authors use a generalization of quad-trees and oct-trees to compute compatible cellular decompositions for level sets of a given differentiable function. Then, after computing the generators of the associated homology groups at each level of the filtration, together with their inclusions, they obtain the persistent homology of a differentiable function. It is worth mentioning that through the computation of the persistence of inclusions between two fixed level sets, the authors of [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF] obtain an exact persistence result for the chosen levels. Similar idea to compute standard persistent homology was used in [START_REF] Mrozek | Coreduction homology algorithm for inclusions and persistent homology[END_REF]. This, however, requires the successful termination of their algorithm. Unfortunately, even in the examples presented in their paper, their algorithm cannot validate a significant number of levels, which is due to grid alignment issues which were first pointed out in [START_REF] Cochran | A randomized subdivision algorithm for determining the topology of nodal sets[END_REF]. They arise naturally in the setting of [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF], since the authors have to use binary subdivisions for the persistence computations, and since varying the threshold level ensures that the level curves will frequently come close to the binary subdivision lines. In the context of one fixed level set, this effect has been one of the shortcomings of the method in [START_REF] Day | Verified homology computations for nodal domains[END_REF], but it could be removed through randomized subdivisions in [START_REF] Cochran | A randomized subdivision algorithm for determining the topology of nodal sets[END_REF]. Unfortunately, the simple randomization trick cannot be applied in the setting of [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF]. Moreover, the total computational cost of the method in [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF] may in some cases be higher than the one of our method presented below. This is due to the fact that in [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF] the authors need to exactly verify every intermediate level set, and then perform independent homology computations at each filtration level.

The technique presented in this paper is of more general use than just for the computation of persistent homology, which was the only focus of the technique presented in [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF]. Our main goal is to provide a constructive version of the simplical approximation theorem. If, however, we restrict our attention only to computations of persistent homology, the technique presented in the present paper turns out to be orthogonal to the technique introduced in [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF]. The latter one is expected to work very well in the case of functions whose first derivatives only change moderately, especially if the considered functions have a large range. In this case one would expect to obtain the approximation of most of the sublevel sets of the function in a fast and easy way. If, on the other hand, the function exhibits significant variations within a limited range, then using the technique introduced in [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF] seems to be quite challenging. In this latter case, the technique presented in the present paper still works well, since it is based on function's range enclosures rather than the properties of its derivatives.

Other approaches to approximating the persistence of functions, which also take errors into account, were presented in [START_REF] Bendich | Computing robustness and persistence for images[END_REF] and [START_REF] Bendich | Persistent homology under non-uniform error[END_REF]. In the former paper, the authors employ oct-trees to approximate the persistence of a digital image as accurately as possible for a fixed number of simplification steps. In fact, one of our algorithms, the one presented in Section 3.5, is directly inspired by this approach. Finally, the paper [START_REF] Bendich | Persistent homology under non-uniform error[END_REF] is concerned with persistent homology computations in the presence of non-uniform error models.

The remainder of this paper is organized as follows. In Section 2 we give a short introduction to rectangular CW-complexes, homology and (multi-dimensional) persistent homology, and we survey necessary results on the stability of persistent homology. Section 3 contains the main results of the paper. After collecting some basic results from interval arithmetic in Section 3.1, we then turn our attention to the main algorithms. In Section 3.2 we begin with a method for computing a piecewise constant approximation of a continuous function f : D ⊂ R n → R m , and discuss termination criteria as well as complexity issues for a large class of functions f . This is followed in Section 3.3 by a rectangle subdivision algorithm which, during the subdivision process, keeps track of all boundary elements, and which establishes a lower semi-continuous approximation. In the end, this method allows us to construct a cell complex which is suitable for the computation of persistence. This will be outlined in detail in Section 3.4, where we use the subdivision algorithm to first construct a filtered complex, and then its (multi-dimensional) persistent homology. It will be shown that the resulting persistence diagram is close to the persistent homology of the underlying function f . A second approximation algorithm, which is motivated by the approach in [START_REF] Bendich | Computing robustness and persistence for images[END_REF], is the subject of Section 3.5. This time, we aim to approximate the persistent homology of f as effectively as possible using a limited number of subdivisions. Finally, Section 4 contains a number of numerical case studies.

Persistent Homology

In this section we provide a brief introduction to persistent homology, in order to keep the present paper as self-contained and accessible as possible. For a more-detailed introduction to standard algebraic topology we refer the reader to [START_REF] Massey | A Basic Course in Algebraic Topology[END_REF], an extensive introduction to applied and computational topology can be found in [START_REF] Edelsbrunner | Computational Topology[END_REF].

Rectangular CW-Complexes

In the following, the term interval is always used for a compact interval I = [a, b] ⊂ R with a ≤ b. We say that the interval I is degenerate if a = b, otherwise it is called nondegenerate. In order to describe the algebraic boundary of intervals, we say that a degenerate interval has no faces, while the nondegenerate interval I =

. A rectangle in R d is a product Q = Q 1 × Q 2 × . . . × Q d of d intervals Q 1 , Q 2 , . . . , Q d . The dimension of Q, abbreviated as dim Q, is the number of nondegenerate intervals among Q 1 , Q 2 , . . . , Q d . Consider now two rectangles Q = Q 1 × . . . × Q d and P = P 1 × . . . × P d .
Then the rectangle P is called a primary face of Q if P ⊂ Q, the dimensions satisfy dim P = dim Q -1, and if there exists an index j such that P j is a face of Q j , in the sense defined above for intervals. Notice that in this case, the index j is uniquely determined, and we have dim P i = dim Q i for all i = j. We would like to point out, however, that it is not required that the identity P i = Q i holds for any of these indices -only the inclusions P i ⊂ Q i are necessary. Finally, the rectangle P is called a face of Q, if either there is a descending sequence of primary faces joining Q to P , or if P is the empty set.

We now turn to the definition of the boundary of a rectangle. For an interval I = [a, b], its boundary is defined as bdr I = {a, b} if a = b, and as bdr I = ∅ if a = b. The boundary of a higher-dimensional rectangle

Q = Q 1 × . . . × Q d is defined as the union bdr Q = d i=1 Q 1 × . . . × bdr Q i × . . . × Q d . (1)
One can easily see that the boundary bdr Q corresponds to boundary of the dim Q-dimensional manifold Q.

The above-defined rectangles are the building blocks for more general structured sets. For the purposes of this paper, we consider a rectangular CW-complex, which can be defined as follows. A rectangular structure is a finite collection Q of rectangles in R d such that for any P, Q ∈ Q we either have P ∩ Q = ∅, or P ∩ Q is a common face of both P and Q which in addition belongs to Q. Furthermore, we assume that for any rectangle Q ∈ Q, its above-defined boundary bdr Q satisfies Figure 2: An example of a rectangular CW-complex. The first image shows the initial collection Q (0) of 46 zero-dimensional rectangles, the middle image depicts the 70 one-dimensional rectangles contained in Q (1) , and the last image shows the 25 two-dimensional rectangles in Q (2) which fill in the holes. The final rectangular CW-complex is a square. the equality bdr Q = {P ∈ Q : P is a primary face of Q} .

(

) 2
We would like to point out that while (1) defines the boundary of a rectangle, the condition stated in [START_REF] Bendich | Computing robustness and persistence for images[END_REF] ensures that in a rectangular CW-complex, the usually finer decomposition of the boundary still covers the complete boundary of Q through primary faces. For example, consider the top left rectangle in the middle image in Figure 2. While the top, left, and bottom edges are represented by one primary face each, the right edge of the rectangle is the union of three primary faces. Finally, a rectangular CW-complex |Q| is given by the union of some rectangular structure Q, i.e., it is the subset of Euclidean space which is occupied by the rectangles in Q. Any 0-dimensional rectangle in Q is called a vertex , and if Q ∈ Q is an n-dimensional rectangle, then the set Q \ bdr Q is called an n-cell . In order to illustrate this definition further, we note that the rectangles in a rectangular structure Q satisfy the following two properties.

• Let Q (k) denote the rectangles in Q of dimension k. Then for k ∈ N 0 and rectangles P, Q ∈ Q (k) we either have P ∩ Q = ∅, or the intersection P ∩ Q is a face of both P and Q.

• For k ∈ N and P ∈ Q (k) , the union of all rectangles in Q (k-1) which are primary faces of P equals the boundary bdr P .

Notice that the above description can be used to give an algorithmic description of building Q.

Initially, one starts with a collection of zero-dimensional rectangles, i.e., points, in R d . Then one adds one-dimensional rectangles connecting some of the points, adds two-dimensional rectangles filling holes, etc. For an example of a rectangular CW-complex, see Figure 2. The depicted CWcomplex decomposes a square into smaller rectangles through selected dyadic subdivisions. In the remainder of this paper, we will only consider rectangular CW-complexes of this type, even though the above definition is more general. This is due to the fact that the algorithms presented below avoid the use of derivatives, which would allow for different adaptive decompositions.

Persistent Homology of Rectangular CW-Complexes

We now turn our attention to the definition of homology. It was shown in [START_REF] Kaczynski | Coreduction Homology Algorithm for Regular CW-Complexes[END_REF]Theorem 3.9] that every rectangular CW-complex is a CW-complex in the classical sense, see for example [START_REF] Massey | A Basic Course in Algebraic Topology[END_REF]. This implies that if a rectangular CW-complex defines a polyhedron, then its homology is isomorphic to its singular homology. In order to further define and then efficiently compute the homology of such a rectangular CW-complex, we make use of the notion of incidence coefficients or incidence numbers. In the simplest case of homology with coefficients in Z 2 , the incidence number of two rectangles Q, P ∈ Q satisfies κ(Q, P ) = 1 if and only if P is a primary face of Q, and we define κ(Q, P ) = 0 in all other cases. Even for the more involved case of homology with coefficients in a field with at least three elements, the incidence numbers of a rectangular CW-complex can still be computed efficiently, for more details see [START_REF] Kaczynski | Coreduction Homology Algorithm for Regular CW-Complexes[END_REF]Theorem 3.10]. Suppose now that X ⊂ R d denotes a rectangular CW-complex with rectangular structure Q. In order to rigorously define the homology of X or Q, we begin by establishing the notions of chains, cycles, and boundaries. Let F denote an arbitrary coefficient field. Then a k-chain c is a formal sum of the form

c = A∈Q (k) α A A , where α A ∈ F and Q (k) = {Q ∈ Q : dim Q = k} .
The set of all k-chains forms an additive group, which is called the k-th chain group and denoted by

C k (Q). For a rectangle Q ∈ Q (k) , the boundary ∂Q of Q is the (k -1)-chain given by ∂Q = P ∈Q (k-1)
κ(Q, P )P , which in view of the above discussion of the incidence numbers κ(Q, P ) is in fact only a sum over the primary faces of Q, see also [START_REF] Kaczynski | Coreduction Homology Algorithm for Regular CW-Complexes[END_REF]Section 3.1]. Through linear extension, this defines a homomorphism ∂ :

C k (Q) → C k-1 (Q),
which is called the boundary operator. One can show that the boundary operator satisfies

∂ 2 = ∂ • ∂ = 0 : C k (Q) → C k-2 (Q) . A chain c ∈ C k (Q) is called a k-cycle, if ∂c = 0.
The set of all k-dimensional cycles forms an additive group which is denoted by Z k (Q). A k-cycle d is called a k-boundary, if there exists a chain e ∈ C k+1 (Q) such that ∂e = d holds. The set of all k-boundaries is abbreviated by B k (Q), and one can easily see that it is a subgroup of Z k (Q). Finally, the homology group of Q in dimension k is then the quotient group

H k (Q) = Z k (Q)/B k (Q).
From a geometric point of view, the k-th homology group H k (Q) measures the number of k-dimensional holes in Q, and we refer the reader to [START_REF] Edelsbrunner | Computational Topology[END_REF][START_REF] Massey | A Basic Course in Algebraic Topology[END_REF] for more details.

As our last concept of this section, we now consider persistent homology. While we keep our discussion brief, a more extensive introduction can be found in [START_REF] Lotko | Topological microstructure analysis using persistence landscapes[END_REF]. Consider a rectangular CWcomplex given by the rectangular structure

Q. Then a subcollection Q 0 of Q is called a subcomplex of Q, if the set Q 0 defines a rectangular structure in its own right. In addition, a sequence of subcomplexes ∅ ⊂ Q 0 ⊂ . . . ⊂ Q n-1 ⊂ Q n = Q is called
a filtration of the rectangular CW-complex given by Q. Filtrations lie at the heart of persistent homology, which basically tracks topological changes in the CW-complexes of the filtration as the index increases. While filtrations of CWcomplexes can arise in a number of ways, the following situation will be central for the present paper. Consider a function f : Q → R which assigns a real value to every rectangle in Q. Moreover, assume that f is increasing in the sense that f (P ) ≤ f (Q) whenever P ∈ Q is a face of Q ∈ Q. One can easily see that if a 0 ≤ a 1 ≤ . . . ≤ a n are the finitely many real numbers in the range of f , then the sets

Q k = {Q ∈ Q : f (Q) ≤ a k }
are rectangular structures, and therefore

Q 0 ⊂ . . . ⊂ Q n-1 ⊂ Q n =
Q is a filtration of the rectangular CW-complex given by Q. We refer to this setting as a filtered rectangular CW-complex |Q| given by (Q, f ).

Consider now an arbitrary filtration

Q 0 ⊂ . . . ⊂ Q n-1 ⊂ Q n of rectangular CW-complexes. Then for 0 ≤ a ≤ b ≤ n, we have Q a ⊂ Q b , which in turn implies C k (Q a ) ⊂ C k (Q b ) for all k ∈ N 0 .
This inclusion of the chain groups can be used to define homomorphisms of the associated homology groups. More precisely, we define the mapping ι k a,b :

H k (Q a ) → H k (Q b ) by mapping any representa- tive k-chain c of a homology class in H k (Q a ) to the homology class generated by c in H k (Q b ).
One can show that this is well-defined, and that the mappings ι k a,b are group homomorphisms. Then the k-th persistent homology groups are defined as

H a,b k (Q) = Im ι k a,b ,
i.e., as the images of the homomorphisms ι k a,b . As mentioned above, persistent homology tracks topological changes through the filtration, and this leads to the following notions for homology classes. We say that a homology class

γ ∈ H k (Q ) is born in H k (Q ), if γ ∈ Im ι k -1, . The homology class γ dies entering H k (Q m )
, if either it became trivial, or if it merges with an older class in Q m . In this context, a class is called older, if it was born earlier than the homology class γ. The index values and m are called the birth time and death time of the homology class γ, respectively. We refer the reader to [START_REF] Edelsbrunner | Computational Topology[END_REF] for more details.

The collection of all persistent homology groups is usually referred to as the persistence module. Despite its abstract nature, it can be completely encoded via the notions of birth and death of homology classes. For every homology class γ, we record its birth time and its death time as an ordered pair, where the death time is set to

∞ if the class γ is nontrivial in Q n . In this way, a multiset of pairs of numbers {(a i , b i )} i∈I is created, where -∞ < a i ≤ b i ≤ ∞.
This multiset of points naturally lies above the line b = a. Such a representation of the persistence module is called a persistence diagram. Another, more intuitive, demonstration of the persistent homology of a cubical complex can be found in Figure 3.

If the persistence module is for a filtered rectangular CW-complex X = |Q| given by (Q, f ), then we denote the persistence diagram by D(f, X), or shorter by D(f ) if the underlying space X is clear from context. Examples of persistence diagrams can be found in the Figure 4. In the top left image, a continuous function f is shown together with its piecewise constant approximation f . Now let Q (0) denote the endpoints of the interval [0, 1], together with the collection of all points of discontinuities of f . If we denote the points in the set Q (0) by 0 = x 0 < . . . < x N = 1 and let

Q (1) = {[x k-1 , x k ] : k = 1, . . . , N }, then Q = Q (0) ∪ Q (1)
is a rectangular structure, which renders the interval X = [0, 1] a filtered rectangular CW-complex given by (Q, f ). The associated persistence diagram is shown in the upper right panel. The two panels below show persistence diagrams for two more approximations of f , which will be discussed more later.

Stability of Persistent Homology

One of the prominent features of persistent homology is the fact that the collection of persistence diagrams can be turned into a metric space in various ways, and we will present two of those below. In order to set the stage, we have to slightly modify the notion of persistence diagram. So far, a persistence diagram is a finite multiset of points (a, b) with a < b, which represent birth/deathtime pairs for homology generators. For reasons that will become clearer momentarily, we now add infinitely many of copies of all points (a, a) for a ∈ R to every persistence diagram, i.e., in addition 3: From left to right, the images show a filtration on a simple 3 by 3 cubical complex (which is a special case of a rectangular CW complex). At each stage in the filtration, the black cubes are part of the complex. The numbers from 0 to 4 below the images denote the value of the filtering function. Below the filtering function the corresponding persistence intervals in dimension 0 (d0) and dimension 1 (d1) are presented. At level 0, both the complex and the collection of persistence intervals are empty. For complex 1, the four added cubes constitute four connected components, indicated by as of yet unbounded persistence intervals in dimension 0. Upon moving to complex 2 in the filtration, two pairs of the connected components merge, and therefore two persistence intervals which were created at level 1 now terminate. With complex 3, only one connected component remains, and therefore another one of the two unbounded intervals in dimension 0 terminates at time 3. At the same time, a homology class in dimension 1 is created, as indicated by the new and as of yet unbounded persistence interval in dimension 1. When we reach the final complex 4, the homology class in dimension 1 becomes trivial, which is indicated by the termination of the unique persistence interval in dimension 1. Since the filtration does not evolve any further, the only remaining unbounded persistence interval in dimension 0 stays infinite. to the birth/death pairs every persistence diagram contains infinitely many copies of the diagonal b = a. This modification ensures that any two persistence diagrams A and B have the same infinite cardinality, and there always are bijections η : A → B. Notice that without the addition of the diagonal, such bijections would generally not exist. Now we can define the two central notions of distance between A and B.

• The bottleneck distance between two persistence diagrams A and B is defined as

W ∞ (A, B) = inf sup p∈A p -η(p) ∞ : η : A → B is a bijection .
• The degree-q Wasserstein distance between two persistence diagrams A and B is defined as

W q (A, B) = inf        p∈A p -η(p) q ∞   1/q : η : A → B is a bijection      .
In other words, both distances try to identify a transport plan η : A → B which transforms the persistence diagram A into the persistence diagram B, with in some sense minimal movement of the points in the diagrams. In the case of infinite persistence intervals of a form p = (b p , +∞) and q = (b q , d q ), we define p -q ∞ = |b p -b q | if d q = +∞, and p -q ∞ = +∞ otherwise. While the bottleneck distance only considers the largest displacement, the Wasserstein distance takes an average over all displacements. See also Figure 5. Now that we are able to measure changes in persistence diagrams, one can try to assess how changes in a filtration function f affect the associated persistence diagrams. In fact, it is possible to establish a wide-ranging stability property for persistent homology in this context, which is one of the central reasons for its usefulness in the applied sciences. Consider a function f : X → R which is defined on a triangulable topological space X. One can define persistent homology for the infinite filtration created by sublevel sets of f , as long as the function f satisfies an additional property. To describe this in slightly more detail, let X a = f -1 (-∞, a] for all a ∈ R. Then we clearly have the inclusion X a ⊂ X b whenever a ≤ b. As before, one would like to define the persistent homology groups as images of the maps on homology induced by the inclusions X a → X b . Since we are faced with a filtration with infinite index set, certain finiteness assumptions have to be made. For this, one calls a number a ∈ R a homological critical value, if there is no value ε > 0 such that the inclusion X a-ε → X a+ε induces only isomorphisms on the homology level. In other words, at a homological critical value the homology groups change. We then call the function f tame, if it only has finitely many homological critical values, and if all homology groups of all sublevel sets X a have finite rank. For tame functions, it is possible to define persistence diagrams similar to the procedure described above, see for example [18, p. 182]. Moreover, the following stability result was established by Cohen-Steiner, Edelsbrunner, and Harer [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF].

Theorem 2.1. Let X be a triangulable space, and let f, g : X → R denote two continuous tame functions. Then the persistence diagrams D(f ) and D(g) of f and g, respectively, satisfy

W ∞ (D(f ), D(g)) ≤ f -g ∞ .
While the above theorem was a milestone for the applicability of persistent homology, it did suffer from a serious drawback. In its original formulation, the functions f and g have to be at least continuous. This of course would preclude our idea of approximating a continuous functions via a piecewise constant one. Fortunately, this continuity restriction has been removed in the work by Chazal et al. [START_REF] Chazal | Proximity of persistence modules and their diagrams[END_REF][START_REF] Chazal | The Structure and Stability of Persistence Modules[END_REF], which is based on interleavings of persistence modules. As a consequence, we have the following result.

Theorem 2.2. Let X be a rectangular CW-complex, and let f, g : X → R denote two functions, each of which is either continuous and tame, or piecewise constant on the rectangles of the rectangular structure associated with X and lower semi-continuous. Then the persistence diagrams D(f ) and D(g) of f and g, respectively, satisfy

W ∞ (D(f ), D(g)) ≤ f -g ∞ .
In other words, small perturbations in the function f only lead to small perturbations of the persistence diagram D(f ) with respect to the bottleneck distance.

One of the crucial properties of persistent homology is its inherent computability and a number of software packages have been developed for this, for example dionysus [START_REF] Morozov | Dionysus[END_REF], PHAT [START_REF] Bauer | Persistent Homology Algorithm Toolbox)[END_REF], Perseus [START_REF] Nanda | Perseus, the Persistent Homology Software[END_REF] and Gudhi [START_REF]GUDHI: User and Reference Manual[END_REF]. While two of these are based on standard matrix reductions [START_REF] Bauer | Persistent Homology Algorithm Toolbox)[END_REF][START_REF] Morozov | Dionysus[END_REF], one is based on discrete Morse theory [START_REF] Nanda | Perseus, the Persistent Homology Software[END_REF], see also [START_REF] Mischaikow | Morse theory for filtrations and efficient computation of persistent homology[END_REF]. All of the packages are designed to compute the persistence of a finite cell complex and are therefore not capable of solving the problem of computing persistence of a continuous function of a compact set to within a pre-specified error bound. This gap in the literature is filled by the present paper.

Persistent homology can also be defined for vector-valued filtering functions f : X → R n , where X is a rectangular CW-complex or a triangulable space, and this leads to the concept of multi-dimensional persistence. However, for multi-dimensional persistence one can show that there is no finite and complete invariant [START_REF] Carlsson | The theory of multidimensional persistence[END_REF] which is similar to the persistence intervals in the onedimensional case f : X → R. Needless to say, this implies that there is no hope for an algorithm which computes multi-dimensional persistence. To address this issue, a number of research groups have tried to find some descriptors of multi-dimensional persistence, which -though incompletecan still prove to be useful in practice and which are computable. It seems likely that in order to compute such invariants it is necessary to construct a finite cell complex. In fact, our algorithms presented below will work for any continuous vector-valued function which is defined on a rectangular domain, and they produce a piecewise constant lower semi-continuous approximation of the function. Given such an approximation, one can use various stability results for multi-dimensional persistence to show that the multi-dimensional persistence information of the approximation is close to the multi-dimensional persistence of the given continuous function.

To close this section, we present one such stability result for multi-dimensional persistence which is taken from the paper [START_REF] Cerri | Betti numbers in multidimensional persistent homology are stable functions[END_REF]. Their construction is based on finding a matching distance between rank invariants of multi-dimensional persistent homology. For this, the following concepts are introduced in [START_REF] Cerri | Betti numbers in multidimensional persistent homology are stable functions[END_REF]. Let u, v ∈ R n be given, and write u = (u 1 , . . . , u n ) and v = (v 1 , . . . , v n ). We say that u v if for every i ∈ {1, . . . , n} we have

u i ≤ v i . Define the set δ + = {(u, v) ∈ R n × R n : u
v}, and consider the usual maximum norm • ∞ on R n . For a fixed triangulable space X and a function f : X → R n , define the sublevel set X(f, u) = {x ∈ X : f (x) u}. Now consider two vectors u v, and let π u,v k : H k (X(f, u)) → H k (X(f, v)) denote the homomorphism induced by inclusion. Then the image of π u,v k is called the multi-dimensional k-th persistent homology group of (X, f ) at (u, v), where the authors of [START_REF] Cerri | Betti numbers in multidimensional persistent homology are stable functions[END_REF] use Čech homology over a field, see [START_REF] Massey | A Basic Course in Algebraic Topology[END_REF] for more details. Additionally, the k-th rank invariant of a pair (X, f ) over a fixed coefficient field K is defined as the function ρ(X, f, k) :

δ + → N ∪ {∞} given by ρ(X, f, k)(u, v) = rank π (u,v) k
. For the so-defined rank invariant, the following stability theorem was established in [START_REF] Cerri | Betti numbers in multidimensional persistent homology are stable functions[END_REF]Theorem 4.4].

Theorem 2.3. Let X be a triangulable space, let f, g : X → R n be continuous, and let k ∈ Z be arbitrary. Then there exists a distance function D match on the collection of all possible rank invariant functions such that

D match (ρ(X, f, k), ρ(X, g, k)) ≤ max x∈X f (x) -g(x) ∞ .
As before, Theorem 2.3 provides a bound on the change of the rank invariant of multi-dimensional persistent homology in terms of the change of the underlying function f measured in the maximum norm. We would like to point out that the above results is just one example of a number of similar results. In the same paper [START_REF] Cerri | Betti numbers in multidimensional persistent homology are stable functions[END_REF], the authors derive a similar stability result for persistent Betti numbers, and in [START_REF] Cerri | Hausdorff stability of persistence spaces[END_REF] the case of multi-dimensional persistence spaces was considered. In all of these examples, the right-hand side of the estimate is of the form presented in Theorem 2.3. It is therefore clear that bottleneck-type stability can be observed in multi-dimensional persistence, and the discretization process presented below may be useful for computations of some of the above invariants.

Verified Cubical Approximation of Continuous Functions

In this section, we present the main algorithms for approximating continuous functions by piecewise constant lower semi-continuous functions. Our approach makes use of rigorous computations via interval arithmetic, and will allow for the approximation of the persistent homology of a function with rigorous error bounds.

Rigorous Numerics and Interval Arithmetic

We begin by collecting a few definitions and results from interval arithmetic, which are necessary for our approach. In what follows, we will use interval arithmetic to obtain rigorous mathematical statements about approximations of continuous functions. We would like to point out, however, that alternative rigorous computer arithmetic methods could be used just as well, see for example [START_REF] De Figueiredo | Affine arithmetic: Concepts and applications[END_REF][START_REF] Hansen | A generalized interval arithmetic[END_REF]. For a comprehensive introduction to interval arithmetic we refer the reader to [START_REF] Moore | Introduction to Interval Analysis[END_REF].

Throughout this paper our computations work with a fixed set R of real numbers which are representable in a computer, called the representable numbers R. In order to account for rounding errors, rather than working with numbers we make use of intervals. This leads to the concept of interval arithmetic, which represents every real number in R as an interval with endpoints in R. More precisely, if x ∈ R is a representable number, it is represented by the degenerate interval [x, x], while any number x ∈ R \ R is represented as [x, x] with x being the largest representable number which is smaller than x, and similarly for x.

In order to work with functions, one has to define arithmetic operations on intervals. This has to be done in such a way that any potential rounding errors are taken into account. To describe this in more detail, let I(R) denote the set of all intervals with both endpoints in the set of representable numbers. Given two intervals a, b ∈ I(R) with a = [a 1 , a 2 ] and b = [b 1 , b 2 ], and denoting any of the four basic arithmetic operations by , the analogous interval operation a b has to be defined in such a way that {x y : x, y ∈ R and x ∈ a, y ∈ b} ⊂ a b ∈ I(R) .

In practice, this is achieved in computer implementations of interval arithmetic by outward rounding to the next representable numbers. In this way, the resulting interval a b always contains all possible outcomes of the actual operation. In addition to the four basic arithmetic operations, implementations of interval arithmetic also have to provide interval versions of the standard elementary functions. This is usually accomplished via Taylor series expansions, where the resulting interval answer will also include any potential truncation errors. For a more comprehensive treatment of these issues we refer the reader to [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]. Common to all practical computer implementations of interval arithmetic is that they provide a collection of predefined continuous functions, which usually at least contains functions in Φ = {abs, sqrt, sqr, exp, ln, sin, cos, arctan}. For bounds on the complexity of our algorithms, we need to assume later on that the functions of interest are at least Lipschitz continuous. Recall that a function f :

D ⊂ R n → R m is called Lipschitz continuous on the domain D 0 ⊂ D if there exists a constant L ≥ 0 such that f (x 1 ) -f (x 2 ) ∞ ≤ L x 1 -x 2 ∞ for all x 1 , x 2 ∈ D 0 .
Clearly, every elementary function from the set Φ is Lipschitz continuous, at least on compact rectangles over which the function is differentiable. In fact, the same statement is true for any arithmetical expression which is formed using functions from Φ, the four basic arithmetic operations, and composition of functions. For a more formal definition, see [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]Section 1.4].

Despite the fact that arithmetical expressions usually lead to Lipschitz continuous functions, this does not immediately carry over to our machine computations. In order for our results below to hold, we need to make sure that if f is a Lipschitz continuous function, then its actual interval implementation f as a function from a subset of R n into R m is Lipschitz continuous as well. This is the subject of the following result from [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF]Theorem 2.1.1]. Note that this theorem holds for any implementation of interval arithmetic. Theorem 3.1 (Arithmetical Expressions Inherit Lipschitz Continuity). Let f denote an arithmetical expression in n variables, which takes m-dimensional values, i.e., suppose f : D ⊂ R n → R m , and suppose that f is Lipschitz continuous on D 0 ⊂ D with constant L. Furthermore, let f denote its interval implementation as described above, and let D 0 denote the interval vectors in I(R) n which represent the points in D 0 . Then f is Lipschitz continuous on D 0 . More precisely, we have

f (x 1 ) -f (x 2 ) ∞ ≤ α(L) x 1 -x 2 ∞ for all x 1 , x 2 ∈ D 0 ,
where the constant α(L) depends on the arithmetical expression for f , and can be computed explicitly using [29, Table 2.1].

The above Theorem 3.1 implies that interval evaluations of a Lipschitz continuous function determine again a Lipschitz continuous function. Of course, one would expect that the Lipschitz constants of the new function f would increase in size, as the implementation deals with intervals. But as shown in [29, Table 2.1], this effect is minor.

Piecewise Constant Approximations

In this section we show how a continuous function f : R → R m , where R ⊂ R n is a rectangular domain, can be approximated by a piecewise constant function. The approximation will be constructed as a collection of compact rectangles whose union is R, where each of the rectangles in the decomposition is carrying the following information:

• Each rectangle is described by an n-vector of pairs of representable numbers {(b i , e i )} n i=1 which satisfy the inequalities b i ≤ e i . The interval [b i , e i ] is then the projection of the considered rectangle to the i-th coordinate. for every i ∈ {1, . . . , n}.

• Each rectangle has an associated value, which is an m-dimensional vector of representable numbers. For a rectangle Q in the final decomposition of R, its value is denoted by val(R).

The algorithm which leads to the desired decomposition of R is motivated by the classical simplicial approximation theorem. In its original form, this theorem states that every continuous function between triangulable spaces can be approximated by a simplicial map. In order to obtain this simplicial map, one has to subdivide the simplices, possibly many times. As far as we are aware, there is no criterion that allows to check if the approximation of a desired accuracy has been achieved. In this way, the simplicial approximation theorem is not constructive. The constructive algorithm presented below uses a similar idea. By employing rigorous arithmetic, for example interval arithmetic as discussed in the previous section, we can ensure that the final Algorithm 1 Piecewise constant approximation of continuous functions.

Input: Compact rectangle R ⊂ R n , continuous function f : R → R m , error bound ε > 0. Output: Rectangular partition of R which gives an ε-approximation of f on R. Queue Q; Put R to Q; while Q is not empty do Pop rectangle Q from Q; v =componentwise midpoint( f (Q)); if (max( f (Q)) -v ≤ ε) and (v -min( f (Q)) ≤ ε) then val(Q) = v; else
Subdivide Q in all directions. Put all resulting rectangles in Q; return subdivision of the rectangular domain R is such that the value val(Q) of every rectangle Q in the final subdivision of R satisfies the following estimate

f (x) -val(Q) ∞ ≤ ε for all x ∈ Q .
To achieve this property, we repeatedly subdivide R until the above estimate holds for every rectangle in the subdivision. By a minor abuse of notation, we sometimes let ε denote a vector of length n which consists of entries equal to ε. In this case, all inequalities hold componentwise.

To describe the algorithm in more detail, we need to introduce some notation. Consider a rectangle R ⊂ R n . Then we say that the rectangles R 1 , . . . , R k form a partition of the rectangle R if we have R = R 1 ∪ . . . ∪ R k , and if for arbitrary i, j ∈ {1, . . . , k} the n-dimensional measure of the intersection R i ∩R j is zero. The intersection R i ∩R j can either be empty, or a rectangle of dimension smaller than n.

A partition R 1 , . . . , R k is an ε-approximation of a function f : R ⊂ R n → R m , if we have f (x) -val(R i ) ∞ ≤ ε for all x ∈ R i and i = 1, . . . , k .
With each ε-approximation of f , we can associate a lower semi-continuous function f which is defined via f (x) = val(R i ) for all x in the interior of R i , and for any point x which lies in the boundary of at least one of the partition rectangles we define f

(x) = min{val(R i 1 ), . . . , val(R i )}, where x is contained in the rectangles R i 1 , . . . , R i . Finally, for a rectangle R = [a 1 , b 1 ]×. . .×[a n , b n ]
its midpoint is given by midpt(R) = ((a 1 + b 1 )/2, . . . , (a n + b n )/2), and its geometric realization, i.e., the set of points belonging to the rectangle, is denoted by |R|. In terms of the given function f , we assume that it is of a form suitable for evaluation over rectangles using the rigorous computer arithmetic. With these definitions, the algorithm for determining a piecewise constant approximation of a continuous function is given by Algorithm 1. In the algorithm, the rigorous arithmetic implementation of f is denoted by f . The properties of rigorous computer arithmetic discussed in the previous section readily imply that if Algorithm 1 terminates, then the constructed partition is an ε-approximation of f . Of course, there are instances where the algorithm does not terminate. To see this, suppose we are computing on a machine whose smallest strictly positive representable number is

ρ min > 0. Choose δ > 0 such that δ (ρ min + δ) < ρ min 2π , (3)
and consider the function f : R → R defined via

f (x) = sin 1 x + δ for all x ∈ R = [0, 1] . (4)
Any partition of [0, 1] into compact intervals whose endpoints are representable numbers has to contain an interval [0, a] with ρ min ≤ a ≤ 1. One can easily see that (3) is equivalent to

1 0 + δ - 1 a + δ ≥ 1 0 + δ - 1 ρ min + δ > 2π ,
and therefore the range of f over the interval [0, a] is the interval [-1, 1]. This immediately implies that any evaluation using rigorous computer arithmetic yields f ([0, a]) ⊃ [-1, 1], and the algorithm cannot terminate if we set ε < 1.

The above example illustrates that even for simple analytic real-valued functions on a compact one-dimensional domain the presented algorithm may not terminate. Therefore, in order to avoid infinite loops one has to set a maximal number of subdivisions that the algorithm is allowed to perform. If this number is exceeded, the algorithm terminate with the answer "cannot decide," i.e., the algorithm is a partial algorithm. For clarity of presentation, this simple modification is not part of the pseudo-code of Algorithm 1.

Nevertheless, on a purely theoretical level, the algorithm has the potential to terminate for every continuous function f : R → R m defined on a compact rectangle R ⊂ R n . This is a direct consequence of the resulting uniform continuity of f on R. In fact, if we assume more regularity of f , and that the computer implementation has infinite precision, then it is even possible to give an upper bound on the number of rectangles in an ε-approximation of f . This is the subject of the following result. Theorem 3.2 (Algorithm Complexity for Lipschitz Functions). Let R ⊂ R n denote a compact rectangle, and assume that the function f : R → R m is Lipschitz continuous with Lipschitz constant L. Suppose further that Algorithm 1 was successfully executed for an arithmetical expression of f , for the algorithm parameter ε > 0, and that the algorithm was run on a computer with infinite precision. Then the algorithm creates a subdivision of the rectangle R with

at most α(L)diam(R) 2ε n rectangles,
where the constant α(L) is taken from Theorem 3.1, and diam(R) denotes the diameter of R with respect to the maximum norm.

Proof. In infinite precision, the if-condition in Algorithm 1 is satisfied as soon as we have

max( f (R)) -min( f (R)) ≤ 2ε ,
where f denotes the interval implementation of f . In other words, the algorithm terminates on a rectangle A ⊂ R, as soon as the evaluation of f has diameter at most 2ε. Moreover, Theorem 3.1 implies that the interval representation has the Lipschitz constant α(L). Thus, for diam(A)

≤ 2ε α(L)
the rigorous evaluation of f on A has diameter at most 2ε. In the worst case, this means that the rectangle R has to be subdivided (α(L)diam(R))/(2ε) times in each coordinate direction, and this will lead to at most ((α(L)diam(R))/(2ε)) n rectangles in the final subdivision.

The above result was formulated for infinite precision in order to exclude examples such as (4). In actual computer implementations, the number of rectangles in the final subdivision of the rectangle is usually very close to the bound in the theorem, but it might be slightly larger. This indicates that the techniques which will be presented below can be expected to work extremely well in low dimensions. On average, however, we expect the complexity of the procedures developed in this paper to grow exponentially with the dimension of the ambient space, in accordance with the above theorem. Of course, the main factors in the complexity are the topological features of the considered function, and there are many functions for which this exponential growth is not of practical importance. Yet, as soon as a function is Lipschitz-continuous with optimal Lipschitz constant L ≥ 1, the computational complexity grows exponentially with the dimension. We would like to point out that this is yet another aspect of the so-called curse of dimensionality. When constructing cubical complexes of an n-dimensional box, we are in fact experiencing the same obstacles as the ones encountered in the construction of Delaunay triangulations or Voronoi diagrams of k points in R n , which requires a computation time of the order O(k (n+1)/2 ), see for example [START_REF] Chazelle | An optimal convex hull algorithm in any fixed dimension[END_REF].

The Subdivision Algorithm

The algorithm presented in Section 3.2 operates entirely on top-dimensional cells. This is sufficient if one's only interest is the computation of an ε-approximation of a continuous function. If, however, the goal of the computation is the (multi-dimensional) persistent homology of the underlying continuous function, then it is necessary to create a complex of some sort, together with an associated filtration. In our situation, this complex will be a rectangular CW-complex, and in addition to the top-dimensional cells one needs to efficiently store also the necessary lower-dimensional cells. The creation of this complex will be described in the current section.

In the following, we assume that the initial rectangle R used in Algorithm 1 is given as a rectangular CW-complex. More precisely, we assume that it is stored in a structure which contains all the boundary elements, and such that each rectangle in the structure can access both its boundary and coboundary rectangles. For this type of input, we describe a straightforward subdivision scheme which allows one to subdivide R and incorporate the resulting new rectangles, while keeping the structure of the rectangular CW-complex updated after each subdivision. In pseudo-code notation, this procedure can be implemented as in Algorithm 2. This algorithm subdivides a rectangle R in the direction of the i-th coordinate into two sub-rectangles R 1 and R 2 which have the same dimension as R, and which share a common face R 3 . The algorithm uses references in R to the boundary and coboundary elements of R to maintain the global structure of the rectangular CW-complex. See also Figure 6 for an illustration. While the above algorithm only subdivides the rectangle R along one coordinate direction, in some cases one might like to perform subdivisions in all directions at once. In this case, one can use Algorithm 3.

To close this brief section, we would like to point out that the use of rectangular CW-complexes as described above for the computation of persistent homology is only one possible approach. Alternatively, one could construct a simplicial complex directly from the collection of top-dimensional cells obtained by the approximation algorithm from the last section. Since all rectangles considered in our approach are convex sets, one can construct a weighted nerve complex which inherits the (multidimensional) persistent homology from the original rectangular CW-complex associated with the lower semi-continuous approximation of the underlying continuous function. This is a consequence of the celebrated nerve theorem, see [START_REF] Edelsbrunner | Computational Topology[END_REF] for more details. The nerve complex can be constructed as follows. Its vertices are given by the top-dimensional rectangles in the ε-approximation. Moreover,

R A 1 B 2 B 3 B 1 R 3 R 2 R 1 A 3 2 A
Figure 6: Illustration of Algorithm 2. Given a two-dimensional rectangle R as shown on the left, as well as a coordinate direction i, the rectangle is first divided into two two-dimensional rectangles R 1 and R 2 by subdividing the edges of R parallel to the i-th coordinate direction. This also creates their common face R 3 . Subsequently, the boundary elements of R which are parallel to the i-th coordinate direction and which contain the subdivision coordinate in their interior are processed by recursive calls of the algorithm. This creates the new one-dimensional rectangles A 1 , A 2 , B 1 , and B 2 , and the new zero-dimensional rectangles A 3 and B 3 . Once all of these have been created, the boundaries and coboundaries are updated accordingly.

as transition elements between one function value and another. In addition, our approach makes it easier to determine persistent homology in the case of periodic boundary conditions. For this, one only has to impose the periodic boundary conditions on the initial CW-complex representation of the domain -and the presented subdivision algorithms automatically lead to the correct persistent homology, as long as this initial CW-complex configuration contains at least two top-dimensional rectangles in each coordinate direction.

Rigorous Persistence Approximation

In this section, we combine the ε-approximation algorithm presented in Section 3.2 with the rectangle subdivision scheme presented in Section 3.3. This allows us to construct a filtered regular CW-complex which rigorously approximates a given continuous function f on a rectangular domain R, and at the same time provides a rectangular CW-complex representation of R compatible with the filtering function. The construction will be performed in such a way that the (multidimensional) persistence of the final rectangular CW-complex can easily be computed, and that the corresponding persistence diagram deviates from the one for the underlying function f by less than ε with respect to the bottleneck distance. The detailed procedure is listed in Algorithm 4, where the initial rectangle R is subdivided according to either of the algorithms contained in Section 3.3 until the rigorous enclosure f (Q) of the range of f over a rectangle Q deviates from the midpoint value by at most ε. When this criterion is met for all rectangles Q in the subdivision of R, the final CW-complex decomposition is returned together with the cell values set as the function values of f at the cell midpoints. See Algorithm 4 for more details.

It is clear from our construction, that at the end of the while-loop in Algorithm 4 one has constructed a lower semi-continuous piecewise constant function f of the given continuous function f as described earlier in Section 3.2. We would like to recall that this approximation satisfies the estimate |f (x) -f (x)| ≤ ε for all x ∈ R , mination, the algorithm returns a rigorous upper bound on the approximation error which is exactly the largest of the range enclosure radii. The precise implementation of the algorithm can be found in Algorithm 5, and it furnishes a piecewise constant approximation f for the given continuous function f defined on the rectangle R, together with the rigorous error bound Error of Approximation. This algorithm can trivially be extended to the case of vector-valued functions f : R → R m .

Numerical Case Studies

In this final section of the paper we present a few case studies which demonstrate how the proposed algorithms perform in practice. The numerical experiments include a one-dimensional continuous, but not everywhere differentiable function, a two-dimensional sample function which serves as a standard benchmark in nonlinear optimization, as well as an example from mathematical materials sciences. An implementation of the algorithms used to perform these computations will be made available upon acceptance of this paper. It can be found at https://github.com/pdlotko/FunTop, and is available under the GPLv3 license.

Two Benchmarks from Optimization

As a first example of the algorithms in this paper we consider the function

f (x) = |sin(6πx)| 1 + x 2 + 3 cos(2πx) 10 for x ∈ R = [0, 1] .
This function is clearly continuous, but not everywhere differentiable. While it is not possible to use the methods of [START_REF] Jaquette | Rigorous computation of persistent homology, accepted to Mathematics of Computation[END_REF] to approximate the persistence of f , our Algorithm 4 can be applied. The results of this method for ε = 0.18, ε = 0.06, and ε = 0.02 were already presented in Figure 4. In these persistence diagrams, the dashed red line indicates all points above the diagonal whose maximum norm distance to the diagonal is exactly equal to ε. In view of Remark 3.4, this means that only in the upper right persistence diagram, which corresponds to ε = 0.18, one of the generators cannot be established rigorously. However, for the two smaller ε-values all generators have been proven to exist. They indicate that the function has seven local minima, that these minima start appearing in a stratified way with respect to increasing function values, first one, and then in three pairs of two, and that between any two consecutive minima a local maximum annihilates homology components. In this sense, our algorithms can be easily used to rigorously solve one-dimensional optimization problems. But the methods are by no means restricted to one-dimensional spaces. Consider for example the two-dimensional function

f (x, y) = g(45x -15, 45y -15) , (6)
where the function g is the standard two-dimensional Ackley function given by g(x, y) = a + e -ae -b √ (x 2 +y 2 )/2 -e (cos(cx)+cos(cy))/2 .

For our simulations, we choose the parameters a = 20, b = 1/5, and c = π/2, and we consider the function f on the rectangle R = [0, 1] 2 ⊂ R 2 . A partially constant approximation of this function has already been shown in Figure 1, and the actual function f is shown in Figure 7. The figure shows the graph of the function f defined in [START_REF] Cerri | Hausdorff stability of persistence spaces[END_REF], which is based on the standard twodimensional Ackley function g given in [START_REF] Cerri | Betti numbers in multidimensional persistent homology are stable functions[END_REF]. The function g is one of the basic benchmark functions in nonlinear optimization due to its unique global minimum at (0, 0), combined with a large number of local minima. Recall that a rigorous piecewise constant approximation f for ε = 2 is shown in the right panel of Figure 1 7, for a variety of values of the parameter ε.

For the above function f we determined rigorous ε-approximations of the associated persistence diagram, for a variety of ε-values. The results from these computations can be found in Table 1.

i.e., the radius of the interval is exactly half its width. However, when using finite arithmetic on a computer, either or both of these numbers might not be a representable number, which leads to interval answers. In this case, the definition in [START_REF] Carlsson | The theory of multidimensional persistence[END_REF] ensures that rad(I) contains a rigorous and tight upper bound on the radius of the interval. These computational times seem to exhibit quadratic growth in 1/ε, as would be expected from our Theorem 3.2. In fact, simple quadratic regression implies a near perfect fit of the computational data with the polynomial T (s) = 71.94s 2 -25.36s + 1.11, and with correlation coefficient 0.9869.

The Diblock-Copolymer Model

Topological methods have seen a wide variety of applications in the physical sciences, see for example the articles in the special issue [START_REF] Day | Topology in dynamics, differential equations, and data[END_REF]. We therefore close this paper with a related application of our methods to the rigorous approximation of persistence in an example from materials science. We consider the so-called diblock copolymer model which is the fourth-order parabolic partial differential equation given by which is an evolution equation for the unknown function u : R + 0 × Ω → R, where Ω ⊂ R d denotes an arbitrary domain with sufficiently smooth boundary. In this formulation, the function F is the negative derivative of a double-well potential, and one usually considers F (u) = u -u 3 . Moreover, the directional derivative ∂u/∂ν denotes the derivative of u in the direction of the outward unit normal vector at a point on the boundary ∂Ω.

u t = -∆ ε 2 ∆u + F (u) -σ(u -µ) in Ω , (8)
While a more complete introduction to the model can be found in [START_REF] Wanner | Topological analysis of the diblock copolymer equation[END_REF], we only focus on a few important aspects. The partial differential equation ( 8) models phase separation phenomena in polymer materials which are composed of two chemically incompatible monomers, say of type A and type B, which are arranged in long polymer chains. The value of the phase variable u describes the local material composition in the following way. Function values of the solution u(t, x) which are close to +1 are interpreted as only monomer A being present near point x ∈ Ω and at time t ≥ 0, and the value -1 indicates that only monomer B is present. Finally, values in between correspond to mixtures of the two components, with zero representing an equal mixture. The parameter µ denotes the average mass of the mixture using the same convention, and the two remaining parameters ε > 0 and σ ≥ 0 are dimensionless interaction lengths. Informally, ε > 0 being small corresponds to short range repulsions being strong, inducing a strong compulsion to separate, while σ being large represents strong long range chain elasticity forces, inducing a strong compulsion to hold together. Notice that for σ = 0 the diblock copolymer model [START_REF] Chazal | Proximity of persistence modules and their diagrams[END_REF] reduces to the celebrated Cahn-Hilliard equation, which serves as a basic model for the phase separation phenomena.

If one considers solutions u which start at a small random perturbation u 0 of the constant initial state µ, and if one assumes that the average value of u 0 equals µ, then solutions quickly start to grow in amplitude, until the function values reach values close to ±1. At this point, the phases have separated and one can observe a wide variety of patterns, see for example Figure 8. One can clearly see from these images that the topology of the phase separated state depends on the total mass µ. In fact, even though different initial states lead to different microstructures, the observed patterns for fixed µ share the same characteristic features. In other words, the parameter µ clearly affects the observed topology of the pattern.

In our recent paper [START_REF] Lotko | Topological microstructure analysis using persistence landscapes[END_REF] we studied a stochastic version of the Cahn-Hilliard equation to show that in fact the converse is true as well, i.e., that the topology of the observed patterns can be used to determine the underlying total mass µ to surprising accuracy. This was accomplished in the following way. For every mass value, we performed a number of Monte-Carlo type simulations, to obtain a collection of persistence diagrams for the resulting microstructures. By equivalently reformulating the diagrams as persistence landscapes (see [START_REF] Lotko | Topological microstructure analysis using persistence landscapes[END_REF][START_REF] Bubenik | A persistence landscapes toolbox for topological statistics[END_REF]), it is then possible to determine an averaged persistence landscape, which can be thought of as the "typical" pattern topology for a given µ-value. If one then considers an arbitrary microstructure, one can determine the associated mass value by finding the closest averaged landscape. This classification scheme exhibits surprising accuracy in determining the parameter µ.

Rather than performing an analogous study of the diblock copolymer equation for σ = 0 in full detail, we use the methods developed in this paper to determine the distance matrices between the averaged persistence landscapes for each considered mass value. We consider the same parameter values as in Figure 8, i.e., we let ε = 0.025 and σ = 10, and consider the two-dimensional square domain Ω = (0, 1) 2 . Finally, we simulated the equation at the 51 mass values µ = k • 0.01, with multiplier k = 0, . . . , 50. In each case, one hundred Monte-Carlo simulations were performed to determine the averaged persistence landscape as in [START_REF] Bubenik | A persistence landscapes toolbox for topological statistics[END_REF], from the associated 100 persistence diagrams of phase separated patterns. All of the necessary persistence diagrams were computed using Algorithm 4 with the error bound set to 0.2. The simulations of the partial differential equation employed a semi-linear spectral method with 128 2 eigenmodes for the solution representation, and 256 2 modes for the nonlinearity computations. Note that, in contrast to [START_REF] Lotko | Topological microstructure analysis using persistence landscapes[END_REF], in the present paper the computations of persistence diagrams of the Fourier series approximation of the solution of the partial differential equation are rigorous computations.

Figure 1 :

 1 Figure 1: Sample approximating functions f . The left image is for the one-dimensional function f (x) = 2 -25x + 108x 2 -162x 3 + 81x 4 on the interval D = [0, 1], and with approximation parameter ε = 1/2. The function f is shown in blue, while its approximation f is depicted in red. Red circles indicate the function values at the points in D which are contained in more than one intervals of the final decomposition. The right panel shows an approximation for the Ackley function with ε = 2.

  [a, b] has the two faces [a] := [a, a] and [b] := [b, b]. Intervals as defined above are the building blocks of higher-dimensional sets in the following sense

  (1,),(1,),(1,),(1,)} d1: {} d0: {(1,2),(1,2),(1,),(1,)} d1: {} d0: {(1,2),(1,2),(1,3),(1,)} d1: {(3,)} d0: {(1,2),(1,2),(1,3),(1,∞)} d1: {(3,4)} d0: {} d1: {}

Figure

  Figure3: From left to right, the images show a filtration on a simple 3 by 3 cubical complex (which is a special case of a rectangular CW complex). At each stage in the filtration, the black cubes are part of the complex. The numbers from 0 to 4 below the images denote the value of the filtering function. Below the filtering function the corresponding persistence intervals in dimension 0 (d0) and dimension 1 (d1) are presented. At level 0, both the complex and the collection of persistence intervals are empty. For complex 1, the four added cubes constitute four connected components, indicated by as of yet unbounded persistence intervals in dimension 0. Upon moving to complex 2 in the filtration, two pairs of the connected components merge, and therefore two persistence intervals which were created at level 1 now terminate. With complex 3, only one connected component remains, and therefore another one of the two unbounded intervals in dimension 0 terminates at time 3. At the same time, a homology class in dimension 1 is created, as indicated by the new and as of yet unbounded persistence interval in dimension 1. When we reach the final complex 4, the homology class in dimension 1 becomes trivial, which is indicated by the termination of the unique persistence interval in dimension 1. Since the filtration does not evolve any further, the only remaining unbounded persistence interval in dimension 0 stays infinite.

Figure 4 :

 4 Figure 4: Sample persistence diagrams in dimension zero for piecewise constant approximations of the non-differentiable function f (x) = | sin(6πx)|/(1+x 2 )+ 3 cos(2πx)/10. The function f , together with its approximation f for ε = 0.18, is shown in the upper left image, while the upper right image contains the persistence diagram for f . The two lower panels show the persistence diagrams for f if we use instead ε = 0.06 (left) or ε = 0.02 (right).

Figure 5 :

 5 Figure 5: Sample persistence diagram matching. The left two images show persistence diagrams with seven and four birth/death pairs, respectively. The third diagram shows a sample transport plan η between the two persistence diagrams.

Figure 7 :

 7 Figure 7: Visualization of the scaled Ackley function in two dimensions on the domain [0, 1] × [0, 1].The figure shows the graph of the function f defined in[START_REF] Cerri | Hausdorff stability of persistence spaces[END_REF], which is based on the standard twodimensional Ackley function g given in[START_REF] Cerri | Betti numbers in multidimensional persistent homology are stable functions[END_REF]. The function g is one of the basic benchmark functions in nonlinear optimization due to its unique global minimum at (0, 0), combined with a large number of local minima. Recall that a rigorous piecewise constant approximation f for ε = 2 is shown in the right panel of Figure1.

Figure 8 :

 8 Figure 8: Sample patterns produced by the diblock copolymer model (8) on the two-dimensional square domain Ω = (0, 1) 2 . The images are all for the parameter values ε = 0.025 and σ = 10, and for end times which correspond to the formation of distinct phases. From top left to bottom right the images are for mass values µ = 0, 0.1, . . . , 0.5.

Table 1 :

 1 . Computation times for rigorously determining the persistent homology of the function f shown in Figure

	ε	1	0.75	0.5	0.4	0.3	0.25
	time (s) 0.003487 0.008145 0.018178 0.027688 0.029342 0.132159
	ε	0.2	0.1	0.075	0.05	0.025	0.02
	time (s) 0.203003 2.17143	6.62386	38.8436	632.568	1704.78

Notice that in infinite precision, the two members of the set on the right-hand side of (5) have to be the same,

Acknowledgments

The first author of the paper would like to thank Michael Kerber for suggesting the algorithm presented in the Section 3.5, and Steve Oudot for all the helpful comments. The first author was supported by the Advanced Grant of the European Research Council GUDHI (Geometric Understanding in Higher Dimensions). The second author was partially supported by NSF grants DMS-1114923 and DMS-1407087.

Algorithm 2 Rectangle subdivision in one coordinate direction.

Input: R = [r 1,0 , r 1,1 ] × . . . × [r n,0 , r n,1 ], i ∈ {1, . . . , n}, x ∈ (r i,0 , r i,1 ); Output: Rectangles R 1 , R 2 , R 3 which are added to the complex such that |R 1 | ∪ |R 2 | = |R| and R 3 is the intersection of R 1 and R 2 ; List of rectangles L; for every S in bdr R such that x is in the interior of the projection of S onto the i-th coordinate direction. do Put S into L; Create a rectangle R 1 = [r 1,0 , r 1,1 ] × . . . × [r i,0 , x] × . . . × [r n,0 , r n,1 ]; Create a rectangle R 2 = [r 1,0 , r 1,1 ] × . . . × [x, r i,1 ] × . . . × [r n,0 , r n,1 ]; Create a rectangle R 3 = [r 1,0 , r 1,1 ] × . . . × [x, x] × . . . × [r n,0 , r n,1 ]; Put R 3 into the boundary of R 1 and R 2 ; for every A in the boundary of R such that A ∈ L do Remove R from the coboundary of A and replace it by either R 1 or R 2 ; Put A into the boundary of R 1 or R 2 wherever it belongs; for every B in the coboundary of R do

Remove R from boundary of B and replace it by both R 1 and R 2 ; Put B into the coboundary of both R 1 and R 2 ; Remove R from the complex; for every S ∈ L do (S 1 , S 2 , S 3 ) = SubdivideRectangle(S, i, x); Put S into the boundary of R , and R into the coboundary of S for ∈ {1, 2, 3};

Algorithm 3 Rectangle subdivision in all coordinate directions.

) for i ∈ {1, . . . , n}; Output: Collection of rectangles obtained after subdividing R in all coordinate directions; List of rectangles L; Add R to L; for dim = 1 to n do List of rectangles L ; for Every rectangle R in the list L do Call Algorithm 2 for R, dim and x dim , which returns the top-dimensional rectangles R 1 and R 2 ; Add R 1 , R 2 to L ; L = L ; Return all the rectangles obtained in the last iteration of the loop; a set of k + 1 of these vertices forms a k-dimensional simplex, if the associated top-dimensional cells have nonempty intersection. Clearly, in the situation of this paper which is based on the subdivision of rectangles, for a planar set only 2 2 = 4 rectangles can have nonempty intersections. For rectangular sets embedded in R 3 one can have at most 2 3 = 8 intersections, and in the general case of rectangles in R d the number of intersections is bounded by 2 d , i.e., it is of exponential complexity with respect to the embedding dimension d. We would like to point out, however, that it suffices to include only simplices up to the embedding dimension d, which leads to a significant reduction in the size of the nerve complex without changing its associated persistent homology. From a computational perspective, it is not immediately clear whether the nerve complex leads to a more memory efficient representation than the rectangular CW-complex constructed via the subdivision algorithm. Nevertheless, we have decided to use the rectangular CW-complex approach, as it reflects the function approximation. In a sense, the lower-dimensional rectangles can be thought of Algorithm 4 Rigorous approximation of the persistence of a continuous function.

Subdivide Q in all directions. Put all resulting rectangles in Q; Determine the lower star filtration of the final complex and compute its persistent homology; return Persistence of the resulting rectangular CW-complex; which is guaranteed through the use of rigorous interval arithmetic. From a computational perspective, the algorithm presented above returns a rectangular CW-complex X in which every topdimensional rectangle Q has been assigned a value val(Q). For the computation of persistence, however, we need a fast way to determine the filtration on X which is induced by the approximating function f . This can be accomplished as follows. According to our construction, the value of a top-dimensional rectangle Q is given by val(Q). If, on the other hand, we consider a lowerdimensional rectangle S in X , and if S 1 , . . . , S k denote all top-dimensional rectangles in X which contain S in their boundary, then we define the value of S via val(S) = min{val(S 1 ), . . . , val(S k )}. These definitions give rise to a filtration of all rectangles in X , and it is not difficult to verify that for every threshold α ∈ R the collection of all rectangles Q with val(Q) ≤ α is a CW-subcomplex of X , whose geometric representation consists of all x ∈ R for which f (x) ≤ α. In other words, this filtration of rectangles provides a fast way to determine the filtration on X induced by f . The above construction is referred to as the lower star filtration on X induced by the values of the top-dimensional rectangles.

Using the lower star filtration on all rectangles one can now use standard persistence software to determine the associated persistence diagram. In the case of one-dimensional persistence, it follows from Theorem 2.2 that this persistence diagram has bottleneck distance at most ε from the persistence diagram of f over the rectangle R. In addition, we would like to point out that Algorithm 4 also works -without significant changes -for a vector-valued function f : R → R m . In this case one obtains an approximation f which can be used to approximate multi-dimensional persistence. In view of Theorem 2.3 and similar stability theorems, one retains the same error bound ε for the approximation results.

Remark 3.3. For the above presentation, we have always assumed that the underlying function f is at least continuous. The main reason for this is our requirement that f has to be amenable to rigorous arithmetic evaluations. Most available software packages for rigorous arithmetic, such as the interval arithmetic implementations used throughout this paper, consider only standard arithmetic operations and standard functions, all of which are continuous on their domains of definition. Nevertheless, in principle the above approach can be extended to general tame functions, which ensure that the associated persistence diagrams are finite. See [START_REF] Chazal | Proximity of persistence modules and their diagrams[END_REF][START_REF] Chazal | The Structure and Stability of Persistence Modules[END_REF] for more details. Note, however, that in general it is not easily possible to verify whether a given function is tame or not. Moreover, it is possible that Algorithm 4 terminates for a non-tame function, provided there are only finitely many persistence Algorithm 5 Greedy approximation of the persistence of a function. intervals which are longer than ε. Needless to say, this latter fact is a necessary condition for the termination of Algorithm 4.

Remark 3.4. In the collection of persistence intervals produced by Algorithm 4, there generally are intervals of length at most ε. Intervals of this type do not have to correspond to nontrivial intervals in the persistence diagram of the given function f . For this reason, such intervals will typically be ignored.

Greedy Approximation

It was already mentioned earlier that the algorithm presented in Section 3.4 does not necessarily terminate, even if the underlying function f is Lipschitz continuous. On the one hand, this can be a consequence of the finite precision inherent in rigorous computer arithmetic. When a certain subdivision depth is reached, the range enclosures of f over a rectangle may not become smaller even if further subdivisions are performed. On the other hand, a given continuous function f does not necessarily have to be tame, i.e., it will not give rise to a finite persistence diagram.

To accommodate such situations, we close this section with an algorithm which removes the input parameter ε, inspired by the method in [START_REF] Bendich | Computing robustness and persistence for images[END_REF]. This leads to a procedure which iteratively subdivides the rectangles in the subdivision of R in a greedy way until the process is stopped by the userwhich could be accomplished either by an interactive interrupt, or by specifying a fixed number of subdivisions ahead of time. The overall goal of the algorithm is to increase the precision of the piecewise constant function approximation f at every subdivision. For this, the algorithm always subdivides the subrectangle Q of R for which the radius of the range enclosure f (Q) is maximal, i.e., the regions over which the range of f has the largest variation are subdivided first. In this context, the radius of an interval I ⊂ R is defined by

and this definition does provide an upper bound on the actual radius in computer arithmetic implementations. 1 The whole procedure can easily be implemented using a priority queue, and upon ter- The results from these computations can be found in Figure 9. This figure contains four heat plots for the distance matrices of the averaged persistence landscapes. In each of the plots, the horizontal and vertical axes correspond to the scaled mass parameter 100µ, i.e., they reflect total mass values µ between µ = 0 and µ = 0.5 in steps of 0.01. The top row contains distance matrices for the 0-dimensional persistence landscapes, while the bottom row is for dimension 1. Moreover, the left column uses the L 1 -norm to measure the distance between the landscapes, and the right column employs the L 2 -norm. These heat maps indicate that depending on the dimension of the persistence landscape, different regions of µ-values can be differentiated. Since we are using nonnegative mass values, for large µ one expects a larger number of persistence intervals, which correspond to the isolated (negative) droplets in the (positive) background matrix. Thus, the heat maps in the first row show greater variability for µ ≈ 0.5. In contrast, in the droplet regime, the sublevel sets of the function u are very unlikely to contain any homology generators in dimension one, and this is reflected in the large region of distance close to zero in the two bottom panels of Figure 9. Nevertheless, for mass values µ close to zero nontrivial homology can be observed in dimension one, and this leads to the greater variability of the heat maps in the lower left corner. These observations are very much in line with the ones in [START_REF] Lotko | Topological microstructure analysis using persistence landscapes[END_REF]. However, we refrain from a more detailed analysis, as it is not the main topic of the present paper. Nevertheless, these preliminary results can be considered as a first step towards a quantitative topological understanding of the topology of microstructures generated by the diblock copolymer equation.