Skeletonisation Algorithms with Theoretical Guarantees for Unorganised Point Clouds with High Levels of Noise

TL;DR

This paper evaluates algorithms for skeletonising noisy, unorganised point clouds, providing theoretical guarantees for shape reconstruction, including a universal method called HoPeS that adapts across scales.

Contribution

It introduces a homologically persistent skeleton (HoPeS) algorithm with guarantees for accurate shape reconstruction from noisy point clouds, without requiring extra parameters.

Findings

HoPeS is effective across various noise levels.

Algorithms can reconstruct underlying shapes with correct homotopy type.

Maximum noise levels for successful reconstruction are identified.

Abstract

Data Science aims to extract meaningful knowledge from unorganised data. Real datasets usually come in the form of a cloud of points with only pairwise distances. Numerous applications require to visualise an overall shape of a noisy cloud of points sampled from a non-linear object that is more complicated than a union of disjoint clusters. The skeletonisation problem in its hardest form is to find a 1-dimensional skeleton that correctly represents a shape of the cloud. This paper compares several algorithms that solve the above skeletonisation problem for any point cloud and guarantee a successful reconstruction. For example, given a highly noisy point sample of an unknown underlying graph, a reconstructed skeleton should be geometrically close and homotopy equivalent to (has the same number of independent cycles as) the underlying graph. One of these algorithm produces a Homologically Persistent Skeleton (HoPeS) for any cloud without extra parameters. This universal skeleton contains sub-graphs that provably represent the 1-dimensional shape of the cloud at any scale. Other subgraphs of HoPeS reconstruct an unknown graph from its noisy point sample with a correct homotopy type and within a small offset of the sample. The extensive experiments on synthetic and real data reveal for the first time the maximum level of noise that allows successful graph reconstructions.

Equations

$$f:H_1((\text{MST}(C;\alpha )\cup E(\alpha ))/\text{MST}(C;\alpha ))\to H_1(Q(C;\alpha )/\text{MST}(C;\alpha ))$$

$$(\text{MST}(C;\alpha )\cup E(\alpha ))/\text{MST}(C;\alpha )\to Q(C;\alpha )/\text{MST}(C;\alpha ).$$

$${\Delta }^k=\{(t_1,\dots ,t_k)\in {\mathbb{R}}^k\mid t_1+\cdots +t_k\le 1,t_i\ge 0\}\subset {\mathbb{R}}^k.$$

Full text

This paper compares different algorithms that solve the above skeletonisation problem for any point cloud and guarantee a successful reconstruction. For example, given a highly noisy point sample of an unknown underlying graph, a reconstructed skeleton should be geometrically close and homotopy equivalent to (has the same number of independent cycles as) the underlying graph.

One of these algorithms produces a Homologically Persistent Skeleton (HoPeS) for any cloud without extra parameters. This universal skeleton contains subgraphs that provably represent the 1-dimensional shape of the cloud at any scale. Other subgraphs of HoPeS reconstruct an unknown graph from its noisy point sample with a correct homotopy type and within a small offset of the sample. The extensive experiments on synthetic and real data reveal for the first time the maximum level of noise that allows successful graph reconstructions.

††journal: Pattern Recognition

1 Introduction: Motivations, Problem Statement and Contributions

The central problem in Data Science is to represent unorganised data in a simple and meaningful form. Our input data will be any cloud of points given by pairwise distances (a distance matrix) or by a smaller distance graph that spans all given points and specifies distances between endpoints of any edge.

Real data rarely splits into well-defined clusters. A meaningful representation of such data is a 1D skeleton consisting of straight-line edges that should provably approximate a given cloud. These 1D skeletons are relevant for curve recognition for surfaces [1], topological shapes of micelles [2] and posture identification [3].

Branches of a skeleton may reveal new classes of data points. This paper compares the Mapper algorithm with the most relevant methods ($\alpha$-Reeb [4] and HoPeS [5]), because they start from the same input (an unorganised point cloud) and solve the problem below with topological and geometric guarantees. All required concepts such as homology are introduced in Appendix A.

Data skeletonisation problem. Given a noisy cloud $C$ of points sampled from a graph $G$ in a metric space, find conditions on $G$ and $C$ such that a reconstructed graph $G^{\prime}$ has the same homology $H_{1}(G^{\prime})=H_{1}(G)$ and geometrically approximates $G$ in the sense that $G^{\prime}\subset G^{\alpha}$ and $G\subset(G^{\prime})^{\alpha}$ for a suitable parameter $\alpha$ depending on $G$ and $C$, see $\alpha$-offsets defined in Appendix A.

Informally, the homology condition above implies that the reconstructed graph $G^{\prime}$ can be continuously deformed to the original graph $G$. The first two graphs in Fig. 2 have the same homology $H_{1}$, because the pair of adjacent edges at each corner vertex in the grid graph $G(2,2)$ can be merged and continuously deformed into one corresponding diagonal edge in the wheel graph $W(4)$. Geometric similarity means that $G$ and $G^{\prime}$ are close to each other with respect to a distance between graphs, e.g. $G^{\prime}$ is in a small neighbourhood of $G$ and vice versa.

The contributions of the paper to data skeletonisation are the following:

$\bullet$ Section 4 gives a much simpler proof of Optimality Theorem 4.8 than for the higher-dimensional version of a Homologically Persistent Skeleton (HoPeS) by Kalisnik et al. [6]. This result shows that $\mathrm{HoPeS}$ extends a classical Minimum Spanning Tree to an optimal graph with cycles representing the 1-dimensional shape of a point cloud $C$ by correctly capturing ‘cycles’ of $C$ at any scale $\alpha$. Briefly, the reduced skeleton $\mathrm{HoPeS}(C,\alpha)$ has the minimum length among all graphs that span $C$ and have the same $H_{0}$ and $H_{1}$ homology as the offset $C^{\alpha}$.

$\bullet$ Key Stability Theorem 3.6 of Topological Data Analysis is extended to Graph Reconstruction Theorems 5.4 and 5.9, which are proved in section 5 after being announced in [5]. Corollary 5.10 proves a global stability of derived subgraphs of HoPeS for the first time, which justifies its applicability to noisy data [7, 8].

$\bullet$ Different algorithms are extensively compared in section 8 by their runtime and topological and geometric measures on the same datasets of real and synthetic clouds with high noise whose generation is explained in section 6. Figs. 23 and 25 show the maximum noise parameters that allow correct reconstructions.

2 Review of Related Past Work on Data Skeletonisation Algorithms

Among many skeletonisation algorithms we review the most relevant ones that can accept any point cloud as in Definition A.1 and provide guarantees with respect to the skeletonisation problem in section 1, see Appendix A.

Past algorithms without guarantees. Singh et al. [9] approximated a cloud $C\subset\mathbb{R}^{d}$ by a subgraph of a Delaunay triangulation, which requires $O(n^{\lceil{d/2}\rceil})$ time for $n$ points of $C$ and three thresholds: a minimum number $K$ of edges in a cycle and $\delta_{min},\delta_{max}$ for inserting/merging 2nd order Voronoi regions. Similar parameters are needed forprincipal curves [10], which were later extended to iteratively computed elastic maps [11]. Since it is hard to estimate a rate of convergence for iterative algorithms, we discuss below non-iterative methods.

Skeletonisation via Reeb graphs. Starting from a noisy sample $C$ of an unknown graph $G$ with a scale parameter $\alpha$, X. Ge et al. [12] considered the Reeb graph of the Vietoris-Rips complex on $C$ at a scale $\alpha$, see Definition A.5.

The $\alpha$-Reeb graph $G$ was introduced by F. Chazal et al. [4] for a metric space $X$ at a user-defined scale $\alpha$, see Definition B.1. If $X$ is $\epsilon$-close to an unknown graph with edges of minimum length $8\epsilon$, the output $G$ is $34(\beta(G)+1)\epsilon$-close to the input $X$, where $\beta(G)$ is the first Betti number of $G$, see [4, Theorem 4.9]. The $\alpha$-Reeb graph has a metric, but isn’t embedded into any space even if $C\subset\mathbb{R}^{2}$. The algorithm has the fast time $O(n\log n)$ for $n$ input points of $X$.

Mapper [13] extends any clustering and outputs a network of interlinked clusters by using a user-defined filter function $f:C\to\mathbb{R}$, which helps to link different clusters of a cloud $C$. M. Carriére et al. [14, Theorem 5.2] found a connection between Mapper and the Reeb graph via MultiNerve Mapper.

Metric graph reconstruction. M. Aanjaneya et al. [15] studied a related problem approximating a metric on a large input graph $Y$ by a metric on a small output graph $\hat{X}$. Let $Y$ be a good $\epsilon$-approximation to an unknown graph $X$. Then [15, Theorem 2] is the first guarantee for the existence of a homeomorphism $X\to\hat{X}$ that distorts the metrics on $X$ and $\hat{X}$ with a multiplicative factor $1+c\epsilon$ for $c>\frac{30}{b}$, where $b>14.5\epsilon$ is the length of the shortest edge of $X$.

Graph reconstruction by discrete Morse theory. A Homological Spanning Forest [16] uses a given pixel grid of 2D images, hence cannot be applied to an arbitrary point cloud $C$. Similarly, the recent algorithm by T. Dey et al. [17] in addition to a cloud $C$ requires a density field $\rho:\mathbb{R}^{d}\to\mathbb{R}$ (usually on a regular grid), which ‘concentrates’ around a hidden graph. Section 8 compares the three algorithms that accept an input cloud $C$ without any extra structure.

3 $\mathrm{HoPeS}(C)$: a Homologically Persistent Skeleton of a cloud $C$

The Homologically Persistent Skeleton (HoPeS) was introduced in 2015 [5] and extended to higher dimensions by S. Kalisnik et al. [6, Definition 4.5], which has crucially clarified potential choices of homology representatives. The 1-dimensional HoPeS was motivated by extending a reconstruction of hole boundaries in unorganised clouds of edge pixels in images to a 1D skeleton [18, 19].

3.1 A Minimum Spanning Tree of any point cloud in a metric space

Definition 3.1** (MST).**

Fix a filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$. A Minimum Spanning Tree $\text{MST}(C)$ is a connected graph with vertex set $C$ and the minimum total length of edges. The length of an edge $e$ is the minimum $\alpha$ such that $e\subset Q(C;\alpha)$. For any scale $\alpha\geq 0$, a forest $\text{MST}(C;\alpha)$ is obtained from $\text{MST}(C)$ by removing interiors of all edges that are longer than $2\alpha$. $\blacksquare$

$\text{MST}(C)$ may not be unique if different edges enter the filtration $\{Q(C;\alpha)\}$ at the same scale $\alpha$, e.g. $\text{MST}(C)$ in Fig. 4 may include the upper horizontal edge instead of the lower one. For all these choices, $\text{MST}(C;\alpha)$ enjoys the optimality in Lemma 3.2. A graph $G$ spans a possibly disconnected complex $Q$ on a cloud $C$ if $G$ has the vertex set $C$, every edge of $G$ belongs to $Q$ and the inclusion $G\subset Q$ induces a 1-1 correspondence between connected components.

Lemma 3.2**.**

Fix a filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$ in a metric space. For any fixed scale $\alpha\geq 0$, the forest $\text{MST}(C;\alpha)$ has the minimum total length of edges among all graphs that span $Q(C;\alpha)$ at the same scale $\alpha$.

Proof.

Let $e_{1},\dots,e_{m}\subset\text{MST}(C)$ be all edges longer than $2\alpha$, so $\text{MST}(C)=e_{1}\cup\dots\cup e_{m}\cup\text{MST}(C;\alpha)$. Assume that there is a graph $G$ that spans $Q(C;\alpha)$ and is shorter than $\text{MST}(C;\alpha)$. Since both inclusions $G\subset Q(C;\alpha)$ and $\text{MST}(C;\alpha)\subset Q(C;\alpha)$ induce 1-1 maps on connected components, the graph $G\cup e_{1}\cup\dots\cup e_{m}$ is connected and is shorter than $\text{MST}(C)$, which contradicts Definition 3.1. ∎

