Source-linked AI summary
Backpropagation Algorithms and Reservoir Computing in Recurrent Neural Networks for the Forecasting of Complex Spatiotemporal Dynamics
Pantelis R. Vlachas, Jaideep Pathak, Brian R. Hunt, Themistoklis P. Sapsis, Michelle Girvan, Edward Ott, Petros Koumoutsakos
TL;DR
The paper asks how Reservoir Computing and backpropagation-trained RNNs compare for forecasting full-order and reduced-order chaotic spatiotemporal dynamics. It evaluates these methods on Lorenz-96 and Kuramoto-Sivashinsky systems, including parallel architectures and Lyapunov spectra. RC performs strongly for full-state and parallel forecasting with lower training time, while gated RNNs are more robust for reduced-order data; both capture positive KS Lyapunov exponents but miss two zero exponents.
Problem
Comparative evidence is limited on the accuracy, efficiency, and scalability of RNNs trained by BPTT versus RC for full-order and reduced-order chaotic systems.
Method
The study benchmarks RC and gated or unitary RNNs on Lorenz-96 and Kuramoto-Sivashinsky forecasting, using parallel architectures and Lyapunov-spectrum evaluation.
Results
RC is competitive or superior for full-state and parallel forecasting with lower training time, whereas LSTM and GRU forecast reduced-order dynamics more stably; both match positive KS Lyapunov exponents.
Takeaways & Limitations
RNNs provide a computational framework for learning and forecasting complex spatiotemporal dynamics across full-state, reduced-order, and parallel settings.
Takeaways & Limitations
Neither RC nor GRU reproduces the KS zero Lyapunov exponents Λ7 and Λ8, although augmenting the spectra with them yields excellent agreement.
Abstract
from arXiv · showhide
We examine the efficiency of Recurrent Neural Networks in forecasting the spatiotemporal dynamics of high dimensional and reduced order complex systems using Reservoir Computing (RC) and Backpropagation through time (BPTT) for gated network architectures. We highlight advantages and limitations of each method and discuss their implementation for parallel computing architectures. We quantify the relative prediction accuracy of these algorithms for the longterm forecasting of chaotic systems using as benchmarks the Lorenz-96 and the Kuramoto-Sivashinsky (KS) equations. We find that, when the full state dynamics are available for training, RC outperforms BPTT approaches in terms of predictive performance and in capturing of the long-term statistics, while at the same time requiring much less training time. However, in the case of reduced order data, large scale RC models can be unstable and more likely than the BPTT algorithms to diverge. In contrast, RNNs trained via BPTT show superior forecasting abilities and capture well the dynamics of reduced order systems. Furthermore, the present study quantifies for the first time the Lyapunov Spectrum of the KS equation with BPTT, achieving similar accuracy as RC. This study establishes that RNNs are a potent computational framework for the learning and forecasting of complex spatiotemporal systems.
1. Introduction
The paper addresses limited comparative evidence on RNNs and Reservoir Computing for model-agnostic forecasting of chaotic dynamical systems. It compares their accuracy, computational efficiency, and scalability for full-order and reduced-order dynamics.
- Earlier neural-network studies often focused on intrinsically low-order systems or used dimensionality reduction tools.
- The study compares backpropagation-trained RNNs and Reservoir Computing for forecasting fully observed and partially observed chaotic systems.
- Limited comparative studies have evaluated these two frameworks on challenging physical problems.
- RNNs can capture temporal dependencies and implicitly identify embeddings needed to forecast dynamical systems.
- The paper also examines parallel prediction of high-dimensional dynamics, model efficiency, and implementation limitations.
2. Methods - Sequence Modeling
The methods use recurrent hidden states to forecast observables, with architectures differing in recurrent mappings and training procedures. The study combines gated and unitary RNNs trained by BPTT with RC, which trains only output weights.
- Each model maps the current observable and hidden state to the next observable forecast through a recurrent hidden-state representation.
- RC randomly selects recurrent weights to satisfy the echo state property and trains only the output weights, for example with ridge regression.
- LSTM and GRU cells use gates to retain, forget, and update information, addressing vanishing-gradient problems in plain Elman RNNs.
- Unitary RNNs constrain the recurrent matrix to complex unitary matrices with spectral radius one to alleviate vanishing gradients.
- BPTT computes loss derivatives through repeated chain-rule applications and updates network parameters with first-order optimization.
- BPTT training uses forward horizon κ1, truncation length κ2, and teacher-forcing interval κ3; stateless models cannot learn dependencies longer than κ2.
3. Comparison Metrics
The paper evaluates forecasting accuracy, long-term statistical fidelity, and chaotic dynamics using normalized error, valid prediction time, power spectra, and Lyapunov exponents. These measures separate short-term prediction quality from long-term dynamical agreement.
- Model evaluations use a single Nvidia Tesla P100 GPU and extensive hyperparameter grid searches.
- Normalized root mean square error compares forecasts with targets after scaling each state component by its temporal standard deviation.
- The maximal Lyapunov exponent provides a characteristic timescale for judging prediction-error growth in chaotic systems.
- Valid prediction time is the largest forecast time whose NRMSE remains below ε=0.5, normalized by the system’s maximal Lyapunov exponent.
- Mean power spectral density across observable components measures whether forecasts reproduce long-term frequency statistics.
4. Forecasting Reduced Order Observable Dynamics in the Lorenz-96
The Lorenz-96 experiments compare RC and BPTT-trained recurrent architectures for forecasting full-state and reduced-order observables. RC performs best with full-state information, whereas GRU and LSTM are more stable and accurate when temporal dependencies must be inferred from reduced-order data.
- Experimental setting: Accurate long-term forecasting is difficult because small initial errors can grow exponentially in chaotic systems, with the maximal Lyapunov exponent setting a characteristic propagation timescale.The study therefore evaluates forecasting of observable dynamics rather than assuming direct access to the full state.
- Experimental setting: The observable is formed by transforming the Lorenz-96 state into SVD modes and retaining either the 35 most energetic modes or all 40 modes.The 40-dimensional case corresponds to full-state information, while the 35-dimensional case is reduced order.
- Reduced-order forecasting: At F=10, LSTM and GRU predict up to 2 Lyapunov times before NRMSE reaches 1, compared with 1 Lyapunov time for RC and Unitary RNNs.Using the stricter NRMSE threshold of 0.5, all models remain useful for only one Lyapunov time on reduced-order data.
- Full-state forecasting: RC achieves the highest full-state VPT, reaching approximately 2.31 for F=8 and 2.35 for F=10.For F=8, GRU, LSTM, and Unitary RNNs reach approximately 1.34, 0.97, and 0.58, respectively; for F=10, RC remains highest among all models.
- Reduced-order forecasting: In reduced-order forecasting, GRU and LSTM outperform RC, with GRU reaching approximately 0.98 VPT at F=8 and both gated models reaching about 1.22 and 1.17 at F=10.RC reaches only approximately 0.55 at F=8 and 0.60 at F=10, while Unitary RNNs show the lowest forecasting ability.
- Long-term statistics: Full-state power spectra match ground truth across models, whereas reduced-order RC fails to match the statistics and GRU/LSTM retain superior forecasting and long-term-statistics performance.The full-state result supports RC as a strong predictive tool when the complete state is available.
- Generalization and limitations: Reduced-order RC models overfit more readily and can diverge because reservoir size may be insufficient for high-dimensional, non-local reduced-order dynamics.GRU and LSTM use validation-based early stopping, while RC performance remains inferior despite tested Tikhonov regularization values.
5. Parallel Forecasting Leveraging Local Interactions
The paper uses ensembles of locally interacting RNNs to forecast high-dimensional spatial dynamics in parallel. In Lorenz-96 and KS experiments, RC performs strongly for full-state forecasting, while gated BPTT networks can match or exceed it with smaller models in some settings.
- Parallel architecture: Local parallel models split the observable into spatial groups, with each RNN forecasting one group using neighboring inputs.Networks train independently, while communication between neighboring groups is required during long-term forecasting.
- Lorenz-96 setup: The parallel architecture replaces one global Lorenz-96 model with 20 RNNs, each forecasting 2 nodes while using 4 neighboring nodes on each side.
- Lorenz-96 results: RC models with 3000–6000 nodes require much less training time than BPTT networks, while achieving slightly lower VPT in the full-state Lorenz-96 experiment.RC models with 12000 nodes require more than 3 GB of RAM per rank, and RC sizes above 12000 were infeasible because of memory requirements.
- Lorenz-96 results: RC with 6000 or 12000 nodes is needed to match the predictive performance of an LSTM with 100 nodes in full-state Lorenz-96 forecasting.All evaluated models accurately capture the reconstructed power spectrum.
- Kuramoto-Sivashinsky results: In KS forecasting, an RC with 3000 hidden units reaches VPT ≈4.8, while LSTM models reach ≈4 and GRU models reach ≈3.5.The 3000-node RC reaches or surpasses comparable RNN performance with similar RAM usage and much less training time.
- Kuramoto-Sivashinsky results: The RC, GRU, and LSTM models successfully reproduce KS power-spectrum statistics, with a small 80-unit LSTM matching a 1000-unit RC.
6. Calculation of Lyapunov Exponents in the Kuramoto-Sivashinsky Equation
The paper evaluates whether recurrent surrogates can recover the Lyapunov spectrum of the Kuramoto-Sivashinsky equation. GRU and RC reproduce the positive exponents and Kaplan-Yorke dimension accurately, although neither captures the two zero exponents directly.
- Method: The KS Lyapunov spectrum is computed from the reference equations and surrogate hidden-state Jacobians using QR or Gram-Schmidt orthonormalization.A GRU is used instead of an LSTM because the LSTM has two coupled hidden states, making spectrum computation more involved and expensive.
- Results: The largest true Lyapunov exponent is Λ1 ≈0.08844; RC estimates Λ1 ≈0.08378 with 5.3% relative error, while GRU estimates Λ1 ≈0.0849 with ≈4% error.
- Results: Both surrogate models identify the correct Kaplan-Yorke dimension, KY ≈15.The dimension is defined using the largest index for which the cumulative Lyapunov-exponent sum remains positive.
- Results: GRU and RC match the positive KS Lyapunov exponents well, but neither surrogate captures the two zero exponents Λ7 and Λ8.After augmenting the spectra with these two zero exponents, agreement with the true spectrum becomes excellent.
- Results: After zero-exponent augmentation, the mean absolute errors are 0.012 for RC and 0.008 for GRU, while mean relative errors are 0.23 and 0.22, respectively.
7. Conclusions
The study compares RNNs trained by BPTT with RC for forecasting chaotic systems, including reduced-order and high-dimensional settings. Performance, scalability, training cost, Lyapunov-spectrum accuracy, and unresolved limitations differ across methods and representations.
- Comparative framework: The study compares LSTM, GRU, Unitary RNNs, and RC for forecasting chaotic systems while evaluating temporal dependencies, scalability, overfitting, and applicability to spatiotemporal dynamics.RC trains only output weights by least-squares regression, whereas the gated RNNs are trained with BPTT.
- Reduced-order forecasting: LSTM and GRU outperform RC on reduced-order Lorenz-96 forecasting and usually remain stable during iterative prediction, unlike RC and Unitary forecasts.The reduced-order observable retains the most energetic SVD modes under forcing regimes F={8,10}.
- Reduced-order forecasting: Large RC models overfit more easily than LSTM and GRU models because validation-based early stopping and regularization are not directly applicable to RC.RC validation instead involves tuning the Tikhonov regularization parameter.
- Parallel architectures: Parallelization exploiting local state interactions addresses scalability problems in high-dimensional systems for both RNNs and RC.The approach is evaluated on Lorenz-96 and the Kuramoto-Sivashinsky equation.
- Parallel architectures: Parallel LSTM and GRU achieve predictive performance similar to parallel RC on Lorenz-96 and Kuramoto-Sivashinsky, while RC training is significantly faster.RC uses 1000–6000 reservoir nodes per member, compared with a few hundred nodes for parallel GRU/LSTM.
- Lyapunov-spectrum analysis: Both RC and GRU accurately capture positive Lyapunov exponents of the Kuramoto-Sivashinsky equation, but neither reproduces the two zero exponents Λ7 and Λ8.After discarding those two exponents, the methods have comparable relative and absolute spectrum errors.
- Limitations and future work: The study identifies unresolved limitations involving zero Lyapunov exponents, RC memory and eigenvalue-decomposition costs, model invertibility, and initialization of complex gated architectures.The authors also propose hybrid RC–LSTM training and RC-informed BPTT initialization as future directions.
- Overall conclusion: The authors conclude that RNN-based surrogate models can efficiently represent high-dimensional chaotic systems and can be parallelized to alleviate scaling problems.They describe this as a promising research subject requiring further analysis.
8. Data and Code
The paper states that its code and data will be released to support reproducibility, with software implemented in Python using TensorFlow and PyTorch.
- Reproducibility: Code and data will be available upon publication through the authors’ GitHub repository to assist reproducibility.The repository link is https://github.com/pvlachas/RNN.
- Implementation: The software uses Python, TensorFlow for automatic differentiation, and PyTorch for neural-network architecture design.
CRediT authorship contribution statement
The contribution statement assigns work across conceptualization, methodology, software, data curation, visualization, writing, supervision, and review.
- Contributions: P. R. Vlachas led conceptualization, methodology, software, data curation, visualization, and the original draft.
- Contributions: J. Pathak contributed conceptualization, methodology, software, data curation, and review and editing.
- Contributions: B. R. Hunt, T. P. Sapsis, M. Girvan, and E. Ott contributed conceptualization and supervision, with some also contributing to review and editing.
A. Memory Efficient Implementation of RC Training
The memory-efficient RC implementation accumulates batched hidden-state and target products instead of storing the full training matrices, while comparing alternative linear solvers.
- Matrix construction: Teacher forcing produces hidden-state matrix H and one-step-shifted target matrix Y from d_N training time steps.H has dimensions d_N×d_h, while Y has dimensions d_N×d_o.
- Linear solve: The Moore–Penrose pseudoinverse is used as a classical solver for the RC output-weight system.
- Memory reduction: When d_N is about 10^5 and much larger than d_h, the implementation forms H^T H and Y^T H with dimensions d_h×d_h and d_o×d_h.The accumulated products reduce the memory burden relative to retaining the full matrices.
- Batched accumulation: Time batches of d_n≪d_N steps form H_b and Y_b, whose contributions are added to the accumulated products after each batch.
- Solver comparison: Among the tested RC solvers, the pseudoinverse produced the most accurate Lorenz-96 results and was treated as a model hyperparameter alongside alternative routines.The alternatives included regularized least squares with LSQR and stochastic gradient descent.
B. Regularizing Training with Noise
Adding modest noise to training data improves RC forecasting efficiency and reduces divergent iterative predictions, while the section also describes the SVD-based reduction and Lyapunov-spectrum workflow.
- Noise regularization: 5–10 ‰ noise improves RC generalization and forecasting efficiency on both training and testing data.The same regularization also significantly alleviates divergence during iterative forecasts.
- Dimensionality reduction: SVD reduces dataset dimensionality by projecting centered state data onto modes ordered by explained variance.The squared singular values quantify each mode’s energy contribution.
- Dimensionality reduction: The Lorenz-96 reduced-order observable uses the 35 modes with the highest cumulative energy contribution.
- RNN training: GRU and LSTM BPTT models share a common observable-to-hidden-to-output architecture while differing in their internal mappings.
- Lyapunov-spectrum calculation: Lyapunov exponents are estimated from hidden-state Jacobian products along surrogate-RNN trajectories, with periodic orthonormalization to prevent numerical divergence.The algorithm warms up the hidden state using true data before evolving deviation vectors and checking convergence.
- Model evaluation: The study evaluates RC, GRU, LSTM, and Unitary models across Lorenz-96 settings using separately specified hyperparameters.
F. Additional Results - Lorenz-96 - Divergence of Unitary and RC RNNs
The Lorenz-96 divergence analysis compares iterative forecast failures across architectures, showing different stability patterns for fully observed and reduced-order data.
- Divergence patterns: Fully observed Lorenz-96 models often produce zero divergent predictions for RC, GRU, and LSTM across both forcing regimes.
- Divergence patterns: Most Unitary hyperparameter sets diverge from the attractor during iterative forecasting in both reduced-order and full-order scenarios.
- Divergence patterns: Reduced-order divergence is more prominent for RC and Unitary networks than for GRU and LSTM models.
- Forecast comparison: Figure 20 compares spatiotemporal forecasts, NRSE contours, and component-average NRSE across architectures for reduced-order and full-state observables.
G. Additional Results - Lorenz-96 - 퐹= 10
Additional Lorenz-96 experiments at F=10 examine forecast behavior across observable dimensions and provide contour-based comparisons of model errors.
- Additional results: F=10 results agree with the main conclusions established for F=8.
- Forecast visualization: Figure 22 compares GRU, LSTM, RC, and Unitary forecasts against true Lorenz-96 evolution using NRSE contours and component-average NRSE.
- Forecast visualization: For observable dimension d_o=40, forecasts diverge more slowly than for the reduced-order observable.
- Temporal dependencies: In reduced-order Lorenz-96 forecasting, larger truncated sequence and prediction lengths improve the efficiency of GRU and LSTM models.
H. Temporal Dependencies and Backpropagation
The study examines temporal-dependency settings for BPTT models and tests generalization in parallel Lorenz-96 and Kuramoto-Sivashinsky forecasting.
- BPTT parameters: κ1 denotes truncated backpropagation length, while κ2 denotes the number of future steps used for loss computation and gradient propagation.The hyperparameter study varies κ1∈{8,16} and κ2∈{1,8}.
- Temporal dependencies: Large κ1 and κ2 are pivotal for high forecasting efficiency in reduced-order observables, indicating relevant temporal correlations.
- Parallel generalization: Parallel models show similar VPT on training and testing data for both Lorenz-96 and Kuramoto-Sivashinsky forecasts.This agreement indicates no observed overfitting in the parallel setting.