capbtabboxtable[][]

# Edge of Stability Echo State Network

Andrea Ceni<EMAIL_ADDRESS>Claudio Gallicchio
<EMAIL_ADDRESS>Department of Computer Science, University of
Pisa, Largo Bruno Pontecorvo, 3 - 56127, IT

###### Abstract

In this paper, we propose a new Reservoir Computing (RC) architecture, called
the Edge of Stability Echo State Network (ES2N). The introduced ES2N model is
based on defining the reservoir layer as a convex combination of a nonlinear
reservoir (as in the standard ESN), and a linear reservoir that implements an
orthogonal transformation. We provide a thorough mathematical analysis of the
introduced model, proving that the whole eigenspectrum of the Jacobian of the
ES2N map can be contained in an annular neighbourhood of a complex circle of
controllable radius, and exploit this property to demonstrate that the ES2N’s
forward dynamics evolves close to the edge-of-chaos regime by design.
Remarkably, our experimental analysis shows that the newly introduced
reservoir model is able to reach the theoretical maximum short-term memory
capacity. At the same time, in comparison to standard ESN, ES2N is shown to
offer a favorable trade-off between memory and nonlinearity, as well as a
significant improvement of performance in autoregressive nonlinear modeling.

###### keywords:

Echo state networks , Reservoir computing , Memory capacity , Recurrent neural
networks , Input-driven systems , Edge of chaos.

##  1 Introduction

Recurrent neural networks (RNNs) [17] are computational models designed to
extract features from data with temporal structures. Applications ranges from
speech recognition to classification of time series. The most common way of
training RNNs is via stochastic gradient descent methods, usually via the
backpropagation through-time algorithm [36]. Unfortunately, these methods come
with a significant computational effort. Modern hardware unleashed the power
of parallel computing techniques, allowing to reduce the computational time of
training deep learning models as RNNs. However, the price to pay is a massive
energy consumption. Moreover, a fundamental limitation of theoretical nature
prevents RNNs to be fully exploited, namely the vanishing/exploding (V/E)
gradient issue [2]. In this regard, an appealing alternative is represented by
Reservoir Computing (RC) [23, 33], a different paradigm of training RNNs
dodging the V/E while being computationally fast, and energy efficient. The
key idea of RC is to inject the input signal into a large random untrained
recurrent layer, the _reservoir_ , from which a readout layer is optimised to
fit the desired target signal. The Echo State Network (ESN) [13, 15] provides
a popular discrete-time class of RC machines. The name ESN recalls the imagery
of the input signal echoing and reverberating within the pool of neuronal
activations in the reservoir, which in turn serves as an high-dimensional
representation of the past history of the input. Although gradients are not
backpropagated in the RC paradigm, the ESN’s forward dynamics are ruled by the
very same equation of conventional RNNs. Thus, ESNs inherit a problem which
closely relates to the V/E gradient problem of plain RNNs, namely the
degradation of memory [12].

In this paper, we propose and analyse a novel RC model, called _Edge of
Stability Echo State Networks (ES 2N)_, that copes with the degradation of
memory of nonlinear reservoir systems. We tackle the degradation of memory of
RNNs from a dynamical system perspective, framing the problem within the _edge
of chaos_ hypothesis. The proposal of this hypothesis can be traced back at
least to the late eighties [19, 27], where it has been observed that extensive
computational capabilities are achieved by adaptive systems whose dynamics are
neither chaotic nor ordered but somewhere in between order and chaos. In this
paper we link these ideas rooted in the study of adaptive systems to the case
of partially randomised RNNs, as ESNs are.

ESNs work “properly” provided with the so called Echo State Property (ESP)
[13, 38]. In coarse terms, the ESP guarantees the ESN to possess a unique
input-driven solution such that all the trajectories originating from
different initial conditions (in the infinite past) synchronise with it (in
present time). Following the imagery, such a unique input-driven solution
would represent the echo of the input signal from the infinite past. The
simplest known criterion to ensure the ESP is to impose the maximum singular
value of the reservoir’s connections matrix to be less than one. This
condition implies straight contraction in the phase space at each time step.
Therefore, any two different internal states of the RNN, when driven by the
same input sequence, will get closer and closer to each other as time flows
ahead. Although stable, such a dynamical system would have truly little margin
to exploit the transient dynamics for computational purposes, due to the
straightforward contraction. An ideal situation would be for the RNN to stay
in a regime of balance between stable contractive dynamics and unstable
chaotic dynamics, i.e. along the _edge of chaos_ [3, 20]. This led the RC
community to adopt the rule of thumb of setting the reservoir matrix to have
spectral radius approximately one. However, this rule of thumb inevitably
overweights the contractive dynamics, especially when considering the action
of the input driving neurons towards the saturation regime. This overweighting
of the contractive dynamics reflects in the degradation of memory in the
forward dynamics of ESNs. A stronger condition would be to have the entire
eigenspectrum of the reservoir matrix to be “around” the complex unit circle.
We study the eigenspectrum of ES2N, and prove in Theorem 3.3 that all the
eigenvalues of an ES2N lie within an annular neighbourhood of the complex unit
circle. The radius of such neighbourhood can be tuned via a specific
hyperparameter of the ES2N that we called the _proximity_ hyperparameter.
Moreover, in the limit of small values of the proximity hyperparameter, we
prove in Theorem 3.5 that the ES2N dynamics narrowly hover over the edge of
chaos while being on average over time in the stable regime, hence the name
Edge of Stability Echo State Network. Through experiments, we show the
empirical advantages of the proposed approach in terms of short-term memory
capacity, memory-nonlinearity tradeoff, and autoregressive time-series
modeling.
The rest of this paper is organized as follows. In Section 2 we introduce the
reader to the RC fundamentals and the classical ESN model. In Section 3 we
propose our ES2N model and provide a theoretical analysis of its dynamics.
Section 4 is dedicated to the experimental results, focusing on the short-term
memory capacity, the trade-off between nonlinearity and memory, and
autoregressive time series modeling. Finally, in Section 5, we discuss our
findings, draw the conclusions, and point out interesting research directions
to explore further.

##  2 The ESN model

RC neural networks [23, 33] identify a class of fastly trainable RNNs, in
which a non-linear untrained dynamical layer is followed by a linear trainable
readout component. In this contribution, we focus on the Echo State Networks
(ESNs) [15, 13] approach within RC, and we recall the well-established [16,
23, 22] formulation of the leaky ESN model, which consists of a nonlinear
reservoir layer made up of leaky recurrent neurons followed by a linear
readout. The equations read as follows:

$\displaystyle x[t]=$
$\displaystyle\alpha\tanh(\rho\mathbf{W}_{r}x[t-1]+\omega\mathbf{W}_{in}u[t])+(1-\alpha)x[t-1],$
(1) $\displaystyle z[t]=$ $\displaystyle\mathbf{W}_{o}x[t].$ (2)

The internal state $x[t]$, input $u[t]$, and output $z[t]$ are, respectively,
$N_{r}$-dimensional, $N_{i}$-dimensional and $N_{o}$-dimensional vectors of
real values. The nonlinearity is expressed by the element-wise applied
hyperbolic tangent function $\tanh$, and the system is typically initialised
in the origin, i.e., $x[0]=0$. Matrices $\mathbf{W}_{r},\mathbf{W}_{in}$ are,
respectively, the recurrent reservoir weight matrix and the input weight
matrix, both randomly instantiated and left untouched. In this paper, we
initialise $\mathbf{W}_{in}$ with i.i.d. random uniformly distributed entries
in $(-1,1)$, and $\mathbf{W}_{r}$ with i.i.d. normally distributed entries
with zero mean and standard deviation $\dfrac{1}{\sqrt{N_{r}}}$. This
initialisation scheme for $\mathbf{W}_{r}$ ensures that, for large $N_{r}$,
the spectral radius of $\mathbf{W}_{r}$ is approximately 1, thanks to the
circular law from random matrix theory [25]. This allows us to interpret the
hyperparameter $\rho$ as the spectral radius (i.e., the largest eigenvalue in
modulus) of the effective recurrent matrix, i.e. of $\rho\mathbf{W}_{r}$.
While, the input scaling $\omega$ is an hyperparameter entitled to rescale the
weight of the current input into the reservoir dynamics. The leaky ESN owes
its name to the presence of the hyperparameter $\alpha\in(0,1]$, the _leak
rate_ , that is designated to modify the time scale of the ESN dynamics
according to the input at hand. In this paper we always consider the recurrent
reservoir matrix $\mathbf{W}_{r}$ a fully connected weight matrix.

