1 Introduction↩︎

Reservoir computing (RC) has seen broad success in the domain of time series forecasting [1]–[3]. Part of its appeal stems from the fact that the classical echo-state-network (ESN) architecture is realizable in various physical systems [4]. Next-generation reservoir computing (NG-RC) [5] has been introduced as an alternative to ESN-type RC: it can be more efficiently trained, entails better control over hyperparameters, and is easier to interpret due to an explicit feature library. Recent applications include, in particular, complex spatio-temporal and non-stationary chaotic systems [6]–[11]. Historically, NG-RC has been obtained by approximating ESN-type RCs as a nonlinear vector autoregressive process (NVAR) [12]. NVARs have a long history in time series forecasting [13]–[15]. In the context of static data processing, NG-RCs are also known as extreme learning machines [16]–[19], which are closely related to random feature methods [20], i.e., single-layer feed-forward networks with random weights. Evolving a hand-chosen dictionary (e.g., instantaneous and time-delay coordinates, and their low-order products) without necessarily projecting back to the physical state is practiced in the Extended Dynamic Mode Decomposition (EDMD) method [21], [22]. Since the training objective enforces that the dictionary evolves within itself, EDMD renders a finite-dimensional approximation to the Koopman operator [23]–[29]. The Koopman operator approach constitutes an overarching theoretical framework for modeling complex dynamical systems and has found numerous applications [30]–[34]. Choosing a monomial basis for the Koopman operator leads to the Carleman linearization method [35].

Characterizing the representations that a model learns from its training data is a central issue of system identification and inverse problem theory [15], [36]–[48]. Integrable dynamical systems can in principle be described by many equivalent surrogate models as a consequence of the underlying symmetries [49], [50]. In teacher–student setups with well-specified (bias-free) models, the student can learn the data-generating process up to a similarity transform [51], resulting in near perfect generalization [52]–[55]. Physics-informed learning supports generalization by incorporating additional structure into a model [56], [57].

A straightforward approach to system identification is provided by the SINDy method [38] and its refinements [58]–[60]. However, these methods typically require approximating the time derivative of the state variables, which can limit robustness for noisy data [58]. We instead focus on the discrete-time setting, where the dynamics is specified by the flow map. Moreover, the presence of delays can obscure the direct relationship between the readout coefficients and the flow map [5], [61], [62]. Various studies have discussed this in the context of numerical instability [63], noise robustness [64], and empirical time-stepping formulations [62], [65]. For stationary signals, viable surrogate models are harmonic time series [24], [66], [67]. These can be exactly represented by linear delay models (linear VARs) [68]–[70]. If the underlying system is linear or follows a limit cycle (periodic motion), the approximation error is mainly limited by the number of frequencies included in the model. Nonlinearities in the underlying system are reflected in coupled phases and nontrivial polyspectra (i.e., Fourier-transformed higher-order cumulants of the process) of the harmonic surrogate.

The discrepancy between one-step training error and autonomous forecast accuracy is well known in system identification as the distinction between prediction-error and simulation-error fitting [71]–[73], and in sequence modeling as the teacher-forcing or exposure-bias problem [74]–[76]. One-step-ahead minimization can select models with poor autonomous rollout behavior, especially for highly sampled data where adjacent lags are strongly correlated [62], [77]. In Koopman-based approaches such as EDMD, prediction accuracy is known to depend on how close the dictionary is to being Koopman-invariant, and that residual error alone does not encode the quality of the feature space [78]–[80].

In the present work, we consider the problem of learning flow maps of fully observed Markovian dynamical systems using NVAR/NG-RC models with instantaneous and delay-augmented feature libraries. We formulate a training error scaling theory based on the Lie-series expansion of the flow map and specialize it to polynomial and Fourier feature libraries. Closely related to the present work are also error analyses of the Koopman and Carleman linearization approaches [26], [46], [79]–[86], including polynomial and delay-augmented embeddings [87], [88], and of learned multi-step methods [89], [90]. The present work hopes to contribute a useful addition to these studies in the context of NVAR/NG-RC modeling.

2 NVAR/NG-RC Model↩︎

In the Koopman/EDMD formulation [22], [27], [29], the next-step prediction problem is given by \[\boldsymbol{\psi}(t+\Delta t) = K \boldsymbol{\psi}(t), \label{eq95nextstep95task}\tag{1}\] where the feature vector \(\boldsymbol{\psi}(t)=\mathbf{G}(\mathbf{x}(t),\mathbf{x}(t-\tau),\ldots, \mathbf{x}(t)^2,\ldots )\in \mathbb{R}^n\) is a dictionary of nonlinear or delayed observables constructed from some time series \(\mathbf{x}(t)\in \mathbb{R}^d\). The matrix \(K\in \mathbb{R}^{n\times n}\) approximates the Koopman operator on the space defined by the observable set \(\boldsymbol{\psi}\). It is estimated from \(P\) observations via linear regression, \[\hat{K}=\arg\min_{K} \sum_t \|\boldsymbol{\psi}(t+\Delta t)- K\boldsymbol{\psi}(t)\|^2 = \arg\min_K \| \boldsymbol{\Psi}_{t+\Delta t} - K \boldsymbol{\Psi}_t \|^2, \label{eq95nvar95optprob}\tag{2}\] where \[\boldsymbol{\Psi}_t=[\boldsymbol{\psi}(t_1), \boldsymbol{\psi}(t_2),\ldots, \boldsymbol{\psi}(t_P)]\in \mathbb{R}^{n\times P} \label{eq95design95matrix}\tag{3}\] is a design matrix with columns given by the time-sampled observables \(\boldsymbol{\psi}(t_i)\). If \(P>n\) and \(\boldsymbol{\Psi}_t\) has full row rank, the optimal solution is given in terms of the pseudoinverse by \[\hat{K} = \boldsymbol{\Psi}_{t+\Delta t} \boldsymbol{\Psi}_t^+ = \boldsymbol{\Psi}_{t+\Delta t}\boldsymbol{\Psi}_t^T (\boldsymbol{\Psi}_t \boldsymbol{\Psi}_t^T)^{-1}. \label{eq95nvar95optsol}\tag{4}\] In order to improve conditioning when features are highly correlated, the pseudoinverse can be replaced by a Tikhonov-regularized form \(\boldsymbol{\Psi}_t^T (\boldsymbol{\Psi}_t \boldsymbol{\Psi}_t^T + \lambda \mathbb{I})^{-1}\) with regularization parameter \(\lambda\). In particular, when \(\boldsymbol{\psi}\) consists of delay observables with small lag \(\tau\simeq {\Delta t}\), scaling \(\lambda \propto P\) mitigates the growth of the condition number [62]. 4 can be compactly expressed in terms of empirical correlation matrices \(\hat{C}(\Delta t) = \frac{1}{P}\sum_{p=1}^P \boldsymbol{\psi}(t_p+\Delta t)\boldsymbol{\psi}(t_p)^T\): \[\hat{K} = \hat{C}(\Delta t) \hat{C}(0)^{-1} \simeq \mathbb{I}+ \Delta t \hat{C}'(0) \hat{C}(0)^{-1}, \label{eq95linreg95correl}\tag{5}\] where the last approximation is valid for small \(\Delta t\).

In the NG-RC method, instead of 1 , one uses only the state-part of \(\boldsymbol{\psi}\) as the target and thus considers the regression problem \[\mathbf{x}(t+\Delta t)=Q \boldsymbol{\psi}(t+\Delta t) \simeq H \boldsymbol{\psi}(t),\qquad Q=\begin{bmatrix}I_d & \mathbf{0}_{d\times(n-d)}\end{bmatrix}\in\mathbb{R}^{d\times n} \label{eq95ngrc95reduct}\tag{6}\] where \(H\in \mathbb{R}^{d\times n}\) denotes the NG-RC next-step operator and the projector \(Q\) selects \(\mathbf{x}\) from the observable set \(\boldsymbol{\psi}\). Thus, the minimization problem of 2 is replaced by \(\hat{H} = \arg\min_H \|X_{t+{\Delta t}} - H \boldsymbol{\Psi}_t\|^2\). Using \(X_{t+{\Delta t}} = Q \boldsymbol{\Psi}_{t+{\Delta t}}\) gives \[\hat{H} = X_{t+{\Delta t}} \boldsymbol{\Psi}_t^+ = Q \hat{K},\] showing that 1 contains the solution of 6 . Predicting beyond the next time step involves recursively lifting the state to the latent space and applying the learned next-step operator: \[\hat{\mathbf{x}}(t+{\Delta t}) = \hat{H} \mathbf{G}(\hat{\mathbf{x}}(t),\ldots),\qquad \hat{\mathbf{x}}(0)=\mathbf{x}_0. \label{eq95ngrc95autonom}\tag{7}\] The result differs, in general, from \(Q\hat{K}^j \mathbf{G}(\mathbf{x}(0),\ldots)\), because the selected observables \(\boldsymbol{\psi}\) typically do not form a Koopman-invariant subspace on which \(\mathbf{G}(\mathbf{x}(t+{\Delta t}),\ldots)=\hat{K}\mathbf{G}(\mathbf{x}(t),\ldots)\) would hold.

3 Learning Markovian dynamical systems↩︎

We now discuss learning of flow maps of fully observed Markovian dynamical systems using instantaneous and delay-augmented feature libraries. Consider a Markovian (and possibly nonlinear) dynamical system described by \[\dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t)) \label{eq95dynsys95ODE}\tag{8}\] with a state vector \(\mathbf{x}\) and a vector field \(\mathbf{f}:\Omega\subset \mathbb{R}^d \to \mathbb{R}^d\). The discrete time evolution is described by the associated flow map (propagator) \(\Phi_{\Delta t}(\mathbf{x}(t))\), which, for small \({\Delta t}\), can be expressed as a Taylor/Lie series: \[\begin{align} \mathbf{x}(t+\Delta t) = \Phi_{\Delta t}(\mathbf{x}(t)) &= \mathbf{x}(t) + \mathbf{f}(\mathbf{x}(t)) \Delta t + \sum_{l=1}^d \frac{\partial \mathbf{f}}{\partial x_l} f_l(\mathbf{x}(t)) \frac{(\Delta t)^2}{2!} + \sum_{q=3}^{\infty} \frac{(\Delta t)^q}{q!} \frac{d^q}{d t^q} \mathbf{x}(t) \\ &= \mathbf{x}(t) + \sum_{q=1}^\infty \frac{(\Delta t)^q}{q!}\mathbf{F}_q(\mathbf{x}(t)) , \end{align}\label{eq95flowmap}\tag{9}\] where \(\mathbf{F}_q\) denotes the \(q\)-th Lie-series coefficient, defined recursively by \(\mathbf{F}_1 = \mathbf{f}\) and \(\mathbf{F}_{q+1} = L_\mathbf{f}\mathbf{F}_q = (\nabla \mathbf{F}_q)\mathbf{f}\), where \(L_\mathbf{f}\mathbf{h}:= \mathbf{f}\cdot \nabla \mathbf{h}\) denotes the Lie derivative along the vector field \(\mathbf{f}\). Equivalently, along trajectories one has \(\mathbf{F}_q(\mathbf{x}(t)) = \frac{d^q}{dt^q}\mathbf{x}(t)\).

