Brain-inspired photonic signal processor for periodic
pattern generation and chaotic system emulation
Piotr Antonik,1, ∗ Marc Haelterman,2 and Serge Massar1

arXiv:1802.02026v1 [] 6 Feb 2018

1

Laboratoire d’Information Quantique, Université libre de Bruxelles,
50 Avenue F. D. Roosevelt, CP 224, B-1050 Brussels, Belgium
2
Service OPERA-Photonique, Université libre de Bruxelles,
50 Avenue F. D. Roosevelt, CP 194/5, B-1050 Brussels, Belgium

Reservoir computing is a bio-inspired computing paradigm for processing time-dependent signals.
Its hardware implementations have received much attention because of their simplicity and remarkable performance on a series of benchmark tasks. In previous experiments the output was uncoupled
from the system and in most cases simply computed offline on a post-processing computer. However,
numerical investigations have shown that feeding the output back into the reservoir would open the
possibility of long-horizon time series forecasting. Here we present a photonic reservoir computer
with output feedback, and demonstrate its capacity to generate periodic time series and to emulate
chaotic systems. We study in detail the effect of experimental noise on system performance. In
the case of chaotic systems, this leads us to introduce several metrics, based on standard signal
processing techniques, to evaluate the quality of the emulation. Our work significantly enlarges
the range of tasks that can be solved by hardware reservoir computers, and therefore the range of
applications they could potentially tackle. It also raises novel questions in nonlinear dynamics and
chaos theory.
I.

INTRODUCTION

Reservoir Computing (RC) is a set of machine learning
methods for designing and training artificial neural networks, introduced independently in [1] and in [2]. The
idea behind these techniques is that one can exploit the
dynamics of a recurrent nonlinear network to process
time series without training the network itself, but simply adding a general linear readout layer and only training the latter. This results in a system that is significantly easier to train (the learning is reduced to solving
a system of linear equations [3]), yet powerful enough to
match other algorithms on a series of benchmark tasks.
RC has been successfully applied to, for instance, channel
equalisation [1], phoneme recognition [4] and won an international competition on prediction of future evolution
of financial time series [5].
Reservoir Computing allows to efficiently implement
simplified recurrent neural networks in hardware, such
as e.g. optical components. Optical computing has been
investigated for decades as photons propagate faster than
electrons, without generating heat or magnetic interference, and thus promise higher bandwidth than conventional computers [6]. RC would thus allow to build highspeed and energy efficient photonic computing devices.
Several important steps have been taken towards this
goal with electronic [7], opto-electronic [8–11], all-optical
[12–14] and integrated [15] experimental RC implementation reported since 2012.
Forecasting is one of the central problems in science:
how can we predict the future from the past? Over the
past few decades, artificial neural networks have gained a

∗

piotr.antonik@ulb.ac.be

significant popularity in the time series forecasting community. Compared to previously employed statisticsbased techniques, they are both data driven and nonlinear, more flexible and do not require an explicit model
of the underlying process. A review of artificial neural
networks models for time series forecasting can be found
in [16]. Reservoir computing can be directly applied to
short-term prediction tasks, that focus on generating a
few future timesteps. As for long-horizon forecasting,
that involves predicting the time series for as long as
possible, it is possible with a small modification of the architecture, namely by feeding the RC output signal back
into the reservoir. This additional feedback significantly
enriches the internal dynamics of the system, enabling
it to generate time series autonomously, that is without
receiving any input signal. With this modification, reservoir computing can be used for long-term prediction of
chaotic series [1, 5, 17–19]. In fact, this approach holds,
to the best of our knowledge, the record for such chaotic
time series prediction [1, 5]. A reservoir computer with
output feedback can also achieve the easier task of generating periodic signals [20–22], and of producing a tunable
frequency [23–25].
The aim of the present work is to explore these novel
applications experimentally. Indeed, they have been
widely studied numerically, but no experimental implementation has been reported so far. There are multiple
motivations for this investigation. First of all, reservoir
computing is a biologically inspired algorithm. Indeed
one of the main motivations of the seminal paper [2] was
to propose how microcircuits in the neocortex could process information. More recently it has been realised that
the cerebellum has a structure very similar to that of a
reservoir computer [26, 27]. Generating time series with
specific attributes is an important property of biological
neural and chemical circuits (for e.g. movement control,

2
biological rhythms, etc). Are biological circuits similar to
reservoir computers used to generate trainable time series? Investigating this process experimentally can shed
light on this tantalising question, for instance by clarifying which kinds of time series, and what training processes are robust to experimental imperfections.
Second, generation of time series with specific properties is an important task in signal generation and processing. Given the possibility that photonic reservoir computing could carry out ultra-fast and low energy optical
signal processing, this is an important area to explore,
again with the aim of understanding which tasks are robust to experimental imperfections.
Finally, this investigation raises a new fundamental
question in nonlinear dynamics: given a system that emulates a known chaotic time series, how does one quantify
the quality of the emulation. Answering this question
becomes vital in case of experimental implementations,
as physical systems are affected by noise, and thus can
not output anything better than an approximate, noisy
emulation of the target chaotic time series. The comparison techniques previously used in numerical investigations fail in such situations, and one needs to develop
new evaluation metrics.
Experimental implementation of these ideas requires,
in principle, a fast readout layer capable of generating
and feeding back the output signal in real-time. Several analogue solutions have been under investigation recently [28–30], but none are as yet capable of realising
this application. In fact, to successfully train an analogue readout layer with offline learning methods, used
in most experimental RC setups up to now, a very precise
model of the readout setup is required, which is hardly
achievable experimentally, as shown in [29]: it is virtually
impossible to characterise each hardware component of
the setup with sufficient level of accuracy. The reason for
this sensitivity is that the output is a weighted sum with
positive and negative coefficients of the internal states of
the reservoir. Therefore errors in the coefficients quickly
build up and become comparable to the value of the desired output. For this reason, we chose the approach of
a fast real-time digital readout layer implemented on a
Field-Programmable Gate Array (FPGA) chip. The use
of high-speed dedicated electronics makes it possible to
compute the output signal in real time and feed it back
into the reservoir. In order to keep the experiment simple, we used as reservoir the opto-electronic delay system introduced in [7–9], that has shown state-of-the-art
results on several benchmark tasks and is fairly easy to
operate. The coupling of an FPGA board to an optoelectronic reservoir was already reported in [31], where
the capacity to compute the output in real time was
used to solve tasks that change in time. Note however
that the FPGA design, i.e. the program implemented on
the chip, differs greatly from the one used in the present
work. Here we use the FPGA to feed the output of the
reservoir back into itself, and thereby to generate periodic time signals and emulate chaotic systems. This is a

completely different problem.
Our experiments show that the system successfully
solves two periodic time series generation tasks: frequency and random pattern generation, that have been
previously investigated numerically [21, 23, 24]. The first
task allows to demonstrate different timescales within the
neural network, and the second can be used as a memory metric. The photonic computer manages to generate both sine waves and random patterns with unlimited
stability. Furthermore, we apply the RC to emulation
of two chaotic attractors: Mackey-Glass [32] and Lorenz
[33] systems. In the literature, the performance on these
tasks is quantified by computing the prediction horizon,
i.e. the duration for which the RC can accurately follow a given trajectory on the chaotic attractor [1]. However, this method fails in the presence of a relatively high
level of experimental noise, with a Signal-to-Noise Ratio
(SNR) of roughly 40 dB, as will be discussed in section
V A. This noise was not problematic in previous experiments using the same opto-electronic reservoir [8, 31],
but turns out to be intolerable for a system with output
feedback. This raises the question of how to evaluate a
system that emulates a known chaotic time series in the
presence of noise. In this study, we introduce several new
approaches, such as frequency spectrum comparison and
randomness tests. These approaches are based on wellknown signal analysis techniques, but they are employed
for the first time here for the evaluation of a chaotic signal generated by a reservoir computer. Our results show
that, although the RC struggles at following the target
trajectory on the chaotic attractor, its output accurately
reproduces the core characteristics of the target time series.
The paper is structured as follows. Sections II and III
introduce the basic principles of the reservoir computing
and the time series generation tasks investigated in this
work. The experimental setup, FPGA design and numerical simulations are outlined in section IV. All experimental and numerical results are presented and discussed in
section V, and section VI concludes the paper.

II.

BASIC PRINCIPLES OF RESERVOIR
COMPUTING

A typical reservoir computer contains a large number
N of internal variables xi (n) evolving in discrete time
n ∈ Z, as given by


N
−1
X
xi (n + 1) = f 
aij xj (n) + bi u(n) ,
(1)
j=0

where f is a nonlinear function, u(n) is some external
signal that is injected into the system, and aij and bi are
time-independent coefficients, drawn from some random
distribution with zero mean, that determine the dynamics of the reservoir. The variances of these distributions

3
are adjusted to obtain the best performances on the task
considered.
The nonlinear function used here is f = sin(x), as in
[8, 9]. To simplify the interconnection matrix aij , we exploit the ring topology, proposed in [7, 34], so that only
the first neighbour nodes are connected. This architecture provides performances comparable to those obtained
with complex interconnection matrices, as demonstrated
numerically in [3, 34] and experimentally in [7–9, 12, 13].
Under these conditions Eq. (1) becomes

x5

x0
Teacher input
signal u(n)

x1

Input layer

y(n) =

N
−1
X

wi xi (n),

(3)

i=0

where wi are the readout weights, trained either offline
(using standard linear regression methods, such as the
ridge regression algorithm [35] used here), or online [31],
in order to minimise the Mean Square Error (MSE) between the output signal y(n) and the target signal d(n),
given by
E
D
2
MSE = (y(n) − d(n)) .

(4)

The introduction of the output feedback requires a minor change of notations. Since the RC can now receive
two different signals as input, we shall denote I(n) the input signal, which can be either the external input signal
I(n) = u(n), or its own output, delayed by one timestep
I(n) = y(n − 1).
The reservoir computer is operated in two stages, depicted in Fig. 1: a training phase and an autonomous
run. During the training phase, the reservoir computer is
driven by a time-multiplexed teacher signal I(n) = u(n),
and the resulting states of the internal variables xi (n)
are recorded. The teacher signal depends on the task
under investigation (which will be introduced in section
III). The system is trained to predict the next value of
the teacher time series from the current one, that is, the
readout weights wi are optimised so as to get as close as
possible to y(n) = u(n + 1). Then, the reservoir input
is switched from the teacher sequence to the reservoir
output signal I(n) = y(n − 1), and the system is left running autonomously. The reservoir output y(n) is used to
evaluate the performance of the experiment.

wi xi (n)

Output
signal y(n)

x2

Output layer

Reservoir

(a)Training stage
(a)Training
stage
(a)Training stage

x4



Mi y(n − 1)

(2a)
(2b)

with i = 1, . . . , N − 1; α and β parameters that are used
to adjust the feedback and the input signals, respectively;
and Mi the input mask, drawn from a uniform distribution over the the interval [−1, +1], as in [8, 12, 34].
The reservoir computer produces an output signal
y(n), given by a linear combination of the states of its
internal variables


x3

x5

x0 (n + 1) = sin (αxN −1 (n − 1) + βM0 u(n)) ,
xi (n + 1) = sin (αxi−1 (n) + βMi u(n)) ,

x4

Mi u(n)

x0
Delayed output
y(n − 1)

x3

x1

Input layer

wi xi (n)

x2

Reservoir

Output
signal y(n)

Output layer

(b)Autonomous run
(b)Autonomous
run

Figure 1. Schematic representation of the training stage (a)
and the autonomous run (b) of our reservoir computer (here
with N = 6 nodes). During the training phase, the reservoir is
driven by a teacher input signal u(n), and the readout weights
wi are optimised for the output to be as close as possible to
u(n + 1). During the autonomous run, the teacher signal is
switched off and the reservoir is driven by its own output
signal. The readout weights wi are kept constant and the
performance of the system is measured in terms of how long
or how well it can generate the desired output.

III.

TIME SERIES GENERATION TASKS

Feeding the output back into the reservoir allows the
computer to autonomously (i.e. without any external
input) generate time series. We tested the capacity of the
experiment to generate both periodic and chaotic signals,
with two tasks in each category.

A.

Frequency generation

Frequency generation is the simplest time series generation task considered here. The system is trained to
generate a sine wave given by
u(n) = sin (νn) ,

(5)

where ν is a real-valued relative frequency. The physical
frequency f of the sine wave depends on the experimental
roundtrip time T (see section IV) as follows
f=

ν
.
2πT

(6)

This task allows to measure the bandwidth of the system
and investigate different timescales within the neural network.

4
B.

III B), during the autonomous run and use this figure
to evaluate the performance of the system.

Random pattern generation

Random pattern generation is a natural step forward
from the frequency generation task to a more complex
problem – instead of a regularly-shaped continuous function, the system is trained to generate an arbitrarilyshaped discontinuous function (that remains periodic,
though). Specifically, a pattern is a short sequence of
L randomly chosen real numbers (here within the interval [−0.5, 0.5]) that is repeated periodically to form
an infinite time series [18]. Similarly to the physical frequency in section III A, the physical period of the pattern
is given by τpattern = L · T . The aim is to obtain a stable pattern generator, that reproduces precisely the pattern and doesn’t deviate to another periodic behaviour.
To evaluate the performance of the generator, we compute the MSE between the reservoir output signal and
the target pattern signal during the training phase and
the autonomous run, with a maximal threshold set to
10−3 . This value is somewhat arbitrary, and one could
have picked a different threshold. As will be illustrated
in Figs. 11 and 15 in Sec. V, the 10−3 level corresponds
to the point where the RC strongly deviates from the
starting trajectory on the chaotic attractor. For consistency, we have used this threshold in all our experiments,
for all tasks. If the error doesn’t grow above the threshold during the autonomous run, the system is considered
to accurately generate the target pattern. We also tested
the long-term stability on several patterns by running the
system for several hours, as will be described in section
V.

D.

Lorenz chaotic series prediction

The Lorenz equations, a system of three ordinary differential equations
dx
= σ (y − x) ,
dt
dy
= −xz + rx − y,
dt
dz
= xy − bz,
dt

(8b)
(8c)

with σ, r, b > 0, was introduced as a simple model for atmospheric convection [33]. The system exhibits chaotic
behaviour for σ = 10, b = 8/3 and r = 28 [37], that
we used in this study. This yields a chaotic attractor
with the highest Lyapunov exponent of λ = 0.906 [1].
The system was solved using Matlab’s ode45 solver and
a stepsize of 0.02, as in [1]. We used all computed points,
meaning that one timestep of the reservoir computer corresponds to a step of 0.02 in the Lorenz time scale. To
avoid unnecessary computations and save time we pregenerated a sequence of 105 samples that we used for
all numerical and experimental investigations. Following
[1], we used the x-coordinate trajectory for training and
testing, that we scaled by a factor of 0.01.

IV.
C.

(8a)

EXPERIMENTAL SETUP

Mackey-Glass chaotic series prediction

Our experimental setup, schematised in Fig. 2, consists of two main components: the opto-electronic reservoir and the FPGA board.

The Mackey-Glass delay differential equation
dx
x(t − τ )
=β
− γx
dt
1 + xn (t − τ )

(7)

