Bifunction and Interlevel Delaunay Trifiltrations

Given a finite point cloud $X\subset {ℝ}^d$, the offset filtration of $X$ is the nested sequence of spaces $𝒪(X)={\left(𝒪{(X)}_r\right)}_{r\in {ℝ}^{+}}$, where ${ℝ}^{+}=[0,\mathrm{\infty })$ and $𝒪{(X)}_r$ is the union of closed balls of radius $r$ centered the points of $X$ (Figure 1). In topological data analysis (TDA), one often computes the persistent homology of $𝒪(X)$. However, directly computing $𝒪(X)$ is difficult, due to its non-combinatorial nature. Therefore, one usually instead computes the Delaunay filtration $𝒟(X)$ (also called the $\alpha$-filtration) [30, 28] (Figure 1), a filtration of the Delaunay triangulation $\mathrm{Del}(X)$ of $X$ that is topologically equivalent (i.e., weakly equivalent, see Definition 19) to $𝒪(X)$, and therefore has the same persistent homology. For low-dimensional Euclidean data, $𝒟(X)$ is the standard choice of filtration in TDA, because it is far more computationally efficient than alternatives like the Rips and Čech filtration. Bauer and Edelsbrunner [5] showed that $𝒪(X)$ is also weakly equivalent to a variant of $𝒟(X)$ called the Delaunay-Čech filtration and denoted $𝒟𝒞(X)$, where the birth index of each simplex in the filtration is the radius of its minimum enclosing ball.

In many TDA applications, the point cloud $X$ comes equipped with a function $\delta :X\to ℝ$ [20, 18]. It is then natural to seek a refinement of the persistent homology of $X$ which is sensitive to the structure of this function. Functions considered in the TDA literature include (co)density, eccentricity (centrality), discrete curvature [18], partial charge [17, 16], and time [46]. To extend the persistence pipeline to point clouds equipped with multiple functions, Carlsson and Zomorodian [18] proposed multiparameter persistence, and specifically, the sublevel offset bifiltration $𝒪(\delta )={\left(𝒪{(\delta )}_{s,r}\right)}_{(s,r)\in ℝ\times {ℝ}^{+}}$, where

Sublevel Čech and sublevel offset bifiltrations are defined analogously, and can be computed in essentially the same way as the usual Rips and Čech filtrations. Defining an analogous sublevel Delaunay bifiltration is more subtle because of the non-monotonicity of Delaunay complexes with respect to insertion of points. However, recent work [3] introduced sublevel Delaunay and sublevel Delaunay-Čech bifiltrations that are both weakly equivalent to the sublevel offset bifiltration and amenable to efficient computation. Experiments with a publicly available implementation showed that, in practice, computing the sublevel Delaunay-Čech bifiltration is only modestly more expensive than computing the ordinary Delaunay filtration. The work [3] is part of a large body of recent work aiming to realize multiparameter persistence as a computationally viable data analysis tool; for an introduction, see [14].

In this paper, we extend the recent results from [3] on sublevel-Delaunay bifiltrations to functions $\gamma :X\to {ℝ}^2$, where $X\subset {ℝ}^d$ is in general position. Specifically, we introduce the sublevel Delaunay and sublevel Delaunay-Čech trifiltrations of $\gamma$, denoted $𝒟(\gamma )$ and $𝒟𝒞(\gamma )$, respectively. We show that both trifiltrations have size $O({|X|}^{\lceil (d+1)/2\rceil +1})$, can be computed in time $O({|X|}^{\lceil d/2\rceil +2})$, and are weakly equivalent to the sublevel offset trifiltration, defined as follows:

Our approach to defining and computing the trifiltrations $𝒟(\gamma )$, $𝒟𝒞(\gamma )$ parallels the approach of [3] for $ℝ$-valued functions, but is necessarily more involved. To illustrate the problem, consider the four point clouds $X_a\subset X_b,X_c\subset X_d$ of Figure 2. While for any $r\in {ℝ}^{+}$, the offsets satisfy $𝒪{(X_a)}_r\subset 𝒪{(X_b)}_r,𝒪{(X_c)}_r\subset 𝒪{(X_d)}_r$, these containment relations do not hold for the corresponding Delaunay triangulations, nor for the Delaunay complexes of radius $r$. The construction of [3] applies when the point clouds are totally ordered by containment, but does not apply here because neither $X_b$ nor $X_c$ contains the other. In this paper, we adapt the approach of [3] to address this case.

