# A computationally efficient framework for vector representation of
persistence diagrams

Kit C Chan
Department of Mathematics and Statistics
Bowling Green State University
Bowling Green, OH 43403
<EMAIL_ADDRESS>
&Umar Islambekov
Department of Mathematics and Statistics
Bowling Green State University
Bowling Green, OH 43403
<EMAIL_ADDRESS>
&Alexey Luchinsky
Department of Mathematics and Statistics
Bowling Green State University
Bowling Green, OH 43403
<EMAIL_ADDRESS>
&Rebecca Sanders
Department of Mathematical and Statistical Sciences
Marquette University
Milwaukee, WI 53233
<EMAIL_ADDRESS>

###### Abstract

In Topological Data Analysis, a common way of quantifying the shape of data is
to use a _persistence diagram_ (PD). PDs are multisets of points in
$\mathbb{R}^{2}$ computed using tools of algebraic topology. However, this
multi-set structure limits the utility of PDs in applications. Therefore, in
recent years efforts have been directed towards extracting informative and
efficient summaries from PDs to broaden the scope of their use for machine
learning tasks. We propose a computationally efficient framework to convert a
PD into a vector in $\mathbb{R}^{n}$, called a _vectorized persistence block_
(VPB). We show that our representation possesses many of the desired
properties of vector-based summaries such as stability with respect to input
noise, low computational cost and flexibility. Through simulation studies, we
demonstrate the effectiveness of VPBs in terms of performance and
computational cost within various learning tasks, namely clustering,
classification and change point detection.

_K_ eywords Topological data analysis $\cdot$ Persistence diagram $\cdot$
Persistence block

###### Contents

  1. 1 Introduction
  2. 2 Related Work
  3. 3 Theory
    1. 3.1 Persistence Block
    2. 3.2 Identification Map
    3. 3.3 Vectorization
  4. 4 Simulations
    1. 4.1 Clustering using K-medoids on a synthetic dataset
    2. 4.2 Clustering using K-medoids on a real dataset
    3. 4.3 Classification of parameter values of a discrete dynamical system
  5. 5 Conclusion and future work

##  1 Introduction

In recent years, topological data analysis (TDA) has become a popular approach
to study data that have inherent shape structure [1, 2]. One of the primary
tools of TDA is _persistent homology_ (PH) [3, 4]. The idea behind PH is to
construct a sequence of nested topological spaces, called _simplicial
complexes_ (which are indexed by some scale parameter), on top of data points
and keep a record of various topological features (such as connected
components, one dimensional holes, two dimensional holes, etc.) that emerge
and later disappear as the scale parameter increases. A simplicial complex
acts as a bridge between discrete data points and the underlying object they
are sampled from, and its combinatorial nature makes it suitable for
computational purposes. The topological features are detected using the
concept of homology from algebraic topology (see e.g. [5, 6] for a detailed
treatment of the topic). The output of applying PH is a _persistence diagram_
(PD) which is a multiset of points in $\mathbb{R}^{2}$ where each point
$(b,d)$ represents a topological feature that is "born" at scale value $b$ and
"dies" at scale value $d$. A PD can be viewed as a multiscale topological
signature of the data from which it is computed.

In applications, PDs are often used within statistical and machine learning
frameworks for further analysis of the data. A direct way in this regard is to
endow the space of PDs with the bottleneck or Wasserstein metrics [7, 8] and
use them in distance-based learning settings such as k-medoids or k-nearest
neighbors (kNN). Though stability results exist for PDs with respect to these
metrics [8, 9, 10], this approach has two limitations: the scope of
applicability and the computational cost. Firstly, it is not suitable for
machine learning methods requiring input features to be elements of a Hilbert
space (the space of PDs can not be turned into a Hilbert space [11, 7]).
Secondly, as the number of points in PDs grows, so does the cost of computing
the distance between them. Therefore, in order to facilitate the use of PDs
for a broader class of machine learning methods, much research in TDA centers
around the problem of extracting useful and computationally efficient
summaries from PDs that are identified as vectors of a Hilbert space. There
are two common ways to achieve this. The first method to map PDs into a
Hilbert space is by defining a kernel function which encodes similarity
between PDs and using it in kernel-based machine learning settings [12, 13,
14, 15]. The second approach is to extract a vector from a PD which is an
element of $\mathbb{R}^{n}$ [16, 17, 18, 19, 20, 21, 22]. According to [18],
the desired properties of vectorization-based summaries are stability with
respect to input noise, computational efficiency, interpretable connection to
PDs and flexibility that allows different treatments of points of PDs based on
their relative importance.

In the present paper, we introduce a computationally efficient framework to
extract vectorizaton-based summaries that have all of the above
characteristics. First, we switch from the birth-death to birth-persistence
coordinates using the map $(b,d)\mapsto(b,p)$, where $p=d-b$. Then, for a
given persistence diagram $D$ in the new coordinates, we construct a
_persistence block_ (PB) which is a function
$\tilde{f}(D):\mathbb{R}^{2}\to\mathbb{R}$ defined as a linear combination of
characteristic functions of squares centered at the points of $D$. The width
of a square centered at point $(b,p)$ depends on the persistence value $p$ and
the scale-free parameter $\tau\in[0,1]$. The value of $\tau$ can be set in
advance or optimally tuned within the machine learning method in question (see
Section 4 for more details). Finally, we superimpose a rectangular grid on the
domain of the PB $\tilde{f}(D)$ and compute the integral of$\tilde{f}(D)$ with
respect to some chosen weight function $w$ over each grid cell. By arranging
those integral values into a vector in $\mathbb{R}^{n}$, where $n$ in the
number of grid cells, we obtain a _vectorized persistence block_ (VPB). The
idea behind a VPB is motivated by the work on _persistence images_ (PI) in
[18]. In Section 2, we discuss in more detail the similarities and differences
between the two approaches.

The rest of the paper is organized as follows. In Section 3, we formally
introduce PBs and VPBs and their continuity and stability results with respect
to the Wasserstein distance. The results of the simulation studies are
presented in Section 4. Finally, in Section 5, we summarize the contributions
of our work and highlight potential directions for future research. Our R code
to compute VPBs can be found at https://github.com/uislambekov/TDA-VPB.git.

##  2 Related Work

There are several vectorization-based summaries found in the TDA literature
such as persistence landscapes [16], persistence images [18], persistence
curves [22], persistence indicator functions [17] and persistent entropy [20],
to name a few. Of these different summaries, the work on persistence images
(PI) serves as a motivation for our method. In this section we primarily focus
on presenting a detailed discussion about the similarities and differences of
the two approaches.

A PI is obtained by integrating a persistence surface over each cell of a
superimposed grid. A persistence surface is a weighted sum of bivariate
Gaussian distributions centered at the points of the PD with the covariance
matrix $\sigma^{2}I_{2\times 2}$. The variance $\sigma^{2}$ of the Gaussians
is the same for all points in the PD and is selected by a user. The weights,
one for each point, encode a relative importance of the points in the PD. If
the birth values are all zero in a PD, as it often happens for the homological
dimension $H_{0}$, then one uses univariate Gaussian distributions instead. In
the simulation studies of [18], the weights are computed from a piecewise
linear function depending only on the persistence value and all the grid cells
are of equal size.

In contrast, for the purposes of our simulation studies, we construct a PB as
an unweighted sum of characteristic functions of the squares centered at the
points of the PD (see Section 3.2 for more details). The sizes of the squares
are not the same but rather change according to their center locations and the
parameter $\tau$. The $\tau$ parameter partially corresponds to $\sigma^{2}$,
which controls the spread of the Gaussians in a persistence surface. However,
there is an important difference: while $\sigma^{2}$ is a scale-dependent
parameter, $\tau$ is a scale-free parameter with value in $[0,1]$ which makes
it easier to select and tune in applications. Our weight function
$w:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is introduced in the vectorization
step when we compute the integral of a PB with respect to the weighted area
measure $wdA$ over each grid cell. Due to simplicity of construction, one can
derive a closed-form formula to compute these integrals for a large class of
weight functions which facilitates a fast computation of PVBs in practice.
Finally, unlike the grid cells in [18], in our simulations they are of
different sizes and selected in a data-driven way based on sample percentiles
of the birth and persistence values.

##  3 Theory

In this section, we provide a mathematical background for our analysis of
persistence diagrams. We set forth a foundation and establish continuity and
stability results for representing a persistence diagram as a function in an
infinite dimensional Hilbert space $L^{2}$ and as a vector in
$\mathbb{R}^{n}$.

###  3.1 Persistence Block

An important aspect of topological data analysis is its ability to track
various topological features appearing and disappearing within a filtration
built upon data. One approach to code this information is to record the birth
$b$ which is the time a topological feature appears, and the death $d$ which
is time the exact same feature disappears. In this way, every topological
feature corresponds to exactly one point $(b,d)$ in $\mathbb{R}^{2}$. Since a
topological feature disappears after it appears, we have $0\leq b\leq d$, and
consequently the point $(b,d)$ lies on or above the line $y=x$. An alternative
approach is to record the birth $b$ of a topological feature and its
persistence $p$, or in other words the life span (i.e. $p=d-b$), of the same
feature. Again, each topological feature corresponds to exactly one point
$(b,p)$ in $\mathbb{R}^{2}$. Since the persistence $p$ can be any nonnegative
number and if necessary by shifting the birth coordinate values by an
appropriate constant, we can assume that the point $(b,p)$ lies in the first
quadrant of $\mathbb{R}^{2}$. For the framework of the present paper, we take
the second approach and assume that a persistence diagram (PD) is a collection
of points in the first quadrant of $\mathbb{R}^{2}$, with the interpretation
that a point $(b,p)$ in a PD represents the birth and persistence of some
topological feature.

Keeping in mind that in practice we only deal with a finite number of PDs and
topological features, we can assume all persistence diagrams are finite and
there are positive real numbers $\alpha,\beta$ such that all persistence
diagrams are subsets of the compact rectangle
$\Omega=[0,\alpha]\times[0,\beta]$. The set of all finite persistence diagrams
is denoted by $P(\Omega)$ and the set of all persistence diagrams with at most
$N$ points is denoted by $P_{N}(\Omega)$. Technically speaking we can assume
all persistence diagrams are in $P_{N}(\Omega)$ for some integer $N$. However,
in a careful treatment of our theory dealing with persistence diagrams with
different numbers of points, we consider $P_{N}(\Omega)$ as a subset of
$P(\Omega)$ and $P(\Omega)=\bigcup_{N}P_{N}(\Omega)$.

To identify a point in a persistence diagram $D$ in $P(\Omega)$ with a
function in some Hilbert space $L^{2}$, we define a length function
$\lambda:\Omega\to[0,\infty)$, which is a continuous function such that
$\lambda(x,0)=0$ whenever $(x,0)\in\Omega$, and

$\displaystyle 0<\lambda(x,y)<2y\ \text{ whenever $(x,y)\in\Omega$ with
$y>0$.}$ (3.1.1)

For any point $(x,y)\in\Omega$, let $E(x,y)$ be the square centered at $(x,y)$
with side length $\lambda(x,y)$. That is,

$\displaystyle
E(x,y)=\left[x-\frac{\lambda}{2},\,x+\frac{\lambda}{2}\right]\times\left[y-\frac{\lambda}{2},\,y+\frac{\lambda}{2}\right],\
\text{where $\lambda=\lambda(x,y).$}$ (3.1.2)

Note part of the square $E(x,y)$ may be in the second quadrant, although that
part of the square may or may not be included in our integration in Lemma 3.2
depending on how we apply our theory in the applications. Regardless, we first
study the properties of the whole square $E(x,y)$ in the following lemma,
which explains the relationship between the area of overlapping squares and
the length function $\lambda$.

######  Lemma 3.1.

For any two points $(x,y),(a,b)\in\Omega$, let $E=E(x,y)$ and $E_{0}=E(a,b)$
be two squares with side lengths $\lambda=\lambda(x,y)$ and
$\lambda_{0}=\lambda(a,b)$, centered at $(x,y)$ and $(a,b)$ respectively. We
have the following statements:

  1. (1)

Area$(E\cap E_{0})>0$ if and only if $\|(x,y)-(a,b)\|_{\infty}\ <\
\frac{\lambda+\lambda_{0}}{2};$

  2. (2)

Area$(E\cup E_{0})=\lambda^{2}+{\lambda}_{0}^{2}$ whenever
$\|(x,y)-(a,b)\|_{\infty}\ \geq\ \frac{\lambda+\lambda_{0}}{2};$

  3. (3)

