Robust Topological Inference: Distance To a Measure and Kernel Distance

Introduction

Figure 1 shows three complex point clouds, based on a model used for simulating cosmology data. Visually, the three samples look very similar. Below the data plots are the persistence diagrams, which are summaries of topological features defined in Section 2. The persistence diagrams make it clearer that the third data set is from a different data generating process than the first two.

This is an example of how topological features can summarize structure in point clouds. The field of topological data analysis (TDA) is concerned with defining such topological features; see Carlsson (2009). When performing TDA, it is important to use topological measures that are robust to noise. This paper explores some of these robust topological measures.

The sublevel sets $L_{t}=\{x:\ \Delta_{S}(x)\leq t\}$ provide multiscale topological information about $S$. As $t$ varies from zero to $\infty$, topological features — connected components, loops, voids — are born and die. Persistent homology quantifies the evolution of these topological features as a function of $t$. See Figure 2. Each point on the persistence diagram represents the birth and death time of a topological feature.

Given a sample $X_{1},\ldots,X_{n}\sim P$, the empirical distance function is defined by

If $P$ is supported on $S$, and has a density bounded away from zero and infinity, then $\widehat{\Delta}$ is a consistent estimator of $\Delta_{S}$, i.e., $\sup_{x}|\widehat{\Delta}(x)-\Delta_{S}(x)|\stackrel{{\scriptstyle P}}{{\to}}0.$ However, if there are outliers, or noise, then $\widehat{\Delta}(x)$ is no longer consistent. Figure 3 (bottom) shows that a few outliers completely change the distance function. In the language of robust statistics, the empirical distance function has breakdown point zero.

A more robust approach is to estimate the persistent homology of the super-level sets of the density $p$ of $P$. As long as $P$ is concentrated near $S$, we expect the level sets of $p$ to provide useful topological information about $S$. Specifically, some level sets of $p$ are homotopic to $S$ under weak conditions, and this implies that we can estimate the homology of $S$. Note that, in this case, we are using the persistent homology of the super-level sets of $p$, to estimate the homology of $S$. This is the approach suggested by Bubenik (2012), Fasy et al. (2014b) and Bobrowski et al. (2014). A related idea is to use persistent homology based on a kernel distance (Phillips et al., 2014). In fact, the sublevel sets of the kernel distance are a rescaling of the super-level sets of $p$, so these two ideas are essentially equivalent. We discuss this approach in Section 5.

A different approach, more closely related to the distance function, but robust to noise, is to use the distance-to-a-measure (DTM), $\delta\equiv\delta_{P,m}$, from Chazal et al. (2011); see Section 2. An estimate $\widehat{\delta}$ of $\delta$ is obtained by replacing the true probability measure with the empirical probability measure $P_{n}$, or with a deconvolved version of the observed measure Caillerie et al. (2011). One then constructs a persistence diagram based on the sub-level sets of the DTM. See Figure 1. This approach is aimed at estimating the persistent homology of $S$. (The DTM also suggests new approaches to density estimation; see Biau et al. (2011).)

The density estimation approach and the DTM are both trying to probe the topology of $S$. But the former is using persistent homology to estimate the homology of $S$, while the DTM is directly trying to estimate the persistent homology of $S$. We discuss this point in detail in Section 9.1.

In this paper, we explore some statistical properties of these methods. In particular:

We show that $\sqrt{n}(\widehat{\delta}^{2}(x)-\delta^{2}(x))$ converges to a Gaussian process. (Theorem 5).

We show that the bootstrap provides asymptotically valid confidence bands for $\delta$. This allows us to identify significant topological features. (Theorem 18).

We find the limiting distribution of a key topological quantity called the bottleneck distance. (Section 4.1).

We also show that, under additional assumptions, there is another version of the bootstrap — which we call the bottleneck bootstrap — that provides more precise inferences. (Section 6).

We show similar results for the kernel distance. (Section 5).

We propose a method for choosing the tuning parameter $m$ for DTM and the bandwidth $h$ for the kernel distance. (Section 7.1).

We show that the DTM and the KDE both suffer from boundary bias and we suggest a method for reducing the bias. (Section 7.2).

Notation. $B(x,\epsilon)$ is a Euclidean ball of radius $\epsilon$, centered at $x$. We define $A\oplus\epsilon=\bigcup_{x\in A}B(x,\epsilon)$, the union of $\epsilon$-balls centered at points in $A$. If $x$ is a vector then $||x||_{\infty}=\max_{j}|x_{j}|$. Similarly, if $f$ is a real-valued fiction then $||f||_{\infty}=\sup_{x}|f(x)|$. We write $X_{n}\rightsquigarrow X$ to mean that $X_{n}$ converges in distribution to $X$, and we use symbols like $c,C,\ldots,$ as generic positive constants.

Remark: The computing for the examples in this paper were done using the R package TDA. See Fasy et al. (2014a). The package can be downloaded from http://cran.r-project.org/web/packages/TDA/index.html.

Remark: In this paper, we discuss the DTM which uses a smoothing parameter $m$ and the kernel density estimator which uses a smoothing bandwidth $h$. Unlike in traditional function estimation, we do not send these parameters to zero as $n$ increases. In TDA, the topological features created with a fixed smoothing parameter are of interest. Thus, all the theory in this paper treats the smoothing parameters as being bounded away from 0. See also Section 4.4 in Fasy et al. (2014b). In Section 7.1, we discuss the choice of these smoothing parameters.

Background

In this section, we define several distance functions and distance-like functions, and we introduce the relevant concepts from computational topology. For more detail, we refer the reader to Edelsbrunner and Harer (2010).

Let $L_{t}=\{x:\ \Delta_{S}(x)\leq t\}$. We will refer to the parameter $t$ as “time.”

Given the nested family of the sublevel sets of $\Delta_{S}$, the topology of $L_{t}$ changes as $t$ increases: new connected components can appear, existing connected components can merge, cycles and cavities can appear or be filled, etc. Persistent homology tracks these changes, identifies features and associates an interval or lifetime (from $t_{\textrm{birth}}$ to $t_{\textrm{death}}$) to them. For instance, a connected component is a feature that is born at the smallest $t$ such that the component is present in $L_{t}$, and dies when it merges with an older connected component. Intuitively, the longer a feature persists, the more relevant it is.

