Computing the Multicover Bifiltration

Given a finite set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\subset {\mathbb {R}}^d$$\end{document}A⊂Rd, let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Cov}_{r,k}$$\end{document}Covr,k denote the set of all points within distance r to at least k points of A. Allowing r and k to vary, we obtain a 2-parameter family of spaces that grow larger when r increases or k decreases, called the multicover bifiltration. Motivated by the problem of computing the homology of this bifiltration, we introduce two closely related combinatorial bifiltrations, one polyhedral and the other simplicial, which are both topologically equivalent to the multicover bifiltration and far smaller than a Čech-based model considered in prior work of Sheehy. Our polyhedral construction is a bifiltration of the rhomboid tiling of Edelsbrunner and Osang, and can be efficiently computed using a variant of an algorithm given by these authors. Using an implementation for dimension 2 and 3, we provide experimental results. Our simplicial construction is useful for understanding the polyhedral construction and proving its correctness.


Introduction
Let A be a finite subset of R d , whose points we call sites.For r ∈ [0, ∞) and an integer k ∈ N = {0, 1, 2, . ..}, we define Cov r,k := b ∈ R d | ||b − a|| ≤ r for at least k sites a ∈ A .
Thus, Cov r,k is the union of all k-wise intersections of closed balls of radius r centered at the sites; see Figure 1.Define a bifiltration to be a collection of sets such that C r,k ⊆ C r ,k whenever r ≤ r and k ≥ k .Clearly, the sets Cov := (Cov r,k ) (r,k)∈[0,∞)×N form a bifiltration.This is known as the multicover bifiltration.It arises naturally in topological data analysis (TDA), and specifically, in the topological analysis of data with outliers or non-uniform density [19,28,46].Figure 1: The 2-and 3-fold cover of a few points with respect to a certain radius.The first homology of the 2-fold cover is trivial, while the first homology of the 3-fold cover is non-trivial.
We wish to study the topological structure of the bifiltration Cov algorithmically in practical applications, via 2-parameter persistent homology [11].For this, the natural first step is to compute a combinatorial model of Cov, that is, a purely combinatorial bifiltration C which is topologically equivalent to Cov.This step is the focus of the present paper.For computational efficiency, C should not be too large.
In fact, we propose two closely related combinatorial models C, one polyhedral and one simplicial.The polyhedral model is a bifiltration of the rhomboid tiling, a polyhedral cell complex in R d+1 recently introduced by Edelsbrunner and Osang to study the multicover bifiltration [28].Edelsbrunner and Osang have given an efficient algorithm for computing the rhomboid tiling [29], and this adapts readily to compute our bifiltration.We use the simplicial model to prove that the polyhedral model is topologically equivalent to Cov.

Motivation and prior work
For k = 1 fixed, (Cov r,1 ) r∈[0,∞) is the well-known offset filtration (also known as the union of balls filtration), a standard construction for analyzing the topology of a finite point sample across scales [27].It is a central object in the field of persistent homology.While the persistent homology of this filtration is stable to small geometric perturbations of the sites [21], it is not robust with respect to outliers, and it can be insensitive to topological structure in high density regions of the data.
Within the framework of 1-parameter persistent homology, there have been many proposals for alternative constructions which address these issues.These approaches include the removal of low density outliers [10], filtering by a density function [7,16,17], distance to measure constructions [1,9,19,31], kernel density functions [44], and subsampling [4].A detailed overview of these approaches can be found in [6].
Several of these constructions have good stability properties or good asymptotic behavior.However, as explained in [6], all of the known 1-parameter persistence strategies for handling outliers or variations in density share certain disadvantages: First, they all depend on a choice of a parameter.Typically, this parameter specifies a fixed spatial scale or a density threshold at which the construction is carried out.In the absence of a priori knowledge about the structure of the data, it may be unclear how to select such a parameter.And if the data exhibits topological features at multiple spatial or density scales, it may be that no single parameter choice allows us to capture all the structure present in the data.Second, constructions that fix a scale parameter are unable to distinguish between small spatial features and large ones, and constructions that fix a density or measure parameter are unable to distinguish features in densely sampled regions of the data from features involving sparse regions.
A natural way to circumvent these limitations is to consider a 2-parameter approach, where one constructs a bifiltration from the data, rather than a 1parameter filtration [11].The multicover bifiltration is one natural option for this.Alternatives include the density bifiltrations of Carlsson and Zomorodian [11], and the degree bifiltrations of Lesnick and Wright [41]; again, we refer the reader to [6] for a more detailed discussion.Among these three options, the multicover bifiltration has two attractive features which together distinguish it from the others.First, its construction does not depend on any additional parameters.Second, the multicover bifiltration satisfies a strong stability property, which in particular guarantees robustness to outliers [6].
There is a substantial and growing literature on the use of bifiltrations in data analysis.Most approaches begin by applying homology with coefficients to the bifiltration, to obtain an algebraic object called bipersistence module.In contrast to the 1-parameter case, where the algebraic structure of a persistence module is completely described by a barcode [51], it is well-known that defining the barcode of bipersistence modules is problematic [11].Nevertheless, one can compute invariants of a bipersistence module which serve as useful surrogates for a barcode, and a number of ideas for this have been proposed [11,14,32,41,49].
Regardless of which invariants of the multicover bifiltration we wish to consider, to work with this bifiltration computationally, the natural first step is to find a reasonably sized combinatorial (i.e., simplicial or polyhedral) model for the bifiltration.With such a model, recently developed algorithms such as those described in [42] and [36] can efficiently compute minimal presentations and standard invariants of the homology modules of the bifiltration.
In the 1-parameter setting, there are two well-known simplicial models of the offset filtration.The Čech filtration, is given at each scale by the nerve of the balls; the equivalence of the offset and Čech filtrations follows from the Persistent Nerve Theorem.For large point sets, the full Čech filtration is too large to be used in practical computations.However, the alpha filtration (also known as the Delaunay filtration) [25,27] is a much smaller subfiltration of the Čech filtration which is also simplicial model for the offset filtration.It is given at each scale by intersecting each ball with the Voronoi cell [50] of its center, and then taking the nerve of the resulting regions.For d small (say d ≤ 3), the Delaunay filtration is readily computed in practice for many thousands of points.
It is implicit in the work of Sheehy [46] that the multicover bifiltration has an elegant simplicial model, the subdivision-Čech bifiltration, obtained via a natural filtration on the barycentric subdivision of each Čech complex; see also [12, Appendix B] and [6,Section 4].However, the subdivision-Čech bifiltration has exponentially many vertices in the size of the data, making it even more unsuitable for computations than the ordinary Čech filtration.
Edelsbrunner and Osang [28] therefore seek to develop the computational theory of the multicover bifiltration using higher-order Delaunay complexes [26,38], taking the alpha filtration as inspiration.Assuming the sites are in general position, they define a polyhedral cell complex in R d+1 called the rhomboid tiling, which contains all higher-order Delaunay complexes as planar sections.Using the rhomboid tiling, they present a polynomial time algorithm to compute the barcodes of a horizontal or vertical slice of the multicover bifiltration (i.e., of a one-parameter filtration obtained by fixing either one of the two parameters r, k).The case of fixed r and varying k is more challenging because the order-k Delaunay complexes do not form a filtration.The authors construct a zigzag filtration for this case.The problem of efficiently computing 2-parameter persistent homology of the multicover bifiltration is not addressed by [28].
We note that the subdivision-Čech bifiltration is more general than the Rhomboid tiling: the rhomboid tiling is defined only for Euclidean data, whereas the topological equivalence of the subdivision-Čech and multicover bifiltrations extends to data in any metric space where finite intersections of balls are contractible.

