# Echo State Networks trained by Tikhonov least squares are $L^{2}(\mu)$
approximators of ergodic dynamical systems

Allen G. Hart<EMAIL_ADDRESS>
University of Bath, UK  James L. Hook<EMAIL_ADDRESS>
University of Bath, UK  Jonathan H. P. Dawes<EMAIL_ADDRESS>
University of Bath, UK

###### Abstract

Echo State Networks (ESNs) are a class of single-layer recurrent neural
networks with randomly generated internal weights, and a single layer of
tuneable outer weights, which are usually trained by regularised linear least
squares regression. Remarkably, ESNs still enjoy the universal approximation
property despite the training procedure being entirely linear. In this paper,
we prove that an ESN trained on a sequence of observations from an ergodic
dynamical system (with invariant measure $\mu$) using Tikhonov least squares
regression against a set of targets, will approximate the target function in
the $L^{2}(\mu)$ norm. In the special case that the targets are future
observations, the ESN is learning the next step map, which allows time series
forecasting. We demonstrate the theory numerically by training an ESN using
Tikhonov least squares on a sequence of scalar observations of the Lorenz
system.

###### keywords:

Reservoir computing; liquid state machine; time series analysis; Lorenz
equations; dynamical system; delay embedding; Ergodic theory; recurrent neural
networks.

††journal: Physica D

##  1 Introduction

Echo state networks (ESNs) are a class of single layer recurrent neural
networks introduced at the turn on the millennium independently by Jaeger
[2001] and Maass et al. [2002]. These relatively simple neural networks have
been used to solve a range of machine problems where the input data is a time
series, including speech recognition [Skowronski and Harris, 2007], learning
the rules of grammar [Tong et al., 2007], financial time series prediction
[Ilies et al., 2007], [Lin et al., 2009], short term traffic forecasting [Ser
et al., 2020], placing UAV base stations [Peng et al., 2019] and learning
about the behaviour of seals [Ser et al., 2020]. ESNs are also a plausible
model for the information processing of biological neurons [Gürel and Egert,
2010]. In this paper, we will present just enough definitions and theory to
make sense of our results, but encourage the interested reader to read the
recent review paper by Tanaka et al. [2019] who cover recent developments and
open questions in the field of _reservoir computing_ , a field of which ESN
comprise a subset. The ESN is defined by the recursion relation

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

where the $x_{k}$ are $T$ dimensional state vectors,
$\sigma:\mathbb{R}^{T}\to\mathbb{R}^{T}$ is the activation function, $A$ is
the $T\times T$ _reservoir matrix_ , representing the connection weights
between neurons, $C$ is the $T\times d$ _input matrix_ connecting the
$d$-dimensional inputs $z_{k}$ to the reservoir matrix $A$, and
$b\in\mathbb{R}^{T}$ is a bias vector. The reservoir matrix $A$, input matrix
$C$ and bias vector $b$ are initialised randomly and remain unchanged. The ESN
can be trained to approximate a sequence of target scalars $u_{k}$ by solving
the regularised linear squares problem

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

where $\lambda>0$ is the Tikhonov regularisation parameter. If the target
scalars $u_{k}$ are equal to the observations $z_{k}$, then the ESN is being
trained to predict the future. To see this, we can set up a sequence of
scalars $v_{k}$ defined by the recurrent relation

$\displaystyle v_{k+1}=W^{\top}s_{k}$ $\displaystyle
s_{k+1}=\sigma(As_{k}+Cv_{k+1}+b)$ (1)

and we then hope that $v_{k}\approx u_{k}$ for sufficiently many future values
of $k$. We can view $s_{k}$ as the state of a discrete time autonomous
dynamical system which we will call the ESN autonomous phase. In this paper,
we will suppose $z_{k}$ are a sequence of sequential observations from an
ergodic dynamical system, with invariant measure $\mu$. We will go on to prove
that ESNs trained by least squares can approximate arbitrary target functions
(including one that returns future observations) of the ergodic dynamical
system in the $L^{2}(\mu)$ norm. This theorem is closely related to recent
work by Verzelli et al. [2020] discussing the connection between ergodic
dynamical systems and feasible learning. The result also explains the
remarkable success of ESNs trained on dynamical systems explored numerically
by, for example, Jaeger [2001], Xi et al. [2005], Schrauwen et al. [2007], Shi
and Han [2007], Yong Song et al. [2010], Yildiz et al. [2012], Pathak et al.
[2017], Løkse et al. [2017], Yeo [2019], Chattopadhyay et al. [2019], Vlachas
et al. [2019], Hart et al. [2019].

The remainder of the paper is organised as follows. In section 2 we define an
ergodic dynamical system and present Birkhoff’s ergodic theorem. Next, in
section 3, we introduce the major result of this paper (Theorem 3.6), stating
that an ESN trained on a sequence of observations from an ergodic dynamical
system using Tikhonov least squares will $L^{2}(\mu)$ approximate an arbitrary
target function. This arbitrary target function could be the _next step map_
used for forecast the future of the time series. Furthermore, we discuss the
_central limit theorem for ergodic dynamical systems_ in connection with the
number of data points that are required for a good approximation.

In section 4 we present the work of Luzzatto et al. [2005] culminating in a
proof that the Lorenz attractor is stably mixing, hence ergodic - explaining
the success of so many authors using an ESN to forecast the trajectory of the
Lorenz system.

In section 5 we numerically simulate a trajectory of the Lorenz system. We
observed the $x$-component of the system (which we called $\xi$ to avoid
notational clash) and assigned the $z$ components (which we denote $\zeta$) as
targets. We explored how the approximation of the targets $\zeta_{k}$ given
the observations $\xi_{k}$ improved as the number of data points
$(\xi_{k},\zeta_{k})$ grew. Finally, in section 6 we summarise the paper and
discuss ideas for future work.

##  2 Preliminaries on Ergodic Theory

We require that the underlying dynamical system is ergodic so that minimising
the mean square differences between observations and targets does not create a
bias toward areas with lots of training data. The ergodicity ensures that that
training data generated from a trajectory initialised at almost any point
$m_{0}\in M$ will represent all dynamics on $M$. To make this formal, we will
introduce the definition of ergodicty and the celebrated ergodic theorem.

######  Definition 2.1.

(Generic Point) Suppose $\phi:M\to M$ is a measure preserving map with respect
to the measure space $(M,\Sigma,\mu)$. Then $m_{0}\in M$ is called a generic
point if the orbit of $m_{0}$ is uniformly distributed over $M$ according to
the measure $\mu$.

######  Proposition 2.2.

Suppose $\phi:M\to M$ is a measure preserving map with respect to the
probability space $(M,\Sigma,\mu)$ and $s\in L^{1}(\mu)$. Suppose $m_{0}$ is a
generic point in $M$ then

$\displaystyle\lim_{\ell\to\infty}\frac{1}{\ell}\sum_{k=0}^{\ell-1}s\circ\phi^{k}(m_{0})=\int_{M}s\
d\mu.$

######  Definition 2.3.

(Ergodic) Let $\phi:M\to M$ be a measure preserving transformation on the
probability space $(M,\Sigma,\mu)$. Then $\phi$ is ergodic if for every
$\sigma\in\Sigma$ with $\phi^{-1}(\sigma)=\sigma$ either $\mu(\sigma)=0$ or
$\mu(\sigma)=1$.

######  Theorem 2.4.

(Ergodic Theorem [Birkhoff, 1931]) Suppose $\phi:M\to M$ is ergodic with
respect to the probability space $(M,\Sigma,\mu)$ and $s\in L^{1}(\mu)$. Then
$\mu$-almost all $m_{0}\in M$ are generic hence for $\mu$-almost all $m_{0}\in
M$

$\displaystyle\lim_{\ell\to\infty}\frac{1}{\ell}\sum_{k=0}^{\ell-1}s\circ\phi^{k}(m_{0})=\int_{M}s\
d\mu.$ (2)

The left hand side of (2) is called the _time average_ taken from initial
point point $m_{0}\in M$, and the right hand side called the _space average_.
The ergodic theorem then states that the time average taken from almost all
initial points equals the space average.