Lemma 3.2 implies that, for any scale $\alpha$, the complex $Q(C;\alpha)$ can be replaced by the smaller graph $\text{MST}(C;\alpha)$ that has the same connected components and the minimal total length among all graphs with this connectivity property.

Visually, the 0-dimensional shape (connected components) of each offset $C^{\alpha}$ on a cloud $C\subset\mathbb{R}^{d}$ can be represented by the forest $\text{MST}(C;\alpha)$ instead of the larger $\alpha$-complex $C(\alpha)$ homotopically equivalent to $C^{\alpha}$, see Lemma A.6.

3.2 Persistent homology and its stability under perturbations

A Homologically Persistent Skeleton ($\mathrm{HoPeS}$) will extend the optimality of MST in Lemma 3.2 to $H_{1}$ in Theorem 4.8. $\mathrm{HoPeS}$ is based on the persistent homology that summarises the evolution of homology as formalised below.

Definition 3.3** (births and deaths).**

For any filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$ in a metric space, a homology class $\gamma\in H_{1}(Q(C;\alpha_{i}))$ is born at $\alpha_{i}=\text{birth}(\gamma)$ if $\gamma$ is not in the full image under the induced homomorphism $H_{1}(Q(C;\alpha))\to H_{1}(Q(C;\alpha_{i}))$ for any $\alpha<\alpha_{i}$. The class $\gamma$ dies at $\alpha_{j}=\text{death}(\gamma)\geq\alpha_{i}$ when the image of $\gamma$ under $H_{1}(Q(C;\alpha_{i}))\to H_{1}(Q(C;\alpha_{j}))$ merges into the image under $H_{1}(Q(C;\alpha))\to H_{1}(Q(C;\alpha_{j}))$ for some $\alpha<\alpha_{i}$. $\blacksquare$

The birth-death pairs from Definition 3.3 can be similarly defined for higher-dimensional homology groups $H_{k}$ with $k>1$ and coefficients in a field, though coefficients in $\mathbb{Z}_{2}=\{0,1\}$ are used in practice and in our paper. For $k=0$ and the filtration of offsets $C^{\alpha}$, the homology group $H_{0}(C^{\alpha})$ is generated by connected components of $C^{\alpha}$. Then all homology classes of $H_{0}(C^{\alpha})$ are born at $\alpha=0$ (isolated dots) and die at $\alpha$ equal to half-lengths of edges in $\text{MST}(C)$.

Definition 3.4** (persistence diagram).**

Fix a filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$. Let $\alpha_{1},\dots,\alpha_{m}$ be all scales when a homology class is born or dies in $H_{1}(Q(C;\alpha))$. Let $\mu_{ij}$ be the number of independent classes in $H_{1}(Q(C;\alpha))$ that are born at $\alpha_{i}$ and die at $\alpha_{j}$. The persistence diagram $\text{PD}\{Q(C;\alpha)\}\subset\mathbb{R}^{2}$ is the multi-set consisting of the birth-death pairs $(\alpha_{i},\alpha_{j})$ with integer multiplicities $\mu_{ij}\geq 1$ and all diagonal dots $(x,x)$ with infinite multiplicity. $\blacksquare$

For the cloud $C$ in Fig. 3 and the filtration $\{C^{\alpha}\}$, the homology $H_{1}$ has 2 classes that persist over $1.5\leq\alpha<R$ and $2\leq\alpha<R$, where $R=\frac{15}{8}\sqrt{17}$ is the circumradius of the largest Delaunay triangle in $\mathrm{Del}(C)$ with sides $4,\sqrt{17},5$. We use the value $R\approx 2.577$. Hence the persistence diagram $\text{PD}\{C^{\alpha}\}$ has 2 off-diagonal red dots $(1.5,2.577)$ and $(2,2.577)$ in the last picture of Fig. 3.

The persistence diagram $\text{PD}\{C^{\alpha}\}$ is invariant under any isometry of $C$, i.e. any transformation that preserves distances between points. The key advantage of $\text{PD}\{C^{\alpha}\}$ over other geometric invariants is Stability Theorem 3.6, which requires the bottleneck distance between persistence diagrams defined below.

Definition 3.5** (bottleneck distance $d_{B}$).**

For points $p=(x_{1},y_{1})$, $q=(x_{2},y_{2})$ in $\mathbb{R}^{2}$, the $L_{\infty}$-distance is $||p-q||_{\infty}=\max\{|x_{1}-x_{2}|,|y_{1}-y_{2}|\}$. The bottleneck distance between persistence diagrams $\text{PD},\text{PD}^{\prime}$ is defined as $d_{B}=\inf\limits_{\psi}\sup\limits_{q\in\text{PD}}||q-\psi(q)||_{\infty}$ for all bijections $\psi:\text{PD}\to\text{PD}^{\prime}$ between multi-sets. $\blacksquare$

Stability Theorem 3.6 below informally says that any small perturbation of original data leads to a small perturbation of the persistence diagram. A metric space $M$ is totally bounded if $M$ has a finite $\epsilon$-sample $C\subset M$ for any $\epsilon>0$.

Theorem 3.6** (stability of persistence).**

[20, simplified Theorem 5.6] Let $C$ be any $\epsilon$-sample of a graph $G$ in a totally bounded metric space $M$. Then the persistence diagrams of C̆ech filtrations on $G$ and $C$ are $\epsilon$-close, namely $d_{B}(\text{PD}\{\text{\emph{\u{C}h}}(G,M;\alpha)\},\text{PD}\{\text{\emph{\u{C}h}}(C,M;\alpha)\})\leq\epsilon$. This inequality holds for the filtrations of $\alpha$-offsets by Nerve Lemma A.6, namely $d_{B}(\text{PD}\{G^{\alpha}\},\text{PD}\{C^{\alpha}\})\leq\epsilon$.

3.3 Gap search method for stable subdiagrams of a persistence diagram

Stability Theorem 3.6 extends to important subdiagrams with highly persistent dots representing features that persist over a long interval of a scale.

Definition 3.7** (diagonal gaps and subdiagrams).**

For any metric space $C$, a diagonal gap of $\text{PD}\{C^{\alpha}\}$ is a strip $\{0\leq a<y-x<b\}$ that has dots in both boundary lines $\{y-x=a\}$ and $\{y-x=b\}$, but not inside the strip. For any $k\geq 1$, the $k$-th widest diagonal gap $\text{dgap}_{k}(C)$ has the $k$-th largest vertical width $|\text{dgap}_{k}(C)|=b-a$. The subdiagram $\text{DS}_{k}(C)\subset\text{PD}\{C^{\alpha}\}$ consists of only the dots above the lowest of the first $k$ widest $\text{dgap}_{i}(C)$, $i=1,\dots,k$. Each $\text{DS}_{k}(C)$ is bounded below by $y-x=a$ and has the diagonal scale $\text{ds}_{k}(C)=a$. $\blacksquare$

In Definition 3.7, if $\text{PD}\{C^{\alpha}\}$ has different gaps of the same width, we say that a lower diagonal gap has a larger width. If $\text{PD}\{C^{\alpha}\}$ has dots only in $m$ different lines $\{y-x=a_{i}>0\}$, $i=1,\dots,m$, we have $m$ diagonal gaps $\text{dgap}_{i}(C)$. We ignore the highest gap $\{y-x>\max a_{i}\}$, so we set $\text{dgap}_{i}(C)=\emptyset$ for $i>m$.

The cloud $C$ in Fig. 3 has the persistence diagram $\text{PD}\{C^{\alpha}\}$ with only 2 off-diagonal dots $(1.5,2.577)$ and $(2,2.577)$. Then $\text{dgap}_{1}(C)=\{0<y-x<0.577\}$ is below $\text{dgap}_{2}(C)=\{0.577<y-x<1.577\}$. So $\text{DS}_{1}(C)=\text{DS}_{2}(C)$ consists of $(1.5,2.577),(2,2.577)$, $\text{ds}_{1}(C)=\text{ds}_{2}(C)=0.577$ in Fig. 6. Definition 3.7 makes sense in any persistence diagram, say for $\alpha$-offsets $G^{\alpha}$ of a graph $G\subset M$. The graph $G$ in Fig. 7 has the 1D diagram $\text{PD}\{G^{\alpha}\}$ containing only one off-diagonal dot $(0,2.577)$ of multiplicity 2. Then $\text{dgap}_{1}(G)=\{0<y-x<2.577\}$, $|\text{dgap}_{2}(G)|=0$, $\text{DS}_{1}(G)=\{(0,2.577),(0,2.577)\}$ and $\text{ds}_{1}(G)=2.577$.

If we need to reconstruct cycles of $G$ from a noisy sample $C$ using the perturbed diagram $\text{PD}\{C^{\alpha}\}$, we should look for birth-death pairs having a high persistence $\text{death}-\text{birth}$ and small birth. Hence we search for a vertical gap to separate the dots $(\text{birth},\text{death})$ with a small birth from all other dots.

Definition 3.8** (vertical gap $\text{vgap}_{k,l}$, subdiagram $\text{VS}_{k,l}$, scale $\text{vs}_{k,l}$).**

In the diagonal subdiagram $\text{DS}_{k}(C)$ from Definition 3.7, a vertical gap is the widest vertical strip $\{a<x<b\}$ whose boundary contains a dot from $\text{DS}_{k}(C)$ in the line $\{x=a\}$, but not inside the strip. For $l\geq 1$, the $l$-th widest vertical gap $\text{vgap}_{k,l}(C)$ has the $l$-th largest horizontal width $|\text{vgap}_{k,l}(C)|=b-a$. The vertical subdiagram $\text{VS}_{k,l}(C)\subset\text{DS}_{k}(C)$ consists of only the dots to the left of the first $l$ widest vertical gaps $\text{vgap}_{k,j}(C)$, $j=1,\dots,l$. So each $\text{VS}_{k,l}(C)$ is bounded on the right by $x=b$ and has the vertical scale $\text{vs}_{k,l}(C)=a$. $\blacksquare$

In Definition 3.8, if there are different vertical gaps with the same horizontal width, we say that the leftmost vertical gap has a larger width. We prefer the leftmost vertical gap, while in Definition 3.7 we prefer the lowest diagonal gap.

We allow the case $b=+\infty$, so the widest vertical gap $\text{vgap}_{k,1}(C)$ always has the form $\{x>a\}$, $VS_{k,1}(C)=\text{DS}_{k}(C)$ and we set $|\text{vgap}_{k,1}(C)|=+\infty$.

For the cloud $C$ in Fig. 3, the diagonal subdiagram $\text{DS}_{1}(C)$ has $\text{vgap}_{1,1}(C)=\{x>2\}$, $\text{vgap}_{1,2}(C)=\{1.5<x<2\}$, $\text{VS}_{1,1}(C)=\{(1.5,2.577),(2,2.577)\}$ and $\text{VS}_{1,2}(C)=\{(1.5,2.577)\}$, so $\text{vs}_{1,1}(C)=2$, $\text{vs}_{1,2}(C)=1.5$, see Fig. 6. Any cloud $C$ has no cycles at $\alpha=0$, hence $\text{PD}\{C^{\alpha}\}$ has no dots in the vertical axis, so any $\text{vgap}_{k,l}(C)$ has a left boundary line $\{x=a>0\}$ and all $\text{vs}_{k,l}(C)>0$.

Definition 3.8 makes sense for any persistence diagram, say for $\alpha$-offsets $G^{\alpha}$ of a graph $G$ in a metric space. The diagonal subdiagram $\text{DS}_{1}(\theta)$ for the graph $\theta$ in Fig. 7 consists of the doubled dot $(0,2.577)$. The only vertical subdiagram $\text{VS}_{1,1}(\theta)$ consists of the same doubled dot $(0,2.577)$, so $\text{vs}_{1,1}(\theta)=0$.

3.4 $\mathrm{HoPeS}(C)$ is the persistence-based extension of $\text{MST}(C)$ for a cloud $C$