Contributions
We introduce the first efficiently computable combinatorial models of the multicover bifiltration Cov.First, we introduce a simplicial model, whose construction is based on two main ideas: In order to connect the higher-order alpha complex constructions for (r, k) and (r, k + 1), we simply overlay their underlying covers to a "double-cover", whose nerve is a simplicial complex that contains both alpha complexes.This yields a zigzag of simplicial filtrations.The second main idea is that this zigzag can be "straightened out" to a (non-zigzaging) bifiltration, simply by taking unions of prefixes in the zigzag sequence.This straightening technique has previously been used by Sheehy to construct sparse approximations of Vietoris-Rips complexes [47].Together, these two ideas give rise to a bifiltration S-Del of simplicial complexes.
The bifiltration S-Del can also be obtained directly as the persistent nerve of a "thickening" of Cov constructed via mapping telescopes.This observation leads to a simple proof of topological equivalence (i.e., weak equivalence; see Section 2) of Cov and S-Del via the Nerve Theorem.It follows that the persistent homology modules of Cov and S-Del are isomorphic.
Our second contribution is to show that the rhomboid tiling as defined in [28] also gives rise to a (non-zigzaging) bifiltration of polyhedral complexes that is topologically equivalent to the multicover.We proceed in two steps: First, we slice every rhomboid at each integer value k (slightly increasing the number of cells) and adapt the straightening trick used to construct S-Del.We prove the topological equivalence of this construction with the multicover bifiltration by relating the slice rhomboid filtration with S-Del.The main observation is a one-to-one correspondence of maximal-dimensional cells in both constructions, which leads to a proof via the Nerve Theorem.Second, we relate the sliced and unsliced rhomboid tilings at every scale via a deformation retraction.
We give size bounds for both of the bifiltrations we introduce.For n points in R d , we show that their size is O(n d+1 ).This is a decisive improvement over Sheehy's Čech-based construction, which has exponential dependence on n.
An efficient algorithm for computing rhomboid tilings has recently been presented in [29]; hence, using the accompanying implementation rhomboidtiling of this algorithm, our second contribution gives us an efficient software to compute a bifiltration of cell complexes equivalent to Cov, currently for points in R 2 and R 3 .We combine this implementation with the libraries mpfree and rivet to demonstrate that minimal presentations of multicover persistent homology modules can now be efficiently computed, often within seconds, as can invariants such as the Hilbert function.

Filtrations
For P a poset, define a (P -indexed) filtration to be a collection of topological spaces X = (X p ) p∈P indexed by P , such that X p ⊆ X q whenever p ≤ q ∈ P .For example, an N-indexed filtration X is a diagram of spaces and inclusions of the following form: A morphism ϕ : X → Y of P -indexed filtrations is a collection of continuous maps (ϕ p : X p → Y p ) p∈P which commute with the inclusions in X and Y .In the language of category theory, a P -indexed filtration is a functor P → Top whose internal maps are inclusions, and a morphism is a natural transformation.Recall that the product poset P × Q of posets P and Q is defined by taking (p, q) ≤ (p , q ) if and only if p ≤ p and q ≤ q .When P is the product of two totally ordered sets, we call a P -indexed filtration a bifiltration.
In the classical homotopy theory of diagrams of spaces, there is a standard analogue of the notion of homotopy equivalence for diagrams of spaces, called weak equivalence.We now define a version of this for P -indexed filtrations: A morphism of filtrations ϕ : X → Y is called an objectwise homotopy equivalence if each ϕ p : X p → Y p is a homotopy equivalence.If there exists a finite zigzag diagram of objectwise homotopy equivalences connecting X and Y , then we say that X and Y are weakly equivalent.The terminology originates from the theory of model categories [23,34].See [5,39,45] for discussions of weak equivalence of diagrams in the context of TDA.
Remark 1.To motivate the consideration of zigzags in the definition above, we note that for X and Y a pair of weakly equivalent P -indexed filtrations, there is not necessarily an objectwise homotopy equivalence f : X → Y .For example, let P = R, X be the offset filtration on {0, 1} ⊂ R, and Y be the nerve filtration of X.It is easy to check that there is no objectwise homotopy equivalence f : X → Y .On the other hand, it follows from the Persistent Nerve Theorem (Theorem 3 below) that X and Y are weakly equivalent.Moreover, one can construct a similar example of weakly equivalent filtrations for which there is no objectwise homotopy equivalence in either direction. (1,1) Figure 2: Left: A bifiltration of good covers over {1 < 2 < 3} × {3 < 2 < 1}.Right: A bifiltration consisting of the nerves of the covers.The Persistent Nerve Theorem ensures that not only the individual spaces at scales (n, m) are homotopy equivalent, but also that the two bifiltrations are weakly equivalent.
An objectwise homotopy equivalence ϕ : X → Y induces isomorphisms on the persistent homology modules of X and Y .Hence, weakly equivalent filtrations have isomorphic persistent homology modules.
We say a P -indexed filtration X is Euclidean if X p ⊂ R n for some n and all p ∈ P .