##  3 A Training Theorem for Echo State Networks

###  3.1 Preliminaries

Suppose we have have an ergodic dynamical system $\phi:M\to M$, and we can
observe the dynamics via an observation map $g:M\to\mathbb{R}^{T}$ and target
map $u:M\to\mathbb{R}$. A trajectory originating from a generic point
$m_{0}\in M$ will ergodically explore the space $M$ and yield a sequence of
observations $g\circ\phi^{k}(m_{0})$ and targets $u\circ\phi^{k}(m_{0})$ for
$k=0,1,2,...,\ell$.

Suppose we compute the vectors $W_{\ell}\in\mathbb{R}^{T}$ minimising the
regularised least squares difference between the mapping of the observations
$W^{\top}g\circ\phi^{k}(m_{0})$ and the targets $u\circ\phi^{k}(m_{0})$. We
prove in the next lemma that as the number of data points $\ell$ grows large,
the least squares solution $W_{\ell}$ minimises the ergodic average difference
between the mapping on the observations $W^{\top}g\circ\phi^{k}(m_{0})$ and
the targets $u\circ\phi^{k}(m_{0})$.

######  Lemma 3.1.

Let $(M,\Sigma)$ be a measurable space, and suppose that $\phi:M\to M$ is
ergodic with invariant measure $\mu$. Let $m_{0}$ be a generic point in $M$.
Let $g\in L^{2}(\mu)(M,\mathbb{R}^{T})$ be an observation function and suppose
that $u\in L^{2}(\mu)(M,\mathbb{R})$ is a target function we wish to
approximate.

Let $\lambda>0$. Define the sequence $(W_{\ell})_{\ell\in\mathbb{N}}$ such
that, for each $\ell\in\mathbb{N}$, the vector $W_{\ell}\in\mathbb{R}^{T}$ is
the unique minimiser of the regularised least squares difference

$\displaystyle\frac{1}{\ell}\bigg{(}\sum_{k=0}^{\ell-1}\rVert
W^{\top}g\circ\phi^{k}(m_{0})-u\circ\phi^{k}(m_{0})\lVert^{2}+\lambda\lVert
W\rVert^{2}\bigg{)}.$

Then, the sequence $(W_{\ell})_{\ell\in\mathbb{N}}$ converges to

$\displaystyle W_{\infty}$ $\displaystyle=\bigg{(}\int_{M}g(m)g(m)^{\top}\
d\mu(m)+\lambda I\bigg{)}^{-1}$ $\displaystyle\times\int_{M}u(m)g(m)\ d\mu(m)$

which is the unique minimiser of

$\displaystyle\lVert W^{\top}g-u\rVert_{L^{2}(\mu)}^{2}+\lambda\lVert
W\rVert^{2}.$

###### Proof.

Consider the map $\Psi:\mathbb{R}^{T}\to\mathbb{R}$ defined

$\displaystyle\Psi(W)$ $\displaystyle=\lVert
W^{\top}g-u\rVert_{L^{2}(\mu)}^{2}+\lambda\lVert W\rVert^{2}$
$\displaystyle=\int_{M}\lVert W^{\top}g(m)-u(m)\rVert^{2}\
d\mu(m)+\lambda\lVert W\rVert^{2}.$

The minimiser of $\Psi$ satisfies $D\Psi=0$ where $D$ is the derivative
operator, so we consider

$\displaystyle 0$ $\displaystyle=(D\Psi)(W)$
$\displaystyle=D\bigg{(}\int_{M}\lVert W^{\top}g(m)-u(m)\rVert^{2}\
d\mu(m)+\lambda\lVert W\rVert^{2}\bigg{)}$ $\displaystyle=\int_{M}D\lVert
W^{\top}g(m)-u(m)\rVert^{2}\ d\mu(m)+\lambda D\lVert W\rVert^{2}$
$\displaystyle=\int_{M}2(W^{\top}g(m)-u(m))g(m)^{\top}\ d\mu(m)+2\lambda
W^{\top}$ $\displaystyle=\int_{M}(W^{\top}g(m)-u(m))g(m)^{\top}\
d\mu(m)+\lambda W^{\top}$ $\displaystyle=W^{\top}\int_{M}g(m)g(m)^{\top}\
d\mu(m)-\int_{M}u(m)g(m)^{\top}\ d\mu(m)$ $\displaystyle+\lambda W^{\top}I$
$\displaystyle=W^{\top}\bigg{(}\int_{M}g(m)g(m)^{\top}\ d\mu(m)+\lambda
I\bigg{)}$ $\displaystyle-\int_{M}u(m)g(m)^{\top}\ d\mu(m),$

which upon rearrangement yields

$\displaystyle W$ $\displaystyle=\bigg{(}\int_{M}g(m)g(m)^{\top}\
d\mu(m)+\lambda I\bigg{)}^{-1}$ $\displaystyle\times\int_{M}u(m)g(m)\
d\mu(m).$

Since this is the unique solution to $0=D\Psi(W)$, this stationary point is
unique, and we will denote it $W_{\infty}$. We can see it is a minimum because
the Hessian $H\Psi$ is positive definite. Next, define the map

$\displaystyle\Phi:\\{y\in C^{1}(\mathbb{R}^{T},\mathbb{R})\ |\ y\text{ has a
unique minimum}\\}\to\mathbb{R}^{T}$

as the mapping on the $C^{1}$ functions with a unique minumum that returns their unique minimum. We can see that $\Phi$ is continuous with respect to the $C^{1}$ topology and standard topology on $\mathbb{R}$ respectively. We consider the family of functions $y_{\ell}\in\\{y\in C^{1}(\mathbb{R}^{T},\mathbb{R})$ | $y$ has a unique minimum}

$\displaystyle y_{\ell}(W)=\frac{1}{\ell}\bigg{(}\sum_{k=0}^{\ell-1}\rVert
W^{\top}g\circ\phi^{k}(m_{0})-u\circ\phi^{k}(m_{0})\lVert^{2}+\lambda\lVert
W\rVert^{2}\bigg{)},$

so that by definition $W_{\ell}=\Phi(y_{\ell}(W))$ and hence

$\displaystyle\lim_{\ell\to\infty}W_{\ell}$
$\displaystyle=\lim_{\ell\to\infty}\Phi(y_{\ell}(W))$
$\displaystyle=\Phi\bigg{(}\lim_{\ell\to\infty}y_{\ell}(W)\bigg{)}$
$\displaystyle=\Phi\bigg{(}\lVert
W^{\top}g-u\rVert_{L^{2}(\mu)}^{2}+\lambda\lVert W\rVert^{2}\bigg{)}$
$\displaystyle=W_{\infty}.$

where we have used, respectively, continuity of $\Phi$ and the Ergodic
Theorem. ∎

###  3.2 Echo State Networks

An Echo State Network is a special case of a more general system called a
state space system, or reservoir system. These are maps of the form
$F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$, which admit an ESN as
a special case when

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

If a state space system is contracting in the state variable, i.e there exists
a $c\in[0,1)$ such that

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

and the inputs $u_{k}$ are the observations of a dynamical system i.e
$u_{k}=\omega\circ\phi^{k}(m_{0})$ then there is a continuous map $f\in
C^{0}(M,\mathbb{R}^{N})$ synchronising the dynamics of $\phi$ on $M$ to the
dynamics of the reservoir states $x_{k}$. The map $f$ is called a state
synchronisation map (SSM) and is a generalised synchronisation in the sense
described by Kocarev and Parlitz [1996]. We can guarantee that an ESN is state
contracting by bounding the 2-norm of the reservoir matrix $\lVert
A\rVert_{2}<1$. An important existence result for SSMs is the following
theorem, due to Grigoryeva et al. [2020].

######  Theorem 3.2.

[Grigoryeva et al., 2020] Let $M$ be a topological space,
$\phi\in\text{Hom}(M)$ be a dynamical system, and $\omega\in
C^{0}(M,\mathbb{R}^{d})$ an observation function. Suppose that the state space
system $F:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ is state
contracting, i.e there exists a $c\in[0,1)$ such that

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