$\mathrm{HoPeS}(C)$ will be obtained from $\text{MST}(C)$ by adding critical edges below.

Definition 3.9** (critical edges).**

[6, modified from Definition 4.5] Fix a filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$. An edge $e$ is critical if the addition of $e$ to the filtration $\{Q(C;\alpha)\}$ creates a new homology class that does not immediately die, see two red critical edges in Fig. 5. The birth of a critical edge $e$ is the scale $\alpha$ when $e$ enters the filtration. Define $E(\alpha)=\{e_{1},...,e_{s}\}$ to be the set of critical edges with $\text{birth}(e_{i})\leq\alpha$ that have not yet been assigned a death value. Then $[e_{1}],...,[e_{s}]$ form a basis of $H_{1}((\text{MST}(C;\alpha)\cup E(\alpha))\>/\>\text{MST}(C;\alpha))$. Denote

[TABLE]

to be the map on homology induced by the inclusion

[TABLE]

Let $\{b_{1},...,b_{r}\}$ be a basis of $\text{ker}(f)$, where $r=\text{dim ker}(f)$ is equal to the number of critical edges that die at exactly scale $\alpha$, so we have $r\leq s$. We can expand each basis element as $b_{i}=\sum_{j=1}^{s}c_{ij}[e_{j}]$, with $c_{ij}\in\mathbb{Z}_{2}$ (coefficients for $H_{1}$). We consider (in $\mathbb{Z}_{2}$) the system of equations $\sum_{j=1}^{s}c_{ij}x_{j}=0$ for $i\in\{1,...,r\}$. Since basis elements are linearly independent, so are these equations. Then this system of equations can be solved with $r$ leading variables, each of which is expressed in terms of the remaining $s-r$ free variables. Let $\mathit{I}\subseteq\{1,...,s\}$ be the set of all indices of leading variables ($\mathit{I}$ could be empty if $f$ is trivial). For each $i\in\mathit{I}$, we declare the death value of the critical edge $e_{i}$ to be $\alpha$. $\blacksquare$

The elder rule selects those leading variables whose corresponding set of critical edges has the greatest combined birth value. If there are distinct sets with the greatest combined birth value, then we have a genuine choice. However, Optimality Theorem 4.8 will be valid regardless of the choice.

Definition 3.10** (HoPeS).**

For any filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$, a Homologically Persistent Skeleton $\mathrm{HoPeS}(C)$ is the union of a minimum spanning tree $\text{MST}(C)$ and all critical edges with labels $(\text{birth},\text{death})$. For any $\alpha\geq 0$, the reduced skeleton $\mathrm{HoPeS}(C;\alpha)$ is obtained from $\mathrm{HoPeS}(C)$ by removing all edges longer than $2\alpha$ and all critical edges $e$ with $\text{death}(e)\leq\alpha$. $\blacksquare$

4 Optimality guarantees for reduced subgraphs of $\mathrm{HoPeS}(C)$

Optimality Theorem 4.8 will prove that, for any scale $\alpha$, the reduced subgraphs $\mathrm{HoPeS}(C;\alpha)$ are minimal graphs that span a cloud $C$ and have the same components and cycles as the offset $C^{\alpha}$. We fix a filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$, but use the simpler notation $\mathrm{HoPeS}(C)$ for the skeleton.

Lemma 4.1**.**

Let a class $\gamma\in H_{1}(Q(C;\alpha))$ be born due to a critical edge $e$ added to $Q(C;\alpha)$. Then $2\text{birth}(\gamma)$ is the length $|e|$ relative to the filtration $\{Q(C;\alpha)\}$.

Proof.

By Definition 3.9 the critical edge $e$ is the last edge added to a cycle $L\subset Q(C;\alpha)$ giving birth to the homology class $\gamma$ at $\alpha=\text{birth}(\gamma)$. The length $|e|$ equals the doubled scale $2\alpha$ when $e$ enters $Q(C;\alpha)$, so $|e|=2\cdot\text{birth}(\gamma)$. ∎

Lemma 4.2**.**

The reduced skeleton $\mathrm{HoPeS}(C;\alpha)$ is a subgraph of the complex $Q(C;\alpha)$ for any $\alpha$, e.g. $\mathrm{HoPeS}(C;\alpha)\subset\mathrm{Del}(C)$ for the filtration of $\alpha$-complexes.

Proof.

By Definition 3.10, the skeleton $\mathrm{HoPeS}(C;\alpha)$ for any $\alpha$ consists of the forest $\text{MST}(C;\alpha)$ and all critical edges with lengths $\text{death}(e)\leq|e|\leq 2\alpha$. Lemma 3.2 implies that $\text{MST}(C;\alpha)\subset Q(C;\alpha)$. Any critical edge $e\subset\mathrm{HoPeS}(C)$ belongs to $Q(C;\alpha)$ for $2\alpha\geq|e|$. Hence $\mathrm{HoPeS}(C;\alpha)\subset Q(C;\alpha)$. ∎

The inclusion $\mathrm{HoPeS}(C;\alpha)\subset Q(C;\alpha)$ in Lemma 4.2 induces a homomorphism $f_{*}$ in $H_{1}$. Lemmas 4.3 and 4.4 analyse what happens with $f_{*}$ when a critical edge $e$ is added to (or deleted from) $Q=Q(C;\alpha)$ and $G=\mathrm{HoPeS}(C;\alpha)$.

Lemma 4.3** (addition of a critical edge).**

Let an inclusion $f:G\to Q$ of a graph $G$ into a simplicial complex $Q$ induce an isomorphism $f_{*}:H_{1}(G)\to H_{1}(Q)$. Between some vertices $u,v\in G$ let us add an edge $e$ to $G$ and $Q$ that creates $\gamma\in H_{1}(Q\cup e)$. Then $f_{*}$ extends to an isomorphism $H_{1}(G\cup e)\to H_{1}(Q\cup e)$.

Proof.

Let $L\subset G\cup e$ be a cycle containing the added edge $e$. Then $f$ extends to the inclusion $G\cup e\to Q\cup e$ and induces the isomorphism $H_{1}(G\cup e)\cong H_{1}(G)\oplus\langle[L]\rangle$. Now $f(L)$ is a cycle in $Q\cup e$ and $H_{1}(Q\cup e)\cong H_{1}(Q)\oplus\langle[f(L)]\rangle$. Mapping $[L]$ to $[f(L)]\in H_{1}(Q\cup e)$, we extend $f_{*}$ to a required isomorphism $H_{1}(G)\oplus\langle[L]\rangle\to H_{1}(Q)\oplus\langle[f(L)]\rangle$. ∎

Lemma 4.4** (deletion of a critical edge).**

Let an inclusion $f:G\to Q$ of a graph $G$ into a simplicial complex $Q$ induce an isomorphism $f_{*}:H_{1}(G)\to H_{1}(Q)$. Let $\gamma\in H_{1}(Q)$ die after adding a triangle $T$ to $Q$. Let $e$ be a longest open edge of the triangle $T$. Then $f_{*}$ descends to an isomorphism $H_{1}(G-e)\to H_{1}(Q\cup T)$.

Proof.

Adding $T$ to $Q$ kills the homology class $\partial T$ of the boundary $\partial T\subset G$. Then $H_{1}(Q\cup T)\cong H_{1}(Q)/\langle[\partial T]\rangle$. Deleting the open edge $e$ from the boundary $\partial T\subset G$ makes the homology smaller: $H_{1}(G-e)\cong H_{1}(G)/\langle[\partial T]\rangle$. The isomorphism $f_{*}$ descends to the isomorphism $H_{1}(G)/\langle[\partial T]\rangle\to H_{1}(Q)/\langle[\partial T]\rangle$. ∎

Proposition 4.5**.**

The inclusion $\mathrm{HoPeS}(C;\alpha)\to Q(C;\alpha)$ from Lemma 4.2 induces an isomorphism of 1D homology groups: $H_{1}(\mathrm{HoPeS}(C;\alpha))\to H_{1}(Q(C;\alpha))$.

Proof.

$\mathrm{HoPeS}(C;0)$ and the complex $Q(C;0)$ coincide with the cloud $C$, so their 1-dimensional homology groups are trivial. Each time a homology class is born or dies in $H_{1}(Q(C;\alpha))$, the isomorphism in $H_{1}$ induced by the inclusion $\mathrm{HoPeS}(C;\alpha)\to Q(C;\alpha)$ is preserved by Lemmas 4.3 and 4.4. ∎

Lemma 4.6**.**

Let a cycle $L\subset Q(C;\alpha)$ represent a homology class $\gamma\in H_{1}(Q(C;\alpha))$ in $\text{PD}\{Q(C;\alpha)\}$. Then any longest edge $e\subset L$ has length $|e|\geq 2\cdot\text{birth}(\gamma)$.

Proof.

Let a longest edge $e$ of a cycle $L\subset Q(C;\alpha)$ representing the class $\gamma$ have a half-length $\alpha<\text{birth}(\gamma)$. Then $L$ enters $Q(C;\alpha)$ earlier than $\text{birth}(\gamma)$ and can not represent the class $\gamma$ that starts living from $\alpha=\text{birth}(\gamma)$. ∎

Recall that the forest $\text{MST}(C;\alpha)$ on a cloud $C$ at a scale $\alpha$ is obtained from a minimum spanning tree $\text{MST}(C)$ by removing all open edges longer than $2\alpha$. An edge $e$ is splitting a graph $G$ if removing the open edge $e$ makes $G$ disconnected. Otherwise the edge $e$ is called non-splitting and belongs to a cycle of $G$.

Proposition 4.7**.**

Let a graph $G\subset Q(C;\alpha)$ span $Q(C;\alpha)$ and $H_{1}(G)\to H_{1}(Q(C;\alpha))$ be an isomorphism induced by the inclusion. Let $(b_{i},d_{i})$, $i=1,\dots,m$, be all dots in the persistence diagram $\text{PD}\{Q(C;\alpha)\}$ such that $\text{birth}\leq\alpha<\text{death}$. Then the total length of $G$ is bounded below by the total length of edges of the forest $\text{MST}(C;\alpha)$ plus $2\sum\limits_{i=1}^{m}b_{i}$.

Proof.

Let the subgraph $G_{1}\subset G$ consist of all non-splitting edges of $G$ and $e_{1}\subset G_{1}$ be a longest open edge. Removing $e_{1}$ from $G$ makes $H_{1}(G)$ smaller. Hence there is a cycle $L_{1}\subset G_{1}$ containing $e_{1}$ and representing a class $\gamma_{1}\in H_{1}(Q(C;\alpha))$ that corresponds to some off-diagonal dot in $\text{PD}\{Q(C;\alpha)\}$, say $(b_{1},d_{1})$. Then $\gamma_{1}$ lives over the interval $b_{1}=\text{birth}(\gamma_{1})\leq\alpha<d_{1}=\text{death}(\gamma_{1})$. Lemma 4.6 implies that $|e_{1}|\geq 2b_{1}$. Let the graph $G_{2}\subset G-e_{1}$ consist of all non-splitting edges and $e_{2}\subset G_{2}$ be a longest open edge. Similarly, find the corresponding point $(b_{2},d_{2})$, conclude that $|e_{2}|\geq 2b_{2}$ and so on until we get $\sum\limits_{i=1}^{m}|e_{i}|\geq 2\sum\limits_{i=1}^{m}b_{i}$. After removing all open edges $e_{1},\dots,e_{m}$, the remaining graph $G-(e_{1}\cup\dots\cup e_{m})$ still spans the (possibly disconnected) complex $Q(\alpha)$. Indeed, each time we removed a non-splitting edge. So the total length of $G-(e_{1}\cup\dots\cup e_{m})$ is not smaller than the total length of $\text{MST}(C;\alpha)$ by Lemma 3.2. ∎

Theorem 4.8** (optimality of reduced skeletons).**

