THE THREE GAP THEOREM IN PERSISTENT HOMOLOGY

By

Luis Suarez Salas

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Applied Mathematics—Doctor of Philosophy
Computational Mathematics, Science, and Engineering—Dual Major

2024

ABSTRACT

The time delay (or Sliding Window) embedding is a technique from dynamical systems to recon-

struct attractors from time series data. Recently, descriptors from Topological Data Analysis (TDA)

— specifically, persistence diagrams — have been used to measure the shape of said attractors in

applications including periodicity and quasiperiodicity quantification. Despite their utility, the fast

computation of persistence diagrams of sliding windows embeddings is still poorly understood.

In this work, we present theoretical and computational schemes to approximate the persistence

diagrams of sliding window embeddings from quasiperiodic functions. We do so by combining

the Three Gap Theorem from number theory with the Persistent Künneth formula from TDA, and

derive fast and provably correct persistent homology approximations. The input to our procedure

is the spectrum of the signal, and we provide numerical as well as theoretical evidence of its utility

to capture the shape of toroidal attractors.

We begin our efforts by documenting the stability of using the Three Gap Theorem to compute

persistence diagrams. The theorem relies on the continued fraction expansion of an irrational

parameter. In turn, this expansion yields the k-th convergent which lie at the core of the result.

This relation is leveraged to show that the number of matching continued fraction expansion terms

on two parameter values can be used to bound the bottleneck distance of their corresponding

persistence diagrams. This stability is then extended to the number of matching terms in their

decimal expression. This is valuable since in practice we extract our parameters using the Fast

Fourier Transform (FFT). Our results show exponential decay in bottleneck distance with respect to

matching decimal terms. Ultimately, they validate the reliability of algorithms relying on continued

fraction information, like our method 3G.

In the experiments presented here, our method took less than a second to run, on average. In

linear computational time, it approximates persistence diagrams, which usually take exponential

computational time. We demonstrated its performance by applying it to dynamical systems from a

wide range of scientific disciplines. We are able to successfully approximate persistence diagrams

within known error bounds and our work contributes to the implementation of TDA on larger data

sets.

Copyright by
LUIS SUAREZ SALAS
2024

TABLE OF CONTENTS

CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

1

CHAPTER 2

BACKGROUND . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 11

CHAPTER 3

STABILITY OF THE THREE GAP THEOREM IN PERSISTENCE
DIAGRAMS .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 33

CHAPTER 4

ESTIMATION OF PERSISTENCE DIAGRAMS VIA THE THREE
GAP THEOREM . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 45

CHAPTER 5

APPLICATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 57

CHAPTER 6

CONTRIBUTION AND FUTURE WORK . . . . . . . . . . . . . . . . 69

BIBLIOGRAPHY .

.

.

.

.

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

v

CHAPTER 1

INTRODUCTION

Scientific observation is at the heart of models and theories. Measurements obtained from the real

world guide and validate a framework. For instance, in neuroscience, the Wilson-Cowan equations

were able to recreate oscillatory behavior that had been independently observed in spiking activities

of motor neurons [1]. The model consists of two driving forces: excitatory and inhibitory neurons,

depicted in Figure 1.1. The spiking of these two types of neurons can be directly measured in the

Figure 1.1 We illustrate the oscillatory behavior the Wilson-Cowan equations can replicate. On the
left, we plot different trajectories in the 𝐼 𝐸-plane. On the right, we plot the solution corresponding
to the red trajectory.

real world. This generates a time series which is a measurement of a variable across time. A time

series is capable of capturing critical information that is difficult to retrieve. How much information

can it contain? A task not worth considering, being ill-posed, is aiming to use observed spiking

activities of motor neurons to infer the number of neurons and their size. Rather, reconstructing

qualitative information is within reach. This type of information is what the model aimed for in the

first place as it is better suited for the study of higher functions in the brain: The number of neurons

discards the possibility of measuring them individually, pattern recognition is a global process,

and there is an abundance of local randomness that gains precision on long-range interactions

[1]. Other fields of science care about qualitative information for analogues reasons. Hence, it is

imperative to be able to determine qualitative information of a system using a time series as input.

Topological data analysis (TDA) offers a framework to achieve this that can be applied across

1

different fields of science. However, a current limitation is the computational complexity it takes.

This motivated our work which provides a much more feasible alternative. It relies on computing

an approximation leveraging results in Fourier analysis and number theory. Ultimately, we hope to

alleviate computational constraints, and in turn facilitate implementation to a wider class of data

sets.

Mathematically, we denote a time series as 𝑓

: R → C. Conceptually, we think of it as

the information obtained from measuring a mechanism that changes over time. The spiking of

excitatory and inhibitory neurons, 𝐸 (𝑡) and 𝐼 (𝑡), respectively, are examples of time series data.

They provide a glimpse into a more complex structure: the nervous system. This is an example

of a dynamical system, which pertains to a theory fit to model non-linear processes. A dynamical

system is composed of [2]:

1- A phase space, which consists of all the possible states of the system. This takes the form

of an N-dimensional manifold 𝑀.

2- The dynamical part which captures how the system changes across time. This is a

continuous map Φ : R × 𝑀 → 𝑀 such that Φ(0, 𝑝) = 𝑝 and Φ(𝑠, Φ(𝑡, 𝑝)) = Φ(𝑠 + 𝑡, 𝑝)

for 𝑝 ∈ 𝑀 and 𝑡, 𝑠 ∈ R.

A dynamical system, like the ones shown in Figure 1.2, can be obtained from a set of differential

Figure 1.2 We illustrate two dynamical systems. By looking at their phase space, without the
need of any coordinate axis information, we can conclude they exhibit periodicity.
Indeed, the
convergence of trajectories to the red trajectory indicates the presence of a toroidal attractor.

2

equations like the Wilson-Cowan equations. The equations represent the constraints of the system.

For example in physics, they could be the ones governing the motion of an object. In neuroscience,

they could be the ones describing the action potential interaction of neurons. The end result is

the evolution of a system across time described entirely by equations. Nevertheless, for many

dynamical systems of interest, the governing equations are unknown. All that can be extracted are

measurements of the system. This leads to the task of reconstructing the underlying dynamical

system using time series data. The endeavour of reconstructing an unknown dynamical system

using observational data has been treated by Takens’ Theorem [3]. It allows for the use of signals

such as 𝐸 (𝑡) or 𝐼 (𝑡) to recover behavior in phase space such as the red trajectories shown in Figure

1.1 and 1.2. Concretely, it demonstrates that almost any observation from an unknown dynamical

Figure 1.3 Top row: On the left we plot the butterfly wings generated by the Lorenz system. On
the right we show a solution of the system. Bottom row: We can reconstruct the Lorenz attractor
𝑥(𝑡), shown on the left, and 𝑦(𝑡), shown on the right, thanks to Takens’ Theorem.

3

system can be used to recover important information from the system: the topology of attractors.

In other words, we start with an unknown topological metric space, sampled by a trajectory in the

phase space pf the dynamical system, and are able to reconstruct a topologically equivalent copy

in Euclidean Space, like shown in Figure 1.3. This is of great significance as the shape traced in

the phase space provides qualitative information of the system [4]. In particular, it can indicate the

presence of recurrent behavior such as periodicity or quasiperiodicity [4, 5, 6, 7]. Indeed, as shown

in Figure 1.1, the red trajectory corresponds to oscillations of 𝐸 (𝑡) and 𝐼 (𝑡).

Quasiperiodicity is a repetitive behavior generalizing the notion of periodicity, see Figure 1.4.

It does so by incorporating at least two Q-linearly independent frequencies, i.e. incommensurate.

This type of signal is amply documented in the scientific literature. It has been observed in crystal

formations, studies of mammal vocalizations, fMRI scans obtained from mice, and electrocardio-

grams [8, 9, 10, 11]. The list goes on. Furthermore, it is known to be an intermediary state between

chaos and stability [12, 13]. Thus, its presence is of great interest across scientific fields. Fortu-

nately, we are able to better understand this behavior using tools from topological data analysis.

Indeed, persistence diagrams have been successfully used to detect the presence of quasiperiodicity

in signals [5, 6, 7].

The method used to detect quasiperiodic signals, depicted in Figure 1.4, is the Sliding Window

Embedding Technique (SW). This takes a time series 𝑓 and creates an embedding

𝑆𝑊𝑑,𝜏 𝑓 (𝑡) =















𝑓 (𝑡)

𝑓 (𝑡 + 𝜏)
...

𝑓 (𝑡 + 𝑑𝜏)















∈ C𝑑+1.

By Takens’ Theorem, this map can preserve the topology of attractors of the original dynamical

space. Furthermore, it has been shown that the embedding will be dense in a space homeomorphic to

an 𝑁-dimensional torus if and only if 𝑓 is a quasiperiodic signal with 𝑁 incommensurate frequencies

[5]. Thus, quasiperiodicity detection reduces to showing that the shape of the reconstruction is

topologically an 𝑁-dimensional torus. This is an ideal task for persistent homology as it creates

4

Figure 1.4 The SW starts with a signal, followed by a reconstruction of the phase space, and
concludes with the computation of persistent homology. Top row: SW done with a periodic signal.
Bottom row: SW done with a quasiperiodic signal generated by two incommensurate frequencies.

an abstract simplicial complex from a discrete set of points. It then computes the homology of

the complex across a distance parameter. This provides a robust way of estimating the underlying

topology of a set of points. For our purposes, we look for the characterizing homology of the

𝑁-torus. In the periodic case, such as the Wilson-Cowan dynamical system and the ones shown in

Figure 1.2, we have the topology of a 1-dimensional torus, namely a circle.

Much work has been done proving the usefulness of SW [6, 7, 14]. Furthermore, there is a

well defined way of optimizing the embedding parameters 𝑑 and 𝜏 [5]. A current limitation is the

time it takes to compute persistent homology. The first persistence algorithm is of the order 𝑂 (𝑛3)

where 𝑛 is the number of simplices [15]. In practice, one uses Ripser [16], which improves running

time but is still exponential [17]. This proves to be computationally taxing for most data sets, and

motivated our project of creating an approximation method that would provide a much faster and

computationally accessible alternative.

Our method relies on The Three Gap Theorem which in turn benefits from the theory of

5

continued fractions. The latter studies the expansion of an irrational number, say 𝜔, of the form

𝜔 = 𝑎1 +

1

𝑎2 +

1
𝑎3+ 1
𝑎

:= [𝑎1, 𝑎2, 𝑎3, · · · ],

4+···
where 𝑎1 ∈ Z and 𝑎𝑖 ∈ N. There are infinite non-zero terms 𝑎𝑖 when 𝜔 is irrational [18]. Fur-

thermore, these terms capture important information.

Indeed, consider the continued fraction

expansion of the golden ratio:

𝜑 = 1 +

1

1 +

1
1+ 1
1+···

.

One can readily use this expression to argue 𝜑 is the most irrational number. There is much more

than can be said by looking at continued fraction expansions. Their connection with geometry is

deeply rooted [19].

We leverage continued fractions to better understand the distribution of points in the set

𝑆𝜔,𝑇 := {[𝑖𝜔]}𝑇
𝑖=0

,

where [𝑥] = 𝑥 mod 1. By the Three Gap Theorem, the distance between adjacent points in this set

can take at most three different values [20]. Remarkably, this result is independent of 𝜔 and 𝑇.

By identifying the endpoints of [0, 1], we can think of 𝑆𝜔,𝑇 as a sampling of the circle. Thus,

knowing the gaps in this set can directly transfer to knowing its 0-th dimensional persistence

homology as depicted in Figure 1.5. We then make the jump to the 𝑁-torus by taking a Cartesian

product of the form

𝑁
(cid:214)

𝑖=0

𝑆𝜔𝑖,𝑇 .

We can also compute persistence homology in this case thanks to the Künneth Formula in persistence

homology [14]. For the study of quasiperiodicity, this provides us with an alternative to using Ripser

to compute persistence diagrams. Indeed, we can use the fast fourier transform of time series data

to estimate the incommensurate frequencies. That this will work for general quasiperiodic signals

is thanks to the results in [5]. We are thus in a position to leverage frequency information thanks to

The Three Gap Theorem. This makes it possible to compute an approximation that includes error

bounds without the use of Ripser at any step, resulting in linear time complexity, see Figure 4.5.

6

Figure 1.5 We illustrate the result of the three gap theorem and how they translates directly to results
in persistence homology. Left column: Depiction of the sets 𝑆𝜋−1,17 (top) and 𝑆√
5,17 (bottom).
Middle column: Corresponding persistence barcodes. Right column: Corresponding persistence
diagrams. We note 𝑑2(𝜋 − 1,

5) = 0.094.

√

1.1 Outline and Overview

We move on to detailing the structure of this dissertation. The next section will cover background

definitions followed by our stability results for the use of The Three Gap Theorem in persistence

computations, then we will present our approximation method, 3G, and we will conclude with

applications of it to the study of dynamical systems.

1.1.1 Background

Persistent homology is the main tool used for our work. Thus, we begin by presenting a

complete treatment of definitions and concepts. We also include important theorems and results

we use in our method. Ultimately, persistent homology is intended to investigate reconstructed

dynamical systems. This leads to the Sliding Window Embedding, which is how we take a time

series and create a point cloud in a higher dimension. We then compute its persistence homology to

analyze its shape. This is followed by a brief subsection on Fourier series, in particular we present

the multidimensional case and relate it to persistence diagrams. We then introduce the theory of

7

continued fraction expansions. This leads naturally to The Three Gap Theorem. We end with a

schematic of our method.

1.1.2 Stability of The Three Gap Theorem in Persistence Diagrams

Our work begins by presenting theoretical justifications for the use of The Three Gap Theorem

in persistence homology. We do so by looking at the stability of the bottleneck distance of the

corresponding persistence diagrams from sampling parameters. Our work shows conditions that

justify the use of continued fraction expansion of decimals obtained from numerical approximations.

In particular, we show that for two parameters with the same first 𝑗 decimals:

Proposition 3.2.8 Suppose 𝜔 and ¯𝜔 are two irrational numbers having the same first 𝑗 decimals,

denoted as

𝑥 𝑗 = ¯𝑥 𝑗 ,

for some 𝑗 ∈ N. If 0 < 𝑇 < 𝑞𝑘 𝑗 (𝜔)−1, where 𝑞𝑛 (𝜔) denotes the n-th convergent of 𝜔, then for 𝑖 ≥ 0

the bottleneck distance 𝑑𝐵 is bounded by

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

1
𝑞𝑘 𝑗 (𝜔)

.

Furthermore, we can show the error decays exponentially in probability

Theorem 3.2.10 For 0 < 𝜖 < 𝑧0 there exists positive constants 𝐶, 𝜆 (depending on 𝜖) with 0 < 𝜆 < 1

such that for all integers 𝑗 ≥ 1 for which

𝑥 𝑗 = ¯𝑥 𝑗 ,

and 𝑇 < 𝑞𝑘 𝑗 (𝜔)−1, if 𝑘 𝑗 (𝜔) ≥ 𝑗 𝑧0, then for 𝑖 ≥ 0

(cid:16)

𝑃

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

(cid:17)

1
𝑗 (𝜖 + 𝑧0)

≤ 𝐶𝜆 𝑗 ,

otherwise

(cid:16)

𝑃

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

(cid:17)

1
𝑗 (𝑧0 − 𝜖)

≥ 1 − 𝐶𝜆 𝑗 .

These results are corroborated by our applications at the end of the section. They illustrate

numerical stability that carries over to our approximation method based on The Three Gap Theorem.

8

1.1.3 Estimation of Persistence Diagrams Via The Three Gap Theorem

In this section we present our approximation method 3G. It is based on theoretical guarantees

of The Three Gap Theorem described above and the Künneth formula in persistence homology.

Specifically, we are able to take a quasiperiodic signal and perform computations using each

incommensurate frequency in parallel. We then put everything together to obtain an approximation

to the persistence diagram of the Sliding Window point cloud of interest. We obtain error bounds

with our approximation and demonstrate it computes much faster. An outline of our contribution

is as follows:

1. Start with an observation signal 𝑓 .

2. Reconstruct the phase space by computing

{𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇

(cid:0)(cid:0)(cid:64)(cid:64)3. Compute dgm𝑅

𝑡=0.
𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2) using Ripser.

3. (a) Use the FFT to retrieve the frequencies of 𝑓 and

then use them to compute continued fraction expansions.

(b) Use them as shown in Section 4.2 and then apply

the results from Section 4.3 to approximate

dgm𝑅

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2).

Figure 1.6 Steps involved in SW. Our method provides an alternative to number 3.

1.1.4 Applications

We conclude by considering different quasiperiodic dynamical systems. Using the SW we

are able to verify this behavior. We then apply our 3G method to approximate the persistence

diagrams obtained from SW. We show our method computes much faster and plot the resulting

persistence diagrams with their corresponding error bounds. Shown below is the result obtained

when considering the pendulum attached to a sliding block. Understanding this system translates

to improvement in earthquake damping technology for skyscrapers [21, 22].

9

Figure 1.7 Top row: Depiction of anti-earthquake technology taken from [23]. The system on
the right models the dynamics of this technology. Middle row: Plot of (𝑥, (cid:164)𝑥) from the pendulum
attached to a sliding block. The solution 𝑥(𝑡) we used for SW is plotted in the middle followed
by the modulus of its Fast Fourier Transform. Bottom row: Persistence diagrams showing our
proposed approximation 𝐾3𝐺. The rectangle depicts the theoretical approximation bound.

10

CHAPTER 2

BACKGROUND

In this section we cover results and definitions needed to present our results. For the most part, our

work pertains to persistent homology. Moreover, our contribution involves theoretical guarantees

from number theory and Fourier analysis. We thus attempt to present these separate elements in

a concise but complete fashion. We begin by thoroughly defining persistent homology and the

important results we rely on. This is followed by a definition of dynamical systems. Our main

aim in this section is to highlight the reconstruction theorem of Takens and its connection to the

Sliding Window Embedding technique (SW). We then move to recall the Fourier transform and

detail its connection to SW. This is followed by the presentation of continued fraction expansions

and their properties. This will naturally lead to The Three Gap Theorem, the pillar of our work.

We conclude with a schematic of our method 3G that highlights the integration of the different

components involved.

2.1 Persistent Homology

This section covers the abstract notion of shape we can assign to discrete data. We begin by

presenting the main object to which we can assign a shape to. Most definitions were taken from

[24, 5].

Definition 2.1.1. Given a set X, an abstract simplicial complex K with vertices in X, is a collection

of subsets of P(X) such that

1. Every 𝜎 ∈ 𝐾 is finite,

2. For any 𝜎 ∈ 𝐾, if 𝜏 ≠ ∅ and 𝜏 ⊂ 𝜎 then 𝜏 ∈ 𝐾.

The sets in 𝐾 are called simplices and the dimension of a simplex is defined as dim(𝜎) = |𝜎| −1.