A feature, or more precisely its lifetime, can be represented as a segment whose extremities have abscissae $t_{\textrm{birth}}$ and $t_{\textrm{death}}$; the set of these segments is called the barcode of $\Delta_{S}$. An interval can also be represented as a point in the plane with coordinates $(u,v)=(t_{\textrm{birth}},t_{\textrm{death}})$. The set of points (with multiplicity) representing the intervals is called the persistence diagram of $\Delta_{S}$. Note that the diagram is entirely contained in the half-plane above the diagonal defined by $u=v$, since death always occurs after birth. This diagram is well-defined for any compact set $S$ (Chazal et al. (2012), Theorem 2.22). The most persistent features (supposedly the most important) are those represented by the points furthest from the diagonal in the diagram, whereas points close to the diagonal can be interpreted as (topological) noise.

Let $S_{1}$ and $S_{2}$ be compact sets with distance functions $\Delta_{1}$ and $\Delta_{2}$ and diagrams $D_{1}$ and $D_{2}$. The bottleneck distance between $D_{1}$ and $D_{2}$ is defined by

where the minimum is over all bijections between $D_{1}$ and $D_{2}$. In words, the bottleneck distance is the maximum distance between the points of the two diagrams, after minimizing over all possible pairings of the points (including the points on the diagonals).

A fundamental property of persistence diagrams is their stability. According to the Persistence Stability Theorem (Cohen-Steiner et al. (2005); Chazal et al. (2012))

Here, $H$ is the Hausdorff distance, namely,

Given a sample $X_{1},\ldots,X_{n}\sim P$, the empirical distance function is defined by

Suppose that $P$ is supported on $S$, and has a density bounded away from zero and infinity. Then

See also Cuevas and Rodríguez-Casal (2004). The previous lemma justifies using $\widehat{\Delta}$ to estimate the persistent homology of sublevel sets of $\Delta_{S}$. In fact, the sublevel sets of $\widehat{\Delta}$ are just unions of balls around the observed data. That is,

The persistent homology of the union of the balls as $t$ increases may be computed by creating a combinatorial representation (called a Cech complex) of the union of balls, and then applying basic operations from linear algebra (Edelsbrunner and Harer, 2010, Sections VI.2 and VII.1).

However, as soon as there is noise or outliers, the empirical distance function becomes useless, as illustrated in Figure 3. More specifically, suppose that

where $\pi\in$, $R$ is an outlier distribution (such as a uniform on a large set), $Q$ is supported on $S$, $\star$ denotes convolution, and $\Phi_{\sigma}$ is a compactly supported noise distribution with scale parameter $\sigma$.

Recovering the persistent homology of $\Delta_{S}$ exactly (or even the homology of $S$) is not possible in general since the problem is under-identified. But we would still like to find a function that is similar to the distance function for $S$. The empirical distance function fails miserably even when $\pi$ and $\sigma$ are small. Instead, we now turn to the DTM.

2 Distance to a Measure

Given a probability measure $P$, for $0<m<1$, the distance-to-measure (DTM) at resolution $m$ (Chazal et al., 2011) is defined by

where $G_{x}(t)=P(\|X-x\|\leq t)$. Alternatively, the DTM can be defined using the cdf of the squared distances, as in the following lemma:

Let $F_{x}(t)=P(\|X-x\|^{2}\leq t)$. Then

Given a sample $X_{1},\ldots,X_{n}\sim P$, let $P_{n}$ be the probability measure that puts mass $1/n$ on each $X_{i}$. It is easy to see that the distance to the measure $P_{n}$ at resolution $m$ is

where $k=\lceil mn\rceil$ and $N_{k}(x)$ is the set containing the $k$ nearest neighbors of $x$ among $X_{1},\ldots,X_{n}$. We will use $\widehat{\delta}$ to estimate $\delta$.

Now we summarize some important properties of the DTM, all of which are proved in Chazal et al. (2011). First, recall that the Wasserstein distance of order $p$ between two probability measures $P$ and $Q$ is given by

where the infimum is over all joint distributions $J$ for $(X,Y)$ such that $X\sim P$ and $Y\sim Q$. We say that $P$ satisfies the $(a,b)$-condition if there exist $a,b>0$ such that, for every $x$ in the support of $P$ and every $\epsilon>0$,

The next theorem summarizes results from Chazal et al. (2011) and Buchet et al. (2013).

If $P$ satisfies (11) and is supported on a compact set $S$, then

In particular, $\sup_{x}|\delta_{P,m}(x)-\Delta_{S}(x)|\to 0$ as $m\to 0$.

If $P$ and $Q$ are two distributions, then

If $Q$ satisfies (11) and is supported on a compact set $S$ and $P$ is another distribution (not necessarily supported on $S$), then

Hence, if $m\asymp W_{2}(P,Q)^{2b/(2+b)}$, then $\sup_{x}|\delta_{Q,m}(x)-\Delta_{S}(x)|=O(W_{2}(P,Q)^{2/(2+b)})$.

Let $D_{P}$ be the diagram from $\delta_{P,m}$ and let $D_{Q}$ be the diagram from $\delta_{Q,m}$. The bottleneck distance is bounded by

We conclude this section by bounding the distance between the diagrams $D_{\delta_{P,m}}$ and $D_{\Delta_{S}}$.

We first apply the stability theorem and part 3 in the previous result:

The term $W_{2}(P,Q)$ can be upper bounded as follows:

These two terms can be bounded with simple transport plans. Let $Z$ be a Bernoulli random variable with parameter $\pi$. Let $X$ and $Y$ be random variables with distributions $R$ and $Q\star\Phi_{\sigma}$. We take these three random variables independent. Then, the random variable $V$ defined by $V=ZX+(1-Z)Y$ has for distribution the mixture distribution $P$. By definition of $W_{2}$, one has

It can be checked in a similar way that $W_{2}\left(Q\star\Phi_{\sigma},Q\right)\leq\sigma$ (see for instance the proof of Proposition 1 in Caillerie et al. (2011)) and the Lemma is proved. ∎

Limiting Distribution of the Empirical DTM

In this section, we find the limiting distribution of $\widehat{\delta}$ and we use this to find confidence bands for $\delta(x)$. We start with the pointwise limit.

Let $\delta(x)\equiv\delta_{P,m}(x)$ and $\widehat{\delta}(x)\equiv\delta_{P_{n},m}(x)$, as defined in the previous section.