Fix a filtration $\{Q(C;\alpha)\}$ of complexes on a cloud $C$ in a metric space. For any $\alpha>0$, the reduced skeleton $\mathrm{HoPeS}(C;\alpha)$ has the minimum total length of edges over all graphs $G\subset Q(C;\alpha)$ that span $Q(C;\alpha)$ and induce an isomorphism $H_{1}(G)\to H_{1}(Q(C;\alpha))$.

Proof.

For any $\alpha>0$, the inclusion $\mathrm{HoPeS}(C;\alpha)\to Q(C;\alpha)$ induces an isomorphism in $H_{1}$ by Proposition 4.5. Let classes $\gamma_{1},\dots,\gamma_{m}$ represent all $m$ dots counted with multiplicities in the ‘rectangle’ $\{\text{birth}\leq\alpha<\text{death}\}\subset\text{PD}\{Q(C;\alpha)\}$. Then $\gamma_{1},\dots,\gamma_{m}$ form a basis of $H_{1}(Q(C;\alpha))\cong H_{1}(\mathrm{HoPeS}(C;\alpha))$ by Definition 3.4. The total length of $\mathrm{HoPeS}(C;\alpha)$ equals the length of $\text{MST}(C;\alpha)$ plus $2\sum\limits_{i=1}^{m}\text{birth}(\gamma_{i})$ by Lemma 4.1. By Proposition 4.7 this length is minimal over all graphs $G\subset Q(C;\alpha)$ that span $Q(C;\alpha)$ and have the same $H_{1}$ as $Q(C;\alpha)$. ∎

5 Guarantees for graph reconstructions using derived $\mathrm{HoPeS}$

This section shows that other subgraphs (derived skeletons) of HoPeS in Definition 5.1 guarantee correct reconstructions in Theorems 5.4, 5.9, see Fig. 8.

Definition 5.1** (derived skeletons $\mathrm{HoPeS}_{k,l}$).**

Let $C$ be a finite cloud in a metric space. For $k,l\geq 1$, the derived skeleton $\mathrm{HoPeS}_{k,l}(C)$ is obtained from $\mathrm{HoPeS}(C)$ by removing all edges that are longer than $2\text{vs}_{k,l}(C)$ and by keeping only critical edges with $(\text{birth},\text{death})\in\text{VS}_{k,l}(C)$ and with $\text{death}>\text{vs}_{k,l}(C)$, see Fig. 8. $\blacksquare$

Definition 5.2** (radius $\rho$ of a cycle and thickness $\theta$ of a graph).**

For a graph $G$ in a metric space $M$, the radius $\rho$ of a non-self-intersecting cycle $L\subset G^{\alpha}$ is the persistence of its corresponding homology class in the filtration, see Fig. 9. When $\alpha$ is increasing, new holes can be born in $G^{\alpha}$. The thickness $\theta(G)$ is the maximum persistence of any hole born after $\alpha=0$ in the diagram $\text{PD}\{G^{\alpha}\}$. $\blacksquare$

If a cycle $L\subset\mathbb{R}^{2}$ encloses a convex region, the only hole of $L^{\alpha}$ completely dies when $\alpha$ equals the radius $\rho(L)$, so $\theta(L)=0$. The ‘heart’ cycle $L$ in Fig. 9 encloses a non-convex region, however no new holes are born in $L^{\alpha}$, so $\theta(L)=0$. The figure-eight-shaped graph $G$ in Fig. 9 has a positive thickness $\theta(G)$ equal to the radius $\rho(L_{1})$ of the largest cycle born in $G^{\alpha}$ when $\alpha$ is increasing.

Lemma 5.3**.**

For a graph $G$ in a metric space $M$, the 1D diagram $\text{PD}\{G^{\alpha}\}$ has $m$ dots (with multiplicities) on the death axis, where $m=\text{rank}\;H_{1}(G)$.

Proof.

By Definition 3.3 all $m$ dots $(0,d_{i})$ in $\text{PD}\{G^{\alpha}\}$ correspond to homology classes born at $\alpha=0$, hence present in $H_{1}(G)$ of the initial graph $G^{0}=G$. ∎

Theorem 5.4**.**

Let $M$ be a metric space, and let $C$ be any $\epsilon$-sample of a connected graph $G\subset M$ with $\text{rank}\;H_{1}(G)=m$. Let $m$ dots on the vertical axis of $\text{PD}\{G^{\alpha}\}$ from Lemma 5.3 have ordered deaths $y_{1}\leq\cdots\leq y_{m}$. If $y_{1}>7\epsilon+2\theta(G)+\underset{i=1,\ldots,m-1}{\text{max}}\{y_{i+1}-y_{i}\}$, we get the lower bound for noise $\text{vs}_{1,1}(C)\leq\epsilon$. The derived skeleton $\mathrm{HoPeS}_{1,1}(C)\subset G^{2\epsilon}$ has the same homology $H_{1}$ as $G$.

Proof.

Apart from the dots on the $y$-axis, all other dots in the 1D persistence diagram $\text{PD}\{G^{\alpha}\}$ come from cycles born later in the $\alpha$-offsets $G^{\alpha}$. The maximum persistence of these dots is bounded above by the thickness $\theta(G)$, see Fig. 10.

The inequality $y_{1}>7\epsilon+2\theta(G)+\underset{i=1,...,m-1}{\text{max}}\{y_{i+1}-y_{i}\}$ guarantees that the gap $\{\theta(G)<y-x<y_{1}\}$ in $\text{PD}\{G^{\alpha}\}$ is wider than any other gaps including the higher gaps $y_{i+1}-y_{i}$ between the dots $(0,y_{i})\in\text{PD}\{G^{\alpha}\},i=1,\dots,m-1$, since the inequality implies $y_{1}-\theta(G)>\underset{i=1,...,m-1}{\text{max}}\{y_{i+1}-y_{i}\}$ and $y_{1}-\theta(G)>\theta(G)$.

By Stability Theorem 3.6, the perturbed diagram $\text{PD}\{C^{\alpha}\}$ is in the $\epsilon$-offset of $\text{PD}\{G^{\alpha}\}\subset\cup_{i=1}^{m}(0,y_{i})\cup\{y-x<\theta(G)\}$ with respect to the $L_{\infty}$ metric on $\mathbb{R}^{2}$. All noisy dots near the diagonal in $\text{PD}\{C^{\alpha}\}$ can not be higher that $\theta(G)+2\epsilon$ after projecting along the diagonal $\{x=y\}$ to the vertical axis $\{x=0\}$. The remaining dots can not be lower than $y_{1}-2\epsilon$ after the same projection $(x,y)\to y-x$. So the smaller diagonal strip $\{\theta(G)+2\epsilon<y-x<y_{1}-2\epsilon\}$ of the vertical width $y_{1}-4\epsilon-\theta(G)$ is empty in the perturbed diagram $\text{PD}\{C^{\alpha}\}$.

By Stability Theorem 3.6, any dot $(0,y_{i})\in\text{PD}\{G^{\alpha}\},i\geq 1$, can not jump lower than the line $y-x=y_{i}-2\epsilon$ or higher than $y-x=y_{i}+\epsilon$ (not $y-x=y_{i}+2\epsilon$ because the dot is on the $y$-axis and hence can not move in the negative $x$ direction). Then the widest diagonal gap between these perturbed dots has a vertical width at most $\underset{i=1,...,m-1}{\text{max}}\{y_{i+1}-y_{i}\}+3\epsilon$, see the last picture of Fig. 10. All dots near the diagonal have diagonal gaps not wider than $\theta(G)+2\epsilon$. The given inequality $y_{1}-4\epsilon-\theta(G)>\theta(G)+\underset{i=1,\ldots,m-1}{\text{max}}\{y_{i+1}-y_{i}\}+3\epsilon$ implies that all other diagonal gaps have a vertical width smaller than $y_{1}-4\epsilon-\theta(G)$.

Hence the first widest gap in $\text{PD}\{C^{\alpha}\}$ covers the diagonal strip $\{\theta(G)+2\epsilon<y-x<y_{1}-2\epsilon\}$ within the first widest gap $\{\theta(G)<y-x<y_{1}\}$ in the original diagram $\text{PD}\{G^{\alpha}\}$. Then the subdiagram $\text{DS}_{1}(C^{\alpha})$ above the line $y-x=y_{1}-2\epsilon$ contains exactly $m$ perturbations of the original dots $(0,y_{i})$ in the vertical strip $\{0\leq x\leq\epsilon\}$. Hence the reduced skeleton $\mathrm{HoPeS}_{1,1}(C)$ contains exactly m critical edges corresponding to all m dots of the subdiagram $\text{DS}_{1}\{C^{\alpha}\}$.

It remains to prove that $\mathrm{HoPeS}_{1,1}(C)$ is $2\epsilon$-close to $G$. Let the critical scale $\alpha(C)$ be the maximum birth over all dots in $\text{DS}_{1}\{C^{\alpha}\}$. These dots are at most $\epsilon$ away from their corresponding points $(0,y_{i})$ in the vertical axis, so critical scale $\alpha(C)$ is at most $\epsilon$. A longest edge $e$ of any cycle in $\mathrm{HoPeS}_{1,1}(C)$ persisting over birth $\leq\alpha<$ death has the half-length equal to its birth, because adding $e$ gave birth to the homology class. Then all edges in $\mathrm{HoPeS}_{1,1}(C)$ have half-lengths at most $\alpha(C)\leq\epsilon$. Hence $\mathrm{HoPeS}_{1,1}(C)$ is covered by the disks that have the radius $\epsilon$ and centres at all points of $C$, so $\mathrm{HoPeS}_{1,1}(C)\subset C^{\epsilon}\subset G^{2\epsilon}$. ∎

Reconstruction Theorem 5.9 for other derived subgraphs $\mathrm{HoPeS}_{k,l}(C)$ will follow from Lemmas 5.5, 5.6, 5.7, 5.8 and will also imply Corollary 5.10.

Lemma 5.5**.**

Let $C$ be any $\epsilon$-sample of a connected graph $G$ in a metric space. If $|\text{dgap}_{k}(G)|-|\text{dgap}_{k+1}(G)|>8\epsilon$, then there is a bijection $\psi:\text{DS}_{k}(G)\to\text{DS}_{k}(C)$ so that $||q-\psi(q)||_{\infty}\leq\epsilon$ for all $q\in\text{DS}_{k}(G)$.

Proof.

By Stability Theorem 3.6 there is a bijection $\psi:\text{PD}\{G^{\alpha}\}\to\text{PD}\{C^{\alpha}\}$ such that $q,\psi(q)$ are $\epsilon$-close in the $L_{\infty}$-distance on $\mathbb{R}^{2}$ for all $q\in\text{PD}\{G^{\alpha}\}$. The $\epsilon$-neighbourhood of a dot $q=(x,y)$ in the $L_{\infty}$ distance is the square $[x-\epsilon,x+\epsilon]\times[y-\epsilon,y+\epsilon]$. Under the diagonal projection $\text{pr}(x,y)=y-x$, this square maps to the interval $[y-x-2\epsilon,y-x+2\epsilon]$. Hence any diagonal gap $\{a<y-x<b\}$ in $\text{PD}\{G^{\alpha}\}$ can become thinner or wider in $\text{PD}\{C^{\alpha}\}$ by at most $4\epsilon$ due to dots $q$ ‘jumping’ to $\psi(q)$ by at most $2\epsilon$ under the projection pr.

By the given inequality, the first $k$ widest gaps $\text{dgap}_{i}(G)$ in $\text{PD}\{G^{\alpha}\}$ for $i=1,\dots,k$ are wider by at least $8\epsilon$ than all other $\text{dgap}_{i}(G)$ for $i>k$. Hence all dots between any two successive gaps from the first $k$ widest can not ‘jump’ over these wide gaps and remain ‘trapped’ between corresponding diagonal gaps in the perturbed diagram $\text{PD}\{C^{\alpha}\}$. Although the order of the first $k$ widest gaps $\text{dgap}_{i}(G)$, $i=1,\dots,k$, may not be preserved under $\psi$, the lowest $\{a<y-x<b\}$ of these $k$ diagonal gaps is respected by $\psi$ as follows. The thinner strip $S=\{a+2\epsilon<y-x<b-2\epsilon\}$ has no dots from $\text{PD}\{C^{\alpha}\}$ and has the vertical width $|S|\geq|\text{dgap}_{k}(G)|-4\epsilon>|\text{dgap}_{k+1}(G)|+4\epsilon\geq|\text{dgap}_{k+1}(C)|$.