If Area$(E\cap E_{0})>0$, then

$\displaystyle{\rm Area}(E\cap E_{0})$
$\displaystyle=\bigg{(}\frac{\lambda+\lambda_{0}}{2}-|x-a|\bigg{)}\bigg{(}\frac{\lambda+\lambda_{0}}{2}-|y-b|\bigg{)}\leq\left(\frac{\lambda+\lambda_{0}}{2}\right)^{2};$

  4. (4)

If Area$(E\cap E_{0})>0$, then

$\displaystyle{\rm Area}(E\setminus E_{0})+{\rm Area}(E_{0}\setminus E)$
$\displaystyle\leq\frac{|\lambda-\lambda_{0}|^{2}}{2}+2(\lambda+\lambda_{0})\|(x,y)-(a,b)\|_{\infty}.$

###### Proof.

The proofs of Statements (1), (2) and (3) follow directly from visualizing two
overlapping and two non-overlapping squares. To establish Statement (4),
observe that whenever $\mbox{Area}(E\cap E_{0})>0$, we have

$\displaystyle\mbox{Area}(E\setminus E_{0})$
$\displaystyle=\mbox{Area}(E)-\mbox{Area}(E\cap E_{0})$
$\displaystyle=\lambda^{2}-\bigg{(}\frac{\lambda+\lambda_{0}}{2}-|x-a|\bigg{)}\bigg{(}\frac{\lambda+\lambda_{0}}{2}-|y-b|\bigg{)}$
$\displaystyle\leq\lambda^{2}-\left(\frac{\lambda+\lambda_{0}}{2}\right)^{2}+\frac{\lambda+\lambda_{0}}{2}(|x-a|+|y-b|)$
$\displaystyle\leq\lambda^{2}-\left(\frac{\lambda+\lambda_{0}}{2}\right)^{2}+(\lambda+\lambda_{0})\|(x,y)-(a,b)\|_{\infty}.$

By interchanging the roles of $E$ and $E_{0}$ in the above inequality, we have

$\displaystyle\mbox{Area}(E_{0}\setminus E)$
$\displaystyle\leq\lambda_{0}^{2}-\left(\frac{\lambda+\lambda_{0}}{2}\right)^{2}+(\lambda+\lambda_{0})\|(x,y)-(a,b)\|_{\infty}.$

Therefore,

$\displaystyle\mbox{Area}(E\setminus E_{0})+\mbox{Area}(E_{0}\setminus E)$
$\displaystyle\leq\lambda^{2}+\lambda_{0}^{2}-2\left(\frac{\lambda+\lambda_{0}}{2}\right)^{2}+2(\lambda+\lambda_{0})\|(x,y)-(a,b)\|_{\infty}$
$\displaystyle=\frac{(\lambda-\lambda_{0})^{2}}{2}+2(\lambda+\lambda_{0})\|(x,y)-(a,b)\|_{\infty},$

which concludes the proof of Lemma 3.1. ∎

Based on the definition of the length function $\lambda(x,y)$, each square
$E(x,y)$ is contained in $[-\beta,\alpha+\beta]\times[0,2\beta]$, and hence
the union
$\bigcup_{(x,y)\in\Omega}E(x,y)\subseteq[-\beta,\alpha+\beta]\times[0,2\beta]$.
By letting $\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ and letting
$w:\Omega^{\prime}\longrightarrow[0,\infty)$ be a continuous nonnegative
weight function such that

$\displaystyle w(x,y)>0,\mbox{ whenever $y>0$,}$ (3.1.3)

we consider the weighted area measure $wdA$ on $\Omega^{\prime}$ and its
corresponding Hilbert space $L^{2}(wdA)$. Any function
$f:\Omega\longrightarrow\mathbb{R}$ induces a map
$\widetilde{f}:P(\Omega)\longrightarrow L^{2}(wdA)$ by

$\displaystyle\widetilde{f}(D)=\sum_{(a,b)\in D}f(a,b)\chi_{E(a,b)}$ (3.1.4)

where $\chi_{E(a,b)}$ is the characteristic function of the square $E(a,b)$
centered at $(a,b)$ as defined in (3.1.2). Since
$P(\Omega)=\bigcup_{N=1}^{\infty}P_{N}(\Omega)$, the induced map
$\widetilde{f}$ can be restricted to $P_{N}(\Omega)$ for any integer $N\geq
1$. The vector $\widetilde{f}(D)$ in the Hilbert space $L^{2}(wdA)$ is called
a persistence block, as we 3-dimensionally visualize a linear combination of
$\chi_{E(a,b)}$ in (3.1.4) as a block. The map $\widetilde{f}$ is called the
persistence block transformation induced by $f$, or simply the induced
persistence block map if the function $f$ is understood. Taking
$\widetilde{f}(D)$ as a function in $L^{2}(wdA)$, we generally integrate the
persistence block $\widetilde{f}(D)$ over the rectangle $\Omega^{\prime}$.
However, there are instances within the theory or within the applications
where one may choose to integrate the persistence block $\widetilde{f}(D)$
over the smaller rectangle $\Omega$. In these instances, we only require the
nonnegative weight function $w$ to be defined on the rectangle $\Omega$ with
$w$ continuously extended to $\Omega^{\prime}$ as necessary.

In Theorem 3.3 we establish a necessary and sufficient condition for an
induced persistence block map $\widetilde{f}$ to be continuous, in the case
when $f$ is bounded. For that we need the following lemma, noting that
$\Omega$ and $\Omega^{\prime}$ are compact and
$\lambda:\Omega\longrightarrow[0,\infty)$ and
$w:\Omega^{\prime}\longrightarrow[0,\infty)$ are continuous, and hence the
functions $\lambda$ and $w$ are necessarily bounded.

######  Lemma 3.2.

Let $\Omega=[0,\alpha]\times[0,\beta]$ and
$\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ with $\alpha,\beta>0$.
Suppose $f:\Omega\longrightarrow\mathbb{R}$ is bounded and $M$ is a positive
scalar such that

$\displaystyle\|f\|_{\infty}\leq M,\ \|\lambda\|_{\infty}\leq M,\mbox{ and
}\|w\|_{\infty}\leq M.$ (3.1.5)

Then we have the following estimations for any two points $(x,y),(a,b)$ in
$\Omega$.

  1. (1)

If ${\rm Area}(E(x,y)\cap E(a,b))>0$, then

$\displaystyle\ \int_{\Omega}|f(x,y)\chi_{E(x,y)}-f(a,b)\chi_{E(a,b)}|^{2}wdA$
$\displaystyle\leq$ $\displaystyle\
\int_{\Omega^{\prime}}|f(x,y)\chi_{E(x,y)}-f(a,b)\chi_{E(a,b)}|^{2}wdA$
$\displaystyle\leq$ $\displaystyle\
4M^{4}\bigg{(}|f(x,y)-f(a,b)|+|\lambda(x,y)-\lambda(a,b)|+\|(x,y)-(a,b)\|_{\infty}\bigg{)}.$

  2. (2)

If ${\rm Area}(E(x,y)\cap E(a,b))=0$, then

$\displaystyle\int_{\Omega}|f(x,y)\chi_{E(x,y)}-f(a,b)\chi_{E(a,b)}|^{2}wdA$
$\displaystyle\leq\int_{\Omega^{\prime}}|f(x,y)\chi_{E(x,y)}-f(a,b)\chi_{E(a,b)}|^{2}wdA$
$\displaystyle\leq 2M^{4}\|(x,y)-(a,b)\|_{\infty}.$

###### Proof.

For convenience, let $E=E(x,y)$, $\lambda=\lambda(x,y)$ and $E_{0}=E(a,b)$,
$\lambda_{0}=\lambda(a,b)$. To establish the first inequalities in Statements
(1) and (2), note that $\Omega\subseteq\Omega^{\prime}$ and the integrand for
the two integrals is nonnegative.

Focusing on the second half of the inequalities, observe that

$\displaystyle\left|f(x,y)-f(a,b)\right|^{2}$
$\displaystyle\leq(|f(x,y)|+|f(a,b)|)\left|f(x,y)-f(a,b)\right|$ (3.1.6)
$\displaystyle\leq 2M\left|f(x,y)-f(a,b)\right|.$

By repeated applications of the inequalities in (3.1.5), we have

$\displaystyle\int_{\Omega^{\prime}}\left|f(x,y)\chi_{E}-f(a,b)\chi_{E_{0}}\right|^{2}wdA$
(3.1.7) $\displaystyle=$ $\displaystyle\int_{E\cap
E_{0}}\left|f(x,y)-f(a,b)\right|^{2}wdA+\int_{E\setminus
E_{0}}\left|f(x,y)\right|^{2}wdA+\int_{E_{0}\setminus
E}\left|f(a,b)\right|^{2}wdA$ $\displaystyle\leq$ $\displaystyle\
\left|f(x,y)-f(a,b)\right|^{2}\int_{E\cap E_{0}}wdA+M^{2}\int_{E\setminus
E_{0}}wdA+M^{2}\int_{E_{0}\setminus E}wdA$ $\displaystyle\leq$ $\displaystyle\
M\left|f(x,y)-f(a,b)\right|^{2}\mbox{Area}(E\cap
E_{0})+M^{3}\mbox{Area}(E\setminus E_{0})+M^{3}\mbox{Area}(E_{0}\setminus E)$
$\displaystyle\leq$ $\displaystyle\
2M^{2}\left|f(x,y)-f(a,b)\right|\mbox{Area}(E\cap
E_{0})+M^{3}\bigg{(}\mbox{Area}(E\setminus E_{0})+\mbox{Area}(E_{0}\setminus
E)\bigg{)},\ \text{by (\ref{0617.A.3}).}$

If $\mbox{Area}(E\cap E_{0})>0$, by Statement (3) of Lemma 3.1, we have

$\displaystyle\mbox{Area}(E\cap E_{0})$
$\displaystyle\leq\bigg{(}\frac{\lambda+\lambda_{0}}{2}\bigg{)}^{2}\leq
M^{2},\mbox{ by (\ref{0617.A.1})}.$ (3.1.8)

Similarly, by Statement (4) of Lemma 3.1, we have

$\displaystyle\mbox{Area}(E\setminus E_{0})+\mbox{Area}(E_{0}\setminus E)$
$\displaystyle\leq\frac{|\lambda-\lambda_{0}|^{2}}{2}+2(\lambda+\lambda_{0})\|(x,y)-(a,b)\|_{\infty}$
(3.1.9)
$\displaystyle\leq\frac{(\lambda+\lambda_{0})|\lambda-\lambda_{0}|}{2}+2(\lambda+\lambda_{0})\|(x,y)-(a,b)\|_{\infty}$
$\displaystyle\leq M|\lambda-\lambda_{0}|+4M\|(x,y)-(a,b)\|_{\infty},\mbox{ by
(\ref{0617.A.1})}$
$\displaystyle=4M\bigg{(}|\lambda-\lambda_{0}|+\|(x,y)-(a,b)\|_{\infty}\bigg{)}.$

Substituting (3.1.8) and (3.1.9) into (3.1.7) yields

$\displaystyle\
\int_{\Omega^{\prime}}\left|f(x,y)\chi_{E}-f(a,b)\chi_{E_{0}}\right|^{2}wdA$
(3.1.10) $\displaystyle\leq$ $\displaystyle\
2M^{4}\left|f(x,y)-f(a,b)\right|+4M^{4}\bigg{(}|\lambda-\lambda_{0}|+\|(x,y)-(a,b)\|_{\infty}\bigg{)}$
$\displaystyle=$ $\displaystyle\
4M^{4}\bigg{(}\left|f(x,y)-f(a,b)\right|+|\lambda-\lambda_{0}|+\|(x,y)-(a,b)\|_{\infty}\bigg{)}.$

If $\mbox{Area}(E\cap E_{0})=0$, then

$\displaystyle\mbox{Area}(E\setminus E_{0})+\mbox{Area}(E_{0}\setminus E)$
$\displaystyle=\mbox{Area}(E)+\mbox{Area}(E_{0})$ (3.1.11) $\displaystyle=\
\lambda^{2}+\lambda_{0}^{2}$ $\displaystyle\leq M\lambda+M\lambda_{0},\mbox{
by (\ref{0617.A.1})}$ $\displaystyle=M(\lambda+\lambda_{0})$
$\displaystyle\leq 2M\|(x,y)-(a,b)\|_{\infty},\text{ by Statement (2) of Lemma
\ref{0511.A}}.$