Then there exists a unique $f\in C^{0}(M,\mathbb{R}^{N})$ called the state
synchronisation map (SSM) such that, for any $m_{0}\in M$ and
$x_{0}\in\mathbb{R}^{N}$ the sequence

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

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

In order to approximate the arbitrary dynamics of $\phi$ via the observation
function $\omega$ using state space systems, we require that the state space
maps $F$ possess some sort of universal approximation property. Thus, we will
define a class of _linear universal approximators_ with respect to an
arbitrary complete norm $\lVert\cdot\rVert$. Every class of linear universal
approximators contains maps, which after composition with another suitable
map, forms a state map.

######  Definition 3.3.

Let $\mathcal{F}$ be a sequence of maps
$\\{F_{T}\\}:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{T}$. Let
$C\subset\mathbb{R}^{N}$ and $K\subset\mathbb{R}^{d}$ be vectors and let
$\Omega(C\times K,\mathbb{R})$ be a Banach space of real valued functions on
$C\times K$, with norm denoted $\lVert\cdot\rVert_{\Omega}$. If, for any
$g\in\Omega(C\times K,\mathbb{R})$ and any $\epsilon>0$ there exists an
$T_{0}\in\mathbb{N}$ such that for any $T>T_{0}$ there exists a
$W_{*}\in\mathbb{R}^{T}$ such that

$\displaystyle\rVert W_{*}^{\top}F_{T}-g\lVert_{\Omega}<\epsilon$

then we say that $\mathcal{F}$ is a class of _linear universal approximators_
on $\Omega(C\times K,\mathbb{R}^{N}\times\mathbb{R}^{d})$.

A widely used class of linear universal approximators is the class of Echo
State Networks with randomly initialised internal weights, as shown by the
following result.

######  Theorem 3.4.

Let $\mathcal{F}$ denote the sequence of maps
$\\{F_{T}\\}:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{T}$ defined by

$\displaystyle F_{T}(x,z)=\sigma(Ax+Cz+b)$

where

  * 1.

$\sigma\in C^{1}(\mathbb{R})$ is $1$-finite (see Hornik et al. [1990] for the
definition of $\ell$-finite)

  * 2.

$A$ is a $T\times N$ random matrix, where $T>N$ and the first $N$ rows of $A$
form an $N\times N$ random submatrix with 2-norm less than 1 almost surely.
The $j^{\mathrm{th}}$ row of $A$ (where $j>N$), denoted $A_{j}$, is a random
variable with full support on $(\mathbb{R}^{N})^{\top}$

  * 3.

$C$ is a $T\times d$ random matrix with $j^{\mathrm{th}}$ row $C_{j}$, a
random variable with full support on $(\mathbb{R}^{d})^{\top}$

  * 4.

$b$ is a random $T$-vector with $j\mathrm{th}$ entry $b_{j}$, a random
variable with full support on $\mathbb{R}$.

Let $C\times K$ be an arbitrary compact subset of
$\mathbb{R}^{N}\times\mathbb{R}^{d}$. Then, almost surely, $\mathcal{F}$ is a
class of linear universal approximators on $L^{2}(C\times K,\mathbb{R})$.

###### Proof.

Fix $g\in L^{2}(C\times K,\mathbb{R})$ and $\epsilon>0$. Then for any
$\alpha\in(0,1)$, it follows from the Random Universal Approximation Theorem
[Hart et al., 2019, Theorem 2.4.5.]) that there exists a $T_{0}\in\mathbb{N}$
such that for any $T>T_{0}$, with probability at least $\alpha$,

$\displaystyle\rVert W^{\top}F_{T}-g\lVert_{L^{2}}<\epsilon,$

hence $\mathcal{F}$ is a class of linear universal approximators. Since
$\mathcal{F}$ is a class of linear universal approximators for any
$\alpha\in(0,1)$, $\mathcal{F}$ is almost surely a class of linear universal
approximators. ∎

To construct such an ESN in practice, we create a reservoir system
$F:\mathbb{R}^{T}\times\mathbb{R}^{d}\to\mathbb{R}^{T}$ by defining

$\displaystyle F(x,z)=\sigma\big{(}[A,0]x+Cz+b\big{)}$

where $[A,0]$ is the $T\times T$ matrix where the first $N$ columns form the
matrix $A$ and the remaining columns are 0. Suppose we truncate at $N$ the
state vectors $x\in\mathbb{R}^{T}$ by applying the canonical projection
$\pi:\mathbb{R}^{T}\to\mathbb{R}^{N}$, and denote the truncation
$\pi(x)=\bar{x}\in\mathbb{R}^{N}$. The dynamics of the truncated vectors
$\bar{x}$ are given by the (state contracting) state space system $\pi\circ
F_{T}:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$, which is also an
ESN as is defined by

$\displaystyle\pi\circ
F_{T}(\bar{x},z)=\sigma(\bar{A}\bar{x}+\bar{C}z+\bar{b}).$

Here, the $N\times N$ reservoir matrix $\bar{A}$ is created by truncating at
$N$ the rows and columns of $A$. The $N\times d$ input matrix $\bar{C}$ is
created by truncating at $N$ the rows of $C$. The $N$-vector $\bar{b}$ is
created by truncating at $N$ the entries of $b$. We conclude that Echo State
Networks with (appropriately chosen) randomly generated internal weights are a
class of linear universal approximators that each give rise to a state
synchronisation map.

We demanded that the $T\times T$ reservoir matrix take the form $[A,0]$,
whereas in practice, the reservoir matrix does not have this structure. We
imposed this condition to simplify the proofs, but we believe, based on
numerical evidence in the literature, that this choice of shape is not
necessary.

There is one more technical lemma we will include here before presenting the
main theorem (Theorem 3.6) of the paper. Recall that topological spaces have a
natural Borel sigma algebra and are therefore measurable spaces. On such
spaces we can integrate real valued functions. If $A$ and $B$ are homeomorphic
topological spaces, then integration on $A$ is essentially the same as
integration on $B$. We use this observation in Theorem 3.6 to move between
integration on the topological space $M$ to integration on the image $f(M)$.
This demands the highly non-trivial assumption that the SSM $f$ is a
homeomorphism. The observation is made formal in the following lemma.

######  Lemma 3.5.

(Change of variables) Let $A,B$ be homeomorphic topological spaces and suppose
$y\in\text{Hom}(A,B)$. The topologies on $A,B$ induce Borel Sigma algebras
$\mathscr{A},\mathscr{B}$ on $A,B$ respectively. Let $\mu_{A}$ be a measure on
$A$ and $\mu_{B}$ a measure on $B$ (called the pushforward measure) defined
$\mu_{B}(b)=\mu_{A}(y^{-1}(b))$ for all $b\in\mathscr{B}$. Then for any
$\mu_{B}$ measurable function $g:B\to\mathbb{R}$

$\displaystyle\int_{A}g\circ y\ d\mu_{A}=\int_{B}g\ d\mu_{B}.$

###### Proof.

This is a special case of Theorem 3.6.1 in Bogachev [2007]. ∎

###  3.3 A Training Theorem For ESNs

Before we finally plunge into the statement and proof of the main theorem, we
will describe the result in words. Suppose we have an ergodic dynamical system
$\phi:M\to M$, which we observe via the function $\omega:M\to\mathbb{R}^{d}$
and that our goal is to approximate a target function $u:M\to\mathbb{R}$.
Suppose we have at our disposal a class $\mathcal{F}$ of linear universal
approximating state maps. For example, $\mathcal{F}$ could be a collection of
arbitrarily high dimensional ESNs. Make the additional (and non trivial)
assumption that the state maps give rise to an SSM that is homeomorphic onto
its image. Suppose then that the state map $F$ is driven with observations of
a trajectory $z_{k}=\omega\circ\phi^{k}(m_{0})$ originating from a generic
point $m_{0}$. This creates a sequence of reservoir states $x_{k}$ that
satisfy

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

We also assemble a sequence of scalar targets $u\circ\phi^{k}(m_{0})$.