We focus henceforth on the case where the vector field \(\mathbf{f}\) has polynomial dependence on \(x_j\) with a maximum degree \(p\). The degree \(M'\) of monomials occurring in the expanded flow map \(\Phi_{\Delta t}\), truncated at the \(r\)-th time derivative, can be estimated by noting that \(\frac{d^r}{dt^r} x(t)\) has a maximum total polynomial degree of \[M'= r(p-1) + 1, \label{eq95flowmap95polydeg}\tag{10}\] leaving an error of \(\mathcal{O}((\Delta t)^{r+1})\) in the Taylor series ¹.

3.1 Learning with instantaneous feature libraries↩︎

We consider now the specific NG-RC regression problem, where one learns an empirical flow map by fitting a given feature library \(\mathbf{g}(\mathbf{x})\in\mathbb{R}^n\) \[\mathbf{x}(t+\Delta t) \simeq H \mathbf{g}(\mathbf{x}(t)), \label{eq95NGRC}\tag{11}\] where \(H\in \mathbb{R}^{d\times n}\) is a matrix of trainable weights. The associated (mean-square) training error is given by \[E_{\text{train}} \equiv \varepsilon_{\text{train}}^2 = \frac{1}{P}\sum_{j=1}^P \|\mathbf{x}(t_j+{\Delta t}) - \hat{H} \mathbf{g}(\mathbf{x}(t_j)) \|^2 , \label{eq95ngrc95trainerr}\tag{12}\] with \(\varepsilon\) denoting the root-mean-square error (RMSE). Analogously, a test error \(E_{\text{test}} = \mathbb{E}\|\mathbf{x}(t+{\Delta t}) - \hat{H} \mathbf{g}(\mathbf{x}(t)) \|^2\) can be defined by replacing the sample average by an average over the statistical distribution of \(\mathbf{x}\). In the large sample number limit (\(P\to\infty\)), and assuming ergodicity/mixing, the training error \(E_{\text{train}}\) converges to the test error \(E_{\text{test}}\).

We assume in the following that the state \(\mathbf{x}\) can be exactly represented by \(\mathbf{g}\). Using 9 and expanding the learned flow map as \(\hat{H} \mathbf{g}(\mathbf{x}) = \mathbf{x}+ {\Delta t}\hat{\mathbf{f}}(\mathbf{x}) + \mathcal{O}({\Delta t}^2)\), the training error becomes \(E_{\text{train}} \simeq {\Delta t}^2 \mathbb{E}\| \mathbf{f}(\mathbf{x}) - \hat{\mathbf{f}}(\mathbf{x}) \|^2 + \mathcal{O}({\Delta t}^3)\) (large \(P\) limit). This motivates defining the truncation error [91] \[E_{\text{trunc}}\equiv \varepsilon^2_{\text{trunc}}\equiv \frac{E_{\text{train}}}{{\Delta t}^2} \label{eq95ngrc95truncerr}\tag{13}\] as a measure indicating how well the flow map is learned. Accordingly, \(E_{\text{trunc}}=\mathrm{const}\) as \({\Delta t}\to 0\) indicates that already the vector field \(\mathbf{f}\) is not exactly representable by the feature library \(\mathbf{g}\). This heuristic is made precise in 6.

Focusing now on the case of polynomial vector fields, we can distinguish two cases:

1. Polynomial feature library. If \(\mathbf{g}\) consists of instantaneous monomials up to total degree \(M\) in \(x_j(t)\), the largest representable Lie order \(r^\ast\) in 9 follows from \(M\ge M'\) as \[r^\ast=\left\lfloor \frac{M-1}{p-1} \right\rfloor, \label{eq95dtexp95poly}\tag{14}\] where \(\lfloor \cdot \rfloor\) is the largest integer less than or equal to the argument. The trivial linear case \(p=1\) can be treated separately. Accordingly, the training and truncation errors behave as [see 45 ] \[E_{\text{train}} \propto {\Delta t}^{2(r^\ast+1)},\qquad E_{\text{trunc}} \propto {\Delta t}^{2 r^\ast}. \label{eq95trainerr95poly}\tag{15}\] For \(M<p\), the model is misspecified and \(E_{\text{trunc}}>0\) as \({\Delta t}\to 0\).

2. Non-polynomial feature library.

   1. Generic case. If the vector field is not exactly representable by the feature span, the trivial asymptotic scaling [41 ] \[E_{\text{train}} \propto {\Delta t}^2,\qquad E_{\text{trunc}}=\mathrm{const}.,\qquad {\Delta t}\ll {\Delta t}_\times \label{eq95Etrunc95norep}\tag{16}\] applies. An estimate for the crossover time scale \({\Delta t}_\times\) is given in 43 . For larger \({\Delta t}\), a pre-asymptotic scaling law of the form \[E_{\text{trunc}} \propto {\Delta t}^{2r_*},\qquad {\Delta t}\gg {\Delta t}_\times \label{eq95Etrunc95preasymp}\tag{17}\] can emerge provided the approximation error of Lie orders earlier than \(r_*+1\) is sufficiently small [see 42 ].

   2. Fourier feature library. A particular relevant case is a Fourier feature library, \[\mathbf{g}(\mathbf{x})=\left\{\mathbf{x}, \sin\!\bigl(\omega_0 \, \mathbf{k}\cdot \mathbf{x}\bigr),\cos\!\bigl(\omega_0 \, \mathbf{k}\cdot \mathbf{x}\bigr):\mathbf{k}\in \{0,\dots,N\}^d\right\}, \label{eq95Fourier95library}\tag{18}\] with \(N\) modes (per dimension) and base frequency \(\omega_0\) chosen such that \(\omega_0 R\ll 1\) on the data domain \(\Omega = [-R,R]^d\) (see 6.4 for details; the state \(\mathbf{x}\) is included here to transfer those results to the full flow-map prediction). Spectral models are widely used in forecasting [92]. While also here the asymptotic scaling 16 holds for small \({\Delta t}\), a pre-asymptotic regime described by 17 with an exponent \[r_* \gtrsim \left\lfloor \frac{N-1}{p-1} \right\rfloor \label{eq95Fourier95scal95exp}\tag{19}\] is expected [see 58 ].

The above statements apply also to non-polynomial dynamical systems and generic feature libraries (see 6), except for the specific predictions involving the polynomial degree \(p\) [14 19 ]. Effects of ill-conditioning [62], [63] and noise give rise to a lower bound \(E_{\text{floor}}(M,\lambda)\) to \(E_{\text{train}}\). When the feature library can reproduce the flow map to some nontrivial order, the training error also controls the forecast horizon (see 8).

As will be discussed in detail in 3.2, delay terms do not by themselves remove a structural mismatch to the nonlinear Markovian vector field, but they can change the one-step training-error scaling through multistep extrapolation. Thus the truncation-error diagnostic is cleanest for instantaneous libraries, or for delay models only after separating genuine flow-map representation from delay-induced interpolation.

The scaling laws derived here describe the approximation-bias component of the one-step prediction risk (see 6.5). A feature library is called well specified only if the flow map belongs to the span of the library. However, this is usually not fulfilled for finite libraries, since even for polynomial vector fields the exact flow map is typically of non-polynomial form. We therefore use the weaker notion of being well specified to Lie order \(r\), meaning that the first \(r\) Lie-series coefficients of the flow map are contained in the feature span. In this sense, a polynomial library with \(M\geq p\) is bias-free at the vector-field level, whereas polynomial libraries with \(M<p\), generic finite Fourier libraries, and linear delay-only models for nonlinear Markovian systems are misspecified.

3.2 Learning Markovian systems with delay-feature libraries↩︎

We now extend the previous setting by including time-lagged states \(\mathbf{x}(t-k\tau)\), \(k=0,\dots,\Gamma-1\) in the feature library (where \(\Gamma\) denotes the total number of taps), and analyze how the resulting multistep predictor represents the Markovian flow map 9 . We assume \(\tau=\Delta t\) for simplicity. The learning problem is then given by \[\mathbf{x}(t+\Delta t) = K \mathbf{G}\big(\{\mathbf{x}(t-k\Delta t)\}_{k=0}^{\Gamma-1}\big), \label{eq95delay95pred}\tag{20}\] where the delay feature library \(\mathbf{G}\in\mathbb{R}^{D_D}\) is constructed from all monomials in the delay coordinates \(\{\mathbf{x}(t-k\Delta t)\}\) up to some maximal total degree \(M\).

According to 9 , one has \[\mathbf{x}(t-k\Delta t) = \Phi_{-k\Delta t}(\mathbf{x}(t)) = \sum_{l=0}^{\infty} \frac{(-k\Delta t)^l}{l!} \frac{\mathrm{d}^l}{\mathrm{d}t^l} \mathbf{x}(t) = \mathbf{R}_k(\mathbf{x}(t)) + \mathcal{O}\big((k\Delta t)^{r+1}\big), \label{eq95delay95taylor}\tag{21}\] where, by truncating the expansion at order \(l=r\), we can express each delay state as a vector \(\mathbf{R}_k\) of polynomials (of degree at most \(M' = r(p-1)+1\)) in \(\mathbf{x}(t)\). Let the feature vector \(\boldsymbol{\psi}(\mathbf{x}) \in \mathbb{R}^{D_I}\) consist of all instantaneous monomials in \(\mathbf{x}\) up to some maximal degree \(M_I\ge M'\), such that all components of \(\mathbf{R}_k(\mathbf{x})\) can in principle be written as linear combinations of the entries of \(\boldsymbol{\psi}(\mathbf{x})\).

By 21 , every component of the library \(\mathbf{G}\) is itself a polynomial in \(\mathbf{x}(t)\) of degree at most \(M M'\). Hence, provided \(M_I \ge M M'\), we can express \(\mathbf{G}\) as a linear combination of the entries of \(\boldsymbol{\psi}\): \[\mathbf{G}\big(\{\mathbf{x}(t-k\Delta t)\}_{k=0}^{\Gamma-1}\big) \approx A \boldsymbol{\psi}(\mathbf{x}(t)), \label{eq95delay95change95basis}\tag{22}\] where \(A\in\mathbb{R}^{D_D\times D_I}\) follows from the expansion in 21 . Assume that the truncated delay library spans the chosen instantaneous polynomial basis used by the truncated flow map, so that \(A\) has full column rank and thus \(A^+ A = \mathbb{I}_{D_I}\). As in 11 , we express the Markovian flow map in the instantaneous basis, \(\mathbf{x}(t+\Delta t) = \Phi_{\Delta t}(\mathbf{x}(t)) \approx H \boldsymbol{\psi}(\mathbf{x}(t))\) where the entries of \(H\in\mathbb{R}^{d\times D_I}\) follow from the Taylor expansion 9 . We then invert 22 to express \(\boldsymbol{\psi}\) in terms of the delay library: \[\mathbf{x}(t+\Delta t) \approx H A^+ \mathbf{G}\big(\{\mathbf{x}(t-k\Delta t)\}_{k=0}^{\Gamma-1}\big) + \mathcal{O}(\|H A^+\|{\Delta t}^{r+1}). \label{eq95delay95flow95representation}\tag{23}\] This shows that a library of delay coordinates can be linearly transformed into a next-step predictor of the form \(\hat{K} = H A^+\). Note, however, that for small \(\Delta t\) the columns of \(A\) [22 ] become nearly linearly dependent because the different delays share the same time derivatives in 21 , leading to an ill-conditioned pseudoinverse \(A^+\) and to the strong correlation and regularization effects observed when using dense delay libraries in NG-RC [62]. This can also cause the error in 23 to diverge for \({\Delta t}\to 0\).

Depending on the choice of the library features, the regressor can realize several predictor types via fitting \(\hat{K}\). If \(\mathbf{G}\) contains only linear delay coordinates, the learned update is linear in the delay state. It may nevertheless achieve small one-step error by exploiting temporal interpolation/extrapolation along the observed trajectory. However, it does not represent a state-only nonlinear Markovian flow map and generally does not yield a consistent autonomous rollout for nonlinear systems (see also 7). This is consistent with the need for lifted observables in Koopman/EDMD approaches [22], [27], [29]. Nevertheless, one can use such a linear multistep predictor to understand the general impact of delay features on the training error. A canonical linear predictor is the extrapolation stencil \[\tilde{\mathbf{x}}(t+\Delta t)= \sum_{k=0}^{\Gamma-1} \zeta_k \mathbf{x}(t-k\Delta t) + \mathcal{O}(\Delta t^\Gamma),\qquad \zeta_k=(-1)^k \binom{\Gamma}{k+1}. \label{eq95delay95extrap95predictor}\tag{24}\] The coefficients \(\zeta_k\) follow straightforwardly by inserting the Taylor expansion of \(\mathbf{x}(t-j\Delta t)\) into the left and right hand sides and equating the coefficients at equal powers of \(\Delta t\). For generic data and at fixed \(M\), this canonical predictor achieves a training error of ² \[\begin{align} E_{\mathrm{train}}(\Gamma) &\le \mathbb{E}\big\|\mathbf{x}(t+\Delta t)-\tilde{\mathbf{x}}(t+\Delta t)\big\|^2 \simeq {\Delta t}^{2 \Gamma} \mathbb{E}\|\mathbf{x}^{(\Gamma)}(t)\|^2 \sim \mathcal{O}(\Delta t^{2 \Gamma}). \label{eq95delay95trainerr95asy} \end{align}\tag{25}\] Since the predictor 24 is contained in the feature class for \(M\ge 1\), the optimal training error admits the combined bound following from 25 15 : \[E_\mathrm{train}(\Gamma) \sim \mathcal{O}\bigl(\Delta t^{2\max\{r^\ast(M)+1,\Gamma\}}\bigr). \label{eq95delay95trainerr95asy95maxform}\tag{26}\] Accordingly, delay coordinates generally decrease the (unregularized) training error monotonically with the number of taps \(\Gamma\). Note, however, that these scaling laws can be affected by \(\Gamma\)-dependent prefactors. By allowing the \(\zeta_k\) to pick up frequency information, the predictor 24 can represent a harmonic signal, as shown in 7. In numerical experiments, predictor forms specific to the dataset are often observed [62], [65].

4 Numerical experiments↩︎

We now discuss numerical experiments on flow map reconstruction of nonlinear Markovian dynamical systems, using polynomial as well as Fourier feature libraries. To improve numerical conditioning and stability, we typically train on increments \(\mathbf{x}(t+{\Delta t})-\mathbf{x}(t)\) instead of \(\mathbf{x}(t+{\Delta t})\) [5], [65].

Metrics.

To assess the quality of the flow map reconstruction, we consider, besides the training error [12 ], the forecast horizon \(T_{\text{fch}}\), defined as the smallest time where the autonomously predicted trajectory \(\hat{\mathbf{y}}(t)\equiv \hat{\mathbf{x}}(t+{\Delta t})\) [7 ] deviates more than a certain threshold from the true trajectory \(\mathbf{y}(t)\equiv \mathbf{x}(t+{\Delta t})\): \[T_{\text{fch}}= \inf \left\{ t : \| \mathbf{y}(t) - \hat{\mathbf{y}}(t) \| > \kappa \sigma \right\}, \label{eq95fc95hor}\tag{27}\] Here, \(\sigma^2= \frac{1}{P}\sum_{i=1}^P \| \bar \mathbf{y}- \mathbf{y}(t_i)\|^2\) is the estimated variance, \(\bar \mathbf{y}= \frac{1}{P}\sum_{i=1}^P \mathbf{y}(t_i)\), and we typically take \(\kappa=1\). In the noiseless, overdetermined regime (\(P\gg n\)) considered here, the train-test gap is typically small, so we focus on training error. We use low solver tolerances \(\sim 10^{-13}\) [93], for which the training error can reach \(\varepsilon\sim 10^{-13}\).

Model systems.

As primary benchmarks, we use the Halvorsen and Lorenz-63 models [94], which have polynomial vector fields of order 2. The Halvorsen model is given by: \[\begin{align} \dot{x} &= -a x - b y - b z - y^2, \\ \dot{y} &= -a y - b z - b x - z^2, \\ \dot{z} &= -a z - b x - b y - x^2, \end{align} \label{eq95Halvorsen}\tag{28}\] with standard parameters \(a=1.4\) and \(b=4\). Reported values for the largest Lyapunov exponent range between \(0.72\) [95] and \(0.81\) [96], depending slightly on the method used. The Lorenz-63 model is given by: \[\begin{align} \dot{x} &= \sigma (y - x) ,\\ \dot{y} &= x (\rho - z) - y ,\\ \dot{z} &= xy - \beta z , \end{align} \label{eq95Lorenz63}\tag{29}\] with standard parameters \(\sigma = 10\), \(\rho = 28\), \(\beta = 8/3\). The largest Lyapunov exponent is reported between \(0.86\) [95] and \(0.91\) [97]. We also consider two other chaotic systems with cubic nonlinearities: the Sprott cubic jerk system [98], [99], defined by \[\begin{align} \dot{x} &= y, \\ \dot{y} &= z, \\ \dot{z} &= -a z + x y^2 - x^3 \end{align} \label{eq95sprott95sys}\tag{30}\] with \(a=3.6\) and largest Lyapunov exponent \(\lambda_1\approx 0.14\), and the Rabinovich-Fabrikant system [100], defined by \[\begin{align} \dot{x} &= y(z - 1 + x^2) + \gamma x, \\ \dot{y} &= x(3z + 1 - x^2) + \gamma y, \\ \dot{z} &= -2z(\alpha + xy), \end{align} \label{eq95RabFab95sys}\tag{31}\] with \(\alpha=1.1\), \(\gamma=0.87\) and largest Lyapunov exponent \(\lambda_1\approx 0.12\).

4.1 Polynomial feature library↩︎

Here, the NG-RC feature library consists of all monomials of the state variables and their delay terms up to depth \(\Gamma\) (\(\Gamma=1\) corresponding to no delay terms) and nonlinear degree \(M\). We typically whiten the features by subtracting the empirical mean and normalize by the empirical standard deviation. Specifically, given the design matrix \(\Psi \in \mathbb{R}^{n \times P}\) (3 , excluding the bias row if present), we use in the NG-RC fit the standardized quantities \[\tilde{\Psi}_{jp} = \frac{\Psi_{jp} - \mu_j}{\sigma_j},\qquad \text{where}\quad \mu_j = \frac{1}{P} \sum_{p=1}^{P} \Psi_{jp}, \quad \sigma_j^2 = \frac{1}{P} \sum_{p=1}^{P} (\Psi_{jp} - \mu_j)^2. \label{eq95whitened95design95matrix}\tag{32}\]

4.1.1 Effect of monomial degree↩︎

1 illustrates the training error and forecast horizon (obtained for the Halvorsen model) in dependence of the monomial degree \(M\) of the feature library and the Ridge regularization \(\lambda\). For raw monomial libraries, the RMS training error behaves approximately as \(\varepsilon\sim a ({\Delta t})^{b M}+ \varepsilon_{\text{floor}}(M,\lambda)\), with constants \(a,b\). The first term is the Taylor-remainder term, which monotonically shrinks as \(M\) increases [see 15 ]. The second term can dominate for ill-conditioned features, which arises when features have strong disparities in their magnitude or are nearly collinear. By contrast, for whitened features, the term \(\varepsilon_{\text{floor}}(M,\lambda)\) vanishes as demonstrated by the bottom row plotted in 1 (a). Hence, in this case it is advantageous to train without Ridge regularization (\(\lambda=0\)). Including delays can improve the training error to a certain extent at fixed monomial order, as shown in 1 (b). We will return to a more detailed analysis of delays below.

4.1.2 Effect of temporal resolution↩︎

2 (a,b) shows in more detail the behavior of the training error, following closely 15 until a floor is reached. The forecast horizon \(T_{\text{fch}}\) grows according to 75 with decreasing \(\varepsilon\). The saturation of the training error and forecast horizon is due to the integrator accuracy used to generate the trajectories and can be further increased by using a higher floating point precision [11] ³.

2 (c) demonstrates that the scaling behavior with \({\Delta t}\) of the truncation error \(\varepsilon_{\text{trunc}}=\varepsilon/{\Delta t}\) [13 ] can distinguish a misspecified model (\(\varepsilon_{\text{trunc}}\sim\mathrm{const}\)) from a model that is able to learn the vector field (\(\varepsilon_{\text{trunc}}\sim{\Delta t}^{r^\ast}\) with a \(r^\ast>0\)). The former case corresponds to \(M<p\), where \(p\) is the nonlinear degree of the vector field. The latter case corresponds to \(M\geq p\) and leads to a consistent autoregressive evolution with a nonzero \(T_{\text{fch}}\).

Notably, this characteristic scaling behavior also occurs for an effective \({\Delta t}= k{\Delta t}_0\) (with \(k\in\mathbb{N}\)) obtained via subsampling the time series of a higher temporal resolution \({\Delta t}_0\) (see 2 (d)). This is relevant for practical situations, where one may not be able to change the temporal resolution \({\Delta t}\) of the data.

In 3, we repeat the same analysis for the Lorenz-63 model, the Sprott cubic jerk, and the Rabinovich-Fabrikant system [30 31 ]. The latter two systems have cubic nonlinearities, such that the truncation error vanishes with \({\Delta t}\) only for feature maps \(M\geq 3\), as confirmed by the numerics (panels c–f).

4.1.3 Semigroup consistency↩︎

The semigroup consistency of a learned flow map \(\widehat\Phi_{\Delta t}\) requires \[\widehat\Phi_{\ell{\Delta t}}\approx \widehat\Phi_{{\Delta t}}^\ell \label{eq95semigroup95flow}\tag{33}\] where \(\widehat\Phi_{{\Delta t}}^\ell = \widehat\Phi_{{\Delta t}} \circ \widehat\Phi_{{\Delta t}} \circ \cdots \circ \widehat\Phi_{{\Delta t}}\) (\(\ell\) times). When the model can reproduce the early Lie derivatives of the true flow map, the error in the approximation in 33 is proportional to the training error [see 66 ], \[\mathcal{S}_m:= \|\widehat\Phi_{\ell{\Delta t}}-\widehat \Phi_{{\Delta t}}^\ell\|\sim \varepsilon_{\text{train}}. \label{eq95semigroup95consistency}\tag{34}\] which is illustrated in 4. Note that, in general, the semigroup consistency only ensures that the flow map generates autonomous dynamics, not that the model has learned the ground truth.

4.1.4 Effect of delay terms↩︎

5 demonstrates the influence of the number of time taps \(\Gamma\) on training error and forecast horizon at fixed monomial degree \(M\) (Halvorsen model with \({\Delta t}=0.01\)). When using only linear delays (\(M=1\)), the NG-RC model becomes a linear autoregressive process and learns a harmonic approximation, often with a slowly decaying amplitude envelope (see 6) [68]–[70]. While arbitrarily low training errors can be achieved up to noise and conditioning effects [see 26 ], the model is unable to produce accurate forecasts for nonlinear systems.

As shown in 7 in more detail, the empirical training error follows the theoretical scaling prediction 26 only for small delay numbers \(\Gamma\) and sufficiently large monomial orders \(M\). For larger \(\Gamma\), the training error eventually reaches the noise floor controlled by accuracy of the generated trajectories. (This happens earlier for large \(M\) since training error \(\varepsilon_0\) without delays is already small.) For small monomial orders \(M\), the observed scaling of \(\varepsilon\) with \(\Gamma\) differs significantly from the asymptotic prediction 26 . Besides a possible \(\Gamma\)-dependence of prefactors, we attribute this discrepancy to the finite temporal resolution.

7 shows moreover that adding delay terms does not automatically improve the forecast horizon: this in general requires the model to learn some finite-\(\Delta t\) approximation of the flow map, which is only possible when the monomial degree \(M\) is at least as large as the maximum degree \(p\) of the vector field of the system (\(p=2\) for the Halvorsen and Lorenz-63 systems, and \(p=3\) for the Sprott and Rabinovich-Fabrikant systems). However, delay terms, especially when combined with low-order nonlinearities, can support a harmonic approximation on top of a rough flow map representation. This can result in nonzero forecast horizons even for \(M<p\), as observed for the Rabinovich-Fabrikant system in 7 (d).

4.1.5 Readout weights↩︎

In 8, we compare the learned readout weights of a NG-RC trained on the Halvorsen system to the theoretical predictions in 3.2. A feature map consisting only of linear delay terms results in the generic weights of 24 (panels a,b). Conversely, if the feature map consists only of instantaneous monomials (panels c,d), the true flow map 9 is recovered to that order (see also [61], [62], [65]). A purely Markovian flow map can also be approximated by a linear combination of instantaneous and delay features. In this case, the learned weights (panel e) are approximated by the expression \(H A^+\) [23 , panel f].

4.2 Fourier feature library↩︎

We now consider a NG-RC model with Fourier features [18 ] as an example of reconstructing a truncated flow map with a non-polynomial feature library. As illustrated in 9, while the overall prediction accuracy is comparable to the polynomial case (depending slightly on the considered dynamical system), the most distinct difference is the presence of a plateau of the truncation error \(\varepsilon/{\Delta t}\) for small \({\Delta t}\). This is due to the functional mismatch between the Fourier features and the low-order Lie coefficients generated by a polynomial vector field [see 16 ]. Since these coefficients are of polynomial form, a finite Fourier feature library can only approximate them, generically leading to nonzero first defect [see 37 for details].

Actual power law behavior of the truncation error can appear only in a pre-asymptotic regime. This is approximately described by 17 with an empirical exponent \(r_*\propto N\) for small \(N\), which is consistent with the conservative lower bound in 19 . For larger \(N\gtrsim 4\), we find the exponent to saturate. Due to the highly nonlinear nature of the Fourier feature library, the truncation error does not show a clear threshold behavior as in the polynomial case that would indicate exact flow map representation.

10 illustrates the semigroup defect for the learned Fourier-feature flow map. Although the Fourier features cannot exactly reproduce the Lie derivatives, the deviations are still small enough for the semigroup defect to closely follow the training error. The learned representation of the flow map can be understood by Taylor expanding the Fourier feature map, thereby bringing it into a polynomial form. 11 reveals close agreement with the theoretical expectations, demonstrating that the model linearly combines the Fourier features in order to represent the monomial terms occurring in the true flow map.

5 Summary↩︎

We analyzed here NVAR/NG-RC models trained on time series generated by Markovian nonlinear dynamical systems and showed that a structural mismatch between the vector field and the feature map results in characteristic scaling behaviors of the (RMS) training error \(\varepsilon_{\text{train}}\sim{\Delta t}^{1+r^\ast}\) with temporal resolution \({\Delta t}\). The exact asymptotic exponent \(r^\ast\) is determined by the first order in the Taylor/Lie expansion of the flow map that can not be represented exactly within the feature space. Libraries with exact reproduction result in a non-zero asymptotic exponent \(r^\ast\) because they annihilate the first few orders of the Lie derivatives exactly. By contrast, libraries that can merely approximate the Lie coefficients generally lead to a “trivial” asymptotic law \(\varepsilon_{\text{train}}\sim {\Delta t}\) (\(r^\ast=0\)) for \({\Delta t}\to 0\). However, they may exhibit steeper pre-asymptotic scaling when the low-order Lie coefficients are approximated with a sufficiently small error, resulting in an effective exponent \(r_*>0\) at larger time resolutions \({\Delta t}\).

In order to illustrate this generic mechanism, we focused on dynamical systems with polynomial vector fields. Since the Lie coefficients are themselves polynomials in this case, a NVAR/NG-RC equipped with a delay-free monomial library can exactly reproduce the flow map to a given order in \({\Delta t}\). Accordingly, this is reflected by a superlinear scaling of the RMS training error, \(\varepsilon_{\text{train}}\sim{\Delta t}^{1+r^\ast}\) with an \(r^\ast>0\) that directly depends on the nonlinear order of the vector field [see 15 and 2 3]. By contrast, a finite Fourier feature library can only approximately reproduce the polynomial Lie coefficients, leading to an asymptotic plateau of the truncation error (see 9). Note that an investigation of the forecast horizon alone would not immediately reveal these qualitative differences of the feature maps: although the forecast horizon is primarily controlled by the training error (see 8), it can be additionally affected by model instabilities and the intrinsic Lyapunov divergence.

Under Markovian dynamics, delay features are deterministic functions of the current state and can thus in principle be represented as linear combinations of instantaneous features [see 23 ]. They are thus not crucial for learning Markovian systems, but still generally reduce the one-step training error, since they allow the regressor to form better interpolants [see, e.g., [62] and 25 ]. However, they do not generically lead to better generalization behavior, i.e., longer forecast horizons, unless also sufficiently high-order nonlinearities are present in the feature map (see 5 7). This constitutes an example of a model failing to generalize despite having a small training error. Delay terms are essential for learning non-Markovian systems, which will be interesting to analyze in future work.

In summary, one can identify three distinct learning regimes of an NVAR/NG-RC model when trained on time series of nonlinear dynamical systems: (i) a linear delay model learns a harmonic approximation (see 7), which can render a small training error but inconsistent forecasts [68]–[70]; (ii) a nonlinear model with a structural mismatch to the underlying vector field can achieve long forecast horizons, but the RMS training error will scale suboptimally (linearly) as \({\Delta t}\to 0\); (iii) a nonlinear model that structurally matches the vector field achieves long forecast horizons and an optimal (superlinear) scaling of the RMS training error with \({\Delta t}\). In situations where the underlying data-generating process is unknown, the truncation error scaling can be used via subsampling to assess the structural properties of the underlying system and to guide the design of the feature library. It will be useful to analyze this further in the future. The semigroup property of the flow map is a standard diagnostic for flow map learning [29], [43], [46]. For a well-specified model, the semigroup defect closely follows the training error (see 4 10). We finally remark that the present results can also help to construct optimal encodings in quantum machine learning models [101], including quantum RCs [102], for which the feature library can often be explicitly stated [103]. In particular, common quantum encodings lead to a finite Fourier feature library [104].

This project was made possible by the DLR Quantum Computing Initiative and the Federal Ministry for Research, Technology and Space; https://qci.dlr.de/NeMoQC.

6 Feature-library structure and asymptotic error scaling↩︎

As in 3, we consider a dynamical system \[\dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t)), \qquad \mathbf{x}(t) \in \mathbb{R}^d,\] with flow map \(\Phi_{\Delta t}\). We study predicting the increment \[\mathbf{F}_{\Delta t}(\mathbf{x}) := \Phi_{\Delta t}(\mathbf{x})-\mathbf{x}\] from a prescribed feature library \(\mathbf{g}(\mathbf{x}) = (g_1(\mathbf{x}),\dots,g_n(\mathbf{x})) \in \mathbb{R}^n\) by a linear readout. A central question is how the structure of the feature library influences the small-\({\Delta t}\) behavior of the mean-square error and gives rise to possible pre-asymptotic scaling laws. The results apply to the full state prediction problem in the main text since we assume that \(\mathbf{x}\) is included in \(\mathbf{g}\) if necessary.

