# Reducing hyperparameter dependence by external timescale tailoring

Lina Jaurigue and Kathy Lüdge  Technische Universität Ilmenau
<EMAIL_ADDRESS>
###### Abstract

Task specific hyperparameter tuning in reservoir computing is an open issue,
and is of particular relevance for hardware implemented reservoirs. We
investigate the influence of directly including externally controllable task
specific timescales on the performance and hyperparameter sensitivity of
reservoir computing approaches. We show that the need for hyperparameter
optimisation can be reduced if timescales of the reservoir are tailored to the
specific task. Our results are mainly relevant for temporal tasks requiring
memory of past inputs, for example chaotic timeseries prediciton. We consider
various methods of including task specific timescales in the reservoir
computing approach and demonstrate the universality of our message by looking
at both time-multiplexed and spatially multiplexed reservoir computing.

##  1 Introduction

Reservoir computing is a machine learning method that has received a lot of
research attention recently. There are several reasons for this. Firstly, the
learning method is very simple; only the output layer is trained via linear
regression [1, 2]. This makes the implementation of the reservoir computing
relatively straightforward and therefore accessible for researcher form a wide
range of backgrounds. Secondly, reservoir computing has been shown to perform
well on predicting complex dynamics [3, 4, 5, 6]. And, thirdly, due to the
simplicity of the training method it is well suited for hardware
implementation [7, 8]. In particular, time-multiplexed reservoir computing of
delay-based reservoirs has been a popular method of hardware implementation.
The idea to utilize a feedback delay loop instead of a network was first
published in [9] and has since then been implemented and tested for opto-
electronic [9, 10, 11, 12], optical [13, 14, 15, 16, 17] and electromechanical
[18] realizations. We also refer the reader to recent reviews on this topic
[19, 20, 21, 8, 3].

An issue that remains is the lack of an efficient method for optimising the
parameters (hyperparamters) of the reservoir in a task specific manner. This
is particularly relevant for hardware implemented reservoir computing, where
parameters can be difficult to tune. There have been many studies that have
looked into various approaches including grid searches, Bayesian optimisation,
genetic algorithms, or new network measures that correlate well with the
reservoir performance [22, 23, 24, 25, 26]. However, these approaches can
greatly increase the computational cost of the training phase, or are not
suitable for hardware implementation. In this paper we aim to contribute to
solving the issue of hyperparameter optimisation, with a particular focus on
hardware implementation. Specifically, we investigate the influence of
augmenting the memory of the reservoir at either the input or the output level
such that task relevant timescales are present in the system.

For the case of delay-based reservoirs it is known that the feedback delay-
time can be an important tuning parameter, as it directly influences the
memory of the reservoir [27, 28]. In the delay-based reservoir computing
community, the delay-time is mostly chosen to be either resonant with the
time-scale of the input clock-cycle or offset from it by one input mask
interval. However, depending on the task, other values of the delay-time are
beneficial [29, 30]. We therefore also address the influence of tuning the
memory of the reservoir directly via the internal delay as one possible
optimization step.

The fact that task relevant timescales play an important role in timeseries
prediction tasks is well established and in the context of reservoir computing
is also closely related to Takens embedding [31, 32]. In this paper we
highlight how we can externally add task relevant timescales by adding
additional delay lines at the input or output layer. These approaches are
universal for all reservoir computing realizations, ranging from electrical,
optical to digital implementations, and should not significantly decrease the
speed or the energy efficiency in hardware implementation. We show that by
externally adding the required memory to the reservoir, the performance can be
improved over the entire parameter space, thereby reducing the need for
hyperparameter optimisation.

The authors of [33] have performed a comprehensive investigation of of closed-
loop timeseries prediction using ESNs. They consider the influence of the
discretisation and look at the difference between providing full and partial
information while the authors of [34] consider the relationship between
discretisiation and dynamic timescales. We also discuss the influence of the
discretisation here as one optimization possibility, but with a focus on the
resulting memory properties. The latter depend strongly on the information fed
into the reservoir, as more memory is needed if not the full information about
the system state is provided.

There have been several other recent studies into various methods of
augmenting the reservoir computing approach. In [35] and [36] the authors
concatenation the state matrix with randomly selected delayed system states.
In [37] trainable state delays are investigated. In [38] and [39]
concatenation of the state matrix with past states is study in the context of
reservoir size reduction. In the context of hardware implemented delay-based
reservoir computing, recent efforts toward performance enhancement include
using ensembles of delay-based reservoirs [40], sequentially coupled delay-
based reservoirs [41, 42] or delay-coupled networks [29]. Our contribution is
to highlight the importance of task specific timescales, and to show how these
can easily and tunably be incorporated in hardware implemented reservoir
computing.

The paper is structured as follows. The concepts of reservoir computing (RC),
external memory augmentation and the used RC models and tasks are introduced
in Sec. 2, 3, and 4, respectively. The results for improving our delay-based
RC with adding additional memory are presented in Sec. 5 for the Mackey-Glass
and the Lorenz 63 tasks before we conclude in Sec.6.

##  2 Reservoir Computing

Reservoir computing is a scheme wherein the input response of a dynamical
system (the reservoir) is sampled a number of times and these sampled system
responses are then linearly combined to construct an output. During the
training phase, the weights of the linear combinations of the system responses
(output weights) are trained via linear regression. This is done by feeding
the reservoir a sequence of $K_{T}$ inputs and collecting $S$ responses of the
reservoir to each input into a state matrix ${\bf
S}\in\mathbb{R}^{K_{T}}\times\mathbb{R}^{S+1}$, adding a bias term of one to
each row, and then finding the vector of weights ${\bf w}\in\mathbb{R}^{S+1}$
that minimises the difference between the target output ${\bf
y}\in\mathbb{R}^{K_{T}}$ and the output ${\bf\hat{y}}={\bf S}{\bf
w}\in\mathbb{R}^{K_{T}}$. The solution to this problem is given by

${\bf w}=({\bf S}^{\textrm{T}}{\bf S}+\lambda{\bf I})^{-1}{\bf
S}^{\textrm{T}}{\bf y},$ (1)

using the Moore-Penrose pseudoinverse. Equation (1) includes Tikhonov
regularisation, with the regularisation parameter $\lambda$ and identity
matrix ${\bf I}\in\mathbb{R}^{S+1}\times\mathbb{R}^{S+1}$.

###  2.1 Spatially-multiplexed Reservoir Computing

The responses of the reservoir can either be multiplexed in space or in time
(or combinations whereof). Echo-state networks are an example of spatially
multiplexed reservoirs [1]. In this scenario, there are a number of spatially
separated dynamical nodes which evolve discretely in time. These nodes are fed
with a weighted input signal and the responses of some (or all) of the nodes
are collected into the state matrix. The evolution of such a reservoir can be
described by an equation of the form

${\bf x}\left(k+1\right)=f\left({\bf W}_{int}{\bf x}\left(k\right)+{\bf
w}_{in}I\left(k+1\right)\right),$ (2)