We also define dim(𝐾) = max{dim(𝜎)}. A face of 𝜎 is any non empty proper set of it. The 𝑛-th

skeleton of 𝐾 is denoted by 𝐾 (𝑛) := {𝜎 ∈ 𝐾 : dim(𝜎) ≤ 𝑛}. In particular, 𝐾 (0) is called the set of

vertices of 𝐾 and 𝐾 (1)/𝐾 (0) is called the set of edges of 𝐾. A subcomplex 𝐿 of 𝐾, which we denote

𝐿 ⊂ 𝐾, is an abstract simplicial complex for which 𝐿 (𝑛) ⊂ 𝐾 (𝑛) for all 𝑛 ∈ Z.

11

Example 2.1.2. Let X be a finite set. One can readily check that the power set, 𝑃(𝑋), generates an

abstract simplicial complex.

Example 2.1.3. Any partially ordered set X can form an abstract simplical complex. Indeed, the

faces can be defined as the totally ordered subsets of X. This abstract simplicial complex is known

as the order complex of X.

We can construct an abstract simplicial complex from any discrete set of points in a metric

space. Indeed, this is achieved by constructing the Rips complex which treats each data point in

a metric space, say 𝑥𝑖 ∈ (𝑋, 𝑑), as a vertex. Edges {𝑥0, 𝑥1} are created between two vertices if

their distance is less than some fixed 𝜖 > 0, i.e. 𝑑 (𝑥0, 𝑥1) < 𝜖. Connecting three vertices creates a

triangle {𝑥0, 𝑥1, 𝑥2}, four a tetrahedron {𝑥0, 𝑥1, 𝑥2, 𝑥3}, and so on. Once can think of the result as a

triangulated space.

Definition 2.1.4. Given a finite metric space (X,d) and a real number 𝜖 > 0 we define the Rips

complex 𝑅𝜖 (𝑋, 𝑑) to be the simplicial complex defined as

(cid:8) {𝑥0, . . . , 𝑥𝑛} ⊂ 𝑋 | 𝑑 (𝑥𝑖, 𝑥 𝑗 ) < 𝜖 ∀0≤𝑖, 𝑗 ≤𝑛 (cid:9).

With this construction we can begin using a powerful shape descriptor: homology. This tool

assigns a topological invariant to an abstract simplicial complex. For instance, the 0-th dimensional

homology corresponds to connected components. The 1-st dimensional homology detects loops

and the 2-nd dimensional homology closed surfaces that are hollow. In general the n-th dimensional

homology detects n-dimensional holes. We first present the algebraic structures required.

Let 𝐾 be a simplicial complex and 𝐺 and abelian group.

Definition 2.1.5. The n-th chain group 𝐶𝑛 (𝐾; 𝐺) of K with coefficients in G is

𝐶𝑛 (𝐾; 𝐺) :=

(cid:26) ∑︁

𝑖∈𝐼

𝑔𝑖𝜎𝑖 : 𝜎𝑖 ∈ 𝐾 (𝑛) \ 𝐾 (𝑛−1), 𝑔𝑖 ∈ 𝐺 and 𝑔𝑖 ≠ 0 for only finitely many 𝑖 ∈ 𝐼

(cid:27)

.

We refer to a 𝜏 ∈ 𝐶𝑛 (𝐾; 𝐺) as a 𝑛-chain.

12

Definition 2.1.6. The n-th boundary map 𝜕𝑛 : 𝐶𝑛 (𝐾; 𝐺) −→ 𝐶𝑛−1(𝐾; 𝐺) is a group homomorphism
defined for any 𝜎 = [𝑣0, . . . , 𝑣𝑛] ∈ 𝐾 (𝑛)/𝐾 (𝑛)−1 as

𝜕𝑛 (𝜎) =

𝑛
∑︁

𝑖=0

(−1)𝑛 [𝑣0, . . . , ˆ𝑣𝑖, . . . , 𝑣𝑛]

where [𝑣0, . . . , ˆ𝑣𝑖, . . . , 𝑣𝑛] denotes the (𝑛 − 1)-th face of 𝜎 obtained from deleting the vertex 𝑣𝑖

from the set {𝑣0, . . . , 𝑣𝑛}.

Definition 2.1.7. A chain complex is a collection of 𝐶∗ = {𝐶𝑖, 𝑓𝑖}𝑖 of groups 𝐶𝑖 and morphisms

𝑓𝑖 : 𝐶𝑖 −→ 𝐶𝑖−1 such that 𝑓𝑖−1 𝑓𝑖 = 0.

Definition 2.1.8. The i-th homology group of a chain complex 𝐶∗ with coefficients in a group G is

defined as

𝐻𝑖 (𝐶∗; 𝐺) :=

ker( 𝑓𝑖)
𝑖𝑚𝑔( 𝑓𝑖+1)

.

Our work restricts to the chain complex 𝐶∗ = {𝐶𝑖 (𝐾; 𝐺), 𝜕𝑖}𝑖∈N, thus we adopt the notation

𝐻𝑖 (𝐾; 𝐺). Furthermore, our abstract simplicial complex is obtained from the Rips complex of

discrete data in a metric space (𝑋, 𝑑), i.e. 𝐾 = 𝑅𝜖 (𝑋, 𝑑) and we work with coefficient in a field,

i.e. 𝐺 = F. We denote the resulting i-th homology in this case as 𝐻𝑖 (𝑅𝜖 (𝑋, 𝑑); F). We note the

result is a vector space whose basis elements represent 𝑖-th dimensional holes in the Rips complex

𝑅𝜖 (𝑋, 𝑑).

Persistent homology is obtained from the computation of homology for all 𝜖 ∈ [0, ∞). Since

𝑅𝜖 (𝑋, 𝑑) ⊂ 𝑅𝜖 ′ (𝑋, 𝑑) for 𝜖 ≤ 𝜖 ′, we can obtain linear maps

𝑇𝜖,𝜖 ′ : 𝐻𝑖 (𝑅𝜖 (𝑋, 𝑑); F) −→ 𝐻𝑖 (𝑅𝜖 ′ (𝑋, 𝑑); F)

between vector spaces. Thus, we can keep track of the number of 𝑖-th dimensional holes as 𝜖

changes. This corresponds to changes in the topology of the Rips complex. For example, three

connected components can become one, or a hollow closed surface can be filled. We detail this

notion.

13

Definition 2.1.9. A filtered complex K is a collection of simplicial complex {𝐾𝜖 }𝜖 ≥0 such that

𝐾𝜖 ⊂ 𝐾𝜖 ′

when 𝜖 ≤ 𝜖 ′.

We refer to K as a filtration.

Definition 2.1.10. Let

𝐻𝑖 (K; F) = (cid:8) 𝑇𝜖,𝜖 ′ : 𝐻𝑖 (𝐾𝜖 ; F) −→ 𝐻𝑖 (𝐾𝜖 ′; F), 𝜖 ≤ 𝜖 ′ (cid:9)

denote the family of F-vector spaces and linear transformations 𝑇𝜖,𝜖 ′ induced by the inclusion maps

𝐾𝜖 ↩→ 𝐾𝜖 ′, for 𝜖 ≤ 𝜖 ′. The 𝑖-th persistent homology groups are

𝐻𝜖,𝜖 ′
𝑖

(K; F) := 𝐼𝑚𝑔(𝑇𝜖,𝜖 ′)

and their dimension over F are the persistent Betti numbers

𝛽𝜖,𝜖 ′
𝑖

(K) := 𝑟𝑎𝑛𝑘 (𝑇𝜖,𝜖 ′) = 𝑑𝑖𝑚F (cid:0)𝐻𝜖,𝜖 ′

𝑖

(K; F)(cid:1).

For our purposes, 𝐻𝑖 (K; F) will satisfy the pointwise-finite condition, i.e. 𝛽𝜖,𝜖
𝑖

(K) < ∞ for

every 𝜖 . This will imply that the isomorphism type of 𝐻𝑖 (K; F) is uniquely determined by a multiset

of intervals 𝐼 ⊂ [0, ∞] [25]. These intervals are called the barcode of 𝐻𝑖 (K; F) and are denoted

bcd𝑖 (K). We use the notation bcd𝑅

𝑖 (𝑋, 𝑑) to highlight we are working with the Rips complex of

a metric space. A bar in the barcode, with endpoints 𝑎, 𝑏 ∈ R≥0 and 𝑎 < 𝑏, indicates when a 𝑖-th

dimensional hole first appears and disappears. Indeed, it represents a new basis element at 𝜖 = 𝑎

(birth) and that this basis element will persist until it is lost at 𝜖 = 𝑏 (death). We refer to 𝑏 − 𝑎 as

the lifetime of a basis element.

Plotting bcd𝑅

𝑖 (𝑋, 𝑑) in the plane, now as (𝑎, 𝑏) ∈ R2, gives the persistence diagram. On it,

points that are far away from the line 𝑦 = 𝑥 represent topological features in data that are more
reliable. Namely, the ones having longer lifetime. We denote them analogously as dgm𝑅

𝑖 (𝑋, 𝑑).

The space of persistence diagrams can be given a metric, the bottleneck distance 𝑑𝐵 [26]. In the

following definition we consider the points in the diagonal 𝑦 = 𝑥 as part of the persistence diagram.

14

Definition 2.1.11. Let K1 and K2 be two filtrations. We define the bottleneck distance between

their 𝑖-th persistence diagram as

𝑑𝐵 (𝑑𝑔𝑚𝑖 (K1), 𝑑𝑔𝑚𝑖 (K2)) =

inf
𝜙:𝑑𝑔𝑚𝑖 (K1)→𝑑𝑔𝑚𝑖 (K2)

{

sup
𝑦∈𝑑𝑔𝑚𝑖 (K1)

{∥𝑦 − 𝜙(𝑦) ∥∞}},

where 𝜙 is a bijection of multisets.

Conceptually, given two persistence diagrams we consider a map that pairs points between

them. This forms a pairing to which we can assign a value using the infinity norm of each pair.

We then look at all the combinations of pairings and find the smallest value possible. In essence,

this captures how similar two persistence diagrams are to each other in the plane. We note that

including the infinite points on the diagonal in each persistence diagram allows for a pairing every

time.

When working with metric spaces, we can connect the bottleneck distance with the Gromov-

Hausdorff distance. The latter measures the similarity between bounded metric spaces [5].

Definition 2.1.12. Given two non-empty subsets 𝐴 and 𝐵 of a metric space (𝑋, 𝑑), the Hausdorff

distance 𝑑𝐻 ( 𝐴, 𝐵) is defined as:

𝑑𝐻 ( 𝐴, 𝐵) = max

(cid:26)

sup
𝑎∈𝐴

inf
𝑏∈𝐵

𝑑 (𝑎, 𝑏), sup
𝑏∈𝐵

inf
𝑎∈𝐴

𝑑 (𝑏, 𝑎)

(cid:27)

,

where 𝑑 (𝑎, 𝑏) denotes the distance between points 𝑎 and 𝑏 in 𝑋.

Definition 2.1.13. For two metric spaces (𝑋1, 𝑑1) and (𝑋2, 𝑑2), the Gromov-Hausdorff distance

𝑑𝐺𝐻 (𝑋1, 𝑋2) is defined as the infimum of the Hausdorff distances between the images of 𝑋1 and 𝑋2

in any common metric space 𝑍, over all isometric embeddings 𝜙1 : 𝑋1 → 𝑍 and 𝜙2 : 𝑋2 → 𝑍:

𝑑𝐺𝐻 (𝑋1, 𝑋2) = inf
𝑍,𝜙1,𝜙2

𝑑𝐻 (𝜙1(𝑋1), 𝜙2(𝑋2)).

The connection is illustrated in the following result.

Theorem 2.1.14. Let 𝑋1 and 𝑋2 be totally bounded metric spaces. Then

𝑑𝐵 (bcd𝑅

𝑖 (𝑋1, 𝑑), bcd𝑅

𝑖 (𝑋2, 𝑑)) ≤ 2𝑑𝐺𝐻 (𝑋1, 𝑋2).

15

We conclude this section by presenting a result that applies to the barcodes of a cross product

space.

Definition 2.1.15. Let (𝑋, 𝑑𝑋) and (𝑌 , 𝑑𝑌 ) be metric spaces. The maximum metric 𝑑𝑀 is given by

𝑑𝑀 (cid:0)(𝑥, 𝑦), (𝑥′, 𝑦′)(cid:1) := max{𝑑𝑋 (𝑥, 𝑥′), 𝑑𝑌 (𝑦, 𝑦′)},

where (𝑥, 𝑦), (𝑥′, 𝑦′) ∈ 𝑋 × 𝑌 .

We note (𝑋 × 𝑌 , 𝑑𝑀) is a metric space. Furthermore, its barcodes are given by [14]:

Theorem 2.1.16 (Persistent Künneth Formula ). Let (𝑋, 𝑑𝑋) and (𝑌 , 𝑑𝑌 ) be metric spaces. Then,

(cid:8)𝐼 ∩ 𝐽 | 𝐼 ∈ 𝑏𝑐𝑑 𝑅

𝑖 (𝑋, 𝑑𝑋), 𝐽 ∈ 𝑏𝑐𝑑 𝑅

𝑗 (𝑌 , 𝑑𝑌 )(cid:9),

𝑏𝑐𝑑 𝑅

𝑛 (𝑋 × 𝑌 , 𝑑𝑀) =

(cid:216)

𝑖+ 𝑗=𝑛

for all 𝑛 ∈ N.

2.2 Dynamical Systems

Now that we presented persistent homology as the tool that allow us to assign shape to discrete

data, we now present the framework that models the source of data. Indeed, real world scientific

observations depict an underlying mechanism. Although this mechanism may be unknown to us or

of immense complexity, the theory of dynamical systems enable us to describe it mathematically.

The framework consists of a phase space 𝑀 that represents all of the relevant states of the system

and a function Φ that keeps track of the evolution of the system. Definitions were taken from [2].

Definition 2.2.1. A global continuous time dynamical system is a pair (𝑀, Φ), where 𝑀 is a

topological space and Φ : R × 𝑀 −→ 𝑀 is a continuous map so that Φ(0, 𝑝) = 𝑝, and

Φ(𝑠, Φ(𝑡, 𝑝)) = Φ(𝑠 + 𝑡, 𝑝) for all 𝑝 ∈ 𝑀 and all 𝑡, 𝑠 ∈ R.

We note that a system of ordinary differential equations can give rise to a dynamical system

in which 𝑀 is a smooth manifold. In this case, the dynamics Φ are given by the integral curves

obtained from the system of equations. One can then see that systems modeled by differential

equations are dynamical systems. The evolution of the system is completely determined by the

constrains captured by the equations and the starting state, i.e. the initial condition.

16

Example 2.2.2. Consider the motion of a pendulum without damping or external driving forces.

In this case, the equation governing motion is given by

𝑑2𝜃
𝑑2𝑡

+

𝑔
𝐿 sin(𝜃) = 0,

where 𝜃 is the angle made from the downward vertical, 𝑔 is the acceleration due to gravity, and

L is the length of the pendulum. One can show the phase space 𝑀 is two dimensional depending

only on 𝜃 and (cid:164)𝜃. This means that if we known their value at any point in time, we can then fully

describe the evolution of the system. We illustrate different trajectories in the 𝜃 (cid:164)𝜃-plane representing

the dynamics obtained when placing the hanging mass 𝜃𝑜 from the vertical and then letting it fall,

see Figure 2.1.

Figure 2.1 The phase space of the simple pendulum lies in the 𝜃 (cid:164)𝜃-plane (Example 2.2.2). We
illustrate different trajectories with initial condition (𝜃𝑜, 0). The motion is well known to be
periodic, this is captured by the shape of circular trajectories in 𝑀.

Example 2.2.3. The Wilson-Cowan equations were introduced in 1972 [1]. They have found

great success accounting for behavior observed in neurons [27, 28]. The phase space 𝑀 is two

dimensional. It lies in the 𝐼 𝐸-plane, where 𝐼 corresponds to spiking inhibitory neurons and 𝐸 for

spiking excitatory neurons. The equations involve a sigmoid function S, we present the one in [1],

However, as noted in their work, any sigmoid function is valid.

𝜏1

𝑑𝐸
𝑑𝑡

= −𝐸 + (𝑘1 − 𝑟1𝐸)S1(𝑐1𝐸 − 𝑐2𝐼 + 𝑃(𝑡)),

17

where

𝜏2

𝑑𝐼
𝑑𝑡

= −𝐼 + (𝑘2 − 𝑟2𝐼)S2(𝑐3𝐸 − 𝑐4𝐼 + 𝑄(𝑡)),

S𝑖 (𝑥) =

1
1 + 𝑒−𝑎𝑖 (𝑥−𝜃𝑖)

−

1
1 + 𝑒𝑎𝑖 𝜃𝑖

,

and 𝜏𝑖, 𝑘𝑖, 𝑟𝑖, 𝑐𝑖 are parameters. 𝑃(𝑡) accounts for external input to the excitatory subpopulation and

𝑄(𝑡) external input to the inhibitory subpopulation of neurons. The equations are able to recreate

periodic behavior observed in [27, 28]. We plot trajectories with different initial conditions in this

state, see Figure 2.2. We also plot the solution corresponding to the red trajectory.

Figure 2.2 Dynamical system resulting from the Wilson-Cowan equations (Example 2.2.3). Left:
The phase space 𝑀 is illustrated in the 𝐼 𝐸-plane. We note how different trajectories tend to the red
trajectory. Right: We plot the solution to the red trajectory. It clearly exhibits periodic behavior.

Representing a model in this way allows for a topological understanding of the system. This is

done by looking at the trajectory in the phase space 𝑀. Although most systems can only be solved

numerically, a general understanding of the shape formed in 𝑀 provides qualitative information

of the system [4]. In particular, knowing the topology of an attractor is of great significance. An

attractor is a subset of 𝑀 that pulls nearby trajectories into it.

Definition 2.2.4. A set 𝐴 ⊂ 𝑀 is called an attractor if

1. it is compact,

2. it is an invariant set, i.e. if 𝑎 ∈ 𝐴 then Φ(𝑡, 𝑎) ∈ 𝐴 for all 𝑡 ≥ 0,
3. there is an invariant open neighborhood 𝑈 of 𝐴, so that 𝐴 = (cid:209)𝑡≥0{ Φ(𝑡, 𝑝) : 𝑝 ∈ 𝑈 }.

18

Example 2.2.5. Consider the dynamical system given by the radial equation

(cid:164)𝑟 = 𝑟 (1 − 𝑟 2),

(cid:164)𝜃 = 1

where 𝑟 ≥ 0. This system has a circle as an attractor in the 𝑥𝑦-plane, shown in red in Figure 2.3.

Figure 2.3 We depict trajectories of Example 2.2.5 in the 𝑥𝑦-plane. The attractor of the system is
shown in red.