The Persistent nerve theorem
We say that the cover is good if it is finite and consists of closed, convex sets [27].
One version of the Nerve Theorem asserts that X and Nrv (X) are homotopy equivalent whenever X is a good cover of X [27,40].It is TDA folklore that this version of the Nerve Theorem can be extended to a persistent version; a proof appears in [3]; see also [18] for formulation of the Persistent Nerve Theorem in terms of open covers.
In order to state the Persistent Nerve Theorem for closed, convex covers, we first extend the definition of a cover.
Definition 2 (Cover of a filtration).Let P be a poset and X a P -indexed Euclidean filtration.A cover of X is a collection X = {X i } i∈I of P -indexed filtrations such that for each p ∈ P , The definition of the nerve above extends immediately to yield a nerve filtration Nrv (X) associated to any cover of a filtration.Theorem 3 (Persistent Nerve Theorem [3]).Let P be a poset, X a P -indexed Euclidean filtration, and X a good cover of X.There exists a diagram of objectwise homotopy equivalences As shown in [3], the intermediate filtration ∆X in the statement of the theorem can be taken to be a homotopy colimit of a diagram constructed from X, just as in the proof of the Persistent Nerve Theorem for open covers [18].Note that if P is a singleton set, the Persistent Nerve Theorem specializes to the classical version of the Nerve Theorem mentioned above.

Multicovers
As indicated in the introduction, the multicover bifiltration of a finite set where as above, Given this, one may wonder why we allow k to take the value 0 in our definition of Cov.As we will see in Section 4, this turns out to be convenient for comparing Cov to the rhomboid bifiltration.
As a first step towards constructing a simplicial model of Cov, we identify a good cover for Cov r,k for fixed r and k.The set of all order-k Voronoi regions yield a decomposition of R d into closed convex subsets having the same k closest points of A in common.This decomposition is called the order-k Voronoi diagram [2,30].We denote it as Vor k .
For any r ∈ R and k ∈ N, intersecting each order-k Voronoi region with the corresponding multicovered region of fixed radius r yields the following good cover of Cov r,k .
For an illustration, see Figure 3.The nerve of Vor r,k , which we will denote Del r,k , is called an order-k Delaunay complex.By the Nerve Theorem, Del r,k and Cov r,k are homotopy equivalent.Note that Del r,1 is the alpha complex of radius r [25,27], whereas Del r,0 is a single point for all r ∈ [0, ∞).A different but related concept is the order-k Delaunay mosaic, which is the geometric dual of the order-k Voronoi diagram [28]; see Section 4.

A simplicial Delaunay bifiltration
For fixed r ≥ 0, we have inclusions but there are no analogous inclusions Del r,k → Del r,k−1 between the higherorder Delaunay complexes.Indeed, we do not even have inclusions of the vertex sets.Consequently, (Vor r,k ) (r,k)∈[0,∞)×N is not a bifiltered cover of Cov.
To correct for this, we first define Del r,k := Nrv (Vor r,k ∪ Vor r,k+1 ) .
See Figure 4  For an illustration, see Figure 5.Note that S-Del r,k is generally not equal to the nerve of the union of all Vor r,i , i ≥ k, which is a much larger object.

Weak equivalence of Cov and S-Del
We next prove that Cov and S-Del are weakly equivalent.The proof will use the following version of the usual mapping telescope construction [33, Section 3.F]: given any sequence of continuous maps let M C , the mapping telescope of C, be the quotient of the disjoint union Theorem 4. The multicover bifiltration Cov is weakly equivalent to S-Del.
Proof.Our proof strategy is similar to ones used for sparse filtrations [13,48].
We will observe that S-Del is isomorphic to the nerve of a good cover of a bifiltration X which is weakly equivalent to Cov.The result then follows from the Persistent Nerve Theorem.Let X r,k be the mapping telescope of the sequence Letting r and k vary, the spaces X r,k assemble into a ([0, ∞) × N op )-indexed bifiltration X, and the deformation retractions X r,k → Cov r,k assemble into an objectwise homotopy equivalence X → Cov.Letting denote the quotient map, we have a good cover , and the collection of all such covers as r and k varies assembles into a good cover X of X.
To finish the proof, it suffices to observe that Nrv(X) is isomorphic to S-Del.First, note that vertices of Nrv(X r,k ) correspond bijectively to non-empty elements of n i=k Vor r,i , as do vertices of S-Del r,k .This gives us a bijection from the vertex set of Nrv(X r,k ) to the vertex set of S-Del r,k .In fact, this bijection is a simplicial isomorphism ϕ r,k : Nrv(X r,k ) → S-Del r,k , because for all if and only if both of the following conditions hold: The isomorphisms ϕ r,k are natural in r and k, so assemble into an isomorphism from Nrv(X) to S-Del.

Truncations of S-Del
When the point cloud is large, it may be computationally difficult to construct the full bifiltration S-Del.We instead consider, for K ∈ N, the bifiltration S-Del ≤K constructed in the same way, but disregarding all order-k Voronoi cells with k > K, S-Del ≤K r,k := Del r,i .

