Samuele Mazzi

David Zarzoso

Parametric Validation of the Reservoir-Computing-Based Machine Learning Algorithm Applied to Lorenz System Reconstructed Dynamics

Keywords: Reservoir Computing, Lorenz system, Hyperparameters, Error quantification, Machine Learning

A detailed parametric analysis is presented, where the recent method based on the Reservoir Computing paradigm, including its statistical robustness, is studied. It is observed that the prediction capabilities of the Reservoir Computing approach strongly depend on the random initialisation of both the input and the reservoir layers. Special emphasis is put on finding the region in the hyper-parameter space where the ensemble-averaged training and generalization errors together with their variance are minimized. The statistical analysis presented here is based on the Projection on Proper Elements method.

Introduction

The numerical modelling of multi-scale and complex systems has been continuously growing in the last decade, reaching very high levels of fidelity and accuracy. Some examples can be found in the modelling of turbulence in nuclear fusion plasmas [START_REF] Fasoli | Computational challenges in magnetic-confinement fusion physics[END_REF], in the prediction of weather and climate in the atmosphere [START_REF] Lynch | The origins of computer weather prediction and climate modeling[END_REF], in the multi-scale methods in systems biology [START_REF] Dada | Multi-scale modelling and simulation in systems biology[END_REF] and in the analysis of the financial market [START_REF] Xiao | A multiscale modeling approach incorporating ARIMA and ANNs for financial market volatility forecasting[END_REF]. This growing phase has inevitably led to the requirement of extremely efficient numerical resources, with the current development of exa-scale high-performance-computing systems, in order to run such very demanding simulations. In some practical applications, however, the elapsed computational time for the exploitation of these numerical modellings represent the bottleneck of the aforementioned validation and/or prediction analyses. Furthermore, in parallel with the improvement of the employed numerical techniques, the incorporation of newly discovered physical processes or improved complex multi spatiotemporal scale interactions into the models makes the development of classical numerical methods even more challenging. For these reasons, the application of the rapidly advancing Machine Learning (ML) algorithms has gained a lot of attraction in the framework of the datadriven modelling of complex dynamical systems. The advantages of using such an approach goes from the acceleration of the simulations due to the highly reduced computational cost, to the quality improvement of the physical assessment of the numerical outcomes.