Combining the above inequality with (3.1.7), we get

$\displaystyle\
\int_{\Omega^{\prime}}\left|f(x,y)\chi_{E}-f(a,b)\chi_{E_{0}}\right|^{2}wdA$
$\displaystyle\leq$ $\displaystyle\ M^{3}(\mbox{Area}(E\setminus
E_{0})+\mbox{Area}(E_{0}\setminus E))$ $\displaystyle\leq$ $\displaystyle\
2M^{4}\|(x,y)-(a,b)\|_{\infty},$

which concludes the proof of the lemma. ∎

In order to establish a necessary and sufficient condition to ensure the
continuity of an induced persistence block map, we utilize the $p$-Wasserstein
distance for the metric on the space of persistence diagrams. For $1\leq
p<\infty$, the $p$-Wasserstein distance $W_{p}(D,D_{0})$ between two
persistence diagrams $D,D_{0}$ in $P_{N}(\Omega)$ is defined by

$\displaystyle W_{p}(D,D_{0})=\inf\left\\{\biggl{(}\sum_{(a,b)\in
D}\|\gamma(a,b)-(a,b)\|_{\infty}^{p}\biggr{)}^{1/p}\bigg{|}\mbox{
$\gamma:D\longrightarrow D_{0}$ is a bijection}\right\\}.$ (3.1.12)

where additional points along the $x$-axis may be added to $D$ or $D_{0}$ in
order to define a bijection $\gamma$. For $0<p<1$, we similarly define the
$p$-Wasserstein distance $W_{p}(D,D_{0})$ by

$\displaystyle W_{p}(D,D_{0})=\inf\left\\{\sum_{(a,b)\in
D}\|\gamma(a,b)-(a,b)\|_{\infty}^{p}\bigg{|}\mbox{ $\gamma:D\longrightarrow
D_{0}$ is a bijection}\right\\},$ (3.1.13)

To explain the difference between the two formulas in (3.1.12) and (3.1.13),
let $x=(x_{1},x_{2},\dots,x_{n})$ and $y=(y_{1},y_{2},\dots,y_{n})$ be two
vectors in $\mathbb{R}^{n}$. The formula
$d_{p}(x,y)=\left(\sum_{j=1}^{n}|x_{j}-y_{j}|^{p}\right)^{1/p}$ defines a
metric on $\mathbb{R}^{n}$ only when $1\leq p<\infty$. However, in the case
when $0<p<1$, the above formula fails to define a metric. On the positive
side, by modifying the formula to
$d_{p}(x,y)=\sum_{j=1}^{n}|x_{j}-y_{j}|^{p}$, we generate a metric on
$\mathbb{R}^{n}$; see Rudin [23, Section 1.47, page 36]. Thus in both (3.1.12)
and (3.1.13), we take the infimum of a metric.

Before we state the main theorem of the current subsection, we remark that the
weight function $w$ is defined on $\Omega^{\prime}$, and hence the norm of the
Hilbert space $L^{2}(wdA)$ is given by

$\displaystyle\|f\|_{2}=\left(\int_{\Omega^{\prime}}|f|^{2}wdA\right)^{1/2}\mbox{
for any $f\in L^{2}(wdA)$.}$ (3.1.14)

######  Theorem 3.3.

Let $0<p<\infty$ and $N\geq 1$. For $\alpha,\beta>0$, let
$\Omega=[0,\alpha]\times[0,\beta]$. A bounded function
$f:\Omega\longrightarrow\mathbb{R}$ is continuous on
$\Omega\setminus\\{(x,0):0\leq x\leq\alpha\\}$ if and only if the induced
persistence block map $\widetilde{f}:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$
is continuous with respect to the p-Wasserstein distance $W_{p}$ on
$P_{N}(\Omega)$ and the norm $\|\cdot\|_{2}$ on the Hilbert space
$L^{2}(wdA)$.

###### Proof.

Suppose the induced persistence block map
$\widetilde{f}:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ is continuous. To show
the underlying function $f:\Omega\longrightarrow\mathbb{R}$ is continuous on
$\Omega\setminus\\{(x,0):0\leq x\leq\alpha\\}$, let $(a,b)\in\Omega$ with
$b>0$ and let $\epsilon>0$. We are to show there is a $\delta>0$ such that

$\displaystyle|f(x,y)-f(a,b)|<\epsilon,\mbox{ \ whenever $(x,y)\in\Omega$ with
$\|(x,y)-(a,b)\|_{\infty}<\delta$.}$

To this end, by the continuity of the length function $\lambda$ and the
assumption that $b>0$, there exists a real number $\delta_{1}$ with
$0<\delta_{1}<\lambda(a,b)/4$ such that

$\displaystyle\mbox{Area}(E(x,y)\cap
E(a,b))\geq\frac{1}{2}\mbox{Area}(E(a,b))>0,\ \text{whenever
$\|(x,y)-(a,b)\|_{\infty}<\delta_{1}$. }$ (3.1.15)

Consider the compact square $K$ centered at $(a,b)$ given by

$\displaystyle K=\\{(x,y)\in\Omega:\|(x,y)-(a,b)\|_{\infty}\leq\delta_{1}\\}.$

Since $0<\delta_{1}<\lambda(a,b)/4$, the compact square $K$ lies entirely
above the $x$-axis, and hence $0<\lambda(x,y)<2y$ for all $(x,y)\in K$, by our
definition of the length function $\lambda$. Furthermore, since $\lambda$ is
continuous, the function $(x,y)\mapsto y-\frac{\lambda(x,y)}{2}$ is continuous
and stays positive on the compact set $K$. Thus, the set $\bigcup_{(x,y)\in
K}E(x,y)$ is contained in a compact set that lies entirely above the $x$-axis.
Thus, by the continuity the functions $w$ and $\lambda$, there is a real
number $\eta>0$ such that

$\displaystyle\eta$ $\displaystyle<\min\\{\lambda(x,y):(x,y)\in K\\},\mbox{
and }$ (3.1.16) $\displaystyle\eta$
$\displaystyle<\min\left\\{w(u,v):(u,v)\in\bigcup_{(x,y)\in
K}E(x,y)\right\\}.$ (3.1.17)

By assumption, the induced persistence block map
$\widetilde{f}:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ is continuous at the
persistence diagram $D_{0}=\\{(a,b)\\}$ in $P_{N}(\Omega)$ that consists of
the single point $(a,b)$. Thus, there exists a $\delta_{2}>0$ such that

$\displaystyle\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}<\frac{\eta^{3/2}\epsilon}{\sqrt{2}},\
\text{whenever $W_{p}(D,D_{0})<\delta_{2}$.}$ (3.1.18)

Select a $\delta>0$ such that
$0<\max\\{\delta,\delta^{p}\\}<\min\\{\delta_{1},\delta_{2}\\}$. For any
$(x,y)\in\Omega$ with

$\displaystyle\|(x,y)-(a,b)\|_{\infty}<\delta,$ (3.1.19)

consider the persistence diagram $D=\\{(x,y)\\}\in P_{N}(\Omega)$. Observe
that

$\displaystyle W_{p}(D,D_{0})$
$\displaystyle\leq\begin{cases}\|(x,y)-(a,b)\|_{\infty},&\text{ if $1\leq
p<\infty$}\\\ \|(x,y)-(a,b)\|_{\infty}^{p},&\text{ if $0<p<1$}\end{cases}$
$\displaystyle<\max\\{\delta,\delta^{p}\\}$ $\displaystyle<\delta_{2}.$

Hence, the upper estimation for
$\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}$ in (3.1.18) holds for any
$(x,y)\in\Omega$ satisfying (3.1.19). For a lower estimation, observe that

$\displaystyle\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}^{2}$
$\displaystyle=\int_{\Omega^{\prime}}\left|f(x,y)\chi_{E(x,y)}-f(a,b)\chi_{E(a,b)}\right|^{2}wdA$
(3.1.20) $\displaystyle\geq\int_{E(x,y)\cap
E(a,b)}\left|f(x,y)-f(a,b)\right|^{2}wdA$
$\displaystyle\geq\eta\mbox{Area}(E(x,y)\cap
E(a,b))\left|f(x,y)-f(a,b)\right|^{2},\mbox{ by (\ref{0517.A.3})}$
$\displaystyle\geq\frac{\eta}{2}\mbox{Area}(E(a,b))\left|f(x,y)-f(a,b)\right|^{2},\mbox{
by (\ref{0517.A.1.1}) }$
$\displaystyle\geq\frac{\eta^{3}}{2}\left|f(x,y)-f(a,b)\right|^{2},\mbox{ by
(\ref{0517.A.2}). }$

Combining estimations in (3.1.18) and (3.1.20) yields

$\displaystyle\left|f(x,y)-f(a,b)\right|$
$\displaystyle\leq\frac{\sqrt{2}}{\eta^{3/2}}\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}<\frac{\sqrt{2}}{\eta^{3/2}}\cdot\frac{\epsilon\eta^{3/2}}{\sqrt{2}}=\epsilon,$

which completes the proof of the continuity of the function $f$.

Conversely, suppose the bounded function $f:\Omega\longrightarrow\mathbb{R}$
is continuous on $\Omega\setminus\\{(x,0):0\leq x\leq\alpha\\}$. To prove the
induced persistence block map $\widetilde{f}:P_{N}(\Omega)\longrightarrow
L^{2}(wdA)$ is continuous, select a persistence diagram $D_{0}\in
P_{N}(\Omega)$ and an $\epsilon>0$. We are to show there is a real number
$\delta>0$ such that

$\displaystyle\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}<\epsilon,\mbox{ \
whenever $D\in P_{N}(\Omega)$ with $W_{p}(D,D_{0})<\delta$.}$ (3.1.21)

To this end, since the function $f$ is bounded and the functions $\lambda,w$
are continuous on the compact sets $\Omega=[0,\alpha]\times[0,\beta]$ and
$\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ respectively, we may
select a constant $M>0$ such that

$\displaystyle\|f\|_{\infty}$
$\displaystyle=\sup\\{|f(x,y)|:(x,y)\in\Omega\\}\leq M,$ (3.1.22)
$\displaystyle\|\lambda\|_{\infty}$
$\displaystyle=\sup\\{\lambda(x,y):(x,y)\in\Omega\\}\leq M,$ (3.1.23)
$\displaystyle\|w\|_{\infty}$
$\displaystyle=\sup\left\\{w(x,y):(x,y)\in\Omega^{\prime}\right\\}\leq M.$
(3.1.24)

Furthermore, by the continuity of the length function $\lambda$ on the compact
set $\Omega$, the function $\lambda$ is uniformly continuous on $\Omega$.
Thus, there exists a real number $\delta_{1}>0$ such that

$\displaystyle|\lambda(x,y)-\lambda(a,b)|<\frac{\epsilon^{2}}{192M^{4}N^{2}},\
\text{whenever $\|(x,y)-(a,b)\|_{\infty}<\delta_{1}$. }$ (3.1.25)

Next, by the definition of the length function in (3.1.1), select a real
number $c$ such that

$\displaystyle 0<c<\min\\{\lambda(a,b):(a,b)\in D_{0}\mbox{ with
}b>0\\}<2\beta.$ (3.1.26)

Likewise by the continuity on a compact set, the function $f$ is uniformly
continuous on the compact set $[0,\alpha]\times[\frac{c}{2},\beta]$, and so
there exists a real number $\delta_{2}>0$ such that

$\displaystyle|f(x,y)-f(a,b)|<\frac{\epsilon^{2}}{192M^{4}N^{2}}.$ (3.1.27)

whenever $(x,y),(a,b)\in[0,\alpha]\times[\frac{c}{2},\beta]$ with
$\|(x,y)-(a,b)\|_{\infty}<\delta_{2}$. Select a real number $\delta>0$ such
that

$\displaystyle
0<\max\\{\delta,\delta^{1/p}\\}<\min\left\\{\frac{c}{2},\,\delta_{1},\,\delta_{2},\,\frac{\epsilon^{2}}{192M^{4}N^{2}}\right\\}.$
(3.1.28)

To show this is the desired $\delta$ for (3.1.21) to hold, let $D$ be any
persistence diagram in $P_{N}(\Omega)$ with

$\displaystyle W_{p}(D,D_{0})<\delta.$ (3.1.29)