We assume that the data are distributed according to a probability measure \(\mu\) on a domain \(\Omega \subset \mathbb{R}^d\), and we identify the idealized mean-square training error with the projection error in \(L^2(\mu;\mathbb{R}^d)\). This setting captures the regime in which the empirical training error is dominated by approximation rather than by finite-sample, optimization, or numerical effects.

6.1 A library-agnostic projection framework↩︎

Given the feature map \(\mathbf{g}\), the associated model class consists of all vector-valued functions that can be represented by a linear readout, \[V := \left\{ H \mathbf{g}(\cdot) : H \in \mathbb{R}^{d \times n} \right\} \subset L^2(\mu;\mathbb{R}^d).\] Here \(H\) is only a generic coefficient matrix parameterizing the elements of \(V\). For a fixed target function \(\mathbf{h}\in L^2(\mu;\mathbb{R}^d)\), the projection \(\mathcal{P}\mathbf{h}\) is the best approximation to \(\mathbf{h}\) among all functions in \(V\). We denote by \(H_{\mathbf{h}}\) a coefficient matrix that represents this projected function, \[(\mathcal{P}\mathbf{h})(\mathbf{x}) = H_{\mathbf{h}}\mathbf{g}(\mathbf{x}).\] This target-dependent matrix is determined by the normal equations for the least-squares projection: defining the population feature covariance and cross-covariance \[G := \int_\Omega \mathbf{g}(\mathbf{x})\mathbf{g}(\mathbf{x})^T \, d\mu(\mathbf{x}), \qquad C_{\mathbf{h}} := \int_\Omega \mathbf{h}(\mathbf{x})\mathbf{g}(\mathbf{x})^T \, d\mu(\mathbf{x}),\] one obtains \[H_{\mathbf{h}} G = C_{\mathbf{h}},\] so that \(H_{\mathbf{h}}=C_{\mathbf{h}}G^{-1}\) when \(G\) is invertible. If \(G\) is singular, the representation of the projected function \(\mathcal{P}\mathbf{h}\) by coefficients need not be unique. In this case one may either remove redundant features or choose the Moore–Penrose solution \(H_{\mathbf{h}}=C_{\mathbf{h}}G^+\). Equivalently, the projector is obtained by Gram–Schmidt orthonormalization of \(g_1,\dots,g_n\) in \(L^2(\mu)\).