Size of S-Del
By the size of a bifiltration X, we mean the number of simplices in the largest simplicial complex in X.If the sites A ⊂ R d are not in general position, the size of S-Del ≤K can be huge; indeed, if all points of A lie on a circle in R 2 and r is at least the radius of this circle, then Del r,k has 2 ( n k ) simplices, so S-Del ≤K is at least as large.However, if A is in general position, then the situation is far better: Proposition 5. Let A ⊂ R d be a set of n sites in general position, with a constant dimension d.Then S-Del ≤K has size O(n ).In particular, S-Del has size O(n d+1 ).
In contrast to Proposition 5, the number of vertices in Sheehy's subdivision-Čech model [46] grows exponentially with n, regardless of whether A is in general position.
In brief, the idea of the proof is to bound the number of maximal simplices in S-Del ≤K using a bound of Clarkson and Shor on the number of Voronoi vertices at levels ≤ k [20].The result then follows by observing that the dimension of the complex is a constant that only depends on d.However, this dependence on d is doubly exponential, so the O-notation hides a large factor if d is large.We now give details of the proof.
For k ∈ N, and a bifiltration X, we let X ∞,k denote the largest simplicial complex of the form X r,k , provided this exists.
In the following lemmas, A ⊂ R d is assumed to be in general position.
Lemma 6.For all k ∈ N, we have Proof.A simplex σ ∈ Del ∞,k is a set of order-k Voronoi regions with a point of common intersection.Let x be such a point.Each order-k Voronoi region R ∈ σ is indexed by a subset of A with k elements, which we denote as sites(R).Let Lemma 7.For all K ∈ N, we have Proof.We first observe that To finish the proof, we will show that the vertex set of Vor k+1/2 is the union of the vertex sets of Vor k and Vor k+1 .It suffices to show that for any cell σ ⊂ R d of Vor k with dim(σ) < d, σ is contained in a cell of Vor k+1 ; then Vor k+1/2 does not subdivide the (d − 1)-skeleton Vor k , implying that no new vertex is created.Let x ∈ σ.As in Lemma 6, we let and we note that |A in | < k.Moreover, since dim(σ) < d, we have |A on | ≥ 2. The set of order-k Voronoi regions containing x (and hence any y ∈ σ) is Therefore A on and A in are independent of our choice of x ∈ σ.Hence, for any x ∈ σ, the set of order-(k + 1) Voronoi regions containing x is The intersection of these Voronoi regions is the closure of a cell of Vor k+1 which contains σ.(In fact, similar reasoning shows that if |A in | + |A on | ≥ k + 2, then σ is a cell in Vor k+1 , though this is not needed for the argument.) Del ∞,k , Lemma 8 implies that the number of its maximal simplices is at most where V ≤K is the number of Voronoi vertices at levels ≤ k.Since the number of simplices of a d-dimensional simplicial complex with m maximal simplices is at most 2 d+1 m, Lemma 7 implies that the size of S-Del ≤K is bounded above by 2 2 Since the dimension bound only depends on the constant d, it is therefore enough to bound V ≤K .By a result of Clarkson and Shor [20], we have that e ≤K , the number of Voronoi regions (i.e., d-dimensional Voronoi cells) at levels ≤ K is ).
This bound hides a constant that depends doubly exponentially in d.Moreover, by a counting argument appearing in [24, Theorem 3.3], we have that the total size of the Voronoi diagram at level k is bounded by with e i the number of Voronoi regions at level i.The same bound applies to V k as well.A simple calculation shows that V ≤K is then bounded by (2d−2)e ≤K = O(e ≤K ). 4 The rhomboid bifiltration

The rhomboid tiling
Let A ⊆ R d be a set of n sites in general position, and S an arbitrary (d − 1)sphere in R d .Then S yields a decomposition A = A in A on A out with A in the sites in the interior of S, A on the sites on the sphere, and A out the sites in the exterior of S. We define the combinatorial rhomboid of S to be the collection of sets We call elements of ρ S combinatorial vertices, and call the (combinatorial) rhomboid tiling of A. Elements of Rhomb(A) are called rhomboids.Since A is fixed throughout, we write Rhomb instead of Rhomb(A).As observed in [28], the combinatorial rhomboid tiling can be geometrically realized as a polyhedral cell complex [37,Def 2.38].For that, a combinatorial vertex {a 1 , . . ., a k } (where a 1 , . . ., a k are sites in R d ) is embedded as ( k i=1 a i , −k) in R d+1 .We call k the depth of the vertex.Embedding a combinatorial rhomboid as the convex hull of its embedded vertices yields an actual rhomboid in R d+1 whose dimension equals the cardinality of A on in the corresponding partition of A. The collection of these rhomboids is the (geometric) rhomboid tiling for A. We illustrate the construction in Figure 6.In what follows, we identify vertices and rhomboids with their combinatorial description.In particular, we will use Rhomb both for the combinatorial and the geometric rhomboid tiling.(d+1)!(n + 1) d+1 ≤ 2(n + 1) d+1 .
For any rhomboid ρ ∈ Rhomb, we let r ρ = inf {r | there exists a sphere S of radius r such that ρ S = ρ}.
It is easily checked that if ρ is a subset of ρ, then r ρ ≤ r ρ , and therefore, for any r ≥ 0, the sublevel set Rhomb r = {ρ ∈ Rhomb | r ρ ≤ r} is also a polyhedral complex.

