# Generalised Synchronisations, Embeddings, and Approximations for Continuous
Time Reservoir Computers

Allen G Hart

###### Abstract

We establish conditions under which a continuous time reservoir computer, such
as a leaky integrator echo state network, admits a generalised synchronisation
$f$ between between the source dynamics and reservoir dynamics. We show that
multiple generalised synchronisations can exist simultaneously, and connect
this to the multi-Echo-State-Property (multi-ESP). In the special case of a
linear reservoir computer, we derive a closed form expression for the
generalised synchronisation $f$. Furthermore, we establish conditions under
which $f$ is of class $C^{1}$, and conditions under which $f$ is a topological
embedding on the fixed points of the source system. This embedding result is
closely related to Takens’ embedding Theorem.

We also prove that the embedding of fixed points occurs almost surely for
randomly generated linear reservoir systems. With an embedding achieved, we
discuss how the universal approximation theorem makes it possible to forecast
the future dynamics of the source system and replicate its topological
properties. We illustrate the theory by embedding a fixed point of the
Lorenz-63 system into the reservoir space using numerical methods. Finally, we
show that if the observations are perturbed by white noise, the GS is
preserved up to a perturbation by an Ornstein-Uhlenbeck process.

###### keywords:

Generalised Synchronisation , Reservoir Computing

###### PACS:

0000 , 1111

###### MSC:

0000 , 1111

††journal: Physica D: Nonlinear Phenomena

[inst1]organization=University of Bath,addressline=Claverton Down, city=Bath,
postcode=BA2 7AY, country=UK

##  1 Introduction

###  1.1 Reservoir Computing

After the seminal papers of [18] and [26] at the turn of the millennium,
reservoir computers have grown in popularity, and are now widely studied in
mathematics [12, 13, 8, 4] physics [17, 38], computer science and robotics
[31, 38]. Authors often use a reservoir map
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ to analyse a discrete
time series $\\{z_{k}\in\mathbb{R}^{d}\\}_{k\in\mathbb{Z}}$ by creating
reservoir states $\\{x_{k}\in\mathbb{R}^{N}\\}_{k\in\mathbb{Z}}$ via the
iteration

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

for some initial state $x_{0}\in\mathbb{R}^{N}$. Having created the reservoir
states $\\{x_{k}\in\mathbb{R}^{N}\\}_{k\in\mathbb{Z}}$ the practitioner may
train the reservoir system to approximate a series of targets
$\\{u_{k}\in\mathbb{R}^{s}\\}_{k\in\mathbb{Z}}$ by optimising a set of
parameters $W$ such that a map (often a neural network) parametrised by $W$
approximately maps the reservoir states $x_{k}$ to the targets $u_{k}$.

A popular choice of reservoir map is the Echo State Network (ESN) [18, 19] of
the form

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

where

  * 1.

$\sigma:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is an activation function

  * 2.

$A$ is a random real square $N\times N$ matrix, called the reservoir matrix

  * 3.

$C$ is a random real $N\times d$ matrix, called the input matrix

  * 4.

$b$ is a random real $N$-vector, called the bias.

The targets $u_{k}$ are often approximated by minimising over the $N\times s$
matrices $W$ the Tikhonov regularised least squares

$\displaystyle L(W)=\sum_{k=1}^{\ell}\lVert
W^{\top}x_{k}-u_{k}\rVert^{2}+\lambda\lVert W\rVert^{2}$

for $L$ the loss function, $\lambda>0$ the regularisation parameter, and
$\ell$ the finite number of training points.

ESNs possess the universal approximation property [12, 9, 14], which allows
them to approximate arbitrary relationships between the reservoir states
$x_{k}$ and targets $u_{k}$. ESNs have proved themselves competitive in
forecasting chaotic time series [20] where the targets $u_{k}=z_{k}$ are the
observations. A close relative of the ESN is the leaky integrator ESN [21, 25,
24]

$\displaystyle F(x,z):=-\alpha x+\sigma(Ax+Cz+b)$

where $\alpha>0$ represents the information ‘leak’. The leaky integrator ESN
is usually applied to continuous time series $z:\mathbb{R}\to\mathbb{R}^{d}$
instead of discrete series $\\{z_{k}\\}_{k\in\mathbb{Z}}$. When dealing with a
continuous time series, we generate a continuous stream of reservoir states by
integrating the non-autonomous ODE

$\displaystyle\dot{x}(t)=F(x(t),z(t))$ (1)

from initial condition $x_{0}\in\mathbb{R}^{N}$. The continuous time case is
especially interesting when modelling physical reservoir computers; which are
implemented on physical systems other than ordinary computers. Such systems
include photonic node arrays [38] and _origami_ [2] and have been studied in
the context of robotic locomotion [2]. The physics of these exotic reservoir
systems are generally understood in terms of ODEs, in contrast to discrete
time maps which are usually preferred when the reservoir computer is
implemented digitally on an ordinary computer. These continuous time reservoir
systems described by ODEs are the central object of this paper, which focuses
on the particular setting where $z(t)$ is given by taking scalar (or low
dimensional) observations of a higher dimensional source system.

To be more specific, we imagine that there is vector field $\mathcal{V}$ on a
manifold $M$, called the source system, which is hidden from view. We are able
only to observe a trajectory of the source system via a scalar observation
function $\omega:M\to\mathbb{R}$. Our goal is to feed this continuous scalar
observation into a reservoir system such as (1) in the hope that we can
replicate the hidden source system in the reservoir space. If we are
successful, then there is a map from the source dynamics to the reservoir
space $\mathbb{R}^{N}$ called a generalised synchronisation (GS)
$f_{(\omega,\mathcal{V},F)}:M\to\mathbb{R}^{N}$. Our goal is then to learn
about the source dynamics from the image of $f_{(\omega,\mathcal{V},F)}$, and
consider whether the GS $f_{(\omega,\mathcal{V},F)}$ allows us to train the
reservoir system for tasks including forecasting future trajectories.

We formalise these ideas throughout the remainder of this paper, which is set
out as follows. In section 2 we establish conditions under which a reservoir
system admits a generalised synchronisation in the sense defined by [22] and
studied in [33, 30, 36, 3, 5]. In section 2 we connect the existence of
multiple synchronisation manifolds to the multi-ESP introduced in [4]. In
section 4 we consider the special case of a linear reservoir system

$\displaystyle\dot{x}(t)=F(x(t),z(t))=-Ax(t)+Cz(t)$

for $A$ a square $N\times N$ reservoir matrix and $C$ a rectangular $N\times
d$ input matrix. In this case, we derive a closed form expression for the
associated generalised synchronisation $f_{(\omega,\mathcal{V},F)}$, establish
conditions under which $f_{(\omega,\mathcal{V},F)}\in C^{1}$ is continuously
differentiable, and show that for randomly generated $A,C$ the generalised
synchronisation $f_{(\omega,\mathcal{V},F)}$ is a topological embedding on the
fixed points of the source dynamics. These results are continuous time
analogues of very recent results [11] that hold for discrete time reservoir
systems.

Our approach is similar in many respects to the approach taken in [11] except
we frame our results in terms of ODEs while the authors of [11] frame theirs
in terms of discrete time maps. We believe proving results in a continuous
time setting may further the understanding of popular continuous time
reservoir systems, like leaky integrator ESNs, which are often used to model
physical reservoir systems.

##  2 GS for Nonlinear Reservoir Systems

We will begin with the definition of a $V$-invariant ODE. Roughly speaking, a
solution to the ODE which originates in the set $V$ will stay in $V$ for all
future time.

######  Definition 2.1.

($V$-invariant) Let $z\in C^{0}(\mathbb{R},\mathbb{R}^{d})$ be a bounded
function and $F\in C^{0}(\mathbb{R}^{N}\times\mathbb{R}^{d},\mathbb{R}^{N})$
be Lipschitz continuous. Then the non-autonomous ODE

$\displaystyle\dot{x}(t)=F(x(t),z(t))$

is called $V\subset\mathbb{R}^{N}$ invariant with respect to $z$ if for any
initial point $x_{0}\in V$ the solution $x:\mathbb{R}\to\mathbb{R}^{N}$
originating at $x_{0}\in V$ remains in $V$ i.e $x(t)\in V$ for all
$t\in\mathbb{R}^{+}$.

Next, we define the concept of uniform $V$-asymptotic stability. Roughly we
call a non-autonomous ODE uniformly $V$-asymptotically stable if, given any
two initial points in $V$ the solutions originating from those points converge
toward each other at a rate that is uniform over driving inputs $z$.

######  Definition 2.2.

(Uniform $V$-asymptotic stability) Let $M$ be a set, and for each $m\in M$ let
$z_{m}\in C^{0}(\mathbb{R},\mathbb{R}^{d})$. Then the non-autonomous ODEs

$\displaystyle\dot{x}(t)=F(x(t),z_{m}(t))$ (2)

are uniformly $V$-asymptotically stable with respect to $\\{z_{m}\\}_{m\in M}$
if there exists a positive valued function $\psi:\mathbb{R}\to\mathbb{R}^{+}$
converging to $0$ such that for all $m\in M$ and initial points
$x_{0},y_{0}\in V\subset\mathbb{R}^{N}$ the solutions $x_{m}(t)$ and
$y_{m}(t)$ to ODE (2) originating from $x_{0},y_{0}$ satisfy

$\displaystyle\lVert x_{m}(t)-y_{m}(t)\rVert\leq\psi(t)$

for all $t\geq 0$.

######  Remark 2.3.

The notion of uniform $V$-asymptotic stability is closely related to the Echo
State Property given in Definition 1 of [21], and global input related
stability discussed in [27].

We are especially interested in ODEs given by a reservoir map $F$ driven by
observations from source dynamics evolving on a manifold $M$. This system is
called a reservoir system.

######  Definition 2.4.

(Reservoir system) Suppose that $\mathcal{V}$ is a smooth vector field on a
smooth manifold $M$ with associated evolution operators
$\\{\phi^{t}\in\text{Diff}^{1}(M)\ |\ t\in\mathbb{R}\\}$ that form a group
under composition such that $\phi^{t_{1}+t_{2}}=\phi^{t_{1}}\phi^{t_{2}}$. Let
$\omega\in C^{0}(M,\mathbb{R}^{d})$ be a continuous observation function, and
let $F\in C^{1}(\mathbb{R}^{N}\times\mathbb{R}^{d},\mathbb{R}^{N})$ be
Lipschitz continuous. Then we can define the continuous time reservoir system

$\displaystyle\dot{x}(t)=F(x(t),\omega\phi^{t}(m)).$ (3)

Associated to the manifold $M$ and vector field $\mathcal{V}$ of the source
system is the Lie derivative $\mathcal{L}_{\mathcal{V}}$. We introduce the Lie
derivative here because it will arise later in an expression involving the
generalised synchronisation $f_{(\omega,\mathcal{V},F)}$.

######  Definition 2.5.

(Lie derivative) Suppose that $\mathcal{V}$ is a smooth vector field on a
smooth manifold $M$ with associated evolution operators
$\\{\phi^{t}\in\text{Diff}^{1}(M)\ |\ t\in\mathbb{R}\\}$. Then the Lie
derivative of a map $u:M\to\mathbb{R}^{N}$ is defined by

$\displaystyle\mathcal{L}_{\mathcal{V}}u(m)=\frac{d}{dt}u\phi^{t}(m)\rvert_{t=0}$

when the derivative on the right hand side exists.

We are now ready to introduce the concept of $V$-generalised synchronisation
studied by [22].

######  Definition 2.6.

($V$-Generalised Synchronisation, [22]) Reservoir system (3) admits a
$V$-generalised synchronisation
$f_{(\omega,\mathcal{V},F)}:M\to\mathbb{R}^{N}$ if, for any initial $r\in
V\subset\mathbb{R}^{N}$ and $m\in M$, the solution $x^{r}_{m}(t)$ of system
(3) satisfies

$\displaystyle\lim_{t\to\infty}\lVert
f_{(\omega,\mathcal{V},F)}\phi^{t}(m)-x^{r}_{m}(t)\rVert=0.$