As in the previous work on $ℝ$-valued functions, the most important part of our approach is the definition and computation of the maximum complex (i.e., colimit) of the trifiltrations $𝒟(\gamma )$ and $𝒟𝒞(\gamma )$. This complex is the same for both trifiltrations; we call it the incremental complex and denote it as $ℐ$ (Definition 3). The complex $ℐ$ has vertex set $X$, and contains the union of the Delaunay triangulations of all sublevel sets of $\gamma$ as a subcomplex. We give a characterization of $ℐ$ in terms of conflict pairs (see Definition 6), and use this to obtain a simple algorithm for computing $ℐ$, based on the standard Bowyer-Watson algorithm for incremental computation of Delaunay triangulations. We then introduce optimizations to this algorithm, exploiting the fact that Delaunay triangulations are local, in a sense made precise by Lemma 16.

For $\sigma \subset X$ non-empty, let $\gamma (\sigma )={\bigvee }_{x\in \sigma }\gamma (x)$, where $\bigvee$ denotes the join, i.e., coordinate-wise maximum. Given $ℐ$, the trifiltration $𝒟(\gamma )$ is defined by assigning each simplex $\sigma \in ℐ$ the birth index $(\gamma (\sigma ),{\omega }_{\sigma })$, where ${\omega }_{\sigma }$ is the incremental Delaunay radius, which we will define in Definition 3. $𝒟𝒞(\gamma )$ is defined by assigning each $\sigma \in ℐ$ the birth index $(\gamma (\sigma ),m_{\sigma })$, where $m_{\sigma }$ is the radius of the minimal enclosing ball of $\sigma$. To show that $𝒟(\gamma )$, $𝒟𝒞(\gamma )$, and $𝒪(\gamma )$ are topologically equivalent, we use an argument based on two-parameter mapping telescopes and functorial version of the nerve theorem [6]. We introduce a practical algorithm for computing $𝒟(\gamma )$ and $𝒟𝒞(\gamma )$, and implement the algorithm to compute $𝒟𝒞(\gamma )$. We perform experiments on point clouds in ${ℝ}^2$ and ${ℝ}^3$ with up to 16000 points. In these experiments, the memory usage scales nearly linearly, while run time scales nearly quadratically, suggesting that substantially larger computations are feasible. Our implementation is publicly available [37].

While there is a substantial recent literature on computing bifiltrations of point cloud and metric data in TDA [23, 3, 9, 44, 2, 4], to the best of our knowledge, our work is the first to seriously engage with the algorithmic aspects of trifiltrations. We hope that our work will lower the barrier to practical data analysis with 3-parameter persistence, which has seen very limited use so far.

Our specific motivation for introducing computable Delaunay-type models of sublevel-offset trifiltrations is threefold: First, we wish to enable persistence analysis which simultaneously considers a pair of functions on a point cloud, e.g., density and eccentricity. Indeed, if one is interested in doing a topological analysis of a function on a point cloud with noise or non-uniformities in density, it becomes important to take density into account, and one simple way to do this is to treat density as an additional filtration parameter.

Second, the definition of the sublevel-offset bifiltration involves a choice of filtering direction (i.e., whether we filter by sublevel or superlevel sets). While for some functions $\delta :X\to ℝ$, like codensity, it is natural to fix a direction, for other functions like time or partial charge, it can be unnatural to fix a filtering direction, and one would like a more symmetric construction. To this end, one can consider the interlevel-offset trifiltration of $\delta$, defined as follows:

Note that ${𝒪}^{\updownarrow }(\delta )$ is a special case of the sublevel-offset trifiltration, namely ${𝒪}^{\updownarrow }(\delta )=𝒪(\gamma )$, where $\gamma =(\delta ,-\delta )$. We are especially interested in computing the homology of the interlevel-offset filtration in the case that $\delta$ is a time function, as this could potentially be useful in the analysis of dynamic data such as that arising in the study of collective motion of animals [1], in the study of dynamical systems, or in the sliding-window analysis of time series [50]. While ${𝒪}^{\updownarrow }(\delta )$ is only one of several possible choices of a multifiltration that can be built from time-varying data (e.g., see [41, 55]), it is notable for its flexibility – unlike some alternatives, the definition does not require a consistent labeling of the points at each time step, and it accommodates different numbers of points at different time steps.