First suppose that $\widehat{F}_{x}^{-1}(m)>F_{x}^{-1}(m)$.

Then, by integrating “horizontally” rather than “vertically”, we can split the integral into two parts, as illustrated in Figure 4:

Next, it can be easily checked that (18) is also true when $\widehat{F}_{x}^{-1}(m)<F_{x}^{-1}(m)$ if we take $\int_{a}^{b}f(u)du:=-\int_{b}^{a}f(u)du$ when $a>b$. Now, since $F_{x}$ is differentiable at $m$, we have that \Bigl{|}F_{x}^{-1}(m)-\widehat{F}_{x}^{-1}(m)\Bigr{|}=O_{P}(1/\sqrt{n}), see for instance Corollary 21.5 in van der Vaart (2000). According to the DKW inequality we have that $\sup_{t}\left|F_{x}(t)-\widehat{F}_{x}(t)\right|=O_{P}(\sqrt{1/n})$ and thus

with $\lim_{u\rightarrow 0}\omega_{\mathcal{X}}(u)=\omega_{\mathcal{X}}(0)=0$. When such modulus of continuity $\omega$ exists, note that it always can be chosen non decreasing and this allows us to consider its generalized inverse $\omega^{-1}$.

One may ask if the existence of the uniform modulus of continuity over a compact domain $\cal X$ is a strong assumption or not. To answer this issue, let us introduce the following assumption:

Note that Assumption $\left(H_{\omega,\mathcal{X}}\right)$ is not very strong. For instance it is satisfied for a measure $P$ supported on a compact and connected manifold, with $P_{x}$ absolutely continuous for the Hausdorff measure on $P$. The following Lemma derives from general results on quantile functions given in Bobkov and Ledoux (2014) (see their Appendix A); the lemma shows that a uniform modulus of continuity for the quantiles exists under Assumption $\left(H_{\omega,\mathcal{X}}\right)$.

Let $x\in\mathcal{X}$. According to Proposition A.17 in Bobkov and Ledoux (2014), Assumption $\left(H_{\omega,\mathcal{X}}\right)$ is equivalent to the absolute continuity of $F_{x}^{-1}$ in $[0,1)$. We can then define a modulus of continuity of $F_{x}^{-1}$ by

According to Lemma 8, we have that for any $(x,x^{\prime})\in\mathcal{X}^{2}$:

where $C$ only depends on $P$ and $\cal X$. According to (19), for any $(m,m^{\prime})\in(0,1)^{2}$, and for any $(x,x^{\prime})\in\mathcal{X}^{2}$:

By taking the supremum over the $m$ and the $m^{\prime}$ such that $|m^{\prime}-m|<u$, it yields:

and $x\mapsto\omega_{x}(u)$ is thus Lipschitz at any $u$. For any $u\in(0,1)$, let

which gives that $\omega_{\cal X}(u_{\phi(n)})$ and $\omega_{\cal X}(u_{n})$ both tend to zero because $\omega_{\bar{x}}$ is continuous at zero. Thus $\omega_{\cal X}$ is continuous at zero and the Lemma is proved. ∎

Therefore $\sqrt{F_{x+a}^{-1}(m)}\leq\sqrt{F_{x}^{-1}(m)}+\|a\|$. Similarly,

which implies $\sqrt{F_{x}^{-1}(m)}\leq\sqrt{F_{x+a}^{-1}(m)}+\|a\|$.

We are now in position to state the functional limit of the distance to measure to the empirical measure.

Note that the functional limit is valid for any value of $m\in(0,1)$. A local version of this result could be also proposed by considering the (local) modulii of continuity of the quantile functions at $m$. For the sake of clarity, we prefer to give a global version.

In the proof of Theorem 5 we showed that $\sqrt{n}(\widehat{\delta}^{2}(x)-\delta^{2}(x))=A_{n}(x)+R_{n}(x)$ where

First, we show that $\sup_{x\in\mathcal{X}}|R_{n}(x)|=o_{P}(1)$. Then we prove that $A_{n}(x)$ converges to a Gaussian process.

Note that $|R_{n}(x)|\leq\frac{\sqrt{n}}{m}|S_{n}(x)||T_{n}(x)|$ where

Let $\xi_{i}\sim$Uniform (0,1), for $i=1,\dots,n$ and let $H_{n}$ be their empirical distribution function. Define $k=mn$. Then $\widehat{F}_{x}^{-1}(m)\stackrel{{\scriptstyle d}}{{=}}F_{x}^{-1}(\xi_{(k)})=F_{x}^{-1}\left(H_{n}^{-1}(m)\right)$, where $\xi_{(k)}$ is the $k$th order statistic. Thus, for any $m>0$ and any $x\in\mathcal{X}$:

In the last line we used inequality 1 page 453 and Point (12) of Proposition 1 page 455 of Shorack and Wellner (2009). Note that $\omega_{\mathcal{X}}^{-1}(\epsilon)>0$ for any $\varepsilon>0$ because $\omega_{\mathcal{X}}$ is assumed to be continuous at zero by definition.

Fix $\varepsilon>0$. There exists an absolute constant $C_{\mathcal{X}}$ such that there exists an integer $N\leq C_{\mathcal{X}}\varepsilon^{-d}$ and $N$ points $(x_{1},\dots,x_{N})$ laying in $\mathcal{X}$ such that $\bigcup_{j=1\dots N}B_{j}\supseteq\mathcal{X}$, where $B_{j}=B(x_{j},\varepsilon)$. Now, we apply Lemma 8 with $P$, and with $P_{n}$ and we find that for any $x\in B_{j}$:

where $C$ is a positive constant which only depends on $\cal X$ and $P$. Using a union bound together with (20) , we find that

Thus, $\sup_{x\in\mathcal{X}}|S_{n}(x)|=o_{P}(1)$. Then

Since $\sup_{x\in\mathcal{X}}|R_{n}(x)|=o_{P}(1)$, it only remains to prove that the process $A_{n}$ converges to a Gaussian process.

Now, we consider the process $A_{n}$ on $\mathcal{X}$. Let us denote $\nu_{n}:=\sqrt{n}(P_{n}-P)$ the empirical process. Note that

Hadamard Differentiability and The Bootstrap

