# Reservoir Computing with Dynamical Systems

Allen G Hart

###### Abstract

A reservoir computer is a special type of neural network, where most of the
weights are randomly fixed and only a subset are trained.

In this thesis we prove results about reservoir computers trained on
deterministic dynamical systems, and stochastic processes. We focus mostly on
a special type of reservoir computer called an Echo State Network (ESN).

In the deterministic case, we prove (under some assumptions) that if a
reservoir computer has the Echo State Property (ESP), then there is a $C^{1}$
generalised synchronisation between the input dynamical system and the
dynamics in the reservoir space. Furthermore, we prove that a reservoir
computer with the local ESP in several disjoint subsets of the reservoir space
will admit several distinct generalised synchronisations. In the special case
that the reservoir map is linear, and has the ESP, we prove that the
generalised synchronisation is generically an embedding. This result admits
Takens’ embedding Theorem as a special case.

We go to show that ESNs trained on scalar observations of an ergodic dynamical
system can approximate an arbitrary target function, including the next step
map used in time series forecasting. This universal approximation property
holds despite the training process being entirely linear.

We prove analogous results for ESNs trained on observations of a stochastic
process, which are not be Markovian in general. We use these results to
develop supervised learning, and reinforcement learning algorithms supported
by an ESN.

In the penultimate chapter of this thesis, we use a reservoir computer to
numerically solve linear PDEs. In the final chapter, we conclude and discuss
directions for future work.

## Acknowledgements

Thank you everyone who supported me over the last four years, and for all the
conversations about reservoir computing and dynamical systems.

## List of Publications and Preprints

This thesis is largely derived from the following publications and preprints.

  1. 1.

Hart, A., Hook, J. & Dawes, J. H. P. (2020), ‘Embedding and approximation
theorems for echo state networks’, Neural Networks 128, 234–247.

Results from this paper appear in Chapters 2, 3, and 5. A G Hart is the lead
author with J L Hook and J H P Dawes having supervisory roles.

  2. 2.

Grigoryeva, L., Hart, A. & Ortega, J.-P. (2021), ‘Chaos on compact manifolds:
Differentiable synchronizations beyond the Takens theorem’, Phys. Rev. E103,
062204.

The paper forms the basis of Chapter 3. All authors contributed equally.

  3. 3.

Grigoryeva, L., Hart, A. & Ortega, J.-P. (2021), ‘Learning Strange Attractors
with reservoir maps’, _arXiv:2108.05024_. [PREPRINT]

This paper forms the basis of Chapter 4. All authors contributed equally.

  4. 4.

Hart, A. G., Hook, J. L. & Dawes, J. H. P. (2021), ‘Echo state networks
trained by Tikhonov least squares are $L^{2}(\mu)$ approximators of ergodic
dynamical systems’, Physica D: Nonlinear Phenomena p. 132882.

Results from this paper appear in Chapter 5. A G Hart is the lead author with
J L Hook and J H P Dawes having supervisory roles.

  5. 5.

Hart, A. G., Olding, K. R., Cox, A. M. G., Isupova, O. & Dawes, J. H. P.
(2021), ‘Using Echo State Networks to Approximate Value Functions for
Control’, _arXiv:2004.08170_. [PREPRINT]

This paper forms the basis of Chapter 6. A G Hart was primarily responsible
for the ESN approximation results, K R Olding was primarily responsible for
derivations surrounding the market making problem, while J H P Dawes, A M G
Cox, and O Isupova had supervisory roles.

###### Contents

  1. 1 Introduction
    1. 1.1 Preamble
    2. 1.2 ESNs and dynamical systems
    3. 1.3 Training an ESN on the Lorenz system
      1. 1.3.1 Creating the reservoir states
      2. 1.3.2 Learning the targets $\zeta$ from the observations $\xi$
      3. 1.3.3 Predicting the future
    4. 1.4 The remainder of the thesis
  2. 2 Computational Methods
    1. 2.1 Fixed Points
    2. 2.2 Lyapunov Spectra
    3. 2.3 Persistent homology
      1. 2.3.1 Introduction
      2. 2.3.2 Theory
      3. 2.3.3 Persistent homology of the Lorenz attractor
  3. 3 Generalised Synchronisation (GS)
    1. 3.1 Background
    2. 3.2 Echo State Property
    3. 3.3 A theorem for global GS
    4. 3.4 Motivating example
    5. 3.5 Theorems for local GS
    6. 3.6 Curious examples
    7. 3.7 Relationship between GS and the Echo Index
    8. 3.8 Biological interpretation
  4. 4 Embeddings
    1. 4.1 Background
    2. 4.2 Summary of novel results
    3. 4.3 An embedding allows for learning
    4. 4.4 Novel results in detail
      1. 4.4.1 Existence of $C^{2}$ GS
      2. 4.4.2 Immersive GS
      3. 4.4.3 Embedding GS
      4. 4.4.4 An embedding almost surely
  5. 5 Universal Approximation
    1. 5.1 Universal Approximation Theorems
    2. 5.2 Ergodic Theory
    3. 5.3 A Training Theorem for Echo State Networks
      1. 5.3.1 Regularised Least Squares Regression
      2. 5.3.2 Linear Universal Approximators
      3. 5.3.3 Linear Universal Approximators Trained by Least Squares
      4. 5.3.4 Convergence rate of the time average to the space average
    4. 5.4 The Lorenz attractor is stably mixing
    5. 5.5 Online Learning
      1. 5.5.1 Fixed learning parameter $\alpha$
      2. 5.5.2 Learning rate $\alpha_{k}=1/k$
    6. 5.6 Value Functions and Control
    7. 5.7 Biological interpretation
    8. 5.8 Learning Diffeomorphic Dynamics
  6. 6 Stochastic Dynamical Systems and Reinforcement Learning
    1. 6.1 Filters and Functionals
    2. 6.2 System isomorphism
    3. 6.3 Supervised learning on stationary ergodic processes
    4. 6.4 Reinforcement learning on stationary ergodic processes
    5. 6.5 Online learning on stationary ergodic processes
  7. 7 Partial Differential Equations
    1. 7.1 Introduction
    2. 7.2 Theory
    3. 7.3 Laplace’s equation on the disc
      1. 7.3.1 Analytic solution
      2. 7.3.2 Reservoir computer solution
    4. 7.4 Future Work and Open Questions
  8. 8 Conclusions and Future Work

## Nomenclature

  * Stochastic Dynamical Systems
  * $(\mathbb{R}^{N})^{\mathbb{Z}}$

The set of bi-infinite $\mathbb{R}^{N}$ valued sequences, i.e. the set of maps
$\mathbb{Z}\mapsto\mathbb{R}^{N}$

  * $\gamma\in[0,1)$

Discount factor

  * $\mathcal{R}:(\mathbb{R}^{d})^{\mathbb{Z}}\to\mathbb{R}$

Reward functional

  * $T:(\mathbb{R}^{N})^{\mathbb{Z}}\to(\mathbb{R}^{N})^{\mathbb{Z}}$

Time shift operator $T(z)_{k}:=z_{k+1}$

  * $V:(\mathbb{R}^{d})^{\mathbb{Z}}\to\mathbb{R}$

Value functional

  * Deterministic Dynamical Systems
  * $\omega:M\to\mathbb{R}^{d}$

Observation function

  * $\phi:M\to M$

Evolution operator on $M$

  * $\text{Diff}^{n}(M)$

Diffeomorphism of order $n\in\mathbb{N}$ on $M$

  * $C^{n}(M,\mathbb{R}^{N})$

The $n$-times differentiable maps from $M$ to $\mathbb{R}^{N}$

  * $D\omega(m):T_{m}\to\mathbb{R}^{d}$

The differential of $\omega$ at $m\in M$

  * $m\in M$

Point on $M$

  * $M$

$q$-dimensional manifold

  * $q\in\mathbb{N}$

Dimension of the manifold $M$

  * $T_{m}$

Tangent space at $m\in M$

  * $T_{m}\phi:T_{m}\to T_{\phi(m)}$

Tangent map for $\phi$ at $m\in M$.

  * $u:M\to\mathbb{R}$

Target function

  * Echo State Networks and reservoir maps
  * $\ell\in\mathbb{N}$

Number of sample points

  * $\mathcal{F}$

Class of linear universal approximators

  * $\sigma:\mathbb{R}^{N}\to\mathbb{R}^{N}$

Activation function

  * $A\in\mathbb{M}_{N\times N}(\mathbb{R})$

Reservoir matrix

  * $b\in\mathbb{R}^{N}$

Bias vector

  * $C\in\mathbb{M}_{N\times d}(\mathbb{R})$

Input matrix

  * $d\in\mathbb{N}$

Number of input channels / dimension of the input

  * $F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$

reservoir map

  * $f_{(\phi,\omega,F)}:M\to\mathbb{R}^{N}$

Generalised Synchronisation

  * $N\in\mathbb{N}$

Number of neurons / dimension of the reservoir

  * $W\in\mathbb{R}^{N}$

Output vector / linear readout layer

  * $x_{k}\in\mathbb{R}^{N}$

Reservoir state at time $k\in\mathbb{Z}$

  * $z_{k}\in\mathbb{R}^{d}$

Input at time $k\in\mathbb{Z}$

  * Partial Differential Equations
  * $\Delta$

The Laplace operator $\Delta:C^{2}(\mathbb{R}^{d},\mathbb{R})\to
C^{0}(\mathbb{R}^{d},\mathbb{R})$

  * $\phi:\Omega\to\mathbb{R}$

Solution to the boundary value problem

  * $h:\partial\Omega\to\mathbb{R}$

Boundary data

  * $Z\in\Omega$

Random interior point

  * $Z^{\prime}\in\partial\Omega$

Random boundary point

##  Chapter 1 Introduction

###  1.1 Preamble

A single layer feedforward neural network $f:\mathbb{R}^{d}\to\mathbb{R}$ is a
map the form

$\displaystyle f(z)=\sum_{k=1}^{N}W_{k}\sigma(C_{k}^{\top}z+b_{k})$

paramatrised by weights $C_{k}\in\mathbb{R}^{d}$, $W_{k}\in\mathbb{R}$ and
biases $b_{k}\in\mathbb{R}$ (Higham & Higham, 2019). The map
$\sigma:\mathbb{R}\to\mathbb{R}$ is called the activation function. Usual
choices of the activation function include

  * •

the rectified linear unit (ReLU) defined coordinate wise
$\sigma(z)_{i}=\max(z_{i},0)$

  * •

the hyperbolic tangent also defined coordinate wise
$\sigma(z)_{i}=\tanh(z_{i})$.

A multilayer (or deep) neural network can be created by composing several
feedforward layers together. A composition of $n$ feedforward neural networks
is called an $n$-layer neural network (Higham & Higham, 2019). Both deep and
shallow (single-layer) neural networks are dense in some appropriate function
space so that for a function $g$ in the appropriate space, and given
sufficiently many neurons $N$, we can choose weights $W_{k},C_{k}$ and biases
$b_{k}$ such that the feedforward neural network $f$ approximates $g$. This is
made formal in the universal approximation theorem for feedforward neural
networks (discussed in Chapter 5).

We can use feedforward neural networks to solve supervised learning problems.
These problems involve a supervisor, an agent and a collection of data. The
data are labelled by the supervisor, and the agent is tasked with learning the
relationship between the data and its labels, in such a way that the
relationship can be generalised to unlabelled data. Learning this relationship
is called training. The success of the training is measured by how accurately
the agent can correctly label a set of unlabelled data.

There are 2 major categories of supervised learning: regression and
classification. In the former, the agent seeks a map from the data
(observations) to the labels (targets) such that any unlabelled datum
(observation) on a continuum can be assigned a label (target) also on a
continuum. The problem of time series forecasting is often formulated as a
regression problem. Classification is similar, except that the labels do not
lie on a continuum but instead lie in a finite set of categories. A
classification algorithm may for example classify a set of images as either
being cats or dogs.

When we use a feedforward neural network to solve a supervised learning
problem, the data are denoted $\\{z_{i}\in\mathbb{R}^{d}\\}_{i\in I}$ for some
finite index $I$. We assume the existence of a relationship
$g:\mathbb{R}^{d}\to\mathbb{R}$ from the data to the targets
$\\{g(z_{i})\\}_{i\in I}$. Then our goal is to find weights $W_{k},C_{k}$ and
biases $b_{k}$ such that the feedforward neural network $f$ approximates $g$
by finding $W_{k},C_{k},b_{k}$ that minimise the so called loss $\mathcal{L}$
(sometimes) defined

$\displaystyle\mathcal{L}=\sum_{i=1}^{|I|}\lVert f(z_{i})-g(z_{i})\rVert^{2}.$

We can view the loss $\mathcal{L}$ as a smooth function of the weights
$W_{k},C_{k}$ and biases $b_{k}$, so we can try to minimise the loss with
smooth optimisation techniques like gradient descent or stochastic gradient
descent. This optimisation problem is generally non-convex, so finding a
global minimum is generally difficult.

Many supervised learning problems involve data drawn from a time series where
the data has an important temporal structure. Such problems encompass speech
recognition, the forecasting of chaotic systems, and decision making in
robotics. For such problems, it is often better to replace the feedforward
neural network with a _reservoir computer_.

In this thesis we focus on mathematical results that hold for reservoir
computers. We focus on reservoir computers trained on a discrete, equally
spaced time series that is either:

  1. 1.

A (partial) observation of a system of ODEs, or

  2. 2.

A realisation of a stochastic process.

A reservoir computer is a pair of maps $(F,h)$ where
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ is called the
reservoir map and takes as input a reservoir state vector $x\in\mathbb{R}^{N}$
and observation vector $z\in\mathbb{R}^{d}$, then returns a new reservoir
state vector $x\in\mathbb{R}^{N}$. Given an initial state
$x_{0}\in\mathbb{R}^{N}$ and a time series of observations
$z_{1},z_{2},\ldots\in\mathbb{R}^{d}$ the reservoir computer creates a
sequence of states $x_{1},x_{2},\ldots\in\mathbb{R}^{N}$ like so

$\displaystyle x_{k+1}=F(x_{k},z_{k}).$

The state $x_{k+1}$ then depends on all observations
$\ldots,z_{k-2},z_{k-1},z_{k}$ prior to and including timestep $k$. The second
part of the reservoir computer is a map $h:\mathbb{R}^{N}\to\mathbb{R}^{s}$
which takes a reservoir state $x\in\mathbb{R}^{N}$ and returns an output. The
map $h$ could be a feedforward neural network or simply a row vector
$W^{\top}$. The maps $F$ and $h$ together form a reservoir computer, which
takes a time series $z_{1},z_{2},\ldots\in\mathbb{R}^{d}$ of observations and
computes an output. We can define a reservoir computer formally as follows:

######  Definition 1.1.1.

(Reservoir Computer) We call
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ a reservoir map,
$h:\mathbb{R}^{N}\to\mathbb{R}^{s}$ a readout map, and the pair $(F,h)$ a
reservoir computer.

Reservoir computing (RC) exploded in popularity at the start of the early
2000s after the seminal papers of Jaeger2001 and
doi:10.1162/089976602760407955. RC has since been applied to both signal
processing and machine learning, and has captured the interest of researchers
in computer science and robotics (10.1007/978-3-540-25940-4_14), physics
(Inubushi2017) (chaos_on_compacta), mathematics (GRIGORYEVA2018495),
(CENI2020132609), and electrical engineering (TANAKA2019100). Furthermore,
there is some biology and neuroscience literature suggesting that a sensory
input stimulating a nervous system can be described as a reservoir computer
driven by a time series (izhikevich2007dynamical), (TANAKA2019100). Rather
ambitious applications of reservoir computing include “building by 2050, a
team of fully autonomous humanoid robots to beat the human winning team of the
FIFA Soccer World Cup” (10.1007/978-3-540-25940-4_14), predicting magnetic
storms (kataokareconstructing), and forecasting electricity sales
(li2020research).

Reservoir computing can be performed on an ordinary digital computer, which is
often referred to reservoir computing _in silico_. Alternatively, reservoir
computing can be performed on purpose build hardware; often known as reservoir
computing _in materio_. Reservoirs have been physically built from magnets,
photonic node arrays (TANAKA2019100), and field programmable gate arrays
(patent) (Apostel2021), and (possibly) occur naturally in the form of
biological brains. In this thesis, all reservoir computing experiments are
performed in silico, though the mathematical results hold more generally. Some
authors have argued that there is potential to develop extremely efficient
reservoir computers in materio, which could eventually outperform ordinary
computers (9106624).

The goal in this thesis is to further develop a rigorous mathematical
description of reservoir computing on time series, which we hope will serve
the highly multidisciplinary reservoir computing community.

###  1.2 ESNs and dynamical systems

In this section we will go through a numerical experiment which was undertaken
at the start of the PhD. The outcome of the experiment is illuminating, and
motivates many of the results that appear in subsequent chapters in this
thesis.

Suppose we have a dynamical system, but we do not know its governing
equations, and instead have access only to a low dimensional observation of a
trajectory. We could have for instance the celebrated Lorenz system
(doi:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2)

$\displaystyle\dot{\xi}$ $\displaystyle=10(\upsilon-\xi)$
$\displaystyle\dot{\upsilon}$ $\displaystyle=\xi(28-\zeta)-\upsilon$ (1.1)
$\displaystyle\dot{\zeta}$ $\displaystyle=\xi\upsilon-(8/3)\zeta$

(where we have fixed the parameters to their ‘usual values’) and have access
to only the $\xi$ component, at regularly spaced time points
$t_{0},t_{1},\ldots t_{\ell-1}$. Given this finite time series, the key
question is: can we predict the future trajectory without knowledge of the
underlying equations? We note that the observation at time $t_{\ell}$ is not
determined from the observation at time $t_{\ell-1}$ alone, so we will need to
use many observations from the past $\ldots,t_{\ell-2},t_{\ell-1}$ to have any
hope of estimating the observation at $t_{\ell}$.

