∎

11institutetext: A. Myers 22institutetext: Michigan State University, College
of Engineering
22email<EMAIL_ADDRESS>33institutetext: F. Khasawneh 44institutetext:
Michigan State University, College of Engineering
44email<EMAIL_ADDRESS>55institutetext: B. Fasy 66institutetext:
Montana State University, School of Computing & Dept. of Mathematical Sciences
66email<EMAIL_ADDRESS>
# ANAPT: Additive Noise Analysis for Persistence Thresholding

Audun D. Myers  Firas A. Khasawneh  Brittany T. Fasy

(Preprint as of )

###### Abstract

We introduce a novel method for Additive Noise Analysis for Persistence
Thresholding (ANAPT) which separates significant features in the sublevel set
persistence diagram of a time series based on a statistics analysis of the
persistence of a noise distribution. Specifically, we consider an additive
noise model and leverage the statistical analysis to provide a noise cutoff or
confidence interval in the persistence diagram for the observed time series.
This analysis is done for several common noise models including Gaussian,
uniform, exponential and Rayleigh distributions. ANAPT is computationally
efficient, does not require any signal pre-filtering, is widely applicable,
and has open-source software available. We demonstrate the functionality ANAPT
with both numerically simulated examples and an experimental data set.
Additionally, we provide an efficient $\Theta(n\log(n))$ algorithm for
calculating the zero-dimensional sublevel set persistence homology.

###### Keywords:

Topological Signal Processing Statistics Topological Data Analysis Sublevel
Set Persistence Sublevel Sets Confidence Intervals Persistent Homology Signal
Processing Extrema Detection Time Series Analysis Cutoff

##  1 Introduction

In this work, we are interested in separating noise from signal in time series
data. We approach this task through studying the statistics of sublevel set
persistence for several common noise models. Specifically, we are using
persistence as a tool from Topological Data Analysis (TDA), which is a method
of data analysis Edelsbrunner2002  used to study datasets from the perspective
of topology or shape. This intersection between TDA and time series analysis
is known as Topological Signal Processing (TSP). TSP includes mathematical
tools for analyzing signals from dynamical systems through TDA, which includes
applications such as dimensionality reduction, detecting bifurcations in the
system behavior, or dynamic state detection (regular or chaotic). TSP methods
help study these applications by leveraging the mathematics of dynamical
systems theory, algebraic topology, information theory, and graph theory.
Because of the novel tools used in TSP, it can be used to reveal information
that may not be currently accessible through standard dynamic systems methods.

In general, there are two common methods for using TSP to study single
variable time series. The first is to study the sublevel sets of a time series
and how they persist as the height of the sublevel set is increased. The
second is to transform the time series into a shape in Euclidean space (e.g.,
a sliding window embedding Takens1981  creates an shape known as the _phase
space_) and study the shape through the Vietoris–Rips filtration. Studying the
shape of the phase space has been successfully used for periodicity
quantification with applications in signal analysis Perea2014 , detecting
quasiperidoicity in video data Tralie2018b , chatter detection in turning
processes Khasawneh2018a ; Yesilli2019 , stability of delayed equations
Khasawneh2017 , and complex network analysis Myers2019 . While the persistent
homology of point clouds has been shown to be a promising tool for time series
analysis, in this work, we take the first approach and study the statistical
analysis of persistence of the sublevel sets of the original time series.

Through persistence, we can capture the significance of sublevel sets in the
time series, which can be extended to higher dimensional data-sets (a more
detailed introduction follows in Section 2). Sublevel set persistence has been
successfully used for a wide range of time series analysis applications
including local extrema detection Gholizadeh2018 , true step detection in
signal processing Khasawneh2018b , fourier spectrum analysis Myers2020 ,
arrythmia detection Dindin2020 , and cancer histology Lawson2019 .

One of the most attractive features of sublevel set persistence as a time
series analysis tool is its stability to noise perturbations (see, e.g.,
Cohen-Steiner2006 ). However, even with the stability of persistence diagrams,
there can be many points in the diagram which we wish to distinguish as
associated to either noise or true signal. One method of doing this is by
estimating a confidence interval for each point in the persistence diagram,
which has recently been of interest in the field of TDA Bubenik2015 ;
Chazal2013 ; Fasy2014 . While methods based on bootstrap resamping
chazal2017robust ; Chazal2013 ; Fasy2014  and entropy Atienza2017  have been
developed for separating noise from signal in the persistence diagrams, they
were not strictly designed for sublevel set persistence. This can lead to
problems in both implementation and computational demand when applied to
sublevel set persistence. These issues are discussed in Section 2 and then
demonstrated in Section 7.1.

This manuscript is organized as follows. In Section 2, we first provide an
introduction of persistence for time series data and the effects of additive
noise on the persistence diagram as well as commonly used methods for
seperating noise from signal. In Section 3, we introduce a novel analysis of
the statistics of the lifetimes of the sublevel set persistence diagram. In
Section 4, we realize the statistical analysis for cutoff expressions for
Gaussian, uniform, Rayleigh, and exponential additive noise distribution. In
Section 5 and 6, we provide a method to approximate the distribution parameter
of the additive noise as well as a signal compensation term, which is required
in the cutoff equation. Lastly, in Section 7, we apply our ANAPT method to
both numeric and experimental time series to validate its functionality and
make comparisons to existing methods.

##  2 Sublevel Set Persistence

We now provide an introduction to persistence, as it applies to computing
persistence of the sublevel sets of a time series. Let us begin with the
single variable function $f:\mathbb{R}\rightarrow\mathbb{R}$. Given
$r\in\mathbb{R}$, we define the _sublevel set below $r$_ as
$f^{-1}(-\infty,r]$. As the filtration parameter $r$ increases, the sublevel
sets may grow but remain the same (up to homology) until a local extrema
(i.e., a local minimum or maximum) is reached.111Here, we assume that the
function $f$ satisfies some mild “niceness” conditions, such as being q-tame
chazal2016structure . If the extrema is a local minimum, then a new connected
component or set is “born” at height $r_{B}$; we label that set with the value
$r_{B}$. On the other hand, if the extrema is a local maximum, two previously
existing sets are combined. If the two sets were labeled $r_{B}$ and
$r_{B}^{\prime}$, with $r_{B}\leq r_{B}^{\prime}$ and the maximum attained at
$r_{D}$, then, by the Elder Rule (Edelsbrunner2010, , p. 150), we say that the
component born at $r_{B}^{\prime}$ dies going into $r_{D}$, and the resulting
set assumes the label $r_{B}$. The pair $(r_{B}^{\prime},r_{D})$ is called a
persistence pair. As $r$ ranges from $-\infty$ to $\infty$, the _persistence
diagram_ is the collection of all $n$ such pairs,
$\mathcal{D}(f)=\\{(b_{i},d_{i})\\}_{i=1}^{n}$. Any unpaired births are called
_essential classes_ and are paired with a death coordinate of $\infty$; thus,
$\mathcal{D}(f)$ is embedded in the _extended plane_
$\overline{\mathbb{R}}^{2}$. The _lifetime_ or _persistence_ of a point
$(b_{i},d_{i})\in\mathcal{D}(f)$ is defined as $\ell_{i}=d_{i}-b_{i}$. In this
paper, our functions are only sampled on a finite domain, with the first
sample at time $t_{a}$ and the last sample at time $t_{b}$. We obtain a
continuous function over $[t_{a},t_{b}]$ by using a piecewise linear
interpolation between consecutive samples, and extending the function to
$\pm\infty$ by extending the first (resp., last) edges to rays. Doing so
allows us to define a persistence diagram that does not have critical points
on the boundary of our time series. As such, we study the persistence points
where both coordinates are finite, and omit persistence points that contain an
unbounded coordinate.

To give a concrete example of a persistence diagram, we demonstrate a simple
example for the function shown in Fig. 1. This function has thirteen sample
points, two local minima, and two local maxima.