Suppose we use regularised least squares regression to minimise the difference
between the linear mapping on the observations $W^{\top}x_{k}$ and the targets
$u\circ\phi^{k}(m_{0})$. Then we can conclude that the ergodic average
difference between the mapping on the data and the target map $u$ can be made
smaller than the arbitrary threshold $\epsilon$. This requires that the
trajectory length $\ell$ and state map dimension $T$ are sufficiently large,
while ensuring the regularisation parameter $\lambda>0$ is sufficiently small.

We remark that a notable weakness of Theorem 3.6 is its non-constructive
natue, because the actual values for $\ell$, $T$ and $\lambda$ are not
computed in terms of $\epsilon$.

######  Theorem 3.6.

Let $M$ be a topological space, and suppose that $\phi\in\text{Hom}(M)$ is
ergodic with invariant measure $\mu$. Let $m_{0}$ be a generic point in $M$.
Let $\omega\in C^{0}(M,\mathbb{R}^{d})$ be the observation function and
suppose that $u\in L^{2}(\mu)(M,\mathbb{R})$ is a target function we wish to
approximate.

Suppose that $\mathcal{F}$ is a class of linear universal approximators on
$L^{2}(C\times K,\mathbb{R})$ on every compact
$C\subset\mathbb{R}^{N},K\subset\mathbb{R}^{d}$. Let
$(s_{T})_{T\in\mathbb{N}}:\mathbb{R}^{T}\to\mathbb{R}^{N}$ be a sequence of
maps. Suppose (for each large enough $T$) the state map $s_{T}\circ
F_{T}:\mathbb{R}^{N}\times\mathbb{R}^{d}\to\mathbb{R}^{N}$ admits an SSM
$f\in\text{Hom}(M,f(M))$. For each $T,\ell\in\mathbb{N}$, and $\lambda>0$ let
$W_{\ell}\in\mathbb{R}^{T}$ be the vector obtained by minimising the
regularised least squares difference

$\displaystyle\sum_{k=0}^{\ell}\rVert
W^{\top}F_{T}(f\circ\phi^{k-1}(m_{0}),\omega\circ\phi^{k}(m_{0}))-u\circ\phi^{k}(m_{0})\lVert^{2}$
$\displaystyle+\lambda\lVert W\rVert^{2}$ .

Then, for any $\epsilon>0$, there exists $\lambda^{*}>0$ and
$\ell_{0},T_{0}\in\mathbb{N}$ such that for all $\lambda\in(0,\lambda^{*})$
and $\ell>\ell_{0},T>T_{0}$

$\displaystyle\lVert
W_{\ell}^{\top}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}<\epsilon.$

###### Proof.

Let $y:M\to y(M)\subset(\mathbb{R}^{N}\times\mathbb{R}^{d})$ be defined by

$\displaystyle y(m)=(f\circ\phi^{-1}(m),\omega(m))\ \forall m\in M$

and note that $F_{T}(f\circ\phi^{-1},\omega)=F_{T}\circ y$ and that
$y\in\text{Hom}(M,y(M))$ because $f\in\text{Hom}(M,f(M))$. Now fix
$\epsilon>0$. Let $\mu^{\prime}$ be a measure defined on
$y(M)\subset(\mathbb{R}^{N}\times\mathbb{R}^{d})$ by
$\mu^{\prime}(\sigma)=\mu(y^{-1}(\sigma))$ for all measurable subsets $\sigma$
of $f(M)$. Using the assumption that $\mathcal{F}$ is a class of linear
universal approximators, we can choose $T_{0}$ sufficiently large that for any
$T>T_{0}$ there exists $W_{*}\in\mathbb{R}^{T}$ such that

$\displaystyle\lVert W_{*}^{\top}F_{T}-u\circ
y^{-1}\rVert_{L^{2}(\mu^{\prime})}^{2}<\frac{\epsilon}{3},$

hence (by lemma 3.5)

$\displaystyle\lVert W_{*}^{\top}F_{T}\circ y-u\rVert_{L^{2}(\mu)}^{2}=\lVert
W_{*}^{\top}F_{T}-u\circ
y^{-1}\rVert_{L^{2}(\mu^{\prime})}^{2}<\frac{\epsilon}{3}.$

Now let

$\displaystyle\lambda^{*}=\frac{\epsilon}{3\lVert W_{*}\rVert^{2}}$

and $\lambda\in(0,\lambda^{*})$. Define the sequence
$(W_{\ell})_{\ell\in\mathbb{N}}$ such that, for each $\ell\in\mathbb{N}$, the
vector $W_{\ell}\in\mathbb{R}^{T}$ is the unique minimiser of the regularised
least squares difference

$\displaystyle\frac{1}{\ell}\bigg{(}\sum_{k=0}^{\ell-1}\rVert
W^{\top}F_{T}(f\circ\phi^{k-1}(m_{0}),\omega\circ\phi^{k}(m_{0}))-u\circ\phi^{k}(m_{0})\lVert^{2}$
$\displaystyle+\lambda\lVert W\rVert^{2}$ $\displaystyle\bigg{)}.$

By lemma 3.1, $(W_{\ell})_{\ell\in\mathbb{N}}$ converges as $\ell\to\infty$ to
$W_{\infty}$ which minimises

$\displaystyle\lVert
W^{\top}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}+\lambda\lVert
W\rVert^{2}.$

Now we choose $\ell_{0}$ such that for all $\ell>\ell_{0}$

$\displaystyle\rVert
W_{\ell}^{T}F_{T}(f\circ\phi^{-1},\omega)-W_{\infty}^{T}F_{T}(f\circ\phi^{-1},\omega)\lVert_{L^{2}(\mu)}^{2}<\frac{\epsilon}{3}.$

Now the proof proceeds directly

$\displaystyle\lVert
W_{\ell}^{T}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}$
$\displaystyle=$ $\displaystyle\lVert
W_{\ell}^{T}F_{T}(f\circ\phi^{-1},\omega)-W_{\infty}^{T}F_{T}(f\circ\phi^{-1},\omega)$
$\displaystyle+$ $\displaystyle
W_{\infty}^{T}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}$
$\displaystyle\leq$ $\displaystyle\lVert
W_{\ell}^{T}F_{T}(f\circ\phi^{-1},\omega)-W_{\infty}^{T}F_{T}(f\circ\phi^{-1},\omega)\rVert_{L^{2}(\mu)}^{2}$
$\displaystyle+$ $\displaystyle\lVert
W_{\infty}^{T}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}$
$\displaystyle<$ $\displaystyle\frac{\epsilon}{3}+\lVert
W_{\infty}^{T}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}$
$\displaystyle\leq$ $\displaystyle\frac{\epsilon}{3}+\lVert
W_{\infty}^{T}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}+\lambda\lVert
W_{\infty}\rVert^{2}$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{3}+\lVert
W_{*}^{T}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}+\lambda\lVert
W_{*}\rVert^{2}$ $\displaystyle<$
$\displaystyle\frac{\epsilon}{3}+\frac{\epsilon}{3}+\lVert
W_{*}^{T}F_{T}(f\circ\phi^{-1},\omega)-u\rVert_{L^{2}(\mu)}^{2}$
$\displaystyle=$ $\displaystyle\frac{\epsilon}{3}+\frac{\epsilon}{3}+\lVert
W_{*}^{T}F_{T}\circ y-u\rVert_{L^{2}(\mu)}^{2}$ $\displaystyle<$
$\displaystyle\frac{\epsilon}{3}+\frac{\epsilon}{3}+\frac{\epsilon}{3}=\epsilon.$

∎

Theorem 3.6 guarantees an approximation in the $L^{2}(\mu)$ norm, which is
sadly weaker than the $C^{1}$ norm. That is to say, a sequence which converges
in $C^{1}$ also converges in $L^{2}(\mu)$, but the converse does not hold in
general. This distinction is particularly relevant when the problem is chaotic
time series forecasting. In this case, the target function is the _next step
map_ $u=\omega\circ\phi$, and we recursively feed predictions into the state
space map to create a trajectory into the future. An example is the ESN
autonomous phase (equation (1)). A weakness of using ESN autonomous dynamics
for time series forecasting is that small approximation errors accumulate
resulting in a predicted trajectory that diverges from the true trajectory in
the far future. That said, Hart et al. [2019] show that under certain
conditions (crucially that the next step map $u=\omega\circ\phi$ is well
approximated in the $C^{1}$ norm) the ESN autonomous phase will adopt dynamics
that are topologically conjugate to the original dynamical system.

