Inference of Unknown Dynamical Components Using Next Generation Reservoir Computing: From Chaotic Systems to Climate Data Jule Budnicka , Andrew Keanea , Serhiy Yanchuka,b,∗ a School of Mathematical Sciences, University College Cork, Cork, Ireland
arXiv:2609.24754v1 [cs.LG] 21 Sep 2026
b Potsdam Institute for Climate Impact Research, Potsdam, Germany
Abstract We investigate next generation reservoir computing (NGRC) as a data-driven approach for inferring unseen components of dynamical systems. We compare NGRC with traditional reservoir computing (RC) using the Lorenz and Rössler system, where two unknown components are inferred from one given component. For both systems, NGRC achieves accurate results while requiring fewer training data and less computational time than RC. We identified an inverse proportional behavior between the number of time-delayed steps needed for NGRC and the temporal resolution, indicating that the physical time span covered by the delay interval is an important factor in determining the required number of delayed steps. Finally, we apply NGRC to the observational climate data of ENSO (El Niño–Southern Oscillation) and infer one observable from the remaining variables. Despite the noise and complexity of the real-world data, the NGRC shows promising results. Our findings demonstrate the potential of NGRC for efficient inference of unseen components in both controlled dynamical systems and real-world data. Keywords: NGRC, RC, Unseen components, Chaotic dynamical systems, ENSO
1. Introduction Data-driven approaches, especially machine learning (ML) algorithms, have become the state-of-the-art framework to predict and reproduce chaotic motions of dynamical systems in their various application areas [1–5]. A special ML paradigm is reservoir computing (RC), whose time-series data may be successfully analyzed through the lens of dynamical systems theory [6–8]. A reservoir computer is a recurrent neural network (RNN) characterized by interconnected neurons that are arranged in a so-called reservoir (a pool of nodes) instead of layers. In contrast to standard RNN algorithms, only the output weights of the RC are trained with a simple and efficient least squares method. All internal network parameters, such as the strengths of the node-to-node connections or the input coefficients, are chosen randomly. With these simplifications, the RC algorithm enables shorter training times with smaller training data sets, compared to standard RNNs, even for high-dimensional time-dependent tasks [9]. Despite its simplicity, the RC performs as well as other ML techniques for certain tasks [10]. It is widely used for predicting dynamical systems output [11, 12], learning climate data [13, 14], and inferring unseen variables [15]. Furthermore, a class of RCs adapted to nonlinear systems with delayed feedback, known as delay-based RCs, have shown excellent performance with efficient and fast computation [16–22]. RC is closely related to the vector autoregressive (VAR) method for modeling and forecasting multiple interrelated time series. In the specific case of a linear RC (i.e., the activation function is linear), [23] shows the equivalence of a RC with linear readout to a VAR as well as of an quadratic readout reservoir to a ∗ Corresponding author
Email addresses: [email protected] (Jule Budnick), [email protected] (Andrew Keane), [email protected] (Serhiy Yanchuk)
nonlinear VAR. This may explain the surprising success and efficiency of RCs [24]. Further work has also shown that under certain conditions a RC is a universal approximator [25]. When implementing the RC approach an optimization of hyperparameters is necessary. These hyperparameters characterize the general properties of a RC and strongly influence its performance. However, the corresponding parameter space can be large, including the number of reservoir nodes, the spectral radius of the network, the sparsity or the input amplitude [26]. There are some explicitly known results, for example, in relation to the spectral radius of the network and the need for recurrent connections [6, 26]. Additionally, there are more general results regarding the sparsity or low-connectivity of the reservoir [12, 27], or the leaking rate [26, 28]. Nonetheless, in practice, an extensive parameter search is still required. Possible optimization algorithms are a grid search [29], Bayesian [30, 31] or gradient descent methods [28], all of which are computationally expensive or only work on a continuous parameter set. As a further development of RC, next generation reservoir computing (NGRC) [32] greatly simplifies the optimization procedure. The NGRC algorithm depends on only two hyperparameters: the number of time delayed steps and a regression parameter. At each time step i, the network receives the input data Xj from the previous k time (delayed) steps (i ≥ j > i−k). These delayed inputs constitute the feature vector, so the choice of k determines both its size and the associated computational expense. The regression parameter results from a least squares method used to train the algorithm. Recent studies have demonstrated the advantages of NGRC over traditional RC when only limited training data are available [33]. According to [32], there are three common benchmark tasks which a ML algorithm should accomplish for a dynamical system: (i) short-term forecasting, (ii) reconstructing the attractor (long-term forecasting), and (iii) inferring unseen data. We investigate the performance of NGRC and compare it to the standard RC for the inference of unknown dynamical components of chaotic systems and climate data. Inference tasks provide a means to recover unobserved variables of a dynamical system from a subset of variables that can be measured. This can be achieved if a dynamical system exhibits sufficiently strong internal couplings between variables and the observed variables carry enough information to reconstruct the unseen components. Mathematical and statistical techniques including state-space reconstruction or projections into latent spaces using autoencoders or proper orthogonal decomposition, allow applications such as inferring neural activity from partial electrode recordings [34] or reconstructing surface pressure fields from sparse measurements [35]. The same logic is highly valuable in climate science, where observational constraints vary widely across variables. Identifying which climate variables can be reliably inferred from more accessible measurements reduces the burden on observing networks and increases resilience to missing or intermittent data [36–38]. In this paper we investigate how the NGRC performs against the traditional RC for the task of inferring unseen variables/observables. We consider the Lorenz [39] and Rössler [40] chaotic systems, as well as climate data measurements related to the El Niño–Southern Oscillations (ENSO) system [41, 42]. Table 1 gives an overview of some known applications of NGRC and RC to these three systems, as well as new contributions from this study. In the following two sections we provide a brief summary of the NGRC methodology and the three dynamical systems we consider. In Section 4 we demonstrate the applicability of NGRC to the inference of unknown components of chaotic dynamical systems in the previously untested setups: two unknown components of the Lorenz system, two unknown components of the Rössler system, and one unknown data series from ENSO measurements. The latter shows the usefulness of NGRC for complex and noisy real-world data. We find that neither algorithm consistently outperforms the other regarding different metrics, however, the NGRC has the advantage of requiring less training data and having a significant reduced computational cost. Both approaches perform surprisingly well with the noisy, real-world ENSO data by capturing prominent trends. However, it is also clear that there is much room for improvement. Additionally, we demonstrate how the optimal choice of the number of time delayed steps depends on the temporal resolution of the data and the desired accuracy of inference.
2
Previous results
This study
Inference of one component for Lorenz system [32]
Inference of two components for Lorenz system
–
Inference of two components for Rössler system
–
Inference of one component for ENSO
Inference of two components for Lorenz system [15]
Comparing NGRC to RC for Lorenz system
Inference of two components for Rössler system [15]
Comparing NGRC to RC for Rössler system
–
Inference of one component for ENSO
NGRC
RC
Comparing NGRC to RC for ENSO Table 1: A summary of known previous results on inference of unknown components and new results from this work. 2. Next Generation Reservoir Computing The architecture of a traditional RC is described in [7] or [11]. The NGRC framework was first described by [32] and we now present the main ideas of their approach. The first step is the construction of a feature vector directly from the input data. Definition 1. (Feature vector) For a time discretization ti , i = 1, .., n, let Xi ∈ RNin be the input data at time ti . The feature vector of the NGRC is given by (p)
Ototal,i := c ⊕ Olin,i ⊕ Ononlin,i ∈ RNtotal ,
(1)
(p)
where c ∈ R is a constant, Olin,i represents the linear and Ononlin,i the nonlinear part, and ⊕ denotes vector concatenation. The linear part of the feature vector is given by Olin,i := Xi ⊕ Xi−1 ⊕ · · · ⊕ Xi−(k−1) ∈ RkNin , where k < i is the number of time delayed steps for which the data Xj , i − k < j ≤ i is used. The nonlinear part of the feature vector is a function of the linear part defined as (p)
Ononlin,i := Olin,i ⌈⊗⌉ Olin,i ⌈⊗⌉ . . . ⌈⊗⌉ Olin,i , where the operator ⌈⊗⌉ denotes the calculation of the outer product, then takes the upper triangular matrix (the unique monomials) and flattens them into a vector. The product is taken p-times, hence, the parameter p describes the maximal polynomial power appearing in the feature vector. The total dimension Ntotal depends on the input dimension Nin , the time delayed steps k, the maximal power p and is equal to Nin k+p . p For example, the feature vector for a system with Nin = 3, described by X = [x, y, z], at time step ti 3
with k = 2 and p = 2 is given by Ototal,i = [ci , xi , yi , zi , x2i , xi yi , xi zi , yi2 , yi zi , zi2 , . . . 2 2 xi−1 , yi−1 , zi−1 , x2i−1 , xi−1 yi−1 , xi−1 zi−1 , yi−1 , yi−1 zi−1 , zi−1 ,...,
xi xi−1 , xi yi−1 , xi zi−1 , yi xi−1 , yi yi−1 , yi zi−1 , zi xi−1 , zi yi−1 , zi zi−1 ]. The size of the feature vector will grow very quickly with larger choices of k and p. [32] has shown that feature vectors containing only low-order monomials may be sufficient to perform accurate predictions. Therefore, we choose p = 2 for this study. The training of a NGRC focuses on the accurate prediction from the state Xi to Xi+1 . This is described in the following definition. Definition 2. (NGRC setup and training) Let Xi ∈ RNin be the input data at time ti and let Ototal,i Xi , Xi−1 , . . . , Xi−(k−1) ∈ RNtotal the feature vector according to Definition 1. The NGRC model is given by for i = 1, .., n, (2) Xi+1 = Xi + Wout Ototal,i Xi , Xi−1 , . . . , Xi−(k−1) where the weight matrix Wout is trained analogously to the RC approach. In matrix notation, Equation (2) leads to ∆X = Wout Ototal ,
(3)
where ∆X = (X2 − X1 , . . . , Xn+1 − Xn ) ∈ RNin ×n , Wout ∈ RNin ×Ntotal and Ototal = (Ototal,1 , . . . , Ototal,n ) ∈ RNtotal ×n . Therefore, Tikhonov regularization (also known as ridge regression) can be applied to find Wout 2 2 f f Wout = argminW fout ∥∆Xd − Wout Ototal ∥F + α∥Wout ∥F −1 = ∆Xd OTtotal Ototal OTtotal + αI , (4) where ∆Xd ∈ RNin ×n represents the desired difference between two states of the system at consecutive time steps (given by the training data set), α is the regularization parameter and v uX m u n X aij for A ∈ Rm×n ∥A∥F := t i=1 j=1
the Frobenius-Norm of a matrix. For the inference of unknown components, an adjustment of the NGRC algorithm is necessary. Consider a Nsys -dimensional system. We now want to predict Ninfer unknown components from the other Nin given components, such that Nsys = Nin +Ninfer . Therefore, the flow equation (2) needs the following modification: Yi+1 = Wout Ototal,i , where Yi+1 ∈ RNinfer represents the vector of the Ninfer components to predict. For p = 2, the nonlinear part of the feature vector at time ti is given by (p)
(2)
Ononlin,i = Ononlin,i = Olin,i ⌈⊗⌉ Olin,i . Hence, the outer product is a symmetric (kNin × kNin ) matrix and the dimension of the nonlinear part is
4
kNin (kNin +1) . Therefore, Ototal,i has 2
1 + kNin +
kNin (kNin + 1) 2
components. It is useful to mention the following two main advantages of NGRC over the RC, as discussed in [32]. First, the randomly chosen matrices in the RC approach are replaced by a fixed construction of the feature vectors. Therefore, the number of possible hyperparameters is reduced. This significantly reduces the number of required training data points. Secondly, there are fewer warm-up points (i.e., the points which are necessary to initialize the algorithm) and training points needed in comparison to RC. For a traditional RC, the warm-up period can last from 103 to 105 data points [12–14] since longer warm-up times are needed to guarantee that the RC becomes independent of the RC initial conditions. The necessary number of warm-up points in the NGRC algorithm is merely the number of time delayed steps k which are needed to create the feature vector at the first time step. 3. Three systems used in this study The three systems considered are the Lorenz system, the Rössler system and as a climate data application the El Niño–Southern Oscillation (ENSO). The Lorenz and Rössler systems are three-dimensional models that exhibit chaos and contain a strange attractor for a suitable choice of parameters. The datasets used to train the NGRCs are generated by numerical integration of these models using the fourth-order Runge-Kutta method. The application to the ENSO phenomenon, on the other hand, is based on observational data only. In the following we describe each system. 3.1. Lorenz system The fundamental work on chaotic behavior goes back to Edward Lorenz who developed a system to model a simplification of atmospheric convection [39]. This rather simple system exhibits chaos and contains a strange attractor. The three-dimensional Lorenz system is given by ẋ(t) = σ (y(t) − x(t)) , ẏ(t) = rx(t) − y(t) − x(t)z(t),
(5)
ż(t) = x(t)y(t) − bz(t), with (x(t), y(t), z(t)) ∈ R3 and parameters σ, r, b ∈ R [39]. We will consider the original parameter choice of Lorenz: σ = 10, b = 38 , and r = 28, for which the system exhibits a strange attractor [43, 44]. Fig. 1(a) shows a trajectory of the Lorenz system, in terms of time t and the phase space, calculated by numerically integrating Eqs. (5) and removing transients. There are already several studies of short- and long-term forecasting for the Lorenz system. For example, a forecasting task has been performed by [10] or [12] for the traditional RC. Similarly, [32] showed forecasting of the Lorenz attractor using NGRC and achieved accurate short-term predictions and nearly exact reconstruction of the attractor. For the inference of unknown components, [15] showed accurate predictions of two unseen components with RC. For the NGRC algorithm, [32] performed the inference of the z component from both the x and y components. In this work, we perform an inference of two unknown components using NGRC. 3.2. Rössler system The Rössler system is given by ẋ(t) = − (y(t) + z(t)) , ẏ(t) = x(t) + ay(t), ż(t) = b + z(t) (x(t) − c) , 5
(6)
Figure 1: a) A chaotic trajectory of the Lorenz system for σ = 10, b = 83 , and r = 28. b) A chaotic trajectory of the Rössler system (6) for a = 0.2, b = 0.2 and c = 5.7. c) ENSO data sets plotted against time (in months) for sea surface temperature (SST), southern oscillation index (SOI), precipitation index (ESPI) and zonal winds. with (x(t), y(t), z(t)) ∈ R3 and a, b, c ∈ R. We consider the parameter values a = 0.2, b = 0.2 and c = 5.7 as in [40]. Fig. 1(b) shows an example of a chaotic trajectory as time series and also in the (x, y, z)-space. For the Rössler system, [10] demonstrates long- and short-term forecasting of a chaotic trajectory. The inference of two unknown components was performed in [15] using RC. However, no application of the NGRC to the inference of unknown components exists for the Rössler system (6). 3.3. El Niño–Southern Oscillation data series The El Niño–Southern Oscillation (ENSO) describes a climate phenomenon in the equatorial Pacific. ENSO is one of the most important “weather-makers” with worldwide changes in temperature and rainfall [42]. An extreme ENSO can affect agriculture [45], ecosystems [46, 47], power generation [48] and extreme weather events worldwide such as wildfires, floods and tropical cyclones [49–51]. Therefore, considering the massive impact of the ENSO on the global climate and economy, an ability to understand and precisely forecast its pattern is crucial. ENSO constitutes two interacting components in the ocean and atmosphere of the equatorial Pacific. El Niño events are associated with large sea-surface temperatures (relative to the long-term mean) and the Southern Oscillations that describe fluctuations in the surface air pressure. Additionally, other weather components play a role, such as winds and precipitation [42]. 6
Name
Description
Time Coverage
Dataset
Source
SST
Sea surface temperature anomalies in the Niño 3.4 region (5°N to 5°S; 170°W to 120°W) [52]
monthly, since 1950
HadISST1.1
NOAA CPC
SOI
Southern Oscillation Index. Atmospheric components of ENSO. Sea level pressure differences between Darwin and Tahiti, normalized [53]
monthly, since 1951
CRU
NOAA PSL
ESPI
Precipitation anomalies in two areas: eastern tropical Pacific (10°S to 10°N, 160°E to 100°W) and Maritime Continent (10°S to 10°N, 90°E to 150°E), normalized [54]
monthly, since 1979
GPCP
NOAA PSL
200mb Zonal Winds
Zonal Winds equator anomalies (2.5°S to 2.5°N; 165°W to 110W°)
monthly, since 1979
NCEP Reanalysis
NOAA PSL
Table 2: Description of the time series data with corresponding sources. The data sources are either NOAA (National Oceanic and Atmospheric Administration) PSL (Physical Sciences Laboratory) or CPC (Climate Prediction Center). For our machine learning task, we consider four different data sets described in Table 2: sea-surface temperature anomolies (SST), southern oscillation index (SOI), ENSO Precipitation Index (ESPI) and zonal wind anomalies. Fig. 1(c) shows time series of these data sets where time t is in terms of months. The maxima of SST [55] and ESPI [54] as well as the minima of SOI [56] and winds [57] (above or below a certain threshold) correspond to an El Niño event. The real-world ENSO system is clearly higher dimensional than the three-dimensional Lorenz and Rössler systems. While progress has been made in the use of RC for forecasting El Niño events (for example, see [58, 59]), we do not focus on forecasting here. Instead, we consider the task of using only some observed data sets to predict another “unseen” observable of the ENSO system. For example, can the ESPI be accurately predicted using only the SST, SOI and wind data sets? 4. Inference of unknown components We focus on the inference of two unseen data components for the Lorenz and the Rössler chaotic systems and one component for ENSO. In each case we find the optimal performance of the NGRC algorithm, which depends on two hyperparameters: the number of time delayed data points fed into the network and the ridge regression parameter. Note that we fix the nonlinearity degree to p = 2. We examine the influence of the hyperparameters on the performance of the NGRC as well as the influence of the discretization in time. For this purpose, the NGRC is optimized with respect to its performance in the testing phase. We compare the performance of the NGRC with traditional RC regarding its accuracy, the number of training points and the training time. We use the same data sets for RC and NGRC. The data set is divided into a training and a testing part. The training data for RC and NGRC may differ since RC needs, in general, more training points for a good performance [32]. To measure the performance, the normalized root mean-squared error (NRMSE) of the resulting predicted time series and the normalized maximal distance 7
is calculated via r NRMSE :=
dmax :=
1 Pn ||X̂i − Xi ||2 n i=1 , ||Xmax − Xmin ||
maxi=1,...,n ||x̂i − xi || xmax − xmin
(7) (8)
where X̂i is the (one- or two-dimensional) predicted state, Xi is the observed state, n is the length of the time series. The maximal distance is calculated for each component of the predicted state.
Figure 2: NGRC applied to Lorenz system (5). (a) The NRMSE of the testing phase plotted on a semilogarithmic scale against parameters α ∈ [10−9 , 0.1] and k ∈ [1, 30]. The plotted NRMSE is the median over 100 randomly chosen initial conditions for each parameter pair (α, k). The dashed lines indicate the optimal values for α and k. (b) Lorenz attractor with optimal NGRC performance for α = 0.001 and k = 15, ∆t = 0.05 and initial condition (x0 , y0 , z0 ) = (5.7, −1.5, 3.1). The blue curve shows the target data, the dashed orange curve the prediction. The corresponding testing error is NRMSE ≈ 3.92 × 10−3 . (c) NRMSE for the training and testing phases for different time delayed steps k with fixed α = 0.001 plotted on a logarithmic scale. The dashed line shows a decrease NRMSE ∼ k −1 . Parameters and initial conditions: ∆t = 0.05, (x0 , y0 , z0 ) = (−12, 5, 22). (d) Distribution of the optimal k values found for 200 random initial conditions, and the median of the NRMSE for each optimal k. The bars indicate the first and third quartile.
8
4.1. Inference of two unknown components of the Lorenz system with NGRC 4.1.1. Performance and hyperparameter optimization In the following we analyze the performance of the NGRC at inferring the y and z components of the Lorenz system (5), given its x component. First, we find the optimal hyperparameters k (time delayed steps), α (ridge regression parameter) and ∆t (timestep). For this purpose, 50 warm-up points, 400 training points, and 800 testing points are selected successively from a given (numerically simulated) data set, in a similar fashion to [32]. A grid search is performed to find the optimal k and α, for which the NRMSE of the testing phase is minimal. The NRMSE is calculated in Fig. 2(a) for the 30 × 100 grid of parameters k ∈ [1, 30], α ∈ [10−9 , 0.1] and fixed ∆t = 0.05. For each parameter pair, k and α, the NRMSE is calculated for 100 randomly chosen initial conditions ((x0 , y0 ) ∈ [−20, 20]2 , z0 ∈ [0, 50]) with the plotted value representing the median. We observe a relatively large NRMSE for small values of α, as well as for large and small values of k. The minimum is located at k = 15 and α = 0.001, as indicated by the two dashed lines in panel (a). Panel (b) demonstrates the inference of variables y and z based on x, using the initial condition (x0 , y0 , z0 ) = (−12, 5, 22). The inference does appear very successful, since the prediction of y and z in orange overlaps the target data in blue. Figure 2(c) shows the NRMSE for the training and testing phases as a function of the time delayed steps k with fixed α = 0.001. For comparison, the dashed line shows a decrease NRMSE ∼ k −1 . The minimum of the NRMSE for the testing phase corresponds to the optimal k = 15. Overfitting can be observed for larger k, as the training error decreases while the testing error increases. To demonstrate the dependence on initial conditions, panel (d) shows the distribution of optimal k values found for 200 randomly chosen initial conditions and a range of k ∈ [1, 20]. This distribution peaks between 14 and 17. The optimal NRMSE lies in most of the cases near 3 × 10−3 showing accurate prediction for different initial conditions. Note that we take the median instead of the mean since the 200 initial conditions are chosen randomly and the algorithm fails in some special cases, where the training data is not representative of the chaotic attractor. We now examine the influence of the timestep discretization, ∆t, on the optimal value of time delayed steps k. For this, we consider 20 equally distributed values ∆t ∈ [0.01, 0.1], 30 values k ∈ [1, 30], and 50 different initial conditions for each ∆t. Figure 3(a) displays the optimal k for each ∆t, while panel (b) shows the corresponding NRMSE in terms of the median and the first and third quartile. The horizontal axis is cut off for small values of ∆t in panel (b), since we set an upper limit of k = 30, resulting in large MSE values in panel (b) and the plateau in panel (a) for small time steps. In panel (a) we observe a clear trend for ∆t > 0.02, indicating
Figure 3: NGRC applied to Lorenz attractor. (a) The optimal number of time-delayed steps k and (b) the corresponding NRMSE as a function of the discretization timestep ∆t. For each ∆t, the median and the first and third quartiles are calculated over 50 randomly fixed initial conditions each time step. The orange curve shows k∆t = 0.77. 9
how the optimal time-delayed steps is inversely proportional to the timestep size. By assuming an inverse proportional dependence, we obtain the mean value with standard deviation of k∆t as (9)
k∆t = 0.77 ± 0.06,
plotted in orange in Fig. 3 (a). This value represents the amount of memory that the NGRC needs in order to correctly infer the unknown components. 4.1.2. Comparison of NGRC and RC First, we describe the training and testing setups for NGRC and RC. The same 800 data points were used for testing the RC and NGRC, however the training sets are different. Since hyperparameter optimization for RC has already been performed by [15], where the y and z components were predicted for a given x input, we use the hyperparameters for RC from [15]. For the ridge regression parameter, an optimal value of 10−9 was estimated by considering 5 different initial conditions for 15 regression values and 30 trials each. The RC training was conducted using 5200 [15] data points to facilitate accurate training, and also using 400 training points for a more direct comparison with the NGRC. Table 3 summarizes the results. The median and the 95% confidence interval are calculated over 100 initial conditions. The computation time is the training and prediction time using data from a single initial condition 100 times, then divided by 100. This is done 7 times to calculate a mean and standard deviation. Generally, the RC with 5200 training points performs slightly better than the NGRC regarding both the NRMSE and the maximal distance. It should be noted that the training errors are less comparable since the training sets are different. However, the NGRC requires only 1/13 of the training points compared to the RC algorithm resulting in much larger computational times for the RC, even considering the same training set of 400 points. NRMSEtest ×10−4
NRMSEtrain ×10−4
ymax ×10−3
zmax ×10−3
comp. time 17 ±0.1ms
NGRC [400]
median
27
17
22
18
95%
[21, 99]
[10, 22]
[12, 78]
[4, 62]
RC [5200]
median
2.8
2.2
0.2
4.0
95%
[1.2, 17]
[1.0, 5.2]
[0.06, 0.5]
[1.2, 37]
RC [400]
median
13.5
0.5
1.1
26
95%
[3.8, 166]
[0.2, 3.5]
[0.2, 6.5]
[6.4, 186]
635 ±17ms 284 ±11ms
Table 3: Comparison of NRMSE, normalized maximal distance (8) and computation time for the Lorenz system using the NGRC with 400 training points and RC with 400 and 5200 training points. The number of testing points is 800 in all cases. The medians and the 95% confidence intervals are calculated over 100 different initial conditions. The average computation time is calculated for the initial condition (x0 , y0 , z0 ) = (−12, 5, 22). Timestep is ∆t = 0.05.
4.2. Inference of two unknown components of the Rössler system with NGRC 4.2.1. Performance and hyperparameter optimization We apply a similar framework for the analysis of the NGRC performance, with respect to the parameters k, α and ∆t, for the Rössler system (6) as we did for the Lorenz system. We choose 50 warm-up points, 600 training points and 1000 testing points successively from a given dataset. We consider a given x component and use the trained NGRC to infer y and z components. A grid search is performed to find the optimal α and k, where the time step size is fixed as ∆t = 0.1. The NRMSE is calculated in Fig. 4(a) for a 100 × 30 grid of parameters α ∈ [10−8 , 0.1] and k ∈ [1, 30]. 10
For each pair of parameters, the median of the NRMSE is shown for 100 randomly chosen initial conditions ((x0 , y0 , z0 ) ∈ [−1, 1]3 ). In contrast to the Lorenz system, the NRMSE decreases for further increasing k without an overfitting effect, at least until k = 30. In other words, there is no interior minimum of NRMSE in Fig. 4(a), as there is in Fig. 2(a), instead there is only a boundary minimum at k = 30. Hence, instead of selecting 30 as the optimal k value, we select the (α, k) pair that satisfy a certain threshold NRMSEtest ≤ ε. while minimising k. We call this value ε the stopping threshold. For Fig. 4, we choose a stopping threshold of ε = 0.001 which is satisfactory for the NRMSE in our setting. Under this assumption, the optimal hyperparameter setting is α = 10−7 and k = 15, as indicated by the dashed lines, which provides the minimal value of k with a NRMSE equal to or less than 0.001. Panel (b) illustrates the performance of the NGRC with the optimal hyperparameter setting and the initial condition (x0 , y0 , z0 ) = (−0.6, 0.2, −0.1). Again, the prediction of the inferred variables in orange overlap with the target data, indicating that 0.001
Figure 4: (a) NGRC applied to Rössler system (6) with stopping threshold ε = 0.001. The NRMSE of the testing phase is plotted on a semi-logarithmic scale against the parameters α ∈ [10−8 , 0.1] and k ∈ [1, 30]. The plotted NRMSE is the median over 100 randomly chosen initial conditions for each parameter pair (α, k). (b) Inference of y and z with optimal NGRC performance for α = 10−7 and k = 15, ∆t = 0.1, and initial condition (x0 , y0 , z0 ) = (−0.6, 0.2, −0.1). The target data is shown in blue, the prediction with dashed orange. The testing error is NRMSE ≈ 5.9 × 10−4 . (c) NRMSE for the training and testing phase plotted on a logarithmic scale versus the time delayed steps k ∈ [1, 60] with α = 10−7 fixed. The dashed line shows NRMSE ∼ k −2 . Parameters and initial condition: (x0 , y0 , z0 ) = (−0.6, 0.2, −0.1), ∆t = 0.1. (d) Distribution of the optimal k values found for 200 random initial conditions for two different stopping thresholds ε = 0.005 and ε = 0.001, α = 10−7 . 11
is a sufficient stopping threshold. Fig. 4(c) shows the influence of varying the number of time delayed steps, k ∈ [1, 60], on the NRMSE. For comparison, the dashed line represents NRMSE ∼ k −2 . Additionally, panel (d) shows the distribution of optimal k values for α = 10−7 and two different stopping thresholds ε = 0.005 and ε = 0.001. An optimal k is found for 200 randomly chosen initial conditions. The optimal k peaks at k = 6 for the larger threshold and lies mostly between 12 and 15 for the smaller one. This suggests that, since the size of the feature vector scales quadratically with k, the computational effort can be significantly reduced by slightly lowering the stopping threshold.
Figure 5: NGRC applied to Rössler attractor depicting the influence of the step size ∆t on the optimal k. (a) Stopping threshold (a) ε = 0.005 and (b) ε = 0.001. For each ∆t, the median is calculated over 50 randomly chosen initial conditions, and the bars indicate the first and third quartile. The orange curve shows k∆t = 0.66 for ∆t ∈ [0.05, 0.1] in (a) and k∆t = 1.44 in (b). To examine the influence of ∆t on the optimal k, we take 20 equidistantly distributed values for ∆t ∈ [0.05, 0.15], 30 values for k ∈ [1, 30] and 50 initial conditions for each time step. Fig. 5 shows the results for two different stopping thresholds, ε = 0.005 in panel (a) and 0.001 in panel (b). We find the inverse proportional dependence k∆t = 1.44 ± 0.11 for the stopping threshold ε = 0.001. Increasing the stopping threshold to ε = 0.005 reduces the value to k∆t = 0.66 ± 0.10 for ∆t ∈ [0.05, 0.1]. Hence, fewer delayed time steps k are sufficient when a lower prediction accuracy is tolerated. 4.2.2. Comparison of NGRC and RC Analogously to the Lorenz system, the data points are chosen such that the testing phases of both algorithms take the same points to enable a comparison between the two algorithms. However, the training sets may differ. We take the parameter values from [15] and optimize the RC with respect to the regression parameter by choosing the same stopping threshold and considering 5 different initial conditions with 30 trails each. The training was conducted using 2600 training points, and a further comparison to NGRC was made by using 600 training points. Table 4 summarizes the results. The NRMSE for both algorithms (NGRC with 600 and RC with 2600 training points) lie in the same range, since both are implemented to stop at this threshold. However, the NGRC algorithm shows larger deviations in the 90% confidence interval. Nevertheless, to achieve this accurate prediction, a much smaller training data set is necessary for the NGRC algorithm. As for the Lorenz attractor, the RC is performed on the same training data set as the NGRC resulting in a much larger NRMSE in the testing phase. In this setting, the testing NRMSE does not achieve the desired threshold of 1 × 10−3 . The normalized maximal distance highlights the better performance of the NGRC algorithm. For 12
both components, the NGRC has the smallest normalized maximal distance compared to the RC method with 2600 and 600 training points. Moreover, Table 4 highlights the average computation time needed to calculate the prediction. The computation time is much smaller for the NGRC as for the RC algorithm, even if the training data set is the same. NRMSEtest ×10−4
NRMSEtrain ×10−4
ymax ×10−3
zmax ×10−3
comp. time 21 ±0.4ms
NGRC [600]
median
9.0
5.1
6.8
5.4
95%
[5.7, 80]
[3.2, 6.1]
[2.6, 82]
[2.0, 67]
RC [5200]
median
10
8.9
8.6
6.7
95%
[4.4, 15]
[5.6, 13]
[2.5, 14]
[1.9, 10]
RC [600]
median
17
10
14
5.4
95%
[6.9, 40]
[5.4, 18]
[4.0, 29]
[3.1, 22]
491 ±85ms 294 ±0.8ms
Table 4: Comparison of NRMSE, normalized maximal distance 8 and computation time for the Rössler system using the NGRC with 600 training points and RC with 600 and 2500 training points. The number of testing points is in all cases 1000. The median and the 95% confidence interval are calculated over 100 different initial conditions. The time average is calculated with the standard deviation over 7 runs, 100 loops each. Timestep ∆t = 0.1 and initial condition (x0 , y0 , z0 ) = (−0.6, 0.2, −0.1).
4.3. Inference of one Component of the El Niño–Southern Oscillation In this section we consider the ENSO as a climate application with real-world data sets. Given that the data points are observed on a monthly basis, there are only 516 available data points for the shortest time series which dates back to 1979. We take the same number of warm-up, training and testing points for both algorithms: 60 warm-up points, 180 training points and 276 testing points. To simplify the algorithms, the reservoir parameters for the RC algorithm are kept the same as for the Lorenz system. Since there is no visible improvement for higher k, only values of k between 1 and 10 are chosen. Additionally, 50 values of α ∈ [0.1, 10−8 ] are considered. Table 5 summarizes the different data sets used for training and testing, along with the resulting NRMSE values. Although the results are less accurate than for the Lorenz and Rössler systems, the model still shows a promising performance, especially given the complexity of real-world ENSO data. Fig. 6 and A.7–A.9 illustrate the inference of the different observables with the predictions in orange and the true target data in blue. In all cases, the predictions reflect key trends of the underlying dynamics quite well. Smaller-scale features are not yet fully reproduced. This is reflected in the fluctuations of the normalized absolute error, which is calculated for each time step as eabs :=
x̂ − x , xmax − xmin
(10)
where x is the observed data series and x̂ is the predicted time series. This error is usually between ±0.2 and indicates room for improvement. Comparing the algorithms, we note that the NGRC algorithm performs similar to the RC algorithm in the testing phase for all experiments. Again, the time needed to calculate the results is consistently and significantly smaller for the NGRC than for the RC algorithm.
13
input data SOI, SST, winds SST, ESPI, winds SOI, ESPI, winds SOI, ESPI, SST
target data
algorithm
NRMSEtest ×10−1
NRMSEtrain ×10−1
opt. k
time
dmax ×10−1
NGRC
0.9
1.0
1
383 µs ± 9 µs
3.6
RC
1.1
1.0
–
240 ms ± 2 ms
4.7
NGRC
1.1
1.2
1
389 µs ± 3 µs
4.1
RC
1.3
1.2
–
241 ms ± 0.5 ms
4.1
NGRC
1.1
0.8
2
1.1 ms ± 11 µs
3.5
RC
1.0
0.8
–
242 ms ± 0.6 ms
3.4
NGRC
1.5
1.2
1
390 µs ± 6 µs
4.4
RC
1.4
1.1
–
240 ms ± 0.5 ms
5.4
ESPI
SOI
SST
winds
Table 5: Comparison of NRMSE, normalized maximal distance 8 and computational time for the inference of different ENSO observables. We take 180 training points and 276 testing points. The optimal value for α is 0.1 in all cases. The average computation time with standard deviation is calculated over 7 runs, each 100 loops. 5. Discussion In this study we have demonstrated the accurate inference of two unknown components of the Lorenz and Rössler systems using the NGRC algorithm. The NRMSE is very small and there are no visible differences between the inferred values and the target data. Compared to the traditional RC algorithm, the NGRC algorithm required less training data and computation time. Given the simpler architecture of the NGRC algorithm, a shorter training time was expected. We considered different data sets when calculating the optimal number of time delayed steps k, sampled for 100 different initial conditions. For the Lorenz system, the optimal number of time delayed steps k was in most cases between 14 and 17, and for the Rössler system between 12 and 15 for a stopping threshold ε = 0.001. This showed that the NGRC algorithm shows an accurate performance also for different initial conditions. However, the values for the number of time delayed steps was unexpected as [32] only considered k ≤ 4. Our findings revealed a decrease of order −1 (Lorenz) and of order −2 (Rössler) in the NRMSE as k increases, emphasizing the importance of considering a larger number of time delayed steps for more accurate predictions. Interestingly, these trends begin to falter for very large k where the NRMSE actually begins to increase. We also showed the expected inversely proportional dependence between k and ∆t, which demonstrates that the physical time span covered by the time delayed steps, rather than the number of steps itself, is crucial for the prediction accuracy of the NGRC model. For the Lorenz system, we find k∆t = 0.76 ± 0.05. With coarser time step sizes, fewer time delayed steps are required since a larger segment of the curve is covered. However, to achieve a lower NRMSE, finer discretization is necessary for a more detailed representation of the dynamics, leading to an increase in the number of required time delayed steps to adequately capture the system’s dynamical behavior. However, a large number of time delayed steps results in higher computational 14
(a) Optimal NGRC performance, NRMSENGRC = 9.6 × 10−2 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
(b) Optimal RC performance, NRMSERC = 2.1 × 10−1 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
Figure 6: ENSO system inferring ESPI. ∆t = 1 represents monthly time steps. The left side shows the training, the right side the testing phase which are the same for both algorithms. costs, since the dimension of the feature vector increases significantly. Therefore, we found that selecting an appropriate time step size with a corresponding value of k is crucial, based on the desired task’s objectives. For the Rössler system with a stopping threshold ε = 10−4 we find k∆t = 1.44 ± 0.11. With coarser time step sizes, fewer time delayed steps are required, since a larger time span of the dynamics is captured. Even for very coarse time step sizes, an NRMSE below 10−4 is still reached. Therefore, to lower the computational costs, coarser time step sizes and therefore smaller k needs to be chosen. 15
Overall, neither algorithm consistently outperforms the other considering all metrics (NRMSE, maximal distance and computational time). The performance of NGRC and RC depends on the specific dynamical system and the evaluation metric. While RC achieves slightly lower prediction errors for the Lorenz system, the NGRC reaches slightly better prediction errors for the Rössler system. However, the simpler NGRC algorithm needs less training points and computational time for the results in all experiments. For the ENSO system, we achieved promising results, with the NGRC performing slightly better than the traditional RC. The same level of accuracy as the other two dynamical systems could not be achieved. Arguably, this was to be expected as the ENSO system is a real-world application and not a system described by a known ODE. Therefore, the underlying dynamics are mathematically not as well understood as for the other two and the measured data sets include noise, which generally requires tailored solutions to ensure accurate forecasting and inference (e.g. [60]). In conclusion, we have shown that the NGRC is a powerful tool for the inference of unknown parameters in chaotic systems. The results indicate that the choice between NGRC and RC involves a trade-off between accuracy and computational efficiency. As complex and chaotic behavior is a common feature in many dynamical systems, the NGRC algorithm has the potential to be an effective, accurate and fast solution for inferring observables across a wide range of applications. Acknowledgments This work emanated from the research funded by Taighde Éireann-Research Ireland (Grant No. FFPA/ 12066). We acknowledge Florian Stelzer for helpful discussions and for providing his reservoir computing code. The NGRC code used in this work is based on the implementation provided by [32]. References [1] R. Bakker, J. C. Schouten, C. L. Giles, F. Takens, C. M. v. d. Bleek, Learning chaotic attractors by neural networks, Neural Computation 12 (10) (2000) 2355–2383. doi:10.1162/089976600300014971. [2] S. Shahi, F. H. Fenton, E. M. Cherry, Prediction of chaotic time series using recurrent neural networks and reservoir computing techniques: A comparative study, Machine Learning with Applications 8 (2022) 100300. doi:10.1016/j.mlwa.2022.100300. [3] S. Dewitte, J. P. Cornelis, R. Müller, A. Munteanu, Artificial intelligence revolutionises weather forecast, climate monitoring and decadal prediction, Remote Sensing 13 (16) (2021) 3209. doi: 10.3390/rs13163209. [4] A. Dingli, , K. S. Fournier, Financial time series forecasting - a deep learning approach, International Journal of Machine Learning and Computing 7 (5) (2017) 118–122. doi:10.18178/ijmlc.2017.7.5. 632. [5] H. Zhao, A chaotic time series prediction based on neural network: Evidence from the shanghai composite index in china, in: 2009 International Conference on Test and Measurement, Vol. 2, 2009, pp. 382–385. doi:10.1109/ICTM.2009.5413024. [6] H. Jaeger, The" echo state" approach to analysing and training recurrent neural networks-with an erratum note’, Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148 (01 2001). [7] W. Maass, T. Natschläger, H. Markram, Real-time computing without stable states: A new framework for neural computation based on perturbations, Neural Computation 14 (11) (2002) 2531–2560. doi: 10.1162/089976602760407955. [8] M. Yan, C. Huang, P. Bienstman, P. Tino, W. Lin, J. Sun, Emerging opportunities and challenges for the future of reservoir computing, Nature Communications 15 (1) (2024) 2056. 16
[9] P. Vlachas, J. Pathak, B. Hunt, T. Sapsis, M. Girvan, E. Ott, P. Koumoutsakos, Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics, Neural Networks 126 (2020) 191–217. doi:10.1016/j.neunet.2020.02.016. [10] S. Bompas, B. Georgeot, D. Guéry-Odelin, Accuracy of neural networks for the simulation of chaotic dynamics: Precision of training data vs precision of the algorithm, Chaos 30 (113118) (2020). doi: 10.1063/5.0021264. [11] H. Jaeger, H. Haas, Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication, Science 304 (5667) (2004) 78–80. doi:10.1126/science.1091277. [12] A. Griffith, A. Pomerance, D. J. Gauthier, Forecasting chaotic systems with very low connectivity reservoir computers, Chaos 29 (123108) (2019). doi:10.1063/1.5120710. [13] Z. Lu, B. R. Hunt, E. Ott, Attractor reconstruction by machine learning, Chaos 28 (061104) (2018). doi:10.1063/1.5039508. [14] A. Röhm, D. J. Gauthier, I. Fischer, Model-free inference of unseen attractors: Reconstructing phase space features from a single noisy trajectory using reservoir computing, Chaos 31 (103127) (2021). doi:10.1063/5.0065813. [15] Z. Lu, J. Pathak, B. Hunt, M. Girvna, R. Brockett, E. Ott, Reservoir observers: Model-free inference of unmeasured variables in chaotic systems, Chaos 27 (041102) (2017). doi:10.1063/1.4979665. [16] L. Appeltant, M. Soriano, G. V. der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. Mirasso, I. Fischer, Information processing using a single dynamical node as complex system, Nature Communications 2 (1) (Sep. 2011). doi:10.1038/ncomms1476. [17] G. Van der Sande, D. Brunner, M. C. Soriano, Advances in photonic reservoir computing, Nanophotonics 6 (3) (2017) 561–576. doi:10.1515/nanoph-2016-0132. [18] L. Larger, A. Baylón-Fuentes, R. Martinenghi, V. S. Udaltsov, Y. K. Chembo, M. Jacquot, High-Speed Photonic Reservoir Computing Using a Time-Delay-Based Architecture: Million Words per Second Classification, Physical Review X 7 (1) (2017) 011015. doi:10.1103/PhysRevX.7.011015. [19] J. D. Hart, L. Larger, T. E. Murphy, R. Roy, Delayed dynamical systems: Networks, chimeras and reservoir computing, Phil. Trans. A. 377 (2153) (2019) 20180123. doi:10.1098/rsta.2018.0123. [20] F. Stelzer, A. Röhm, K. Lüdge, S. Yanchuk, Performance boost of time-delay reservoir computing by non-resonant clock cycle, Neural Networks 124 (2020) 158–169. arXiv:1905.02534, doi:10.1016/j. neunet.2020.01.010. [21] M. Goldmann, F. Köster, K. Lüdge, S. Yanchuk, Deep time-delay reservoir computing: Dynamics and memory capacity, Chaos 30 (9) (2020) 093124. doi:10.1063/5.0017974. [22] F. Koster, S. Yanchuk, K. Ludge, Master memory function for delay-based reservoir computers with single-variable dynamics, IEEE Transactions on Neural Networks and Learning Systems (2022) 1– 14doi:10.1109/tnnls.2022.3220532. URL https://doi.org/10.1109/tnnls.2022.3220532 [23] E. Bollt, On explaining the surprising success of reservoir computing forecaster of chaos? the universal machine learning dynamical system with contrast to var and dmd, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (1) (2021) 013108. doi:10.1063/5.0024890. [24] L. Jaurigue, K. Lüdge, Connecting reservoir computing with statistical forecasting and deep neural networks, Nature Communications 13 (1) (Jan. 2022). doi:10.1038/s41467-021-27715-5. 17
[25] S. Sugiura, R. Ariizumi, T. Asai, S.-i. Azuma, Necessary and sufficient reservoir condition for universal reservoir computing, Mathematics 13 (21) (2025) 3440. [26] M. Lukoševičius, A practical guide to applying echo state networks, in: Lecture Notes in Computer Science, Springer Berlin Heidelberg, 2012, pp. 659–686. doi:10.1007/978-3-642-35289-8_36. [27] Q. Song, Z. Feng, Effects of connectivity structure of complex echo state network on its prediction performance for nonlinear time series, Neurocomputing 73 (10-12) (2010) 2177–2185. doi:10.1016/j. neucom.2010.01.015. [28] H. Jaeger, M. Lukoševičius, D. Popovici, U. Siewert, Optimization and applications of echo state networks with leaky- integrator neurons, Neural Networks 20 (3) (2007) 335–352. doi:10.1016/j. neunet.2007.04.016. [29] A. Rodan, P. Tino, Minimum complexity echo state network, IEEE Transactions on Neural Networks 22 (1) (2011) 131–144. doi:10.1109/TNN.2010.2089641. [30] J. Yperman, T. Becker, Bayesian optimization of hyper-parameters in reservoir computing (2016). doi:10.48550/ARXIV.1611.05193. [31] J. R. Maat, N. Gianniotis, P. Protopapas, Efficient optimization of echo state networks for time series datasets, in: 2018 International Joint Conference on Neural Networks (IJCNN), 2018, pp. 1–7. doi: 10.1109/IJCNN.2018.8489094. [32] D. J. Gauthier, E. Bollt, A. Griffith, W. A. S. Barbosa, Next generation reservoir computing, Nature Communications 12 (5564) (2021). doi:10.1038/s41467-021-25801-2. [33] A. Haluszczynski, D. Koeglmayr, C. Räth, Controlling dynamical systems to complex target states using machine learning: next-generation vs. classical reservoir computing, in: 2023 International Joint Conference on Neural Networks (IJCNN), 2023, pp. 1–7. doi:10.1109/IJCNN54540.2023.10191257. [34] S. Talukder, J. J. Sun, M. Leonard, B. W. Brunton, Y. Yue, Deep neural imputation: A framework for recovering incomplete brain recordings, arXiv preprint arXiv:2206.08094 (2022). [35] W. Qiu, Z. Sun, J. Wang, Y. Bai, S. Yao, D. Guo, G. Yang, A compressed sensing based framework for surface pressure field reconstruction from sparse measurement, Physics of Fluids 37 (4) (2025). [36] J. Chao, B. Pan, Q. Chen, S. Yang, J. Wang, C. Nai, Y. Zheng, X. Li, H. Yuan, X. Chen, et al., Learning to infer weather states using partial observations, Journal of Geophysical Research: Machine Learning and Computation 2 (1) (2025) e2024JH000260. [37] C. Kadow, D. M. Hall, U. Ulbrich, Artificial intelligence reconstructs missing climate information, Nature Geoscience 13 (6) (2020) 408–413. [38] W. Konrad, G. Katul, A. Roth-Nebelsick, M. Grein, A reduced order model to analytically infer atmospheric co2 concentration from stomatal and climate data, Advances in Water Resources 104 (2017) 145–157. [39] E. N. Lorenz, Deterministic nonperiodic flow, Journal of the Atmospheric Sciences 20 (2) (1963) 130– 141. doi:10.1175/1520-0469(1963)020<0130:dnf>2.0.co;2. [40] O. Rössler, An equation for continuous chaos, Physics Letters A 57 (5) (1976) 397–398. doi:10.1016/ 0375-9601(76)90101-8. [41] G. WALKER, WORLD WEATHER, Monthly Weather Review 56 (5) (1928) 167–170. doi:10.1175/ 1520-0493(1928)56<167:ww>2.0.co;2. 18
[42] E. S. Sarachik, M. A. Cane, The El Niño-Southern Oscillation Phenomenon, Cambridge University Press, 2010. doi:10.1017/CBO9780511817496. [43] R. F. Williams, The structure of lorenz attractors, Publications Mathematiques de l’IHES 50 (1979) 73–99. [44] W. Tucker, The lorenz attractor exists, Comptes Rendus de l’Académie des Sciences-Series IMathematics 328 (12) (1999) 1197–1202. [45] D. Wilhite, D. Wood, S. Meyer, M. Glantz, R. Katz, Climate crisis (1987). [46] R. B. Aronson, W. F. Precht, I. G. Macintyre, T. J. T. Murdoch, Coral bleach-out in belize, Nature 405 (6782) (2000) 36–36. doi:10.1038/35011132. [47] P. W. Glynn, W. H. de Weerdt, Elimination of two reef-building hydrocorals following the 1982-83 el niño warming event, Science 253 (5015) (1991) 69–71. doi:10.1126/science.253.5015.69. [48] J. Y. Ng, S. W. D. Turner, S. Galelli, Influence of el niño southern oscillation on global hydropower production, Environmental Research Letters 12 (3) (2017) 034010. doi:10.1088/1748-9326/aa5ef8. [49] W. Cai, S. Borlace, M. Lengaigne, P. van Rensch, M. Collins, G. Vecchi, A. Timmermann, A. Santoso, M. J. McPhaden, L. Wu, M. H. England, G. Wang, E. Guilyardi, F.-F. Jin, Increasing frequency of extreme el niño events due to greenhouse warming, Nature Climate Change 4 (2) (2014) 111–116. doi:10.1038/nclimate2100. [50] E. M. Vincent, M. Lengaigne, C. E. Menkes, N. C. Jourdain, P. Marchesiello, G. Madec, Interannual variability of the south pacific convergence zone and implications for tropical cyclone genesis, Climate Dynamics 36 (9-10) (2009) 1881–1896. doi:10.1007/s00382-009-0716-3. [51] S. Philander, Meteorology: Anomalous el niño of 1982–83, Nature 305 (5929) (1983) 16–16. doi: 10.1038/305016a0. [52] R. W. Reynolds, N. A. Rayner, T. M. Smith, D. C. Stokes, W. Wang, An improved in situ and satellite SST analysis for climate, Journal of Climate 15 (13) (2002) 1609–1625. doi:10.1175/1520-0442(2002) 015<1609:aiisas>2.0.co;2. [53] C. F. Ropelewski, P. D. Jones, An extension of the tahiti-darwin southern oscillation index, Monthly Weather Review 115 (9) (1987) 2161–2165. doi:10.1175/1520-0493(1987)115<2161:aeotts>2.0. co;2. [54] S. Curtis, R. Adler, ENSO indices based on patterns of satellite-derived precipitation, Journal of Climate 13 (15) (2000) 2786–2793. doi:10.1175/1520-0442(2000)013<2786:eibopo>2.0.co;2. [55] K. E. Trenberth, The definition of el niño, Bull. Am. Meteorol. Soc. 78 (12) (1997) 2771–2777. [56] E. M. Rasmusson, J. M. Wallace, Meteorological aspects of the el niño/southern oscillation, Science 222 (4629) (1983) 1195–1202. [57] J. Bjerknes, Atmospheric teleconnections from the equatorial pacific1, Mon. Weather Rev. 97 (3) (1969) 163–172. [58] T. Jinno, T. Mitsui, K. Nakai, Y. Saiki, T. Yoneda, Long-term prediction of el niño-southern oscillation using reservoir computing with data-driven realtime filter, Chaos: An Interdisciplinary Journal of Nonlinear Science 35 (5) (2025). [59] F. Guardamagna, C. Wieners, H. A. Dijkstra, Explaining the high skill of reservoir computing methods in el niño prediction, Nonlinear Processes in Geophysics 32 (2) (2025) 201–224. 19
[60] G. A. Gottwald, S. Reich, Combining machine learning and data assimilation to forecast dynamical systems from noisy partial observations, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (10) (2021).
20
Appendix A. Inferring components of ENSO
(a) Optimal NGRC performance, NRMSENGRC = 1.1 × 10−1 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
(b) Optimal RC performance, NRMSERC = 2.0 × 10−1 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
Figure A.7: ENSO system inferring SOI. ∆t = 1 represents monthly time steps. The left side shows the training, the right side the testing phase which are the same for both algorithms.
21
(a) Optimal NGRC performance, NRMSENGRC = 1.1 × 10−1 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
(b) Optimal RC performance, NRMSERC = 2.7 × 10−1 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
Figure A.8: ENSO system inferring SST. ∆t = 1 represents monthly time steps. The left side shows the training, the right side the testing phase which are the same for both algorithms.
22
(a) Optimal NGRC performance, NRMSENGRC = 1.5 × 10−1 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
(b) Optimal RC performance, NRMSERC = 2.6 × 10−1 . Top panel: The target is shown in blue, the prediction with dashed orange. Bottom panel: normalized absolute error in each time step.
Figure A.9: ENSO system inferring zonal winds. ∆t = 1 represents monthly time steps. The left side shows the training, the right side the testing phase which are the same for both algorithms.
23