Figure 1: Persistence diagram summarizing the sublevel sets of a function
$f(t)$ over finite domain $t\in[t_{a},t_{b}]$. This function has two local
minima (blue squares) and two local maxima (red circles.

The lowest finite critical value of the function occurs at height $v_{0}$. For
all $r<v_{0}$, $f^{-1}(-\infty,r]$ is the ray $[f^{-1}(r),\infty)$, where we
observe that $f^{-1}(r)$ is a well-defined and greater than $t_{b}$. This
connected component is labeled with $-\infty$, since it is “born” at
$-\infty$. Then, at height $r=v_{0}$, a second connected component is born.
The next topological change occurs at height $r=v_{1}$, where a third
connected component is born. The next extrema is reached when $r=p_{0}$. At
this extrema, the sublevel set that was born at $r=v_{1}$ dies, while the
sublevel set born at $r=v_{0}$ persists based on the Elder Rule. This pair
$(v_{1},p_{0})$ is recorded in the persistence diagram. From here, the next
change happens at $r=p_{1}$, where the second sublevel set dies and is
recorded in the persistence diagram as $(v_{0},p_{1})$. Then, no further
topological changes occur, but this sublevel set continues to grow as $r$
grows. This essential class is recorded in the persistence diagram as
$(-\infty,\infty)$ and is not studied in our analysis. As shown in the
persistence diagram, the point $(v_{1},p_{0})$ is close to the diagonal (the
line $y=x$), which signifies that the sublevel set only persisted for a short
range of heights ($r$); on the other hand, the point $(v_{0},p_{1})$ is far
from the diagonal, suggesting it was from a significant sublevel set. The
algorithm used to calculate the this persistence diagram is detailed in
Appendix A.2.

The idea of persistence can be extended to higher dimensions allowing for the
analysis of the shape of high-dimensional data sets. However, for our work, we
only need to analyze the zero-dimensional features (i.e., connected
components) of a one-dimensional function. A more thorough background on TDA,
and persistent homology specifically, can be found in Edelsbrunner2008 ;
Munch2017 ; Perea2018 . Other common ways for studying time series with a
similar perspective is through merge trees or dendograms Carlsson2012 ;
Chowdhury2017 ; Lee2012 .

_Sublevel Set Persistence with Additive Noise_.

We now investigate the stability of sublevel set persistence diagrams to
additive noise for single variable functions. To illustrate the stability, we
first take an example time series with additive noise as $x(t)+\epsilon$,
where $x(t)$ is sampled at a uniform rate $f_{s}$ and $\epsilon$ is additive
noise from the noise model $\mathcal{N}$. An example of a persistence diagram
from the time series with additive noise $\mathcal{D}(x+\epsilon)$ is shown in
Fig. 2, along with the diagram without the additive noise $\mathcal{D}(x)$.
This example also demonstrates how a cutoff $C_{\alpha}$ can be used to
separate the significant points in the persistence diagram and those
associated to the additive noise.

Figure 2: Sublevel set persistence applied to $x(t)$ of a single variable
function or time series with and without additive noise $\epsilon$ from
$\mathcal{N}$, shown in red and blue, respectively. This demonstrates the
stability of persistent homology with the time series (left) with and without
additive noise and the small effect on the resulting persistence diagrams
(right). In addition, the light red region separates the significant features
from those associated to additive noise.

This example demonstrates that the addition of (small) noise does not have a
large effect on the position of significant sublevel sets in the persistence
diagram with the distances between significant points ($d_{1}$ and $d_{2}$)
all being relatively small. This is no surprise due to the stability theorem
of the bottleneck distance for persistence diagrams Cohen-Steiner2006 ;
chazal2016structure , where the bottleneck distance is defined as the minimum
distance to match two persistence diagrams. For example, if we assume
$d_{1}>d_{2}>d_{3}>d_{4}$, then the bottleneck distance would be $d_{1}$.
However, additive noise does introduce several points in the persistence
diagram located near the diagonal with relatively small lifetimes. These
noise-artifact persistence pairs are formed from the peak-valley pairs in the
additive noise. This work focuses on a statistical analysis of these lifetimes
to develop a method for separating the significant persistence diagram points
from those of additive noise, shown in as light red region in example
persistence diagram of Fig. 2, through a cutoff $C_{\alpha}$ with
$\alpha\in[0,1]$ as the given confidence level.

_Current Methods for Separating Noise from Signal in the Persistence Diagram_.
The methods for developing confidence sets and separating noise from signal
for persistence diagrams using resampling techniques Chazal2013 ; Fasy2014
were not strictly developed for persistence of sublevel sets of real-valued
time series. This means that using the bottleneck bootstrap as presented in
chazal2017robust  gives conservative confidence sets. One solution to this
problem could be to tailor these methods for sublevel set persistence by
applying resampling techniques for time series (e.g., stationary bootstrap
Politis1994  and the sieve bootstrap Buehlmann1997 ). However, these methods
have general requirements of weakly dependent stationary observations or
exceptionally high sampling rates such that $f_{s}>>f^{\prime}$, where
$f^{\prime}$ the highest dominant frequency of the time series. An alternative
approach to applying a bootstrap resampling is to use a function fitting
procedure that is suitable for strongly dependent time series. Using the
function fit of the signal, the residuals can be resampled to the time series
to apply the bottleneck bootstrap technique. Unfortunately, this method is
plagued with additionally computational steps including the function fitting
initial step, repeated sublevel set persistence calculations, and the
bottleneck distance calculation. The possibility of signal distortion when
fitting to the curve can also cause for inaccurate cutoff estimations. Another
method for separating noise from significant features in a persistence diagram
is persistent entropy Atienza2017 . However, this method may not properly
distinguish between noise and significant features in the persistence diagram
if the number of significant points in the persistence diagram is relatively
large compared to the amount of noise.

Because of the significant drawbacks of the current methods for developing
confidence sets and associated cutoffs in the persistence diagram for sublevel
set persistence, we have developed a new method that directly analyzes the
statistics of additive noise in the sublevel set persistence diagram. We have
also made the Python code for implementing the $\Theta(n\log(n))$ algorithm
(see Section A.2) discussed throughout the manuscript publicly available
through the Python package teaspoon222https://lizliz.github.io/teaspoon/.

##  3 Statistics of Additive Noise in the Persistence Diagram

Before studying a time series with additive noise,
$x+\epsilon\colon\mathbb{R}\to\mathbb{R}$, we analyze the statistics of
sublevel set persistence diagrams of the noise alone. Our goal is to leverage
this analysis in order to generate a cutoff in the persistence diagram to
separate out these noise-artifact points in the persistence diagram for
$\mathcal{D}(x+\epsilon)$ from the points that capture true features of $x$.

###  3.1 Sublevel Set Persistence Diagram Statistics Background

We start with the noise, which can be thought of as a (sampled) function
$\epsilon\colon\mathbb{R}\to\mathbb{R}$, where, for each $t\in\mathbb{R}$, the
value $\epsilon(t)$ is a random independently and identically distributed
(iid) variable sampled from a predefined noise distribution $\mathcal{N}$. In
our noise model, there is no covariance structure between these
$\mathbb{R}$-indexed family of random variables. Let $\mathcal{N}(\mathbb{R})$
denote the induced probability distribution over functions from $\mathbb{R}$
to $\mathbb{R}$, so $\epsilon\sim\mathcal{N}(\mathbb{R})$.

The first step in developing a cutoff based on the persistence diagram
statistics of additive noise $\mathcal{D}(\epsilon)$ is to determine a
relationship between the descriptive additive noise distribution parameters
and the distribution of the lifetimes. To do this, we develop an expression
for the expected lifetime of points in $\mathcal{D}(\epsilon)$. Let
$f\colon\mathbb{R}\to\mathbb{R}$ and $F\colon\mathbb{R}\to\mathbb{R}$ be the
probability density function and cumulative density function of $\mathcal{N}$,
respectively. Let $f_{B}\colon\mathbb{R}\to\mathbb{R}$ and
$f_{D}\colon\mathbb{R}\to\mathbb{R}$ be the probability density functions for
the local minima (corresponding to births) and maxima (corresponding to
deaths) of the sublevel sets from a random noise function
$\varepsilon\sim\mathcal{N}(\mathbb{R})$. Let $F_{B}$ and $F_{D}$ be the
corresponding cumulative density functions. By the commutative property of
addition and the definition of a lifetime being the difference between the
death and birth times, the expected or mean lifetime $\mu_{L}$ is the
difference between the expected birth times $\mu_{B}:=\mathbb{E}(B)$ and death
times $\mu_{D}:=\mathbb{E}(D)$, where $B$ and $D$ are the sets of birth and
death values, as

$\mu_{L}:=\mu_{D}-\mu_{B}=\int_{-\infty}^{\infty}z\left[f_{D}(z)-f_{B}(z)\right]dz.$
(1)

A formal proof of this relationship is provided in Theorem A.1 of the
appendix. From Equation (1), we move forward knowing that $\mu_{L}$ can be
defined using only expressions for $f_{B}$ and $f_{D}$. In other words, we
only need the distribution of birth and death times, not of the lifetimes,
which would require knowing how the births and deaths are paired.

Let $z\in\mathbb{R}$ and $\epsilon_{L},\epsilon_{R}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}$. Because $\epsilon_{L}$ and $\epsilon_{R}$ are
independent, by the multiplication rule for probability, we can state that the
probability that $\max\\{\epsilon_{L},\epsilon_{R}\\}<z$ is

$p(\epsilon_{L}<z)p(\epsilon_{R}<z).$ (2)

Recalling that $F$ is the cumulative distribution function of $\mathcal{N}$,
we can equate these probabilities as

$p(\epsilon_{L}<z)=p(\epsilon_{R}<z)=F(z).$

Combining with Equation (2), we define the probability distribution function
of maximas $f_{\max}\colon\mathbb{R}\to\mathbb{R}$ by

$f_{\max}(z)=F^{2}(z).$ (3)

Similarly, we can show that the local minima probability distribution function
$f_{\min}\colon\mathbb{R}\to\mathbb{R}$ is defined by

$f_{\min}(z)=[1-F(z)]^{2},$ (4)

since $p(\epsilon_{L}>z)=p(\epsilon_{R}>z)=1-F(z)$.

Using these distributions, we calculate the local maxima probability density
function $f_{D}$, which represents the distribution of death times. We proceed
in the discrete setting and consider $\epsilon\sim\mathcal{N}(\mathbb{R})$
represented by a finite sequence of values; that is, we have
$\epsilon=(\epsilon_{1},\epsilon_{2},\ldots,\epsilon_{n})$ where
$\epsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}$. Fixing
$z=\epsilon_{i}$ $\epsilon_{L}=\epsilon_{i-1}$, and
$\epsilon_{R}=\epsilon_{i+1}$, we know the probability that $\epsilon_{i}$ is
a maximum (or a minimum) from above. To calculate $f_{D}$ (resp., $f_{B}$), we
take the convolution of $f$ and $f_{\max}$ (resp., of $f$ and $f_{\min}$) and
normalize (to ensure the density integrates to unity):

$\begin{split}{f}_{D}(z)&:=\frac{f(z)F^{2}(z)}{N_{D}},\\\
{f}_{B}(z)&:=\frac{f(z)[1-F(z)]^{2}}{N_{B}},\end{split}$ (5)

where $N_{B}=\int_{-\infty}^{\infty}f(z)[1-F(z)]^{2}dz$ and
$N_{D}=\int_{-\infty}^{\infty}f(z)F^{2}(x)dz$. We calculate $N_{B}$ and
$N_{D}$ from Equation (5) by substituting $f(z)=F^{\prime}(z)$, which reduces
the $N_{D}$ equation to