We approach this time series forecasting problem with a reservoir computer
called an Echo State Network (ESN). An ESN takes a time series of observations
$z_{0},z_{1},\ldots\in\mathbb{R}^{d}$ and produces a sequence of reservoir
state vectors $x_{0},x_{1},\ldots\in\mathbb{R}^{N}$ defined like so:

$\displaystyle x_{k+1}=F(x_{k},z_{k}):=\sigma(Ax_{k}+Cz_{k}+b)$ (1.2)

where

  * •

$\sigma:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is the activation function defined to
act coordinate-wise: $\sigma_{i}(x)=\tanh(x_{i})$ for $i=1,\ldots N$

  * •

$A\in\mathbb{M}_{N\times N}(\mathbb{R})$ is the $N\times N$ real reservoir
matrix

  * •

$C\in\mathbb{M}_{N\times d}(\mathbb{R})$ is the $N\times d$ real input matrix

  * •

$b\in\mathbb{R}^{N}$ is the real bias vector.

We will sometimes call the Lorenz system the _drive system_ and the dynamics
of the reservoir states the _response system_.

The recursive nature of (1.2) ensures that the reservoir state $x_{k+1}$
depends on all past observations $\ldots,z_{k-1},z_{k}$. The ESN therefore has
a sort of memory, because its representation of the world at time step $k$ is
$x_{k}$ which depends on everything it has previously seen. Equation (1.2) is
sometimes called a single layer (or shallow) ESN in contrast to multilayer (or
deep) ESNs studied by gallicchio2018design and others.

It seems natural that observations far in the past should have less impact on
the reservoir state $x_{k}$ than those observations in the near past. This is
called the Fading Memory Property (FMP) defined in GRIGORYEVA2018495. The FMP
is closely related to the Echo State Property (ESP) presented by Jaeger2001
and is introduced and discussed in Chapters 3 and 6. We can ensure that an ESN
has the ESP by demanding that the activation function $\sigma$ and reservoir
matrix $A$ are contracting. The matrix $A$ is called the reservoir matrix (and
the field called _reservoir computing_) because the entries of $A$ are fixed
and represent a sort of computational reservoir that is capable of encoding a
rich variety of inputs with a reservoir state vector $x$.

ESNs are often used to approximate the relationship between a time series of
observations $z_{0},z_{1},\ldots\in\mathbb{R}^{d}$ and a time series of
targets $u_{0},u_{1},\ldots\in\mathbb{R}$. In the special case that $d=1$ and
the targets $u_{k}=z_{k}$ then the problem is time series forecasting. There
is substantial numerical evidence in the literature suggesting that ESNs are
very good at this. They have performed remarkably well on problems ranging
from seizure detection, to robot control, handwriting recognition, and
financial forecasting, where ESNs have won competitions (LUKOSEVICIUS2009127),
(LukoseviciusMantas2012RCT; RodanA2011MCES; NIPS2010_4056). Impressively, ESNs
outperformed recurrent neural networks (RNNs) and long short term memory
networks (LSTMs) at a chaotic time series prediction task by a factor of over
2400 (Jaeger78). ESNs have also proved themselves competitive in reinforcement
learning (10.1007/11840817_86) and control (Pietz_2021).

To approximate the relationship between the observations and the targets, we
seek a readout vector $W^{*}$ that minimises (over $W$) the regularised least
squares difference between the reservoir states (which represent observations)
and targets, so

$\displaystyle
W^{*}=\operatorname*{arg\,min}_{W\in\mathbb{R}^{N}}\bigg{(}\sum_{k=0}^{\ell-1}\rVert
W^{\top}x_{k}-u_{k}\lVert^{2}+\lambda\lVert W\rVert^{2}\bigg{)}$

where $\lambda>0$ is the Tikhonov regularisation parameter. We can find the
optimal $W^{*}$ easily using linear regression.

###  1.3 Training an ESN on the Lorenz system

####  1.3.1 Creating the reservoir states

To demonstrate everything we have discussed in section 1.2, we will present
two numerical experiments. In section 1.3.2, we use an ESN to learn a mapping
from the discrete time series of values of the $\xi$ component of a trajectory
in the Lorenz attractor (observations) to the $\zeta$ component (targets). In
section 1.3.3, we learn a mapping from the $\xi$ component of the same
discrete trajectory in the Lorenz attractor (observations) to the next value
of the discretely-sampled $\xi$ component (targets). We then use this next
step map to generate a future trajectory for the $\xi$ components. To this
end, let $\phi:\mathbb{R}^{3}\to\mathbb{R}^{3}$ denote a discretisation of the
Lorenz system (1.1) with time step $\tau$ i.e. a discrete-time map of the form

$\displaystyle\phi(\xi,\upsilon,\zeta)=(\xi,\upsilon,\zeta)+\int_{0}^{\tau}(\dot{\xi},\dot{\upsilon},\dot{\zeta})\
dt.$

We set the timestep $\tau=0.01$ and initial condition
$(\xi_{0},\upsilon_{0},\zeta_{0})=(0,1.0,1.05)$. The timestep is small in
comparison to the maximal Lyapunov exponent (see Section 2.2) of the Lorenz
system, $\lambda_{\text{max}}\approx 0.9056$ (Sprott_2003) and the initial
condition is close to the attractor. For these initial conditions and the
parameter values as in (1.1), we computed a trajectory for a 200 time units
(i.e. $\ell=20000$ timesteps), illustrated in Figure 1-1.

Figure 1-1: A typical trajectory of the Lorenz system (1.1) computed for 4000
timesteps, represented by the individual dots at time intervals $\tau=0.01$.
Colour indicates the direction of travel along the trajectory: darkest colours
(blue) at the earliest times and lightest colours (yellow) at the most recent
times.

The observation and target functions are the first and third components of the
Lorenz system $\xi$ and $\zeta$ i.e. we have an observation function
$\omega(\xi,\upsilon,\zeta)=\xi$ so that the observations $z_{k}$ are the
$\xi$ components of the trajectory at the sampled time points $k\tau$, so that

$\displaystyle z_{k}=\omega\circ\phi^{k}(\xi_{0},\upsilon_{0},\zeta_{0}).$

The target function is $u(\xi,\upsilon,\zeta)=\zeta$ so the targets $u_{k}$
are the $\zeta$ components of the trajectory:

$\displaystyle u_{k}=u\circ\phi^{k}(\xi_{0},\upsilon_{0},\zeta_{0}).$

The trajectories of these two components of observations and targets are shown
in Figure 1-2(a) and (b), respectively.

(a) The $\xi$-component of the Lorenz trajectory (vertical axis) plotted
against time (horizontal axis).

(b) The $\zeta$-component of the Lorenz trajectory (vertical axis) plotted
against time (horizontal axis). The black line at the $k^{\mathrm{th}}$
timestep indicates the approximation to this target time series given by
$W^{*\top}x_{k}$.

Figure 1-2: Observations $z_{k}$ and targets $u_{k}$ drawn from the Lorenz
trajectory.

Our goal is to use an ESN to predict the targets based on the observations.
So, we set up an ESN with the following parameters:

  * •

Reservoir size: $N=300$,

  * •

Input matrix $C$ and bias vector $\zeta$: i.i.d uniform random variables $\sim
U[-0.05,0.05]$,

  * •

Reservoir matrix $A$: i.i.d uniform random variables rescaled so that the
matrix 2-norm $\lVert A\rVert_{2}=1$,

  * •

Regularisation parameter $\lambda=10^{-9}$,

which we obtained by hand tuning so that the experiment went smoothly. We have
(so far) little theory for choosing the parameters.

Iterating the ESN with observations $z_{k}$ creates a discrete-time sequence
of reservoir states $x_{k}$, illustrated in Figure 1-3, which shows a
projection of the reservoir states onto their first 3 principal components. To
compute the components, we created a matrix
$X^{\top}\in\mathbb{M}_{N\times\ell}(\mathbb{R})$ with $k$th column $x_{k}$,
and took the singular value decomposition (SVD) $X=U\Sigma V^{\top}$. The
first 3 principal components are the first 3 columns of $V$ which we denote
$V_{1},V_{2},V_{3}$. Then $(V_{1}^{\top}X,V_{2}^{\top}X,V_{3}^{\top}X)$ is the
projection of the reservoir states onto the first 3 principal components.

Figure 1-3: Illustration of the reservoir states of the ESN driven by inputs
$z_{k}$ being the discrete-time samples observed from a trajectory of the
Lorenz system. The figure shows the projection of the reservoir states onto
their first 3 principal components.

The Lorenz system in Figure 1-1 and reservoir dynamics in Figure 1-3 look
rather similar. It appears that the ESN has mapped the Lorenz attractor
$\mathcal{A}$ directly into the reservoir space via a map
$f:(\mathcal{A}\subset\mathbb{R}^{3})\to\mathbb{R}^{N}$. The smoothness of the
Lorenz system, ESN, and reservoir dynamics together suggest that $f$ is at
least continuously differentiable ($C^{1}$). We show in Chapter 3 that $f$ is
a continuously differentiable mapping of the Lorenz system into the reservoir
space. Consequently, the Lorenz dynamics are synchronised to the reservoir
dynamics, via the map $f$, which is called a generalised synchronisation (GS),
in the sense of PhysRevLett.76.1816. The map $f$ is essentially ‘learned’
during training, and this is discussed in detail in Chapter 3.

The image of the GS $f$ appears to smoothly and injectively reproduce the
Lorenz dynamics. This suggests that properties such as Lyapunov spectra,
eigenvalues of fixed points, and homology groups may be preserved under $f$
and replicated in the reservoir dynamics. We explore these claims numerically
using computational methods detailed in Chapter 2. From a purely mathematical
perspective, we know these properties are preserved under diffeomorphism. A
map which is diffeomorphic to its image is called an embedding. Since the GS
$f$ appears to preserve these properties we propose that $f$ is an embedding,
and analyse this claim in Chapter 4.

These observations appear to hold over many different realisations of the ESN
(which has random weights and biases) suggesting the properties hold with
probability close to 1 or perhaps exactly 1. This possibility is studied in
Chapter 4.

####  1.3.2 Learning the targets $\zeta$ from the observations $\xi$

Having made our observations about the reservoir dynamics, we proceed to solve
the least squares problem

$\displaystyle\operatorname*{arg\,min}_{W\in\mathbb{R}^{N}}\sum_{k=0}^{\ell-1}\lVert
W^{\top}x_{k}-u_{k}\rVert^{2}+\lambda\lVert W\rVert^{2}$ (1.3)

to determine the output layer $W^{*}$ using the SVD once again. In particular,
if we denote the $k$th singular value by $\sigma_{k}$, the $k$th column of $U$
by $U_{k}$ and define a vector $Y\in\mathbb{R}^{\ell}$ with $\ell$th component
$u_{\ell}$ then it is a standard calculation (Deblurring_2006) that

$\displaystyle
W^{*}=\sum_{k=1}^{N}\frac{\sigma_{k}U_{k}^{\top}Y}{\sigma_{k}^{2}+\lambda}.$
(1.4)

With this $W^{*}$, we found a good approximation of a mapping from the $\xi$
component to the $\zeta$ component, and the results are shown in Figure 1-2.

The good fit of the nonlinear relationship between the observations and the
targets suggests that the ESN has a sort of universal approximation property,
despite the fitting process being entirely linear. This remarkable feature of
ESNs is analysed in Chapter 5.

####  1.3.3 Predicting the future

Having found a mapping from $\xi$ to $\zeta$, we are now interested in finding
the next step map for time series forecasting. We can feed our predicted next
step back into the ESN to recursively generate a trajectory into the future.
To make this more explicit we define the ESN autonomous phase

$\displaystyle x_{k+1}=\sigma(Ax_{k}+C(W^{*\top}x_{k})+b)$ (1.5)

and recognise that if the ESN is ‘well trained’ (i.e. we have made a good
choice of $W^{*}$) then $W^{*\top}x_{k}\approx u_{k}=z_{k}$ for some
sufficiently long future time, yielding a good forecast for the time series.
So we define the target function

$\displaystyle u(\xi,\upsilon,\zeta):=\omega\circ\phi(\xi,\upsilon,\zeta)$

where the observation function remains $\omega(\xi,\upsilon,\zeta)=\xi$. We
solve the least squares problem (1.3) using the SVD (1.4) once again, then
generated a 4000 timestep (40 time units) future trajectory, starting at
$t=40$, using the ESN autonomous phase (1.5). The future $\xi$ components and
reservoir states generated by the autonomous ESN are shown in Figures 1-4(b)
and 1-4(a) respectively.

(a) $\xi$ component of the Lorenz system (vertical axis) evolves in time
(horizontal axis) up until time $t=40$. After this point, the autonomous ESN
predicts a future trajectory (multicoloured), which is compared to the true
future trajectory (black).

(b) Illustration of the reservoir dynamics of the autonomous ESN, projected
onto the first 3 principal components of the (non-autonomous) reservoir
dynamics.

Figure 1-4: ESN autonomous phase.

The ESN autonomous dynamics appear to have replicated a trajectory that is
diffeomorphic to the Lorenz dynamics, even as the trajectories diverge in
absolute value.

This diffeomorphism is not guaranteed by a universal approximation result,
which does not rule out the possibility of small errors accumulating until the
autonomous dynamics no longer resemble the Lorenz dynamics at all. We explore
in Chapter 5 how the _structural stability_ of the Lorenz system ensures that
qualitative features of the dynamics are preserved even as the trajectories
diverge.

###  1.4 The remainder of the thesis

Predicting the future trajectory of the Lorenz system is a supervised learning
problem. In almost all real applications of supervised learning, the data is
noisy and there is uncertainty. This is in contrast to the deterministic
evolution of the Lorenz system we considered in the previous sections. Thus,
in Chapter 6 we develop stochastic analogues for many of the deterministic
results that appear in Chapters 3-5 and that were touched on here.

In Chapter 6 we extend the stochastic set-up to encompass reinforcement
learning problems. In the reinforcement learning paradigm an agent explores
its environment by executing actions and collects rewards influenced by the
environment and its actions. The goal of the agent is to learn the impact of
its actions on the environment, and develop a policy that maximises its
expected sum of future rewards. Learning about the environment and refining
the policy often happen at the same time. We prove that ESNs are especially
well suited to reinforcement learning problems where the environment evolves
randomly, in a generally non-Markovian manner, and the ESN’s memory of the
past is crucial for making good decisions.

In Chapter 7, we use reservoir computers to solve linear PDEs. The assumptions
made in this Chapter are rather different to those made in the preceding
Chapters, where we imagine that our agent does not know the equations that
govern the environment. In previous chapters, the agent must learn a model
from scratch using a reservoir computer trained on data. However in Chapter 7
we allow the agent access to the PDE explicitly, and use a reservoir computer
to approximate the solution.

##  Chapter 2 Computational Methods

The numerical results in Chapter 1 suggests that under appropriate conditions
the ESN autonomous dynamics are diffeomorphic to the dynamics of the drive
system. In this case, the autonomous dynamics should inherit qualitative
properties of the drive system, for example the existence of fixed points,
eigenvalues of linearisations around them, and for chaotic attractors, their
Lyapunov spectra and homology groups. Properties like these can be constructed
from a time series using an assortment of techniques. In particular the
computation of homology groups is a central idea in the field known loosely as
_computational topology_ also known as _algorithmic topology_ or _applied
topology_. We will demonstrate some of these methods by recovering the
topological properties of the Lorenz system using an ESN. The presentation in
this chapter closely follows the numerical experiments in
embedding_and_approximation_theorems. To demonstrate these computational
methods, we set up an ESN with the following parameters:

  * •

Reservoir size: $N=300$,

  * •

Elements of the input matrix $C$ and bias vector $\zeta$: i.i.d uniform
Gaussian variables $\sim\mathcal{N}(0,0.1^{2})$,

  * •

Reservoir matrix $A$: Erdős-Rényi matrix with mean $6$ and connection weights
(where they are non-zero) i.i.d Gaussian, re-scaled such that the spectral
radius $\rho=1$.

  * •

Regularisation parameter $\lambda=10^{-6}$.

###  2.1 Fixed Points

If the Generalised Synchronisation (GS) $f$ is an embedding, then $f$ will
embed the fixed points of the Lorenz system into the reservoir space. Moreover
if the autonomous ESN approximates the embedded Lorenz dynamics sufficiently
well, the autonomous dynamics will have fixed points with the same
linearisation as the Lorenz system. To verify this, we searched for the
autonomous ESN’s fixed points using Newton’s method.

In particular, we found the fixed point of

$\displaystyle\psi(x):=\sigma(Ax+C(W^{\top}x)+b)$ (2.1)

using Newton iterations

$\displaystyle x_{k+1}=x_{k}-(J(x_{k})-I)^{-1}(\psi(x_{k})-x_{k})$

where $J(x_{k})$ is the Jacobian of $\psi$ evaluated at $x_{k}$. To reduce
numerical instability and computational time we do not find
$(J(x_{k})-I)^{-1}$ explicitly but instead solve

$\displaystyle(J(x_{k})-I)(x_{k+1}-x_{k})=-\psi(x_{k})+x_{k}$

for $(x_{k+1}-x_{k})$. The iterates are shown in Figure 2-1, and it appears by
eye the algorithm was successful and a fixed point was found.

Figure 2-1: The driven reservoir dynamics are plotted in blue and autonomous
dynamics are plotted in red. Both were projected onto the first three
principal components of the driven dynamics, then the axes are rotated such
that the projection appears on the first 2 components. The black line
indicates the iterates of Newton’s method, used to locate a fixed point - the
method eventually converges to a fixed point in the middle of the right wing
of the figure. We can see by eye that the reservoir dynamics appear by eye to
be diffeomorphic to the Lorenz system.

