Improved Topological Approximations by Digitization

TL;DR

This paper introduces an efficient approximation scheme for bech complexes that significantly reduces computational complexity for large Euclidean datasets by strategic sampling, enabling scalable topological data analysis.

Contribution

It presents a novel bech complex approximation method with size bounds independent of data spread, improving existing schemes in Euclidean space.

Findings

Achieves bech complex approximation with size (n(1/psilon))^{O(d)}

Adds additional sample points to improve approximation accuracy

Bound is independent of the data spread by pre-identifying critical scales

Abstract

\v{C}ech complexes are useful simplicial complexes for computing and analyzing topological features of data that lies in Euclidean space. Unfortunately, computing these complexes becomes prohibitively expensive for large-sized data sets even for medium-to-low dimensional data. We present an approximation scheme for $(1+\epsilon)$-approximating the topological information of the \v{C}ech complexes for $n$ points in $\mathbb{R}^d$, for $\epsilon\in(0,1]$. Our approximation has a total size of $n\left(\frac{1}{\epsilon}\right)^{O(d)}$ for constant dimension $d$, improving all the currently available $(1+\epsilon)$-approximation schemes of simplicial filtrations in Euclidean space. Perhaps counter-intuitively, we arrive at our result by adding additional $n\left(\frac{1}{\epsilon}\right)^{O(d)}$ sample points to the input. We achieve a bound that is independent of the spread of the point set by pre-identifying the scales at which the \v{C}ech complexes changes and sampling accordingly.

Equations

$$n{2}^{O(d\log d+dk)}{\left(\frac{1}{\epsilon }\right)}^{O(d)}$$

$$\mathrm{nerve}(U):=\{V\subseteq U\mid \bigcap _{U_i\in V}{U}_i\ne \varnothing \}.$$

$$V_{\frac{\alpha }{c}}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

$${\mathcal{B}}_{\alpha }:={\left(B(p,\alpha )\right)}_{p\in P},$$

$$I:=\{{\alpha }_k:=(1+\epsilon {)}^k\mid k\in \mathbb{Z}\}.$$

$$|{\mathcal{S}}_{\alpha }|\subseteq {\mathcal{B}}_{(1+\epsilon /2)\alpha }\subseteq |{\mathcal{S}}_{(1+\epsilon )\alpha }|.$$

$$(|{\mathcal{S}}_{\alpha }|{)}_{\alpha >0}.$$

$$\left(I_1^{(0)}\times \dots \times I_d^{(0)}\right),\dots ,\left(I_1^{(k)}\times \dots \times I_d^{(k)}\right)$$

$$\left(I_1^{(0)}\cap \dots \cap I_1^{(k)}\right)\times \dots \times \left(I_d^{(0)}\cap \dots \cap I_d^{(k)}\right),$$

$$({\mathcal{X}}_{\alpha }{)}_{\alpha >0}.$$

$$n{\left(\frac{1}{\epsilon }\right)}^d{2}^{O(d\log d+dk)}.$$

$${\left(\frac{\alpha }{\frac{\epsilon \beta }{8\sqrt{d}}}\right)}^d\le {\left(\frac{16\sqrt{d}}{\epsilon }\right)}^d={\left(\frac{1}{\epsilon }\right)}^d{2}^{O(d\log d)}$$

$$n{\left(\frac{1}{\epsilon }\right)}^d{2}^{O(d\log d)}{2}^{O(dk)}=n{\left(\frac{1}{\epsilon }\right)}^d{2}^{O(d\log d+dk)},$$

$$n{\left(\frac{1}{\epsilon }\right)}^d{2}^{O(d\log d+dk)}.$$

$$diam(A_i),diam(B_i)\le \delta \cdot \underset{p\in A_i,q\in B_i}{\min }\Vert p-q\Vert ,$$

$$R_{(A,B)}:=[d(A,B)/8,8d(A,B)]$$

$$R_p:=\bigcup \{R_{(A,B)}\mid (A,B)\in W\text{ and p=repA or p=repB}\}\subset [0,\infty )$$

$$r_p(\alpha ):=\max \{r\le \alpha \mid r\in R_p\}$$

$${\mathcal{LB}}_{\alpha }:=\bigcup _{p\in P}B(p,r_p(\alpha )).$$

$${\tilde{R}}_{(A,B)}:=[d(A,B)/4,4d(A,B)]$$

$${\tilde{R}}_p:=R_p\cup \bigcup \{{\tilde{R}}_{(A,B)}\mid (A,B)\in W\text{ and p∈A or p∈B}\}\subset [0,\infty ),$$

$${\tilde{r}}_p(\alpha ):=\max \{r\le \alpha \mid r\in {\tilde{R}}_p\},$$

$${\widetilde{\mathcal{LB}}}_{\alpha }:=\bigcup _{p\in P}B(p,{\tilde{r}}_p(\alpha )).$$

$${\tilde{B}}_{\alpha }:=\{B(p,{\tilde{r}}_p(\alpha ))\mid p\in P\}.$$

$$d(A,B)\ge \Vert p_i-p_j\Vert -diam(A)-diam(B)\ge \rho -2\delta d(A,B),$$

$$\rho \le (1+2\delta )d(A,B)\le 2d(A,B)$$

$$\rho \ge \frac{(1-2\delta )}{2}d(A,B)\ge d(A,B)/4$$

$$H({\mathcal{LB}}_{\alpha })\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

$$H({\mathcal{LB}}_{c\alpha })\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}\text{ignorespaces}$$

$$n\log n{2}^{O(d)}+n{\left(\frac{1}{\delta }\right)}^{O(d)}.$$

Full text

1 Introduction

Context and Motivation.

Topological data analysis attempts to extract relevant information out of a data set by interpreting data as a shape and understanding the connectivity of this shape. Connectivity clearly includes the study of connected components, but also the presence of tunnels, voids, and other high-dimensional topological features.

How do we interpret data as a shape? We consider the common scenario that our data is a set of $n$ points $P$ in a possibly high-dimensional Euclidean space $\mathbb{R}^{d}$. Fixing a positive real scale parameter $\alpha$, we define $\mathcal{B}_{\alpha}(P)$ (or simply $\mathcal{B}_{\alpha}$) as the union of balls of radius $\alpha$ centered at each point in $P$. The shape $\mathcal{B}_{\alpha}$ can be seen as an approximation of the shape formed by $P$ for a parameter $\alpha$. As seen in Figure 1, for different values of $\alpha$, the union $\mathcal{B}_{\alpha}$ may form different topological configurations raising the question of the “right” $\alpha$-value that accurately captures the shape of $P$.

One answer, formalized by the theory of persistent homology, is to generate a summary of every topological feature formed at any scale. Note that $\mathcal{B}_{\alpha}\subseteq\mathcal{B}_{\alpha^{\prime}}$ for $\alpha\leq\alpha^{\prime}$. The collection $(\mathcal{B}_{\alpha})_{\alpha\geq 0}$ together with inclusion maps $\mathcal{B}_{\alpha}\to\mathcal{B}_{\alpha^{\prime}}$ whenever $\alpha\leq\alpha^{\prime}$ forms a filtration. The topological properties of this filtration can be summarized by a finite multi-set of points in the plane, called the persistence diagram. For every topological feature, a point in the diagram denotes the range of scales for which the feature is present. This range is called the persistence of the feature.

Importantly, there are several notions of a distance between two diagrams. In this paper, we use the multiplicative bottleneck distance to compare persistence diagrams. Intuitively, diagrams with small distance have similar features that have long persistence; note that this does not imply that their homology is equal on any scale (see Section 2 for details).

Another important aspect of the above theory is the possibility of efficiently computing persistence diagrams. The Čech complex $\mathcal{C}_{\alpha}$ is the intersection/nerve complex of the balls of radius $\alpha$ and has the same homological properties as the union of balls. The derived Čech filtration is a combinatorial object yielding the same persistence diagram as the union of balls. The problem with this standard approach is the sheer size of the produced filtration. It is common practice to only look at the $k$-skeleton of the Čech filtration, that is, ignore intersection of $(k+2)$ or more balls, where $k\leq d$ is another constant parameter. That filtration consists of $O(n^{k+1})$ simplices, which is too large for large data sets, even for relatively small values of $k$. In this paper, we will present a compact and efficiently computable approximation of the Čech filtration.

Problem statement and prior work.

We study the question: can we efficiently compute a discrete filtration $(A_{\alpha_{i}})_{i=1,\ldots,\ell}$, called the approximate filtration, that is significantly smaller than the Čech filtration with the persistence diagrams of Čech and approximate filtration being close to each other? Note that the question defines three quality criteria: the size, computation time and the quality of the approximation.

Sheehy [27] showed that an $\varepsilon$-approximate filtration of size $n(\frac{1}{\varepsilon})^{O(dk)}$ can be computed efficiently for the $k$-skeleton. For $k=d$, the size of this filtration is $n(1/\varepsilon)^{O(d^{2})}$. His result, formulated for the closely related Vietoris-Rips complex, also extends to the Čech complex. Note that the size is linear in the number of points, which is a significant improvement over $O(n^{k+1})$ especially when $n$ is large and $d$, $k$, and $1/\varepsilon$ are of moderate size. The technical idea of the approximation scheme is to compute a net-tree based hierarchical clustering and using it to combine balls whose centers are in the same cluster into a single ball. As a result of this, a large number of ball intersections can be avoided. Similar ideas with the same guarantees, have appeared in [5, 16, 24]. The dimension $d$ in the complexity bounds can also be replaced with the doubling dimension of the metric space, which makes the result especially useful for spaces with low doubling dimension.

