In reservoir computing, an input signal drives a dynamical system, the reservoir, and the state variables of this system (evolving in time) are used for generating some desired output, like a temporal forecast of the input signal or a cross-prediction of some other signal that is assumed to be related to the input signal.

The concept of reservoir computing can be implemented and studied in continuous or discrete time. The discrete time framework is used, for example, if a time series sampled in time is input of a dynamical system simulated on a digital computer (e.g., a recurrent network). With discrete time, the reservoir system reads²

where r[n]∈ℝNres is the state of the reservoir system and u[n]∈ℝNin the input time series, that is often given by sampling a time continuous signal u[n] = u(nT_smp) where T_smp denotes the sampling time.

If reservoir computing is implemented on an analog computer or using any driven real-world (physical, biological, etc.) system, a formal description using continuous time systems [i.e., ordinary differential equations (ODEs)] is in general more appropriate. If the input u(t)∈ℝNin is given by N_in signals in time, the reservoir system reads

where r∈ℝNres is the state of the reservoir system and N_res its dimension.

The output v∈ℝNout of the reservoir system is given by a parameterized readout function g(u, r) of the reservoir state r and the input u.³

During training (or learning), the parameters of g are the only quantities to be estimated and adapted. Since the functional form of the function g should enable efficient training, functions that are linear in their parameters are preferred for reservoir computing design, such as linear superpositions of K_j basis functions bjk to define the component v_j of the output v as

where the coefficients w_kj constitute the parameters to be estimated using training data (see Section 2.3). The basis functions bjk can be any linear or non-linear function [e.g., polynomials [43] and radial basis functions [46]], but in most applications affine linear or quadratic polynomials have been used as readout functions of reservoir systems. Affine linear output functions

are given by an N_esv = 1 + N_in + N_res dimensional extended state vector x=(1;u;r)∈ℝNesv and an N_esv×N_out output matrix W^out whose elements are the parameters to be estimated during training [here (·;·;·) denotes concatenation and x is a row vector]. Equation (4) provides the N_out dimensional output signal (or time series) v where the j-th column of W^out contains superposition coefficients w_kj for computing the element v_j of the row vector v and the number of basis functions K_j equals N_esv for all j = 1, …, N_out components of the output v.

In general, the number of basis functions K_j can be different for each component v_j and can be larger than the dimension N_esv of the extended state vector x, for example, if (additional) non-linear basis functions like polynomials [43] or radial basis functions are used. Such non-linear output functions not only extend the possibilities of approximating complicated functional relations, but can also be used to break unwanted symmetries of the reservoir system [47].

In this and the following sections, the reservoir systems are assumed to be finite dimensional. However, also dynamical systems with an infinite dimensional state space, like delay differential equations or extended systems described by partial differential equations may act as dynamical reservoirs. For example, reservoir computing using feedback loops (e.g., ring lasers) has been demonstrated experimentally and in numerical simulations. The dynamics of these delay-based reservoir systems⁴ is governed by delay differential equations⁵ ṙ(t) = f(r(t), r(t−τ), u(t)) where τ equals the time delay due to the internal feedback loop (round-trip time of the light in case of a ring laser, for example). The output v(t)=∑n=1Nwnr(t-τn-NN) is obtained by sampling the solution r(t_n) at discrete times tn=t-τN-nN along the delay line. For a more detailed description of this type of delay-based reservoir computing, see the pioneering work of Appeltant et al. [30, 48], a recent review by Chembo [49], references [31, 32, 50, 51] for implementations and extensions, and Hülser et al. [52] for a detailed discussion of the influence of (multiple) delay times as successfully used by Stelzer et al. [53]. We will not go further into delay-based reservoir computing, but most of the aspects and features discussed below also apply to this approach.

Reservoir computing can be used for essentially all machine learning tasks including:

In all these cases, the reservoir system is trained such that it provides a task-dependent target value (the training is described in more detail in Section 2.3). Tasks (i) and (ii) can also be combined to learn and provide a signal, which is then used to control another dynamical process. Using different readout functions, a single reservoir system can perform different tasks like (i), (ii), and (iii) in parallel.