In this section, we use the bootstrap to get a confidence band for $\delta$. Define $c_{\alpha}$ by

Let $X_{1}^{*},\ldots,X_{n}^{*}$ be a sample from the empirical measure $P_{n}$ and let $\widehat{\delta}^{*}$ be the corresponding empirical DTM. The bootstrap estimate $\widehat{c}_{\alpha}$ is defined by

As usual, $\widehat{c}_{\alpha}$ can be approximated by Monte Carlo. Below we show that this bootstrap is valid. It then follows that

A different approach to the bootstrap is considered in Section 6.

conditionally given $X_{1},X_{2},\ldots$, in probability.

We will establish the above result using the functional delta method, which entails showing that the distance to measure function is Hadamard differentiable at $P$. In fact, the proof further shows that the process

This result is consistent with the result established in Theorem 9, but in order to establish Hadamard differentiability, we use a slightly different assumption. Theorem 9 is proved by assuming an uniform modulus of continuity on the quantile functions $F_{x}^{-1}$ whereas in Theorem 11 an uniform lower bound on the derivatives is required. These two assumptions are consistent: they both say that $F^{-1}_{x}$ is well behaved in a neighborhood of $m$ for all $x$. However, (24) is stronger.

Let us first give the definition of Hadamard differentiability, for which we refer the reader to, e.g., Section 3.9 of van der Vaart and Wellner (1996). A map $\phi$ from a normed space $(\mathcal{D},\|\cdot\|_{\mathcal{D}})$ to a normed space $(\mathcal{E},\|\cdot\|_{\mathcal{E}})$ is Hadamard differentiable at the point $x\in\mathcal{D}$ if there exists a continuous linear map $\phi_{x}^{\prime}:\mathcal{D}\to\mathcal{E}$ such that

whenever $\|h_{t}-h\|_{\mathcal{D}}\to 0$ as $t\to 0$.

We also recall the functional delta method (see, e.g. van der Vaart and Wellner, 1996, Theorem 3.9.4): suppose that $T_{n}$ takes values in $\mathcal{D}$, $r_{n}\to\infty$, $r_{n}(T_{n}-\theta)\rightsquigarrow T$, and suppose that $\phi$ is Hadamard differentiable at $\theta$. Then $r_{n}(\phi(T_{n})-\phi(\theta))\rightsquigarrow\phi_{\theta}^{\prime}(T)$. Moreover, by Theorem 3.9.11 of van der Vaart and Wellner (1996) the bootstrap has the same limit. More precisely, given $X_{1},X_{2},\ldots$, we have that $r_{n}(\phi(T_{n}^{*})-\phi(T_{n}))$ converges conditionally in distribution to $\phi_{\theta}^{\prime}(T)$, in probability. This implies the validity of the bootstrap confidence sets.

The pair $(\mathcal{M},\|\cdot\|_{\mathcal{B}})$ is a normed space.

where the infimum over the empty set is define to be $\infty$. If $P$ is a probability measure and $m\in(0,1)$ then $F^{-1}_{P,x}(m)$ is just the $m$-th quantile of the random variable $\|X-x\|^{2}$, $X\sim P$.

Fix a $m\in(0,1)$ and let $\mathcal{M}_{m}=\mathcal{M}_{m}(\mathcal{X})$ denote the subset of $\mathcal{M}$ consisting of all finite signed measure $\mu$ such that, there exists a value of $r>0$ for which $\inf_{x\in\mathcal{X}}\mu\left(B(x,\sqrt{r})\right)\geq m$. Thus, for any $\mu\in\mathcal{M}_{m}$ and $x\in\mathcal{X}$, $F^{-1}_{\mu,x}(m)<\infty$. Let $\mathcal{D}_{m}$ be the image of $\mathcal{M}_{m}$ by the mapping (26).

Let $\mathcal{E}$ the set of bounded, real-valued function on $\mathcal{X}$, a normed space with respect to the sup norm. Finally, we define $\phi\colon\mathcal{D}_{m}\rightarrow\mathcal{E}$ to be the mapping

Notice that if $P$ is a probability measure, simple algebra shows that $\phi(P)(x)$ is the square value of the distance to measure of $P$ at the point $x$, i.e. $\delta^{2}_{p}(x)$; see Figure 5.

Below we will show that, for any probability measure $P$, the mapping (27) is Hadamard differentiable at $P$.

For an arbitrary $Q\in\mathcal{Q}$, let $\{Q_{t}\}_{t>0}\subset\mathcal{D}$ be a sequence of signed measure such that $\lim_{t\rightarrow 0}\|Q_{t}-Q\|_{\mathcal{B}}=0$ and such that $P+tQ_{t}\in\mathcal{D}_{m}$ for all $t$. Sequences of this form exist: since $\|tQ_{t}\|_{\mathcal{B}}\rightarrow 0$ as $t\rightarrow 0$, for any arbitrary $0<\epsilon<1-m$ and all $t$ small enough,

By the boundedness of $\mathcal{X}$ and compactness of $S$, this implies that there exists a number $r>0$ such that

so the image of $P+tQ_{t}$ by (27) is an element of $\mathcal{E}$ (i.e. it is a bounded function).

To demonstrate Hadamard differentiability (see 25), we will prove that, as $t\to 0$, the expression in (28), as a bounded function of $x\in\mathcal{X}$, will converge in $\mathcal{E}$ to the bounded function

Towards that end, we have, for all $t$ and any $x\in\mathcal{X}$,

where, for $a<b$, we write $\int_{b}^{a}=-\int_{a}^{b}$.

To handle the three terms appearing in the last display, we first state and prove two useful results.

We first prove (31). Let $\mathcal{A}_{x,t}=\{y\colon\|y-x\|^{2}\leq F^{-1}_{x,t}(m)\}$. Then,

Since $\sup_{x\in\mathcal{X}}tQ_{t}(\mathcal{A}_{x,t})\rightarrow 0$ as $t\rightarrow 0$ we obtain (30). The claim (31) follows from (30) using the facts that $F_{x}$ is monotone for each $x\in\mathcal{X}$ and that $\inf_{x\in\mathcal{X}}F_{x}^{\prime}(F^{-1}_{x}(m))>0$.

To show (29), the relation (32) combined with the fact that $m=F_{x}\left(F^{-1}_{x}(m)\right)$ for all $x$ yields that