Example 2.2.6. Consider the van der Pol equation

(cid:164)(cid:164)𝑥 + 𝜇(𝑥2 − 1) (cid:164)𝑥 + 𝑥 = 0

where 𝜇 ≥ 0 is a parameter. This system also has an attractor in the 𝑥 (cid:164)𝑥-plane as shown in Figure

2.4. In this case, it is not a round unit circle, yet it is topologically equivalent to it.

An attractor that is homeomorphic to a circle correspond to a system exhibiting periodicity.

In general, an attractor that is homeomporphic to an 𝑁-dimensional torus, a toroidal attractor,

comes from a quasiperiodic system [5]. Quasiperiodicity has been widely documented in the

scientific literature such as in fMRI of mice, electrocardiograms, crystal formation, and mammalian

vocalization [10, 11, 8, 9]. Its presence is abundant in the real world. Moreover, it has also been

shown to be a transition step between a chaos and stability [12, 13]. Mathematically:

Definition 2.2.7. Let {𝜔𝑖}𝑁

𝑖=1 be positive incommensurate real numbers, i.e. 𝜔𝑖 ∈ R>0 and they
are linearly independent over Q. We say that 𝑓 : R → C is quasiperiodic with frequency vector

19

Figure 2.4 Trajectories of the van der Pol system in the 𝑥 (cid:164)𝑥-plane. The attractor of the system is
shown in red.

𝜔 = (𝜔1, . . . , 𝜔𝑁 ), if

𝑓 (𝑡) = 𝐹 (𝜔1𝑡, . . . , 𝜔𝑁𝑡),

where 𝐹 : T𝑁 → C is a complex-valued continuous function on the 𝑁-torus T𝑁 = R𝑁 /Z𝑁 , called

the parent function of 𝑓 .

In the case of a single frequency, 𝑓 is just a periodic function. Thus, quasiperiodicity generalizes

this notion.

It is a recurrent behavior that can be detected in a time series, see Figure 2.5.

In

practice, one only has an observation of the dynamical system of interest in the form of a time

series. The underlying equations governing the system are unknown i.e. we don’t know 𝑀 and

in turn are unclear about Φ. Nevertheless, there is a way of reconstructing this unknown system

while preserving qualitative information. Concretely, the remarkable 1981 result of Floris Takens

assures us that a time series observation can be enough to reconstruct the topology of the original

attractors [3].

Theorem 2.2.8 (Taken’s Embedding ). Let M be a m-dimensional compact Riemannian manifold.

For pairs (𝜙, ¯𝑓 ), where 𝜙 ∈ 𝐶2(𝑀, 𝑀) and ¯𝑓 ∈ 𝐶2(𝑀, R), it is a generic property that the map

Φ𝜙, ¯𝑓 : 𝑀 → R2𝑀+1 defined as

Φ𝜙, ¯𝑓 ( 𝑝) = (cid:0) ¯𝑓 ( 𝑝), ¯𝑓 (𝜙( 𝑝)), ¯𝑓 (𝜙2( 𝑝)), . . . , ¯𝑓 (𝜙2𝑚 ( 𝑝))(cid:1)

is an embedding.

20

Figure 2.5 Top row: We illustrate a periodic signal on the left and show its sliding window
embedding point cloud in three dimensions. Clearly, it lies on a circle. Bottom row: In this case
we have a quasiperiodic signal generated with two incommensurate frequencies. This time, the
embedding lies on a 2-dimensional torus.

Thus, our reconstruction may well look different than the original but we are assured it will be

topologically equivalent. This is a powerful guarantee that validates the Sliding Window embedding

technique (SW). This is a pipeline enabling the reconstruction of dynamical systems from a time

series. In particular, it has been successfully used to detect quasiperiodicity [5, 6, 7].

Definition 2.2.9. For a function 𝑓 : R → C, an integer 𝑑 > 0 called the embedding dimension,

and a real number 𝜏 > 0 called the time delay, the sliding window embedding of 𝑓 at 𝑡 is given by:

𝑆𝑊𝑑,𝜏 𝑓 (𝑡) =















𝑓 (𝑡)

𝑓 (𝑡 + 𝜏)
...

𝑓 (𝑡 + 𝑑𝜏)















∈ C𝑑+1.

For 𝑇 ∈ N, we denote {𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇

𝑡=0 as the sliding window point cloud. It follows by Takens’

Theorem that for adequate choice of embedding dimension the sliding window point cloud will be

21

an embedding of 𝑀. Indeed, we interpret the time series 𝑓 as the composition

𝑓 (𝑡) = ¯𝑓 (𝜙(0, 𝑝𝑡)),

where ¯𝑓 : 𝑀 → C is an observation function, 𝜙 : R × 𝑀 → 𝑀 is the flow associated to Φ, and

𝑝𝑡 = 𝜙(𝑡, 𝑝0) for some 𝑝0 ∈ 𝑀. Noting

𝑓 (𝑡 + 𝑖𝜏) = ¯𝑓 (𝜙(𝑡, 𝜙(𝑖𝜏, 𝑝0))) = ¯𝑓 (𝜙(𝑡, 𝜙𝑖 (𝜏, 𝑝0))),

for 𝑖 ∈ N, justifies the claim.

Example 2.2.10. Consider the Lorenz equations

(cid:164)𝑥 = 𝜎(𝑦 − 𝑥)

(cid:164)𝑦 = 𝑟𝑥 − 𝑦 − 𝑥𝑧

(cid:164)𝑧 = 𝑥𝑦 − 𝑏𝑧.

Figure 2.6 Top row: We depict the Lorenz system, Example 2.2.10, with parameter values 𝜎 =
10, 𝑟 = 28, and 𝑏 = 8/3. We note 𝑀 is a bounded subset of R3. On the right we show the solution
of the system for a particular initial condition. Bottom row: We use the solution of the system
to reconstruct the phase space. Namely, on the left we have the result when using the embedding
𝑆𝑊𝑑,𝜏𝑥(𝑡) and on the right when using 𝑆𝑊𝑑,𝜏 𝑦(𝑡).

22

This dynamical system exhibits fractal behavior in three dimensional space when 𝜎 = 10, 𝑏 =

8/3, 𝑟 = 28. We reconstruct this behavior using the sliding window embedding map. We illustrate

two different reconstructions by using 𝑥(𝑡) and 𝑦(𝑡), see Figure 2.6.

For the reconstruction of quasiperiodic systems, an adequate choice of parameters 𝑑 and 𝜏 is

established in [5]. We are thus in the position of probing for recurrent behavior in a dynamical

system by analyzing the sliding window point cloud of observed time series. The tool capable for

the task is persistent homology as it provides a convenient geometrical representation at a multi-

scale level. This leads to a robust approach capable of discerning noise and synthetic artifacts.

In particular, for toroidal attractors, we are looking for homological features of an 𝑁-dimensional

torus, as depicted in Figure 2.7. These steps constitute SW:

Figure 2.7 Top row: We perform SW to the radial system. The reconstruction is done using the 𝑥(𝑡)
solution belonging to the attractor. Bottom row: SW to the van der Pol system. The reconstruction
is also done using 𝑥(𝑡) of the trajectory corresponding to the attractor.

1. Start with a time series pertaining to a measurement of a system,

2. Construct the sliding window point cloud from it,

23

3. Compute their persistence diagrams to discern for the presence of an 𝑁-torus.

Example 2.2.11. We illustrate SW by considering the radial equation and the van der Pol equation

illustrated in Figure 2.3 and Figure 2.4, respectively. They both have an attractor that is topologi-

cally equivalent to the circle which translates to periodic behavior. Let us consider the solution of the

attractor. Using 𝑥(𝑡) as input, we perform SW and reconstruct the attractor in three dimensions. As

illustrated in Figure 2.7, both persistence diagrams show a persistent 1-dimensional feature and no

significant 2-dimensional features. This is corresponds to the homology characterizing topologies

equivalent to the circle.

Although SW has been well established and found multiple applications [6, 7, 14], computing

persistence diagrams is computationally taxing. This motivated our work of developing a faster

alternative. We move on to cover the tools that helped us achieve this effort.

2.3 Fourier Series

For our purposes, we leverage the fast fourier transform (FFT) of time series data to infer the

frequency vector. This is validated by theoretical results we present in this section. They can be

found in [5].

Definition 2.3.1. Let T be the quotient space R/(2𝜋Z). For a single variable periodic function

𝑔 ∈ 𝐿2(T) and 𝑍 ∈ N, its Fourier series and its 𝑍-truncated Fourier polynomial are written as

𝑔(𝑡) =

∞
∑︁

𝑘=−∞

ˆ𝑔(𝑘)𝑒𝑖𝑘𝑡

and

𝑆𝑍 𝑔(𝑡) =

𝑍
∑︁

𝑘=−𝑍

ˆ𝑔(𝑘)𝑒𝑖𝑘𝑡

where the Fourier coefficients of 𝑔 are

ˆ𝑔(𝑘) =

1
2𝜋

∫ 𝜋

−𝜋

𝑔(𝑘)𝑒−𝑖𝑘𝑡 .

If 𝑔 is continuously differentiable, one can show pointwise converge [29]. Let us now introduce

the 𝑁-dimensional Fourier series.

Definition 2.3.2. For 𝐾 ∈ N, let

𝐾 = {k ∈ Z𝑁 | ∥k∥∞ ≤ 𝐾 },
𝐼 𝑁

24

where ∥k∥∞ = max

𝑖

|𝑘𝑖 |. The 𝐾-truncated Fourier polynomial of 𝐹 ∈ 𝐿2(T𝑁 ) is the function

𝑆𝐾 𝐹 (t) =

ˆ𝐹 (k)𝑒𝑖⟨k,t⟩

∑︁

k∈𝐼 𝑁
𝐾

where t ∈ R𝑁 , ⟨·, ·⟩ is the standard inner product in R𝑁 , and

ˆ𝐹 (k) =

1
(2𝜋) 𝑁

∫ 2𝜋

∫ 2𝜋

· · ·

0

0

𝐹 (t)𝑒−𝑖⟨k,t⟩𝑑𝑡1 · · · 𝑑𝑡𝑁 = (cid:10)𝐹, 𝑒𝑖⟨k,·⟩(cid:11)

𝐿2

is the k-Fourier coefficient of 𝐹, for k ∈ Z𝑁 .

Analogously to the one dimensional case, the sequence {𝑆𝐾 𝐹}𝐾∈N coverges to 𝐹 in 𝐿2(T𝑁 )

as 𝐾 → ∞. This analytical tool has been used to show the persistence diagrams obtained from

a quasiperiodic function can be approximated in bottleneck distance with a truncated sum of

exponentials. This result is paramount to the treatment of general quasiperiodic functions.

Theorem 2.3.3. [5] Let 𝑓 : R → C be a quasiperiodic function with frequency vector 𝜔 and parent

function 𝐹. For k ∈ Z𝑁 and 𝐾 ∈ N

ˆ𝐹 (k) = lim
𝜆→∞

1
𝜆

∫ 𝜆

0

𝑓 (𝑡)𝑒−𝑖⟨k,𝑡𝜔⟩𝑑𝑡.

Furthermore, for 𝑗 ∈ N, dgm𝑅

𝑗 (𝑆𝑊𝑑,𝜏 𝑓 ) can be approximated in bottleneck distance by

dgm𝑅

𝑗 (𝑆𝑊𝑑,𝜏𝑆𝐾 𝑓 )

as

𝐾 → ∞,

where

𝑆𝐾 𝑓 (𝑡) =

∑︁

ˆ𝐹 (k)𝑒𝑖⟨k,𝑡𝜔⟩.

The approximation is of order 𝑂 (cid:0)𝐾 𝑁

2 −𝑟 (cid:1) when | ˆ𝐹 (k)| = 𝑂 (∥k∥−𝑟

2 ) and 𝑟 > 𝑁/2.

∥k∥∞≤𝐾

This result allow us to focus on quasiperiodic functions that look like a sum of exponentials

𝑓 (𝑡) =

∑︁

𝑐 𝑗 𝑒𝑖𝜔 𝑗 𝑡,

𝑗

25

where 𝑐 𝑗 ∈ C and 𝜔 𝑗 ∈ R. In this case the sliding window embedding can be expressed as

𝑆𝑊𝑑,𝜏 𝑓 (𝑡) =

√

where ¯𝑐𝑖 =

· · ·
. . .
. . .

1

1

√

1
𝑑 + 1

𝑒𝑖𝜔1𝜏
...

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)
𝑑 + 1 𝑐𝑖, 𝑑 ∈ N, and 𝜏 ∈ R. The vector 𝑣𝑡 lies in a space homeomorphic to the 𝑁-torus,

¯𝑐1𝑒𝑖𝜔1𝑡
...

𝑒𝑖𝜔𝑀 𝜏
...

· · · 𝑒𝑖𝜔𝑀 𝑑𝜏

¯𝑐𝑀 𝑒𝑖𝜔𝑀 𝑡

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

= 𝐴𝑣𝑡,

𝑒𝑖𝜔1𝑑𝜏

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

for some 𝑁 ≤ 𝑀; furthermore 𝑑 and 𝜏 can be chosen for the same to hold for 𝑆𝑊𝑑,𝜏 𝑓 (𝑡). Indeed,

this was shown in [5] where it was concluded that 𝑁 corresponds to the number of incommensurate

frequencies.

The Cartesian product

𝐺𝑇 =

(cid:214)

𝑗

{ ¯𝑐 𝑗 𝑒𝑖𝜔 𝑗 𝑡 }𝑇 ′
𝑡=0

offers a way of approximating {𝑣𝑡 }𝑇

𝑡=0 in bottleneck distance via Theorem 2.1.14. This is useful since
the persistence diagrams of 𝐺𝑇 can be readily computed using the Persistent Künneth Formula.

That this approach can also be used to approximate the sliding window embedding of 𝑓 is presented

in Theorem 4.3.1.

2.4 Continued Fraction Expansion

We now introduce a notion that is at the center of our method. It allow us to infer information

from our frequency parameter and compute an approximation of the target persistence diagram.

The idea is to express a number as a nested sequence of fractions. This can be achieved by a

continued application of the division algorithm. For a rational number, this process will end after

finitely many steps. We are primarily concerned with irrational numbers which have an infinite

number of nested fractions. Definitions and results were retrieved from [18].

Definition 2.4.1. Let 𝜔 be a positive irrational number. The continued fraction of 𝜔 is given by

𝜔 = 𝑎1 +

1

𝑎2 +

1
𝑎3+ 1
𝑎

4+···

:= [𝑎1, 𝑎2, 𝑎3, · · · ],

where 𝑎1 ∈ Z and 𝑎𝑖 ∈ N, for 𝑖 > 1.

26

As we alluded to, the expansion is infinite if and only if 𝜔 is irrational. Let us go over the

iterative process to compute the expansion, namely to get the 𝑎𝑖’s from 𝜔. We denote by ⌊𝑥⌋ the

floor function, i.e. the greatest integer less than x. We start by letting 𝑎1 = ⌊𝜔⌋ and writing

where 𝑥2 is the irrational number given by

𝜔 = 𝑎1 +

,

1
𝑥2

𝑥2 =

1
𝜔 − 𝑎1

.

To calculate 𝑎2 we note

where

𝜔 = 𝑎1 +

1
⌊𝑥2⌋ + 1
𝑥3

,

𝑥3 =

1
𝑥2 − ⌊𝑥2⌋

.

Thus, 𝑎2 = ⌊𝑥2⌋. Noticing the pattern, we can conclude that 𝑎𝑖 = ⌊𝑥𝑖⌋ where

𝑥𝑖 =

1
𝑥𝑖−1 − ⌊𝑥𝑖−1⌋

and 𝑥1 = 𝜔.

Example 2.4.2. Let us consider the case 𝜔 =

√

2. Since

𝑥2 =

√

1
2 − 1

= 2.41421 . . . ,

𝑎2 = ⌊𝑥2⌋ = 2. Similarly, 𝑎3 = 2 since

𝑥3 =

1
2.41421 · · · − 2

= 2.41421 . . . .

In fact, this shows that for all 𝑖, 𝑎𝑖 = 2, i.e.

√

2 = 1 +

1

2 +

1
2+ 1
2+···

= [1, 2, 2, · · · ].

27

Example 2.4.3. The case 𝜔 =

√

2 is an example of a quadratic irrational, namely a number of the

form

√

𝑃 ±
𝑄

𝐷

,

where 𝑃, 𝑄, 𝐷 are integers and 𝐷 is positive and not a perfect square. Lagrange proved that these

irrational numbers have periodic continued fractions, i.e. there exists a 𝑗𝑜 such that 𝑖 > 𝑗𝑜 implies

𝑎𝑖 = 𝑎𝑖+𝑟𝑜

for some fixed 𝑟𝑜 ≥ 1. Indeed, as show for

√

expansions with a bar above the repeating terms. With this notation,

2, this is true for 𝑗𝑜 = 1 and 𝑟𝑜 = 1. We denote periodic
√

2 = [1, ¯2]. Similarly, one

can readily show that

in this case 𝑗𝑜 = 0 and 𝑟𝑜 = 3.

√

10

1 +

3

= [1, 2, 1],

Example 2.4.4. Gauss used a generalized notion of continued fractions to express the series [30]

𝐹 (𝑎, 𝑏, 𝛾, 𝑥) = 1 +

∞
∑︁

𝑖=0

(cid:206)𝑖
(𝑖 + 1)! (cid:206)𝑖

𝑗=0(𝑎 + 𝑗) (𝑏 + 𝑗)
𝑗=0(𝛾 + 𝑗)

𝑥𝑖+1.

Using this result, one can readily obtain closed formulas for the continued fraction expansion of

elementary functions. Indeed, one can show that for the positive integer 𝑛, the rational number

𝑝
𝑞 ,

where 𝑔𝑐𝑑 ( 𝑝, 𝑞) = 1, and 𝐼𝑛 (𝑥) denoting the hyperbolic Bessel function of the first kind,

𝑒1/𝑛 = [1, 𝑛 − 1, 1, 1, 3𝑛 − 1, 1, 1, 5𝑛 − 1, 1, 1, . . . ],

tanh(1/𝑛) = [0, 𝑛, 3𝑛, 5𝑛, 7𝑛, 9𝑛, 11𝑛, . . . ],

and

𝐼𝑝/𝑞 (2/𝑞)
𝐼1+𝑝/𝑞 (2/𝑞)

= [ 𝑝 + 𝑞, 𝑝 + 2𝑞, 𝑝 + 3𝑞, 𝑝 + 4𝑞, . . . ].

From here onwards, we let 𝜔 denote a positive irrational number.

Definition 2.4.5. Let 𝑖 ∈ N. The 𝑖-th convergent of 𝜔 is given by

𝑝𝑖
𝑞𝑖

= [𝑎1, 𝑎2, · · · , 𝑎𝑖].

28