We will introduce a slightly stronger concept of synchronisation which we call
uniform $V$-generalised synchronisation, which ensures that the map
$f_{(\omega,\mathcal{V},F)}$ is continuous.

######  Definition 2.7.

(Uniform $V$-Generalised Synchronisation) Reservoir system (3) admits a
uniform $V$-generalised synchronisation $f_{(\omega,\mathcal{V},F)}\in
C^{0}(M,V)$ if there exists a positive function
$\psi:\mathbb{R}\to\mathbb{R}^{+}$ converging to $0$ such that for any initial
$r\in V\subset\mathbb{R}^{N}$ and $m\in M$ the solution $x^{r}_{m}(t)$ of
system (3) satisfies

$\displaystyle\lVert
f_{(\omega,\mathcal{V},F)}\phi^{t}(m)-x^{r}_{m}(t)\rVert\leq\psi(t)$

for all $t\geq 0$.

The seminal paper by [22] proves that asymptotic stability is equivalent to
the existence of a generalised synchronisation. We will prove a very similar
result - stating that uniform $V$-asymptotic stability is equivalent to the
existence of a uniform $V$-generalised synchronisation. The uniform
$V$-generalised synchronisation $f_{(\omega,\mathcal{V},F)}$ is continuous,
and the image $f_{(\omega,\mathcal{V},F)}(M)$ is therefore a manifold. These
two properties of $f_{(\omega,\mathcal{V},F)}$ are convenient and not
established explicitly in [22].

######  Theorem 2.8.

[22] Reservoir system (3) admits a uniform $V$-generalised synchronisation
$f_{(\omega,\mathcal{V},F)}\in C^{0}(M,\mathbb{R}^{N})$ if and only if (3) is
uniformly $V$-asymptotically stable with respect to
$\\{\omega\phi^{t}(m)\\}_{m\in M}$. Furthermore, the uniform $V$-generalised
synchronisation $f_{(\omega,\mathcal{V},F)}$ solves the quasilinear PDE

$\displaystyle\mathcal{L}_{\mathcal{V}}f_{(\omega,\mathcal{V},F)}=F(f_{(\omega,\mathcal{V},F)},\omega).$

###### Proof.

Suppose first of all that (3) is uniformly $V$ asymptotically stable with
respect to $\\{\omega\phi^{t}(m)\\}_{m\in M}$. Let $x^{r}_{m}(t)$ denote the
solution of ODE (3) originating from the initial point $r\in
V\subset\mathbb{R}^{N}$. Then by uniform $V$-asymptotic stability there exists
a positive function $\psi:\mathbb{R}\to\mathbb{R}^{+}$ converging to $0$ such
that, for any $s,\sigma>0$ that satisfy $0\leq s\leq\sigma$, the solutions
$x^{r}_{\phi^{-s}(m)}(t)$ and
$x^{x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)}_{\phi^{-s}(m)}(t)$ originating from
$r\in V$ and $x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)\in V$ satisfy

$\displaystyle\lVert
x^{r}_{\phi^{-s}(m)}(t)-x^{x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)}_{\phi^{-s}(m)}(t)\rVert\leq\psi(t)$

for all $t\geq 0$. Then at the particular time $t=s$

$\displaystyle\lVert
x^{r}_{\phi^{-s}(m)}(s)-x^{x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)}_{\phi^{-s}(m)}(s)\rVert\leq\psi(s).$

Then by Lemma 8.1

$\displaystyle
x^{r}_{\phi^{-\sigma}(m)}(\sigma)=x^{x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)}_{\phi^{-s}(m)}(s)$

so

$\displaystyle\lVert
x^{r}_{\phi^{-s}(m)}(s)-x^{r}_{\phi^{-\sigma}(m)}(\sigma)\rVert\leq\psi(s).$

Now $\psi$ is a positive function independent of $m\in M$ that converges to
$0$ so the convergence of

$\displaystyle\lim_{s\to\infty}x^{r}_{\phi^{-s}(m)}(s)=:f^{r}(m)$

is uniform over $m\in M$. The uniform convergence ensures $f^{r}\in
C^{0}(M,\mathbb{R}^{N})$ is continuous. Furthermore for any other initial
point $\rho\in V\subset\mathbb{R}^{N}$ it follows from $V$-asymptotic
stability that

$\displaystyle\lVert
x^{\rho}_{\phi^{-s}(m)}(s)-x^{r}_{\phi^{-s}(m)}(s)\rVert\leq\psi(s)$

hence

$\displaystyle
f^{\rho}(m)=\lim_{s\to\infty}x^{\rho}_{\phi^{-s}(m)}(s)=\lim_{s\to\infty}x^{r}_{\phi^{-s}(m)}(s)=f^{r}(m)$

so $f^{\rho}(m)=:f(m)$ does not depend on $\rho\in V$. Now

$\displaystyle
f\phi^{t}(m)=\lim_{s\to\infty}x^{r}_{\phi^{-s+t}(m)}(s)=\lim_{s\to\infty}x^{x^{r}_{\phi^{-s}(m)}(s)}_{m}(t)=x_{m}^{f(m)}(t)$

so $f\phi^{t}(m)$ is the solution to the reservoir system originating from the
initial point $f(m)\in V\subset\mathbb{R}^{N}$. Then for any initial $\rho\in
V$ and $m\in M$ the solution $x^{\rho}_{m}(t)$ satisfies

$\displaystyle\lVert x^{\rho}_{m}(t)-f\phi^{t}(m)\rVert\leq\psi(t)$

which establishes that $f$ is a uniform $V$-generalised synchronisation. Now
to prove the converse observe that for any $m\in M$ and $r,\rho\in
V\subset\mathbb{R}^{N}$

$\displaystyle\lVert x^{r}_{m}(t)-x^{\rho}_{m}(t)\rVert$
$\displaystyle\leq\lVert
x^{r}_{m}(t)-f_{(\omega,\mathcal{V},F)}\phi^{t}(m)\rVert+\lVert
f_{(\omega,\mathcal{V},F)}\phi^{t}(m)-x^{\rho}_{m}(t)\rVert$
$\displaystyle\leq 2\psi(t)$

so the existence of the uniform $V$-GS implies uniform $V$ asymptotic
stability. Now we can take the Lie derivative of $f\phi^{t}(m)$ to see that

$\displaystyle\mathcal{L}_{\mathcal{V}}f(m)$
$\displaystyle=\frac{d}{dt}f\phi^{t}(m)\rvert_{t=0}$
$\displaystyle=\frac{d}{dt}x_{m}^{f(m)}(t)\rvert_{t=0}$
$\displaystyle=\dot{x}_{m}^{f(m)}(t)\rvert_{t=0}$
$\displaystyle=F(x_{m}^{f(m)}(0),\omega(m))$ $\displaystyle=F(f(m),\omega(m))$

hence in general

$\displaystyle\mathcal{L}_{\mathcal{V}}f=F(f,\omega).$

∎

It is useful to establish conditions on the reservoir map $F$ that ensure that
the associated reservoir system (3) is uniformly $V$-asymptotically stable
with respect to $\\{\omega\phi^{t}(m)\\}_{m\in M}$, as this ensures the
existence of a uniform $V$-GS $f_{(\omega,\mathcal{V},F)}$.

######  Theorem 2.9.

Let $V\subset\mathbb{R}^{N}$ be a bounded convex set and suppose reservoir
system (3) is $V$-invariant. Suppose further there exists a $\delta>0$ such
that for any $w,v\in V$ and $z\in\omega(M)$

$\displaystyle\frac{v^{\top}D_{x}F(w,z)v}{\lVert v\rVert^{2}}<-\delta,$ (4)

where $DF_{x}(x,z)$ denotes the derivative of $F(x,z)$ with respect to $x$.
Then reservoir system 3 is uniformly $V$-asymptotically stable with respect to
$\\{\omega\phi^{t}(m)\\}_{m\in M}$ hence admits a uniform $V$-generalised
synchronisation $f_{(\omega,\mathcal{V},F)}\in C^{0}(M,V)$.

###### Proof.

$V\subset\mathbb{R}^{N}$ is bounded so there exists a $K>0$ such that for any
$x_{0},y_{0}\in V$ the following bound

$\displaystyle\lVert x_{0}-y_{0}\rVert\leq K$ (5)

holds. Now let $x,y:\mathbb{R}\to\mathbb{R}^{N}$ be solutions to (3)
originating from initial points $x_{0},y_{0}\in V\subset\mathbb{R}^{N}$
respectively. Reservoir system (3) is $V$-invariant so the solutions
$x(t),y(t)\in V$ for all $t>0$. Then for any $m\in M$ it follows from the
convexity of $V$ and the Mean Value Theorem (MVT) that there exists a curve
$s:\mathbb{R}^{+}\to\mathbb{R}^{N}$ such that

$\displaystyle
D_{x}F(s(t),\omega\phi^{t}(m))(x(t)-y(t))=F(x(t),\omega\phi^{t}(m))-F(y(t),\omega\phi^{t}(m)).$

Now we observe that

$\displaystyle\frac{d}{dt}\lVert x(t)-y(t)\rVert^{2}$
$\displaystyle=2(x(t)-y(t))^{\top}(\dot{x}(t)-\dot{y}(t))$
$\displaystyle=2(x(t)-y(t))^{\top}(F(x(t),\omega\phi^{t}(m))-F(y(t),\omega\phi^{t}(m)))$
$\displaystyle=2(x(t)-y(t))^{\top}D_{x}F(s(t),\omega\phi^{t}(m))(x(t)-y(t)),\
\text{(by MVT)}$
$\displaystyle=\frac{2(x(t)-y(t))^{\top}D_{x}F(s(t),\omega\phi^{t}(m))(x(t)-y(t))}{\lVert
x(t)-y(t)\rVert^{2}}\lVert x(t)-y(t)\rVert^{2}$
$\displaystyle\leq-2\delta\lVert x(t)-y(t)\rVert^{2}$

Then using separation of variables, and then bound (5), it follows that

$\displaystyle\lVert x(t)-y(t)\rVert^{2}\leq e^{-2\delta t}\lVert
x_{0}-y_{0}\rVert^{2}\leq e^{-2\delta t}K^{2}$

so

$\displaystyle\lVert x(t)-y(t)\rVert\leq e^{-\delta t}K=:\psi(t)$

where $\psi(t)$ is a positive function converging to 0. This establishes that
(3) is uniformly $V$-asymptotically stable with respect to
$\\{\omega\phi^{t}(m)\\}_{m\in M}$ which implies (3) admits a uniform
$V$-generalised synchronisation $f_{(\omega,\mathcal{V},F)}\in
C^{0}(M,\mathbb{R}^{N})$ by Theorem 2.8.

∎

It is possible for reservoir system (3) to admit several distinct generalised
synchronisations simultaneously. To see this, suppose that
$V_{1},\ldots,V_{n}$ are pairwise disjoint subsets of $\mathbb{R}^{N}$ and
that for each $i=1,\ldots,n$ (3) admits a $V_{i}$-generalised synchronisation
$f_{i(\omega,\mathcal{V},F)}\in C^{0}(M,V_{i})$. The existence of multiple
generalised synchronisations is closely related to the multi-ESP studied in
[4], and the discrete synchronisation results in [10]. A numerical experiment
which produces multiple generalised synchronisations is described and
performed in Section 3.

##  3 Numerical Illustration of Multiple GS

To demonstrate the existence of multi-GS we consider a hidden source system
defined by the vector field $\mathcal{V}$ on the manifold $M=\mathbb{R}^{2}$
defined by the ODE

$\displaystyle\dot{u}$ $\displaystyle=-v$ $\displaystyle\dot{v}$
$\displaystyle=u$ (6)

which has an associated group of evolution operators
$\\{\phi^{t}\in\text{Diff}^{\infty}(M)\ |\ t\in\mathbb{R}\\}$ defined by

$\displaystyle\phi^{t}(u_{0},v_{0})=(u_{0},v_{0})+\int_{0}^{t}(\dot{u}(\tau),\dot{v}(\tau))\
d\tau$