Within the several possible ML techniques applied to the forecasting of complex systems, the Reservoir Computing (RC) has shown appealing features, such as the simplicity and the speed of the training process, and revealed to be able to achieve state-of-the-art results, compared to other Artificial Intelligence Systems (AISs). The RC method takes inspiration from the paradigm of the reservoir computers, which were developed in two independent studies of the early current century [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF][START_REF] Maass | Real-time computing without stable states: A new framework for neural computation based on perturbations[END_REF]. Due to the difficulty of designing efficient recurrent neural network (RNN) architectures and the subsequent inaccuracy of simplistic RNNs [START_REF] Xue | Decoupled echo state networks with lateral inhibition[END_REF][START_REF] Ma | Multistep-prediction of chaotic time series based on co-evolutionary recurrent neural network[END_REF], Jaeger and Maass proposed the approach of the reservoir computing, which differs from the standard RNN in being essentially split into two levels. The internal layer is composed by the so-called reservoir, which is a randomly initialised RNN. The output level, usually called readout, is a feed-forward layer which realizes an optimized output function in order to obtain the predicted vector. In our case, as in the case of the RC proposed by Jaeger (the Echo State Network (ESN) [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF][START_REF] Jaeger | Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication[END_REF][START_REF] Schrauwen | An overview of reservoir computing: theory, applications and implementations[END_REF][START_REF] Lukoševičius | Reservoir computing approaches to recurrent neural network training[END_REF][START_REF] Lukoševičius | A practical guide to applying echo state networks[END_REF][START_REF] Gauthier | Reservoir computing: Harnessing a universal dynamical system[END_REF][START_REF] Verstraeten | An experimental unification of reservoir computing methods[END_REF]), the output function is obtained from a linear regression. The weights of the connected nodes of the reservoir are fixed randomly, whereas the readout layer is trained to obtain its weights. For more details on the application of the RC technique on the forecasting of dynamical systems one can refer to recent contributions such as, e.g., Refs. [START_REF] Xi | Analyzing the state space property of echo state networks for chaotic system prediction[END_REF][START_REF] Song | Multi-steps prediction of chaotic time series based on echo state network[END_REF][START_REF] Bianchi | An overview and comparative analysis of recurrent neural networks for short term load forecasting[END_REF][START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Lu | Reservoir observers: Model-free inference of unmeasured variables in chaotic systems[END_REF][START_REF] Mcdermott | An ensemble quadratic echo state network for non-linear spatio-temporal forecasting[END_REF][START_REF] Lu | Attractor reconstruction by machine learning[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF][START_REF] Zimmermann | Observing spatio-temporal dynamics of excitable media using reservoir computing[END_REF][START_REF] Tanaka | Recent advances in physical reservoir computing: A review[END_REF][START_REF] Chattopadhyay | Data-driven predictions of a multiscale lorenz 96 chaotic system using machinelearning methods: Reservoir computing, artificial neural network, and long short-term memory network[END_REF][START_REF] Bianchi | Reservoir computing approaches for representation and classification of multivariate time series[END_REF]. Furthermore, theoretical efforts have been spent throughout the last decade to study the underlying mechanism of the RC algorithm [START_REF] Verstraeten | Memory versus non-linearity in reservoirs[END_REF][START_REF] Marzen | Difference between memory and prediction in linear recurrent networks[END_REF][START_REF] Massar | Mean-field theory of echo state networks[END_REF][START_REF] Goudarzi | Memory and information processing in recurrent neural networks[END_REF][START_REF] Grigoryeva | Echo state networks are universal[END_REF][START_REF] Hart | Embedding and approximation theorems for echo state networks[END_REF][START_REF] Tino | Dynamical systems as temporal feature spaces[END_REF][START_REF] Hart | Echo state networks trained by Tikhonov least squares are l2 (µ) approximators of ergodic dynamical systems[END_REF]. Among those results, detailed insights prove that a configuration of the weights of the readout layer W out exists for which the RC-ESN output dynamics is topologically connected to the input chaotic system [START_REF] Hart | Embedding and approximation theorems for echo state networks[END_REF][START_REF] Hart | Echo state networks trained by Tikhonov least squares are l2 (µ) approximators of ergodic dynamical systems[END_REF]. However, the capability of the RC algorithm to find those weights, for which the difference with the input data in the training phase is minimized, still remains elusive.

The present study is motivated by recent contributions [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF] in which RC paradigm is used to predict the dynamics of the well-known Lorenz chaotic system [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]. The Lorenz system is given by the equations: These equations were derived for the first time by Lorenz, after simplifying a set of partial differential equations describing the motion of a fluid in between two layers of the atmosphere. Therefore, the parameters a, b and c have a real physical meaning, especially b, which correspond to the convection rate. The purpose of this study is to validate and evaluate the accuracy of the RC method applied for reconstructing the dynamics of the Lorenz system. The reported statistical analyses, based on the Projection of Proper elements method [START_REF] Cartier-Michaud | Code verification and model reduction: Application to transport in turbulent plasmas[END_REF][START_REF] Cartier-Michaud | Projection on proper elements for code control: Verification, numerical convergence, and reduced models. Application to plasma turbulence simulations[END_REF][START_REF] Cartier-Michaud | Verification of turbulent simulations using PoPe: quantifying model precision and numerical error with data mining of simulation output[END_REF], demonstrate that the large randomness of the RC approach is strongly related to the random initialisation of the reservoir network components and is highly affected by the chosen configuration parameters. Large variation of the error in the prediction phase is obtained for a broad set of parameters of the reservoir, leading to strong inaccuracy in predicting the time evolution of the Lorenz trajectories for a significant sample of realisations. The results presented in recent contributions [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF] are shown to fall within a narrow region of the scanned parameter series. Nevertheless, the prediction error of such a narrow favorable region is affected by large statistical variation as well, making the reconstructed dynamics by this specific ML technique frequently not reliable.

     ẋ = F x (x,
The remainder of the paper is organized as follows: the specific RC approach that has been applied on the Lorenz system [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] and the evaluation techniques, including the Projection on Proper elements method [START_REF] Cartier-Michaud | Code verification and model reduction: Application to transport in turbulent plasmas[END_REF][START_REF] Cartier-Michaud | Projection on proper elements for code control: Verification, numerical convergence, and reduced models. Application to plasma turbulence simulations[END_REF][START_REF] Cartier-Michaud | Verification of turbulent simulations using PoPe: quantifying model precision and numerical error with data mining of simulation output[END_REF], are described in Section 2; Section 3 is divided into five subsections, illustrating the principal results that are achieved. Different hyper-parameters of the reservoir network and of the model are considered and deeply analyzed, with statistical insight of the probability distribution function of the errors; dedicated analysis on the time evolution of the Lorenz coefficients is carried out; the main conclusions are summarized in Section 4.

Methods

In this paper we produce the data set by solving the Lorenz system given by Equations 1 using a 4 th order Runge-Kutta method (RK4). This integration will produce the data set {x i , y i , z i }, with N the maximum number of time steps. We will use ∆t = 0.02 and N = 6250. Such data set will be divided into a training set and a prediction set, with respective lengths n train = 5000 and n pred = 1250, such that n tot = n train + n pred . The training set will be used to compute the parameters of the AIS and the prediction set will be used to test the predictive capabilities of the AIS. Here we follow exactly the methodology employed by Pathak et al. in Refs. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF] with the same set of parameters they used to build the RC. The architecture of the AIS is shown in figure 1 and will be briefly explained in the following. The input vector is called u(t) ∈ R 3 and contains the three components of the Lorenz system, x, y and z. It is fed into the reservoir, represented by a square matrix whose element (i, j) indicates the connection between the neuron i and the neuron j. Such matrix is called A ∈ R Dr×Dr , with D r the number of neurons. If the element A i,j = 0, the neurons i and j are not connected. Between the input and the reservoir, one needs an input layer W in ∈ R 3×Dr , mathematically representing an application from R 3 to R Dr . Both the reservoir matrix and the input layer are randomly initialised. The elements of both W in and A are drawn from a uniform distribution function in [-1,1], with an additional multiplying scaling factor γ scal = 0.1 applied only to the input layer elements. The elements of the reservoir are re-scaled so that the largest magnitude of the eigenvalues (the so-called spectral radius) is equal to the desired value ρ. In the present study, based on the example of Refs. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF], the reservoir matrix A is given by a sparse Erdős-Rènyi network [START_REF] Erdős | On random graphs I[END_REF][START_REF] Gilbert | Random graphs[END_REF], with an average degree of d = 6. It is to be noted that, since the degree of the Erdős-Rènyi network is fixed to d = 6, the rewiring probability p = d/D r is decreased (resp. increased) when D r is increased (resp. decreased) and therefore the connectivity matrix of the neural network becomes sparser (resp. less sparse). The results can be reproduced also with less weakly connected network [START_REF] Hart | Embedding and approximation theorems for echo state networks[END_REF], and therefore the choice of the sparse Erdős-Rènyi network is made only to be consistent with the aforementioned publications. As it will be presented in the following, the number of reservoir neurons and the spectral radius of the adjacency matrix are widely scanned and the impact of each couple of input parameters on the measured errors is then evaluated. It must be stressed that the Echo State Property [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF], which was introduced as a necessary feature for the effectiveness of the RC technique, was originally linked to the rule of thumb ρ < 1 for the reservoir network. However, this latter consideration has been proven not to be a necessary criterion for the ESP to hold [START_REF] Yildiz | Re-visiting the echo state property[END_REF][START_REF] Basterrech | Empirical analysis of the necessary and sufficient conditions of the echo state property[END_REF]. Hence, detailed studies in this framework [START_REF] Manjunath | Echo state property linked to an input: Exploring a fundamental characteristic of recurrent neural networks[END_REF][START_REF] Ceni | The echo index and multistability in input-driven recurrent neural networks[END_REF][START_REF] Verzelli | Learn to synchronize, synchronize to learn[END_REF] proposed less restrictive conditions, thereby allowing to explore the hyper-parameter space also for the spectral radius ρ > 1 of the adjacency matrix A. Justified by this, the analyzed range of ρ values is extended also beyond the unity, as it will be presented in the following.