Further, if the GS $f$ is a $C^{1}$ embedding of the original dynamics, we
expect $f$ to preserve the stability of fixed points, i.e. we expect the
eigenvalues of the linearisation of the autonomous phase to be preserved at
every fixed point. Now, comparing the eigenvalues of the linearisation of the
Lorenz system and autonomous phase at the respective fixed points requires
some subtlety, because the Lorenz system is a continuous time flow, while the
autonomous phase is a discrete time map. So, we began by considering one of
the known fixed points found in the Lorenz attractor’s wings

$\displaystyle m^{*}=(6\sqrt{2},6\sqrt{2},27)$

and noted the Jacobian $J$ of the continuous time Lorenz system evaluated at
the fixed point $m^{*}$ is therefore

$\displaystyle J\Big{|}_{m^{*}}=\begin{bmatrix}-10&10&0\\\ 1&-1&-6\sqrt{2}\\\
6\sqrt{2}&6\sqrt{2}&-8/3\end{bmatrix}.$

So the discretised Lorenz system

$\displaystyle\phi(\xi,\upsilon,\zeta)=(\xi,\upsilon,\zeta)+\int_{0}^{\tau}(\dot{\xi},\dot{\upsilon},\dot{\zeta})\
dt.$

admits a linearisation $\phi^{*}$ about the fixed point $m^{*}$

$\displaystyle\phi^{*}(\xi,\upsilon,\zeta)=\exp\bigg{(}J\Big{|}_{m^{*}}\tau\bigg{)}\begin{bmatrix}\xi\\\
\upsilon\\\ \zeta\end{bmatrix},$

where $\exp:\mathbb{M}_{3\times 3}(\mathbb{R})\to\mathbb{M}_{3\times
3}(\mathbb{R})$ denotes the matrix exponential defined

$\displaystyle\exp(Q)=\sum_{n=0}^{\infty}\frac{Q^{n}}{n!}.$

The matrix $\exp(J\rvert_{m^{*}}\tau)$ has 3 eigenvalues, which we have
compared with the ESN autonomous eigenvalues in Figure 2-2.

Figure 2-2: Comparison of the eigenvalues of the discrete time evolution
operator for the Lorenz system and the autonomous ESN (2.1). The three
eigenvalues of the linearisation of the Lorenz system on the fixed point
inside one of the Lorenz attractor’s wings are represented by blue crosses.
The 300 eigenvalues of the linearisation of the ESN autonomous system at the
fixed point found with Newton’s method are represented by red dots. The
concentric grey circles have radius 1/2 and 1 respectively.

If the GS $f$ is indeed a $C^{1}$ embedding, the dynamics of the autonomous
phase are diffeomorphic to the discrete time Lorenz system on some three
dimensional submanifold. This manifold (linearised at the fixed point $m^{*}$)
is spanned by three eigenvectors, each with an associated eigenvalue, which
will coincide with the eigenvalues of the linearisation of the Lorenz system
at the fixed point. Figure 2-2 appears to show three overlapping eigenvalues,
suggesting that the autonomous phase is diffeomorphic to the Lorenz system (at
least in a neighbourhood of $m^{*}$). This is particularly remarkable because
$m^{*}$ is distant from the training data. It appears that the ESN has
inferred the existence, position and eigenvalues of a fixed point from
training data which contains no fixed points. In machine learning parlance, we
might say that the ESN has generalised patterns in the training data to an
unseen region of the phase space. The remaining 297 eigenvalues of the ESN
autonomous system have absolute value less than 1, suggesting that the fixed
point is attractive with respect to them.

###  2.2 Lyapunov Spectra

Another invariant of the Lorenz system preserved under diffeomorphism is the
Lyapunov spectrum, which captures how quickly very close trajectories diverge
from eachother, and is used as a measure of chaos. To define the spectrum, let
$J$ be the Jacobian of the evolution operator of a continuous time dynamical
system. Let $Y$ be the solution of the ODE $\dot{Y}=JY$ with initial condition
$Y(0)=m_{0}$. Then the Lyapunov spectrum of the invariant set containing
$m_{0}$ is the spectrum of the matrix $\Lambda$ defined

$\displaystyle\Lambda=\lim_{t\to\infty}\frac{1}{2t}YY^{\top}.$

Each eigenvalue in the spectrum is called a _Lyapunov exponent_ to signify
that two initially close trajectories diverge or converge exponentially fast
with exponentiation constant in the direction of each eigenvector of $J$ given
by a Lyapunov exponent. Further details and discussion is given in
DARBYSHIRE1996287. The Lyapunov spectrum for the Lorenz system was estimated
by Sprott_2003 as $\\{0.9056,0,-14.5723\\}$. In order to compare the Lorenz
spectrum to the spectrum of the autonomous ESN, we computed the autonomous
system’s spectrum using the discrete time $QR$ method discussed in
DARBYSHIRE1996287 and plotted each Lyapunov exponent against the known
exponents of the Lorenz system in Figure 2-3. The method is described by the
following algorithm

Algorithm 1 Compute the Lyapunov spectrum

1:Let $Q_{0}$ be an arbitrary $N\times N$ orthogonal matrix

2:Let $J\lvert_{m_{0}}$ be the Jacobi matrix of the autonomous ESN evaluated
at $m_{0}$

3:for each $j$ from $0$ to $n$

4: Take the $QR$ decomposition $Q_{j+1}R_{j+1}=Q_{j}J\lvert_{\phi^{j}(m_{0})}$

5: Compute $\Lambda_{j+1}=\log(\text{diag}(R_{j+1}))$ where $\text{diag}(A)$
creates a vector from the diagonal entries of $A$)

6:Compute $\Lambda=\frac{1}{n}\sum_{j=1}^{n}\Lambda_{j}$. The exact spectrum
is obtained as $n\to\infty$

We found that the largest 2 exponents of the Lorenz system and the ESN map are
in good agreement with each other while there was a significant difference
between the next-largest in each case. This problem was also noted and
encountered by Pathak2017, and we do not have a satisfactory explanation for
this.

Figure 2-3: The Lyapunov spectrum of the autonomous phase as the iterates
increases is shown. The true Lyapunov exponents of the autonomous phase is
given by the limit of these exponents as the iterations tend to infinity.
These autonomous exponents are compared to the black dotted lines representing
the 3 exponents of the Lorenz system.

In the next section we will consider a somewhat different invariant that is
preserved under the embedding $f$: the homology groups of the attractor (for a
finite length of time intergration). We note that the true attractor has a
fractal structure, and that the closure of the true attractor contains the
fixed points at the origin and those in the middle of each wing. We therefore
proceed carefully in the upcoming sections. We will compute the these groups
using _persistent homology_.

###  2.3 Persistent homology

####  2.3.1 Introduction

The objective of an unsupervised learning algorithm is to reveal the hidden
relations and structures within a dataset, without the use of a training set.
A classic unsupervised learning algorithm is principal component analysis
(PCA), which reveals linear relations in high dimensional data. Another is
$k$-means clustering which partitions a data set into $k$ clusters. In this
section, we will explore a completely different kind of structure that a data
set might have: homology.

Roughly speaking, homology is the study of objects under continuous
deformation, which preserves their connected components, holes and voids. The
sphere is made up of 1 connected component, and encloses a 3 dimensional void.
The torus is also made up of 1 connected component, and also encloses a 3
dimensional void, but unlike the sphere contains a 2 dimensional hole, or
tunnel. It is this distinction between the sphere and the torus that
distinguish their homology. Given a constellation of data points then, what
does it mean to say it has a hole or encloses a void? Astronomers have for
thousands of years looked into the night sky, and in their minds’ eye,
connected the stars to create shapes with holes and voids. The emerging field
of persistent homology offers a mathematical formalism for obtaining similar
results. This subsection will outline the main results of persistent homology
and explain how to practically compute the homology of a real data set. We
draw heavily on the work of Ghrist2007 and Ghrist2014.

####  2.3.2 Theory

##### Simplicial complexes

A data set can be transformed into a geometric object called a simplicial
complex with computable homology. Every singleton set containing a single
point is a point, also known as a 0-simplex. Any pair of points in the data
set can be connected together to form a line, or 1-simplex. Any triple of
points can be collected to form a triangular face, also known as a 2-simplex.
Any quadruple of points can together form a tetrahedron, (a 3-simplex) and so
on. A set of simplices form a simplicial complex. To make this idea rigorous
we will introduce a few definitions.

######  Definition 2.3.1.

(Closed under restriction) A collection of sets $S$ is closed under
restriction if any subset of a set in $S$ is also in $S$.

######  Definition 2.3.2.

(Simplicial Complex) Let $X$ be a finite set. A simplicial complex $S$ is a
set of subsets of $X$ that is closed under restriction.

######  Definition 2.3.3.

(Simplex) Let $S$ be a simplicial complex. Then a non-empty $\sigma\in S$
containing $k$ elements is called a $k-1$ simplex.

Given a data set, how do we decide whether a pair of points ought to be
connected to form a 1-simplex? Or indeed how do we decide whether some set of
$k$ points ought to be brought together to form a $k-1$ simplex? One method
developed by Leopold Vietoris1927 and Eliyahu Rips supposes every data point
is the centre of a ball of radius $\epsilon$. If a set of $k$ points are each
contained in the ball of every other, then they form a $k-1$ simplex,
otherwise they do not.

######  Definition 2.3.4.

(Data set) A data set $X$ is a finite list of elements in $\mathbb{R}^{n}$.

######  Definition 2.3.5.

(Vietoris-Rips Complex) The Vietoris-Rips (Rips for short) complex
$\mathcal{R}_{\epsilon}$ of a data set $X$ is the set of subsets of $X$
defined as those subsets $\sigma\subset X$ with members whose pairwise
distance (in the Euclidean metric) is less than $\epsilon$.

Another method of connecting the dots is called the C̆ech Complex, and
supposes every data point is the centre of a ball of radius $\epsilon/2$. The
division by $2$ is a useful convention for reasons that will soon be made
clear. This time, we say a set of $k$ data points form a $k-1$ simplex if and
only if their balls have a common intersection.

######  Definition 2.3.6.

(C̆ech Complex) The C̆ech complex $\mathcal{C}_{\epsilon}$ of a data set $X$
is a set of subsets of $X$ defined by the statement that $\sigma\subset X$ is
in $\mathcal{C}_{\epsilon}$ if and only if

$\displaystyle\bigcap\limits_{x\in\sigma}B_{\epsilon/2}(x)\neq\emptyset.$

######  Remark 2.3.7.

The Rips complex and C̆ech Complex are clearly closed under restriction and
are therefore simplicial complexes.

The C̆ech Complex has the elegant feature of being homotopy equivalent to the
union of $\epsilon/2$ balls centred at each data point.

######  Definition 2.3.8.

(Homotopy) Let $X,Y$ be topological spaces and $f,g:X\to Y$ continuous maps. A
homotopy between $f,g$ is a continuous function $H:X\times[0,1]\to Y$ such
that $H(x,0)=f(x)$ and $H(x,1)=g(x)$.

######  Definition 2.3.9.

(Homotopy Equivalent) Two topological spaces $X$, $Y$ are homotopy equivalent
if there is a pair of continuous maps $f:X\to Y$ and $g:Y\to X$ such that
$f\circ g$ is homotopic to the identity $\text{id}_{Y}$ and $g\circ f$ is
homotopic to the identity $\text{id}_{X}$.

One could say the C̆ech Complex is the nervous system found inside a muscular
subspace of $\mathbb{R}^{n}$ composed of $\epsilon$-balls. To put this more
formally we introduce homotopy and the Nerve Lemma. An example is shown in
Figure 2-4.

######  Lemma 2.3.10.

(Nerve Lemma) The C̆ech Complex $\mathcal{C}_{\epsilon}$ of a data set $X$ is
homotopic to

$\displaystyle\bigcup\limits_{x\in X}B_{\epsilon/2}(x)\subset\mathbb{R}^{n}.$

###### Proof.

Alexandroff1928. ∎

Doctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of Philosophy Figure 2-4: The C̆ech
complex is homotopic to the union of $\epsilon$-balls, which in this figure,
is a genus 2 surface.

##### Persistence complexes

Given a data set and a particular $\epsilon$, we can create a Rips or C̆ech
complex with holes and voids, and therefore a homology we can scrutinise. But
the question remains as to how to chose $\epsilon$. A sufficiently small
$\epsilon$ will not connect any points and result in a boring list of $0$
dimensional simplices, while a sufficiently large $\epsilon$ will connect
every point to every other forming a single gigantic simplex. In persistent
homology we consider all values of $\epsilon$, and observe which homological
features of the data persist over long intervals of $\epsilon$ and which
flutter rapidly in and out of existence. Those features which persist are
considered real, while their more transient companions are dismissed as noise.

To make this idea concrete, we consider a data set $X$ and create a sequence
of complexes called a _persistence complex_ as follows. The first complex
$\mathcal{R}_{\epsilon_{0}}$ is a set of $0$-simplices; no data point is
connected to any other. Starting from $\epsilon=0$ we let $\epsilon$ grow
until some pair of points become connected to form a $1$-simplex, giving rise
to the next complex in the sequence $\mathcal{R}_{\epsilon_{1}}$. We let
$\epsilon$ grow some until a new pair of points become connected, giving rise
to the third complex in the sequence $\mathcal{R}_{\epsilon_{2}}$. We continue
to add terms to the sequence until $\epsilon$ is large enough that every point
in the data set belongs to the same simplex. Whether some set of $k$ points
are connected and form a $k-1$ simplex depends of course on whether we are
considering a sequence of Rips complexes or C̆ech complexes. The following
definitions formalise these ideas

######  Definition 2.3.11.

(Filtration) A filtration is a collection of sets $X_{i}$ with the property
that if $i<j$, then $X_{i}\subset X_{j}$.

######  Definition 2.3.12.

(Rips Persistence complex) Let $X$ be a data set and notice that the subsets
of $\\{\mathcal{R}_{\epsilon}|\epsilon>0\\}$ form a filtration
$\mathcal{R}_{\epsilon_{0}}\subset\mathcal{R}_{\epsilon_{1}}\subset\mathcal{R}_{\epsilon_{2}}\subset\ldots\subset\mathcal{R}_{\epsilon_{m}}$
for $\epsilon_{0}<\epsilon_{1}<\ldots<\epsilon_{m}$. The sequence of inclusion
maps

$\displaystyle\mathcal{R}_{\epsilon_{0}}\xhookrightarrow{\iota_{1}}\mathcal{R}_{\epsilon_{1}}\xhookrightarrow{\iota_{2}}\mathcal{R}_{\epsilon_{2}}\xhookrightarrow{\iota_{3}}\ldots\xhookrightarrow{\iota_{m}}\mathcal{R}_{\epsilon_{m}}$

is called the Rips persistence complex.

######  Definition 2.3.13.

(C̆ech Persistence complex) Let $X$ be a data set and consider that the
subsets of $\\{\mathcal{C}_{\epsilon}|\epsilon>0\\}$ form a filtration
$\mathcal{C}_{\epsilon_{0}}\subset\mathcal{C}_{\epsilon_{1}}\subset\mathcal{C}_{\epsilon_{2}}\subset\ldots\subset\mathcal{C}_{\epsilon_{m}}$.
The sequence of inclusion maps

$\displaystyle\mathcal{C}_{\epsilon_{0}}\xhookrightarrow{\iota_{1}}\mathcal{C}_{\epsilon_{1}}\xhookrightarrow{\iota_{2}}\mathcal{C}_{\epsilon_{2}}\xhookrightarrow{\iota_{3}}\ldots\xhookrightarrow{\iota_{m}}\mathcal{C}_{\epsilon_{m}}$

is called the C̆ech persistence complex.

For a given data set $X$ the Rips and C̆ech chain complexes are not the same
in general, but we can always squeeze a C̆ech complex between two Rips
complexes as follows.

######  Lemma 2.3.14.

(Squeezing Lemma) For any $\epsilon>0$ there is a chain of inclusion maps

$\displaystyle\mathcal{R}_{\epsilon}\xhookrightarrow{}\mathcal{C}_{\epsilon\sqrt{2}}\xhookrightarrow{}\mathcal{R}_{\epsilon\sqrt{2}}.$

###### Proof.

diSilva ∎

It follows from this Lemma that any homological feature of the Rips complex
that persists from $\epsilon$ to $\epsilon\sqrt{2}$ is also a feature of the
C̆ech complex $\mathcal{C}_{\epsilon\sqrt{2}}$. An illustration of this
squeezing phenomenon is shown in Figure 2-5.

Doctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of Philosophy Figure 2-5: A filtration of
C̆ech complexes for increasing $\epsilon$. We can see that for $\epsilon<1$
the C̆ech complex is a collection of $3$ points, for $1<\epsilon<\sqrt{2}$ the
C̆ech complex forms a 1 dimensional loop, and for $\epsilon>1$ the complex is
a triangle. The Rips complex for the same 3 points is $3$ isolated points when
$\epsilon<1$ and a triangle otherwise. The one dimensional loop in the
sequence of C̆ech complexes is squeezed between 2 Rips complexes; 3 isolated
points and a triangle.

##### Simplicial homology

We first seek the homology of a simplicial complex defined at a particular
value of $\epsilon$. To this end, we note that a simplicial complex is
composed of some $m$-simplices, $(m-1)$-simplices, $(m-2)$-simplices and so
on, until a final set of $0$-simplices. Each $k$ simplex has a boundary, which
is composed of $(k-1)$-simplices. The boundary of a 1-simplex (line) is
composed of the two $0$-simplices (points) with one at each end. The boundary
of a $2$-simplex (triangle) is composed of the three $1$-simplices (lines)
that enclose it. The boundary of a $3$-simplex (a tetrahedron) is composed of
the four $2$-simplices (triangles) that enclose it and so on. The boundary of
a $0$-simplex is the empty set.