The terms 𝑝𝑖 and 𝑞𝑖 play a special role in the theory. They can be obtained recursively thanks

to the following result.

Proposition 2.4.6. The numerators 𝑝𝑖 and the denominators 𝑞𝑖 of the 𝑖-th convergent of 𝜔 satisfy

the equations

for 𝑖 ≥ 1, where

𝑝𝑖 = 𝑎𝑖 𝑝𝑖−1 + 𝑝𝑖−2,

𝑞𝑖 = 𝑎𝑖𝑞𝑖−1 + 𝑞𝑖−2.

𝑝0 = 1,

𝑞0 = 0,

𝑝−1 = 0,

𝑞−1 = 1.

This provides a convenient computational avenue for obtaining the 𝑖-th convergent. It can also

be used to show the following result.

Proposition 2.4.7. 𝑝𝑖 and 𝑞𝑖 have no common divisor other than 1 or −1.

We can also assert how good of an approximation they are to 𝜔.

Proposition 2.4.8. Let 𝑖 ≥ 1. Then

1
2𝑞𝑖𝑞𝑖+1

<

(cid:12)
(cid:12)
(cid:12)

𝜔 −

𝑝𝑖
𝑞𝑖

(cid:12)
(cid:12)
(cid:12)

<

1
𝑞𝑖𝑞1+1

< 1
𝑞2
𝑖

.

Furthermore, this result allow us to establish the existence and uniqueness of infinite continued

fraction expansions. Indeed, this follows by noting 𝑞𝑖 < 𝑞𝑖+1 for all 𝑖 ≥ 1. This justifies the equality

sign in the definition we presented. The next result illustrates more of the type of convergence. It

characterizes subsequences converging from below and from above.

Proposition 2.4.9. Let 𝑖 ≥ 1 be odd. Then

𝑝𝑖
𝑞𝑖

< 𝜔 <

𝑝𝑖+1
𝑞𝑖+1

.

29

Example 2.4.10. The results presented here can be used to argue the golden ration 𝜑 is the ‘most’

irrational number. Indeed, by noting

𝜑 = 1 +

1

1 +

1
1+ 1
1+···

= [1, 1, 1, · · · ],

one can show this expression corresponds to the slowest rate of convergence possible. Furthermore,

the Klein diagram allow us to obtain a geometrical picture of the results presented, see Figure 2.8.

Figure 2.8 Klein diagram of the golden ratio.

We detailed how exponential functions play a critical role in the study of quasiperiodic functions.

Let us consider a periodic function of the form

𝑓 (𝑡) = 𝑒𝑖𝜔𝑡 .

One can topologically describe { 𝑓 (𝑡)}𝑇

𝑡=0 as a sampling of the circle. We can consider the analogues

sampling using ‘flat’ coordinates.

Definition 2.4.11. Let 𝑁 ∈ N and [𝑥] = 𝑥 mod 1. We define the set

𝑆𝜔,𝑇 := {[0], [𝜔], . . . , [𝑇𝜔]}.

30

By identifying the endpoints of [0, 1], this set can also be considered a sampling of the circle

as shown in Figure 2.9. Furthermore, the function

Figure 2.9 On the left we depict 𝑆𝜋−1,17. Each different gap size is shown in a different color. We
note in this case there is only two gaps. On the right we show 𝑆√
5,17. This case does have three
different gaps.

¯𝑑 (𝑥, 𝑦) = min{|𝑥 − 𝑦|, |1 − 𝑦 + 𝑥|, |1 − 𝑥 + 𝑦|}

defines a metric on this quotient space aligned with the notion of considering the space as a circle.

It makes (𝑆𝜔,𝑇 , ¯𝑑) a metric space. As we will illustrate in this work, understanding 𝑆𝜔,𝑇 is sufficient

for the treatment of sampling of the form { 𝑓 (𝑡)}𝑇

𝑡=0. The set 𝑆𝜔,𝑇 is known in the literature as the
Kronecker sequence. Its asymptotic behavior has been well documented [31]. Moreover, there is a

remarkable result about this set: The Thee Gap Theorem, also known as the Steinhaus conjecture,

it states [20]:

Theorem 2.4.12 (Three Gap Theorem). Let 𝑇 ≥ 1. The points in 𝑆𝜔,𝑇 partition [0, 1] into 𝑇 + 1

intervals, such that their lengths take at most 3 different values 𝛿 𝐴, 𝛿𝐵 and 𝛿𝐶, with 𝛿𝐶 = 𝛿 𝐴 + 𝛿𝐵.

We note the result is independent of 𝜔 or 𝑇. Furthermore, the gaps from the theorem are closely

related to the convergents of 𝜔 [32]. This forms the basis of our method, The Three Gap Theorem

Method (3G), which consists in the application of the results of Theorem 4.2.1 and Theorem 4.3.2.

We depict it schematically in Figure 2.10. It uses code, Algorithm 2.1, validated by The Three

Gap Theorem to compute persistence diagrams that we then put back together using the Persistent

31

Figure 2.10 Schematic of our method 3G. We start with a sum of exponentials approximating a
general quasiperiodic functions. The frequency parameters are retrieved using FFT. We apply
our results to each frequency independently. We then put back the result using the Persistent
Künneth formula in persistence homology. The result approximates a cross product space that is
topologically similar to the sliding window point cloud of interest.

Algorithm 2.1 3 Gap Code
1: Input frequency 𝜔𝑖
2: Scale frequency 𝜔𝑖 = 𝜔𝑖/(2𝜋)
3: Obtain continued fraction expansion, C.F.E., of 𝜔𝑖
4: Use C.F.E. as shown in Theorem 4.2.1 to obtain gaps & their multiplicity
5: Scale gaps as detailed in Theorem 4.3.2 [𝑎, 𝑏) → ¯𝑐𝑖 [ ¯𝑎, ¯𝑏)
6: Output gaps

Künneth Formula. The resulting diagrams pertain to 𝐺𝑇 which can approximate the sliding window

point cloud of the quasiperiodic function of interest, see Theorem 4.3.1.

32

CHAPTER 3

STABILITY OF THE THREE GAP THEOREM IN PERSISTENCE DIAGRAMS

We begin by establishing the stability of The Three Gap Theorem in persistence diagrams. Our

results validate the usage of the theorem in a computational setting. As a result we will establish the

The Three Gap Theorem Method (3G) in the next chapter. Our stability results consist in looking

at the number of matching continued fraction terms of two frequency values and the bottleneck

distance of the corresponding sets. Furthermore, we present this stability in terms of the number

of matching decimals of the two parameter values. Our results quantify the error one could expect

when working with numerical approximations in this context.

3.1 Preliminaries

When a dynamical system transitions from a stable to a chaotic state, it undergoes an interme-

diary transition that exhibits repetitive like motions: quasiperiodicity [12]. The Sliding Window

Embedding (SW) method has been successfully used for quasiperiodicity detection [5]. It consists

in reconstructing a time series using the sliding window map and then computing the persistence

homology of the reconstruction. Since the reconstruction of a quasiperiodic signal with 𝑁 incom-

mensurable frequencies will span a 𝑁 dimensional torus, the geometric characterisation provided

by looking at the persistence diagram is of most value. However, obtaining this information is

computationally taxing and is thus limited to small sample sizes. This was the motivation for an

approximation alternative: The Three Gap Theorem Method (3G). This avenue provides a much

quicker and computationally feasible approach allowing for scaling of the sample size [33]. To

achieve this, we leverage the continued fraction expansion of the frequency parameters. This is

possible due to the Three Gap Theorem [20]. The theorem is depicted in Figure 3.1. Being the

basis of our approximation method, we begin by validating its stability. As shown in Figure 3.1, two

different sampling parameters close in norm can have very different outcomes. Our work addresses

this issue by providing a better stability condition, namely decimal precision. Indeed, denoting
𝑥 𝑗 = ⌊10 𝑗 𝜔⌋
10 𝑗

, the first 𝑗 decimal terms of 𝜔, we prove

33

Theorem 3.2.9 Suppose

𝑥 𝑗 = ¯𝑥 𝑗

for some 𝑗 ∈ N. If 0 < 𝑇 < 𝑞𝑘 𝑗 (𝜔)−1, then for 𝑖 ≥ 0

𝑑𝐵 (𝑑𝑔𝑚 𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) = 𝑂 (𝑘 𝑗 (𝜔)−1)

as 𝑗 → ∞.

The function 𝑘 𝑗 (𝜔) is a non-decreasing function that allow us to connect continued fraction

expansion precision with decimal precision. Furthermore, we can also state the rate of decay of the

bottleneck distance of two persistence diagrams as we increase decimal precision

Theorem 3.2.10 For 0 < 𝜖 < 𝑧0 there exists positive constants 𝐶, 𝜆 (depending on 𝜖) with 0 < 𝜆 < 1

such that for all integers 𝑗 ≥ 1 for which

𝑥 𝑗 = ¯𝑥 𝑗 ,

and 𝑇 < 𝑞𝑘 𝑗 (𝜔)−1, if 𝑘 𝑗 (𝜔) ≥ 𝑗 𝑧0, then for 𝑖 ≥ 0

(cid:16)

𝑃

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

(cid:17)

1
𝑗 (𝜖 + 𝑧0)

≤ 𝐶𝜆 𝑗 ,

otherwise

(cid:16)

𝑃

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

(cid:17)

1
𝑗 (𝑧0 − 𝜖)

≥ 1 − 𝐶𝜆 𝑗 .

This exponential decay is exhibited in Figure 3.3 and Figure 3.4. It is a valuable guaranteed

since 3G entails using the fast Fourier transform to retrieve frequencies. We can now proceed

knowing the error we could expect given a sample size.

3.2 Stability Results

3.2.1 Stability with respect to C.F.E.

We begin by working with the continued fraction expansion of two parameter values. For what

follows, we let 𝜔 and ¯𝜔 denote two positive irrational numbers. We denote their 𝑖-th continued

fraction terms and 𝑗-th convergents by 𝑎𝑖, ¯𝑎𝑖 and 𝑝 𝑗 /𝑞 𝑗 , ¯𝑝 𝑗 / ¯𝑞 𝑗 , respectively.

34

Figure 3.1 We illustrate the result of the three gap theorem and how they translates directly to
results in persistence homology. Top row: The set 𝑆𝜋−1,17 is illustrated with the two gaps it creates.
Bottow row: The set 𝑆√
5,17 creates three gaps. Due to clustering, it is a poor covering of the circle.
We note 𝑑2(𝜋 − 1,

5) = 0.094.

√

Proposition 3.2.1. Suppose there exists a 𝑗 ∈ N such that 𝑎𝑖 = ¯𝑎𝑖 for all 𝑖 ≤ 𝑗. If 0 < 𝑇 < 𝑞 𝑗−1,

𝑇 |𝜔 − ¯𝜔| < 1
𝑞 𝑗

.

Proof. W.L.O.G we assume 𝜔 < ¯𝜔. If 𝑗 − 1 is even,

𝑝 𝑗 −1
𝑞 𝑗 −1

≤ 𝜔 < ¯𝜔. Since | ¯𝜔 − ¯𝑝 𝑗 −1
¯𝑞 𝑗 −1

| <

1

¯𝑞 𝑗 −1 ¯𝑞 𝑗 and

noting 𝑞𝑖 = ¯𝑞𝑖, 𝑝𝑖 = ¯𝑝𝑖 for all 𝑖 ≤ 𝑗, we conclude

|𝜔 − ¯𝜔| = ¯𝜔 − 𝜔 < ¯𝜔 −

¯𝑝 𝑗−1
¯𝑞 𝑗−1

<

1
𝑞 𝑗−1𝑞 𝑗

.

Similarly, if 𝑗 − 1 is odd we note 𝜔 < ¯𝜔 ≤ ¯𝑝 𝑗 −1
¯𝑞 𝑗 −1

=

𝑝 𝑗 −1
𝑞 𝑗 −1

. Using |𝜔 −

𝑝 𝑗 −1
𝑞 𝑗 −1

| <

1

𝑞 𝑗 −1𝑞 𝑗 , we conclude

Thus,

|𝜔 − ¯𝜔| <

𝑝 𝑗−1
𝑞 𝑗−1

− 𝜔 <

1
𝑞 𝑗−1𝑞 𝑗

.

𝑇 |𝜔 − ¯𝜔| <

𝑞 𝑗−1
𝑞 𝑗 𝑞 𝑗−1

=

.

1
𝑞 𝑗

35

□

Corollary 3.2.2. Suppose there exists a 𝑗 ∈ N such that 𝑎𝑖 = ¯𝑎𝑖 for all 𝑖 ≤ 𝑗. If 0 < 𝑇 < 𝑞 𝑗−1,

¯𝑑 ( [𝑟𝜔], [𝑟 ¯𝜔]) < 1
𝑞 𝑗

where 𝑟 ∈ N and 𝑟 ≤ 𝑇.

Proof. By the previous lemma, we note

𝑇 |𝜔 − ¯𝜔| < 1
𝑞 𝑗

≤

.

1
2

This implies that for 𝑟 ∈ N and 𝑟 ≤ 𝑇

¯𝑑 ([𝑟𝜔], [𝑟 ¯𝜔]) = 𝑟 |𝜔 − ¯𝜔|.

The result follows.

□

Corollary 3.2.2 provides us with a matching of 𝑆𝜔,𝑇 and 𝑆 ¯𝜔,𝑇 such that pairs have a distance

less than 1/𝑞 𝑗 . This translates to an upper bound of 𝑑𝐻 (𝑆𝜔,𝑇 , 𝑆 ¯𝜔,𝑇 ). Moreover, by translating one

of the sets by a factor of 𝑇 |𝜔 − ¯𝜔|/2 we obtain an analogues bound for 𝑑𝐺𝐻 (𝑆𝜔,𝑇 , 𝑆 ¯𝜔,𝑇 ):

Proposition 3.2.3. Suppose there exists a 𝑗 ∈ N such that 𝜔𝑖 = ¯𝜔𝑖 for all 𝑖 ≤ 𝑗. If 0 < 𝑇 < 𝑞𝑘−1,

then

𝑑𝐺𝐻 (𝑆𝜔,𝑇 , 𝑆 ¯𝜔,𝑇 ) ≤

1
2𝑞 𝑗

.

Proof. W.L.O.G. we assume 𝜔 < ¯𝜔. Consider the function 𝑓 : [0, 1) → [0, 1) given by

𝑓 (𝑥) =

(cid:104)

𝑥 −

𝑇 |𝜔 − ¯𝜔|
2

(cid:105)

.

Clearly, 𝑓 is an isometry. Furthermore, since

𝑇 |𝜔− ¯𝜔|
2

< 1

2, we note that for 𝑟 ∈ N and 𝑟 ≤ 𝑇

¯𝑑 ( 𝑓 ([𝑟 ¯𝜔]), [𝑟𝜔]) =

(cid:12)
(cid:12)
(cid:12)

¯𝑑 ( [𝑟 ¯𝜔], [𝑟𝜔]) −

𝑇 |𝜔 − ¯𝜔|
2

(cid:12)
(cid:12)
(cid:12)

≤

𝑇 |𝜔 − ¯𝜔|
2

.

Thus,

This shows the result.

𝑑𝐻 ( 𝑓 (𝑆 ¯𝜔,𝑇 ), 𝑆𝜔,𝑇 ) ≤

𝑇 |𝜔 − ¯𝜔|
2

< 1
2𝑞 𝑗

.

□

36

We can then apply Theorem 2.1.14 to conclude the next result.

Corollary 3.2.4. Suppose there exists a 𝑗 ∈ N such that 𝜔𝑟 = ¯𝜔𝑟 for all 𝑟 ≤ 𝑗. If 0 < 𝑇 < 𝑞𝑘−1,

then for 𝑖 ≥ 0

𝑑𝐵 (𝑑𝑔𝑚 𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

1
𝑞𝑘

.

Our results up to now have been with respect to continued fraction expansions. We are able

able to relate them to decimal precision using the function 𝑘𝑖 (𝜔) which takes the first 𝑖-th terms in

the decimal expansion of 𝜔 and returns the number of matching partial quotients of the continued

fraction expansion of 𝜔.

3.2.2 Stability with respect to Decimal Precision

We introduce the function 𝑘𝑖 (𝜔) by presenting examples on how to compute it. To do so, we

need two decimals

𝑥𝑖 =

⌊10𝑖𝜔⌋
10𝑖

and

𝑦𝑖 = 𝑥𝑖 +

1
10𝑖

,

where ⌊𝜔⌋ denotes the floor of 𝜔. We obtain 𝑘𝑖 (𝜔) by tracking the number of the first matching

quotients in the continued fraction expansions of 𝑥𝑖 and 𝑦𝑖.
√
For example, consider 𝜔 = 3

7 = 1.91293118.... To compute 𝑘3(𝜔) we note

𝑥3 = 1.912

and

𝑦3 = 1.913,

Their continued fraction expansions are [1, 1, 10, 2, 1, 3] and [1, 1, 10, 2, 43], respectively. Thus,

√
𝑘3( 3

7) = 3. Similarly, for 𝑘6(𝜔) we note

𝑥6 = 1.912931

and

𝑦6 = 1.912932,

Their respective continued fraction expansions are [1, 1, 10, 2, 16, 3, 21, 4, 2] and

√
[1, 1, 10, 2, 16, 2, 11, 2, 1, 2, 3]. Thus, 𝑘6( 3

7) = 4.

For a fixed 𝜔, the function 𝑘𝑖 (𝜔) can always be computed using these steps. Furthermore, this

function has been well studied in the past [34]:

Theorem 3.2.5. For almost all irrationals 𝜔, with respect to Lebesgue measure, we have

𝑘𝑖 (𝜔)
𝑖

=

6 𝑙𝑜𝑔(10) 𝑙𝑜𝑔(2)
𝜋2

.

lim
𝑖→∞

37

We denote the limit above by 𝑧0 and apply a more recent result to our work [35]:

Theorem 3.2.6. For all 𝜖 > 0, there exists positive constants 𝐶, 𝜆 (depending on 𝜖) with 0 < 𝜆 < 1

such that

for all integers 𝑖 ≥ 1.

𝑃

(cid:16)(cid:12)
(cid:12)
(cid:12)

𝑘𝑖 (𝜔)
𝑖

− 𝑧0

(cid:12)
(cid:12)
(cid:12)

(cid:17)

≥ 𝜖

≤ 𝐶𝜆𝑖

Figure 3.2 We illustrate the result of Theorem 3.2.6. We compute the 𝑘𝑖 function for 5,000 values
of 𝜔. These are obtained from sampling the interval (0,1) uniformly.

We now relate the 𝑘 function to our work. The following lemma is obtained by construction.

Let us denote by ¯𝑥 𝑗 the decimal obtained from ¯𝜔.

Lemma 3.2.7. Suppose

𝑥 𝑗 = ¯𝑥 𝑗

for some 𝑗 ∈ N. Then 𝑎𝑟 = ¯𝑎𝑟 for all 𝑟 ≤ 𝑘 𝑗 (𝜔).

We can now combined our previous results to obtain a relation in terms of decimal precision:

Proposition 3.2.8. Suppose

𝑥 𝑗 = ¯𝑥 𝑗

for some 𝑗 ∈ N. If 0 < 𝑇 < 𝑞𝑘 𝑗 (𝜔)−1, then for 𝑖 ≥ 0

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

1
𝑞𝑘 𝑗 (𝜔)

.

38

This provides us with a rate of decay for the bottleneck distance. We further simplify this

observation by noting that 𝑞𝑘 𝑗 (𝜔) ≥ 𝑘 𝑗 (𝜔).

Theorem 3.2.9. Suppose

for some 𝑗 ∈ N. If 0 < 𝑇 < 𝑞𝑘 𝑗 (𝜔)−1, then for 𝑖 ≥ 0

𝑥 𝑗 = ¯𝑥 𝑗

𝑑𝐵 (𝑑𝑔𝑚 𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) = 𝑂 (𝑘 𝑗 (𝜔)−1)

𝑎𝑠 𝑗 → ∞.

We can further obtain a relation illustrating the rate at which the bottleneck distance decays.

Theorem 3.2.10. For 0 < 𝜖 < 𝑧0 there exists positive constants 𝐶, 𝜆 (depending on 𝜖) with

0 < 𝜆 < 1 such that for all integers 𝑗 ≥ 1 for which

𝑥 𝑗 = ¯𝑥 𝑗

and 𝑇 < 𝑞𝑘 𝑗 (𝜔)−1, if 𝑘 𝑗 (𝜔) ≥ 𝑗 𝑧0, then for 𝑖 ≥ 0

(cid:16)

𝑃

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

(cid:17)

1
𝑗 (𝜖 + 𝑧0)

≤ 𝐶𝜆 𝑗 ,

otherwise

(cid:16)

𝑃

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,𝑇 )) ≤