Reservoir Computing Setup

Each neuron of the reservoir is characterised by a state r (t), which represents the activation of the neuron (-1: de-activated, +1: activated), computed at each time step by a tanh activation function:

r(t + ∆t) = tanh(A • r(t) + W in • u(t)) (2)
One can note that the recurrence in the RC paradigm is essentially considered in Equation 2, where the neuron states at the i th time step r i are linked to the neuron states at the previous time steps (r i-1 , r i-2 , r i-3 , and so on). This is also illustrated in the schematic view of the RC architecture of Fig. 1. The particular choice of using the activate function of Equation 2 is made in order to follow the example of Refs. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Lu | Reservoir observers: Model-free inference of unmeasured variables in chaotic systems[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF], nonetheless diverse RC architectures have been already exploited successfully [START_REF] Haynes | Reservoir computing with a single time-delay autonomous boolean node[END_REF][START_REF] Antonik | Brain-inspired photonic signal processor for generating periodic patterns and emulating chaotic systems[END_REF][START_REF] Larger | High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification[END_REF]. The underlying idea of the reservoir computing is to predict the dynamics of the system at the time step n + 1, i.e. u n+1 by means of the relation:

u n+1 = W out • r n + c out (3)
The weights of the readout layer W out ∈ R 3×n can be subsequently computed by minimizing the difference between the actual u(t) and the predicted v(t) Lorenz trajectories for each time step of the training phase. Therefore, in order to train the reservoir computing, one has to find W out and c ∈ R 3 minimizing the following quadratic form:

||W out • r + c -v|| 2 (4)
for all the time steps of the training phase, where ||q|| 2 = q T • q.

The so-called "ridge regression parameter" is set to zero, i.e. β = 0, since the over-fitting is already avoided [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF]. Moreover, it is to be noted that the assumption of symmetry in computing the first two components x → -x and y → -y of the Lorenz equations [START_REF] Lu | Reservoir observers: Model-free inference of unmeasured variables in chaotic systems[END_REF][START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Chattopadhyay | Data-driven predictions of a multiscale lorenz 96 chaotic system using machinelearning methods: Reservoir computing, artificial neural network, and long short-term memory network[END_REF] in the predicted vector v is retained, in order to be consistent with the assumption made in Refs. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF][START_REF] Lu | Reservoir observers: Model-free inference of unmeasured variables in chaotic systems[END_REF]. For this reason, Equation 4includes r instead of r. The modified neuron state matrix r presents the second half of the third row given by the nonlinear combination ri = r 2 i (for i > D r /2). However, the same analysis has also been carried out relaxing such symmetry assumption with very similar results. Then, the minimization problem described in Equation 4 reduces to a linear regression [START_REF] Jaeger | The "echo state" approach to analysing and training recurrent neural networks[END_REF], whose solution in the training phase is:

W * out = (U • R T ) • (R • R T ) -1 (5)
where U ∈ R 3×n is the array containing the actual dynamics of the Lorenz system, R ∈ R 3×n the array containing the neuron states (with the nonlinear combination as described above and in Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF]) and W * out is the particular solution of the minimization problem. As a result, the readout layer W out is trained and thereby fixed. Hence, the same procedure can be applied in the prediction phase. Indeed, the neuron states are now computed using the predicted vector by:

r(t + ∆t) = tanh(A • r(t) + W in • v(t)) (6)
and subsequently applied as:

v n+1 = W out • v n (7)
for the entire prediction phase. The vector v n+1 = [x, ỹ, z] n+1 contains the components of the Lorenz dynamics predicted by the RC approach at the time step n + 1.

With respect to traditional RNNs, the RC architecture benefits from the possibility of training the readout layer at once, which reduces significantly the elapsed time for the training phase [START_REF] Vlachas | Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics[END_REF], and from the absence of the vanishing gradient issue [START_REF] Pascanu | On the difficulty of training recurrent neural networks[END_REF][START_REF] Chattopadhyay | Data-driven predictions of a multiscale lorenz 96 chaotic system using machinelearning methods: Reservoir computing, artificial neural network, and long short-term memory network[END_REF].

Table 1 summarizes the principal parameters of the analyzed RC configurations.

Error Quantification through the PoPe Method

The AIS predictive capabilities were evaluated in Refs. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF][START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF] based on the ability to predict the climate, i.e. the ergodic properties of the chaotic system quantified by the Lyapunov exponents. In the present Letter, we employed a different criterion, based on the Projection on Proper elements (PoPe) using the Euclidean distance as the measure to quantify the error of the system to reproduce the Lorenz dynamics. The PoPe method is described in detail in Refs. [START_REF] Cartier-Michaud | Code verification and model reduction: Application to transport in turbulent plasmas[END_REF][START_REF] Cartier-Michaud | Projection on proper elements for code control: Verification, numerical convergence, and reduced models. Application to plasma turbulence simulations[END_REF][START_REF] Cartier-Michaud | Verification of turbulent simulations using PoPe: quantifying model precision and numerical error with data mining of simulation output[END_REF]] and here we provide only a brief explanation for our purposes.

Let us assume that we have a set of data d = {x i , y i , z i } 1≤i≤n pred provided by the RC. The PoPe method is employed to determine to Table 1: Reservoir parameters that has been used for the standard simulations, if not otherwise clearly noted. As explained in the main body of the Letter, the parameters are: β the ridge regression parameter, γ scal the scaling factor applied to the input layer W in , D r the reservoir number of neurons, ρ the spectral radius of the reservoir network, d the average degree of the Erdős-Rènyi network, ∆ t the time step of the considered Lorenz dynamics, n train and n pred the number of time step in the training and prediction phase respectively, and δ the period of the moving time window for the error calculation. Very similar parameters have been used also in Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF].

Parameter Value β 0 γ scal 0.1 D r [50, 610] ρ [0.5, 2.5] d 6 ∆ t 0.02 n train 5000 n pred 1250 δ 1.25
what extent the set of data d, produced by the the RC technique, has been generated by the Lorenz system given by Equation 1. This is done by computing the coefficients (a, b, c) that better fit d.