where ${\bf x}\left(k\right)\in\mathbb{R}^{M}$ is the vector describing the
state of the $M$ nodes at time $k$, ${\bf
W}_{int}\in\mathbb{R}^{M}\times\mathbb{R}^{M}$ is the matrix of internal
reservoir node couplings (adjacency matrix), $I\left(k\right)$ is the input
signal at time $k$, and ${\bf w}_{in}\in\mathbb{R}^{M}$ is the vector of input
weights. The function $f\left(\cdot\right)$ is a nonlinear function which is
applied element wise and commonly referred to as the activation function. In
Eq. (2) the input is one dimensional, i.e. there is a single temporal input
signal. To generalise to $p$ input signals, $I\left(k\right)$ becomes a vector
${\bf i}_{in}\left(k\right)\in\mathbb{R}^{p}$ and the vector of input weigths
${\bf w}_{in}$ becomes a matrix ${\bf
W}_{in}\in\mathbb{R}^{M}\times\mathbb{R}^{p}$.

The $k^{\textrm{th}}$ row of the state matrix is filled by sampling $S$ of the
reservoir states ${\bf x}\left(k\right)$, with $S\leq M$. For $S=M$ the
elements of the state matrix are $s_{k,j}=x_{j}\left(k\right)$ with
$x_{j}\left(k\right)$ ($j\in\left[0,..S\right)$) being the elements of the
reservoir state vector ${\bf x}\left(k\right)$ .

Spatially continuous reservoirs are also possible. In such a case the
reservoir evolution would be described by a suitable system of partial
differential equations, and the state matrix would be filled by sampling
discretely in space (and in time).

###  2.2 Time-multiplexed Reservoir Computing

In the case of time-multiplexed reservoir computing, typically, a masked input
signal is sequentially fed into the reservoir and the $S$ reservoir responses
needed to fill the state matrix are sampled sequentially in time. The input
mask is analogous to the input weights in the spatially multiplexed case, but
rather than simultaneously feeding the input signal into spatially separated
nodes with varying weights, the input signal is fed into the same node
multiple times with varying weights (mask values). If $I\left(k\right)$ is the
task-dependent input signal of length $K_{T}$, then the masked reservoir input
is given by

$J\left(i\right)=J\left(k,j\right)=I\left(k\right)m_{1}\left(j\right),$ (3)

where $i=Sk+j$, $k\in\left[0,..K_{T}\right)$, $j\in\left[0,...,S\right)$ and
$m_{1}\left(j\right)$ are the mask values. The state matrix entry
corresponding to the input $J\left(i\right)$ is given by
$s_{k,j}=x\left(i\right)=x\left(k,j\right)$, where $x\left(i\right)$ is the
sampled reservoir response. If the reservoir is a time-continuous system, then
the discrete input sequence can be transformed into a piecewise constant
function $\tilde{J}\left(t\right)$ by holding each input constant for a fixed
time $\theta$: In such a case the time for one input mask sequence, $\theta
S$, is typically referred to as the clockcycle $T$ [27]. In this study we
consider a time discrete reservoir, see Sec. 4.1.1, in which case the
dimensionless clockcycle $\tilde{T}$ is equal to the number of input mask
values $S$.

Here we have taken the number of mask values to be equal to the number of
sampled reservoir responses $S$. This is not a necessary condition, but is the
typically choice. We adhere to this convention for simplicity. We also note
that, in delay-based reservoir computing the number of sampled reservoir
responses $S$ is often referred to as the number of virtual nodes.

##  3 Reservoir Computing Memory Augmentation Methods

The memory requirements of timeseries prediction tasks are very task specific
and can vary greatly, even for tasks using the same input signal but trained
on different target signals, e.g. for different prediction horizons. To
perform well, the reservoir must have memory on the relevant timescales. This
means that there must be information of past values of the input signal in the
$k^{\textrm{th}}$ row of the state matrix in order to accurately predict the
$k^{\textrm{th}}$ target $y\left(k\right)$. Note that here we are using the
term ”memory” to refer to any influences of past inputs on the state of the
reservoir, i.e. also nonlinear transforms of past inputs. In the typical
reservoir computing scheme, this memory comes from the dynamics of the
reservoir. In this study we compare methods of augmenting the memory of a
reservoir. There are many ways in which this can be done and the choice of
method comes down to the objective one aims to achieve, which means there is
always a tradeoff between computational cost and absolute performance. The
motivation of this study is to reduce the need for hyperparameter optimisation
at a low computational cost, as well as allowing for easy hardware
implementation. We investigate two simple schemes of augmenting the memory
externally at the levels of the input and the output layer of the reservoir
(described below) and analyse how variations of the internal delay and the
discretisation change the memory.

Figure 1: Sketch of (a) a possible implementation for a delay-based reservoir
with additional delayed input and state matrix concatenation (delayed output),
and (b) the input data processed in those external delay schemes to predict a
target (orange crosses and bended arrows indicate the situation for
$n_{S}=6,n_{I}=2,\Delta s=1$).

###  3.1 Delayed-input

We consider the influence of one additional time-delayed input implemented as
in [43, 44]. The input function Eq. (3) is modified to

$\tilde{J}\left(i\right)=I\left(k\right)m_{1}\left(j\right)+I\left(k-n_{I}\right)m_{2}\left(j\right),$
(4)

where $n_{I}$ is the input-delay and $m_{2}\left(j\right)$ is the mask applied
to the delayed version of the input sequence, $I\left(k-n_{I}\right)$. This
delayed input scheme is illustrated in Fig. 1b.

In terms of a photonic setup, such a delayed input could be implemented by
adding a delay-loop before the reservoir, as depicted in Fig. 1a. Since
$n_{I}$ describes the input shift via the number of k-steps, the actual length
of the delay line depends on the sampling interval.

###  3.2 Delayed State Matrix Concatenation (Output delay)

On the level of the output layer, a method of including a task-dependent
delay, is to concatenate the current and delayed versions of the reservoir
responses, i.e. to concatenate the state matrix with a delayed version of
itself:

$\tilde{{\bf
S}}=\begin{pmatrix}s_{1,1}&\cdots&s_{1,S}&s_{1-n_{S},1}&\cdots&s_{1-n_{S},S}&1\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots&1\\\
s_{K_{T},1}&\cdots&s_{K_{T},S}&s_{K_{T}-n_{S},1}&\cdots&s_{K_{T}-n_{S},S}&1\end{pmatrix},$
(5)

where $n_{S}$ is the output delay that described how many steps (index $k$) we
go back to collect the input for the state matrix, see Fig. 1b. This approach
has recently been investigated in [39] and [38] in the context of reducing the
size of a recurrent neural network reservoir.

In order to make fair comparisons with other methods it is important to
consider whether one should compare results using the same reservoir sampling
dimension $S$ (number of virtual nodes), or whether this should be halved such
that the final state matrix dimensions are the same. The latter case has the
advantage that the dynamical system can be smaller, or that the input
clockcycle can be reduced in the case of time multiplexed reservoir computing.
In this study we choose to keep the reservoir sampling dimension $S$ the same
because in the context of hardware implemented reservoirs we believe comparing
the same sized physical reservoirs is fairer.