(cid:17)

1
𝑗 (𝑧0 − 𝜖)

≥ 1 − 𝐶𝜆 𝑗 .

Proof. Let 𝜖 > 0. By Theorem 3.2.6 there exists positive constants 𝐶, 𝜆 (depending on 𝜖) with

0 < 𝜆 < 1 such that

If 𝑘 𝑗 (𝜔) ≥ 𝑗 𝑧0 this implies

𝑃

(cid:16)(cid:12)
(cid:12)
(cid:12)

𝑘 𝑗 (𝜔) − 𝑗 𝑧0

(cid:17)

≥ 𝑗𝜖

≤ 𝐶𝜆 𝑗 .

(cid:12)
(cid:12)
(cid:12)

(cid:16)

𝑃

1
𝑘 𝑗 (𝜔)

≤

(cid:17)

1
𝑗𝜖 + 𝑗 𝑧0

≤ 𝐶𝜆 𝑗 ,

combined with Theorem 3.2.8 give the first half of the result.

If 𝑘 𝑗 (𝜔) < 𝑗 𝑧0 we can use the

inequality, following from Theorem 3.2.6 as well,

𝑃

(cid:16)(cid:12)
(cid:12)
(cid:12)

𝑘 𝑗 (𝜔) − 𝑗 𝑧0

(cid:17)

< 𝑗𝜖

(cid:12)
(cid:12)
(cid:12)

≥ 1 − 𝐶𝜆 𝑗 ,

39

to conclude

1
𝑗 𝑧0 − 𝑗𝜖
since 𝑗 𝑧0 − 𝑗𝜖 > 0. Similarly, applying Theorem 3.2.8 shows the second part of the result and

≥ 1 − 𝐶𝜆 𝑗

1
𝑘 𝑗 (𝜔)

≤

𝑃

(cid:16)

(cid:17)

completes the proof.

3.3 Applications

□

We now analyze the error in computations using our main results. We detail their applica-

tion when retrieving frequencies using the Fast Fourier Transform (FFT) and when dealing with

quasiperiodic signals.

3.3.1 Fast Fourier Transform

We illustrate the approach by considering a signal of the form

𝑓𝜔 (𝑡) = 𝑒𝑖𝑡𝜔.

As shown in [5], this type of signal is related to a periodic dynamical system.

In practice, the

frequency 𝜔 is approximated using FFT, say ¯𝜔, and thus it is subject to error. Since 3G relies on

this information to approximate the SW obtained from the original signal, we can apply our work

to keep track of the stability of the method.

Let’s consider the case

√

3 = 1.7320508075688772...

𝜔1 =

and that we are able to recover four decimals of precision, say

𝜔2 = 1.7320641266595873... .

The Three Gap Theorem method will use the frequencies, as shown in 4.3.2,

and

𝜔 =

𝜔1
2𝜋

= 0.27566444771089604...

¯𝜔 =

𝜔2
2𝜋

= 0.27566656751002... .

40

Figure 3.3 We depict the results from 3.3.1 where 𝑇 = {0, 1, ..., 43}. Top Row: We plot the change
of the bottleneck distance as we add more decimals of precision to ¯𝜔. Bottom Row: Depicted are
the pairs achieving the bottleneck distance. We note our bound is almost identical to the bottleneck
distance.

In this case, 𝑥5 = ¯𝑥5, 𝑘5(𝑎) = 7, 𝑞6 = 43, and 𝑞7 = 51. Thus, by Theorem 3.2.8, for 𝑖 ≥ 0

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 (𝑆𝜔,43), 𝑑𝑔𝑚 𝑅

𝑖 (𝑆 ¯𝜔,43)) ≤

1
51

= 0.01960....

Furthermore, we can compute the exact distance in dimension 0 and 1, say 𝜖0, 𝜖1, respectively,

using the results from 4.2.2. We can then obtain the bound