Then the diagonal strip $S=\{a+2\epsilon<y-x<b-2\epsilon\}$ is within the lowest gap of the first $k$ widest gaps $\text{dgap}_{i}(C)$, $i=1,\dots,k$. Hence all dots above $S$ remain above $S$ under the bijection $\psi$. By Definition 3.7, all these dots above $S$ in $\text{PD}\{G^{\alpha}\}$ and in $\text{PD}\{C^{\alpha}\}$ form the diagonal subdiagrams $\text{DS}_{k}(G)$ and $\text{DS}_{k}(C)$ respectively. So $\psi$ descends to a bijection $\text{DS}_{k}(G)\to\text{DS}_{k}(C)$. ∎

Lemma 5.6**.**

Let $\psi:\text{DS}_{k}(G)\to\text{DS}_{k}(C)$ be a bijection such that $||q-\psi(q)||_{\infty}\leq\epsilon$ holds for all dots $q\in\text{DS}_{k}(G)$ as in Lemma 5.5. If $|\text{vgap}_{k,l}(G)|-|\text{vgap}_{k,l+1}(G)|>4\epsilon$ for some $l\geq 1$, then $\psi$ restricts to a bijection $\text{VS}_{k,l}(G)\to\text{VS}_{k,l}(C)$.

Proof.

The $x$-coordinate of any dot $q\in\text{DS}_{k}(G)$ changes under the given bijection $\psi$ by at most $\epsilon$. Similarly to the proof of Lemma 5.5 each vertical gap $\text{vgap}_{k,l}(G)$ becomes thinner or wider by at most $2\epsilon$. By the given inequality the first $l$ gaps $\text{vgap}_{k,j}(G)$ in $\text{PD}_{k}\{G^{\alpha}\}$ for $j=1,\dots,l$ are wider by at least $4\epsilon$ than all other $\text{vgap}_{k,j}(G)$ for $j>l$. Hence all dots between any two successive gaps from the first $k$ widest can not ‘jump’ over these wide gaps and remain ‘trapped’ between corresponding vertical gaps in the perturbed diagram $\text{PD}\{C^{\alpha}\}$.

Although the order of the first $l$ widest $\text{vgap}_{k,j}(G)$, $j=1,\dots,l$, may not be preserved under $\psi$, the leftmost $\{a<x<b\}$ of these $l$ vertical gaps is respected by $\psi$ in the following sense. The thinner strip $S=\{a+\epsilon<x<b-\epsilon\}$ contains no dots from $\text{DS}_{k}(C)$ and has the horizontal width $|S|\geq|\text{vgap}_{k,l}(G)|-2\epsilon>|\text{dgap}_{k,l+1}(G)|+2\epsilon\geq|\text{dgap}_{k,l+1}(C)|$. Then the vertical strip $S=\{a+\epsilon<x<b-\epsilon\}$ is within the leftmost of the first $l$ widest $\text{vgap}_{k,j}(C)$, $j=1,\dots,l$. Hence all dots to the left of $S$ remain to the left of $S$ under the bijection $\psi$. By Definition 3.8, all these dots to the left of $S$ in $\text{PD}\{G^{\alpha}\}$ and in $\text{PD}\{C^{\alpha}\}$ form the vertical subdiagrams $\text{VS}_{k,l}(G)$ and $\text{VS}_{k,l}(C)$ respectively. So $\psi$ descends to a bijection $\text{VS}_{k,l}(G)\to\text{VS}_{k,l}(C)$ between smaller vertical subdiagrams. ∎

Lemma 5.7**.**

The derived skeleton $\mathrm{HoPeS}_{k,l}(C)$ is within $\mathrm{HoPeS}(C;\text{vs}_{k,l}(C))$.

Proof.

By Definition 3.10 all edges of the reduced skeleton $\mathrm{HoPeS}(C;\text{vs}_{k,l}(C))$ have a half-length at most $\text{vs}_{k,l}(C)$ and $\text{death}>\text{vs}_{k,l}(C)$ for all critical edges. Definition 5.1 also imposes the extra restriction on critical edges of $\mathrm{HoPeS}_{k,l}(C)$, namely each dot $(\text{birth},\text{death})$ is in the vertical subdiagram $\text{VS}_{k,l}(C)$. So $\mathrm{HoPeS}_{k,l}(C)\subset\mathrm{HoPeS}(C;\text{vs}_{k,l}(C))$. ∎

Lemma 5.8** (approximation by reduced $\mathrm{HoPeS}$).**

Let $C$ be any finite $\epsilon$-sample of a subspace $G$ in a metric space $M$. Then the reduced skeleton $\mathrm{HoPeS}(C;\alpha)$ for any scale $\alpha>0$ is contained within the $(\epsilon+\alpha)$-offset $G^{\epsilon+\alpha}\subset M$.

Proof.

Any edge $e\subset\mathrm{HoPeS}(C;\alpha)$ has a length at most $2\alpha$ by Definition 3.10, hence is covered by the balls with the radius $\alpha$ and centres at the endpoints of $e$. Then $\mathrm{HoPeS}(C;\alpha)\subset C^{\alpha}\subset G^{\epsilon+\alpha}$ since $C\subset G^{\epsilon}$ is an $\epsilon$-sample of $G$. ∎

Theorem 5.9 guarantees a reconstruction with a correct homology $H_{1}$ and within a $2\epsilon$-offset for any $\epsilon$-sample of a graph $G$ with $k$ cycles whose $\text{PD}\{G^{\alpha}\}$ has $k$ dots separable by the $k$-th widest diagonal gap from noisy artefacts. This is formalised below in the 4 conditions, which fail if the noise level $\epsilon$ is too high.

Theorem 5.9** (reconstruction by derived skeletons).**

Let $C$ be any $\epsilon$-sample of an unknown graph $G$ in a metric space. Let $G$ satisfy the conditions below.

(1) All cycles $L\subset G$ are ‘persistent’, i.e. $\text{death}(L)\geq\text{ds}_{k}(G)$ for some $k\geq 1$.

(2) $|\text{dgap}_{k}(G)|$ ‘jumps’, i.e. $|\text{dgap}_{k}(G)|-|\text{dgap}_{k+1}(G)|>8\epsilon$ for $k$ from (1).

(3) No cycles are born in offsets $G^{\alpha}$ for small $\alpha$, i.e. $\text{vs}_{k,l}(G)=0$ for some $l$.

(4) $|\text{vgap}_{k,l}(G)|$ ‘jumps’: $|\text{vgap}_{k,l}(G)|-|\text{vgap}_{k,l+1}(G)|>4\epsilon$ for $k,l$ in (1),(3).

Then we get the lower bound for noise $\epsilon\geq\text{vs}_{k,l}(C)$ and the derived skeleton $\mathrm{HoPeS}_{k,l}(C)\subset G^{2\epsilon}$ has the same homology $H_{1}$ as the underlying graph $G$.

Proof.

Due to condition (2), Lemma 5.5 implies that there is a bijection $\psi:\text{DS}_{k}(G)\to\text{DS}_{k}(C)$ so that $||q-\psi(q)||_{\infty}\leq\epsilon$ for all $q\in\text{DS}_{k}(G)$. Lemma 5.6 due to condition (4) implies that $\psi$ descends to a bijection $\text{VS}_{k,l}(G)\to\text{VS}_{k,l}(C)$. All cycles in a graph $G$ give birth to corresponding homology classes in $H_{1}(G^{\alpha})$ at the scale $\alpha=0$. These classes may split later at $\alpha>0$, but will eventually die and always give dots $(0,\text{death})\in\text{PD}\{G^{\alpha}\}$ in the vertical death axis. For any cycle $L\subset G$, let $\text{death}(L)$ be the maximum $\alpha$ such that $H_{1}(G^{\alpha})$ has the class $[L]$. The graph $\theta$ in Fig. 7 has 3 cycles with $\text{death}(L)\approx 2.577$.

Let $L_{1},\dots,L_{m}\subset G$ be all $m$ cycles generating $H_{1}(G)$. Then the 1D diagram $\text{PD}\{G^{\alpha}\}$ contains $m$ dots $(0,\text{death}(L_{i}))$, $i=1,\dots,m$, because each class $[L_{i}]$ persists over $0\leq\alpha<\text{death}(L_{i})$ by Definition 3.4. Condition (1) implies that all dots $(0,\text{death}(L_{i}))$ belong to the subdiagram $\text{DS}_{k}(G)$, hence to $\text{VS}_{k,l}(C)$.

Condition (3) $\text{vs}_{k,l}(G)=0$ means that the leftmost of the first $l$ widest $\text{vgap}_{k,j}(G)$, $j=1,\dots,l$, is attached to the vertical death axis in $\text{PD}\{G^{\alpha}\}$, which should contain the vertical subdiagram $\text{VS}_{k,l}(G)$. So $\text{VS}_{k,l}(G)$ consists of the $m$ dots $(0,\text{death}(L_{i}))$, $i=1,\dots,m$. Then the vertical subdiagram $\text{VS}_{k,l}(C)$ for the cloud $C$ also has exactly $m$ dots, which are ‘noisy’ images $\psi(0,\text{death}(L_{i}))$, $i=1,\dots,m$. Moreover, all these dots in the vertical subdiagram $\text{VS}_{k,l}(C)$ are at most $\epsilon$ away from the vertical axis, so the vertical scale $\text{vs}_{k,l}(C)$ is at most $\epsilon$.

The lowest of the points $\psi(0,\text{death}(L_{i}))$ has a death with the lower bound $\min\limits_{i=1,\dots,m}\text{death}(L_{i})-\epsilon\geq\text{DS}_{k}(G)-\epsilon>|\text{dgap}_{k}(G)|-\epsilon>7\epsilon>\text{vs}_{k,l}(C)$.

So the condition $\text{death}>\text{vs}_{k,l}(C)$ from Definition 5.1 is satisfied. Hence the reduced skeleton $\mathrm{HoPeS}_{k,l}(C)$ contains exactly $m$ critical edges corresponding to all $m$ dots in the subdiagram $\text{VS}_{k,l}(C)$. All these $m$ critical edges of $\mathrm{HoPeS}_{k,l}(C)$ generate $H_{1}$ of the required rank $m$. The geometric approximation $\mathrm{HoPeS}_{k,l}(C)\subset G^{2\epsilon}$ follows from Lemmas 5.7, 5.8 for $\alpha=\text{vs}_{k,l}(C)\leq\epsilon$. ∎

The first attempt in a graph reconstruction from a cloud $C$ by Theorem 5.9 can be the derived skeleton $\mathrm{HoPeS}_{1,1}(C)$. If the conditions fail for $k=l=1$, Theorem 5.9 might work for other pairs of the parameters $k,l$, see Fig. 46.

Corollary 5.10**.**

In the conditions of Theorem 5.9, if another cloud $\tilde{C}$ is $\delta$-close to $C$, then the perturbed skeleton $\mathrm{HoPeS}_{k,l}(\tilde{C})$ is $(2\delta+4\epsilon)$-close to $\mathrm{HoPeS}_{k,l}(C)$.

Proof.

The condition that the perturbed cloud $\tilde{C}$ is $\delta$-close to the original cloud $C$, which is $\epsilon$-close to the graph $G$, implies that $\tilde{C}$ is $(\delta+\epsilon)$-close to $G$.