In a photonically implemented reservoir, this approach could be incorporated
by including an additional delay-loop after the reservoir, as sketched in Fig.
1a. Again, the physical length of the delay line is given by the product of
$n_{S}$ and the sampling time.

###  3.3 Input Discretisation

Another approach to improve the reservoir performance by incorporating task
specific memory requirements is to modify the task itself. When constructing a
time series model from data, the discretisation $\Delta t$ (sampling step) of
the input timeseries can have a significant influence on the model and on the
resulting performance [45, 34]. In the context of reservoir computing, the
input discretisation influences the memory required of the reservoir and is
related to Takens (delay) embedding [31]. Recent studies have shown that under
certain conditions reservoirs perform a delay embedding of the input
timeseries [32, 46]. For this to occur, and to yield accurate results, the
reservoir must have memory on the appropriate time scales, which are related
to the dynamics of the task [47]. This means that if the discretisation of the
input timeseries is finer than is necessary, the reservoir will need to be
able to remember input steps which lie further in the past in order for it to
perform the delay embedding. The impact of different $\Delta t$ is discussed
in Sec. 5.2.

##  4 Tasks and Methods

###  4.1 Reservoir models

We investigate the influence of the memory augmentation methods focusing
mainly on a time-multiplexed delay-based reservoir. The methods are, however,
not specific to such reservoirs. We therefore also show an example of spatial
multiplexed reservoir (i.e. a random recurrent network reservoir). Both
reservoir are described below.

####  4.1.1 Time-delayed Reservoir Model

The time-delayed reservoir (TDR) model we consider is the following iterative
map:

$x\left(i\right)=\mathcal{G}\left(Kx\left(i-N\right)+J(i)\right),$ (6)

where $\mathcal{G}\left(\cdot\right)$ is the nonlinear activation function,
$J\left(i\right)$ is the input sequence given by Eq. (3), $K$ is the strength
with which the current iteration is coupled with a past state of the system
and $N$ is the feedback delay (see Fig. 1a). Remember that the index $i$ runs
over all mask steps and input points ($i=Sk+j$).

Motivated by photonic reservoir implementations, we use a nonlinear function
which models a semiconductor optical amplifier (SOA) for the activation
function. An SOA can be described by a dynamical differential equation model
which describes the light-matter interaction and the complex carrier dynamics
[48, 49]. In our case, only the static characteristic is important, which is
why we use the simple input-response of an SOA that is given by [50]:

$\mathcal{G}\left(x\right)=\frac{g_{0}x}{1+x/a}.$ (7)

Using this activation function, Eq. (6) describes the input-response of an SOA
with weak time-delayed feedback in the limit where the time between mask steps
is much longer than the characteristic time-scale of the carrier dynamics in
the SOA, i.e. the limit that the reservoir reaches a steady state before it is
sampled at the end of each mask step. That means, we do not include any
dynamic response of the SOA and all the memory of the setup comes from the
delay-loop. A sketch of such a setup is depicted in Fig. 1a.

For the input masks in Eqs. (3) and (4) we use a set of $S$ randomly drawn
values from a uniform distribution between zero and one, $U\left(0,1\right)$,
and multiply these by a scaling factor $G_{\alpha}$, i.e.
$m_{\alpha}\left(j\right)=G_{\alpha}U^{\alpha}_{j}$ for $\alpha=1,2$ and
$U^{\alpha}_{j}$ sampled from $U\left(0,1\right)$, where $G_{1}$ is the input
scaling of the normal input and $G_{2}$ the scaling factor for the delayed
input in Eq.(4).

Unless stated otherwise, the parameters for the time delay based reservoir are
given in Tab. 1.

Table 1: Dimensionless time-delay reservoir and input parameters. Time-delayed
Reservoir Parameters
---
Parameter | Value | Parameter | Value
SOA gain parameter $g_{0}$ | 40 | node number $S$ | 250
SOA nonl. parameter $a$ | 1 | clock cycle $\tilde{T}$ | 250
regression parameter $\lambda$ | 5e-6 | input gain $G_{1}$ | 0.2
number training steps | 10000 | internal delay $N$ | 251
number of testing steps $K_{T}$ | 5000 | feedback strength $K$ | 0.03
number of realizations $Z_{R}$ | 10 |  |

####  4.1.2 Random Network Reservoir Model

We use a random network as described by Eq. (2) with

$f\left(x\right)=\Theta\left(x-1\right),$ (8)

where $\Theta\left(\cdot\right)$ is the Heaviside step function. The internal
node coupling weights are given by ${\bf W}_{int}=K_{r}{\bf W}_{int}^{u}$
where the entries of ${\bf W}_{int}^{u}$ are chosen randomly from a uniform
distribution between zero and one, and $K_{r}$ is a scaling factor that
determines and coupling strength and the spectral radius $\rho$. Analogous to
the time delay based reservoir, the input weights are also sampled from
$U\left(0,1\right)$ and scaled by a factor $G_{1}$.

The motivation of our choice of random network is not to achieve the best
performance possible, but rather to demonstrate the universality of our
results on a reservoir with very different properties to the time delay based
reservoir described above.

###  4.2 Error Measure

We quantify the reservoir computing performance using the normalised root mean
squared error (NRMSE), defined as

$\textrm{NRMSE}=\sqrt{\frac{\sum_{k=1}^{K_{T}}\left(y_{k}-\hat{y}_{k}\right)^{2}}{K_{T}\textrm{var}\left({\bf
y}\right)}},$ (9)

where $y_{k}$ are the target values, $\hat{y}_{k}$ are the outputs produced by
the reservoir computer, the number of testing steps $K_{T}$ is the length of
the vector ${\bf y}$ and $\textrm{var}\left({\bf y}\right)$ is the variance of
the target sequence.

###  4.3 Mackey-Glass Time-series Prediction

Figure 2: Mackey-Glass task: (a) Timeseries as given by Eq. (10) showing
chaotic dynamics for two different discretisation steps $\Delta t=12$ (orange
dots) and $\Delta t=1$ (green dots), (b) corresponding autocorrelation
function as a function of the time lag.

Multi-step-ahead prediction of timeseries produced by the Mackey-Glass
equation is a common benchmarking task [51, 52, 42]. The Mackey-Glass delay-
differential equation is given by [53]

$\frac{dx}{dt}=\beta\frac{x\left(t-\tau\right)}{1+x\left(t-\tau\right)^{n}}-\gamma
x.$ (10)

Using the standard parameters, $\tau=17$, $n=10$, $\beta=0.2$ and
$\gamma=0.1$, this system exhibits chaotic dynamics.

In this study we focus mainly on $\Delta s=12$ step-ahead prediction of this
system with an input discretisation of $\Delta t=1$ (open-loop). To create the
input sequence $I\left(k\right)$, we generate a timeseries by numerically
integrating Eq.(10) using a Runge-Kutta fourth order method with Hermitian
interpolation for the mid-steps of the delayed terms and with a time step of
$h=10^{-2}$. This timeseries is then sampled with a time step of $\Delta t=1$.
The corresponding target sequence is given by $I\left(k+\Delta s\right)$. We
have chosen $\Delta t=1$ as our default sampling, as this is a commonly used
discretisation and it is small enough that finer features of the time series
can be resolved, as can be seen by the green dots in Fig.2a. It is however not
the optimal discretisation for delay-embedding [45].