𝑑𝐵 (𝑑𝑔𝑚𝑅

0 ( 𝑓𝜔 (𝑇)), 𝑑𝑔𝑚 𝑅

0 ( 𝑓 ¯𝜔 (𝑇)) = 5.3089 × 10−4 ≤ 2𝜋𝜖0 = 5.3277 × 10−4

and

𝑑𝐵 (𝑑𝑔𝑚𝑅

1 ( 𝑓𝜔 (𝑇)), 𝑑𝑔𝑚 𝑅

1 ( 𝑓 ¯𝜔 (𝑇)) = 1.4574 × 10−4 ≤ 2𝜋𝜖1 = 1.4651 × 10−4,

41

where 𝑇 = {0, 1, ..., 43}. This follows by noting that to recover distances in the complex plane we

can use the map 𝑥 → 2 sin(𝜋𝑥) as illustrated in 4.3. Furthermore, since

2 sin(𝜋𝑥) − 2 sin(𝜋𝑦) ≤ 4 sin(

𝜋(𝑥 − 𝑦)
2

) ≤ 2𝜋(𝑥 − 𝑦),

using (𝑥 − 𝑦) << 1, the upper bounds follow. Moreover, by combining these observations with

Theorem 3.2.8, we note

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑖 ( 𝑓𝜔 (𝑇)), 𝑑𝑔𝑚 𝑅

𝑖 ( 𝑓 ¯𝜔 (𝑇))) ≤

2𝜋
𝑞7

.

3.3.2 Quasiperiodic Functions

We now apply a similar analysis to the more complicated quasiperiodic signal

𝑓𝜔 (𝑡) = 𝑒𝑖𝑡 + 𝑒𝑖𝑡𝜔.

In this case, 3G consists in constructing a Cartesian product space to approximate (cid:8)𝑆𝑊𝑑,𝜏 ( 𝑓𝜔 (𝑡))(cid:9)𝑇

𝑡=0,

see 4.3. For this example, we denote the Cartesian product space as

𝐺𝜔,𝑇 = {𝑒𝑖𝑡 }𝑇

𝑡=0 × {𝑒𝑖𝑡𝜔}𝑇
𝑡=0

The persistence of this Cartesian product space can be computed using the persistence Künneth

formula [14] and the results in 4.3. We highlight that doing so does not require Ripser at any point.

Our error analysis consists in varying the decimal precision of 𝜔 to obtain a ¯𝜔. We then compare

the bottleneck distance as before but now of the sets 𝐺𝜔,𝑇 and 𝐺 ¯𝜔,𝑇 . Furthermore, as a direct

consequence of the Künneth formula, we can obtain analogue bounds to the ones in Section 3.3.1:

𝑑𝐵 (𝑑𝑔𝑚 𝑅

1 (𝐺𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

1 (𝐺 ¯𝜔,𝑇 )) ≤ 2𝜋𝜆1

and

𝑑𝐵 (𝑑𝑔𝑚 𝑅

2 (𝐺𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

2 (𝐺 ¯𝜔,𝑇 )) ≤ 2𝜋𝜆2,

where 𝜆 𝑗 = 𝑑𝐵 (𝑑𝑔𝑚 𝑅

𝑗 (𝑆 1

2 𝜋 ,𝑇 × 𝑆𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑗 (𝑆 1

2 𝜋 ,𝑇 × 𝑆 ¯𝜔,𝑇 )), which can be computed directly using

3G.

42

Let us consider the case 𝜔 =

√

5 and say we approximate with precision of five decimals

¯𝜔 = 2.23606. Computing 𝜔 and ¯𝜔 as before, we obtain 𝑥5 = ¯𝑥5, 𝑘5(𝜔) = 7, 𝑞6 = 100, and

𝑞7 = 121. To apply our results from Theorem 3.2.8 to the bound above, we let 𝑇 = 100. Noting the

bound, for 𝑗 = 1, 2,

where 𝜖𝑖 is as defined in Section 3.3.1, we conclude that

𝜆 𝑗 ≤ max{𝜖0, 𝜖1},

𝑑𝐵 (𝑑𝑔𝑚𝑅

𝑗 (𝐺𝜔,𝑇 ), 𝑑𝑔𝑚 𝑅

𝑗 (𝐺 ¯𝜔,𝑇 )) ≤

2𝜋
𝑞7

.

Figure 3.4 We depict the results from 3.3.2. Top Row: We plot the change of the bottleneck distance
as we add more decimals of precision to ¯𝜔. Bottom Row: Depicted are the pairs achieving the
bottleneck distance. We note our bound is still close to the bottleneck distance.

43

3.4 Concluding Remarks

We have successfully presented stability results that validate the use of The Three Gap Theorem

method under an adequate sample size. Furthermore, we are able to quantify the divergence in

the outcome of the approximation. This is relative to the number of matching continued fraction

exponents. It can also be related to decimal precision thanks to the introduction of the 𝑘 𝑗 function.

Overall, we hope this work helps in the further implementation of 3G and on its role of better

understanding and detecting quasiperiodic signals.

44

CHAPTER 4

ESTIMATION OF PERSISTENCE DIAGRAMS VIA THE THREE GAP THEOREM

We present theoretical and computational schemes to approximate the persistence diagrams of

sliding window embeddings from quasiperiodic functions. We do so by combining the Three Gap

Theorem from number theory, the Persistent Künneth formula from TDA, and our stability results

to and derive fast and provably correct persistent homology approximations. The input to our

procedure is the spectrum of the signal, and we provide numerical evidence of its utility to capture

the shape of toroidal attractors.

4.1 Preliminaries

Figure 4.1 depicts a system with a behavior more complex than periodicity. By using any

Figure 4.1 Pendulum attached to a sliding block.

measurement of the system, say the horizontal position 𝑥(𝑡) or the angular position 𝜃 (𝑡), we can

determine the qualitative aspects of the system and conclude it exhibits an oscillatory pattern.

Currently, a significant limitation is the computational time it takes to analyze the reconstructed

signal. Our work addresses this issue by detailing an estimation method of the standard approach

capable of proving results orders of magnitude faster. Explicitly, we provide an alternative to step

3 below, depicted in Figure 4.2:

45

1. Start with an observation signal 𝑓 .

2. Reconstruct the phase space by computing

{𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇

(cid:0)(cid:0)(cid:64)(cid:64)3. Compute dgm𝑅

𝑡=0.
𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2) using Ripser.

3. (a) Use the FFT to retrieve the frequencies of 𝑓 and

then use them to compute continued fraction expansions.

(b) Use them as shown in Section 4.2 and then apply

the results from Section 4.3 to approximate

dgm𝑅

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2).

Figure 4.2 The SW starts with a signal, followed by a reconstruction of the phase space, and
concludes with the computation of persistent homology. Top row: Signal obtained from the
periodic motion of an ideal pendulum. The sliding window point cloud recovers a circle, hence its
persistence diagram only has one hole in dimension 1. Bottom row: A quasiperiodic function with
two incommensurate frequencies (left), a PCA projection of its sliding window embedding (center)
recovers a 2-torus. The persistence diagrams verify this topology: two 1-dim holes, and one 2-dim
hole.

46

4.2 Main Results

Let us consider the Kronecker Sequence with one parameter 𝑆𝜔,𝑁 . We note (𝑆𝜔,𝑇 , ¯𝑑) is a metric

space and that the distribution of points in 𝑆𝜔,𝑇 can behave drastically different for parameters close

in 𝑑2.

4.2.1 Persistence Diagrams

We are able to compute dgm𝑅

0 (𝑆𝜔,𝑇 , ¯𝑑) applying the information from The Three Gap Theorem

which involves the 𝑖-th convergents.

Figure 4.3 We illustrate above how the 0-th dimensional persistent homology can be computed using
the Three Gap Theorem. At 𝜖 = 0 all of the points in (𝑆𝜔,𝑇 , ¯𝑑) count as a basis in 𝐻0(𝑅𝜖 (𝑆𝜔,𝑇 , ¯𝑑); F)
(birth time). At 𝜖 = 𝛿 𝐴 is the first instance points get connected, in fact the three possible gaps
are the only time when components are lost (death time). Top row: The set 𝑆𝜋−1,17 is illustrated
with the two gaps it creates. These can be traced in the barcodes and persistent diagram as shown.
Bottow row: The set 𝑆√
5,17 creates three gaps. Due to clustering, it is a poor covering of the circle.
We note 𝑑2(𝜋 − 1,

5) = 0.094.

√

Theorem 4.2.1. Let 𝑇 ∈ N and 𝜔 ∈ R\Q with continued fraction [𝑎1, 𝑎2, 𝑎3, · · · ] and 𝑖-th
𝑝𝑖
𝑞𝑖 . We assume there are three different intervals generated by 𝑆𝜔,𝑇 and that 𝛿𝐶 < 1/2.

convergent

If 𝑘 ≥ 0 is the unique integer for which

𝑞𝑘 + 𝑞𝑘−1 ≤ 𝑇 < 𝑞𝑘 + 𝑞𝑘+1,

47

𝐷𝑖 = 𝑞𝑖𝜔 − 𝑝𝑖, and r,s are the unique integers satisfying

𝑇 = 𝑟𝑞𝑘 + 𝑞𝑘−1 + 𝑠,

1 ≤ 𝑟 ≤ 𝑎𝑘+1,

0 ≤ 𝑠 ≤ 𝑞𝑘 − 1,

then 𝐻0(𝑅𝜖 (𝑆𝜔,𝑇 , ¯𝑑); F) =

F𝑇+1

0 ≤ 𝜖 < |𝐷 𝑘 |

F𝑞𝑘

|𝐷 𝑘 | ≤ 𝜖 < |𝐷 𝑘+1| + (𝑎𝑘+1 − 𝑟)|𝐷 𝑘 |

F𝑞𝑘 −𝑠−1

|𝐷 𝑘+1| + (𝑎𝑘+1 − 𝑟)|𝐷 𝑘 | ≤ 𝜖

< |𝐷 𝑘+1| + (𝑎𝑘+1 − 𝑟 + 1)|𝐷 𝑘 |

.




F

|𝐷 𝑘+1 | + (𝑎𝑘+1 − 𝑟 + 1)|𝐷 𝑘 | ≤ 𝜖


We also compute dgm𝑅

1 (𝑆𝜔,𝑇 , ¯𝑑) using information from the Three Gap Theorem. We assume
𝑆𝜔,𝑇 has three different gaps and the notation from before. The case in which it has only two gaps

can be handle analogously.

Theorem 4.2.2. Suppose 𝛿𝐶 < 1/3. Let Γ be the set containing the pairs (𝑥, 𝑦), 𝑥, 𝑦 ∈ 𝑆𝜔,𝑇 and 𝑥 <

𝑦, for which there exists a 𝑧 ∈ ([0, 𝑥) ∪ (𝑦, 1)) ∩ 𝑆𝜔,𝑇 such that 1 − ¯𝑑 (𝑥, 𝑦) ≤ 2 ¯𝑑 (𝑦, 𝑧) ≤ 2 ¯𝑑 (𝑥, 𝑦).

If

then 𝐻1(𝑅𝜖 (𝑆𝜔,𝑇 , ¯𝑑); F) =

𝜆 = min{ ¯𝑑 (𝑥, 𝑦)| (𝑥, 𝑦) ∈ Γ},

F |𝐷 𝑘+1 | + (𝑎𝑘+1 − 𝑟 + 1)|𝐷 𝑘 | ≤ 𝜖 < 𝜆

.

∅ 𝑒𝑙𝑠𝑒






Proof. Let {𝑥𝑖}𝑇
(cid:205)𝑇

𝑖=0 denote the points of 𝑆𝜔,𝑇 in ascending order. Consider the 1-simplex 𝜎 =
𝑖=0 {𝑥𝑖, 𝑥𝑖⊕1} where ⊕ is mod(T+1) addition. We note 𝜕 (𝜎) = 0 and 𝜎 ∈ 𝑅𝜖 (𝑆𝜔,𝑇 , ¯𝑑) if and only

48

if 𝜖 ≥ 𝛿𝐶 = |𝐷 𝑘+1 | + (𝑎𝑘+1 − 𝑟 + 1)|𝐷 𝑘 | by the Three Gap Theorem. Furthermore, any 1-simplex

𝜎0 for which 𝜕 (𝜎0) = 0 is homologous to 𝜎.

Now we demonstrate there are no 2-simplex 𝜏 ∈ 𝑅𝛿𝐶 (𝑆𝜔,𝑇 , ¯𝑑) such that 𝜕 (𝜏) = 𝜎. We

assume this to be true and arrive at a contradiction. By reindexing and translating, we can assume

¯𝑑 (0, 𝑥1) = 𝛿𝐶. By assumption, there needs to exist a {0, 𝑥1, 𝑥𝑖} ∈ 𝜏2 for some 𝑖 > 1. If 𝑥𝑖 ≤ 1/2

then ¯𝑑 (0, 𝑥𝑖) = ¯𝑑 (0, 𝑥1) + ¯𝑑 (𝑥1, 𝑥𝑖) > 𝛿𝐶 and if 𝑥𝑖 > 1/2, then ¯𝑑 (𝑥1, 𝑥𝑖) =

min{ ¯𝑑 (0, 𝑥1) + ¯𝑑 (0, 𝑥𝑖), 1 − ( ¯𝑑 (0, 𝑥1) + ¯𝑑 (0, 𝑥𝑖))} > 1/3,

(𝛿𝐶 < 1/3). Thus, {0, 𝑥1, 𝑥𝑖} ∉ 𝜏2 for all 𝑖 > 1, a contradiction.

To conclude the proof, it suffices to show that for 𝜖 ∈ (𝛿𝐶, 𝜆] there exists a 2-simplex 𝜏 ∈

𝑅𝜖 (𝑆𝜔,𝑇 , ¯𝑑) such that 𝜕 (𝜏) = 𝜎 if and only if 𝜖 = 𝜆. We first show the forward direction.

, 𝑥𝑛2) ∈ Γ be such that 𝜆 = ¯𝑑 (𝑥𝑛1

Let (𝑥𝑛1
described in the statement of the theorem. We note ¯𝑑 (𝑥𝑛2
, 𝑥𝑛3) ≤ ¯𝑑 (𝑥𝑛2
1 − ¯𝑑 (𝑥𝑛1
, 𝑥𝑛3) ≤ 𝜆. Thus, {𝑥𝑛1

, 𝑥𝑛2) − ¯𝑑 (𝑥𝑛2

, 𝑥𝑛3) ≤ ¯𝑑 (𝑥𝑛1
, 𝑥𝑛2

, 𝑥𝑛2) and 𝑥𝑛3 be a point associated with the pair as
, 𝑥𝑛2) = 𝜆 and ¯𝑑 (𝑥𝑛21, 𝑥𝑛3) ≤

, 𝑥𝑛3 } ∈ 𝑅𝜆 (𝑆𝜔,𝑇 , ¯𝑑). Let

∑︁

𝜏1 =

𝑛1≤ 𝑖 <𝑛2−1

{𝑥𝑖, 𝑥𝑖+1, 𝑥𝑛2 }.

By keeping track of the index at zero, one can define analogues 𝜏2, 𝜏3 corresponding to the pairs
, 𝑥𝑛3 }. One can verify 𝜏 ∈ 𝑅𝜆 (𝑆𝜔,𝑇 , ¯𝑑)

, 𝑥𝑛1), respectively. Let 𝜏 = 𝜏1+𝜏2+𝜏3+{𝑥𝑛1

, 𝑥𝑛2

, 𝑥𝑛3), (𝑥𝑛3

(𝑥𝑛2
and that 𝜕 (𝜏) = 𝜎 as desired.

For the reverse direction, let 𝜖 ∈ (𝛿𝐶, 𝜆).

If one assumes the existence of such 2-simplex,

say 𝜏0 ∈ 𝑅𝜖 (𝑆𝜔,𝑇 , ¯𝑑), one can show this will contradict the minimality of 𝜆 by considering

max{ ¯𝑑 (𝑥, 𝑦) | there exists {𝑥, 𝑦, 𝑧} ∈ 𝜏2

0 }.

□

4.2.2 Bottleneck Distance

We can apply our result to obtain explicit formulas for the the stability with respect to the

bottleneck distance as detailed in Section 2. These results provide a way of obtaining the 𝜖0 and

𝜖1 shown in 3.3.1 without using Ripser or Persim [16, 36]. We present our results following the

49

notation used in Section 2. Furthermore, we assume there are three gaps and that the biggest

one has multiplicity of at least two. We also impose assumptions in the statement for the sake of

simplifying the presentation. The proof for the other cases is analogues.

Theorem 4.2.3. Suppose

for some 𝑗 ∈ N and

𝑥 𝑗 = ¯𝑥 𝑗

𝑞𝑘 𝑗 (𝛼)−1 + 𝑞𝑘 𝑗 (𝛼)−2 ≤ 𝑁 < 𝑞𝑘 𝑗 (𝛼)−1 + 𝑞𝑘 𝑗 (𝛼).

Let 𝑟, 𝑠 be the unique integers satisfying

𝑁 = 𝑟𝑞𝑘 𝑗 (𝛼)−1 + 𝑞𝑘 𝑗 (𝛼)−2 + 𝑠,

1 ≤ 𝑟 ≤ 𝑎𝑘 𝑗 (𝛼),

0 ≤ 𝑠 ≤ 𝑞𝑘 𝑗 (𝛼)−1 − 1.

We assume 𝑟 < 𝑎𝑘 𝑗 (𝛼), 𝑁 − 𝑞𝑘 𝑗 (𝛼)−1 > 𝑠, 𝛼 > ¯𝛼, and 𝑘 𝑗 (𝛼) is even. Let 𝐷𝑖 (𝑥) = |𝑞𝑖𝑥 − 𝑝𝑖 |,

𝐾1(𝑥) = 𝐷 𝑘 𝑗 (𝛼) (𝑥) + (𝑎𝑘 − 𝑟)𝐷 𝑘 𝑗 (𝛼)−1(𝑥), 𝐾2(𝑥) = 𝐷 𝑘 𝑗 (𝛼) (𝑥) + (𝑎𝑘 − 𝑟 + 1)𝐷 𝑘 𝑗 (𝛼)−1(𝑥)

and

Δ1 = (𝛼 − ¯𝛼)(𝑞𝑘 𝑗 (𝛼) − 𝑞𝑘 𝑗 (𝛼)−1(𝑎𝑘 𝑗 (𝛼) − 𝑟)), Δ2 = (𝛼 − ¯𝛼) (𝑞𝑘 𝑗 (𝛼) − 𝑞𝑘 𝑗 (𝛼)−1(𝑎𝑘 𝑗 (𝛼) − 𝑟 + 1)).

Then,

𝑑𝐵 (𝑑𝑔𝑚 𝑅

0 (𝑆𝛼,𝑁 ), 𝑑𝑔𝑚 𝑅

0 (𝑆 ¯𝛼,𝑁 )) = max{min{Δ1,

𝐾1( ¯𝛼)
2

}, min{Δ2,

𝐾2( ¯𝛼)
2

}}.

Proof. For 1 ≤ 𝑖 ≤ 3, let 𝛿𝑖 and ¯𝛿𝑖 denote the length of the three gaps in increasing size of 𝑆𝛼,𝑁

and 𝑆 ¯𝛼,𝑁 respectively. The assumptions imply ¯𝛿3 > 𝛿3, ¯𝛿2 > 𝛿2, and that we can match each point

(0, 𝛿𝑖) with a (0, ¯𝛿𝑖). This leads to a ( ¯𝛿2 − 𝛿2)-matching since

𝛿1 − ¯𝛿1 = (𝛼 − ¯𝛼)𝑞𝑘 𝑗 (𝛼)−1 ≤ (𝛼 − ¯𝛼) (𝑞𝑘 𝑗 (𝛼)−2 + 𝑞𝑘 𝑗 (𝛼)−1𝑟) = 𝛿2 − ¯𝛿2

and ( ¯𝛿3 − 𝛿3) = ( ¯𝛿2 − 𝛿2) − (𝛿1 − ¯𝛿1). Now, we argue matching (0, ¯𝛿3) with (0, 𝛿2) or (0, 𝛿1) leads

to a worse matching. Indeed, this follows by noting

¯𝛿3 − 𝛿2 = ( ¯𝛿2 − 𝛿2) + ¯𝛿1

50

and ¯𝛿3 − 𝛿1 > ¯𝛿3 − 𝛿2. Thus, the only possibly better matches can be obtain from matching to
(0, 𝛿3) or to ( ¯𝛿3
, ¯𝛿3
2 ),. The best matching obtained from the latter scenario is a
¯𝛿3 > 𝛿3. For the former case, the situation can only improve by matching (0, ¯𝛿2) with ( ¯𝛿2
with (0, 𝛿1), but

¯𝛿3
2 -matching since
2 ) or

, ¯𝛿2

2

2

¯𝛿2 − 𝛿1 = ( ¯𝛿2 − 𝛿2) + (𝛿2 − 𝛿1).

We conclude the best matching is a r-matching, where

𝑟 = max{min{( ¯𝛿2 − 𝛿2),

¯𝛿2
2

}, min{( ¯𝛿3 − 𝛿3),

¯𝛿3
2

}}.

The result now follows by substituting the expressions of the gaps using The Three Gap Theorem

and noting that for 0 ≤ 𝑖 ≤ 𝑘 𝑗 (𝛼)

𝛼 >

𝑝𝑖
𝑞𝑖

iff ¯𝛼 >

𝑝𝑖
𝑞𝑖

.

□

Similarly, we can apply the results of The Three Gap Theorem to obtain the bottleneck distance

in the case of the 1-st dimensional persistence diagrams. We make the assumptions in the previous

results and use the notation 𝜆𝛼 to denote the 𝜆 from Theorem 4.2.2 corresponding to the sampling

parameter 𝛼. As before, the assumptions made are to simplify the presentation of the statement.

Theorem 4.2.4. Let

then

𝑇 (𝑥) =

1
2

(𝜆𝑥 − 𝐾2(𝑥)),

𝑑𝐵 (𝑑𝑔𝑚𝑅

1 (𝑆𝑎,𝑁 ), 𝑑𝑔𝑚 𝑅

1 (𝑆 ¯𝑎,𝑁 )) = min{max{Δ2, |𝜆𝑎 − 𝜆 ¯𝑎 |}, max{𝑇 (𝑎), 𝑇 ( ¯𝑎)}}.

4.3 Approximation method

We present our method in the next theorem in which we denote by 𝜎𝑚𝑖𝑛 and 𝜎𝑚𝑎𝑥 the smallest

and biggest, respectively, eigenvalue of 𝐴. We also denote by 𝑑𝐺𝐻 the Gromov-Hausdorff distance

computed in (C𝑁 , 𝑑∞), where 𝑑∞ denotes the supremum metric.

51

Figure 4.4 Depiction of 𝜙(13) (red triangles) and 𝐺 (13) (blue dots) in the case of two incommen-
surate frequencies. Left: Points are in the torus. Right: Points are plotted using the rectangular
representation of the torus. The pair marked with an x achieve the Hausdorff distance between the
two sets.

Theorem 4.3.1. Let

𝑓 (𝑡) =

𝑀
∑︁

𝑗=1

𝑐 𝑗 𝑒𝑖𝜔 𝑗 𝑡,

𝑐 𝑗 ∈ C, 𝜔 𝑗 ∈ R.

Suppose 𝑁 ≤ 𝑀 is the maximum number of incommensurate frequencies, say {𝜔𝑙 }𝑁

𝑙=1. For

𝑁 < 𝑙 ≤ 𝑀 and ¯𝑡 ∈ N𝑁 let

¯𝜔𝑙 (¯𝑡) =

𝑁
∑︁

𝑗=1

¯𝑡 𝑗𝜆𝑙

𝑗 𝜔 𝑗 ,

where 𝜆𝑖

𝑗 ∈ Q and ¯𝜔𝑙 ((1, 1, . . . , 1)) = 𝜔𝑙 . Given 𝑇 ∈ N and ¯𝑐𝑖 =

√

𝑑 + 1 𝑐𝑖 we define 𝜙(𝑇) =

{( ¯𝑐1𝑒𝑖𝜔1𝑡, . . . , ¯𝑐𝑀 𝑒𝑖𝜔𝑀 𝑡)⊺}0≤𝑡≤𝑇 ,

𝐺𝑇 =

𝑘 = max{𝜎−1

𝑚𝑖𝑛, 𝜎𝑚𝑎𝑥

√

satisfies

{( ¯𝑐1𝑒𝑖𝜔1 ¯𝑡1, . . . , ¯𝑐𝑁+1𝑒𝑖 ¯𝜔𝑁 +1 (¯𝑡), . . . , ¯𝑐𝑀 𝑒𝑖 ¯𝜔𝑀 (¯𝑡))⊺}0≤¯𝑡𝑖 ≤𝑇 ,

𝑀 } and 𝜆 = 𝑑𝐺𝐻 (𝜙(𝑇), 𝐺 (𝑇 ′)), where 𝑇 ′ ≤ 𝑇. If (𝑎1, 𝑏1) ∈ dgm𝑅

𝑗 (𝐺 (𝑇 ′), 𝑑∞)

then there exists a unique (𝑎2, 𝑏2) ∈ dgm𝑅

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2) satisfying

𝑏1 − 2𝜆
𝑎1 + 2𝜆

> max{𝑘 2, 1},

max

(cid:26)

0,

(cid:27)

𝑏1 − 2𝜆
𝑘

≤ 𝑏2 ≤ 𝑘 (𝑏1 + 2𝜆)

max

(cid:26)

0,

(cid:27)

𝑎1 − 2𝜆
𝑘

≤ 𝑎2 ≤ 𝑘 (𝑎1 + 2𝜆).

52

Proof. Our first task will be to apply the Stability Theorem [37] to relate (𝑎1, 𝑏1) with a unique
(𝑎0, 𝑏0) ∈ dgm𝑅
with a unique (𝑎2, 𝑏2) ∈ dgm𝑅

𝑗 (𝜙(𝑇), 𝑑∞). In turn, by using the Isometry Theorem [38] we will match (𝑎0, 𝑏0)

, 𝑑2). We conclude the result by keeping track of

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

the relations throughout.

By Theorem 2.1.14, if 𝑏1 − 𝑎1 > 4𝜆 then there exists a unique (𝑎0, 𝑏0) ∈ dgm𝑅

𝑗 (𝜙(𝑇), 𝑑∞) such

that |𝑏1 − 𝑏0| < 2𝜆 and |𝑎1 − 𝑎0| < 2𝜆. To relate (𝑎0, 𝑏0) to a (𝑎2, 𝑏2) ∈ dgm𝑅

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2)

we note

and

𝜎𝑚𝑖𝑛𝑑∞(𝑣𝑡1

, 𝑣𝑡2) ≤ 𝑑2( 𝐴𝑣𝑡1

, 𝐴𝑣𝑡2)

𝑑2( 𝐴𝑣𝑡1

, 𝐴𝑣𝑡2) ≤

√

𝑀𝜎𝑚𝑎𝑥 𝑑∞(𝑣𝑡1

, 𝑣𝑡2),

where 𝑆𝑊𝑑,𝜏 𝑓 (𝑡) = 𝐴𝑣𝑡. This allow us to find an interleaving between the two persistent modules

using a logarithmic scale, namely a 𝑙𝑛(𝑘) interleaving. Thus, by the Isometry Theorem we can say

that if

𝑏0
𝑎0

> 𝑘 2 then there exists a unique (𝑎2, 𝑏2) ∈ dgm𝑅

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2) such that

1
𝑘

< max{𝑏0, 𝑏2}
min{𝑏0, 𝑏2}

, max{𝑎0, 𝑎2}
min{𝑎0, 𝑎2}

< 𝑘.

Putting it all together, we can say that if

𝑏1 − 2𝜆
𝑎1 + 2𝜆

> max{𝑘 2, 1}

then there exist a unique (𝑎0, 𝑏0) ∈ dgm𝑅

𝑗 (𝜙(𝑇), 𝑑∞) such that 𝑎0 < 𝑎1 + 2𝜆 and 𝑏0 > 𝑏1 − 2𝜆.

Thus,

implying there exists a unique

𝑏0
𝑎0

>

𝑏1 − 2𝜆
𝑎1 + 2𝜆

> 𝑘 2

(𝑎2, 𝑏2) ∈ dgm𝑅

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2)

such that, say if 𝑏0 = max{𝑏0, 𝑏2},

1
𝑘

<

𝑏0
𝑏2

< 𝑘

53

and thus

𝑏1 − 2𝜆
𝑘

< 𝑏2 < 𝑘 (𝑏1 + 2𝜆).

The same result holds when 𝑏2 = max{𝑏0, 𝑏2}. One can show an analogues inequality for the birth

time as well. The result follows.

□

We note that for the case 𝑁 = 𝑀, dgm𝑅

𝑗 (𝐺 (𝑇 ′), 𝑑∞) can be computed for 𝑗 = 1, 2 using our

main results. This is achieved by using the persistent Künneth formula [14]. We are able to leverage

our results since

𝑑2(𝑐𝑙 𝑒𝑖𝜔𝑙𝑡1, 𝑐𝑙 𝑒𝑖𝜔𝑙𝑡2) = 2|𝑐𝑙 |sin(𝜋 ¯𝑑 ( [𝜔′

𝑙𝑡1], [𝜔′

𝑙𝑡2])),

for 𝑡1, 𝑡2 ∈ N and 𝑤′

𝑙 = 𝜔𝑙

2𝜋 . Furthermore, by noting

𝑑2(𝑐𝑙 𝑒𝑖𝜔𝑙𝑡1, 𝑐𝑙 𝑒𝑖𝜔𝑙𝑡2) = 𝑑2(𝑐𝑙 𝑒−𝑖𝜔𝑙𝑡1, 𝑐𝑙 𝑒−𝑖𝜔𝑙𝑡2)

we can obtain the following result:

Theorem 4.3.2. Let

𝑓 (𝑡) =

𝑁
∑︁

𝑗=1

𝑐 𝑗 (𝑒𝑖𝜔 𝑗 𝑡 + 𝑒−𝑖𝜔 𝑗 𝑡),

𝑐 𝑗 ∈ C, 𝜔 𝑗 ∈ R,

where {𝜔 𝑗 }𝑁
2𝜋 , and for 𝐼 = [𝑎, 𝑏)
we let ¯𝐼 = [ ¯𝑎, ¯𝑏) and write 𝑠𝐼 to denote [𝑠𝑎, 𝑠𝑏). Given 𝑇 ∈ N and ¯𝑐𝑙 = √︁(4𝑁 + 2)|𝑐𝑙 |2, we define

𝑗=1 is incommensurate. For 𝑎, 𝑏, 𝑠 ∈ R, let ¯𝑎 = 2 sin(𝜋𝑎), 𝜔′

𝑙 = 𝜔𝑙

𝜙(𝑇) =

𝐺𝑇 =

{( ¯𝑐1𝑒𝑖𝜔1𝑡, . . . , ¯𝑐𝑁 𝑒𝑖𝜔𝑁 𝑡)⊺}0≤𝑡≤𝑇 ,

{( ¯𝑐1𝑒𝑖𝜔1 ¯𝑡1, . . . , ¯𝑐𝑁 𝑒𝑖𝜔𝑁 ¯𝑡 𝑁 )⊺}0≤¯𝑡𝑖 ≤𝑇 ,

𝑘 = max{𝜎−1

𝑚𝑖𝑛, 𝜎𝑚𝑎𝑥

√

𝑁 } and 𝜆 = 𝑑𝐺𝐻 (𝜙(𝑇), 𝐺 (𝑇 ′)), where 𝑇 ′ ≤ 𝑇. If [𝑎1, 𝑏1) ∈

(cid:216)

(cid:205)𝑙 𝑟𝑙= 𝑗