$\begin{split}N_{D}&=\int_{-\infty}^{\infty}F^{\prime}(z)F^{2}(z)dz=\int_{-\infty}^{\infty}\frac{1}{3}{(F^{3}(z)})^{\prime}dz\\\
&=\frac{1}{3}F^{3}(z)|_{-\infty}^{\infty}=\frac{1}{3},\end{split}$ (6)

since $F(\infty)=1$ and $F(-\infty)=0$. Similarly,

$\begin{split}N_{B}&=\int_{-\infty}^{\infty}f(z)[1-F(z)]^{2}dz=\int_{-\infty}^{\infty}f(z)[1-2F(z)+F^{2}(z)]dz=N_{D}+\int_{-\infty}^{\infty}f(z)[1-2F(z)]dx\\\
&=N_{D}+\int_{-\infty}^{\infty}F^{\prime}(z)dz-\int_{-\infty}^{\infty}(F^{2}(z))^{\prime}dz=N_{D}+\left[F(z)-F^{2}(z)\right]|_{-\infty}^{\infty}=N_{D}=\frac{1}{3}.\end{split}$
(7)

We can now reduce Equation (5) to

$\begin{split}{f}_{B}(z)&=3{f(z)[1-F(z)]^{2}},\\\
{f}_{D}(z)&=3{f(z)F^{2}(z)},\end{split}$ (8)

For demonstrative purposes we now assume $f(z)$ is of a Gaussian distribution
to validate our expressions in Equation (5). Specifically, the Gaussian
probability density function is defined as

$f(z)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(z-\mu)^{2}}{2\sigma^{2}}},$
(9)

with a cumulative distribution

$F(z)=\frac{1}{2}\left[1+{\rm
erf}\left(\frac{z-\mu}{\sigma\sqrt{2}}\right)\right].$ (10)

To validate the resulting expressions for ${f}_{B}(z)$ and ${f}_{D}(z)$ in
Equation (8), a numerical simulation of a normal distribution
$\mathcal{N}_{n}(\mu=0,\sigma^{2}=1)$ of length $n=10E5$ is used (see Fig. 3).
This analysis shows a very similar result between the histograms $h(*)$ and
distributions.

Figure 3: Histograms $h(*)$ of the zero mean normal distriubtion
$\mathcal{N}(0,\sigma^{2}=1)$ and the resulting birth times $B$ and death
times $D$, which are compared to the density distributions from Equation (5).

Now that we have shown that our expressions for the probability distribution
of the minima and maxima are correct, we proceed to correlate the mean
lifetime $\mu_{L}$ to the additive noise distribution parameters. From our
results in Equation (8) applied to Equation (1) we calculate the mean lifetime
as

$\mu_{L}=3\int_{-\infty}^{\infty}zf(z)\left[F^{2}(z)-(1-F(z))^{2}\right]dz=3\int_{-\infty}^{\infty}z\left[(F^{2}(z))^{\prime}-F^{\prime}(z)\right]dz,$
(11)

which is simplified using integration by parts as

$\boxed{\mu_{L}=3\int_{-\infty}^{\infty}F(z)\left[1-F(z)\right]dz.}$ (12)

From the numerical simulation in Fig. 3 we also found the $\bar{L}\approx
1.6921$, where $\bar{L}$ is the sample mean of the sublevel set persistence
lifetimes from $\mathcal{N}(0,\sigma^{2}=1)$. This is very similar compared to
the numerically estimated theoretical $\mu_{L}\approx 1.6925$ from Equation
(12). These results suggest that Equation (12) is correct.

###  3.2 Cutoff Background

To determine a suitable cutoff using our ANAPT method, we again start by
assuming we have $n$ random samples from our noise distribution:
${\epsilon}=\\{\epsilon_{1},\epsilon_{2},\ldots,\epsilon_{n}\\}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}$ with a cumulative probability function $F$. The
probability that the minimum of ${\bf\epsilon}$ is less than the value $a$ is
equivalent to

$P(\min({\bf\epsilon})<a)=1-P(\epsilon_{1}>a,\epsilon_{2}>a,\ldots,\epsilon_{n}>a),$
(13)

where $P(\epsilon_{i}>a)=1-F(a)$. If we extend this relationship to all $n$
realizations, we can express this probability as

$P(\min({\bf\epsilon})<a)=1-(1-F(a))^{n}.$ (14)

Similarly, for $b>a$, an expression for the probability that an element of
${\epsilon}$ is greater than $b$ is

$P\left(\max({\epsilon})>b\right)=1-(F(b))^{n}.$ (15)

If we now take both of these probabilities, we can extend them to the maximum
finite lifetime as $\max(L)\lessapprox\max({\epsilon})-\min({\epsilon})$. we
can use this to generate a probability of a lifetime being greater than $b-a$
as

$\alpha=P(\max(L)>b-a)\gtrapprox
P\big{(}\max({\epsilon})>b,\min({\epsilon})<a\big{)}=(1-[F(b)]^{n})(1-[1-F(a)]^{n}),$
(16)

where $\alpha$ is the confidence of this event occurring. If the $f(z)$
associated to $F(z)$ of Equation (16) is symmetric about some mean $\mu$ such
that $c=b-\mu=\mu-a$, we can reduce Equation (16) to

$\alpha={(1-[F(c)]^{n})}^{2}$ (17)

since $F(b)=1-F(a)$ for the symmetric case. Equation (17) can be then solved
for $c$ as

$c=F^{-1}\left[\left(1-\sqrt{\alpha}\right)^{1/n}\right].$ (18)

Additionally, we know that a cutoff should be set such that
$C_{\alpha}=b-a=2c$ resulting in a cutoff equation as

$\boxed{C_{\alpha}=2F^{-1}\left[{\left(1-\sqrt{\alpha}\right)}^{1/n}\right].}$
(19)

If there is no symmetry in the distribution then we need a new cutoff
equation. To do this, we return to our probability equation as

$\alpha=P(\max(L)>b-a)\gtrapprox
P(\min({\epsilon})<a,\max({\epsilon})>b)=(1-{[1-F(a)]}^{n})(1-[F(b)]^{n}),$
(20)

However, unlike Equation (19), we can not solve Equation (16) for a parameter
$c$ due to their being no symmetry between $a$ and $b$ about a mean $\mu$
which means we must simplify Equation (20). To do this, we assume that
$P(\min({\epsilon})<a)=P(\max({\epsilon})>b)$ or
$1-{[1-F(a)]}^{n}=1-[F(b)]^{n}=\sqrt{\alpha}$. We then solve for $a$ and $b$
separately as

$a=F^{-1}\left[1-{(1-\sqrt{\alpha})}^{1/n}\right]$ (21)

and

$b=F^{-1}\left[{(1-\sqrt{\alpha})}^{1/n}\right].$ (22)

With $C_{\alpha}=b-a$ and the values of $a$ and $b$ from Equation (21) and
Equation (22), respectively, we can solve for our general cutoff expression as

$\boxed{C_{\alpha}=F^{-1}\left[{(1-\sqrt{\alpha})}^{1/n}\right]-F^{-1}\left[1-{(1-\sqrt{\alpha})}^{1/n}\right].}$
(23)

For our application we want to have a high confidence level that no
persistence pairs associated to noise occur with a lifetime greater than the
cutoff. We suggest using a confidence level of $\alpha=0.001$ or $0.1\%$ for a
low probability of this occurring.

As a summary, the ANAPT cutoff equations in Equation (19) and Equation (23)
are only dependent on the desired confidence $\alpha$, the signal length $n$
and the cumulative probability distribution $F(z)$. In Section 4 we
demonstrate how to apply Equation (19) and Equation (23) for the Gaussian,
uniform, Rayleigh, and exponential distribution.

##  4 Cutoff for Noise Models

For applying noise models to the confidence levels in Equations (16) and (17)
for our ANAPT method, we need to be either given the additive noise
parameters, or estimate them from the lifetimes. However, before this can be
done, we need to understand which parameters are needed given the additive
noise distribution $f(z)$. We do this analysis for Gaussian (normal), Uniform,
Rayleigh, and exponential distributions as shown in Fig. 4.

Figure 4: Additive noise probability distributions $f(z)$ for the four models
realized in this work: uniform, Gaussian, Rayleigh, and exponential.

###  4.1 Cutoff for Gaussian Noise

We start our analysis with the commonly used Gaussian distribution model. The
Gaussian (normal) probability distribution function is defined as

$f(z)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(z-\mu)^{2}}{2\sigma^{2}}},$
(24)

with a cumulative distribution function

$F(z)=\frac{1}{2}\left[1+{\rm
erf}\left(\frac{z-\mu}{\sigma\sqrt{2}}\right)\right].$ (25)

We start by solving for the inverse of Equation (25) as

$F^{-1}(u)=\sqrt{2}\sigma{\rm erf}^{-1}\left(2u-1\right)+\mu.$ (26)

Since the mean shift $\mu$ has no effect on the sublevel set lifetimes we can
ignore it and apply Equation (26) with $\mu=0$ to solve for the cutoff from
Equation (19) as

$C_{\alpha}=2^{3/2}\sigma\>{\rm
erf}^{-1}\left[2(1-\sqrt{\alpha})^{1/n}-1\right].$ (27)

With a full development of the statistics of sublevel set persistence for
Gaussian (normal) additive noise we are able to determine a suitable cutoff
for iid noise with only the distribution parameter $\sigma$ needing estimated.

###  4.2 Cutoff for Uniform Noise

Let $a<b\in\mathbb{R}$. The uniform distribution over the interval $[a,b]$ has
a probability density function defined as

