Next generation reservoir computing

The study addresses how to efficiently learn and forecast the behavior of dynamical systems—often chaotic or spatiotemporal—using limited observed time-series data. While reservoir computing (RC) is effective and data-efficient compared to many ML approaches, it relies on randomly sampled connectivity matrices and numerous metaparameters, creating variability and requiring costly optimization. The paper builds on theoretical results showing that an RC with linear nodes and nonlinear readout is equivalent to nonlinear vector autoregression (NVAR), implying that a reservoir is not explicitly required. The authors pose the question of whether an explicitly constructed NVAR—termed the next generation reservoir computer (NG-RC)—can match or exceed traditional RC performance on core benchmarks while reducing training data, computation time, and hyperparameter tuning. The purpose is to demonstrate NG-RC’s efficacy on three benchmarks: short-term forecasting, long-term climate (attractor) reproduction, and inference of unseen variables, and to highlight its advantages and interpretability. The work is significant because it simplifies RC deployment, improves computational efficiency, and provides theoretical clarity by linking RC to NVAR/Volterra representations.

The paper situates NG-RC within RC and nonlinear system identification literature. Classical RC uses fixed random recurrent networks with trained linear readout and has shown strong performance on chaotic and high-dimensional systems but suffers from metaparameter sensitivity and random matrix selection challenges. Recent theory establishes universality of RC both with nonlinear reservoir activations and with linear reservoirs plus nonlinear readouts, with NVAR shown to be equivalent to a linear RC with polynomial readout. Prior efforts have used synchronization criteria to guide RC design, but no general design rules exist. The NG-RC connects to NARX (nonlinear autoregression with exogenous inputs) methods developed since the 1980s, but distinguishes itself via Tikhonov regularization and a theoretical equivalence to an implicit RC. Related data-driven linearization and sparse system identification (e.g., SINDy, LASSO, information-theoretic methods) aim to recover vector fields, while NG-RC learns the flow (finite-time mapping), which can differ from the vector field. Prior RC studies on large spatiotemporal systems suggest scalability, and recent domain-driven and simplified NG-RC-like models show promising results but without the broader theoretical framework presented here.

The authors first review traditional RC. A recurrent network with N nodes receives input X via fixed random input weights W and evolves according to r_{i+1} = (1−γ) r_i + γ f(A r_i + W X_i + b), with nonlinear activation f (commonly tanh). The output Y_{i+1} = W_out O_total,i+1 is a linear transformation of features O_total constructed from the reservoir state. Training uses ridge-regularized least squares: W_out = Y_d O_total^T (O_total O_total^T + α I)^{-1}. They then consider RC variants with linear reservoir nodes (f(r) = r) and nonlinear output features, e.g., O_total = r ⊕ (r ⊙ r), maintaining universal approximation capability. NG-RC (NVAR) constructs features directly from the observed input time series without an explicit reservoir. The linear feature block O_lin,i concatenates k time-delay vectors of X spaced by s steps: O_lin,i = [X_i, X_{i−s}, ..., X_{i−(k−1)s}]^T (dimension d k). Nonlinear features O_nonlinear are formed as polynomial monomials in O_lin (Volterra series perspective), e.g., via upper-triangular elements of outer products for quadratic or higher-order tensor products for higher degrees p. The total feature vector is O_total = c O_lin ⊕ O_nonlinear. Training uses the same ridge-regularized least squares as RC. Warm-up requires only s k steps to populate delays—much shorter than traditional RC warm-up. For driven systems or parameter-varying systems, O_lin(t) includes the drive signal and accessible parameters. Polynomial nonlinearities are used (low-order sufficient in practice). The authors note the proven equivalence between NVAR and linear RC with polynomial readout, implying NG-RC implicitly defines an RC but avoids explicit network construction and metaparameter optimization. Tasks and models:

• Lorenz63 system: ẋ = 10(y − x), ẏ = x(28 − z) − y, ż = x y − (8/3) z; Lyapunov time ≈ 1.1.
• Double-scroll electronic circuit (dimensionless): V̇1 = V1/R − ΔV/R2 − 2 I_s sinh(β ΔV), V̇2 = ΔV/R2 + 2 I_s sinh(β ΔV) − I, Ī = V2 − R1 I, with ΔV=V1−V2 and parameters R1=1.2, R2=3.44, R1 (typo likely R3)=0.193, β=11.6, I_s=2.25×10^{-5}; Lyapunov time 7.81. Chosen to test nonpolynomial vector fields; polynomial features still effective due to Volterra theory. Forecasting setup (tasks 1 and 2): Train NG-RC to perform one-step-ahead prediction using X(t+Δt) = W_out O_total(t) with ridge regularization. In testing, close the loop by feeding predictions back as inputs, making NG-RC autonomous. Feature choices:
• Lorenz63 forecasting: O_total = O_lin ⊕ O_nonlinear^(3) (cubic polynomial terms).
• Double-scroll forecasting: To respect odd symmetry and zero mean, O_total = O_lin ⊕ O_nonlinear^(2) (quadratic terms) as defined in text (the counts reflect monomial sets per chosen degree and embedding size). To emphasize learning the flow, they use a residual/Euler-like formulation: X_{t+1} = X_t + W_out O_total,t. Training uses short datasets and minimal warm-up. Inference of unseen data (task 3): Train with full Lorenz63 state [x,y,z] to predict one-step-ahead z from [x,y,z]. In testing, provide only [x,y] and infer z, using time-delay embedding guided by Takens’ theorem. Hyperparameters include Δt, k (number of delays), s (skip between delays), polynomial degree, and ridge α. Implementation details: Python 3.7.9, NumPy 1.20.2, SciPy 1.6.2 on x86-64 Windows 10; code is publicly available.