Doctor of PhilosophyDoctor of PhilosophyDoctor of Philosophy1Doctor of
Philosophy2Doctor of Philosophy3Doctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of Philosophy4Doctor of Philosophy5Doctor of Philosophy6
Figure 2-6: The boundary of this 2-complex is the set of 1-simplices
$\\{(12),(23),(34),(45),(56),(61)\\}$

More generally, a set of $k$-simplices has a boundary composed of
$(k-1)$-simplices. For example, a pair of triangles welded together along one
side has a boundary composed of four $1$-simplices. A pair of squares welded
together along one edge is shown in Figure 2-6. A pair of tetrahedra welded
together along one face has a boundary comprised of $6$ triangular faces. This
operation of welding simplices together will be denoted with a $+$. Under the
operation $+$ the $k$ simplices form a group.

######  Definition 2.3.15.

(Simplex Group) Let $s_{k}$ be the set of all $k$-simplices in a simplicial
complex. Let $S_{k}=P(s_{k})$: the powerset of $s_{k}$. Let $+$ be a binary
operation on the elements of $S_{k}$ defined

$\displaystyle\sigma+\tau=(\sigma\cup\tau)/(\sigma\cap\tau)$

for $H,G\in S_{k}$. Under the operation $+$ the set $S_{k}$ forms an Abelian
group with identity $\emptyset$ called the simplex group.

It therefore makes sense for each $k$ to define a map $\partial_{k}$ which
sends a set of $k$ simplices to the set of $(k-1)$ simplices which comprise
its boundary. It turns out this map is a group homomorphism.

######  Lemma 2.3.16.

Let $s_{k}$ be the set of all $k$-simplices in a simplicial complex. Let
$S_{k}=P(s_{k})$. Let $\partial_{k}:S_{k}\to S_{k-1}$ be defined

$\displaystyle\partial_{k}(\sigma)=s_{k-1}\cap\sum_{\varsigma\in\sigma}\varsigma$

then $\partial_{k}$ is a homomorphism.

###### Proof.

The proof proceeds directly

$\displaystyle\partial_{k}(\sigma)+\partial_{k}(\tau)=\big{(}\partial_{k}(\sigma)\cup\partial_{k}(\tau)\big{)}/\big{(}\partial_{k}(\sigma)\cap\partial_{k}(\tau)\big{)}$
$\displaystyle=\Bigg{(}s_{k-1}\cap\sum_{\varsigma\in(\sigma\cup\tau)}\varsigma\Bigg{)}\Bigg{/}\Bigg{(}s_{k-1}\cap\sum_{\varsigma\in(\sigma\cap\tau)}\varsigma\Bigg{)}$
$\displaystyle=s_{k-1}\cap\sum_{\varsigma\in(\sigma\cup\tau)/(\sigma\cap\tau)}\varsigma$
$\displaystyle=s_{k-1}\cap\sum_{\varsigma\in(\sigma+\tau)}\varsigma$
$\displaystyle=\partial_{k}(\sigma+\tau).$

∎

A set of $k$ simplices that have an empty boundary is called a cycle. For
example a set of $1$-simplices arranged in a pentagon form a loop, which has
no $0$ dimensional boundary. Similarly a set of $2$-simplices welded together
to form a cube form a closed surface, which has no $1D$ boundary. Formally, a
set of $k$ simplices is a cycle if the boundary map $\partial_{k}$ sends the
set to the empty set $\emptyset$. This motivates the following definition.

######  Definition 2.3.17.

(Cycle Group) Let $S$ be a simplicial complex. The group of $k$-cycles of $S$
is

$\displaystyle Z_{k}=\text{ker}(\partial_{k}).$

Since the boundary map returns boundaries, the image of the $(k+1)$ boundary
map is the group of boundaries among the group of $k$-simplices. This is
formalised in the following definition.

######  Definition 2.3.18.

(Boundary Group) Let $S$ be a simplicial complex. The group of $k$-boundaries
of $S$ is

$\displaystyle B_{k}=\text{Im}(\partial_{k+1})$

We say that 2 cycles are equivalent if we can weld them together to form a
boundary. Under this equivalence relation, the equivalence class of $k$-cycles
can be interpreted as a $(k+1)$-dimensional hole. For example an equivalence
class of $1$-cycles represents a $2$ dimensional hole - exactly the type of
hole that distinguishes the torus from the sphere. An equivalence class of
$2$-cycles can be interpreted as a $3$ dimensional void, like the volume
enclosed by the torus, or the volume enclosed by the sphere. The equivalence
classes of cycles form a group as defined below.

######  Definition 2.3.19.

(Homology Group) The $k$th homology group is the quotient of the cycles group
by the boundary group:

$\displaystyle
H_{k}=\text{ker}(\partial_{k})/\text{Im}(\partial_{k+1})=Z_{k}/B_{k}.$

######  Definition 2.3.20.

(Betti Numbers) The $k$th Betti number is

$\displaystyle\beta_{k}=\text{rank}(H_{k}).$

The Betti numbers $\beta_{0},\beta_{1},\beta_{2},...$ of a simplicial complex
characterise its homology. $\beta_{0}$ is the number of connected components,
$\beta_{1}$ is the number of $2$ dimensional holes, $\beta_{2}$ is the number
of $3$ dimensional voids, and the remaining $\beta_{k}$ represent higher
dimensional voids. The Betti numbers of a simplicial complex are the
coefficients of a polynomial called the Poincaré Polynomial.

######  Definition 2.3.21.

(Poincaré polynomial) The Poincaré polynomial of a simplicial complex $S$ is
defined

$\displaystyle p[X]=\sum_{i}\beta_{i}X^{i}.$

##### Persistent homology

In the previous section we introduced the homology group of a single
simplicial complex - but we can do better! We will now present the persistent
homology group which characterises the homology of a persistence complex that
persists over a range of $\epsilon$.

######  Definition 2.3.22.

(Induced homomorphism) Let $i<j$ and $\mathcal{R}_{\epsilon_{i}}$,
$\mathcal{R}_{\epsilon_{j}}$ be a pair of Rips complexes from the same Rips
persistence complex. Let $Z^{i}_{k}$ and $Z^{j}_{k}$ be the $k$-cycle groups
of the complex $\mathcal{R}_{\epsilon_{i}}$ and $\mathcal{R}_{\epsilon_{j}}$
respectively. Let $\iota:Z^{i}_{k}\to Z^{j}_{k}$ be the inclusion map. Then
$\iota_{*}:H^{i}_{k}\to H^{j}_{k}$ is called the induced homomorphism of
$\iota$.

######  Lemma 2.3.23.

Let $x$ be a $k$-cycle in $\mathcal{R}_{\epsilon_{i}}$ so $x\in Z^{i}_{k}$.
Let $[x]_{i}$ denote the equivalence class of $x$ under the quotient relation

$\displaystyle H^{i}_{k}=Z^{i}_{k}/B^{i}_{k}.$

Then the induced homomorphism $\iota_{*}$ satisfies

$\displaystyle[\iota(x)]_{j}=\iota_{*}([x]_{i}).$

###### Proof.

Ghrist2014. ∎

######  Definition 2.3.24.

(Persistent Homology Group) Let $\mathcal{R}$ be a Rips persistence complex.
The $(i,j)$ persistent $k$ homology group of $\mathcal{R}$ is

$\displaystyle H^{i\to j}_{k}=\text{Im}(\iota_{*}).$

We should take a moment to celebrate. The $(i,j)$ persistent $k$ homology
group $H^{i\to j}_{k}$ represents exactly those $(k+1)$ dimensional holes that
persist over the interval $[\epsilon_{i},\epsilon_{j}]$. For example $H^{1\to
2}_{2}$ contains all $3$-dimensional voids that persist from
$\epsilon=\epsilon_{1}$ to $\epsilon=\epsilon_{2}$. This group is exactly what
we want.

##### Computing the persistent homology

We have presented an algebraic derivation of the persistent homology group.
This section will provide some details about how to determine the group in
practice. The key observation is that we can convert the group theory problem
into a linear algebra problem. Let $S$ be an $n$-simplex (and therefore a
simplical complex). Let the $k$th simplex in $S$ be the $k$th canonical unit
vector of the $n$ dimensional vector space over the field of $2$ elements
$\mathbb{F}^{n}_{2}$. i.e. the first simplex is $(1,0,0,...,0)$, the second is
$(0,1,0,...,0)$ and so on. It turns out the $k$th simplex group $S_{k}$ is
isomorphic to a subspace of $\mathbb{F}^{n}_{2}$. The boundary map
$\partial_{k}:S_{k}\to S_{k-1}$ turns out to be a linear map between
subspaces. The cycle group $Z_{k}$ is isomorphic to a boundary space
determined by computing $\text{ker}(\partial_{k})$ using the standard
techniques of linear algebra. The boundary group $B_{k}$ is just the column
span of the matrix representation of $\partial_{k}$. The quotient

$\displaystyle H_{k}=Z_{k}/B_{k}$

is isomorphic to the column span of the vectors that are not in the span of
the $B_{k}$. Using this vector space representation, we are ready to outline
an algorithm for computing the persistent homology group. Suppose we have an
$n$ dimensional data set $X$ comprising $m$ points. Then Algorithm 2 is
pseudocode for computing the Rips persistent homology groups of the data set
$X$.

Algorithm 2 Compute $H^{i\to j}_{k}$

1:assign to every simplex in $P(X)$ a unique natural number

2:compute the distance between every pair of points in $X$

3:list the distances in ascending order

4:create a filtration of Rips complexes

5:for each $j$ indexing the $j$th Rips complex in the filtration

6: partition the simplices of the Rips complex into simplices of equal
dimension $k$

7: for all dimensions $k>0$

8: for each simplex in the set

9: retrieve the unique natural number $\ell$ of the simplex

10: set the $\ell$th column of matrix representing $\partial_{k}$ as the
vector representing the boundary of the $\ell$th simplex

11: find a basis for the vector space isomorphic to
$Z_{k}=\text{ker}(\partial_{k})$

12: find a basis for the vector space isomorphic to
$B_{k-1}=\text{Im}(\partial_{k})$

13: if $k>1$

14: find a basis for the vector space isomorphic to
$H^{j}_{k-1}=Z_{k-1}/B_{k-1}$

15: for all $i<j$

16: for each vector $v$ in the basis of $H^{i}_{k}$

17: if $v$ is in the span of $B^{j}_{k}$

18: $v$ is not in the $(i,j)$ persistent $k$ homology group

19: else

20: $v$ is in the $(i,j)$ persistent $k$ homology group

This algorithm is naïve, and if properly implemented would surely run slower
than more efficient algorithms like Ripser developed by ctralie2018ripser.

Computing the persistent homology group in practice is a confusing ordeal, so
we will offer a hopefully illustrative example. Consider $6$ nodes arranged in
a regular hexagon, and contemplate the resultant Rips filtration. This is
illustrated in Figure 2-7. Tables 2.1 and 2.2 show the matrix representations
$D_{2}$ and $D_{1}$ of the boundary maps $\partial_{2}$ and $\partial_{1}$. We
can find a basis for the column span of $D_{2}$ and null space of $D_{1}$ and
determine that the homology group

$\displaystyle H^{2}_{1}=\text{ker}(\partial_{1})/\text{Im}(\partial_{2})$

is isomorphic to the space spanned by those vectors in the null space of
$D_{1}$ not in the column span of $D_{2}$. We shall call this list of vectors
$\alpha$. We can similarly compute a basis for a vector space isomorphic to
$H^{3}_{1}$, which we will call $\beta$. The persistent homology group
$H^{2\to 3}_{1}$ is isomorphic to the span of those vectors $\alpha$ not in
the span of the vectors in $\beta$.

Doctor of PhilosophyDoctor of Philosophy1Doctor of Philosophy6Doctor of
Philosophy5Doctor of Philosophy4Doctor of Philosophy3Doctor of
Philosophy2Doctor of Philosophy$\epsilon=0$

Doctor of PhilosophyDoctor of Philosophy1Doctor of Philosophy6Doctor of
Philosophy5Doctor of Philosophy4Doctor of Philosophy3Doctor of
Philosophy2Doctor of Philosophy$\epsilon=1$

Doctor of PhilosophyDoctor of Philosophy1Doctor of Philosophy6Doctor of
Philosophy5Doctor of Philosophy4Doctor of Philosophy3Doctor of
Philosophy2Doctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of Philosophy$\epsilon=\sqrt{3}$

Doctor of PhilosophyDoctor of Philosophy1Doctor of Philosophy6Doctor of
Philosophy5Doctor of Philosophy4Doctor of Philosophy3Doctor of
Philosophy2Doctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyDoctor of Philosophy$\epsilon=2$

Figure 2-7: A filtration of Rips complexes $\mathcal{R}_{0}\xhookrightarrow{}\mathcal{R}_{1}\xhookrightarrow{}\mathcal{R}_{\sqrt{3}}\xhookrightarrow{}\mathcal{R}_{2}.$ | (123) | (234) | (345) | (456) | (561) | (612)
---|---|---|---|---|---|---
(12) | 1 | 0 | 0 | 0 | 0 | 1
(23) | 1 | 1 | 0 | 0 | 0 | 0
(34) | 0 | 1 | 1 | 0 | 0 | 0
(45) | 0 | 0 | 1 | 1 | 0 | 0
(56) | 0 | 0 | 0 | 1 | 1 | 0
(61) | 0 | 0 | 0 | 0 | 1 | 1
(13) | 1 | 0 | 0 | 0 | 0 | 0
(35) | 0 | 0 | 1 | 0 | 0 | 0
(51) | 0 | 0 | 0 | 0 | 1 | 0
(62) | 0 | 0 | 0 | 0 | 0 | 1
(24) | 0 | 1 | 0 | 0 | 0 | 0
(46) | 0 | 0 | 0 | 1 | 0 | 0
Table 2.1: For $\sqrt{3}<\epsilon<2$, the Rips complex in Figure 2-7 has $6$ traingular faces (2-simplices) each denoted by listing the 3 nodes they contain, e.g. (123). The boundary map $\partial_{2}$ sends a face to a set of 3 edges ($1$-simplices). For example $\partial_{2}$ sends the face (123) to the set of edges { (12),(23),(13) }. We can therefore represent $\partial_{2}$ as a matrix with $ij$th entry the $ij$th entry of this table. | (12) | (23) | (34) | (45) | (56) | (61) | (13) | (35) | (51) | (62) | (24) | (46)
---|---|---|---|---|---|---|---|---|---|---|---|---
(1) | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0
(2) | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0
(3) | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0
(4) | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1
(5) | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0
(6) | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1
Table 2.2: For $\sqrt{3}<\epsilon<2$, the Rips complex in Figure 2-7 has $12$
edges (1-simplices) each denoted by listing the 2 nodes they contain, e.g.
(12). The boundary map $\partial_{1}$ sends an edge to a set of 2 nodes
($0$-simplices). For example $\partial_{1}$ sends the face (12) to the set of
nodes { (1),(2) }. We can therefore represent $\partial_{1}$ as a matrix with
$ij$th entry the $ij$th entry of this table.

At any Rips complex in the filtration, a basis element that has not featured
in any previous complex may appear in an associated homology group. This basis
element represents a hole or void. At some later point in the filtration, the
element representing the hole or void will disappear. We can record the
ordered pair of $\epsilon$-values at which the void is born and dies. Then,
for a particular $k$, we can take all birth-death pairs for persistent
homology groups of order $k$ and plot them on a 2-dimensional plane, called a
persistence diagram.

####  2.3.3 Persistent homology of the Lorenz attractor

We compared the homology groups of the Lorenz attractor to the persistent
homology groups of the autonomous and driven attractors. We followed the lead
of GARLAND201649 who computed the persistent homology of the Lorenz system
reconstructed from a sequence of 1D observations of a Lorenz trajectory using
the delay observation map described in Takens’ Theorem. The authors used the
open source software Javaplex created by Javaplex to find the Witness Complex
(outside the scope of this thesis) for the delay embedded Lorenz attractor and
computed the homology of the complex. They discuss a few subtleties that
arise, in particular that the Lorenz attractor is a fractal, whose structure
cannot be reconstructed exactly from any finite number of sample points. The
authors therefore satisfied themselves by approximating the Lorenz attractor
with a branched manifold model presented by PMIHES_1979__50__73_0 which has
the homology of the _double loop_ (which is topologically a drawing of the
number 8). We made the same approximation, and expected to find that the
application of persistent homology to the Lorenz system, driven ESN dynamics,
and autonomous ESN dynamics would reveal that all three have the double loop
homology groups. In particular the persistence diagrams of these three systems
would exhibit a pair of $H_{1}$ persistent homology groups floating well above
the diagonal. To verify this, we produced persistence diagrams using the open
source software Ripser produced by ctralie2018ripser and plotted the results
in Figure 2-8.

In the Figure we plot the $H_{1}$ persistence diagrams of the driven ESN
dynamics, autonomous ESN dynamics, and Lorenz dynamics as blue circles, red
downward triangles, and purple upward triangles. Each symbol represents a
homology group. The $x$-coordinate of the symbol represents the length
$\epsilon$ at which a homology group appears (is born) and the $y$-coordinate
represents the length at which the homology group vanishes (dies). All symbols
are above the diagonal because all homology groups die after they are born.
Symbols close to the diagonal represent homology groups that die shortly after
they are born, and are the result of noise, or an artefact of finite data.
Symbols far above the diagonal represent homology groups that persist long
after they are born and therefore represent persistent homology groups.

We can see a pair of blue circles, red downward triangles, and purple upward
triangles floating well above the diagonal, suggesting the existence of 2
persistent homology groups for the driven ESN dynamics, autonomous ESN
dynamics, and Lorenz dynamics. This is consistent with our expectation that
the driven ESN dynamics, autonomous ESN dynamics, and Lorenz dynamics all
adopt the topology of the double loop, which of course has 2 holes. The blue
circles, red downward triangles, and purple upward triangles are not exactly
in the same location because the driven ESN dynamics, autonomous ESN dynamics,
and Lorenz dynamics are not isometric.