with τ , γ, β, n > 0 was introduced to illustrate the appearance of complex dynamics in physiological control
systems [32]. To obtain chaotic dynamics, we set the parameters as in [1]: β = 0.2, γ = 0.1, τ = 17 and n = 10.
With these characteristics, the Kaplan-Yorke dimension
of the chaotic attractor is 2.1 [36].
The equation was solved using the Runge-Kutta 4
method with a stepsize of 1.0. To avoid unnecessary
computations and save time, both in simulations and experiments, we pre-generated a sequence of 106 samples
that we used for all numerical and experimental investigations.
During autonomous run, without the correct teacher
signal, the system slowly deviates from the desired trajectory. The MSE is used to evaluate both the training
phase and the autonomous run. We then compute the
number of correct prediction steps, i.e. steps for which
the MSE stays below the 10−3 threshold (see section

A.

Opto-electronic reservoir

The opto-electronic reservoir is based on previously
published works [8, 9]. The reservoir states are encoded into the intensity of an incoherent light signal, produced by a superluminiscent diode (Thorlabs SLD1550PA40). The Mach-Zehnder (MZ) intensity modulator
(EOSPACE AX-2X2-0MSS-12) implements the nonlinear function, its operating point is adjusted by applying a
bias voltage, produced by a Hameg HMP4040 power supply. A fraction (10%) of the signal is extracted from the
loop and sent to the readout photodiode (TTI TIA-525I)
and the resulting voltage signal is sent to the FPGA. An
optical attenuator (Agilent 81571A) is used to set the
feedback gain α of the system (see Eqs. (2a) and (2b)).
The resistive combiner sums the electrical feedback signal, produced by the feedback photodiode (TTI TIA525I), with the input signal from the FPGA to drive the
MZ modulator, with an additional amplification stage of

5
xi (n)

SLD

Clock

Pr

Att

MZ

ADC

Amp

Spool

90/10

Comb

Matlab

xi (n)

Out

Mask

y(n)

wi

Mi

Mi

DAC

Pf

OR

In

u(n)

u(n)
Mi × [u(n)

OR

y(n − 1)]

FMC151

ML605

Opto-electronic reservoir

Figure2.2.Schematic
Schematicrepresentation
representation of
of the
the experimental
experimental setup.
areare
Figure
setup. Optical
Optical and
andelectronic
electroniccomponents
componentsofofthe
thephotonic
photonicreservoir
reservoir
showniningrey
greyand
andblack,
black,respectively.
respectively. It
It contains
contains an
an incoherent
(MZ),
shown
incoherent light
light source
source(SLD),
(SLD),aaMach-Zehnder
Mach-Zehnderintensity
intensitymodulator
modulator
(MZ),
90 10beam
beamsplitter,
splitter, an
an optical
optical attenuator
attenuator (Att),
(Att), aa fibre
a a90/10
fibre spool
spool (Spool),
(Spool), two
twophotodiodes
photodiodes(P
(Prand
andPP),f ),a aresistive
resistivecombiner
combiner
(Comb)
and
an
amplifier
(Amp).
The
FPGA
board
implements
the
readout
layer
and
computes
the
output
signal
) in in
real
(Comb) and an amplifier (Amp). The FPGA board implements the readout layer and computes the output signal y(n)
real
). The computer, running Matlab,
time. It also generates the analogue input signal
time.
It also generates the analogue input signal I(n)) and
and acquires
acquires the
the reservoir
reservoirstates
states xi (n).
The computer, running Matlab,
controls the devices, performs the offline training and uploads all the data (
),
and
) on the FPGA.
controls the devices, performs the offline training and uploads all the data (u(n), wi and Mi ) on the FPGA.

(Spool) and the sampling frequency of the Analogue-to+27 dB (ZHL-32A+ coaxial amplifier) to span the entire
Digital converter (ADC). While increasing the latter inVπ interval of the modulator.
volves relatively high costs, one can lengthen the delay
As fairly
the neurons
arethis
time-multiplexed,
the spools
maximal
line
easily. In
work, we used two
of
reservoir
size fibre
depends
on the1.6
delay
the
fibre
spool
single mode
of lengths
km from
and 10
km,
approxi(Spool)
the sAnalogue-tomately.and
Thethe
firstsampling
producedfrequency
a delay ofof
7 93
and allowed
Digital
(ADC). into
While
latter
into fit converter
= 100 neurons
theincreasing
reservoir. the
The
second
volves
costs,theone
canuplengthen
delay
spool relatively
was used tohigh
increase
delay
to 49 the
s and
the
line
fairly size
easily.
work,
we used
of
reservoir
up toIn this
= 600.
In both
cases,two
the spools
reservoir
single
of lengths
1.6 km and
10MHz
km, (the
approxistatesmode
were fibre
sampled
at approximately
200
precise frequency,
in section
IV B,
depends
on allowed
the demately.
The firstgiven
produced
a delay
of 7.93
µs and
state into
was averaged
over 16The
samples
in
tolay
fitloop)
N =and
100each
neurons
the reservoir.
second
orderwas
to decrease
the noisethe
and
remove
the
transients
inspool
used to increase
delay
up to
49.2
µs and the
duced bysize
theup
finite
bandwidth
the Digital-to-Analogue
reservoir
to N
= 600. Inofboth
cases, the reservoir
converter
states
were (DAC).
sampled at approximately 200 MHz (the preciseThe
frequency,
givenisinoperated
section IV
depends
on the
experiment
as B,
follows.
First,
the deinlay
loop)
and
each
state
was
averaged
over
16
samples
in
put mask
and the teacher signal
), generated in
order
to decrease
the noise
andFPGA
removeboard,
the transients
inMatlab,
are uploaded
on the
which then
duced
by the
finite
bandwidth
of the Digital-to-Analogue
generates
the
masked
input signal
), sent to
converter
(DAC).
the reservoir
via the DAC. The resulting reservoir states
sampled is
and
averagedasbyfollows.
the FPGA,
and
The) are
experiment
operated
First,
thethen
insentmask
to the
real time.
That
is, the
design alput
Mcomputer
teacher
signal
u(n),
generated
in
i and thein
lows to are
capture
the reservoir
any desired
Matlab,
uploaded
on the states
FPGAfor
board,
which time
then
interval. After
traininginput
of the signal
reservoir,
optimal
readgenerates
the masked
Mithe
× u(n),
sent
to
out
weights
are
uploaded
on
the
board.
Because
of
the
the reservoir via the DAC. The resulting reservoir states
long delay
for by
thethe
offline
training,
the
xirelatively
(n) are sampled
andneeded
averaged
FPGA,
and then
reservoir
to beinreinitialised
order
sent
to theneeds
computer
real time. in
That
is, to
therestore
designthe
aldesired
dynamics
the internal
states
running
it
lows
to capture
theofreservoir
states
forprior
any to
desired
time
autonomously.
For
this
reason,
we
drive
the
system
with
interval. After training of the reservoir, the optimal readan weights
initialisation
ofon
128the
timesteps
(as illustrated
out
wi aresequence
uploaded
board. Because
of the
in
Fig.
9),
before
coupling
the
output
signal
with the the
inrelatively long delay needed for the offline training,
put and letting the reservoir computer run autonomously.
reservoir needs to be reinitialised in order to restore the
In this stage, the FPGA computes the output signal
desired dynamics of the internal states prior to running it
in real time, then creates a masked version
) and
autonomously. For this reason, we drive the system with
sends it to the reservoir via the DAC.
an initialisation sequence of 128 timesteps (as illustrated
As the
are processed
sequentially,
propin Fig.
9), neurons
before coupling
the output
signal due
withtothe
inagation
delay between
the intensity
modulator
(MZ) and
put
and letting
the reservoir
computer
run autonomously.
In this stage, the FPGA computes the output signal y(n)

the ADC, the output signal
) can only be computed
in real time, then creates a masked version M ×y(n) and
in time to update the 24-th neuron 23 + 1).i For this
sends it to the reservoir via the DAC.
reason, we set the first 23 elements of the input mask
As
neurons
are processed
sequentially,
to propto the
zero.
That way,
all neurons
contribute due
to solving
agation
delay
between
the
intensity
modulator
(MZ)
and
the task, but the first 23 do not “see” the input signal
the).ADC,
the
output
signal
y(n)
can
only
be
computed
Note that this reflects an aspect that is inherent to
in time
to updateimplementation
the 24-th neuron
x23 (n +computing
1). For this
any
experimental
of reservoir
reason,
we set
the first there
23 elements
of the input
with
output
feedback:
will necessarily
be amask
delayMi
to zero. the
That
way,
all neurons
contribute
to solving
the
between
time
when
the output
is available
and the
task,when
but the
23 isdosent
notback
“see”
thethe
input
signal
time
the first
output
into
system.
As I(n).
we
Noteseethat
thisthis
reflects
aspectonthat
is inherent
to any
will
below
has anan
impact
system
performance.
experimental implementation of time-multiplexed reservoir computing with feedback. In principle, the output
y(n) has to be computed
after the acquisition of the last
B. FPGA board
neuron xN −1 (n) at timestep n, but before the first neuron
x0 (n + 1) of the following timestep. However, in timeIn this work
weimplementations,
use a Xilinx ML605
evaluation
multiplexed
RC
these
states areboard,
consecpowered
by
a
Virtex
6
XC6VLX240T
FPGA
utive, and the experiment cannot be pausedchip.
to letThe
y(n)
board is paired with a 4DSP FMC151 daughter card, conbe computed. There may therefore be a delay (whose
taining one two-channel ADC and one two-channel DAC.
duration depends on the hardware used) before y(n) is
The ADC’s maximum sampling frequency is 250 MHz
computed and can be fed back into the reservoir. In the
with 14-bit resolution, while the DAC can sample at up
present experiment, this delay is approximately 115 ns,
to 800 MHz with 16-bit precision.
which corresponds to 23 neuron durations. As we will
The FPGA design is written in standard IEEE 1076see below this has an impact on system performance.
1993 VHDL language [38, 39] and compiled with Xilinx
ISE Design Suite 14.7, provided with the board. We also
used Xilinx ChipScope Pro Analyser to monitor signals
B.for FPGA
board
on the board, mostly
debugging
and testing.
The simplified schematics of the design is depicted in
In 3.
this
work we use
a Xilinx
ML605
evaluation
board,
Fig.
Rectangular
boxes
represent
modules
(entities),
and
rounded
for electronic
components
powered
by arectangles
Virtex 6 stand
XC6VLX240T
FPGA
chip. The
on
the is
ML605
the FMC151
daughtercard
board
pairedboard,
with anamely
4DSP FMC151
daughter
card, conand
the one
onboard
Marvell ADC
Alaska
PHY
(88E1111)
taining
two-channel
and
onedevice
two-channel
DAC.
for
(ETH).
TheEthernet
ADC’s communications
maximum sampling
frequency is 250 MHz
The14-bit
operation
of thewhile
FPGA
controlled
from at
theup
with
resolution,
theisDAC
can sample
computer,
running
Matlab,
via a simple custom prototo 800 MHz
with 16-bit
precision.
colThe
through
a Gbit
Ethernet
connection.
DataIEEE
and varFPGA
design
is written
in standard
1076ious
such as
or state
1993commands,
VHDL language
[38,memory
39] andread/write,
compiled with
Xilinx
ISE Design Suite 14.7, provided with the board. We also

6
used Xilinx ChipScope Pro Analyser to monitor signals
on the board, mostly for debugging and testing.
The simplified schematics of the design is depicted in
Fig. 3. Rectangular boxes represent modules (entities),
and rounded rectangles stand for electronic components
on the ML605 board, namely the FMC151 daughtercard
and the onboard Marvell Alaska PHY device (88E1111)
for Ethernet communications (ETH).
The operation of the FPGA is controlled from the
computer, running Matlab, via a simple custom protocol through a Gbit Ethernet connection. Data and various commands, such as memory read/write, or state
change, are encapsulated into standard UDP packages.
The etherphy module interfaces the FPGA design with
the Marvell Ethernet PHY chip, and the etherctrl module receives the UDP packets (frames) and decodes the
commands and the data. It also creates the frames for
sending data from FPGA to the computer.
Blocks of Random-Access Memory (BRAM) are used
to store data, such as teacher inputs u(n), input masks
Mi and readout weights wi , that are generated on the
computer and uploaded on the FPGA. Each type of data
is assigned a specific module, since they vary in size (e.g.
600 values for the input mask and up to 3000 for the
teacher signal) and resolution, as will be explained below. The bramRec is a buffer-like module, designed to
transfer the signal recorded from the experiment directly
to the computer through Ethernet, without permanently
storing it in memory. It consists of two blocks of RAM
of 2048 bytes each, that are used as follows: while the
recorded signal is written into the first block, bram2ether
reads the contents of the second, that is then encapsulated into four UDP frames sent to the computer. When
the first block is full, the blocks are switched and the
process continues.
The FMC151 daughtercard is interfaced with the rest
of the design through the fmc151 module, that outputs
two 14-bit signals from the ADCs and receives two 16-bit
signals from the DACs. The FMC151 card is also used
to deliver a clock signal from an external clock generator, that produces a high-precision signal, allowing to
synchronise the FPGA with the delay loop of the experimental setup. This clock signal was generated by
the Hewlett Packard 8648A signal generator. In our experiments with two delay loops (see section IV A), we
fine-tuned the clock frequency so as to fit 16 samples per
neuron into the roundtrip time T . Specifically, we sampled the N = 600 reservoir states at 195.4472 MHz with
a large fibre spool, and at 203.7831 MHz with a small
spool and N = 100.
The fpga2exp module controls the signal sent to the
opto-electronic reservoir through the DAC. At the training phase, it generates the masked input signal Mi ×u(n)
by multiplying the inputs u(n) by the mask Mi , both
being read from the BRAMs. During the autonomous
run, it receives the reservoir output signal y(n), computed by the comprcout module, masks it and transfers
to the DAC.