Let \(\mathcal{Q}:= I-\mathcal{P}\), so that \(\mathcal{Q}\mathbf{h}\) is the component of \(\mathbf{h}\) orthogonal to \(V\), i.e., the part of \(\mathbf{h}\) that is not exactly representable by the feature library. The population mean-square error for increment prediction is then given by \[\mathcal{E}({\Delta t})=\inf_{\mathbf{h}\in V} \|\mathbf{F}_{\Delta t}-\mathbf{h}\|_{L^2(\mu)}^2 =\|\mathcal{Q}\mathbf{F}_{\Delta t}\|_{L^2(\mu)}^2. \label{eq:population95projection95error}\tag{35}\] This expression depends only on the span of \(V\). Therefore any two feature libraries with the same span, including monomial and orthogonal polynomial bases of the same polynomial subspace, have identical projection error for every \({\Delta t}\) and thus identical asymptotic scaling.

The increment map admits the usual Taylor/Lie-series expansion, using the same Lie-series coefficients \(\mathbf{F}_m\) as in 9 : \[\mathbf{F}_{\Delta t}(\mathbf{x}) =\sum_{m=1}^\infty \frac{{\Delta t}^m}{m!} \mathbf{F}_m(\mathbf{x}), \qquad \mathbf{F}_1 = \mathbf{f}, \qquad \mathbf{F}_{m+1} = L_\mathbf{f}\mathbf{F}_m, \label{eq:Lie95series95increment}\tag{36}\] where \(L_\mathbf{f}\mathbf{h}:= \mathbf{f}\cdot \nabla \mathbf{h}\) is the Lie derivative along the vector field. Defining the \(m\)-th library defect by \[\mathbf{q}_m := \mathcal{Q}\mathbf{F}_m, \qquad \varepsilon_m := \|\mathbf{q}_m\|_{L^2(\mu)}, \label{eq:defect95profile}\tag{37}\] we obtain the asymptotic expansion \[\mathcal{E}({\Delta t}) = \left\| \sum_{m=1}^\infty \frac{{\Delta t}^m}{m!} \mathbf{q}_m \right\|_{L^2(\mu)}^2 = \sum_{m,n \ge 1} \frac{{\Delta t}^{m+n}}{m! n!}\langle \mathbf{q}_m, \mathbf{q}_n \rangle_{L^2(\mu)}. \label{eq:exact95error95expansion}\tag{38}\] It is useful to introduce \(r^\ast\) as the largest order to which the library represents the successive Lie derivatives exactly, i.e., \[r^\ast := \max \left\{ r \ge 0 : \mathbf{F}_1,\dots,\mathbf{F}_r \in V \right\},\] with the convention \(r^\ast=0\) if already \(\mathbf{F}_1=\mathbf{f}\notin V\). The first non-vanishing defect order is then \[m_0 := \min \{ m \ge 1 : \varepsilon_m > 0 \} = r^\ast+1, \label{eq:defect95order95nonv}\tag{39}\] where the last equality holds provided \(r^\ast<+\infty\). This defines the leading term of the error as \[\mathcal{E}({\Delta t})=\frac{\varepsilon_{m_0}^2}{(m_0!)^2} {\Delta t}^{2m_0}+\mathcal{O}\!\left({\Delta t}^{2m_0+1}\right), \label{eq:true95asymptotic95scaling}\tag{40}\] implying the exact asymptotic exponent \(2m_0\). The asymptotic exponent for the (mean-square) truncation error [13 ] is, correspondingly, \(2r^\ast = 2(m_0-1)\). This shows that the small-\({\Delta t}\) behavior of the error is determined by the structural alignment between the feature map and the successive Lie derivatives \(\mathbf{f}, L_\mathbf{f}\mathbf{f}, L_\mathbf{f}^2 \mathbf{f},\dots\).

6.2 Asymptotic versus pre-asymptotic regimes↩︎

If the vector field is not exactly representable by the feature library, \(\mathbf{f}\notin V\), then \(\varepsilon_1 > 0\) and the true asymptotic law is necessarily \[\mathcal{E}({\Delta t}) = \varepsilon_1^2 {\Delta t}^2 + \mathcal{O}({\Delta t}^3), \label{eq:asymptotic95scaling95norep}\tag{41}\] If \(\varepsilon_1\) is much smaller than \(\varepsilon_m\) for some \(m>1\), different pre-asymptotic scalings are possible. Then, over a window of \({\Delta t}\) in which the \(m\)-th term dominates in 38 , the error behaves effectively as \[\mathcal{E}({\Delta t})\approx \frac{\varepsilon_m^2}{(m!)^2} {\Delta t}^{2m}. \label{eq:preasymptotic95scaling95generic}\tag{42}\] Accordingly, one should distinguish between the exact asymptotic exponent \(2m_0\) [40 ], which is determined by the first nonzero defect, and the effective pre-asymptotic exponent \(2 m\). This is the main difference between libraries that reproduce the relevant Lie derivatives exactly and libraries that only approximate them well.

The crossover from an intermediate \({\Delta t}^{2m}\) regime back to the asymptotic \({\Delta t}^2\) regime is determined by a balance of defect amplitudes. If \(m>1\) is the dominant pre-asymptotic order, then the crossover scale is obtained by equating the first- and \(m\)-th defect contributions, \(\varepsilon_1 {\Delta t}\sim (\varepsilon_m/m!) {\Delta t}^m\), which gives \[{\Delta t}_{\times}^{(1,m)} \asymp \left( \frac{m! \, \varepsilon_1}{\varepsilon_m} \right)^{1/(m-1)}. \label{eq:crossover95scale}\tag{43}\]

6.3 Polynomial vector fields↩︎

In order to make the previous discussion more concrete, consider the case when \(\mathbf{f}\) is polynomial of total degree \(p\). In that case each Lie coefficient \(\mathbf{F}_m = L_\mathbf{f}^{m-1} \mathbf{f}\) is again polynomial, with the degree bound \[\mathrm{deg}\,\mathbf{F}_m \le 1 + m(p-1), \qquad p>1. \label{eq:degree95growth}\tag{44}\] Hence the question reduces to how well the library reproduces polynomials of the degrees that are produced by the successive Lie derivatives.

6.3.1 Libraries with exact polynomial reproduction↩︎

