A Sparse Multicover Bifiltration of Linear Size

For $k=1$, the $1$-cover is the union of balls of radius $r$ around points of $X$. The union of balls is commonly used in reconstructing submanifolds from samples [46, 4, 3] and is a cornerstone of persistent homology methods in topological data analysis, being homotopy equivalent to the Čech complex and the alpha complex [31, 32], and related to the (Vietoris-)Rips complex. The union of balls assembles into a $1$-parameter filtration: by increasing the scale parameter $r$ we obtain inclusions ${\mathrm{Cov}}_{r,1}\hookrightarrow {\mathrm{Cov}}_{r^{\prime },1}$, for $r\le r^{\prime }$. Taking the homology of each ${\mathrm{Cov}}_{r,1}$, we obtain a persistence module, an algebraic descriptor of the topology of the filtration across the different scales. Such a persistence module is stable to perturbations of the points [26, 25], but is not robust to outliers (already appreciated in the first column of Figure 1) and insensitive to differences in density in the point cloud.

There are multiple methods that address the lack of robustness to outliers and changes in density within the $1$-parameter framework. These include density estimation [47, 20, 8], distance to a measure [18, 38, 39, 19, 2, 1, 11, 10], and subsampling [5]. However, they all depend on choosing a parameter. The problem is that it is not clear how to choose such a parameter for all cases, and such a choice might focus on a specific range of scales or densities. We refer to [7, Section 1.7] for a complete overview of the methods and their limitations.

It is then natural to consider constructions over two parameters, scale and density. Examples include the density bifiltration [13], degree bifiltration [44], and the multicover, closely related to both the distance to a measure [18] and $k$-nearest neighbors [53]. The advantage of the multicover is that it does not depend on any further choices (like choosing a density estimation function) and that it is robust to outliers, as a consequence of its strong stability properties [7].

However, one currently problematic aspect of the multicover is its computation. Sheehy [53] introduced an influential simplicial model of the multicover called the subdivision Čech bifiltration, based on the barycentric subdivision of the Čech filtration. It has exponentially many vertices on the number of input points, making its computation infeasible. A crucial ingredient in the theory is the multicover nerve theorem [53, 14, 7] that establishes the topological equivalence (see Section 3.1 for a precise definition) of the subdivision Čech and multicover bifiltrations. Such a multicover nerve theorem has its analogue for the sparse multicover, the sparse multicover nerve theorem (Theorem 14).

Dually to an hyperplane arrangement in ${ℝ}^{d+1}$, Edelsbrunner and Osang [34, 33] define the rhomboid tiling and use it to compute, up to homotopy, slices (fixing $k$ and varying $r$ or viceversa) of the multicover bifiltration. Corbet, Kerber, Lesnick and Osang [28, 27] build on it to define two bifiltrations topologically equivalent to the multicover, of size $O({|X|}^{d+1})$.

Recently, Buchet, Dornelas and Kerber [12] introduced an elegant $(1+\epsilon )$-approximation of the multicover whose $m$-skeleton has size that is linear on $|X|$ but that incurs an exponential dependency on the maximum order $k$ we want to compute. A crucial difference between our work and the Buchet-Dornelas-Kerber (BDK) sparsification is that BDK work at the level of intersection of balls, while we work directly at the level of balls. This starting point is what allows us to ultimately obtain a construction of linear size, independently of $k$.

In addition, BDK work by freezing the lenses: at a certain scale they stop growing. It is not clear (or is technically challenging, in the words of BDK) how to compute exactly the first scale at which freezing lenses first intersect – a problem they sidestep by discretizing the scale parameter (at no complexity cost). In contrast, by letting sparse balls grow slowly at a certain rate, rather than freezing them completely, we can compute their first intersection time exactly, as we explain in Section 5.3.

Persistent nets, defined below, are a reinterpretation of farthest point sampling, or greedy permutations, as known in the sparsification literature [16, 55, 12], and as such are guaranteed to exist by Gonzalez algorithm [37]. We will revisit these concepts in Section 5.

Let $(X,\partial )$ be a finite metric space. We say that a subset $S\subset X$ is a $r$-net for a radius $r\ge 0$, if it is

1. 1.

   a $r$-covering: for every $x\in X\setminus S$ there is an $a\in S$ with , and

2. 2.

   a $r$-packing: for every two $a,b\in S$ we have $\partial (a,b)\ge r$.

In the spirit of persistence, a persistent net $𝒮$ of $X$ is a collection of $r$-nets $𝒮(r)\subset X$, one for each $r\ge 0$, such that for any two $r^{\prime }\ge r$ we have $𝒮(r^{\prime })\subset 𝒮(r)$.