There are two main options for predicting the evolution of a target value over time: (a) performing the desired prediction in a single step or (b) decomposing the prediction interval into many small steps, resulting in an iterated prediction for discrete systems and an autonomous ODE in the case of a continuous reservoir system. Both approaches are now described in more detail.

The only parameters to be determined during training are the elements of the output matrix W^out. If the output function g is chosen to be linear in these parameters, training consists essentially in solving a linear set of equations as will be shown in the following. First, a training set {(v[m], y[m])} consisting of m = 1, …, M sampled output values v[m] and target values y[m]∈ℝNout has to be generated, where v[m] = v(mT_smp) and y[m] = y(mT_smp) in the case of a continuous reservoir system sampled with a sampling time T_smp. To make sure that the dynamical evolution of the reservoir system is (asymptotically) independent from its (in general unknown) initial conditions, this sampling has to be preceded by a sufficiently long transient (or washout) time T_trs = M_trsT_smp of M_trs time steps, a requirement that will be discussed in more detail in Section 2.4.

The output of the reservoir system matches the target if W^out can be chosen such that the linear equation

holds, where X is an M × N_esv matrix whose rows are extended states x[m] and Y is an M × N_out matrix whose rows are the corresponding target values y[m]. If M > N_esv, this linear set of equations has more equations (M·N_out) than unknowns (N_esv·N_out), while for M < N_esv, it is underdetermined. In the latter case, regularization is required to avoid solutions W^out with very large matrix elements, but also for the over-determined case M > N_esv, the solution may suffer from instabilities in the case that some columns of X are almost collinear, i.e., that elements of u and/or r are highly correlated. This problem can be addressed by means of Tikhonov regularization (also known as Ridge Regression) where a cost function

with

is minimized that contains a regularization term punishing too large matrix elements Wijout. The regularization parameter β>0 controls the impact of the regularization and has to be chosen suitably. The minimum of this cost function (Equation 8) is given by

where X^T denotes the transposed of X. Alternatively, the solution W^out (Equation 9) minimizing the cost function (Equation 8) can be computed using the singular value decomposition (SVD) of X = USV^T, where S is a diagonal matrix with S_ij ≥ 0 for i = j and S_ij = 0 for i ≠ j. In this case, the minimum of the cost function (Equation 8) is given by

with

Instead of Tikhonov regularization, one can also use the Moore-Penrose-pseudo inverse X ⁺ of X, such that

Using the SVD of X = USV^T, the pseudo-inverse X ⁺ can be written as X ⁺ = VEU^T with

where ϵ>0 represents the numerical resolution but can also be chosen larger and then serves as a regularization parameter.⁸ The difference between the definitions of the matrices D and E is illustrated in Figure 1. In both cases, the impact of very small or vanishing singular values S_ii is limited by bounding their inverse values 1/S_ii (see Equations 10 and 11). For sufficiently long training data, i.e., with M > N_esv, the pseudo-inverse based on Equation (11) may provide better results than Tikhonov regularization (Equation 10), because the latter has a stronger bias for small but not too small singular values.

An alternative to minimizing the quadratic cost function (Equation 8) is L1-optimization, where the Euclidean norm ∥·∥₂ is replaced by the L1-norm ∥·∥₁, and the solution of the minimum of the resulting cost function can be computed using the Lasso algorithm [66–68]. In this case, the computation of the optimal coefficients is more time consuming, but since many of them are zero after training, the calculation of the readout consists of fewer operations and is therefore more efficient.

An important feature of a trained reservoir system is its capability for generalization on similar but previously unseen data. This is notoriously difficult to determine in advance, but for reservoir systems consisting of (a certain class of) recurrent neural networks, Han et al. [69] derived a bound for the generalization error.