Let \(\Pi_M^d\) denote the space of \(\mathbb{R}^d\)-valued polynomials of total degree at most \(M\). We define the polynomial reproduction order of the library as the largest integer \(M_\ast\) such that \[\Pi_{M_\ast}^d \subset V.\] Any library whose span \(V\) contains \(\Pi_{M_\ast}^d\) thus yields the same asymptotic scaling for polynomial dynamics, since all Lie coefficients satisfying \(1+m(p-1) \le M_\ast\) are represented exactly. Specifically, one has \(\mathbf{F}_1,\dots,\mathbf{F}_{r^\ast}\in V\) with \[r^\ast = \left\lfloor \frac{M_\ast-1}{p-1} \right\rfloor. \label{eq:poly95scaling95generic}\tag{45}\] For a generic polynomial vector field of degree \(p\), the first unresolved Lie coefficient thus appears at order \(r^\ast+1= m_0\), which provides the asymptotic exponent in 40 .

Thus polynomial libraries change the true asymptotic exponent because they annihilate the first defects exactly. This is in sharp contrast with the generic behavior of finite non-polynomial libraries, for which the same early defects are typically small but not exactly zero.

6.3.2 Libraries without exact polynomial reproduction↩︎

For a general non-polynomial feature library the most informative quantity is not a single degree cutoff, but rather the defect profile against polynomial test functions. Consider the (worst-case) error in reconstructing a degree-\(D\) polynomial by the feature library, \[\delta(D) := \sup_{\substack{\mathbf{h}\in \Pi_D^d \\ \|\mathbf{h}\| = 1}} \|\mathcal{Q}\mathbf{h}\|_{L^2(\mu)}. \label{eq:polynomial95defect95profile}\tag{46}\] If \(\mathbf{f}\) is polynomial of degree \(p\), then the Lie coefficient \(\mathbf{F}_m \in \Pi_{1+m(p-1)}^d\), so 46 yields \[\varepsilon_m = \|\mathcal{Q}\mathbf{F}_m\|_{L^2(\mu)} \le \delta(1+m(p-1)) \, \|\mathbf{F}_m\|_{L^2(\mu)}. \label{eq:defect95bound95via95delta}\tag{47}\] If \(\delta(D)=0\) up to a certain degree, then the corresponding Lie coefficients are represented exactly, yielding an asymptotic exponent larger than 2. If \(\delta(D)\) is merely small, then the true asymptotic exponent remains \(2\) (unless the first defects vanish exactly), but pre-asymptotic windows with larger effective exponents may occur.

6.4 Fourier feature libraries↩︎

We now consider flow map approximation with a \(N\)-frequency Fourier feature library: \[\mathbf{g}^{(N)}(\mathbf{x})=\left\{\sin\!\bigl(\omega_0 \, \mathbf{k}\cdot \mathbf{x}\bigr),\cos\!\bigl(\omega_0 \, \mathbf{k}\cdot \mathbf{x}\bigr):\mathbf{k}\in \{0,\dots,N\}^d\right\}, \label{eq:Fourier95library}\tag{48}\] with base frequency \(\omega_0\). Such a library can show extended pre-asymptotic error scaling before crossing over to the asymptotic law \(\propto {\Delta t}^2\). We assume the data reside in a box \(\Omega \subset [-R,R]^d\), and define \[\rho := \omega_0 R.\] The relevant regime for approximating polynomial dynamics by low-frequency Fourier features is given by \(\rho \ll 1\).

Monomials are generated by frequency derivatives: \[\partial_{\mathbf{k}}^\alpha e^{\mathrm{i}\omega_0 \mathbf{k}\cdot \mathbf{x}} = (\mathrm{i}\omega_0)^{|\alpha|} \mathbf{x}^\alpha e^{\mathrm{i}\omega_0 \mathbf{k}\cdot \mathbf{x}}, \qquad \mathbf{x}^\alpha = x_1^{\alpha_1} \cdots x_d^{\alpha_d}, \label{eq:monomial95from95frequency95derivative}\tag{49}\] which gives at \(\mathbf{k}=0\): \[\mathbf{x}^\alpha = (\mathrm{i}\omega_0)^{-|\alpha|} \partial_{\mathbf{k}}^\alpha e^{\mathrm{i}\omega_0 \mathbf{k}\cdot \mathbf{x}}\big|_{\mathbf{k}=0}.\] A finite Fourier library may thus be viewed as a finite-difference approximation to these derivatives in frequency space. More explicitly, let \(a_{\alpha_j,\nu_j}\), \(\nu_j=0,\dots,N\), be the one-dimensional finite-difference weights for the \(\alpha_j\)-th derivative at \(0\) using the nodes \(0,\dots,N\), so that \[\sum_{\nu_j=0}^N a_{\alpha_j,\nu_j} \nu_j^m = \alpha_j! \, \delta_{m,\alpha_j},\qquad m=0,\dots,N.\] Setting \[c_{\alpha,\nu} := \prod_{j=1}^d a_{\alpha_j,\nu_j},\qquad \nu=(\nu_1,\dots,\nu_d) \in \{0,\dots,N\}^d,\] one obtains the Fourier approximant for the monomial \(\mathbf{x}^\alpha\), \[Q_{\alpha}^{\mathrm{Four}}(\mathbf{x}):= \frac{1}{(\mathrm{i}\omega_0)^{|\alpha|}} \sum_{\nu \in \{0,\dots,N\}^d} c_{\alpha,\nu} e^{\mathrm{i}\omega_0 \nu \cdot \mathbf{x}}. \label{eq:Fourier95monomial95approximant}\tag{50}\] A Taylor expansion of the exponential then yields the uniform bound \[\|Q_{\alpha}^{\mathrm{Four}} - \mathbf{x}^\alpha\|_{L^\infty(\Omega)} \le C_{\alpha,d,N} R^{|\alpha|} \rho^{\,N+1-\|\alpha\|_\infty}, \qquad \|\alpha\|_\infty \le N, \label{eq:universal95Fourier95bound}\tag{51}\] for a constant \(C_{\alpha,d,N}\) depending only on the indicated parameters. Thus every mixed monomial \(\mathbf{x}^\alpha\) whose largest coordinate exponent does not exceed \(N\) is approximated with an error suppressed by a power of \(\rho\). In particular, the approximation quality improves rapidly as either \(N\) increases or \(\omega_0 R\) decreases.

This mechanism is readily illustrated by the following low-order examples: \[x_j=\frac{\sin(\omega_0 x_j)}{\omega_0}+\mathcal{O}(\omega_0^2 x_j^3),\qquad x_j^2=\frac{2(1-\cos(\omega_0 x_j))}{\omega_0^2}+\mathcal{O}(\omega_0^2 x_j^4), \label{eq:low95order951d95fourier95examples}\tag{52}\] as well as, \[x_i x_j=\frac{\cos(\omega_0 x_i) + \cos(\omega_0 x_j) - 1 - \cos\!\bigl(\omega_0(x_i+x_j)\bigr)}{\omega_0^2}+\mathcal{O}(\omega_0^2 R^4). \label{eq:bilinear95fourier95example}\tag{53}\] Moreover, for even \(m\) and \(N \ge m/2\), there exist coefficients \(a_n^{(m)}\) such that \[x_j^m=\omega_0^{-m}\sum_{n=0}^N a_n^{(m)} \cos(n \omega_0 x_j)+\mathcal{O}\!\left( R^m \rho^{\,2N+2-m} \right),\] whereas, for odd \(m\) and \(N \ge (m+1)/2\), there exist coefficients \(b_n^{(m)}\) such that \[x_j^m=\omega_0^{-m}\sum_{n=1}^N b_n^{(m)} \sin(n \omega_0 x_j)+\mathcal{O}\!\left( R^m \rho^{\,2N+1-m} \right).\] These formulas show that pure powers may be easier to approximate than general mixed monomials.

For a polynomial vector field, the Lie coefficients \(\mathbf{F}_m\) are finite sums of monomials. Let \[\eta_m := \max_{\alpha \in \mathrm{supp}(\mathbf{F}_m)} \|\alpha\|_\infty \label{eq:max95coord95exponent}\tag{54}\] denote the maximum coordinate exponent of the monomials in \(\mathbf{F}_m\). Then 51 implies that, whenever \(\eta_m \le N\), \[\varepsilon_{m} = \|\mathcal{Q}\mathbf{F}_m\|_{L^2(\mu)} \le C_{m,N} \rho^{\,N+1-\eta_m},\qquad (\eta_m \le N) \label{eq:Fourier95defect95bound}\tag{55}\] for a suitable constant \(C_{m,N}\). This shows that low-order Lie coefficients are not represented exactly, but their defects are strongly suppressed by a positive power of \(\rho\). Consequently, for fixed \(N\), a finite Fourier library generically satisfies \(\varepsilon_{1}>0\), so the true asymptotic law remains \[\mathcal{E}({\Delta t}) \sim \varepsilon_1^2 {\Delta t}^2.\qquad \text{(asymptotically {\Delta t}\to 0)}\]

The pre-asymptotic regime is determined by the first Lie coefficient whose degree leads an unsuppressed defect (or, at least, where 55 does not guarantee suppression). This corresponds to a non-positive power of \(\rho\) in 55 , implying a conservative bound \[m_*(N) := \min \{ m \ge 1 : \eta_m > N \}. \label{eq:Fourier95effective95mstar}\tag{56}\] Over the \({\Delta t}\)-range in which the defects of orders \(m < m_*(N)\) remain suppressed relative to the \(m_*(N)\)-th defect, one thus expects \[\mathcal{E}_N({\Delta t}) \approx A_N {\Delta t}^{2m_*(N)}, \label{eq:Fourier95preasymptotic}\tag{57}\] with a constant \(A_N\sim \varepsilon_{m_*(N)}^2/(m_*(N)!)^2\). For generic degree-\(p\) polynomial dynamics one has \(\eta_m \le 1+m(p-1)\), which yields the conservative estimate \[m_*(N) \geq 1+ r_*,\qquad r_* = \left\lfloor \frac{N-1}{p-1} \right\rfloor. \label{eq:Fourier95preasy95exp}\tag{58}\] Because \(\varepsilon_1\) is not exactly zero, this pre-asymptotic regime must ultimately cross over to the true \({\Delta t}^2\) asymptotics at a scale \[{\Delta t}_{\times} \asymp \left( \frac{m_*(N)! \, \varepsilon_{1}}{\varepsilon_{m_*(N)}} \right)^{1/(m_*(N)-1)}.\label{eq:Fourier95crossover}\tag{59}\] Note that in practice, empirically observed values of \(m_*(N)\) can be larger since bounds like 51 are not necessarily tight.