###  2.1 Training ESNs via ridge regression

Given a training set of input-target samples,
$\\{u[t],y[t]\\}_{t=1,\ldots,T}$, we train a leaky ESN by means of optimising
the readout matrix $\mathbf{W}_{o}$ in order to solve the linear regression
problem $y[t]=\mathbf{W}_{o}x[t]$. Usually, this is achieved via _ridge
regression (or Tichonov regularisation)_ [23] by means of the following
formula:

$\mathbf{W}_{o}=\mathbf{Y}\mathbf{X}^{T}(\mathbf{X}\mathbf{X}^{T}+\mu\mathbf{I})^{-1},$
(3)

where $\mathbf{X}$ is the matrix of dimension $N_{r}\times T$, containing all
the internal states $x[t]$ of the ESN driven by the input $u[t]$ for
$k=1,\ldots,T$, $\mathbf{Y}$ the matrix of dimension $N_{o}\times T$,
containing all the target values $y[t]$, $\mathbf{I}$ is the identity matrix
of dimension $N_{r}\times N_{r}$, and $\mu$ is the regularisation parameter.
Throughout the paper we denote the transposed of a matrix $\mathbf{X}$ with
$\mathbf{X}^{T}$. Usually, when dealing with regression tasks, a number of
initial time steps are discarded for the computation of eq. (3) to allow the
reservoir to “warm up”. This is to ensure that the ESN transient dynamics are
washed out, so that the internal dynamics of the ESN get linked with the
driving input.

###  2.2 Echo state property

ESNs work under the fundamental assumption of the echo state property (ESP), a
condition ensuring a unique stable input-driven response [13]. Roughly
speaking, the ESP guarantees that the internal state $x[t]$ is uniquely
determined by the entire past history of the input signal. The easiest
condition to ensure the ESP is to set a reservoir such that

$\lVert\rho\mathbf{W}_{r}\rVert<1,$ (4)

where, for a given matrix $\mathbf{M}\in\mathbb{R}^{N\times N}$,
$\lVert\mathbf{M}\rVert$ denotes the matrix norm induced by the Euclidean norm
in $\mathbb{R}^{N}$, or equivalently the maximum singular value of
$\mathbf{M}$. For the rest of the paper, we will always use
$\lVert\cdot\rVert$ to denote either the Euclidean norm on $\mathbb{R}^{N}$,
or the matrix norm induced by that, depending on whether the argument of the
norm is respectively a vector or a matrix.
The condition expressed by eq. (4) implies contraction of the internal states
at each time step for whatever input; that is a strong stability condition
implying a Markovian state space organisation [9]. However, stability is not
the only property we crave from an ideal recurrent neural system. We strive
for a stable dynamical system that is able to recall the information conveyed
from the external input for as long as possible. In other words, we want to
keep the ESP while staying close to the border of its domain of existence.
This led the RC community to adopt the less restrictive rule of thumb of
setting the spectral radius of $\rho\mathbf{W}_{r}$ approximately (or often
slightly less than) $1$, e.g. initialising $\mathbf{W}_{r}$ to have spectral
radius $1$ and then setting setting the hyperparameter $\rho\approx 1$.
Intuitively, $\rho$ controls the amount of nonlinearity into the reservoir and
the contribution of the past activations. A small $\rho$ promotes stable
dynamics at the expenses of forgetting faster the past internal activations,
thus amplifying the fading memory property. However, as previous works
demonstrated [38], the only constraint of $\rho<1$ is not generally sufficient
to guarantee the ESP in an input-driven ESN; and, due to the presence of the
external input driving the dynamics, it is not even necessary. In fact, the
ESP is not a property of the reservoir alone, but rather of the reservoir plus
the forcing input. Some efforts to comprehend the input-dependence in the
analysis of the ESP in RC have been exerted in [24, 35, 8, 5]. As demonstrated
in the literature, the ESP might hold even with large values of $\rho$, as
long as the amplitude of the inputs (or the hyperparameter $\omega$) are large
enough to counterbalance the effect. As a consequence, $\rho$ should be
optimised in synergy with the input scaling $\omega$.

Although ESNs present many advantages, it remains unclear how to tune a priori
an ESN close to the edge of chaos keeping the benefits of a stable nonlinear
computational system. Driven by the need to reconcile the properties of
stability, nonlinearity, long-term memory, and ease of computation, we
introduce in the next section a new RC architecture.

##  3 The ES2N model

We propose a variant of ESN, called Edge of Stability ESN (ES2N), whose
equations with linear readout read as follows:

$\displaystyle x[t]=$
$\displaystyle\beta\;\phi(\rho\mathbf{W}_{r}x[t-1]+\omega\mathbf{W}_{in}u[t])+(1-\beta)\;\mathbf{O}x[t-1],$
(5) $\displaystyle z[t]=$ $\displaystyle\mathbf{W}_{o}x[t],$ (6)

Where $x[t]$, $u[t]$, $z[t]$, $\mathbf{W}_{r}$, $\mathbf{W}_{in}$, and
$\mathbf{W}_{o}$ are as in eq. (1) and eq. (2). Here, in addition, the matrix
$\mathbf{O}$ is a randomly generated _orthogonal_ matrix, and $\beta\in(0,1]$
a real valued hyperparameter that we call _proximity_. All the random
orthogonal matrices in the experiments of this paper have been obtained by
generating first a random matrix $\mathbf{D}$ of the desired dimension
$N_{r}\times N_{r}$ with i.i.d. uniformly random entries in $(-1,1)$, hence
performing a QR decomposition of $\mathbf{D}=\mathbf{QR}$, and taking the
resulting $N_{r}\times N_{r}$ orthogonal matrix $\mathbf{Q}$ as the random
orthogonal matrix $\mathbf{O}$ in eq. (5). Although apparently the
hyperparameters $\alpha$ of a leaky ESN of eq. (1) and $\beta$ of an ES2N of
eq. (5) share the same position in the equation, they play quite a different
role in the dynamics of the RNNs, as it will be evident after our theoretical
and experimental analysis. In fact, while the value of $\alpha$ in eq. (1) is
intended to slow down the speed of the reservoir dynamics relative to those of
the input signal [16], as we will see later in this section, the role of
$\beta$ in eq. (5) is to determine the proximity of the reservoir dynamics to
the edge of chaos. With the setting in eq. (5), the reservoir of the ES2N
results into a convex combination of a standard nonlinear input-driven
reservoir (first term on the right hand side) and a linear orthogonal input-
free reservoir (second term on the right hand side). Note that for
$\beta\rightarrow 0$, and $\alpha\rightarrow 0$, the ES2N model of eq. (5) and
the leaky ESN model of eq. (1) both coincide to the standard (non-leaky) ESN
model.

###  3.1 Edge of chaos in ES2N

In this section we present a detailed mathematical analysis of the structure
of the eigenspectrum for the newly introduced ES2N model, and its maximum
local Lyapunov exponent. We start introducing our notation. Denoting the ES2N
reservoir map as
$G(u,x)=\beta\phi(\rho\mathbf{W}_{r}x+\omega\mathbf{W}_{in}u)+(1-\beta)\mathbf{O}x$,
where we use $\phi$ to generally indicate the nonlinear activation function,
then the corresponding Jacobian reads as follows:

$\dfrac{\partial G}{\partial
x}(u,x)=\beta\,\,\mathbf{D}_{(u,x)}\,\,\rho\mathbf{W}_{r}\,\,+\,\,(1-\beta)\,\,\mathbf{O},$
(7)

where we defined the following diagonal matrix

$\mathbf{D}_{(u,x)}:=\operatorname*{diag}(\phi^{\prime}(\rho\mathbf{W}_{r}x+\omega\mathbf{W}_{in}u)),$
(8)

whose supremum of its norm we denote as follows

$\gamma:=\sup_{u,x}\lVert\mathbf{D}_{(u,x)}\rVert.$ (9)

Note that whenever $|\phi^{\prime}|\leq 1$, then it holds that
$\gamma\in[0,1]$. Considering a specific input-driven trajectory
$\\{(u[t+1],x[t])\\}_{t=0,\ldots,T-1}$, we denote as

$\mathbf{J}[t]:=\dfrac{\partial G}{\partial x}(u[t+1],x[t])$ (10)

the Jacobian map of eq. (7) evaluated along the input-driven trajectory.
Finally, we will denote the maximum singular value of the matrix
$\rho\mathbf{W}_{r}$ as follows

$\sigma:=\lVert\rho\mathbf{W}_{r}\rVert.$ (11)

We start our analysis providing a sufficient condition for the ESP to hold for
the ES2N model, as stated in the proposition below.

###### Proposition 3.1

Let us assume that $|\phi^{\prime}|\leq 1$, e.g., $\phi=\tanh$. If $\sigma<1$
then the ES2N has the ESP for all inputs.

Proof.
Following the same proof of [13, Proposition 3], for an ES2N to have the ESP
for all input, it is sufficient to demonstrate that $\lVert\dfrac{\partial
G}{\partial x}(u,x)\rVert<1$ for all $x,u$. The following holds

$\lVert\dfrac{\partial G}{\partial
x}(u,x)\rVert\leq\beta\lVert\mathbf{D}_{(u,x)}\rVert\lVert\rho\mathbf{W}_{r}\rVert+(1-\beta)\lVert\mathbf{O}\rVert\leq\beta\sigma+(1-\beta)=1-\beta(1-\sigma).$

The first inequality is the triangle inequality, the second holds since
$\lVert\mathbf{D}_{(u,x)}\rVert\leq 1$ for all activation functions such that
$|\phi^{\prime}|\leq 1$, and because of the isometric property of orthogonal
matrices which implies $\lVert\mathbf{O}\rVert=1$. Now the thesis holds by
hypothesis since $\sigma<1\Longrightarrow 1-\beta(1-\sigma)<1$.

###### Remark 3.1

Proposition 3.1 informs us that, regardless of the particular value of
$\beta\in(0,1]$, the ES2N model owns the ESP under the same contractivity
condition of a standard ESN.

In an ES2N with a relatively low value of $\beta$, the dependence from the
hyperparameters $\rho$ and $\omega$ is attenuated. In fact, the forward
dynamics of such an ES2N evolve by design close to the edge of chaos. Since
the Jacobian of the ES2N map in eq. (7) is itself a convex combination, we can
weight more the orthogonal part of the ES2N approaching the hyperparameter
$\beta$ to $0$, which by a continuity argument tune the eigenvalues of the
Jacobian to stay closer to the unitary circle, pretty much regardless of the
spectral radius of the reservoir and the scaling of the input matrix (provided
that the spectral radius is bounded). We formalise this intuition exploiting
the Bauer-Fike theorem [1], that we report here for ease of comprehension.

######  Theorem 3.2 (Bauer-Fike)

Let $\mathbf{A}$ be a diagonalisable matrix, and let $\mathbf{V}$ be the
eigenvector matrix such that
$\mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}$ where
$\mathbf{\Lambda}$ is the diagonal matrix of the eigenvalues of $\mathbf{A}$.
Let $\mathbf{E}$ be an arbitrary matrix of the same dimension of $\mathbf{A}$.
Then, for all $\mu$ eigenvalues of $\mathbf{A}+\mathbf{E}$, there exists an
eigenvalue $\lambda$ of $\mathbf{A}$ such that

$|\lambda-\mu|\leq\lVert\mathbf{V}\rVert\lVert\mathbf{V}^{-1}\rVert\lVert\mathbf{E}\rVert.$
(12)

The Bauer-Fike theorem allows us to derive the following characterisation of
the eigenspectrum of the Jacobian of the ES2N map.

###### Theorem 3.3

Let us consider an ES2N model whose state update equation is given by eq. (5),
and recall the definitions of $\sigma$ and $\gamma$ of eq. (11) and eq. (9).
Then, the eigenspectrum of the Jacobian of the ES2N map is confined in the
annular neighbourhood of radius $\beta\gamma\sigma$ of the circle centered in
the origin of radius $1-\beta$. In formulas, for each eigenvalue $\mu$ of the
Jacobian matrix of eq. (7) there exists a $\theta\in[0,2\pi)$ such that

$|(1-\beta)e^{i\theta}-\mu|\leq\beta\gamma\sigma.$ (13)

Proof.
Define $\mathbf{A}=(1-\beta)\mathbf{O}$, and
$\mathbf{E}=\beta\mathbf{D}_{(u,x)}\rho\mathbf{W}_{r}$, so that the Jacobian
of the ES2N model is $\dfrac{\partial G}{\partial
x}(u,x)=\mathbf{A}+\mathbf{E}$. The matrix $\mathbf{O}$ is orthogonal, hence
there exists an unitary matrix $\mathbf{V}$ such that
$\mathbf{O}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}$, where
$\mathbf{\Lambda}$ is the diagonal matrix of the eigenvalues of $\mathbf{O}$.
In particular, each eigenvalue of $\mathbf{O}$ is of the kind $e^{i\theta}$,
for some argument $\theta\in[0,2\pi)$, due to the orthogonality of
$\mathbf{O}$. Therefore,
$\mathbf{A}=(1-\beta)\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}$, and all
eigenvalues of $\mathbf{A}$ have the form $(1-\beta)e^{i\theta}$, for some
argument $\theta\in[0,2\pi)$. In other words, all the eigenvalues of
$\mathbf{A}$ lie on the complex circle centered in the origin with radius
$(1-\beta)$. Now, since $\mathbf{V}$ is unitary, we have that
$\lVert\mathbf{V}\rVert=\lVert\mathbf{V}^{-1}\rVert=1$. Therefore, eq. (12)
tells us that for each $\mu$ eigenvalue of $\dfrac{\partial G}{\partial
x}(u,x)$ there exists a complex number $\lambda=(1-\beta)e^{i\theta}$ such
that

$|(1-\beta)e^{i\theta}-\mu|\leq\lVert\beta\mathbf{D}_{(u,x)}\rho\mathbf{W}_{r}\rVert\leq\beta\gamma\lVert\rho\mathbf{W}_{r}\rVert=\beta\gamma\sigma.$
(14)

This implies that each eigenvalue of $\dfrac{\partial G}{\partial x}(u,x)$
must be inside the circle of radius $1-\beta+\beta\gamma\sigma$, and outside
the circle of radius $1-\beta-\beta\gamma\sigma$, both centered in zero, which
is the thesis.

A similar argument leads to the following characterisation of a leaky ESN’s
eigenspectrum.

###### Corollary 3.4

Let us consider a leaky ESN model whose state update equation is given by eq.
(1). Then, the eigenspectrum of the Jacobian of the leaky ESN map is confined
in a neighbourhood of radius $\alpha\gamma\sigma$ of the complex number
$(1-\alpha,0)$. In formulas, for each eigenvalue $\mu$ of the Jacobian matrix
of a leaky ESN’s map it holds

$|(1-\alpha)-\mu|\leq\alpha\gamma\sigma.$ (15)