We must conclude that least squares regression does not guarantee a
topologically conjugate autonomous phase, but we note that real data sets are
contaminated by noise and finite precision arithmetic where an $L^{2}(\mu)$
approximation may be most suitable. Moreover, computing the (regularised)
least squares solution using the SVD decomposition, or some other algorithm,
is much faster than minimising the maximal pointwise distance, which may be
necessary to yield a good $C^{1}$ approximation. Indeed, despite the
theoretical limitations of the regularised least squares approach it seems to
work well in practice. In fact we can interpret bad $C^{1}$ approximations in
the parlance of machine learning as overfitted solutions, as they fit the
training data well, in exactly the terms that we define a good fit, but may
fail to make good predictions about the unseen future.

###  3.4 Convergence rate of the time average to the space average

Theorem 3.6 guarantees, under appropriate conditions, that with sufficiently
many neurons $T$ and a sufficiently many training data $\ell$ we can obtain an
arbitrarily good $L^{2}(\mu)$ approximation of a target function $u$. It is
natural to wonder how many training data is required to achieve a given
$L^{2}(\mu)$ approximation. To answer this, we turn our attention to the
convergence rate of the time average to the space average

$\displaystyle\lim_{\ell\to\infty}\frac{1}{\ell}\sum_{k=0}^{\ell-1}s\circ\phi^{k}(m_{0})=\int_{M}s\
d\mu$ (2)

as the timespan over which training data is collected grows. We want a uniform
estimate for the rate of convergence for $s$ over all ergodic maps $\phi$.
Unfortunately, no such estimate can possibly exist. Kachurovskii [1996]
presents negative results that (in the author’s words) _leave no hope that
estimates of the rate of convergence depending only on the averaged function
$s$ can be obtained in ergodic theorems_. The negative results presented by
Kachurovskii [1996] prove that the amount of training data required is
strictly dependant on the dynamical system.

Though we cannot say exactly how many data points we need for a good
$L^{2}(\mu)$ approximation, the central limit theorem for ergodic dynamical
systems suggests that for an initial point chosen uniformly over the invariant
measure of $\phi$, the difference between the finite time average and space
average converges to a mean $0$ normal distribution with standard deviation
$1/\sqrt{\ell}$. This is made precise by the central limit theorem for ergodic
dynamical systems. Before we state the theorem, we recall the definition of
Hölder continuity.

######  Definition 3.7.

(Hölder continuous) Let $(M,d)$ be a metric space. A map $s:M\to\mathbb{R}$ is
called _Hölder continuous_ if there exist constants $p\in(0,1]$ and $K>0$ such
that

$\displaystyle\lVert s(m)-s(m^{\prime})\rVert\leq Kd(m,m^{\prime})^{p}$

for all $m,m^{\prime}\in M$.

######  Theorem 3.8.

(Central limit theorem for ergodic dynamical systems) Let $\phi:M\to M$ be
ergodic with respect to the probability space $(M,\Sigma,\mu)$. Let $X_{0}$ be
a uniform random variable with respect to the space $(M,\Sigma,\mu)$. Let
$s\in L^{1}(\mu)(M,\mathbb{R})$ be Hölder continuous and denote the space
average of $s$ by

$\displaystyle\mathbb{E}[s]:=\int_{M}s\ d\mu.$

Let the random variables $X_{j}:=s\circ\phi^{j}(X_{0})$ for
$j=0,\ldots,\ell-1$ and denote the partial sum
$S_{\ell}=X_{0}+\cdots+X_{\ell-1}$. Then, for some $\sigma>0$, the partial sum
$S_{\ell}$ satisfies the central limit theorem:

$\displaystyle\lim_{\ell\to\infty}\mu\bigg{(}\bigg{\\{}\frac{S_{\ell}-\ell\mathbb{E}[s]}{\sqrt{\ell}}\leq
z\bigg{\\}}\bigg{)}=\frac{1}{2\pi\sigma}\int_{-\infty}^{z}\mathrm{e}^{-\frac{\tau^{2}}{2\sigma^{2}}}\
d\tau$

almost surely, or in other words $(S_{\ell}-\ell\mathbb{E}[s])/\sqrt{\ell}$
converges in law to $\mathcal{N}(0,\sigma^{2})$.

###### Proof.

Camí [2010]. ∎

To see the connection between the central limit theorem and the work in this
paper, suppose we choose a map $s$ that returns the matrix vector pair

$\displaystyle s(m_{0})$
$\displaystyle=\bigg{(}\big{[}f(m_{0})f^{\top}(m_{0})+I\lambda\big{]},f(m_{0})u(m_{0})\bigg{)}$
$\displaystyle=:(\Sigma_{0},v_{0}),$

and define a sequence of pairs with $\ell$th pair

$\displaystyle(\Sigma_{\ell},v_{\ell}):=\frac{1}{\ell}\sum_{k=0}^{\ell-1}s\circ\phi^{k}(m_{0}).$

Then it follows that

$\displaystyle W_{\ell}=\Sigma_{\ell}^{-1}v_{\ell}$

is the linear readout layer obtained by regularised least squares regression
using $\ell$ data points. Furthermore, it follows from the central limit
theorem that for random initial points $m_{0}$ (distributed uniformly with
respect to the invariant measure $\mu$) the sequence
$(\Sigma_{\ell},v_{\ell})_{\ell\in\mathbb{N}}$ converges in law to a
(multivariate) normal distribution, with variance converging with order
$1/\ell$, and mean $(\Sigma,v)$ which satisfies

$\displaystyle W_{\infty}=\Sigma^{-1}v.$

We note that the convergence of $(\Sigma_{\ell},v_{\ell})_{\ell\in\mathbb{N}}$
to $(\Sigma,v)$ with order $1/\sqrt{\ell}$ does not necessarily imply that
$(W_{\ell})_{\ell\in\mathbb{N}}$ converges to $W_{\infty}$ at the same rate.

##  4 The Lorenz attractor is stably mixing

We have shown that we can approximate, in the $L^{2}(\mu)$ sense, any target
function on an ergodic dynamical system using an ESN and Tikhonov least
squares. This partially explains the success enjoyed by Jaeger [2001], Xi et
al. [2005], Schrauwen et al. [2007], Shi and Han [2007], Yong Song et al.
[2010], Pathak et al. [2017], Løkse et al. [2017], Yeo [2019], Chattopadhyay
et al. [2019], Vlachas et al. [2019], and Hart et al. [2019]. Many authors
including Chattopadhyay et al. [2019] successfully predict the future
observations of the Lorenz system, while Pathak et al. [2017], Vlachas et al.
[2019], and Hart et al. [2019] additionally recover topological invariants
including Lyapunov exponents, fixed point eigenvalues and homology groups. The
authors are successful in their numerical experiments because the Lorenz
attractor is _mixing_ which implies it is ergodic, suggesting the conditions
Theorem 3.6 hold and we can $L^{2}(\mu)$ approximate target functions on the
Lorenz attractor.

Proving that the Lorenz attractor is mixing was a tremendous achievement,
built upon the works of Afraimovich et al. [1977], Guckenheimer and Williams
[1979], Pesin [1992], Williams [1979], and Tucker [1999] culminating with the
seminal paper by Tucker [2002], which resolved Smale’s $14^{\mathrm{th}}$
problem _‘Is the dynamics of the ordinary differential equations of Lorenz
(1963) that of the geometric Lorenz attractor of Williams, Guckenheimer and
Yorke? ’_ [Smale, 1998]. To formalise some of these ideas, we will begin with
the definition of a mixing dynamical system.

######  Definition 4.1.