Since $D$, $D_{0}$ are in $P_{N}(\Omega)$, the $p$-Wasserstein distance
$W_{p}(D,D_{0})$ is attained for a specific bijection $\gamma$. Furthermore,
in order to have a bijection $\gamma$ between $D$ and $D_{0}$, we may need to
add at most $N$ points along the $x$-axis to each persistence diagram. Thus,
we may assume the persistence diagrams $D$, $D_{0}$ each has $m$ points with
$1\leq m\leq 2N$, and without loss of generality, we can assume the points in
$D$ and $D_{0}$ are arranged as

$\displaystyle D_{0}=\\{(a_{j},b_{j}):1\leq j\leq m\\}\mbox{ and
}D=\\{(x_{j},y_{j}):1\leq j\leq m\\}$ (3.1.30)

so that we use (3.1.12) and (3.1.13) to rewrite (3.1.29) as

$\displaystyle W_{p}(D,D_{0})$
$\displaystyle=\begin{cases}\left(\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{p}\right)^{1/p},&\text{
if $1\leq p<\infty$}\\\
\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{p},&\text{ if
$0<p<1$}\end{cases}$ (3.1.31) $\displaystyle<\delta.$

Moreover, we can assume, by the symmetry of the roles of $D$ and $D_{0}$ in
(3.1.31), there exists an integer $m_{0}$ with $1\leq m_{0}\leq m$ and

$\displaystyle b_{j}>0\mbox{ if and only if }1\leq j\leq m_{0}.$

In order to show $\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}<\epsilon$, we
first establish the following claim.

Claim. For integers $j$ with $1\leq j\leq m$, we have

$\displaystyle\int_{\Omega^{\prime}}\left|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}\right|^{2}wdA<\frac{\epsilon^{2}}{16N^{2}}.$

Proof of Claim. In the case that $1\leq m_{0}<m$, we see that if an integer
$j$ satisfies $m_{0}<j\leq m$, then we have $b_{j}=0$ and
$\lambda(a_{j},b_{j})=0$, which in turn implies
$\mbox{Area}(E(x_{j},y_{j})\cap E(a_{j},b_{j}))=0$. Thus, if $m_{0}<j\leq m$,
by Lemma 3.2, Statement (2), we have

$\displaystyle\int_{\Omega^{\prime}}\left|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}\right|^{2}wdA$
$\displaystyle\leq 2M^{4}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}$
$\displaystyle\leq 2M^{4}\max\\{\delta,\delta^{1/p}\\},\mbox{ by
(\ref{0517.A.31})}$
$\displaystyle<2M^{4}\cdot\frac{\epsilon^{2}}{192M^{4}N^{2}},\ \text{by
(\ref{0517.A.14})}$ $\displaystyle<\frac{\epsilon^{2}}{16N^{2}},$

which establishes the claim for the integers $j$ satisfying $m_{0}<j<m$.

For an integer $j$ with $1\leq j\leq m_{0}$, we have $b_{j}>0$ and by
(3.1.26), we have $\lambda(a_{j},b_{j})\geq c>0$. It follows from (3.1.1) that
$(a_{j},b_{j})\in[0,\alpha]\times[\frac{c}{2},\beta]$. Also, observe that by
(3.1.31)

$\displaystyle\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}$
$\displaystyle\leq\begin{cases}W_{p}(D,D_{0}),&\text{ if $1\leq p<\infty$}\\\
W_{p}(D,D_{0})^{1/p},&\text{ if $0<p<1$}\end{cases}$
$\displaystyle<\max\\{\delta,\delta^{1/p}\\}$
$\displaystyle\leq\frac{c}{2},\mbox{ by (\ref{0517.A.14})}$
$\displaystyle\leq\frac{\lambda(a_{j},b_{j})}{2},\mbox{ by (\ref{0517.A.9})}$
$\displaystyle\leq\frac{\lambda(x_{j},y_{j})+\lambda(a_{j},b_{j})}{2},\mbox{
because $\lambda(x_{j},y_{j})\geq 0$.}$

Hence, we have
$(x_{j},y_{j}),(a_{j},b_{j})\in[0,\alpha]\times[\frac{c}{2},\beta]$ with
$\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}<\max\\{\delta,\delta^{1/p}\\}$.
Furthermore, by Lemma 3.1, we have $\mbox{Area}(E(x_{j},y_{j})\cap
E(a_{j},b_{j}))>0$.

Returning to the integral in our claim, observe that Statement (1) of Lemma
3.2 gives us

$\displaystyle\
\int_{\Omega^{\prime}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA$
(3.1.32) $\displaystyle\leq$ $\displaystyle\
4M^{4}\bigg{(}|f(x_{j},y_{j})-f(a_{j},b_{j})|+|\lambda(x_{j},y_{j})-\lambda(a_{j},b_{j})|+\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}\bigg{)}$
$\displaystyle<$ $\displaystyle\
4M^{4}\left(\frac{\epsilon^{2}}{192M^{4}N^{2}}+\frac{\epsilon^{2}}{192M^{4}N^{2}}+\frac{\epsilon^{2}}{192M^{4}N^{2}}\right),\mbox{
by (\ref{0517.A.13}), (\ref{0517.A.12.1}), (\ref{0517.A.14})}$
$\displaystyle=$ $\displaystyle\frac{\epsilon^{2}}{16N^{2}},$

which concludes the proof of our claim. $\Box$

Using the Claim, we now finish the proof of the theorem by showing
$\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}<\epsilon$. Observe that

$\displaystyle\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}$
$\displaystyle=\left\|\sum_{j=1}^{m}f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-\sum_{j=1}^{m}f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}\right\|_{2}$
$\displaystyle\leq\sum_{j=1}^{m}\|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}\|_{2}$
$\displaystyle=\sum_{j=1}^{m}\left(\int_{\Omega^{\prime}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA\right)^{1/2}.$

In the case that $m_{0}<m$, we use our claim to get

$\displaystyle\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}$
$\displaystyle<\sum_{j=1}^{m_{0}}\frac{\epsilon}{4N}+\sum_{j=m_{0}+1}^{m}\frac{\epsilon}{4N}$
(3.1.33) $\displaystyle=m_{0}\frac{\epsilon}{4N}+(m-m_{0})\frac{\epsilon}{4N}$
$\displaystyle\leq 2N\frac{\epsilon}{4N}+2N\frac{\epsilon}{4N}$
$\displaystyle=\epsilon.$

Lastly, in the case that $m_{0}=m$, we do not have the second summand on the
right-hand side of (3.1.33), but the same argument gives
$\|\widetilde{f}(D)-\widetilde{f}(D_{0})\|_{2}<\epsilon$. ∎

A continuous function $f:[0,\infty)\times[0,\infty)\longrightarrow\mathbb{R}$
is always bounded on the compact set $\Omega=[0,\alpha]\times[0,\beta]$ with
$\alpha,\beta>0$ and continuous on the set $\Omega\setminus\\{(x,0):0\leq
x\leq\alpha\\}$. Thus, Theorem 3.3 yields the following corollary.

######  Corollary 3.4.

Let $0<p<\infty$ and let
$f:[0,\infty)\times[0,\infty)\longrightarrow\mathbb{R}$ be a continuous
function. For any integer $N\geq 1$ and $\Omega=[0,\alpha]\times[0,\beta]$
with $\alpha,\beta>0$, the induced persistence block map
$\widetilde{f}:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ is continuous with
respect to the $p$-Wasserstein distance $W_{p}$ on $P_{N}(\Omega)$ and the
norm $\|\cdot\|_{2}$ of $L^{2}(wdA)$.

The corollary provides a theoretical result on the continuity of the induced
persistence block map. On a practical note, we infer from the corollary that
the computation of the vector $\widetilde{f}(D)$ for any given persistence
diagram $D$ is robust with respect to the noises that incur in $D$.

###  3.2 Identification Map

We continue the study of induced persistence block maps
$\widetilde{f}:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ in the special case
that the underlying function $f$ is simply the constant function $1$. In this
situation, the induced persistence block map becomes what we call the
identification map $\rho:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ given by

$\displaystyle\rho(D)=\sum_{(a,b)\in D}\chi_{E(a,b)}.$ (3.2.1)

Since the constant function is continuous on $\Omega$, by Theorem 3.3, the
identification map $\rho:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ is
continuous with respect to the p-Wasserstein distance $W_{p}$ for any $p$ with
$0<p<\infty$. When the length function $\lambda$ is Lipschitz, the
identification map $\rho$ is stable with respect to the $1/2$-Wasserstein
distance.

While one may be more familiar with the $2$-Wasserstein distance, the
$1/2$-Wasserstein distance also naturally arises due to the square root in the
definition of the Hilbert space norm on $L^{2}(wdA)$. In some situations, the
$1/2$-Wasserstein distance provides a tighter bound than the $2$-Wasserstein
distance. For an illustration, take two persistence diagrams
$D=\\{(1,1),(20,5)\\}$ and $D_{0}=\\{(2,2),(20,1)\\}$. Observe that

$\displaystyle W_{1/2}(D,D_{0})$
$\displaystyle=\|(1,1)-(2,2)\|_{\infty}^{1/2}+\|(20,5)-(20,1)\|_{\infty}^{1/2}$
$\displaystyle=\sqrt{1}+\sqrt{4}=3,$ $\displaystyle W_{2}(D,D_{0})$
$\displaystyle=\left(\|(1,1)-(2,2)\|_{\infty}^{2}+\|(20,5)-(20,1)\|_{\infty}^{2}\right)^{1/2}$
$\displaystyle=\sqrt{1^{2}+4^{2}}=\sqrt{17}\approx 4.12,$

and so $W_{1/2}(D,D_{0})<W_{2}(D,D_{0})$. Using the norm notation
$\|\cdot\|_{2}$ of $L^{2}(wdA)$ given by (3.1.14), we state the main theorem
of this subsection.

######  Theorem 3.5.

Let $N\geq 1$, and let $\Omega=[0,\alpha]\times[0,\beta]$ and
$\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ for $\alpha,\beta>0$.
Suppose the length function $\lambda:\Omega\longrightarrow[0,\infty)$ is
Lipschitz with respect to the sup-norm on $\Omega$ and
$w:\Omega^{\prime}\longrightarrow[0,\infty)$ satisfying (3.1.3). Then, the
identification map $\rho:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ is stable
with respect to the $1/2$-Wasserstein distance. In particular, for any
positive real number $M$ which satisfies

$\displaystyle|\lambda(x,y)-\lambda(a,b)|\leq M\|(x,y)-(a,b)\|_{\infty}\
\text{whenever $(x,y),(a,b)\in\Omega$, and }$ (3.2.2)
$\displaystyle\|\lambda\|_{\infty}\leq M\ \text{and}\ \|w\|_{\infty}\leq M,$
(3.2.3)

we have

$\displaystyle\left(\int_{\Omega}|\rho(D)-\rho(D_{0})|^{2}wdA\right)^{1/2}\leq\|\rho(D)-\rho(D_{0})\|_{2}\leq
2M\sqrt{M+1}\ W_{1/2}(D,D_{0})$ (3.2.4)

for any $D,D_{0}\in P_{N}(\Omega)$.

###### Proof.

The first inequality in (3.2.4) is obvious because
$\Omega\subseteq\Omega^{\prime}$. To establish the second inequality (3.2.4),
let $D$, $D_{0}$ be any two persistence diagrams in $P_{N}(\Omega)$. Following
the same reasons as in the proof of Theorem 3.3, the $1/2$-Wasserstein
distance $W_{1/2}(D,D_{0})$ is attained for a specific bijection
$\gamma:D\longrightarrow D_{0}$. For this particular bijection $\gamma$ to
exist, we may need to add at most $N$ points along the $x$-axis to each
persistence diagram. Hence, we may assume the persistence diagrams $D$,
$D_{0}$ each has $m$ points with $1\leq m\leq 2N$, and without loss of
generality, the points in $D$ and $D_{0}$ are arranged such that

$\displaystyle D=\\{(x_{j},y_{j}):1\leq j\leq m\\}\mbox{ and
}D_{0}=\\{(a_{j},b_{j}):1\leq j\leq m\\},$

and

$\displaystyle
W_{1/2}(D,D_{0})=\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{1/2}.$

In order to estimate $\|\rho(D)-\rho(D_{0})\|_{2}$, observe that