The chaotic Mackey-Glass timeseries is shown in Fig. 2a, along with its
corresponding autocorrelation function in Fig. 2b. The dynamics is chaotic,
but exhibits clear oscillations on a scale of approximately 50 time units. For
the reservoir to reconstruct the dynamics of this timeseries, it would need to
have memory on the order of 0.1-0.5 of these features [47]. Furthermore, the
autocorrelation shows that points separated by one time unit are highly
correlated. This indicates that, with an input discretisation of $\Delta t=1$,
consecutive inputs will contain redundant information and the reservoir will
need to be able to retain information from about 10 steps into the past for a
good delay embedding. We therefore laster also consider other combinations of
$\Delta t$ and $\Delta s$ to compare the influence of the discretisation of
the input signal on the memory requirements of the reservoir and the
performance.

###  4.4 Lorenz 63 Time-series Prediction

The Lorenz system [54] is given by

$\displaystyle\frac{dX}{dt}=c_{1}Y-c_{1}X,\quad\frac{dY}{dt}=X(c_{2}-Z)-Y,\quad\textrm{and}\quad\frac{dZ}{dt}=XY-
c_{3}Z.$ (11)

With $c_{1}=10$, $c_{2}=28$ and $c_{3}=8/3$ this system exhibits chaotic
dynamics as can be seen for the X and Z coordinate in Fig. 3a and b. We use
the $X$ variable, sampled with a step size of $\Delta t=0.1$ (orange dots), as
the input $I\left(k\right)$ for two time series prediction tasks. The first
task is one step-ahead ($\Delta s=1$) prediction of the $X$ variable. The
second task is cross prediction of the $Z$ variable ($\Delta s=0$). For the
cross prediction task we also consider a discretisation of $\Delta t=0.02$
(gren dots). We generate the Lorenz timeseries using Runge-Kutta fourth order
numerical integration with a time step of $h=10^{-3}$.

Figure 3: Lorenz 63 task: Timeseries for X (a) and Z (b) as given by Eq. (11)
(chaotic dynamics) for two different discretisation steps $\Delta t=0.1$
(orange dots) and $\Delta t=0.02$ (green dots). (c) The corresponding
normalized autocorrelation $\langle X(t)X(t-\tau)\rangle$ and cross
correlation function $\langle X(t)Z(t-\tau)\rangle$.

When doing timeseries prediction using the an ODE system like the Lorenz 63
system, there is a qualitative difference between the tasks where all
dynamical variables are provided as input and tasks where this is not the
case. Using the Lorenz 63 system, a comparitive study of these two cases was
made in [33]. When all variables are provided, no memory of past steps is
needed to determine the evolution of the dynamics. However, if only partial
information is provided, then the reservoir must perform a delay embedding
[32]. For our study, only the latter case is relevant as we only provide one
of the three dynamical variables for the Lorenz63 task and not the complete
history for the Mackey-Glass task. The sampling step for an optimal embedding
for the Lorenz 63 system is about 0.1 [45, 34]. Figure 3c shows the
autocorrelation of the $X$ variable (green line) and the cross-correlation
between the $X$ and $Z$ variables (orange line) for the parameters given above
that lead to the chaotic dynamics (Fig. 3a,b). We will see later on in Sec.
5.3 that the differences in the correlation between the input and the target
for the $X\left(t\right)\rightarrow X\left(t+\Delta t\right)$ and
$X\left(t\right)\rightarrow Z\left(t\right)$ tasks have important consequences
for the optimal discretisation.

###  4.5 Simulation Parameters

For all tasks we rescale the input and target sequences to the range between
zero and one. Before beginning the training and testing phases for these
tasks, we initialise the reservoirs by feeding in 10000 inputs
$I\left(k\right)$. In the training phase we use 10000 inputs while in the
testing phase 5000 inputs are used. All results presented in the next section
show the testing error averaged over $Z_{R}$ realisations of the random input
masks. For the NRMSE results the associated error bars are the standard
deviation of the $Z_{R}$ realisations.

##  5 Results and Discussion

###  5.1 Mackey-Glass task performance comparison

We are interested in the influence of including task specific timescales in
the reservoir computing scheme. The delayed-input and delayed-state matrix
methods described above can be applied to all reservoir computing schemes.
However, when using a delay-based reservoir, we also have the option of
directly tuning the memory properties by changing the feedback delay time. For
the reservoir used here (Eq. (6)), this corresponds to the parameter $N$. In
most delay-based reservoir computing studies, the authors typically choose
$N=\tilde{T}$ (delay-time $\tau=T$ in the time-continuous case) or
$N=\tilde{T}+1$ ($\tau=T+\theta$) [9, 11, 13, 14, 55, 40, 56, 57], where
$\tilde{T}$ is the input clockcycle. It has been shown that resonances between
the feedback delay time and clock-cycle are generally detrimental to the
reservoir computing performance [30, 28, 27], we therefore choose the
desynchronised configuration $N=\tilde{T}+1$ ($N=S+1$) as our default and
compare its impact with our suggested memory augmentation methods. The other
parameters of the time-delay based reservoir are the feedback strength $K$,
the input scaling $G_{1}$ and the parameters of SOA nonlinearity $g_{0}$ and
$a$. The latter two are determined by properties of the SOA and are therefore
not freely tunable. Since we are motivated by hardware implementability, in
the following, we only consider the influence of $K$ and $G_{1}$, while
keeping $g_{0}$ and $a$ fixed.

Figure 4: Time-delay based reservoir \- NRMSE for the Mackey-Glass 12-step-
ahead prediction task ($\Delta t=1$, $\Delta s=12$) color coded in dependence
of the feedback strength $K$ and the input scaling $G_{1}$ for two reservoir
sizes, $S=10$ (upper row (a-d) with $Z_{R}=50$) and $S=500$ (lower row (e-h)
with $Z_{R}=10$) and 4 different delay schemes (columns): (b,f) delayed input
with $G_{2}=0.6$, $n_{I}=8$, (c,g) delayed output with $n_{S}=8$, and (d,h)
optimal internal delay $N=5S+1$. Parameters as in Tab. 1.

Figure 4a shows the NRMSE for the Mackey-Glass 12-step-ahead prediction task
as a function of the feedback strength $K$ and the input scaling $G_{1}$ for a
reservoir consisting of $S=10$ virtual nodes. The parameter ranges are chosen
such that the characteristics of the nonlinearity span from linear to highly
nonlinear and the feedback strength covers experimentally realisable values.
Over the entire parameter plane, relatively poor performance is achieved
(orange colors in Fig. 4a). Increasing the number of virtual nodes to $S=500$
(Fig. 4e), the minimum NRMSE is decreased from $0.37\pm 0.05$ to $0.064\pm
0.002$. However, this comes at a cost of a factor 50 longer clockcycle
$\tilde{T}=S$ and requires sensitive tuning of the reservoir parameters.