(Mixing) Let $\phi:M\to M$ be a measure preserving transformation on the
measure space $(M,\Sigma,\mu)$ with $\mu(M)=1$. Then $\phi$ is mixing if for
any $A,B\in\Sigma$

$\displaystyle\lim_{\ell\to\infty}\mu\big{(}A\cap\phi^{-\ell}(B)\big{)}=\mu(A)\mu(B).$

######  Lemma 4.2.

(Mixing implies ergodic) Let $\phi:M\to M$ be a measure preserving
transformation on the measure space $(M,\Sigma,\mu)$ with $\mu(M)=1$. Suppose
$\phi$ is mixing, then $\phi$ is ergodic.

###### Proof.

Suppose $\phi$ is mixing and $A,B\in\Sigma$. Then

$\displaystyle\lim_{\ell\to\infty}\mu\big{(}A\cap\phi^{-\ell}(B)\big{)}$
$\displaystyle=\mu(A)\mu(B)$
$\displaystyle\implies\lim_{\ell\to\infty}\frac{1}{\ell}\sum_{k=0}^{\ell-1}\mu\big{(}A\cap\phi^{-k}(B)\big{)}$
$\displaystyle=\mu(A)\mu(B)$
$\displaystyle\implies\lim_{\ell\to\infty}\frac{1}{\ell}\sum_{k=0}^{\ell-1}\mu\big{(}A\cap\phi^{-k}(A)\big{)}$
$\displaystyle=\mu(A)^{2}.$ (3)

Now suppose $\mu(A)=\mu\big{(}\phi^{-1}(A)\big{)}$. Then (3) reduces to
$\mu(A)=\mu(A)^{2}$ hence $\mu(A)=1$ or $\mu(A)=0$, so $\phi$ is ergodic. ∎

######  Definition 4.3.

(Stably mixing) Let $\phi:M\to M$ be a measure preserving transformation on
the measure space $(M,\Sigma,\mu)$ with $\mu(M)=1$. Then $\phi$ is stably
mixing if sufficiently small $C^{1}$ perturbations of $\phi$ are mixing.

######  Theorem 4.4.

The Lorenz [1963] system

$\begin{array}[]{ll}\dot{\xi}&=\sigma(\upsilon-\xi)\\\
\dot{\upsilon}&=\xi(\rho-\zeta)-\upsilon\\\
\dot{\zeta}&=\xi\upsilon-\beta\zeta\end{array}$ (4)

with parameters $\sigma=10$, $\beta=8/3$, $\rho=28$ admits a robust attractor
that is stably mixing.

###### Proof.

Luzzatto et al. [2005] ∎

Since the Lorenz attractor is stably mixing, so is any sufficiently good
$C^{1}$ approximation to the evolution operator $\phi$, obtained by numerical
methods. Consequently, a numerically approximated Lorenz system is ergodic, by
Lemma 4.2. Thus, we expect that an ESN, trained using Tikhonov least squares,
on a sequence of observations of a numerically integrated trajectory of the
Lorenz attractor will $L^{2}(\mu)$ approximate arbitrary target functions on
the attractor.

##  5 Numerical experiments

Our goal is to use an ESN to learn a mapping from the $\xi$ component of the
Lorenz attractor to the $\zeta$ component. We will sample data from a single
trajectory of the Lorenz attractor. To this end, let
$\phi:\mathbb{R}^{3}\to\mathbb{R}^{3}$ denote a discretisation of the Lorenz
system (4) with time step $\tau$ i.e effectively 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)$. For these initial conditions
and the parameter values as in 4.4, we computed a trajectory for a 40 time
units (i.e. 4000 timesteps), illustrated in Figure 1.

Figure 1: A typical trajectory of the Lorenz system (4) 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 lighest colours at the most recent times
(yellow).

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

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

We select the target function to be $\omega(\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 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}_{\infty}x_{k}$.

Figure 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:

  * 1.

Reservoir size: $T=300$,

  * 2.

Activation function: $\sigma=\tanh$,

  * 3.

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

  * 4.

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

  * 5.

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

Iterating the ESN with observations $z_{k}$ creates a discrete-time sequence
of reservoir states $x_{k}$, illustrated in Figure 3, which shows a projection
of the reservoir states onto their first the principal components.

Figure 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.

We then solved the least squares problem

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

to determine the output layer $W$ using the SVD. This offline learning method
is described by Hansen et al. [2006]. Our aim here is to understand how
increasing the number of data points $\ell$ improves our approximation of the
target function $u$. So we repeated this process with _fewer_ observation-
target pairs, from 300 in increments of 100 up to 4000. For each value of
$\ell$, we compute the best-fit readout layer $W$. We repeated this process
once more for a 20,000 time step (i.e. 200 time unit) trajectory and computed
the readout layer which for this case we denote by $W_{\infty}$, assuming that
it is extremely close to the readout layer we would obtain in the limit of
infinitely many time steps. For each readout layer $W$ obtained using fewer
data points ($300\leq\ell\leq 4000$) we estimated the error on the readout
layer which we denote by WE:

$\displaystyle\text{WE}=\frac{\rVert W-W_{\infty}\lVert}{\lVert
W_{\infty}\rVert},$

and the root mean square error (RMSE) between the targets and the
approximation _for the entire 20,000 point trajectory_

$\displaystyle\text{RMSE}=\sqrt{\frac{1}{20000}\sum_{k=0}^{20000-1}\rVert
W^{\top}x_{k}-u_{k}\lVert^{2}}.$

We expect that as the the number of data points $\ell$ grows the WE and RMSE
will converge. The central limit theorem suggests that the matrix vector pairs
$(\Sigma_{\ell},v_{\ell})_{{}_{\ell\in\mathbb{N}}}$ (which satisfy the Gauss
normal equations $\Sigma_{\ell}W_{\ell}=v_{\ell}$) will converge in law to a
multivariate normal distribution, with standard deviation converging with
order $1/\sqrt{\ell}$; as the number of data points $\ell$ tends to infinity.
This suggests (but does not strictly imply) that the WE and RMSE might
converge at a similar rate. We have been unable to derive expressions for the
convergence of the RMSE and WE and remark that the need to compute $W$ via a
least-squares fit means it is not obvious that these would share the
convergence rate of $\Sigma_{\ell}$ and $v_{\ell}$. Typical numerical results
for the convergence of the RMSE and WE are illustrated in Figure 4.

(a) Log–log plot of the error on the linear readout layer $W$ (vertical axis)
against number of data points (horizontal axis) used to train the readout
layer $W$. The line $y=45/\sqrt{\ell}$ is plotted in black as a guide to the
eye.

(b) Log–log plot of the root mean square error (RMSE) (vertical axis) against
number of data points (horizontal axis) used to train the readout layer $W$.
The line $y=150/\sqrt{\ell}$ is plotted in black as a guide to the eye.

Figure 4: Convergence of the error on the readout layer (WE) and convergence
of the root mean square error (RMSE) displayed in log–log plots. Black lines
indicate convergence with order $1/\sqrt{\ell}$ and are shown in order to
compare the convergence to what might be expected if a central limit theorem
applied.

The figures reveal that the convergence of the RMSE and WE is complicated. We
observe sudden jumps which appear when the Lorenz trajectory switches to a
different wing in the attractor, at such times presumably the ESN rapidly
acquires new independent information which improves the fit. Furthermore, the
convergence at least over this range of trajectory lengths does not
(convincingly) converge with order $1/\sqrt{\ell}$. Since the sudden jumps
occur on a timescale intrinsic to the dynamical system $\phi$ we conclue that
the internal structure of the attractor and its dynamics plays an important
role in the evolution of the error; appealing to the asymptotic behavior may
not always be useful.

We pushed the numerics further, hoping to detect an asymptotic regime by
repeating the numerical experiments with a much longer trajectory. We computed
$W_{\infty}$ for a $100\ 000$ point trajectory and compared this to the $W$
obtained for shorter time series of lengths $\ell=1000,2000,\ldots,98000$. For
each $\ell$ we computed the WE with 10 randomly generated realisations of the
ESN. The results are shown in Figure 5 and are also (sadly) inconclusive;
there is no obvious regime over which the error decreases as a power law.
Sudden decreases as the trajectory switches lobes on the attractor are still
visible, and the rate of convergence remains complicated.