The neuron states xi (n) from the photonic reservoir are
sampled and averaged by the exp2fpga module. At the
training phase, these are buffered in bramRec, then processed by the bram2ether module and sent to the computer. During the autonomous run, the reservoir states
are used by the comprcout, together with the readout
weights wi , read from the bramWts memory, to compute
the reservoir output y(n). It is then injected back into
the reservoir through the fpga2exp, and also transferred
to the computer through the bramRec and bram2ether
modules.
The design is driven by two clocks: the experimental
clock expclk (around 200 MHz, depending on the loop
delay T ) that operates data acquisition and generation
modules and allows to synchronise the FPGA with the
experiment, and the 125 MHz Ethernet clock ethclk.
Both clocks have to be managed properly within the design, as several signals, such as inputs or weights, coexist in both clock domains. That is, data to BRAMs
comes from the Ethernet modules, and is thus driven by
the ethclk clock. On the other hand, this same data
is used by the fpga2exp module, and has to appear in
the expclk clock domain. To this end, we exploit the
dual-port capability of Xilinx block RAMs. That is, data
is written into memory blocks through port A at clock
ethclk and read from port B at clock expclk (and vice
versa for the bramRec). This allows for smooth transition of data between clock domains. The two clock domains are depicted in Fig. 3 as follows: signals running
at expclk are shown in solid lines, and those clocking at
ethclk are drawn with dashed lines.
The arithmetic operations computed by the FPGA are
performed on real numbers. However, the chip is a logic
device, designed to operate bits. The performance of the
design thus highly depends on the bit-representation of
real numbers, i.e. the precision. The main constraint
comes from the ADC and DAC, limited to 14 and 16
bits, respectively. Numerical simulations, reported in
[18], show that such precision is sufficient for all tasks
studied in this work. It was also shown in [18] that the
precision of the readout weights wi has a significant impact on the performance of the system. For this reason we designed the experiment for optimal utilisation of
the resolution available. The reservoir states were tuned
to lie within a ] − 1, +1[ interval. They are thus represented as 16-bit integers, with 1 bit for the sign and
15 bits for the decimal part. Another limitation comes
from DSP48E slices, used to multiply the states xi (n)
by the readout weights wi . These blocks are designed to
multiply a 25-bit integer by a 18-bit integer. To meet
these requirements, we also keep the readout weights wi
within the ] − 1, 1[ interval and represent them as 25-bit
integers, with 1 sign bit and 24 decimal bits. To ensure
that wi ∈]− 1, 1[, we amplify the reservoir states digitally
inside the FPGA. That is, the xi (n) are multiplied by 8
after acquisition, prior to computing the output signal
y(n).

7

ETH

Matlab

etherphy

etherctrl

bram2ether

ethclk
Mi × I(n)

fpga2exp

DAC

Experiment

bramInp
y(n)

fmc151
xi (n)

BRAMs

u(n)

xi (n)

exp2fpga

ADC

Mi

comprcout

bramMsk
wi

CLK

Clock Gen

expclk

FMC151

bramRec

bramWts

ML605

Figure 3. Simplified schematics of the FPGA design. Modules (entities) are represented by rectangular boxes, onboard electronic
components are shown with rounded rectangles. External hardware, such as the computer, running Matlab, the opto-electronic
reservoir and the external clock generator are shown in grey. The design is driven by two clocks: the experimental clock expclk
and the Ethernet clock ethclk. Signals from these two clock domains are drawn in solid and dashed lines, respectively.

C.

Numerical simulations

In addition to the physical experiments, we investigated the proposed setup in numerical simulations, to
have a point of comparison and identify possible malfunctions. To this end, we developed three models that
simulate the experiment to different degrees of accuracy.
Our custom Matlab scripts are based on [8, 18].
Idealised model: It incorporates the core characteristics of our reservoir computer, i.e. the ring-like
architecture, the sine nonlinearity and the linear
readout layer (as described by equations 2 and 3),
disregarding all experimental considerations. We
use this model to define the maximal performance
achievable in each configuration.
Noiseless experimental model: This model emulates
the most influential features of the experimental
setup, such as the high-pass filter of the amplifier,
the finite resolution of the ADC and DAC, and precise input and feedback gains. This model allows to
cross-check the experimental results and to easily
identify the problematic points.
Noisy experimental model: Contrary to the noiseless
numerical models introduced above, our experimental implementation is noisy, which, as will be
explained below, has a significant impact on performance. In order to compare our experimental
results to a more realistic model, we estimated the
level of noise present in the experimental system
(see section V A), and incorporated this noise into
the noisy version of the experimental model.
V.

RESULTS

In this section we present the experimental results,
compare them to numerical simulations and discuss the

performance of the reservoir computer on each task introduced in section III.
The two periodic signal generation tasks were solved
using a small reservoir with N = 100 and a fibre spool
of approximately 1.6 km. The chaotic signal generation
tasks, being more complex, required a large reservoir of
N = 600 for decent results, that we fit in a delay line of
roughly 10 km.

A.

Noisy reservoir

For most tasks studied here, we found the experimental
noise to be the major source of performance degradation
in comparison to numerical investigations. In fact, previously reported simulations [18] considered an ideal noiseless reservoir, while our experiment is noisy. This noise
is generated by the active components of the setup: the
amplifier, which has a relatively high gain and is therefore very sensitive to small parasitic signals on the input,
the DAC and the photodiodes. In-depth experimental
investigations (not presented in this paper) have shown
that, in fact, each component contributes more or less
equally to the overall noise level. Thus, it can not be reduced by replacing one “faulty” component. Neither can
it be averaged out, as the output value has to be computed at each timestep. This noise was found to have
a significant impact on the results, as will be shown in
the following sections. For this reason we estimated the
level of noise present in the experimental system and incorporated it to the numerical models. This allows us to
“switch off” the noise in simulations, which is impossible
experimentally.
Fig. 4 shows numerical and experimental reservoir
states of a 100-neuron reservoir, as received by the readout photodiode. That is, the curves depict the timemultiplexed neurons: each point represents a reservoir
state x0...99 (n) at times n = 1 and n = 2. The system
does not receive any input signal I(n) = 0. The exper-

8

Noise amplitude (arb. units)

0.01

B. Frequency
Frequency generation
generation
B.
Experiment
Simulation with noise

0.005

0

−0.005
x0 (1)

x50 (1)

x0 (2)

x50 (2)

x99 (2)

Time-multiplexed reservoir states xi (n)

Figure 4.
4. Illustration
the experimental
experimental
Figure
Illustration of
of the
the noise inside the
reservoir. Experimental
Experimental (solid
(solid grey
grey line)
line) and
and numerical
numerical (dot(dotreservoir.
ted black
black curve)
curve) reservoir
reservoir states
states xi (n)) are
are shown
shown in
in the
the case
case
ted
when the
the input
input signal
signal is
is null
null I(n)) =
= 0,
0, scaled
scaled so
so that
thatin
innormal
normal
when
experimental conditions
conditions (non-zero
(non-zero input)
input) they
they would
would lie
lie in
in aa
experimental
1] interval.
interval. Although
Although the
the input
input signal
signal isis null
null I(n)) =
= 0,
0,
[−1, 1]
actual neurons
neurons are
are non-zero
non-zero because
because of
of noise.
noise. NumeriNumerithe actual
cal noise
noise was
was generated
generated with
with aa Gaussian
Gaussian random
random distribution
distribution
cal
with standard
standard deviation
deviation of
of 11 × 10
10−3 so
so that
that to
to reproduce
reproduce the
the
with
noise level
level of
of the
the experiment.
experiment.
noise

signal is
is plotted
plotted with
with aa solid
solid grey
grey line.
line. We
We use
use
imental signal
compute the
the experimental
experimental noise
noise level
level by
by taking
taking the
the
it to compute
deviation of
of the
the signal,
signal, which
which gives
gives 22 × 10
10−3 .
standard deviation
then employed
employed this
this noise
noise level
level in
in the
the noisy
noisy experexperWe then
model to
to compare
compare experimental
experimental results
results to
to nunuimental model
simulations. The
The dotted
dotted black
black curve
curve in
in Fig.
Fig. 44
merical simulations.
the response
response of
of the
the noisy
noisy experimental
experimental model,
model, with
with
shows the
same amount
amount of
of Gaussian
Gaussian noise
noise (standard
(standard deviation
deviation
the same
of 2.0
The choice
choice of
of aa
as in
in the
the experiment.
experiment. The
2 × 10
10−3 )) as
Gaussian noise
noise distribution
distribution was
was validated
validated by
by experimenexperimental measurements.
measurements.
The level of
of the
the experimental
experimental noise
noise can
can also
also be
be charcharacterised by
by the
the Signal-to-Noise
Signal-to-Noise Ratio
Ratio (SNR),
(SNR), defined
defined as
as
[40]
!
2
RMS
RMSsignal
signal
SNR
,
SNR =
= 10
10 log
log10
2
10
RMS
RMSnoise
noise
where RMS
RMS is
is the
the Root
Root Mean
Mean Square
Square value,
value, given
given by
by
v
u
N
u1 X
t
x2i .
RMS(x
)
=
RMS( i ) = N
i=1
=1

We
= 0.2468 and RMS noise =
We measured
measured RMS
RMSsignal
signal = 0 2468 and RMSnoise
0.0023,
0023, so
so the
the SNR
SNR is
is equal
equal to
to approximately
approximately 40
40 dB
dB in
in
this
case.
Note
that
this
figure
is
given
as
an
indicator
this case. Note that this figure is given as an indicator
of
as the
the RMS
RMS of
of the
the reservoir
reservoir
of order
order of
of magnitude
magnitude only
only as
states
depends
on
the
gain
parameters
(α
and
states depends on the gain parameters ( and β,, see
see secsection
II)
and
varies
from
one
experiment
to
another.
tion II) and varies from one experiment to another.

We found
found the
the frequency
frequencygeneration
generationtask
taskto
tobe
bethe
theonly
only
We
one not
not affected
affected by
by noise.
noise. That
That is,
is, our
our experimental
experimentalrereone
sults matched
matched accurately
accurately the
the numerical
numerical predictions
predictions reresults
ported in
in [24].
[24]. Concretely,
Concretely, we
we expected
expected aa bandwidth
bandwidth ofof
ported
ν ∈ [0.06,
[0 06,π]
π] with
with aa 100-neuron
100-neuron reservoir.
reservoir. The
The upper
upper
limit
limit isis aa signal
signaloscillating
oscillatingbetween
between−11and
and11and
andisisgiven
given
by
by half
half of
of the
the sampling
sampling rate
rate of
of the
the system
system (the
(the Nyquist
Nyquist
frequency).
frequency). The
The lower
lower limit
limit isis caused
caused by
by the
the memory
memory
limitation
limitation of
of the
the reservoir.
reservoir. In
In fact,
fact, low-frequency
low-frequency ososcillations
cillations correspond
correspond to
to longer
longer periods,
periods, and
and the
the neural
neural
network
network can
can no
no longer
longer “remember”
“remember” aa sufficiently
sufficiently long
long
segment
segment of
of the
the sine
sine wave
wave so
so as
as to
to keep
keep generating
generating aa sisinusoidal
nusoidal output.
output. These
These numerical
numerical results
results are
are confirmed
confirmed
experimentally
experimentally here.
here.
We
rangingfrom
from0.01
0 01
We tested
testedour
oursetup
setupon
onfrequencies
frequenciesν ranging
to
to π,, and
andfound
foundthat
thatfrequencies
frequencieswithin
within[0.1,
[0 ,π]
π]are
aregenergenerate
ate accurately
accurately with
with any
any random
random input
input mask.
mask. Lower
Lowerfrefrequencies
quencies between
between 0.01
0 01 and
and 0.1,
0 1, however,
however, were
were produced
produced
properly
properly with
with some
some random
random masks,
masks, but
but not
notall.
all. For
Forthis
this
reason,
reason, we
we investigated
investigated the
the [0.01,
[0 01 0.1]
1] interval
interval more
more preprecisely,
cisely,since
sincethis
thisisiswhere
wherethe
thelower
lowerlimit
limitofofthe
thebandwidth
bandwidth
lies.
lies. For
For each
each frequency,
frequency, we
we ran
ran the
theexperiment
experiment10
10times
times
for
for 10k
10k timesteps
timestepswith
withdifferent
differentrandom
randominput
inputmasks
masksand
and
counted
counted the
the number
number of
of times
times the
the reservoir
reservoir produced
produced aa
sine
see
sine wave
wave with
with the
the desired
desired frequency
frequency (MSE
(MSE< 10
10−3 , , see
section
section III
IIIB)
B) and
and amplitude
amplitude of
of 1.1. The
Theresults
resultsare
areshown
shown
in
in Fig.
Fig. 5.
5. Frequencies
Frequencies below
below 0.05
0.05 are
are not
not generated
generatedcorcorrectly
rectly with
with most
most input
input masks.
masks. At
At ν == 0.7
0 7 the
the output
output
isis correct
correct most
most of
of the
the times,
times, and
and for
for ν ==0.08
0 08and
andabove
above
the
theoutput
outputsine
sinewave
waveisiscorrect
correctwith
withany
anyinput
inputmask.
mask. The
The
bandwidth
bandwidth of
of this
this experimental
experimental RC
RC isis thus
thus ν ∈ [0.08,
[0 08,π].
π].
Given
Given the
the roundtrip
roundtrip time
time T == 7.93
7 93 µs,s, this
this results
results inin
aa physical
physical bandwidth
bandwidth of
of 1.5
1 5 –– 63
63 kHz.
kHz. Note
Note that
that frefrequencies
quencies within
within this
this interval
interval can
can be
be generated
generated with
with any
any
random
Lowerfrequencies,
frequencies,down
downto
to0.02,
0
randominput
inputmask
maskMi .. Lower
could
could also
also be
be generated,
generated, but
but only
only with
with aa suitable
suitable input
input
mask.
mask.
Fig.
Fig. 66 shows
shows an
an example
example of
of the
the output
output signal
signal during
during
the
the autonomous
autonomous run.
run. The
The system
system was
was trained
trained for
for 1000
1000
timesteps
timestepsto
togenerate
generateaafrequency
frequencyof
ofν ==0.1,
0 1,and
andsuccesssuccess−9
fully
fully accomplished
accomplished this
this task
task with
with aa MSE
MSEofof5.6
5 × 10
10 .
The
The above
above results
results were
were obtained
obtained by
by scanning
scanning the
the ininput
put gain
gain β and
and the
the feedback
feedback gain
gain α to
to obtain
obtain the
the best
best
results.
results. It
It was
was found
found that
that β has
haslittle
littleimpact
impacton
onthe
thesyssystem
tem performance
performance so
so long
long as
as itit isis chosen
chosen in
in the
the interval
interval
β ∈ [0.02,
[0 02 0.5],
5], while
while the
the feedback
feedback gain
gain α, , on
on the
the concontrary,
trary, has
has to
to lie
lie within
within aa narrow
narrow interval,
interval, approximately
approximately
between
between 0.9
0 9 and
and 1.0,
1 0, otherwise
otherwise the
the reservoir
reservoir yields
yields very
very
poor
results.
poor results.
C.
C. Random
Random pattern
pattern generation
generation

The
The random
random pattern
pattern generation
generation task
task isis more
more complex
complex
than
frequency
generation
and
is
slightly
affected
than frequency generation and is slightly affectedby
bythe
the
experimental
noise.
The
goal
of
this
task
is
two-fold:
experimental noise. The goal of this task is two-fold: “re“re-

9
1616
0.16

1
1515
0.15
Reservoir
Reservoirstate
statex0 (n)
Reservoir state

Success
Success
rate
Successrate
rate

0.8
0.6
0.4

00
0

02
000.02
02

04
06
000.04
04
000.06
06
Relative
frequency
Relative
Relativefrequency
frequency ν

08
000.08
08

000.1

Figure 5.
5. Determination
Determination
of
the
lower
limit
of
the
reservoir
Figure
5.
Determination
of
the
lower
limit
of
the
reservoir
5.
Determination of
ofthe
thelower
lowerlimit
limitof
ofthe
thereservoir
reservoir
Figure
computer bandwidth.
bandwidth.
Frequencies
above
08
are
generated
computer
bandwidth.
Frequencies
above
00008
are
generated
computer
bandwidth. Frequencies
Frequenciesabove
above0.08
08are
aregenerated
generated
computer
very well
well with
with
any
of
the
10
random
input
mask,
and
are
very
well
with
any
of
the
10
random
input
mask,
and
are
very
very
well
with any
any of
of the
the 10
10random
randominput
inputmask,
mask,and
andare
are
therefore
not
shown
on
the
plot.
Frequencies
below
0
05
fail
therefore
not
shown
on
the
plot.
Frequencies
below
0
05
fail
therefore
therefore not
not shown
shownon
onthe
theplot.
plot. Frequencies
Frequenciesbelow
below0.05
0 05fail
fail
with most
most input
input
masks.
We
thus
consider
08
as
the
lower
most
masks.
We
thus
consider
00008
as
the
lower
with
with
most
input masks.
masks. We
Wethus
thusconsider
consider0.08
08as
asthe
thelower
lower
limit of
of the
the bandwidth,
bandwidth, but
but keep
keep in
inmind
mindthat
thatfrequencies
frequenciesas
as
limit
bandwidth,
but
keep
in
mind
that
frequencies
limit of the bandwidth, but keep in mind that frequenciesasas
low as
as 0.02
0 02
02 could
could also
also be
be generated,
generated,but
butonly
onlywith
withaaacarefully
carefully
0
could
also
be
generated,
but
only
with
carefully
low
low as 0 02 could also be generated, but only with a carefully
picked input mask.
mask.
picked
picked input
input mask.
mask.
0.8

Output
signals
(arb.
units)
Output
signals
(arb.
Output
signals
(arb.units)
units)

1313
0.13
1212
0.12
1111
0.11

0.2
0

1414
0.14

0.6

RC output
output
RC
RC output
output
Target
Target
output
Target output

0.4
0.2
0
−0.2
−0.4
1000
1000
1000

1100
1200
1100
1200
1100
1200
Discrete time
time n
Discrete
Discrete time

1300
1300
1300

Figure 6.
6. Example
Example of
of an
an autonomous
autonomous run
run output
outputsignal
signalfor
for
Figure
Example
of
an
autonomous
run
output
signal
for
Figure
6. generation
Example of
anwith
autonomous
run
output
signal
for
frequency
generation
task
with
=0.1.
1. The
Theexperiment
experiment
congeneration
task
with
=
00 1.
The
experiment
confrequency
task
ν=
confrequency
generation
task
with
=
0
1.
The
experiment
continues beyond
beyond the
the range
range of
of the
the figure.
figure.
the
range
of
the
figure.
tinues
tinues beyond the range of the figure.