Slicing
Next, we slice the rhomboid tiling by cutting every rhomboid along the hyperplanes {x ∈ R d+1 | −x d+1 = k} with k = 0, . . ., n.In this way, a rhomboid decomposes into its intersections with these hyperplanes and with slabs of the form {x ∈ R d+1 | k ≤ −x d+1 ≤ k + 1}.The resulting polyhedra again form a polyhedral complex that we call the sliced rhomboid tiling S-Rhomb.We refer to its cells as sliced rhomboids.
For a sliced rhomboid ρ, we define k ρ as the minimum depth among its vertices.Moreover, there is a unique (unsliced) rhomboid ρ of smallest dimension that contains ρ, and we define r ρ := r ρ .Define S-Rhomb r,k := {ρ ∈ S-Rhomb | r ρ ≤ r, k ρ ≥ k} and observe that for r ≤ r and k ≥ k , we have S-Rhomb r,k ⊆ S-Rhomb r ,k .Hence, (S-Rhomb r,k ) (r,k)∈[0,∞)×N op is a bifiltration of combinatorial cell complexes.Again, we will abuse notation and use the symbol S-Rhomb both for the sliced rhomboid tiling and the bifiltration (S-Rhomb r,k ) (r,k)∈[0,∞)×N op .As shown in [28], the restriction of S-Rhomb to cells in the hyperplane −x d+1 = k is the order-k Delaunay mosaic, i.e., the geometric dual of the order-k Voronoi diagram.