for all $x\in\mathcal{X}$. By (31) and (24),

uniformly in $x\in\mathcal{X}$ and for all $t$ small enough. The relation (29) follows from the fact that $|Q_{t}(S)|=|Q(S)|+o(1)<\infty$. ∎

We now analyze the terms $A_{1}(x,t)$, $A_{2}(x,t)$ and $A_{3}(x,t)$ separately.

Term $A_{1}(x,t)$. As $t\rightarrow 0$, $Q_{t}\rightarrow Q$ and, uniformly in $x\in\mathcal{X}$ and $z>0$, $|Q_{t}(\mathcal{A}_{x,z})|\leq|Q_{t}(S)|=|Q(S)|+o(1)<\infty$. Furthermore, $\sup_{x\in\mathcal{X}}F^{-1}_{x}(m)<\infty$ by compactness of $\mathcal{X}$ and $S$. Therefore, using the dominated convergence theorem,

Term $A_{2}(x,t)$. Since $P(\mathcal{A}_{x,u})$ is non-decreasing in $u$ for all $x$, we have

Term $A_{3}(x,t)$. Finally, since $|Q_{t}(S)|\leq|Q(S)|+o(1)$ as $t\rightarrow 0$ and using (31), we obtain

Therefore, from (28), (33), (34), and (35),

is the Hadamard derivative of $\delta^{2}$ at $P$.

Fasy et al (2014) showed how to use the bootstrap to test the significance of a topological feature. They did this for distance functions and density estimators but the same idea works for DTM as we now explain.

Given a feature with birth and death time $(u,v)$, we will say that the feature is significant if $|v-u|>2c_{\alpha}/\sqrt{n}$ where $c_{\alpha}$ is defined by

In particular, $c_{\alpha}$ can be estimated from the bootstrap as we showed in the previous section. Specifically, define $\widehat{c}_{\alpha}$ by

Then $\widehat{c}_{\alpha}$ is a consistent estimate of $c_{\alpha}$.

To see why this makes sense, let ${\cal D}$ be the set of persistence diagrams. Let $D\equiv D_{\delta}$ be the true diagram and let $\widehat{D}\equiv D_{\widehat{\delta}}$ be the estimated diagram. Let

as $n\to\infty$. Now $|v-u|>2\widehat{c}_{\alpha}/\sqrt{n}$ if and only if the feature cannot be matched to the diagonal for any diagram in ${\cal C}$. (Recall that the diagonal corresponds to features with zero lifetime.)

We can visualize the significant features by putting a band of size $2c_{\alpha}/\sqrt{n}$ around the diagonal of $\widehat{D}$. See Figure 6.

Theory for Kernels

In this section, we consider an alternative to the DTM, namely, kernel based methods. This includes the kernel distance and the kernel density estimator.

Phillips et al. (2014) suggest using the kernel distance for topological inference. Given a kernel $K(x,y)$, the kernel distance between two probability measures $P$ and $Q$ is

It can be shown that $D_{K}(P,Q)=\|\mu_{P}-\mu_{Q}\|$ for vectors $\mu_{P}$ and $\mu_{Q}$ in an appropriate reproducing kernel Hilbert space (RKHS). Such distances are popular in machine learning; see Sriperumbudur et al. (2009), for example.

Given a sample $X_{1},\ldots,X_{n}\sim P$, let $P_{n}$ be the probability measure that put mass $1/n$ on each $X_{i}$. Let $\vartheta_{x}$ be the Dirac measure that puts mass one on $x$. Phillips et al. (2014) suggest using the discrete kernel distance

for topological inference. This is an estimate of the population quantity

The most common choice of kernel is the Gaussian kernel $K(x,y)\equiv K_{h}(x,y)=\exp\left(-\frac{\|x-y\|^{2}}{2h^{2}}\right)$, which has one tuning parameter $h$. We recall that, in topological inference, we generally do not let $h$ tend to zero. See the related discussion in Section 4.4 of Fasy et al. (2014b).

Recall that the kernel density estimator is defined by

Here, we used the fact that $n^{-1}\sum_{i=1}^{n}p(X_{i})=1+o_{P}(1)$ and $||\widehat{p}_{h}-p||_{\infty}=O_{P}(\sqrt{\log n/n})$.

We see that up to small order terms, the sublevel sets of $D_{K}(x)$ are a rescaled version of the super-level sets of the kernel density estimator. Hence, the kernel distance approach and the density estimator approach are essentially the same, up to a rescaling. However, $D^{2}_{K}$ has some nice properties; see Phillips et al. (2014).

The limiting properties of $\widehat{D}_{K}^{2}(x)$ follow immediately from well-known properties of kernel density estimators. In fact, the conditions needed for $\widehat{D}_{K}^{2}$ are weaker than for the DTM.

The Bottleneck Bootstrap

More precise inferences can be obtained by directly bootstrapping the persistence diagram. Define $\widehat{t}_{\alpha}$ by

The quantile $\widehat{t}_{\alpha}$ can be estimated by Monte Carlo. We then use a band of size $2\widehat{t}_{\alpha}$ on the diagram $D$.

In the following, we show that $\widehat{t}_{\alpha}$ consistently estimates the population value $t_{\alpha}$ defined by

The reason why the bottleneck bootstrap can lead to more precise inferences than the functional bootstrap from the previous section is that the functional bootstrap uses the fact that $W_{\infty}(\widehat{D},D)\leq||\widehat{\delta}-\delta||_{\infty}$ and finds an upper bound for $||\widehat{\delta}-\delta||_{\infty}$. But in many cases the inequality is not sharp so the confidence set can be very conservative. Moreover, we can obtain different critical values for different dimensions (connected components, loops, voids, …) and so the inferences are tuned to the specific features we are estimating. See Figure 7.

Although the bottleneck bootstrap can be used with either the DTM or the KDE, we shall only prove its validity for the KDE. First, we need the following result. For any function $p$, let $g=\nabla p$ denote its gradient and let $H=\nabla^{2}p$ denotes its Hessian. We say that $x$ is a critical point if $g(x)=(0,\ldots,0)^{T}$. We then call $p(x)$ a critical value. A function is Morse if the Hessian is non-degenerate at each critical point. The More index of a critical point $x$ is the number of negative eigenvalues of $H(x)$.