member”
member” aa pattern
pattern of
of aaagiven
givenlength
lengthL and
andbe
beable
ableto
torerepattern
of
given
length
and
be
able
to
remember”
pattern
of a given
length We
beshown
able to
reproduce
an
duration.
have
nuproduce ititafor
for
an unlimited
unlimited
duration.
Weand
have
shown
nuunlimited
duration.
We
have
shown
nuproduce
for an
unlimited
duration.
We
have
shown
numerically
aa noiseless
51-neuron
reservoir
is
capable
mericallyitthat
that
noiseless
51-neuron
reservoir
is
capable
noiseless 51-neuron reservoir is capable
merically
thatpatterns
a noiseless
51-neuron
reservoir
is capable
of
up
long
This
of generating
generating
patterns
up to
to 51-element
51-element
long[18].
[18].
This
patterns
up
to
51-element
long
[18].
This
is
logical
as,
each
the
of
patterns
up to 51-element
longof
This
is aagenerating
logical result,
result,
as, intuitively,
intuitively,
eachneuron
neuron
of[18].
thesyssysresult,
as,
intuitively,
each
neuron
of
the
system
is
to
“memorise”
one
value
of
is
a logical
result,
intuitively,
each
neuron
the system
is expected
expected
to as,
“memorise”
one
value
ofthe
theofpattern.
pattern.
to
“memorise”
one
value
of
the
pattern.
Simulations
of
noisy
100-neuron
similar
to
tem
is expected
“memorise”
onereservoir,
value of the
pattern.
Simulations
of aato
noisy
100-neuron
reservoir,
similar
to
noisy
100-neuron
reservoir,
similar
to
the
that
maximum
pattern
Simulations
of asetup,
noisy show
100-neuron
to
the experimental
experimental
setup,
show
thatthe
thereservoir,
maximumsimilar
pattern
setup,
show
that
the
maximum
pattern
length
isis reduced
L=
noise
length
reduced down
down
toshow
=13.
13. This
This
meansthat
thatpattern
noise
the
experimental
setup,to
that
the means
maximum
down
to
=
13.
This
means
that
noise
significantly
reduces
the
effective
of
system.
significantly
reduces
theto
effective
memory
ofthe
thethat
system.
length
is reduced
down
= 13.memory
This means
noise
the
effective
memory
of
the
system.
In
fact,
the
noisy
neural
network
has
to
take
into
account
In fact, the noisy
neural
network
has
to take
takeofinto
into
account
neural
to
significantly
reduces
thenetwork
effectivehas
memory
theaccount
system.
the
slight
of
output
from
the
patthefact,
slight
deviations
of the
the
output
from
the target
target
patdeviations
of
the
output
from
the
target
patIn
thedeviations
noisy neural
network
has
to take
into
account
tern
so
as
to
be
able
to
follow
the
pattern
disregarding
tern
so
as
to
be
able
to
follow
the
pattern
disregarding
be
able
to
follow
the
pattern
disregarding
the slight deviations of the output from the target pattern so as to be able to follow the pattern disregarding

0.1
1250
1250
1250

1300
1350
1300
1350
1300
1350
Discrete
time
Discrete
time
n
Discrete
time

1400
1400
1400

Figure
7.7.Example
Example
ofof
behaviour
ofone
one
neuron
inain
Figure
7.
Example
of
behaviour
of
one
neuron
aanoisy
noisy
Figure
7.
of
behaviour
of
neuron
inin
Figure
Example
behaviour
of
one
neuron
anoisy
noisy
experimental
reservoir.
For
clarity,
the
range
the
experimental
reservoir.
For
clarity,
the
range
the
axis
experimental
reservoir.
For
clarity,
the
range
ofofofthe
YYYaxis
experimental
reservoir.
For
clarity,
the
range
of
the
Yaxis
axis
limited
toto
the
area
ofinterest.
interest.
Because
ofnoise,
noise,
despite
limited
to
the
area
of
interest.
Because
noise,
despite
isisisis
limited
to
the
area
of
Because
ofof
despite
aaa a
limited
the
area
of
interest.
Because
of
noise,
despite
periodic
input
signal
the
reservoir
state
takes
similar,
periodic
input
signal
the
reservoir
state
takes
similar,
periodic
input
signal
u(n),),),the
reservoir
state
takes
similar,
periodic
input
signal
),
the
reservoir
state
takes
similar,
but
not
identical
values.
but
not
identical
values.
but
not
identical
values.
but
not
identical
values.

theseimperfections.
imperfections. Fig.
Fig. 777illustrates
illustratesthis
thisissue.
issue. PeriPerithese
this
issue.
these
theseimperfections.
imperfections.Fig.
Fig. 7illustrates
illustrates
this
issue.PeriPeriodicoscillations
oscillationsof
ofone
oneneuron
neuronof
ofthe
thereservoir
reservoirare
areshown,
shown,
odic
oscillations
of
one
neuron
of
the
reservoir
are
shown,
odic
odic oscillations of one neuron of the reservoir are shown,
withintended
intendedfocus
focuson
onthe
theupper
uppervalues
valuesand
andan
anadequate
adequate
with
on
the
upper
values
and
an
adequate
with
withintended
intendedfocus
focus
on
the
upper
values
and
an
adequate
magnification
so
as
to
see
the
small
variations.
shows
magnification
soso
asas
toto
seesee
thethe
small
variations.
shows
magnification
so
as
to
see
the
small
variations.
ItItItshows
magnification
small
variations.
It shows
that
the
neuron
oscillates
between
similar,
but
not
identhat
the
neuron
oscillates
between
similar,
butbut
notnot
identhat
the
neuron
oscillates
between
similar,
but
not
identhat
the
neuron
oscillates
between
similar,
identical
values.
This
makes
the
generation
task
much
more
tical
values.
This
makes
the
generation
task
much
more
tical
values.
makes
the
generation
task
much
more
tical
values.This
This
makes
the
generation
task
much
more
complex,and
andrequires
requiresmore
morememory,
memory,hence
hencethe
themaximal
maximal
complex,
more
memory,
hence
the
maximal
complex,
complex,and
andrequires
requires
more
memory,
hence
the
maximal
patternlength
lengthisisisshorter.
shorter.
pattern
length
shorter.
pattern
pattern length is shorter.
Weobtained
obtainedsimilar
similarresults
resultsin
inthe
theexperiments.
experiments.Fig.
Fig.888
We
obtained
similar
results
in
the
experiments.
Fig.
We
Wethe
obtained
similar
results
in the experiments.
Fig.
shows
evolution
ofthe
the
MSEmeasured
measured
duringthe
thefirst
first 8
shows
the
of
MSE
during
shows
theevolution
evolution
of
the
MSE
measured
during
the
first
shows
the
evolution
of
the
MSE
measured
during
the
first
1ktimesteps
timestepsof
of10k-timestep
10k-timestepautonomous
autonomousruns
runswith
withdifdif1k
1k
timesteps
of
10k-timestep
autonomous
runs
with
dif1k
timesteps
of
10k-timestep
autonomous
runs
with
differentpattern
patternlengths.
lengths.Plotted
Plottedcurves
curvesare
areaveraged
averagedover
over
ferent
ferent
pattern
lengths.
Plotted
curves
are
averaged
over
ferent
pattern
lengths.
Plotted
curves
are
averaged
over
100
runs
of
the
experiment,
with
5
random
input
masks
100
100 runs
runs of
of the
the experiment,
experiment, with
with 55 random
randominput
inputmasks
masks
10020
runs
of the
experiment,
with
5 random
inputminmasks
and
20
random
patterns
foreach
each
length.
Theinitial
initial
minand
random
patterns
for
length.
The
and
20
random
patterns
for
each
length.
The
initial
minand
20
random
patterns
for
each
length.
The
initial
minimum(at
(atn =
=128)
128)corresponds
correspondsto
tothe
theinitialisation
initialisationof
imum
imum
(at
=
128)
corresponds
to
the
initialisation
ofof
imum
(at
=
128)
corresponds
to
the
initialisation
thereservoir
reservoir(see
(seesection
sectionIV
IVA),
A),then
thenthe
theoutput
outputisisiscoucou- of
the
the
reservoir
(see
section
IV
A),
then
the
output
coutheback
reservoir
(see
sectionruns
IV A),
then the output
is coupled
backand
andthe
thesystem
system
runs
autonomously.
Patterns
pled
autonomously.
Patterns
pled
back
and
the
system
runs
autonomously.
Patterns
pled
back
and
the
system
runs
autonomously.
Patterns
withL =
=12
12or
orless
lessare
aregenerated
generatedvery
verywell
welland
andthe
theererwith
with
=
12
or
less
are
generated
very
well
and
the
erwith
=
12Patterns
or
less are
generated
very
well
and the
ror
stayslow.
low.
Patterns
oflength
length13
13show
showan
anincrease
increase
inerror
stays
of
inin
ror
stays
low.
Patterns
of
length
13
show
an
increase
ror
stays
low.
Patterns
of
length
13
show
an
increase
MSE,but
butthey
theyare
arestill
stillgenerated
generatedreasonably
reasonablywell.
well.For
Forin
MSE,
MSE,
but
they
are
still
generated
reasonably
well.
For
MSE,
but
they
are
still
generated
reasonably
well.
For
longerpatterns,
patterns,the
thesystem
systemdeviates
deviatesto
toaaadifferent
differentperiperilonger
longer
patterns,
the
system
deviates
to
different
peri−3
longer
patterns,
the
system
deviates
to
a
different
periodicbehaviour,
behaviour,and
andthe
theerror
errorgrows
growsabove
above10
10 .
odic
odic
behaviour,
and
the
error
grows
above
10
odic
and
the error
grows
above
10 during
Fig.
999shows
example
of
output
signal
Fig.behaviour,
showsan
an
example
ofthe
the
output
signal
during
Fig.
shows
an
example
of
the
output
signal
during
Fig.
9
shows
an
example
of
the
output
signal
during
the
autonomous
run.
The
system
was
trained
for
the
autonomous
run.
The
system
was
trained
for1000
1000
the autonomous run. The system was trained
for
1000
the
autonomous
run.
The
system
was
trained
for
1000
timesteps
to
generate
a
pattern
of
length
10.
The
resertimesteps
to
generate
a
pattern
of
length
10.
The
resertimesteps to generate a pattern of length 10. The resertimesteps
to generate
a pattern
ofthe
length
10. pattern
The
reservoir
computer
successfully
learned
desired
voir
computer
successfully
learned
the
desired
pattern
voir
computer
successfully
learned
the
desired
pattern
voir
computer
successfully
learned
the
desired
and
the
output
matches
the
signal.
Fig.
and
the
outputperfectly
perfectly
matches
thetarget
target
signal.pattern
Fig.
and
the
output
perfectly
matches
the
target
signal.
Fig.
and
the output
perfectly
matches
the target
signal.
Fig.
10
illustrates
aaacase
with
(L
10
illustrates
case
withaaalonger
longerpatter
patter
14),that
that
10
illustrates
case
with
longer
patter
(( ===14),
14),
that
could
not
by
the
can
from
10 illustrates
a case
with
asystem.
longerAs
patter
(beseen
=
14),
that
could
notbe
belearned
learned
by
thesystem.
As
canbe
seen
from
could
not
be
learned
by
the
system.
As
can
be
seen
from
the
plot,
the
RC
captured
the
general
patcould
not
be
learned
by the
As
canofof
be
seen
from
the
plot,
the
RC
captured
thesystem.
generalshape
shape
ofthe
the
patthe
plot,
the
RC
captured
the
general
shape
the
pattern,
individual
the but
plot,
thenot
RCaccurately
captured generate
the
general
shape ofpoints.
the pattern,
butcan
can
not
accurately
generate
individual
points.
tern,
but
can
not
accurately
generate
individual
points.
−3
The
MSE
this
run
isisis5.2
,,,which
isisisabove
the
The
MSE
of
this
run
10
which
abovepoints.
the
tern,
butof
can
not
accurately
generate
individual
The
MSE
of
this
run
55 × 10
10
which
above
the
−3
acceptable
10
threshold.
threshold.
acceptable
10
The MSE10
of this
run is 5
10 , which is above the
acceptable
threshold.
We
acceptable
10 the
threshold.
Wealso
alsotested
tested
thestability
stabilityof
ofthe
thegenerator
generatorby
byrunrunWe
also
tested
the
stability
of
the
generator
by
run9
ning
itititfor
several
hours
(∼
10
timesteps)
with
random
ningWe
for
several
hours
(
10
timesteps)
with
random
ning
for
several
hours
(
10
timesteps)
with
random
also tested the stability of the generator by runpatterns
11
was
patterns
oflengths
lengths
10,
11and
and
12.The
Theoutput
outputsignal
signal
was
patterns
of
lengths
10,
11
and
12.
The
output
signal
was
ning it of
for
several10,
hours
( 12.
10
timesteps)
with
random
patterns of lengths 10, 11 and 12. The output signal was

10
10
10
10

Autonomous
run
MSE
Autonomous run
run MSE
MSE
Autonomous

10
10
10−1

Autonomous
run
starts
Autonomous
Autonomousrun
runstarts
starts

10
10
10−2

L

10
10
10−3

10
10
10−4
000

200
200
200

400
400
400

600
600
600

800
800
800

L

· 16
· ······16
· ········
· 16

1616
16

· 15
· ······15
· ········
· 15
· 14
· ······14
· ········
· 14

1515
15
1414
14

· 13
· ······13
· ········
· 13

1313
13

· 12
· ······12
· ········
· 12
· 11
· ······11
· ········
· 11
· 10
· ······10
· ········
· 10

1212
12
1111
11
1010
10

1000
9000
1000
1000· ······ 9000
9000
Discrete
time
Discrete
time
Discrete time n

9200
9200
9200

9400
9400
9400

9600
9600
9600

9800
9800
9800

10000
10000
10000

Figure 8.
8.
Evolution
of
MSE
during
experimental
autonomous
generation
periodic
random
patterns
lengths
10,
, ,16
Figure
MSE
during
experimental
autonomous
generation
ofof
periodic
random
patterns
ofof
lengths
===
10,
Figure
8. Evolution
Evolutionof
ofMSE
MSEduring
duringexperimental
experimentalautonomous
autonomousgeneration
generationof
ofperiodic
periodic
random
patterns
of
lengths
10,
.16
16
Figure
8.
Evolution
of
random
patterns
of
lengths
L=
10,
.......,.,.16
as
function
of
discrete
time
.
The
autonomous
run
starts
at
=
128,
as
indicated
by
the
arrow.
Patterns
shorter
than
13
are
as
function
of
discrete
time
.
The
autonomous
run
starts
at
=
128,
as
indicated
by
the
arrow.
Patterns
shorter
than
13
are
as function
function of
of discrete
discretetime
timen. . The
Theautonomous
autonomousrun
runstarts
startsatatn ==128,
128,asasindicated
indicatedbyby
the
arrow.Patterns
Patterns
shorter
than
are
as
the
arrow.
shorter
than
13 13
are
−3 . .In
reproduced with
with
low
MSE
10
while
patterns
longer
than
14
are
not
generated
correctly
with
MSE
10
the
latter
reproduced
10
while
patterns
longer
than
14
are
not
generated
correctly
with
MSE
10
the
latter
reproduced
with low
lowMSE
MSE< 10
10−3 ,, ,,while
whilepatterns
patternslonger
longerthan
than14
14are
arenot
notgenerated
generatedcorrectly
correctly
with
MSE
10
. In
In
the
latter
reproduced
with
low
MSE
with
MSE
> 10
. In
the
latter
cases, the
the
reservoir
dynamics
remains
stable
and
periodic,
but
the
output
only
remotely
resembles
the
target
pattern.
cases,
remains
stable
and
periodic,
but
the
output
only
remotely
resembles
the
target
pattern.
cases,
the reservoir
reservoirdynamics
dynamicsremains
remainsstable
stableand
andperiodic,
periodic,but
butthe
theoutput
outputonly
onlyremotely
remotely
resembles
the
target
pattern.
cases,
the
reservoir
dynamics
resembles
the
target
pattern.

0.6