For this purpose, let us note that the Lorenz system can be rewritten as follows

     ẋ = a (y -x) + α x xz + y + ẏ = bx + α y -xy + ż = -cz + α z (8a) (8b) (8c)
where we have included (α x , α y , α z ) for the sake of generality. Such expression highlights the linear dependence on the parameters (a, b, c, α x , α y , α z ). Therefore, the calculation of the coefficients which better fit the RC predictions is reduced to a linear regression. To do this, the time derivatives on the LHS of Equations 8 are computed from d using an arbitrarily high precision method (in the present work we use a 4 th order finite difference scheme). These derivatives are noted x, ỹ, z . The Lorenz equations can then be rewritten in a matrix form, combining the RC-predicted trajectories d i with the expression of the computed time derivatives

M i = N i • Γ (9)
with

M i = ( xi , ỹi + x i z i + y i , zi -x i y i ), N i = [(y i -x i , x i , -z i ), (1, 1, 1)] and Γ = [(a, b, c), (α x , α y , α z )]. The linear regression consists in calcu- lating Γ = [(a , b , c ), (α x , α y , α z )] which minimizes ε = n i=1 (M i -N i • Γ ) 2 (10)
Hence, it translates into the condition

∂ε/∂Γ = - n i=1 2N T i • (M i -N i • Γ ) = 0 (11)
Further mathematical derivations lead to:

Γ = (N T • N ) -1 • (N T • M ) ( 12
)
where M and N are the matrices containing the considered data points, i.e. M = {M 1 , M 2 , ..., M n } and N = {N 1 , N 2 , ..., N n }. Therefore, the coefficients (a , b , c , α x , α y , α z ) better fitting the RC-predicted trajectories can be compared to the ones that were used for the numerical integration of the Lorenz system, namely (a, b, c, 0, 0, 0). Such comparison provides the error pred between the numerical integration of 1 (what we call the exact solution) and the prediction made by the AIS. The error is finally expressed as the euclidean distance between the two set of coefficients, i.e.:

pred = (a -a ) 2 + (b -b ) 2 + (c -c ) 2 + α 2 x + α 2 y + α 2 z ( 13
)
Note also that another way to measure the prediction error can be used, where the distance between the exact and the predicted trajectory are computed. However, for chaotic time series this is not a useful quantity, since any small error get amplified. Therefore, even though the systems are physically close to each other, their time series can be very disparate. For this reason, we better quantify the prediction error through Relation 13. Let us stress also that the prediction error is computed based on the coefficients (a , b , c , α x , α y , α z ) via a linear regression using 3n pred data points. Therefore, the prediction error may depend on the length of the prediction phase n pred .

It is worth to highlight that the PoPe method requires the prior knowledge of the system governing equations. For this reason, using the PoPe method for evaluating the prediction accuracy of the modelfree RC technique is not generally allowed. Nevertheless, the purpose of the present study is the parametric validation of the RC algorithm in reconstructing the Lorenz dynamics, whose equations are fully available. For this reason, the use of the PoPe method to evaluate the quality of the RC prediction is justified.

In this paper, we quantify the accuracy of the RC approach to predict Lorenz dynamics through a scan on the number of neurons D r and the spectral radius ρ of the Erdős-Rènyi network. In addition, we explore the impact of the length of the prediction phase on the prediction error, Equation 13. As already stated, the RC approach is based on the fact that all the parameters of the system are initialised randomly, except those of the output layer. This implies that the training phase is reduced to a simple linear regression. But this also implies that the results may depend on the random initialisation of the parameters. To overcome this difficulty, for each point in our scan (D r , ρ), we perform N = 500 realisations. Each of them will be different, since the initialisation of the network and of the input layer is random. For each realisation s we calculate pred , where the subscript rls stands for realisation.

The prediction error can be calculated in two different ways: (1) for time windows whose left endpoint is set to 0 and the right endpoint is upper-bounded by t pred = 25 and increasing for each time window as κδ, with δ = 1.25 and κ ∈ [1, 2, ..., n] denoting the number of the time window; (2) for moving time windows of the same length and whose endpoints are increasing as [κδ, (κ + 1)δ], with δ = 1.25 and κ ∈ [0, 1, 2, ..., n]. In other words, in the first way the error is computed over time windows whose lengths are increasing, considering therefore the dynamics from t pred = 0 up to κδ; and in the second way, the prediction error pred is calculated over moving time windows of fixed length δ = 1.25. In the following, we use the first method to compute the errors. We call this method increasing-time-window procedure. In this way the dynamics of the entire considered prediction phase is taken into account in calculating the deviation from the actual Lorenz trajectories. Nonetheless, it could be useful to compare these results with the moving average error, i.e. the second procedure described here above, as it is done for instance in section 3.3. The choice of the time window length of δ = 1.25 yields a sufficient number of points to calculate the linear regression and the subsequent predicted coefficients (a , b , c , α x , α y , α z ) in each time window of the moving error procedure. However, if not otherwise noted, the standard procedure for computing the errors is the increasing-time-window procedure.

Finally, it is also instructive to quantify the training error train , which is computed using the same Euclidean distance as in equation 13, but applied only to the training phase.

Results