Choudhary et al. [15] presented an alternative approximation scheme, where the approximation quality is only $O(d)$, but the size of the approximation is significantly smaller at $n2^{O(d\log k)}$. This result is obtained by tiling the space into permutahedra (a generalization of the hexagonal grid in $\mathbb{R}^{2}$) where the diameter of the permutahedra is controlled by the scale parameter $\alpha$. Every permutahedron containing at least one input point is selected and the nerve of every selected permutahedron is reported as the approximation. The improvement in size stems from the fact that at most $(d+1)$-permutahedra can intersect in $d$ dimensions which upper bounds the number of simplices in the nerve. The approximation quality has been improved to $O(\sqrt[4]{d})$ in subsequent work, with size $n2^{O(d\log k)}$ [14]. In [15], it is also shown that, for any $\varepsilon<1/\log^{1+c}n$ with $c\in(0,1)$, any $\varepsilon$-approximate filtration must have a size of at least $n^{\Omega(\log\log n)}$.

Our contribution.

We derive an approximation scheme of size

[TABLE]

for an $\varepsilon$-approximation of the Čech complex for $\varepsilon\in(0,1]$. For constant dimension $d$, this simplifies to $n\left(\frac{1}{\varepsilon}\right)^{O(d)}$, improving all previous results for the Euclidean case. We achieve this approximation based on the techniques devised in [14, 15]. On a fixed scale $\alpha$, we tile the space with a cubical grid, carefully select a subset of them and take their nerve as our approximation complex. The novelty lies in the resolution of the tiling. Unlike in [15] where the diameter of a permutahedron was in the order of $\alpha$, we use a much smaller diameter of roughly $\varepsilon\alpha$, and we select a hypercube in the approximation if its center is $\alpha$-close to an input point. This is equivalent to approximating the union of $\alpha$-balls $\mathcal{B}_{\alpha}$ with cubical “pixels” at resolution $\varepsilon$ (see Figure 4). The number of selected pixels can be bounded by $n\left(\frac{d}{\varepsilon}\right)^{d}$ which is significantly larger than the number of points. Perhaps counter-intuitively, its nerve is still smaller than in previous approaches, because each vertex of the nerve is only incident to $2^{O(dk)}$ simplices in the $k$-skeleton. This technique resembles previous work in computational geometry where adding points (referred to Steiner points) helps to reduce the size of triangulations [1]. Hudson et al. [22] used similar techniques to compute an $\varepsilon$-approximation of the Čech complex of size $n/\varepsilon^{O(d^{2})}$.

Our results pave the way for important extensions whose proofs we omit for the sake of brevity of the presentation. We only announce them to underline the importance of our techniques and postpone the technical discussion to an extended version of this manuscript:

• •

Our approximation from above is not a filtration, but a simplicial tower, that is a sequence of simplicial complexes connected by simplicial maps. Our technique allows, however, for a definition of an actual approximate filtration of the same size. This filtration has the additional property of being a flag complex at each scale, which has computational benefits as its persistence diagram can be computed using faster algorithms.111e.g., Ripser by U.Bauer: https://github.com/Ripser/ripser

• •

The factor $2^{O(d\log d+dk)}$ in our size bound can be further reduced to $2^{O(d\log d)}$ by replacing the cubical grid by a permutahedral grid [15]. A consequence of this result is that when $d=\Theta(\log n)$ and $\varepsilon=1/\log^{1+c}n$, for some $c\in(0,1)$, our approximation with the permutahedral grid has a size of $n^{O(\log\log n)}$ matching the lower bound of [15].

The result of this work are proved in two steps. We first devise a simple approximation scheme which approximates the union of balls by pixels of a carefully chosen size on each scale. It follows with standard techniques that the approximation quality is $(1+\varepsilon)$. However, the total size of the approximation has an additional factor that depends on the spread of the point set. To achieve spread-independence, we use the observation that when we increase the scale from $\alpha$ to $(1+\varepsilon)\alpha$, the growth in radius of a ball does only affect the topology of the union of balls if the enlarged balls gives rise to a new intersection pattern with other balls. In other words, we can delay the growth of a ball (and hence, avoiding resampling the balls by pixels) if $\alpha$ is not in a critical range of scales. Moreover, we can bound the critical range of all balls efficiently through a well-separated pair decomposition (WSPD) of the point set. This suffices to reduce the total number of pixels used in the approximation to $n(\frac{1}{\varepsilon})^{O(d)}$. One technical difficulty arising from this approach is that on a fixed scale, the approximation complex consists of cubes of different sizes.

Outline of the paper

We review the topological concepts needed for our result in Section 2. We present the simple, spread-dependent approximation in Section 3, and extend it to a spread-independent approximation in Section 4. We conclude in Section 5.

2 Preliminaries

We give a short introduction to topological concepts that are essential for our results. For more details, we refer to standard references such as [6, 12, 17, 20, 25].

Simplicial complexes and simplicial maps

A simplicial complex $K$ on a finite set of elements $S$ is a collection of subsets $\{\sigma\subseteq S\}$ called simplices such that each subset $\tau\subset\sigma$ is also in $K$. In other words, simplicial complexes are closed under taking subsets. We say that a simplex $\sigma\in K$ has dimension $k:=|\sigma|-1$, in which case $\sigma$ is called a $k$-simplex. A simplex $\tau$ is a face of $\sigma$ if $\tau\subseteq\sigma$. The $k$-skeleton of $K$ consists of all simplices of $K$ of dimension $k$ or lower. For instance, the $1$-skeleton of $K$ is a graph defined by its [math]-simplices and $1$-simplices.

For a point set $P\subset\mathbb{R}^{d}$ and a real number $\alpha\geq 0$, the Čech complex $\mathcal{C}_{\alpha}$ on $P$ at scale $\alpha$ is the collection of all simplices $\sigma\subseteq P$ such that there is a non-empty intersection between the Euclidean balls of radius $\alpha$ centered at the points of $\sigma$. Equivalently, $\sigma\in\mathcal{C}_{\alpha}$ if the minimum enclosing ball of $\sigma$ has radius at most $\alpha$. We denote the union of $\alpha$-balls as $\mathcal{B}_{\alpha}$. See Figure 2 for an example.

A simplicial complex is called a flag complex, if it has the following property: whenever a set $\{p_{0},\ldots,p_{k}\}\subseteq P$ has the property that every $1$-simplex $\{p_{i},p_{j}\}$ is in the complex, then the $k$-simplex $\{p_{0},\ldots,p_{k}\}$ is also in the complex.

A simplicial complex $K^{\prime}$ is a subcomplex of $K$ if $K^{\prime}\subseteq K$. For instance, $\mathcal{C}_{\alpha}$ is a subcomplex of $\mathcal{C}_{\alpha^{\prime}}$ for $0\leq\alpha\leq\alpha^{\prime}$. Let $L$ be another simplicial complex. Let $\hat{\varphi}$ be any map that assigns to each vertex of $K$ a vertex of $L$. A simplicial map is a map $\varphi:K\rightarrow L$ defined using some vertex map $\hat{\varphi}$, such that for every simplex $\sigma=\{p_{0},\ldots,p_{k}\}$ in $K$, the set of vertices $\varphi(\sigma):=\{\hat{\varphi}(p_{0}),\ldots,\hat{\varphi}(p_{k})\}$ is a simplex of $L$. For $K^{\prime}\subseteq K$, the inclusion map $\mathrm{inc}:K^{\prime}\rightarrow K$ is an elementary example of a simplicial map. Note that any simplicial map is completely specified by its action on the vertices of the domain.

Homotopy and Nerve theorem

Let $f_{1},f_{2}:X\rightarrow Y$ be two continuous maps between topological spaces $X,Y$. A continuous function $H:X\times[0,1]\rightarrow Y$ is said to be a homotopy between $f_{1}$ and $f_{2}$, if $H(x,0)=f_{1}(x)$ and $H(x,1)=f_{2}(x)$. In this case, $f_{1}$ and $f_{2}$ are said to be homotopic to each other, and we record this relation as $f_{1}\overset{h}{\simeq}f_{2}$. Informally, the second parameter of $H$ can be interpreted as time, so that $H$ is a continuous deformation of $f_{1}$ into $f_{2}$, as time varies from 0 to 1.

$X$ and $Y$ are said to be homotopy equivalent if there exist continuous maps $f:X\rightarrow Y$ and $g:Y\rightarrow X$ such that $f\circ g\overset{h}{\simeq}\mathrm{id}_{Y}$ and $g\circ f\overset{h}{\simeq}\mathrm{id}_{X}$, where $\mathrm{id}_{X}$ and $\mathrm{id}_{Y}$ are the identity functions on $X$ and $Y$, respectively. We record this relation as $X\overset{h}{\simeq}Y$. Intuitively, two spaces are homotopy equivalent if they can be continuously transformed into one another.

Let $U:=\{U_{1},\ldots,U_{m}\}$ denote a finite collection of sets. The nerve of $U$ is an abstract simplicial complex $\mathrm{nerve}(U)$ consisting of simplices corresponding to non-empty intersections of elements of $U$, that is,

[TABLE]

For instance, the Čech complex at scale $\alpha$ is $\mathcal{C}_{\alpha}=\mathrm{nerve}(\mathcal{B}_{\alpha})$.

Theorem 2.1** (nerve theorem).**

Let us denote by $U:=\{U_{1},\ldots,U_{m}\}$ a finite collection of closed sets in Euclidean space, such that all non-empty intersections of the form $\{\cap_{U_{i}\in V}U_{i}\neq\varnothing\mid V\subseteq U\}$ are contractible. Then, $\mathrm{nerve}(U)\overset{h}{\simeq}\bigcup_{i=1}^{m}U_{i}$, that is, they are homotopy equivalent (see [2, 3, 4, 28] for more general versions of the theorem).

In particular, if $U$ is a collection of convex sets, then all non-empty common intersections are contractible, so that the nerve theorem applies to $U$. As an example, $\mathcal{C}_{\alpha}\overset{h}{\simeq}\mathcal{B}_{\alpha}$.

Filtrations and simplicial towers