To look at the influence of a delayed-input, as described in Sec.3.1, we
select the parameters corresponding to the lowest error in Fig. 4a (Fig. 4e)
and perform a grid search in the input-delay $n_{I}$ and input scaling
$G_{2}$. Please see Fig. 9 in A for the corresponding scans. The $n_{I}$ and
$G_{2}$ parameters for which the lowest error is achieved
($n_{I}=8,G_{2}=0.6$) are then used for the $G_{1}$-$K$ scan shown in Fig. 4b
(Fig. 4f). The delayed-input parameters are therefore not optimised for each
pair of values in the $G_{1}$-$K$ parameter plane, however, we have checked
that the optimal $n_{I}$ does not vary significantly. With the delayed input
included, not only is the minimum error reduced, the variation of the error in
dependence of $G_{1}$ and $K$ is also greatly reduced (mostly green colors in
Fig. 4f). For the delayed-state matrix concatenation method described in
Sec.3.2, we use the same approach as above, i.e. we optimise the state matrix
delay $n_{S}$ for the optimal $G_{1}$,$K$ pair from Fig. 4a (Fig. 4e) and then
use this value $n_{S}=8$ to produce Fig. 4c (Fig. 4g). The results of
optimising the internal delay of the reservoir, again using the same approach
as above, are shown is Fig. 4d (Fig. 4h). The optimised delay values for the
Mackey-Glass 12 step ahead prediction task ($\Delta t=1$) are $n_{I}=8$,
$n_{S}=8$ and an internal delay length of about 5 times the clock cycle
$N=5S+1$.

Figure 5: Impact of reservoir size \- (a) Minimum NRMSE and (b) parameter
sensitivity $P_{0.1}$ for the Mackey-Glass 12-step-ahead prediction task
($\Delta t=1$, $\Delta s=12$) in dependence of the virtual node number $S$ for
the default delay based setup (black), with delayed input $G_{2}=0.6$ and
$n_{I}=8$ (orange), with delayed output $n_{S}=8$ (green) and for optimal
internal delay $N=5S+1$ (grey). Parameters as in Table 1.

From inspecting the color distribution in the three rightmost columns of Fig.
4, we see that all three methods of augmenting or changing the memory of the
reservoir to optimise the performance for this Mackey-Glass task lead to a
large decrease in the error over most of the depicted $G_{1}$-$K$ parameter
plane (see minimum error achieved is always depicted in the inset of each
panel). To summarise the results of Fig. 4, and similar calculations for
virtual node numbers between $S=10$ and $S=500$, we show the minimum error as
a function of the number of virtual nodes $S$ in Fig. 5a. To quantify the
sensitivity of the error to the reservoir parameters $G_{1}$ and $K$ we
introduce a quantity $P_{0.1}$, which is the percentage of the depicted
$G_{1}$-$K$ plane for which the NRMSE is under a threshold value of 0.1.
Figure. 5b shows $P_{0.1}$ as a function of the number of virtual nodes $S$
(the results for another threshold leading to $P_{0.05}$ are given in Fig.10
in B). Figure 5 shows that, regardless of the method, including a task
relevant time scale vastly improves the performance over a wide range of
reservoir parameters. The smallest error in Fig. 5 is achieved when the
internal delay $N$ of the reservoir is optimised (for $S=500$). Whereas, the
least parameter sensitivity is achieved using the delayed state matrix
concatenation method (dark green points in Fig. 5b).

It should be noted that which of the memory augmentation methods produces the
best results is strongly task dependent. This is because all methods lead to a
different nonlinear mixing of the input signal. An example where the delayed-
input method out performs the delayed state matrix concatenation method in
terms of the sensitivity measure $P_{0.1}$ is shown in Fig. 12 in B where a
Mackey Glass task one step prediction with $\Delta t=12$, $\Delta s=1$ was
used. The strongest difference compared to Fig. 5 are seen for the default
case with $N=S+1$. The larger discretisation $\Delta t=12$ is a much better
delay embedding for the Mackey Glass timeseries and thus already the non-
optimized internal delay configuration leads to good prediction results. The
next section will focus on this issue.

###  5.2 Relationship between optimal delays, task and reservoir

Figure 6: Impact of prediction horizon on output delay \- NRMSE for the
Mackey-Glass prediction task ($\Delta t=1$) color coded as a function of the
prediction step $\Delta s$ and the externally added state matrix delay
$n_{S}$. (a) extreme learning machine ($K=0$), (b) default internal delay
setup, and (c) optimal internal delay setup. Parameters: $S=50$ and remaining
parameters as in Tab. 1.

The optimal augmentation delays ($n_{I}$ and $n_{S}$) depend on the
requirements of the task and on the memory properties of the reservoir. We
demonstrate this by altering the Mackey-Glass task. Firstly, by changing
$\Delta s$ for a fixed discretisation $\Delta t$ (Fig. 6) and secondly, by
altering $\Delta s$ and $\Delta t$ under the constraint that $\Delta
s\times\Delta t$ is fixed, i.e. that the time span of the prediction step is
constant (Fig. 7) . We then consider three implementations of the reservoir
which have very different recall abilities. Results for these cases are shown
in Fig. 6 and Fig. 7. In Fig. 6a and Fig. 7a the feedback strength $K$ is set
to zero. In this limit the reservoir itself is a single layer feedforward
network, i.e. an extreme learning machine [58, 52], and has no memory on its
own while memory is only added via the delayed output. Figure 6b and Fig. 7b
show results for the default configuration, where the internal delay is one
time step longer than the input clockcycle. In this configuration the system
behaves similar to a simple ring [27]; from one clockcycle to the next
neighbouring nodes couple and the rate of information decay is determined by
$K$. The output delay $n_{S}$ adds additional memory. In Fig. 6c and Fig. 7c
the internal delay is about five times as long as the input clockcycle,
$N=5S+1$. This means that the reservoir has zero ability to recall information
between one and four steps into the past (remember that the optimal sampling
step for Mackey Glass are 12 steps for $\Delta t=1$).

For $K=0$, Fig. 6a clearly shows that as $\Delta s$ is changed, so does the
optimal input delay, ($n_{s}$ corresponding to the darkest green). The sum
between optimal output delay and predicted step $\Delta s$, however stays
constant until $\Delta s=19$ lading to a z-like behaviour of the optimal
output delay. For the default internal delay configuration (Fig. 6b), the
optimal state matrix delay is shifted slightly to lower values compared with
the $K=0$ case. This is because the state matrix delay compensates for the
memory that is missing in the reservoir and for the default case there is
memory due to the coupling between adjacent nodes. When the internal delay is
$N=5S+1$ (Fig. 6c) the reservoir is already well suited for delay embedding
the Mackey-Glass input signal, therefore adding the state matrix concatenation
has only a very small influence. Only slight improvements can be seen for
$n_{S}$ values that are multiples of 5.