6.5 Connection to bias–variance decomposition↩︎

The projection error 35 is the approximation-bias component of the usual bias–variance decomposition. In the increment-prediction setting of 6.1, the target is \(\mathbf{F}_{\Delta t}\in L^2(\mu;\mathbb{R}^d)\) and the model class is \[V = \{ H \mathbf{g}(\cdot) : H \in \mathbb{R}^{d \times n} \}.\] Let \(\mathbf{h}_\star := \mathcal{P}\mathbf{F}_{\Delta t}\) denote the \(L^2(\mu)\)-orthogonal projection of \(\mathbf{F}_{\Delta t}\) onto \(V\). Then, for any learned predictor \(\hat{\mathbf{h}}\in V\), \[\|\mathbf{F}_{\Delta t}-\hat{\mathbf{h}}\|_{L^2(\mu)}^2 = \underbrace{\|\mathbf{F}_{\Delta t}-\mathbf{h}_\star\|_{L^2(\mu)}^2}_{=\mathcal{E}({\Delta t})=\text{model bias}^2} + \underbrace{\|\hat{\mathbf{h}}-\mathbf{h}_\star\|_{L^2(\mu)}^2}_{\text{estimation error}},\] because \(\mathbf{F}_{\Delta t}-\mathbf{h}_\star = \mathcal{Q}\mathbf{F}_{\Delta t}\perp V\). Averaging over training sets decomposes the second term into the usual finite-sample estimator bias and variance. For the noiseless, overdetermined, fixed-library large-sample regime considered here, this estimation term vanishes, leaving only the approximation bias. The Lie-defect expansion in 38 40 is therefore an asymptotic expansion of this model-bias term.

6.6 Truncated flow-map error↩︎

Since the exact expression for the flow map \(\mathbf{F}_{\Delta t}\) is usually not available, consider its approximation to order \(r\) \[\mathbf{F}_{\Delta t}^{[r]}(\mathbf{x}):=\Phi_{\Delta t}^{[r]}(\mathbf{x})-\mathbf{x}= \sum_{m=1}^{r}\frac{\Delta t^m}{m!}\mathbf{F}_m(\mathbf{x}),\] and define the difference to the true flow map as \(\mathbf{R}_{\Delta t}^{[r]}:=\mathbf{F}_{\Delta t}-\mathbf{F}_{\Delta t}^{[r]}\). Consider now the difference between \(\mathbf{F}_{\Delta t}^{[r]}\) and the approximation of the true flow map by the feature library, \(\mathcal{P}\mathbf{F}_{\Delta t} \approx \hat{H} \mathbf{g}\): in the same approximation-dominated regime as in 35 , one has \[\begin{align} \mathcal{D}_{\Phi,r}(\Delta t) &= \left\| \mathcal{P}\mathbf{F}_{\Delta t} - \mathbf{F}_{\Delta t}^{[r]}\right\|_{L^2(\mu)}^2 \\ &= \left\|-\mathcal{Q}\mathbf{F}_{\Delta t}^{[r]}+\mathcal{P}\mathbf{R}_{\Delta t}^{[r]} \right\|_{L^2(\mu)}^2 \\ &= \left\|\sum_{m=1}^{r}\frac{\Delta t^m}{m!}\mathbf{q}_m\right\|_{L^2(\mu)}^2+\left\| \sum_{m=r+1}^{\infty} \frac{\Delta t^m}{m!}\mathcal{P}\mathbf{F}_m \right\|_{L^2(\mu)}^2 , \label{eq:truncated95flowmap95decomp} \end{align}\tag{60}\] where the last equality uses the orthogonality of \(V\) and \(V^\perp\). Thus the flow-map approximation error consists of an unrepresented part of the retained Lie coefficients and the represented part of the omitted Taylor remainder.

Let \(m_0\) be the first non-vanishing defect order defined in 39 , and set \[s_r:=\min\left\{m>r: \|\mathcal{P}\mathbf{F}_m\|_{L^2(\mu)}>0 \right\}.\] Then, generically, \[\mathcal{D}_{\Phi,r}(\Delta t) \sim \begin{cases} \dfrac{\varepsilon_{m_0}^2}{(m_0!)^2}\Delta t^{2m_0}, & m_0\le r,\\[1.2em] \dfrac{\|\mathcal{P}\mathbf{F}_{s_r}\|_{L^2(\mu)}^2}{(s_r!)^2} \Delta t^{2s_r}, & m_0>r . \end{cases} \label{eq:truncated95flowmap95scaling}\tag{61}\] The relation to the training error [35 ] follows from the triangle inequality: \[\left| \mathcal{D}_{\Phi,r}^{1/2}(\Delta t) - \mathcal{E}^{1/2}(\Delta t) \right| \le \|\mathbf{R}_{\Delta t}^{[r]}\|_{L^2(\mu)} =\mathcal{O}(\Delta t^{r+1}). \label{eq:truncated95flowmap95training95bound}\tag{62}\] Consequently, if the truncation order includes the first unresolved Lie coefficient, \(r\ge m_0\), then the truncated flow-map error is asymptotically equivalent to the training error: \[\mathcal{D}_{\Phi,r}(\Delta t) \simeq \mathcal{E}(\Delta t) .\] If \(r<m_0\), however, the error is controlled by the Taylor terms omitted from the reference map.

6.7 Semigroup consistency↩︎

The exact flow map satisfies the semigroup property [29] \[\Phi_{s+t}=\Phi_s\circ \Phi_t .\] This can be tested for the learned flow map \(\widehat{\Phi}_h\), and we show that it essentially follows the same \({\Delta t}\)-scaling as the training error. In the population limit, \[\widehat\Phi_h(\mathbf{x}):=\mathbf{x}+\mathcal{P}\mathbf{F}_h(\mathbf{x}),\] the \(\ell\)-step semigroup defect is given by \[\mathcal{S}_\ell(h)^2:=\left\|\widehat\Phi_{\ell h}-\widehat\Phi_h^\ell\right\|_{L^2(\mu)}^2 , \label{eq:semigroup95defect}\tag{63}\] where \(\widehat\Phi_h^\ell\) denotes \(\ell\)-fold composition. Since \(\Phi_{\ell h}=\Phi_h^\ell\) for the exact flow, it follows that \(\widehat\Phi_{\ell h}-\widehat\Phi_h^\ell = (\widehat\Phi_{\ell h}-\Phi_{\ell h}) + (\Phi_h^\ell-\widehat\Phi_h^\ell)\) and the semigroup defect can be bounded as \[\mathcal{S}_\ell(h) \leq \sqrt{\mathcal{E}(\ell h)}+C_\ell \sqrt{\mathcal{E}(h)} \label{eq:semigroup95bound95training}\tag{64}\] in terms of the increment-prediction training error [35 ], \(\mathcal{E}(h)=\|\mathcal{Q}\mathbf{F}_h(\mathbf{x})\|_{L^2(\mu)}^2 =\|\Phi_h(\mathbf{x})-\widehat\Phi_h(\mathbf{x})\|^2\). Here, \(C_\ell\) depends on Lipschitz constants of the flow, assuming stability along intermediate trajectories.

If the first unresolved Lie coefficient occurs at order \(m_0\geq 2\), i.e. \(\mathbf{q}_1=\cdots=\mathbf{q}_{m_0-1}=0\), \(\mathbf{q}_{m_0}\neq 0\), then, for fixed \(\ell\), \[\widehat\Phi_{\ell h}-\widehat\Phi_h^\ell=-\frac{\ell^{m_0}-\ell}{m_0!} h^{m_0}\mathbf{q}_{m_0} + \mathcal{O}(h^{m_0+1}) \label{eq:multistep95semigroup95leading}\tag{65}\] Combining this with 40 , i.e., \(\mathcal{E}(h)=\frac{h^{2m_0}}{(m_0!)^2} \|\mathbf{q}_{m_0}\|_{L^2(\mu)}^2 + \mathcal{O}(h^{2m_0+1})\), gives \[\mathcal{S}_\ell(h)^2=(\ell^{m_0}-\ell)^2\mathcal{E}(h)+\mathcal{O}(h^{2m_0+1}). \label{eq:semigroup95training95relation}\tag{66}\] Thus, whenever the vector field and the early Lie coefficients are exactly represented, the semigroup defect has the same asymptotic exponent as the training error. In particular, for a polynomial vector field and a polynomial library, \(m_0=r^\ast+1\) with \(r^\ast\) given by 45 .

By contrast, if the vector field is not exactly representable, \(\mathbf{q}_1\neq 0\), then the \(\mathcal{O}(h)\) terms cancel in \(\widehat\Phi_{\ell h}-\widehat\Phi_h^\ell\), such that the semigroup defect starts at order \(h^2\). Consequently, semigroup consistency is in this case at best an internal diagnostic of whether the learned map is has flow-map character.

7 Remarks on the linear delay extrapolation stencil↩︎

A linear delay model is defined by \[\boldsymbol{\psi}((n+1){\Delta t}) = H_\mathrm{ld}\boldsymbol{\psi}(n{\Delta t}), \label{eq:app95delay95statespace}\tag{67}\] with the delay state \[\boldsymbol{\psi}(n {\Delta t}) = \begin{bmatrix} \mathbf{x}(n{\Delta t})\\ \mathbf{x}((n-1){\Delta t})\\ \vdots\\ \mathbf{x}((n-\Gamma+1){\Delta t}) \end{bmatrix} \in \mathbb{R}^{d \Gamma}\] and the companion (delay-shift) matrix \[H_\mathrm{ld} = \begin{bmatrix} \zeta_0 \mathbb{I}_d & \zeta_1 \mathbb{I}_d & \cdots & \zeta_{\Gamma-1} \mathbb{I}_d\\ \mathbb{I}_d & 0 & \cdots & 0\\ 0 & \mathbb{I}_d & \cdots & 0\\ \vdots & & \ddots & \vdots\\ 0 & \cdots & \mathbb{I}_d & 0 \end{bmatrix} \in \mathbb{R}^{d \Gamma\times d \Gamma}. \label{eq:app95Hst}\tag{68}\] For general coefficients \(\zeta_j\), the block companion matrix has characteristic polynomial \[\chi_{H_\mathrm{ld}}(\lambda) = \left( \lambda^\Gamma-\sum_{j=0}^{\Gamma-1}\zeta_j\lambda^{\Gamma-1-j} \right)^d .\] For the polynomial extrapolation stencil [24 ] \(\zeta_j=(-1)^j\binom{\Gamma}{j+1}\), this reduces to \(\chi_{H_\mathrm{ld}}(\lambda)=(\lambda-1)^{d\Gamma}\), which implies a single eigenvalue \(\lambda=1\) with algebraic multiplicity \(d \Gamma\). A \(d=1\) harmonic signal \(x(t) = \sum_{k=1}^N (a_k \cos(k\omega_0 t) + b_k \sin(k\omega_0 t))\) can be modeled with \(\Gamma= 2N\) time taps [68]–[70]. The \(\zeta_k\) are then given by the coefficients of the characteristic polynomial \(\prod_{k=1}^N(\lambda^2-2\cos(k\omega_0\Delta t)\lambda+1)\) expanded in \(\lambda\) [105], [106].