where $u(t),v(t)$ are the solutions of the ODEs (6) under the initial
condition $u(0)=u_{0},v(0)=v_{0}$. The vector field $\mathcal{V}$ represents
circular motion and the image of the trajectory $\\{\phi^{t}(0,1)\ |\
t\in\mathbb{R}\\}$ is a circle. The trajectory $\phi^{t}(0,1)$ is observed via
the function $\omega(u,v)=u$. The sequence of observations is shown in Figure
1, along with a so-called washout period to give the dynamics time to
synchronise.

Figure 1: The sequence of observations $\\{\omega\phi^{t}(0,1)\ |\
t\in(0,100)\\}$ with a washout period $(0,35)$ coloured blue.

The reservoir map $F:\mathbb{R}^{2}\times\mathbb{R}\to\mathbb{R}^{2}$ is
defined by

$\displaystyle F(x(t),y(t);z(t))=\begin{bmatrix}\sin(2\pi
x(t))+\lambda\sin(z(t))\\\ \sin(2\pi y(t))+\lambda\cos(z(t))\end{bmatrix}.$
(7)

Then we fix $\lambda=1$ and solve the reservoir system

$\displaystyle(\dot{x}(t),\dot{y}(t))=F(x(t),y(t);\omega\phi^{t}(0,1))$

under 4 initial conditions for the reservoir states $(x(0),y(0))=$
$(1/2,1/2)$, $(-1/2,1/2)$, $(1/2,-1/2)$, $(-1/2,-1/2)$. Each of these initial
points belong to one of four disjoint subsets $V_{1},V_{2},V_{3},V_{4}$ for
which the ODEs have a $V_{i}$-generalised synchronisation. The synchronised
dynamics are plotted in Figure 2.

Figure 2: The vector field defined by reservoir map (7) with $\lambda=0$ along
with the reservoir states $x(t),y(t)$ for times $t\in(35,100)$ originating
from 4 initial conditions $(x(0),y(0))=$ $(1/2,1/2)$, $(-1/2,1/2)$,
$(1/2,-1/2)$, $(-1/2,-1/2)$. This reveals the image of four distinct
generalised synchronisations each mapping the circle into a different region
$V_{i}$ of the reservoir space. The change in colour is indicative of the
change in time.

##  4 GS for Linear Reservoir Systems

###  4.1 Differentiable generalised synchronisation

In this paper, we have so far considered reservoir maps $F$ that are in
general nonlinear. In this section we will explore in detail the special case
where $F$ is a linear map. In particular, we define the linear reservoir
system by the ODE

$\displaystyle\dot{x}(t)=F(x(t),\omega\phi^{t}(m)):=-Ax(t)+C\omega\phi^{t}(m)$
(8)

where $A\in\mathbb{M}_{N\times N}(\mathbb{R})$ is symmetric positive definite,
which ensures that condition (4) in Theorem 2.9 is satisfied, and
$C\in\mathbb{M}_{N\times d}(\mathbb{R})$.

The linear reservoir system has been called a _next-generation reservoir
computer_ in a recent influential paper [6], and is much more amenable to
mathematical analysis than the generally nonlinear case studied so far.
Despite its simplicity, the linear reservoir computer has interesting
properties we will explore in the remainder of this paper. In the following
theorem, we will obtain a closed form expression for the reservoir states
$x(t)$ and associated generalised synchronisation $f_{(\omega,\mathcal{V},F)}$
for the linear system.

######  Proposition 4.1.

Suppose that $\omega\in C^{0}(M,\mathbb{R})$ is bounded. Then for any $m\in M$
and initial state $x(0)=x_{0}\in\mathbb{R}^{N}$ reservoir system (8) admits
the unique solution $x:\mathbb{R}\to\mathbb{R}^{N}$ defined by

$\displaystyle x(t)=\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau+t}(m)\
d\tau+e^{-At}\bigg{(}x_{0}-\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau}(m)\bigg{)},$
(9)

and associated generalised synchronisation

$\displaystyle
f_{(\omega,\mathcal{V},F)}(m)=\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau}(m)\
d\tau.$ (10)

###### Proof.

We start by expressing (8) as

$\displaystyle\dot{x}+Ax=C\omega\phi^{t}(m)$

and consider first the homogeneous equation

$\displaystyle\dot{x}+Ax=0.$

We solve this using separation of variables, and obtain the general solution
$x_{h}(t)=e^{-At}v$ where $v\in\mathbb{R}^{N}$ is a constant of integration.
We will now show that

$\displaystyle x_{i}(t)=\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau+t}(m)\
d\tau$

is an inhomogeneous solution to (3) by taking the derivative of $x_{i}$ and
verifying that the result satisfies (3).

$\displaystyle x_{i}(t)$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau+t}(m)\ d\tau$
$\displaystyle=-\int_{-t}^{\infty}e^{-A(t-u)}C\omega\phi^{u}(m)\
du\qquad(u=-\tau+t)$ (11)

hence

$\displaystyle\dot{x_{i}}(t)$
$\displaystyle=-\frac{d}{dt}\int_{-t}^{\infty}e^{-A(t-u)}C\omega\phi^{u}(m)\
du$ $\displaystyle=A\int^{\infty}_{-t}e^{-A(t-u)}C\omega\phi^{u}(m)\
du+C\omega\phi^{t}(m)\qquad\text{(Leibniz Integral Rule)}$
$\displaystyle=-Ax_{i}(t)+C\omega\phi^{t}(m)\qquad\text{by \eqref{xt}}.$

Then the full solution

$\displaystyle
x(t)=e^{-At}v+\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau+t}(m)\ d\tau$

must satisfy the initial condition

$\displaystyle x_{0}=x(0)=v+\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau}(m)\
d\tau$

hence

$\displaystyle v=x_{0}-\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau}(m)\
d\tau.$

Now we define the map $f_{(\omega,\mathcal{V},F)}:M\to\mathbb{R}^{N}$ by

$\displaystyle
f_{(\omega,\mathcal{V},F)}(m)=\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau}(m)\
d\tau$ (12)

for each $m\in M$ and observe that for any initial reservoir state
$x_{0}\in\mathbb{R}^{N}$ the solution $x(t)$ converges to
$f_{(\omega,\mathcal{V},F)}\phi^{t}(m)$. This establishes that
$f_{(\omega,\mathcal{V},F)}$ is a generalised synchronisation.

∎

In the following proposition, we will state conditions under which
$f_{(\omega,\mathcal{V},F)}\in C^{1}(M,\mathbb{R}^{N})$. This is especially
interesting because a differentiable synchronisation preserves geometric
information including the eigenvalues of fixed points and Lyapunov exponents
of the underlying system. This results in a higher quality reconstruction [32,
28].

######  Proposition 4.2.

Let $\omega\in C^{1}(M,\mathbb{R})$, $\phi\in\text{Diff}^{1}(M)$, and
$T_{m}\phi$ denote the differential of $\phi$ at $m\in M$. Let
$\sigma[A]_{\min}$ denote the smallest eigenvalue of $A$. Suppose there exists
$c\in(0,1)$ and $K>0$ such that

$\displaystyle\sup_{m\in M}\lVert T_{m}\phi^{-t}\rVert\leq
Ke^{c\sigma[A]_{\min}t},\qquad\sup_{m\in M}\lVert
D\omega(m)\rVert<\infty,\qquad\sup_{m\in M}\lVert\omega(m)\rVert<\infty$

then the $\mathbb{R}^{N}$-generalised synchronisation
$f_{(\omega,\mathcal{V},F)}$ is of class $C^{1}(M,\mathbb{R}^{N})$.

###### Proof.

Our goal is to prove that both $f_{(\omega,\mathcal{V},F)}$ and the derivative
$Df_{(\omega,\mathcal{V},F)}$ exist and are continuous. First of all,
$f_{(\omega,\mathcal{V},F)}$ exists and is continuous if

$\displaystyle\int_{0}^{\infty}\sup_{m\in M}\lVert
e^{-A\tau}C\omega\phi^{-\tau}(m)\rVert\ d\tau$

is finite, by lemma 8.3. To prove the integral is finite we observe that

$\displaystyle\int_{0}^{\infty}\sup_{m\in M}\lVert
e^{-A\tau}C(\omega\phi^{-\tau})(m)\rVert\ d\tau$
$\displaystyle\leq\int_{0}^{\infty}\lVert e^{-A\tau}C\rVert\ d\tau\ \sup_{m\in
M}\lVert\omega(m)\rVert$
$\displaystyle\leq\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\ d\tau\ \lVert
C\rVert\sup_{m\in M}\lVert\omega(m)\rVert\text{ by lemma \ref{eig_lemma}}$
$\displaystyle\leq\frac{\lVert C\rVert}{\sigma[A]_{\min}}\sup_{m\in
M}\lVert\omega(m)\rVert.$

The derivative $Df_{(\omega,\mathcal{V},F)}$ exists and is continuous if

$\displaystyle\int_{0}^{\infty}\sup_{m\in M}\lVert
e^{-A\tau}CD(\omega\phi^{-\tau})(m)\rVert\ d\tau$

is finite, using lemma 8.3 once again. To prove this second integral is finite
we observe that

$\displaystyle\int_{0}^{\infty}\sup_{m\in M}\lVert
e^{-A\tau}CD(\omega\phi^{-\tau})\rVert\ d\tau(m)$ $\displaystyle\leq$
$\displaystyle\int_{0}^{\infty}\lVert e^{-A\tau}\rVert\lVert C\rVert\sup_{m\in
M}\lVert D(\omega\phi^{-\tau})(m)\rVert\ d\tau$ $\displaystyle\leq$
$\displaystyle\lVert
C\rVert\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\sup_{m\in M}\lVert
D(\omega\phi^{-\tau})(m)\rVert\ d\tau\ \text{ by lemma \ref{eig_lemma}}$
$\displaystyle=$ $\displaystyle\lVert
C\rVert\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\sup_{m\in M}\lVert
D\omega\phi^{-\tau}(m)T_{m}\phi^{-\tau}\rVert\ d\tau$ $\displaystyle\leq$
$\displaystyle\lVert
C\rVert\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\sup_{m\in M}\lVert
D\omega(m)\rVert\sup_{m\in M}\lVert T_{m}\phi^{-\tau}\rVert\ d\tau$
$\displaystyle\leq$ $\displaystyle\lVert
C\rVert\int_{0}^{\infty}Ke^{-(1-c)\sigma[A]_{\min}\tau}\ d\tau\ \sup_{m\in
M}\lVert D\omega(m)\rVert$ $\displaystyle=$ $\displaystyle\frac{K\lVert
C\rVert}{(1-c)\sigma[A]_{\min}}\sup_{m\in M}\lVert D\omega(m)\rVert.$

∎

###  4.2 Embedding the fixed points

In a highly influential paper written in 1981, Floris Takens [37] proved that
the delay observation map is an embedding for a generic pair $(\omega,\phi)$
of observation functions $\omega$ and discrete time source system $\phi$. This
launched the field of embedology [34] which explores novel and useful ways to
embed a source system into $\mathbb{R}^{N}$ using only a time series of scalar
(or low dimensional) observations of the source system. Many authors [14] [11]
[40] [41] [35] have observed that reservoir computing with observations from a
dynamical system is a special type of embedology; as long as the generalised
synchronisation $f_{(\omega,\mathcal{V},F)}$ is an embedding.

It is especially desirable that the generalised synchronisation
$f_{(\omega,\mathcal{V},F)}$ is an embedding because this property allows a
topologically faithful reconstruction of the source dynamics, which allows for
future forecasting and other forms of learning. In fact an embedding creates a
vector field in the reservoir space that is diffeomorphic to the vector field
of the source system. It is very challenging to prove that a generalised
synchronisation is embedding globally on $M$, so in this paper we prove a
simpler result: that for generic observation functions $\omega$ the
generalised synchronisation $f_{(\omega,\mathcal{V},F)}$ is an embedding on
the isolated fixed points of $\mathcal{V}$. First of all this simpler result
is a necessary first step in the quest to prove a global embedding. Moreover,
isolated fixed points are of particular interest in their own right because
they are each contained by a small neighbourhood of $M$ on which the vector
field $\mathcal{V}$ is approximately linear. Furthermore, when a trajectory of
the source dynamics is at a point $m$ close to a fixed point $m^{*}$, a
nonlinear reservoir system converges to the linear reservoir system