To use a trained reservoir system in any application, e.g., prediction, its output has to be reproducible and unique if the same input signal is presented again, independently from the (often unknown) initial state of the reservoir at the beginning. In reservoir computing, this feature is known as echo state property (ESP) or fading memory, and it is formally defined as follows [4, 70, 71].

For all input sequences u[0],u[1],u[2],...∈U⊂ℝNin from a compact set U and any pair of initial conditions of the reservoir system, r[0] and r~[0], the iterated application of Equation (1) results in limn→∞∥r[n]-r~[n]∥=0 (or correspondingly limt→∞∥r(t)-r~(t)∥=0 for continuous time systems (Equation 2) with input signal u(t)).

If the reservoir system fulfills this criterion, it asymptotically “forgets” its initial state after some transient or washout time T_trs, and its (current) state is given by a so-called echo state function ξ of the semi-infinite input sequence r[n] = ξ(…, u[n−1], u[n]) [4].

If the input signal was generated by a dynamical system z.=fsrc(z) or z[n + 1] = f_src(z[n]), this system and the reservoir system can be seen as a pair of uni-directionally coupled dynamical systems (i.e., a drive-response configuration), where the echo state property occurs due to generalized synchronization [72–80]. With generalized synchronization, the state of the driven system, i.e., the reservoir, r(t) is (asymptotically) a function ψ(z(t)) of the state z(t) of the driving system generating the input signal, i.e., limt→∞∥r(t)-ψ(z(t))∥=0 (or limn→∞∥r[n]-ψ(z[n])∥=0 for discrete systems). In this way, each state variable r_i(t) of the reservoir system represents after a sufficiently long (synchronization) transient time T_trs, a dynamically emerging function ψ_i(z) of the state of the driving system that is used in the output Equation (4) as a non-linear basis function for approximating the target v [39]. This function ψ_i(z) corresponds to the echo state function r[n] = ξ(…, u[n−1], u[n]), if the input u is generated by the dynamical system and given by an observation function u = h(z). In this case, (previous) values of the input time series are given by functions of the current state z[n] of the driving system, because u[n−k] = h(z[n−k]) = h(ϕ^−k(z[n])), where ϕ^−k denotes the inverse flow of the driving system. The same holds for continuous reservoir systems driven by an input signal from a continuous dynamical system.

So in principle, any driven dynamical system may be used as a reservoir system if it exhibits generalized synchronization with the dynamical system generating the input signal. This includes analog computers or even real physical processes. The occurrence of generalized synchronization can by checked, for example, using conditional Lyapunov exponents or the auxiliary system approach [74, 76, 79]. However, generalized synchronization is only a necessary condition for fulfilling the echo state property but is not sufficient to guarantee the same unique response for all initial conditions of the reservoir system. There may be different coexisting attractors that lead to different reproducible asymptotic responses depending on the initial conditions of the driven system. This form of multistability has to be taken in to account [63–65] or excluded to guarantee a unique response and thus the echo state property.

The emerging functions ψ_i are continuous, and their smoothness properties depend on the (transversal) contraction properties of the response system. In general, the stronger this contraction, i.e., the faster two neighboring trajectories r(t) and r~(t) converge with limt→∞∥r(t)-r~(t)∥=0, the smoother are the functions ψ_i. A mathematically more rigorous treatment of this feature can be found, for example, in a series of publications by Stark et al. on forced dynamical systems [81–83]. A detailed mathematical analysis of the differentiability properties of the echo state functions and their importance for successful reservoir computing was carried out by Grigoryeva and Ortega [84]. Hart et al. [85] proved under some mild assumption that ψ almot surely provides an embedding of the dynamics of the driving system into the state space of the reservoir system, and that there exists a linear readout function φ(r) that predicts arbitrarily well the input values provided that the dimension of the reservoir system is high enough. If this prediction is used to substitute the input u of the reservoir system an autonomous system r[n + 1] = f(r[n], φ(r[n])) or r.=f(r,φ(r)) is obtained (see Section 2.2.2). A well-trained autonomous reservoir system reproduces the Lyapunov exponents of the driving system, i.e., the Lypunov exponents of the input signal are a subset of the Lyapunov exponents of the autonomous reservoir [42]. In this context, Platt et al. [86] found that the prediction time for autonomous reservoir systems with a low Kaplan–Yorke dimension is the longest. Carroll [87] obtained a similar result for the Kaplan–Yorke dimension calculated from the full Lyapunov spectrum of the coupled drive-response system.