Our third motivation is the following: The degree-Rips bifiltration has emerged in recent years as a popular choice of density-sensitive bifiltration on metric data [52, 10, 45]. The construction is appealing in part because it requires no parameter choice, whereas the sublevel-Rips filtration of a (co)density function requires a choice of bandwidth parameter. In the same way, one can define parameter-free degree-offset and degree-Čech filtrations. We are interested in the problem of computing an analogous degree-Delaunay filtration which is topologically equivalent to the degree-offset bifiltration and more computable than degree-Rips or degree-Čech for low-dimensional data. In ongoing work on this problem, we noticed that the problems of computing the sublevel Delaunay(-Čech) trifiltration and the degree-Delaunay bifiltrations are very closely related. As computing the former is somewhat simpler but already rather involved, we view the results of this paper as a key step towards the development of computationally efficient degree-Delaunay bifiltrations. However, degree-Delaunay bifiltrations will not be studied in this paper.

Our work and the prior work [3] solve particular instances of a general computational problem: Given a filtered point cloud $Z={(Z_p)}_{p\in P}$ indexed by a poset $P$, one has an associated offset filtration indexed by the product poset $P\times [0,\mathrm{\infty })$; the problem is to define and compute a generalized Delaunay filtration of reasonable size which is weakly equivalent to this offset filtration. The case that $Z$ is of the form $X\to X\cup Y\leftarrow Y$ was studied by Reani and Bobrowski [51], who call the resulting generalized Delaunay filtration the coupled $\alpha$-filtration. A more general case, where instead of considering two point clouds and their union, one considers $n\ge 2$ point clouds and their union, was treated in [27, 7]. The authors call the resulting Delaunay objects chromatic $\alpha$-complexes.

Many works have studied time-varying point cloud or metric data through the lens of persistence. To avoid the difficulties of multiparameter persistence, such work most often either considers persistence with respect to a scale parameter and ignores persistence across the time parameter [49, 1, 59, 35, 48, 22], or fixes a scale parameter and considers extended or zig-zag persistence across time [56, 43, 42]. In particular, the problem of extending Delaunay complexes to time-varying data has been studied in [38, 29] using 1-parameter (extended) persistent homology with a fixed radius parameter. It is natural to simultaneously consider persistence with respect to both scale and time, working in the multiparameter setting. This was previously explored by Kim and Mémoli [41], who introduce and study a Rips trifiltration for time-varying metric data with consistent labels across time, and by Saunders [53] and Dey and Samaga [26], who study a version of multiparameter persistence indexed over a product of a zigzag poset and a totally ordered set.

Finally, we note that a construction similar to the one used in [3] was previously used by Sheehy to construct sparse Delaunay filtrations [54].

Section 2 covers preliminaries that are needed throughout the paper. In Section 3, we define the incremental Delaunay complex $ℐ$ of a function $\gamma :X\to {ℝ}^2$, thereby completing the definition of the trifiltrations $𝒟(\gamma )$ and $𝒟𝒞(\gamma )$. Section 4 introduces conflict pairs and presents our characterization of $ℐ$ in terms of these. In Section 5, we apply this characterization to obtain an algorithm for computing $ℐ$. We also discuss optimizations to this algorithm. In Section 6, we establish that $𝒟(\gamma )$ and $𝒟𝒞(\gamma )$ are weakly equivalent to the sublevel offset trifiltration $𝒪(\gamma )$. Finally, in Section 7, we describe our implementation and report the results of our experiments.

To arrive at a contradiction, suppose , and $(\sigma ,y)$ is a conflict pair at $p$. Then . Thus, by injectivity of ${\gamma }_1$ and ${\gamma }_2$, $\gamma (x)\le p$. Since $(\sigma ,x)$ is a conflict pair, we have $\sigma \notin \mathrm{Del}(X_p)$, which contradicts that $(\sigma ,y)$ is a conflict pair at $p$. $◀$

Lemma 7 implies that for any $d$-simplex $\sigma$ in ${\bigcup }_p\mathrm{Del}(X_p)$, the set

is totally ordered by ${\gamma }_1$-coordinate. Henceforth, we assume that $T_{\sigma }$ is given this order.

If $x,y$ are ordered consecutively, then it is easily checked that $(\sigma ,x,y)$ is a conflict triple at $p$, where $p\nearrow =\gamma (x)\vee \gamma (y)$. Conversely, if there exists $z\in T_{\sigma }$ with in the order on $T_{\sigma }$, then for any $p$ with $x\in X_{p\uparrow }$ and $y\in X_{p\to }$, we have $z\in X_p$, so $\sigma \notin \mathrm{Del}(X_p)$. Hence $(\sigma ,x,y)$ is not a conflict triple. $◀$

The next lemma follows immediately from the definition of $ℐ$.