$\displaystyle\dot{x}(t)$
$\displaystyle=D_{x}F(f_{(\omega,\mathcal{V},F)}(m^{*}),\omega(m^{*}))x(t)+D_{z}F(f_{(\omega,\mathcal{V},F)}(m^{*}),\omega(m^{*}))\omega\phi^{t}(m)$
$\displaystyle=-Ax(t)+C\omega\phi^{t}(m)$

where

$\displaystyle A$
$\displaystyle=D_{x}F(f_{(\omega,\mathcal{V},F)}(m^{*}),\omega(m^{*})),$
$\displaystyle C$
$\displaystyle=D_{z}F(f_{(\omega,\mathcal{V},F)}(m^{*}),\omega(m^{*})).$

We can easily evaluate the GS evaluated at a fixed point $m^{*}$

$\displaystyle f_{(\omega,\mathcal{V},F)}(m^{*})$
$\displaystyle=\int^{\infty}_{0}e^{-A\tau}C\omega(m^{*})d\tau=A^{-1}C\omega(m^{*}).$

Now, before we prove that $f_{(\omega,\mathcal{V},F)}$ is an embedding of the
fixed points for generic $\omega\in C^{1}(M,\mathbb{R})$, we will define an
immersion, an embedding, and the word generic.

######  Definition 4.3.

Let $T_{m}M$ denote the tangent space of $M$ at $m$. A map
$f:M\to\mathbb{R}^{N}$ is an immersion at $m\in M$ if the derivative
$Df(m):T_{m}M\to\mathbb{R}^{N}$ exists and is full rank.

######  Definition 4.4.

Let $X$ be a compact subset of $M$. Then a map $f:M\to\mathbb{R}^{N}$ is an
embedding on $X$ if $f$ is an immersion for each $m\in X$ and the restriction
$f\rvert_{X}$ is injective.

######  Remark 4.5.

If $X$ is a compact subset of $M$ and $f:M\to\mathbb{R}^{N}$ is an embedding
on $X$ then there is an open subset $\Omega\subset M$ such that
$X\subset\Omega\subset M$ and the restriction $f\rvert_{\Omega}$ is an
embedding. Furthermore, $f\rvert_{\Omega}$ is a diffeomorphism onto its image.

######  Definition 4.6.

Let $X$ be compact subset of $M$. A property that holds on a dense open subset
of $X$ is a generic property of $X$.

Now to prove that $f_{(\omega,\mathcal{V},F)}$ is an embedding on the fixed
points for generic $\omega\in C^{1}(M,\mathbb{R})$ we will prove:

  1. 1.

In lemma 4.7 that $f_{(\omega,\mathcal{V},F)}$ is an embedding on the fixed
points for observation functions on an open (possibly empty) subset of
$C^{1}(M,\mathbb{R})$.

  2. 2.

In theorem 4.8 that $f_{(\omega,\mathcal{V},F)}$ is an embedding on the fixed
points for observation functions in a dense subset of $C^{1}(M,\mathbb{R})$.
This combined with lemma 4.7 completes the proof.

Our strategy is very similar to Huke’s [16] strategy to prove Takens’ Theorem,
and the strategy appearing in [11] to prove a similar result in the discrete
time case.

######  Lemma 4.7.

Let $\omega\in C^{1}(M,\mathbb{R})$, $\phi\in\text{Diff}^{1}(M)$, and
$T_{m}\phi$ denote the differential of $\phi$ at $m\in M$. Let
$\sigma[A]_{\min}$ denote the smallest eigenvalue of $A$. Suppose there exists
$c\in(0,1)$ and $K>0$ such that

$\displaystyle\sup_{m\in M}\lVert T_{m}\phi^{-t}\rVert\leq
Ke^{c\sigma[A]_{\min}t},\qquad\sup_{m\in M}\lVert
D\omega(m)\rVert<\infty,\qquad\sup_{m\in M}\lVert\omega(m)\rVert<\infty$

Let $X$ be a compact subset of $M$. Then for an open (possibly empty) subset
of $C^{1}(M,\mathbb{R})$ the GS $f_{(\omega,\mathcal{V},F)}$ is an embedding
on $X$.

###### Proof.

Embeddings on $X$ are an open set in $C^{1}(M,\mathbb{R}^{N})$ by Theorem 1.4
in [15]. Furthermore, the inverse image of any open set under a continuous map
is open. Thus, we will show that the map $\Psi:C^{1}(M,\mathbb{R})\to
C^{1}(M,\mathbb{R}^{N})$ defined by

$\displaystyle\Psi(\omega)=f_{(\omega,\mathcal{V},F)}$

is continuous, and this will complete the proof. To this end, consider a
sequence $\\{\omega_{n}\\}_{n\in\mathbb{N}}$ that converges to $\omega_{*}$ in
the $C^{1}(M,\mathbb{R})$ topology. Then for any $\epsilon>0$ there exists
$n^{\prime}>0$ such that for any $n>n^{\prime}$

$\displaystyle\sup_{m\in
M}\lVert\omega_{n}(m)-\omega_{*}(m)\rVert<\frac{\epsilon\sigma[A]_{\min}}{2\lVert
C\rVert},$ $\displaystyle\sup_{m\in M}\lVert
D\omega_{n}(m)-D\omega_{*}(m)\rVert<\frac{\epsilon(1-c)\sigma[A]_{\min}}{2K\lVert
C\rVert}.$

Then our goal is to establish that $\Psi(\omega_{n})$ converges to
$\Psi(\omega_{*})$ in the $C^{1}(M,\mathbb{R}^{N})$ topology, i.e. show that

$\displaystyle\lVert\Psi(\omega_{n})-\Psi(\omega_{*})\rVert_{C^{1}(M,\mathbb{R})}<\epsilon.$

Now

$\displaystyle\lVert\Psi(\omega_{n})-\Psi(\omega_{*})\rVert_{C^{1}(M,\mathbb{R})}$
$\displaystyle=$ $\displaystyle\lVert
f_{(\omega_{n},\mathcal{V},F)}-f_{(\omega^{*},\mathcal{V},F)}\rVert_{C^{1}(M,\mathbb{R}^{N})}$
$\displaystyle=$
$\displaystyle\bigg{\lVert}\int_{0}^{\infty}e^{-A\tau}C(\omega_{n}\phi^{-\tau}-\omega^{*}\phi^{-\tau})\
d\tau\bigg{\rVert}_{C^{1}(M,\mathbb{R}^{N})}$ $\displaystyle\leq$
$\displaystyle\int_{0}^{\infty}\lVert
e^{-A\tau}C(\omega_{n}\phi^{-\tau}-\omega^{*}\phi^{-\tau})\rVert_{C^{1}(M,\mathbb{R}^{N})}\
d\tau$ $\displaystyle\leq$ $\displaystyle\int_{0}^{\infty}\sup_{m\in M}\lVert
e^{-A\tau}C(\omega_{n}(m)-\omega^{*}(m))\rVert\ d\tau$
$\displaystyle\qquad+\int_{0}^{\infty}\sup_{m\in M}\lVert
e^{-A\tau}C(D(\omega_{n}\phi^{-\tau})(m)-D(\omega^{*}\phi^{-\tau})(m))\rVert\
d\tau$ $\displaystyle\leq$ $\displaystyle\int_{0}^{\infty}\lVert
e^{-A\tau}\rVert\ d\tau\ \lVert C\rVert\sup_{m\in
M}\lVert\omega_{n}(m)-\omega^{*}(m)\rVert$ $\displaystyle\qquad+\lVert
C\rVert\int_{0}^{\infty}\lVert e^{-A\tau}\rVert\sup_{m\in
M}\lVert(D(\omega_{n}\phi^{-\tau})(m)-D(\omega^{*}\phi^{-\tau})(m))\rVert\
d\tau$ $\displaystyle\leq\lVert
C\rVert\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\ d\tau\ \sup_{m\in
M}\lVert\omega_{n}(m)-\omega^{*}(m)\rVert$ $\displaystyle\qquad+\lVert
C\rVert\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\sup_{m\in
M}\lVert(D(\omega_{n}\phi^{-\tau})(m)-D(\omega^{*}\phi^{-\tau})(m))\rVert\
d\tau$ by lemma 8.4 $\displaystyle\leq\frac{\lVert
C\rVert}{\sigma[A]_{\min}}\sup_{m\in
M}\lVert\omega_{n}(m)-\omega^{*}(m)\rVert$ $\displaystyle\qquad+\lVert
C\rVert\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\sup_{m\in
M}\lVert(D\omega_{n}\phi^{-\tau}(m)-D\omega^{*}\phi^{-\tau}(m))T_{m}\phi^{-\tau}\rVert\
d\tau$ $\displaystyle\leq\frac{\lVert C\rVert}{\sigma[A]_{\min}}\sup_{m\in
M}\lVert\omega_{n}(m)-\omega^{*}(m)\rVert$ $\displaystyle\qquad+\lVert
C\rVert\int_{0}^{\infty}e^{-\sigma[A]_{\min}\tau}\sup_{m\in M}\lVert
T_{m}\phi^{-\tau}\rVert\ d\tau\ \sup_{m\in
M}\lVert(D\omega_{n}(m)-D\omega^{*}(m))$ $\displaystyle\leq\frac{\lVert
C\rVert}{\sigma[A]_{\min}}\sup_{m\in
M}\lVert\omega_{n}(m)-\omega^{*}(m)\rVert$ $\displaystyle\qquad+\lVert
C\rVert\int_{0}^{\infty}Ke^{-(1-c)\sigma[A]_{\min}\tau}\ d\tau\ \sup_{m\in
M}\lVert(D\omega_{n}(m)-D\omega^{*}(m))$ $\displaystyle\leq\frac{\lVert
C\rVert}{\sigma[A]_{\min}}\sup_{m\in
M}\lVert\omega_{n}(m)-\omega^{*}(m)\rVert$ $\displaystyle\qquad+\frac{K\lVert
C\rVert}{(1-c)\sigma[A]_{\min}}\ \sup_{m\in
M}\lVert(D\omega_{n}(m)-D\omega^{*}(m))$
$\displaystyle<\frac{\epsilon}{2}+\frac{\epsilon}{2}=\epsilon.$

∎

To show that $f_{(\omega,\mathcal{V},F)}$ is embedding about the fixed points
for a dense set of observation functions, we will prove that for an arbitrary
observation function $\omega$, we can always make an arbitrarily small
perturbation $\omega^{\prime}$ such that the resulting generalised
synchronisation $f_{(\omega^{\prime},\mathcal{V},F)}$ is an embedding on the
fixed points

######  Theorem 4.8.

Let $\omega\in C^{1}(M,\mathbb{R})$, $\phi\in\text{Diff}^{1}(M)$, and
$T_{m}\phi$ denote the differential of $\phi$ at $m\in M$. Let
$\sigma[A]_{\min}$ denote the smallest eigenvalue of $A$. Suppose there exists
$c\in(0,1)$ and $K>0$ such that

$\displaystyle\sup_{m\in M}\lVert T_{m}\phi^{-t}\rVert\leq
Ke^{c\sigma[A]_{\min}t},\qquad\sup_{m\in M}\lVert
D\omega(m)\rVert<\infty,\qquad\sup_{m\in M}\lVert\omega(m)\rVert<\infty$

Let the dimension of $M$ be $q$ and dimension of the reservoir space $N\geq
q$. Suppose that the smooth vector field $\mathcal{V}$ on $M$ admits a finite
number of fixed points. For each fixed point $m\in M$ let $J_{m}:T_{m}M\to
T_{m}M$ denote the Jacobian at $m$. Suppose that for each fixed point $m$

  1. 1.

The eigenvalues of $J_{m}$ denoted $\lambda_{1},\ldots\lambda_{q}$ are
distinct.

  2. 2.

The vectors

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

are linearly independent.