If the response system exhibits no contraction, generalized synchronization and the echo state property cannot occur. If, on the other hand, contraction is too strong, the emerging functions ψ_i are very (too) smooth and, therefore, not very well-suited to be used for approximating functional relations of target signals y, whose dependence on the state z is given by a less smooth and more complicated function. It is therefore important that the contraction properties of the response system can be adapted to generate the appropriate level of smoothness of the output for a given modeling task. In case of echo state networks (see Section 2.5), contraction is controlled by the spectral radius of the reservoir matrix W^res, and in many studies, it has been observed that weak contraction provides best performance of the reservoir system. Weak contractions lead to long synchronization transients, and an increased memory capacity of the reservoir system quantified by the ability to recover past input signals u[n−k] from the current state r[n] of the reservoir system [88]. A comprehensive discussion of the role of generalized synchronization and the question under which conditions iterated prediction of a reservoir reproduces the attractor (i.e., the climate) of the driving system was given by Lu et al. [44]. Mathematical criteria and statements for the latter can also be found in [85].

Dambre et al. [89] pointed out that the information processing capacity of a driven dynamical system with fading memory (i.e., echo state property) and linear readout is bounded by the number of linearly independent functions of the input the system can compute (or generate). In general, good performance of any reservoir system can be expected if the time series r_i[n] given by the state variables of the reservoir system are linearly independent, because in this case the dimension of the linear subspace spanned by the columns of the matrix X is large. This feature can be quantified by the rank of the covariance matrix X^TX, also called covariance rank [90–93]. A large covariance rank will not guarantee good performance, but in general, a reservoir with a large rank will be superior to a reservoir generating output with a lower covariance rank.

An extension of generalized synchronization is consistency [94], which quantifies the degree of functional dependency of a driven dynamical system on its input in terms of the linear correlation of the reservoir variables when repeatedly driven by the same input. If the echo property is fulfilled, the consistency is equal to one. Consistency has also been used for studying reservoir computing, in particular with experimental set-ups [31, 32]. Lymburn et al. [95] found that even in cases of non-complete consistency, such inconsistent echo state networks may exhibit very good performance because there are subspaces in the state space of the reservoir system that represent signals with high consistency. This observation may lead to an enhanced understanding of the computational capacity of chaotic or noisy response systems and novel unsupervised optimization procedures. Extending the basic concept of consistency, Jüngling et al. [96] have developed an approach to investigate the propagation and distribution of information-carrying signals in a reservoir system.

Last but not least, reservoir computing using echo state networks has also been used to approximate the emergent functions that arise from the generalized synchronization of different coupled systems [97].

An echo state network [4, 98] is a recurrent neural network acting as a dynamical reservoir that can be given by a set of ordinary differential equations (ODEs)

or a set of updating rules for echo state networks that operate in discrete time

In both cases, q(·) denotes an activation function (applied component-wise), which in most applications is a sigmoid function like q(·) = tanh(·). b∈ℝNres is a vector of different bias constants b_i∈ℝ affecting the dynamics of the state variable r_i of the i-th node of the network and may be used to break the symmetry of the dynamical equations [47]. The constants γ > 0 and α ≥ 0 are called leaking rates and can be used to adjust the temporal scale of the echo state network to that of the input signal and also to control the (fading) memory of the network dynamics. The parameter σ > 0 controls the amplitude of the input signal. Both matrices, Wⁱⁿ and W^res, are chosen randomly. Since the impact of the input can be controlled by the parameter σ, the elements of the input matrix Wⁱⁿ can, for example, be uniformly distributed in the interval [−1, 1]. The internal coupling matrix of the reservoir W^res and the leaking rates γ or α have to be chosen in a way that the echo state network possesses the echo state property [99].