The training and prediction errors averaged over all the realisations are plotted in figures 2(a) and (b), respectively, fixing the length of the prediction phase to t pred = 25, i.e. n pred = 1250. The isocontours of the averaged training and prediction errors are represented in logarithmic scale as a function of the number of neurons D r and the spectral radius ρ of the reservoir network. The error in the prediction phase is calculated for the entire prediction time, i.e. t pred = [0, 25]. It is observed that the training error is minimum in a rather narrow region in the (D r , ρ) plane. Moreover, mild variations of (D r , ρ) may lead to an abrupt increase of the training error. This means that there is a strong restriction in the allowed values of D r and ρ. It is particularly important to realise that increasing the number of neurons does not necessarily reduce the averaged training error. Moreover, we can observe that the averaged prediction error is significant, increased by several orders of magnitude with respect to the training error. The values selected in [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF] were D r = 300 and ρ ∈ {1.2, 1.45}, which fall into the beneficial narrow region where the training error is minimal, as the markers in figure 2(a) show. Yet, it is observed that the averaged prediction error is far from being negligible, considering the entire prediction phase t pred = 25. Within the set of scanned parameters in figure 2, a region where pred rls is minimized is noted for D r > 400 and ρ > 1.5. In the remainder of the paper, such a region, and the configuration included therein, will be called favorable. It is stressed again that, although the ESP originally required the spectral radius to be ρ < 1, the range of values of the spectral radius leading to optimal forecasting is obtained for ρ > 1, consistently with the results of Ref. [START_REF] Jiang | Model-free prediction of spatiotemporal dynamical systems with recurrent neural networks: Role of network spectral radius[END_REF]. Furthermore, the region where D r > 300 and ρ < 1.2 presents the largest measured averaged errors of the analyzed range of parameters in the prediction phase. For this reason, it will be termed as the unfavorable region of the (D r , ρ) plane. It must be noted that, although in a different framework from the present parametric validation study, very similar dependencies of the errors on the spectral radius have also been reported in a recent publication [START_REF] Verzelli | Learn to synchronize, synchronize to learn[END_REF]. It is also observed that the averaged errors in both training and prediction phases for each couple (D r , ρ) may be accompanied by a large variation within the set of performed realisations. In figure 3, the standard deviation σ of the averaged errors within the set of N = 500 realisations is plotted as a function of the number of neurons D r and the spectral radius ρ for the training and the prediction phases in (a) and (b), respectively. The set of realisations is the same displayed in figure 2, and the errors shown in panel (b) for the prediction phase are calculated for t pred = [0, 25], consistently with figure 2. The large dispersion that can be appreciated by measuring σ demonstrates the strong impact of the randomness on the prediction of Lorenz chaotic trajectories by the RC technique. The plot in figure 3(a) provides additional evidence that the possible efficient configurations of this AIS Complex Systems, Volume (year) 1-1+ based on the RC approach for the prediction of chaotic trajectories are limited to a narrow region in the (D r , ρ) plane. Furthermore, similarly to what has been observed in figure 2, the difference of the standard deviation for the set of realisations between the training phase and the prediction phase is significant, especially in the favorable region of the scan. Indeed, this proves the more frequent occurrence of large error events in the prediction phase, which are then responsible also for the increase of the averaged prediction error within the set of the realisations with respect to the training phase. Insightful information on the randomness of the RC prediction can be found by further analysing the two specific study-cases described in Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF]. From the statistically-relevant averaged results illustrated in figure 2, one can infer that the case with (D r , ρ) = (300, 1.45) performs better than the case with (D r , ρ) = (300, 1.2) (respectively, red diamond and red star in figure 2), as the pred rls is minor for the former than for the latter case. This result might not be in full agreement with what is reported in Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF], where the case with (D r , ρ) = (300, 1.45) is shown to be less performing to predict chaotic trajectories, with respect to the case with (D r , ρ) = (300, 1.2). This might be due to the random initialization of the input layer and the reservoir array. Indeed, for one single realization, the case with (D r , ρ) = (300, 1.2) can perform better than that with (D r , ρ) = (300, 1.45), consistently with Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF]. This is illustrated in figure 4, where the time traces of x, y and z are plotted for D r = 300 and the two values ρ = 1.2 (left column panels) and ρ = 1.45 (right column panels). In each panel, the realisation with minimum error is plotted by a red curve. Another randomly chosen realization is plotted by the blue curve. As a reference, the exact solution is given by the dashed black curve. It is observed that depending on the realization, the case ρ = 1.2 can perform better or worse than the case ρ = 1.45.

Statistical Analysis on the Error Distribution

Given the large statistical variation exhibited by the measured errors, we analyse the relevance of using the mean value as an indicator of the error for each single (D r , ρ) configuration within the set of realisations. The histograms of log 10 ( pred ) for four couples (D r , ρ) configurations are shown in figure 5. The same analysis on the error in the training phase log 10 ( train ) revealed analogous results. These four different cases are chosen to represent: (a) a favorable case, (b) and (c) the same cases studied by Pathak et al. in Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF] and (d) a case in the unfavorable region of the (D r , ρ) plane. In this particular analysis, the number of realisations has been increased up to N = 5000 to improve the statistical relevance of the data set, and the histograms are captured in n bins = 100 bins evenly spaced. As can be seen, the histograms are shifted to large error values going from panel (a) to panel (d). In the four panels, the Gaussian probability density function (PDF), whose mean and standard deviation are respectively equal to the statistical mean and standard deviation of the corresponding set of realisations, is over-plotted with black dashed lines. In this way, it is straightforward to compare the histograms and the Gaussian PDFs. Whereas for panels (b) and (c) the distribution of the error over the set of realisations follows a Gaussian, for panels (a) and (d) the histograms exhibit a departure from the normal PDF. In particular for the configuration (D r , ρ) = (550, 0.9), an inset with the rightmost tail of the distribution function is displayed in a log-log plot. In this way, it is possible to appreciate the deviation for the large-error events of such configuration from the tail of the Gaussian distribution. This deviation is present in all the analyzed configurations of the unfavorable region of the (D r , ρ) plane. Therefore, it is demonstrated the accuracy of using the statistical mean, i.e. the first moment of the distribution function, to quantify the error in both training and prediction phases for a significant range of analyzed configurations. Yet, the unfavorable region of the (D r , ρ) plane presents a non-negligible deviation from the Gaussian distribution, especially concerning the large error events. Additionally, it could be observed that the slight deviation of the histograms from the normal distribution function, especially in panel (d), can be explained by studying the third and fourth standardized moments of the PDF, i.e. the skewness µ 3 and the kurtosis µ 4 respectively. The skewness measures the loss of asymmetry of a distribution function. For a symmetric PDF the skewness is null (µ 3 = 0), hence a positive (negative) value indicates that the data-set is skewed to the right (left) with respect to the principal mode -the mean value for a unimodal Gaussian PDF. On the other hand, the kurtosis property measures the flatness of the PDF. Generally, the kurtosis is compared to the kurtosis of a Gaussian distribution function, which is equal to 3 (µ 4 = 3). Therefore, it is commonly used the excess kurtosis μ4 , which is the kurtosis re-scaled to 3 (μ 4 = µ 4 -3). A large excess kurtosis thus indicates that the distribution has strongly populated tails. Therefore, by inspecting the loss of asymmetry (skewness) and strength of the tails (excess kurtosis) together, it is possible to evaluate the behaviour of the measured PDFs. The measured distribution function in figure 5(a) and (d) reveal a positive skewness, denoting a major broadening towards large error values, as can be also inferred from the plots. Whereas for the configuration (D r , ρ) = (470, 2.2) the skewness is mildly above zero, the RC with (D r , ρ) = (550, 0.9) shows a four-time larger value of the skewness, which evidences the more frequent occurring of large error events. For this latter configuration the re-scaled fourth moment, the excess kurtosis, is positive and large indicating that the distribution of the data-set is strongly peaked around the mean value. In conclusion, the statistical properties of the cases in panels (a), (b) and (c) indicate that the distribution functions could be well approximated with a Gaussian around the corresponding mean values. On the other hand, for the unfavorable case in panel (d), the statistical mean value is strongly affected by a more frequent occurring of large-error events.

