On the emergence of numerical instabilities in Next Generation Reservoir Computing

TL;DR

This work analyzes how the numerical conditioning of NGRC’s feature matrix, built from polynomial evaluations on time‑delay coordinates, governs long‑term forecasting stability. By exposing the Vandermonde‑ and Hankel‑like structure of Ψ and leveraging ergodic theory, it links hyperparameters (polynomial degree p, delay dimension k, time lag τ, and training length N_train) to conditioning and numerical sensitivity. The study shows κ(Ψ) grows exponentially with p and can blow up when k>1 or τ is small, but increasing τ or using SVD for training can yield accurate forecasts without regularization, guiding practical hyperparameter choices. In addition to Lorenz dynamics, results on the Double‑Scroll circuit and partial measurements demonstrate broad applicability of conditioning‑aware NGRC design, with concrete recommendations such as preferring SVD and tuning τ to balance stability and fidelity. The findings offer a path to reduced computational cost in NGRC by avoiding unnecessary regularization while maintaining dynamical stability, with implications for other data‑driven linear readout models and reservoir computing variants.

Abstract

Next Generation Reservoir Computing (NGRC) is a low-cost machine learning method for forecasting chaotic time series from data. Computational efficiency is crucial for scalable reservoir computing, requiring better strategies to reduce training cost. In this work, we uncover a connection between the numerical conditioning of the NGRC feature matrix -- formed by polynomial evaluations on time-delay coordinates -- and the long-term NGRC dynamics. We show that NGRC can be trained without regularization, reducing computational time. Our contributions are twofold. First, merging tools from numerical linear algebra and ergodic theory of dynamical systems, we systematically study how the feature matrix conditioning varies across hyperparameters. We demonstrate that the NGRC feature matrix tends to be ill-conditioned for short time lags, high-degree polynomials, and short length of training data. Second, we evaluate the impact of different numerical algorithms (Cholesky, singular value decomposition (SVD), and lower-upper (LU) decomposition) for solving the regularized least-squares problem. Our results reveal that SVD-based training achieves accurate forecasts without regularization, being preferable when compared against the other algorithms.

Figures (14)

• Figure 1: NGRC model accurately reproduces the Lorenz attractor. The top panel shows the NGRC reconstruction (red) of the Lorenz attractor (blue). The bottom panel shows the NGRC model's performance over the training and testing phases. The horizontal axis shows time in Lyapunov time $\frac{1}{\Lambda}$ where $\Lambda = 0.9056$ is the maximum Lyapunov exponent of the Lorenz system. The parameters are $h = 0.01$, delay dimension $k = 1$, maximum degree $p = 2$, time lag $\tau = 1$, $N_{\mathrm{train}} = 500$, $N_{\mathrm{test}} = 10000$ (but a smaller time window is shown), and regularizer parameter $\beta = 0$.
• Figure 2: NGRC model captures topological and statistical features of the original dynamics. The left panel displays the Poincaré return map of the successive local maxima of the original dynamics (in blue) and the reconstructed one (in red). The right panel depicts the power spectrum density of the original dynamics (in blue) and the reconstructed one (in red). Both curves mostly overlap for all frequencies, indicating that the NGRC model captures correctly the long-term statistics on the attractor. The parameters are the same used for \ref{['fig:NGRC_k_1_dt_001']}.
• Figure 3: Regularization reduces the difference among numerical algorithms during NGRC training. The top panel shows the maximum closeness of fit for Cholesky (dark gray), SVD (gray), and LU (light gray) solvers plotted together with the regularized condition number of $\Psi$ (red) with respect to the regularizer parameter $\beta$. All numerical algorithms attain similar fitting accuracy at $\beta = 10^{-10}$. While the fitting performance gets worse as the regularizer parameter is increased, the bottom panel shows that the maximum pairwise difference among the solutions decays monotonically. The dots are the median, and bars are $75\%$ and $25\%$ quantiles over 50 different initial conditions. The parameters are $h = 0.01$, delay dimension $k = 2$, time lag $\tau = 1$, the maximum degree $p = 2$, and $N_{\mathrm{train}} = 5000$.
• Figure 4: Regularization also reduces the difference among numerical algorithms during NGRC testing. Comparison between three different numerical algorithms: Cholesky, SVD, and LU as the regularizer parameter $\beta$ is increased under the testing phase. From top to bottom, panels show box plots with respect to bounded models over 50 different initial conditions: valid prediction time (VPT), the distance between the induced successive local maxima maps, and the error $E$ between the power spectrum density, respectively. The colored hashed areas correspond to: absent/small (green), large (purple) regularization, and level of regularization in which all algorithms perform similarly within statistical confidence (orange). Blank spaces correspond to unbounded NGRC models over all initial conditions.. The parameters are $h = 0.01$, delay dimension $k = 2$, time lag $\tau = 1$, the maximum degree $p = 2$, and $N_{\mathrm{train}} = 5000$.
• Figure 5: Exponential growth with respect to the maximum degree. Box plot of the condition number $\kappa(\hat{\Psi})$ over 25 distinct initial conditions for increasing maximum degree $p$. The exponential growth is confirmed by the exponential curve (black dashed line) plotted for reference. The numerical integrator is the explicit forward Euler method with step size $h = 0.01$, and the number of training data points $N_{\mathrm{train}} = 5000$.

Theorems & Definitions (3)

• proof
• Example 4.2: $x$-coordinate case
• Remark 4.3: Scaled regularizer parameter $\beta$ by the length of training data $N_{\mathrm{train}}$