8 Training error and forecast horizon↩︎

We review here how the training error of an NG-RC-type model trained on dynamical system trajectories \(x(t)\) affect the forecast horizon. This is established theory in numerical analysis of time series [107]–[109], recast here in the language of NG-RC.

8.1 Residual error and model bias↩︎

Consider the flow map \(\Phi_{\Delta t}(x)\), acting as \[x(t+\Delta t)=\Phi_{\Delta t}\bigl(x(t)\bigr),\] and its approximation by a learned next-step map \[\hat{x}(t+\Delta t)= \hat{H}\,g\bigl(\hat{x}(t)\bigr),\] where \(g=(x_1,\ldots,x_d,\psi_1,\ldots,\psi_m)\in \mathbb{R}^n\) is a reservoir vector. Consider the residual \[R(x) \equiv \Phi_{\Delta t}(x) - \hat{H}\,g(x),\] which we assume to be Lipschitz continuous with some constant \(L_R\), i.e., \(\|R(x)-R(y)\| \le L_R \|x-y\|\). On the training data, we assume the uniform bound \[\|R(x)\| \le \kappa \varepsilon_{\mathrm{train}}, \qquad (\kappa\sim \mathcal{O}(1)) \label{eq95trainerr95bound}\tag{69}\] where we assume the residual is sufficiently concentrated such that it can be bounded by the training error up to a factor of \(\mathcal{O}(1)\) (which we neglect henceforth). Denoting by \(\hat{x}\) a point outside the training data (generated, e.g., by the autonomous forecast), it follows from the triangle inequality that \[\|R(\hat{x})\| \le \|R(x)\| + L_R \|x-\hat{x}\| \le \varepsilon_{\mathrm{train}} + L_R \|x-\hat{x}\|. \label{eq95residual95bound}\tag{70}\] The second term on the r.h.s.captures the effect of model misspecification (bias), which happens if the trained model is structurally different from the true model.

8.2 Error recurrence in autonomous forecasting↩︎

We study the error evolution during autonomous forecasting. Let \(x_{n+1}=\Phi_{\Delta t}(x_n)\) be the true dynamics and \(\hat{x}_{n+1}=\hat{H}\,g(\hat{x}_n)\) be the autonomous prediction obtained from the trained model. Denoting the deviation between the two by \(e_n=\|x_n-\hat{x}_n\|\), we have \[e_{n+1} =\|\Phi_{\Delta t}(x_n)-\hat{H}\,g(\hat{x}_n)\| \le \|\Phi_{\Delta t}(x_n)-\Phi_{\Delta t}(\hat{x}_n)\| +\|\Phi_{\Delta t}(\hat{x}_n)-\hat{H}\,g(\hat{x}_n)\|. \label{eq95trainerr95evol}\tag{71}\] Lipschitz continuity of the flow map implies \(\|\Phi_{\Delta t}(x_n)-\Phi_{\Delta t}(\hat{x}_n)\|\le L_\Phi e_n\). In general, for small \(\Delta t\), one can express the Lipschitz constant as \(L_\Phi\approx 1+ \Delta t \|\mathcal{D}f\|\) in terms of the Jacobian \(\mathcal{D}f\) of the dynamical system. We now use the residual bound of 70 in the second term in 71 , and assume the model generalizes well, such that the one-step error \(\|R(x_n)\|\) on the true trajectory is also bounded by \(\varepsilon_{\mathrm{train}}\). This yields the recurrence relation \[e_{n+1}\le \underbrace{(L_\Phi+L_R)}_{L} e_n+\varepsilon_{\mathrm{train}}, \label{eq95err95recur}\tag{72}\] where we introduced an effective Lipschitz constant \(L\) that accounts for both the flow map approximation error and model misspecification.

8.2.1 Stable dynamics↩︎

For non-chaotic (precisely, non-expansive) systems, we take \(L_\Phi \le 1\), while the error may still grow exponentially if \(L_R\) is large. If the system is absolutely stable (\(L < 1\)), iterating 72 leads to the error bound \(e_n \le L^n e_0 + \varepsilon_{\text{train}} (1 - L^n)/(1 - L)\). As \(n \to \infty\), the error converges to a finite value \(e_\infty = \varepsilon_{\text{train}} / (1 - L)\). If this asymptotic error is smaller than some given tolerance \(E_{\text{tol}}\), the forecast horizon can be considered to be effectively infinite. If \(e_\infty > E_{\text{tol}}\), the horizon is finite but typically determined by the initial transient rather than the accumulation of \(\varepsilon_{\text{train}}\).

In the case of marginal stability (\(L = 1\), e.g., for conservative or periodic systems with \(L_\Phi=1\) and \(L_R=0\)), the recurrence bound becomes \(e_{n+1} \le e_n + \varepsilon_{\text{train}}\). Iterating this bound suggests a worst-case scenario where errors accumulate coherently, leading to linear error growth: \(e_n \le e_0 + n \varepsilon_{\text{train}}\). If the forecast horizon \(T_{\text{fch}} = n \Delta t\) is defined by when \(e_n\) reaches \(E_{\text{tol}}\), this implies \[T_{\text{fch}} \approx \Delta t (E_{\text{tol}} - e_0) / \varepsilon_{\text{train}}.\] By contrast, assume now that the one-step prediction errors \(\delta_n = \Phi_{\Delta t}(\hat{x}_n) - \hat{H} g(\hat{x}_n)\) (with \(\|\delta_n\| \le \varepsilon_{\text{train}}\)) accumulate incoherently, akin to steps in a random walk—a useful heuristic if \(\varepsilon_{\text{train}}\) originates from observation noise that decorrelates sufficiently fast. Then the expected squared error grows linearly with time: \(\mathbb{E}[\|e_n\|^2] \approx \|e_0\|^2 + n \mathbb{E}[\|\delta_k\|^2] \approx \|e_0\|^2 + n \varepsilon_{\text{train}}^2\) (interpreting \(\varepsilon_{\text{train}}^2\) as the mean squared one-step error). In this scenario, setting the error magnitude to the tolerance \(E_{\text{tol}}\) yields a forecast horizon \[T_{\text{fch}} \approx \Delta t (E_{\text{tol}}^2 - \|e_0\|^2) / \varepsilon_{\text{train}}^2, \label{eq95pred95hor95incoh}\tag{73}\] which scales as \(\varepsilon_{\text{train}}^{-2}\). The relevant scaling depends on the nature of the approximation error and how it propagates in the specific neutrally stable system.

8.2.2 Chaotic dynamics forecast horizons↩︎

We finally consider the case \(L>1\), which not only includes chaotic dynamics, but also exponential error growth due to large \(L_R\). For chaotic dynamics, one has \(L_\Phi \simeq e^{\lambda_1 \Delta t}\) with \(\lambda_1>0\) the largest Lyapunov exponent. Accordingly, iterating 72 over \(n\) steps gives \[e_n \lesssim L^n e_0 + \varepsilon_{\text{train}} \sum_{k=0}^{n-1} L^k, \label{eq95errprop95sum}\tag{74}\] where \(e_0\) is introduced as the error representing the mismatch with the exact initial condition. Evaluating the geometric series in 74 as \(\sum_{k=0}^{n-1} L^k = \frac{1 - L^n}{1 - L} \approx \frac{L^n}{L - 1}\) for large \(n\), we obtain \[e_n \lesssim L^n \left( e_0 + \frac{\varepsilon_{\text{train}}}{L - 1} \right).\] Both an imperfect initial condition and a nonzero training error feed into the exponential divergence between true and predicted dynamics.

Defining a tolerance \(\varepsilon_{\mathrm{tol}}\) for the accumulated error, such that \(e_n \lesssim \varepsilon_{\mathrm{tol}}\) signifies good forecasting accuracy, one obtains the prediction horizon \[T_{\text{fch}}=n \Delta t \simeq \frac{1}{\lambda} \ln \left(\frac{\varepsilon_{\text{tol}}}{e_0+ c\varepsilon_{\text{train}}}\right), \label{eq95forec95hor}\tag{75}\] with the rate constant (effective Lyapunov exponent) \(\lambda\equiv {\Delta t}^{-1} \ln L\) and the constant \(c\equiv 1/(L-1)\). As an example, if \(L\simeq e^{\lambda_1{\Delta t}}\) with \(\lambda_1\sim\mathcal{O}(1)\) and \({\Delta t}\sim\mathcal{O}(10^{-1}-10^{-3})\), then \(c=(L-1)^{-1}\sim\mathcal{O}(10-1000)\). Taking \(\varepsilon_{\text{tol}}\sim\mathcal{O}(10^{-1})\) and \(\varepsilon_{\text{train}}\sim\mathcal{O}(10^{-2}-10^{-7})\), the dimensionless horizon \(\lambda T_{\text{fch}}\) ranges from values close to zero, when \(c\varepsilon_{\text{train}}\gtrsim\varepsilon_{\text{tol}}\), up to \(\mathcal{O}(10)\) for the smallest training errors, neglecting the influence of the initial condition, i.e., \(e_0\approx 0\). 12 illustrates the relation between the training error and the forecast horizon for a well-specified (bias-free) model.

References↩︎