Figure 5: The error on the readout layer (WE) (vertical axis) shown against
the number of data points $\ell$ (horizontal axis). The black line has
equation $y=45/\sqrt{\ell}$ as a guide to the eye. Results for 10 separate
realisations of an ESN are shown.

##  6 Conclusions and future work

The main result of this paper (Theorem 3.6) states that an ESN trained on a
sequence of observations from an ergodic dynamical system (with invariant
measure $\mu$) using Tikhonov least squares will $L^{2}(\mu)$ approximate any
target function $u$. We then summarised the result by Luzzatto et al. [2005]
which implies the Lorenz attractor exists and is mixing, hence ergodic. This
allowed us to conclude that an ESN trained on a sequence of scalar
observations taken from the Lorenz system using Tikhonov least squares should
$L^{2}(\mu)$ approximate the dynamics in the attractor. In section 5 we
simulated the Lorenz system ourselves and designated the $\xi$ and $\zeta$
components observations and targets respectively. We confirmed that as the
number of data points $(\xi_{k},\zeta_{k})$ grew, the approximation of the
target function improved. A good approximation was reached before the number
of data points was large enough for the central limit theorem to (perhaps)
become relevant. This suggests that (perhaps unfortunately) this asymptotic
result may have limited practical use.

We discussed in section 2 that the $L^{2}(\mu)$ norm is weaker than the
$C^{1}$ norm, in the sense that convergence in $C^{1}$ implies convergence in
$L^{2}(\mu)$, while the converse does not hold. This is somewhat unsatisfying,
because (topologically conjugate) time series forecasting requires the
autonomous phase of the ESN to be a $C^{1}$ approximator of the embedded
(structurally stable) dynamics.

It may be a fruitful to develop a training method beyond Tikhonov least
squares that guarantees a $C^{1}$ approximation. Alternatively, it may be
intriguing to explore under what conditions Tikhonov least squares does
provide a sufficiently good $C^{1}$ approximation, which appears to happen
frequently in simulations. Authors including Pathak et al. [2017], Vlachas et
al. [2019], and Hart et al. [2019] have demonstrated that an ESNs trained with
Tikhonov least squares can replicate topological invariants of dynamical
systems like Lyapunov exponents, fixed point eigenvalues, and homology groups,
suggesting a sufficiently good $C^{1}$ approximation was achieved.

Though the $L^{2}(\mu)$ approximation may not be sufficient for topological
results, it may be powerful enough to prove interesting results about ESNs
applied to control problems. We can view a control system as a dynamical
system, for which we have at every state $x\in M$ a set of actions
$a\in\mathcal{A}$ available to us. Then we seek a map $\pi:M\to\mathcal{A}$,
called an optimal controller (in control theory), or an optimal policy (in
reinforcement learning), which maximises some reward function. To determine
the value of a policy $\pi$ it suffices to determine the value function
$u:M\to\mathbb{R}$ which, we can in principal approximate with an ESN from
only partial observations of the control system. Developing algorithms to find
the optimal controller/policy may be a rewarding direction of future work.

We also believe much of the theory presented here could be generalised or
modified for other recurrent neural networks such as long short term memory
networks (LSTMs). LSTMs are used extensively in industry and perform very well
at _context dependant_ time series problems. These are problems where events
that happened a long time in the past may suddenly become important in the
present. The ESN is not well suited to such problems, because the importance
of events necessarily decays (at least) exponentially quickly as we move
further into the past, while the structure of an LSTM sidesteps this problem.
A detailed explanation of the architecture is provided by Gers [1999].
Equations for a peephole LSTMS are listed below

$\displaystyle f_{k}$
$\displaystyle=\varphi_{g}(A_{f}c_{k}+W^{\text{in}}_{f}u_{k})$ $\displaystyle
i_{k}$ $\displaystyle=\varphi_{g}(A_{i}c_{k}+W^{\text{in}}_{i}u_{k})$
$\displaystyle o_{k}$
$\displaystyle=\varphi_{g}(A_{o}c_{k}+W^{\text{in}}_{o}u_{k})$ $\displaystyle
c_{k}$ $\displaystyle=f_{k}\odot
c_{k-1}+i_{k}\odot\varphi_{c}(W^{\text{in}}_{c}u_{k})$ $\displaystyle h_{k}$
$\displaystyle=\varphi_{h}(o_{k}\odot c_{k})$

where $f_{k},i_{k},o_{k},c_{k},h_{k}\in\mathbb{R}^{n}$ are the vectors of the
forget gate, input gate, output gate, cell state, and hidden state (also known
as the output state) associated to the LSTM at time $k$. Next,
$u_{k}\in\mathbb{R}$ is the scalar input of the LSTM at time $k$ and
$\varphi_{g}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a componentwise sigmoid
function, $\varphi_{c}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the componentwise
$\tanh$ function, and $\varphi_{h}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is some
function that is usually the identity map. $A_{f},A_{i},A_{o},A_{c}$ are
$n\times n$ matrices and
$W^{\text{in}}_{f},W^{\text{in}}_{i},W^{\text{in}}_{o},W^{\text{in}}_{c}$ are
$1\times n$ matrices. Finally, the symbol $\odot$ here represents the Hadamard
product (taking componentwise product of 2 vectors).

We can see that LSTMs admit ESNs as a special case by fixing $A_{f}=0$,
$W^{\text{in}}_{f}=0$, $b_{f}=0$, $b_{i}=0$, $b_{c}=\text{arc}\tanh(1/2)$,
$W^{\text{in}}_{c}=0$. It may therefore interest the academic community
studying LSTMs, as well as those with industrial applications in mind, to
generalise the theory of ESNs presented here and elsewhere to LSTMs.

One shortcoming of Echo State Networks (that is typical for a machine learning
paradigm) is that physical information about the underlying dynamical system
is typically ignored. The question of how one might integrate some basic
knowledge of the underlying dynamical system into the ESN architecture was
recently explored numerically by Huhn and Magri [2020] and Doan et al. [2020].
Developing their ideas further may be an intriguing direction of future work.

## Acknowledgements

We thank the examiners of A.G. Hart’s PhD confirmation viva, Alastair Spence
and Chris Guiver, who offered helpful criticism of much of the material which
formed the basis of this paper. We are grateful to an anonymous reviewer for
their comments and careful reading of the manuscript which have significantly
helped to improve its presentation. A.G. Hart is supported by a scholarship
from the EPSRC Centre for Doctoral Training in Statistical Applied Mathematics
at Bath (SAMBa), under the project EP/L015684/1.