In a way that the following proposition makes precise, every simplex is a face of a $(d+2)$-simplex or $(d+1)$-simplex coming from conflicts as above. For this, we must assume that $X$ satisfies the $\mathrm{\Delta }$-property, i.e., the total orders and coincide in the first $d+1$ points and $X$ is contained in the convex hull of these first $d+1$ points. This technical condition can be easily enforced in practice by simply adding $d+1$ artificial points “at infinity.” The analogous condition for a single total order (or a slight variant thereof) is commonly exploited in the incremental computation of Delaunay triangulations [31, 11, 12], and is also used in the previous work on sublevel-Delaunay bifiltrations [3].

We apply the above results, together with a result about the incremental construction of Delaunay triangulations [31], to bound the size of $ℐ$:

We make essential use of the Bowyer-Watson algorithm [15, 57], a standard incremental algorithm for computing Delaunay triangulations. Given the Delaunay triangulation $\mathrm{Del}(Y)$ of a point cloud $Y\subset {ℝ}^d$ and a point $y\in {ℝ}^d$, the Bowyer-Watson algorithm computes $\mathrm{Del}(Y\cup \{y\})$. To do so, it first identifies each $d$-simplex $\sigma$ such that $y$ lies inside the circumsphere of $\sigma$. Each such simplex $\sigma$ is then removed from $\mathrm{Del}(Y)$, resulting in a star-convex region $R$ centered at $y$. Finally, to obtain $\mathrm{Del}(Y\cup \{y\})$, the region $R$ is retriangulated by forming the simplicial cone with base the boundary of $R$ and apex $y$ (See Figure 4). For $\sigma$ a $d$-simplex as above, we call the pair $(\sigma ,y)$ a BW-conflict (for $Y$).

To compute all conflict pairs, a natural strategy is to traverse each horizontal and each vertical line in $\mathrm{Grid}(\gamma )$ in increasing order, iteratively applying the Bowyer-Watson algorithm. To explain this, we focus on the traversal of vertical lines, as the case of horizontal lines is symmetric. A naive version of the strategy, which already satisfies our main complexity bound (Theorem 13), is as follows: For each $a\in \mathrm{im}{\gamma }_1$, iterate through $\mathrm{im}{\gamma }_2=\{b_1,\mathrm{\dots },b_n\}$ in increasing order. For each , let $p=(a,b_i)$. Given $\mathrm{Del}(X_p)$, if $X_{p\uparrow }\setminus X_p$ is non-empty, hence equal to $\{x\}$ for some $x\in X$, then we apply Bowyer-Watson to compute both $\mathrm{Del}(X_{p\uparrow })$ and all conflict pairs of the form $(\sigma ,x)$ at $p$.

Once $ℐ$ is computed, completing the computation of either $𝒟(\gamma )$ or $𝒟𝒞(\gamma )$ amounts to computing the birth index of each simplex $\sigma \in ℐ$ in the trifiltration. Recall that the birth index of $\sigma$ in $𝒟𝒞(\gamma )$ is $(\gamma (\sigma ),m_{\sigma })$, where $m_{\sigma }$ is the radius of the smallest enclosing ball of $\sigma$. Computation of the first two components of this triple is trivial. In general, smallest enclosing balls can be computed via standard algorithms with efficient implementations [58, 34, 33, 36, 32]. Assuming $d$ is constant, computing $m_{\sigma }$ also requires constant time. Further, recall that the birth index of $\sigma$ in $𝒟(\gamma )$ is $(\gamma (\sigma ),{\omega }_{\sigma })$, where ${\omega }_{\sigma }$ is the incremental Delaunay radius of Definition 3. We discuss the computation of ${\omega }_{\sigma }$ in the full version, where we observe that ${\omega }_{\sigma }$ can be computed for all $\sigma \in ℐ$ in amortized constant time per simplex. From this and Theorem 13, we obtain the following:

A cover of a set $X\subset {ℝ}^d$ is a set $C={\{C(i)\}}_{i\in I}$ of subsets of $X$ whose union is $X$. The nerve of $C$ is the simplicial complex

We say $C$ is a good cover if every intersection of cover elements is either empty or contractible.

These definitions extend naturally to $𝒫$-indexed filtrations, as follows: For $F$ a $𝒫$-indexed filtration, a cover of $F$ is a collection of $𝒫$-indexed filtrations $𝒞={\{𝒞(i)\}}_{i\in I}$, such that for each $p\in 𝒫$, ${\{𝒞{(i)}_p\}}_{i\in I}$ is a cover of $F_p$. We say that $𝒞$ is good if for each $p$, ${\{𝒞{(i)}_p\}}_{i\in I}$ is a good cover of $F_p$. We define $\mathrm{Nrv}(𝒞)$, the nerve of $𝒞$, to be the $𝒫$-indexed filtration given by $\mathrm{Nrv}{(𝒞)}_p=\mathrm{Nrv}({\{𝒞{(i)}_p\}}_{i\in I})$, with structure maps the inclusions.