$q$ is a Morse function with exactly $k$ critical points $c_{1}^{\prime},\ldots,c_{k}^{\prime}$ say, and, after a suitable re-labeling of indices,

Moreover, $c_{j}$ and $c_{j}^{\prime}$ have the same Morse index.

Now let $\epsilon>0$ be small enough such that $2\epsilon<\min_{i\not=j}\|c_{i}-c_{j}\|$, and for any $i\not=j$, $p(B(c_{i},\epsilon))\cap p(B(c_{j},\epsilon))=\emptyset$. Then $\eta_{1}=\eta_{1}(\epsilon)=\min_{i\not=j}d(p(B(c_{i},\epsilon)),p(B(c_{j},\epsilon)))$ where $d(A,B)=\min_{a\in A,b\in B}|a-b|$ and $\eta_{2}=\eta_{2}(\epsilon)=\inf\{\|\nabla p(x)\|:x\in S\setminus\cup_{i=1}^{k}B(c_{i},\epsilon)\}$ are both positive. If $q$ satisfies the assumptions of the lemma for any $0<\eta\leq\min(\eta_{1},\eta_{2})$, then the critical values of $q$ have to be in $\cup_{i}p(B(c_{i},\epsilon))$ and the critical points $c_{i}^{\prime}$ have to be in $\cup_{i}B(c_{i},\epsilon)$.

More precisely, notice that since $p$ is a Morse function, for $\epsilon$ small enough, $\eta_{2}=O(\epsilon)$, and, for any $i\in\{1,\cdots,k\}$, the Taylor series of $\nabla p$ about $c_{j}$ yields

where $r(z)\rightarrow 0$ as $\|z\|\rightarrow 0$ and $H_{i}$ is the Hessian of $p$ at $c_{i}$. Let $\lambda_{\min}$ be the smallest absolute eigenvalue of the Hessians at all the critical points. Since $p$ is a Morse function, the matrix $H_{i}$ is full rank and $\lambda_{\min}$ is positive. As a consequence, for all $x\in S\setminus\cup_{i=1}^{k}B(c_{i},\epsilon)$ and $\epsilon$ small enough, $\|\nabla p(x)\|\geq\frac{\lambda_{\min}}{2}\epsilon$. Since $\eta_{1}$ is a non-incresing function of $\epsilon$, we have that, for $\epsilon$ small enough, $\eta=\eta_{2}\geq\frac{\lambda_{\min}}{2}\epsilon$.

To conclude the proof of the lemma, we need to prove that each ball $B(c_{i},\epsilon)$ contains exactly one critical point of $q$. Indeed, for $t\in$, the functions $q_{t}(x)=p(x)+t(q(x)-p(x))$ are Morse functions satisfying the same properties as $q$. Now, since each $c_{i}$ is a non-degenerate point of $p$, it follows from the continuity of the critical points (see, e.g. Prop. 4.6.1 in Demazure (2013)) that, restricting $\epsilon$ if necessary, there exist smooth functions $c_{i}:\to S$, $c_{i}(0)=c_{i},c_{i}(1)=c_{i}^{\prime}$ such that $c_{i}(t)$ is the unique critical point of $q_{t}$ in $B(c_{i},\epsilon)$. Moreover, since all the $q_{t}$ are Morse functions and since the Hessian of $q_{t}$ at $c_{i}(t)$ is a continuous function of $t$, then for any $t\in$, $c_{i}(t)$ is a non-degenerate critical point of $q_{t}$ with same index as $c_{i}$. ∎

Consider now two smooth functions such that the critical points are close, as illustrated in Figure 8. Next we show that, in this circumstance, the bottleneck distance takes a simple form.

Let $p$ and $q$ be two Morse functions as in Lemma 16, with finitely many critical points $C=\{c_{1},\ldots,c_{k}\}$ and $C^{\prime}=\{c_{1}^{\prime},\ldots,c_{k}^{\prime}\}$ respectively. Let $D_{p}$ and $D_{q}$ be the persistence diagrams from the upper level set filtrations of $p$ and $q$ respectively and let $a=\min_{i\neq j}|p(c_{i})-p(c_{j})|$ and $b=\max_{j}|p(c_{j})-q(c_{j}^{\prime})|$. If $b\leq a/2-||p-q||_{\infty}$ and $a/2>2||p-q||_{\infty}$, then $W_{\infty}(D_{p},D_{q})=b$.

The topology of the upper level sets of the Morse functions $p$ and $q$ only changes at critical values (Theorem 3.1 in Milnor (1963)). As a consequence the non diagonal points of $D_{p}$ (resp. $D_{q}$) have their coordinates among the set $\{p(c_{1}),\ldots,p(c_{k})\}$ (resp. $\{q(c_{1}^{\prime}),\ldots,p(c_{k}^{\prime})\}$) and each $p(c_{i})$ is the coordinate of exactly one point in $D_{p}$. Moreover, the pairwise distances between the points of $D_{p}$ are lower bounded by $a$ and all non diagonal points of $D_{p}$ are at distance at least $a$ from the diagonal. From the persistence stability theorem Cohen-Steiner et al. (2005); Chazal et al. (2012), $W_{\infty}(D_{p},D_{q})\leq||p-q||_{\infty}$. Since $a>4||p-q||_{\infty}$ and $a\geq 2b+2||p-q||_{\infty}$, the (unique) optimal matching realizing the bottleneck distance $W_{\infty}(D_{p},D_{q})$ is such that if $(p(c_{i}),p(c_{j}))\in D_{p}$ then it is matched to the point $(q(c_{i}^{\prime}),q(c_{j}^{\prime}))$ which thus have to be in $D_{q}$. It follows that $W_{\infty}(D_{p},D_{q})=b$. ∎

Now we establish the limiting distribution of $\sqrt{n}W_{\infty}(\widehat{D},D)$.

where $Z=(Z_{1},\ldots,Z_{k})\sim N(0,\Sigma)$ and

Let $\widehat{c}=\{\widehat{c}_{1},\widehat{c}_{2},\dots\}$ be the set of critical points of $\widehat{p}_{h}$. Let $g$ and $H$ be the gradient and Hessian of $p_{h}$. Let $\widehat{g}$ and $\widehat{H}$ be the gradient and Hessian of $\widehat{p}_{h}$. By a standard concentration of measure argument (and recalling that the support is compact), for any $\eta>0$ there is an event $A_{n,\eta}$ such that, on $A_{n,\eta}$,

