Backpropagation Algorithms and Reservoir Computing in Recurrent Neural Networks for the Forecasting of Complex Spatiotemporal Dynamics

1. Introduction

• The introduction frames complex-system applications as a growing opportunity for machine learning, spanning weather, climate, energy, medicine, ocean dynamics, and turbulent flows.It notes that machine-learning deployment on complex systems remains in its infancy.Source: §1 · p. 1
• Complex systems’ interacting spatiotemporal scales and chaotic dynamics challenge classical prediction methods, while data-driven RNNs can help when governing equations are unavailable or costly to evaluate.RNN hidden states encode information about the history of system dynamics.Source: §1 · p. 2
• RC has been used for high-dimensional chaotic systems, including parallel ensembles that exploit local state-space interactions to scale forecasting (Pathak et al., 2018a).In that approach, each network forecasts a group of modes and receives interacting modes as inputs.Source: §1 · p. 2
• LSTM gates alleviate vanishing gradients by selectively remembering and forgetting information, and prior work applied LSTMs to physical-system forecasting.RNNs had also been used for physical dynamics before gated architectures, though learning could suffer vanishing or exploding gradients.Source: §1 · p. 2
• The study compares RNNs trained with backpropagation and RC on fully observed and reduced-order systems, assessing forecasting, efficiency, Lyapunov exponents, and frequency spectra.It also considers Unitary, GRU, and LSTM architectures and parallel RNN ensembles for high-dimensional forecasting.Source: §1 · pp. 2–3
• The comparison addresses model-agnostic forecasting when system equations are unknown, including reduced-order observations that omit some degrees of freedom.The study also examines parallel prediction of full-state dynamics and discusses model implementation, RAM requirements, training time, and limitations.Source: §1 · pp. 2–3

2. Methods - Sequence Modeling

• RC and BPTT-trained networks share a hidden-state architecture that forecasts the next observable from the current input and encoded history.The models differ in their recurrent mapping and in how their parameters are learned.Source: §2 · p. 3

2.1. Long Short-Term Memory

• LSTMs use forget, input, and output gates with a cell state to control information flow and alleviate vanishing gradients.The hidden-state dimension determines the cell’s capacity to encode history; the forget-gate bias is initialized to one to accelerate training.Source: §2.1 · p. 4

2.2. Gated Recurrent Unit

• GRUs use update and reset gates but omit the LSTM output gate, yielding fewer parameters while retaining a similar gating mechanism.The paper notes comparable LSTM and GRU performance on polyphonic music and speech datasets.Source: §2.2 · p. 4

2.3. Unitary Evolution

• Unitary RNNs constrain the recurrent weights to complex unitary matrices with spectral radius one, aiming to prevent vanishing gradients without gating.The recurrent matrix is optimized on the subspace of complex unitary matrices.Source: §2.3 · p. 5

2.4. Back-Propagation Through Time

• BPTT computes loss gradients by applying the chain rule backward through time and updates network parameters with a first-order optimizer.The study uses BPTT for LSTM, GRU, and Unitary RNNs; its tunable parameters include the forward-pass accumulation length κ1 and gradient truncation length κ2, which must capture relevant temporal dependencies.Source: §2.4 · p. 5
• The study trains only stateful BPTT networks, carrying hidden states between consecutive batches and teacher-forcing between updates for κ3 = κ2 + κ1 − 1 steps.Unlike stateless models, stateful models can propagate dependencies across batches; teacher forcing supplies ground-truth inputs without back-propagation to prepare the next batch’s hidden state.Source: §2.4 · p. 6, Figure 2 · p. 6
• Under multi-round learning-rate reduction, the study found little result sensitivity to Adam’s β1 and β2 settings, but did not explore them further because of limited compute time.The authors report sensitivity when using a single round and treating the learning rate as an additional hyperparameter.Source: §2.4 · p. 7

2.5. Reservoir Computing

• RC fixes randomly selected recurrent weights to create a reservoir state and trains the hidden-to-output mapping, with the reservoir scaled to satisfy the echo state property.The method uses regularized least-squares output training and adds Gaussian noise to training data; the noise level is tuned, and squared hidden-state features are used for half the reservoir nodes.Source: §2.5 · p. 7

Comparison Metrics