Then for generic $\omega\in C^{1}(M,\mathbb{R})$ the GS
$f_{(\omega,\mathcal{V},F)}$ is an embedding on the fixed points.

###### Proof.

Let $m\in M$ be a fixed point. Suppose that $v_{1},\ldots,v_{q}$ denote the
eigenvectors of $J_{m}$ associated to the eigenvalues
$\lambda_{1},\ldots,\lambda_{q}$. Let $\psi\in C^{\infty}(M,\mathbb{R})$ be a
smooth bump function with support that contains $m$ and no other fixed point.
Furthermore suppose that

$\displaystyle D\psi(m)=w^{\top}$

where $w$ is the unique solution to

$\displaystyle\begin{bmatrix}v_{1}^{\top}\\\ v_{2}^{\top}\\\ \vdots\\\
v_{q}^{\top}\end{bmatrix}w=\begin{bmatrix}1\\\ 1\\\ \vdots\\\ 1\end{bmatrix}$

Now define for $\epsilon>0$ the perturbed observation function

$\displaystyle\omega^{\prime}(m)=\omega+\epsilon\psi(m).$

Then the perturbed GS has the form

$\displaystyle
f_{(\omega^{\prime},\mathcal{V},F)}(m)=\int_{0}^{\infty}e^{-A\tau}C\omega^{\prime}\phi^{-\tau}(m)\
d\tau$

so

$\displaystyle Df_{(\omega^{\prime},\mathcal{V},F)}(m)$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}CD(\omega^{\prime}\phi^{-\tau})(m)\
d\tau$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}CD\omega^{\prime}(m)T_{m}\phi^{-\tau}\
d\tau$ $\displaystyle=\int_{0}^{\infty}e^{-A\tau}CD\omega(m)T_{m}\phi^{-\tau}\
d\tau+\epsilon\int_{0}^{\infty}e^{-A\tau}Cw^{\top}T_{m}\phi^{-\tau}\ d\tau.$

Now $m$ is a fixed point so

$\displaystyle T_{m}\phi^{-\tau}=e^{-J_{m}\tau}.$

Now the eignevectors $v_{j}$ of the Jacobian $J_{m}$ are linearly independent
so it suffices to find arbitrarily small $\epsilon>0$ such that

$\displaystyle\\{Df_{(\omega^{\prime},\mathcal{V},F)}v_{j}\\}_{j=1,\ldots,q}$

are linearly independent. Now

$\displaystyle Df_{(\omega^{\prime},\mathcal{V},F)}(m)v_{j}$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}CD\omega(m)T_{m}\phi^{-\tau}v_{j}\
d\tau+\epsilon\int_{0}^{\infty}e^{-A\tau}Cw^{\top}T_{m}\phi^{-\tau}v_{j}\
d\tau$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}CD\omega(m)e^{-J_{m}\tau}v_{j}\
d\tau+\epsilon\int_{0}^{\infty}e^{-A\tau}Cw^{\top}e^{-J_{m}\tau}v_{j}\ d\tau$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}CD\omega(m)e^{-\lambda_{j}\tau}v_{j}\
d\tau+\epsilon\int_{0}^{\infty}e^{-A\tau}Cw^{\top}e^{-\lambda_{j}\tau}v_{j}\
d\tau$
$\displaystyle=\int_{0}^{\infty}e^{-(A+\lambda_{j}\mathbb{I})\tau}CD\omega(m)v_{j}\
d\tau+\epsilon\int_{0}^{\infty}e^{-(A+\lambda_{j}\mathbb{I})\tau}Cw^{\top}v_{j}\
d\tau$
$\displaystyle=\int_{0}^{\infty}e^{-(A+\lambda_{j}\mathbb{I})\tau}CD\omega(m)v_{j}\
d\tau+\epsilon\int_{0}^{\infty}e^{-(A+\lambda_{j}\mathbb{I})\tau}C\ d\tau$
$\displaystyle=\int_{0}^{\infty}e^{-(A+\lambda_{j}\mathbb{I})\tau}CD\omega(m)v_{j}\
d\tau+\epsilon(A+\lambda_{j}\mathbb{I})^{-1}C.$

Now the vectors

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

are linearly independent by assumption. Hence we can choose a sufficiently
small $\epsilon>0$ so that

$\displaystyle\bigg{\\{}\int_{0}^{\infty}e^{-(A+\lambda_{j}\mathbb{I})\tau}CD\omega(m)v_{j}\
d\tau+\epsilon(A+\lambda_{j}\mathbb{I})^{-1}C\bigg{\\}}_{j=1,\ldots,q}$

are linearly independent, by lemma 8.5. Thus we have an immersion of the point
$m$. This immersion has no effect on the other fixed points, we so can repeat
this procedure on each of the finitely many fixed points in turn without
spoiling the immersion on any of the previous points. By Theorem 1.1 in [15]
the immersions of the fixed points form an open set, so any sufficiently small
perturbation of the observation function will preserve the immersion on the
fixed points. Thus, we will construct an arbitrarily small perturbation
$\omega^{\prime}$ of an arbitrary observation function $\omega$ which ensures
that $f_{(\omega^{\prime},\mathcal{V},F)}$ restricted to the fixed points is
injective.

For each fixed point $m_{i}\in M$ define a smooth bump function
$\varphi_{i}\in C^{\infty}(M,\mathbb{R})$ such that $\varphi_{i}(m_{i})=1$ and
$\text{supp}(\varphi_{i})$ are disjoint. Then define the perturbed observation
function

$\displaystyle\omega^{\prime}(m)=\omega(m)+\sum_{i}\epsilon_{i}\varphi_{i}(m)$

for $\epsilon_{i}>0$. Now
$\omega^{\prime}\phi^{-\tau}(m_{i})=\omega^{\prime}(m_{i})$ for all $\tau$
because $m_{i}$ is a fixed point so

$\displaystyle f_{(\omega^{\prime},\varphi,F)}(m_{i})$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}C\omega^{\prime}\phi^{-\tau}(m_{i})\
d\tau$ $\displaystyle=\omega^{\prime}(m_{i})\int_{0}^{\infty}e^{-A\tau}\ d\tau
C$ $\displaystyle=\omega^{\prime}(m_{i})A^{-1}C$
$\displaystyle=(\omega(m_{i})+\epsilon_{i})A^{-1}C.$

Then we construct arbitrarily small $\\{\epsilon_{i}\\}$ such that

$\displaystyle\omega(m_{i})+\epsilon_{i}=\omega(m_{j})+\epsilon_{j}\iff i=j$

and this proves that $f_{(\omega,\mathcal{V},F)}$ is an embedding of the fixed
points for a dense set of observation functions. Then lemma 4.7 gives an
embedding of the fixed points for generic $\omega$ and the proof is complete.
∎

One fundamental feature of reservoir computers that distinguish them from
ordinary neural networks is that the reservoir neurons are randomly generated
prior to training, and subsequently untouched by training. We are therefore
interested in the likelihood that conditions 1. and 2. of Theorem 4.8 hold for
randomly generated $A$ and $C$. It is shown in lemma 8.2 that the conditions
hold almost surely for suitably randomly generated $A,C$.

##  5 A Central Limit Theorem for Reservoir Computing

If $f_{(\omega,\mathcal{V},F)}$ is injective, we can define the autonomous
reservoir system

$\displaystyle\dot{x}=F(x,\omega f_{(\omega,\mathcal{V},F)}^{-1}(x))$ (13)

which in the linear case is

$\displaystyle\dot{x}=-Ax+C\omega f_{(\omega,\mathcal{V},F)}^{-1}(x).$

If $f_{(\omega,\mathcal{V},F)}$ is an embedding on some open subset
$\Omega\subset M$ then the autonomous system (13) defines a vector field on
$f(\Omega)\subset\mathbb{R}^{N}$ which is diffeomorphic to the source vector
field $\mathcal{V}$ on $\Omega$.

We cannot evaluate the map $\omega f^{-1}_{(\omega,\mathcal{V},F)}$ directly,
but we can approximate it with a random neural network
$h:\mathbb{R}^{N}\times\Theta\to\mathbb{R}$, where $\Theta$ is a parameter
space from which we draw random weights and biases. A common choice of
‘neuron’ $h$ comprising a random neural network is

$\displaystyle h(x,\theta)=h(x,\alpha,\beta)=\tanh(\alpha^{\top}x+\beta)$

where $\alpha\in\mathbb{R}^{N}$ and $\beta\in\mathbb{R}$ are random weights
and biases. Instead of using $\tanh$ we could use other activation functions
including the relu activation or radial basis activations. To improve
generality we will approximate an arbitrary target function
$u:\mathbb{R}^{N}\to\mathbb{R}$, motivated by the special case $u=\omega
f^{-1}_{(\omega,\mathcal{V},F)}$. The strategy is to first randomly generate
weights and biases $\\{\theta_{i}\\}_{1,\ldots,D}$, then take samples
$x_{j}=f(m_{j})$ from a trajectory in the reservoir space. Then minimise over
$\\{w_{i}\in\mathbb{R}\\}_{i,\ldots,D}$ the loss

$\displaystyle
L(w_{1},\ldots,w_{D})=\frac{1}{\ell}\sum_{j=1}^{\ell}\bigg{\lVert}u(x_{j})-\frac{1}{D}\sum_{i=1}^{D}w_{i}h(x_{j},\theta_{i})\bigg{\rVert}^{2}$

then we hope

$\displaystyle u(x)\approx\frac{1}{D}\sum_{i=1}^{D}w_{i}h(x,\theta_{i})$ (14)

for all $x\in f(M)$. In this section, will study the error on the
approximation (14) as we increase the number of neurons $D$. This is closely
related to the work in [8] and [7]. We do not study the error on (14) as the
number of sample points $\ell$ increases, which is addressed in [13] for
discrete ergodic source systems.

We show in Theorem 5.2 that if the target $u(x)$ can be expressed as a
weighted integral over $\theta\in\Theta$ of $h(x,\theta)$ then $h$ can
approximate $u$ with asymptotically normal approximation error, with variance
proportional to $1/D$. This follows directly from the central limit theorem
(CLT), which we will state here.

######  Theorem 5.1.

(Central Limit Theorem) Let $(\Theta,\mathcal{F},\mathbb{P})$ be a probability
space and $\\{\boldsymbol{Y}_{i}\\}_{i\in\mathbb{N}}$ a collection of IID
$\Theta$-valued random variables uniformly distributed with respect to
$\mathbb{P}$. Suppose that

$\displaystyle\mu:=\mathbb{E}[\boldsymbol{Y}_{i}]\qquad\text{and}\qquad\sigma^{2}:=\text{Var}[\boldsymbol{Y}_{i}]$

are finite. Then

$\displaystyle\mu-\frac{1}{D}\sum^{N}_{i=1}\boldsymbol{Y}_{i}\xrightarrow[D\to\infty]{}\mathcal{N}(0,\sigma^{2}/D)$

in the sense that

$\displaystyle\lim_{N\to\infty}\mathbb{P}\bigg{[}\mu-\frac{1}{D}\sum_{i=1}^{D}\boldsymbol{Y}_{i}\in(a,b)\bigg{]}=\frac{\sqrt{D}}{\sqrt{2\pi}\sigma}\int_{a}^{b}\exp\bigg{(}-\frac{z^{2}D}{2\sigma^{2}}\bigg{)}\
dz$

for any real interval $(a,b)\subset\mathbb{R}$.

We are now are ready to state the result.

######  Theorem 5.2.

(CLT for random neural networks) Let $(\Theta,\mathcal{F},\mathbb{P})$ be a
probability space and $x\in\mathbb{R}^{N}$. Let
$h:\mathbb{R}^{N}\times\Theta\to\mathbb{R}$ be a map where $h(x,\cdot)$ is
$\mathbb{P}$-measurable. Let $\boldsymbol{\theta}$ and
$\\{\boldsymbol{\theta}_{i}\\}_{i\in\mathbb{N}}$ be real valued IID
$\Theta$-valued random variables uniformly distributed with respect to
$\mathbb{P}$. Suppose that the target function $u:\mathbb{R}^{N}\to\mathbb{R}$
is of the form