Reconstruction Theorem 5.9 for the $\epsilon$-sample $C$ and $(\delta+\epsilon)$-sample $\tilde{C}$ of $G$ says that $\mathrm{HoPeS}_{k,l}(C)$ is $2\epsilon$-close to $G$ and $\mathrm{HoPeS}_{k,l}(\tilde{C})$ is $(2\delta+2\epsilon)$-close to $G$. Hence $\mathrm{HoPeS}_{k,l}(C)$ and $\mathrm{HoPeS}_{k,l}(\tilde{C})$ are $(2\delta+4\epsilon)$-close as required. ∎

6 The 79K Dataset of Random Noisy Point Samples of Planar Graphs

Mapper, $\alpha$-Reeb and $\mathrm{HoPeS}$ have guarantees on the homotopy type (or the number of independent cycles in homology $H_{1}$) of reconstructed graphs, not on a more advanced homeomorphism type that captures branches or filaments outside closed cycles. Hence, for experimental comparisons, it is enough to generate noisy samples of planar graphs without degree 1 vertices.

Subsection 6.1 introduces three patterns of planar graphs around which point clouds are randomly sampled: wheel, grid and hexagonal, see examples in Fig. 2. Subsection 6.3 discusses two noise models: uniform and Gaussian.

6.1 Patterns: families of connected planar graphs for generating noisy samples

$\bullet$ The $k$-wheel graph $W(k)\subset\mathbb{R}^{2}$ has $k\geq 3$ circumference vertices equally distributed along the unit circle centred at the central vertex at the origin. $W(k)$ has edges between the central vertex and all circumference vertices, and also edges between successive circumference vertices, see $W(4)$ in Fig. 2.

$\bullet$ For $k,l\geq 1$, the $(k,l)$-grid graph $G(k,l)$ has vertices at the points $(i,j)$ with integer coordinates $0\leq i\leq k$, $0\leq j\leq l$. Each vertex $(i,j)$ is connected to up to 4 neighbours $(i\pm 1,j)$, $(i,j\pm 1)$ in the rectangle $[0,k]\times[0,l]\subset\mathbb{R}^{2}$, see Fig. 2.

$\bullet$ The $k$-hexagons graph $H(k)$ consists of the boundaries of $k$ regular hexagons with edges of unit length, e.g. $H(1)$ is one regular hexagon. The $(k+1)$-hexagons graph is obtained from the $k$-hexagons graph by adding the boundary of a new hexagon. The last picture of Fig. 2 shows the order of added hexagons for $k\leq 7$. For $k>7$, more hexagons are similarly added in circular layers.

6.2 Generating random points in edges of a graph according to length

For a fixed graph $G\subset\mathbb{R}^{2}$ from one of the patterns above, each random point is uniformly generated along edges of $G$ as follows. Fix any order of edges of $G$. Then the continuous graph $G$ is parameterised by a single variable $t$ that takes values in the interval $[0,\sum\limits_{i=1}^{k}l_{i}]$, where $l_{1},\dots,l_{k}$ are the edge-lengths of $G$.

If the uniform random variable $t$ belongs to the $j$-th interval $[\sum\limits_{i=1}^{j-1}l_{i},\sum\limits_{i=1}^{j}l_{i}]$, a point is chosen in the $j$-th edge as the weighted combination $w\vec{u}+(1-w)\vec{v}$ of the endpoints $\vec{u},\vec{v}\in\mathbb{R}^{2}$ with $w=(t-\sum\limits_{i=1}^{j-1}l_{i})/l_{j}$. We sample 100 points per unit length of the embedded graph $G\subset\mathbb{R}^{2}$, e.g. 400 points for $G_{1,1}$.

6.3 Noise models: uniform and Gaussian perturbations of points

For each sampled point $p$ lying on an edge $e$ of a graph, we generate two independent random shifts $d_{e},d_{\perp}$ with the same distribution as follows.

$\bullet$ Uniform noise with an upper bound $\mu$: $d_{e},d_{\perp}$ are uniform in $[-\mu,\mu]$.

$\bullet$ Gaussian noise with mean 0 and a standard deviation $\sigma$: $d_{e},d_{\perp}\in\mathbb{R}$ have the Gaussian density $f(t,\sigma)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{t^{2}}{2\sigma^{2}}}$ for $t\in\mathbb{R}$.

Then the point $p\in e$ is shifted by $d_{e}$ units parallel to $e$ (in one of two fixed directions), and by $d_{\perp}$ units in the straight line perpendicular to $e$, see Fig. 11.

6.4 Types of randomly sampled clouds obtained by varying parameters

The Planar Graph Cloud (PGC) dataset contains 79K clouds, 200 clouds for each of the following 395 types of clouds counted over 3 graph patterns.

$\bullet$ 70 wheel types: $k$-wheel graphs $W_{k}$ have 7 values of $k$ from $3$ to $9$; the bound $\mu$ of uniform noise has 5 values from $0.05$ to $0.25$ in $0.05$ intervals; the deviation $\sigma$ of Gaussian noise has 5 values from $0.02$ to $0.1$ in $0.02$ intervals.

$\bullet$ 108 grid types: grids have 6 pairs $(k,l)=(1,1),(2,1),(3,1),(2,2),(3,2),(3,3)$; the bound $\mu$ of uniform noise has 8 values from $0.05$ to $0.4$ in $0.05$ intervals; the deviation $\sigma$ of Gaussian noise has 10 values from $0.02$ to $0.2$ in $0.02$ intervals.

$\bullet$ 217 hexagonal types: the $k$-triangles graph $H_{k}$ has 7 values of $k$ from 1 to 7; the bound $\mu$ of uniform noise has 15 values from $0.05$ to $0.75$ in $0.05$ intervals; the deviation $\sigma$ of Gaussian noise has 16 values from $0.02$ to $0.32$ in $0.02$ intervals.

The above intervals for noise parameters were chosen to report maximum noise when the skeletonisation algorithms produce correct reconstructions (for the 1st Betti number) with high success rates, see Fig. 23, 25 in section 8.

7 Drawing and simplifying skeletons of point clouds

By Lemma 4.2, for any cloud $C\subset\mathbb{R}^{d}$ and the filtration $\{C^{\alpha}\}$ of $\alpha$-offsets, $\mathrm{HoPeS}(C)$ is a subgraph of a Delaunay triangulation $\mathrm{Del}(C)\subset\mathbb{R}^{d}$. Hence the skeleton $\mathrm{HoPeS}(C)$ is embedded into $\mathbb{R}^{d}$ (has no intersection of edges by definition), i.e. visualises the shape of $C$ directly in the space where $C$ lives.

7.1 Drawing the abstract Mapper and $\alpha$-Reeb graphs in the plane

Both Mapper and $\alpha$-Reeb outputs are abstract graphs, and so do not naturally embed in $\mathbb{R}^{d}$. Therefore, in order to draw the graphs in the plane when $C\subset\mathbb{R}^{2}$, it is necessary for us to project them to $\mathbb{R}^{2}$ as naturally as possible.

Each vertex $v$ of the Mapper graph corresponds to a cluster $C_{v}$ in a given cloud $C\subset\mathbb{R}^{d}$ and is naturally placed at the geometric centre $\frac{1}{|C_{v}|}\sum\limits_{p\in C_{v}}p$.

For $\alpha$-Reeb graphs, it is far less natural to appropriately embed them in $\mathbb{R}^{2}$, because final stage 5 in section B.3 identifies several disjoint intervals. Every central point $v$ of such an interval $J_{v}$ is mapped to the geometric centre of the corresponding connected subgraph similarly to the Mapper graph.

In all cases the edges between vertices are drawn as straight-line segments.

Since the algorithms are compared on noisy samples around planar graphs without degree 1 vertices, the outputs of the 3 algorithms Mapper, $\alpha$-Reeb and $\mathrm{HoPeS}$ before comparison will be simplified by removing all filamentary branches in subsection 7.2. Then the $\alpha$-Reeb graph in Fig. 32 will become a topological circle that better resembles the original noisy sample of an ellipse.

7.2 Simplification of output skeletons by pruning branches

After running all algorithms, we iteratively remove all degree 1 vertices (with their single edge), to get a graph without degree 1 vertices, see Fig. 8, 12, 13.

By Definition 3.10, $\mathrm{HoPeS}(C)$ contains $\text{MST}(C)$, which spans all points of a cloud $C$. Both $\text{MST}(C)$ and $\mathrm{HoPeS}(C)$ approximate $C$ with zero geometric error, i.e. any point of $C$ is at distance 0 from $\mathrm{HoPeS}(C)$. Hence pruning all branches in $\mathrm{HoPeS}(C)$ above will make the experimental comparison fairer.

An ideal output should have much fewer vertices than a given point cloud $C$. We have decided to simplify $\mathrm{HoPeS}$ by Algorithm 7.1 below, which preserves the homology by Corollary 7.2. This extra simplification is applied only to $\mathrm{HoPeS}$ and makes a geometric approximation worse, while Mapper and $\alpha$-Reeb graphs keep their geometric errors with many more intermediate vertices.

Algorithm 7.1**.**

For a given threshold $\epsilon$, we iteratively collapse all edges of $\mathrm{HoPeS}(C)$ with lengths less than $\epsilon$, starting from shortest edges as follows.

Step (7.1a). Any edge shorter than $\epsilon$ and contained in a triangular cycle will not be collapsed by steps below in order to preserve the homology group $H_{1}$.

Step (7.1b). To collapse an edge $e$ with endpoints $v_{1},v_{2}$, we first remove $v_{1},v_{2}$ and all edges incident on them. We now add a new vertex $v$ as follows.

Step (7.1c). If $\text{deg}(v_{1})=\text{deg}(v_{2})=2$, or $\text{deg}(v_{1})\neq 2$ and $\text{deg}(v_{2})\neq 2$, the new vertex $v$ is placed at the midpoint of the straight-line edge between $v_{1},v_{2}$.

Step (7.1d). If (say) $\deg(v_{1})\neq 2$ but $\deg(v_{2})=2$, then $v$ is put at the original location of $v_{1}$, to best preserve the geometry of the skeleton $\mathrm{HoPeS}(C)$.

Step (7.1e). We add an edge from the new vertex $v$ to each of the original neighbours of the old vertices $v_{1}$ and $v_{2}$, see the final picture in Fig. 8.

Step (7.1f). If any edges become intersecting, the collapse is reversed.

The output of Algorithm 7.1 can be denoted by $\mathrm{simHoPeS}_{\epsilon}(C)$. We always use $\epsilon$ equal to the maximum death of any dot above the first widest diagonal gap in the persistence diagram $\text{PD}\{C^{\alpha}\}$ and denote the result by $\mathrm{simHoPeS}(C)$ for simplicity. This choice of $\epsilon$ is justified by the fact that maximum death is the scale $\alpha$ when all (most persistent) cycles of $\mathrm{HoPeS}(C)$ become contractible in the offset $C^{\alpha}$, so there is no sense to collapse edges longer than $\alpha$.

Corollary 7.2**.**

Under the conditions of Theorems 5.4 and 5.9, $\mathrm{simHoPeS}(C)$ from Algorithm 7.1 has the same homology $H_{1}$ as the underlying graph $G$.

Proof.

Removing degree 1 vertices does not change the homology group $H_{1}$ of a graph. Collapsing a short edge in Algorithm 7.1 would only change $H_{1}$ if this edge was in a triangular cycle or, if considering the output as an embedded graph in $\mathbb{R}^{2}$, we get intersections of edges. Algorithm 7.1 prevents both operations by Steps (7.1a) and (7.1f). Hence $H_{1}$ remains unchanged under all simplifications and is isomorphic to $H_{1}(G)$ in both cases specified by Theorems 5.4, 5.9. ∎

Fig. 14, 15, 16 show typical outputs. Since Mapper and $\alpha$-Reeb graphs are abstract, their projections to $\mathbb{R}^{2}$ can have self-intersections. The full skeleton $\mathrm{HoPeS}(C)$ and simplified skeleton $\mathrm{simHoPeS}(C)$ have no self-intersections.

8 Comparison of the algorithms on synthetic and real data

8.1 Running over wide ranges of parameters of Mapper and $\alpha$-Reeb graphs