• NRMSE compares forecasts with targets using each state component’s temporal standard deviation, and its time evolution is averaged over 100 randomly sampled attractor initial conditions.The state-space average is used in the NRMSE definition.Source: §3 · p. 8
• Lyapunov exponents quantify average exponential separation of nearby trajectories; their ordered spectrum typically has as many values as the state-space dimension.The spectrum is evaluated in Section 6.Source: §3 · p. 8
• The maximal Lyapunov exponent (MLE) is the largest exponent, governs separation for randomly oriented nearby trajectories, and defines the Lyapunov time as its inverse.A positive MLE characterizes chaotic motion of a bounded trajectory, and Lyapunov time provides a scale for judging prediction-error growth.Source: §3 · p. 8
• Valid prediction time is the longest forecast duration with NRMSE below 0.5, normalized by the system’s MLE.It serves as a single metric of predictive performance.Source: §3 · p. 8
• Mean power spectral density across all state components is used to assess how well models capture the dynamical system’s long-term statistics.The power spectrum is the complex Fourier spectrum of each component’s state evolution.Source: §3 · p. 8

4. Forecasting Reduced Order Observable Dynamics in the Lorenz-96

• RNN hidden states can capture delayed-observation embedding information without manually tuning the time lag and embedding dimension.Reduced observables are harder to forecast because information loss makes their dynamics more irregular, and even small initial errors grow exponentially in chaotic systems.Source: §4 · pp. 8–9

4.1. Lorenz-96 Model

• The Lorenz-96 experiments use J=40 with F=8 or F=10; the 35 leading SVD modes explain approximately 98% of system energy, while 40 modes provide the full state.Forecasting is evaluated on reduced observables because incomplete state information requires networks to capture temporal dependencies.Source: §4.1 · p. 9, §4 · p. 9

4.2. Results on the Lorenz-96 Model

• RC leads with full-state information, whereas GRU and LSTM are more stable and outperform RC in reduced-order Lorenz-96 forecasts.For F=10 reduced-order data, GRU and LSTM forecast up to 2 Lyapunov times before NRMSE reaches 1, versus 1 Lyapunov time for RC and Unitary RNNs; for the stricter NRMSE threshold 0.5, all models are limited to one Lyapunov time on reduced-order data and two on full-state data.Source: §4.2 · pp. 9–10, Figure 3 · p. 10
• For F=8, full-state RC reaches a maximum VPT of ≈2.31, compared with ≈1.34 for GRU, ≈0.97 for LSTM, and ≈0.58 for Unitary RNNs.For F=10, full-state RC has the highest maximum VPT at ≈2.35; LSTM reaches ≈1.73 and GRU ≈1.59.Source: §4.2 · pp. 10–11
• With reduced-order data, GRU has the highest maximum VPT for F=8 at ≈0.98, while GRU and LSTM lead for F=10 at ≈1.22 and ≈1.17, respectively.Reduced-order RC reaches maximum VPTs of ≈0.55 for F=8 and ≈0.60 for F=10, below both gated networks in each regime; high-VPT hyperparameter sets are rarer than in the full-state setting.Source: §4.2 · p. 11
• All models match the true power spectra in full-state forecasts, but reduced-order RC fails to match the statistics that GRU and LSTM capture.Unitary RNN forecasts generally diverge from the attractor and are omitted from the long-term statistics comparison.Source: §4.2 · pp. 11–12, Figure 5 · p. 12
• For reduced-order forecasts, BPTT-trained GRU and LSTM achieve higher VPT than RC despite RC receiving an order of magnitude more memory.RC's maximum reservoir size may not capture high-dimensional dynamics with reduced-order information and non-local interactions; its training time also grows cubically with reservoir size for the Moore-Penrose solve.Source: §4.2 · p. 12
• GRU and LSTM have lower reduced-order generalization errors than RC, in part because validation-based early stopping guards against overfitting.The tested Tikhonov regularization values did not eliminate RC's higher generalization error; in full-state forecasting, RC instead achieves superior accuracy and generalization, aided by Gaussian-noise regularization.Source: §4.2 · p. 13, Figure 6 · p. 13

5. Parallel Forecasting Leveraging Local Interactions

• For high-dimensional systems, model size, training time, and RAM requirements limit single-network forecasting; local interactions enable parallel networks to forecast local regions.When dynamics are translation invariant, the local networks can be identical, allowing training a single network and copying its weights.Source: §5 · p. 15

5.1. Parallel Architecture

• Each of N_g RNNs predicts G neighboring state elements from those elements plus I neighboring inputs on each side; networks train independently and exchange predictions for long-term forecasting.The implementation uses MPI for communication, while translation-invariant systems can use one trained RNN with copied weights.Source: §5.1 · pp. 15–16