$\displaystyle\|\rho(D)-\rho(D_{0})\|_{2}=$ $\displaystyle\
\left\|\sum_{j=1}^{m}\chi_{E(x_{j},y_{j})}-\sum_{j=1}^{m}\chi_{E(a_{j},b_{j})}\right\|_{2}$
(3.2.5) $\displaystyle\leq$ $\displaystyle\
\sum_{j=1}^{m}\|\chi_{E(x_{j},y_{j})}-\chi_{E(a_{j},b_{j})}\|_{2}$
$\displaystyle=$ $\displaystyle\
\sum_{j=1}^{m}\bigg{(}\int_{\Omega^{\prime}}|\chi_{E(x_{j},y_{j})}-\chi_{E(a_{j},b_{j})}|^{2}wdA\bigg{)}^{1/2}$
$\displaystyle=$ $\displaystyle\
\sum_{j=1}^{m}\bigg{(}\int_{E(x_{j},y_{j})\setminus
E(a_{j},b_{j})}\,wdA+\int_{E(a_{j},b_{j})\setminus
E(x_{j},y_{j})}\,wdA\bigg{)}^{1/2}$ $\displaystyle\leq$ $\displaystyle\
\sum_{j=1}^{m}\bigg{(}M\mbox{Area}(E(x_{j},y_{j})\setminus
E(a_{j},b_{j}))+M\mbox{Area}(E(a_{j},b_{j})\setminus
E(x_{j},y_{j}))\bigg{)}^{1/2}$ $\displaystyle=$ $\displaystyle\
\sum_{j=1}^{m}\sqrt{M}\bigg{(}\mbox{Area}(E(x_{j},y_{j})\setminus
E(a_{j},b_{j}))+\mbox{Area}(E(a_{j},b_{j})\setminus
E(x_{j},y_{j}))\bigg{)}^{1/2}.$

For each integer $j$ in the above summation, we either have
$\mbox{Area}(E(x_{j},y_{j})\cap E(a_{j},b_{j}))>0$ or
$\mbox{Area}(E(x_{j},y_{j})\cap E(a_{j},b_{j}))=0$. For the case when
$\mbox{Area}(E(x_{j},y_{j})\cap E(a_{j},b_{j}))>0$, by inequality (3.1.9) in
the proof of Lemma 3.2, we have

$\displaystyle\ \mbox{Area}(E(x_{j},y_{j})\setminus
E(a_{j},b_{j}))+\mbox{Area}(E(a_{j},b_{j})\setminus E(x_{j},y_{j}))$ (3.2.6)
$\displaystyle\leq$ $\displaystyle\
4M\bigg{(}|\lambda(x_{j},y_{j})-\lambda(a_{j},b_{j})|+\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}\bigg{)}$
$\displaystyle\leq$ $\displaystyle\
4M\bigg{(}M\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}+\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}\bigg{)},\mbox{
by (\ref{0526.A.1})}$ $\displaystyle=$ $\displaystyle\
4M(M+1)\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}.$

For the case when $\mbox{Area}(E(x_{j},y_{j})\cap E(a_{j},b_{j}))=0$, by
inequality (3.1.11) in the proof of Lemma 3.2, we have

$\displaystyle\ \mbox{Area}(E(x_{j},y_{j})\setminus
E(a_{j},b_{j}))+\mbox{Area}(E(a_{j},b_{j})\setminus E(x_{j},y_{j}))$ (3.2.7)
$\displaystyle\leq$ $\displaystyle\
2M\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}$ $\displaystyle\leq$
$\displaystyle\ 4M(M+1)\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}.$

Returning to the estimation of $\|\rho(D)-\rho(D_{0})\|_{2}$ in (3.2.5), the
above inequalities (3.2.6) and (3.2.7) give us

$\displaystyle\|\rho(D)-\rho(D_{0})\|_{2}$
$\displaystyle\leq\sum_{j=1}^{m}\sqrt{M}\bigg{(}\mbox{Area}(E(x_{j},y_{j})\setminus
E(a_{j},b_{j}))+\mbox{Area}(E(a_{j},b_{j})\setminus
E(x_{j},y_{j}))\bigg{)}^{1/2}$
$\displaystyle\leq\sum_{j=1}^{m}\sqrt{M}\bigg{(}4M(M+1)\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}\bigg{)}^{1/2}$
$\displaystyle=2M\sqrt{M+1}\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{1/2}$
$\displaystyle=2M\sqrt{M+1}\ W_{1/2}(D,D_{0}),$

which finishes the whole proof of the Theorem. ∎

The identification map $\rho$ is defined on $P_{N}(\Omega)$, but the
estimation in inequality (3.2.4) in Theorem 3.5 is independent of the integer
$N$. From this observation, we can say the identification map
$\rho:P(\Omega)\longrightarrow L^{2}(wdA)$ is stable with respect to the
$1/2$-Wasserstein distance. While we initially establish the stability of the
identification map $\rho$ relative to the $1/2$-Wasserstein distance, we can
extend our stability result using any $p$-Wasserstein distance for
$\frac{1}{2}<p<\infty$ when we restrict to $P_{N}(\Omega)$ for any fixed
integer $N\geq 1$. To achieve the estimation, we require Hölder’s Inequality:
for any vectors $x=(x_{1},x_{2},\dots,x_{n})$ and
$y=(y_{1},y_{2},\dots,y_{n})$ in $\mathbb{R}^{n}$, we have

$\displaystyle\sum_{j=1}^{n}|x_{j}y_{j}|\leq\left(\sum_{j=1}^{n}|x_{j}|^{q_{0}}\right)^{1/q_{0}}\left(\sum_{j=1}^{n}|y_{j}|^{p_{0}}\right)^{1/p_{0}},$
(3.2.8)

where $p_{0},q_{0}$ satisfy $1<p_{0},q_{0}<\infty$ and
$\frac{1}{p_{0}}+\frac{1}{q_{0}}=1$. Keeping Hölder’s inequality in mind, we
have the following corollary to Theorem 3.5.

######  Corollary 3.6.

Let $N\geq 1$, and let $\Omega=[0,\alpha]\times[0,\beta]$ and
$\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ for $\alpha,\beta>0$.
Suppose the length function $\lambda:\Omega\longrightarrow[0,\infty)$ and the
weight function $w:\Omega^{\prime}\longrightarrow(0,\infty)$ satisfy the
hypotheses of Theorem 3.5. The identification map
$\rho:P_{N}(\Omega)\longrightarrow L^{2}(wdA)$ satisfies

$\displaystyle\left(\int_{\Omega}|\rho(D)-\rho(D_{0})|^{2}wdA\right)^{1/2}\leq\|\rho(D)-\rho(D_{0})\|_{2}\leq\begin{cases}4NM\sqrt{M+1}\
W_{p}(D,D_{0})^{1/2},&\text{ if $1\leq p<\infty$}\\\ 4NM\sqrt{M+1}\
W_{p}(D,D_{0})^{1/2p},&\text{ if $\frac{1}{2}<p<1$}\end{cases},$

for any $D,D_{0}\in P_{N}(\Omega)$. In particular for $p=2$, we have

$\displaystyle\left(\int_{\Omega}|\rho(D)-\rho(D_{0})|^{2}wdA\right)^{1/2}\leq\|\rho(D)-\rho(D_{0})\|_{2}\leq
4NM\sqrt{M+1}\ W_{2}(D,D_{0})^{1/2}.$

###### Proof.

Let the real number $p$ satisfy $\frac{1}{2}<p<\infty$. As in the proof of
Theorem 3.5, for any two persistence diagrams $D,D_{0}\in P_{N}(\Omega)$,
select the bijection $\gamma:D\longrightarrow D_{0}$ and the integer $m$ with
$1\leq m\leq 2N$ so that we arrange the points in $D$ and $D_{0}$ as

$\displaystyle D=\\{(x_{j},y_{j}):1\leq j\leq m\\}\mbox{ and
}D_{0}=\\{(a_{j},b_{j}):1\leq j\leq m\\},$

for which

$\displaystyle
W_{p}(D,D_{0})=\begin{cases}\left(\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{p}\right)^{1/p},&\text{
if $1\leq p<\infty$}\\\
\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{p},&\text{ if
$\frac{1}{2}<p<1.$}\end{cases}$ (3.2.9)

It follows from (3.2.4) that

$\displaystyle\|\rho(D)-\rho(D_{0})\|_{2}$ $\displaystyle\leq
2M\sqrt{M+1}\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{1/2}$
(3.2.10)
$\displaystyle=\sum_{j=1}^{m}2M\sqrt{M+1}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{1/2}.$

Applying Hölder’s inequality to (3.2.10) with $p_{0}=2p>1$ and
$q_{0}=\frac{1}{1-\frac{1}{p_{0}}}=\frac{2p}{2p-1}>1$ yields

$\displaystyle\ \|\rho(D)-\rho(D_{0})\|_{2}$ $\displaystyle\leq$
$\displaystyle\
\left(\sum_{j=1}^{m}(2M\sqrt{M+1})^{q_{0}}\right)^{1/q_{0}}\left(\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{p_{0}/2}\right)^{1/p_{0}}$
$\displaystyle=$ $\displaystyle\
2m^{1/q_{0}}M\sqrt{M+1}\left(\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{p}\right)^{1/2p},\mbox{
because $p_{0}=2p$}$ $\displaystyle\leq$ $\displaystyle\
4NM\sqrt{M+1}\left(\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{p}\right)^{1/2p},\
\mbox{ because $m^{1/q_{0}}\leq(2N)^{1/q_{0}}\leq 2N$}$ $\displaystyle=$
$\displaystyle\begin{cases}4NM\sqrt{M+1}\ W_{p}(D,D_{0})^{1/2},&\text{ if
$1\leq p<\infty$}\\\ 4NM\sqrt{M+1}\ W_{p}(D,D_{0})^{1/2p},&\text{ if
$\frac{1}{2}<p<1,$}\end{cases}\ \mbox{ by (\ref{0527.A.2})}.$

∎

It follows from the above theorem and corollary that with respect to
$W_{p}(D,D_{0})$ for any $p\geq 1/2$, the computation of $\rho(D)$ for any
persistence diagram $D$ is robust with respect to the noises that incur in
$D$, under an assumption on the length function $\lambda$ and the weight
function $w$.

###  3.3 Vectorization

Now, for a persistence diagram $D$, we would like to associate every
persistence block $\widetilde{f}(D)$, which is some vector in the Hilbert
space $L^{2}(wdA)$, with a vector in $\mathbb{R}^{n}$. There are two ways to
do that, depending whether one wants to focus on
$\Omega=[0,\alpha]\times[0,\beta]$ or
$\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ as discussed in the
previous two subsections. Since both ways share the same mathematical theory,
we develop them together as one in this subsection. To this end, select a
partition $\mathcal{P}$ of either $\Omega$ or $\Omega^{\prime}$, whichever one
chooses to focus, into $n$ grid cells $S_{1},S_{2},\ldots,S_{n}$. For any
persistence diagram $D$ in $P(\Omega)$, let

$\displaystyle I_{f,S_{i}}(D)=\int_{S_{i}}\,\widetilde{f}(D)\,wdA,$

and define the vector

$\displaystyle\widetilde{f}_{\mathcal{P}}(D)=(I_{f,S_{1}}(D),I_{f,S_{2}}(D),\ldots,I_{f,S_{n}}(D))$

in $\mathbb{R}^{n}$ to be the vectorized persistence block (VPB). Thus we can
think of $\widetilde{f}_{\mathcal{P}}$ as a map
$\widetilde{f}_{\mathcal{P}}:P(\Omega)\longrightarrow\mathbb{R}^{n}$.

To ensure the theory accounts for both types of partitions, let

$\displaystyle G=\begin{cases}\Omega,&\text{ if $\mathcal{P}$ is a partition
of $\Omega$}\\\ \bigcup_{(x,y)\in\Omega}E(x,y),&\text{ if $\mathcal{P}$ is a
partition of $\Omega^{\prime}$}\end{cases}.$ (3.3.1)

Note that $\Omega\subseteq G\subseteq\Omega^{\prime}$, and the persistence
block $\widetilde{f}(D)$ is a linear combination of characteristic functions
$\chi_{E(x,y)}$ of the square $E(x,y)$ centered at $(x,y)\in\Omega$. Thus,
whether $S_{1},S_{2},\dots,S_{n}$ form a partition $\mathcal{P}$ of $\Omega$
or $\Omega^{\prime}$, we have

$\displaystyle I_{f,S_{i}}(D)$ $\displaystyle=\int_{S_{i}}\widetilde{f}(D)wdA$
(3.3.2) $\displaystyle=\int_{S_{i}\cap G}\widetilde{f}(D)wdA.$

The VPB $\widetilde{f}_{\mathcal{P}}(D)$ is a vector in $\mathbb{R}^{n}$, and
so we may use the standard norms on $\mathbb{R}^{n}$ to measure the distance
between two VPBs. In particular, for any two persistence diagrams $D,D_{0}\in
P_{N}(\Omega)$ and for $1\leq p<\infty$, we have