Dependence of the Error on the Prediction Phase Length

To shed some light on the validity of the RC to reconstruct the Lorenz dynamics, we additionally explore the dependence of the prediction error on the maximum length of the prediction phase. As illustrated in Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF], the AIS based on the RC approach correctly predicts the short-term trajectories, while significantly deviating from the actual Lorenz trajectories in the long-term phase. The considered length of the prediction phase has thereby a strong impact on the computed error. Therefore, we analyse the dependence of the averaged prediction error pred on the length of the prediction phase in the two different procedures that have been described at the end of section 2.2. The results are illustrated in figure 6. Hence, in panels (a)-(c) the prediction error is calculated for time windows whose period is [0, κδ], with δ = 1.25 and κ = [1, 2, ..., n], up to t pred = 25. Indeed, panels (d)-(f) represent the prediction error calculated in the second fashion, for which the period of the time windows is [κδ, (κ + 1)δ], with δ = 1.25 and κ = [0, 1, 2, ..., n]. The two procedures are labelled in the figure [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF] and eventually (c) and (f) cases from the unfavorable region of the (D r , ρ) plane explored in figure 2.

as increasing-time-window average error and moving average error, respectively.

As it has already been noted, even for the favorable cases in panel 6(a) in previous sections, the log 10 ( pred ) rls indicator is very large compared to the error in the corresponding training phase. It is also observed that there is a strong dependence on the length of the prediction phase in the short-term prediction. Focusing on the increasing-timewindow average procedure, after the prompt increase in the short-term prediction, the error remains constant (after an elbow) in the longterm prediction, consistently with Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF]. Such a behaviour is actually present for the favorable configurations also when the error is calculated with the moving average procedure, as panels (d) and (e) show.

Furthermore, smoother curves in the time evolution of the calculated error appear when the increasing-time-window average procedure is employed. This is to be understood as, in this latter procedure, the calculation has a kind of 'memory' from previous time windows, since the analyzed time period increases by including also the dynamics of the previous time windows. On the other hand, there is no such 'memory effect' in the moving average procedure, as the solution coefficients are calculated for each time window independently from the previous ones.

Another striking observation comes from the value of the error at the end of the first time period of the prediction phase (t pred ∈ [0, 1.25]). Indeed, the average prediction error pred rls already presents value much larger than the error evaluated for training phase for the same configuration. The difference between these two values can also be inferred from the plots in figure 2. Especially for the favorable configurations, the difference could be also more than 4 orders of magnitude. Therefore, the RC approach could show a strong disparity from the expected dynamics also for t pred < 2.5, i.e. for only 3 ∼ 4 oscillations of the Lorenz system. This could essentially be due to two reasons: either the reservoir layers are failed to be trained and the error is thereby propagated in the prediction phase inevitably, or the RC technique predicts strongly different dynamics in the short-term phase. The former condition, for which the training phase already produces a large error and thus the reservoir layers are not well-trained, is less frequent, as the lower mean values of the error in figure 2(a) and the standard deviations in figure 3(a) illustrate. Nevertheless, if a large error event occur in the training phase, this inevitably produces a large error event also in the prediction phase, since the RC is ill-trained and so unable to recover the correct Lorenz dynamics. The latter condition, for which the RC network is well-trained but fails in predicting the short-term dynamics, is indeed the most frequent. As a result, the average error in the first time window of the prediction phase is much higher than the one measured in the training phase.

Dependence of the Error on the Train Phase Length

It is now considered the effect of the length of the training phase on the error in predicting the Lorenz trajectories by means of the RC technique. In figure 7, the dependence of both train rls and pred rls on the length of the training phase is displayed for various chosen reservoir configurations. The plots illustrate the result for n train > 1000, as for smaller number of time steps the RC definitely fails in being trained and subsequently in predicting the chaotic dynamic of the system. The parameters of the various configurations are the same analyzed in panels (a) and (b) (and also (d) and (e)) of figure 6. The error in the prediction phase is considered for the entire prediction phase t pred = [0, 25]. A striking observation about the training phase for the configurations with (D r , ρ) = (430, 2.3) and (D r , ρ) = (470, 2.2) is the non-monotonic behaviour of the function. Indeed, for a narrow range of the training phase length, the train rls curves presents a minimum. Such a behaviour has been found also in other configurations within the favorable region of the (D r , ρ) plane. Instead, for the configuration setups with (D r , ρ) = (300, 1.2) and (D r , ρ) = (300, 1.45), the same analyzed by Pathak et al. in [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF], the error in the training phase is basically decreasing with the increasing number of considered time steps for training the system until a plateau is reached. This plateau could be expected due to the quasi-periodicity of the chaotic oscillations of the Lorenz system. Regarding the prediction phase, the behaviours of all the detailed configurations resemble. After a sharp decrease for n train > 1000, the curves mildly continue to decrease until n train = 4000 ∼ 5000, where they become averagely constant. Therefore, the increase of the number of time steps in the training phase beyond n train = 10000 leads to enhanced elapse time and complexity of the system, without improving the prediction of the chaotic trajectories significantly. Consistent results have been recovered also in Ref. [START_REF] Chattopadhyay | Data-driven predictions of a multiscale lorenz 96 chaotic system using machinelearning methods: Reservoir computing, artificial neural network, and long short-term memory network[END_REF], although a modified version of the Lorenz 96 chaotic model [START_REF] Thornes | On the use of scale-dependent precision in earth system modelling[END_REF] was employed as test-bed system instead of the original Lorenz system here studied. Furthermore, it is to be noted the jagged behaviour of pred rls in function of the number of time steps n train . Consecutive points of this particular scan could lead to very different results. Thus, increasing n train does not automatically lead to an improvement of the results, but it could be indeed detrimental in average.