Proof. The proof follows the same steps of the proof of Theorem (3.3),
replacing $\beta$ with $\alpha$, and $\mathbf{O}$ with the identity matrix
$\mathbf{I}$. Now, since $\mathbf{A}=(1-\alpha)\mathbf{I}$ has a unique
eigenvalue $1-\alpha$ with multiplicity the dimension of $\mathbf{A}$, then
eq. (13) reads $|(1-\alpha)-\mu|\leq\alpha\gamma\sigma$, which is the thesis.

###### Remark 3.2

The key feature that differentiates an ES2N model from a leaky ESN model is
that the latter has an eigenspectrum that shrinks towards $1$ (for small
values of $\alpha$), while the former has an eigenspectrum that tends to
dispose along the unitary circle (for small values of $\beta$). This spread of
the eigenspectrum can add richness and diversity to the resulting reservoir
dynamics, while the “collapse” of the eigenspectrum towards $1$ might harm the
expressiveness of the recurrent neural dynamics in some tasks like retrieving
memory. In Fig. 1 the eigenspectrum of the Jacobian of the ES2N model and
leaky ESN model are plotted for various combinations of the hyperparameters
$\rho,\omega,\beta$ and $\alpha$.

Figure 1: Eigenspectrum of the Jacobian of the proposed ES2N model of eq. (5)
(in red) and the leaky ESN model of eq. (1) (in green) for various
combinations of $\rho,\omega,\beta$, and $\alpha$. In black the unitary
complex circle. The input driving the reservoirs is set to the constant $1$,
which is then scaled by the value of the hyperparmameter $\omega$.

In the context of autonomous dynamical systems, the notion of Maximum Lyapunov
Exponent (MLE) is widely adopted to detect whether a system is sensitive to
initial conditions [34]. The idea is to consider two infinitesimally close
initial conditions and measure the average (over time) maximum expansion rate
of the distance between those two initial conditions. Although there exists a
spectrum of Lyapunov Exponents, exactly one for each dimension of the system,
the maximum among them is the most important. If the MLE is less than zero, it
means that any perturbation of an initial condition gets damped on average
over time. On the contrary if the MLE is greater than zero, then there exists
at least one direction in tangent space along which the perturbation gets
magnified on average over time; this characteristic expansion behaviour of the
linearised system is often one of the basic ingredients for the definition of
chaotic dynamics [26, 32]. Therefore, in the literature the edge of chaos is
often defined as the locus of parameters where the MLE is exactly zero [4].
However, when we allow an external input to drive the dynamics, as usual in
RNNs, the definition of the Lyapunov exponents becomes input-dependent. In the
literature [33, 21], one way to assess the degree of regularity of the
dynamics of an input-driven system is to compute its Maximum Local Lyapunov
Exponent (MLLE) on a given input-driven trajectory, whose definition is
provided below.

###### Definition 3.1

Let us consider an input-driven system of equation $x[t]=G(u[t],x[t-1])$, e.g.
the one defined by the ES2N’s map whose state update equation is given by eq.
(5). Let be given an initial internal state condition $x[0]$ at time $t=0$,
and a sequence of $T$ inputs $u[1],\ldots,u[T]$. Then the input-driven
trajectory $\\{(u[t+1],x[t])\\}_{t=0,\ldots,T-1}$ is well defined via eq. (5).
Then, the MLLE of the ES2N on such input-driven trajectory is defined as
follows

$\Lambda:=\max_{n=1,\ldots,N_{r}}\dfrac{1}{T}\sum_{t=0}^{T-1}\log(r_{n}[t])$
(16)

where $N_{r}$ is the dimension of the reservoir matrix $\mathbf{W}_{r}$, and
$r_{n}[t]$ is the square root of the modulus of the $n$th eigenvalue of the
symmetric real matrix $\mathbf{J}[t]\mathbf{J}[t]^{T}$, with $\mathbf{J}[t]$
defined as in eq. (10).

Roughly speaking, eq. (16) gives us an estimation of the maximum expansion
rate, averaged over the window of time $[0,T]$, locally to a given input-
driven trajectory. On the same vein of autonomous dynamical systems, a
$\Lambda<0$ denotes local contractive stable dynamics, $\Lambda>0$ is a
blueprint of chaotic dynamics since it implies local exponential divergence of
trajectories, while $\Lambda=0$ characterises the edge of chaos.
Below we provide an estimation of the MLLE for the ES2N model.

###### Theorem 3.5

Let us consider an ES2N model whose state update equation is given by eq. (5).
Then, for all time lenghts $T$, initial internal state $x[0]$, and inputs
$u[1],\ldots,u[T]$, the MLLE defined by eq. (16) is bounded as follows

$\log\bigl{(}1-\beta(\gamma\sigma+1)\bigl{)}\leq\Lambda\leq\log\bigl{(}1+\beta(\gamma\sigma-1)\bigl{)},$
(17)

In particular, in the first order approximation of small values of $\beta$ it
holds that

$\Lambda\approx-\beta.$ (18)

Proof. Thanks to Theorem 3.3 we know that the modulus of each eigenvalue $\mu$
of the Jacobian $\mathbf{J}[t]$ is bounded as
$1-\beta-\beta\gamma\sigma\leq\,\,|\mu|\,\,\leq 1-\beta+\beta\gamma\sigma$,
regardless of $t$. Each eigenvalue $\nu$ of $\mathbf{J}[t]\mathbf{J}[t]^{T}$
is bounded as
$(1-\beta-\beta\gamma\sigma)^{2}\leq\,\,|\nu|\,\,\leq(1-\beta+\beta\gamma\sigma)^{2}$,
regardless of $t$. In particular, the square root $r_{n}[t]$ of the modulus of
the $n$th eigenvalue of the matrix $\mathbf{J}[t]\mathbf{J}[t]^{T}$ is bounded
as $1-\beta-\beta\gamma\sigma\leq\,\,r_{n}[t]\,\,\leq
1-\beta+\beta\gamma\sigma$, regardless of $n\in\\{1,\ldots,N_{r}\\}$, and $t$.
Therefore, it follows (17) from Definition 3.1. In particular, for small
values of $\beta\approx 0$ we have a tight squeeze that justifies the
estimation of $\Lambda$ as the arithmetic mean of the first order
approximation of the left bound,
$\log\bigl{(}1-\beta(\gamma\sigma+1)\bigl{)}\approx-\beta(\gamma\sigma+1)$,
and the right bound,
$\log\bigl{(}1+\beta(\gamma\sigma-1)\bigl{)}\approx\beta(\gamma\sigma-1)$,
which results in the approximation $\Lambda\approx-\beta$.

###### Remark 3.3

Theorem 3.5 implies that, in an ES2N, we can tune the recurrent neural
dynamics towards the edge of chaos via tuning the proximity hyperparameter
$\beta$, regardless of the input. More precisely, for decreasing values of
$\beta$ the bounds of eq. (17) become tighter, and in the approximation of
small values of $\beta$ the ES2N model gets close to the edge of chaos
($\Lambda\approx 0$) while being on average over time in the stable regime
characterised by $\Lambda<0$.

##  4 Experiments

In this section we present our experimental analysis on the ES2N model, in
comparison with well-established approaches from the ESN literature.
Specifically, the short-term memory capacity is analysed in Section 4.1, while
in Section 4.2 we investigate the trade-off between nonlinearity and memory.
Finally, in Section 4.3 we focus on the autoregressive time series modeling,
using the multiple superimposed oscillators case as a reference task.

###  4.1 Memory Capacity

The Memory Capacity (MC) task was introduced by Jaeger in [14] to measure the
short-term memory ability of an ESN to retrieve from the pool of internal
neuronal activations the past input history. The task later became very
popular for analysing the computational properties of RC-based models [28]. A
systematic analysis of MC varying various hyperparameters of ESNs can be found
in [7], in which the authors also propose gradient descent based
orthogonalization procedures to increase the MC. In the following, we set an
MC experiment similarly to [10].