0.9

Target
output
Target
Targetoutput
output
RC
output
RC
RCoutput
output

Output
signals
(arb.
units)
Outputsignals
signals(arb.
(arb.units)
units)
Output

Output
signals
(arb.
units)
Outputsignals
signals(arb.
(arb.units)
units)
Output

0.9

0.3
0

−0.3

Target
output
Target
output
Target
output
RC
output
RC
output
RC
output

0.6
0.3
0

−0.3
−0.6

−0.6
Autonomous
Autonomousrun
runstarts
starts
Autonomous run starts

50
100
50
100
50
100
Discrete
Discretetime
time n
Discrete time

150
150
150

Figure
Figure 9.
9. Example
Example of
of an
an output
output signal
signal for
for random
randompattern
pattern
Figure 9. Example of an output signal for random pattern
generation
reservoir
isis
generation task,
task, with
with aapattern
patternofoflength
length10.
10. The
Thereservoir
reservoiris
generation task, with a pattern of length 10. The reservoir is
first
section
firstdriven
drivenby
bythe
thedesired
desiredsignal
signalfor
for128
128timesteps
timesteps(see
(seesection
section
first driven by the desired signal for 128 timesteps (see section
IV
IVB),
B), and
and then
then the
the input
input isis connected
connectedto
tothe
theoutput.
output. Note
Note
IV B), and then the input is connected to the output. Note
that
that in
in this
this example
example the
the reservoir
reservoir output
output requires
requiresabout
about50
50
that in this example the reservoir output requires about 50
timesteps
timesteps to
to match
match the
the driver
driver signal.
signal. The
The autonomous
autonomousrun
run
timesteps to match the driver signal. The autonomous run
continues
continues beyond
beyond the
the scope
scopeofofthe
thefigure.
figure.
continues beyond the scope of the figure.

visualised
visualised on
on aa scope
scope and
and remained
remained stable
stableand
andaccurate
accurate
visualised
on a scope
and remained stable and accurate
through
through the
the whole
whole test.
test.
through the whole test.
The
above
results
were
The
obtained
by
scanning
the
input
The above
aboveresults
resultswere
wereobtained
obtainedby
byscanning
scanningthe
theinput
input
The
above
results
were
obtained
by scanning
the
input
gain
β
and
the
feedback
gain
α
to
obtain
the
best
results.
gain
and
the
feedback
gain
to
obtain
the
best
results.
gain and the feedback gain to obtain the best results.
gain
and
the
feedback
gain
to
obtain
the
best
results.
As
for
frequency
generation,
it
was
found
that
β
has
little
As
Asfor
forfrequency
frequencygeneration,
generation,ititwas
wasfound
foundthat
that has
haslittle
little
As for frequency
generation,
it was found
that
little
impact
on
the
system
performance
is
impact
so
long
as
itithas
isischochoimpact on
on the
the system
system performance
performance so
so long
long as
asit
choimpact
on
the
system
performance
so
long
as
it
is
chosen
in
the
interval
β
∈
[0.1,
1],
while
the
feedback
gain
sen
[0
sen in
in the
the interval
interval
[0 1],
1], while
while the
the feedback
feedback gain
gain
sen
in
the
interval
[0
1],
while
the
feedback
gain
α,
on
the
contrary,
has
to
lie
within
a
narrow
interval,
,, on
on the
the contrary,
contrary, has
has to
to lie
lie within
within aa narrow
narrow interval,
interval,
, on the contrary,
has0 to
lie within
approximately
between
approximately
and
110.0. a narrow interval,
approximately between
between 0.9
099 and
and1.0.
approximately between 0 9 and 1 0.

1950
1950
1950

1960
1960
1960

1970
1980
1970
1980
1970
1980
Discrete
Discretetime
timen
Discrete time

1990
1990
1990

2000
2000
2000

Figure
10.
ofofan
an
autonomous
run
output
after
1950
Figure10.
10. Example
Exampleof
anautonomous
autonomousrun
runoutput
outputafter
after1950
1950
Figure
Figure 10.Example
Example of
an autonomous
run output
after 1950
timesteps,
with
pattern
ofoflength
length
==14.
14.
The
RC
outputs
timesteps,with
withaaapattern
patternof
lengthL =
14.The
TheRC
RCoutputs
outputsaaa
timesteps,
timesteps, with a pattern of length = 14. The RC outputs a
periodic
signal
clearly
does
not
match
the
target
pattern
periodicsignal
signalthat
thatclearly
clearlydoes
doesnot
notmatch
matchthe
thetarget
targetpattern
pattern
periodic
periodic signalthat
that clearly
does not
match the
target pattern
−3
(MSE
==5.2
55 × 10
10
).).
(MSE=
10
(MSE
).
(MSE = 5
10 ).
D.
D. Mackey-Glass
Mackey-Glassseries
seriesprediction
prediction
Mackey-Glass
series
prediction
D.D. Mackey-Glass
series
prediction

Chaotic
Chaotictime
timeseries
seriesgeneration
generationtasks
taskswere
werethe
themost
mostafafChaotictime
timeseries
seriesgeneration
generation
tasks
were
the
most
Chaotic
tasks
thesurprising,
most
af-affected
by
noise.
This
isisnot
fected
bythe
theexperimental
experimental
noise.
Thiswere
not
surprising,
fectedbybythe
theexperimental
experimentalnoise.
noise.This
This
is
not
surprising,
fected
isvery
not
surprising,
since
sensitive
sincechaotic
chaoticsystems
systemsare,
are,by
bydefinition,
definition,
very
sensitivetoto
sincechaotic
chaoticsystems
systemsare,
are,byby
definition,
very
sensitive
since
definition,
very
sensitive
to to
noise.
noise. Reservoir
Reservoircomputing
computingwas
wasfirst
firstapplied
appliedtotothis
thisclass
class
noise.
Reservoir
computing
was
first
applied
to
this
class
noise.
Reservoir
computing
was firstwork,
applied
toauthors
this class
ofoftasks
inin[1].
the
intasks
[1]. In
Intheir
theirnumerical
numerical
work,
the
authors
inof
tasks
in
[1].
In
their
numerical
work,
the
authors
of
tasks
in
[1].
In
their
numerical
work,
the
authors
in-investigated
vestigatedthe
thecapacity
capacityofofthe
thecomputer
computertotofollow
followaagiven
given
vestigatedinthe
capacityofspace
of the
computer
follow
a given
vestigated
capacity
computer
toto
follow
a given
trajectory
of
attractor.
trajectorythe
in the
the phase
phase the
space
of the
the chaotic
chaotic
attractor.
trajectory
in
the
phase
space
of
the
chaotic
attractor.
trajectory
in the
phase
space but
of
chaotic
attractor.
We
this
approach,
since
our
We also
also tried
tried
this
approach,
butthe
since
ourexperimenexperimenWealso
alsotried
tried
thisapproach,
approach,
but
sinceour
our
experimenWe
this
but
since
tal
performs
as
emulator
ofofexperimenthe
talsystem
system
performs
asaa“noisy”
“noisy”
emulator
thechaotic
chaotic
tal
system
performs
as
a
“noisy”
emulator
of
the
chaotic
tal
system
performs
as
a
“noisy”
emulator
of
the
chaotic
attractor,
the
attractor, its
its trajectory
trajectory deviates
deviatesvery
veryquickly
quicklyfrom
from
the
attractor,
its
trajectory
deviates
very
quickly
from
the
attractor,
its
trajectory
deviates
very
quickly
from
the
target
target one,
one, especially
especiallywith
withaaSNR
SNRasaslow
lowasas40
40dB
dB(see
(see
target
one,
especially
with
a
SNR
as
low
as
40
dB
(see
target
one,
especially
with
a
SNR
as
low
as
40
dB
(see
section
sectionVVA).
A). For
Forthis
thisreason,
reason,we
weconsidered
considereddifferent
differentapapsectionVVA).
A).For
For
this
reason,
considered
different
apsection
this
reason,
wewe
considered
different
approaches
to
the
ofofthe
system,
as
proaches
to evaluate
evaluate
theperformance
performance
the
system,
as
proaches to evaluate the performance of the system, as

11

series
MGMG
series
RCRC
output
output

0

−0.2

Autonomous run MSE

0.2

10
Autonomous run MSE

Outputsignals
signals (arb.
(arb. units)
units)
Output

0.4

11
11

10−1

10

10−2

10

10−3

−0.4
10

run starts
MSE exceeds threshold
−0.6 Autonomous
Autonomous run starts
MSE exceeds threshold
0
100
200
300
400
500
600
0
100
200
300
400
500
600
Discrete time n
Discrete time

Figure11.
11.Example
Exampleof
of reservoir
reservoir computer
Figure
computeroutput
outputsignal
signal y(n)
Figure 11.
Example
of reservoir
computer
output
signal
(dotted
black
line) during
during
autonomous
run
(dotted
black
line)
autonomous
runon
onthe
theMackeyMackey(dotted
black
line)
during
autonomous
run
on
the
MackeyGlass task. The system was driven by the target signal
Glass
task.
The
system
was
driven
by
the
target
signal
Glassgrey
task.
system
was driven
byleft
therunning
target ausignal
(solid
line)The
for 128
timesteps
and then
(solid
grey
line)
for
128
timesteps
and
then
left
running
(solid grey line)
timesteps
then left
running
autonomously,
with for )128
coupled
to theand
input
) (see
Mainautonomously,
y(n)
coupled
(see
text).
The with
MSE
threshold
was to
set
to
10input
. I(n)
The )photonic
tonomously,
with
) coupled
tothe
theinput
(seeMain
Main
−3
text).
MSE
threshold
reservoir
computer
with
=was
600
was
to .generate
up to
text).The
The
MSE
threshold
wasset
settoable
to1010
. The
Thephotonic
photonic
reservoir
computer
N ==
600
was
435
correct
values. with
reservoir
computer
with
600
wasable
abletotogenerate
generateup
uptoto
435435
correct
values.
correct values.

will be described below.
Thebesystem
was below.
trained
over 1500 input
samples
andas
proaches
to
evaluate
the performance
of the
system,
will
described
wasbe
running
autonomously
6001500
timesteps.
In partic-and
will
described
below.
The
system
was
trainedfor
over
input samples
ular,
we
prepared
2100
steps
of
was running
autonomously
forthe
600Mackey-Glass
timesteps.
Inseries
particThe
system
was
trained
over
1500
input samples
and
for
each
run
of
the
experiment
and
used
the
first
1500
we prepared
2100 steps
of the
Mackey-Glass
series
wasular,
running
autonomously
for 600
timesteps.
In particas
teacher
) to train the
system
last
forawe
each
runsignal
of the
and
usedand
the the
first
1500
ular,
prepared
2100experiment
steps of the
Mackey-Glass
series
600 both as an initialisation sequence (see section IV B)
aseach
a teacher
signal
) to train
theused
system
and
the
last
forand
run
of
the
experiment
and
the
first
1500
as a target signal
) to compute the MSE of the
600
both as
an initialisation
sequence
(see and
section
IV
B)
asoutput
a teacher
signal
u(n)
to
train
the
system
the
last
signal
). These 2100 samples were taken from
and
as
a
target
signal
)
to
compute
the
MSE
of
the
600
both starting
as an initialisation
several
points (seesequence
equation (see
(7)) section
in orderIVtoB)
output
signal
). These
2100
samples
were
taken
from
and
asthe
a target
signal
d(n) on
to
compute
the
MSE
the
test
reservoir
computer
different
instances
ofofthe
several
starting
points
(see
equation
(7))
in
order
output
signal
y(n).
These
2100
samples
were
taken
from
Mackey-Glass series. We scanned the input gain and the to
test the
reservoir
computer
instances
of to
the
several
starting
points
(see on
equation
(7))
in order
feedback
attenuation
( t and
indifferent
equations
(1))
to find
Mackey-Glass
series.
scanned
the reservoir
input
gain
the
optimal
dynamics
of theWe
opto-electronic
forand
this
test
the reservoir
computer
on different
instances
of
the
feedback
and
equations
(1))
to
find
task.
We attenuation
used
[0 ( scanned
3] and in
tuned
the gain
optical
atMackey-Glass
series. We
the
input
and
the
optimal
dynamics
of
the
opto-electronic
reservoir
for
this
tenuator
in
the
range
[4
25
25]
dB,
which
corresponds
feedback attenuation (β and α in equations (1)) to find
task.
We
used
3]
and
tuned
the
optical
atapproximately
to
[0
85
95],
with
slightly
different
optimal dynamics of the opto-electronic reservoir for this
tenuator
in
the
range
[4
25
25]
dB,
which
corresponds
values
for
different
instances
of
the
Mackey-Glass
series.
task. We used β ∈ [0.1, 0.3] and tuned the optical atapproximately
[0 85of the
95],reservoir
with slightly different
Fig. 11 shows to
an example
tenuator
in the range
[4.25, 5.25]
dB, whichoutput
corresponds
(dotted
black
line) during
theofautonomous
run. The
values for
different
instances
the Mackey-Glass
series.
approximately to α ∈ [0.85, 0.95], with slightly different
target
series is shown
grey. The
MSE
Fig.Mackey-Glass
11 shows an example
of the in
reservoir
output
values
for different
instances
of the
Mackey-Glass
series.
reservoir computer
threshold
was set
to 10
(dotted black
line)
duringand
thetheautonomous
run. The
Fig.
11
shows
an
example
of
the
reservoir
output
predicted
435 correct values
this example.
12
target Mackey-Glass
series in
is shown
in grey. Fig.
The y(n)
MSE
(dotted
black
line)
during
the
autonomous
run.
The
displays
thewas
evolution
MSEthe
recorded
during
the
reservoir
computer
threshold
set to of10the and
target
Mackey-Glass
series
isplotted
shown
inexample.
grey.
MSE
same
autonomous
run.
The
error
curveThe
was
av- 12
Fig.
predicted
435 correct
values
in this
−3
threshold
was
set
to
10
and
the
reservoir
computer
eraged
over
200-timestep
intervals.
It
exceeds
the
10
displays the evolution of the MSE recorded during the
threshold
within
[500The
600]
and
a constant
predicted
435
correctrun.
values
inplotted
this reaches
example.
Fig.
same autonomous
error curve
was 12
avvalue
ofthe
approximately
10
timesteps.
displays
evolution
of1 the
MSEafter
recorded
during
the
eraged
over
200-timestep
intervals.
It2500
exceeds
the 10
At this
point, the run.
generated
time series
is completely
offavsame
autonomous
The
error
curvea was
threshold
within
[500 plotted
600] and
reaches
constant
−3
the target
(see
Fig. 13 forintervals.
illustration).
eraged
over
200-timestep
It
exceeds
the
10
value of approximately 1
10 after 2500 timesteps.
The noise
inside
the
opto-electronic
reservoir,
disthreshold
within
n
∈
[500,
600]
and
reaches
a
constant
At thisinpoint,
the generated the
time series isofcompletely
off
cussed
section V A, makes
expervalue
approximately
10−1outcome
after 2500an
timesteps.
the of
target
(see Fig. 131.1
for×illustration).
At this
generated
time series is completely
Thepoint,
noisethe
inside
the opto-electronic
reservoir, off
disthecussed
targetin(see
Fig. V
13A,
formakes
illustration).
section
the outcome of an experThe noise inside the opto-electronic reservoir, discussed in section V A, makes the outcome of an exper-

10−4
0

0

1000

2000
3000
4000
1000
2000
3000
Discrete time
Discrete time n

5000
4000

5000