$\displaystyle
u(x)=\mathbb{E}[w(\boldsymbol{\theta})h(x,\boldsymbol{\theta})]=\int_{\Theta}w(\theta)h(x,\theta)\
d\mathbb{P}(\theta)$

for some $\mathbb{P}$-measurable $w:\Theta\to\mathbb{\mathbb{R}}$, and that

$\displaystyle\sigma^{2}(x):=\text{Var}[w(\boldsymbol{\theta})h(x,\boldsymbol{\theta})]=\int_{\Theta}\lvert
w(\theta)h(x,\theta)-u(x)\rvert^{2}\ d\mathbb{P}(\theta)$

is finite for all $x\in\mathbb{R}^{N}$. Then

$\displaystyle
u(x)-\frac{1}{D}\sum_{i=1}^{D}w(\boldsymbol{\theta}_{i})h(x,\boldsymbol{\theta}_{i})\xrightarrow[D\to\infty]{}\mathcal{N}(0,\sigma^{2}(x)/D)$

in the sense that

$\displaystyle\lim_{D\to\infty}\mathbb{P}\bigg{[}u(x)-\frac{1}{D}\sum_{i=1}^{D}w(\boldsymbol{\theta}_{i})h(x,\boldsymbol{\theta}_{i})\in(a,b)\bigg{]}=\frac{\sqrt{D}}{\sqrt{2\pi}\sigma(x)}\int_{a}^{b}\exp\bigg{(}-\frac{z^{2}D}{2\sigma^{2}(x)}\bigg{)}\
dz$

for any real interval $(a,b)\subset\mathbb{R}$.

###### Proof.

We let

$\displaystyle\boldsymbol{Y}_{i}(x):=w(\boldsymbol{\theta}_{i})h(x,\boldsymbol{\theta}_{i})$

then the result follows immediately from the central limit theorem (CLT). ∎

Theorem 5.2 (roughly) states that when using a random neural network composed
of $D$ neurons to approximate a target function $u$, the error is
asymptotically normal with a variance $\sigma^{2}(x)/D$, so the reciprocal
error is polynomial in $D$. Furthermore, the variance $\sigma^{2}(x)/D$ does
not depend explicitly on the dimension $N$ of the reservoir dynamics, so the
random neural network does not suffer the curse of dimensionality. It is
important to remark that the result holds under a nontrivial assumption that
the target function $u$ is a weighted integral of the random neural network
$h$. It is shown in [8] and [7] that a broad class of target functions $u$ can
be expressed as weighted integrals of random neural networks with relu
activations.

##  6 Numerical Illustration of an Embedding

The theoretical portion of this paper covers two major ideas

  1. 1.

Given scalar observations of a source system and a linear reservoir system, it
is possible to embed the source dynamics (about a fixed point) into the
reservoir space.

  2. 2.

If an embedding is achieved it follows from the CLT that we can train a
reservoir computer to approximate the future trajectory of the observations.
Furthermore the autonomous dynamics of the trained reservoir computer will be
diffeomorphic to the source dynamics, hence the reservoir dynamics will
inherit some geometrical properties of the source dynamics, including the
eigenvalues of the linearisation about the fixed point.

In this section we will illustrate both of these these ideas with a numerical
experiment. In particular we will take the Lorenz-63 system [23]

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

as the source system. Suppose we can observe a trajectory of the Lorenz system
via the observation function $\omega(m)=\omega(\xi,\upsilon,\zeta)=\xi$. Then
our first goal is to use a linear reservoir system to embed the fixed point
$m^{*}=(6\sqrt{2},6\sqrt{2},27)$ into the reservoir space using only an
observed trajectory in a neighbourhood of $m^{*}$. This illustrates idea 1.

To achieve this embedding we generate random reservoir matrices $A,C$ using
the following Python script

⬇

import numpy as np

from scipy.stats import ortho_group

N = 7

D = 300

# Q is a random orthogonal matrix drawn from the Haar

# distribution.

Q = ortho_group.rvs(dim=N)

# generate random positive definite A

A = np.random.rand(N,)

A = np.diag(A)

A = Q @ A @ np.transpose(Q)*30

# generate random C

C = np.random.rand(N,)-0.5

Then we integrate the Lorenz equations together with the equations of the
linear reservoir system

$\displaystyle\dot{\xi}$ $\displaystyle=10(\upsilon-\xi)$ (15)
$\displaystyle\dot{\upsilon}$ $\displaystyle=\xi(28-\zeta)-\upsilon$
$\displaystyle\dot{\zeta}$ $\displaystyle=\xi\upsilon-(8/3)\zeta$
$\displaystyle\dot{x}$ $\displaystyle=Ax+C\xi$

using a numerical integrator implemented in the Python Scipy library

⬇

scipy.integrate.RK45(fun=system, t0=0,

y0=y0, t_bound=200, rtol=1e-9)

with system (15) initial time t0=0, final time t_bound=200, relative error
tolerance rtol=1e-9, and initial point y0 defined by

$\displaystyle\xi^{*}$ $\displaystyle=6\sqrt{2}$ $\displaystyle\upsilon^{*}$
$\displaystyle=6\sqrt{2}$ $\displaystyle\zeta^{*}$ $\displaystyle=27$
$\displaystyle x^{*}$ $\displaystyle=A^{-1}C\xi$

so that the Lorenz system and reservoir system is initialised, up to machine
precision, at $m^{*}$ and $f(m^{*})$ respectively. This returns a finite set
of observations $\xi_{1},\ldots,\xi_{\ell}$ and reservoir states
$x_{1},\ldots,x_{\ell}$. Because the fixed point is unstable, the error on the
initial point due to machine precision causes the trajectory to spiral out
from the fixed point and fill the attractor. The observations
$\xi_{1},\ldots,\xi_{\ell}$ are plotted in Figure 3, the full Lorenz system is
plotted in Figure 4, and the reservoir states $x_{1},\ldots,x_{\ell}$
projected onto the first 3 principal components are plotted in Figure 5.

Figure 3: The observations $\xi_{1},\ldots,\xi_{\ell}$ of the Lorenz system
over the time interval $(0,200)$. The coloring of the trajectory is consistent
with Figures 4 and 5. Figure 4: The Lorenz-63 system is the source system of
the experiment. The initial point of the trajectory is
$m^{*}=(6\sqrt{2},6\sqrt{2},27)$ up to machine precision. Because the fixed
point is unstable, the error on the initial point due to machine precision
cause the trajectory to spiral out from the fixed points and fill the
attractor. The coloring of the trajectory is consistent with Figures 3 and 5.
Figure 5: The reservoir states projected onto their first three principal
components in the reservoir space. No axis labels or tickmarks are included
because these metric quantities are not preserved by the embedding, and are
(for our purposes) meaningless. The coloring of the trajectory is consistent
with Figures 3 and 4. The reservoir dynamics appear to be diffeomorphic to the
Lorenz dynamics, suggesting the GS $f_{(\omega,\mathcal{V},F)}$ exists, and is
an embedding. It is interesting to note that it appears we have achieved a
global embedding, though the theory we have established in this paper is
limited to a local embedding about a fixed point.

If an embedding is achieved on $m^{*}=(6\sqrt{2},6\sqrt{2},27)$ then it
follows from the CLT that we can train the reservoir computer on a
neighbourhood of $f(m^{*})=f(6\sqrt{2},6\sqrt{2},27)$ to predict the future
trajectory of the observations. Furthermore, the autonomous dynamics of the
trained reservoir computer may be diffeomorphic to the Lorenz system in the
vicinity of the fixed point. In this case the geometrical properties of the
fixed point $m^{*}$, including the eigenvalues of the linearisation
$J_{m^{*}}$ about $m^{*}$, will be preserved by the reservoir system. Our goal
is therefore to compute the eigenvalues of the linearisation about $f(m^{*})$
of the reservoir dynamics and compare them to those of the linearisation about
$m^{*}$ of the Lorenz dynamics. This illustrates idea 2.

The linearisation of the Lorenz system about the fixed point
$m^{*}=(6\sqrt{2},6\sqrt{2},27)$ is the matrix

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

and the eigenvalues of the matrix are plotted in Figure 6.

Our goal now is to train a linear reservoir system to approximate the future
trajectory, find the eigenvalues of the linearisation about $f(m^{*})$ of the
autonomous reservoir system, and compare them to the eigenvalues of the Lorenz
system.

For each $i=1,\ldots,D$ we randomly generate independent components of
$\alpha_{i}\in\mathbb{R}^{N}$ and $\beta_{i}\in\mathbb{R}$ from the uniform
distribution $U\sim[-0.5,0.5]$. Then we find real numbers
$w_{1},\ldots,w_{\ell}$ that minimise the loss

$\displaystyle
L(w_{1},\ldots,w_{D})=\frac{1}{\ell}\sum_{j=1}^{\ell}\bigg{\lVert}\xi_{j}-\frac{1}{D}\sum_{i=1}^{D}w_{i}\tanh(\alpha_{i}^{\top}x_{j}+\beta_{i})\bigg{\rVert}^{2}$

using scipy.sparse.linalg.lsqr with damp=0. Then we have the autonomous system

$\displaystyle\dot{x}=-Ax+C\bigg{(}\frac{1}{D}\sum_{i=1}^{D}w_{i}\tanh(\alpha_{i}^{\top}x+\beta_{i})\bigg{)}$

which has Jacobian

$\displaystyle
J^{\prime}_{x}=-A+C\bigg{(}\frac{1}{D}\sum_{i=1}^{D}w_{i}\text{sech}^{2}(\alpha_{i}^{\top}x+\beta_{i})\alpha_{i}^{\top}\bigg{)}.$

Then we compute the eigenvalues of the Jacobian at the fixed point
$x^{*}=A^{-1}C(6\sqrt{2})$. We compare these eigenvalues to the eigenvalues of
the Lorenz system in Figure 6.

Figure 6: The eigenvalues of the autonomous reservoir system are red discs and
the eigenvalues of the Lorenz system are blue crosses. If an embedding is
achieved at the fixed point $m^{*}$, then after training we expect the
Jacobian $J^{\prime}_{f(m^{*})}$ to have 3 linearly independent eigenvectors
that lie in the tangent space of $f_{(\omega,\mathcal{V},F)}(M)$ with
associated eigenvalues that approximate those of the Lorenz system. We can see
3 eigenvalues of $J^{\prime}_{f(m^{*})}$ clearly approximating the 3
eigenvalues of the Lorenz system, suggesting the experiment was a success.

##  7 Generalised Synchronisations with Noisy Observations

Our analysis has so far has made the rather ideal assumption that the
observations are unperturbed by any noise. Reservoir computers appear to be
fairly resilient to noise in time series forecasting problems [29], which
suggests that under noisy observations an (injective) GS is preserved in some
form. With this inspiration we will assume that the observations are perturbed
by Gaussian white noise with an amplitude $\sigma\in C^{0}(M,\mathbb{R}^{+})$.
The reservoir states associated to a linear reservoir system are no longer
deterministic, but instead an $\mathbb{R}^{N}$-valued stochastic process
$X_{t}$ that satisfies the SDE

$\displaystyle
dX_{t}=-AX_{t}dt+C\omega\phi^{t}(m)dt+C\sigma\phi^{t}(m)dW_{t}.$ (16)

We will show in the following theorem, that for any initial reservoir state
$X_{0}$ and point $m\in M$ of the source system, the reservoir states converge
to a the image of a deterministic GS $f_{(\omega,\mathcal{V},F)}$ perturbed by
an error process with explicit form.

######  Theorem 7.1.

Suppose that $\omega,\sigma\in C^{0}(M,\mathbb{R})$ are bounded. Then for any
$m\in M$ and initial state $X_{0}\in\mathbb{R}^{N}$ the stochastic reservoir
system (16) admits the unique solution $X_{t}$ which converges in distribution
to

$\displaystyle
X_{t}=f_{(\omega,\mathcal{V},F)}\phi^{t}(m)+\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)\
dW_{\tau}$

as $t\to\infty$.

###### Proof.