We consider reservoirs of $N_{r}=100$ neurons, with a linear readout trained
by ridge regression. The input $u[t]$ is an i.i.d. signal uniform in
$[-0.8,0.8]$ of discrete-time length $T=6000$. The first 5000 time steps are
exploited for training (excluding the very first 100 time steps to warm up the
reservoir), and the last 1000 time steps are left for test. The task is to
reproduce in output a signal $z_{k}[t]$ that is a delayed version of $k$ time
steps of the input signal, i.e., to have $z_{k}[t]$ as close as possible to
$u[t-k]$. The MC score is defined as follows:

$MC=\sum_{k=1}^{\infty}MC_{k},$ (19)

where $MC_{k}$ is the squared correlation coefficient between the output
$z_{k}[t]$ and the target $u[t-k]$, defined as follows:

$MC_{k}=\dfrac{\Bigl{(}<\bigl{(}z_{k}[t]-<z_{k}[t]>_{t}\bigl{)}\bigl{(}u[t-k]-<u[t-k]>_{t}\bigl{)}>_{t}\Bigl{)}^{2}}{<\Bigl{(}z_{k}[t]-<z_{k}[t]>_{t}\Bigl{)}^{2}>_{t}\,\,<\Bigl{(}u[t-k]-<u[t-k]>_{t}\Bigl{)}^{2}>_{t}}.$
(20)

Angular bracket in eq. (20) denotes average over time, and are calculated with
regard to the test session, i.e. for $t=5001,\ldots,6000$. Moreover, the
calculation of the sum in eq. (19) has been truncated to $k=200$, i.e. twice
the reservoir size. This choice of truncating the sum in eq. (19), widely used
in the RC community, makes sense considering that the maximum MC achievable by
an $N_{r}$-dimensional reservoir is $N_{r}$; for a proof of this fact see [14,
Proposition 2].

For the calculations, both ESNs and ES2N have been set with a spectral radius
of $0.9$, and an input scaling of $0.1$. This setting has been tested as good
for ESNs in previous works [10, 31]. Keeping fixed those hyperparameters, the
leak rate $\alpha$ for the ESN model, and the proximity hyperparameter $\beta$
for the ES2N model, have been varied in $(10^{-3},1)$. Precisely, we used the
same grid of 50 values for $\alpha$ and $\beta$, generated via the formula
$a10^{-s}$, with $a$ random uniform in $(0.1,1)$, and $s$ random uniform in
$\\{0,1,2\\}$. This methodology, while covering a large range of values for
the leak rate of ESNs, also ensures to properly explore values close to zero
where the ES2N model presents an interesting behaviour.

Figure 2: Left: MC values of eq. (19) averaged over 10 trials for various
values of $\alpha$ (for the ESN), and $\beta$ (for the ES2N). Right: Output
signals of ESN with $\alpha=1$ (green), and ES2N with $\beta=0.05$ (red), over
the ground truth signal (dashed), for the two cases of delay $k=30$ and
$k=98$.

Both ESNs and ES2Ns have been run for 10 different initialisations for each
delay $k$, and the computed MC has been averaged over these trials in order to
have statistical significance. In the left picture of Figure 2 the computed MC
values for both ESN and ES2N are plotted. The MC of ES2N exhibits a peculiar
nonlinear dependence on the proximity hyperparameter $\beta$, peaking around
$\beta=0.05$. On the contrary, the MC of ESN appears monotonic in $\alpha$,
reaching the highest value at $\alpha=1$. As expected, for $\beta$ values
approaching to 1, the ES2N’s MC and the ESN’s MC curves overlap on the right
part of the graph. In the right plots of Figure 2, two examples with delay
$k=30$ and $k=98$ are reported for leaky ESN and ES2N (each with its best
hyperparameter setting found) highlighting the supremacy of the ES2N model
over a standard ESN. Note that recalling the input signal up to 98 time steps
in the past is challenging for a reservoir of 100 units. The computed mean MC
values for the optimal $\alpha$ and $\beta$ are reported in Table 4 along with
their empirical standard deviations.

Additionally, we computed the MC of linear ESN, i.e. the model of eq. (1) with
$\alpha=1$ and the identity function as $\phi$ (called linearESN in Table 4);
the MC of an ESN with an orthogonal reservoir, i.e. the model of eq. (1) with
$\alpha=1$ and a randomly generated orthogonal matrix $\mathbf{W}_{r}$ (called
orthoESN in Table 4); and the MC of a linear ESN with a specific orthogonal
structure that realises a circular shift (called linearSCR in Table 4), i.e.
the model of eq. (1) with $\alpha=1$ and $\mathbf{W}_{r}$ with nonzero
elements in the lower subdiagonal and the upper-right corner, all filled with
$1$.111The linearSCR model implements the same reservoir topology used in the
_Simple Cycle Reservoir_ in [28], from which the name linearSCR. All the other
hyperparameters of linearESN, orthoESN, and linearSCR have been set
identically to those specified previously. Linear ESNs are known to perform
better than nonlinear ESNs in the MC task. However, as evident from the large
standard deviation of linearESN in Table 4, linear models occasionally give
very poor performance. One way to stabilise linear ESN’s performance is to
employ an orthogonal matrix as reservoir. We used linearSCR as benchmark
because theoretical results (see [14, Proposition 4] and [28, Theorem 1])
guarantee for it optimal performance on the MC task. Remarkably, the ES2N
model can get very close to the MC optimal value (of 100, for a reservoir of
100 neurons) with a noticeably narrow standard deviation, almost matching the
MC of linearSCR, see Table 4.

Finally, we plot in the left picture of Figure 4 the squared correlation
coefficient $MC_{k}$ for all $k=1,\ldots,200$, for all of the five considered
models: leaky ESN (with $\alpha=1$, the best found), linearESN, orthoESN,
linearSCR, and ES2N (with $\beta=0.05$, the best found). These $MC_{k}$ values
are not averaged over more trials. However, apart from linearESN (which
sometimes fails), the $MC_{k}$ values of all the other four models are
relatively insensitive to the random initialisation. As evident from the left
plot of Figure 4, ES2N presents an MC curve (red) particularly close to the
optimal one of linearSCR (black), while we note that despite the large MC of
the orthoESN model, the orthoESN never excels in reconstructing the delayed
version of the input, even for the simplest case of $k=1$. In the right plot
of Figure 4 a comparison between orthoESN and ES2N outputs is plotted, for the
case of delay $k=1$. The ES2N (red) is able to perfectly reconstruct the
target (black dashed line), while the orthoESN (blue) seems to struggle.
Interestingly, the orthoESN trades its poor reconstruction of the delayed
signal with the ability to mildly correlate its output with the target for
very large delays, and given the definition of $MC=\sum_{k}MC_{k}$, this
results in an overall large value of memory capacity.

Model | MC
---|---
leaky ESN | $30.40\pm 3.76$
linearESN | $49.35\pm 17.13$
orthoESN | $89.42\pm 1.50$
linearSCR | $99.09\pm 0.01$
ES2N | $98.43\pm 0.11$

Figure 3: Mean and standard deviation of the MC computed over 10 different
initialisations of reservoir models of 100 neurons. Leaky ESN is with
$\alpha=1$, ES2N is with $\beta=0.05$

Figure 4: Left: $MC_{k}$ values for delays $k=1,\ldots,200$, for all the 5
models.
Right: orthoESN versus ES2N in the simplest case of $k=1$.

###  4.2 Memory-nonlinearity trade-off

From previous RC literature, it is well known the existence of a trade-off
between the strength of nonlinearity of a dynamical reservoir system and its
short-term memory abilities [6]. In [12] the authors propose a task with the
aim of measuring this trade-off. The task consists of extracting from an
i.i.d. uniform input signal $u[t]$ in $[-1,1]$ a target signal of the form
$y[t]=\sin(\nu*u[t-\tau])$; here $\nu$ quantifies the _nonlinearity strength_
, while $\tau$ measures the _memory depth_. We use this task to benchmark the
memory-nonlinearity trade-off for various combinations of $\tau$ and $\nu$
comparing ES2N against leaky ESN, and linearSCR. For this experiment, an input
signal of length 6000 has been generated, of which the first 5000 for training
(excluding the very first 100 steps), and the remaining 1000 for testing. The
metric used to evaluate the task is the NRMSE between target $y(t)$ and output
$z(t)$ in the test session, defined as follows:

$NRMSE(y,z)=\sqrt{\dfrac{<\lVert y[t]-z[t]\rVert^{2}>_{t}}{<\lVert
y[t]-<y[t]>_{t}\rVert^{2}>_{t}}},$ (21)

where $\lVert v\rVert$ represents the Euclidean norm of a vector $v$, and
angular brackets means the average over time. According to this metric, the
lower the better.

In the following experiment, we considered a grid of values of
$(\log(\nu),\tau)$, with $\tau\in[1,20]$, with a step of $1$, and
$\log(\nu)\in[-1.6,1.6]$ with a step of $0.1$. For each pair
$(\log(\nu),\tau)$, we ran 100 instantiations of leaky ESN, linearSCR, and
ES2N. For each run we generated randomly the following hyperparameters

  * 1.

uniformly random input scaling in $[0.2,6]$,

  * 2.

uniformly random spectral radii in $[0.1,3]$,

  * 3.

$\alpha$ (for leaky ESN) and $\beta$ (for ES2N) generated via the formula
$a10^{-s}$, with $a$ uniformly random in $(0.1,1)$, and $s$ uniformly random
in $\\{0,1\\}$, so that they vary in $(10^{-2},1)$.

In Figure 5 are plotted the best NRMSE found on test on a coloured scale from
black ($NRMSE=0)$ to yellow ($NRMSE=1$ or greater), the lower the better.
Results in Figure 5 show that ES2N significantly outperforms both leaky ESN
and linearSCR. Note that, linearSCR starts to increasingly underperform as
soon as $\log(\nu)>0$, i.e. in the region where nonlinearity is needed, while
ES2N is able to retrieve the information for much stronger nonlinearly
transformed input signals. This indicates that the ES2N model can truly
exploit nonlinearity in the computation. On the other hand, leaky ESN is able
to reconstruct the input signal in the strong nonlinear regime, but only for
very small delays. In particular, for the challenging case of $\log(\nu)>1$
(strong nonlinearity), leaky ESN’s performance significantly degrades already
at $\tau=4$, with NRMSE values always above $0.5$. On the contrary, the ES2N
model obtains NRMSE values always below of $0.5$ up to delays of $\tau=16$ in
the strong nonlinearity regime of $\log(\nu)>1$.
In Figure 6 are plotted the output signals in the test session for the best
hyperparameters found on the three models ES2N, linearSCR, and leaky ESN, in
comparison with the ground truth target signal for the case of strong
nonlinearity with $\log(\nu)=1.3$ and $\tau=10$. These results demonstrate how
ES2N can conciliate the two contrasting properties of having a large memory
capacity and the ability to perform nonlinear computations.

Figure 5: Results of the memory-nonlinearity trade-off task explained in
Section 4.2. NRMSE values ranging from black (NRMSE$=0$) to yellow (NRMSE$=1$
or greater) are plotted for various combinations of delay $\tau$ and
nonlinearity strength $\nu$. Left: best NRMSE values after a random search of
100 trials on a leaky ESN. Centre: best NRMSE values after a random search of
100 trials on a linearSCR. Right: best NRMSE values after a random search of
100 trials on our proposed ES2N model. Figure 6: Input $u[t]$ is i.i.d.
uniformly sampled in $[-1,1]$. Target is a strong nonlinear transformation of
a 10-delayed version of the input given by the function $\sin(\nu u[t-10])$
with $\nu\approx 3.67$ represented in dashed line. The best hyperparameters
found for the setting of $\log(\nu)=1.3$ and $\tau=10$ are the following in
the form of (input scaling, spectral radius, $\alpha$ or $\beta$):
($0.43,1.03,0.98$) for leaky ESN, ($0.35,0.49,1.0$) for linearSCR, and
($1.44,0.10,0.18$) for ES2N.

###  4.3 Multiple superimposed oscillators in auto-regressive mode

The Multiple Superimposed Oscillators pattern generation (MSO) is a popular
benchmark task where an RNN is trained to generate autonomously (i.e. without
the input driving the dynamics) a one-dimensional signal made from the
superpositions of a few incommensurate sines. This is achieved providing in
input the target signal during the training phase, then, once the readout
matrix is trained, closing the loop in the testing phase via $u[t]=z[t]$ in
eq. (5), i.e. self-driving the ESN dynamics with its own generated output. In
these auto-regressive tasks, the ridge regression training of the readout
matrix can be regarded as a _teacher forcing_ training strategy. In the
literature [37, 30, 11, 29, 18] the MSO task with different numbers of sine
waves has been inspected. In all of the mentioned studies the frequencies of
the sine waves were taken from the same set:
$\nu_{1}=0.2,\nu_{2}=0.311,\nu_{3}=0.42,\nu_{4}=0.51,\nu_{5}=0.63,\nu_{6}=0.74,\nu_{7}=0.85$,
and $\nu_{8}=0.97$. In this section, we only consider the case of 8
frequencies, the most challenging one among those mentioned above, which we
will denote concisely with MSO8. Thus, the target signal takes the following
form:

$y[t]=\sum_{i=1}^{8}\sin(\nu_{i}t).$ (22)

The MSO8 signal of eq. (22) is then rescaled in order to have zero mean and be
bounded in $(-1,1)$. From now on, we will refer to $y[t]$ as the normalised
signal. The function of eq. (22) has an almost-period of $\tau_{P}=6283$ time
steps.222More precisely, for the normalised signal, it turns out that
$<\bigl{|}y[t+\tau_{P}]-y[t]\bigl{|}>_{t=0,\ldots,\tau_{P}}\approx 0.024\pm
0.02$. We use 6383 time steps for training, i.e. a whole “period” of 6283
steps excluding the first 100 time steps to wash out transient dynamics. In
this phase, the target signal (the teacher) is injected into the reservoir as
the input signal, i.e. $u[t]=y[t]$ for $t=1,\ldots,6383$. During the training
phase, a small Gaussian noise of zero mean and standard deviation of $10^{-4}$
has been added to the argument of the $\tanh$ in order stabilise the dynamics.
Therefore, the ES2N’ state-update equation during training was

$x[t]=\beta\phi(\mathbf{W}_{r}x[t-1]+\mathbf{W}_{in}u[t]+\eta[t])+(1-\beta)\mathbf{O}x[t-1],$
(23)

with $\eta\in N(0,10^{-4})$ the source of Gaussian noise. We set zero
regularisation, i.e. $\mu=0$ in eq. (3). Thus, we train the linear readout to
reproduce the target signal. Then, we close the loop feeding back the output
into the reservoir in place of the external input, i.e. $u[t]=z[t]$ for $t\geq
6384$. From this moment on, the RNN runs autonomously (without noise).

####  4.3.1 Random hyperparameter search

First of all, we accomplished a random search to optimise the hyperparameters
of both leaky ESN and ES2N. Varying the reservoir size revealed that, as
expected, larger models leads to better results, for both leaky ESNs and
ES2Ns. We considered leaky ESNs of 600 reservoir neurons, and ES2Ns of 100
reservoir neurons. We considered leaky ESNs with reservoir size 6 times larger
of the ES2N’s, in order for the leaky ESN to get competitive performance on
the challenging MSO8 task. We ran 10000 different initialisations for both
leaky ESN and ES2N with uniformly random generated hyperparameters as follows:

  * 1.

uniformly random spectral radius $\rho\in[0.8,1.2]$

  * 2.

uniformly random input scaling $\omega_{in}\in[0,0.4]$

  * 3.

(for leaky ESN) uniformly random $\alpha\in(0.1,1)$

  * 4.