Figure 12. Evolution of MSE during experimental auFiguregeneration
12. Evolution
Evolution
of MSE
MSE during
during
experimental
auFigure
12.
of
autonomous
of the Mackey-Glass
chaoticexperimental
time series
tonomous
generation
of the
the
Mackey-Glass
(same
run as in
Fig. (11). of
The
error
curve, averaged
overtime
200 series
tonomous
generation
Mackey-Glass
chaotic
time
series
(samerun
run
asin
in
Fig.
(11).
The error
error
curve, averaged
averaged
over
timesteps,
crosses
the
10 (11).
threshold
approximately
between
(same
as
Fig.
The
curve,
over 200
200
timesteps,
crosses
the 10
10−3 threshold
threshold approximately
approximately between
=
500 and crosses
=
600. the
timesteps,
between
500and
andn ==600.
600.
n ==500

iment inconsistent. That is, several repetitions of the
experiment
with same parameters
result
in signifi-of
iment inconsistent.
inconsistent.
That is,
is, may
several
iment
That
several
repetitions
of the
the
cantly
differentwith
prediction
In fact,
impact
experiment
with
same lengths.
parameters
may the
result
in
experiment
same
parameters
may
result
in signifisignifiof cantly
noise
varies
from prediction
one
trial tolengths.
another. In
In fact,
some the
casesimpact
cantly
different
prediction
lengths.
In
fact,
different
the
impact
it does
not disturb
the system
much.
But in most
cases cases
of
noise
varies
from
one
trial
to
another.
In
some
of noise varies from one trial to another. In some cases
it induces
anot
significant
error
on themuch.
outputBut
value most
), cases
doesnot
disturb the
the system
system
ititdoes
disturb
much.
But in
in most
cases
so that
the
neural
network
deviates
very
quickly
from
itinduces
induces aa significant
significant error
error on
on the
the output
output value
),
value
y(n),
theitso
target
trajectory.
To
estimate
the
variability
of
the
that the
the neural
neural network
network deviates
deviates very
very quickly
from
so
that
quickly
from
results,
we performed
50 consecutive
runs of the
thetarget
target
trajectory.
To estimate
estimateautonomous
the variability
variability
the
To
with
the sametrajectory.
readout weights
and thethe
same optimal of the
results, we
we performed
performed 50 consecutive
consecutive autonomous
autonomous runs
results,
experimental
parameters. 50
While the system
produced runs
with
the
same
readout
weights
and
the
optimal
with very
the good
samepredictions
readout weights
the same
same
optimal
several
(of orderand
of 400),
most of
experimental
parameters.
While
the
system
produced
system
produced
theexperimental
outcomes wereparameters.
rather poor, While
with anthe
average
predicseveral
very
good
predictions
(of order of 400),
most of
tion
lengthvery
of 63good
7 andpredictions
a standard deviation
2. We
several
(of order ofof65400),
most of
the
outcomes
were
rather
poor,
with
an
average
predicobtained
similar behaviour
with
the noisy
experimental
the outcomes
were rather
poor,
with an
average prediction
length
of
63
7
and
a
standard
deviation
of
65
2. We
model,
as in the
experiments.
tion using
lengththeofsame
63.7level
and of
a noise
standard
deviation
of 65.2. We
obtained
similar
behaviour
with theinnoisy
experimental
Changing
the
ridge
regression
parameter
the
training
obtained similar behaviour with the noisy experimental
model,
using
theII)
same
level
of noise as in
the experiments.
process
(see
section
didlevel
not improve
This
model,
using
theridge
same
of noise the
as inresults.
the experiments.
Changing
the reservoir
regression
the training
suggests
that the
computerparameter
emulates ain“noisy”
Changing
thesection
ridge regression
inresults.
the training
process (see
II) therefore,
did notparameter
improve
Mackey-Glass
system, and
using itthe
to follow This
process
(see
section
II)
did
not
improve
the
results.
This
suggests
that the
reservoir
a “noisy”
a given
trajectory
doesn’t
make computer
much senseemulates
with such
a
suggests
that
the
reservoir
computer
emulates
a
“noisy”
system, and
it to follow
highMackey-Glass
noise level. Nevertheless,
thetherefore,
noise doesusing
not prevent
system,
and
therefore,
it to–such
follow
theMackey-Glass
from emulating
themake
Mackey-Glass
system
asystem
given trajectory
doesn’t
much using
sense
with
a
ahigh
given
trajectory
doesn’t
make
much
sense
with
such a
even
if the
output
deviates
from
the
target,
noise
level.quickly
Nevertheless,
the
noise
does
not it
prevent
high
noise
level.
Nevertheless,
theMackey-Glass
noise
does notsystem
prevent
still
resembles
the
original
time the
series.
Therefore,
we
the
system
from
emulating
–
the
system
from
emulating
the
Mackey-Glass
system
tried
a
few
distinct
methods
of
comparing
the
output
of
even if the output quickly deviates from the target, it–
ifresembles
the
quickly
from Therefore,
the target,we
it
theeven
system
withoutput
thethe
target
time deviates
series.
still
original
time series.
still
resembles
time
series.
Therefore,
We
performed
a the
new original
set
of experiments,
where,
tried
a few distinct
methods
of comparing
theafter
output we
of
tried
a few
distinct
of the
comparing
therunoutput of
a training
phase
of 1500
timesteps,
system was
the system
with
themethods
target time
series.
theautonomously
system
with
the
target
time
series.
in
order
to
collect
ning
for
10k
timesteps
We performed a new set of experiments, where, after
enough
points
for
dataaofanalysis.
We
then proceeded
with
performed
new
of experiments,
where,
a We
training
phase
1500set
timesteps,
the system
was after
runa simple
visual
inspection
of
the
generated
time
series,
to runaning
training
phase of 1500
timesteps,
theinsystem
was
autonomously
for 10k
timesteps
order to
collect
check
whether
it still looks
similar
to the Mackey-Glass
ning
autonomously
for analysis.
10k
timesteps
in order
to collect
enough
points for data
We then
proceeded
with
time series, and does not settle down to simple periodic
enough
points
for
data analysis.
We
then proceeded
with
a
simple
visual
inspection
of
the
generated
time
series,
to
oscillations. Fig. 13 shows the output of the experimenacheck
simplewhether
visual inspection
of similar
the generated
series, to
it still looks
to the time
Mackey-Glass
check
whether
it does
still looks
similar
to to
thesimple
Mackey-Glass
time series,
and
not settle
down
periodic
time
series, and
not settle
down to
periodic
oscillations.
Fig.does
13 shows
the output
of simple
the experimen-

oscillations. Fig. 13 shows the output of the experimental reservoir computer at the end of the 10k-timestep

12

0.2
0
−0.2

timestep
output is
irregular
ations.
me series
spectra.
the origie output
ack line).
very ac, with its
es.
nt of the
ed in the
0 01 for
lue comass series
unov ex-

0

0 05

Lorenz series prediction
01

0 15

02

0 25

0

This task was investigated in a similar way to the
Frequency
9700
9800
9900
10000
previous one. The reservoir computer was trained over
Discrete time n
3000 input samples and was run autonomously for 1000
Figure 14. Comparison of Fast Fourier Transforms of the
timesteps.
The 4000
were
from
original
Mackey-Glass
series samples
(solid grey
line)taken
and the
timean interFigure
13.
Output
of
the
experimental
reservoir
computer
val
with
even
distribution
of
transitions
between
Figure 13. Output of the experimental reservoir computer series generated by the photonic reservoir computer (dashed the two
(dotted black line) at the end of a long run of 10k timesteps.
“wings”
ofplot
theis Lorenz
In as
fact,
we have noline). The
limited toattractor.
low frequencies
the power
(dotted black line) at the end of a long run of 10k timesteps. black
Although the system does not follow the starting trajectory
at higher
frequencies
is almost
null.
Dominant
frequencies
Although
the
system
does
not
follow
the
starting
trajectory
ticed
that
the
first
1000
samples
of
the
sequence
gener(solid grey line), its output still resembles visually the target
multiples
of 1/τ
(seesection
section III
(solid
line), its output still resembles visually the target correspond
ated byto the
ode45
solver 06
(see
IIIC).
D)The
contained
timegrey
series.
reproduces the target spectrum notably well.
12experiment
time series.
more oscillations above zero than below, that is, a transient from the starting point to the actual chaotic attractal reservoir
computer at the end of the 10k-timestep
tor.
uneven
distribution
the reservoir
8
tor.
ThisThis
uneven
distribution
forced forced
the reservoir
com- comseries output is
autonomous run. It shows that theMG
reservoir
puter
to
generate
a
biased
output.
We
thus
discarded
the
puter to generate a biased output. We thus discarded the
RC output
still similar to the target time series,
that is, irregular
first
1000
values
and
trained
the
system
over
the
interval
first 1000 values and trained the system over the interval
and consisting of the same kind of uneven oscillations.
[1000, 4000]
(these
transients
wereremoved
also removed
6
[1000
(these
initialinitial
transients
were also
A more thorough way of comparing two time series
in [1]).
For For
optimal
performance
of the of
opto-electronic
in [1]).
optimal
performance
the opto-electronic
that “look similar” is to compare their frequency spectra.
reservoir,
we set
gain to
theand
feed-the feedreservoir,
wethe
setinput
the input
gain=to0 5β and
= 0.5
Fig. 14 shows the Fast Fourier Transforms of the origiback
attenuation
to =
6α1 =
dB.
4
back
attenuation
to
6.1
dB.
nal Mackey-Glass series (solid grey line) and the output
Fig.Fig.
15 shows
an example
of the of
reservoir
output output y(n)
15 shows
an example
the reservoir
of the experiment after a long run (dotted black line).
(dotted
black
line)
during
the
autonomous
run.
The
tar-The tar(dotted black line) during the autonomous
run.
Remarkably, the reservoir computer reproduces very acget Lorenz series is shown in grey. With the MSE thresh2
curately
the spectrum of the chaotic time series, with its
get Lorenz series is shown in grey. With the MSE threshold set to 10 , the
−3 system predicted 122 correct steps,
main frequency and several secondary frequencies.
old settwo
to 10
, the system
correct
including
transitions
betweenpredicted
the wings 122
of the
at- steps,
Finally, we estimated the Lyapunov exponent of the
including
two
transitions
between
the
wings
of the attractor. As in the Mackey-Glass study, we performed
generated
time series, using the method described in the
0
tractor.
As
in
the
Mackey-Glass
study,
we
performed
50 autonomous runs with identical readout weights and
0
0.05Material
0.1 of [1].
0.15 We 0.2
0.250 01 0.3
Supplementary
obtained
for
50optimal
autonomous
runs with
identicalanreadout
and
same
parameters,
and obtained
average weights
pre−1
our experimental implementation,
Frequency t while the value comsamehorizon
optimal
parameters,
an deviaverage prediction
of 46
0 timestepsand
withobtained
a standard
monly found in the literature for the Mackey-Glass series
ation
of 19 5.
Takingofinto
the higher
of devidiction
horizon
46.0account
timesteps
with adegree
standard
is 0 006. The slightly higher value of the Lyapunov exchaos
of the
Lorenz
attractor,
and
given the
same
probFigure
14.
Comparison
of
Fast
Fourier
Transforms
of
the
ation
of
19.5.
Taking
into
account
the
higher
degree of
Figure
14.
Comparison
of
Fast
Fourier
Transforms
of
the
ponent may simply reflect the presence of noise in the
lemschaos
related
to
noise,
it is attractor,
hard to expect
agiven
betterthe
perfororiginal
Mackey-Glass
series
(solid
grey
line)
and
the
time
of
the
Lorenz
and
same
proboriginal
Mackey-Glass
series
(solid
grey
line)
and
the
time
emulator.
the reservoir
computer
at following
target perforseries generated
generated by
by the
thephotonic
photonicreservoir
reservoircomputer
computer(dashed
(dashedmance
series
lemsofrelated
to noise,
it is hard
to expectthe
a better
depicts computer
the evolution
of the MSE
black line).
line). The
Theplot
plotisislimited
limitedtotolow
lowfrequencies
frequenciesasasthe
thepower
powertrajectory.
black
mance ofFig.
the 16
reservoir
at following
the target
the autonomous run. The error curve was avat higher
higher frequencies
frequencies isis almost
almostnull.
null. Dominant
Dominantfrequencies
frequenciesduring
at
Fig.
16
depicts
the
evolution
of
the MSE
trajectory.
E. Lorenz series prediction
over 100-timestep intervals. The initial dip corcorrespond to
to multiples
multiplesofof1/τ
1/τ ≈ 0.06
06(see
(seesection
sectionIII
IIIC).
C).The
Theeraged
correspond
during
the
autonomous
run.
The
error
curve
resonds to the teacher-forcing of the reservoir computer was avexperiment reproduces
reproducesthe
thetarget
targetspectrum
spectrumnotably
notablywell.
well.
experiment
over signal
100-timestep
intervals.asThe
initialindip corwitheraged
the target
for 128 timesteps,
discussed
This task was investigated in a similar way to the
to the
the reservoir
computer
around
Sec.resonds
IV A. The
errorteacher-forcing
exceeds the 10 of threshold
previous one. The reservoir computer was trained over
signal
for 128
timesteps,
in
the with
= the
250 target
mark and
reaches
a constant
valueasofdiscussed
ap3000 input samples and was run autonomously for 1000
−3
tor.
This
uneven
distribution
forced
the
reservoir
comautonomous
run.
It
shows
that
the
reservoir
output
is
Sec.
IV
A.
The
error
exceeds
the
10
threshold
around
proximately
1
10
after
less
than
1000
timesteps.
timesteps. The 4000 samples were taken from an interputer
to generate
biasedoftime
output.
We between
thus discarded
still
to distribution
thea target
series,
that
is, the
irregular
this npoint,
themark
reservoir
lost thevalue
tar- of ap= 250
andcomputer
reaches ahas
constant
valsimilar
with
even
transitions
two theAt the
first
1000
values
and
trained
the
system
over
the
interval
and
consisting
of
the
same
kind
of
uneven
oscillations.
get proximately
trajectory, but1.5
keeps
on−2
generating
a time
“wings” of the Lorenz attractor. In fact, we have no× 10
after less
thanseries
1000with
timesteps.
(these
initial
removed
[1000
thereservoir
Lorenz series
(see Fig.
forthe tarticed
thatthorough
the
first 1000
of thewere
sequence
generA more
waysamples
oftransients
comparing
twoalso
time
series properties
At thissimilar
point,tothe
computer
has 17
lost
byFor
the
(see section
III D)
contained
inated
[1]).
optimal
performance
of frequency
the
opto-electronic
that
“look
similar”
issolver
to
compare
their
spectra. illustration).
get trajectory, but keeps on generating a time series with
more
above
zero
thanTransforms
below,
is, the
a the
tranto thesimilar
Mackey-Glass
task, we series
performed
vi- 17 for
reservoir,
we set
input
gain
to
= that
0 5 and
feed- Similar
Fig.
14oscillations
shows
thethe
Fast
Fourier
of
origiproperties
to the Lorenz
(see aFig.
sient
from the starting
the
attracinspection of the generated Lorenz series after a long
back
attenuation
to
=(solid
6 1todB.
nal
Mackey-Glass
seriespoint
greyactual
line) chaotic
and the
output sualillustration).
15 shows anafter
example
of run
the reservoir
output
of Fig.
the experiment
a long
(dotted black
line).
Similar to the Mackey-Glass task, we performed a viRemarkably,
reservoir
reproduces
(dotted blackthe
line)
during computer
the autonomous
run. very
The actar- sual inspection of the generated Lorenz series after a long
curately
theseries
spectrum
of the
timethe
series,
with
its
get Lorenz
is shown
in chaotic
grey. With
MSE
threshrun, and compared the frequency spectra. Fig. 17 shows
main
frequency
theseveral
systemsecondary
predictedfrequencies.
122 correct steps, the output of the experiment near the end of a 95k auold set
to 10 ,and
including
transitionsthebetween
the exponent
wings of of
thethe
at- tonomous run. Although the system is quite far from the
Finally, two
we estimated
Lyapunov
tractor. As in the Mackey-Glass study, we performed
50 autonomous runs with identical readout weights and
same optimal parameters, and obtained an average prediction horizon of 46 0 timesteps with a standard deviation of 19 5. Taking into account the higher degree of
chaos of the Lorenz attractor, and given the same probPower (arb. units)

computer
imesteps.
trajectory
the target

E.

−0.4
−0.6
9600

10000

generated time series, using MG
theseries
method described in the
output
Supplementary Material ofRC[1].
We obtained 0.01 for
our experimental implementation, while the value commonly found in the literature for the Mackey-Glass series
is 0.006. The slightly higher value of the Lyapunov exponent may simply reflect the presence of noise in the
emulator.

Power (arb. units)

Output signals (arb. units)

0.4

MG series
RC output

12

0.4
0.2
0

−0.2
Autonomous
Autonomousrun
runstarts
starts MSE
MSEexceeds
exceedsthreshold
threshold
Autonomous
run
starts
MSE
exceeds
threshold
Autonomous 100
run
MSE exceeds
300
400
500
100starts200
200
300 threshold
400
500

−0.4

100
100

200 Discrete
300
time
Discrete
time
200
300

400
400

RC output

0.2

0

−0.2

−0.4

500
500

93600
93600

93600
93600

Discrete time
time n
Discrete

Figure
Figure15.
15.Example
Exampleofofreservoir
reservoircomputer
computeroutput
output signal
signal
(dotted
(dottedblack
blackline)
line)during
duringautonomous
autonomous runs
runs on
on the
the Lorenz
Lorenz
Figure15.
15. Example of
of reservoir computer
computer output
output signal
signal y(n)
Figure
Figure
15.
task.
The
was
task.Example
Thesystem
systemreservoir
wasdriven
drivenby
bythe
thetarget
target signal
signal (solid
(solid grey
grey
(dotted
black
line)
during
autonomous
runs
on
the
Lorenz
(dotted
black
line)
during
runs
on the Lorenz
(dotted
black
line)
during autonomous
autonomous
runs
line)
for128
128timesteps
timesteps
beforerunning
running
autonomously
(see secsecline)
for
before
autonomously
(see
task.The
The
system
was
driven
bythe
thewas
target
signal
(solid
grey
task.
system
was
driven
by
target
signal
task.
The
system
was
driven
by
the
target
signal
grey
tion
IV
A).The
The
MSE
threshold
was
setto
to
10 ..(solid
The photonic
photonic
tion
IV
A).
MSE
threshold
set
10
The
line)for
for
128timesteps
timesteps
before
running
autonomously
(see
sec-exline)
for
128
timesteps
before
running
autonomously
secline)
128
before
running
autonomously
(see
system
with ==
600generated
generated
122
correctvalues
values
in this
this
exsystem
with
600
122
correct
in
−3 . The photonic
tion
IV
A).
The
MSE
threshold
was
set
to
10
ample,
and
predicted
twoswitches
switches
thetrajectory
trajectory
from one
one
tionIV
IV
A).The
The
MSE
threshold
was
setto
to
10
tion
A).
MSE
threshold
was
set
.. The photonic
ample,
and
predicted
two
ofof10
the
from
system
with
=attractor
600generated
generated
122
correct values
values in
in this
this exexlobe
the
attractor
the122
other.
system
with
600
generated
122
correct
values
lobe
ofof
totothe
other.
system
with
Nthe
==
600
correct
ample,and
andpredicted
predictedtwo
twoswitches
switchesof
ofthe
the trajectory
trajectory from
from one
one
ample,
and
predicted
two
switches
of
the
trajectory
ample,
lobeofof
ofthe
theattractor
attractorto
tothe
theother.
other.
lobe
the
attractor
to
the
other.
lobe