$\displaystyle\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{p}$
$\displaystyle=\left(\sum_{i=1}^{n}|I_{f,S_{i}}(D)-I_{f,S_{i}}(D_{0})|^{p}\right)^{1/p},$
(3.3.3)

and

$\displaystyle\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{\infty}$
$\displaystyle=\sup_{1\leq i\leq n}|I_{f,S_{i}}(D)-I_{f,S_{i}}(D_{0})|.$
(3.3.4)

The induced persistence block map $\widetilde{f}:P_{N}(\Omega)\longrightarrow
L^{2}(wdA)$ is continuous whenever the underlying function
$f:[0,\infty)\times[0,\infty)\longrightarrow\mathbb{R}$ is bounded and
continuous, by Theorem 3.3. In order to explore a similar continuity result
for the VPB relative to the Hilbert space norm $\|\cdot\|_{2}$ given in
(3.3.3) above, we establish an upper estimate for
$\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{2}$,
which involves the measure $\mu$ defined by

$\displaystyle\mu(E)=\int_{E}d\mu=\int_{E}wdA.$

######  Lemma 3.7.

Let $\Omega=[0,\alpha]\times[0,\beta]$ and
$\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ with $\alpha,\beta>0$,
let $\mathcal{P}$ be a partition of $\Omega$ or $\Omega^{\prime}$ into $n$
grid cells $S_{1},S_{2},\dots,S_{n}$, and let $G$ be the set defined in
(3.3.1). For persistence diagrams $D=\\{(x_{j},y_{j}):1\leq j\leq m\\}$ and
$D_{0}=\\{(a_{j},b_{j}):1\leq j\leq m\\}$, we have

$\displaystyle\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{2}^{2}$
$\displaystyle\leq\sum_{i=1}^{n}\mu(S_{i}\cap
G)\left[\sum_{j=1}^{m}\left(\int_{S_{i}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA\right)^{1/2}\right]^{2}.$

###### Proof.

Observe that for integers $i$ with $1\leq i\leq n$, we have

$\displaystyle\ \left|I_{f,S_{i}}(D)-I_{f,S_{i}}(D_{0})\right|^{2}$
$\displaystyle=$ $\displaystyle\ \left|\int_{S_{i}\cap
G}\widetilde{f}(D)wdA-\int_{S_{i}\cap
G}\widetilde{f}(D_{0})wdA\right|^{2},\mbox{ by (\ref{0803.2})}$
$\displaystyle\leq$ $\displaystyle\ \left(\int_{S_{i}\cap
G}|\widetilde{f}(D)-\widetilde{f}(D_{0})|wdA\right)^{2}$ $\displaystyle\leq$
$\displaystyle\ \left(\int_{S_{i}\cap G}1^{2}wdA\right)\left(\int_{S_{i}\cap
G}|\widetilde{f}(D)-\widetilde{f}(D_{0})|^{2}wdA\right),\mbox{ by the Cauchy-
Schwarz Inequality}$ $\displaystyle\leq$ $\displaystyle\ \mu(S_{i}\cap
G)\int_{S_{i}}|\widetilde{f}(D)-\widetilde{f}(D_{0})|^{2}wdA,\mbox{ because
$S_{i}\cap G\subseteq S_{i}$ and the definition of $\mu$}$ $\displaystyle\leq$
$\displaystyle\ \mu(S_{i}\cap
G)\left(\sum_{j=1}^{m}\|(f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})})\chi_{S_{i}}\|_{2}\right)^{2},\
\text{by the Triangle Inequality}$ $\displaystyle=$ $\displaystyle\
\mu(S_{i}\cap
G)\left(\sum_{j=1}^{m}\left(\int_{S_{i}}|(f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA\right)^{1/2}\right)^{2}.$

Summing over all the grid cells $S_{i}$ in the partition $\mathcal{P}$ leads
to the desired inequality. ∎

Similar to the identification map $\rho$, under some basic assumptions, we now
show the vectorized persistence block map
$\widetilde{f}_{\mathcal{P}}:P(\Omega)\longrightarrow\mathbb{R}^{n}$ is stable
relative to the $1/2$-Wasserstein distance.

######  Theorem 3.8.

Let $\Omega=[0,\alpha]\times[0,\beta]$ and
$\Omega^{\prime}=[-\beta,\alpha+\beta]\times[0,2\beta]$ with $\alpha,\beta>0$,
let $\mathcal{P}$ be a partition of $\Omega$ or $\Omega^{\prime}$ into $n$
grid cells $S_{1},S_{2},\dots,S_{n}$, and let $G$ be the set defined in
(3.3.1). Suppose the function $f:\Omega\longrightarrow\mathbb{R}$ and the
length function $\lambda:\Omega\longrightarrow[0,\infty)$ are Lipschitz. Then
the vectorized persistence block map
$\widetilde{f}_{\mathcal{P}}:P(\Omega)\longrightarrow\mathbb{R}^{n}$ is stable
relative to the $1/2$-Wasserstein distance. In particular, for any positive
real number $M$ which satisfies

$\displaystyle|f(x,y)-f(a,b)|\leq M\|(x,y)-(a,b)\|_{\infty}\ \text{for all
$(x,y),(a,b)\in\Omega$; }$ (3.3.5)
$\displaystyle|\lambda(x,y)-\lambda(a,b)|\leq M\|(x,y)-(a,b)\|_{\infty}\
\text{for all $(x,y),(a,b)\in\Omega$; }$ (3.3.6)
$\displaystyle\|f\|_{\infty}\leq M,\ \|\lambda\|_{\infty}\leq M,\ \text{and}\
\|w\|_{\infty}\leq M,$ (3.3.7)

we have

$\displaystyle\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{2}\leq
2M^{2}\sqrt{\mu(G)(2M+1)}\ W_{1/2}(D,D_{0}),$ (3.3.8)

for any $D,D_{0}\in P(\Omega)$.

###### Proof.

Assume the positive real number $M$ satisfies inequalities (3.3.5), (3.3.6)
and (3.3.7). For any two persistence diagrams $D,D_{0}\in P(\Omega)$, let $N$
be the maximum cardinality of $D$ and $D_{0}$. As in the proof of previous
results, we may assume $D,D_{0}$ each has $m$ points with $1\leq m\leq 2N$,
and the points in $D,D_{0}$ are arranged such that

$\displaystyle D=\\{(x_{j},y_{j}):1\leq j\leq m\\}\mbox{ and
}D_{0}=\\{(a_{j},b_{j}):1\leq j\leq m\\}$

and

$\displaystyle
W_{1/2}(D,D_{0})=\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{1/2}.$

As in the proof of Theorem 3.5, we have two cases based on the area of
$E(x_{j},y_{j})\cap E(a_{j},b_{j})$. For the case when
$\mbox{Area}(E(x_{j},y_{j})\cap E(a_{j},b_{j}))>0$, Lemma 3.2 and assumptions
(3.3.5), (3.3.6) give us

$\displaystyle\
\int_{S_{i}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA$
$\displaystyle\leq$ $\displaystyle\
\int_{\Omega^{\prime}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA,\mbox{
because $S_{i}\subseteq\Omega^{\prime}$}$ $\displaystyle\leq$ $\displaystyle\
4M^{4}\bigg{(}|f(x_{j},y_{j})-f(a_{j},b_{j})|+|\lambda(x_{j},y_{j})-\lambda(a_{j},b_{j})|+\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}\bigg{)}$
$\displaystyle\leq$ $\displaystyle\
4M^{4}\bigg{(}M\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}+M\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}+\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}\bigg{)}$
$\displaystyle=$ $\displaystyle\
4M^{4}(2M+1)\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}.$

For the case when $\mbox{Area}(E(x_{j},y_{j})\cap E(a_{j},b_{j}))=0$, Lemma
3.2 gives us

$\displaystyle\
\int_{S_{i}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA$
$\displaystyle\leq$ $\displaystyle\
\int_{\Omega^{\prime}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA$
$\displaystyle\leq$ $\displaystyle\ 2M^{4}\
\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}$ $\displaystyle\leq$ $\displaystyle\
4M^{4}(2M+1)\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}.$

To finish the estimation of
$\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{2}$
using Lemma 3.7, observe that

$\displaystyle\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{2}^{2}$
$\displaystyle\leq\sum_{i=1}^{n}\mu(S_{i}\cap
G)\left[\sum_{j=1}^{m}\left(\int_{S_{i}}|f(x_{j},y_{j})\chi_{E(x_{j},y_{j})}-f(a_{j},b_{j})\chi_{E(a_{j},b_{j})}|^{2}wdA\right)^{1/2}\right]^{2}$
$\displaystyle\leq\sum_{i=1}^{n}\mu(S_{i}\cap
G)\left[\sum_{j=1}^{m}\left(4M^{4}(2M+1)\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}\right)^{1/2}\right]^{2}$
$\displaystyle=4M^{4}(2M+1)\sum_{i=1}^{n}\mu(S_{i}\cap
G)\left[\sum_{j=1}^{m}\|(x_{j},y_{j})-(a_{j},b_{j})\|_{\infty}^{1/2}\right]^{2}$
$\displaystyle=4M^{4}(2M+1)\ W_{1/2}(D,D_{0})^{2}\sum_{i=1}^{n}\mu(S_{i}\cap
G)$ $\displaystyle=4M^{4}\mu(G)(2M+1)\ W_{1/2}(D,D_{0})^{2},$

which concludes the proof. ∎

Similar to the identification map $\rho$ discussed in the previous subsection,
the stability of the vectorized persistence block map
$\widetilde{f}_{\mathcal{P}}$ relative to the $1/2$-Wasserstein distance is
independent of the number of points in the persistence diagrams $D,D_{0}\in
P(\Omega)$. If we restrict our attention to the persistence diagrams in
$P_{N}(\Omega)$, the following corollary shows the vectorized persistence
block map $\widetilde{f}_{\mathcal{P}}$ is stable relative to the
$p$-Wasserstein distance for $\frac{1}{2}<p<\infty$. The proof of the
corollary utilizes Hölder’s inequality, following the exact same argument as
in the proof of Corollary 3.6, and we need not repeat the argument here.

######  Corollary 3.9.

Assuming the hypotheses of Theorem 3.8, the VPBs
$\widetilde{f}_{\mathcal{P}}(D)$ and $\widetilde{f}_{\mathcal{P}}(D_{0})$
satisfy

$\displaystyle\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{2}\leq\begin{cases}4NM^{2}\sqrt{\mu(G)(2M+1)}\
W_{p}(D,D_{0})^{1/2},&\text{ if $1\leq p<\infty$}\\\
4NM^{2}\sqrt{\mu(G)(2M+1)}\ W_{p}(D,D_{0})^{1/2p},&\text{ if
$\frac{1}{2}<p<1$}\end{cases},$

for any $D,D_{0}\in P_{N}(\Omega)$. In particular for $p=2$, we have

$\displaystyle\|\widetilde{f}_{\mathcal{P}}(D)-\widetilde{f}_{\mathcal{P}}(D_{0})\|_{2}\leq
4NM^{2}\sqrt{\mu(G)(2M+1)}\ W_{2}(D,D_{0})^{1/2}.$

The vectorization map $\widetilde{f}_{\mathcal{P}}(D)$ is stable with respect
to the $p-$Wasserstein distance, when $1/2\leq p<\infty$. The estimations
provided in Theorem 3.8 and Corollary 3.9, do not even depend on the number
$n$ of grid cells in the partition $\mathcal{P}$ of $\Omega$ or
$\Omega^{\prime}$. In practice, we can take $\mu(G)$ to be a fixed constant
and we can have a maximum number $N$ of points in a persistence diagram. Thus
there is a natural upper bound for the integer $N$ in our estimations in the
conclusion of Corollary 3.9. However, in the case for $p=1/2$, the estimation
provided in Theorem 3.8 does not depend on the integer $N$.

##  4 Simulations

In this section we present our simulation results to illustrate the utility of
vectorized persistence blocks (VPB) for classification, clustering and change
point detection tasks. In all of the simulations, our results are compared
with the corresponding results produced from persistence images (PI). Two of
the simulations closely mimic those of Sections 6.1 and 6.4 in [18]. The
simulations involve both synthetic and real datasets. On 2D datasets, we
extract topological information using the Vietoris-Rips filtration [24]. For
large 3D datasets, we resort to the $\alpha$-shape complexes [25, 1] to avoid
computational challenges that arise from using the Vietoris-Rips complexes in
that situation.