(for ES2N) uniformly random $\beta\in(0.01,0.1)$.

We also tried to vary the reservoir connectivity, but it did not influence the
performance. Thus, we used fully connected reservoirs. In the training phase,
a small Gaussian noise with zero mean and standard deviation of $10^{-4}$ has
been introduced in the state-update equation as in eq. (23), for both leaky
ESN and ES2N. To evaluate the performance we compute the NRMSE as defined in
eq. (21) for 300 time steps in the testing phase, i.e. for
$t=6384,\ldots,6684$. We refer the reader to Appendix A for the visualisation
of the data obtained from this random search. In summary, from this search
emerged that leaky ESNs need spectral radii strictly around 1 for good results
on the MSO8 task. On the contrary, ES2Ns are quite insensitive on the choice
of the spectral radius. However, by design ES2N reflects the insensitivity on
the spectral radius on a dependence on the proximity hyperparameter $\beta$,
which for this task reaches its optimum around $\beta=0.03$. Apart from that,
the most influential hyperparameter turned out to be the input scaling, in
line with previous works. In particular, we found the optimal combination of
$\rho=0.99$, $\omega_{in}=0.05$, and $\alpha=0.9$ for leaky ESNs. While for
ES2N $\omega_{in}=0.11$, and $\beta=0.03$ (regardless of $\rho$). All the
NRMSEs computed in this random search are reported in the histograms in the
left plots of Figure 7. The mean and standard deviation of NRMSE for leaky ESN
is $1.44\pm 13.34$, while for the ES2N is $0.05\pm 0.11$. The difference of
mean NRMSE values highlights the supremacy of ES2N over the leaky ESN on the
MSO8 task, despite a 6 times smaller reservoir size. As a side result, this
random search revealed how ES2Ns are characterised by a wider “good”
hyperparameter region compared to leaky ESNs, testified via the more than 120
times larger standard deviation of leaky ESN over ES2N.333There are few
hundreds NRMSE values of leaky ESN greater than 2 which do not appear in the
histogram in Figure 7, some of them exceeding NRMSE of 100.

####  4.3.2 Stability in the long run

In this section, we compare the quality of the learned signal in the long run.
For the purpose, we fixed the hyperparameters in their respective optimal
setting, precisely we set $(\rho=1,\omega_{in}=0.11,\beta=0.03)$ for ES2Ns,
and $(\rho=0.99,\omega_{in}=0.05,\alpha=0.9)$ for leaky ESNs.444The choice of
$\rho=1$ for ES2N was arbitrary since it does not influence the outcome of
this experiment, see also Appendix A. Therefore, we trained a large leaky ESN
of 3000 neurons, and a relatively small ES2N of 300 neurons, on 6383 training
steps as explained above (i.e. 100 for transient, and 6283 for actual
training). The resulting output signals in the test session of these trained
models are reported in the centre plots of Figure 7. In the beginning both
models follow tightly the target (dashed line). After 15000 time steps the
leaky ESN output (green) already deviates significantly from the target, while
the ES2N output (red) still performs very well. After 20000 time steps the
leaky ESN output is completely decorrelated with the target. On the contrary,
the ES2N manages to generate a meaningful output signal even after 50000 time
steps of autonomous run. Remarkably, the ES2N model is able to substantially
outperform the leaky ESN model while having a number of neurons which is an
order of magnitude smaller.

Figure 7: MSO task with 8 frequencies. In all plots the colour green
corresponds to the leaky ESN model, while red to the ES2N model. Left:
hystograms of NRMSE of 10000 initialisations of leaky ESN (600 reservoir
neurons) and ES2N (100 reservoir neurons) with hyperparameters uniformly
generated as explained in Section 4.3.1. In the upper plot, there are almost
300 cases with NRMSE over 2 which have been cut off from the picture (the
maximum reaching a NRMSE of 942). In the bottom plot both leaky ESN and ES2N’s
histograms are plotted together and magnified around low NRMSE values of
$0.05$. Centre and right: A large fine tuned leaky ESN of 3000 neurons and a
relatively small ES2N of 300 neurons have been trained to reproduce with the
feedback of the output the MSO8 signal (dashed line). The output signals
generated by leaky ESN (green) and ES2N (red) after various time intervals
$\Delta t$ of running in auto-generation mode are plotted, from top to bottom,
$\Delta t=0,15000,20000,50000$. On the right plots, the eigenvalues of the
Jacobians of the corresponding trained leaky ESN model, and ES2N model, on a
randomly selected time step in the testing session.

As explained in [16], the MSO is a relatively easy task for a linear ESN, as
long as one has at least two reservoir neurons per frequency. However, the
apparent perfection of trained linear ESNs hides an intrinsic unstable phase-
coupling dynamics, which manifests itself as soon as little perturbations are
applied to the system. On the other hand, nonlinear ESNs are more resilient to
perturbations, but they struggle to learn functions composed of even just two
superimposed oscillators, especially if one wishes to maintain the learned
oscillations for long time. The ES2N is a nonlinear model lying in between
these two extremes, since it can be trained to self-sustain complex
oscillatory dynamics for very long time, as shown in the centre plots of
Figure 7. This can be attributed to the peculiar form of the ES2N’s Jacobian
spectrum which tends to dispose its eigenvalues along the unitary circle
promoting the dynamics to take place at the edge of stability as evident from
the right plots in Figure 7.

##  5 Conclusions

In this paper, we have developed a new RC architecture known as the Edge of
Stability Echo State Network (ES2N), which has the unique feature of being
controllably tunable to the edge of chaos dynamical behavior. Our mathematical
analysis first showed that the proposed model has dynamic behavior that is
able to exhibit contracting dynamics and the Echo State Property in a manner
similar to a standard ESN. Furthermore, and relevantly, we provided precise
analytical bounds for the entire eigenspectrum of the Jacobian of the forward
map that hold on each input-driven trajectory and for all time steps. As a
result, the reservoir exhibits a behavior whose quality is determined by a
specific proximity hyper-parameter. By architectural construction, smaller
values of this hyper-parameter result in dynamics that are progressively
closer to the edge of chaos. Overall, the introduced model takes advantage of
both the benefits of having a linear orthogonal reservoir and a nonlinear
contracting dynamics.

We empirically showed that ES2N can reach the maximum memory capacity
obtainable within a given reservoir size, showing significant advantages over
standard nonlinear and linear ESN alternatives. Furthermore, we tested the
trade-off between nonlinear computation and long short-term memory, and found
that ES2N can reconstruct strongly nonlinear transformations of relatively
large delayed input signals, where both standard (nonlinear) ESN and linear
ESN fail. Finally, we empirically demonstrated the superiority of ES2N in the
generation of complex oscillatory patterns. Remarkably, the recurrent network
driven with its own output signal can produce meaningful oscillations for 50
thousands time steps of autonomous run on the MSO8 task.

The analytical and experimental results presented in this paper have shown,
already in this form, the advantages of combining nonlinear dynamics and
orthogonal transformations in the state space of a recurrent network. Building
upon these findings, our future research aims to investigate alternative forms
of reservoir construction in a ES2N. For instance, we plan to explore
reservoir construction methods based on permutation or circular shift
matrices, which on the one hand allow to further improve the computational
efficiency of the approach while maintaining its computational properties, and
on the other hand are prone to implementation in edge or neuromorphic devices.
Moreover, in forthcoming studies, we will explore the performance of ES2N in
various applications, such as time series forecasting, attractor
reconstruction, and classification tasks.

## Acknowledgements

This work is partially supported by the EC H2020 programme under project
TEACHING (grant n. 871385), and by the EU Horizon research and innovation
programme under project EMERGE (grant n. 101070918).