For $\eta$ smaller than a fixed value $\eta_{0}$, we can apply Lemma 16, we get that on $A_{n,\eta}$, $\widehat{c}$ and $c$ have the same number of elements and can be indexed so that

where $C$ is the same constant is in Lemma 16. We then take $\eta_{n}:=\sqrt{\frac{\log n}{n}}$ and we consider the events $A_{n}:=A_{n,\eta_{n}}$. Then, for $n$ large enough, on $A_{n}$ we get

whereas $P\left(A_{n}^{c}\right)=o(1)$. In the following, we thus can restrict to $A_{n}$.

The critical values of $p_{h}$ are $v=(v_{1}\equiv p_{h}(c_{1}),\ldots,v_{k}\equiv p_{h}(c_{k}))$ and the critical values of $\widehat{p}_{h}$ are $\widehat{v}=(\widehat{v}_{1}\equiv\widehat{p}_{h}(\widehat{c}_{1}),\ldots,\widehat{v}_{k}\equiv\widehat{p}_{h}(\widehat{c}_{k}))$. Now we use Lemma 17 to conclude that $W_{\infty}(\widehat{D},D)=\max_{j}\|\widehat{v}_{j}-v_{j}\|_{\infty}$ for $n$ large enough. Hence,

Then, using a Taylor expansion, for each $j$,

Since $g(c_{j})=(0,\ldots,0)$ we can write the last equation as

For the second term, note that $\sqrt{n}(\widehat{c}_{j}-c_{j})=O(\log n)$ and $(\widehat{g}(c_{j})-g(c_{j}))=O_{P}(1/\sqrt{n})$. So

By the multivariate Berry-Esseen theorem (Bentkus, 2003),

By Lemma 17, $W_{\infty}(\widehat{D},D)=\max_{j}|\widehat{v}_{j}-v_{j}|$. The result follows. ∎

Let $X_{1}^{*},\ldots,X_{n}^{*}\sim P_{n}$ where $P_{n}$ is the empirical distribution. Let $\widehat{D}^{*}$ be the diagram from $\widehat{p}_{h}^{*}$ and let

be the bootstrap approximation to $F_{n}$.

Next we show that the bootstrap quantity $F_{n}(t)$ converges to the same limit as $F_{n}(t)$.

Assume the same conditions as the last theorem. Then,

The proof is essentially the same as the proof of Theorem 18 except that $\widehat{p}_{h}$ replaces $p_{h}$ and $\widehat{p}_{h}^{*}$ replaces $\widehat{p}_{h}$. Using the same notations as in the proof of Theorem 18, we note that on the set $A_{n}$, for $n$ larger than a fixed value $n_{0}$, the function $\widehat{p}_{h}$ is a Morse function with two uniformly bounded continuous derivatives and finitely many critical points $\widehat{c}=\{\widehat{c}_{1},\ldots,\widehat{c}_{k}\}$. We can restrict the analysis to sequence of events $A_{n}$ since $P(A_{n})$ tends to zero. Assuming that $A_{n}$ is satisfied, using the same argument as in Theorem 18, we get that:

where $\widetilde{Z}\sim N(0,\widehat{\Sigma})$ with

and $C_{2}^{*}$ depends on the empirical third moments of $h^{-d}K((x-X^{*})/h)$. There exists an upper bound $C_{2}$ on $C_{2}^{*}$, which only depends on $K$ and $P$. Since $\max_{j,k}|\widehat{\Sigma}_{j,k}-\Sigma_{j,k}|=O_{P}(\log n/\sqrt{n})$ and $\max_{j}\|\widehat{c}_{j}-c_{j}\|=O_{P}(\log n/\sqrt{n})$, we conclude that

Extensions

In this section, we discuss how to deal with three issues that can arise: choosing the parameters, correcting for boundary bias, and dealing with noisy data.

An unsolved problem in topological inference is how to choose the smoothing parameter $m$ (or $h$). Guibas et al. (2013) suggested tracking the evolution of the persistence of the homological features as the tuning parameter varies. Here we make this method more formal, by selecting the parameter that maximizes the total amount of significant persistence.

These measures are small when $m$ is small since $c_{\alpha}(m)$ is large. On the other hand, they are small when $m$ is large since then all the features are smoothed out. Thus we have a kind of topological bias-variance trade-off. We choose $m$ to maximize $N(m)$ or $S(m)$. The same idea can be applied to the kernel distance and kernel density estimator. See the example in Figure 9.

2 Boundary Bias

It is well known that kernel density estimators suffer from boundary bias. For topological inference, this bias manifests itself in a particular form and the same problem affects the DTM. Consider Figure 10. Because of the bounding box, many of the loops are incomplete. The result is that, using either the DTM or the KDE we will miss many of the loops.

There is a large literature on reducing boundary bias in the kernel density estimation literature. Perhaps the simplest approach is to reflect the data around the boundaries (see for example Schuster (1958)). But there is a simpler fix for topological inference: we merely need to close the loops at the boundary. This can be done by adding points uniformly around the boundary.

3 Two Methods for Improving Performance

We can improve the performance of all the methods if we cam mitigate the outliers and noise. Here we suggest two methods to do this. We focus on the kernel density estimator.

First, a simple method to reduce the number of outliers is to truncate the density, that is, we eliminate $\{X_{i}:\ \widehat{p}(X_{i})<t\}$ for some threshold $t$. Then we re-estimate the density.

Secondly, we sharpen the data as described in Choi and Hall (1999) and Hall and Minnotte (2002). The idea of sharpening is to move each data point $X_{i}$ slightly in the direction of the gradient $\nabla\widehat{p}(X_{i})$ and then re-estimate the density. The authors show that this reduces the bias at peaks in the density which should make it easier to find topological features. It can be seen that the sharpening method amounts to running one or more steps if the mean-shift algorithm. This is a gradient ascent which is intended to find modes of the density estimator. Given a point $x$, we move $x$ to

which is simply the local average centered at $x$. For data sharpening, we do one (or a few) iterations of this to each data point $X_{i}$. Then the density is re-estimated. In fact, we could also use the subspace constrained mean shift algorithm (SCMS) which moves points towards ridges of the density; see Ozertem and Erdogmus (2011). Figure 11 shows these methods applied to a simple example.