We can solve SDE (16) exactly. First, let

$\displaystyle Z_{t}=e^{At}X_{t}$

then

$\displaystyle dZ_{t}$ $\displaystyle=e^{At}dX_{t}+e^{At}AX_{t}dt$
$\displaystyle=e^{At}[-AX_{t}dt+C\omega\phi^{t}(m)dt+C\sigma\phi^{t}(m)dW_{t}]+e^{At}AX_{t}dt$
$\displaystyle=e^{At}C\omega\phi^{t}(m)dt+e^{At}C\sigma\phi^{t}(m)dW_{t}$

hence

$\displaystyle
Z_{t}-Z_{-s}=\int_{-s}^{t}e^{Au}C\omega\phi^{u}(m)du+\int_{-s}^{t}e^{Au}C\sigma\phi^{u}(m)dW_{u}.$

Now we substitute $X_{t}$ back in

$\displaystyle
e^{At}X_{t}-e^{-As}X_{-s}=\int_{-s}^{0}e^{Au}C\omega\phi^{u}(m)du+\int_{-s}^{t}e^{Au}C\sigma\phi^{u}(m)dW_{u}$

so

$\displaystyle
X_{t}-e^{-A(t+s)}X_{-s}=\int_{-s}^{t}e^{A(u-t)}C\omega\phi^{u}(m)du+\int_{-s}^{t}e^{A(u-t)}C\sigma\phi^{u}(m)dW_{u}$

Then making the substitution $\tau=t-u$

$\displaystyle
X_{t}-e^{-A(t+s)}X_{-s}=\int_{0}^{t+s}e^{-A\tau}C\omega\phi^{-\tau+t}(m)d\tau+\int_{0}^{t+s}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)dW_{\tau}$
(17)

Now choosing $t=0$

$\displaystyle
X_{0}-e^{-As}X_{-s}=\int_{0}^{s}e^{-A\tau}C\omega\phi^{-\tau}(m)d\tau+\int_{0}^{s}e^{-\alpha\tau}C\sigma\phi^{-\tau}(m)dW_{\tau}$

so

$\displaystyle
e^{-As}X_{-s}=X_{0}-\int_{0}^{s}e^{-A\tau}C\omega\phi^{-\tau}(m)d\tau-\int_{0}^{s}e^{-A\tau}C\sigma\phi^{-\tau}(m)dW_{\tau}$

can now be substituted into equation (17)

$\displaystyle X_{t}$
$\displaystyle=\int_{0}^{t+s}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)d\tau+\int_{0}^{t+s}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)dW_{\tau}$
$\displaystyle\qquad+e^{-At}\bigg{[}X_{0}-\int_{0}^{s}e^{-A\tau}C\omega\phi^{-\tau}(m)d\tau-\int_{0}^{s}e^{-A\tau}C\sigma\phi^{-\tau}(m)dW_{\tau}\bigg{]}.$

Now let $s\to\infty$

$\displaystyle X_{t}$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)d\tau+\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)dW_{\tau}$
$\displaystyle\qquad+e^{-At}\bigg{[}X_{0}-\int_{0}^{\infty}e^{-A\tau}C\omega\phi^{-\tau}(m)d\tau-\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau}(m)dW_{\tau}\bigg{]}.$

The process $X_{t}$ converges in distribution to

$\displaystyle X_{t}$
$\displaystyle=\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)d\tau+\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)dW_{\tau}$
$\displaystyle=f_{\omega,\mathcal{V},F}\phi^{t}(m)+\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)dW_{\tau}$

as $t\to\infty$. ∎

The stochastic reservoir states converge to the image of a GS, perturbed by
the stochastic error term

$\displaystyle\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)dW_{\tau}.$

In the special case that $\sigma(m)=\sigma>0$ is a constant function the
stochastic error term

$\displaystyle\int_{0}^{\infty}e^{-A\tau}C\sigma\phi^{-\tau+t}(m)dW_{\tau}=\sigma\int_{0}^{\infty}e^{-A\tau}CdW_{\tau}=\mathcal{N}(0,\sigma
A^{-1})$

is a multivariate normal with mean $0$ and covariance matrix
$\sigma^{2}A^{-1}$ [39]. This is the stationary distribution of the error
process, which is a multivariate Ornstein-Uhlenbeck process [39]

$\displaystyle dY_{t}=-AY_{t}dt+\sigma CdW_{t}.$

The errors $Y_{t}$ given by this process are normal and identically
distributed, but they are not independent. This explicit form of the error may
be important when analysing the error or uncertainty on the reconstruction of
the source system.

##  8 Conclusions

In the paper, we established conditions under which a continuous time
reservoir system admits a generalised synchronisation
$f_{(\omega,\mathcal{V},F)}$, and showed that $f_{(\omega,\mathcal{V},F)}$
solves a PDE involving the Lie derivative $\mathcal{L}_{\mathcal{V}}$. We
discuss how it is possible for a reservoir system to admit several distinct
generalised synchronisations simultaneously, and how this is related to the
multi-ESP studied by [4] and the discrete time results in [10]. In the special
case of a linear reservoir system, we derived a closed form expression for
$f_{(\omega,\mathcal{V},F)}$ and establish conditions that ensure that
$f_{(\omega,\mathcal{V},F)}$ is continuously differentiable, i.e
$f_{(\omega,\mathcal{V},F)}\in C^{1}(M,\mathbb{R}^{N})$. Furthermore, it was
shown that for randomly generated weights the linear reservoir system is
almost surely an embedding of the fixed points of the source system. Having
achieved an embedding, we discuss how the central limit theorem applied to
reservoir systems makes it possible to use the reservoir states to learn a
vector field that is diffeomorphic to the source dynamics. We demonstrated
this by embedding the fixed points of the Lorenz-63 system into the reservoir
space, and training the reservoir computer to learn the eigenvalues of the
linearisation about one of the fixed points belonging to the Lorenz system.
Finally, it was shown that when the observations are perturbed by noise the GS
is preserved, but the reservoir states are perturbed by an error process. When
the errors on the observations are IID Brownian increments, the error process
on the reservoir states is an Ornstein-Uhlumbeck process.

A possible direction of future work is to establish conditions under which the
GS $f_{(\omega,\mathcal{V},F)}\in C^{1}(M,\mathbb{R}^{N})$ is a global
embedding. Very similar results have been proved in discrete time by [11],
using the same tools and methods as the celebrated Takens embedding theorem.
It has been conjectured [14] that for a larger class of nonlinear discrete
time reservoir systems the associated GS is an embedding, and this may also
hold in continuous time.

There may also be other ways to extend the results for linear systems to
nonlinear systems. The best empirical results in forecasting and
classification problems are typically achieved with a nonlinear reservoir map
$F$, but the mathematical analysis is much more challenging. Furthermore, real
data is often better modelled as a noisy trajectory of a deterministic source
system, or a trajectory of a stochastic differential equation (SDE), which
further complicates the analysis. A study of reservoir systems in these
contexts may be shed more light on the performance of reservoir computing in
practical applications.

## Appendix

######  Lemma 8.1.

Let $x^{r}_{m}(t)$ denote the solution of ODE (3) originating from the initial
point $r\in V\subset\mathbb{R}^{N}$. Let $0\leq s\leq\sigma$ Then

$\displaystyle
x^{r}_{\phi^{-\sigma}(m)}(\sigma)=x^{x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)}_{\phi^{-s}(m)}(s).$

###### Proof.

Consider ODE (3) initialised from the point $r\in V\subset\mathbb{R}^{N}$ in
the reservoir space and $\phi^{-\sigma}(m)$ in the source system. Suppose we
integrate the trajectory of (3) under these initial conditions over a time
interval of length $\sigma-s$. We label the solution at this time point:

$\displaystyle x^{r}_{\phi^{-\sigma}(m)}(\sigma-s),$

and note the value of the source system is $\phi^{-s}(m)$. Now starting at the
initial points $x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)$ and $\phi^{-s}(m)$ we
integrate the trajectory forward by a time interval of length $s$. We label
the solution:

$\displaystyle x^{x^{r}_{\phi^{-\sigma}(m)}(\sigma-s)}_{\phi^{-s}(m)}(s),$

and note that this solution was achieved by integrating (3) from initial
points $r$ and $\phi^{-\sigma}$ over a time interval $\sigma-s+s=\sigma$ which
we could equivalently label

$\displaystyle x^{r}_{\phi^{-\sigma}(m)}(\sigma).$

∎

######  Lemma 8.2.

Let $A$ be a real $N\times N$ symmetric positive definite matrix and $C$ a
random vector in $\mathbb{R}^{N}$ and assume the entries of $C$, and upper
triangular entries of $A$ are stochastically independent and each of them is
non-singular. Let $\lambda_{1}\ldots\lambda_{q}\in\mathbb{C}$ be distinct
complex numbers. Then the vectors

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

are linearly independent almost surely.

###### Proof.

It follows from Lemma 8.6 that none of the complex numbers
$\lambda_{1},\ldots,\lambda_{q}$ lie in the spectrum of $-A$ almost surely. On
this event it follows that for each $j=1,\ldots,q$ the matrix
$(A+\lambda_{j}\mathbb{I})$ is invertible, hence the product

$\displaystyle\prod_{i=1}^{q}(A+\lambda_{i}\mathbb{I})$

is invertible. Now it follows that the vectors

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

are linearly independent if and only if the vectors

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

are linearly independent. All that remains is therefore to show that the
vectors

$\displaystyle\\{p_{j}(A)C\\}_{j=1,\ldots,q}$

are linearly independent, where $p_{1},\ldots p_{q}$ are the polynomials

$\displaystyle
p_{j}(X)=\prod_{i=1}^{q}(X+\lambda_{i}\mathbb{I})(X+\lambda_{j}\mathbb{I})^{-1}=\prod_{i\neq
j}^{q}(X+\lambda_{i}\mathbb{I}).$

By lemma 8.7 is suffices to show that the polynomials $p_{1},\ldots,p_{q}$ are
linearly independent. To show this, suppose that for scalars $c_{1}\ldots
c_{q}$ we have

$\displaystyle 0=\sum_{i=1}^{q}c_{i}p_{i}(X)$

then for $X=-\lambda_{k}\mathbb{I}$

$\displaystyle
0=\sum_{i=1}^{q}c_{i}p_{i}(-\lambda_{k}\mathbb{I})=c_{k}\prod_{i\neq
k}^{q}(\lambda_{k}\mathbb{I}-\lambda_{i}\mathbb{I})=\sum_{i=1}^{q}c_{i}p_{i}(-\lambda_{k})=c_{k}\prod_{i\neq
k}^{q}(\lambda_{k}-\lambda_{i})\mathbb{I}.$

Now the complex numbers $\lambda_{1},\ldots,\lambda_{q}$ are distinct by
assumption so $c_{k}=0$. ∎

######  Lemma 8.3.

Let $M$ be a topological space and suppose
$g:C^{0}(M\times\mathbb{R},\mathbb{R}^{N})$ satisfies

$\displaystyle\int_{0}^{\infty}\sup_{m\in M}\lVert g(m,\tau)\rVert\
d\tau<\infty.$

Then

$\displaystyle h(m)=\int_{0}^{\infty}g(m,\tau)\ d\tau$

is continuous on $M$.

###### Proof.

For each $n\in\mathbb{N}$ let

$\displaystyle\varphi_{n}(m):=\int_{n}^{n+1}g(m,\tau)\ d\tau.$

Now

$\displaystyle\sum_{n=0}^{\infty}\sup_{m\in M}\lVert\varphi_{n}(m)\rVert$
$\displaystyle=\sum_{n=0}^{\infty}\sup_{m\in
M}\bigg{\rVert}\int_{n}^{n+1}g(m,\tau)\ d\tau\bigg{\rVert}$
$\displaystyle\leq\sum_{n=0}^{\infty}\int_{n}^{n+1}\sup_{m\in M}\rVert
g(m,\tau)\rVert\ d\tau$ $\displaystyle=\int_{0}^{\infty}\sup_{m\in M}\lVert
g(m,\tau)\rVert\ d\tau<\infty$