Let $I\subseteq\mathbb{R}$ be a set of real values which we refer to as scales. A filtration is a collection of simplicial complexes indexed by $I$, $(K_{\alpha})_{\alpha\in I}$ such that $K_{\alpha}\subseteq K_{\alpha}^{\prime}$ for all $\alpha\leq\alpha^{\prime}\in I$. For instance, $(\mathcal{C}_{\alpha})_{\alpha\geq 0}$ is the Čech filtration. A (simplicial) tower is a sequence $(K_{\alpha})_{\alpha\in J}$ of simplicial complexes with $J$ being a discrete set (for instance $J=\{2^{k}\mid k\in\mathbb{Z}\}$), together with simplicial maps $\varphi_{\alpha}:K_{\alpha}\rightarrow K_{\alpha^{\prime}}$ between complexes at consecutive scales.

A tower of the form $(K_{\alpha})_{\alpha\in J}$ on simplicial maps $g^{.,.}$ with index set $J=\{0\leq j_{1}<j_{2}<\ldots\}$ can always extend the tower to $(K_{\alpha})_{\alpha\geq 0}$, by setting $K_{\alpha}=K_{j_{i}}$ and $g^{\alpha,j_{i+1}}=g^{j_{i},j_{i+1}}$ for all $\alpha\in[j_{i},j_{i+1})$. We call this technique of extending towers the standard filling technique in our paper for simplicity. The approximation constructed in this paper will be an example of such a tower.

We say that a simplex $\sigma$ is included in the tower at scale $\alpha^{\prime}$ if $\sigma\in K_{\alpha^{\prime}}$ is not in the image of the map $\varphi_{\alpha}:K_{\alpha}\rightarrow K_{\alpha^{\prime}}$, where $\alpha$ is the scale preceding $\alpha^{\prime}$ in the tower. The size of a tower is the number of simplices included over all scales. If a tower arises from a filtration, its size is simply the size of the largest complex in the filtration (or infinite, if no such complex exists). However, this is not true in general for simplicial towers, since simplices can collapse in the tower and the size of the complex at a given scale may not take into account the collapsed simplices which were included at earlier scales in the tower.

Persistence diagram and Interleavings

A collection of vector spaces $(V_{\alpha})_{\alpha\in I}$ connected with homomorphisms (linear maps) $\lambda_{\alpha_{1},\alpha_{2}}:V_{\alpha_{1}}\rightarrow V_{\alpha_{2}}$ is called a persistence module, if $\lambda_{\alpha,\alpha}$ is identity on $V_{\alpha}$ for all $\alpha\in I$, and $\lambda_{\alpha_{2},\alpha_{3}}\circ\lambda_{\alpha_{1},\alpha_{2}}=\lambda_{\alpha_{1},\alpha_{3}}$ for all $\alpha_{1}\leq\alpha_{2}\leq\alpha_{3}\in I$ for the index set $I$.

We can generate persistence modules using simplicial complexes. Given a simplicial tower $(K_{\alpha})_{\alpha\in I}$, we can fix some base field $\mathcal{F}$ to obtain a sequence $(H(K_{\alpha}))_{\alpha\in I}$ of vector spaces with linear maps $\overline{\varphi}^{\ast}$, that forms a persistence module. This is true because of the functorial properties of homology [25]. The same construction is also applicable to filtrations.

Under certain tameness conditions, persistence modules admit a decomposition into a collection of intervals of the form $[\alpha,\beta]$ (with $\alpha,\beta\in I$). These intervals can be represented as a set of points in the plane, called the persistence diagram, where each interval $[\alpha,\beta]$ is simply represented as the point $(\alpha,\beta)$. See Figure 3 for an example. The persistence diagram of a persistence module characterizes the module uniquely up to isomorphism. If the persistence module is generated by a simplicial complex, each point $(\alpha,\beta)$ in the diagram corresponds to a homological feature (a “hole”) that comes into existence at complex $K_{\alpha}$ and persists until it disappears at $K_{\beta}$.

Two persistence modules $(V_{\alpha})_{\alpha\geq 0}$ and $(W_{\alpha})_{\alpha\geq 0}$ with respective linear maps $\phi_{\cdot,\cdot}$ and $\psi_{\cdot,\cdot}$ are said to be (multiplicatively) strongly $c$-interleaved if there exist a pair of families of linear maps $\gamma_{\alpha}:V_{\alpha}\rightarrow W_{c\alpha}$ and $\delta_{\alpha}:W_{\alpha}\rightarrow V_{c\alpha}$ for $c>0$, such that Diagram (1) commutes for all $0\leq\alpha\leq\alpha^{\prime}$ (the subscripts of the maps are excluded for readability). In such a case, the persistence diagrams of the two modules are said to be $c$-approximations of each other in the sense of [10]. More precisely, there is a partial matching between the points of the two diagrams, after a suitable logarithmic scaling of the diagrams. Further details can be found in [7]. See Figure 3 for an example.

[TABLE]

Next, we discuss a special case that relates to the equivalence of persistence modules [9, 18]. Two persistence modules $\mathbb{V}=(V_{\alpha})_{\alpha\in I}$ and $\mathbb{W}=(W_{\alpha})_{\alpha\in I}$ on linear maps $\phi,\psi$ respectively are isomorphic, if there exists an isomorphism $f_{\alpha}:V_{\alpha}\rightarrow W_{\alpha}$ for each $\alpha\in I$ such that the diagram

[TABLE]

commutes for all $\alpha\leq\beta\in I$. Isomorphic persistence modules have identical persistence barcodes.

3 A simple digitization scheme

In this section we describe our first approximation scheme which is conceptually simple and lays down the foundation for the more technical approximation scheme that we present in Section 4. We describe our scheme for a set of $n$ points $P$ in $d$-dimensional Euclidean space. We assume throughout that $\varepsilon\in(0,\frac{1}{5}]$.

A cubical subdivision

Let $J:=\{2^{i}\mid i\in\mathbb{Z}\}$. For any $\alpha\in J$, we define the lattice $L_{\alpha}:=\left(\frac{\varepsilon\alpha}{4\sqrt{d}}\right)\mathbb{Z}^{d}$ as the $\mathbb{Z}^{d}$ lattice whose basis vectors have been scaled by $\frac{\varepsilon\alpha}{4\sqrt{d}}$. We define the pixels of this lattice as the smallest $d$-dimensional cubes whose vertices are the lattice points. The diameter of the pixels of $L_{\alpha}$ is at most $\varepsilon\alpha/4$. Our approximation complexes are built as nerves of the union of pixels of $L_{\alpha}$.222It seems simpler at first to define the approximation as a cubical complex directly instead of taking the nerve, but it is more complicated to construct the chain maps connecting different scales in this cubical setup.

Each pixel of lattice $L_{\alpha}$ is fully contained in some pixel of $L_{\beta}$ when $\beta>\alpha$ with $\beta,\alpha$ both in $J$. This also implies that each pixel center of $L_{\alpha}$ has a unique nearest-neighbor among the pixel centers in $L_{\beta}$.

3.1 Approximation complex

Let $B(p,\alpha)$ denote the $d$-dimensional Euclidean ball of radius $\alpha$, centered at any input point $p\in P$. We denote by

[TABLE]

the set of $\alpha$-balls centered at the input points. Naturally, the Čech complex on $P$ at scale $\alpha$ is the nerve of $\mathcal{B}_{\alpha}$, that is, $\mathcal{C}_{\alpha}=\mathrm{nerve}(\mathcal{B}_{\alpha})$.

Let $I$ denote the set of scales

[TABLE]

We now define our approximation complex $\mathcal{X}_{\alpha}$ for $\alpha\in I$. Consider the maximal $\beta\in J$ such that $\beta\leq\alpha$. As discussed, there exists a cubical subdivision based on the lattice $L_{\beta}=(\frac{\varepsilon\beta}{4\sqrt{d}})\mathbb{Z}^{d}$. Let $\mathcal{S}_{\alpha}$ denote the set of pixels of that subdivision whose center lies in $\mathcal{B}_{\alpha}$, that is, which are in distance at most $\alpha$ to a point in $P$. We define $\mathcal{X}_{\alpha}:=\mathrm{nerve}(\mathcal{S}_{\alpha})$. Also, we write $|\mathcal{S}_{\alpha}|$ for the union of all pixels in $\mathcal{S}_{\alpha}$. See Figure 4 for an illustration of $\mathcal{B}_{\alpha}$ and $|\mathcal{S}_{\alpha}|$.

Lemma 3.1** (Sandwich lemma).**

For any $\alpha\in I$,

[TABLE]

Proof.

For $\alpha\in I$, let $\beta\in J$ be as in the definition of $|\mathcal{S}_{\alpha}|$. For the first inclusion, if $x\in|\mathcal{S}_{\alpha}|$, $x$ lies in some pixel $C$ of the grid $L_{\beta}$. Let $c$ denote the center of this pixel. There is some $p\in P$ with $\|c-p\|\leq\alpha$. The diameter of $C$ is $\varepsilon\beta/4<\varepsilon\alpha$, hence $\|c-x\|\leq\varepsilon\alpha/2$, and by the triangle inequality, $\|x-p\|\leq(1+\varepsilon/2)\alpha$.

For the second inclusion, fix $x\in\mathcal{B}_{(1+\varepsilon/2)\alpha}$ and let $p\in P$ such that $\|x-p\|\leq(1+\varepsilon/2)\alpha$. Note that $(1+\varepsilon)\alpha\in I$ as well, and the approximation complex is either constructed using the same $\beta$ as for $\alpha$, or using $2\beta$ if $2\beta=(1+\varepsilon)\alpha$. In both cases, $x$ lies in a pixel $C$ with diameter at most $\varepsilon\beta/2<\varepsilon\alpha$, so the distance of $x$ to the center $c$ of $C$ is at most $\varepsilon\alpha/2$. By the triangle inequality, $\|x-p\|\leq(1+\varepsilon)\alpha$, so the pixel belongs to $|\mathcal{S}_{(1+\varepsilon)\alpha}|$. ∎

