When are dynamical systems learned from time series data statistically accurate?

Results require full Jacobian but it can be approximated with available derivative information in practice

You are absolutely right in that the Jacobian is computationally hard to calculate (through AD or finite difference) in high dimensions; thank you and Reviewer vAp3 for pointing this out! This in turn would slow the training process and might make it prohibitive for learning high-dimensional chaotic systems. We emphasize that our primary contribution is toward the foundational question of understanding generalization in the context of learning chaotic dynamics. Our main result deconstructs the effect of the ergodic properties of the true system on the generalization ability of the neural network. To isolate this effect, we consider only a minimal illustrative setting of regression using the common MSE loss. When some Jacobian information (local information about the short-term dynamics and its linear perturbation structure) is added to the loss, surprisingly we find that the physical measure (global statistics) is learned. Our results are aimed at explaining this observation via the mechanism of shadowing. As a minimal setting to verify our results in practice without any confounders, we consider systems in which the full Jacobian can be calculated and used in training. Having said that however, our results lend themselves to many practical training strategies where the Jacobian information is only approximated. For instance, in scenarios where adjoints or other surrogate models for the Jacobian are available, which are common in many scientific applications in the geosciences and aerospace turbulent flows, these can be used as proxies for the full Jacobian matrix. As we mention in response to Reviewer vAp3, who raised the same concern, we are interested in analyzing whether random Jacobian matrix-vector products can yield informative generalization bounds such as Theorem 2. Since the novelty and focus of our present paper is in the theoretical connection between using Jacobian information and statistical accuracy and why this obviates the need for performing generative modeling, we defer such analyses to future work.

Focus on autonomous systems