Examples

The data in Figure 12 are 10,000 data points on a 2D grid. We add Gaussian noise plus 1,000 outliers and compute the persistence diagrams of Kernel Density Estimator, Kernel distance, and Distance to Measure. The pink bands show 95% confidence sets obtained by bootstrapping the corresponding functions. The black lines show 95% confidence bands obtained with the bottleneck bootstrap for dimension 0, while the red lines show 95% confidence bands obtained with the bottleneck bootstrap for dimension 1. The Distance to Measure, which is less sensitive to the density of the points, correctly captures the topology of the data. The Kernel Distance and KDE find some extra significant connected component, corresponding to high density regions at the intersection of the grid.

Figure 13 shows the field position of two soccer players. The data come from body-sensor traces collected during a professional soccer game in late 2013 at the Alfheim Stadium in Tromso, Norway. The data are sampled at 20 Hz. See Pettersen et al. (2014). Although the data is a function observed over time, we treat it as a point cloud. Points on the boundary of the field have been added to avoid boundary bias. The DTM captures the difference between the two players: the defender leaves one big portion of the filed uncovered (1 significant loop in the persistence diagram), while the midfielder does not cover the 4 corners (4 significant loops). Nonetheless, the Kernel distance, which is more sensible to the density of these points, fails to detect significant topological features.

We will sample points around the the nodes, lines and faces that are formed at the intersection of the Voronoi regions. A Voronoi wall model is a sampling scheme that returns points within or around the Voronoi faces. Similarly, by sampling points exclusively around the lines or exclusively around the nodes, we can construct Voronoi filament models and Voronoi cluster models.

These models were introduced by Icke and van de Weygaert (1991) to mimic key features of cosmological data; see also van de Weygaert et al. (2011).

In this example we generate data from filament models and wall models using the basic definition of Voronoi diagram, computed on a fine grid in $^{3}$. We also add random Gaussian noise. See Figure 14: the first two rows show 100K particles concentrated around the filaments of 8 and 64 Voronoi cells, respectively. The last two rows show 100K particles concentrated around the walls of 8 and 64 Voronoi cells. 60K points on the boundary of the boxes have been added to mitigate boundary bias. For each model we present the persistent diagrams of the distance function, distance to measure and kernel density estimator. We chose the smoothing parameters by maximizing the quantity $S(\cdot)$, defined in Section 7.1.

The diagrams illustrate the evolution of the filtrations for the three different functions: at first, the connected components appear (black points in the diagrams); then they merge forming loops (red triangle), that eventually evolve into 3D voids (blue squares).

The persistence diagrams of the three functions allow us to distinguish the different models (see Figure 1 for a less trivial example) and the confidence bands, generated using the bootstrap method of Section 4.1, allow us to separate the topological signal from the topological noise. In general, the DTM performs better than the KDE, which is more affected by the high density of points around the nodes and filaments. For instance, this is very clear in the third row of Figure 14. The DTM diagram correctly captures the topology of the Voronoi wall model with 8 nuclei: one connected component and 8 voids are significant, while the remaining homological features fall into the band and are classified as noise.

Discussion

In this paper, we showed how the DTM and KDE can be used for robust topological inference. Further, we showed how to use the bootstrap to identify topological features that are distinguishable from noise. We conclude by discussing two issues: comparing DTM and KDE, and using persistent homology versus selecting a single level set.

The DTM and the KDE have the same broad aim: to provide a means for extracting topological features from data. However, these two methods are really focused on different goals. Consider again the model $P=\pi R+(1-\pi)(Q\star\Phi_{\sigma})$ and let $S$ be the support of $Q$. As before, we assume that $S$ is a “small set” meaning that either it has dimension $k<d$ or that it is full dimensional but has small Lebesgue measure. When $\pi$ and $\sigma$ are small, the persistent homology of the upper level sets of the density $p$ will be dominated by features corresponding to the homology of $S$. In other words, we are using the persistent homology of $\{p>t\}$ to learn about the homology of $S$. In contrast, the DTM is aimed at estimating the persistent homology of $S$. Both are useful, but they have slightly different goals.

This also raises the intriguing idea of extracting more information from both the KDE and DTM by varying more parameters. For example, if we look at the sets $\{p_{h}>t\}$ for fixed $t$ but varying $h$, we get information very similar to that of the DTM. Conversely, for the DTM, we can vary the tuning parameter $m$. There are many possibilities here which we will investigate in future work.

2 Persistent Homology Versus Choosing One Level Set

We have used the persistent homology of the upper level sets $\{\widehat{p}_{h}>t\}$ to probe the homology of $S$. This is the approach used in Bubenik (2012) and Phillips et al. (2014).

Bobrowski et al (2104) suggest a different approach. They select a particular level set $\{p>t\}$ and they form a robust estimate of the homology of this one level set. They have a data-driven method for selecting $t$. (This approach is only one part of the paper. They also consider persistent homology.)

They make two key assumptions. The first is that there exists $A<B$ such that $\{p>t\}$ is homotopic to $S$ for all $A<t<B$. (If two sets are homotopic, then they have the same homology.) This is a very reasonable assumption. In the mixture model $P=\pi R+(1-\pi)(Q\star\Phi_{\sigma})$ this assumption will be satisfied when $S$ is a small set and when $\pi$ and $\sigma$ are small. In this case, persistent homology will also work well: the dominant features in the persistence diagram will correspond to the homology of $S$.

Bobrowski et al (2104) make an additional assumption. They assume that the dimension $k$ of $S$ is known and that the rank of the $k^{\rm th}$ homology group is 0 for all $t>B$. This assumption is critical for their approach to choosing a single level set. Currently, it is not clear how strong this assumption is. In future work, we plan to compare the robustness of the single-level approach versus persistent homology.

3 Future Work

Lastly, we would like to mention that several issues deserve future attention. In particular, the methods we discussed for choosing the tuning parameters, for mitigating boundary bias and for sharpening the data, all deserve further investigation.

In a companion paper we will show how the ideas presented in this work can be used to develop hypothesis tests for comparing point clouds.

Acknowledgements

The authors are grateful to Jérome Dedecker for pointing out the key decomposition (18) of the DTM. The authors also would like to thank Jessi Cisewski and Jisu Kim for their comments.

References