The reader may wonder why we would use persistent homology to show that the
Lorenz system, driven ESN dynamics, and autonomous ESN dynamics all have the
homology of the double loop when this can clearly be seen in Figure 2-1. The
homology of a 3D system is usually apparent from a plot, but persistent
homology can reveal the holes, voids and higher dimensional hypervoids of high
dimensional systems that cannot be easily visualised. For example MULDOON19931
computed the homology of a delay embedded time series from a fluid dynamics
experiment, which could in general be of much higher dimension.

Figure 2-8: We have plotted the $H_{1}$ persistence diagrams of the driven ESN
dynamics, autonomous ESN dynamics, and Lorenz dynamics as blue circles, red
downward triangles, and purple upward triangles. If a homology group appears
(is born) at time $\epsilon_{0}>0$ and then vanishes (dies) at time
$\epsilon_{1}>\epsilon_{0}>0$ then we plot this homology group at the point
$(\epsilon_{0},\epsilon_{1})\in\mathbb{R}^{2}$ on a so-called persistence
diagram. Points that are close to the diagonal die shortly after they are
born, and are therefore transient homology groups, which we view as artefacts.
Points far from the diagonal die long after they are born and are therefore
persistent homology groups which we interpret as truly describing the topology
of the underlying object. We can see that each of these 3 objects has a pair
of points floating well above the diagonal, suggesting each has 2 holes. This
is consistent with our expectation that all three adopt the topology of the
double loop.

##  Chapter 3 Generalised Synchronisation (GS)

###  3.1 Background

In this Chapter we will examine the generalised synchronisation (GS) i.e. the
map $f$ from an underlying dynamical system $(M,\phi)$ to the reservoir space
$\mathbb{R}^{N}$. GSs appear in the seminal papers by PhysRevE.51.980,
PhysRevLett.76.1816, doi:10.1063/1.166278, stark_1999, BOCCALETTI20021, and a
recent paper by doi:10.1080/00107514.2017.1345844. Work on GS appears broadly
in mathematics and physics (PhysRevLett.64.821), (PhysRevA.44.2374),
(doi:10.1080/00107514.2017.1345844), and extensively in reservoir computing
(doi:10.1063/1.5022276), (Verzelli), (chaos_on_compacta),
(doi:10.1063/1.5128898).

For the purpose of this thesis we will define the GS (in discrete time) in
terms of a drive and response system.

######  Definition 3.1.1.

(Drive and response systems) Let $M$ be a smooth manifold and
$\phi\in\text{Diff}^{1}(M)$ a diffeomorphism. This is called the _drive
system_. Let $G:\mathbb{R}^{N}\times M\to\mathbb{R}^{N}$, and define the
system

$\displaystyle x_{k+1}=G(x_{k},\phi^{k}(m_{0})).$

This is called the _response system_.

A reservoir computer is essentially a drive-response system where the
underlying process $(M,\phi)$ is the drive system, and the evolution of the
reservoir states (the reservoir dynamics) is the response system. We are
interested in finding a map $f:M\to\mathbb{R}^{N}$ relating the drive to the
response, and this is called a generalised synchronisation (GS).

######  Definition 3.1.2.

(Generalised Synchronisation) Let $V\subseteq\mathbb{R}^{N}$. A generalised
synchronisation (GS) between the drive system and response system is a map
$f:M\to\mathbb{R}^{N}$ such that for any $x_{0}\in V\subseteq\mathbb{R}^{N}$
and any $m_{0}\in M$ we have $f\circ\phi^{k}(m_{0})=G(x_{k},\phi^{k}(m_{0}))$
for all $k\in\mathbb{Z}$. If $V$ is a proper subset of $\mathbb{R}^{N}$ then
$f$ is called a _local GS_ , and if $V=\mathbb{R}^{N}$ then $f$ is called a
_global GS_.

We will see in section 3.2 that the concept of GS is closely related the _Echo
State Property_ (ESP), which was coined by (Jaeger2001), and is discussed
extensively in the reservoir computing world. The name evokes the image of
past observations influencing the present moment with a fading intensity -
just as the intensity of an echo fades with time.

###  3.2 Echo State Property

The existence of the GS $f$ is closely related to the Echo State Property
(ESP) first introduced in the seminal paper by Jaeger2001. There are several
different formulations of the ESP, including in Jaeger2001, YILDIZ20121, and
GRIGORYEVA2018495. We will proceed with one inspired by YILDIZ20121 for now.
To be as general as possible we introduce an arbitrary reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ which reduces to an
ESN as a special case when

$\displaystyle F(x,z)=\sigma(Ax+Cz+b).$

######  Definition 3.2.1.

(Global Echo State Property) We say that a reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ has the _global ESP_
if, given any sequence of inputs $(z_{k}\in\mathbb{R}^{d})_{k\in\mathbb{N}}$,
and initial reservoir states $x_{0},y_{0}\in\mathbb{R}^{N}$ the sequences

$\displaystyle x_{k+1}=F(x_{k},z_{k})\qquad\text{and}\qquad
y_{k+1}=F(y_{k},z_{k})$

satisfy $\lVert x_{k+1}-y_{k+1}\rVert\to 0$ as $k\to\infty$.

We can ensure that a reservoir map $F$ has the global ESP if $F$ is globally
state contracting.

######  Definition 3.2.2.

A reservoir map $F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ is
called globally state contracting if there exists a $c\in(0,1)$ such that for
any $x,y\in\mathbb{R}^{N}$ and $z\in\mathbb{R}^{d}$ it follows

$\displaystyle\lVert F(x,z)-F(y,z)\rVert\leq c\lVert x-y\rVert.$

######  Theorem 3.2.3.

If a reservoir map $F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ is
globally state contracting, then $F$ has the global ESP.

###### Proof.

Consider that

$\displaystyle\lVert x_{k+1}-y_{k+1}\rVert=\lVert
F(x_{k},z_{k})-F(y_{k},z_{k})\rVert\leq c\lVert x_{k}-y_{k}\rVert$

so we can show by recursion that

$\displaystyle\lVert x_{k+1}-y_{k+1}\rVert\leq c^{k}\lVert x_{0}-y_{0}\rVert$

so $\lVert x_{k+1}-y_{k+1}\rVert\to 0$ as $k\to\infty$. ∎

With this, we are ready to introduce our first major result. Suppose the
inputs are a sequence of scalar observations from a dynamical system i.e.
$z_{k}=\omega\circ\phi^{k}(m_{0})$ where

  * •

$m_{0}\in M$ is an initial point on a smooth compact $q$-manifold $M$.

  * •

$\phi\in\text{Diff}^{1}(M)$ is a smooth diffeomorphism on $M$

  * •

$\omega\in C^{1}(M,\mathbb{R})$ is a smooth observation function.

Suppose further that $F$ is globally state contracting. Then, there is a
continuously differentiable GS $f\in C^{1}(M,\mathbb{R}^{N})$ synchronising
the drive system system (represented by $\phi$) to dynamics in the reservoir
space $\mathbb{R}^{N}$.

###  3.3 A theorem for global GS

Before we state the theorem, we will need some formalities. First we insist
that the manifold $M$ is compact, connected, second countable and Hausdorff so
that it can be endowed with a Riemannian metric $g$ (Riemannian_Geometry).
Then for any $f\in C^{1}(M,\mathbb{R}^{N})$ we define

$\displaystyle\lVert Df\rVert_{\infty}=\sup_{m\in M}\lVert
Df(m)\rVert_{2},\qquad\lVert Df(m)\rVert_{2}=\sup_{v\in
T_{m}M/\\{0\\}}\frac{\lVert Df(m)v\rVert_{2}}{\sqrt{g(m)(v,v)}}$

with $T_{m}M$ the tangent space of $M$ at $m$ and $Df(m)$ the differential of
$f$ at $m$. Analogously for $\phi\in\text{Diff}^{1}(M)$ we define

$\displaystyle\lVert T\phi\rVert_{\infty}=\sup_{m\in M}\lVert
T_{m}\phi\rVert,\qquad\lVert T_{m}\phi\rVert=\sup_{v\in
T_{m}M/\\{0\\}}\sqrt{\frac{g(\phi(m))(T_{m}\phi v,T_{m}\phi v)}{g(m)(v,v)}}$

where $T_{m}\phi$ is the tangent map for $\phi$ at $m$, and $(v,w)$ is the
inner product on $T_{m}$.

######  Theorem 3.3.1.

(Grigoryeva et al. (2021a)) Let $M$ be a compact, smooth, $q$-manifold (i.e a
$q$-dimensional manifold) and $\phi\in\text{Diff}^{1}(M)$. Let $\omega\in
C^{1}(M,\mathbb{R}^{d})$ be the observation function on $M$. Suppose that the
reservoir map $F\in C^{2}(\mathbb{R}^{N}\times\mathbb{R}^{d},\mathbb{R}^{N})$
is globally state contracting. Furthermore, suppose that

$\displaystyle c_{x}:=\sup_{(x,z)\in\mathbb{R}^{N}\times\mathbb{R}^{d}}\lVert
D_{x}F(x,z)\rVert_{2}$

(where $D_{x}$ is the partial derivative with respect to $x$) satisfies

$\displaystyle c_{x}<1,\qquad c_{x}\lVert T\phi^{-1}\rVert_{\infty}<1.$

Then, there exists a unique solution $f\in C^{1}(M,\mathbb{R}^{N})$ of the
equation

$\displaystyle f=F(f\circ\phi^{-1},\omega)$

such that for all initial $x_{0}\in\mathbb{R}^{N}$ and $m_{0}\in M$ the
sequence

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))$

converges to $f\circ\phi^{k}(m_{0})$ as $k\to\infty$.

The proof of Theorem 3.3.1 is omitted because the theorem is a special case of
Theorem 3.5.5 (obtained by setting $V=\mathbb{R}^{N}$) which we will state and
prove in section 3.5. The more general result proved in section 3.5 concerns
the case where the reservoir map $F$ admits multiple generalised
synchronisations. To build some intuition for this, we will proceed with an
example. Our example is closely related to the work by GREBOGI1984261 which
studies _strange attractors that are not chaotic_.

###  3.4 Motivating example

Let $F:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ be defined by

$\displaystyle F(x,z)=\tanh(2x+z).$

Under the constant driving input $z_{k}=0$ for all $k$, we have the autonomous
system $F(x,0)=\tanh(2x)$, which has 3 fixed points at $-x^{*},0,x^{*}$ where
$x^{*}\approx 0.958$ illustrated in Figure 3-1 below. $-x^{*}$ and $x^{*}$ are
stable nodes, while $0$ is an unstable node.

Doctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyUntable NodeDoctor of PhilosophyStable
NodeDoctor of PhilosophyStable NodeDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of Philosophy Figure 3-1: A graph of $y=\tanh(2x)$ intersecting $y=x$ at 3
points. These 3 points are fixed points of the map $F(x,0)=\tanh(2x)$.

Multiple equilibria prove that $F$ does not have the global ESP for this
input. Now, let us consider the behavior of $F$ under a more interesting
driving input. Let $\phi\in\text{Diff}^{1}(S^{1})$ be a diffeomorphism on the
circle defined by $\phi(m)=m+\epsilon$ for some fixed $\epsilon>0$. Let
$\omega\in C^{1}(S^{1},\mathbb{R})$ be defined to be
$\omega(m)=\frac{1}{2}\sin(m)$. Now suppose the driving input
$u_{k}=\omega\circ\phi^{k}(m_{0})$ for some initial $m_{0}$. Now, any GS
$f:S^{1}\to\mathbb{R}$ satisfies, $\forall m\in S^{1}$,

$\displaystyle f(m)$ $\displaystyle=F(f\circ\phi^{-1}(m),\omega(m))$
$\displaystyle=\tanh\bigg{(}2f(m+\epsilon)+\frac{1}{2}\sin(m)\bigg{)}.$ (3.1)

Suppose we choose $\epsilon=2\pi/100$, and take an initial
$x_{0}\in\mathbb{R}$, then compute $x_{k+1}=F(x_{k},z_{k})$, and plot the
pairs $(\phi^{k}(m_{0}),x_{k})\in\mathbb{R}^{2}$. It appears from numerical
experiments, that for any initial $x_{0}$ in the vicinity of the fixed point
$-x^{*}$, after a short transient period, we arrive at the same graph $G_{-}$
of the pairs $(\phi^{k}(m_{0}),x_{k})\in\mathbb{R}^{2}$. Furthermore, the
graph plausibly satisfies the relation (3.1) of being a GS. Intriguingly, if
we take any $x_{0}$ in the vicinity of $x^{*}$ this time, then after a short
transient period, we arrive at a different graph $G_{+}$ of the pairs
$(\phi^{k}(m_{0}),x_{k})\in\mathbb{R}^{2}$. This suggests the existence of (at
least) two distinct (locally) stable GSs with respect to the input $z_{k}$.
The graphs $G_{-}$ and $G_{+}$ are shown in Figure 3-2.

(a) $G_{-}$

(b) $G_{+}$

Figure 3-2: The graphs $G_{-}$ and $G_{+}$ of the points
$(\phi^{k}(m_{0}),x_{k})$.

###  3.5 Theorems for local GS

This numerical experiment suggests that there exist reservoir maps $F$ that do
not have global ESP, and admit multiple GSs. To start to make sense of this,
we will introduce the concept of local ESP with admits the global ESP as a
special case.

######  Definition 3.5.1.

(Local ESP) Let $V\subset\mathbb{R}^{N}$ be closed under the Euclidean
topology. We say that a reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ has the $(V,W)$-local
ESP if, given any sequence of inputs $(z_{k}\in
W\subset\mathbb{R}^{d})_{k\in\mathbb{N}}$, and initial reservoir states states
$x_{0},y_{0}\in V$ the sequences

$\displaystyle x_{k+1}=F(x_{k},z_{k})\qquad\text{and}\qquad
y_{k+1}=F(y_{k},z_{k})$

are contained in $V$ and $\lVert x_{k+1}-y_{k+1}\rVert\to 0$ as $k\to\infty$.

We observe first of all that it is possible for $F$ to have $(V_{i},W)$-local
ESP on disjoint sets $V_{i}$ with $i=1,...,n$. We can also ensure that a
reservoir map $F$ has the $(V,W)$-local ESP by demanding that $F$ has the
$(V,W)$-local state contracting property.

######  Definition 3.5.2.

(Local state contracting) A reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ is called
$(V,W)$-locally state contracting if there exists a $c\in(0,1)$ such that for
any $x,y\in V$ and $z\in W\subset\mathbb{R}^{d}$ it follows that

  * •

$\lVert F(x,z)-F(y,z)\rVert\leq c\lVert x-y\rVert$

  * •

$F(x,z)\in V$.

######  Theorem 3.5.3.

If a reservoir map $F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ is
$(V,W)$-locally state contracting, then $F$ has the $(V,W)$-local ESP.

###### Proof.

First of all, since $F(x,z)\in V$, it follows that the sequence

$\displaystyle x_{k+1}=F(x_{k},z_{k})$

is contained by $V$. Next, for any $x_{0},y_{0}\in V$ we have

$\displaystyle\lVert x_{k+1}-y_{k+1}\rVert=\lVert
F(x_{k},z_{k})-F(y_{k},z_{k})\rVert\leq c\lVert x_{k}-y_{k}\rVert$

so we can show by recursion that

$\displaystyle\lVert x_{k+1}-y_{k+1}\rVert\leq c^{k}\lVert x_{0}-y_{0}\rVert$

so $\lVert x_{k+1}-y_{k+1}\rVert\to 0$ as $k\to\infty$. ∎

Now we are ready to state two theorems for the existence of local GS. The
first result establishes the existence of a continuous GS $f$ on a topological
space, while the second establishes a differentiable GS $f$ on a smooth
manifold. A less general version of these results (with a less rigorous proof)
appears in embedding_and_approximation_theorems (Theorem 2.2.2) and a fully
rigorous more general version appears in Grigoryeva et al. (2021a) (Theorem
III).

######  Theorem 3.5.4.

(Grigoryeva et al. (2021a)) Let $M$ be a compact topological space and
$\phi\in\text{Hom}(M)$. Let $\omega\in C^{0}(M,\mathbb{R}^{d})$ be the
observation function on $M$. Let $V\subset\mathbb{R}^{N}$ be closed and
suppose that the reservoir map $F\in C^{0}(V\times\omega(M),\mathbb{R}^{N})$
is $(V,\omega(M))$-locally state contracting with contraction coefficient $c$.
Then, there exists a unique solution $f\in C^{0}(M,V)$ of the equation

$\displaystyle f=F(f\circ\phi^{-1},\omega)$

such that for all initial $x_{0}\in V$ and $m_{0}\in M$ the sequence

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))$

converges to $f\circ\phi^{k}(m_{0})$ as $k\to\infty$.

###### Proof.

Define the map $\Psi:C^{0}(M,V)\to C^{0}(M,V)$ by

$\displaystyle\Psi(h)=F(h\circ\phi^{-1},\omega).$

Notice that

$\displaystyle\lVert\Psi(h)-\Psi(g)\rVert_{\infty}$ $\displaystyle=\lVert
F(h\circ\phi^{-1},\omega)-F(g\circ\phi^{-1},\omega)\rVert_{\infty}$
$\displaystyle\leq c\lVert
h\circ\phi^{-1}-g\circ\phi^{-1}\rVert_{\infty}=c\lVert h-g\rVert_{\infty}$

by the definition of $c$. So $\Psi$ is a contraction mapping on the Banach
space $C^{0}(M,V)$ and therefore (by the Banach1922 Fixed Point Theorem)
admits a unique fixed point $f\in C^{0}(M,V)$ such that

$\displaystyle\Psi(f)=F(f\circ\phi^{-1},\omega)=f.$