5.2. Results on the Lorenz-96

• Exploiting local interactions raised Lorenz-96 VPT from ≈1.4 for nonparallelized models to ≈3.9 with parallel networks.BPTT networks achieved predictions comparable to RC with 100–500 hidden nodes versus 6000–12000 for RC; RC with 3000 or 6000 nodes had slightly lower VPT than GRU and LSTM but significantly lower training times. RC with 12000 nodes required more than 3 GB per rank, and all models reproduced the system’s power spectrum.Source: §5.2 · p. 16, Figure 9 · p. 17, Figure 10 · p. 17

5.3. Kuramoto-Sivashinsky

• The Kuramoto–Sivashinsky benchmark used a full-state observable of d_o=512 grid nodes.The discretized trajectory was split into 100,000 training and 100,000 testing samples after truncating initial transients; parallel models used G=8, I=8, and N_g=64, with 24 inputs per network.Source: §5.3 · p. 17, §5.4 · p. 18

5.4. Results on the Kuramoto-Sivashinsky Equation

• On Kuramoto–Sivashinsky, VPT was ≈4 for LSTM, ≈3.5 for GRU, ≈3.2 for RC with d_h=500, ≈3.9 for RC with d_h=1000, and ≈4.8 for RC with d_h=3000.The d_h=3000 RC slightly outperformed GRU and LSTM while using a similar amount of RAM and much less training time; the authors attribute GRU’s lower performance partly to the risk of poor local minima across parallel networks.Source: §5.4 · pp. 18–19, Figure 12 · p. 19
• For Kuramoto–Sivashinsky at state dimension d₀=512, averaged over 100 initial conditions, an LSTM with 80 hidden units matched a 1000-hidden-unit RC in NRMSE, and all models captured system statistics.Source: §5.4 · p. 19, Figure 13 · p. 19

• The KS reference spectrum used Fourier-space ETDRK4 and QR; surrogate hidden-state Jacobians used Gram–Schmidt, and GRU was chosen over the more costly two-state LSTM.The GRU spectrum required an iterative rollout of N=10^4 time steps to converge.Source: §6 · p. 21
• Λ1≈0.0849 from GRU (≈4% error) improved on RC’s Λ1≈0.08378 (5.3% error) against the true Λ1≈0.08844; both surrogates identified KY≈15.The comparison used a GRU with 2,000 hidden units and an RC with 9,000 nodes.Source: §6 · p. 21, Figure 15 · p. 22
• Neither GRU nor RC captured zero exponents Λ7 and Λ8; adding them to each surrogate spectrum produced excellent agreement with the true spectrum.The exponents Λ7 and Λ8 correspond to the zero crossings in the KS spectrum.Source: §6 · p. 21, Figure 15 · p. 22
• 0.008 mean absolute error for GRU versus 0.012 for RC, and 0.22 mean relative error for GRU versus 0.23 for RC, after augmenting with zero exponents.Source: §6 · p. 21, Figure 16 · p. 22

Conclusions

• The study benchmarks LSTM, GRU, and Unitary cells trained with BPTT against RC with random recurrent weights and least-squares-trained output weights.These RNN cells are considered as approaches to mitigating vanishing and exploding gradients in BPTT.Source: §7 · p. 23
• For reduced-order Lorenz-96, LSTM and GRU forecast more stably and accurately than RC and Unitary, while full-state RC outperforms all other models and reproduces the frequency spectrum.The reduced-order comparisons use forcing regimes F = 8 and 10; LSTM and GRU have similar memory requirements to RC but take longer to train.Source: §7 · p. 23
• Parallel GRU, LSTM, and RC reproduce long-term attractor climates and power spectra for Lorenz-96 and Kuramoto-Sivashinsky; LSTM and GRU match RC’s predictive performance with fewer nodes but longer training.Memory requirements are comparable, but RC needs reservoirs of 1000–6000 nodes per member to match parallel GRU/LSTM models with a few hundred nodes.Source: §7 · p. 23
• Both GRU and RC capture positive Kuramoto-Sivashinsky Lyapunov exponents but miss zero exponents Λ7 and Λ8; excluding those yields comparable spectrum errors.The comparison considers relative and absolute error of the Lyapunov spectrum.Source: §7 · p. 23
• Future work includes investigating memory capacity, data requirements, RC-based RNN initialization, and hybrid RC–LSTM training to improve predictive performance.In the proposed hybrid, backpropagation would optimize the RC output and LSTM weights, with the LSTM correcting RC errors; training cost and suitable echo-state-property initialization for gated cells remain challenges.Source: §7 · pp. 23–24
• Batched RC training reduces memory requirements, but matrix inversion and sparse-weight eigenvalue decomposition still pose memory and numerical challenges.The authors identify further work on alleviating the memory burden of matrix inversion and numerical problems with eigenvalue decomposition.Source: §7 · p. 23
• The study leaves model invertibility unexamined; noninvertible RNNs may produce spurious dynamics, while invertible RNNs are proposed as a possible alternative.The cited concern applies when forecasting full-state dynamics.Source: §7 · p. 24