Comparison of S-Rhomb and S-Del
The next lemma establishes a close relationship between the bifiltrations S-Rhomb and S-Del.It is closely related to [28,  Proof.For the first part, note that a set of sites v = {a 1 , . . ., a k } is a vertex in S-Rhomb r,k if and only if k ≥ k and for all r > r there is a sphere S of radius at most r whose associated decomposition A = A in A on A out satisfies Note that (3) holds if and only if the center of S has v among its k closest sites, which is equivalent to the condition that the order-k Voronoi region Vor(v) contains the center of S. Thus, such a sphere S exists if and only v is a vertex of Del r ,k .Thus, v ∈ S-Rhomb r,k if and only if k ≥ k and v ∈ Del r ,k for all r > r.But the latter holds if and only if v ∈ S-Del r,k because, first, the vertex sets of S-Del r ,k and k ≥k Del r ,k are equal and, second, v ∈ S-Del r ,k for all r > r implies v ∈ S-Del r,k .Thus, the vertex sets of S-Rhomb r,k and S-Del r,k are equal, as claimed.
For the second part, let v 1 , . . ., v m denote the vertices of a sliced rhomboid in S-Rhomb r,k .Let ρ denote the smallest rhomboid of Rhomb containing this sliced rhomboid.For any r > r, there exists a sphere S of radius at most r with ρ S = ρ.Let x denote the center of S and A = A in A on A out be the decomposition with respect to S. Now, each v i is the union of A in with a subset of A on , and hence, the point x belongs to the Voronoi region of v i .Since i is arbitrary, it follows that all Voronoi regions intersect, and x has distance ≤ r to each site in each v i , so v 1 , . . ., v m span a simplex in S-Del r ,k .Since this holds for all r > r, this simplex is also contained in S-Del r,k .
For the third part, consider vertices v 1 , . . ., v m that span a simplex in S-Del r,k .Assume that some v i has order k ≥ k and that the remaining vertices are of order k or k − 1.Since v 1 , . . ., v m span a simplex, the corresponding higherorder Voronoi regions, intersected with balls of radius r around the involved sites, have non-empty intersection.Let x be a point in this intersection.
If k = 0, then m = 1 and v 1 = ∅.As {∅} is a combinatorial rhomboid of dimension 0, the desired result holds in this case.If k = 1 and x is a site, then either m = 1 and v 1 = {x}, or else m = 2, in which case {v 1 , v 2 } = {{x}, ∅}.In either case, the desired result again holds.
Otherwise, there is smallest sphere S centered at x that includes k sites (either on the sphere or in its interior).S induces a partition A = A in A on A out and a rhomboid ρ.Each vertex v i of order k must contain all sites of A in , and some subset of the sites of A on , meaning that v i is in ρ.Furthermore, at least one site lies on S; otherwise, there would be a smaller sphere.This implies that each vertex v j of order k − 1 also has to contain all sites of A in and some subset of A on .Thus, all vertices lie in ρ, and in particular in its (k − 1, k )-slice.Finally, observe that because one of the sites lies on S, the radius of S is the distance of that site to x, which is at most r.
Remark 11.Parts 1 and 2 of Lemma 10 establish that we have a vertexpreserving injection J from the cells of S-Rhomb r,k to the simplices of S-Del r,k .Moreover, the third part of Lemma 10 implies that J restricts to a bijection from the maximal cells of S-Rhomb r,k to the maximal simplices of S-Del r,k .Figure 7 illustrates that J itself needn't be a bijection.
We note that J does not preserve dimension: For ν a cell of S-Rhomb spanned by vertices of cardinality k (i.e., a cell in the order-k Delaunay mosaic), we have dim(ν) ≤ d.But if d ≥ 3, it can be that dim(J (ν)) > d, even when the sites are in general position.For a cell ν of S-Rhomb r,k spanned by vertices of cardinality k and k + 1, we may have dim(ν) = dim(J (ν)) even for d = 2, see Figure 7.
Theorem 12.The bifiltrations S-Rhomb and S-Del are weakly equivalent.
Proof.We define good covers of both bifiltrations: First, for S-Rhomb r,k , we choose the cover that consists of all its cells.This is a good cover because the for some intermediate bifiltration ∆U.Moreover, we obtain a cover V r,k of S-Del r,k whose elements are the simplices spanned by the vertices of the sliced rhomboids in S-Rhomb r,k .By the second part of Lemma 10, these cover elements indeed exist.Moreover, in view of Remark 11, every maximal simplex of S-Del r,k is an element of V r,k .We thus obtain a good cover V of S-Del.Applying the Persistent Nerve Theorem again, we obtain objectwise homotopy equivalences S-Del ∆V Nrv(V).
Finally, Nrv(U) and Nrv(V) are isomorphic: The elements of U and of V are in 1-to-1 correspondence, with corresponding cover elements having the same vertex set.In both cases, an intersection of cover elements is non-empty if and only if the elements share a vertex, which is determined by their vertex sets.Hence, we have objectwise homotopy equivalences

Unslicing the rhomboid
Next, we define a bifiltration on the (unsliced) rhomboid tiling.Recall that for a rhomboid ρ, we already defined r ρ as the radius of the smallest sphere that gives rise to that rhomboid.As in the sliced version, we define k ρ as the minimal depth among the vertices of ρ, and This yields a bifiltration (Rhomb r,k ) (r,k)∈[0,∞)×N op , which we denote by Rhomb.
Lemma 13.The bifiltrations Rhomb and S-Rhomb are weakly equivalent.
Proof.For a rhomboid ρ in Rhomb, set k min as the minimum depth and k max as the maximum depth among the vertices in ρ.Note that k ρ = k min .r and k fixed, we say ρ is dangling if r ρ ≤ r and k min < k < k max .If ρ is dangling then ρ ∈ Rhomb r,k , but some of the slices of ρ are contained in S-Rhomb r,k .
In fact, all cells of S-Rhomb r,k not contained in (the geometric realization of) Rhomb r,k are of this form.For example, taking k = 2 and r very large, the shaded rhomboid {c, bc, cd, bcd} of Figure 6 is dangling.S-Rhomb r,2 contains the cell {bc, cd, bcd} but Rhomb r,2 does not.
Observe that there is a deformation retraction of S-Rhomb r,k onto Rhomb r,k which, for each dangling rhomboid ρ, "pushes" onto the boundary of ρ; for instance, in the example above {bc, cd, bcd} is pushed onto {bc, bcd} ∪ {cd, bcd}.Thus, for every choice of r and k, the inclusion Rhomb r,k → S-Rhomb r,k is a homotopy equivalence.Moreover, these inclusions commute with the inclusion maps in Rhomb and S-Rhomb, hence define an objectwise homotopy equivalence.
Combining the previous lemma with Theorem 12 and Theorem 4 yields the following result: Theorem 14.The bifiltrations Rhomb and Cov are weakly equivalent.
Remark 15 (Size of the Rhomboid Bifiltration).In view of Proposition 9, Rhomb has at most 2(n + 1) d+1 = O(n d+1 ) cells.One can also bound the size of a truncated version of Rhomb, defined analogously to the truncation of S-Del considered in Proposition 5. Indeed, Rhomb is clearly smaller (in terms of number of cells) than S-Rhomb, and by Remark 11, S-Rhomb is at least as small as S-Del.Moreover, this extends to truncations of these bifiltrations.Thus, the size bound of Proposition 5 also holds for truncations of Rhomb.

Computation
In [29,43], a relatively simple algorithm is given for computing the rhomboid bifiltration.(In fact, [29] explicitly considers only the computation of the rhomboid tiling and the radius r ρ of each rhomboid ρ; the rhomboid bifiltration is not mentioned.But it is trivial to extend the algorithm to compute the depth k ρ of each rhomboid ρ, thus computing the rhomboid bifiltration.) We briefly outline the approach.Given a (d + 1)-dimensional rhomboid ρ, let σ denote the intersection of ρ with the hyperplane k = k ρ + 1.We call σ the generation-1 slice of ρ.Note that σ is a d-simplex.Given the combinatorial vertices of σ, we can easily recover ρ [29, Lemma 2].
The algorithm of [29] computes the rhomboid tiling by computing the generation-1 slices of all (d+1)-dimensional rhomboids.For each k in increasing order, a weighted Delaunay triangulation W k is computed which triangulates the order-k Delaunay mosaic and has the same vertex set.Given the vertex set, W k can be computed via any algorithm for weighted Delaunay triangulation computation, e.g., via a d + 1-dimensional convex hull computation [8,Section 4.4.4].We explain below how the vertex set is computed.
A simple combinatorial criterion [29, Lemma 3] tells us whether a d-simplex in W k is a generation-1 slice of a rhomboid.Thus, one can efficiently identify all generation-1 slices in the triangulation by iterating through all the d-simplices of W k .In this way, we identify all (d + 1)-dimensional rhomboids ρ with It remains to explain how the vertex set of W k is computed.The vertices of W 1 are just the sites A. For k ≥ 2, [29, Lemma 3] establishes that every vertex v in W k appears in a rhomboid ρ with k ρ ≤ k − 2. We thus discover ρ, and hence v, by the time we finish processing W k−1 .
Complexity The complexity of this algorithm is discussed in [29,Section 4].While an explicit runtime bound is not given, it is easy to extract naive bounds from the discussion; we now do so.We distinguish between two contributions to the runtime: 1. computing W k at all levels k, given the vertices, 2. checking, for each k, whether each d-simplex in W k is a generation-1 slice and if so, storing the corresponding rhomboid and its faces.The complexity of computing the triangulations W k depends on a choice of algorithm for computing weighted Delaunay triangulations.Some well-known algorithms have output-sensitive complexity bounds.For example, in the case d = 3, a weighted Delaunay triangulation of p points can be computed by the algorithm of [15] in time O((p + m) log 2 m), where m is the size of the output.In our setting, the total size of all the W k is O(n d+1 ) because the size of each W k differs from the size of the order-k Delaunay mosaic by at most a constant factor.Hence for d = 3, which is arguably the case of primarily interest, the total time to compute all of the triangulations W k is O(n 4 log 2 n).Therefore, the total cost of computing the rhomboid bifiltration is O(n 5 + n 4 log 2 n) = O(n 5 ).
For arbitrary d, the approach of [8, Section 4. ).In our setting, each vertex of each W k is a vertex of the rhomboid tiling, so there are a total of O(n d+1 ) vertices among all W k .Thus, for d ≥ 3 the time required to compute all W k is O(n (d+1) d 2 ), and the runtime of the full algorithm satisfies the same asymptotic bound.This bound is rather large, but it seems likely that it could be improved via a more careful analysis.
Implementation The above algorithm has been implemented in the software package rhomboidtiling1 [29].The code computes the sliced and unsliced bifiltrations S-Rhomb and Rhomb as well as their free implicit representations (FIREPs), i.e., chain complexes [42].rhomboidtiling is written in C++, using the Cgal library2 for geometric primitives.The current version accepts only 2-and 3-dimensional inputs, but all steps readily generalize to higher dimensions; adding support for higher-dimensional inputs is a matter of software design rather than algorithm development.That said, handing higherdimension inputs of practical size is likely to be computationally expensive.

Experiments
We performed experiments on point sets in R 2 and R3 .We provide a brief summary here; for detailed results, see Appendix A. We sampled points uniformly at random from [0, 1] 2 and [0, 1] 3 , from a disk, from an annulus, and from an annulus with noise added.We computed the rhomboid bifiltrations Rhomb ≤K and S-Rhomb ≤K .We then used mpfree 3 to compute minimal presentations of 2-parameter persistent homology of our bifiltrations.
In one set of experiments, we found that Rhomb ≤K is up to 47% smaller than S-Rhomb ≤K , and can be computed more than 20% faster.The experiments suggest that the relative performance of Rhomb ≤K improves with increasing K.
We investigated the size of Rhomb ≤K , varying the sample size and the threshold K.For d = 2, our experiments show a clear subquadratic growth of the size of Rhomb ≤K and its FIREP with respect to increasing K.For d = 3, the growth is clearly subcubic.These observations also extend to time complexity.Letting the number of points increase, the size of Rhomb ≤K and its FIREP shows roughly linear growth for both space dimensions, with a slight superlinear tendency.Again, we observed the same behavior for the computation time.
We conclude this section with a data visualization enabled by the ideas of this paper: For i ≥ 0, the i-th Hilbert function assigns to each parameter (r, k) ∈ R × N the rank of i-th homology module of Cov r,k (with coefficients in some fixed field).The Hilbert functions are well known to be unstable invariants.Nevertheless, their visualization can give us a feel for how the Lipschitz stability property of the multicover bifiltration established in [6] manifests itself in random data.Figure 8 shows a few examples, plotted using rivet4 .

Conclusion
We have introduced a simplicial model for the multicover bifiltration, as well as a polyhedral model based on the rhomboid tiling of [28].For a data set of size n in R d with d constant, the size of both constructions is O(n d+1 ).The size can be controlled by thresholding the parameter k of the multicover bifiltration.An algorithm of [29] computes the rhomboid bifiltration, and an implementation is available.In our experimental results, this approach scales well enough to suggest that practical applications could soon be within reach.A natural next step is to begin exploring the use of the multicover bifiltration on real world data.
To obtain our combinatorial models of the multicover bifiltration, we begin with a zigzag of filtrations, and then straighten it out by taking unions of prefixes.Notably, one could in principle compute the persistent homology modules of the multicover bifiltration without straightening out the zigzag, by inverting the isomorphisms on homology induced by the inclusions Del r,k → Del r,k .It seems plausible that this approach could be computationally useful.
We are curious to learn which indecomposables typically arise in the persistent homology modules of multicover bifiltration, and our approach could be used in conjunction with existing algorithms [22,35] to study this.It would also be interesting to investigate whether there is an interplay between the geometry of a space and the multicover bifiltration of a noisy sample of this space; we wonder if invariants of the bifiltration encode additional information about geometric properties, such as the reach or differentiability.
Our experiments show a significant increase in the size of our models of multicover bifiltration for increasing K.This suggests the need for refinements to our algorithmic approach in order to handle large values of K. Aside from the truncations considered in this paper, there are a couple of promising ways forward: One could construct a coarsened bifiltration where some values of k are skipped.Alternatively, one could make use of the inductive nature of our constructions: for the step from k to k + 1, one does not need information about the bifiltrations at indices j < k.Therefore, one could provide the bifiltration as an output stream without storing it completely in memory.Subsequent algorithmic steps would then have to be implemented as streaming algorithms as well.

A Details on experiments A.1 Implementation
Let us mention an important technicality in our pipeline for computing minimal presentations: In order to limit the size of the minimal presentations, we "snap" the radius values of all generators and relations onto a set of 100 evenly spaced points in R. The values of the parameter k are left unchanged.This snapping is done only after the minimal presentation is computed.The snapping process in fact can make a minimal presentation non-minimal, so after snapping, we re-minimize the presentation.All reported results below are for minimal presentations computed using this pipeline.The Hilbert functions shown in Figure 8 were also computed from such "snapped" presentations.

A.2 Experimental output
We present the concrete outcome of some of our experiments.All results are averaged over 5 runs with independently generated data sets.The sizes reflect the number of elements of the corresponding set, and the times were measured in seconds.
We give a brief overview of the experiments, referring to the tables for further details.We were curious about the practical improvements from S-Rhomb to Rhomb.We documented these for a few values of n and K in the plane.The sizes of their truncated versions both grow linearly in the number of points.See Table 1 for more refined results.In further experiments, we only used Rhomb.
Table 2 shows the behavior of fixed K and an increasing number of uniformly sampled points.Conversely, Table 3 documents the behavior of Rhomb in an experiment with a fixed number of points and increasing K.Both in dimension 2 and 3, we investigated the size of the bifiltration, the size of the FIREPs, and the size of minimal presentations thereof.We also kept track of the time needed for the computations.
Finally, we wondered how the measurements change when data sets are sampled from a particular shape.As an example, we sampled points from an annulus with random but bounded perturbations and added uniform background noise to it.Letting both the range of the perturbations and the portion of the background noise vary, we tracked the size of the minimal presentations in Table 4  We sampled points of an annulus around a circle of radius 0.25.In total, we have 10,000 points whereas p% of these points are uniform noise in the surrounding box [0, 1] 2 .The other points are sampled with a random perturbation per coordinate bounded by a number err.We considered K = 8 as maximal value for k.The size of the snapped minimal presentation increases when adding more uniform noise.It may increase more drastically within a certain range of p, i.e., for p ∈ {1, 4}, starting at about at 0.12.We also observed a considerable variance of the individual results in such areas.In particular, the size of the snapped minimal presentations is neither a linear, nor a sub-or superlinear process.We regard this process mostly as a property of the snapping technique.Furthermore, when p is not too big, the perturbations around each sampled point can be quite high, i.e., for p ≤ 4 and err = 0.16, the snapped minimal presentations are still of relative size 5.7% and 8.3%, respectively.Note that the samples only stay inside the surrounding box [0, 1] 2 if err ≤ 0.25.

Figure 3 :
Figure 3: Left: The 2-fold cover of some points with respect to a certain radius overlapped with its Voronoi diagram of order 2. Vor is combinatorially simpler than Cov.Right: The corresponding 3-fold cover overlapped with its Voronoi diagram of order 3.

Figure 4 :
Figure4: Left: The Delaunay complexes of order 2 and 3 of our running example.Right: The construction of the simplicial complex Del r,2 .Del r,2 consists of the Delaunay complexes Del r,2 and Del r,3 , and additional mixed simplices connecting these.This connection arises from intersections of the 2-, and 3-fold covers restricted to their Voronoi diagrams of order 2 and 3, respectively.

Figure 5 :
Figure5: S-Del r,2 is the union of all Del r,i , i ≥ 2. We have that S-Del r,k = ∅ for k greater than the number of sites in A. Thus, in this case, Del r,k is empty for k ≥ 5.
for an illustration.Note that Vor r,k ∪ Vor r,k+1 covers the same space as Vor r,k and that we have inclusions Del r,k+1 Del r,k Del r,k .Letting n := |A|, we see that Del r,n = Del r,n and Del r,k = ∅ for all k > n.We now define a [0, ∞) × N op -indexed simplicial bifiltration S-Del by S-Del r,k := i≥k Del r,i = n i=k Del r,i .

d+1 2 , 1 d+1 2 ,Lemma 8 .
by noting that the dimension on the left is upper bounded by the sum of the maximal number of order-k Voronoi regions and order-(k + 1) Voronoi regions that meet at a fixed point in R d .Since these two summands equal dim(Del ∞,k ) − 1 and dim(Del ∞,k+1 ) − 1, we have the first inequality, and the second one follows from the previous Lemma 6.Since each simplicial complex in S-Del ≤K is a union of simplicial complexes Del r,k , it follows that dim(S-Del ≤K ) ≤ 2 d + If n > d, then for any k ∈ N, the number of maximal simplices in Del ∞,k is at most V k + V k+1 , where V k , V k+1 denote the number of Voronoi vertices at level k or k + 1, respectively.Proof.Recall that Vor k and Vor k+1 denote the order-k and order-(k+1) Voronoi diagrams, respectively.Let Vor k+1/2 denote their overlay, i.e., the polyhedral subdivision of R n induced by the closed polyhedra of the form R ∩ S, where R and S are Voronoi regions of order k such that dim(R ∩ S) = d.The subdivision Vor k+1/2 is called the Voronoi diagram of degree-(k +1)[28].It has been studied in[30].Maximal simplices of Del ∞,k correspond bijectively to cells in Vor k+1/2 with no proper faces.We claim that the only such cells are vertices.If k ∈ {0, n − 1, n}, then either Vor k or Vor k+1 contains only the region R d , and it is clear that the claim holds.For 0 < k < n − 1, we observe that since A is in general position and n > d, each cell of either Vor k or Vor k+1 is bounded by a vertex, which implies that each cell of Vor k+1/2 is bounded by a vertex.The claim follows.

Figure 6 :
Figure6: The rhomboid tiling of 5 points on the real line.The highlighted 2-rhomboid ρ defined by A in (ρ) = {c} and A on (ρ) = {b, d} is the convex hull of the points c, bc, cd, and bcd, simplifying the labels here and, e.g., writing bcd instead of the cell complex associated to {b, c, d}.The horizontal line at depth k intersects the tiling in the order-k Delaunay mosaic.

3 .
Theorem 1 and Lemma 2].Lemma 10.For all (r, k) ∈ [0, ∞) × N, 1.The vertex sets of S-Rhomb r,k and S-Del r,k are equal.2. The vertices of each sliced rhomboid in S-Rhomb r,k span a simplex in S-Del r,k .The vertices of each simplex in S-Del r,k are contained in a sliced rhomboid of S-Rhomb r,k .

