Stabilizing chaotic dynamical system reproduction in reservoir computing

TL;DR

This work tackles the instability of autonomous reservoir computing when reproducing chaotic dynamics by identifying spurious unstable and neutral modes as the primary cause. It introduces a simple, deterministic input-layer design that constructs the input weights $\mathbf{W}_{\mathrm{in}}$ from a $D$-dimensional invariant subspace of the reservoir matrix $\mathbf{A}$, thereby constraining the closed-loop spectrum so that only the slow, learnable modes are adjustable and the rest remain fast and stable. Empirically, this approach yields dramatically improved robustness to initialization and noise, longer valid prediction times, and faithful estimation of the full Lyapunov spectrum across seeds and across 135 chaotic systems, indicating topological semi-conjugacy to the target dynamics. The method reduces hyperparameter tuning, offers practical guidelines (e.g., using a symmetric $\mathbf{A}$ and small spectral radius $\rho$), and has potential for stable, physics-based RC implementations and broader applicability to other state-space models. Overall, the paper provides a principled design principle that stabilizes RC attractor reconstruction and enhances reliability for data-driven chaos modeling.

Abstract

Reservoir Computing (RC), a type of recurrent random neural network, is a powerful framework for modeling complex and chaotic dynamics. However, its autonomous (closed-loop) operation is often plagued by inherent instability. Moreover, performance is highly sensitive to the reservoir's random initialization, leading to vulnerability to noise and/or behaviour that bears no resemblance whatsoever to the target dynamical system. Here we identify a primary cause of this unreliability: the emergence of excessive, spurious unstable or neutral modes in the closed-loop dynamics. We introduce a simple deterministic input layer design principle that directly addresses this vulnerability by structurally suppressing the emergence of these spurious modes a priori (before training). Our approach dramatically improves robustness to both initialization sensitivity and internal noise, doubling the prediction horizon. Furthermore, we demonstrate on chaotic dynamical systems that this design enables robust estimation of the full Lyapunov spectrum (100\% success rate across 50 seeds), signifying that the autonomous RC faithfully emulates the essential properties of the target dynamical system. This work provides a systematic explanation for a common RC failure mode and offers a concrete design guideline, advancing RCs from heuristic trial-and-error tuning toward a reliable tool for modeling complex systems.

Figures (9)