To compute VPBs, we experiment with several choices of the weight function $w$
and the side length function $\lambda$. For example, for $\lambda$ we consider
functions of the form

$\lambda=F(y;\tau,n,m)=2\tau
y\left(\frac{y}{y_{max}}\right)^{n}\left(1-\frac{y}{y_{max}}\right)^{m},$

where $\tau\in(0,1]$, $n$, $m$ are non-negative integers and $y_{max}$ is the
maximum of persistence values $y$ for the PD in question. All our results are
based on the weight function $w(x,y)=x+y$ and the side length
$\lambda=F(y;\tau,0,0)=2\tau y$, where $0<\tau\leq 1$. Note that with the
selected functions for weight and side length, VPBs can be computed
analytically. PDs are computed using the R packages TDA [26] and TDAstats
[27]. To produce PIs we utilize the R package kernelTDA [28].

###  4.1 Clustering using K-medoids on a synthetic dataset

This simulation study considers the synthetic dataset used in Sections 6.1 of
[18] which consists of 25 point clouds of size 500 sampled uniformly from each
of the following six shapes: a unit cube, a circle of radius 1/2, a sphere of
radius 1/2, three clusters with random centers in the unit cube, three
clusters each containing another set of three clusters, a torus with inner and
outer diameters of 1/2 and 1 respectively. Furthermore, Gaussian noise is
added to each point cloud at two levels: $\eta=0.05$ and $\eta=0.1$.
Topological signatures (in the form of VPBs and PIs) are computed from PDs for
homological dimensions $H_{0}$ and $H_{1}$. In this study, the goal is to
cluster the given 150 point clouds into 6 respective shape classes.

For PIs, we use the parameter settings adopted by the authors of [18]:
standard deviation $\sigma=0.1$, a $20\times 20$ grid and the linear weighting
function. VPBs are computed with the candidate values of $\tau$ ranging in
$\\{0.1,0.3,0.5,0.7,0.9\\}$ on a grid of size $6\times 6$. The final value of
$\tau$ is selected based on the internal measure of the Davies-Bouldin Index
[29]. More specifically, among the five different clustering performances (one
for each $\tau$), we choose that for which the Davies-Bouldin Index is
minimum.

Following the approach in [18], clustering is performed employing K-medoids
method [30, 31], which is a robust version of the well-known K-means algorithm
[32]. It uses actual data points as cluster centers (called medoids) instead
of centroids as in K-means and can be applied with any dissimilarity measure.
To implement K-medoids in R, we use cluster::pam() function with its default
parameter settings [33].

A $150\times 150$ matrix of pairwise distances (also called a dissimilarity
matrix) between topological signatures (VPBs or PIs) is computed and passed to
K-medoids (specifying the number of clusters to be six) which returns a list
of cluster assignments. The performance is evaluated based on the accuracy of
cluster predictions against the true cluster labels.

Comparison of the clustering accuracies is given in Figure 2. The results show
that K-medoids produces slightly better accuracies on VPBs. Also, we observe
that the computational cost associated with VPBs (0.09 and 0.188 seconds for
$H_{0}$ and $H_{1}$) is lower than that of PIs (0.615 and 0.714 seconds for
$H_{0}$ and $H_{1}$) despite running K-medoids five times for different values
of $\tau$. Here, the cost is the combined time to compute the topological
signatures and to run the K-medoids algorithm. Here, however, the observed
difference in computational cost is partially due to the grid sizes being
unequal (i.e. $6\times 6$ for VPBs and $20\times 20$ for PIs). So next, we
compare the computational times to produce VPBs and PIs computed over the same
grid. We generate 100 PDs of size $n\in\\{1000,2000,\ldots,19000,20000\\}$,
where the birth and persistence values are randomly sampled from the beta
distributions with the shape parameters $\alpha_{b}=4$, $\beta_{b}=6$ and
$\alpha_{p}=1$, $\beta_{p}=5$ respectively. We superimpose a grid of size
$10\times 10$ to compute VPS and PIs. The plot of the total computational cost
(in seconds) versus the size of PDs is given in Figure 2, which shows that
VBPs are 6-7 times faster to compute than PIs.

Cluster visualizations of VPBs reduced via principal component analysis (PCA)
are given in Figure 3.

Metric | Dim | Noise | PI | VPB
---|---|---|---|---
$L_{1}$ | $H_{0}$ | 0.05 | 0.933 | 0.947
$L_{2}$ | $H_{0}$ | 0.05 | 0.927 | 0.940
$L_{\infty}$ | $H_{0}$ | 0.05 | 0.940 | 0.947
$L_{1}$ | $H_{0}$ | 0.10 | 0.953 | 0.960
$L_{2}$ | $H_{0}$ | 0.10 | 0.953 | 0.960
$L_{\infty}$ | $H_{0}$ | 0.10 | 0.960 | 0.960
$L_{1}$ | $H_{1}$ | 0.05 | 1.000 | 1.000
$L_{2}$ | $H_{1}$ | 0.05 | 1.000 | 1.000
$L_{\infty}$ | $H_{1}$ | 0.05 | 1.000 | 1.000
$L_{1}$ | $H_{1}$ | 0.10 | 0.960 | 0.973
$L_{2}$ | $H_{1}$ | 0.10 | 0.960 | 0.973
$L_{\infty}$ | $H_{1}$ | 0.10 | 0.960 | 0.960
Figure 1: K-medoids accuracy results.

Figure 2: Comparison of computational cost.

Figure 3: Visualization of VPB-based clusters using PCA.

Finally, we explore the sensitivity of accuracies for VPBs with respect to
changes in tuning parameters. From Figures 5 and 5, we can observe that the
accuracies appear to be fairly stable with respect to $\tau$ and grid size
when the noise level $\eta=0.05$.

Figure 4: VPB accuracy as a function of $\tau$. Grid size $6\times 6$, metric
$L_{2}$.

Figure 5: VPB accuracy as a function of grid size. $\tau=0.75$, metric
$L_{2}$.

###  4.2 Clustering using K-medoids on a real dataset

In our second simulation study, we perform clustering and shape retrieval
tasks on a 3D dataset. The dataset consists of 100 point clouds derived from
3D scans of ten non-rigid toys, each having approximately 4000 points [34].
All toys are scanned in ten different poses by moving various body parts
around their joints (see Figures 7 and 7). The dataset contains only vertex
information and no triangle mesh is provided. The PDs are computed via
$\alpha$-shape filtrations built on top of the point clouds.

Figure 6: Toys used for SHREC17 dataset [34].

Figure 7: Sample point clouds from SHREC17 dataset. The points are colored
based on the $y$ and $z$ coordinates [34].

As in Section 4.1, we perform clustering using the K-medoids algorithm. For
PIs, we select a $20\times 20$ grid and the default value of $\sigma$ which is
1/2 of the maximum persistence across all PDs divided by the grid size. VPBs
are computed on a $12\times 12$ grid. All the other configurations carry over
from the previous section. In particular, similar to the previous experiment,
we perform clustering for each homological dimension. We find that VPBs and
PIs are most discriminative for dimension $H_{2}$ and therefore we report our
results only for this dimension. Figure 8 shows the accuracies achieved by the
two methods.

Figure 8: Homological dimension used is $H_{2}$.

The $100\times 100$ dissimilarity matrix used as an input for K-medoids, can
also be viewed as a result of a 3D shape retrieval problem. In shape
retrieval, as in clustering, no class information is allowed to be used to
compute this matrix. We evaluate the retrieval performance using five standard
statistics (Nearest Neighbor (NN), First-Tier (FT), Second-Tier (ST),
E-Measure (E) and Discounted Cumulative Gain (DCG) [35]) computed from the
dissimilarity matrix.

The selected dataset for our experiment was initially used for the 2017 Shape
Retrieval Contest (SHREC) track: Point-Cloud Shape Retrieval of Non-Rigid
Toys. This track invited researchers to contribute to the development of new
discriminative and efficient descriptors of 3D non-rigid shapes. A total of
eight groups participated in this content. Each group was allowed to submit up
to 6 dissimilarity matrices. We compare our results with not only the
performance of PIs but also with the best results of all participants.

Method | NN | FT | ST | E | DCG
---|---|---|---|---|---
BoW-RoPS-DMF-3 | 1.0000 | 0.9778 | 0.9978 | 0.4390 | 0.9979
BPHAPT | 0.9800 | 0.9111 | 0.9544 | 0.4273 | 0.9743
PI | 0.9700 | 0.8989 | 0.9767 | 0.4376 | 0.9740
VPB | 0.9800 | 0.8400 | 0.9567 | 0.4351 | 0.9573
MFLO-FV-IWKS | 0.8900 | 0.7911 | 0.8589 | 0.4024 | 0.9038
CDSPF | 0.9200 | 0.6744 | 0.8156 | 0.4005 | 0.8851
SnapNet | 0.8800 | 0.6633 | 0.8011 | 0.3985 | 0.8663
AlphaVol1 | 0.7900 | 0.5878 | 0.7578 | 0.3980 | 0.8145
SQFD(WKS) | 0.5400 | 0.3111 | 0.4467 | 0.2507 | 0.6032
m3DSH-2 | 0.4400 | 0.1867 | 0.2856 | 0.1932 | 0.4997
Figure 9: Retrieval results on SHREC17 [35]. The methods are ranked according
to overall performance.

We use $L_{1}$ metric to compute the dissimilarity matrices for VPBs and PIs
in homological dimension $H_{2}$. For VPBs, $\tau=0.1$ gives a slightly better
result than that of $\tau=0.3$ for which the Davies-Bouldin Index is minimum.
From Figure 9, we see that PIs slightly outperform VPBs. BPHAPT and our method
show comparable results - we have the same NN score and surpass in terms of ST
and EE scores.

Furthermore, VPBs exhibit substantial computational advantages over both BoW-
RoPS-DMF-3 and BPHAPT. For example, for BoW-RoPS-DMF-3 which utilizes the Bag-
of-Words scheme, it takes around 30 minutes for the codebook module to run,
whereas our method needs only about 48 seconds to compute all topological
summaries (46 seconds for PDs and 2 seconds for VPBs). In terms of data
preprocessing, BPHAPT has an extra step of constructing triangular meshes
using the Poisson reconstruction method, while our approach does not
preprocess the data and only uses the vertex information. Finally, our method
has simpler complexity with regard to the number of tuning parameters involved
- we only need to choose the parameter $\tau$ and the grid size.

###  4.3 Classification of parameter values of a discrete dynamical system

In this section, we explore the utility of VPBs in the context of
classification and change point detection problems. These studies are also
motivated by [18] and involve data generated by the linked twist map defined
by the dynamical system:

$\displaystyle x_{n+1}$ $\displaystyle=x_{n}+r\cdot y_{n}(1-y_{n})\hbox{ mod
1}$ $\displaystyle y_{n+1}$ $\displaystyle=y_{n}+r\cdot
x_{n+1}(1-x_{n+1})\hbox{ mod 1,}$

where $r>0$ [36]. As in [18], we consider values of $r$ ranging in
$\\{2.0,3.5,4.0,4.1,4.3\\}$, and take $\sigma=0.005$ and a $20\times 20$ grid
for the computation of PIs.

For classification, we generate 100 truncated orbits
$\\{(x_{n},y_{n}):n=0,1,\ldots,1000\\}$ for each value of $r$, where the
initial points $(x_{0},y_{0})$ are selected at random from the unit square
(see Figures 11, 13 and 13; a typical scatter plot for $r=3.5$, not shown in
the figures due to space limitations, does not exhibit any pattern and appears
to be uniformly distributed over the unit square). The figures show that the
scatter plots for $r=4.0$ and $r=4.1$ both have one hole in the center but
they slightly differ in size and that scatter plots for the same value of $r$
can show noticeable within-class variation in shape making the classification
task more challenging.

Vietoris-Rips filtrations are used to compute PDs. To curb the computation
cost to produce the PDs, we terminate the filtrations at scale value of 0.45.
In other words, simplices with at least one pairwise distance greater than
0.45 do not enter any complex in the filtration. At this scale value, the
holes in the point clouds (for $r=4.0,4.1,4.3$) are already filled by
2-simplices (triangles). Going beyond that scale value does not essentially
produce any additional useful topological information.