When the input discretisation $\Delta t$ is changed then the required memory
can vary substantially even if the prediction horizon $\Delta s\cdot\Delta t$
stays constant. This can be seen most clearly in Fig. 7a by the varying
position of the minimum error along the $n_{S}$ axis. However, $n_{S}$ is
given in units of $\Delta t$ and the optimal value of $n_{S}$ in units of the
original Mackey-Glass time series stays constant at about
$n_{S}^{opt}\cdot\Delta t=9$ for the extreme learning limit $K=0$. For the
default delay configuration in Fig. 7b additional minima are seen. Inspecting
the minima of the curves in Fig. 7b, it can also be seen that different input
discretisations result in different lowest errors ($\Delta t=3$ is the best
choice here), which was not the case in the extreme learning limit in Fig. 7a.
Thus, the optimal discretisation depends on the memory of the reservoir.

In Fig. 7c, an internal delay configuration of $N=5S+1$ was chosen which was
optimal for the discretisation of $\Delta t=1$ (vyan curve). If this value is
changed, $N$ is not optimized and adding an output delay with delay $n_{S}$
can now provide for the missing memory. In general, a clear predictive
relationship between the optimal augmentation delays, the task requirements
and the properties of the reservoir are non-trivial. We thus recommend
treating the delays as tuning parameters, guided by Takens embedding.

Figure 7: Optimal output delay \- Impact of discretisation on the NRMSE for
the Mackey-Glass prediction task for constant prediction horizon $\Delta
s\cdot\Delta t$ as a function of the state matrix delay $n_{S}$. (a) K=0, (b)
default internal delay of $S=N+1$, and (c) internal delay of $5S+1$.
Parameters: $S=50$, remaining parameters as in Tab. 1.

###  5.3 Lorenz results

So far we have used various versions of the Mackey-Glass timeseries prediction
task. In this section we present a brief summary of results obtained using the
Lorenz 63 system. We do this to make clear that the influence of the memory
augmentation methods is not specific to the Mackey-Glass system and to show
that the benefits of such methods are related to the memory requirements of
the task. In Table 2 we shown the minimum NRMSE, $P_{0.1}$ and $P_{0.05}$,
which we have calculated using the same optimisation procedure as in Sec. 5.1.
For the $X(t)\rightarrow X\left(t+\Delta t\right)$ task the benefits of the
memory augmentation methods are marginal. The state matrix concatenation
method shows the best performance, both in terms of the absolute error and the
parameter sensitivity. However, we remind the reader that the effective output
dimension is doubled when the state matrix concatenation is used (see Sec.
3.2).

Table 2: Time delay based reservoir - Lorenz 63 timeseries prediction
performance for 3 different tasks (rows) and 3 different delay configurations
(columns). Parameters as in Tab. 1.

  * | Task: Lorenz X, one step ahead with $\Delta t$=0.1
---|---
| $X(t)\rightarrow X(t+\Delta t)$
| Default | Delayed Input | Delayed Output
|  |  $n_{I}=2$, $G_{2}=0.2$ | $n_{s}=1$
min NRMSE | $0.0193\pm 0.0007$ | $0.0176\pm 0.0007$ | ${\bf 0.0140\pm 0.0006}$
$P_{0.1}$ | $0.790\pm 0.009$ | $0.830\pm 0.006$ | ${\bf 0.882\pm 0.003}$
$P_{0.05}$ | $0.47\pm 0.01$ | $0.53\pm 0.01$ | ${\bf 0.597\pm 0.008}$
| Task: Lorenz cross prediction with $\Delta t$=0.1
| $X(t)\rightarrow Z(t)$
| Default | Delayed Input | Delayed Output
|  |  $n_{I}=2$, $G_{2}=0.2$ | $n_{s}=2$
min NRMSE | $0.0058\pm 0.0009$ | $0.0046\pm 0.0004$ | ${\bf 0.0030\pm 0.0002}$
$P_{0.1}$ | $0.41\pm 0.02$ | $0.54\pm 0.02$ | ${\bf 0.898\pm 0.008}$
$P_{0.05}$ | $0.13\pm 0.02$ | $0.24\pm 0.03$ | ${\bf 0.55\pm 0.01}$
| Task: Lorenz cross prediction with $\Delta t$=0.02
| $X(t)\rightarrow Z(t)$
| Default | Delayed Input | Delayed Output
|  |  $n_{I}=17$, $G_{2}=0.2$ | $n_{s}=20$
min NRMSE | $0.099\pm 0.006$ | $0.030\pm 0.002$ | ${\bf 0.0167\pm 0.0005}$
$P_{0.1}$ | $0.0003\pm 0.0003$ | ${\bf 0.09\pm 0.02}$ | $0.044\pm 0.006$
$P_{0.05}$ | $0\pm 0$ | $0.007\pm 0.001$ | ${\bf 0.012\pm 0.001}$

For the cross-prediction task, $X\left(t\right)\rightarrow Z\left(t\right)$,
the memory augmentation methods greatly improve the results when $\Delta
t=0.02$, but lead to only a slight improvement for $\Delta t=0.1$. However,
better performance is achieved for the $\Delta t=0.1$ case. This is because
$\Delta t=0.1$ is the optimal delay for delay embedding the Lorenz 63 system
[33, 45]. This means that with $\Delta t=0.02$ the reservoir needs to retain
information from further in the past and, hence, augmenting the memory has a
greater impact.

Although not shown here, in contrast to the cross-prediction task, the
performance for the $X(t)\rightarrow X\left(t+\Delta t\right)$ task improves
when $\Delta t$ is smaller than 0.1. This can be understood by looking at the
autocorrelation plot in Fig. 3c. For small lags the autocorrelation of
$X\left(t\right)$ is high, whereas the cross-correlation with
$Z\left(t\right)$ is close to zero even for small lags. This means that the
$X(t)\rightarrow X\left(t+\Delta t\right)$ task becomes more linear and
therefore easier as $\Delta t$ approaches zero, however, for the
$X\left(t\right)\rightarrow Z\left(t\right)$ task this is not the case.

###  5.4 Random Network

Delay based reservoirs and the more commonly used ESNs, or random recurrent
networks, have very different memory properties. In the case of the ESN the
memory is determined indirectly by the random topology of the reservoir and
will mostly be monotonically decreasing. However, for a delay-based reservoir,
such as the one used here, the memory can be directly tuned by changing the
ratio of the delay and the input clocktime, and will not decrease
monotonically if the delay is sufficiently large compared with the clocktime
and the memory of the solitary nonlinear node [29, 28, 27]. These differences,
however, have no bearing on the influence of the memory augmentation methods.
We demonstrate this in Fig. 8 which shows one example of the performance
improvement that can be achieved with the delayed state matrix concatenation
method. Compared with the time delay based reservoir (see Fig. 4), the overall
performance is worse because our chosen nonlinearity produces a binary output
(Sec. 4.1.2). We also note, that for Fig. 8 we did not optimise the delay
$n_{S}$, we simply set it to $n_{S}=8$, it is therefore possible that the
results could be improved further. However, despite this a stark improvement
is obtained over the entire parameter range when the delayed state matrix
concatenation is included.

Figure 8: Random network \- NRMSE for the Mackey-Glass 12-step-ahead
prediction task in dependence of the spectral radius $\rho$ and the input
scaling $G_{1}$. Parameters: $n_{S}=8$, $S=500$ and $Z_{R}=10$.