The lemma shows that $|\mathcal{S}_{\alpha}|\subset|\mathcal{S}_{(1+\varepsilon)\alpha}|$. Using the standard filling technique from Section 2, we extend this discrete space to a continuous filtration:

[TABLE]

Moreover, the Sandwich lemma together with the strong-interleaving diagram (Diagram (1)) implies at once:

Theorem 3.2**.**

$(H(|\mathcal{S}_{\alpha}|))_{\alpha>0}$* $(1+\varepsilon)^{2}$-approximates the persistence module $(H(\mathcal{B}_{\alpha}))_{\alpha\geq 0}$.*

Note that we obtain $(1+\varepsilon)^{2}$ instead of $(1+\varepsilon)$ because we consider the continuous filtration instead of the discrete filtration of $|\mathcal{S}_{\alpha}|$.

The theorem also implies that $(H(|\mathcal{S}_{\alpha}|))_{\alpha>0}$ is a $(1+\varepsilon)^{2}$-approximation of the Čech filtration since the Čech filtration is dual to $(\mathcal{B}_{\alpha})_{\alpha\geq 0}$ and has the same persistence diagram.

3.2 Connecting the scales

We now turn our attention to the approximation complex $\mathcal{X}_{\alpha}$. While at each scale $\alpha$, $\mathcal{X}_{\alpha}$ is the nerve of $\mathcal{S}_{\alpha}$ and $|\mathcal{S}_{\beta}|$ forms a filtration, it is not true that $\mathcal{X}_{\alpha}\subseteq\mathcal{X}_{\beta}$ for all $\alpha\leq\beta$ since $\mathcal{S}_{\alpha}\not\subset\mathcal{S}_{\beta}$. Therefore, it is not sufficient to apply the persistent nerve lemma of [13] directly, but it requires a more involved analysis, which we describe next. We show that the complexes are connected by simplicial maps.

A first useful property is that $\mathcal{X}_{\alpha}$ is a flag complex. This follows from the following statement, which we prove in more general form for later use. We define an axis-aligned cuboid in $\mathbb{R}^{d}$ as the Cartesian product of $d$ intervals $I_{1}\times\ldots\times I_{d}$ (where the degenerate case is allowed that $I_{j}$ consists of only one point). For instance, all pixels of any lattice $L_{\beta}$ are (non-degenerate) cuboids.

Lemma 3.3**.**

The nerve complex of any finite collection of cuboids is a flag complex.

Proof.

We show that if there is a set of $(k+1)$ cuboids that pairwise intersect, then all of them have a common intersection. The intersection of the set of $(k+1)$ cuboids

[TABLE]

is the cuboid

[TABLE]

and it suffices to show that these intersection are non-empty coordinate-wise. By assumption, $I^{(\ell_{1})}_{1}\cap I^{(\ell_{2})}_{1}$ is non-empty for all $0\leq\ell_{1},\ell_{2}\leq k$. Helly’s theorem [21] for the case of intervals implies that all $I^{(\cdot)}_{1}$ intersect commonly. ∎

Next, we consider $\alpha\in I$ and define a simplicial map $g:\mathcal{X}_{\alpha}\to\mathcal{X}_{(1+\varepsilon)\alpha}$. This is simple if $\mathcal{X}_{\alpha}$ and $\mathcal{X}_{(1+\varepsilon)\alpha}$ are constructed using the same grid, as in that case $\mathcal{S}_{\alpha}\subseteq\mathcal{S}_{(1+\varepsilon)\alpha}$ and consequently, $\mathcal{X}_{\alpha}\subseteq\mathcal{X}_{(1+\varepsilon)\alpha}$. If the grid changes, there is no direct inclusion. There is, however, a natural map $g^{\prime}$ mapping each pixel $C$ in $\mathcal{S}_{\alpha}$ to the unique pixel $g^{\prime}(C)$ in $\mathcal{S}_{(1+\varepsilon)\alpha}$ that contains $C$. That fact that $g^{\prime}(C)$ is indeed in $\mathcal{S}_{(1+\varepsilon)\alpha}$ follows from Sandwich Lemma.

Lemma 3.4**.**

The map $g^{\prime}$ extends to a simplicial map $g:\mathcal{X}_{\alpha}\rightarrow\mathcal{X}_{(1+\varepsilon)\alpha}$.

Proof.

From Lemma 3.3 we know that $\mathcal{X}$ is a flag complex. Therefore it suffices to show that $g$ maps every edge of $\mathcal{X}_{\alpha}$ to either a single vertex or an edge of $\mathcal{X}_{(1+\varepsilon)\alpha}$. But that follows at once, because an edge in $\mathcal{X}_{\alpha}$ corresponds to two pixels on the grid of $\mathcal{X}_{\alpha}$ which are intersecting, and the map $g^{\prime}$ as defined above either maps both of them to the same pixel, or to two pixels which are also intersecting. ∎

By composing the maps above, we obtain maps $g^{\alpha_{1},\alpha_{2}}:\mathcal{X}_{\alpha_{1}}\rightarrow\mathcal{X}_{\alpha_{2}}$ for all $0<\alpha_{1}\leq\alpha_{2}\in I$. Using the standard filling technique, we define the (continuous) simplicial tower

[TABLE]

Interleaving

We next establish a relationship between the approximation tower $(\mathcal{X}_{\alpha})_{\alpha>0}$ and $(|\mathcal{S}_{\alpha}|)_{\alpha>0}$. More precisely, we will show that both towers yield the same persistence diagram. Using Theorem 3.2, this implies that $(\mathcal{X}_{\alpha})_{\alpha>0}$ is a $(1+\varepsilon)^{2}$-approximation of the Čech complex.

Since $\mathcal{S}_{\alpha}$ consists of convex objects, the nerve theorem (Theorem 2.1) asserts that $\mathcal{S}_{\alpha}\overset{h}{\simeq}\mathcal{X}_{\alpha}$, for each scale $\alpha$, so they have isomorphic homology groups. All that is left to show is that the homology map $g^{\ast}$ induced by the simplicial maps $g$ from above commutes with the isomorphisms from the Nerve theorem. We prove this indirectly by introducing an intermediate filtration that is equivalent to both filtrations. This intermediate filtration considers the union of all pixels $\mathcal{S}_{\beta}$ with $\beta\leq\alpha$. We shift the technical details to Appendix A.1 and just state the final result.

Lemma 3.5**.**

The persistence modules $(H(\mathcal{X}_{\alpha}))_{\alpha>0}$ and $(H(|\mathcal{S}_{\alpha}|))_{\alpha>0}$ are isomorphic.

We conclude with the main result of this section. Setting $\varepsilon^{\prime}=\varepsilon/4$ in our approximation, we see that $(1+\frac{\varepsilon}{4})^{2}<1+\varepsilon$. We conclude that

Theorem 3.6**.**

$(H(\mathcal{X}_{\alpha}))_{\alpha>0}$* and $(H(\mathcal{C}_{\alpha}))_{\alpha\geq 0}$ are $(1+\varepsilon)$-approximations of each other, for $\varepsilon\in(0,\frac{1}{20}]$.*

3.3 Size and Computation

Theorem 3.7**.**

For every $\alpha\geq 0$, the $k$-skeleton of the approximation tower $\mathcal{X}_{\alpha}$ has size

[TABLE]

Proof.

Let $\beta\leq\alpha$ be such that the complex $\mathcal{X}_{\alpha}$ is built using the lattice $L_{\beta}$. Note that $\beta\geq\alpha/2$. The sidelength of a pixel of this lattice is $\frac{\varepsilon\beta}{4\sqrt{d}}$. A ball of half this radius is contained inside the pixel. Since the pixels are interior-disjoint, a simple packing argument shows that each ball in $\mathcal{B}_{\alpha}$ is covered by no more than

[TABLE]

pixels. There are $n$ balls of $\mathcal{B}_{\alpha}$, hence, $\mathcal{S}_{\alpha}$ contains at most $n\left(\frac{1}{\varepsilon}\right)^{d}2^{O(d\log d)}$ pixels. Each pixel is incident to at most $3^{d}-1=2^{O(d)}$ other pixels. Each simplex incident to a pixel has vertices among these $2^{O(d)}$ pixels. Therefore, the size of the $k$-skeleton incident to each pixel is then $2^{O(dk)}$. In total, the $k$-skeleton of $\mathcal{X}_{\alpha}$ has size

[TABLE]

independent of $\alpha$. ∎

To compute the complex at a given scale $\alpha\in I$, we first find the pixels at that scale using a simple flooding algorithm that starts at vertices of $P$. Then we inspect the neighborhood of each pixel to compute the $k$-skeleton incident to that vertex. See Appendix A.2 for more details.

Theorem 3.8**.**

At each scale of $I$, the $k$-skeleton of the approximation tower and the simplicial map can be computed in time

[TABLE]

4 Removing the spread

The size and computation bounds from Section 3.3 only hold for a single considered scale. To get a bound on the total size and complexity, one needs to multiply with the considered number of scales. If the entire filtration is of interest, this number can be upper bounded by the logarithm of the spread of the point set (the ratio of the diameter and the closest distance). In this section, we remove that dependence on the spread.

4.1 A short overview of the approximation scheme

We first informally illustrate the idea behind eliminating the dependence on spread. The main exposition starts from Subsection 4.2, where we start with the details of our construction.

The crucial observation guiding our improved approximation scheme is that while the union of balls $\mathcal{B}_{\alpha}$ changes continuously for all $\alpha\geq 0$, the Čech complex changes only at a finite number of scales. For instance, if an edge $(p,q)$ enters the Čech filtration at $\alpha$, then the intersection of the $\alpha$-balls at $p$ and $q$ captures the edge. But at higher scales, it is pointless to allow these balls to grow further to represent this edge. So we allow the balls to grow individually only in a lazy fashion, only when they form new connections. This ensures that the set of edges of the original Čech filtration are correctly captured. Allowing for some slack in growing the balls, that is, if we allow the balls at $p,q$ to grow in the interval $[\alpha,2\alpha]$, we can capture all simplices incident to $p$ and $q$. Naturally, this lazy union of balls is homotopy-equivalent to the original union. See Figure 5 for an example.