Furthermore, in many cases, sparse coupling matrices W^res provided better performance. Therefore, a sparseness parameter δ is introduced quantifying the relative number of non-zero elements of W^res. With a sparse coupling topology, the echo state network decomposes into many loosely coupled subsystems, and this may result in a more diverse dynamics of its nodes [100].

All network parameters mentioned so far, N_res, W^res, Wⁱⁿ, b, σ, γ, α, δ, and the regularization weight β, are hyperparameters, which have to be chosen by the user when designing the echo state network and may have a major impact on its performance. Optimal values of hyperparameters can be found by grid search or using gradient methods [70, 101] or Bayesian optimization [102]. Improved predictive performance and enhanced generalization capabilities have been obtained with first projecting the state of the reservoir system to a lower dimensional space before computing the parameters of the readout function [103]. An approach for reducing the size of the network using proper orthogonal decomposition (or principle component analysis) of its internal dynamics was suggested by Jordanou et al. [104].

An important feature of reservoir systems is their memory capacity, i.e., their ability to reconstruct and output previous input signals from the past. The memory capacity of an echo state network increases with the network size and is also affected by the other hyperparameters [detailed studies using i.i.d. random input time series can, for example, be found in [88, 105, 106]]. For a given input signal and task, there is an optimal amount of memory of the reservoir system maximizing its performance, but it is difficult to determine it in advance without simulating the full reservoir computer [107].

Reservoir systems do not have to be non-linear, because linear system may also generated non-linear output due to their exponential solution functions [108]. Bollt [109] showed that an echo state function with a linear activation function is equivalent to a vector autoregressive (VAR) model and if a quadratic readout function is used this is equivalent to non-linear VAR (NVAR) models. Furthermore, he showed that linear echo state networks are less effective but in combination with quadratic readout prediction results for chaotic time series may be acceptable. Carroll and Pecora also found that linear node dynamics in echo state networks is less effective [90]. In their study of the impact of network structure, they considered the rank of the covariance matrix of the output, symmetries, and the resulting memory capacity and found that sparsity of the internal coupling matrix is beneficial, but much more important for improving the performance may be flipped links in the network. To cope with dynamics on multiple time scales, Tanaka et al. [110] considered echo state networks with heterogeneous leaky integrator (LI) neurons, where the parameter α in Equation (12) is chosen different for each node.

The performance of echo state networks (and reservoir computing in general) also depends on the available input u, more precisely whether full or only partial information about the state of the driving system is provided (e.g., a single component of the state vector of the dynamical system generating the data). Storm at al. [92] demonstrated that if only partial information is available, the reservoir system needs some memory to generated implicitly a delay embedding. This is not necessary if the full state is presented as input, and in this case, even an echo state network without internal dynamics can provide good predictions [62], because static non-linear maps of the input state (given by the activation functions) are sufficient. However, Lymburn et al. [97] have shown that the inclusion of internal dynamics can still be advantageous compared with memoryless reservoir computing, even if the complete state is available as input, if the input–output function to be learned is itself dynamically generated and has a very complex (e.g., fractal) form.

Gauthier et al. [111] compared echo state networks with NVAR modeling using polynomial basis functions and found that these types of NVAR models provide similar or even better predictions, while requiring shorter time series for training and less computational resources.⁹ Jaurigue and Lüdge [115] discussed this approach in the context of other machine learning methods including delay-based reservoir computing with multiple delay loops [53]. An evaluation of the performance of different echo state networks with different network structures compared with other machine learning methods was presented by Shahi et al. [116]. They found that reservoir computing and NVAR models achieve similar prediction results.