Studies on the Time Evolution of the Predicted Lorenz Coefficients

Figure 8: The coefficients (ã, b, c, αx , αy , αz ) corresponding to the best solution of the RC application for the configurations also explored in Ref. [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF] are shown in the first two columns. The horizontal dashed lines represent the values of the relative coefficients in the actual Lorenz system. In the third column, the predicted dynamics of these two configurations are compared to the actual Lorenz dynamics (labeled 'Ground truth').

Additional analyses have been carried out on the evolution in time of the coefficients (ã, b, c, αx , αy , αz ) which better fit the prediction obtained by means of the RC technique. In the Lorenz system, such coefficients are constant throughout the whole time trajectory. This assumption, together with the constraint of having ( αx , αy , αz ) = (0, 0, 0), is crucial for the validity of the Lorenz system. Therefore, in order to ensure the reliability of the predicted dynamics of the chaotic system, the fitted coefficients have been plotted for two different RC configurations, which minimizes the prediction error at the end of the prediction phase pred (t pred = [0, 25]), in figure 8. These two configurations are the same analyzed by Pathak et al. in [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data[END_REF], i.e. (D r , ρ) = (300, 1.2) and (D r , ρ) = Complex Systems, Volume (year) 1-1+ (300, 1.45). As it has already been explained, (ã, b, c, αx , αy , αz ) are the coefficients which better fit the Lorenz equations for the RC predicted dynamics in the prediction phase. In this particular analysis, the coefficients are calculated for consecutive moving time windows of the prediction phase with a period of δ = 1.25. The same procedure has already been described in section 3.3 when the calculation of moving error has been applied to produce the figures 6(d)-(f). We repeat that this choice of δ = 1.25 enables to have sufficient elements in each time window for the calculation of such predicted coefficients by linear regression techniques and to grasp, thereby, their time dependence in the prediction phase. Thus, in the first two columns of figure 8, the coefficients of the solution are plotted against the time of the prediction phase t pred , while in the rightmost column the actual Lorenz dynamics are compared to each component of the AIS predicted systems. However, it is observed that in this optimized case the setup with (D r , ρ) = (300, 1.45) performs better than (D r , ρ) = (300, 1.2), consistently with what is displayed in figures 6(b) and 6(b) and (e). Indeed, the three predicted coefficients for (D r , ρ) = (300, 1.45) oscillates around the actual initial Lorenz coefficients (a, b, c) (highlighted in black horizontal dashed lines in the panels of the first column) throughout the whole prediction phase. The deviations from the exact coefficients of the studied Lorenz system appear almost negligible (±5%). Regarding indeed the case with (D r , ρ) = (300, 1.2), a significant deviation is noted in the time window around t pred = 11. Nevertheless, this deviation does not impair the whole prediction phase, as the actual coefficients are promptly recovered after the deviation occurring. If the third column of figure 8 is inspected for the configuration (D r , ρ) = (300, 1.2) (blue curves), it is possible to see that around t pred = 11 there is a pattern loss, which determines the high spike in the calculation of the predicted coefficients. Yet, such an instantaneous deviation from the exact Lorenz trajectory does not prejudice the overall dynamical properties of the predicted system.

However, as it has already been noted, the solution plotted in figure 8 are the predictions that minimizes pred (t pred = 25) within a statistically relevant set of N = 5000 realisations for each RC configuration. Therefore, the reader must be cautious on evaluating it, since in the same sample of realisations cases with large measured errors were frequently observed, as the average prediction error in figure 2 and the PDF of the prediction error in figure 5 demonstrate. The optimized cases reported in this section, indeed, are located in the tail of the corresponding PDF, and therefore the probability of their occurrence is quite low.

A more detailed analysis on the α coefficients is performed and the results are illustrated in figure 9. Here, the PDF of the α coefficients (α x , αy , αz ) for the configuration with (D r , ρ) = (300, 1.45) [11, 12.5], [23.5, 25]). The PDFs are binned in n bins = 500. The same analysis has been carried out also for other configurations in the favorable region of the (D r , ρ) plane, yielding similar results. What can be inferred here is that the distribution functions of the α coefficients are centered around zero, strongly broadening with increasing time in the prediction phase. Measuring the width of the PDFs, one can notice that it is significantly enhanced going from t pred = 1.25 to t pred = 12.5 for all the three coefficients, whereas no such a difference is measured from t pred = 12.5 to t pred = 25. To corroborate this result, the standard deviation σ of the α coefficients within the set of N = 5000 realisations is plotted in logarithmic scale against t pred . The green vertical dashed lines represent the upper boundary of the analyzed time windows. It is observed that already after the first time window, σ increases of more than two orders of magnitude. Then, it reaches a quasi-stationary phase with jagged behaviours. Therefore, as already proved by figure 9, the frequency of of large value events occurring is higher already after the first time window, but it does not increase additionally in the rest of the prediction phase. However, the standard deviation presents some spikes (especially for αy and αx ), which are symptomatic of a more frequent occurring of the large value events. These large value events, indeed, are correlated with the large error events in the prediction phase, as the relation 13 proves. To conclude, the analysis here performed shows that in the predicted Lorenz system the α coefficients, measured by the PoPe method, are centered around zero. However, quite frequently they differ from the null value, indicating that the actual Lorenz dynamics is definitely impaired. Moreover, the solutions displayed in figure 8 are only representative of the tail of the probability distribution, as the frequency of their occurrence is very low. Hence, even for favorable configurations, the accuracy of such a RC approach must be evaluated cautiously.