##  6 Conclusion

We have shown that including task specific timescales can significantly
improve the performance of a reservoir computer for timeseries prediction
tasks. While the timescale itself matters, the method of including it does not
play such an important role for the overall performance and should be chosen
according to the requirements of a specific task. Our focus was on a simple
means of improving the performance, which reduces the need for hyperparameter
tuning in the context of hardware implementation. There are pros and cons of
all methods. The state matrix concatenation method has the advantage that it
does not influence the reservoir as it comes after the reservoir sampling
step. This is likely the reason why this method showed the best performance in
terms of the hyperparameter sensitivity. It also has the advantage that it
only introduces one tuning parameter, the delay. However, if this method is to
be implemented in hardware, the previous reservoir states need to be stored
and and the state vector is twice as large. The delayed input method, on the
other hand, has the comparative disadvantage that it introduces two tuning
parameters, the delay and the input scaling. However, this method produces
very similar results to the state matrix concatenation without an increased
effective output dimension. In the case of a delay-based reservoir, tuning the
internal delay is also an option, but can be difficult in practise,
particularly in photonic systems which can be very feedback phase sensitivity
[59, 60, 61]. Lastly, one can optimise the discretisation of the input signal,
but here one must also consider the desired output discretisation.

Extensions or variations of methods discussed here are possible and some have
already been study [55, 39, 38, 62, 35, 36]. Whether or not to apply such
methods once again comes down to the objective and a trade-off between
simplicity and absolute performance.

##  Appendix A Delay optimisation for augmented memory

Figure 9: Optimizing delays \- NRMSE for the Mackey-Glass 12-step-ahead
prediction task ($\Delta t=1$) in dependence of (a) state matrix concatenation
delay $n_{S}$, (b) internal delay $N$ (normalized to the number of virtual
nodes $S$ and chosen off-resonant), and (c) input delay $n_{I}$ and input
scaling $G_{2}$. For parameters see Tab. 1.

Figure 9 shows the grid search and parameters scans performed to find the
optimal internal delay $N$, the optimal input delay $n_{I}$ and the best state
matrix concatenation delay $n_{S}$. For the internal delay scan $N$ was
increased in steps of $S$ starting at $N=1$. This was done to avoid the
resonance between the delay and clockcycle, which are known to be detrimental
[43, 30].

##  Appendix B Further Mackey-Glass Performance Results

Figure 10 shows the parameter sensitivity measure $P_{0.05}$ for the Mackey-
Glass 12 step-ahead prediction task with $\Delta t=1$. $P_{0.05}$ is the
percentage of the $G_{1}$-$K$ parameter space which is depicted in Fig. 4 for
which the NRMSE is below 0.05.

Analogous to Figs 5, Fig. 11 shows the minimum NRMSE and the sensitivity
$P_{0.1}$ for an input discretisation of $\Delta t=12$ and one step-ahead
prediction ($\Delta s=1$), i.e. the size of the prediction step is the same as
in the $\Delta t=1$ and $\Delta s=12$ case. The default reservoir now shows
much better performance. This is because less redundant information is being
feed into the reservoir and the memory requirements are reduced. The
performance is still improved when a delay-input or delayed state matrix
concatenation is included, however the impact is reduced and the optimised
delays are smaller (now $n_{I}=4$ and $n_{S}=4$, compared with $n_{I}=8$ and
$n_{S}=8$ for the $\Delta t=1$ and $\Delta s=12$ case). We also note that, in
terms of the sensitivity measure $P_{0.1}$, the delayed input method out
performs the delayed state matrix concatenation method for large node numbers,
even though the effective state matrix dimension is only half as big. This
also holds for $P_{0.05}$, see Fig. 12.

Figure 10: Time-delayed reservoir \- Sensitivity measure $P_{0.05}$ for the
Mackey-Glass 12-step-ahead prediction task $\Delta t=1$ in dependence of the
number of virtual nodes $S$, obtained from the results depicted in Fig. 4.
Parameters: $N=S+1$ (black), $n_{I}=8$ and $G_{2}=0.6$(orange), $n_{S}=8$
(green), $N=5S+1$ (grey) and as in Table 1.

Figure 11: Time-delayed reservoir \- (a) Minimum NRMSE and (b) sensitivity
measure $P_{0.1}$ for the Mackey-Glass $\Delta s=1$ step-ahead prediction task
with $\Delta t=12$ in dependence of the number of virtual nodes $S$ for the
default delay based setup (black), with delayed input $G_{2}=0.6$ and
$n_{I}=8$ (orange), and with delayed output $n_{S}=8$ (green). Parameters as
in Table 1.

Figure 12: Impact of reservoir size \- Sensitivity measure $P_{0.05}$ for the
Mackey-Glass $\Delta s=1$ step-ahead prediction task with $\Delta t=12$ in
dependence of the number of virtual nodes $S$ for the default delay based
setup (black), with delayed input $G_{2}=0.6$ and $n_{I}=8$ (orange), and with
delayed output $n_{S}=8$ (green). Parameters as in Table 1.

## References