Unlike $\mathrm{HoPeS}$, the Mapper and $\alpha$-Reeb algorithms require additional parameters that essentially affect outputs. To get adequate outputs of Mapper and $\alpha$-Reeb, one can manually tune these parameters for every input cloud. However, for different types of synthetic and real clouds, choosing the same parameters would be inadequate. If the parameters such as $\alpha$ of the $\alpha$-Reeb algorithm and $\epsilon$, used for DBSCAN in the Mapper algorithm, are too small, then noisy cycles will also be captured in addition to the dominant cycles. If these parameters are too large, then even the dominant cycles will not be captured. $\mathrm{HoPeS}$ is always run for the two heuristics: the 1st widest diagonal gap to the derived skeleton $\mathrm{HoPeS}_{1,1}(C)$, which is simplified with the threshold $\epsilon$ equal to the max death by Algorithm 7.1. The parameters of the Mapper and $\alpha$-Reeb in the experiments below were searched over wide ranges to separately optimise each measure. Hence only best outputs of Mapper and $\alpha$-Reeb were counted.

For Mapper, the amount of overlap of the intervals in the covering of the range of the filter function was fixed at 50%. For a cloud of $n$ points, the number of intervals was chosen as $\frac{tn}{100}$, rounded to the nearest integer, where $t$ took $10$ values between $1.5$ and $3.3$ in $0.2$ increments. The clustering parameter $\epsilon$ also took $10$ values between $0.05$ and $0.5$ in $0.05$ increments. So in total we have used $100$ configurations of the parameters of Mapper. For the $\alpha$-Reeb algorithm, the scale $\alpha$ took $10$ values between $0.15$ and $0.6$ in $0.05$ increments.

8.2 Four measures for the experimental comparison of the algorithms

For each fixed graph, the measures below are averaged over 200 noisy samples.

The Betti success rate (in percentages) counts the cases when an output contains the expected number of independent cycles, which is called the first Betti number and equal to the rank of the homology group $H_{1}$. The $k$-wheel graph $W(k)$ and $k$-hexagons graph $H(k)$ have the first Betti number $k$.

Homeomorphism success rate (in percentages) is the empirical conditional probability that the output skeleton is homeomorphic to an underlying graph $G$ from which a cloud was sampled assuming the first Betti number was correct.

The root mean square (RMS) distance measures how close geometrically an output is to the cloud $C$. For each point $p\in C$, the distance $d(p,S)$ is computed from $p$ to the closest straight-line edge of the output $S$. Then the RMS error is $\sqrt{\sum_{p\in C}d(p,S)^{2}}$. For each type of cloud, we average the RMS distance over all outputs that are successful on the first Betti number.

Running time (ms) is the time averaged over all clouds and parameters.

Each point in Fig. 17-22 has the $y$-coordinate equal to the average value over 200 clouds (with optimal parameters of Mapper and $\alpha$-Reeb for each cloud) sampled around a graph $G$, the $x$-coordinate is the first Betti number of $G$.

8.3 Comparison of geometric errors on clouds of edge pixels in BSD images

We ran the algorithms over the clouds $C$ of Canny edge pixels extracted from 500 real images in the Berkeley Segmentation Database (BSD500). The last pictures in Fig. 26, 27, 28 show the derived skeleton $\mathrm{HoPeS}_{1,1}$ with critical edges in red. Table 1 shows the RMS distance from the output graphs to $C$ (in pixels, averaged over images and parameters of Mapper and $\alpha$-Reeb).

[TABLE]

The key advantage of the Mapper algorithm is its flexibility due to many parameters, e.g. any clustering algorithm can be combined with Mapper. When a suitable scale parameter $\alpha$ for a given $n$-point cloud $C$ can be quickly guessed, the $\alpha$-Reeb graph is computed with the very fast time $O(n\log n)$. The experiments on 79K clouds and 500 BSD images in section 8 have optimised the Mapper and $\alpha$-Reeb graphs over wide ranges of their important parameters.

Practically, by always using the justified heuristics (the 1st widest diagonal gap for $\mathrm{HoPeS}_{1,1}(C)$ and $\epsilon$ equal to the max death for simplification), the skeleton $\mathrm{simHoPeS}(C)$ has comparable or better results without optimising any parameters in Fig. 23, 25 and Table 1. All algorithms can be improved by a further minimisation of the RMS distance from a skeleton to a cloud $C\subset\mathbb{R}^{d}$.

Theoretically, a Homologically Persistent Skeleton $\mathrm{HoPeS}(C)$ has the advantage of being a parameter-free graph that is also embedded into the space of a cloud $C$. Hence the time-consuming tuning of parameters is avoided and no intersections of edges are guaranteed. At the same time $\mathrm{HoPeS}$ remains flexible in the sense that after the full $\mathrm{HoPeS}(C)$ is found, one can quickly extract reduced subgraphs (for a fixed scale $\alpha$) and derived subgraphs $\mathrm{HoPeS}_{k,l}(C)$.

The paper gives detailed proofs of Optimality Theorem 4.8 and Reconstruction Theorems 5.4, 5.9 for the first time. The last two results are yet to be extended to higher dimensions. Corollaries 5.10 and 7.2 are new results.

The C++ code for all three algorithms is available at https://github.com/Phil-Smith1/Cloud_Skeletonization_3.git. The dataset of 79K point clouds (2GB) for comparison with other skeletonisation algorithms is available by request.

A Basic definitions and results from geometry and topology

Definition A.1** (metric space).**

A metric space $M$ is any set of points with a distance function $d(p,q)$ between any points $p,q$ that satisfies these axioms:

$\bullet$ Positiveness: $d(p,q)\geq 0$ with equality if and only if $p=q$.

$\bullet$ Symmetry: $d(p,q)=d(q,p)$ for any points $p,q\in C$.

$\bullet$ The triangle inequality: $d(p,q)+d(q,r)\geq d(p,r)$ for any $p,q,r\in C$.

For any subset $C$ of a metric space $M$ and a parameter $\alpha\geq 0$, the $\alpha$-offset $C^{\alpha}=\{p\in M\mid d(p,q)\leq\alpha\mbox{ for some }q\in M\}$ is the union of closed balls with the radius $\alpha$ and centres at all points of $C$, e.g. $C^{0}=C$. If a metric space $M$ consists of only finitely many points, this space will be called a point cloud. $\blacksquare$

The Euclidean distance for points $p,q\in\mathbb{R}$ is $d(p,q)=|p-q|$. The function $(p-q)^{2}$ fails the triangle inequality above for the points $p=-1,q=0,r=1$.

Definition A.2** (metric graph and neighbourhood graph).**

A graph is a finite set of vertices with edges being unordered pairs of vertices. A metric graph is a graph that has a length assigned to each edge. The distance between two vertices is the minimum total length of any path (if it exists) from one vertex to the other. An example of a metric graph is the neighbourhood graph $N(C,\epsilon)$ of a point cloud $C$ with a threshold $\epsilon$. The vertex set of the graph $N(C,\epsilon)$ is the cloud $C$. Two vertices $p,q\in C$ are joined by an edge in $N(C,\epsilon)$ if the distance $d(p,q)\leq\epsilon$. The length of this edge is fixed as the distance $d(p,q)$. $\blacksquare$

Definition A.3** (graph homeomorphism).**

A graph $G$ is connected if any two vertices of $G$ are connected by a sequence of edges such that any successive edges share a common vertex. Connected graphs $G,H$ are homeomorphic, $G\simeq H$, if, ignoring all vertices of degree 2, there exists a bijection $\psi$ between the vertices of $G$ and $H$ that respects edges, i.e. any vertices $u,v\in G$ are connected by $m$ edges in $G$ if and only if $\psi(u),\psi(v)\in H$ are connected by $m$ edges in $H$. $\blacksquare$

The graphs $W(4)$ and $G(2,2)$ in Fig. 2 are homeomorphic, because four corner vertices of degree 2 in $G(2,2)$ are topologically trivial so that $G(2,2)$ can be mapped to $W(4)$ via a continuous bijection whose inverse is continuous.

The vertex set $C$ of any connected metric graph $G$ is a point cloud in the sense of Definition A.1, where the distance $d(p,q)$ between any vertices $p,q\in C$ is the total length (sum of edge-lengths) of a shortest path from $p$ to $q$ in $G$.

Definition A.4** (simplicial complex).**

A simplicial complex $Q$ with a finite vertex set $C$ is a collection of subsets $\{v_{0},\dots,v_{k}\}\subset C$ called simplices so that

$\bullet$ any subset (called a face) of a simplex is also a simplex and is included in $Q$;

$\bullet$ any non-empty intersection of simplices is their common face included in $Q$.

Any simplex on $k+1$ vertices has the dimension $k$ and geometric realisation:

[TABLE]

The geometric realisation of $Q$ is obtained by gluing realisations of all its simplices along their common faces. Hence $Q$ inherits the Euclidean topology. The dimension of a complex $Q$ is the maximum dimension of its simplices. $\blacksquare$

A 1-dimensional complex $Q$ is a graph $G$ without loops (edges connecting a vertex to itself) and multiple edges that have the same endpoints. A geometric realisation of $G$ is a continuous function $f:G\to\mathbb{R}^{k}$ mapping edges to straight line segments in $\mathbb{R}^{k}$ that can meet only at common endpoints, see Fig. 2.

Definition A.5** (C̆ech and Vietoris-Rips complexes, Delaunay triangulations and $\alpha$-complexes).**

Let $C$ be any set in a metric space $M$. For abstract data, the space $M$ can coincide with a given cloud $C$. The *C̆ech * complex $\text{\emph{\u{C}h}}(C,M;\alpha)$ has a simplex on $v_{0},\dots,v_{k}\in C$ if the full intersection of $k+1$ closed balls with a radius $\alpha$ and centres $v_{0},\dots,v_{k}$ is not empty, i.e. contains a point from the ambient space $M$, see [20, section 4.2.3] and Fig. 29. The Vietoris-Rips complex $VR(C,M;\alpha)$ has a simplex on vertices $v_{0},\dots,v_{k}\in C$ if each pairwise intersection of the above $k+1$ balls is not empty in $M$. For a cloud $C\subset\mathbb{R}^{d}$, a Delaunay triangulation $\mathrm{Del}(C)$ consists of simplices (with vertices at points of $C$) whose minimal open circumballs contain no points of $C$. For $\alpha\geq 0$, the $\alpha$-complex $C(\alpha)\subset\mathrm{Del}(C)$ consists of all Delaunay simplices whose circumradii are at most $\alpha$, i.e. $C(\alpha)$ grows from $C$ at $\alpha=0$ to $\mathrm{Del}(C)$ when $\alpha\to+\infty$. $\blacksquare$

Lemma A.6 says that the C̆ech complex correctly represents the geometric shape of any union of balls. The Vietoris-Rips complex is determined by its 1-dimensional skeleton, hence requires much less storage than the C̆ech complex. The $\alpha$-complex $C(\alpha)$ is preferred for low $d$, because any $C(\alpha)$ lives in $\mathbb{R}^{d}$, while C̆ech and VR complexes can be only projected to $\mathbb{R}^{d}$ (often with intersections).

Lemma A.6**.**

[21, Theorem 3.2] Let $C$ be a subspace of a metric space $M$. Then any $\alpha$-offset $C^{\alpha}\subset M$ is homotopy equivalent to the C̆ech complex $\text{\emph{\u{C}h}}(C,M;\alpha)$. For a finite cloud of points $C\subset\mathbb{R}^{d}$ in a Euclidean space, any $\alpha$-offset $C^{\alpha}\subset\mathbb{R}^{d}$ is homotopy equivalent to the $\alpha$-complex $C(\alpha)\subset\mathbb{R}^{d}$.

Definition A.7** (Reeb graphs).**