• Lorenz63 forecasting (task 1 and climate reproduction, task 2): Trained on 400 data points per variable (Δt = 0.025, k = 2, s = 1). Feature vector had 28 components; W_out dimension 3×28. Training time < 10 ms on a single-core desktop. During training, NRMSE over one Lyapunov time: (1.06 ± 0.01) × 10^{-4}. In autonomous prediction, NG-RC reproduces the strange attractor and forecasts accurately for ~5 Lyapunov times; testing NRMSE: (2.40 ± 0.53) × 10^{-3}. Many nonzero W_out components differ from explicit vector field terms, reflecting learning of the flow rather than the vector field.
• Double-scroll forecasting: To match Lyapunov-time coverage, training extended from 10 to 100 time units with Δt = 0.25. Using cubic monomials per Eq. 10 with d = 3, k = 2, s = 1 resulted in 62 features; 400 points per variable used for training. NG-RC accurately reconstructs the attractor and forecasts with testing NRMSE = (4.5 ± 1.0) × 10^{-3}.
• Inference of unseen variable (task 3): Using Lorenz63 with k = 4, s = 5, Δt = 0.05 for embedding; total features = 45; W_out dimension 1×45. Training NRMSE for z from [x,y,z]: (9.5 ± 0.1) × 10^{-3} (α = 0.05). In testing, providing only [x,y], inferred z with NRMSE = (1.75 ± 0.3) × 10^{-2}, with predicted traces overlaying ground truth closely.
• Computational efficiency: For Lorenz63 forecasting, NG-RC is estimated to be ~33–162× less costly than a typical efficient traditional RC and >10^6× less costly than a high-accuracy RC for a single metaparameter set. For the double-scroll task (cubic nonlinearity), NG-RC is ~8–41× less costly than a typical efficient RC per metaparameter set.
• Additional advantages: Extremely short warm-up (only s k steps, e.g., 2 points for s=1, k=2) vs 10^3–10^5 for RC; fewer metaparameters; robustness to noise via ridge regularization (Supplementary Note 5 shows noise-driven Lorenz63 performance improves with higher α).

The findings demonstrate that an explicitly constructed NVAR (NG-RC) can match state-of-the-art RC performance on key benchmarks—short-term forecasting, climate reproduction, and inferring unobserved variables—while substantially reducing training data, warm-up, computational cost, and metaparameter burden. This addresses the central challenge of RC design complexity and randomness by removing the need for random reservoirs and offering interpretable time-delay polynomial features. The NG-RC’s success aligns with theory proving equivalence between NVAR and linear RC with a polynomial readout; thus, NG-RC implicitly realizes a traditional RC without constructing it. The results highlight that NG-RC learns the flow map (finite-time evolution) rather than the vector field, explaining why significant W_out components need not mirror the dynamical equations’ terms. Compared to NARX, NG-RC adds principled regularization and a theoretical link to RC, yielding improved performance and robustness. Feature-space parsimony is possible since some weights are small; future reduction can leverage LASSO or information-theoretic criteria, or kernel methods when feature dimension is large relative to data. Regularization confers noise tolerance, enabling recovery of effective noise-free dynamics with appropriate α. Prior RC successes on high-dimensional spatiotemporal systems suggest NG-RC will scale favorably due to its implicit RC and leaner parameterization. Potential applications include digital twins, hybrid data assimilation combining models with observations, adaptive nonlinear control, and accelerating turbulent simulations.

The paper introduces the next generation reservoir computer (NG-RC), an NVAR-based approach that is theoretically equivalent to a linear reservoir computer with a polynomial readout but eliminates the need for explicit random reservoirs and extensive hyperparameter tuning. Empirical results on Lorenz63 and a nonpolynomial double-scroll circuit show state-of-the-art forecasting and attractor reconstruction using only 400 training points per variable, very short warm-up, and orders-of-magnitude lower computational cost than traditional RC. NG-RC also accurately infers unobserved variables using delay embeddings, underscoring its practical utility. Future research directions include systematic feature selection and sparsification (e.g., LASSO, information-theoretic methods), exploring kernelized variants when feature dimensionality becomes large, rigorous noise studies, and extending to complex high-dimensional spatiotemporal systems and real-world datasets. Applications are envisaged in digital twins, data assimilation, nonlinear control, and fast surrogate modeling for multiscale physics.

The study evaluates NG-RC primarily on low-dimensional chaotic systems (Lorenz63 and a double-scroll circuit) using data from noise-free numerical simulations; real-world data and high-dimensional systems are not directly tested in the main text. Performance metrics can vary slightly with software implementations. Although ridge regularization confers noise tolerance (as shown in supplementary results), a comprehensive treatment of measurement and process noise is outside the main scope. Feature selection was not optimized—many small weights suggest potential redundancy—and alternative regularizers or selection schemes were not systematically compared. Finally, hyperparameter choices (k, s, polynomial degree, Δt, α) were selected based on heuristic considerations; automated tuning strategies were not explored.