The performance of a reservoir computing system can be improved by increasing its state-space dimension N_res, which results in a larger number of basis function ψ_i used for generating the output (Equation 4) [89]. This increases, however, the computational load for simulating the reservoir system. An alternative is the use of delayed input and state variables, u(t − τ_u) and r(t − τ_r), respectively, because they provide additional basis functions. Let ϕ^t be the flow of the driving dynamical systems such that z(t-τu)=ϕ-τu(z(t)), and let the input signal of the reservoir system be given by u(t) = h(z(t)) where h denotes the observation function. In this case, after the synchronization transient (washout time) decayed, delayed input and reservoir state variables can be (asymptotically) written as functions of the state z(t) of the driving system at time t

Here, the functions h◦ϕ-τu and ψ◦ϕ-τr, represented by u(t−τ_u) and r(t−τ_r), respectively, can be used to extend the set of basis functions entering in the computation (Equation 4) of the output of the reservoir system.¹⁰

The extension of reservoir computing using past states of the reservoir system has been first studied by Marquez et al. [117], Sakemi et al. [118], and Del Frate et al. [119] who all found significant improvement in the performance of the reservoir system, a result that will be confirmed in the following. Later, Carroll and Hart [91] showed for an opto-electronic delay-based reservoir computer with only a small number of virtual nodes that time-shifted readout increases the rank and memory of the reservoir computer. In general, using many different delay times τ_u and τ_r increases the number of columns of the matrix X and (if suitable) its covariance rank, which leads to an improved performance of the reservoir [91, 93].

Duan et al. [120] showed that linear readout functions including delayed variables may also provide an embedding of the input dynamics, and they discuss the trade-off between the number of delays and the state space dimension of the reservoir system (i.e., the number of nodes in case of a recurrent network).

Jaurigue et al. [121] showed that the performance of an unoptimized reservoir system can significantly be improved by adding a time-delayed input signal. In their study, they simulated time multiplexed delay-based reservoir computing with a single non-linear element that can be implemented using optical hardware. In a recent study, Jaurigue and Lüdge [122] showed that the dependance of reservoir performance on the choice of suitable hyperparameter values can be reduced by including delayed variables.

To demonstrate the effectiveness of reservoir computing using delayed variables, we will consider now a very simple example: a three-dimensional reservoir system based on the Lorenz-63 system [79, 123]

This system fulfills the echo state property for all initial conditions r(0) (see Appendix).

The parameter μ is used to adapt the time scale to that of the driving signal u(t) and controls the transversal contraction properties of the response system (see Appendix). Therefore, it is a hyperparameter that has to be carefully chosen.

Using an affine linear readout function, the output component v_j (j = 1, 2, 3) sampled at times t_n = nT_smp is given by

or using the notation v_j[n] = v_j(t_n)

The minimal reservoir system (Equation 15) will be used to predict the variables z_i of a chaotic Rössler system [79, 124]

from the observations u[m] = u(t_m) = z₁(t_m) of the first variable at times t_m = mT_smp, where T_smp denotes the sampling time.

In our first example, the reservoir system (Equation 15) is used to cross-predict the variables z₂(t_m) and z₃(t_m) of the chaotic Rössler system (Equation 17) from previous observations u[k] = u(t_k) = z₁(t_k) of the first variable z₁ with k ≤ m. This task is known in control theory as state observer. If a training time series u(t_k) of length M is available,¹¹ the computation of the weights w_ij in Equation (16) consists of solving the set of linear equations (Equation 7) with

where

and

with the columns of Y given by the target values z_i[m] = z_i(t_m) of the Rössler system.¹² This set of linear equations (Equation 18) is solved for the weights W^out = X ⁺ Y using the Moore-Penrose pseudo-inverse X ⁺ of X (see Section 2.3) with ϵ=10-7σ1, where σ₁ denotes the largest singular value of X.

Figure 2 shows predictions of z₂ and z₃ obtained with a training time series of length M = 50, 000 and sampling time T_smp = 0.01 corresponding to a time interval of length 500. While the cross-prediction of the z₂ component may be acceptable, the peaks of the z₃ variable are not correctly reconstructed. This relatively poor result should be expected given the fact that only time series of u, r₁, r₂, and r₃ have been used to predict the target signal. The bottom panels of Figures 2C, F show the dependence of the Normalized Root Mean Squares Error (NRMSE) of the test data