Furthermore, $F$ is $(V,\omega(M))$-locally contracting and therefore has the
$(V,\omega(M))$-local ESP, so for all initial $x_{0}\in V$ and $m_{0}\in M$
and $y_{0}=f\circ\phi^{-1}(m_{0})$ the sequences

$\displaystyle
x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))\qquad\text{and}\qquad
y_{k+1}=F(y_{k},\omega\circ\phi^{k}(m_{0}))$

satisfy $\lVert x_{k+1}-y_{k+1}\rVert\to 0$ as $k\to\infty$. Since
$y_{k+1}=f\circ\phi^{k}(m_{0})$ it follows that $x_{k+1}\to
f\circ\phi^{k}(m_{0})$ and this completes the proof. ∎

######  Theorem 3.5.5.

(Grigoryeva et al. (2021a)) Let $M$ be a compact, smooth, $q$-manifold and
$\phi\in\text{Diff}^{1}(M)$. Let $\omega\in C^{1}(M,\mathbb{R}^{d})$ be the
observation function on $M$. Let $V\subset\mathbb{R}^{N}$ be closed and
convex. Suppose that the reservoir map $F\in
C^{2}(V\times\omega(M),\mathbb{R}^{N})$ is $(V,\omega(M))$-locally state
contracting. Furthermore, suppose that

$\displaystyle c_{x}:=\sup_{(x,z)\in V\times\omega(M)}\lVert
D_{x}F(x,z)\rVert_{2}$

(where $D_{x}$ is the partial derivative with respect to $x$) satisfies

$\displaystyle c_{x}<1,\qquad and\qquad c_{x}\lVert
T\phi^{-1}\rVert_{\infty}<1.$

Then, there exists a unique solution $f\in C^{1}(M,V)$ of the equation

$\displaystyle f=F(f\circ\phi^{-1},\omega)$

such that for all initial $x_{0}\in V$ and $m_{0}\in M$ the sequence

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))$

converges to $f\circ\phi^{k}(m_{0})$ as $k\to\infty$.

The proof of the second result is rather more technical, but relies on a
similar argument. We will need to first introduce for each $\delta>0$ a norm
$\lVert\cdot\rVert_{C^{1}(\delta)}$ on $C^{1}(M,V)$ defined by

$\displaystyle\lVert f\rVert_{C^{1}(\delta)}:=\lVert
f\rVert_{\infty}+\delta\lVert Df\rVert_{\infty}.$

It is shown in Chapter 2 of https://doi.org/10.1002/zamm.19690490315 that each
of these endows $C^{1}(M,V)$ with a Banach space structure. Furthermore,
https://doi.org/10.1002/zamm.19690490315 show that all of these norms induce
the same topology, and that these topologies coincide with the $C^{1}$
topology described in Invariant_Manifolds. With this, we will now prove 2
technical lemmas which are used to prove Theorem 3.5.5.

######  Lemma 3.5.6.

Let $r>0$ and $W\subset\mathbb{R}^{N}$ be closed in the Euclidean topology.
Then, for any $\delta>0$, the set

$\displaystyle\Omega(r,W)=\\{f\in C^{1}(M,W)\ |\ \lVert Df\rVert_{\infty}\leq
r\\}$

is closed in the topology induced by the $\lVert\cdot\rVert_{C^{1}(\delta)}$
norm.

###### Proof.

$\Omega(r,W)=\Omega(r)\cap C^{1}(M,W)$ where

$\displaystyle\Omega(r)=\\{f\in C^{1}(M,\mathbb{R}^{N})\ |\ \lVert
Df\rVert_{\infty}\leq r\\}$

so it suffices to the show that both $\Omega(r)$ and $C^{1}(M,W)$ are closed.
We will show $\Omega(r)$ is closed by proving that the complement
$C^{1}(M,\mathbb{R}^{N})\backslash\Omega(r)$ is open. We will show that for an
arbitrary point in the complement $f\in
C^{1}(M,\mathbb{R}^{N})\backslash\Omega(r)$ there exists an $\epsilon>0$ such
that the ball of radius $\epsilon$ centred at $f$ denoted
$B_{C^{1}(\delta)}(f,\epsilon)$ is a subset of the complement
$C^{1}(M,\mathbb{R}^{N})\backslash\Omega(r)$. We fix

$\displaystyle\epsilon=\delta(\lVert Df\rVert_{\infty}-r)$

and then for any $g\in C^{1}(M,\mathbb{R}^{N})/\Omega(r)$

$\displaystyle\lVert Df\rVert_{\infty}$ $\displaystyle=\lVert Df+Dg-
Dg\rVert_{\infty}$ $\displaystyle\leq\lVert Df-Dg\rVert_{\infty}+\lVert
Dg\rVert_{\infty}$ $\displaystyle\leq\lVert Df-
Dg\rVert_{\infty}+\frac{1}{\delta}\lVert f-g\rVert_{\infty}+\lVert
Dg\rVert_{\infty}$ $\displaystyle=\frac{1}{\delta}\lVert
f-g\rVert_{C^{1}(\delta)}+\lVert Dg\rVert_{\infty}$
$\displaystyle<\frac{\epsilon}{\delta}+\lVert Dg\rVert_{\infty}$
$\displaystyle=\lVert Df\rVert_{\infty}-r+\lVert Dg\rVert_{\infty}$

using the definition of the $C^{1}(\delta)$ norm directly. So $\lVert
Dg\rVert_{\infty}>r$ hence $B_{C^{1}(\delta)}(f,\epsilon)\subset
C^{1}(M,\mathbb{R}^{N})\backslash\Omega(r)$. Next, we will show that
$C^{1}(M,W)$ is closed by showing that for an arbitrary sequence $(f_{n}\in
C^{1}(M,W))_{n\in\mathbb{N}}$ that converges to $f^{*}\in
C^{1}(M,\mathbb{R}^{N})$ it follows that $f^{*}\in C^{1}(M,W)$. If $(f_{n}\in
C^{1}(M,W))_{n\in\mathbb{N}}$ converges to $f^{*}\in C^{1}(M,\mathbb{R}^{N})$
then for all $\epsilon>0$ there exists $N(\epsilon)\in\mathbb{N}$ such that,
for any $m\in M$

$\displaystyle\rVert f(m)-f^{*}(m)\lVert$ $\displaystyle\leq\rVert
f-f\lVert_{\infty}$ $\displaystyle\leq\rVert
f-f^{*}\lVert_{\infty}+\delta\lVert Df-Df^{*}\rVert_{\infty}$
$\displaystyle=\lVert f-f^{*}\rVert_{C^{1}(\delta)}<\epsilon$

so $f(m)\to f^{*}(m)$ and $W$ is closed so $f^{*}(m)\in W$. This holds for
arbitrary $m\in M$ so it follows that $f^{*}\in C^{1}(M,W)$. ∎

In what follows we will use the notation $D_{x}F(x,z)$, $D_{z}F(x,z)$ to
denote the partial derivatives of the reservoir map $F$ with respect to $x$
and $z$, respectively.

######  Lemma 3.5.7.

Let $F\in C^{1}(V\times\omega(M),\mathbb{R}^{N})$ be a $(V,\omega(M))$-locally
contracting reservoir map. Let

$\displaystyle c_{x}:=\sup_{(x,z)\in V\times\omega(M)}\lVert
D_{x}F(x,z)\rVert_{2},\qquad c_{z}:=\sup_{(x,z)\in V\times\omega(M)}\lVert
D_{z}F(x,z)\rVert_{2},$

and and choose $r$ such that

$\displaystyle r>\frac{c_{z}\lVert D\omega\rVert_{\infty}}{1-c_{x}\lVert
T\phi^{-1}\rVert_{\infty}}.$ (3.2)

Let $\Psi:\Omega(r,\omega(M))\to C^{1}(M,\mathbb{R}^{N})$ be defined by

$\displaystyle\Psi(h)=F(h\circ\phi^{-1},\omega).$

Then $\Psi(\Omega(r,\omega(M)))\subset\Omega(r,\omega(M))$ and hence the map
$\Psi:\Omega(r,\omega(M))\to\Omega(r,\omega(M))$ is well defined

###### Proof.

Let $h\in\Omega(r,\omega(M))$. It suffices to show that
$\Psi(h)(m)\in\omega(M)$ for all $m\in M$ and that $\lVert
D(\Psi(h))\rVert_{\infty}\leq r$. We easily verify that
$\Psi(h)(m)\in\omega(M)$ for all $m\in M$ because $F$ is
$(V,\omega(M))$-locally contracting. Now, by direct computation of the
derivative

$\displaystyle\lVert D(\Psi(h))\rVert_{\infty}=\lVert
DF(h\circ\phi^{-1},\omega)\rVert_{\infty}$ $\displaystyle=\sup_{m\in M}\lVert
D_{x}F(h\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}+D_{z}F(h\circ\phi^{-1}(m),\omega(m))D\omega(m)\rVert_{2}$
$\displaystyle\leq c_{x}\lVert Dh\rVert_{\infty}\lVert
T\phi^{-1}\rVert_{\infty}+c_{z}\lVert D\omega\rVert_{\infty}$
$\displaystyle=c_{x}r\lVert T\phi^{-1}\rVert_{\infty}+c_{z}\lVert
D\omega\rVert_{\infty}<r.\text{ (by \eqref{r_condition})}$

∎

We are almost ready to prove the result. For what follows we let
$D_{xx}F(x,z)$ and $D_{xz}F(x,z)$ denote the second partial derivatives of $F$
and let

$\displaystyle c_{xx}:=\sup_{(x,z)\in V\times\omega(M)}\lVert
D_{xx}F(x,z)\rVert_{2},\qquad\textrm{and}\qquad c_{xz}:=\sup_{(x,z)\in
V\times\omega(M)}\lVert D_{xz}F(x,z)\rVert_{2}.$

Then it follows from the mean value theorem, which applies because $V$ is
convex, that for any $x,y\in V$ and $z\in\omega(M)$ that

$\displaystyle\lVert D_{x}(F(x,z)-F(y,z))\rVert\leq c_{xx}\lVert
x-y\rVert,\qquad\textrm{and}\qquad\lVert D_{z}(F(x,z)-F(y,z))\rVert\leq
c_{xz}\lVert x-y\rVert.$

###### Proof.

(Proof of Theorem 3.5.5) We have shown that
$\Psi:\Omega(r,\omega(M))\to\Omega(r,\omega(M))$ is a well defined map on a
closed subset of the Banach space $C^{1}(M,\mathbb{R}^{N})$. We will now show
that for sufficiently small $\delta_{0}>0$ the map $\Psi$ is a contraction in
the $\lVert\cdot\rVert_{C^{1}(\delta_{0})}$ norm, and therefore admits a
unique fixed point $f\in\Omega(r,\omega(M))\subset C^{1}(M,\mathbb{R}^{N})$ by
the Banach1922 Fixed Point Theorem. The fixed point $f$ satisfies

$\displaystyle\Psi(f)=F(f\circ\phi^{-1},\omega)=f.$

Furthermore, $F$ is $(V,\omega(M))$-locally contracting and therefore has the
$(V,\omega(M))$-local ESP, so for all initial $x_{0}\in V$ and $m_{0}\in M$
and $y_{0}=f\circ\phi^{-1}(m_{0})$ the sequences

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))\qquad
y_{k+1}=F(y_{k},\omega\circ\phi^{k}(m_{0}))$

satisfy $x_{k+1}-y_{k+1}\to 0$ as $k\to\infty$. Now
$y_{k+1}=f\circ\phi^{k}(m_{0})$ so $x_{k+1}\to f\circ\phi^{k}(m_{0})$. All
that remains to show is that $\Psi$ is a contraction in
$\lVert\cdot\rVert_{C^{1}(\delta_{0})}$ for sufficiently small $\delta_{0}$.
We do this by considering the two parts of the $C^{1}(\delta_{0})$ norm
separately. First, we observe that

$\displaystyle\lVert\Psi(h)-\Psi(g)\rVert_{\infty}$ $\displaystyle=\lVert
F(h\circ\phi^{-1},\omega)-F(g\circ\phi^{-1},\omega)\rVert_{\infty}$
$\displaystyle\leq c\lVert h\circ\phi^{-1}-g\circ\phi^{-1})\rVert_{\infty}$
$\displaystyle=c\lVert h-g\rVert_{\infty}.$

Second, we compute directly that

$\displaystyle\qquad\lVert D\Psi(g)-D\Psi(h)\rVert_{\infty}$
$\displaystyle=\lVert
DF(g\circ\phi^{-1},\omega)-DF(h\circ\phi^{-1},\omega)\rVert_{\infty}$
$\displaystyle=\sup_{m\in M}\lVert
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dg)\circ\phi^{-1}(m)T_{m}\phi^{-1}+D_{z}F(g\circ\phi^{-1}(m),\omega(m))D\omega(m)$
$\displaystyle\qquad-
D_{x}F(h\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi{-1}(m)T_{m}\phi^{-1}-D_{z}F(h\circ\phi^{-1}(m),\omega(m))D\omega(m)\rVert_{2}$
$\displaystyle\leq\sup_{m\in M}\lVert
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dg)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad-
D_{x}F(h\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}\rVert_{2}$
$\displaystyle\qquad+\sup_{m\in M}\lVert
D_{z}F(g\circ\phi^{-1}(m),\omega(m))D\omega(m)-D_{z}F(h\circ\phi^{-1}(m),\omega(m))D\omega(m)\rVert_{2}$
$\displaystyle\leq\sup_{m\in M}\lVert
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dg)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad-
D_{x}F(h\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}\rVert_{2}$
$\displaystyle\qquad+\lVert
D_{z}(F(g\circ\phi^{-1},\omega)-F(h\circ\phi^{-1},\omega))\rVert_{\infty}\lVert
D\omega\rVert_{\infty}$ $\displaystyle\leq\sup_{m\in M}\lVert
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dg)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad-
D_{x}F(h\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}\rVert_{2}$
$\displaystyle\qquad+c_{xz}\lVert g-h\rVert_{\infty}\lVert
D\omega\rVert_{\infty}$ $\displaystyle=\sup_{m\in M}\lVert
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dg)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad-
D_{x}F(h\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad+D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad-
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}\rVert_{2}$
$\displaystyle\qquad+c_{xz}\lVert g-h\rVert_{\infty}\lVert
D\omega\rVert_{\infty}$ $\displaystyle\leq\sup_{m\in M}\lVert
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dg)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad-
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}\rVert_{2}$
$\displaystyle\qquad\sup_{m\in M}\lVert
D_{x}F(h\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}$
$\displaystyle\qquad-
D_{x}F(g\circ\phi^{-1}(m),\omega(m))(Dh)\circ\phi^{-1}(m)T_{m}\phi^{-1}\rVert_{2}$
$\displaystyle\qquad+c_{xz}\lVert g-h\rVert_{\infty}\lVert
D\omega\rVert_{\infty}$ $\displaystyle\leq\lVert
D_{x}F(g\circ\phi^{-1},\omega)\rVert_{\infty}\lVert Dg-Dh\rVert_{\infty}\lVert
T\phi^{-1}\rVert_{\infty}$ $\displaystyle\qquad+\lVert
D_{x}(F(h\circ\phi^{-1},\omega)-F(g\circ\phi^{-1},\omega))\rVert_{\infty}\lVert
Dh\rVert_{\infty}\lVert T\phi^{-1}\rVert_{\infty}$
$\displaystyle\qquad+c_{xz}\lVert g-h\rVert_{\infty}\lVert
D\omega\rVert_{\infty}$ $\displaystyle\leq c_{x}\lVert Dg-
Dh\rVert_{\infty}\lVert T\phi^{-1}\rVert_{\infty}$
$\displaystyle\qquad+c_{xx}\lVert h-g\rVert_{\infty}\lVert
Dh\rVert_{\infty}\lVert T\phi^{-1}\rVert_{\infty}$
$\displaystyle\qquad+c_{xz}\lVert g-h\rVert_{\infty}\lVert
D\omega\rVert_{\infty}$ $\displaystyle\leq c_{x}\lVert Dg-
Dh\rVert_{\infty}\lVert T\phi^{-1}\rVert_{\infty}$
$\displaystyle\qquad+c_{xx}\lVert h-g\rVert_{\infty}r\lVert
T\phi^{-1}\rVert_{\infty}$ $\displaystyle\qquad+c_{xz}\lVert
g-h\rVert_{\infty}\lVert D\omega\rVert_{\infty}$ $\displaystyle=(rc_{xx}\lVert
T\phi^{-1}\rVert_{\infty}+c_{xz}\lVert D\omega\rVert_{\infty})\lVert
g-h\rVert_{\infty}+c_{x}\lVert T\phi^{-1}\rVert_{\infty}\lVert Dg-
Dh\rVert_{\infty}.$

Now we choose

$\displaystyle\delta_{0}<\frac{1-c}{rc_{xx}\lVert
T\phi^{-1}\rVert_{\infty}+c_{xz}\lVert D\omega\rVert_{\infty}}$

so that, combining the two estimates above we obtain

$\displaystyle\lVert\Psi(g)-\Psi(h)\rVert_{C^{1}(\delta_{0})}$
$\displaystyle=\lVert\Psi(g)-\Psi(h)\rVert_{\infty}+\delta_{0}\lVert
D\Psi(g)-D\Psi(h)\rVert_{\infty}$ $\displaystyle=(c+\delta_{0}(rc_{xx}\lVert
T\phi^{-1}\rVert_{\infty}+c_{xz}\lVert D\omega\rVert_{\infty}))\lVert
g-h\rVert_{\infty}$ $\displaystyle\qquad+\delta_{0}c_{x}\lVert
T\phi^{-1}\rVert_{\infty}\lVert Dg-Dh\rVert_{\infty}$
$\displaystyle\leq(c+\delta_{0}(rc_{xx}\lVert
T\phi^{-1}\rVert_{\infty}+c_{xz}\lVert D\omega\rVert_{\infty}))\lVert
g-h\rVert_{C^{1}(\delta_{0})}.$