In this paper, we focus on autonomous systems -- wherein the flow map itself does not depend on time -- which occupy a large class of chaotic systems encountered in nature. As you correctly point out, the analysis and even the background, starting from the existence of a unique physical measure, will completely change when we consider nonautonomous systems such as chaotic systems with control terms and random dynamical systems with chaotic orbits. For these systems, there is a lot of exciting recent development in the theoretical dynamical systems literature (see e.g. this survey article: https://link.springer.com/article/10.1007/s10955-016-1639-0) proving the existence of random physical measures. The next step would be connect the theory of such measures with statistical learning theory. In the present work, we show a path for such a connection in deterministic chaotic systems, deferring random chaotic systems (which are increasingly recognized as being useful models for stochastically parameterized climate systems) to future work. We will add this remark to the appendix in the next revision. Thank you very much for bringing this to our attention!

Related work on invariant measures

Thank you for sharing these references -- we have added them to Related Work, which has spilled into the Appendix! Since the focus is on deriving dynamics-aware generalization bounds, we consider a minimal regression setup for one-step dynamics, which does not require specialized training methods such as the one introduced in [1]. In future work, it will be interesting to derive bounds similar to Theorem 2 for the sparse teacher forcing in [1] and other RNN/ESN-based training methods, as you point out. Many thanks for reference [3] that we have been remiss to overlook. Our setting is that of unique, ergodic, invariant physical measures in this paper, which is less complicated than multiple ergodic basins considered in [3]. It is reassuring that generalization with MSE loss does not imply statistical accuracy for the non-ergodic case considered in Sec 3.2 of [3]. Another interesting connection would be to specialize notions of generalization (Definition 3.2 specifically) introduced in [3] to the unique SRB/Gibbs measure case. While this notion or the W1 distance can be used to derive ERM problems, the sample complexity of these problems may be large. From the theoretical standpoint, we will extend our shadowing-based results to the non-ergodic case by considering weaker notions of finite-time shadowing in every connected component of the attractor. Different from [3], we demonstrate that even when we study ergodic systems, and sample from the entire state space, the traditional MSE does not yield the correct notion for generalization.

Practicality of Jacobian loss and implications of theory

This is indeed an important point, thank you and Reviewer PvwR for pointing this out! The computation of the Jacobian matrix estimated via automatic differentiation or via finite differences scales quadratically with the problem dimension. However, as illustrated in our numerical results, the sample complexity is indeed much smaller than for directly learning the physical measure via a generative model. Moreover, since the training strategy is an elementary regression on any simple architecture (like an MLP), it can be competitive when compared to more sophisticated training algorithms with recurrent architectures or transformers. In practice, whether to choose simple regression or to resort to transformers or recurrent architectures with innovative training stabilization strategies will depend on the dimension and the dynamical complexity (attractor dimension, number of positive Lyapunov exponents, correlation decay time etc). The computational complexity of complicated models such as transformers during inference time is quadratic in the dimension as well, making an elementary architecture appealing as a surrogate model.

Furthermore, our result that adding Jacobian information implies statistical accuracy invites the development of practical supervised learning methods in high dimensions. Even though the theoretical results are derived for $C^1$ generalization, we are currently working on incorporating random Jacobian-vector products into the loss function. Such matrix vector products are easily tractable in several scientific applications where adjoint/tangent equation codes are available (see e.g geophysical turbulence models). One subject of our future work is to evaluate the efficacy of black box Jacobian-vector products in learning the physical measure.

Unrolling dynamics

We have run experiments on the Lorenz '63 system with an unrolled loss function: $\ell_{\rm u}(x) := (1/k) \sum_{t \leq k} \|v_t(x)\|^2,$ where $v_t(x) = F^t_{\rm nn}(x) - F^t(x)$ are vector fields on learned orbits. With $k = 10$ timesteps of unrolling time, we see that the attractor is reproduced well, but atypical orbits are still produced for random initializations. The learned LEs are still incorrect: $[ 0.905, -0.03, -5.6]$ (the true LEs are :$[0.85, 0, -14.5].$) When the unrolling time $k$ is increased, we still obtain the same results, with even the positive LE being overestimated.

Overall, unrolling seems to learn more accurate representations than the MSE model but worse than the JAC model. We also observe that the unrolling time $,k,$ needs to be fine-tuned as a hyperparameter to achieve good generalization; a small perturbation in $k$ can lead to training instabilities. We have these results to Appendix C of our revision.

To understand these results, for short times, $t,$ (when compared to the Lyapunov time, $1/\lambda_\mathrm{max}),$ $v_t$ are in tangent spaces along the NN-orbit, $\{F^t_\mathrm{nn}(x).\}$ Now, we note that $v_t$ is the pushforward through the Jacobian/linearized dynamics of $v_{t-1},$ which yields the recursive relationship, $v_t(x) = dF(F^t(x)) v_{t-1}(x) + v_0(F_{\rm nn}^t(x)).$ Thus, the unrolled loss does contain Jacobian information implicitly although it does not enforce the learned trajectory to be close to the true trajectory in $C^1$-distance. This serves to explain, via Theorem 2, why it works better than the MSE (even in some practical climate emulators, e.g., FourcastNet). Thank you very much for raising this important point, which supports the central argument of this paper: adding Jacobian information leads to statistical accuracy.

W1 distance

We use Python Optimal Transport package, which uses the Sinkhorn algorithm. For our 1D distributions, the optimal transport map is analytically the target inverse CDF $\circ$ CDF of source, which is estimated via projected gradient descent in the package.

Edits

Many thanks for a careful reading and suggesting many improvements (implemented in the revision, including adding limitations to sec 7) to the paper! We have rewritten ``Dynamic generative models'': here, the dynamics $\Phi^t_\phi$ is stochastic, while the encoder and decoder are deterministic. This is needed for greater expressivity of the conditional measures.

When Assumption 1 is expected to hold

A necessary condition for Assumption 1 ($\mathcal{C}^1$ strong generalization) is that the optimization problem with the Jacobian loss is solved "well" -- resulting in low generalization errors. But, as you have carefully observed, this is insufficient to claim this notion of strong generalization. The Jacobian loss must be small along orbits of the learned NN (NN with small Jacobian generalization error). Note that these orbits are $(\epsilon_0, \epsilon_1)$ orbits of the true system. When these orbits are on or near the support of $\mu$ (the physical measure of the true system), then, the Jacobian loss is small at those points, and $\mathcal{C}^1$ strong generalization holds. Thus, a sufficient condition is for the true dynamics to be such that small $\mathcal{C}^1$ perturbations of it lead to small perturbations of the physical measure. Such a condition is called smooth linear response in the parlance of dynamical systems and has been extensively studied in the theoretical literature (e.g., see the review by Baladi here: https://arxiv.org/abs/1408.2937). For uniformly hyperbolic systems, which are considered in this paper, and many chaotic systems observed across physics, such a smooth linear response holds (https://iopscience.iop.org/article/10.1088/0951-7715/22/4/009/meta). But, of course, this is only a sufficient condition, and in practice, Assumption 1 may be satisfied more easily. For instance, when we train by sampling points randomly in a box around the attractor, we automatically sample from points near the support of $\mu$ (the attractor). Thus, in input space regions where orbits of $F_{\rm nn}$ live, the Jacobian loss might evaluate to a small value. That is, whether or not linear response holds, we achieve $\mathcal{C}^1$ strong generalization by choosing training points randomly distributed (according to any density) around the attractor and minimizing the Jacobian loss. We thank you for this excellent question, which has led to an important clarification in the paper. Due to space constraints in the paper, we have added more details in the appendix and only briefly make a clarification after Assumption 1 in the main text.

Why minimizing MSE loss does not yield statistical accuracy

This is an astute observation, thank you. Indeed the fundamental contribution of the paper is to suggest shadowing as a mechanism for learning the physical measure. Shadowing does not hold when we only learn well with respect to the $\mathcal{C}^0$ loss! This is because in order to shadow a chaotic orbit, informally, we need to predict the next step as well as local directions of infinitesimal linear perturbation (linear/Jacobian structure) induced by the next step dynamics. This intuition is implicit in the proof of the shadowing lemma (see e.g., Chapter 18 of Katok and Hasselblatt), and therefore extends to the proof of our high-probability version of the shadowing lemma that leads to our main result. Just learning the short-term dynamics without learning the local directions of growth/decay of infinitesimal perturbations cannot lead to learning shadowing orbits. Since we prove that shadowing is the underlying mechanism for learning the invariant measure, the MSE loss is insufficient. We hope this clarifies our result -- thank you for the great question! Due to space constraints, we have not added this explanation, but we will in the next revision right after Theorem 2.

Neural ODEs are neural parameterizations of the vector field

Neural ODEs (trained with MSE) learn the vector fields, but we need to learn the derivatives of the vector fields (with respect to the state vector) for statistical accuracy. Our results indeed already compare Neural ODEs where the vector field is parameterized by an MLP, MLP with Fourier Layers and ResNet blocks. Our observation is the same: without the Jacobian information in the loss function, even parameterizing the vector field, like a Neural ODE does, is not sufficient to learn physical measure.

Intuitively, one can argue that learning more and more derivatives of a function (like a vector field of a flow or a discrete-time dynamical system, $F$) would lead to more accurate learning the function and hence also give rise to statistical accuracy. But shadowing (which only requires first derivatives) says that higher-order derivative matching is not necessary.

Limitations

We add the discussion on limitations in Section 7 of the revision. Mainly, we illustrate our theoretical results using the full Jacobian, but this does not yield a practical scheme for training. Secondly, our results are derived for mathematically ideal chaotic systems trained with a vanilla regression setup. Training interventions that enforce learning invariant measures are not considered in our analysis. In this regard, please see our responses to Reviewer vAp3 and sRDx.

Implications for training with Sobolev norm

It is indeed interesting to consider errors in the learned dynamics as distributions in a Sobolev space, $W^{k,p}$, as you point out. You are absolutely correct that our MSE loss is a special case for $k = 0, p =2$ and the Jacobian loss for $k= 1, p =2.$ Here, we only consider classical functions (as opposed to distributions) for the true map $F$; the neural network, $F_{\rm nn},$ also does not have singularities since we use a smoothed out ReLU as activation.

Hence, the weak derivatives above are indeed just classical derivatives. Our results (Theorem 2) means the following for the Sobolev loss: we do not gain more statistical accuracy by minimizing errors with larger $k,$ that is, by considering higher-order derivatives. When $p=2,$ we do not gain by including in the error or loss term higher order Fourier coefficients. In this paper, we prove that a sufficient condition for statistical accuracy is the prevalence of shadowing and its typicality (that is, distribution of shadowing orbits according to $\mu$ with high probability). Our high-probability version of shadowing only requires $C^1$ strong generalization: that is, matching derivatives only up to first order. Therefore, when a shadowing based mechanism for generating the physical measure holds in a supervised learning problem, our results imply that matching higher-order derivatives are not necessary. We thank you for making this important suggestion! Due to space constraints, we have not added this to the revision but plan to add a more detailed analysis to the next revision.

Related work

These are valuable additions to our related work, we are greatly indebted to you and Reviewer vAp3 for showing us these references! We have added all of the references in the Related Work and continued the section into an appendix due to lack of space. Regarding the first cited reference, it is very interesting to derive generalization bounds for the dissipativity-enforcing loss proposed by the authors. Without enforcing dissipativity, in our work, we find that incorporating the Jacobian learns an attractor of zero volume (the attractor dimension is related to the LEs and these are correctly obtained) in dissipative systems. We do not consider training interventions as done in [2] and [4] to stabilize the training of chaotic dynamics over longer time horizons. We leave performing theoretical analysis of the statistical accuracy of the empirical approach to train RNNs in [4] for future research. Instead, our focus here is on obtaining provable guarantees for learning of the physical measure; hence, we only consider here the minimal training setting of regression over short time horizons, wherein such stabilization interventions are not needed. The reference [3] is indeed very close in motivation to our work in that invariant measures of chaotic systems are learned by neural network models. However, the approach taken in [3] is markedly different since optimal transport is performed on key statistics or such summary statistics are approximated through contrastive learning. We argue that supervised learning with Jacobian information is less complicated and more tractable in high dimensions even when compared to an efficient OT algorithm (Sinkhorn -- on discrete measures -- is used in the reference). Furthermore, we obtain theoretical guarantees for the error in the learned invariant measure (e.g. in Wasserstein distance), while it would be interesting to derive such guarantees for the contrastive learning approach taken in [3].

Clarification of Figure 1

The y-axis of the 4th and 5th columns of Figure 1 were the empirical PDFs of a random orbit generated by the MSE and Jacobian models respectively. The different scale of the plots show that the empirical distribution of the MSE orbit is incorrect. This exemplifies the thesis of the paper: adding Jacobian information in the regression problem leads to learning the physical measure. To avoid confusion, we have now combined the 4th and 5th plots into one figure, revised the caption and added a Gist (please see the link). Thank you!

Clarifying notation of map

We apologize for the confusion: the h is used as an argument for writing the generalization error as a function of the model $h$. The learned chaotic map is denoted throughout as $F_{\rm nn}$ and the true map by $F$.

Solving OT problems with ergodic measures of chaotic maps

Thank you for the careful observation! Suppose we are trying to minimize the Wasserstein distance, whose dual form we can lower bound, for some $g \in {\rm Lip}^1(M),$
as

for Leb a.e. x. In practice, e.g., when solving OT problems with the Sinkhorn algorithm, we replace the continuous measure with a discrete measure, which in the above case reduces to replacing the ergodic average ($t\to \infty$) with an average over an orbit of finite length. Then, the estimate of the above error is noisy when $F$ and $F_{\rm nn}$ are chaotic, and indeed the variance grows exponentially with $t$ and then saturates (due to both attractors being bounded).

LEs are statistical measures

LEs are indeed long-term evaluation metrics (now added in Section 2). For a vector field $E_i$ in the $i$th Oseledets subspace, we can rewrite the definition of the $i$th LE as an ergodic average: $\lambda_i := \lim_{t \to \infty} \dfrac{1}{t} \sum_{n =1}^t \log \|dF (F^n(x)) E_i(F^n(x))\|.$ We also show empirical distributions of various quantities (like the components of the state vector) estimated over the long orbits directly in Figure 1 and Table 1. Thank you for this suggestion!

Thank you for your response! Your response regarding the Sobolev norm was valuable, and your further clarifications cleared up my concerns. I appreciate your discussion regarding the related works, which helped me better assess the status of your work. I will increase my score to accept.