Interestingly, when generating the orbits for our simulation studies we notice
that for certain values of $r$, some initial points $(x_{0},y_{0})$ lead to
atypical-looking orbits (see Figure 11). Clearly, these orbits would mislead
any classification method, so we devise a simple heuristic to omit such orbits
from our analysis. We compute VPBs on a $7\times 7$ grid and let $\tau$ vary
over $\\{0.1,0.3,0.5,0.7,0.9\\}$. The optimal value of $\tau$ is selected by
10-fold cross-validation. We use standard machine learning methods such as
Random Forest, Support Vector Machines (SVM), GLMNet and Model Averaged Neural
Network to classify the orbits in the test set according to the associated
parameter value of $r$.

Figure 10: Scatter plots of sample truncated orbits corresponding to
$r=2.0,4.0,4.1,4.3$ respectively.

Figure 11: Scatter plots of atypical orbits corresponding to
$r=2.0,4.0,4.1,4.3$ respectively.

Figure 12: Scatter plots of sample truncated orbits corresponding to $r=4.3$
and different $(x_{0},y)$.

Figure 13: Scatter plots of sample truncated orbits corresponding to $r=2.0$
and different $(x_{0},y)$.

We run the classifiers on the $H_{0}$ and $H_{1}$ features as well as when
these features are concatenated. $H_{1}$-based test accuracies for 50/50 and
70/30 splits of the data and averaged over 10 independent trials are given in
Figures 15 and 15. When the $H_{0}$ and $H_{1}$ are concatenated, the
accuracies come out to be only about 1% higher than those for $H_{1}$, while
the $H_{0}$ features yield the lowest results (around 62% accuracy). From the
figures we can see that, VPBs, when compared to PIs, on average produce
slightly higher test accuracies across almost all the machine learning methods
employed.

Figure 14: 50/50 split of data. VPBs and PIs are constructed from $H_{1}$
features.

Figure 15: 70/30 split of data. VPBs and PIs are constructed from $H_{1}$
features.

Next, we present a study on change point detection. Let $r_{1}=2$,
$r_{2}=3.5$, $r_{3}=4.0$, $r_{4}=4.1$, $r_{5}=4.3$. In this study, we generate
50 truncated orbits for each value of $r$ using the linked twist map and
assume they are ordered in time so that the change points occur at times
$t=51,101,151,201$: $X_{t,r}=\\{(x_{n}(r),y_{n}(r)):n=0,1,\ldots,1000\\}$,
where $r=r_{\lceil t/50\rceil}$, $t=1,\ldots,250$ and initial conditions
$(x_{0},y_{0})$ are selected at random. To estimate the true change points, we
use a nonparametric multiple change point detection method called E-Divisive
which determines statistically significant change points based on a
permutation test [37]. Our choice is motivated by the fact that the
distribution of the topological features remains unknown. We input $VPB_{t,r}$
and $PI_{t,r}$ associated with $X_{t,r}$ to e.divisive() function from ecp
package which outputs a list of estimated change points [38]. The change point
results are given in Figures 17 and 17, which contain the mean absolute errors
computed over 10 independent trials. In case the expected change point is not
detected, we set the absolute error equal to 50 - the gap between any two
consecutive true change points.

| 1st CP | 2nd CP | 3rd CP | 4th CP
---|---|---|---|---
VPB | 0.0 | 43.6 | 7.4 | 33.5
PI | 0.0 | 10.5 | 3.6 | 19.2
Figure 16: Mean absolute errors for $H_{0}$.

| 1st CP | 2nd CP | 3rd CP | 4th CP
---|---|---|---|---
VPB | 0.0 | 0.5 | 1.7 | 1.0
PI | 0.0 | 17.7 | 1.6 | 6.6
Figure 17: Mean absolute errors for $H_{1}$.

Finally, we examine the optimal values of the turning parameter $\tau$ that
should have been selected based on the ground truth. We find that if there is
little to no noise in the data, then $\tau=0.3$ was optimal for homological
dimensions $H_{1}$ and $H_{2}$ (when relevant) across all our simulations on
clustering and classification. Moreover, we observe that all the classifiers
identify $\tau=0.3$ as optimal based on 10-fold cross-validation as well.

##  5 Conclusion and future work

The induced persistence block map presents a way to transform a PD into a
vector in the Hilbert space $L^{2}(wdA)$ where the area measure $wdA$ is
induced by a continuous nonnegative weight function $w$. We establish the
necessary and sufficient condition for the induced persistence block map to be
continuous with respect to the $p$-Wasserstein distance. As a special case,
when the persistence block map is induced by the constant function 1, we prove
additional results on stability. A VPB is a finite dimensional vector summary
of the image of the induced persistence block map for a given PD. Being
vectors of $\mathbb{R}^{n}$, VPBs can easily be used as inputs within many
machine learning frameworks. VPBs are faster to compute compared to other
vector-based summaries such as persistence images and persistence landscapes.
The parameters involved in the computation of VPBs are scale-free which makes
their selection or tuning easier in practice. Through several simulation
studies covering various learning task such as clustering, classification and
change point detection, we have shown that VPBs produce improved performance
results over persistence images in most cases.

In our simulation studies, VPBs are derived from persistence blocks induced by
the constant function 1. In the future, we plan to consider a wider class of
such functions and investigate the corresponding stability properties.

## References

  * [1] Herbert Edelsbrunner and John Harer.  Computational topology: an introduction.  American Mathematical Soc., 2010.
  * [2] Gunnar Carlsson.  Topology and data.  Bulletin of the American Mathematical Society, 46(2):255–308, 2009\.
  * [3] Afra Zomorodian and Gunnar Carlsson.  Computing persistent homology.  Discrete & Computational Geometry, 33(2):249–274, 2005.
  * [4] Herbert Edelsbrunner, John Harer, et al.  Persistent homology-a survey.  Contemporary mathematics, 453:257–282, 2008.
  * [5] Joseph J Rotman.  An introduction to algebraic topology, volume 119.  Springer Science & Business Media, 2013.
  * [6] Tomasz Kaczynski, Konstantin Mischaikow, and Marian Mrozek.  Computational homology, volume 157.  Springer Science & Business Media, 2006.
  * [7] Yuriy Mileyko, Sayan Mukherjee, and John Harer.  Probability measures on the space of persistence diagrams.  Inverse Problems, 27(12):124007, 2011.
  * [8] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer.  Stability of persistence diagrams.  Discrete & computational geometry, 37(1):103–120, 2007.
  * [9] Frédéric Chazal, Vin De Silva, and Steve Oudot.  Persistence stability for geometric complexes.  Geometriae Dedicata, 173(1):193–214, 2014.
  * [10] David Cohen-Steiner, Herbert Edelsbrunner, John Harer, and Yuriy Mileyko.  Lipschitz functions have l p-stable persistence.  Foundations of computational mathematics, 10(2):127–139, 2010.
  * [11] Peter Bubenik and Tane Vergili.  Topological spaces of persistence modules and their properties.  Journal of Applied and Computational Topology, 2(3):233–269, 2018\.
  * [12] Yen-Chi Chen, Daren Wang, Alessandro Rinaldo, and Larry Wasserman.  Statistical analysis of persistence intensity functions.  arXiv preprint arXiv:1510.02502, 2015.
  * [13] Genki Kusano, Yasuaki Hiraoka, and Kenji Fukumizu.  Persistence weighted gaussian kernel for topological data analysis.  In International Conference on Machine Learning, pages 2004–2013. PMLR, 2016.
  * [14] Chunyuan Li, Maks Ovsjanikov, and Frederic Chazal.  Persistence-based structural recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1995–2002, 2014.
  * [15] Jan Reininghaus, Stefan Huber, Ulrich Bauer, and Roland Kwitt.  A stable multi-scale kernel for topological machine learning.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4741–4748, 2015.
  * [16] Peter Bubenik.  Statistical topological data analysis using persistence landscapes.  The Journal of Machine Learning Research, 16(1):77–102, 2015.
  * [17] Bastian Rieck, Filip Sadlo, and Heike Leitte.  Topological machine learning with persistence indicator functions.  In Topological Methods in Data Analysis and Visualization, pages 87–101. Springer, 2017.
  * [18] Henry Adams, Tegan Emerson, Michael Kirby, Rachel Neville, Chris Peterson, Patrick Shipman, Sofya Chepushtanova, Eric Hanson, Francis Motta, and Lori Ziegelmeier.  Persistence images: A stable vector representation of persistent homology.  The Journal of Machine Learning Research, 18(1):218–252, 2017.
  * [19] Eric Berry, Yen-Chi Chen, Jessi Cisewski-Kehe, and Brittany Terese Fasy.  Functional summaries of persistence diagrams.  Journal of Applied and Computational Topology, 4(2):211–262, 2020\.
  * [20] Nieves Atienza, Rocío González-Díaz, and Manuel Soriano-Trigueros.  On the stability of persistent entropy and new summary functions for topological data analysis.  Pattern Recognition, 107:107509, 2020.
  * [21] Eitan Richardson and Michael Werman.  Efficient classification using the euler characteristic.  Pattern Recognition Letters, 49:99–106, 2014.
  * [22] Yu-Min Chung and Austin Lawson.  Persistence curves: A canonical framework for summarizing persistence diagrams.  arXiv preprint arXiv:1904.07768, 2019.
  * [23] Walter Rudin.  Functional Analysis, Second Edition.  International Series in Pure and Applied Mathematics. McGraw-Hill, 1991\.
  * [24] Robert Ghrist.  Barcodes: the persistent topology of data.  Bulletin of the American Mathematical Society, 45(1):61–75, 2008\.
  * [25] Herbert Edelsbrunner and Ernst P Mücke.  Three-dimensional alpha shapes.  ACM Transactions on Graphics (TOG), 13(1):43–72, 1994.
  * [26] Brittany T Fasy, Jisu Kim, Fabrizio Lecci, Clement Maria, David L Millman, Vincent Rouvreau, and Maintainer Jisu Kim.  Package ‘tda’, 2021.
  * [27] Raoul R Wadhwa, Drew FK Williamson, Andrew Dhawan, and Jacob G Scott.  Tdastats: R pipeline for computing persistent homology in topological data analysis.  Journal of open source software, 3(28):860, 2018.
  * [28] Tullia Padellini, Francesco Palini, Pierpaolo Brutti, Chih-Chung Chang, Chih-Chen Lin, Michael Kerber, Dmitriy Morozov, Arnur Nigmetov, and Maintainer Tullia Padellini.  Statistical learning with kernel for persistence diagrams, 2020.
  * [29] Maria Halkidi, Yannis Batistakis, and Michalis Vazirgiannis.  Clustering validity checking methods.  ACM Sigmod Record, 31(3):19–27, 2002.
  * [30] LKPJ Rdusseeun and P Kaufman.  Clustering by means of medoids.  In Proceedings of the Statistical Data Analysis Based on the L1 Norm Conference, Neuchatel, Switzerland, pages 405–416, 1987.
  * [31] Hae-Sang Park and Chi-Hyuck Jun.  A simple and fast algorithm for k-medoids clustering.  Expert systems with applications, 36(2):3336–3341, 2009.
  * [32] Stuart Lloyd.  Least squares quantization in pcm.  IEEE transactions on information theory, 28(2):129–137, 1982.
  * [33] Martin Maechler, Peter Rousseeuw, Anja Struyf, Mia Hubert, and Kurt Hornik.  cluster: Cluster Analysis Basics and Extensions, 2021.  R package version 2.1.2 — For new features, see the ’Changelog’ file (in the package source).
  * [34] Frederico A Limberger, Richard C Wilson, M Aono, N Audebert, A Boulch, B Bustos, A Giachetti, A Godil, B Le Saux, B Li, et al.  Shrec’17 track: Point-cloud shape retrieval of non-rigid toys.  In 10th Eurographics workshop on 3D Object retrieval, pages 1–11, 2017.
  * [35] Philip Shilane, Patrick Min, Michael Kazhdan, and Thomas Funkhouser.  The princeton shape benchmark.  In Proceedings Shape Modeling Applications, 2004., pages 167–178. IEEE, 2004.
  * [36] Jan-Martin Hertzsch, Rob Sturman, and Stephen Wiggins.  Dna microarrays: design principles for maximizing ergodic, chaotic mixing.  Small, 3(2):202–218, 2007.
  * [37] David S Matteson and Nicholas A James.  A nonparametric approach for multiple change point analysis of multivariate data.  Journal of the American Statistical Association, 109(505):334–345, 2014.
  * [38] Nicholas A James and David S Matteson.  ecp: An r package for nonparametric multiple change point analysis of multivariate data.  arXiv preprint arXiv:1309.3295, 2013.