so the series

$\displaystyle\sum_{n=0}^{\infty}\varphi_{n}(m)=\sum_{n=0}^{\infty}\int_{n}^{n+1}g(m,\tau)\
d\tau=\int_{0}^{\infty}g(m,\tau)\ d\tau=h(m)$

converges absolutely and uniformly on $M$ by the Weierstrass $M$-test. The
uniform limit of continuous functions is continuous by the uniform limit
theorem. ∎

######  Lemma 8.4.

For any symmetric positive definite matrix $X$ let $\sigma[X]_{\max}$ and
$\sigma[X]_{\min}$ denote the largest and smallest eigenvalues of $X$. then
for any $t\geq 0$ and symmetric positive definite matrix $A$

$\displaystyle\sigma[e^{-At}]_{\max}=e^{\sigma[A]_{\min}t}.$

###### Proof.

Denote the spectrum of $A$ and $e^{-At}$ with $\sigma[A]=\\{\mu_{i}\\}$ and
$\sigma[e^{-At}]=\\{e^{-\mu_{i}t}\\}$ respectively. Now for fixed $t\geq 0$
the function $e^{-xt}$ is decreasing with respect to $x>0$, so the largest
eigenvalue in $\\{e^{-\mu_{i}t}\\}$ is $e^{\sigma[A]_{\min}t}$ ∎

######  Lemma 8.5.

(Lemma 6.1 [11]) Let $A$ and $B$ be two square matrices of the same size such
that $\det(A)=0$ and $\det(B)\neq 0$. Then there exists an $\epsilon>0$ such
that

$\displaystyle\det(A-\epsilon B)\neq 0$

######  Lemma 8.6.

(Lemma 6.2 [11]) Let $X_{1},\ldots,X_{n}$ be independent real random variables
each of which is non singular, and let $p$ be a polynomial in $n$ complex
variables, that is not identically zero. Then

$\displaystyle\mathbb{P}[p(X_{1},\ldots X_{n})=0]=0.$

######  Lemma 8.7.

(Lemma 6.3 [11]) Let $A$ be a real symmetric positive definite $N\times N$
matrix and $C$ a random vector in $\mathbb{R}^{N}$ and assume the entries of
$C$ and the upper triangular entries of $A$ are stochastically independent and
each of them is non-singular. Moreover, let
$p_{1},\ldots,p_{n}\in\mathbb{C}[x]$ be linearly independent polynomials in
one variable of degree at most $n-1$. Then

$\displaystyle\mathbb{P}(\text{det}[p_{1}(A)C,p_{2}(A)C,\ldots,p_{n}(A)C]=0)=0,$

i.e the vectors $p_{1}(A)C,p_{2}(A)C,\ldots,p_{n}(A)C$ are linearly
independent almost surely.

###### Proof.

Lemma 6.3 [11] is stated for real matrices $A$ that are not in general
symmetric positive definite, but the same proof goes through in the symmetric
positive definite case. ∎

## References

  * [1]
  * Bhovad and Li [2021] Bhovad, P. and Li, S. [2021], ‘Physical reservoir computing with origami and its application to robotic crawling’, Scientific Reports 11(1), 13002.
https://doi.org/10.1038/s41598-021-92257-1

  * Boccaletti et al. [2002] Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. and Zhou, C. [2002], ‘The synchronization of chaotic systems’, Physics Reports 366(1), 1–101.
https://www.sciencedirect.com/science/article/pii/S0370157302001370

  * Ceni et al. [2020] Ceni, A., Ashwin, P., Livi, L. and Postlethwaite, C. [2020], ‘The echo index and multistability in input-driven recurrent neural networks’, Physica D: Nonlinear Phenomena 412, 132609.
http://www.sciencedirect.com/science/article/pii/S0167278920300427

  * Eroglu et al. [2017] Eroglu, D., Lamb, J. S. W. and Pereira, T. [2017], ‘Synchronisation of chaos and its applications’, Contemporary Physics 58(3), 207–243.
https://doi.org/10.1080/00107514.2017.1345844

  * Gauthier et al. [2021] Gauthier, D. J., Bollt, E., Griffith, A. and Barbosa, W. A. S. [2021], ‘Next generation reservoir computing’, Nature Communications 12(1), 5564.
https://doi.org/10.1038/s41467-021-25801-2

  * Gonon [2021] Gonon, L. [2021], ‘Random feature neural networks learn black-scholes type pdes without curse of dimensionality’, arXiv:2106.08900 .
  * Gonon et al. [2020] Gonon, L., Grigoryeva, L. and Ortega, J.-P. [2020], ‘Approximation bounds for random neural networks and reservoir systems’, arXiv:2002.05933 .
  * Gonon and Ortega [2021] Gonon, L. and Ortega, J.-P. [2021], ‘Fading memory echo state networks are universal’, Neural Networks 138, 10–13.
https://www.sciencedirect.com/science/article/pii/S0893608021000332

  * Grigoryeva et al. [2021a] Grigoryeva, L., Hart, A. and Ortega, J.-P. [2021a], ‘Chaos on compact manifolds: Differentiable synchronizations beyond the takens theorem’, Phys. Rev. E 103, 062204.
https://link.aps.org/doi/10.1103/PhysRevE.103.062204

  * Grigoryeva et al. [2021b] Grigoryeva, L., Hart, A. and Ortega, J.-P. [2021b], ‘Learning strange attractors with reservoir systems’, arXiv:2108.05024 .
  * Grigoryeva and Ortega [2018] Grigoryeva, L. and Ortega, J. [2018], ‘Echo state networks are universal’, Neural Networks 108, 495 – 508.
  * Hart et al. [2021] Hart, A. G., Hook, J. L. and 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.
  * Hart et al. [2020] Hart, A., Hook, J. and Dawes, J. [2020], ‘Embedding and approximation theorems for echo state networks’, Neural Networks 128, 234–247.
https://www.sciencedirect.com/science/article/pii/S0893608020301830

  * Hirsch [1976] Hirsch, M. W. [1976], Differential Topology, Springer-Verlag.
  * Huke [2006] Huke, J. P. [2006], ‘Embedding nonlinear dynamical systems: A guide to Takens’ theorem’.
http://eprints.maths.manchester.ac.uk/175/1/embed.pdf

  * Inubushi and Yoshimura [2017] Inubushi, M. and Yoshimura, K. [2017], ‘Reservoir computing beyond memory-nonlinearity trade-off’, Scientific Reports 7(1), 10199.
  * Jaeger [2001] Jaeger, H. [2001], ‘The “echo state” approach to analysing and training recurrent neural networks’.
https://www.ai.rug.nl/minds/uploads/EchoStatesTechRep.pdf

  * Jaeger [2002] Jaeger, H. [2002], ‘Tutorial on training recurrent neural networks, covering bppt, rtrl, ekf and the echo state network approach’, 5.
  * Jaeger and Haas [2004] Jaeger, H. and Haas, H. [2004], ‘Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication’, Science 304(5667), 78–80.
  * Jaeger et al. [2007] Jaeger, H., Lukoševičius, M., Popovici, D. and Siewert, U. [2007], ‘Optimization and applications of echo state networks with leaky- integrator neurons’, Neural Networks 20(3), 335 – 352.  Echo State Networks and Liquid State Machines.
  * Kocarev and Parlitz [1996] Kocarev, L. and Parlitz, U. [1996], ‘Generalized synchronization, predictability, and equivalence of unidirectionally coupled dynamical systems’, Phys. Rev. Lett. 76, 1816–1819.
  * Lorenz [1963] Lorenz, E. N. [1963], ‘Deterministic nonperiodic flow’, Journal of the Atmospheric Sciences 20(2), 130–141.
  * Lun et al. [2019] Lun, S.-x., Hu, H.-f. and Yao, X.-s. [2019], ‘The modified sufficient conditions for echo state property and parameter optimization of leaky integrator echo state network’, Applied Soft Computing 77, 750–760.
  * Lun et al. [2015] Lun, S.-X., Yao, X.-S., Qi, H.-Y. and Hu, H.-F. [2015], ‘A novel model of leaky integrator echo state network for time-series prediction’, Neurocomputing 159, 58–66.
  * Maass et al. [2002] Maass, W., Natschläger, T. and Markram, H. [2002], ‘Real-time computing without stable states: A new framework for neural computation based on perturbations’, Neural Computation 14(11), 2531–2560.
  * Manjunath [2022] Manjunath, G. [2022], ‘Embedding information onto a dynamical system’, Nonlinearity 35(3), 1131–1151.
https://doi.org/10.1088/1361-6544/ac4817

  * Mhaskar [1996] Mhaskar, H. N. [1996], ‘Neural networks for optimal approximation of smooth and analytic functions’, Neural Computation 8, 164–177.
  * Park et al. [2022] Park, Y., Gajamannage, K., Jayathilake, D. I. and Bollt, E. M. [2022], ‘Recurrent neural networks for dynamical systems: Applications to ordinary differential equations, collective motion, and hydrological modeling’, arXiv:2202.07022 .
  * Pecora et al. [1997] Pecora, L. M., Carroll, T. L., Johnson, G. A., Mar, D. J. and Heagy, J. F. [1997], ‘Fundamentals of synchronization in chaotic systems, concepts, and applications’, Chaos: An Interdisciplinary Journal of Nonlinear Science 7(4), 520–543.
https://doi.org/10.1063/1.166278

  * Plöger et al. [2004] Plöger, P. G., Arghir, A., Günther, T. and Hosseiny, R. [2004], Echo state networks for mobile robot modeling and control, in D. Polani, B. Browning, A. Bonarini and K. Yoshida, eds, ‘RoboCup 2003: Robot Soccer World Cup VII’, Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 157–168.
  * Poggio et al. [2017] Poggio, T., Mhaskar, H., Rosasco, L., Miranda, B. and Liao, Q. [2017], ‘Why and when can deep-but not shallow-networks avoid the curse of dimensionality: A review’, International Journal of Automation and Computing 14(5), 503–519.
https://doi.org/10.1007/s11633-017-1054-2

  * Rulkov et al. [1995] Rulkov, N. F., Sushchik, M. M., Tsimring, L. S. and Abarbanel, H. D. I. [1995], ‘Generalized synchronization of chaos in directionally coupled chaotic systems’, Phys. Rev. E 51, 980–994.
https://link.aps.org/doi/10.1103/PhysRevE.51.980

  * Sauer et al. [1991] Sauer, T., Yorke, J. A. and Casdagli, M. [1991], ‘Embedology’, Journal of Statistical Physics 65(3), 579–616.
  * Shi and Han [2007] Shi, Z. and Han, M. [2007], ‘Support vector echo-state machine for chaotic time-series prediction’, IEEE Transactions on Neural Networks 18(2), 359–372.
  * Stark [1999] Stark, J. [1999], ‘Regularity of invariant graphs for forced systems’, Ergodic Theory and Dynamical Systems 19(1), 155–199.
  * Takens [1981] Takens, F. [1981], ‘Detecting strange attractors in turbulence’, Lecture Notes in Mathematics, Berlin Springer Verlag 898, 366.
  * Tanaka et al. [2019] Tanaka, G., Yamane, T., Héroux, J. B., Nakane, R., Kanazawa, N., Takeda, S., Numata, H., Nakano, D. and Hirose, A. [2019], ‘Recent advances in physical reservoir computing: A review’, Neural Networks 115, 100 – 123.
  * Vatiwutipong and Phewchean [2019] Vatiwutipong, P. and Phewchean, N. [2019], ‘Alternative way to derive the distribution of the multivariate ornstein–uhlenbeck process’, Advances in Difference Equations 2019(1), 276.
https://doi.org/10.1186/s13662-019-2214-1

  * Verzelli et al. [2020] Verzelli, P., Alippi, C. and Livi, L. [2020], ‘Learn to synchronize, synchronize to learn’, arXiv:2010.02860 .
  * Xi et al. [2005] Xi, J., Shi, Z. and Han, M. [2005], Analyzing the state space property of echo state networks for chaotic system prediction, in ‘Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005.’, Vol. 3, pp. 1412–1417 vol. 3.