## References

  * Bauer and Fike [1960] F. L. Bauer and C. T. Fike.  Norms and exclusion theorems.  _Numerische Mathematik_ , 2(1):137–141, 1960\.
  * Bengio et al. [1994] Y. Bengio, P. Simard, and P. Frasconi.  Learning long-term dependencies with gradient descent is difficult.  _IEEE transactions on neural networks_ , 5(2):157–166, 1994.
  * Bertschinger and Natschläger [2004] N. Bertschinger and T. Natschläger.  Real-time computation at the edge of chaos in recurrent neural networks.  _Neural computation_ , 16(7):1413–1436, 2004\.
  * Boedecker et al. [2012] J. Boedecker, O. Obst, J. T. Lizier, N. M. Mayer, and M. Asada.  Information processing in echo state networks at the edge of chaos.  _Theory in Biosciences_ , 131:205–213, 2012.
  * Ceni et al. [2020] A. Ceni, P. Ashwin, L. Livi, and C. Postlethwaite.  The echo index and multistability in input-driven recurrent neural networks.  _Physica D: Nonlinear Phenomena_ , 412:132609, 2020.
  * Dambre et al. [2012] J. Dambre, D. Verstraeten, B. Schrauwen, and S. Massar.  Information processing capacity of dynamical systems.  _Scientific reports_ , 2(1):1–7, 2012.
  * Farkaš et al. [2016] I. Farkaš, R. Bosák, and P. Gergel’.  Computational analysis of memory capacity in echo state networks.  _Neural Networks_ , 83:109–120, 2016.
  * Gallicchio [2018] C. Gallicchio.  Chasing the echo state property.  _arXiv preprint arXiv:1811.10892_ , 2018.
  * Gallicchio and Micheli [2011] C. Gallicchio and A. Micheli.  Architectural and markovian factors of echo state networks.  _Neural Networks_ , 24(5):440–456, 2011.
  * Gallicchio et al. [2017] C. Gallicchio, A. Micheli, and L. Pedrelli.  Deep reservoir computing: A critical experimental analysis.  _Neurocomputing_ , 268:87–99, 2017.
  * Holzmann and Hauser [2010] G. Holzmann and H. Hauser.  Echo state networks with filter neurons and a delay&sum readout.  _Neural Networks_ , 23(2):244–256, 2010.
  * Inubushi and Yoshimura [2017] M. Inubushi and K. Yoshimura.  Reservoir computing beyond memory-nonlinearity trade-off.  _Scientific reports_ , 7(1):1–10, 2017.
  * Jaeger [2001] H. Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  _Bonn, Germany: German National Research Center for Information Technology GMD Technical Report_ , 148(34):13, 2001.
  * Jaeger [2002] H. Jaeger.  Short term memory in echo state networks. gmd-report 152.  In _GMD-German National Research Institute for Computer Science (2002), http://www. faculty. jacobs-university. de/hjaeger/pubs/STMEchoStatesTechRep. pdf_. Citeseer, 2002.
  * Jaeger and Haas [2004] H. Jaeger and H. Haas.  Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication.  _science_ , 304(5667):78–80, 2004.
  * Jaeger et al. [2007] H. Jaeger, M. Lukoševičius, D. Popovici, and U. Siewert.  Optimization and applications of echo state networks with leaky-integrator neurons.  _Neural networks_ , 20(3):335–352, 2007.
  * Kolen and Kremer [2001] J. F. Kolen and S. C. Kremer.  _A field guide to dynamical recurrent networks_.  John Wiley & Sons, 2001.
  * Koryakin et al. [2012] D. Koryakin, J. Lohmann, and M. V. Butz.  Balanced echo state networks.  _Neural Networks_ , 36:35–45, 2012.
  * Langton [1990] C. G. Langton.  Computation at the edge of chaos: Phase transitions and emergent computation.  _Physica D: nonlinear phenomena_ , 42(1-3):12–37, 1990.
  * Legenstein and Maass [2007] R. Legenstein and W. Maass.  What makes a dynamical system computationally powerful.  _New directions in statistical signal processing: From systems to brain_ , pages 127–154, 2007.
  * Livi et al. [2017] L. Livi, F. M. Bianchi, and C. Alippi.  Determination of the edge of criticality in echo state networks through fisher information maximization.  _IEEE transactions on neural networks and learning systems_ , 29(3):706–717, 2017.
  * Lukoševičius [2012] M. Lukoševičius.  A practical guide to applying echo state networks.  In _Neural networks: Tricks of the trade_ , pages 659–686. Springer, 2012.
  * Lukoševičius and Jaeger [2009] M. Lukoševičius and H. Jaeger.  Reservoir computing approaches to recurrent neural network training.  _Computer Science Review_ , 3(3):127–149, 2009\.
  * Manjunath and Jaeger [2013] G. Manjunath and H. Jaeger.  Echo state property linked to an input: Exploring a fundamental characteristic of recurrent neural networks.  _Neural computation_ , 25(3):671–696, 2013.
  * Meckes [2021] E. Meckes.  The eigenvalues of random matrices.  _arXiv preprint arXiv:2101.02928_ , 2021.
  * Ott [2002] E. Ott.  _Chaos in dynamical systems_.  Cambridge university press, 2002.
  * Packard [1988] N. H. Packard.  Adaptation toward the edge of chaos.  _Dynamic patterns in complex systems_ , 212:293–301, 1988\.
  * Rodan and Tino [2010] A. Rodan and P. Tino.  Minimum complexity echo state network.  _IEEE transactions on neural networks_ , 22(1):131–144, 2010.
  * Roeschies and Igel [2010] B. Roeschies and C. Igel.  Structure optimization of reservoir networks.  _Logic Journal of IGPL_ , 18(5):635–669, 2010\.
  * Schmidhuber et al. [2007] J. Schmidhuber, D. Wierstra, M. Gagliolo, and F. Gomez.  Training recurrent networks by evolino.  _Neural computation_ , 19(3):757–779, 2007.
  * Schrauwen et al. [2008] B. Schrauwen, M. Wardermann, D. Verstraeten, J. J. Steil, and D. Stroobandt.  Improving reservoirs using intrinsic plasticity.  _Neurocomputing_ , 71(7-9):1159–1171, 2008.
  * Strogatz [2018] S. H. Strogatz.  _Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering_.  CRC press, 2018.
  * Verstraeten et al. [2007] D. Verstraeten, B. Schrauwen, M. d’Haene, and D. Stroobandt.  An experimental unification of reservoir computing methods.  _Neural networks_ , 20(3):391–403, 2007.
  * Vulpiani et al. [2009] A. Vulpiani, F. Cecconi, and M. Cencini.  _Chaos: from simple models to complex systems_ , volume 17.  World Scientific, 2009.
  * Wainrib and Galtier [2016] G. Wainrib and M. N. Galtier.  A local echo state property through the largest lyapunov exponent.  _Neural Networks_ , 76:39–45, 2016.
  * Werbos [1990] P. J. Werbos.  Backpropagation through time: what it does and how to do it.  _Proceedings of the IEEE_ , 78(10):1550–1560, 1990.
  * Xue et al. [2007] Y. Xue, L. Yang, and S. Haykin.  Decoupled echo state networks with lateral inhibition.  _Neural Networks_ , 20(3):365–376, 2007.
  * Yildiz et al. [2012] I. B. Yildiz, H. Jaeger, and S. J. Kiebel.  Re-visiting the echo state property.  _Neural networks_ , 35:1–9, 2012.

##  Appendix A Data from the random search of the MSO8 task

Here we report the NRMSE values of the random search of the 10000 runs for the
MSO8 task for both the ES2N model (left plots), and the leaky ESN model (right
plots). Note that a great portion of the 10000 runs of the ES2N model has
NRMSE value less than 2e-2, while the vast majority of runs of the leaky ESN
model has NRMSE value greater than 1e-1.

Figure 8: NRMSEs resulted from the random search of section 4.3.1 plotted
versus each one of the hyperparameter. The horizontal orange lines correspond
to the minimum NRMSE values reached, while the corresponding minimum points
are highlighted in red on the abscissa. Left: data of ES2N. Right: data of
leaky ESN.