The lazy union is non-homogeneous, so it is non-trivial to pixelize it. The ideal approach is to pixelize each ball according to its radius, to get a non-homogeneous union of pixels. As we show later, this does not pose a problem if the pixels are chosen as in Section 3. To construct the pixelization efficiently, we make use of a well-separated pair decomposition. Each pair of points $(p,q)\in P$ has roughly the same distance as the distance between any points of a well-separated pair that covers $(p,q)$. Hence, the lazy balls at the representatives approximate the lazy balls at the points of the pair, and we use the former for pixelization. The use of an $\varepsilon$-WSPD ensures that there are at most $n(1/\varepsilon)^{O(d)}$ lazy balls. See Figure 5 for an example (further details in Subsection 4.2).

We collect the pixels that intersect the union of balls. While the lazy union and its pixelization do not interleave on a space level as in Lemma 3.1, we show that they still form persistence modules that are closely interleaved (Subsection 4.3). The nerves of these pixelizations are our approximation complexes, and they can be connected to form a simplicial tower: for this, we simply use the map that takes a pixel at a lower scale to the pixel at the higher scale that contains it.

To compute the pixelization, we first construct a WSPD and identify the intervals at which the representatives are growing. At each such scale, we identify only those balls that are expanding, pixelize them, and take the nerve. This constructs our simplex inclusions at that scale. More details follow in Subsection 4.5.

4.2 Lazy union of balls

We review the concept of a well-separated pair decomposition (WSPD) from [8], which plays a crucial role in our algorithm. A $\delta$-WSPD of $P$ consists of pairs of the form $(A_{i},B_{i})\subset P\times P$ that satisfy

• •

each $(A_{i},B_{i})$ is a well-separated pair (WSP), that means,

[TABLE]

• •

and for each pair of points $(p,q)\in P$, there exists a WSP $(A_{j},B_{j})$ in the WSPD such that either $(p\in A_{j},q\in B_{j})$ or $(p\in B_{j},q\in A_{j})$. That means, a WSPD covers each pair of points of $P$.

Let $W$ denote a $\delta$-WSPD on $P$, where $\delta\leq 1/10$. For each WSP $(A,B)\in W$ let $\{P_{A},P_{B}\}\subset P$ denote the points of $P$ in $A$ and $B$, respectively. We pick a point $rep_{A}\in P_{A}$ and call it the representative of $A$ (similarly for $B$). We denote the distance between the representatives as $d(A,B):=\|rep_{A}-rep_{B}\|$. From the WSPD-property, $max\{diam(A),diam(B)\}\leq\delta d(A,B)$. Using the triangle inequality, we see that the distance between any two points of $P_{A}$ and $P_{B}$ lies in the interval $[d(A,B)/2,2d(A,B)]$ since $\delta\leq 1/10$. For reasons that becomes apparent later, we scale the interval by a factor of $4$ and call

[TABLE]

the active interval of the pair $(A,B)$. For every point $p\in P$, we define

[TABLE]

as the active interval of $p$. Each scale $\alpha\in R_{p}$ is an active scale. Next, we define for $p\in P$,

[TABLE]

and set

[TABLE]

We can interpret these definitions as follows: $R_{p}$ specifies a range of scales $\alpha$ for which the $\alpha$-ball of $p$ might encounter new intersections. The function $r_{p}(\alpha)$ is monotonously increasing, $r_{p}(0)=0$ for all $p$, and $r_{p}(\alpha)=\alpha$ if $\alpha$ is an active scale; otherwise, the radius of the ball just remains at the last encountered active scale. $\mathcal{LB}_{\alpha}$ is the union of balls with radii given by the $r_{p}$ functions.

Note that $\mathcal{LB}_{\alpha}\subseteq\mathcal{LB}_{\alpha^{\prime}}$ whenever $\alpha\leq\alpha^{\prime}$. Hence, $(\mathcal{LB}_{\alpha})_{\alpha\geq 0}$ is a filtration.

Lemma 4.1**.**

The persistence module $(H(\mathcal{LB}_{\alpha}))_{\alpha\geq 0}$ is a $(1+8\delta)$-approximation of $(H(\mathcal{B}_{\alpha}))_{\alpha\geq 0}$ (and consequently, also of the Čech filtration).

The (somewhat tedious) proof can be summarized as follows: we define an intermediate filtration of balls where not only the balls of representatives, but of all balls participating in a WSP grow. In the first part, we show that this filtration can be sandwiched with that of $\mathcal{LB}_{\alpha}$. In the second part, we show that the intermediate filtration has the same nerve as the union of $\alpha$-balls.

We now define the intermediate filtration. Let

[TABLE]

be a scaled version of the active interval from above. Define

[TABLE]

[TABLE]

and set

[TABLE]

Intuitively, all the balls centered at points in a WSP grow when the WSP is active (but those of the representatives grow for a slightly longer time).

Lemma 4.2**.**

For all $\alpha\geq 0$, $\mathcal{LB}_{\alpha}\subseteq\widetilde{\mathcal{LB}}_{\alpha}\subseteq\mathcal{LB}_{(1+8\delta)\alpha}$.

Proof.

The first inclusion follows at once from the fact that $r_{p}\leq\tilde{r}_{p}$ for all points of $P$ and all scales. For the second inclusion, let $x$ be a point in $\widetilde{\mathcal{LB}}_{\alpha}$ and let $p$ be such that $x\in B(p,\tilde{r}_{p}(\alpha))$. There is some WSP $(A,B)$ such that $\tilde{r}_{p}(\alpha)$ lies in $R_{(A,B)}$ and $p$ is the representative of $A$ or of $B$, or $\tilde{r}_{p}(\alpha)$ lies in some $\tilde{R}_{(A,B)}$ a WSP, and $p\in A$ or $p\in B$. In the first case, it follows that $\tilde{r}_{p}(\alpha)=r_{p}(\alpha)$ and hence $x\in\mathcal{LB}_{\alpha}\subseteq\mathcal{LB}_{(1+8\delta)\alpha}$. For the rest of the proof, we assume without loss of generality that $\tilde{r}_{p}(\alpha)$ lies in some $\tilde{R}_{(A,B)}$ and $p\in A$.

Let $q\neq p$ be the representative of $A$. Then $R_{(A,B)}$ is a subset of $R_{q}$. By our choice of $R_{(A,B)}$ and $\tilde{R}_{(A,B)}$, we know that $R_{(A,B)}$ contains the value $(1+8\delta)\tilde{r}_{p}(\alpha)$ because $(1+8\delta)<2$ for $\delta<1/10$. On the other hand, $\|p-q\|\leq\delta d(A,B)\leq 8\delta\tilde{r}_{p}(\alpha)$, where the last inequality comes from the fact $\tilde{r}_{p}(\alpha)\in R_{(A,B)}$. By triangle inequality, $\|x-q\|\leq(1+8\delta)\tilde{r}_{p}(\alpha)$, which implies the claim. ∎

The lemma implies that the two filtrations $(1+8\delta)$-approximate each other. Next, we consider

[TABLE]