The insertion radius $\mathrm{ins}(x)\in ℝ\cup \{\mathrm{\infty }\}$ of a point $x\in X$, with respect to a persistent net $𝒮$, is the supremum of those $r$ such that $x\in 𝒮(r)$. There is a large enough $r$ so that $𝒮(r)$ is a single point (it has to be a $r$-packing), and so there is only one point $x$ with $\mathrm{ins}(x)=\mathrm{\infty }$.

In what follows, we work over a fixed finite subset $X\subset {ℝ}^d$, with ${ℝ}^d$ equipped with any norm. We also fix an error parameter $\epsilon >0$ and a persistent net $𝒮$ of $X$, that we often drop from the notation. We define the slowing time $\mathrm{slow}(x)$ of a point $x\in X$ by

For example, the radius function of a point $x\in X$ can be

but other options are useful for computation (Section 5.3). In any case, the sparse ball of $x$ has radius $r$ up to scale $\mathrm{slow}(x)$, it then starts to grow slowly with radius at most , until eventually disappearing at scale $(1+3\epsilon )\mathrm{slow}(x)$. We say that the sparse ball around $x\in X$ is slowed at scale $r$ if $r\ge \mathrm{slow}(x)$.

The functions $L_x(r)$ and $U_x(r)$ are defined in this way for technical reasons related to properties of the sparse multicover that we will prove shortly. For now, let us say that whenever a sparse ball disappears (that is, at scale $(1+3\epsilon )\mathrm{slow}(x)$), it is covered by another non-slowed sparse ball:

One of the intuitive principles behind the sparse multicover is that, at every scale $r\ge 0$, the non-empty sparse balls approximately cover the balls of radius $r$ around all the points of $X$. Below, we formally define a map that takes each point $x$ to the non-empty sparse ball that approximately covers the ball of $x$.

We will use the following lemma multiple times. It describes how the covering map inherits covering and packing properties from the persistent net.

We can write the $k$-cover ${\mathrm{Cov}}_{r,k}$ for a scale $r\ge 0$ as

In parallel, the sparse multicover uses sparse balls instead of balls, and the covering weight instead of the number of balls:

The sparse multicover is the bifiltration given by the sparse covers:

We now recall the fundamental definitions of the homotopy theory of filtrations, where filtrations are viewed as diagrams, as is standard; see, e.g., [7]. Let $F$ be a bifiltration of topological spaces indexed, in our case, by radii $r\in [0,\mathrm{\infty })$ and order $k\in \mathbb{N}$; that is, a collection of topological spaces $F_{r,k}$, one for each $r\ge 0$ and $k\in \mathbb{N}$, together with a continuous map $F_{(r,k)\to (r^{\prime },k^{\prime })}:F_{r,k}\to F_{r^{\prime },k^{\prime }}$ for every $r\le r^{\prime }$ and $k\ge k^{\prime }$, such that $F_{(r,k)\to (r,k)}$ is the identity and the following diagram commutes

A pointwise weak equivalence $F\stackrel{\cong }{\to }G$ between two bifiltrations $F$ and $G$, is a collection of weak homotopy equivalences ${\alpha }_{r,k}:F_{r,k}\to G_{r,k}$ such that the following diagram commutes for every $r\le r^{\prime }$ and $k\ge k^{\prime }$:

Pointwise weak equivalence is not a equivalence relation, but the following is. Two bifiltrations $F$ and $G$ are weakly equivalent if there exists a zigzag of pointwise homotopy equivalences that connect them:

We now define the sparse subdivision bifiltration and prove that is weakly equivalent to the multicover bifiltration. As already mentioned, the sparse subdivision bifiltration is an extension of Sheehy’s subdivision Čech bifiltration [53], which we review first.

The rest of this section is dedicated to proving the following theorem. We write $\mathrm{SSub}$ to refer both to the bifiltration and to the largest complex in the bifiltration, $\mathrm{SSub}={\bigcup }_{r,k}{\mathrm{SSub}}_{r,k}$, when no confusion is possible.

The following argument is analogous to the one used in previous sparse filtrations [54, 16]. Consider the map $\min :\mathrm{SSub}\to X$ that takes each simplex $\sigma ={\sigma }_0\subset \mathrm{\cdots }\subset {\sigma }_m$ to the minimum point $x$ in ${\sigma }_m$ with respect to the order given by their insertion radius $\mathrm{ins}(x)$, breaking ties arbitrarily. The points in each preimage of this map $\min :\mathrm{SSub}\to X$ are not too far from each other:

The sparse subdivision bifiltration, unlike the subdivision Čech bifiltration, is multicritical, which means that each simplex may have multiple critical grades. Let us briefly describe the critical grades of a simplex $\sigma ={\sigma }_0\subset \mathrm{\cdots }\subset {\sigma }_m\in \mathrm{SSub}$. Ordering the critical grades by scale $r$, the first critical grade has scale $r$ equal to the first scale at which the sparse balls around the points ${\sigma }_m$ intersect, and order $k$ equal to the sum of the covering weights of the points of ${\sigma }_0$ at such an scale. Further critical grades of $\sigma$ correspond to scales $r^{\prime }$ at which the covering weight of the points of ${\sigma }_0$ increase, until one of the associated sparse balls disappears. Still, there are not many critical grades, as another packing argument shows:

Let $(x_1,\mathrm{\dots },x_n)$ be an ordering of the $n$ points in $X$. Such a sequence is called greedy if for every $x_{i+1}$ and prefix $X_i≔\{x_1,\mathrm{\dots },x_i\}$, one has

where $d(x,X_i)={\min }_{x_j\in X_i}\Vert x-x_j\Vert$. The distance $r_{i+1}≔d(x_{i+1},X_i)$ is the insertion radius $r_{i+1}$ of $x_{i+1}$, where we take $r_1=\mathrm{\infty }$ by convention. A greedy permutation gives a persistent net $𝒮$ by taking ${𝒮}_r≔\{x_i|r\le r_i\}$.

Gonzalez [37] gives a method to compute a greedy permutation. It consists of $n$ phases, divided in an initialization phase and $|X|-1$ update phases. At the end of phase $i$, it maintains a subset $S_i=\{x_1,\mathrm{\dots },x_i\}\subset X$ and for each $x_j\in S_i$ its Voronoi set: the subset of points in $X$ that are closer to $x_j$ than to any other point in $S_i$. We break ties arbitrarily but consistently (say, by minimum index in $S_i$), so that the Voronoi sets form a partition of $X$. In the initialization phase we pick a random point $x_1$ and we set its Voronoi set to be every other point. In the update phase $i$, we obtain the new set of leaders $S_i$ by adding to $S_{i-1}$ the point $x_i$ whose distance to $S_{i-1}$ is maximal, and updating the Voronoi sets. This algorithm can be implemented in $O({|X|}^2)$-time by going over the whole set of points during each phase. An algorithm of Clarkson [23, 24], which also maintains the Voronoi sets, can be shown to take $O(|X|\log \mathrm{\Delta })$-time [40], and has been implemented [24, 51].

As mentioned, at every phase of the algorithm, each point $x\in X$ is assigned a Voronoi set of a point $x_j\in S_i$: its nearest neighbor among those in $S_i$. We call the sequence of such points $x_j$, as the algorithm progresses, the sequence of leaders of a point $x$. Such a sequence necessarily starts with $x_1$ and ends with the point itself.

From the sequence of leaders we can obtain the covering map. The covering sequence of a point is a subsequence of its sequence of leaders, in reverse order. To see this, note that the covering sequence of $x$ starts with $x$ itself, as in the reversed sequence of leaders. Then, when the sparse ball of $x$ disappears the next point in its covering sequence is its nearest neighbor among those points with non-slowed sparse ball. This set of points with non-slowed sparse ball is necessarily equal to one of the subsets $S_i$ we have encountered through the execution of the algorithm. In addition, balls are slowed in the same order as the reversed sequence of leaders. All in all, we conclude the following:

We also need to find those subsets of points whose sparse balls have a non-empty intersection, at any scale. By Lemma 16, any such subset $A$ is contained in a ball of radius $2(1+3\epsilon )\mathrm{slow}(x)$ around $x$, where $x$ is the point of minimal insertion point in $A$. Making use of this observation, and following [12], given a point $x\in X$, we call the points of higher insertion radius than the one of $x$ and within distance $2(1+3\epsilon )\mathrm{slow}(x)$ its friends. Note that the number of friends of a point is $O(1)$, because we are only considering points of higher insertion radius, as in Lemma 17. Once we have computed the friends of $x$, we can go over every possible subset and check whether the associated sparse balls intersect at some scale, which is done in the following section. This procedure is similar to other sparsification schemes [16, 12].

To compute the friends of every point, we can use a proximity search data structure as follows. First, we initialize it empty. In decreasing order of insertion radius (that is, the order given by the greedy permutation), we add the points one by one. After adding a point $x$, we query all the points within $2(1+3\epsilon )\mathrm{slow}(x)$ distance from $x$.

Cavanna, Jahanseir and Sheehy [16] gives such a data structure to compute all friends in such a way in time $O(|X|)$, again for fixed $\epsilon$ and $d$. In practice one can also use the related greedy trees [22], which have been implemented [51].

We now show how to compute the scale at which a subset of sparse balls first intersect, if they do so. We do it only for the Euclidean case; the $l^{\mathrm{\infty }}$-norm, which is also of interest, is an easier case, since it is enough to compute the intersection times pairwise. We show that this problem is of LP-type [45] and can be solved efficiently following the blueprint of the LP-algorithm for computing the smallest enclosing ball of balls [35], which is efficient in practice and implemented in CGAL [36].