8. Data and Code

• Code and data will be available upon publication at the project’s GitHub link to assist reproducibility; Python, TensorFlow, and PyTorch support automatic differentiation and neural-network design.The passage cites Abadi et al. (2016) for TensorFlow and Paszke et al. (2017) for PyTorch.Source: §8 · p. 24

A. Memory Efficient Implementation of RC Training

• Time-batched accumulation computes the quantities needed for the RC output weights without retaining the full training history in memory.The full hidden-state and shifted-target matrices span dN time steps; batches of dn ≪ dN update the accumulated matrices after each batch.Source: A. Memory Efficient Implementation of… · p. 28
• The output weights are obtained by solving a linear system, for which the classical approach is the Moore–Penrose pseudo-inverse with Tikhonov regularization.The section notes that dN is of order 10^5 and much larger than dh, motivating reduced-memory computation of the associated matrices.Source: A. Memory Efficient Implementation of… · p. 28
• In Lorenz-96 tests, the pseudo-inverse was more accurate than iterative regularized least squares or stochastic gradient descent, so it was used for Kuramoto–Sivashinsky to save computational resources.Source: A. Memory Efficient Implementation of… · p. 28

B. Regularizing Training with Noise

• RC models trained with 5–10‰ noise achieve better generalization and forecasting efficiency in both training and testing data.Figure 17 plots testing VPT against training VPT for RC and GRU models trained with different noise levels.Source: B. Regularizing Training with Noise · p. 28, Figure 17 · p. 29
• Added noise significantly alleviates divergent iterative forecasts, whereas it has little impact on GRU performance.Source: B. Regularizing Training with Noise · p. 28, Figure 17 · p. 29

C. Dimensionality Reduction with Singular Value Decomposition

• SVD dimensionality reduction begins by subtracting the temporal mean from the stacked time-evolution data matrix.Source: C. Dimensionality Reduction with… · p. 30
• SVD modes are ordered by descending singular value, and each mode’s squared singular value represents its share of data energy.The percentage of total energy for a mode is calculated by dividing its squared singular value by the sum of all squared singular values.Source: C. Dimensionality Reduction with… · p. 30
• To reduce dimension, the method retains the first r_dim high-energy modes, discards the remaining modes, and projects each state vector onto the retained modes.Source: C. Dimensionality Reduction with… · p. 30
• The Lorenz-96 reduced-order observable uses 35 modes, with its cumulative energy distribution and state and mode evolution illustrated in Figure 18.Source: C. Dimensionality Reduction with… · p. 30, Figure 18 · p. 30

D. Calculation of Lyapunov Spectrum

• The trained GRU surrogate is used to estimate the full Kuramoto-Sivashinsky Lyapunov spectrum, while true exponents are computed in Fourier space using ETDRK4 and QR decomposition.Source: D. Calculation of Lyapunov Spectrum · p. 31
• The method forms the hidden-state Jacobian along the trajectory, evolves its product, and periodically orthonormalizes deviation vectors while tracking convergence of estimated exponents.Orthonormalization uses Gram-Schmidt; the algorithm accepts a warm-up series, number of exponents, unroll limit, normalization interval, and termination threshold.Source: D. Calculation of Lyapunov Spectrum · p. 31
• The procedure applies directly to one-hidden-state RC and GRU models; adapting it to LSTM is left for future research.Source: D. Calculation of Lyapunov Spectrum · p. 31

E. Model Hyperparameters