We will use the following functorial version of the nerve theorem for semi-algebraic sets; see also [6] for a thorough treatment of other variants of the nerve theorem.

This follows immediately from [6, Theorem 5.9], using the fact that an inclusion of closed semi-algebraic sets satisfies the homotopy extension property [24, Theorem 4]. $◀$

For $(p,r)\in {ℝ}^2\times {ℝ}^{+}$, the telescopic offset ${𝒯𝒪}_{p,r}\subset {ℝ}^d\times {ℝ}^2$ is given by

Equivalently,

where $[\gamma (x),p]=\{q\in {ℝ}^2\mid \gamma (x)\le q\le p\}$. Varying $p$ and $r$, we obtain the telescopic offset trifiltration $𝒯𝒪:{ℝ}^2\times {ℝ}^{+}\to \text{Top}$. Let $\psi :𝒯𝒪\to 𝒪(\gamma )$ be the natural transformation induced by the projection of ${ℝ}^d\times {ℝ}^2$ onto the first $d$ coordinates.

We next define our cover of $𝒯𝒪$. For $x\in X$, we define the telescopic Voronoi cell

where $\mathrm{cl}$ denotes the closure. For $(p,r)\in {ℝ}^2\times {ℝ}^{+}$, we define the telescopic Voronoi ball

See Figure 6 for an illustration in the case $d=1$. Allowing $p$ and $r$ to vary, the spaces $𝒯{(x)}_{p,r}$ assemble into a trifiltration $𝒯(x):{ℝ}^2\times {ℝ}^{+}\to \text{Top}$. We let $𝒯={\{𝒯(x)\}}_{x\in X}$.

In what follows, we identify simplices of $\mathrm{Nrv}(𝒯)$ with their corresponding subsets of $X$.

The test suite we use, obtained from [39], consists of point clouds obtained by sampling 1-spheres ($S^1$) and unit squares (${[0,1]}^2$) in ${ℝ}^2$, and 2-spheres ($S^2$), tori ($S^1\times S^1$) and unit cubes (${[0,1]}^3$) in ${ℝ}^3$, with 5% noise drawn from a uniform distribution and perturbation to ensure general position. We consider point clouds of 500, 1000, 2000, 4000, 8000, 16000 points, with four point clouds per type and size.

For each point cloud $X$ we computed the interlevel Delaunay-Čech trifiltration of four different functions $\delta :X\to ℝ$, using the local algorithm. The four functions are:

• $\mathrm{■}$

  codensity, $\delta (x)=-\sum _{y\ne x}\exp \left(-\frac{{\Vert x-y\Vert }^2}{{\sigma }^2}\right)$, with $\sigma$ chosen as the ${0.1}^{\mathrm{th}}$ percentile of the non-zero distances between points in $X$,

• $\mathrm{■}$

  $L_1$-coeccentricity, $\delta (x)=-\sum _{y\in X}\frac{\Vert x-y\Vert }{|X|}$,

• $\mathrm{■}$

  height, $\delta (x)=x_d$, where $x_d$ is the last coordinate of $x$,

• $\mathrm{■}$

  random, where for each $x\in X$, $\delta (x)$ is chosen uniformly at random from $[0,1]$.

All experiments were performed on a computer with an Intel Core i7-5960X CPU @3.00GHz and 64GB of memory, running Ubuntu 20.04.6 LTS. The code was compiled with g++ 9.4.0.

The full results of the experiments are available on Zenodo²²2https://doi.org/10.5281/zenodo.19227801. A representative subset of the results is given in Table 1; this data is plotted in Figure 7. We observe that across all examples, the size of the incremental Delaunay complex and the memory usage grow nearly linearly as a function of the input size. In contrast, the runtime grows nearly quadratically. Thus, our experiments indicate that our approach is memory efficient, and therefore that substantially larger computations should be feasible.

Table 2 compares the performance of the local and non-local algorithms on three types of examples, with $|X|$ up to 4000. While the local algorithm offers no improvement over the non-local algorithm in the first example type (interlevel trifiltrations), for the two other example types, the local algorithm is always faster by a factor of between 3 and 14, and the factor increases with the size of the data set.