• Figure 1: a, Schematic of the Reservoir Computing (RC) framework. During the training phase (switch left, open-loop), the RC is driven by the observed time series $\mathbf{u}_t$. For autonomous prediction (switch down, closed-loop), the RC uses its own one-step-ahead prediction $\hat{\mathbf{u}}_t$ as the input for the next time step $t+1$. b, General formulation of the discrete-time RC dynamics employed in this study. c, Conceptual diagram of topological (semi-)conjugacy. An ideally designed closed-loop RC dynamics $\hat{f}$ becomes topologically (semi-)conjugate to the target dynamical system $f$, enabling faithful emulation of the original dynamics. d--f, Performance of the proposed deterministic input layer design. d, Eigenvalue spectrum of the closed-loop matrix $\mathbf{W}_{\mathrm{cl}}$. By design, only the three eigenvalues corresponding to the learnable modes exhibit large magnitudes, while the remaining non-principal eigenvalues are strictly confined within the stable region defined by the spectral radius of $\mathbf{A}$ (green circle), suppressing spurious modes. e, Schematic of the reservoir state space. The reconstruction of the target system's invariant set is confined to a low-dimensional attractive manifold within the reservoir space, ensuring transversal stability. f, Long-term autonomous prediction trajectory, accurately and stably reproducing the Lorenz attractor. g--i, Instability in the conventional (random) input layer design. g, Eigenvalue spectrum of $\mathbf{W}_{\mathrm{cl}}$, showing the uncontrolled eigenvalue spectrum and the emergence of spurious modes, where some eigenvalues spill over the spectral radius of $\mathbf{A}$ (green circle). h, Instability induced by spurious modes causes trajectories to diverge from the reconstruction. i, Long-term autonomous prediction trajectory, which diverges from the true attractor and exhibits spurious behavior unfaithful to the target dynamics.
• Figure 2: a--c, Eigenvalue distributions of the closed-loop matrix $\mathbf{W}_{\mathrm{cl}}$. a, b, Conventional random designs (using different random seeds) exhibit eigenvalues that spill over the spectral radius of $\mathbf{A}$ (green dashed circle). c, The proposed deterministic design restricts the spectrum such that only the $D=3$ eigenvalues corresponding to the learnable modes lie outside the circle, with all remaining modes strictly confined within. d--f, Closed-loop prediction trajectories (red) for the Rössler system plotted against the ground truth (blue), using the RCs corresponding to a, b, and c, respectively. g--i, Basin stability analysis corresponding to d--f. The plots show the initial points ($\hat{\mathbf{u}}_{t_0} = \mathbf{W}_{\mathrm{out}} \mathbf{r}_{t_0}$, derived from perturbed reservoir states) which subsequently deviated from the original attractor (red) versus those that successfully converged (green). Note that the conventional model in h appears stable in trajectory (e) but is actually fragile to perturbations. j, Dependence of Basin Stability ($S_B$) on the spectral radius $\rho$ of $\mathbf{A}$. The proposed method (green) consistently outperforms the conventional method (orange). Error bars represent the standard deviation over 10 random seeds.
• Figure 3: a--c, Performance evaluation on the Lorenz system. a, Comparison of closed-loop prediction trajectories across 50 random initialization seeds. The plots display the best (longest VPT; red/green dashed lines), worst (shortest VPT; orange/blue dashed lines), and average performance (purple/limegreen solid lines) for the conventional and proposed methods, respectively. b, c, Estimated Lyapunov spectra $\hat{\lambda}_1, \hat{\lambda}_2, \hat{\lambda}_3, \hat{\lambda}_4$ compared to the true values (dashed line). The y-axis uses a symmetric log scale (linear within $[-1, 1]$). b, The proposed method successfully estimates the full spectrum across all 50 seeds. c, The conventional method failed in 36 out of 50 seeds, frequently exhibiting spurious positive exponents. d--f, Corresponding analysis for the Rössler system. d, VPT comparison demonstrating the superior robustness of the proposed method (green/blue/limegreen) over the conventional method (red/orange/purple). e, f, Lyapunov spectrum estimation (linear range $[-0.1, 0.1]$). The proposed method (e) accurately captures the true spectrum across all seeds, whereas the conventional method (f) fails in the majority of trials.
• Figure 4: a--c, Dependence of Valid Prediction Time (VPT) on the invariant subspace selection strategy. The x-axis corresponds to the indices of eigenvalues (and associated eigenvectors) used to construct $\mathbf{W}_{\mathrm{in}}$, sorted in ascending order of their real parts. For instance, the label eigenvec_[1, Re(2), Im(2)] denotes that $\mathbf{W}_{\mathrm{in}}$ was constructed using the eigenvector for the 1st real eigenvalue and the real/imaginary parts of the eigenvector pair for the 2nd complex conjugate eigenvalue. The 'random' bar represents the conventional approach. Results are shown for different reservoir topologies: a, random matrix; b, Erdős–Rényi random graph; and c, symmetric matrix. Error bars indicate the standard deviation across 10 random initialization seeds (averages are calculated over 100 prediction start points). d, Comprehensive parameter sensitivity analysis. Heatmaps display the average VPT as a function of the spectral radius $\rho$ (y-axis) and regularization parameter $\beta$ (x-axis). The grid layout is organized by reservoir topology (columns) and $\mathbf{W}_{\mathrm{in}}$ construction strategy (rows). Strategies include random and the proposed methods selecting eigenvectors with maximal (high), median (intermediate), and minimal (low) real parts.
• Figure S1: Prediction performance evaluated using two metrics: Valid Prediction Time (VPT, left side) and prediction duration based on Symmetric Mean Absolute Percentage Error (SMAPE, right side). Systems are categorized as non-delayed autonomous (114 systems, blue), delayed autonomous (15 systems, orange), and non-autonomous (6 systems, green). Data points show the average over $N_{\mathrm{test}}=100$ prediction start points for each system, further averaged across 10 initial RC seeds; error bars represent the minimum and maximum (min/max) values across these 10 seeds. See Figs. S2–S4 for results under different conditions using the same evaluation metrics and format.

Theorems & Definitions (8)

• Theorem 1: Invariant Eigenvalues
• proof
• Corollary 2: Inheritance of Invariant Subspace
• proof
• Theorem 1: Invariant Eigenvalues
• proof
• Corollary 2: Inheritance of Invariant Subspace
• proof