By our choice of $\delta_{0}$ the constant $(c+\delta_{0}(rc_{xx}\lVert
T\phi^{-1}\rVert_{\infty}+c_{xz}\lVert D\omega\rVert_{\infty}))\in(0,1)$, so
we have shown that $\Psi$ is a contracting map and the proof is complete. ∎

###  3.6 Curious examples

To create some further intuition for Theorem 3.5.5, we present four more
examples of reservoir maps $F$ which are locally state contracting and admit
local GSs. These examples are related to the work by GREBOGI1984261. In the
first example (which appears in Grigoryeva et al. (2021a)) we consider the
Lorenz system (1.1) once again. This time we feed the $\xi$ observations into
a reservoir map that admits 8 distinct GSs.

##### 1\. Finitely many GSs

We simulated a 4000 point (40 time units) trajectory of the Lorenz system
originating from the initial point $m_{0}=(0,1,1.05)$ with time-step $h=0.01$.
Figure 3-3 shows this trajectory for times $t\in(20,40)$.

Figure 3-3: A trajectory of the Lorenz system. $(u,v,w)$ correspond to the
variables $(\xi,\upsilon,\zeta)$ used in equations (1.1).

If we observe only the $\xi$-component of this trajectory, then the
observation function is $\omega(\xi,\upsilon,\zeta)=\xi$. The corresponding
observed trajectory is illustrated in Figure 3-4.

Figure 3-4: $\xi$-component of a trajectory of the Lorenz system. The
horizontal axis is time $t\in(20,40)$.

We now define the reservoir map
$F:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}^{3}$ given by

$\displaystyle
F(u,v,w;z)=(\mathrm{sgn}(u)|u|^{\alpha},\mathrm{sgn}(v)|v|^{\alpha},\mathrm{sgn}(w)|w|^{\alpha})+\lambda(\sin(kz),\cos(kz),\sin^{2}(kz))$

with $\lambda,k\geq 0$ and $\alpha\in(0,1)$. If we set $\lambda=0$ then the
reservoir map is autonomous and has 8 stable fixed points at $(\pm 1,\pm 1,\pm
1)$ taking all sign combinations. A cross section of the phase portrait of
this autonomous system at the $w=1$ plane is shown in Figure 3-5.

Figure 3-5:  A phase portrait of the autonomous system
$F(u,v,w;z)=(u^{\alpha},v^{\alpha},w^{\alpha})$ at the cross section $w=1$.
The four stable fixed points in this cross section are shown.

If we now choose $\lambda=0.009$, $k=0.1$ and $\alpha=0.9$, then we can
construct a box $V_{1}=[0.9,1.1]\times[0.9,1.1]\times[0.9,1.1]$ containing the
fixed point $(1,1,1)$. We can create a similar box about each of the other
seven fixed points and denote them $V_{2},V_{3},\ldots,V_{8}$. We have by
construction for each $i=1,\ldots,8$ that
$F(V_{i},\omega(\mathbb{R}^{3}))\subset V_{i}$ and that $F$ is state
contracting on each box $V_{i}$. Thus, by Theorem 3.5.5 the reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ is
$(V_{i},\omega(M))$-locally state contracting for each $i$, so there exist GSs
$f_{i}\in C^{1}(M,V_{i})$ for each $i$.

To observe the images of the GSs $f_{i}\in C^{1}(M,V_{i})$ we computed the
states

$\displaystyle x_{k+1}=F(x_{k},\omega(\phi^{k}(m_{0}))),\qquad
y_{k+1}=F(y_{k},\omega(\phi^{k}(m_{0}))),$ (3.3)

from two different initial states $x_{0}=(1,1,1)\in V_{1}$ and
$y_{0}=(-1,1,1)\in V_{2}$. For any other initial points $x_{0}^{\prime}\in
V_{1}$ and $y_{0}^{\prime}\in V_{2}$ we expect that the sequences (3.3) will
converge to $f_{1}\circ\phi^{k}(m_{0})$ and $f_{2}\circ\phi^{k}(m_{0})$
respectively. Figure 3-6 illustrates the states (3.3) after a burn-in time
$t=(0,20)$ so that the state sequences very closely approximate
$f_{1}\circ\phi^{k}(m_{0})$ and $f_{2}\circ\phi^{k}(m_{0})$ respectively.

(a) Image of the GS that contains the point $(1,1,1)$ in its image.

(b) Image of the GS that contains the point $(-1,1,1)$ in its image.

Figure 3-6: Image of the Lorenz solution under two different GSs.

##### 2\. Uncountably many GSs

Suppose $F:\mathbb{R}^{2}\times\mathbb{R}\to\mathbb{R}^{2}$ is defined by

$\displaystyle F(\rho,\theta;z)=(\sqrt{\rho}+z,\theta+\delta)$

where $(\rho,\theta)$ is the polar coordinate representation of the state
vector $x\in\mathbb{R}^{2}$, and $\delta\in\mathbb{R}$. Under the constant
driving input $z_{k}=0$ the state space system has a unstable node at the
origin, and admits an attracting invariant set $\\{\rho=1\\}$. A phase
portrait of the dynamics is shown in Figure 3-7(a).

(a) $F(\rho,\theta;z)=(\sqrt{\rho}+z,\theta+\delta)$.

(b) $F(\rho,\theta;z)=(\rho^{2}+z,\theta+\delta)$.

Figure 3-7: Phase portraits of 2 reservoir maps
$F:\mathbb{R}^{2}\times\mathbb{R}\to\mathbb{R}^{2}$ under the constant driving
input $z_{k}=0$ for all $k$.

Once again, suppose the driving input arises from
$\phi\in\text{Diff}^{1}(S^{1})$ a diffeomorphism on the circle defined by
$\phi(m)=m+\epsilon$ for $\epsilon>0$. This time, let the observation function
$\omega\in C^{1}(S^{1},\mathbb{R})$ be defined $\omega(m)=\frac{1}{5}\sin(m)$.
Suppose first of all that $\delta=0$. For each $\theta_{0}\in[-\pi,\pi)$ let

$\displaystyle V(\theta_{0})=\\{(\rho,\theta)\ |\ \rho\geq 1/2\
,\theta=\theta_{0}\\}.$

Note that $V(\theta_{0})\subset\mathbb{R}^{2}$ is closed and that for all
$x,y\in V(\theta_{0})$ and $z\in\omega(S^{1})$

  1. 1.

$\exists\,c\in(0,1)$ such that $\lVert F(x,z)-F(y,z)\rVert\leq c\lVert
x-y\rVert$

  2. 2.

$F(x,z)\in V(\theta_{0})$

so we have the $(V(\theta_{0}),\omega(S^{1}))$-local state contraction
property. Hence, by Theorem 3.5.5, there exists a unique solution
$f_{\theta_{0}}\in C^{1}(S^{1},V(\theta_{0}))$ of the equation

$\displaystyle f_{\theta_{0}}=F(f_{\theta_{0}}\circ\phi^{-1},\omega)$

such that for all initial $x_{0}\in V_{\theta_{0}}$ and $m_{0}\in S^{1}$ the
sequence

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))$

converges to $f_{\theta_{0}}\circ\phi^{k}(m_{0})$ as $k\to\infty$. This holds
for each $\theta_{0}\in[\pi,\pi)$, so there is an uncountable set of GSs
parametrised by $\theta_{0}$.

##### 3\. A bounded trajectory that does not admit a GS

Let $n\in\mathbb{N}$, $\epsilon=2\pi/n$ and $\delta=1$. Notice
$\phi^{n}(m_{0})=m_{0}$. If $x_{0}$ represents the polar angle of an initial
reservoir state $x_{0}$, and proceeding states are defined in the usual manner

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))$

then the polar angle of the $n$th reservoir state $x_{n}$ is
$\theta_{n}=\theta_{0}+n\delta$.

We will show that no GS $f$ exists for initial points on $\\{\rho=1\\}$.
Suppose for contradiction that such a GS does exist. Then
$f(m_{0})=f\circ\phi^{n}(m_{0})$ for all $m_{0}$, so the polar angles
$\theta_{0}$ and $\theta_{n}$ of $f(m_{0})$, and $f\circ\phi^{n}(m_{0})$ must
be equal. This leaves us with a contradiction, because $n\delta$ is not a
multiple of $2\pi$. We can see however that the trajectory
$\\{x_{k}\\}_{k\in\mathbb{N}}$ remains bounded.

##### 4\. Some initial conditions converge to the image of a GS

Suppose $F:\mathbb{R}^{2}\times\mathbb{R}\to\mathbb{R}^{2}$ is defined by

$\displaystyle F(\rho,\theta;z)=(\rho^{2}+z,\theta+\delta)$

where $\delta>0$. A phase portrait of this system under the constant driving
input $z_{k}=0$ for all $k$ is shown in Figure 3-7(b). Once more, suppose the
driving input arises from $\phi\in\text{Diff}^{1}(S^{1})$ a diffeomorphism on
the circle defined by $\phi(m)=m+\epsilon$ for $\epsilon>0$, and let the
observation function $\omega\in C^{1}(S^{1},\mathbb{R})$ be defined by
$\omega(m)=\frac{1}{5}\sin(m)$. Let

$\displaystyle V=\\{(\rho,\theta)\ |\ \rho\leq 1/2\\}.$

Note that $V\subset\mathbb{R}^{2}$ is closed and that for all $x,y\in V$ and
$z\in\omega(S^{1})$

  1. 1.

$\exists\,c\in(0,1)$ such that $\lVert F(x;z)-F(y;z)\rVert\leq c\lVert
x-y\rVert$,

  2. 2.

$F(x;z)\in V$,

so we have the $(V,\omega(M))$-local state contracting property. Hence, by
Theorem 3.5.5, there exists a unique solution $f\in C^{1}(S^{1},V)$ of the
equation

$\displaystyle f=F(f\circ\phi^{-1},\omega)$

such that for all initial $x_{0}\in V$ and $m_{0}\in S^{1}$ the sequence

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))$

converges to $f\circ\phi^{k}(m_{0})$ as $k\to\infty$. This $f$ is a GS. Now
suppose we choose an initial reservoir state $x_{0}$ such that

$\displaystyle\rho_{0}:=\lVert x_{0}\rVert>2.$

Then, for any initial point $m_{0}$, the reservoir states $x_{k}$ that satisfy

$\displaystyle x_{k+1}=F(x_{k},\omega\circ\phi^{k}(m_{0}))$

grow without bound. Thus any GS with image intersecting

$\displaystyle V=\\{(\rho,\theta)\ |\ \rho>2\\}$

must be unbounded. But this is impossible, as the GS is continuous and defined
on the compact manifold $S^{1}$. Thus we must conclude that any GS that exists
cannot have an image that intersects $V$. In other words, a sequence of states
originating in $V$ will not converge to the image of a GS.

###  3.7 Relationship between GS and the Echo Index

We were inspired to generalise the results surrounding global GSs to local GSs
after reading a paper by CENI2020132609. The authors present _uniformly
attracting entire solution_ (UAES) to reservoir maps and introduce the _echo
index_ , which is closely related to local GSs. We will present just enough
theory to link the GS to the Echo Index and recommend the interested reader
consult CENI2020132609 for further details.

CENI2020132609 consider reservoir maps

$\displaystyle F:X\times U\to X$ (3.4)

that satisfy 3 assumptions

  1. 1.

$F$ is continuously differentiable i.e. $F\in C^{1}(X\times U,X)$.

  2. 2.

For all $z\in\mathbb{R}^{d}$ the map $F(\cdot,z):X\to X$ is a local
diffeomorphism onto its image.

  3. 3.

$U\subset\mathbb{R}^{d}$ is compact and $X\subset\mathbb{R}^{N}$ is usually
the compact closure of a $N$-dimensional Cartesian product of real intervals.

The authors introduce the shift map $T:U^{\mathbb{Z}}\to U^{\mathbb{Z}}$
defined $T(z)_{k}=z_{k+1}$ for $k\in\mathbb{Z}$, the canonical projection
$\pi:U^{\mathbb{Z}}\to U$ defined by $\pi(z)=z_{0}$, and the cocycle mapping
which we define below.

######  Definition 3.7.1.

(Cocycle mapping) The reservoir map (3.4) can be described using a cocycle
mapping $\Phi:\mathbb{Z}_{0}^{+}\times U^{\mathbb{Z}}\times X\to X$ defined by

$\displaystyle\Phi(0,z,x_{0})$ $\displaystyle:=x_{0}$ $\displaystyle\forall
z\in U^{\mathbb{Z}},\ x_{0}\in X,$ $\displaystyle\Phi(n,z,x_{0})$
$\displaystyle:=F(\Phi(n-1,z,x_{0}),\pi\circ T^{n}(z))$ $\displaystyle\forall
z\in U^{\mathbb{Z}},\ x_{0}\in X,n\geq 1.$

######  Definition 3.7.2.

(Entire Solutions) An entire solution for the reservoir map (3.4) with input
$z\in\mathbb{Z}$ is a bi-infinite sequence of states $x\in X^{\mathbb{Z}}$
that satisfies (3.4). In other words

$\displaystyle\Phi(s,T^{m}(z),x_{m})=x_{m+s}$ (3.5)

for all $m\in\mathbb{Z}$ and $s\in\mathbb{Z}_{0}^{+}$

######  Definition 3.7.3.

(Positively Invariant Nonautonomous Sets) A family of nonempty compact subsets
$B$ is called a positively invariant nonautonomous set for input $z\in
U^{\mathbb{Z}}$ (or simply a $z$-positively invariant set) if

$\displaystyle\Phi(s,T^{m}(z),B_{m})\subset B_{s+m}$

for all $m\in\mathbb{Z}$ and $s\in\mathbb{Z}_{0}^{+}$.

######  Definition 3.7.4.

(Uniformly Attracting Entire Solution) Consider a fixed input sequence $z\in
U^{\mathbb{Z}}$, and entire solution $x$, and a $z$-positively invariant
nonautonomous set $B$ composed of compact sets.

  1. 1.

If

$\displaystyle\lim_{k\to\infty}\bigg{(}\sup_{j\in\mathbb{Z}}h(\Phi(k,T^{j}(z),B_{j}),x_{j+k})\bigg{)}=0$

(where $h$ is the Hausdorff semi-distance) then we say $B$ is uniformly
attracted to $x$.

  2. 2.

We say that $x$ is a UAES if there is a neighbourhood $B$ of $x$ that is
uniformly attracted to $x$.

######  Definition 3.7.5.

(Decompositions) We say that the reservoir map (3.4) with input $z\in
U^{\mathbb{Z}}$ admits a decomposition into $n\geq 1$ UAESs if there are $n$
UAESs $x^{1},\ldots,x^{n}$ such that for all $\eta>0$ and $i=1,\ldots,n$ there
are neighbourhoods $B^{\eta}_{i}$ uniformly attracted by $x$ and

$\displaystyle\mu\bigg{(}\mathbb{R}^{N}-\bigcup_{i=1}^{n}(B_{i}^{\eta})_{k}\bigg{)}<\eta\qquad\forall
k\in\mathbb{Z}$

where $\mu$ is the Lebesgue measure. We say this is a _proper decomposition_
if in addition

$\displaystyle\inf_{k\in\mathbb{Z}}\lVert
x^{i}_{k}-x^{j}_{k}\rVert>0\qquad\forall i\neq j.$

######  Definition 3.7.6.

(Echo Index) We say that the reservoir map (3.4) driven by input $z\in
U^{\mathbb{Z}}$ has echo index $n\geq 1$ and write

$\displaystyle\mathcal{I}(z)=n$

if it admits a proper decomposition into $n$ UAESs. In this case we say (3.4)
has the $n$-ESP for input $z$.

We can see that the work in CENI2020132609 does not reference an underlying
dynamical system $\phi:M\to M$ and therefore holds for a larger class of input
sequences than those generated by observations of $\phi$. GSs on the other
hand are defined explicitly in terms of the dynamical system $\phi$ and
observation function $\omega$, so we will need to introduce these in order to
connect GSs to the Echo Index.

######  Theorem 3.7.7.

Let $M$ be a compact topological space and $\phi\in\text{Hom}(M)$. Let
$\omega\in C^{0}(M,\mathbb{R}^{d})$ be the observation function on $M$.
Suppose that (3.4) admits a GS $f:M\to V\subset X\subset\mathbb{R}^{N}$. Then
for each $m_{0}\in M$, the sequence
$(x_{k})_{k\in\mathbb{Z}}=(f\circ\phi^{k}(m_{0}))_{k\in\mathbb{Z}}$ is a UAES
attracted by the constant sequence $B_{k}=V$ for all $k$.

###### Proof.

For any $m_{0}\in M$, let
$(z_{k})_{k\in\mathbb{Z}}=(\omega\circ\phi^{k}(m))_{k\in\mathbb{Z}}$. The
constant sequence $B_{k}=V$ for all $k$ is uniformly attracted to
$(x_{k})_{k\in\mathbb{Z}}=(f_{i}^{\eta}\circ\phi^{k}(m))_{k\in\mathbb{Z}}$ by
Theorem 3.5.4. ∎

######  Theorem 3.7.8.

Let $M$ be a compact topological space and $\phi\in\text{Hom}(M)$. Let
$\omega\in C^{0}(M,\mathbb{R}^{d})$ be the observation function on $M$.
Suppose there are exactly $n$ GSs $f_{i}:M\to f_{i}(M)\subset\mathbb{R}^{N}$,
such that, for all $\eta>0$, there exist $n$ sets $V^{\eta}_{i}\supset
f_{i}(M)$ such that $f_{i}^{\eta}:M\to V^{\eta}_{i}$ is a GS and