which is the collection of balls whose union forms $\widetilde{\mathcal{LB}}_{\alpha}$. By the nerve theorem, the nerve of $\widetilde{B}_{\alpha}$ is homotopically equivalent to $\widetilde{\mathcal{LB}}_{\alpha}$, and since the balls have non-decreasing radius, the persistence modules $(H(\widetilde{\mathcal{LB}}_{\alpha}))_{\alpha\geq 0}$ and $(H(\mathrm{nerve}(\widetilde{B}_{\alpha}))_{\alpha\geq 0}$ are persistence-equivalent. The latter, however, is exactly the Čech persistence module, as we show next.

Lemma 4.3**.**

$\mathrm{nerve}(\widetilde{B}_{\alpha})=\mathrm{nerve}(\mathcal{B}_{\alpha})=\mathcal{C}_{\alpha}$* for all scales $\alpha\geq 0$.*

Proof.

Because $\tilde{r}_{p}(\alpha)\leq\alpha$, $\mathrm{nerve}(\widetilde{B}_{\alpha})\subseteq\mathrm{nerve}(\mathcal{B}_{\alpha})$ follows at once. For the other direction, fix any simplex $\sigma=(p_{0},\ldots,p_{k})$ of $\mathcal{C}_{\alpha}$. Let $\rho\leq\alpha$ denote the smallest radius such that the $\rho$-balls around the $p_{i}$ intersect. It suffices to show that $\tilde{r}_{p_{i}}(\alpha)\geq\rho$ for all $i=\{0,\ldots,d\}$. For $p_{i}$ fixed, let $p_{j}$ denote the point of $\sigma$ with maximal distance to $p_{i}$. Clearly, $\|p_{i}-p_{j}\|\geq\rho$. On the other hand, there exists a WSP $(A,B)$ that covers $(p_{i},p_{j})$ giving rise to the interval $\tilde{R}_{(A,B)}=[d(A,B)/4,4d(A,B)]$. With the well-separation property and the triangle inequality,

[TABLE]

which leads to the inequality

[TABLE]

since $\delta\leq 1/10$. Moreover, $\|p_{i}-p_{j}\|\leq 2\rho$, and a similar calculation shows that

[TABLE]

which implies that $\rho\in\tilde{R}_{(A,B)}$. Since $p_{i}\in P_{A}$ or $p_{i}\in P_{B}$, this implies immediately that $\tilde{r}_{p_{i}}(x)\geq\rho$ for $x\geq\rho$. ∎

Combining Lemma 4.3 with Lemma 4.2 completes the proof of Lemma 4.1.

4.3 Pixelization

We next define a collection of pixels $\mathcal{LS}_{\alpha}$ for each $\alpha\in I$. In contrast to Section 3, these pixels will be taken from different grids. For technical reasons, we enlarge the critical intervals $R_{(A,B)}$ by the smallest amount such that their endpoints are points in $I$. The definitions of $R_{p}$, $r_{p}(\alpha)$ and $\mathcal{LB}_{\alpha}$ get adapted accordingly, without affecting the claims of Lemma 4.2 and Lemma 4.3.

We assume again without loss of generality that we start the construction at a minimal scale $\alpha_{0}$ that is sufficiently small and set $\mathcal{LS}_{\alpha_{0}}=P$. We define $\mathcal{LS}_{\alpha}$ for $\alpha>\alpha_{0}$ inductively: let $\beta$ denote the largest power of $2$ that is not larger than $\alpha$ as in Section 3. Call a point $p\in P$ active if $r_{p}(\alpha)=\alpha$, and inactive otherwise. Let $\mathcal{LS}^{\prime}_{\alpha}$ be the collection of all pixels whose centers lie in the any ball of radius $(1+\varepsilon/2)r_{p}(\alpha)$ at an active points $p$. Define $\mathcal{LS}_{\alpha}$ as the inclusion-maximal pixels of the set $\mathcal{LS}_{\alpha/(1+\varepsilon)}\cup\mathcal{LS}^{\prime}_{\alpha}$. We denote by $|\mathcal{LS}_{\alpha}|$ the union of the pixels in $\mathcal{LS}_{\alpha}$. See Figure 6 for an illustration. We use the standard filling technique to define $|\mathcal{LS}_{\alpha}|$ for every $\alpha\geq 0$.

To rephrase our construction: at scale $\alpha$, we resample the pixels of all balls which are growing from $\alpha/(1+\varepsilon)$ to $\alpha$ in the filtration $(\mathcal{LB}_{\alpha})_{\alpha\geq 0}$, and add them to the set of pixels, removing lower-scale pixels which are covered by the new ones.

The idea of the digitization is to approximate $\mathcal{LB}_{\alpha}$ with pixels. Indeed, we have that $\mathcal{LB}_{\alpha}\subseteq|\mathcal{LS}_{\alpha}|$. That follows inductively, observing that the balls that grow from $\mathcal{LB}_{\alpha}$ to $\mathcal{LB}_{(1+\varepsilon)\alpha}$ are covered by pixels. It is also true that $|\mathcal{LS}_{\alpha}|\subseteq\mathcal{LB}_{\alpha}^{1+\varepsilon}$, where $\mathcal{LB}_{\alpha}^{1+\varepsilon}$ is obtained from $\mathcal{LB}_{\alpha}$ by enlarging every ball by a factor of $(1+\varepsilon)$. The proofs are similar to the proof of Lemma 3.1. However, it is not true that $\mathcal{LB}_{\alpha}^{1+\varepsilon}=\mathcal{LB}_{(1+\varepsilon)\alpha}$ in general (more precisely, neither $\mathcal{LB}_{\alpha}^{1+\varepsilon}\subseteq\mathcal{LB}_{(1+\varepsilon)\alpha}$ nor $\mathcal{LB}_{\alpha}^{1+\varepsilon}\supseteq\mathcal{LB}_{(1+\varepsilon)\alpha}$ holds in general). Therefore, the simple sandwiching strategy from Section 3 fails. We need a more refined argument, showing that the interleaving still works on the homology level.

For that, the following observation is crucial:

Lemma 4.4**.**

The spaces $\mathcal{LB}_{\alpha}^{1+\varepsilon}$ and $\mathcal{LB}_{(1+\varepsilon)\alpha}$ are homotopically equivalent, and the maps $f_{1},f_{2}$ realizing this homotopic equivalence can be chosen as deformation retracts to the intersection $\mathcal{LB}_{\alpha}^{1+\varepsilon}\cap\mathcal{LB}_{(1+\varepsilon)\alpha}$.

The statement follows from the following statement on unions of balls, noting that the balls forming $\mathcal{LB}_{\alpha}^{1+\varepsilon}$, $\mathcal{LB}_{(1+\varepsilon)\alpha}$, and $\mathcal{LB}_{\alpha}^{1+\varepsilon}\cap\mathcal{LB}_{(1+\varepsilon)\alpha}$ all have the same nerve by construction.

Lemma 4.5**.**

Let $\mathcal{A}:=\{A_{1},\ldots,A_{n}\}$, $\mathcal{B}:=\{B_{1},\ldots,B_{n}\}$ be collections of closed balls such that $A_{i}\subseteq B_{i}$ and $A_{i}$ and $B_{i}$ have the same center for all $i=1,\ldots,n$. If $\mathrm{nerve}(\mathcal{A})=\mathrm{nerve}(\mathcal{B})$, there is strong deformation retraction from $|\mathcal{B}|$ to $|\mathcal{A}|$.

We only briefly sketch the proof idea and postpone details to an extended version. Note that because the balls are closed, there exists an $\varepsilon>0$, such that we can increase all balls of $\mathcal{B}$ without changing the nerve. This yields a collection of balls $\mathcal{C}:=\{C_{1},\ldots,C_{n}\}$ with the same centers, such that $A_{i}\subsetneq C_{i}$. We can define a collection of distance-like functions $f_{1},\ldots,f_{n}$, such that $A_{i}$ is the sublevel set of $f_{i}$ for value $1$, and $C_{i}$ is the sublevel set of $f_{i}$ for value $2$. Let $f$ denote the lower envelope of these functions. This function is continuous, but not differentiable; however, using techniques similar to the generalized gradient [11], it is possible to define a flow on $|\mathcal{C}|$. The singularities of this flow are outside of $|\mathcal{C}|\setminus|\mathcal{A}|$ because a singular points of the generalized gradient triggers a change in the nerve. That flow hence defines a deformation retract from $|\mathcal{C}|$ to $|\mathcal{A}|$.

We can use the above to define a deformation retract from $|\mathcal{B}|$ to $|\mathcal{A}|$: Simply restricting the retract to $|\mathcal{B}|$ is not sufficient because a point of $|\mathcal{B}|$ may be pushed out of $|\mathcal{B}|$ during the flow. However, note that $\mathcal{C}$ was defined using some parameter $\varepsilon$. Call a point $x\in|\mathcal{B}$ $\varepsilon$-affected if there exists an $i$ such that $x\notin B_{i}$, but $x\in C_{i}$. For every $x\in|\mathcal{B}|$, we can find some $\varepsilon>0$ such that an open neighborhood of $x$ is not $\varepsilon$-affected. If $x$ is not $\varepsilon$-affected, its flow does not change when further decreasing $\varepsilon$. Hence, for every $x\in|\mathcal{B}|$, we can define its flow to be the limiting flow for $\varepsilon\to 0$. Because of the neighborhood property above, the resulting flow is continuous on $|\mathcal{B}|$, and defines the desired deformation retraction.

Lemma 4.6** (Lazy sandwich lemma).**

The two persistence modules $(H(|\mathcal{LS}_{\alpha}|))_{\alpha\geq 0}$ and $(H(\mathcal{LB}_{\alpha}))_{\alpha\geq 0}$ are $(1+\varepsilon)^{2}$-approximations of each other.

Proof.

We need to define interleaving maps such that the four diagrams of  (1) commute. We restrict our attention to scales in $I$ and show a $(1+\varepsilon)$-interleaving for these scales, which implies the result.

From $H(\mathcal{LB}_{\alpha})\to H(|\mathcal{LS}_{(1+\varepsilon)\alpha}|)$, we take the homology map induced by the inclusion $\mathcal{LB}_{\alpha}\subseteq|\mathcal{LS}_{\alpha}|$. With that, the diagram

[TABLE]

commutes, since it already commutes on a space level.

We define the map $\phi:H(|\mathcal{LS}_{\alpha}|)\to H(\mathcal{LB}_{(1+\varepsilon)\alpha})$ as the composition of the homology map $H(|\mathcal{LS}_{\alpha}|)\to H(\mathcal{LB}_{\alpha}^{1+\varepsilon})$ induced by inclusion and the map $f^{\ast}:H(\mathcal{LB}_{\alpha}^{1+\varepsilon})\to H(\mathcal{LB}_{(1+\varepsilon)\alpha})$, where $f$ is the map from above realizing the homotopy equivalence. We consider the diagram

[TABLE]

where $c=(1+\varepsilon)$ and all maps except $\phi$ are induced by inclusion. Now, fixing a (singular) $p$-cycle $z$ over $|\mathcal{LS}_{\alpha}|$, composition with $f$ yields a $p$-cycle $f(z)$ in $H(\mathcal{LB}_{(1+\varepsilon)\alpha})$ which in turn includes as a $p$-cycle in $|\mathcal{LS}_{(1+\varepsilon)^{2}\alpha^{\prime}}|$. Let $f^{\prime}$ denote the map $\mathcal{LB}_{(1+\varepsilon)\alpha}\to\mathcal{LB}_{\alpha}^{1+\varepsilon}$ in the opposite direction of the homotopy equivalence. Since $f^{\prime}$ is realized as a deformation retraction to the intersection, and $f(z)$ lies in the intersection, it follows that $f^{\prime}(f(z))=f(z)$. On the other hand, $f^{\prime}\circ f$ is homotopic to the identity, hence, $f(z)$ is homotopic to $z$. This proves that the $p$-cycles $z$ and $f(z)$ represent the same homology class, proving that the above diagram commutes. We skip the remaining two diagrams, which can be handled with similar methods. ∎

4.4 Approximation tower and equivalence to the lazy pixels filtration

We define our approximation complex as $\mathcal{LX}_{\alpha}:=\mathrm{nerve}(\mathcal{LS}_{\alpha})$. Since $\mathcal{LS}_{\alpha}$ is a collection of pixels, $\mathcal{LX}_{\alpha}$ is a flag complex using Lemma 3.3. In general, $\mathcal{LX}_{\alpha}$ does not include into $\mathcal{LX}_{(1+\varepsilon)\alpha}$ because pixels might be removed when passing to a larger grid.

We define a map $g^{\prime}:\mathcal{LS}_{\alpha}\to\mathcal{LS}_{(1+\varepsilon)\alpha}$ that maps a pixel of $\mathcal{LS}_{\alpha}$ to the unique pixel in $\mathcal{LS}_{(1+\varepsilon)\alpha}$ that contains it. Because $\mathcal{LX}_{\alpha}$ is a flag complex and intersecting pixels are mapped to intersecting pixel under $g^{\prime}$, it follows directly that $g^{\prime}$ induces a simplicial map $g:\mathcal{LX}_{\alpha}\to\mathcal{LX}_{(1+\varepsilon)\alpha}$. Using the same proof strategy as in Appendix A.1, we can show that the towers $(\mathcal{LX}_{\alpha})_{\alpha\geq 0}$ and $(|\mathcal{LS}|_{\alpha})_{\alpha\geq 0}$ are persistence-equivalent.

We summarize the results of the section so far.

Theorem 4.7**.**

Using a $\frac{\varepsilon}{8}$-WSPD, our construction yields a persistence module $(H(\mathcal{LX}_{\alpha}))_{\alpha\geq 0}$ that is a $(1+\varepsilon)^{3}$-approximation of the Čech persistence module.

Proof.

$(H(\mathcal{LS}_{\alpha}))_{\alpha\geq 0}$ is a $(1+\varepsilon)^{2}$-approximation of $(H(\mathcal{LB}_{\alpha}))_{\alpha\geq 0}$ from Lemma 4.6. Further, the persistence module $(H(\mathcal{LB}_{\alpha}))_{\alpha\geq 0}$ is a $(1+\varepsilon)$-approximation of the Čech persistence module from Lemma 4.1. ∎

We can easily get a $(1+\varepsilon)$-interleaving as well by using $\varepsilon^{\prime}=\varepsilon/4$ in our approximation, since $\left(1+\varepsilon^{\prime}\right)^{3}=\left(1+\frac{\varepsilon}{4}\right)^{3}<1+\varepsilon$.

4.5 Size and computation time

Since our construction heavily depends on WSPDs, we state the well-known complexity bound here:

Theorem 4.8** ([8, 19]).**

Given $\delta\in(0,1]$, a $\delta$-WSPD of size $n\left(\frac{1}{\delta}\right)^{O(d)}$ can be computed in time

[TABLE]

Since we require an $(\varepsilon/8)$-WSPD for our construction, the same bounds hold up to a constant factor of the form $2^{O(d)}$ when replacing $\delta$ by $\varepsilon$.

Theorem 4.9**.**

The size of the $k$-skeleton of the tower $(\mathcal{LX}_{\alpha})_{\alpha\geq 0}$ is

[TABLE]

Proof.

Recall that the size of a tower is the number of simplices added in total when increasing the scale. We first bound the number of pixels that are added to the tower by $n\left(\frac{d}{\epsilon}\right)^{O(d)}$. This is done by a charging argument, where we charge the inclusion of a pixel into the tower to one pair in the WSPD and show that each pair is charged at most $\left(\frac{d}{\epsilon}\right)^{O(d)}$ times. The claim follows then with the size bound of the WSPD.

When constructing $\mathcal{LX}_{\alpha}$, we add new pixels to sample a ball around $p\in P$ only if $r_{p}(\alpha)=\alpha$. That, in turn, means that there is a WSP $(A,B)$ such that $\alpha\in R_{(A,B)}$ and $p$ is the representative of $A$ or of $B$. We charge the new pixels for the $\alpha$-ball around $p$ to the WSP $(A,B)$. Since $(A,B)$ has two representatives, it gets charged only by the pixels of two balls. By an analogue packing argument as in Theorem 3.7, an $\alpha$ ball gives rise to not more than $\left(\frac{d}{\epsilon}\right)^{d}$ pixels, so $(A,B)$ gets charged by that many pixels per scale in the worst case.

Finally, the pair $(A,B)$ can only get charged for scales that are included in $R_{(A,B)}$. A simple calculation shows that the number of such scales is the smallest $k$ that satisfies $(1+\varepsilon)^{k}\geq 64$ which can be solved as $k=O(1/\varepsilon)$. That means, multiplying with number of scales on which $(A,B)$ can be charged does not change the bound. This proves the bound on the number of pixels.

We next bound the number of simplex inclusions. Again, we use a charging argument, charging the inclusion of a simplex to one pixels in its boundary that has minimal size among all pixels in the simplex. We show that every pixel gets charged at most $2^{O(dk)}$ times, which completes the proof. A pixel $C$ can only have up to $3^{d}-1$ neighbors that are of same or greater size than $C$. Hence, at the scale $\alpha$ where it gets included, $C$ can only get charged for $3^{dk}=2^{O(dk)}$ simplices of dimension $\leq k$. However, $C$ might also get charged at scales larger than $\alpha$. This, in turn, can only happen if $C$ gets a new neighbor along a face on which it was not intersecting another pixel before (otherwise, the simplex is not included, but originates from a previous scale). Since this can only happen $3^{d}-1$ times, the total number of simplices that $C$ gets charged for is still $2^{O(dk)}$. ∎

Next, we describe an algorithm to compute the tower. By “computing” a tower, we mean the following: the output is a list of tokens of three different types:

• •

“(scale, $s$)”, denoting that the tower passes to scale $s$

• •

“(add, $\sigma$)”, denoting that the simplex $\sigma$ is added to the tower

• •

“(contract, $v_{0},v_{1}$)”, denoting that the vertices $v_{0}$ and $v_{1}$ are contracted in the tower.

These three operations specify the tower completely, because every simplicial map can be written as a sequence of inclusions and contractions [16, 23].

A first step is to compute the $\frac{\varepsilon}{8}$-WSPD $W$. We traverse its pairs and compute the active ranges of each pair. Let

[TABLE]

Clearly, the scales in $I$ that lie outside of $U$ can be disregarded because no new pixels are added to $\mathcal{LS}_{\alpha}$. In total, the algorithm considers only $O(1/\varepsilon)$ scales per WSP, so $n\left(\frac{d}{\epsilon}\right)^{O(d)}$ scales in total (and these scales can be computed in the same complexity).

For each scale, the algorithm determines the new pixels to be sampled. On a scale $\alpha$, we can efficiently determine the points $p\in P$ where $r_{p}(\alpha)=\alpha$, for instance using a priority queue. For each active point, we compute the new pixels of its $\alpha$-ball, using the same flooding algorithm as in Appendix A.2. The complexity is proportional to the number of pixels added, and thus the total complexity is bounded by $n\left(\frac{d}{\epsilon}\right)^{O(d)}$.

Next, the algorithm computes the contraction needed at scale $\alpha$. Note that this problem simply means to find pixels in $\mathcal{LS}_{\alpha/(1+\varepsilon)}$ which are covered by a pixel in $\mathcal{LS}_{\alpha}$. Each such pair of pixels defines one contraction of the old pixel into the new one. To find them efficiently, we store the collection of pixels of $\mathcal{LS}_{\alpha}$ as leaves of a (compressed) quad-tree $Q$ [19]. Whenever a new pixel is added to $\mathcal{LS}_{\alpha}$, it is also added to $Q$ as a new node $v$. All leaves in $Q$ that are children of $v$ are removed from the quad-tree, and the corresponding contraction is enlisted.

Finally, the algorithm finds the $k$-simplices with $k\geq 1$ added to $\mathcal{LX}_{\alpha}$. We only discuss the case $k=1$; the simplices for higher $k$ can be combinatorially enlisted in time proportional to the number of simplices because $\mathcal{LX}_{\alpha}$ is a flag complex (see the end of Appendix A.2). For edges, it suffices to iterate through the newly added vertices once more and query the quad-tree data structure for all neighbors of the corresponding leaf. It is not hard to see that all neighbors can be found in time $O(2^{d}+N)$, where $N$ is the number of neighbors encountered. Hence, the total time to find edges is also bounded by $n\left(\frac{d}{\epsilon}\right)^{O(d)}$. That concludes the description of the algorithm. In total, the running time is dominated by the construction of the WSPD and the enumeration step of all simplices, which is bounded by the size of the tower in Theorem 4.9. The space complexity is dominated by the size of the output to be produced. We summarize the result.

Theorem 4.10**.**

Computing the approximation tower takes time

[TABLE]

and space

[TABLE]

5 Conclusion

We presented an approximation scheme for approximating Čech filtrations in Euclidean space. The main idea is that by introducing additional sample points, we can asymptotically reduce the size of the complex, because the number of close-by points is reduced to a constant independent of $n$ and $\varepsilon$.

We briefly discuss the extensions of our results mentioned in the introduction. While turning our tower into a (flag complex) filtration is mostly a matter of technicalities, employing the permutahedral instead of the cubical grid requires an additional bag of techniques. The major complication comes from the fact that we lose the property of cubical grids that one pixel on the lower scale is always completely contained in a single pixel on the next scale. For instance, this property made the construction of the simplicial maps connecting the scales natural and easy to prove. For the permutahedral case, we need heavier algebraic machinery, including barycentric subdivisions and acyclic carriers to obtain similar results.

While the introduction of sample points yields improved theoretical guarantees, the algorithm in its current form is impractical, because the number of sample points is large even for well-conditioned problem instances. We discuss two possibilities to reduce that number. First of all, instead of covering an $\alpha$-ball uniformly by pixels of a certain scale, we could employ a quad-tree like subdivision and approximate the interior of the ball with larger pixels in practice. Moreover, it might be advisable to reduce the obtained complex further by applying elementary collapses to it which do not change the homotopy type. In very recent work [26], it has been demonstrated how such collapses can be performed efficiently in a simplicial complex. While these extensions would certainly help in practice, the question remains whether they can also lead to a further asymptotic improvements.

Acknowledgements

Aruni Choudhary is partly supported by ERC grant StG 757609. Michael Kerber is supported by Austrian Science Fund (FWF) grant number P 29984-N35. Sharath Raghvendra acknowledges support of NSF CRII grant CCF-1464276.

Appendix A Missing proofs for Section 3

A.1 Proof of Lemma 3.5

In this subsection we prove that the filtration of pixels $(|\mathcal{S}_{\alpha}|)_{\alpha>0}$ is persistence-equivalent to the tower of nerves $(\mathcal{X}_{\alpha})_{\alpha>0}$. For that, we introduce the intermediate object

[TABLE]

where $\beta$ ranges over scales in $I$. We can assume for simplicity that the range of scales is finite, for instance by ignoring all scales which are smaller than a third times the closest point distances among points in $P$, and those scales that are larger than the diameter of $P$. In that way, $S^{\prime}_{\alpha}$ also becomes a finite object.

We define $Y_{\alpha}:=\mathrm{nerve}(S^{\prime}_{\alpha})$ for $\alpha\in I$, and by the filling technique extend it to $(Y_{\alpha})_{\alpha\geq 0}$. This is a filtration, because $S^{\prime}_{\alpha}\subseteq S^{\prime}_{(1+\varepsilon)\alpha}$. By the Sandwich Lemma, we also have that $|S^{\prime}_{\alpha}|=|S_{\alpha}|$. Since the Nerve isomorphism on homology commutes with maps induced by inclusion [13], we obtain directly that $(Y_{\alpha})_{\alpha\geq 0}$ and $(|\mathcal{S}_{\alpha}|)_{\alpha>0}$ are persistence-equivalent. In the rest of the section, we will show that also $(Y_{\alpha})_{\alpha\geq 0}$ and $(\mathcal{X}_{\alpha})_{\alpha\geq 0}$ are persistence-equivalent.

Note that the Nerve Theorem implies that for each $\alpha$, $\mathcal{X}_{\alpha}$ and $Y_{\alpha}$ are homotopically equivalent, hence their homology groups are isomorphic. We will construct such an isomorphism explicitly. Recall from Section 3 that we have constructed a simplicial map $g:\mathcal{X}_{\alpha}\to\mathcal{X}_{(1+\varepsilon)\alpha}$ which maps a pixel of $\mathcal{S}_{\alpha}$ to the pixel of $\mathcal{S}_{(1+\varepsilon)\alpha}$ that contains it. For $\beta\leq\alpha\in I$, we can define the map $g_{\beta\to\alpha}$ as the composition of such $g$-maps, yielding a simplicial map from $\mathcal{X}_{\beta}$ to $\mathcal{X}_{\alpha}$.

Next, we define a vertex map $g_{\to\alpha}^{\prime}:S^{\prime}_{\alpha}\to S_{\alpha}$ as follows: for a pixel $C$ in $S^{\prime}_{\alpha}$, let $\beta\leq\alpha$ be such that $C\in S_{\beta}$. Then, $g_{\to\alpha}^{\prime}(C):=g_{\beta\to\alpha}(C)$. In words, the map assigns to a pixel at scale $\leq\alpha$ the (unique) pixel at scale $\alpha$ that contains it.

Lemma A.1**.**

The map $g_{\to\alpha}^{\prime}$ induces a simplicial map $g_{\to\alpha}:Y_{\alpha}\to\mathcal{X}_{\alpha}$.

Proof.

Note that $Y_{\alpha}$ is a flag complex by Lemma 3.3, so it suffices to consider edges. The statement follows as in the proof of Lemma 3.4 from the fact that if two pixels are intersecting, they either coincide or still intersect when mapped to a larger scale. ∎

By the definition of the $g$-maps, we have that for $\beta\leq\alpha\in I$,

[TABLE]

which implies at once that the following diagram commutes

[TABLE]

To prove the equivalence of $(\mathcal{X}_{\alpha})$ and $(Y_{\alpha})$, it remains to show that the induced map

[TABLE]

is an isomorphism. We do so by showing that the map $\mathrm{inc}^{\ast}:H(\mathcal{X}_{\alpha})\to H(Y_{\alpha})$ induced by inclusion is inverse to $g_{\to\alpha}^{\ast}$. Note that $g_{\to\alpha}$ is the identity on the subcomplex $\mathcal{X}_{\alpha}$, hence $g_{\to\alpha}\circ inc$ is the identity map on $\mathcal{X}_{\alpha}$ as well. It follows that $g_{\to\alpha}^{\ast}\circ inc^{\ast}$ is the identity map on the homology groups.

For the other direction, we make use of the following standard definition [25]: two simplicial maps $f_{1},f_{2}:K\rightarrow L$ are called contiguous if for every simplex $\sigma=(v_{0},\ldots,v_{k})\in K$, the set of vertices $\{f_{1}(v_{0}),\ldots,f_{1}(v_{k}),f_{2}(v_{0}),\ldots,f_{2}(v_{k})\}$ forms a simplex in $L$. Two contiguous maps are homotopic and therefore, the induced maps $f_{1}^{\ast}$, $f_{2}^{\ast}$ on homology are equal.

Lemma A.2**.**

The maps $\mathrm{id}_{Y_{\alpha}}$ and $\mathrm{inc}\circ g_{\to\alpha}:Y_{\alpha}\to Y_{\alpha}$ are contiguous.

Proof.

Fix any simplex $\sigma=(x_{0},\ldots,x_{k})$ in $Y_{\alpha}$ and consider the set of vertices

[TABLE]

We have to show that this set spans a simplex in $Y_{\alpha}$. Since $Y_{\alpha}$ is a flag complex, it suffices to show that all edges formed by two (distinct) vertices are in $Y_{\alpha}$. Since $\sigma\in Y_{\alpha}$, every pair $(x_{i},x_{j})$ is in $Y_{\alpha}$ as well. Since $g_{\to\alpha}$ is simplicial, the same is true for pairs of the form $(g_{\to\alpha}(x_{i}),g_{\to\alpha}(x_{j}))$. For pairs of the form $(x_{i},g_{\to\alpha}(x_{i}))$, note that $g_{\to\alpha}(x_{i})$ contains $x_{i}$, so the edge is also contained in the nerve. Finally, for pairs $(x_{i},g_{\to\alpha}(x_{j}))$ with $i\neq j$, we observe that the pixels $x_{i}$ and $x_{j}$ are intersecting and since $g_{\to\alpha}(x_{j})$ contains $x_{j}$, the pair is also intersecting. ∎

This lemma shows that for all $\alpha$, $g_{\to\alpha}^{\ast}$ is an invertible map and hence an isomorphism. Therefore, the filtrations $(Y_{\alpha})_{\alpha\geq 0}$ and $(\mathcal{X}_{\alpha})_{\alpha\geq 0}$ are persistence equivalent, finishing the proof of Lemma 3.5.

A.2 Algorithm to compute the tower

First, we describe the flooding algorithm to find the pixels in the $\alpha$-ball of each point $p\in P$. The complex $\mathcal{X}_{\alpha}$ is built on the lattice $L_{\beta}$. To find the pixels with centers inside $B_{1}:=B(p,\alpha)$, we first map $p$ to the nearest lattice point $x$ in $L_{\beta}$. Since the diameters of the pixels are upper bounded by $\varepsilon\beta/4$, a simple triangle inequality shows that for each pixel whose center lies in $B_{1}$, the pixel’s center and the centers of its neighboring pixels all lie within $B_{2}:=B(x,\alpha+\varepsilon\beta)$. It suffices to inspect each pixel in $B_{2}$ and to add it to the list of pixels if the distance of its center to $p$ is at most $\alpha$. We simply enqueue the neighbors of $x$ in a queue $Q$. We dequeue an element $y\in Q$ and inspect whether the distance to $p$ is at most $\alpha$; if so, we add it to the list of pixels of $p$. Also, if the distance to $p$ is at most $\alpha+\varepsilon\beta$, we enqueue the neighbors of $y$ into $Q$. We continue dequeuing until $Q$ is empty.

Since the neighborhood of each pixel $y$ whose center is in $B_{1}$ is completely contained inside $B_{2}$, it is easy to see that there is a sequence of pixels $x,c_{1},\ldots,c_{m},y$ such that consecutive pixels intersect. In this way, the search finds $y$, and terminates after a finite number of steps.

To compute the $k$-skeleton, we store the pixels at $\alpha$ in a dictionary. We go over each pixel and inspect whether its $3^{d}-1$ neighbors lie in the dictionary; if so, we put an edge between the pixel and its neighbor. The $k$-skeleton can be obtained from the $1$-skeleton by a simple combinatorial algorithm, traversing the Hasse diagram of the complex (refer to Algorithm 5.8 in [15]).

A.2.1 Proof of Theorem 3.7

Proof.

First, we bound the computation time required to find the pixels associated to each point $p\in P$. We find $x$, the closest point to $p$ in the lattice, and this takes $O(d)$ time. To find the pixels, we use a breadth-first search. A packing argument similar to the one in Theorem 3.7 shows that the number of vertices of this graph is at most $\left(\frac{1}{\varepsilon}\right)^{d}2^{O(d\log d)}$. Each vertex of this graph has degree $2^{O(d)}$. Inspecting the distance of a pixel-center to $p$ takes $O(d)$ time. The breadth-first search then requires $\left(\frac{1}{\varepsilon}\right)^{d}2^{O(d\log d)}$ time. Since there are $n$ points in $P$, the total time required is $n\left(\frac{1}{\varepsilon}\right)^{d}2^{O(d\log d)}$.

Computing the $k$-skeleton requires $2^{O(dk)}$ time per pixel. To compute the simplicial map, it suffices to report the vertex map to the next scale in the tower. When the lattice does not change between scales, the vertex map is simply an inclusion, so there is nothing to compute. Otherwise, we compute the nearest lattice point to each pixel-center of the current scale and this takes $O(d)$ time per pixel-center. This cost is subsumed by the cost of the previous operations. The claim follows. ∎