$f(z)=\left\\{\begin{array}[]{ll}\frac{1}{b-a}&z\in[a,b]\\\
0&\text{otherwise}\end{array}\right.$ (28)

with a cumulative distribution function

$F(x)=\left\\{\begin{array}[]{ll}0&z<a\\\ \frac{z-a}{b-a}&z\in[a,b]\\\
1&z>b.\end{array}\right.$ (29)

By assuming a symmetric distribution about zero (this assumptions does not
influence the resulting cutoff due to the properties sublevel set persistence
lifetime) such that $a=-b$ and $\Delta=b-a$. This changes $F(x)$ to

$F(z)=\left\\{\begin{array}[]{ll}0&z<-\frac{\Delta}{2}\\\
\frac{2z+\Delta}{2\Delta}&z\in[-\frac{\Delta}{2},\frac{\Delta}{2}]\\\
1&z>\frac{\Delta}{2}\end{array}\right.$ (30)

If we now apply Equation (19) to the inverse of the cumulative probability
distribution in Equation (30), we can calculate $C_{\alpha}$ as

$C_{\alpha}={\Delta}\left[{2\left(1-\sqrt{\alpha}\right)}^{1/n}-1\right].$
(31)

Equation (31) only requires an estimate of the distribution parameter $\Delta$
as both $\alpha$ and $n$ are chosen as desired and the length of the time
series, respectively.

###  4.3 Cutoff for Rayleigh Noise

The Rayleigh distribution has a probability density function over the domain
$x\in[0,\infty)$ and is defined as

$f(z)=\frac{z}{\sigma^{2}}e^{-\frac{z^{2}}{2\sigma^{2}}},$ (32)

with a cumulative distribution function

$F(z)=1-e^{-\frac{z^{2}}{2\sigma^{2}}}.$ (33)

Since this distribution is asymmetric we use Equation (23) to calculate
$C_{\alpha}$ as

$C_{\alpha}=\sigma\left(\sqrt{-2\ln\left([1-\sqrt{\alpha}]^{1/n}\right)}-\sqrt{-2\ln\left(1-[1-\sqrt{\alpha}]^{1/n}\right)}\right),$
(34)

where $\sigma$ is the only parameter that needs to be provided to calculate
the cutoff.

###  4.4 Cutoff for Exponential Noise

The exponential distribution has a probability density function over the
domain $z\in[0,\infty)$ and is defined as

$f(z)=\lambda e^{-\lambda z},$ (35)

with a cumulative distribution function

$F(z)=1-e^{-\lambda z},$ (36)

where $\lambda$ is the distribution parameter with $\lambda>0$. This this
distribution is also asymmetric, so we use Equation (23) to calculate
$C_{\alpha}$ as

$C_{\alpha}=-\frac{1}{\lambda}\ln\left([1-\sqrt{\alpha}]^{1/n}-[1-\sqrt{\alpha}]^{2/n}\right),$
(37)

where $\lambda$ is the only parameter that needs to be provided to calculate
the cutoff.

##  5 Cutoff and Distribution Parameter Estimation Method

If the distribution parameter is know ($\sigma$ for Gaussian distributions,
$\Delta$ for uniform distributions, $\sigma$ for Rayleigh distributions, and
$\lambda$ for exponential distributions), then the cutoff $C_{\alpha}$ using
our ANAPT method can be calculated simply with the use of the correct cutoff
equation in Section 4 and the subsequent analysis may be skipped. However, in
most real-world time series it is uncommon to know what this parameter is and
thus it needs to be estimated. While there are some methods for estimating the
additive noise parameters Czesla2018 ; Urbanowicz2003 ; Hu2004 , we introduce
a new method utilizing the relationship between the sublevel set lifetimes
from both the signal and noise and the additive noise distribution parameters.

To generate a theoretical relationship between the mean lifetime $\mu_{L}$ and
the distribution parameters, we recall Equation (12):

${\mu_{L}=3\int_{-\infty}^{\infty}F(z)\left[1-F(z)\right]dz.}$

In the subsequent subsections, we show how this relationship is used for each
of the four noise models analyzed in this work. However, when the signal is
not pure noise, which would be the case for any informative time series, the
mean lifetime is heavily influenced from the lifetimes associated to
significant features. To address this issue, we instead calculate the median
of the lifetimes as it is robust up to $50\%$ outliers and apply a signal
compensation. In this sense outliers are referring to the persistence pairs
associated to signal that are greater than the cutoff. This brings up an
assumption for this distribution parameter estimation method to function
correctly: the number of persistence diagram features associated with noise
$N_{n}$ must be equal to or greater than the number of features from the
signal $N_{s}$. Additionally, when $N_{n}$ approaches $N_{s}$ the cutoff
becomes more conservative due to the robustness limitation of the median. To
minimize this effect, in Section 6, we develop a numeric compensation
multiplier which uses the persistence pairs associated to both additive noise
and signal. In general, the condition for $N_{s}<N_{n}$ is met if the time
series is sampled at a rate sufficiently higher than the Nyquist sampling
criteria $f_{\rm Nyquist}$ and, of course, the time series has some additive
noise. If these conditions are not met, we suggest the use of an alternative
method to estimate the distribution parameter of the additive noise and apply
its associated cutoff equation in Section 4.

For a symmetric distribution of the lifetimes, the median would be an accurate
estimate of the mean. However, for most additive noise distributions (e.g.
Gaussian), the distribution of the resulting sublevel set persistence
lifetimes is not symmetric. Therefore, we resort to approximating the
relationship between the mean and median numerically. While there are methods
to estimate the mean using the median and Inter-Quartile Range (IQR) as
described in Wan2014 . This method is only robust for up to $25\%$ outliers
due to the $Q_{3}$ upper quartile. Therefore, we use the numerically
approximated ratios of $\rho=\bar{L}/\tilde{L}$ as provided in Table 1 for
each of the four distributions investigated, where $\bar{L}$ is the sample
mean lifetime and $\tilde{L}$ is the sample median lifetime. For each of these
numeric estimates a time series of length $10^{5}$ was used. This numeric
experiment was repeated ten times to provide a mean $\rho$ with uncertainty.
This ratio can be used to estimate the mean lifetime as
$\bar{L}\approx\rho\tilde{L}$.

Table 1: Ratios $\rho=\bar{L}/\tilde{L}$ for estimating sample mean from the sample median with uncertainty as three standard deviations Distribution | Gussian | Uniform | Rayleigh | Exponential
---|---|---|---|---
$\rho=\bar{L}/\tilde{L}$ | $1.154\pm 0.012$ | $1.000\pm 0.010$ | $1.136\pm 0.013$ | $1.265\pm 0.016$

###  5.1 Relating the distribution Statistic to the Median Lifetime

We now apply Equation (12) and $\rho$ from Table 1 to find relationships
between the median lifetime $M_{L}$ and the distribution parameter used in
each distribution’s cutoff equation.

_Normal Distribution:_. For estimating $\sigma$ of the Gaussian distribution,
we use Equation (12) and the Gaussian cumulative distribution to estimate
$\mu_{L}$ as a function of $\sigma$. Specifically, by numerically
approximating the integral in Equation (12) using $z\in[-10,10]$ with $10^{6}$
uniformly spaced samples, we found the relationship

$\sigma\approx\frac{\mu_{L}}{1.6925}.$ (38)

We then use $\rho$ to have Eq. Equation (38) as a function of the median
lifetime $M_{L}$ as

$\sigma\approx\frac{\rho M_{L}}{1.692}\approx 0.680M_{L},$ (39)

where $M_{L}$ is the median lifetime. Applying this result to Equation (27)
allows for a cutoff to be calculate as

$C_{\alpha}\approx 1.923\tilde{L}\>{\rm
erf}^{-1}\left[2(1-\sqrt{\alpha})^{1/n}-1\right],$ (40)

where $\tilde{L}$ is the sample median lifetime.

_Uniform Distribution_. Next, we apply Equation (12) to the uniform cumulative
distribution to estimate $\mu_{L}$ as a function of $\Delta$. Substituting
Equation (30) into Equation (12) results in

$\mu_{L}=3\int_{-\Delta/2}^{\Delta/2}\frac{2z+\Delta}{2\Delta}\left[1-\frac{2z+\Delta}{2\Delta}\right]dz.$

Expanding and solving this integral results in the relationship

$\mu_{L}=\frac{\Delta}{2}\implies\Delta=2\mu_{L}=2M_{L}.$ (41)

Applying this result to Equation (31) allows for a cutoff to be calculate as

$C_{\alpha}=2\tilde{L}\left[{2\left(1-\sqrt{\alpha}\right)}^{1/n}-1\right].$
(42)

_Rayleigh Distribution_. For estimating $\sigma$ of the Rayleigh distribution,
we again use Equation (12) with the cumulative Rayleigh distribution to
numerically estimate the relationship between $\mu_{L}$ and $\sigma$ as

$\sigma\approx\frac{\mu_{L}}{1.102}\approx\frac{\rho M_{L}}{1.102}\approx
1.025M_{L},$ (43)

where the integral in Equation (12) was numerically approximated using
$z\in[0,20]$ with $10^{6}$ uniformly spaced samples. Applying this result to
Equation (34) allows for a cutoff to be calculate as

$C_{\alpha}\approx
1.025\tilde{L}\left(\sqrt{-2\ln\left([1-\sqrt{\alpha}]^{1/n}\right)}-\sqrt{-2\ln\left(1-[1-\sqrt{\alpha}]^{1/n}\right)}\right).$
(44)

_Exponential Distribution_. Next, we apply Equation (12) to the exponential
cumulative distribution function to estimate $\mu_{L}$ as a function of
$\lambda$. Substituting Equation (36) into Equation (12) results in

$\mu_{L}=3\int_{0}^{\infty}\left(1-e^{-\lambda z}\right)e^{-\lambda z}dz,$
(45)

which was solved using a $u$-substitution as

$\mu_{L}=\frac{3}{2\lambda}\rightarrow\lambda=\frac{3}{2\mu_{L}}.$ (46)

By then using the appropriate $\rho$ from Table 1 to use $M_{L}$ instead of
$\mu_{L}$, we approximate $\lambda$ from the median lifetime:

$\lambda\approx\frac{1.875}{M_{L}}.$ (47)

Applying this result to Equation (37) allows for a cutoff to be calculate as

$C_{\alpha}\approx-0.533{\tilde{L}}\ln\left([1-\sqrt{\alpha}]^{1/n}-[1-\sqrt{\alpha}]^{2/n}\right).$
(48)

##  6 Signal Compensation for the Cutoff and Distribution Parameter

In this section, we discuss the effects of signal on the cutoff estimation
methods described. In Section 5.1 we assumed that the time series was of the
form
$\epsilon=\\{\epsilon_{1},\epsilon_{2},\ldots,\epsilon_{n}\\}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}$. However, in practice we typically have some
underlying informative signal $s\colon\mathbb{R}\to\mathbb{R}$ and have a time
series of the form $x(t)=s(t)+\epsilon$ with a finite domain as
$t\in[t_{a},t_{b}]$. The resulting sublevel sets from $s(t)+\epsilon$ are
assumed to have some lifetimes from $s(t)$ with the slope of the signal having
an effect on the lifetimes associated with $\mathcal{N}$. Because of these
effects, we attempt to compensate the cutoff calculation and distribution
parameter estimations for these effects for a general signal. Since a general
signal is, in practice, rather subjective, we move away from a theoretical
analysis of the signal and rather analyze the effects of the signal
experimentally. We have partially addressed this issue of signal compensation
by implementing the median lifetime $M_{L}$ instead of the mean lifetime
$\mu_{L}$ with the median being an outlier robust statistic for up to 50%
outliers. Even with the use of the median, we need to further develop a signal
compensation procedure to improve the accuracy of the suggested cutoff.

Figure 5: Example time series showing sample $\delta_{i}$.

To fully understand the effects of signal on estimating the cutoff, we do a
numeric study to develop a method for adjusting the median lifetime such that
$RM_{L}(s(t)+\epsilon)\approx M_{L}(\mathcal{N})$, where $R$ is the
compensation term. This analysis requires a new variable which we term
$\delta$ as the median step size with $\delta_{i}=x(t_{i+1})-x(t_{i})$ as
shown in Fig. 5, where $x(t)$ is a discretely and uniformly sampled signal
with a constant sampling rate $f_{s}$. We now experimentally approximate the
effects of signal on the median lifetime by using three “generic" signals
suggested by Wood1996  as

$f_{1}(t)=t-t^{3}/3,$ (49)

with $t\in[3.1,20.4]$ and sampling rate $f_{s}=20$ Hz,

$f_{2}(t)=\sin(t)+\sin(2t/3),$ (50)

with $t\in[3.1,20.4]$ and sampling rate $f_{s}=20$ Hz, and

$f_{3}(t)=-\sum_{i=1}^{5}\sin((i+1)t+i),$ (51)

with $t\in[-10,10]$ and sampling rate $f_{s}=20$ Hz. Additionally, additive
noise is included in the signal with $s(t)=Af(t)+\epsilon$ with the additive
noise distribution parameter set to one (e.g. $\sigma=1$ for Gaussian) and
signal amplitude $A$ increment by unit steps starting from zero such that the
$\delta$ is also incremented until reaching a value $\delta/\sigma=2$. At each
$\delta$ we calculate the median lifetime $\tilde{L}$ for 100 trials to
provide a mean $\tilde{L}$ with uncertainty $u_{L}$ as one standard deviation
(see Fig. 6 for the Gaussian additive noise example).

Our goal is to find a function to approximately fit this relationship between
$\delta$ and $\tilde{L}$ for each distribution type. By observation of the
median lifetimes in Fig. 6, we experimentally found an approximate template
function:

$\tilde{L}^{*}=\tilde{L}_{0}e^{-c_{1}{\left(\frac{\delta}{\delta+\tilde{L}}\right)}^{c_{2}}},$
(52)

where $\tilde{L}_{0}$ is the median lifetime when $\delta=0$ or when the
signal is just additive noise $\mathcal{N}$.

Figure 6: Numeric function fitting of Equation (52) to the mean of the median
lifetime $\tilde{L}$ of $f_{i}(t)$ for $i\in[1,3]$ where $\mathcal{N}$ is unit
variance Gaussian additive noise with $\delta\in[0,2]$ being incremented to
understand the effects of signal on the median lifetime.

As shown in Fig. 6, the fitted function shows a very similar quality to the
numerically simulated means of the median lifetimes when the two constants in
Equation (52) were set to $c_{1}\approx 0.845$ and $c_{2}\approx 0.809$ for a
Gaussian additive noise, which were chosen using the BFGS minimization of the
$\ell_{2}$ norm cost function on the residuals when fitting to $\tilde{L}$ for
all three generic functions. Another characteristic of these constants is that
they are approximately independent of the additive noise distribution
parameter, sampling frequency, and time series. The two constants from
Equation (52) are provided in Table 2 for the four distributions investigated
in this work.

Table 2: Constants of Equation (52) for each distribution type investigated in this work with associated uncertainty from ten trials. Distribution | Gussian | Uniform | Rayleigh | Exponential
---|---|---|---|---
$c_{1}$ | $0.845\pm 0.029$ | $0.880\pm 0.017$ | $0.726\pm 0.026$ | $0.436\pm 0.036$
$c_{2}$ | $0.809\pm 0.061$ | $0.639\pm 0.026$ | $0.605\pm 0.054$ | $0.393\pm 0.075$

With these constants, we calculate a multiplication compensation term for the
signal as $R$, which is calculated from Equation (52) as

$R=\frac{\tilde{L}_{0}}{\tilde{L}^{*}}=e^{c_{1}{\left(\frac{\delta}{\delta+\tilde{L},}\right)}^{c_{2}}}$
(53)

which is used to compensate for the effects of signal with
$C_{\alpha}^{*}=RC_{\alpha}$ and $\sigma^{*}=R\sigma$.

Unfortunately, when $s(t)$ is unknown, the $\delta$ parameter used in Equation
(53) can no longer be directly calculated from the time series or sublevel set
persistence diagram. To approximate $\delta$ we use the lifetimes greater than
the initial uncompensated cutoff $C_{\alpha}$ as

$\delta\approx\frac{2}{n}\sum L_{C_{\alpha}},$ (54)

where $L_{C_{\alpha}}$ are the lifetimes greater than $C_{\alpha}$.

To validate the accuracy of Equation (53) with $\delta$ approximated from
Equation (54) we estimate $\sigma$ with and without the signal compensation
$R$ from Equation (53), we use a new time series $x(t)=A\sin(\pi t)+\epsilon$
with $\mathcal{N}$ being a Gaussian distribution with unit variance and $A$
incremented to change $\delta\in[0,2]$ for $100$ trials at each $\delta$. As
shown in Fig. 7, the true $\sigma=1$ and the estimated $\sigma$ without
compensation from Equation (39) shows an underestimate as $\delta$ increases
until plateauing around $\delta/\sigma\approx 1$, which would cause for a
cutoff that may not capture all of the lifetimes associated with noise.
However, the signal compensated distribution parameter $\sigma^{*}$ shows an
accurate estimation of $\sigma$ even as $\delta$ becomes significantly large.
This example demonstrates the importance of signal compensation for an
accurate cutoff and distribution parameter estimation.

Figure 7: Demonstration of distribution parameter $\sigma$ estimation of
Gaussian additive noise in $x(t)=A\sin(\pi t)+\mathcal{N}$ using the median
lifetime with and without signal compensation as $\sigma$ and $\sigma^{*}$,
respectively.

##  7 Comparisons and Examples

We first make a comparison between our ANAPT method for estimating a suitable
cutoff in the sublevel set persistence diagram to two common approaches.
Namely, persistent entropy for labeling persistence pairs as signal or noise
based on their relative entropy and bootstrap resampling to generate a cutoff
based on the bottleneck distance between persistence diagrams. In our second
analysis we experimentally validate our ANAPT method for three examples: (1) a
numeric simulation of the noise types with their associated cutoffs
calculated; (2) a numerically generated quasiperiodic time series with
additive noise; and (3) experimental data from a pendulum with natural noise
of an unknown distribution. These three examples demonstrate the wide breadth
of noise distributions for selecting a cutoff, the effects of signal in
estimating an appropriate cutoff, and the ability for the cutoff estimation
method to function for experimental data with an unknown distribution type.

###  7.1 Comparison to Standard Methods

In this section, we analyze the effectiveness of our ANAPT method of the
sublevel sets for estimating a cutoff compared to other standard techniques
for separating topological noise from signal. For our first comparison we
separate topological noise from signal using persistent entropy Rucco2016 ;
Atienza2017 , bootstrap resampling of the bottleneck distance Fasy2014 , and
ANAPT methods for an example signal.

_Persistent Entropy_. Persistent entropy uses the Shannon entropy of the
lifetimes from a persistence diagram to measure the significance of a
persistence pair Rucco2016 . By comparing the entropy of a modified set of
lifetimes to a relative entropy level, the identification of a feature as
either noise or signal is determined. Details on how to implement persistent
entropy for separating signal from topological noise can be found in
Atienza2017 .

_Bootstrap Resampling_. Bootstrap resampling is based on resampling a dataset
with replacement; here, we bootstrap residuals to estimate the additive noise.
Since we are working with dependent time series data, we use a frequency
domain function fitting method to calculate these residuals. The bootstrap
resampling procedure is as follows:

  1. 1.

The data is a sequence of real-valued observations, which we can model as a
function $x\colon\\{t_{1},t_{2},\ldots,t_{n}\\}\to\mathbb{R}$. For each
$t\in\\{t_{1},t_{2},\ldots t_{n}\\}$, this reading is actually a sum of a true
underlying signal and some (zero-centered) additive noise:
$x(t)=s(t)+\epsilon$. However, we do not know exactly what $s(t)$ and
$\epsilon$ are individually.

  2. 2.

Compute $\mathcal{D}(x)$.

  3. 3.

Filter the signal using a Butterworth low-pass filter of the Fast Fourier
Transform (FFT) with a frequency threshold of twice the sampling rate. The
resulting filtered signal $\hat{s}(t)$ is an approximation of $s(t)$.

  4. 4.

For each $t$, set $\hat{\epsilon}(t)=\hat{s}(t)-x(t)$ as an estimate of
$\epsilon$

  5. 5.

Sampling over all $\hat{\epsilon}$ with replacement, we can obtain:
$\epsilon^{*}_{1},\epsilon^{*}_{2},\ldots,\epsilon^{*}_{n}$ (our bootstrapped
noise).

  6. 6.

Our new signal is $x^{*}(t)=\hat{s}(t)+\epsilon^{*}$.

  7. 7.

Compute the persistence diagram $\mathcal{D}(x^{*})$.

  8. 8.

Compute the bottleneck distance between $\mathcal{D}(x)$ and
$\mathcal{D}(x^{*})$.

  9. 9.

Repeat Steps 5 through 8 to get $N$ bootstrapped distances. (In our
experiments, we used $N=1000$.

At the end of this procedure, we have a set of $N$ bootstrapped distances
between the persistence diagram computed from the data, $\mathcal{D}(x)$, and
the diagrams obtained through the bootstrap; let
$\\{d_{1}^{*},d_{2}^{*},\ldots,d_{N}^{*}\\}$ be this set of distances. Letting
$\alpha=0.05$, let $\bar{d}^{*}_{95\%}$ denote the $95\%$-ile of these
distances. Then, with $95\%$ confidence, we can say that
$d(\mathcal{D}(x),\mathcal{D}(s))\leq\bar{d}^{*}_{0.05}$. Hence, for each
point in $\mathcal{D}(x)$, we can classify it as either _true signal_ or
_indistinguishable from noise_ by thresholding the liftetime of the points by
$2\bar{d}^{*}_{0.05}$.333The factor of $2$ appears, since the lifetime of a
persistence point is half of its $L_{\infty}$-distance to the diagonal. It
should be noted that for each resampling the end conditions are held constant
to make the bottleneck distance correspond to additive noise only.

_Data_. The example signal used to analyze the effectiveness of each method is
a chaotic solution to the Lorenz system

$\frac{dx}{dt}=\sigma(y-x),\>\frac{dy}{dt}=x(\rho-z)-y,\>\frac{dz}{dt}=xy-\beta
z.$ (55)

The Lorenz system was sampled at a rate of 200 Hz for 200 seconds with system
parameters $\sigma=10.0$, $\beta=8.0/3.0$, and $\rho=181.0$. Only the last
2500 samples or 12.5 seconds of the signal were used to avoid transients.
Additive Gaussian noise with a Signal to Noise Ratio (SNR) of 23 dB was used
to contaminate the signal. For our analysis, we used the time series $x(t)$ as
shown in the top of Fig. 8.

Figure 8: Comparison between the persistent homology, bootstrap resampling,
and ANAPT methods for separating noise from topological signal in the sublevel
set persistence diagram. The example signal is generated from a chaotic Lorenz
system.

_Comparison_. Fig. 8 shows that the persistent entropy method (bottom left
subfigure) does not accurately separate all of the persistence pairs
associated to noise from those of signal. This is due to its design being for
persistence diagrams with relatively few significant points in the persistence
diagram. However, we empirically observed that the persistent entropy method
has a limitation on the ratio of persistence pairs associated to signal to
those associated to noise as demonstrated in Fig. 8. This reduces its
functionality for persistence diagrams where the number of features associated
to signal is close to the number associated to noise.

The second method, bootstrap resampling (bottom center subfigure in Fig. 8),
has a resulting cutoff ($C\approx 7.70$) that narrowly encapsulates all of the
features associated to noise. We hypothesize that this is due to the
bottleneck distance minimum pairing between persistence diagrams is not
necessarily matching persistence pairs that are associated but rather the
closest persistence pair. Unfortunately, the bootstrap method is also
significantly more computationally demanding than the other techniques. This
is due to the need to calculate the bottleneck distance for each persistence
diagram pair and the need to calculate the sublevel set persistence for the
number of resamplings desired (1000 for our example).

Based on the results in Fig. 8, it is clear that the ANAPT method is optimal
for estimating a cutoff for seperating noise from topological signal. Our
ANAPT method accurately seperates all persistence pairs associated to noise
from signal with a cutoff $C_{\alpha}\approx 12.96$. This is in comparison to
the optimal cutoff of $C_{\alpha}\approx 12.70$ based on knowing $\sigma$ and
the theory in Section 3.2.

Figure 9: Comparison of distribution parameter estimation techniques (low-Pass
filter residuals, cubic spline residuals, and sublevel set persistence) for
estimating $\sigma$ of the additive Gaussian noise contaminating the $x(t)$
solution of the Lorenz system in Eq. 55.

_Estimating Standard Deviation of (Gaussian) Additive Noise_.

To make further conclusions on the comparative performance of our ANAPT
method, we evaluate residual based methods for estimating the standard
deviation of the additive Gaussian noise. The two methods we implement are a
low-pass Butterworth Filter with a frequency cutoff of twice the sampling
frequency and a cubic spline fit on a down sampled version of the signal for
estimating the residuals from the noisy signal. The downsampling rate for the
cubic spline fitting was determined such that the time delay estimated from
the autocorrelation function Box2015  is two, which avoids aliasing and over
fitting. To make comparisons between the residual based methods and the ANAPT
method we again use the chaotic Lorenz signal from Equation (55). However, we
vary the SNR from $0$ to $40$ dB for Gaussian additive noise. An SNR of $0$ dB
represents a noise amplitude equivalent to the amplitude of the signal, while
a value of $40$ dB represents a relatively low amount of additive noise. Our
results are shown in Fig. 9. For each SNR 10 trials are run with the mean and
standard deviation as the error bars. As demonstrated in Fig. 9, the cubic
spline method for function fitting and analyzing the residuals tends to
overestimate $\sigma$. In comparison, both the low-pass Butterworth filter and
the ANAPT methods closely estimate $\sigma$. However the low-pass Butterworth
filter tends to slightly underestimate $\sigma$. This underestimation can
cause for a less conservative cutoff in comparison to the slightly
overestimated $\sigma$ from the ANAPT method. In our opinion, either the low-
pass filter or the ANAPT methods could be used, but due to the additional
computation of filtering the signal, we suggest using our ANAPT method. It
could be beneficial to calculate both $\sigma$ estimates to verify that the
estimate is accurate.

###  7.2 Numerically Simulated Noise Models

The first example is of time series of only additive noise such that
$\epsilon=\\{\epsilon_{1},\epsilon_{2},\ldots,\epsilon_{n}\\}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}$. For this example, we simply calculate the cutoff
from the distributions corresponding cutoff Equation (Equation (27), Equation
(31), Equation (34), and Equation (37)) and show that the cutoff is greater
than all sublevel set lifetimes (see Fig. 10). For this example, the
distribution parameters are $\sigma=\Delta=\lambda=\sigma=1$ for the Gaussian,
uniform, Rayleigh, and exponential distributions, respectively. Fig. 10 shows
the histogram estimates of the probability densities of the additive noise
types on the top row, and on the bottom row are the approximate probability
distributions of the corresponding lifetimes. For each of the lifetime
distributions, the calculated cutoff $C_{\alpha}$ is shown, where
$\alpha=0.001$ and $n=10^{5}$ data points.

Figure 10: Numerical simulations with $10^{5}$ data points for each of the
four investigated noise models with probability densities on the top row and
lifetime probability densities on the bottom row with their associated cutoff
$C_{\alpha}$. The distribution parameters were set as
$\sigma=\Delta=\lambda=\sigma=1$ and were estimated as $\sigma\approx 0.996$,
$\Delta\approx 1.013$, $\lambda\approx 1.009$, and $\sigma\approx 0.998$ for
the Gaussian, uniform, Rayleigh, and exponential distributions, respectively

As shown, the cutoff is greater than the maximum lifetime in the probability
density histograms of the lifetimes for each type of additive noise
investigated, which demonstrates that this method for selecting a cutoff using
the statistics of the sublevel set lifetimes is suitable for selecting a
cutoff if the distribution is known. However, this example does not show what
effects dependency in the time series has on the distribution. In the next
example, we examine the effect of a dependency or signal in the time series on
the cutoff calculation and its accuracy.

###  7.3 Numerically Simulated Signal with Additive Noise

The second example used is a sinusoidal time series with additive noise (see
left column of Fig. 11) calculated as

$x(t)=A(\sin(\pi t)+\sin(t))+\epsilon,$ (56)

where $A=10$ is the signal amplitude, $\epsilon$ is additive noise with a
normal distribution with zero mean $\mu=0$ and standard deviation
$\sigma\in[0.2,1,2,4]$, and $t\in[0,15]$ with a sampling rate of $f_{s}=40$
Hz.

Figure 11: Example cutoff calculation for time series $x(t)=A(\sin(\pi
t)+\sin(t))+\epsilon$, where $\epsilon$ is additive noise from a Gaussian
distribution with zero mean and four different standard deviations as
$\sigma\in[0.2,1,2,4]$. The resulting sublevel set persistence diagrams with
cutoff $C_{\alpha}^{*}$ are shown to the right.

For a signal like this there is a need to multiply by a compensation term to
the cutoff and distribution parameter as discussed in Section 6. This is due
to the slope of the function reducing the lifetimes associated to additive
noise and the effect of lifetimes associated to signal. For the signal with
$\sigma=1$ and calculating the median lifetime and setting $\alpha=0.001$, we
calculate a cutoff $C_{\alpha}\approx 5.36$ and estimated $\sigma\approx
0.751$ for this example. However, if we apply the compensation multiplier $R$
from Equation (53) and approximate delta from Equation (54) as $\Delta\approx
0.953$ we find a cutoff $C_{\alpha}^{*}\approx 7.54$ and $\sigma^{*}\approx
1.05$, which is a much more accurate cutoff and approximate distribution
parameter with the actual $\sigma=1$. Similar improvements were found for the
other noise levels when incorporating signal compensation.

As shown in the example with $\sigma=4$ at the bottom of Fig. 11, the cutoff
will eventually encapsulate persistence pairs associated to signal if the
amplitude of the noise is significant in comparison to the amplitude of the
signal.

This example showed how the method for selecting a cutoff is applicable to a
time series with additive noise and an underlying signal that is not simply
sinuisodal, but quasiperiodic. However, we still knew what the additive noise
distribution type was. In our next example, we look at experimental data,
where the noise does not have a known distribution.

###  7.4 Experimental Example

The experimental data used in this example is from a single pendulum with an
encoder to track the position of the pendulum arm. A full design document for
the pendulum is available through GitHub444https://github.com/Khasawneh-
Lab/simple_pendulum. The data collected is from a free drop response of the
pendulum where the amplitude of the pendulum swing is decreasing due to
natural damping mechanisms. The recorded data from the encoder was translated
from a voltage output to radians, but this calibration and dimensional
analysis are not covered in this manuscript, but rather in the documentation
report with a full uncertainty analysis in Petrushenko2017 . In Fig. 12, a
section of the time series sampled at 100 Hz is shown with a subsection
highlighted in green when the time series no longer has significant
dependencies such that $x(t)\approx\mathcal{N}$. Additionally, in Fig. 12, the
resulting persistence diagram with the corresponding noise cutoff with and
without signal adjustment are shown.

Figure 12: Experimental time series from free drop response of single pendulum
with corresponding cutoffs in the persistence diagram and time-order lifetimes
plot.

To generate this cutoff, an additive noise distribution model is needed. Since
the data was from an experiment this distribution is not exactly known. To
better understand and assume a distribution we used the section of the time
series ($t\in 24,28]$) highlighted in green in the top left of Fig. 12 where
there is approximately no dependency in the data such that
$x(t)\approx\epsilon$. Using this section of data, we generated a histogram of
the additive noise as shown in Fig. 13.

This distribution, qualitatively, seems to best match a Rayleigh distribution,
which was used to determine the cutoff. The fitted Rayleigh distribution is
also shown in Fig. 13, where $\sigma$ was estimated from the median lifetime
using Equation (43) to estimate $\sigma\approx 0.028$. Using this fitted
distribution parameter, a suitable cutoff would be $C_{\alpha}\approx 0.1293$
using Equation (34). However, we are interested in calculating the cutoff
using the entire signal with $t\in[15,28]$ as shown in Fig. 12 since it is not
always possible to have a section of time series that does not have signal
dependencies. If we directly apply Equation (34) we calculate a cutoff as
$C_{\alpha}\approx 0.0905$, which is lower than the desired cutoff
$C_{\alpha}\approx 0.1293$. To account for this, we implement the signal
compensation discussed in Section 6 and calculated the cutoff as
$C_{\alpha}^{*}\approx 0.1365$ which is slightly conservative, but a better
estimate than without signal compensation.

Figure 13: Histogram of the highlighted section of signal in Fig. 12 with
fitted Rayleigh distribution.

This example highlights the importance of signal compensation in the cutoff
calculation and demonstrates the functionality of the ANAPT cutoff estimation
method for experimental data. Additionally, we showed how one of the analyzed
distributions (Rayleigh) can be used to approximate the distribution of the
additive noise and be used to calculate an accurate cutoff.

##  8 Conclusion

In this work we provided a novel statistical analysis of sublevel set
persistence diagrams from a single variable function. Specifically, we show
how to calculate a cutoff $C_{\alpha}$ for separating persistence pairs
associated to noise from those of signal in the persistence diagram. This
method is based on an assumption that the form of the underlying probability
distribution of the additive noise is known (e.g. Gaussian distribution). This
method has several advantages to other methods including a reduced
computational time in comparison to resampling techniques since only a single
persistence diagram needs to be calculated, no signal filtration required,
widely applicable, and functions with significant amounts of additive noise.
Additionally, this method is essentially parameter free with the only
parameter being the confidence level $\alpha$, which is suggested to be
approximately 0.05. These features make ANAPT especially useful for
applications such as extrema detection and peak prominence analysis. While
ANAPT does use TDA, it is easy to implement and use with our open-source
python code available555https://github.com/lizliz/teaspoon.

To determine the effectiveness of the cutoff, we first showed how ANAPT
functions correctly for time series of pure additive noise (i.e.
$x(t)=\mathcal{N}$) with all sublevel set persistence pairs captured by the
cutoff. Next, we tested ANAPT on both an informative, simulated time series
and an experimental time series from a pendulum. For the simulated time series
we added Gaussian additive noise and showed how ANAPT accurately determined a
suitable cutoff by separating the correct sublevel set lifetimes. For the
experimental time series study, we did not know the distribution of the
additive noise and proposed approximating it as the Rayleigh distribution.
This approximation allowed for an accurate calculation of a cutoff and showed
that ANAPT functions even without the exact probability distribution of the
additive noise known.

In the future we plan to investigate a few new directions. The first is to
extend the statistical analysis to sublevel set persistence for multi-variable
functions, where a new theoretical distribution of the birth and death times
needs to be developed. The second improvement would be to develop a new method
to approximate the probability distribution function of the underlying
additive noise to alleviate the need for the distribution to be specified. The
challenge with this new development would be creating a library of
distribution functions to test against, which was out of the scope of this
work.

###### Acknowledgements.

##  9 Acknowledgements

This material is based upon work supported by the National Science Foundation
under grant numbers CMMI-1759823, DMS-1759824, and DMS-1854336.

## Conflict of interest

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## References

  * [1] Nieves Atienza, Rocio Gonzalez-Diaz, and Matteo Rucco.  Persistent entropy for separating topological features from noise in vietoris-rips complexes.  Journal of Intelligent Information Systems, 52(3):637–655, jul 2017\.
  * [2] Peter Bühlmann and Peter Buhlmann.  Sieve bootstrap for time series.  Bernoulli, 3(2):123, jun 1997.
  * [3] George EP Box, Gwilym M Jenkins, Gregory C Reinsel, and Greta M Ljung.  Time series analysis: forecasting and control.  John Wiley & Sons, 2015.
  * [4] Peter Bubenik.  Statistical topological data analysis using persistence landscapes.  Journal of Machine Learning Research, 16:77–102, 2015.
  * [5] Gunnar Carlsson, Jackson Gorham, Matthew Kahle, and Jeremy Mason.  Computational topology for configuration spaces of hard disks.  Physical Review E, 85(1), jan 2012.
  * [6] Frédéric Chazal, Vin De Silva, Marc Glisse, and Steve Oudot.  The Structure and Stability of Persistence Modules.  Springer, 2016.
  * [7] Frédéric Chazal, Brittany Terese Fasy, Fabrizio Lecci, Bertrand Michel, Alessandro Rinaldo, and Larry Wasserman.  Robust topological inference: Distance to a measure and kernel distance.  J. Mach. Learn. Res., 18:5845–84, 2018.
  * [8] Frédéric Chazal, Brittany Terese Fasy, Fabrizio Lecci, Alessandro Rinaldo, Aarti Singh, and Larry Wasserman.  On the bootstrap for persistence diagrams and landscapes.  arXiv preprint arXiv:1311.0376, 2013.
  * [9] S. Chowdhury and F. Mémoli.  Convergence of hierarchical clustering and persistent homology methods on directed networks.  ArXiv, abs/1711.04211, 2017.
  * [10] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer.  Stability of persistence diagrams.  Discrete & Computational Geometry, 37(1):103–120, dec 2006.
  * [11] S. Czesla, T. Molle, and J. H. M. M. Schmitt.  A posteriori noise estimation in variable data sets.  Astronomy & Astrophysics, 609:A39, jan 2018.
  * [12] Cecil Jose A. Delfinado and Herbert Edelsbrunner.  An incremental algorithm for Betti numbers of simplicial complexes on the 3-sphere.  Computer Aided Geometric Design, 12(7):771–784, 1995.
  * [13] Meryll Dindin, Yuhei Umeda, and Frederic Chazal.  Topological data analysis for arrhythmia detection through modular neural networks.  In Advances in Artificial Intelligence, pages 177–188. Springer International Publishing, 2020.
  * [14] Edelsbrunner, Letscher, and Zomorodian.  Topological persistence and simplification.  Discrete & Computational Geometry, 28(4):511–533, 2002.
  * [15] Herbert Edelsbrunner and John Harer.  Persistent homology-a survey.  Contemporary mathematics, 453:257–282, 2008.
  * [16] Herbert Edelsbrunner and John Harer.  Computational Topology - an Introduction. American Mathematical Society, 2010.
  * [17] Brittany Terese Fasy, Fabrizio Lecci, Alessandro Rinaldo, Larry Wasserman, Sivaraman Balakrishnan, Aarti Singh, et al.  Confidence sets for persistence diagrams.  The Annals of Statistics, 42(6):2301–2339, 2014.
  * [18] Shafie Gholizadeh and Wlodek Zadrozny.  A short survey of topological data analysis in time series and systems analysis.  arXiv preprint arXiv:1809.10745, 2018.
  * [19] J. Hu, J.B. Gao, and K.D. White.  Estimating measurement noise in a time series by exploiting nonstationarity.  Chaos, Solitons & Fractals, 22(4):807–819, nov 2004.
  * [20] Firas A. Khasawneh and Elizabeth Munch.  Utilizing Topological Data Analysis for Studying Signals of Time-Delay Systems, pages 93–106.  Springer International Publishing, Cham, 2017.
  * [21] Firas A. Khasawneh and Elizabeth Munch.  Topological data analysis for true step detection in periodic piecewise constant signals.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 474(2218):20180027, oct 2018.
  * [22] Firas A. Khasawneh, Elizabeth Munch, and Jose A. Perea.  Chatter classification in turning using machine learning and topological data analysis.  IFAC-PapersOnLine, 51(14):195–200, 2018.
  * [23] Peter Lawson, Andrew B. Sholl, J. Quincy Brown, Brittany Terese Fasy, and Carola Wenk.  Persistent homology for the quantitative evaluation of architectural features in prostate cancer histology.  Scientific Reports, 9(1), feb 2019.
  * [24] Hyekyoung Lee, Hyejin Kang, M. K. Chung, Bung-Nyun Kim, and Dong Soo Lee.  Persistent brain network homology from the perspective of dendrogram.  IEEE Transactions on Medical Imaging, 31(12):2267–2277, dec 2012\.
  * [25] Andreas Otto Melih C. Yesilli, Firas A. Khasawneh.  Topological feature vectors for chatter detection in turning processes.  arXiv:1905.08671, 2019.
  * [26] Elizabeth Munch.  A user’s guide to topological data analysis.  Journal of Learning Analytics, 4(2), 2017.
  * [27] Audun Myers and Firas A. Khasawneh.  On the automatic parameter selection for permutation entropy.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(3):033130, mar 2020.
  * [28] Audun Myers, Elizabeth Munch, and Firas A Khasawneh.  Persistent homology of complex networks for dynamic state detection.  arXiv preprint arXiv:1904.07403, 2019.
  * [29] Jose A. Perea.  A brief history of persistence.
  * [30] José A. Perea and John Harer.  Sliding windows and persistence: An application of topological methods to signal analysis.  Foundations of Computational Mathematics, pages 1–40, 2015.
  * [31] David Petrushenko and Firas A. Khasawneh.  Uncertainty propagation of system parameters to the dynamic response: An application to a benchtop pendulum.  In Volume 4B: Dynamics, Vibration, and Control. American Society of Mechanical Engineers, nov 2017.
  * [32] Dimitris N. Politis and Joseph P. Romano.  The stationary bootstrap.  Journal of the American Statistical Association, 89(428):1303–1313, dec 1994.
  * [33] Matteo Rucco, Filippo Castiglione, Emanuela Merelli, and Marco Pettini.  Characterisation of the idiotypic immune network through persistent entropy.  In Proceedings of ECCS 2014, pages 117–128. Springer International Publishing, 2016.
  * [34] Floris Takens.  Detecting strange attractors in turbulence.  In David Rand and Lai-Sang Young, editors, Dynamical Systems and Turbulence, Warwick 1980, volume 898 of Lecture Notes in Mathematics, pages 366–381. Springer Berlin Heidelberg, 1981.
  * [35] Christopher J. Tralie and Jose A. Perea.  (quasi)periodicity quantification in video data, using topology.  SIAM Journal on Imaging Sciences, 11(2):1049–1077, jan 2018.
  * [36] Krzysztof Urbanowicz and Janusz A. Hołyst.  Noise-level estimation of time series using coarse-grained entropy.  Physical Review E, 67(4), apr 2003.
  * [37] Xiang Wan, Wenqian Wang, Jiming Liu, and Tiejun Tong.  Estimating the sample mean and standard deviation from the sample size, median, range and/or interquartile range.  BMC Medical Research Methodology, 14(1), dec 2014.
  * [38] G.R. Wood and B.P. Zhang.  Estimation of the lipschitz constant of a function.  Journal of Global Optimization, 8(1), jan 1996.

##  Appendix A Proofs and Algorithms

In the appendix, we provide an omitted proof and algorithm. Specifically, we
have included the theorem showing the relationship between the mean lifetime
and mean birth and death times of a persistence diagram and the algorithm for
calculating the persistence diagram for the sublevel-set filtration of a time
series.

###  A.1 Proof of Expected Lifetime Equation

The following proof supports a claim made in Section 3.1. In what follows, we
use the notation $\mu_{S}$ to denote the expected value of the distribution
over the multi-set $S$.

######  Theorem A.1 (Expected Lifetime)

Let $\mathcal{D}=\\{(b_{i},d_{i})\\}_{i=1}^{n}$ be a persistence diagram. Let
$B$, $D$, and $L$ be the multi-sets of birth times, death times, and
lifetimes, respectively. Then, the average lifetime is:

$\mu_{L}=\mu_{D}-\mu_{B}.$

###### Proof

By definition, $B={\\{b_{i}\\}}_{i=1}^{n}$, $D={\\{d_{i}\\}}_{i=1}^{n}$, and
$L=\\{d_{i}-b_{i}\\}_{i=1}^{n}$. By definition of mean and of $L$, the mean
lifetime is

$\begin{split}\mu_{L}&=\frac{1}{n}\sum_{i=1}^{n}(d_{i}-b_{i}).\end{split}$
(57)

Expanding the sum to two separate sums and using the commutative property of
addition, we get:

$\begin{split}\mu_{L}&=\frac{1}{n}\sum_{i=1}^{n}(d_{i}-b_{i})\\\
&=\frac{1}{n}\sum_{i=1}^{n}d_{i}-\frac{1}{n}\sum_{i=1}^{n}b_{i}\\\
&=\mu_{D}-\mu_{B},\end{split}$ (58)

where the last equality is by definition of $\mu_{D}$ and $\mu_{B}$. Thus, we
conclude $\mu_{L}=\mu_{D}-\mu_{B}$.

###  A.2 Algorithm

In Algorithm 1, we give an $\Theta(n\log n)$ algorithm for computing the
persistence diagram for the sublevel sets of a real-valued function. The time
series is represented as an array $A$ of $n$ real numbers. This array is
interpreted as a continuous function $f\colon[1,n]\to\mathbb{R}$, where
$f(i)=A[i]$ and $f(x)$ is calculated by linearly interpolating between
$f(\left\lfloor x\right\rfloor)$ and $f(\left\lceil x\right\rceil)$. The
algorithm is similar to that presented in [12], but uses a priorty queue
instead of the union-find data structure.

We begin (Line 1) by creating a linked list of local extrema in $A$, ordered
by their index in $A$. This step takes $\Theta(n)$ time. Each $m\in M$ stores
its value or height as $m.height$ as well as a pointer $m.next$ to the next
element of $M$. In order to avoid a boundary effect, if the first value is a
minimum, then we included it as a minimum (resp., maximum) with height
$-\infty$ (resp., $\infty$). We handle the last value similarly. If $p$ is a
pointer to an element of $M$, we use $M[p]$ to denote the corresponding
element of the linked list.

We then create a priority queue $Q$ in order to keep track of the potential
min/max pairs that we could make (where smaller value indicates higher
priority). The priority queue stores a pair $(ptr,v)$, where $ptr$ is a
pointer to an element of $M$ and $v$ is the priority of $M[ptr]$, in
particular, $v=\left|M[ptr].height-M[ptr].next.height\right|$. $Q.pop$ returns
the (value of the) element of $Q$ with lowest priority, and simultaneously
removes it from $Q$. The construction of the priority queue in Line 1 takes
$\Theta(n\log n)$.

The while loop starting on Line 1 maintains the invariant that when entering
the while loop for the $i$-th time the $i-1$ smallest persistence points have
been calculated. The loop is straightforward, other than the update of $Q$ and
$M$ in Line 1. To do so, we first set $m^{\prime}=m.prev$ and remove $m$ and
$m.next$ from $M$ as well as $m^{\prime}$, $m$, and $m.next$ from $Q$.
Finally, we add $m^{\prime}$ to $Q$ with priority
$|m^{\prime}.height-m^{\prime}.next.height|$ The while loop executes
$\Theta(n)$ times; therefore, the total runtime of this algorithm is
$\Theta(n\log n)$.

For our specific python implementation of this algorithm, we use the
SortedList data structure from the Sorted Containers Python package. Our
implementation is available through the python package
teaspoon666https://lizliz.github.io/teaspoon/.

Data: Array $A$ of $n$ Real Numbers

Result: Persistence Diagram

1 begin

2 $M=$ list of local extrema in $A$ (non-endpoint maxima and minima);

3 $Q=$ priority queue of pointers to elements of $M$;

4 while  _$|M| >3$_ do

5 $m\leftarrow Q.pop$;

6 $b\leftarrow\min\\{m.height,m.next.height\\}$;

7 $d\leftarrow\max\\{m.height,m.next.height\\}$;

8 Add $(b,d)$ to $pairs$;

9 Update $Q$ and $M$;

10

11 end while

12 return _pairs_

13 end

14

Algorithm 1 Zero-Dimensional Persistence Algorithm.