Figure 7 :
Figure7: An illustration of the difference between S-Del and S-Rhomb for three points in the plane.The order-1 Voronoi regions of the points {x}, {y}, and {z} intersect in c, as do the order-2 Voronoi regions of {x, y}, {y, z}, and {x, z}.Consequently, S-Del contains a 5-simplex, but the corresponding cell in S-Rhomb on the same vertex set is a 3-dimensional triangular skew prism.Many simplices in S-Del do not correspond to any cell in S-Rhomb, e.g., the 1-simplex {x, yz}.
The latter requires O(k) time per d-simplex in W k .Hence, since the rhomboid tiling has size O(n d+1 ) (seeRemark 15), the total time required over all d-simplices is O(n d+2 ).

4 . 4 ]
computes the weighted Delaunay triangulation of p points in R d in time O(p log p + p d 2

Figure 8 :
Figure 8: An illustration of the first Hilbert function of the multicover bifiltration, using grayscale shading.The instances are samples of an annulus (top), a noisy annulus (middle), and a disk (bottom).The sample size is 100 in the left column, and 200 in the right column.Darkness of the shading is proportional to the value of the Hilbert function, up to some maximum value, above which the shading is taken to be black; the lightest non-white shade of gray corresponds to a Hilbert function value of 1.
By minimality of s, we have |A in | < k, and by the assumption of general position, we have |A on | ≤ d + 1.Moreover, for all R ∈ σ, we have A in ⊂ sites(R) and sites .

Table 3 :
Results on 500 uniformly sampled points in the unit square and unit cube.The sizes of the bifiltration and its FIREP grow clearly subquadratic for d = 2, and clearly subcubic for d = 3.The running times seem close to being quadratic and cubic in dimensions 2 and 3, respectively.We also measured the relation of the sizes of the snapped minimal presentations and FIREPs.