For a simplicial complex $Q$ and a function $f:Q\rightarrow\mathbb{R}$, the level set of $f$ corresponding to a value $t\in\mathbb{R}$ is the set of points $L_{t}(f)=\{x\in Q\mid f(x)=t\}$. We define, for points $x,y\in Q$, the equivalence relation $\sim$ such that $x\sim y$ if and only if $f(x)=f(y)$ and $x$ and $y$ are in the same connected component of the level set of $f(x)$, see Fig. 30. The Reeb graph $\text{Reeb}(Q,f)$ is the quotient space of $Q$ formed by mapping all points that are equivalent to each other under the relation $\sim$ to a single point.

Definition A.7 makes sense when $Q$ is a topological space, which is more general than a simplicial complex. A Homologically Persistent Skeleton ($\mathrm{HoPeS}$) is essentially based on the concept of homology introduced below. A 1-dimensional cycle in a complex $Q$ is a sequence of edges $e_{1},\dots,e_{k}$ such that $e_{i}$ and $e_{i+1}$ have a common endpoint, where $e_{k+1}=e_{1}$. We define the first homology group $H_{1}(Q)$ of a simplicial complex $Q$ only with coefficients in the group $\mathbb{Z}_{2}=\mathbb{Z}/2\mathbb{Z}=\{0,1\}$.

Definition A.8** (homology $H_{1}(Q)$ of a complex).**

Cycles of a complex $Q$ can be algebraically written as finite linear combinations of edges (with coefficients only [math] or $1$) and generate the vector space $C_{1}$ of cycles. The boundaries of all triangles in $Q$ are cycles of 3 edges and generate the subspace $B_{1}\subset C_{1}$. The quotient space $C_{1}/B_{1}$ is the homology group $H_{1}(Q)$. The operation is the addition of cycles, the empty cycle is the zero. We define the $k$-th Betti number of a complex $Q$ to be the rank of the $k$-th homology group. For instance, the first Betti number of a complex $Q$ is the rank of the first homology group, and is equivalent to the number of linearly independent 1-dimensional cycles in $Q$.

The complex $Q=\text{\emph{\u{C}h}}(C,\mathbb{R}^{2};\alpha)$ in Fig. 29 has 5 edges $e_{i}$, $i=1,\dots,5$. The vector space $C_{1}$ is generated by two triangular cycles $e_{1}+e_{2}+e_{5}$ and $e_{3}+e_{4}+e_{5}$. The subspace $B_{1}$ is generated by the boundary cycle $e_{3}+e_{4}+e_{5}$ of the yellow triangle. Hence $H_{1}(Q)=C_{1}/B_{1}=\mathbb{Z}_{2}$ is generated by the cycle $e_{1}+e_{2}+e_{5}$.

B A detailed review of DBSCAN, Mapper and $\alpha$-Reeb algorithms

B.1 DBSCAN: the clustering algorithm used in our Mapper implementation

DBSCAN is a density-based spatial clustering for applications with noise. Mapper in subsection B.2 should work for any number of clusters in subclouds, hence we can not use clustering that needs a specified number of clusters, e.g. $k$-means. We have also tried the single-edge clustering and found that DBSCAN works better for noisy samples from section 6. DBSCAN requires two parameters: (1) $\epsilon$ is the radius around a point within which we search for neighbours; and (2) minPts is the number of points required within a neighbourhood of a point before a cluster can be formed. Here are the stages of DBSCAN.

$\bullet$ Stage 1: A single point $p_{1}$ is selected at random, and we compute the set $\text{Nbhd}(p_{1})$ of all points within a distance $\epsilon$ of $p_{1}$. If $|\text{Nbhd}(p_{1})|<\text{minPts}$, then $p_{1}$ is labelled as noise. If however $|\text{Nbhd}(p_{1})|\geq\text{minPts}$, we label all points in $\text{Nbhd}(p_{1})$ (not already belonging to another cluster, but may have previously been labelled as noise) as belonging to the cluster of $p_{1}$.

$\bullet$ Stage 2: Then, looping over all points $p_{i}$ in $\text{Nbhd}(p_{1})$ (apart from $p_{1}$), we compute $\text{Nbhd}(p_{i})$. All points in $\text{Nbhd}(p_{i})$ (not already belonging to another cluster) are labelled as belonging to the cluster of $p_{1}$, and, if $|\text{Nbhd}(p_{i})|\geq\text{minPts}$, we add all points in $\text{Nbhd}(p_{i})$ to $\text{Nbhd}(p_{1})$. Hence, this loop will finish only when all points in $p_{1}$’s cluster have been identified.

$\bullet$ Stage 3: A new point $p_{2}$, that is not already labelled, is selected, and we continue as in the first two stages for $p_{2}$ instead of $p_{1}$. The process continues until all points are assigned to a cluster or are labelled as noise.

B.2 The Mapper graph: a network of interlinked clusters

The Mapper algorithm is a skeletonisation framework introduced by G. Carlsson et al. [13], which aims to give a simple description of a dataset by a network of interlinked clusters. Mapper takes as input a point cloud $C$ given by pairwise distances with extra parameters below, and outputs a simplicial complex.

A filter function is a function $f:C\to Y$ whose value should be given for all points in a data cloud $C$. The parameter space $Y$ will often be $\mathbb{R}$, but could also be $\mathbb{R}^{2}$ or $S^{1}$. The type of a filter function is up to the user, with common examples being a density estimator, the Euclidean distance from a base point, or the distance from a root point within a neighbourhood graph on $C$.

A covering of the range of $f$. The range of the filter function $f$ must be covered by overlapping regions, e.g. line intervals for $Y=\mathbb{R}$. A cover usually has two parameters – the number of regions and the ratio of overlap – which can be used to control the resolution of the output simplicial complex.

A clustering algorithm. The Mapper algorithm requires subsets of the input point cloud $C$ to be clustered. So a clustering algorithm must be chosen, with a good choice again being dependent on the type of data being analysed.

The main stages of the Mapper algorithm are now described below.

Stage 1. If values of a given filter function $f$ are not yet given explicitly, then the function $f$ is computed on all points of a point cloud $C$, see Fig. 31.

Stage 2. For each region $I$ in the covering of the range of $f$, we cluster the preimage $f^{-1}(I)=\{p\in C\mid f(p)\in I\}$ according to the given clustering algorithm. Each cluster is represented as a vertex in the output complex $Q$.

Stage 3. If any $k\geq 2$ resulting clusters (across all regions in the covering) share a common point of $C$, the corresponding $k$ vertices (representatives of $k$ clusters) span a $(k-1)$-dimensional simplex in the Mapper complex $Q$.

A simple way to visualise the Mapper complex $Q$ is to project $Q$ to the space where the cloud $C$ lives by mapping every vertex to the average point of the corresponding cluster. If edges are drawn as straight-line segments, e.g. in $\mathbb{R}^{2}$, then intersections may be unavoidable, see the 2nd picture in Fig. 1.

Mapper is a versatile tool that if used rightly can be a useful method in visualising large datasets. A significant difficulaty of the method however is that with so many choices to be made, the user often requires existing knowledge of the dataset in order to select parameters that will give meaningful outputs.

The experiments in Section 8 use a filter function mapping a point cloud $C$ to $\mathbb{R}$, and therefore Mapper outputs a one-dimensional complex, which is a graph. The filter function is the Euclidean distance from a base point, which is the most distant from a random point of $C$. The overlap percentage of two adjacent intervals in the covering is $50\%$. The clustering algorithm is DBSCAN, which has two parameters: 1) a radius $\epsilon$ around a point within which we search for neighbours, and 2) the minimum number of points in a cluster. We found it acceptable to fix this minimum at 5 and optimised $\epsilon$ in the experiments.

B.3 The $\alpha$-Reeb graph is a parametric version of the Reeb graph

The Reeb graph, introduced in Definition A.7, is a simplified representation of a simplicial complex. If, instead of having an entire simplicial complex $Q$, we only have a finite number of points sampled from $Q$, even approximating the Reeb graph of $(Q,f)$ is not straightforward. Therefore, to bridge this barrier, F. Chazal et al. introduced a parametric version called the $\alpha$-Reeb graph [4].

Definition B.1**.**

Let $(Q,d)$ be a simplicial complex $Q$ with a continuous function $d:Q\to\mathbb{R}$. Let $\alpha>0$ and let $\mathcal{I}=\{I_{i}\}$ be a covering of the range of $d$, where each $I_{i}$ is a closed interval of the length $\alpha$. Consider the transitive closure of the following relation on points of $Q$: $x\sim_{\alpha}y$ if and only if $d(x)=d(y)$ and $x$ and $y$ are in the same connected component of $d^{-1}(I_{i})=\{x\in Q\>|\>d(x)\in I_{i}\}$ for some interval $I_{i}\in\mathcal{I}$. Then the $\alpha$-Reeb graph associated to the covering $\mathcal{I}$ of a simplicial complex $Q$, is the quotient space formed from $Q$ by mapping all points that are equivalent to each other under $\sim_{\alpha}$ to a single point.

Here are the stages of the $\alpha$-Reeb algorithm [4, section 6].

Stage 0. To get a simplicial complex $Q$ needed for the $\alpha$-Reeb graph in Definition B.1, usually one creates a neighbourhood graph $N(C,\epsilon)$ where all points of $C$ at a distance less than $\epsilon$ are connected by an edge, see Fig. 32.

Stage 1. Select a root vertex in $N(C,\epsilon)$, e.g. as the most distance vertex from a random one. For each vertex, calculate its distance from the root within the graph, which gives the function $d:V\to\mathbb{R}$ on the vertex set $V$ of $N(C,\epsilon)$.

Stage 2. Let $\max(d)$ be the maximal value of $d:V\to\mathbb{R}$. For any $0<\alpha<\max(d)$, the range of $d$ is covered by the overlapped intervals $\mathcal{I}=\{I_{i}\}_{0\leq i\leq m}$, where $I_{i}=\{[\frac{i\alpha}{2},\frac{i\alpha}{2}+\alpha]\}$, $m$ is the smallest integer such that $m\geq\frac{2(\max(d)-\alpha)}{\alpha}$.

Stage 3. For each interval in the covering $\mathcal{I}$, we consider its preimage $d^{-1}(I_{i})$ in the vertex set $V$. We add an edge between two vertices in the preimage if an edge also existed between these vertices in the neighbourhood graph $N(C,\epsilon)$. So each interval $I_{i}$ gives rise to a subgraph of $N(C,\epsilon)$. Every connected component of these subgraphs corresponds to a vertex in an intermediate graph $G$.

Stage 4. Two vertices in the intermediate graph $G$ are connected by an edge if their corresponding components share a common vertex in $N(C,\epsilon)$. Such common vertices are joined by dotted arcs in the last picture of Fig. 32.

Stage 5. For each vertex $v$ in our intermediate graph $G$, we take a copy of the interval $J_{v}$, which are split by $v$ into two halves and partially ordered according to $\mathcal{I}=\{I_{i}\}_{0\leq i\leq m}$. The $\alpha$-Reeb graph is the quotient of the disjoint union of these intervals $J_{v}$. The top half of $J_{v}$ is identified with the bottom half of all higher intervals $J_{u}$ such that the vertices $u,v$ are joined by an edge of $G$.

Similarly to Mapper, the $\alpha$-Reeb graph aims to join close clusters, however clustering in preimages of intervals is replaced by finding connected subgraphs. Informally, the scale $\alpha$ helps to identify points that can be connected by paths of lengths at most $\alpha$ within a neighbourhood graph. In the limit $\alpha\to 0$, the $\alpha$-Reeb graph becomes the Reeb graph from Definition A.7. Hence the $\alpha$-Reeb graph essentially depends on $\alpha$ whose wide range is explored in subsection 8.1.

Appendix C: more experimental results on synthetic and real data

C.1. Comparisons of algorithms on noisy samples of graphs

C.2. Comparisons of algorithms on edge pixels in BSD images

Final Fig. 46 highlights the applicability of $\mathrm{HoPeS}(C)$, which is computed for a cloud $C$ without any parameters. Then a user may explore derived skeletons $\mathrm{HoPeS}_{k,l}(C)$ only by looking at gaps in the persistence diagram $\text{PD}\{C^{\alpha}\}$ and choosing a certain number of persistent cycles without extra computations.