on the parameter μ for both target signals y₂ and y₃, respectively. Here, ȳ_j stands for the mean value of y_j, σyj2 represents its variance, and N_test is the length of the time series used to compute the test error. The NRMSE quantifies the performance of the reservoir compared with a prediction of y_j by its mean value ȳ_j. The lowest errors occur for predictions of z₂ at μ = 0.025 and for z₃ near μ = 1.1. This large difference in optimal μ values may be due to different smoothness properties of the functions ψ_i required to successfully predict z₂ or z₃ (see Section 2.4 and Appendix). Such differences with respect to optimal (hyper) parameters may complicate the simultaneous prediction of different target variables using output from the same reservoir system.

We will now increase the pool of signals (or basis functions) using past values of the input signal (Equation 13) and the state variables (Equation 14) of the reservoir system (Equation 15) in terms of delay vectors u(t) = [u(t − N_uτ_u), …, u(t − τ_u), u(t)] and r_k(t) = [r_k(t − N_rτ_r), …, r_k(t − τ_r), r_k(t)], which can for t = mT_smp, τ_u = L_uT_smp and τ_r = L_rT_smp be written as¹³

To simplify the presentation, we use the same number of N_u = N_r = N_delay delayed variables u and r in the following. Using the row vectors u[m] and r_k[m], the linear set of equations for computing the parameters of the affine linear output function reads

with K = 1 + 4(N_delay + 1) = 5 + 4N_delay basis functions. This set of linear equations can again be solved using the pseudo-inverse of X or any other method for solving a linear set of equations using regularization (see Section 2.3). Note that with delay M + N_delay samples at times t_1−Ndelay, …, t₁, …, t_M−1, t_M are required for constructing the matrix X and that the number of basis functions K = 5 + 4N_delay used for approximating the target function increases linearly with N_delay. Furthermore, the rows of matrix X represent a projection of the state (z, r) of the full coupled system to a K-dimensional space (providing an embedding if K is large enough).

Figure 3 shows prediction results for z₂ and z₃ obtained with N_delay = 100 and delay times (τ_u = 0.32, τ_r = 0.37) for z₂ and (τ_u = 0.03, τ_r = 0.30) for z₃, such that the delay vectors for z₂ cover approximately a period of time of 32 < N_delayτ_{u, r} < 37 corresponding to about six oscillations of z₁ and z₂. As can be seen, not only the z₂-component of the Rössler system is recovered with high precision, but also the peak structure of the input signal and the z₃-component.

The delay times τ_u, τ_r and the number of delays N_u and N_r are additional hyperparameters that have to be carefully chosen for a given prediction problem. One method for optimization is the maximization of the rank correlation [93], which can also be achieved by a targeted use of the zeros or minima of the autocorrelation functions of the input signal and the reservoir states.

So far, only the current values of the target signals z₂(t) and z₃(t) have been reconstructed using input signals z₁(t) until time t. As discussed in Section 2.2, reservoir systems can also be used for predicting the future evolution of the variables of interest. Figure 4 shows prediction results where the output of the reservoir system provides the target signals z₂(t + T_prd) and z₃(t + T_prd) for a prediction time of T_prd = 10 (in units of the time variable of the Rössler system). As indicated by the green horizontal bars in Figures 4A, C, this period of time corresponds to about 1.6 oscillation periods of the z₂ variable. The continuously available future values of the target signals can be used, for example, for implementing a control system or for predicting the occurrence of “extreme events” like the peaks of the z₃ variable.

A prediction scheme based on delayed variables can also be implemented using feedback (see Section 2.2.2). In this case, the output which is fed back as input to the reservoir is a function of delayed state variables v(t) = g(r(t), r(t−τ₁), r(t−τ₂), ….), and the resulting autonomous reservoir system is given by a delay differential equation r.=f(r(t),g(r(t),r(t-τ1),r(t-τ2),....)) with multiple delays.