## References

  * Afraimovich et al. [1977] Afraimovich, V.S., Bykov, V.V., Shilnikov, L.P., 1977.  On the origin and structure of the Lorenz attractor.  Akademiia Nauk SSSR Doklady 234, 336–339.
  * Birkhoff [1931] Birkhoff, G.D., 1931.  Proof of the ergodic theorem.  Proceedings of the National Academy of Sciences 17, 656–660.  doi:10.1073/pnas.17.2.656, arXiv:https://www.pnas.org/content/17/12/656.full.pdf.
  * Bogachev [2007] Bogachev, V., 2007.  Measure Theory. volume 1.  Springer-Verlag.
  * Camí [2010] Camí, J.A.L., 2010.  Ergodic theory.
  * Chattopadhyay et al. [2019] Chattopadhyay, A., Hassanzadeh, P., Palem, K., Subramanian, D., 2019\.  Data-driven prediction of a multi-scale lorenz 96 chaotic system using a hierarchy of deep learning methods: Reservoir computing, ann, and rnn-lstm.  arXiv:1906.08829 .
  * Doan et al. [2020] Doan, N.A.K., Polifke, W., Magri, L., 2020.  Learning hidden states in a chaotic system: A physics-informed echo state network approach.  arXiv:2001.02982 .
  * Gers [1999] Gers, F., 1999.  Learning to forget: continual prediction with lstm.  IET Conference Proceedings , 850--855(5).
  * Grigoryeva et al. [2020] Grigoryeva, L., Hart, A., Ortega, J.P., 2020.  Chaos on compact manifolds: Differentiable synchronizations beyond takens.  arXiv:2010.03218 .
  * Guckenheimer and Williams [1979] Guckenheimer, J., Williams, R.F., 1979\.  Structural stability of lorenz attractors.  Publication Mathématique de l’I.H.É.S. 50, 59--72.
  * Gürel and Egert [2010] Gürel, T., Egert, S.R.U., 2010\.  Functional identification of biological neural networks using reservoir adaptation for point processes.  Journal of Computational Neuroscience , 279–299.
  * Hansen et al. [2006] Hansen, P.C., Nagy, J.G., O’leary, D.P., 2006.  Deblurring Images: Matrices, Spectra, and Filtering.  SIAM.
  * Hart et al. [2019] Hart, A.G., Hook, J.L., Dawes, J.H.P., 2019.  Embedding and approximation theorems for echo state networks.  arXiv:1908.05202 .
  * Hornik et al. [1990] Hornik, K., Stinchcombe, M., White, H., 1990.  Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks.  Neural Networks 3, 551 -- 560.  doi:https://doi.org/10.1016/0893-6080(90)90005-6.
  * Huhn and Magri [2020] Huhn, F., Magri, L., 2020.  Learning ergodic averages in chaotic systems.  arXiv:2001.04027 .
  * Ilies et al. [2007] Ilies, I., Jaeger, H., Kosuchinas, O., Rincon, M., Sakenas, V., Vaskevicius, N., 2007\.  Stepping forward through echoes of the past: forecasting with echo state networks .
  * Jaeger [2001] Jaeger, H., 2001.  The “echo state” approach to analysing and training recurrent neural networks .
  * Kachurovskii [1996] Kachurovskii, A.G., 1996.  The rate of convergence in ergodic theorems.  Russian Mathematical Surveys 51, 653--703.
  * Kocarev and Parlitz [1996] Kocarev, L., Parlitz, U., 1996\.  Generalized synchronization, predictability, and equivalence of unidirectionally coupled dynamical systems.  Phys. Rev. Lett. 76, 1816--1819.
  * Lin et al. [2009] Lin, X., Yang, Z., Song, Y., 2009.  Short-term stock price prediction based on echo state networks.  Expert Systems with Applications 36, 7313 -- 7317.  doi:https://doi.org/10.1016/j.eswa.2008.09.049.
  * Løkse et al. [2017] Løkse, S., Bianchi, F.M., Jenssen, R., 2017.  Training echo state networks with regularization through dimensionality reduction.  Cognitive Computation 9, 364--378.  doi:10.1007/s12559-017-9450-z.
  * Lorenz [1963] Lorenz, E.N., 1963.  Deterministic nonperiodic flow.  Journal of the Atmospheric Sciences 20, 130--141.
  * Luzzatto et al. [2005] Luzzatto, S., Melbourne, I., Paccaut, F., 2005.  The lorenz attractor is mixing.  Communications in Mathematical Physics 260, 393--401.  doi:10.1007/s00220-005-1411-9.
  * Maass et al. [2002] Maass, W., Natschläger, T., Markram, H., 2002.  Real-time computing without stable states: A new framework for neural computation based on perturbations.  Neural Computation 14, 2531--2560.  doi:10.1162/089976602760407955.
  * Pathak et al. [2017] Pathak, J., Lu, Z., Hunt, B.R., Girvan, M., Ott, E., 2017\.  Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data.  Chaos 27.
  * Peng et al. [2019] Peng, H., Chen, C., Lai, C.C., Wang, L.C., Han, Z., 2019\.  A predictive on-demand placement of uav base stations using echo state network.  arXiv:1909.11598 .
  * Pesin [1992] Pesin, Y.B., 1992.  Dynamical systems with generalized hyperbolic attractors: hyperbolic, ergodic and topological properties.  Ergodic Theory and Dynamical Systems 12, 123–151.  doi:10.1017/S0143385700006635.
  * Schrauwen et al. [2007] Schrauwen, B., Verstraeten, D., Van Campenhout, J., 2007.  An overview of reservoir computing: theory, applications and implementations, in: Proceedings of the 15th European Symposium on Artificial Neural Networks. p. 471-482 2007, pp. 471--482.
  * Ser et al. [2020] Ser, J.D., Lana, I., Manibardo, E.L., Oregi, I., Eneko Osaba, J.L.L., Bilbao, M.N., Vlahogianni, E.I., 2020.  Deep Echo State Networks for Short-Term Traffic Forecasting: Performance Comparison and Statistical Assessment.  arXiv:2004.08170 .
  * Shi and Han [2007] Shi, Z., Han, M., 2007.  Support vector echo-state machine for chaotic time-series prediction.  IEEE Transactions on Neural Networks 18, 359--372.
  * Skowronski and Harris [2007] Skowronski, M.D., Harris, J.G., 2007\.  Automatic speech recognition using a predictive echo state network classifier.  Neural Networks 20, 414 -- 423.  doi:https://doi.org/10.1016/j.neunet.2007.04.006. echo State Networks and Liquid State Machines.
  * Smale [1998] Smale, S., 1998.  Mathematical problems for the next century.  The Mathematical Intelligencer 20, 7–15.
  * Tanaka et al. [2019] Tanaka, G., Yamane, T., Héroux, J.B., Nakane, R., Kanazawa, N., Takeda, S., Numata, H., Nakano, D., Hirose, A., 2019.  Recent advances in physical reservoir computing: A review.  Neural Networks 115, 100 -- 123.  doi:https://doi.org/10.1016/j.neunet.2019.03.005.
  * Tong et al. [2007] Tong, M.H., Bickett, A.D., Christiansen, E.M., Cottrell, G.W., 2007\.  Learning grammatical structure with echo state networks.  Neural Networks 20, 424 -- 432.  doi:https://doi.org/10.1016/j.neunet.2007.04.013. echo State Networks and Liquid State Machines.
  * Tucker [1999] Tucker, W., 1999.  The lorenz attractor exists.  Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 328, 1197 -- 1202.
  * Tucker [2002] Tucker, W., 2002.  A rigorous ode solver and smale’s 14th problem.  Foundations of Computational Mathematics 2, 53--117.
  * Verzelli et al. [2020] Verzelli, P., Alippi, C., Livi, L., 2020.  Learn to synchronize, synchronize to learn.  arXiv:2010.02860 .
  * Vlachas et al. [2019] Vlachas, P.R., Pathak, J., Hunt, B.R., Sapsis, T.P., Girvan, M., Ott, E., Koumoutsakos, P., 2019.  Forecasting of spatio-temporal chaotic dynamics with recurrent neural networks: a comparative study of reservoir computing and backpropagation algorithms.  arXiv:1910.05266 .
  * Williams [1979] Williams, R.F., 1979.  The structure of lorenz attractors.  Publications Mathématiques de l’IHÉS 50, 73--99.
  * Xi et al. [2005] Xi, J., Shi, Z., 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., pp. 1412--1417 vol. 3.  doi:10.1109/IJCNN.2005.1556081.
  * Yeo [2019] Yeo, K., 2019.  Data-driven reconstruction of nonlinear dynamics from sparse observation.  Journal of Computational Physics 395, 671 -- 689.  doi:https://doi.org/10.1016/j.jcp.2019.06.039.
  * Yildiz et al. [2012] Yildiz, I.B., Jaeger, H., Kiebel, S.J., 2012.  Re-visiting the echo state property.  Neural Networks 35, 1 -- 9.  doi:https://doi.org/10.1016/j.neunet.2012.07.005.
  * Yong Song et al. [2010] Yong Song, Yibin Li, Qun Wang, Caihong Li, 2010\.  Multi-steps prediction of chaotic time series based on echo state network, in: 2010 IEEE Fifth International Conference on Bio-Inspired Computing: Theories and Applications (BIC-TA), pp. 669--672.  doi:10.1109/BICTA.2010.5645205.