(cid:8) ¯𝑐1 ¯𝐼1 ∩ · · · ∩ ¯𝑐𝑁 ¯𝐼𝑁 | 𝐼𝑙 ∈ 𝑏𝑐𝑑 𝑅

𝑟𝑙 (𝑆𝜔′

𝑙

,𝑇 ′, ¯𝑑)(cid:9),

54

satisfies

𝑏1 − 2𝜆
𝑎1 + 2𝜆

> max{𝑘 2, 1},

then there exists a unique (𝑎2, 𝑏2) ∈ dgm𝑅

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝑡)}𝑇
𝑡=0

, 𝑑2) satisfying

max

(cid:26)

0,

(cid:27)

𝑏1 − 2𝜆
𝑘

≤ 𝑏2 ≤ 𝑘 (𝑏1 + 2𝜆)

max

(cid:26)

0,

(cid:27)

𝑎1 − 2𝜆
𝑘

≤ 𝑎2 ≤ 𝑘 (𝑎1 + 2𝜆).

We note that our results allow us to approximate (𝑎2, 𝑏2) with an error bound. Indeed, when

the right conditions are met, the theorem guarantees the point to lie inside of the rectangle with the

vertices shown below:

(𝑥0, 𝑦1)

(𝑥1, 𝑦1)

(𝑎2, 𝑏2)

(𝑥0, 𝑦0)

(𝑥1, 𝑦0)

where 𝑥0 = max{0, 𝑎1−2𝜆
approximation (𝑎1, 𝑏1) provides the rectangle as a confidence region.

}, 𝑥1 = 𝑘 (𝑎1 + 2𝜆), 𝑦0 = max{0, 𝑏1−2𝜆

𝑘

𝑘

}, and 𝑦1 = 𝑘 (𝑏1 + 2𝜆). Thus, our

4.3.1 Example

We apply our method to the function 𝑓 (𝑡) = 1√
2

𝑒𝑖

√
5𝑡 + 1√
2

√
3𝑡. By appropriate choice of

𝑒𝑖

parameters, 𝑑 = 1, and 𝜏 =

√

𝜋
5−

√
3

, we can make the matrix A be orthogonal. We also compare our

method to two other approximations of 𝑆𝑊 (𝑛): landmarks SW𝐿 (𝑛) and landmarks for the Künneth

formula K𝐿 (𝑛), where 𝑛 is the number of samples used [14].

4.4 Concluding remarks

In this chapter we introduced the approximation method 3G. It approximates the persistence

diagram from SW. We were able to compute much faster results by using frequency information

thanks to the Three Gap Theorem. We then combined results obtained from single frequencies

to treat the quasiperiodic case using the Künneth formula in persistent homology. Under the

right conditions, our approximation can include error bounds obtained from the Isometry Theorem

55

Figure 4.5 Diagrams from example 4.1. We also illustrate the running time it takes when using
Ripser [16] for 𝑆𝑊 and 𝐾𝐿. Our code is denoted as 𝐾3𝐺; we note Ripser is not used at any step in
our code.

and the Stability Theorem. We hope our approach contributes to the use of the sliding window

embedding technique on large data sets.

56

CHAPTER 5

APPLICATIONS

In this section we will apply our pipeline, 3G, to the following examples: Synthetic tremor signals,

in the nonlinear tuned vibration absorber it is an indicator of safe or unsafe operations [39].

Neuroscience, where quasiperiodicity is associated to task-specific functions in certain brain areas

[40]. Celestial mechanics, where quasiperiodicity translates to nice trajectories for a mission [41].

We contribute to this fields by providing a faster alternative (see Figure 4.5 and Table 5.1) to analyze

the reconstructed signal of potential quasiperiodic signals.

5.1 Methodology

We work with dynamical systems for which Φ is known.

In particular, we work with or-

dinary differential equations associated to the system. These equations dictate the evolution of

the system across time. To imitate an observation of the system, we use the solutions to cre-

ate a time series signal

𝑓 . We reconstruct the system by then computing the sliding window

point cloud {𝑆𝑊𝑑,𝜏 𝑓 (𝜖𝑡)}𝑇

𝑡=0, where 𝜖 is chosen following the Nyquist-Shannon sampling theo-
rem [42]. The framework detailed in [5] provides a way of selecting 𝑑 and 𝜏 and approximate
dgm𝑅
𝑗 ({𝑆𝑊𝑑,𝜏𝑆1 𝑓 (𝜖𝑡)}𝑇
𝑡=0

, 𝑑2). Since 𝑆1 𝑓 is a sum of the form

𝑗 ({𝑆𝑊𝑑,𝜏 𝑓 (𝜖𝑡)}𝑇
𝑡=0

, 𝑑2) with dgm𝑅

shown in 4.3, we can apply our results to approximate the latter and as a result get close to the

former. Indeed, in the examples we consider

2
∑︁

𝑆1 𝑓 (𝜖𝑡) =

𝑐 𝑗 (𝑒𝑖𝜖𝜔 𝑗 𝑡 + 𝑒−𝑖𝜖𝜔 𝑗 𝑡).

Thus, 3G can be used to approximate dgm𝑅

𝑗=1
𝑗 ({𝑆𝑊𝑑,𝜏𝑆1 𝑓 (𝜖𝑡)}𝑇
𝑡=0

, 𝑑2) with a known error bound as

detailed in 4.3.

5.2 Results

We consider systems with known quasiperiodic behavior. Specifically, they are systems with two

incommensurate frequencies. Hence, the sliding window point cloud, {𝑆𝑊𝑑,𝜏 𝑓 (𝜖𝑡)}𝑇

𝑡=0, samples
a 2-torus. We label its corresponding persisetence diagrams by 𝑆𝑊 (𝑇) (shown in brown dots).

Similarly, we denote by 𝑆𝑊𝑆1 (𝑇) the ones corresponding to 𝑆1( 𝑓 ) (shown in blue squares). The

57

approximation obtained from 3𝐺 is labeled as 𝐾3𝐺 (𝑇) (shown in orange crosses).

5.2.1 Double Gyre

The driven Double Gyre system provides a model for patterns occurring in geophysical flows

[43]. A topological approach to it has shown great success.

Indeed, as is shown in [44] there

are multiple topological classifications that can be observed in the system which represents the

motion of a fluid particle. In this work they noted the motion of fluid can be sparse, meaning some

initial conditions will imply a particle will be contained in a small region. Indeed, different initial

conditions have shown trajectories with a topology of a standard strip, five-handle structure with a

torsion, torus and a möbius strip. This result was achieved by calculating homologies in a branched

manifold while keeping track of the orientability chains which allowed for the identification of the

branches and the localization of twists or torsions. This approach is called Branched Manifold

analysis through Homologies (BraMAH). At its core, it follows the principle of reconstructing a

dynamical system from a time series. Yet the details at each step are completely different to SW.

Figure 5.1 Top row: Phase space of the driven Double Gyre system starting at 𝑥0. The solution
𝑥1(𝑡) we used for SW is plotted in the middle followed by its fast fourier transform. Bottom row:
Corresponding diagrams.

58

Nevertheless, for the corresponding initial conditions, we independently corroborated the presence

of a torus, i.e. quasiperiodicity in the motion of the fluid particle, see Figure 5.1. This highlight

the robustness offered by a topological approach to dynamical systems.

Concretely, the driven Double Gyre is given by the equation (cid:164)𝑥 = (−

𝜕𝜓
𝜕𝑥2

, 𝜕𝜓
𝜕𝑥1

) where

with

𝜓(𝑥1, 𝑥2, 𝑡) = 𝐴 sin(𝜋𝑔(𝑥1, 𝑡)) sin(𝜋𝑥2)

𝑔(𝑥1, 𝑡) = 𝜂 sin(𝜆𝑡) (𝑥2

1 − 2𝑥1) + 𝑥1,

and parameter values 𝐴, 𝜇 = 0.1 and 𝜆 = 𝜋/5. The system has a toroidal attractor for the initial

condition 𝑥0 = (0.5, 0.625) [44]. We solve the system using 𝑑𝑡 = 0.1 up to 𝑡 = 800. Furthermore,

SW is done using 𝑓 = 𝑥1, 𝑑 = 4, 𝜏 = 119.03, and 𝜖 = 0.1.

5.2.2 Torus in R4

We move on to consider the torus in rectangular coordinates. It is given by the set of equations:

(cid:164)𝑥 = −𝑦 + 𝑥(1 −

√︃

𝑥2 + 𝑦2),

(cid:164)𝑦 = 𝑥 + 𝑦(1 −

√︃

𝑥2 + 𝑦2),

(cid:164)𝑧 = −𝑘𝑟 + 𝑧(4 −

√︁𝑧2 + 𝑟 2),

(cid:164)𝑟 = 𝑘 𝑧 + 𝑟 (4 −

√︁𝑧2 + 𝑟 2).

When 𝑘 is irrational, the solutions span a torus [45]. We considered the case 𝑘 =

√

2 with initial

condition (𝑥0, 𝑦0, 𝑧0, 𝑟0) = (1, 0, 4, 0). We solve the system with 𝑑𝑡 = 0.1 up to 𝑡 = 800. SW is

Figure 5.2 Diagrams obtained for the torus in R4.

59

done using 𝑓 = 𝑥 + 𝑧, 𝑑 = 4, 𝜏 = 87.56, and 𝜖 = 0.1, see Figure 5.2.

5.2.3 Pendulum Attached to a Sliding Block

We now consider a particular type of tune mass damper (TMD). TMD is a device used to

suppress vibration by moving a mass attached to the main structure through springs and dampers

[22]. This model has been successfully used to dampen the effect of long duration earthquake

ground motions [21].

Indeed, one can consider the main structure being a skyscraper and the

incoming wave being generated by the earthquake. Thus, a better understanding of this system

translates to earthquake-resistant technologies.

We consider the case of a pendulum attached to a sliding block, see Figure 5.3. It was shown

to exhibit quasiperiodicity in [22]. The governing equations are:

(cid:165)𝑥 + 𝛼2𝑥 − ¯𝜖𝑔𝜃 − ¯𝜖 𝐿 (cid:164)𝜃2𝜃 = 0

(cid:165)𝜃 + (1 + ¯𝜖) 𝛽2𝜃 − ¯𝜖 ℎ𝛼2𝑥 + ¯𝜖 (cid:164)𝜃2𝜃 = 0

Figure 5.3 Top row: Plot of (𝑥, (cid:164)𝑥) from the pendulum attached to a sliding block. The solution
𝑥(𝑡) we used for SW is plotted in the middle followed by its fast fourier transform. Bottom row:
Persistence diagrams.

60

in which

¯𝜖 =

𝑚
𝑀

, 𝛼 =

√︂ 𝑘
𝑀

, 𝛽 =

√︂ 𝑔
𝐿

, ℎ =

,

1
¯𝜖 𝐿

and 𝑔 is the acceleration of gravity. The parameter values are 𝑚 = 0.5, 𝑀 = 1, 𝐿 = 1, and 𝑘 = 5.

We solve the system with initial condition (𝑥0, (cid:164)𝑥0, 𝜃0, (cid:164)𝜃0) = (0.1, 0, −0.1, 0) and 𝑑𝑡 = 0.27 up to

𝑡 = 540. SW is done using 𝑓 = 𝑥, 𝑑 = 4, 𝜏 = 108.05, and 𝜖 = 0.027.

5.2.4 Generalized Wilson-Cowan Equations

We consider a generalized version of the Wilson-Cowan equations shown in Example 2.2.3.

Traditionally, these equations are derived via a time-coarse graining technique that averages the

response. They also restrict to a weak Gamma distribution of time delays. The extended model has

been treated in [40], it is given by

(cid:164)𝑢(𝑡) = −𝑢(𝑡) + 𝑓1

(cid:164)𝑣(𝑡) = −𝑣(𝑡) + 𝑓2

(cid:16)

(cid:16)

𝜃𝑢 +

𝜃𝑣 +

∫ 𝑡

−∞
∫ 𝑡

−∞

ℎ(𝑡 − 𝑠) (𝑎𝑢(𝑠) + 𝑏𝑣(𝑠)) 𝑑𝑠

(cid:17)

,

ℎ(𝑡 − 𝑠) (𝑐𝑢(𝑠) + 𝑑𝑣(𝑠)) 𝑑𝑠

(cid:17)

,

Figure 5.4 Top row: Plot of (𝑢, 𝑣) from the generalized Wilson-Cowan equations. The solution
𝑢(𝑡) we used for SW is plotted in the middle followed by its fast fourier transform. Bottom row:
Persistence diagrams.

61

where 𝑢(𝑡) and 𝑣(𝑡) model the firing activity in two neuronal populations, a,b,c,d are the connection

weights and 𝜃𝑢, 𝜃𝑣 are background drives. The activation functions 𝑓1, 𝑓2 are smooth and increasing

on the real line. The authors considered three types of delay kernel ℎ : [0, ∞) → [0, ∞), namely,

weak Gamma, strong Gamma, and Dirac kernel. Their analysis indicates a kernel is preferable

based on the function of the model populations.

One can readily verify this system recovers the Wilson-Cowan equations in the case of a weak

Gamma kernel. We consider the case of a Dirac kernel, see Figure 5.4. This system has been

shown to exhibit quasiperiodicity [40, 46]. Furthermore, it has applications to the subthalamic

nucleus - globus pallidus network involved in Parkinson’s Disease (guided by anatomical and

electrophysiological research) [47]. The system of interest becomes

(cid:164)𝑢(𝑡) = −𝑢(𝑡) + 𝑓1(𝜃𝑢 + 𝑎𝑢(𝑡 − 𝜏1) + 𝑏𝑣(𝑡 − 𝜏2)),

(cid:164)𝑣(𝑡) = −𝑣(𝑡) + 𝑓2(𝜃𝑣 + 𝑐𝑢(𝑡 − 𝜏2) + 𝑑𝑣(𝑡 − 𝜏1))

where

𝑓1(𝑥) = 𝑓2(𝑥) =

1
1 + 𝑒−𝛿𝑥

.

The parameter values are 𝜃𝑢 = 0.1, 𝜃𝑣 = 0.2, 𝜏1 = 𝜏2 = 0.152, 𝑎 = 𝑑 = −19, 𝑏 = 𝑐 = 10, and

𝛿 = 10. We solve the system with initial conditions (𝑢0, 𝑣0) = (0.05, 0.05) and 𝑑𝑡 = 0.001 up to

𝑡 = 50. SW is done using 𝑓 = 𝑢, 𝑑 = 4, 𝜏 = 1.712, and 𝜖 = 0.01.

5.2.5 Electromagnetic Radiation on a Wilson Neuron Model

Wilson introduced a simplified model for a neocortical neuron by making assumptions on the

Hodgkin-Huxley model [48]. The biophysics of these neurons is governed by the interplay of about

a dozen ion currents. His model showed the need of only four ion currents to accurately replicate

known spiking behavior. We analyze an extension of his model that incorporates the presence of

electromagnetic radiation (EMR).

Commonplace presence of electronic devices introduces EMR exposure to neurons. The effects

EMR has was imitated with the presence of a flux-controlled memristor in [49]. Their proposed

model was shown to exhibit quasiperiodicity and their results were successfully replicated using

62

Figure 5.5 Top row: Depiction of the phase space in (𝑣, 𝑟, 𝜙). The solution 𝑟 (𝑡) we used for SW is
plotted in the middle followed by its fast fourier transform. Bottom row: Persistence diagrams.

a micro-controller unit based hardware platform. Their mathematical model describe membrane

potential 𝑣, recovery variable 𝑟, and inner state of the memristor 𝜙

𝐶𝑚

𝑑𝑣
𝑑𝑡

= −𝑚∞(𝑣)(𝑣 − 𝐸𝑁𝑎) − 𝑔𝐾𝑟 (𝑣 − 𝐸𝐾) + 𝐼𝑒𝑥𝑡 − 𝑘1𝑊 (𝜙)𝑣,

𝑑𝑟
𝑑𝑡

=

1
𝜏𝑟

(−𝑟 + 𝑟∞(𝑣)),

𝑑𝜙
𝑑𝑡

= 𝑣 − 𝑘2𝜙 + 𝜙𝑒𝑥𝑡

where

and

𝑚∞(𝑣) = 17.8 + 47.6𝑣 + 33.8𝑣2

𝑟∞(𝑣) = 1.24 + 3.7𝑣 + 3.2𝑣2.

We use typical model parameters for the membrane capacitor, 𝐶𝑚 = 1, reverse potential for sodium

and potassium, 𝐸𝑁𝑎 = 0.5 and 𝐸𝐾 = −0.95, respectively, maximal conductance of potassium,

𝑔𝐾 = 26, potassium ion channel activation time constant, 𝜏𝑟 = 5, and external stimulus current

𝐼𝑒𝑥𝑡 = 1 [50]. The EMR external contribution is given by 𝜙𝑒𝑥𝑡 = 𝐴 sin(2𝜋𝐹𝑡), the terms 𝑘1𝑊 (𝜙)𝑣

63

and 𝑘2𝜙 denote the induction current caused by variation of magnetic flux and the leakage of

magnetic flux. The memductance of the memristor is given by 𝑊 (𝜙) = 𝑎 − 𝑏 tanh 𝜙 [51]. The

remaining parameter values are 𝑎 = 1, 𝑏 = 3, 𝑘1 = 1.2, 𝑘2 = 0.5, 𝐴 = 0.35, and 𝐹 = 0.22. We

solve the system using the initial conditions (𝑣0, 𝑟0, 𝜙0) = (0, 0, 0) and 𝑑𝑡 = 0.01 up to 𝑡 = 500.

SW was done with 𝑓 = 𝑟, 𝑑 = 4, 𝜏 = 72.21, and 𝜖 = 0.07, see Figure 5.5.

5.2.6 Competitive Threshold-Linear Network

Threshold-linear networks (TLNs) provide an accessible model acting under the presence of a

threshold nonlinearity [ref]. It provides a refinement to linear approximations of networks. The

latter can be conceptualized as a directed graph, 𝐺, in which nodes interact. The TLN model can

be described by the equations

𝑑𝑥𝑖
𝑑𝑡

= −𝑥𝑖 +

𝑛
∑︁

𝑗=1








𝑊𝑖 𝑗 𝑥 𝑗 + 𝑏𝑖






 +

Figure 5.6 Top row: PCA representation in three dimensions of the competitive TLN solutions.
The solution 𝑥4(𝑡) we used for SW is plotted in the middle followed by its fast fourier transform.
Bottom row: Persistence diagrams.

64

where 𝑖 = 1, . . . , 𝑛 and 𝑥𝑖 represents the activity level of node 𝑖. The term 𝑊𝑖, 𝑗 is the directed

connection strength, 𝑏𝑖 is the external drive, and [𝑦]+ = max{𝑦, 0} is the threshold nonlinearity.

We consider competitive TLNs, which are the case 𝑊𝑖, 𝑗 ≤ 0, 𝑊𝑖,𝑖 = 0, and 𝑏𝑖 ≥ 1. Furthermore,