## References

  * [1] Jaeger H 2001 The ’echo state’ approach to analysing and training recurrent neural networks GMD Report 148 GMD - German National Research Institute for Computer Science
  * [2] Maass W, Natschläger T and Markram H 2002 Neural. Comput. 14 2531 ISSN 0899-7667
  * [3] Nakajima M, Tanaka K and Hashimoto T 2021 Commun. Phys. 4 2399
  * [4] Bueno J, Brunner D, Soriano M C and Fischer I 2017 Opt. Express 25 2401
  * [5] Kuriki Y, Nakayama J, Takano K and Uchida A 2018 Opt. Express 26 5777
  * [6] Shahi S, Fenton F H and Cherry E M 2022 Mach. Learn. Appl. 8 100300 ISSN 2666-8270
  * [7] Van der Sande G, Brunner D and Soriano M C 2017 Nanophotonics 6 561
  * [8] Tanaka G, Yamane T, Héroux J B, Nakane R, Kanazawa N, Takeda S, Numata H, Nakano D and Hirose A 2019 Neural Netw. 115 100–123 ISSN 0893-6080
  * [9] Appeltant L, Soriano M C, Van der Sande G, Danckaert J, Massar S, Dambre J, Schrauwen B, Mirasso C R and Fischer I 2011 Nat. Commun. 2 468
  * [10] Larger L, Soriano M C, Brunner D, Appeltant L, Gutierrez J M, Pesquera L, Mirasso C R and Fischer I 2012 Opt. Express 20 3241–3249
  * [11] Paquot Y, Duport F, Smerieri A, Dambre J, Schrauwen B, Haelterman M and Massar S 2012 Sci. Rep. 2 287
  * [12] Chen Y, Yi L, Ke J, Yang Z, Yang Y, Huang L, Zhuge Q and Hu W 2019 Opt. Express 27 27431
  * [13] Brunner D, Soriano M C, Mirasso C R and Fischer I 2013 Nat. Commun. 4 1364
  * [14] Dejonckheere A, Duport F, Smerieri A, Fang L, Oudar J L, Haelterman M and Massar S 2014 Opt. Express 22 10868
  * [15] Vinckier Q, Duport F, Smerieri A, Vandoorne K, Bienstman P, Haelterman M and Massar S 2015 Optica 2 438
  * [16] Hou Y S, Xia G Q, Yang W Y, Wang D, Jayaprasath E, Jiang Z, Hu C X and Wu Z M 2018 Opt. Express 26 10211
  * [17] Borghi M, Biasi S and Pavesi L 2021 Sci. Rep. 11 15642 ISSN 2045-2322
  * [18] Dion G, Mejaouri S and Sylvestre J 2018 J. Appl. Phys. 124 152132
  * [19] Brunner D, Penkovsky B, Marquez B A, Jacquot M, Fischer I and Larger L 2018 J. Appl. Phys. 124 152004
  * [20] Brunner D, Soriano M C and Van der Sande G 2019 Photonic Reservoir Computing, Optical Recurrent Neural Networks (Berlin, Boston: De Gruyter) ISBN 978-3-11-058200-0
  * [21] Nakajima K 2020 Jpn. J. Appl. Phys. 59 060501
  * [22] Penkovsky B, Larger L and Brunner D 2018 J. Appl. Phys. 124 162101 ISSN 0021-8979
  * [23] Griffith A, Pomerance A and Gauthier D J 2019 Chaos 29 123108
  * [24] Racca A and Magri L 2021 Neural Netw. 142 252–268 ISSN 0893-6080
  * [25] Valencia C H, Vellasco M M B R and Figueiredo K 2023 Neurocomputing 545 126317 ISSN 0925-2312
  * [26] Carroll T L 2020 Chaos 30 083130
  * [27] Hülser T, Köster F, Jaurigue L C and Lüdge K 2022 Opt. Mater. Express 12 1214
  * [28] Köster F, Yanchuk S and Lüdge K 2021 J. Phys. Photonics 3 024011
  * [29] Röhm A, Jaurigue L C and Lüdge K 2019 IEEE J. Sel. Top. Quantum Electron. 26 7700108
  * [30] Stelzer F, Röhm A, Lüdge K and Yanchuk S 2020 Neural Netw. 124 158 ISSN 0893-6080
  * [31] Takens F 1981 Detecting strange attractors in turbulence Dynamical Systems and Turbulence, Warwick 1980 ed Rand D and Young L S (Berlin, Heidelberg: Springer Berlin Heidelberg) pp 366–381 ISBN 978-3-540-38945-3
  * [32] Hart A G, Hook J L and Dawes J H 2020 Neural Netw. 128 234–247 ISSN 0893-6080
  * [33] Storm L, Gustavsson K and Mehlig B 2022 Mach. Learn.: Sci. Technol. 3 045021
  * [34] Tsuchiyama K, Röhm A, Mihana T, Horisaki R and Naruse M 2023 Chaos 33 063145 ISSN 1054-1500 arXiv: 2302.10761 (Preprint 2302.10761) URL https://doi.org/10.1063/5.0143846
  * [35] Del Frate E, Shirin A and Sorrentino F 2021 Chaos 31 121103
  * [36] Carroll T L and Hart J D 2022 Chaos 32 083122 arXiv:2205.02267
  * [37] Holzmann G and Hauser H 2010 Neural Netw. 23 244–256 ISSN 0893-6080
  * [38] Sakemi Y, Morino K, Leleu T and Aihara K 2020 Sci. Rep. 10 21794 ISSN 2045-2322
  * [39] Marquez B A, Suarez-Vargas J and Shastri B J 2019 Phys. Rev. Res. 1 033030
  * [40] Ortín S and Pesquera L 2017 Cognitive Computation 9 327 ISSN 1866-9964
  * [41] Penkovsky B, Porte X, Jacquot M, Larger L and Brunner D 2019 Phys. Rev. Lett. 123 054101
  * [42] Goldmann M, Köster F, Lüdge K and Yanchuk S 2020 Chaos, 30 093124
  * [43] Jaurigue L C, Robertson E, Wolters J and Lüdge K 2021 Entropy 23 1560 ISSN 1099-4300
  * [44] Jaurigue L C, Robertson E, Wolters J and Lüdge K 2022 SPIE Proceedings 12204 1220408
  * [45] Kantz H and Schreiber T 2003 Nonlinear Time Series Analysis 2nd ed (Cambridge University Press)
  * [46] Duan X, Ying X, Leng S, Kurths J, Lin W and Ma H 2023 Phys. Rev. Research 5
  * [47] Fraser A M and Swinney H L 1986 Phys. Rev. A 33(2) 1134–1140
  * [48] Zajnulina M, Lingnau B and Lüdge K 2017 IEEE J. Sel. Top. Quantum Electron. 23 3000112 ISSN 1077-260X
  * [49] Lingnau B and Lüdge K 2017 Quantum-dot semiconductor optical amplifiers Handbook of Optoelectronic Device Modeling and Simulation (Series in Optics and Optoelectronics vol 1) ed Piprek J (Boca Raton, FL, USA: CRC Press) chap 23 ISBN 9781498749466
  * [50] Coldren L A, Corzine S W and Mashanovitch M L 2012 Diode Lasers and Photonic Integrated Circuits (John Wiley & Sons) ISBN 978-0-470-48412-8
  * [51] Shi Z and Han M 2007 IEEE Trans. Neural Netw. Learn. Syst. 18 359–372
  * [52] Ortín S, Soriano M C, Pesquera L, Brunner D, San-Martín D, Fischer I, Mirasso C R and Gutierrez J M 2015 Sci. Rep. 5 14945 ISSN 2045-2322
  * [53] Mackey M C and Glass L 1977 Science 197 287–289
  * [54] Lorenz E N 1963 J. Atmos. Sci. 20 130
  * [55] Takano K, Sugano C, Inubushi M, Yoshimura K, Sunada S, Kanno K and Uchida A 2018 Opt. Express 26 29424
  * [56] Sugano C, Kanno K and Uchida A 2020 IEEE J. Sel. Top. Quantum Electron. 26 1500409 ISSN 1558-4542
  * [57] Argyris A, Cantero J, Galletero M, Pereda E, Mirasso C R, Fischer I and Soriano M C 2020 IEEE J. Sel. Top. Quantum Electron. 26 5100309 ISSN 1558-4542
  * [58] Huang G B, Zhu Q Y and Siew C K 2004 IEEE proceedings 2 985–990
  * [59] Globisch B, Otto C, Schöll E and Lüdge K 2012 Phys. Rev. E 86 046201
  * [60] Lingnau B, Turnwald J and Lüdge K 2019 Phil. Trans. R. Soc. A 377 20180124
  * [61] Lüdge K and Lingnau B 2020 Laser Dynamics and Delayed Feedback, ISBN 978-1-07-160420-5 (New York, NY 10004: Springer Nature) pp Synergetic–18 Encyclopedia of Complexity and Systems Science 1st ed ISBN 978-1-07-160420-5
  * [62] Pauwels J, Van der Sande G, Verschaffelt G and Massar S 2021 Entropy 23 ISSN 1099-4300