Conclusions

In this study, we have analyzed the validity of the reservoir computing ML technique applied to reconstruct the dynamical coefficients of the predicted trajectories of the Lorenz system, a well-known chaotic system, by the Projection on Proper elements method. Such a method is briefly explained in section 2 (more details can be found in Refs. [START_REF] Cartier-Michaud | Projection on proper elements for code control: Verification, numerical convergence, and reduced models. Application to plasma turbulence simulations[END_REF][START_REF] Cartier-Michaud | Verification of turbulent simulations using PoPe: quantifying model precision and numerical error with data mining of simulation output[END_REF]). The validity of the RC technique is statistically measured by running N = 500 realisations for a large set of configuration parameters. Thus, it has been shown that a good accuracy is achieved only for a small range of (D r , ρ) configurations of the reservoir. Yet, the error in the prediction phase is significantly increased with respect to the error in the training phase, with a significantly large variation of the error, as figure 3 shows. This implies that quite frequently the RC approach produces predictions strongly affected by large errors of the Lorenz system. It is shown that this is essentially due to the intrin-sic randomness of the RC technique, which inevitably leads to have large error events frequently also for the favorable configurations of the network. Additional scans on the relevant parameters of the ML technique such as the length of both training and prediction phases are carried out, showing that the range of validity of this RC approach is even further narrowed. Indeed, a minimum number of time steps in the training phase is required to achieve acceptable results, as figure 7 illustrates, but only for predicting the short-term dynamics of the Lorenz system. The deviation from the exact Lorenz dynamics in the long-term prediction phase is shown to become very large. Such results are also analyzed deeper by means of detailed studies on the time evolution of the Lorenz coefficients computed by the AIS. It is shown that the large variation of the Lorenz coefficients computed from the predicted solution can be significant. Therefore, the solution predicted by the AIS based on the RC approach is not necessarily representative of a Lorenz system. This paper, in the context of the RC approach to reconstruction of complex dynamical series, is helpful to establish the range of validity of this AIS technique. It also suggests that further developments of the RC paradigm, as recently proposed in, e.g. , Refs. [START_REF] Løkse | Training echo state networks with regularization through dimensionality reduction[END_REF][START_REF] Yeo | Data-driven reconstruction of nonlinear dynamics from sparse observation[END_REF], are required in order to robustly achieve a good accuracy in predicting chaotic time series.

  y, z; a) = a (y -x)ẏ = F y (x, y, z; b) = x (b -z) -y ż = F z (x, y, z; c) = xy -cz

Figure 1 :

 1 Figure 1: Schematic representation of the reservoir computing paradigm.

  (s) pred and we use the ensemble averaged quantity, pred rls = 1 N s (s)

Figure 2 :

 2 Figure 2: A scan of the error 13, in logarithmic scale, over the input parameters D r and ρ of the reservoir network is displayed for the training and the prediction phase, in panel (a) and (b) respectively. The errors are averaged over the whole set of realisations (N = 500) for each couple (D r , ρ). Highlighted with red markers are the configuration employed in Ref. [18].

3. 1 Figure 3 :

 13 Figure 3: The standard deviation σ of the logarithm in base 10 of the errors, for both training (a) and prediction (b) phases, within the set of N = 500 realisations is displayed for the same scans reported in figure 2. The large variation within the whole set of realisations can be thus appreciated.

Figure 4 :

 4 Figure 4: The dynamics of the three components of the Lorenz system, predicted by the RC, is shown for the two configurations studied in Ref. [18], i.e. (D r , ρ) = (300, 1.2) (left column) and (D r , ρ) = (300, 1.45) (right column). The black dotted curves are the actual chaotic trajectories obtained through numerical integration of the Lorenz equations with the RK4 method (labelled as Ground truth in the legend). Red and blue curves show respectively the best and a non-favorable solution obtained through RC technique over the whole set of N = 500 realisations.

Figure 5 :

 5 Figure 5: The error in logarithmic scale in the prediction phase, calculated at t pred = 25, of an incremented set of N = 5000 realisations is shown in histogram plots for four different RC configurations. The distribution is binned in n bins = 100 samples. The black dashed curves represent the unimodal normal PDF, whose mean value and standard deviation are the mean value and the standard deviation of the error distribution functions for the corresponding configuration.

Figure 6 :

 6 Figure 6: The error in the prediction phase is displayed as a function of the length of the prediction phase for different configurations. The prediction error is calculated at increasing width of the time windows in panels (a)-(c), whereas with moving time window of fixed period δ = 1.25 in panels (d)-(f). The configurations reported are picked from the whole scanned cases, with the purpose of representing: (a) and (d) favorable cases, (b) and (e) the same cases reported in Ref. [18] and eventually (c) and (f) cases from the unfavorable region of the (D r , ρ) plane explored in figure 2.

Figure 7 :

 7 Figure 7: The errors in the training (a) and in the prediction phase (b) are shown as a function of the length of the training phase for four different input configurations of the reservoir network. The configurations are the same analyzed also in figures 5 and 6. The abscissas is expressed both in length of the training phase in normalized time t train and in number of time steps n train in the training phase.

Figure 9 :

 9 Figure 9: The probability distribution function of the coefficients (α x , αy , αz ) for the RC configuration with (D r , ρ) = (300, 1.45) is shown in histograms for three different time windows of an incremented set of N = 5000 realisations. The distribution is binned in n bins = 500 samples. The standard deviation of the PDFs is shown in figure 10 for the entire prediction phase.

Figure 10 :

 10 Figure 10: The evolution in time of the standard deviation σ of the α coefficient PDFs diplayed in figure 9 is shown in logarithmic scale for the prediction phase up to t pred = 25. The green dashed vertical lines represent the time windows displayed in figure 9.

Complex Systems, Volume (year) 1-1+; year Complex Systems Publications, Inc.

Complex Systems, Volume (year) 1-1+

Acknowledgements

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under Grant agreement No 633053. The views and opinions express herein do not necessarily reflect those of the European Commission.

Data Availability Statement

The data and the algorithms that support the findings of this study are available from the corresponding author upon request.

Conflict of interest

The authors have no conflicts to disclose.