we consider the connection strength

𝑊𝑖 𝑗 =






0

if 𝑖 = 𝑗,

−1 + 𝜆

if node 𝑗 is connected to node 𝑖,

−1 − 𝛿

else.

This model has shown complex behavior, including quasiperiodicity [ref]. We replicate the

quasiperiodic behavior of Figure 2 in [ref]. We refer to it for the initial condition and connec-

tion matrix. The rest of the parameter values are 𝑛 = 25, 𝑏𝑖 = 1, 𝜆 = 0.25, and 𝛿 = 0.5.. We solve

up to 𝑡 = 600 and perform SW with 𝑓 = 𝑥4, 𝑑 = 4, 𝜏 = 27.53, and 𝜖 = 0.07, see Figure 5.6.

5.2.7 Restricted Three Body Problem

In celestial mechanics, the restricted three body problem (RTBP) offers an accessible model with

known equilibrium points [52]. The equations describe the progression of three celestial bodies in

which one of them is considered a massless particle. The other two bodies, denotes as primaries,

are assumed to move in circular orbits around their center of mass. By taking a coordinate system

that rotates with the primaries and under the appropriate scaling, it can be assumed the primaries

have masses 1 − 𝜇 and 𝜇, 𝜇 ∈ [0, 1/2], are fixed at (𝜇, 0, 0) and (𝜇 − 1, 0, 0), respectively, and

complete one revolution in 2𝜋 [41]. This framework allows us to express the motion of the massless

particle by the equations

(cid:165)𝑥 − 2 (cid:164)𝑦 = Ω𝑥,

(cid:165)𝑦 + 2 (cid:164)𝑥 = Ω𝑦,

(cid:165)𝑧 = Ω𝑧,

where

Ω =

1
2

(𝑥2 + 𝑦2) +

1 − 𝜇
𝑟1

+

𝜇
𝑟2

+

1
2

𝜇(1 − 𝜇),

65

Figure 5.7 Top row: Phase space depiction of RTBP in (𝑥, 𝑦, 𝑧) space. We include the Earth, shown
in blue and green, and the Moon, shown in grey. The solution 𝑥(𝑡) we used for SW is plotted in the
middle followed by its fast fourier transform. Bottom row: Persistence diagrams.

and

√︃

𝑟1 =

(𝑥 − 𝜇)2 + 𝑦2 + 𝑧2,

√︃

𝑟2 =

(𝑥 − 𝜇 + 1)2 + 𝑦2 + 𝑧2

are the distances from the particles to the primaries. The case of the Earth-Moon system is of

particular interest since it can aid in spacecraft missions interested in the Sun and the magnetosphere

of the Earth [41]. Indeed, near the equilibrium points of the system, quasiperiodicity is present

which translates to nice trajectories for a mission. We replicate this behavior for the Earth-Moon

system, see Figure 5.7. In this case, 𝜇 = 0.0121506. We solve the system with the initial condition

(𝑥0, 𝑦0, 𝑧0, (cid:164)𝑥, (cid:164)𝑦, (cid:164)𝑧) = (−0.5, 0, 0, 0, 0, 0.73) up to 𝑡 = 100. SW was done with 𝑓 = 𝑥, 𝑑 = 4, 𝜏 = 4.37

and 𝜖 = 0.03.

5.2.8 Bicircular Restricted Four Body Problem

Considering the effects of the sun in the RTBP model gives rise to the bicircular restricted four

body problem (BCR4BP). The model is of importance for exploiting the force from the sun in

trajectory designs for lunar missions and has found applications for ballistic lunar transfers to the

lunar region [53]. The new model uses the same coordinate axes as the RTBP and the same circular

66

Figure 5.8 Top row: Phase space depiction of BCR4BP in (𝑥, 𝑦, 𝑧) space. We include the Earth,
shown in blue and green, the Moon, shown in grey, and the Sun, shown in red and yellow. We note
the position of the Earth and Sun are not to scale. The solution 𝑧(𝑡) we used for SW is plotted in
the middle followed by its fast fourier transform. Bottom row: Persistence diagrams.

motion assumptions of the Earth and Moon but it now includes terms pertaining to the influence of

the Sun. The Sun is assumed to lie in the x-y plane at a fixed distance to the origin, 𝑎4, and moving
with a constant angular velocity, (cid:164)𝜃𝑆, in circular motion. Furthermore, the model assumes the Earth

and Moon are not perturbed by solar gravity. Under these assumptions, the equations are [53]

(cid:165)𝑥 = 2 (cid:164)𝑦 +

(cid:165)𝑦 = −2 (cid:164)𝑥 +

,

𝜕Υ
𝜕𝑥
𝜕Υ
𝜕𝑦

,

(cid:165)𝑧 =

𝜕Υ
𝜕𝑧

,

where

Υ =

1 − 𝜇
𝑟13

+

𝜇
𝑟23

+

𝑥2 + 𝑦2
2

+ 𝜆

(cid:32) 𝑚4
𝑟43

−

𝑚4
𝑎3
4

(cid:33)

(𝑥4𝑥 + 𝑦4𝑦 + 𝑧4𝑧)

.

and (𝑥4, 𝑦4, 𝑧4) denotes the position of the sun, 𝑟1,3, 𝑟2,3, 𝑟4,3 the distance of the Earth to the

massless object, the Moon to the massless object, and the Sun to the massless object, respectively,

67

and 𝑚4 is the non-dimensional mass of the Sun. We note the case 𝜆 = 0 reduced to the RTBP.

A study of the system as 𝜆 increases to 1 is done in [53], where they showed the model exhibited

quasiperiodicty. We replicate the quasiperiodic behavior, see Figure 5.8, for parameter values

𝜆 = 1, 𝜇 = 0.012155, 𝑎4 = 388.84, 𝑚4 = 328950.69, and (cid:164)𝜃𝑆 = −0.9251986. The system was

solved with the initial condition (𝑥0, 𝑦0, 𝑧0, (cid:164)𝑥, (cid:164)𝑦, (cid:164)𝑧, 𝜃𝑆) = (1.09, 0, 0, 0, 0.19, 0.05, 2) up to 𝑡 = 200.

SW was done with 𝑓 = 𝑧, 𝑑 = 4, 𝜏 = 39.99 and 𝜖 = 0.009.

5.3 Concluding Remarks

Table 5.1 Running Times

Example
5.2.1
5.2.2
5.2.3
5.2.4
5.2.5
5.2.6
5.2.7
5.2.8

SW
7008.66 sec
4351.15 sec
3126.81 sec
5556.86 sec
5137.54 sec
5805.28 sec
6900.29 sec
6615.33 sec

SW(S1)
7672.85 sec
5351.15 sec
3218.73 sec
7628.48 sec
6505.54 sec
7193.97 sec
7937.03 sec
8231.25 sec

K(3G)
0.87 sec
0.42 sec
0.50 sec
0.98 sec
0.66 sec
0.22 sec
2.09 sec
0.59 sec

In this section we detailed how 3G can be implemented on general quasiperiodic functions.

Our approximation method was used on dynamical systems known to exhibit quasiperiodicity. As

our figures show, we successfully approximated the diagrams obtained from SW. Error bounds can

also be computed for our method. We illustrate them in the plots as rectangles. Furthermore,

our approach significantly reduces computational time as shown in Table 5.1. We hope our work

contributes to the implementation of the sliding window embedding technique on large data sets.

68

CHAPTER 6

CONTRIBUTION AND FUTURE WORK

We detailed how SW is a vital method used for the better understanding of quasiperiodic signals.

The correspondence between a quasiperiodic function of 𝑁 incommensurate frequencies and a

𝑁-dimensional torus make persistent homology an ideal component. However, the exponential

computational cost of the latter limit the datasets that can be impacted by SW. Our contribution is

providing an approximation to the persistence diagrams of interest (Section 1.1). Our method 3G

achieves this within known error bounds. Our computation, on average, took less than 1 second to

run for the examples depicted in Chapter 5. This contrast to an average computation time of 5,563

seconds when using a standard library, Ripser [16], to compute the persistence diagrams of interest

(Table 5.1).

Future work for our project would be to improve our approximation by tightening the error

bounds. A more granular approach would be to inspect the specific transformation matrix instead

of relying on a singular value argument. Another direction is to apply our work to real world data

such as electrocardiogram dataset.

69

BIBLIOGRAPHY

[1] H. R. Wilson and J. D. Cowan, “Excitatory and inhibitory interactions in localized populations

of model neurons,” Biophysical journal, vol. 12, no. 1, pp. 1–24, 1972.

[2]

J. A. Perea, “Topological time series analysis,” Notices of the American Mathematical Society,
vol. 66, no. 5, pp. 686–694, 2019.

[3] F. Takens, “Detecting strange attractors in turbulence,” in Dynamical systems and turbulence,

Warwick 1980, pp. 366–381, Springer, 1981.

[4] P. Skraba, V. De Silva, and M. Vejdemo-Johansson, “Topological analysis of recurrent sys-
tems,” in NIPS 2012 Workshop on Algebraic Topology and Machine Learning, December 8th,
Lake Tahoe, Nevada, pp. 1–5, 2012.

[5] H. Gakhar and J. A. Perea, “Sliding window persistence of quasiperiodic functions,” arXiv

preprint arXiv:2103.04540, 2021.

[6] C. J. Tralie and J. A. Perea, “(quasi) periodicity quantification in video data, using topology,”

SIAM Journal on Imaging Sciences, vol. 11, no. 2, pp. 1049–1077, 2018.

[7]

J. A. Perea, A. Deckard, S. B. Haase, and J. Harer, “Sw1pers: Sliding windows and 1-
persistence scoring; discovering periodicity in gene expression time series data,” BMC bioin-
formatics, vol. 16, no. 1, pp. 1–12, 2015.

[8] D. Levine and P. J. Steinhardt, “Quasicrystals: a new class of ordered structures,” Physical

review letters, vol. 53, no. 26, p. 2477, 1984.

[9]

I. Wilden, H. Herzel, G. Peters, and G. Tembrock, “Subharmonics, biphonation, and deter-
ministic chaos in mammal vocalization,” Bioacoustics, vol. 9, no. 3, pp. 171–196, 1998.

[10] M. Belloy, M. Naeyaert, G. Keliris, A. Abbas, S. Keilholz, A. Van Der Linden, and M. Verhoye,
“Dynamic resting state fmri in mice: detection of quasi-periodic patterns,” Proceeding of the
International Soc. Magn. Reson. Med, vol. 961, 2017.

[11] S. Noorzadeh, M. Niknazar, B. Rivet, J. Fontecave-Jallon, P.-Y. Guméry, and C. Jutten, “Mod-
eling quasi-periodic signals by a non-parametric model: Application on fetal ecg extraction,”
in 2014 36th Annual International Conference of the IEEE Engineering in Medicine and
Biology Society, pp. 1889–1892, IEEE, 2014.

[12] D. Weixing, H. Wei, W. Xiaodong, and C. Yu, “Quasiperiodic transition to chaos in a plasma,”

Physical review letters, vol. 70, no. 2, p. 170, 1993.

[13] M. A. Quiroz-Juárez, R. Vázquez-Medina, E. Ryzhii, M. Ryzhii, and J. L. Aragón, “Quasiperi-
odicity route to chaos in cardiac conduction model,” Communications in Nonlinear Science

70

and Numerical Simulation, vol. 42, pp. 370–378, 2017.

[14] H. Gakhar and J. A. Perea, “Künneth formulae in persistent homology,” arXiv preprint

arXiv:1910.05656, 2019.

[15] J.-D. Boissonnat, T. K. Dey, and C. Maria, “A space and time efficient implementation for

computing persistent homology,” tech. rep., INRIA Research Report 8195, 2012.

[16] U. Bauer, Ripser, vol. 514. URL: https://github. com/Ripser/ripser, 2016.

[17] N. O. Malott, S. Chen, and P. A. Wilsey, “A survey on the high-performance computation
of persistent homology,” IEEE Transactions on Knowledge and Data Engineering, vol. 35,
no. 5, pp. 4466–4484, 2022.

[18] C. Olds, “Continued fractions,” Mathematical Association of America, 1963.

[19] H. Hancock, “Development of the minkowski geometry of numbers,” (No Title), 1939.

[20] T. Van Ravenstein, “The three gap theorem (steinhaus conjecture),” Journal of the Australian

Mathematical Society, vol. 45, no. 3, pp. 360–370, 1988.

[21] M. M. Murudi and S. M. Mane, “Seismic effectiveness of tuned mass damper (tmd) for
different ground motion parameters,” in 13th World Conference on Earthquake Engineering,
vol. 2, pp. 1–8, 2004.

[22] R. Wen, T. Li, and B. Zhen, “Quasi-periodic motions of a pendulum with vibrating suspension

point,” Journal of Vibration Engineering & Technologies, vol. 7, pp. 519–532, 2019.

[23] Phys.org, “Japanese companies look to quake damping pendulums,” August 2013. Accessed:

2024-03-20.

[24] A. Hatcher, Algebraic topology. Cambridge: Cambridge University Press, 2002.

[25] E. J. Hanson and J. D. Rock, “Decomposition of pointwise finite-dimensional sˆ 1 persistence

modules,” arXiv preprint arXiv:2006.13793, 2020.

[26] Y. Cao, A. Monod, A. Vlontzos, L. Schmidtke, and B. Kainz, “Topological information
retrieval with dilation-invariant bottleneck comparative measures,” Information and Inference:
A Journal of the IMA, vol. 12, no. 3, pp. 1964–1996, 2023.

[27] J. D. Cowan, J. Neuman, and W. van Drongelen, “Wilson–cowan equations for neocortical
dynamics,” The Journal of Mathematical Neuroscience, vol. 6, no. 1, pp. 1–24, 2016.

[28] G. F. Poggio and L. J. Viernstein, “Time series analysis of impulse sequences of thalamic

somatic sensory neurons,” Journal of Neurophysiology, vol. 27, no. 4, pp. 517–545, 1964.

71

[29] E. Robert, “Fourier series: a modern introduction,” Springer Science & Business Media,

vol. 2, 2012.

[30] C. F. Gauss, Demonstratio nova theorematis omnem functionem algebraicam. apvd CG

Fleckeisen, 1799.

[31] A. Zhigljavsky and I. Aliev, “Kronecker sequences: Asymptotic distributions of the partition

lengths,”

[32] N. Leong, Sums of Reciprocals and the Three Distance Theorem. PhD thesis, University of

York, 2017.

[33] L. S. Salas and J. A. Perea, “Estimation of persistence diagrams via the three gap theorem.”

Manuscript under review, 2023.

[34] G. Lochs, “Vergleich der genauigkeit von dezimalbruch und kettenbruch,” in Abhandlungen
aus dem Mathematischen Seminar der Universität Hamburg, vol. 27, pp. 142–144, Springer,
1964.

[35] C. Faivre, “On decimal and continued fraction expansions of a real number,” Acta Arithmetica,

vol. 82, no. 2, pp. 119–128, 1997.

[36] N. Saul and C. Tralie, “Persim,” 2019.

[37] D. Cohen-Steiner, H. Edelsbrunner, and J. Harer, “Stability of persistence diagrams,” in
Proceedings of the twenty-first annual symposium on Computational geometry, pp. 263–271,
2005.

[38] F. Chazal, V. De Silva, M. Glisse, and S. Oudot, The structure and stability of persistence

modules, vol. 10. Springer, 2016.

[39] T. Detroux, G. Habib, L. Masset, and G. Kerschen, “Performance, robustness and sensi-
tivity analysis of the nonlinear tuned vibration absorber,” Mechanical Systems and Signal
Processing, vol. 60, pp. 799–809, 2015.

[40] E. Kaslik, E.-A. Kokovics, and A. Rădulescu, “Stability and bifurcations in wilson-cowan
systems with distributed delays, and an application to basal ganglia interactions,” Communi-
cations in Nonlinear Science and Numerical Simulation, vol. 104, p. 105984, 2022.

[41] G. Gómez and J. M. Mondelo, “The dynamics around the collinear equilibrium points of the

rtbp,” Physica D: Nonlinear Phenomena, vol. 157, no. 4, pp. 283–321, 2001.

[42] C. E. Shannon, “Communication in the presence of noise,” Proceedings of the IRE, vol. 37,

no. 1, pp. 10–21, 1949.

72

[43] S. C. Shadden, F. Lekien, and J. E. Marsden, “Definition and properties of lagrangian coherent
structures from finite-time lyapunov exponents in two-dimensional aperiodic flows,” Physica
D: Nonlinear Phenomena, vol. 212, no. 3-4, pp. 271–304, 2005.

[44] G. D. Charó, G. Artana, and D. Sciamarella, “Topology of dynamical reconstructions from

lagrangian data,” Physica D: Nonlinear Phenomena, vol. 405, p. 132371, 2020.

[45] J. Penalva Vadell, “Takens’ theorem: Proof and applications,” 2018.

[46] S. Coombes and C. Laing, “Delays in activity-based neural networks,” Philosophical Trans-
actions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 367,
no. 1891, pp. 1117–1129, 2009.

[47] A. J. N. Holgado, J. R. Terry, and R. Bogacz, “Conditions for the generation of beta oscillations
in the subthalamic nucleus–globus pallidus network,” Journal of Neuroscience, vol. 30, no. 37,
pp. 12340–12352, 2010.

[48] H. R. Wilson, “Simplified dynamics of human and mammalian neocortical neurons,” Journal

of theoretical biology, vol. 200, no. 4, pp. 375–388, 1999.

[49] Z. Ju, Y. Lin, B. Chen, H. Wu, M. Chen, and Q. Xu, “Electromagnetic radiation induced non-
chaotic behaviors in a wilson neuron model,” Chinese Journal of Physics, vol. 77, pp. 214–222,
2022.

[50] Q. Xu, Z. Ju, C. Feng, H. Wu, and M. Chen, “Analogy circuit synthesis and dynamics
confirmation of a bipolar pulse current-forced 2d wilson neuron model,” The European
Physical Journal Special Topics, vol. 230, no. 7, pp. 1989–1997, 2021.

[51] B. Bao, H. Qian, Q. Xu, M. Chen, J. Wang, and Y. Yu, “Coexisting behaviors of asym-
metric attractors in hyperbolic-type memristor based hopfield neural network,” Frontiers in
Computational Neuroscience, vol. 11, p. 81, 2017.

[52] W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, “Heteroclinic connections between
periodic orbits and resonance transitions in celestial mechanics,” Chaos: An Interdisciplinary
Journal of Nonlinear Science, vol. 10, no. 2, pp. 427–469, 2000.

[53] B. P. McCarthy and K. C. Howell, “Four-body cislunar quasi-periodic orbits and their applica-
tion to ballistic lunar transfer design,” Advances in Space Research, vol. 71, no. 1, pp. 556–584,
2023.

73