LZ series
LZ
series
LZoutput
series
LZ
series
RC
RC
RC output
output

0.4

Output
Output signals
signals (arb.
(arb. units)
units)
Outputsignals
signals(arb.
(arb.units)
units)
Output

LZ series
LZ series
LZRC
series
output
LZ
series
RC output
RC
output
RC output

Output signals (arb. units)
Output signals (arb. units)

Output
Outputsignals
signals(arb.
(arb. units)
units)

13
13
1313
13

10−2
10

10−3
10

1010

10
10

4

LZ
LZseries
series
RC
RCoutput
output

0

0

0

500
500

500

1000
1000
Discrete time
Discrete time

1000

1500
1500

1500

2000
2000

2000

0 16.
500
1500 experimental
2000
Evolution
of1000
MSE during
auFigure
Evolutionof
of MSE
experimental auFigure
16. generation
Discrete
time during
tonomous
the
Lorenz
Discrete
time
n chaotic time series (same
tonomous
of the
time series
(same
run as ingeneration
Fig. (15).
TheLorenz
error chaotic
curve, averaged
over
100
run
as in Fig.
(15).
curve,near
averaged
over The
100
timesteps,
crosses
the The
10 error
threshold
= 250.
Evolution
of 10
MSEthreshold
during experimental
experimental
auFigure
16. Evolution
timesteps,
crosses
the
near
250.auThe
Figure
16.
Evolution
of
MSE
during
Figure
16.
of
MSE
initial
dip
corresponds
to theduring
warmupexperimental
of the =
reservoir
(see
tonomous
generation
ofthe
theLorenz
Lorenz
chaotic time
time
series
(same
initial
dipA).
corresponds
to
the chaotic
warmup
of the
reservoir
(see
tonomous
generation
of
the
Lorenz
chaotic
time
tonomous
generation
of
series
(same
Sec.
IV
runasas
as
inFig.
Fig.
(15). The
The error
error curve,
curve, averaged
averaged over
over 100
100
Sec.
IV
A).(15).
run
Fig.
(15).
The
error
curve,
averaged
over
run
inin
100

timesteps,crosses
crossesthe
the 10
10−3 threshold
threshold near
near n =
= 250.
250. The
The
timesteps,
crosses
the
10
threshold
near
=
250.
timesteps,
The
initial dip
dip corresponds
corresponds to
to the
the warmup
warmup of
of the
the reservoir
reservoir (see
(see
initial
initial dip
corresponds to the warmup of the reservoir (see
Sec. IV
IVrun,
A). and compared the frequency spectra. Fig. 17 shows
Sec.
A).
Sec. IVrun,
A).
compared
the frequency
spectra.
Fig.
shows
the and
output
of the experiment
near
the end
of a1795k
authe
output of
theAlthough
experiment
near the
end offara from
95k autonomous
run.
the system
is quite
the
tonomous
run. Although
quite
far from
target trajectory
(plottedthe
in system
grey) atisthis
point,
it is the
aprun, and
and
compared
the
frequency
spectra.
Fig.
17the
shows
target
trajectory
(plotted
in grey)
at this
point,
it is apparent
that(plotted
itthe
has
captured
of17
run,
compared
frequency
spectra.
Fig.
target
trajectory
in grey)the
at dynamics
this point,
it shows
is Lorenz
apthe output
output
of
theitexperiment
experiment
near
the
endthe
ofofaaFast
95kLorenz
auparent
that
has
thedisplays
dynamics
the
system
very
well.captured
Fig. near
18
Fourier
the
of
the
the
end
of
95k
auparentsystem
that itvery
has well.
captured
the
dynamics
of
the
Lorenz
Fig.
18
displays
the
Fast
Fourier
tonomous
run.
Although
the
system
is
quite
far
from
the
Transforms
of
the
generated
time
series
(dotted
black
tonomous
run.
Although
the
system
is quite
far from
the
system
very
well.
Fig.
18
displays
the
Fast
Fourier
Transforms
of
the
generated
time
series
(dotted
black
line)
and
the
computed
Lorenz
series
(solid
grey
line).
target
trajectory
(plotted
in
grey)
at
this
point,
it
is
aptarget trajectory
(plotted
in grey)
atseries
this point,
it is
apTransforms
ofitthe
the
generated
time
(dotted
black
line)
and
the
computed
Lorenz
series
(solid
line).
Unlike
system,
these
specparent
that
hasMackey-Glass
captured the
the
dynamics
offrequency
thegrey
Lorenz
parent
that
it
has
captured
dynamics
of
the
Lorenz
line)
and
thethe
computed
seriesthese
(solid
greyFourier
line).
Unlike
Mackey-Glass
frequency
specsystem
very
well.
Fig. Lorenz
18 system,
displays
the Fast

systemthe
very well. Fig. system,
18 displays
the
Fast Fourier
Unlike
frequency
specTransformsMackey-Glass
of the
the generated
generated time
timethese
series
(dotted black
black
Transforms
of
series
(dotted
tra
do
not
have
any
dominant
frequencies.
That
is,
the
line)
and
the
computed
Lorenz
series
(solid
grey
line).
line) and the computed Lorenz series (solid grey line).
power
distribution
does
not
contain
any
strong
spikes,
Unlike the
the Mackey-Glass
Mackey-Glass system,
system, these
these frequency
frequency specspecUnlike
that could have been used as reference points for comparison. Therefore, comparing the two the spectra is

LZ series
series
LZ
RC output
output
RC

3

2

1

0
10−4
10

94000
94000

Discrete time
time
Discrete

Power
Power (arb.
(arb. units)
units)
Power
Power(arb.
(arb. units)
units)

Autonomous run MSE
Autonomous run MSE

Autonomous
Autonomousrun
runMSE
MSE

1010

94000
94000

93800
Discretetime
time
Discrete
n
93800

Figure17.
17.Output
Outputof
of the
the experiment
experiment (dotted
(dotted
black
line)
Output
of
the
line)
atatat
Figure
experiment
(dottedblack
black
line)
the
end
long
run
of
95k
timesteps
on
the
Lorenz
task.
of
a
long
run
of
95k
timesteps
on
the
Lorenz
task.
the
end17.
of Output
aOutput
long run
of 95k
timesteps
on the
Lorenz
Figure
17.
of the
the
experiment
(dotted
black
line)task.
at
Figure
of
(dotted
black
line)
at
the system
does
not
starting
trajectory
Although
system
doesexperiment
not follow
follow the
the
starting
trajectory
Although
system
does
nottimesteps
follow the
starting
trajectory
the(solid
end grey
ofthe
long
run
of
95k
timesteps
on
the
Lorenz
task.
the
end
of
aa line),
long
run
of
95k
on
the
Lorenz
task.
line), it
it does
does aa fairly
fairly good
good job
jobatatemulating
emulatingthe
the
(solid
greythe
line),
it
does
anot
fairly
good
at emulating
the
Although
the
system
does
not
follow
thejob
starting
trajectory
Although
system
does
follow
the
starting
trajectory
of
the
system.
dynamics
the Lorenz
Lorenz
system.
dynamics
the Lorenz
(solid grey
greyofline),
line),
it does
doessystem.
fairly good
good job
job at
at emulating
emulating the
the
(solid
it
aa fairly
dynamics of
of the
the Lorenz
Lorenz system.
system.
dynamics

1010

10−1
10

93800
93800

0
0

01
0.1

02
0.2

03
04
0.3
0.4
Frequency −1
Frequency t

05
0.5

0

0.6

00 18. Comparison
00 11
00 22 of Fast
00 33 Fourier
00 44Transforms
00 55 of the
00
Figure

Figure
18.
Comparison
of Fast
Fast
Fourier
ofofthe
Figure
Comparison
of
Transforms
the
Frequency
Lorenz18.
series
(solid grey Frequency
line)
and Fourier
the timeTransforms
series generated
Lorenz
series(solid
(solid
grey computer
line) and
and (dotted
the
series
generated
Lorenz
grey
line)
the time
time
series
generated
by theseries
photonic
reservoir
black
line)
during
by
photonicreservoir
reservoir
computer
(dotted
black
line)
95k
timesteps.
Both spectra
are normalised
so as
to during
have
by
thethe
photonic
computer
(dotted
black
line)
during
Figure
18. Comparison
Comparison
of
Fast
Fourier
Transforms
of
the
95k
timesteps.
Both
spectra
areFourier
normalised
so
as
toto
have
Figure
18.
of
Fast
Transforms
of
the
equal
total power.
The
curves
smoothened
by
averaging
95k
timesteps.
Both
spectra
are
normalised
so
as
have
Lorenz
series
(solid
grey
line)
and
the
time
series
generated
equal
total
power.
The
curves
are
smoothened
by
averaging
Lorenz
series
(solid
grey
line)
and
the
time
series
generated
over
50
samples
and
the
plot
is
limited
to
lower
frequencies
equal
power.
The curves
are (dotted
smoothened
by
averaging
byover
thetotal
photonic
reservoir
computer
black
line)
during
50 samples
and the
plot
is limited
toblack
lowerline)
frequencies
by
the
reservoir
computer
(dotted
during
(the
higher
ones
being
close
to zero).
Despite
mismatch,
over
50photonic
samples
and
the
plot
is
limited
to some
lower
frequencies
95k
timesteps.
Both
spectra
are
normalised
so
as
togrey
have
(the
higher
ones
being
close
to
zero).
Despite
some
mismatch,
95k
timesteps.
Both
spectra
are
normalised
so
as
to
have
the
shape
of
the
dotted
curve
is
roughly
similar
to
the
(the
higher
ones
being
to zero).
Despite
some
mismatch,
equal
total power.
power.
Theclose
curves
are
smoothened
bytoaveraging
averaging
the
shape
of the
dotted
curve
is roughly
similar
the grey
line.total
equal
The
curves
are
smoothened
by
the
shape
of the and
dotted
roughlytosimilar
to the grey
line.
over
50 samples
samples
the curve
plot is
isislimited
limited
lower frequencies
frequencies
over
50
and the
plot
to lower
line.
(the
higher
ones
being
close
to
zero).
Despite
some
mismatch,
(the higher ones being close to zero). Despite some mismatch,
the shape
shape of
of the
the dotted
dotted curve
curve is
is roughly
roughly similar
similar to
to the
the grey
grey
the
tra do not have any dominant frequencies. That is, the
line.
line.
That
is,
the
tra
do
not
have
any
dominant
frequencies.
power
distribution
does
any strongthe
spikes,
much
more
subjective
in not
thiscontain
case. Although
curves
power
distribution
does
notas contain
any
strong
spikes,
that
could
have
been
used
reference
points
for
comdo not match, one can still see a certain similitude bethat
could have beencomparing
used as reference
for comparison.
the two points
the spectra
is
tween
them.Therefore,
parison.
Therefore,
comparing
the
two
the
spectra
is
That
is, the
the
tra
do
not
have
any
dominant
frequencies.
much
more
subjective
in
this
case.
Although
the
curves
traIndoaddition
not have
any
dominant
frequencies.
That
is,
to one
those
visual
comparisons,
wethe
performed
much
more
subjective
in
this
case.
Although
curves
do
not
match,
can
still
see
a
certain
similitude
bepower
distribution
does
not
contain
any
strong
spikes,
power
distribution
does
not
contain
any strong
spikes,
athat
specific
randomness
test
ofasthe
generated
series.
exdo
notthem.
match,
one
still
see
a certain
similitude
between
could
have been
beencan
used
reference
points
forWe
comthat
could
have
used
as
reference
points
for
comploited
an
interesting
property
of
the
Lorenz
dynamics.
tween
them.
parison. Therefore, comparing the two the spectra is

parison.
Therefore,
comparing
the two
the spectra
is
Since
basically
switches
between
regions
(thecurves
wings
much it
more
subjective
in this
this
case. two
Although
the
much
more
subjective
in
case.
Although
the
curves
of
with
one tobethe
dothe
notbutterfly),
match, one
one
canrandom
still see
seetransitions
certain from
similitude
do
not
match,
can
still
aa certain
similitude
beother,
one
can
assign
binary
“0”
and
“1”
to
these
regions
tween them.
them.
tween
and thus transform the Lorenz series into a sequence of
random bits. We used this trick to check the randomness

14
of the generated series. We solved the Lorenz equation
and ran the experiment for 95k timesteps and converted
the resulting time series into two sequences of approximately 2400 bits. The two were then analysed with the
ENT program [41] – a well known software for testing
random number sequences – with the results shown in
Tab. I. Their interpretation requires a brief introduction
of the tests performed by the software.
• The first test computes the entropy per byte (8
bits). Since entropy is a measure of disorder, i.e.
randomness, a truly random sequence should have
8 bits of entropy per byte. Both sequences are close
to the maximum value, with the Lorenz series being
slightly more random.
• The compression is a commonly used indirect
method of estimating the randomness of bytes in
a file by compressing it with an efficient compression algorithm (such as e.g. Lempel-Ziv-Renau
algorithm, used by the Zip program). These algorithms basically look for large repeating blocks,
that should not appear in a truly random sequence.
Again, both sequences could only be slightly compressed.
• The mean value is the arithmetic mean of the data
bytes. A random sequence should be evenly distributed, and thus have a random value of 127.5.
The Lorenz series is very close to this value, and
the RC sequence is fairly close.
• The Monte Carlo method of computing the value
of π randomly places points inside a square and
computes the ratio of points located inside an inscribed circle, that is proportional to π. This is a
more complex test that requires a large sequence to
yield accurate results. We note that, nevertheless,
both sequences produce a reasonable estimation of
π.
• Finally, the serial correlation coefficient of a truly
random sequence is zero. Both series present a very
low correlation, yet again the Lorenz series demonstrating a better score.
These results do not prove that the generated sequence
is random. One obviously has to use a much longer sequence of bits for this study, and should also consider
more sophisticated and complete tests, such as Diehard
[42] or NIST Statistical Test Suite [43]. The purpose of
these tests was to show that the output of the RC generator does not consist of trivial oscillations, that only remotly resemble the Lorenz system. The results show that
the randomness of the RC output is similar to the Lorenz
system, which gives reasons to believe in the similarity
between the two time series. This, in turn, indicates, that
our photonic reservoir computer was capable of learning
to effectively emulate the dynamics of the Lorenz chaotic
system.

RC Lorenz
Entropy (/byte)
6.6
7.1
Compression (%) 17
10
Mean (byte)
134.3 125.8
π
2.88
3.00
Correlation
−0.08 −0.02
Table I. Results returned by the ENT program for the bit sequences generated by the experiment (RC) and the integrated
Lorenz system. The Lorenz sequence shows better figures, but
the RC output is not far behind. All these figures are poor
compared to common random series, but this is due to the
very short sequences used here (roughly 300 bytes).

VI.

CONCLUSION

The present work demonstrates the potential of output feedback in hardware reservoir computing and constitutes an important step towards autonomous photonic
implementations of recurrent neural networks. Specifically, we presented a photonic reservoir computer that
is capable of generating both sine waves of different frequencies and short random patterns with substantial stability. Moreover, it could emulate the Mackey-Glass time
series with a remarkably similar frequency spectrum, and
a fairly close highest Lyapunov exponent. Finally, it
could efficiently capture the dynamics of the Lorenz system and generate a time series with similar randomness
properties. To the best of our knowledge, this is the first
report of these task being implemented on an experimental reservoir computer.
The readout of the reservoir computer is carried out
in real time on a fast FPGA chip. This results in a digital output layer in an analogue device, that nevertheless
allows one to investigate many of the issues that will
affect a system with purely analogue feedback. The latter is a much more complicated experiment. Indeed the
only analogue output layers implemented so far on experimental reservoir computers were reported in [28–30].
Using them for output feedback would require adding an
additional electronic circuit consisting of a sample and
hold circuit, amplification, and multiplication by the input mask. The present experiment allows one to investigate the benefits that output feedback has to offer to
experimental reservoir computing, while anticipating the
difficulties and limitations that will affect a fully analogue
implementation. Such a two-step procedure, in which
part of the experiment is analogue and part digital, is
natural, and parallels the development of experimental
reservoir computers in which some of the first experiments were only partially analogue, see e.g. [7, 10].
This work allowed to highlight a critical limitation of
the present opto-electronic setup, namely the relatively
high level of noise generated by the components. While
this was not a concern in previous experiments with offline readout [8, 29], it becomes critical in this study of
the output feedback. Since this noise can not be averaged out, it propagates back into the system with out-

15
put feedback and considerably deteriorates the reservoir
states. This problem does not have a simple solution.
One could rethink the entire experimental setup and rebuild it with new, less noisy components. One could also
switch to a different experimental system, such as the
low-noise passive cavity reported in [14]. In any case,
it will probably be difficult to increase the SNR above
60 dB. There may also be algorithmic solutions, such as
using conceptors [44, 45].
The high level of experimental noise quickly pushed
the Reservoir Computer, emulating a chaotic system,
away from a given trajectory, which initiated the search
for quantitative methods for the evaluation of the experiment performance. In this work, we introduced a
few simple techniques, based on standard signal analysis methods, such as statistics of the prediction length
and visual comparison of the time series and their frequency spectra after a long autonomous run. We have
also proposed case-specific methods, such as the randomness test, that could only be applied to the Lorenz time
series. Overall, these are only the first steps towards the
answer to a very general question: given a noisy emulator of a known chaotic system, how best to evaluate
its performance? It will be interesting to understand the
relationship between the performance obtained on the estimators above, and the properties of the chaotic system,
such as the Lyapunov exponents or the dimension and geometry of the chaotic attractor. These questions should
lead to a rich new direction of enquiry in the theory of
nonlinear dynamics and complex systems.

[1] Herbert Jaeger and Harald Haas, “Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication,” Science 304, 78–80 (2004).
[2] Wolfgang Maass, Thomas Natschläger,
and Henry
Markram, “Real-time computing without stable states:
A new framework for neural computation based on perturbations,” Neural comput. 14, 2531–2560 (2002).
[3] Mantas Lukoševičius and Herbert Jaeger, “Reservoir
computing approaches to recurrent neural network training,” Comp. Sci. Rev. 3, 127–149 (2009).
[4] Fabian Triefenbach, Azarakhsh Jalalvand, Benjamin
Schrauwen, and Jean-Pierre Martens, “Phoneme recognition with large hierarchical reservoirs,” Adv. Neural
Inf. Process. Syst. 23, 2307–2315 (2010).
[5] “The 2006/07 forecasting competition for neural
networks
&
computational
intelligence,”
http://www.neural-forecasting-competition.com/NN3/
(2006).
[6] Henri Arsenault, Optical processing and computing (Elsevier, 2012).
[7] Lennert Appeltant, Miguel Cornelles Soriano, Guy
Van der Sande, Jan Danckaert, Serge Massar, Joni
Dambre, Benjamin Schrauwen, Claudio R Mirasso, and
Ingo Fischer, “Information processing using a single dynamical node as complex system,” Nat. Commun. 2, 468
(2011).

Adding output feedback to experimental reservoir
computers allows them to solve considerably more complex tasks than without output feedback. Future work
could address nonlinear computations that depend on
past information and that require persistent memory
[46], FORCE training [47] (which however requires that
the learning time scale be short compared to the reservoir time scale) and applications such as frequency modulation [23], or implementation of conceptors [44, 45].
Ideally, an entirely analogue feedback should be implemented, like in e.g. [48], rather than the digital feedback
demonstrated here. The present work is therefore just a
first step towards realising these additional applications.
Finally, going back to the question of biological implementation, our work shows that the biologically plausible
structure of reservoir computing [2, 26, 27] can be trained
to generate highly complex temporal patterns, both periodic and chaotic, even in the presence of moderate levels
of noise. Whether nature in fact implements this mechanism remains to be seen, and will depend amongst other
aspects on the amount of noise present in biological implementations of reservoir computing, and whether there
exist biologically plausible training mechanisms for this
kind of signal generation.
We thank Michiel Hermans for insightful discussions.
This work was supported by the Interuniversity Attraction Poles Program (Belgian Science Policy) Project Photonics@be IAP P7-35, by the Fonds de la Recherche Scientifique (FRS-FNRS) and by the Action de Recherche
Concertée of the Fédération Universitaire WallonieBruxelles through Grant No. AUWB-2012-12/17-ULB9.

[8] Yvan Paquot, Francois Duport, Anteo Smerieri, Joni
Dambre, Benjamin Schrauwen, Marc Haelterman, and
Serge Massar, “Optoelectronic reservoir computing,” Sci.
Rep. 2, 287 (2012).
[9] Laurent Larger, MC Soriano, Daniel Brunner, L Appeltant, Jose M Gutiérrez, Luis Pesquera, Claudio R
Mirasso, and Ingo Fischer, “Photonic information processing beyond Turing: an optoelectronic implementation of reservoir computing,” Opt. Express 20, 3241–3249
(2012).
[10] Romain Martinenghi, Sergei Rybalko, Maxime Jacquot,
Yanne Kouomou Chembo, and Laurent Larger, “Photonic nonlinear transient computing with multiple-delay
wavelength dynamics,” Phys. Rev. Let. 108, 244101
(2012).
[11] Laurent Larger, Antonio Baylón-Fuentes, Romain Martinenghi, Vladimir S. Udaltsov, Yanne K. Chembo, and
Maxime Jacquot, “High-speed photonic reservoir computing using a time-delay-based architecture: Million
words per second classification,” Phys. Rev. X 7, 011015
(2017).
[12] François Duport, Bendix Schneider, Anteo Smerieri,
Marc Haelterman, and Serge Massar, “All-optical reservoir computing,” Opt. Express 20, 22783–22795 (2012).
[13] D Brunner, M. C. Soriano, C. R. Mirasso, and I. Fischer,
“Parallel photonic information processing at gigabyte per

16

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]
[23]

[24]

[25]
[26]

[27]

[28]

second data rates using transient states,” Nat. Commun.
4, 1364 (2012).
Quentin Vinckier, François Duport, Anteo Smerieri,
Kristof Vandoorne, Peter Bienstman, Marc Haelterman,
and Serge Massar, “High-performance photonic reservoir
computer based on a coherently driven passive cavity,”
Optica 2, 438–446 (2015).
Kristof
Vandoorne,
Pauline Mechet,
Thomas
Van Vaerenbergh, Martin Fiers, Geert Morthier,
David Verstraeten, Benjamin Schrauwen, Joni Dambre,
and Peter Bienstman, “Experimental demonstration of
reservoir computing on a silicon photonics chip,” Nat.
Commun. 5, 3541 (2014).
G. Peter Zhang, “Neural networks for time-series forecasting,” in Handbook of Natural Computing, edited by
Grzegorz Rozenberg, Thomas Bäck, and Joost N. Kok
(Springer Berlin Heidelberg, Berlin, Heidelberg, 2012)
pp. 461–477.
F. wyffels and B. Schrauwen, “A comparative study of
reservoir computing strategies for monthly time series
prediction,” Neurocomputing 73, 1958 – 1964 (2010).
Piotr Antonik, Michiel Hermans, François Duport, Marc
Haelterman, and Serge Massar, “Towards pattern generation and chaotic series prediction with photonic reservoir computers,” in SPIE’s 2016 Laser Technology and
Industrial Laser Conference, Vol. 9732 (2016) p. 97320B.
Meiling Xu, Min Han, and Shunshoku Kanae, “L1/2
norm regularized echo state network for chaotic time series prediction,” in APNNS’s 23th International Conference on Neural Information Processing, LNCS, Vol. 9886
(2016) pp. 12–19.
Francis wyffels, Benjamin Schrauwen,
and Dirk
Stroobandt, “Stable output feedback in reservoir computing using ridge regression,” in International Conference on Artificial Neural Networks (Springer, 2008) pp.
808–817.
K. Caluwaerts, M. D’Haene, D. Verstraeten,
and
B. Schrauwen, “Locomotion without a brain: Physical
reservoir computing in tensegrity structures,” Artificial
Life 19, 35 – 66 (2013).
Ren Felix Reinhart and Jochen Jakob Steil, “Regularization and stability in reservoir networks with output
feedback,” Neurocomputing 90, 96 – 105 (2012).
Francis wyffels, Jiwen Li, Tim Waegeman, Benjamin
Schrauwen, and Herbert Jaeger, “Frequency modulation
of large oscillatory neural networks,” Biological Cybernetics 108, 145–157 (2014).
Piotr Antonik, Michiel Hermans, Marc Haelterman, and
Serge Massar, “Towards adjustable signal generation
with photonic reservoir computers,” in 25th International Conference on Artificial Neural Networks, Vol.
9886 (2016).
H. Jaeger, “Echo state network,” Scholarpedia 2, 2330
(2007).
Tadashi Yamazaki and Shigeru Tanaka, “The cerebellum
as a liquid state machine,” Neural Networks 20, 290 –
297 (2007).
Christian Rössert, Paul Dean, and John Porrill, “At the
edge of chaos: How cerebellar granular layer network dynamics can provide the basis for temporal filters,” PLOS
Computational Biology 11, 1–28 (2015).
Anteo Smerieri, François Duport, Yvan Paquot, Benjamin Schrauwen, Marc Haelterman, and Serge Massar,
“Analog readout for optical reservoir computers,” (2012)

pp. 944–952.
[29] François Duport, Anteo Smerieri, Akram Akrout, Marc
Haelterman, and Serge Massar, “Fully analogue photonic
reservoir computer,” Sci. Rep. 6, 22381 (2016).
[30] Quentin Vinckier, Arno Bouwens, Marc Haelterman, and
Serge Massar, “Autonomous all-photonic processor based
on reservoir computing paradigm,” (Optical Society of
America, 2016) p. SF1F.1.
[31] Piotr Antonik, Franois Duport, Michiel Hermans, Anteo
Smerieri, Marc Haelterman, and Serge Massar, “Online
training of an opto-electronic reservoir computer applied
to real-time channel equalization,” IEEE Transactions on
Neural Networks and Learning Systems PP, 1–13 (2016).
[32] Michael C Mackey and Leon Glass, “Oscillation and
chaos in physiological control systems,” Science 197,
287–289 (1977).
[33] Edward N Lorenz, “Deterministic nonperiodic flow,”
Journal of the atmospheric sciences 20, 130–141 (1963).
[34] Ali Rodan and Peter Tino, “Minimum complexity echo
state network,” IEEE Trans. Neural Netw. 22, 131–144
(2011).
[35] Andrei Nikolaevich Tikhonov,
AV Goncharsky,
VV Stepanov,
and Anatoly G Yagola, Numerical
methods for the solution of ill-posed problems, Vol. 328
(Springer Netherlands, 1995).
[36] J Doyne Farmer, “Chaotic attractors of an infinitedimensional dynamical system,” Physica D: Nonlinear
Phenomena 4, 366–393 (1982).
[37] Morris W Hirsch, Stephen Smale, and Robert L Devaney,
Differential equations, dynamical systems, and an introduction to chaos (Academic press, Boston, MA, 2003).
[38] “IEEE Standard VHDL Language Reference Manual.”
ANSI/IEEE Std 1076-1993 (1994), 10.1109/IEEESTD.1994.121433.
[39] V.A. Pedroni, Circuit Design with VHDL (MIT Press,
2004).
[40] Paul Horowitz and Winfield Hill, The art of electronics
(Cambridge University Press, 1980).
[41] John
Walker,
“ENT
program,”
http://www.fourmilab.ch/random/.
[42] George Marsaglia, “The Marsaglia Random Number
CDROM including the Diehard Battery of Tests of Randomness,” http://stat.fsu.edu/pub/diehard/.
[43] Andrew Rukhin, Juan Soto, James Nechvatal, Miles
Smid, and Elaine Barker, A statistical test suite for
random and pseudorandom number generators for cryptographic applications, Tech. Rep. (National Institute of
Standards and Technology, 2001).
[44] Herbert Jaeger, “Conceptors: an easy introduction,”
CoRR abs/1406.2671 (2014).
[45] Herbert Jaeger, “Controlling recurrent neural networks
by conceptors,” CoRR abs/1403.3369 (2014).
[46] Andr David Kovac, Maximilian Koall, Gordon Pipa, and
Hazem Toutounji, “Persistent memory in single node
delay-coupled reservoir computing,” PLOS ONE 11, 1–
15 (2016).
[47] David Sussillo and L.F. Abbott, “Generating coherent
patterns of activity from chaotic neural networks,” Neuron 63, 544 – 557 (2009).
[48] Piotr Antonik, Marc Haelterman, and Serge Massar,
“Online training for high-performance analogue readout
layers in photonic reservoir computers,” Cognitive Computation , 1–10 (2017).