$\displaystyle\mu\bigg{(}\mathbb{R}^{N}-\bigcup_{i=1}^{n}V^{\eta}_{i}\bigg{)}<\eta,$
(3.6)

where $\mu$ is the Lebesgue measure. Then for any $m_{0}\in M$, the sequence
$(z_{k})_{k\in\mathbb{Z}}=(\omega\circ\phi^{k}(m_{0}))_{k\in\mathbb{Z}}$ has
echo index $\mathcal{I}[z]=n$.

###### Proof.

Fix $m_{0}\in M$. For each GS $f_{i}^{\eta}:M\to V^{\eta}_{i}$ the sequence
$(x^{i}_{k})_{k\in\mathbb{Z}}=(f^{i}_{\eta}\circ\phi^{k}(m_{0}))_{k\in\mathbb{Z}}$
is a UAES attracted by the constant sequence $(B_{i}^{\eta})_{k}=V_{i}^{\eta}$
for all $k$. This follows from Theorem 3.7.7. Now, with respect to the input
sequence
$(z_{k})_{k\in\mathbb{Z}}=(\omega\circ\phi^{k}(m_{0}))_{k\in\mathbb{Z}}$, the
reservoir map clearly admits a _decomposition_ into $n$ UAESs, so all that
remains is to show that this decomposition is _proper_. The images $f_{i}(M)$
are disjoint compact manifolds and each sequence
$(x^{i}_{k})_{k\in\mathbb{Z}}=(f_{i}\circ\phi^{k}(m_{0}))_{k\in\mathbb{Z}}$
lies in exactly one of these images. Consequently

$\displaystyle\inf_{k\in\mathbb{Z}}\lVert
x^{i}_{k}-x^{j}_{k}\rVert>0\qquad\forall\,i\neq j.$

which is precisely the requirement for the decomposition to be _proper_ in the
sense of definition 3.7.5. ∎

###  3.8 Biological interpretation

These results about GSs admit a biological interpretation. We use the ESN
reservoir sequence

$\displaystyle x_{k+1}=\sigma(Ax_{k}+Cz_{k}+b)$ (3.7)

to represent the dynamics of a brain, composed of neurons, that are connected
together by dendrites and axons. Neuron $n$ has dendrites, which receive
information from neighbouring neurons, causing neuron $n$ to fire a sequence
of action potentials, thus influencing the other neurons
(izhikevich2007dynamical). We (over)simplify this complex interaction by
representing the neurons with nodes, and the connections between them with
arcs. We represent the connection strength between neurons $i$ and $j$ with
the entry $A_{ij}$ of the reservoir matrix $A$. If neurons $i$ and $j$ are not
directly connected then $A_{ij}=0$. Furthermore, the connections need not be
symmetric, so $A_{ij}\neq A_{ji}$ in general.

The $i$th neuron has an _excitation level_ at time $k$ given by the $i$th
component of the reservoir state $x_{k}$. The excitation level of neuron $i$
governs the frequency of action potentials fired by neuron $i$. It therefore
follows from (3.7) that the excitation of each neuron depends on the
excitation of its neighbours. If $\sigma=\tanh$ or some other
saturating/squashing/sigmoidal function, then the contribution of a neuron’s
neighbours obeys a law of diminishing returns.

We view the inputs $z_{k}$ as electrical signals induced in the brain from an
external source. These signals could arise in the body or from external
sensory organs. The inputs interact with the neurons via connections $C$, and
influence the excitation levels of the neurons in the brain. The results we
have presented in this chapter suggest that, as long as the source of the
inputs $z_{k}$ evolves deterministically (and the ESN has ESP) then the state
of the brain (represented by the state vectors $x_{k}$) evolves in synchrony
with the source system. We illustrate this interplay between the input
connections and reservoir connections in Figure 3-8 which illustrates an ESN
with 3 neurons.

Doctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor
of PhilosophyDoctor of PhilosophyDoctor of PhilosophyDoctor of
Philosophy$x_{1}$$x_{1}$Doctor of Philosophy$x_{2}$$x_{2}$Doctor of
Philosophy$x_{3}$$x_{3}$Doctor of PhilosophyDoctor of PhilosophyDoctor of
PhilosophyDoctor of PhilosophyDoctor of Philosophy$z$Doctor of
Philosophy$A_{12}$Doctor of Philosophy$A_{31}$Doctor of
Philosophy$A_{23}$Doctor of Philosophy$A_{21}$Doctor of
Philosophy$A_{13}$Doctor of Philosophy$A_{32}$Doctor of
Philosophy$A_{11}$Doctor of Philosophy$A_{22}$Doctor of
Philosophy$A_{33}$Doctor of Philosophy$C_{1}$Doctor of Philosophy$C_{2}$Doctor
of Philosophy$C_{3}$ Figure 3-8: An ESN driven by input $z$. The input
communicates to three nodes with states $x_{1},x_{2},$ and $x_{3}$ via
connections $C_{1},C_{2},$ and $C_{3}$. The nodes communicate with eachother
via the three by three connectivity matrix $A$ with connection between node
$i$ and $j$ given by $A_{ij}$. The direction of the arrows denotes the
direction of communication.

Interestingly, if the ESN does not have ESP then the vectors $x_{k}$ may
evolve erratically, and fail to synchronise with the source system. Some
authors (manjunath2017evolving) have suggested that this loss of the ESP
describes an epileptic seizure, triggered by an input process $z_{k}$ that is
outside a tolerated range (e.g rapidly flashing lights). We might imagine that
a brain with $(V,W)$-local ESP undergoes an epileptic seizure if exposed to
inputs outside of $W$. The existence of multiple GSs for different sets of
inputs $W$ and sets of initial reservoir states $V$ suggests that the brain
may adopt qualitatively different states under different conditions, that are
each locally stable. This is all very speculative, but in future, it would be
intriguing to explore these ideas in the context of careful neuroscience.

##  Chapter 4 Embeddings

###  4.1 Background

Theorems 3.5.4 and 3.5.5 establish conditions under which an evolution
operator $\phi\in\text{Diff}^{1}(M)$, observation function
$\omega:M\to\mathbb{R}^{d}$ and reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ admit an associated GS
$f_{(\phi,\omega,F)}$. Theorem 3.5.5 establishes conditions under which the GS
$f_{(\phi,\omega,F)}$ is continuously differentiable. In this chapter we will
establish additional conditions, which ensure that $f_{(\phi,\omega,F)}$ is a
continuously differentiable _embedding_. When $f_{(\phi,\omega,F)}$ is an
embedding, the dynamics in the image of $f_{(\phi,\omega,F)}$ (the response
system) replicate the dynamics of the drive system $(M,\phi)$ so that we have
a chance of learning arbitrary target functions $g:M\to\mathbb{R}^{s}$ from
observations alone. Out of the collection of all possible target maps, a
target map in which we are particularly interested is the next step map
$g:=\omega\circ\phi$ which is well defined when $s=d$. For the rest of this
chapter we assume for simplicity, and without loss of generality, that $d=1$
so we have scalar observations. This chapter is based on Grigoryeva et al.
(2021b), which builds on the earlier work in
embedding_and_approximation_theorems.

We will now define what it means for a map to be an embedding, and the closely
related property of being an immersion.

######  Definition 4.1.1.

(Immersion) A map $f:M\to N$ between smooth manifolds $M,N$ is called an
immersion at the point $m\in M$ if the tangent map at $m$ denoted
$T_{m}f:T_{m}M\to T_{f(m)}N$ is injective. If $f$ is an immersion at all
points $m\in M$ then $f$ is called an immersion.

######  Definition 4.1.2.

(Embedding) A map $f:M\to N$ between smooth manifolds $M,N$ is called an
embedding if $f$ is an injective immersion and a diffeomorphism onto its
image.

Notably, if $M$ is compact and $f:M\to N$ is an injective immersion then $f$
is an embedding. Hence we restrict ourselves to compact manifolds and consider
the conditions on $\phi,\omega,F$ ensure that $f_{(\phi,\omega,F)}$ is an
injective immersion. This problem is difficult in general, so we will start by
recalling some fundamental theorems. The first is Whitney’s embedding theorem.

######  Theorem 4.1.3.

(whitney_embedding_thm). Let $M$ be a compact $q$-manifold (i.e a manifold of
dimension $q$) and $N\in\mathbb{N}$ such that $N>2q$. Then a generic $f\in
C^{1}(M,\mathbb{R}^{N})$ is an embedding.

To make sense of this theorem, we need to understand the term generic. Roughly
speaking, a property is generic in some topological space if it holds for most
elements, but excludes a small subset of exceptions. The concept of genericity
is a topological one, which is analogous to (but not identical to) the
concepts of _almost everywhere_ and _prevalence_ that appear in measure
theory.

Ghrist2014 provides some examples of genericity. In the appropriate topology:

  * •

a generic smooth curve $\gamma:\mathbb{R}\to\mathbb{R}^{2}$ self-intersects
countably many times,

  * •

a generic smooth curve $\gamma:\mathbb{R}\to\mathbb{R}^{3}$ does not self-
intersect,

  * •

a generic square matrix is invertible.

We state the formal definition of genericity below.

######  Definition 4.1.4.

(Generic property) Let $(X,\tau)$ be a topological space. A property that
holds on a countable intersection of dense open subsets of $X$ is called a
generic property of $(X,\tau)$.

Notably if $X$ is compact, then a property that holds on a dense open subset
of $X$ is generic. We can now see that Whitney’s embedding theorem implies
that if we set $N>2q$, where $M$ is a $q$ dimensional manifold, and we impose
conditions on $\phi,\omega,F$ such that $f_{(\phi,\omega,F)}\in
C^{1}(M,\mathbb{R}^{N})$, then $f_{(\phi,\omega,F)}$ is either an embedding,
or very close to an embedding. Heuristically, we expect that
$f_{(\phi,\omega,F)}$ is very likely to be an embedding unless $\phi,\omega,F$
are deviously (or unfortunately) chosen such that the resulting GS
$f_{(\phi,\omega,F)}$ is an exceptional map in $C^{1}(M,\mathbb{R}^{N})$ that
is not an embedding. This heuristic is quite compelling, but is not a proof
that $f$ is an embedding for generic choices $\phi,\omega,F$. To progress
toward this result, we will recall Takens’ celebrated embedding theorem. The
specific formulation and proof of the result is by Hukes_thm, and it is this
formulation which we will use in great detail in this Chapter.

######  Theorem 4.1.5.

(TakensThm), (Hukes_thm). Let $M$ be a compact manifold of dimension $q$.
Suppose $\phi\in\text{Diff}^{2}(M)$ has the following two properties:

  * (1)

$\phi$ has only finitely many periodic points with periods less than or equal
to $2q$.

  * (2)

For each period point $m\in M$ with period $n<2q$ then the eigenvalues of the
derivative $T_{m}\phi^{n}$ at $m$ are distinct. $\phi^{n}$ denotes the
$n$-fold composition of $\phi$ with itself.

Then for a generic $C^{2}$ observation function $\omega\in
C^{2}(M,\mathbb{R})$ the $(2q+1)$ delay observation map
$\Phi_{(\phi,\omega)}:M\to\mathbb{R}^{2q+1}$ defined by

$\displaystyle\Phi_{(\phi,\omega)}(m):=(\omega(m),\omega\circ\phi(m),\omega\circ\phi^{2}(m),\ldots,\omega\circ\phi^{2q}(m))$

is an embedding.

Takens Theorem essentially states that: if we take equally spaced (in time)
scalar observations from a dynamical system of dimension $q$, and the system
$(M,\phi)$ satisfies conditions $(1)$ and $(2)$, then a vector formed from
$2q+1$ (or more) sequential observations fully captures the state of the
dynamical system.

###  4.2 Summary of novel results

The connection between reservoir computing and Takens’ Theorem becomes clear
when we observe that the reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}\to\mathbb{R}^{N}$ defined by

$\displaystyle F(x,z)=Ax+Cz$

admits as a GS the Takens delay map

$\displaystyle f_{(\phi,\omega,F)}=\Phi_{(\phi,\omega)}$

in the special case that $A\in\mathbb{M}_{N\times N}(\mathbb{R})$ is the lower
shift matrix

$\displaystyle A=\begin{bmatrix}0&0&0&\ldots&0\\\ 1&0&0&\ldots&0\\\
0&1&0&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&1&0\end{bmatrix}$

and $C\in\mathbb{R}^{N}$ is defined to be $C=(1,0,\ldots,0)$. This has been
observed by 4118282, embedding_and_approximation_theorems,
doi:10.1063/5.0024890, Grigoryeva et al. (2021b) and others. It therefore
seems natural to consider more general conditions on $A$ and $C$ for which the
associated GS $f_{(\phi,\omega,F)}$ is an embedding for generic observation
functions $\omega\in C^{2}(M,\mathbb{R})$. Such conditions are established in
the next result.

######  Theorem 4.2.1.

Let $M$ be a compact manifold of dimension $q$. Suppose
$\phi\in\text{Diff}^{2}(M)$ has the following two properties:

  * (1)

$\phi$ has only finitely many periodic points.

  * (2)

For each periodic point $m\in M$ with period $n$ the eigenvalues of the
derivative $T_{m}\phi^{n}$ at $m$ are distinct.

Let $F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ be defined by

$\displaystyle F(x,z)=Ax+Cz$

for $A\in\mathbb{M}_{N\times N}(\mathbb{R})$ and $C\in\mathbb{R}^{N}$
satisfying the following properties:

  * (A)

$N>\max\\{2q,\ell\\}$ where $\ell\in\mathbb{N}$ is the lowest common multiple
of the periods of all periodic points.

  * (B)

$\lambda_{\text{max}}\,\rho(A^{n_{\text{min}}})<1$ where $n_{\text{min}}$ is
the minimal period over all periodic points and $\lambda_{\text{max}}$ is the
maximal absolute value over all eigenvalues of all derivatives
$T_{m}\phi^{n}$.

  * (C)

For each periodic point $m\in M$ with period $n$ the set of vectors

$\displaystyle\bigg{\\{}(I-\lambda_{j}A^{n})^{-1}(I-A)^{-1}(I-A)^{n}C\bigg{\\}}_{j=1,\ldots,q}$

where $\\{\lambda_{j}\\}_{j=1,\ldots q}$ are the eigenvalues of
$T_{m}\phi^{n}$, are linearly independent.

  * (D)

The vectors $\\{A^{j}C\\}_{j=0,\ldots N-1}$ are linearly independent.

Then there exists an associated GS $f_{(\phi,\omega,F)}\in
C^{2}(M,\mathbb{R}^{N})$ with explicit form

$\displaystyle
f_{(\phi,\omega,F)}(m)=\sum_{k=0}^{\infty}A^{k}C\omega\circ\phi^{-k}(m).$

Furthermore, for generic $\omega\in C^{2}(M,\mathbb{R})$ the GS
$f_{(\phi,\omega,F)}$ is an embedding.

######  Remark 4.2.2.

Condition $(1)$ appearing in Theorem 4.2.1 is stronger than condition $(1)$
appearing in Takens Theorem, so it may be possible to relax the former
condition.

######  Remark 4.2.3.

Condition $(D)$ appears to have interesting connections to controllability
(Katsuhiko, 2010) and Krylov spaces (Nevanlinna, 1993).

The proof proceeds in three major stages, and closely follows the proof of
Takens Theorem presented by Hukes_thm. We prove first of all that

  1. 1.

The associated GS $f_{(\phi,\omega,F)}\in C^{2}(M,\mathbb{R}^{N})$ exists and
has explicit form

$\displaystyle
f_{(\phi,\omega,F)}(m)=\sum_{k=0}^{\infty}A^{k}C\omega\circ\phi^{-k}(m).$

  2. 2.

Next, we show that the set of $\omega\in C^{2}(M,\mathbb{R})$ for which the
associated GS $f_{(\phi,\omega,F)}$ is an embedding is an open subset of
$C^{2}(M,\mathbb{R})$.

  3. 3.

Then all that remains is to show that a dense subset of the observation
functions $C^{2}(M,\mathbb{R})$ yields an embedding. This process is split
into 4 substages.

     * •

We show that for arbitrary $\omega\in C^{2}(M,\mathbb{R})$ we can take an
arbitrarily small perturbation $\omega^{\prime}\in C^{2}(M,\mathbb{R})$ such
that the perturbed GS $f_{(\phi,\omega^{\prime},F)}$ is an immersion on the
periodic points. It immediately follows that $f_{(\phi,\omega^{\prime},F)}$ is
also an immersion on some neighbourhood of the periodic points by the
immersion theorem (Theorem 4.4.6).

     * •

Immersions form an open subset of $C^{1}(M,\mathbb{R}^{N})$, so we can create
a second perturbation $\omega^{\prime\prime}\in C^{2}(M,\mathbb{R})$
sufficently small that the open neighbourhood of the periodic points remains
immersed, while immersing all remaining points. Then we have that
$f_{(\phi,\omega^{\prime\prime},F)}$ is a global immersion. All that remains
is to establish injectivity.

     * •

We take another perturbation $\omega^{\prime\prime\prime}$ sufficiently small
that the immersion is preserved while ensuring that
$f_{(\phi,\omega^{\prime\prime\prime},F)}$ is injective on the periodic
points. By continuity of $f_{(\phi,\omega^{\prime\prime\prime},F)}$, it
follows that there is a neighbourhood of the periodic points on which
$f_{(\phi,\omega^{\prime\prime\prime},F)}$ is injective.

     * •

We take one final perturbation $\omega^{\prime\prime\prime\prime}$
sufficiently small that the immersion is preserved and
$f_{(\phi,\omega^{\prime\prime\prime\prime},F)}$ is globally injective. This
completes the proof.

The full proof is detailed in Grigoryeva et al. (2021b) and section 4.4 of
this thesis. It is also shown in Grigoryeva et al. (2021b) and in section 6.2
of this thesis that conditions $(A)-(D)$ are invariant under _system