• Lorenz-96 hyperparameters are reported for RC, GRU/LSTM, and Unitary models at d_o=35 and d_o=40 in PCA mode.Source: E. Model Hyperparameters · p. 33, Table 2 · p. 33, Table 3 · p. 33
• Parallel RC, GRU/LSTM, and Unitary hyperparameters are tabulated for Lorenz-96 and Kuramoto-Sivashinsky; Adam was used to reduce computational burden.Source: E. Model Hyperparameters · p. 33, Table 5 · p. 34, Table 6 · p. 34
• Most Unitary hyperparameter sets produce divergent iterative forecasts across both Lorenz-96 state sizes and forcing regimes.The comparison covers reduced order d_o=35 and full order d_o=40 for F=8 and F=10.Source: Figure 19 · p. 36
• RC divergence is negligible for full-state Lorenz-96 forecasts but more prominent in reduced order forecasts than for GRU and LSTM.The figure describes results across 100 test initial conditions and all hyperparameter sets.Source: Figure 19 · p. 36
• Stable reduced-order forecasting hyperparameters are easier to identify for GRU and LSTM than for RC and Unitary networks.Source: Figure 19 · p. 36

F. Additional Results - Lorenz-96 - Divergence of Unitary and RC RNNs

• Across both forcing regimes, reduced-order divergence is more prominent in RC and Unitary networks, making hyperparameter tuning for these models more cumbersome than for GRU and LSTM.The comparison covers forcing regimes F∈{8, 10} and 100 initial conditions across tested hyperparameter sets.Source: F. Additional Results - Lorenz-96 -… · p. 36
• In a reduced-order test forecast, RC and Unitary predictions diverge after approximately two Lyapunov times; full-state forecasts diverge more slowly because the dynamics are deterministic.The example uses Lorenz-96 with F=8; reduced-order forecasting is challenging because of initial-condition sensitivity and incomplete state information.Source: F. Additional Results - Lorenz-96 -… · p. 36, Figure 20 · p. 37
• With full-state data, RC achieves superior forecasting accuracy, while many RC, GRU, and LSTM hyperparameter sets avoid divergence; in reduced-order forecasts, GRU and LSTM are stable and reproduce attractor statistics, unlike frequently diverging RC models.Unitary networks also commonly diverge from the attractor in iterative prediction.Source: F. Additional Results - Lorenz-96 -… · p. 36, Figure 20 · p. 37

G. Additional Results - Lorenz-96 - 퐹= 10

• At F=10 and d_o=35, LSTM and GRU remain stable and reproduce long-term attractor statistics, while even RC with d_h=9000 may diverge and miss the attractor climate.Very large RC networks can achieve accurate short-term predictions, but require high memory; at d_o=40, forecasts diverge more slowly because the dynamics are deterministic.Source: G. Additional Results - Lorenz-96 - 퐹=… · p. 38, Figure 22 · p. 39
• At d_o=35, forecasting is challenging because of sensitivity to initial conditions and incomplete state information; at d_o=40, forecasts diverge more slowly because the dynamics are deterministic.The reduced-order observable requires capturing temporal dependencies, unlike the full-state observable.Source: Figure 22 · p. 39
• For reduced-order d_o=35, a large truncated BPTT parameter κ1 is vital for capturing temporal dependencies and forecasting efficiently, whereas a small horizon suffices at d_o=40.Figure 23 reports mean VPT with minimum-to-maximum error bars for stateful LSTM and GRU models trained with different κ1 and κ2 values.Source: Figure 23 · p. 40

H. Temporal Dependencies and Backpropagation

• BPTT training of GRU and LSTM models tunes κ1, the truncated backpropagation or sequence length, and κ2, the number of future steps used for loss and gradient computation per batch.Source: H. Temporal Dependencies and… · p. 40
• The hyperparameter study varied κ1 ∈ {8, 16} and κ2 ∈ {1, 8}.Source: H. Temporal Dependencies and… · p. 40
• For reduced-order Lorenz-96, large sequence length κ1 and prediction length κ2 are pivotal for high forecasting efficiency.The results are averaged over 100 test initial conditions for F ∈ {8, 10}.Source: H. Temporal Dependencies and… · p. 40
• The reduced-order result indicates temporal correlations that models with smaller horizons cannot easily capture.Source: H. Temporal Dependencies and… · p. 40
• In the full-order scenario, models with smaller κ1 perform reasonably well, consistent with less need to capture temporal correlations when the full system state is available.Source: H. Temporal Dependencies and… · p. 40

I. Over-fitting of Parallel Models

• In both Lorenz-96 and Kuramoto-Sivashinsky parallel models, test VPT is very close to training VPT, with no observed overfitting.The authors interpret this as evidence of good generalization.Source: Figure 24 · p. 41, §I · p. 41