When are dynamical systems learned from time series data statistically accurate?
Authors: Jeongjin Park
Nicole Tianjiao Yang
Nisha Chandramoorthy
Venue: NeurIPS 2024 poster (2024)
Keywords: generalization
ergodic theory and dynamical systems
chaotic dynamics
scientific machine learning
Primary area: machine_learning_for_physical_sciences
TL;DR: We combine ergodic theory with generalization to advance the theory and practice of learning chaotic systems using dynamical data.

Abstract: Conventional notions of generalization often fail to describe the ability of learned models to capture meaningful information from dynamical data. A neural network that learns complex dynamics with a small test error may still fail to reproduce its \emph{physical} behavior, including associated statistical moments and Lyapunov exponents. To address this gap, we propose an ergodic theoretic approach to generalization of complex dynamical models learned from time series data. Our main contribution is to define and analyze generalization of a broad suite of neural representations of classes of ergodic systems, including chaotic systems, in a way that captures emulating underlying invariant, physical measures. Our results provide theoretical justification for why regression methods for generators of dynamical systems (Neural ODEs) fail to generalize, and why their statistical accuracy improves upon adding Jacobian information during training. We verify our results on a number of ergodic chaotic systems and neural network parameterizations, including MLPs, ResNets, Fourier Neural layers, and RNNs.

When are dynamical systems learned from time series data statistically accurate?
Jeongjin (Jayjay) Park
Nicole Tianjiao Yang Department of Mathematics Emory University Atlanta, GA 30322 tianjiao, yang@emoryy edu
School of Computational Science and Engineering Georgia Institute of Technology
Atlanta, GA 30332 14ttta1111@aatchhhhcc edu
Nisha Chandramoorthy"
Department of Statisties The University of Chicago Chicago, IL 60637
rishac@uchicago.. edu
Abstract
Conventional notions of generalization often fail describe the ability of learned models to capture meaningful information from dynamical data. A neural network that learns complex dynamics with a small test eтror may still fail reproduce its physical behavior, including associated statistical moments and Lyapunov ex- ponents. To address this gap. we propose an ergodic theoretic approach to gener- alization of complex dynamical models learned from time series data. Our main contribution define and analyze generalization of broad suite of neural repre- sentations of classes of ergodic systems, including chaotic systems, in way that captures emulating underlying invariant, physical measures. Our results provide theoretical justification for why regression methods for generators of dynamical systems (Neural ODEs) fail to generalize, and why their statistical accuracy im- proves upon adding Jacobian information during training. We verify our results on a number of ergodic chaotic systems and neural network parameterizations, including MLPs, ResNets, Fourier Neural layers, and RNNs.
Introduction
Learning a dynamical system from time series data is pervasive challenge across scientific domains. Such data come from expensive experiments and high-fidelity numerical models that simulate the underlying nonlinear, often chaotic, processes. The learning challenge is to train on available data to produce output models that i) provably preserve known symmetries and invariances (e.g., conservation principles); and ii) are inexpensive surrogates for use downstream computations such optimization and uncertainty quantification. The field of physics-guided guided machine learning [vdGSB + 20 FO22 LK22 KKLL21| RPK19 KKL has emerged in response, rapidly integrating neural networks into data-driven modeling and prediction workflows for wide variety of complex dynamics, from geophysical fluid flows to phase transitions in materials (see [KMA 22 CCC+ for surveys). Yet rigorous generalization analyses of neural parameterizations in these applications, wherein the underlying dynamics can exhibit chaotic behavior, have been underexplored.
"Corresponding author
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

TRUE
MSE MLP
JAC MLP
Empirical Density
2000
40
40
True
0.050
MSE
1000
20
20
0.025
JAC
20
X 20 250
500
750
20
20 0.000
20
Figure 1: A random orbit the x-z plane obtained from RK4 integration of the Lorenz vector field ([oor63 (first column), Neural ODE, 'MSE_MLL trained with mean-squared loss (second column) trained with Jacobian loss (third column). The last column shows the empirical PDF of the orbit generated by the true (gray): IMS_MLP'' (red) and 'JAC_MLP' (blue) models. Experimental settings and additional results are in Appendices B and andgeeeettttivef Gist: A model trained well with MSE can produce atypical orbits but when Jacobian information is added to the training, it reproduces the long-term/sttaiitiial behavior accurately.
and Neural ODE,
Here we investigate data-driven neural parameterizations of chaotic ODEs/PDEs and maps (discrete- time dynamical systems) motivated by the following observation: neural representation that is learned well, i.e., with a small generalization error, can still produce unphysical long-term or ensemble behavior. Figure [] shows the classical Lorenz '63 chaotic attractor (Lor63] plotted using long random trajectory/orbi (time integration of the Lorenz equation vector field starting from a random initial condition, which is indicated by '+" sign). The second column shows an orbit from neural network model, 'MSE MLP', which minimizes the mean-sep -rranrd error in the represented vector field at 10,000 training points and shows high accuracy (< 5% average relative error) over 8000 test points on the attractor. Surprisingly, Figure [column 2) shows that the learned 'MSE_ 'LPp' Neural ODE [CRBDT8] model produces an atypical orbit an orbit different from almost every orbit of the true system for the same randomly chosen initial condition. As result, the learned empirical distribution not close to the physical distribution that of almost every true orbit, as shown in Figure I (column 4). On the other hand, the 'JAC_ __LP' model, which is obtained by minimizing the mean--quared error in the vector field and its first derivative (Jacobian matrix), reproduces the Lorenz also
'63 attractor and the physical distribution on the attractor (Figure I column 4). The captures all the Lyapunov exponents (LEs) - measures of asymptotic stability to perturbations, which invariants in ergodic systems accurately, while the 'MSE MLP' only obtains the leading LE accurately.
Naturally, we ask about the wider applicability of these observations. For any ground truth dynamical system, how can we quantify the probability success, including obtaining sample complexity results, of learning its physical or typical behavior? That is, how do we redefine and extend generalization analyses to neural representations of complex dynamical systems? To answer these questions, we start with a deceivingly simple supervised learning setup: given samples from time series, ((r,, 1t, S / S (m - 1), how can we learn a model, Fan, such that, i) 144 R Fnn(2,) for all time L, and ii) the underlying distribution of the states 2t, and other dynamical invariants such as Lyapunov exponents are reproduced by orbits of Fnn? ? It is widely accepted that learning such an Fan involves matching orbits of Fmn with x; over large t during training. However, small errors propagate over orbits, by definition, in a chaotic system, leading to training instabilitics. In response, vast literature has been dedicated developing sophisticated training models based RNNs
[PLH 17
PWF+ 18 RM21), operator learning [LLSK*2 *22 and regularizations LG22 FJNO20
We focus instead on the empirical risk minimization (ERM) for Fm that does not explicitly use the temporal correhationsddnammiaall structure in the data, avoiding training instabilities. Thus, we attempt to characterize when an elementary regression problem F. can still lead to learning a physical representation, leading to a practical theory of learning chaotic systems from data. Our specific contributions are as follows:
Motivated by extensive empirical results, we develop useful notions of generalization that characterize model's ability to reproduce dynamical invariants.
We develop new dynamics-aware generalization bounds for minimization of errors in the C", r = 0, 1, topology.
We leverage shadowing theory from dynamical systems to rigorously characterize failure modes in learning statistically accurate models.

---

2 Generalization of parameterized ergodic dynamics
In this section, we motivate, through illustrative examples, the need for a dynamics-awaan definition of generalization in the context of learning from time series data. We introduce concepts from dynamical systems as needed for a self-contained presentation.
Dynamics: A map, or discrete-timg dynamical system, F, is a function on a closed and bounded (compact) set, M C RR, which is called the state/phase space of the dynamics. We exclude scenarios where the dynamics can be unbounded, and focus on settings where F E C* (M) is differentiable function on M. We denote by Ft, F Z+, the iterates of the dynamics, or the compositions of F with itself times, i.e., F = F F An orbit or trajectory of Ft starting at an initial state x E M is the sequence {F* (x)) In practice, F' may be a numerical ODE solver that approximates the continuous-time solutions, s², E IR_ of the ODE (written in dynamical systems notation ): (x)/dtt = "(pp"(())))) (x)), where u : M + TM is the true vector field describing the governing equations. Fixing some E R. F = "p" We say map F is chaotic if there exists subbundle of TM, called the unstable subbundle, where infinitesimal perturbations grow exponentially under the linearization (Jacobian map), dFF,, of the dynamics. The true map F generates deterministic, autonomous system (see Appendix (A].
Physical measures. A probability measure J M RR+ is physical measure You02 for the dynamies if it a) F-invariant, b) ergodic and c) observable through F. A measure l is observable if time-averages along any orbit starting from randomly chosen initial point converge to constants that are expectations (phase space average) with respect to μ.. That is, for any f f C(M),
(L/T) Eigr ((F(x)) T-x, (z),,, for any initial point x chosen Lebesgue a.c. on set U S M. The orbits starting almost everywhere on the basin of attraction, U, asymptotically enter set, A, called the attractor. In dissipative chaotic systems, the attractor, A, which is the compact support of the physical measure μ. has Lebesgue measure O. Consequently, μ may not have probability density, that is, is not absolutely continuous, or is singular, with respect to Lebesgue measure.
Neural ODE. Introduced in [CRBDI8] Neural ODE, denoted here by, uo - Rd, with parameters, 0, is a vector field represented by a neural network. The vector field can be time integrated to obtain solutions v/b Rd + Rd I the ODE, (x))
(x)) As noted in the introduction, suppose we have n distinct pairs S = {(z;. (xx:)) which could come from a single orbit of F", as our training data. We train the Neural ODE by solving an ERM for the loss, f, that can be chosen to be a square loss, e.g., E(x, s0)) = (x) - F(x)]| That is. we solve for 0 that minimizes the training loss, (/n))zes a(x, 0). Note that the true ODE or vector field is not explicitly used in training, only the solution map at some time intervals, F. We refer to the map, x + TT (z), as Fun, or neural representation of the target map. F.
Data and optimization. Suppose we use an m-length orbit as our training data, i.c., z4+1 = Fxx), then, the data are not, strietly speaking, independent. However, feature of chaotic systems is an exponential decay of correlations, and so training data of the form, {(x;, F (xx)))} starting from an initial condition Ty ~ μt, can be thought of as iid, for large enough w N. In that case, the learned map Fun represents the function FW. During optimization, infinitesimal linear perturbations, will have to be evolved for time w. Since adjoint solutions blow up exponentially, a longer w will lead to training instabilities. To avoid numerical difficulties in training and focus on issues surrounding generalization, choose T 81, small time step, define the target map F and learn neural network representation of this function. When viewed this way, the generalization of Neural ODEs can be analyzed through the conventional lens of supervised learning with the loss
{(z, Fua) = (z) - F(x))1
where Fmn : ht is neural network representing the map. For map h M + M, we define the training and generalization errors in the usual way:
m
(s) = (1/m) Σ ((,,,)).
R(h)= = El(x,,
"We use this notation rather than the customary da(t)/dt u(x(t)). x00) = za, since in dynamical systems, we are interested in all possible orbits, as opposed to a particular orbit/path, for which detining the flow, y², is necessary.
3

---

1.0
0.9
true
0.8
LE3
mse
true
-10
jac
0.7
mse
jac
0.6
15
5000 sample
10000
5000 sample
10000
Figure 2: Learned and true LEs computed over 30,000 time steps using the QR algorithm of Ginelli et al [GCLPI3] starting from 10,000 random initial states. The "MSE" and "JAC" labels indicate computations using the Neural ODE models trained with the loss functions in T and respectively.
where the expectation is over the data distribution, which may be μ or any initial probability distribu- tion of the states.
Example. To illustrate our conceptual findings surrounding generalization, we use the canonical
Lorenz '63 system, which is a 3-variable reduced order model of atmospheric convection [Lor63].
The ODE can be written as (x) dtt = v(yp* (x)), where the vector field v is given by v(x)
z) - y, y,xy 3z] and x = [x. z] are the coordinate functions of the state æ. We use the standard values of the parameters o 10, p = 28, 3 = 8/3, at which the solutions are chaotic. We use the Runge- Kunge-Kutta 4-stage time integrator with a time step size of 0.01 to define the map That is, F(x) is the solution 00.0 (x) approximated by an RK4 time-integrator. Our Neural ODE map. Fun, is learned to approximate F by solving the above optimization with 10, 000 training points along an orbit. We illustrate our numerical results on various Neural ODE models and architectures that have fully connected layers, ResNet blocks with convolutional layers, and Fourier neural operators [KK) Many such models learn accurate representations of the true vector field v, as evidenced by small training and test errors (sample average approximation of the generalization error in 2 over 8,000 points from the data distribution). These are shown in Figure (Appendix C). while other hyperparameter and optimization settings are in Appendix B
Statistical measures and Lyapunov exponents. Our numerical results test the accuracy of models beyond generalization error as defined in (2). In particular, we compute time-averages using the Neural ODE models and compare against expectations with respect to f obtained from the true equations. For the Lorenz system, we note that time averages obtained from the models with small generalization errors can be inaccurate. That is, even the errors in the vector field are small (see also Figure 4 in Appendix C the time-averages can match poorly. This is described for the best performing Neural ODE model in Table 3 The discrepancy in the learned distribution is shown in terms of Wasserstein distance computed using empirical distributions on long orbits (of length 50,000). We find that the Neural ODE model does not learn the ground truth statistics even the training data include transient dynamics off the attractor.
Lyapunov exponents, roughly speaking. measure the asymptotic exponential growth/decay of in- finitesimal perturbations under the Jacobian map dF. In an ergodic system, they are independent of the initial state and can be written as an expectation with respect μ.. In this work, Lyapunov exponents are yet another dynamical invariant (statistical quantity) in ergodic systems that we use to evaluate the statistical fidelity of learned models. In a chaotic system, there is at least one positive Lyapunov exponent, and the number of positive Lyapunov exponents is the dimension of the unstable manifold. The Lorenz system has one positive Lyapunov exponent (LE), one zero LE (corresponding to center direction, the vector field itself) and one negative LE. The ground truth map F, which is a time öt approximation of the flow op has a two-dimensional center-unstable manifold and a one-dimensional stable manifold.
In Figure 2 we plot the Lyapunov exponents obtained using a classical QR iteration- based algorithm [CLPP3]. The Neural ODE model, marked 'MSE', produces a reasonably close approximation of the ground truth value (<x 0.9) for the positive LE but obtains an incorrect approximation of the stable LE (true value x 14.5). We note also that standard deviations in the LE values computed over different orbits quite large, indicating some orbits with atypical behaviors. We show the & error in the computed LEs by the Neural ODE model in Table and full details in Table

---

Table True vs. learned Lorenz system: comparison of statistics. (W Wasserstein-1 Distance, A:
set of LEs, for: empirical distribution of an orbit T. The subscript NN indicates quantities computed
using NN models trained with different loss functions (MSE (T. JAC B
Norm Difference
Model
Loss
W (is00:: (NN,500)
|| ANN
I| (r))00
(x) s00.NN
MLP
MSE
18.9711
9,6950
15.2220
MLP
JAC
0.6800
0.0118
0,6524
ResNet
MSE
1.3567
10.8516
0.7760
Res/ee
JAC
0.1433
0.0106
0,0559
FNO
MSE
10.5409
22.1600
9,4270
FNO
JAC
1.3076
0.0505
0,9748
Jacobian-matching. We now consider Neural ODE models trained with the following loss function.
Fun) [x) F Al|dF an (z) d=F()))1?
(3)
where Fun: .M M is a Neural ODE map, d' (z) T.M - T its (d x d) Jacobian matrix at I, and > is a hyperparameter that scales the relative importance of the two terms. We train in
the usual way, by solving an ERM problem to minimize the training error, Rs,,, which is defined as before, but with the new Jacobian-matching loss:
(h) = (1/m) Σ () Rh)h) = (x
(4)
We train Neural ODE models with similar architectures as above, and find best performing models in terms of the test error (approximation of the generalization error) in (4). As in the training with the loss function f, we find several accurate neural representations of the Lorenz map (see training and test loss plots in Appendix C). Their performance on statistical measures and Lyapunov exponent predictions is remarkably different, however. In Figure [ Neural ODE trained with the Jacobian- matching loss produces an attractor (column 3) that visually similar to the Lorenz attractor. The empirical distributions match closely with the ground truth, shown in column 5 of Figure Tand Table 3 where the model marked with a 'JAC', while one trained with the loss & in ① is indicated with an 'MSE'. The LEs and statistical averages of the state all accurately represented, although the model is not explicitly designed to learn temporal patterns the data. The MSE models learn vector fields with comparable accuracy with the Jacobian lacobian-matching models. That is, they learn accurate representations of F with high probability (over the data distribution). but failing to learn dF accurately leads to statistical inaccuracy. Due to the presence of atypical orbits, the attractor and the physical measure are not reproduced. Even though our formulation of Neural ODEs simply stated as an ERM for regression, the generalization errors R R(Fmn) and (Fnn) cannot determine whether Fm is a physical representation of the dynamics.
When generalization implies statistical accuracy
In the previous section, we observed that adding information about the Jacobian in the training process led to statistically-accure
learning. Does this observation apply more broadly to other ergodic systems? How can one further improve generalization? In this section, we provide answers to these questions by proving iyyaamics-awar generalization bounds for Neural ODEs.
Our ultimate goal is to minimize statistical loss function rather than 1 (3) since this mea- sures how accurately a model can reproduce ergodic averages associated with F. For instance, gstat(r, h) = supseLipp [=[[(z)] - (o /T) /T) f'hh (x))) This is the Wasserstein (W") distance (expressed in its dual form) between the ergodic measure and the ergodic measure associ- ated with the orbit of h. learned model, starting at x, limy-so Unif {2, h(x). (1)xz). ... hh" (x)). As noted, ERMs for this loss function do not have straightforward implementation when the map sought is chaotic. Hence, we seek conditions under which solving ERMs defined with losses (I) and can still minimize gstat, To understand when solving regression problems for F.n lead to statistically accurate physical representations, we first assume that notion of "shadowing" applies to F.
5

---

We say that the shadowing property applies to a map F if any 8 > 0, there exists an e <= = so that for any map G with F|I := suPReM ((=))
Definition 1.
d=())))) - there exists map M + M close to the identity such that ||G (T(z)) F - 5 for all L.
Intuitively, the shadowing property means that an orbit of a nearby dynamical system can closely follow a true orbit called a shadowing orbit- - of F for all time. This kind of miform-in--ime shadowing is a classical result for mathematically ideal class of chaotic systems, called uniformly hyperbolic systems (see Katok and Hasselblatt [KKH95] Ch 18; Appendix ().. For a texthook presentation and extension of shadowing for dynamical systems with some hyperbolicity, see [Pi106]. For uniformly hyperbolic/ F' (see Appendix A), we now assume that neural representation F'n of F trained with n samples generalizes well in terms of C'-distance, That is, an ERM solution for the loss B ('JAC' models in sections and 2 generalizes so that Rx (Fun) is small.
Definition 2 (C* generalization). Given 8 > O, there exist Eo, E1 > O and a function (ō, Eg. E1) - r(ō, Eo, & E N such that Ez~n [F(z)) - Fon (x)] S Eo and ExzulldF((r - ddF nn (x) S for all m 2 T with probability 2 1 - 5 over the randomness of the training data datafrom jC"
We now make an optimistic assumption on learned model, Fun, that satisfies the above definition of Ca generalization. Using Hoeffding's 's inequality, we know that for any óp > 0, with proba- bility at least 1 du over the randomness of x, [xx) - Fna(()) < & + (supxeM [1)(x) Fin (x)(l))) ( (266) and [xdF(x) - dff (x)] S E + (suPZEM [IdF(x) - dm (z)]l) /lol///o//////////// Given 5 > 0, let co = 28, and 61 = 28,. Fixing óg > O, suppose that the trained model Fun is F=)()))) < 12/22/logg2 (6)) and (sup as lIdF(z) dFFun df <
such that (sup as [1=(x)
c //(2//lo (2/5)) Taking a union bound, with probability > 1- (6 + 50), [[ - (1) S co and (ldF(x) - dd' (x)] S Ei. We enhance this inequality to obtain a stronger assumption on Fmn.
Assumption 1 (C* strong generalization), Given 8 >> 0, there exist co. 61 > O and a function (Co, C1, m)) - n((a, C,, (, N such that [[F (1)) F++ (x)]l I co and ||dF (F)n(:)) (x)) dF (F (x))l S Ei for all t n,, m 2 T(Ō,)), and n + oo as m + oo, with probability 2 1 1- d over the initial state x.
That is, we assume that, with high probability, the trained model makes a small error at each time. This stronger notion of generalization can be satisfied when the true model shows a smooth linear response in its statistics [Bal14], or in practice, training is performed with points sampled at random near the attractor, as opposed to with spin-off time to achieve a state on the attractor. Given tuple. (Co. c1), an orbit with initial condition x that satisfies, [FF(i (z)) - Fitt (z)]l S co and [IdF(Ftt un (z)) - d'Fun (i (=)) S é1 for all f S m will be referred to as an (eg. é1) orbit. That is, at each time, the neural representation Fun is (Eo, E1 -close to the true map, F. Under this assumption,
we can follow the proof of the Shadowing lemma (see e.g., Ch 18 of KKH95) for hyperbolic maps to show that a true orbit (of F) shadows every (co, 9) orbit.
Proposition (Shadowing). Let Fun be an approximation of F that satisfies the C strong gen- eralization (Assumption] 7 Given any ō > O, there exist eg; 61,n such that every (Eo.E1) orbit is 5-shadowed by an orbit of F. That is, there is a true orbit, say, {F((x), corresponding to every orbit, {G'(z') such that (()) (z')ll S 8, for all t S n.
Jesn be an n-length orbit of Fun: i.e., xxp = Fnu (=?" We use mm. MM to denote the direct sum @ Tati M of tangent spaces along the orbit, x:" The proof follows Theorem 18.1.3 of KKH95 to apply contraction mapping on a compact ball in M.
See section A.1 for the proof. Let {=:"
The above result defines, for each (ea; 61 )orbit, (^:2") }r, shadowing orbit, 29th == (F(xg" + vt) where o = # is the fixed point of the contraction map in the proof (section A.I. Let pe (I be the empirical measure defined on the shadowing orbit corresponding to an (eo. e Jorbit, {27m Je, i.c. 44th (=g²) = Unif (zgh but, unexpectedly, it may be atypical for the physical measure, Jr. That is, for an atypical shadowing orbit, the time average, (1/n) Eig, [((fh)) does not converge to the expected value (x)) as
.rth.... 2 }.. A shadowing orbit is indeed an orbit of the true map F,
co. This means that the Wasserstein distance, W (4C) (=g"), 1)), does not converge to 0 as n + 0o.
Given an initial condition, -n of an (eo, e1 -orbit, the corresponding shadowing orbit may be typical with some probability (over the distribution of aQ and this probability of finding typical shadowing orbits is property of the true dynamics, F. When μ. the physical measure of F, is highly

---

sensitive to perturbations of F, (see [Rue09, [all4] for surveys on linear response theory, the study of perturbations of statistics) this probability may be small. Although the connection between the sensitivity of statistics and the atypicality of shadowing orbits is not completely known, this can justify the differences in the statistical accuracy of neural parameterizations with good C' generalization 2 A neural model F... which generalizes well in the C* sense (Definition 2 can be thought of as a C² c'-smooth perturbation of the true dynamics F. If smooth perturbations of F can cause large change in μ4, then even when the training size
(see Definition]
oc, and the orbit length + oo, (1 (=)mm, ), 1t), may not converge to zero, resulting in a model that does not preserve the physical behavior of F. On the other hand, for typical shadowing, this possibility is excluded and gives characterization of statistically accurate learning.
Theorem 1 (Statistically accurate learning). Let Fun be a model of F that satisfies C' strong generalization. In addition, let Fun and F be such thatfor any 6 >> 0, there exists an eg > 0 so that (yel lim (I); p) S 62 with probability (over the randomness of x) 2 1 1- 6. Then, for any 5 > such that 1im W (Unifffa, (I)(.).... Fu (z), Fin (z)) p1) S
0, there exists
with probability 2 1 - 5.
Proof: Let p. (z) == (z, F'n .... F'm (z)... ... . (x)} be the empirical measure of an n-length orbit of Fmn starting at Given 8 > O, we choose (eo, t1) (co,E1) so that C strong generalization is satisfied with probability 2 1- 6/2. Thus, an initial condition x is such that (x),(=))ss is an (Co, G1) orbit, where the tuple (co.ci is defined in Proposition with prob- ability 2 8/2. Applying Proposition] [ we have, for any 1-Lipschitz function J: M + R, (L/n) Eign (fin) fFF S (I/n) [nnE een (x)) - F* F'(x))l S 5/2. Taking a supre- mum over f. W (x), (z), (="())) (z)) I 8/2, with probability 2 1- - 6/2. By assumption, there exists some C2 > 0 such that lim,- W (n. (z)) < c2 with probability
8/2. Hence, lim W (n. pi (z)) S lim W (nch (z); p. (x)) - W (u. (=))) S 62 + 8/2 == e, with
probability > 1 - 8, using triangle inequality and taking union bound.
This result explains why training to minimize Jacobian-matching loss B can lead to statistically accurate models, even though, long-time time temporal patterns in the data are not learned explicitly by
regression for the one-time map F. Since C' generalization (i.e., small errors )n is insufficient for learning shadowing orbits, and thus for Proposition and Theorem to hold, the models trained on MSE loss are not expected to learn ergodic/statistical averages with respect j..
When shadowing orbits are atypical with high probability, we observe numerically that C generaliza- tion, i.c., training with Jacobian-matching matching loss, still does not produce statistically accurate dynamics, in line with the above theorem. For instance, for maps with atypical shadowing described in (CW21), we find that Iearned neural representations with
do have good O'l generalization (Figure 6). but do not exhibit good statistical accuracy and learn incorrect Lyapunov exponents (see Table 5 Plucked Tent map).
4 Dynamic generative models
So far, we have focused on understanding the statistical accuracy of supervised learning of dynamical systems. Without any minimization of distances on the space of probability measures, we proved sufficient conditions under which regression with aecobian-matchiin from with high probability. A generative method is an unsupervised learning technique to train on samples from a target distribution to produce more samples (provably) from the target. Naturally, can several popular generative models for our target physical measure here, including score- based methods SE11 SSDK 20), Variational Autoencoders (VAE) [RM15) or normalizing flows [RM13] PNR 21). However, these methods neglect the dynamical relationships in the input samples. In other words, from a vanilla generative model of a physical measure, we cannot also recover the true dynamics. Thus, we focus here on Latent SDEs models from (LWCD20) which combine neural representations of dynamics with generative models. We reinterpret them as dynamic generative models, and analyze their ability faithfully represent both the underlying dynamics as well as the physical measure.
information can yield samples
In a dynamic generative model, we approximate Ft with stochastic map, that can written as, Fis := fe o t go, where the subscript Is stands for "latent SDE". Here, the function Gø : Rd > R4, with Icarnable parameters, o, is a (possibly stochastic) embedding from the data to latent
7

---

space (RR4), ) such that, Pörpi = 96.0. Recall the pushforward notation (f), i... if x μ., then, go(x) ~ qa.0. The dynamical system & Rd 7 Ed acts on the latent space, and defines a sequence of pushforward distributions 4 20 = qs.t. A special case of this setup is "latent ODE", where & is deterministic map instead. With a stochastic latent SDE model instead, we observe, in line with [LWCD20), that the multimodal distribution of the Lorenz '63 attractor is reproduced better. That is, D is a solution map of a Neural SDE: d (2) = (t, & (z))dt + a,,t, (x)))) o dW,- Here the drift term, Wa. and the diffusion term, Gø. are represented as neural networks. The decoder fo IRd: + gd is a deterministic map that defines the conditional, forqs,t = pa (-)z))
This dynamic VAE approach has been found to be expressive for chaotic systems (see Chapters 3 and 5 of Kid22], (KMFL20), wherein the conditional distribution, x + [e,0(2 (=)m)) of Zo given X1-m == {F" (x)()))e (a) is modeled as Gaussian distribution whose learnable parameters are also denoted by 6. Similarly. the conditional (- Ze) is again modeled as a Gaussian with parameters 8. The parameters, & and @ respectively, of the encoder and the decoder are trained by maximizing the following dynamic version of the evidence lower bound (ELBO): l (X Xism; (Ø, 0)) == Di=1 E Pe (-IX | (x<e |+)] + (- 1.mm llpzo)) )], where the prior pz, follows standard Gaussian distribution in the latent dimension. For alternatives to the ELBO objective above, such the Wasserstein- GAN objective, we refer the reader to [Kid22].
We now evaluate both the learned dynamics, Fis, and the learned generative model, pe, which approximates μ. In Figure [T] we present the empirical distribution of generated orbit against that of a typical orbit of the Lorenz system (μ)). We observe that the distributions match well, nevertheless the vector field is not well--apprximmated First, even though the system is deterministic, the learned Ø with minimum generalization error (EEetta) encountered in the hyperparameter search (Appendix C.8) is not, i.e, the diffusion term 0g is not zero. The learned stochastic map, Fis, produces an incorrect stable LE for the Lorenz system (== 11.8), while the leading unstable LE matches reasonably well.
Since the map Fis, or its underlying vector field, on the latent space is not unique, we may obtain maps that do not preserve dynamical structure invariants, even if the generated samples from Po approximately capture As noted in section [ the physical measure ls is often singular, but absolutely continuous on lower-dimensional manifolds, leading to lack of theoretical guarantees for vanilla generative models (Pid22). Finally, the sample complexity (and tight generalization bounds) of generative models, the above variational optimization, are not fully understood theoretically, especially for singular distributions (that satisfy the manifold hypothesis (BCVT3)). We remark that since the minimax rates for approximating distributions have an exponential dependence on the dimension, exploiting the intrinsic dimension (unstable dimension) associated with the support of will be key to tractable generative models for μ in sigh--imensional chaotic systems. Even in the Lorenz '63 system, we require O119)) samples for training reasonably accurate model in Figure I1 while the macchinn (Figure 5) produces smaller Wasserstein distances with fewer samples (10*), and simpler regression problem as opposed to variational optimization above.
5 Numerical Experiments
We conduct experiments with the MSE and JAC losses in (1) and (3) respectively on many canonical chaotic systems: 1D tent maps, 2D Bakers map, 3D Lorenz '63 system and the Kuramoto Sivanshinsky equation (127 dimensional system after discretization) (Appendix C In each system, we identify the Neural ODE model with the lowest generalization errors R and RA in 2 and (4) respectively. by an architecture and optimization hyperparameter search (Appendix B On these models that generalize well, we perform tests of statistical accuracy and LE computations (as described in section 2 for the Lorenz '63 system). Consistent with Theorem I although most of the considered systems are not uniformly hyperbolic, we find that the JAC models are statistically accurate and reproduce the LEs, while the MSE models with R comparable to Rx are not statistically accurate. For the KS system for instance, the MSE models even overpredict the number of positive LEs. Interestingly, the best JAC model can learn more than half of the first 64 LEs, compared to the best MSE model that can learn only 2 out of 64 LEs, with < 10% relative error. The Python code is available at https:/ htte:://giiii http//ggtttt...////// com/ sha- c/ ctsaaCNDDE
8

---

Related work
Neural ODEs and generalization. Introduced as continuous- analogues Weil7
HR17
CRBDI8] of Residual neural networks (ResNets) [HZRS16], Neural ODEs [CRBDI8 HR17 offer a vector field parameterization that can be time-integrated with an ODE solver. Augmented Neural ODEs DDT19 RCDI9 DA22] demonstrate improved expressivity for complex dynamics, and some universal approximation results appear in (ZGUA20) Kid22]. Neural ODEs [CRBD18] as normalizing flows [RMI5] and for density estimation have been tackled in (GCB+19) Training and regularization techniques that allow handling long time series are the subject of numerous works [RMM 20 KJD 21 EJNO20 PER23a PER23b PMSR21
MSKF2T. Our focus is dif- ferent: we characterize when an elementary miccbian-matching regularization serves as an inductive bias toward learning physical representations from irregular, even chaotic data.
Data-driven surrogates of complex systems. The vast and growing literature in the field of physics- informed machine learning [AAK22 LKB22 21 HZ |vdGSB+20 W Y21 RPK19
WY21 KKLL21
LK22] has encouraged the adoption of physical machine learning models that preserve physical properties, symmetries and conservation laws, and are yet applicable in scientific problems [HKU220] 21 HZB 22 RBL22 VAUK22 vdGSB 20
23] At the same time, several purely data-driven methods for complex systems [CNH20 PSH 22),, which do not require an expensive high-fidelity solver, have gained attention for their impressive prediction skill [RCVS+ SLST17
LSGW 23] and both faster training and inference making them suitable
for optimization BB21 RBL22 LP21 ID21 and inverse problems [HV223 ABOV23]. Since the fundamental regression problem study underlies both hybrid and data-driven methods, provide insight into dynamics-awand generalization (e.g.. reproducing ergodic behavior, Lyapunov exponents etc) applicable to different surrogate modeling approaches. Several innovative approaches for ensuring training stability [MMD22 SWP+24 [MBD23 [LOW24)) in chaotic systems when using recurrent architectures have been proposed recently. The failure of generalization notions based on mean -eaan-suared eттor have also been noted in (SWP + 24 [LOW24] and empirical strategies and new definitions of generalization suitable of non-ergodic systems have been introduced in [GHB 24). For ease of theoretical analysis, we do not consider these more sophisticated training approaches, choosing instead learn short-term dynamics which obviates the need for stabilization strategies. Encouragingly, we observe low sample complexity of vanilla regression with Jacobian information to learn physical measures, when compared to the generative modeling approaches (which can be comparable to RNNs as well). A computational analysis of the Jacobian loss training compared to generative modeling/stablized recurrent training is deferred to future work.
Ergodic theory and shadowing. Hyperbolic dynamics and ergodic theory (see e.g. the texthook (KH955) lay the foundation for understanding the -in-tiim/stattiticall physical behavior [You022] of complex systems. The scientific computing community has leveraged ergodic theory and shadowing for rigorous computations that use high-fidelity numerical simulations of
Ano67 B5
Pi106)
chaotic systems GL24
Ni21] and to analyze the correctness of numerical simulations GHYS90). The novelty of our work lies in introducing shadowing as the basis for generalization, thus providing new analysis tools to understand the correctness of learned chaotic systems [LLSK
Wan
HYG87 CW21
Lial7
Sau05
22) An interesting direction for future work is to extend dynamics aware generalization bounds similar to Theorem 1 to operator learning with Sobolev norms introduced in ILLSK+22
7 Conclusion
Our dynamics-aware generalization (Theorem D and empirical results provide a new characterization of statistical accuracy in models learned from dynamical data. These results open many avenues for improving mechanistic understanding and bridging the theory-practics gap in physical neural modeling of complex systems:
Understanding learning attractors. By exposing foundational problems in the elementary and fairly general setting of regression of a dynamical system, our analytical tools in section 3 broaden our conceptual understanding of learning from time series data in more complicated models and paradigms.

---

Learning dynamical representations vs. generative models. We show that physical measures can be produced by accurate neural representations of the dynamics, which can be more data and time-efficient compared with generative modeling for time series generated by chaotic systems.
Dynamics-aware learning of scientific models. Our results imply that generalization does imply statistical accuracy and preservation of dynamical invariants and properties such as Lyapunov exponents, which crucial for trustworthy ML for science. Thus, reinforce the need to move toward context-awarc theory of generalization, organically unifying complex dynamics with learning theory.
Limitations: A limitation of the theoretical results about C' generalization is that we need to assume the typicality of shadowing. Empirically, the Jacobian be expensive to estimate for high-dimensional scientific applications. We leave the study of statistical accuracy for learning atypical shadowing orbits, and the extension of an efficient algorithm for Jacobian information during training as future work.
Acknowledgments We grateful the anonymous reviewers for their useful feedback and their suggestions of additional numerical results, which have strengthened the paper. NC acknowledges support from the James C. Edenfield faculty fellowship provided by the College of Computing Georgia Tech.
10

---

References
[AEOV23] Harbir Antil, Howard C Elman, Akwum Onwunta, and Deepanshu Verma. A deep neu- ral network approach for parameterized pdes and bayesian inverse problems. Machine Learning: Science and Technology, 4(3):035015, 2023.
[Ano67] D. V. Anosov. Geodesic flows on closed riemannian manifolds of negative curvature. Trudy Mat. Inst. Steklov., 90, 1967.
[Bal14] Viviane Baladi. Linear response, or else. arXiv preprint
2014.
[BB21] William Bradley and Fani Boukouvala. Two-stage approach to parameter estimation of differential equations using neural odes. Industrial Engineering Chemistry Research, 60(45):16330-1634, 2021.
[BCV13] Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new perspectives. IEEE iransactions on pattern analysis and machine intelligence, 35(8):1798-1822 2013.
[BJS24] Oscar F Bandtlow, Wolfram Just, and Julia Slipantschuk. A numerical study of rigidity of hyperbolic splittings in simple two-dimensional maps. Nonlinearity, 37(4):045007, mar 2024.
[Bow75] Rufus Bowen, w-limit sets for axiom a diffeomorphisms. Journal of differential equations, 18(2):333333, 1975.
[BW14] Patrick J Blonigan and Qiqi Wang. Least squares shadowing sensitivity analysis of a modified siissityss equation. Chaos, Solitons & Fractals, 64:16-25, 2014.
[CCC+ 19] Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt-Maranto, and Lenka Zdeborová. Machine learning and the physical sciences. Rev. Mod. Phys., 91:045002, Dec 2019.
[CNH20) Ashesh Chattopadhyay, Ebrahim Nabizadeh, and Pedram Hassanzadeh. Analog forecasting of extreme-causing weather patterns using deep learning. Journal of Advances in Modeling Earth Systems, ((2)=e201MMS00
2020.
[CRBD18] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances neural information processing systems, 31, 2018.
[CW21] Nisha Chandramoorthy and Qiqi Wang. On the probability of finding nonphysical solutions through shadowing. Journal of Computational Physics, 440: 10389, 2021.
[CW22] Nisha Chandramoorthy and Qiqi Wang. Efficient computation of linear response of chaotic attractors with dimennsinni unstable manifolds. SIAM Journal on Applied Dynamical Systems, 21(2):735 2022.
[CZH21] Michele Coti Zelati and Martin Hairer. A moise-induced transition in the lorenz system. Communications in Mathematical Physics, 383:2243-221 2274, 2021.
[DA22] Thai Duong and Nikolay Atanasov. Adaptive control of se (3) hamiltonian dynamics with learned disturbance features. IEEE Control Systems Letters, 6:2773-2778 2022.
[DT199] Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural odes. Ad- vances in neural information processing systems, 32, 2019.
[EHL04] Gregory I. Eyink, Tom WN Haine, and Daniel J Lea. Ruelle's linear response formula, ensemble adjoint schemes and lévy flights. Nonlinearity, 17(5):1867, 2004.
[FJNO20] Chris Finlay, Jörn-Henrik Jacobsen, Levon Nurbekyan, and Adam Oberman. How to train your neural ode: the world of jacobian and kinetic regularization. International conference machine learning, pages 3154-3164, PMLR, 2020.

---

[FO22] Alexander G Fletcher and James M Osborne. Seven challenges in the multiscale modeling multicellular tissues. WIREs mechanisms of disease, e4(1) 2022.
[GBD+20] Arnab Ghosh, Harkirat Behl, Emilien Dupont, Philip Torr, and Vinay Namboodiri. Steer: Simple temporal regularization for neural ode. Advances in Neural Information Processing Systems, 111114483848888 2020.
[GC95] Giovanni Gallavotti and Ezechiel Godert David Cohen. Dynamical ensembles in nonequilibrium statistical mechanics. Physical review letters, 74(14):2694, 1995.
[GCB+ 19] Will Grathwohl, Ricky TQ Chen, Jesse Bettencourt, Ilya Sutskever, and David Du- venaud. Ffjord: Free-form continuous dynamics for scalable reversible generative models. 7th International Conference on Learning Representations, ICIR 20/9, 2019.
[GCLP13] Francesco Ginelli, Hugues Chaté, Roberto Livi, and Antonio Politi. Covariant Iyapunov vectors. Journal of Physics A: Mathematical and Theoretical, 46(25):254005, 2013.
(GHB+24) Niclas Göring, Florian Hess, Manuel Brenner, Zahra Monfared, and Daniel Durstewitz. Out-of-domain generalization in dynamical systems reconstruction. arXiv preprint 1881: 2024.
[GHL20] Bernard J Geurts, Darryl D Holm, and Erwin Luesink. Lyapunov exponents of two stochastic lorenz 63 systems. Journal of statistical physics, 179:1343-1365. 2020.
[GHYS90] Celso Grebogi, Stephen M. Hammel, James A. Yorke, and Tim Sauer. Shadowing of physical trajectories in chaotic dynamics: Containment and refinement. Phys. Rev. Lett., 65:1527- Sep 1990.
[GL.2]] Rhys Gilbert and Davide Lasagna. On the application of optimal control techniques to the shadowing approach for time averaged sensitivity analysis of chaotic systems. SIAM Journal on Applied Dynamical Systems, 23(1):505-552 2024.
[KUT20]] Philipp Holl, Vladlen Koltun, Kiwon Um, and Nils Thuerey, phiflow: A differen- tiable pde solving framework for deep learning via physical simulations. In NeurlPS workshop, volume 2, 2020.
[HMBD23] Florian Hess, Zahra Monfared, Manuel Brenner, and Daniel Durstewitz. Generalized teacher forcing for learning chaotic dynamics. arXiv preprint arXiv:2306.04406. 2023.
[HR17] Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse problems, 34(1):01 :00404444 2017.
(HVT23) Benjamin Holzschuh, Simona Vegetti, and Nils Thuerey. Solving inverse physics problems with score matching. Advances in Neural Information Processing Systems, 36, 2023.
[HYG87] Stephen M Hammel, James A Yorke, and Celso Grebogi. Do numerical orbits of chaotic dynamical processes represent true orbits? Journal of Complexity, 3(2): 136- 145, 1987.
(HZB+21)
Changnian Han, Peng Zhang, Danny Bluestein, Guojing Cong, and Yuefan Deng. Artificial intelligence for accelerating time integrations multiscale modeling. Journal of *Coomutational Physics, 12051,, 2021.
[HZRS16] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770-778, 2016.
[JD2]] Elnaz Jamili and Vivek Dua. Parameter estimation of partial differential equations using artificial neural network. Computers & Chemical Engineering, 147:107221, 2021.
[JLOW24] Ruoxi Jiang, Peter Y La, Elena Orlova, and Rebecca Willett. Training neural operators to preserve invariant measures of chaotic attractors. Advances in Neural Information Processing Systems, 36, 2024.
12

---

[Kid22] Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.024: 2022.
[KJD*21] Suyong Kim, Weigi Ji, Sili Deng, Yingbo Ma, and Christopher Rackauckas. Stiff neural ordinary differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, (1))), 2021.
[KKH95] Anatole Katok, AB Katok, and Boris Hasselblatt. Introduction to the modern theory of dynamical systems. Number 54. Cambridge university press, 1995.
[KKL*21] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422- 440, 2021.
[KKLL21] Sung Wook Kim, Iljeok Kim, Jonghwan Lee, and Seungchul Lee. Knowledge integra- tion into deep learning in dynamical systems: an overview and taxonomy. Journal of Mechanical Science and Technology, 35:1331-1342, 2021.
KMA+21) Karthik Kashinath, M Mustafa, Adrian Albert, JL Wu, C Jiang, Soheil Esmaeilzadeh, Kamyar Azizzadenesheli, R Wang. Ashesh Chattopadhyayy A Singh, al. Physics- informed machine learning: case studies for weather and climate modelling. Philo- 2021.
sophical Transactions of the Royal Society A, 379(2199))202000))
[KMFL20] Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differ- ential equations for irregular time series. Advances Neural Information Processing Systems, 33:6696-6707 2020.
[LG22] Alee J. Linot and Michael D. Graham. Data-driven reduced-order modeling of spa- tiotemporal chaos with neural ordinary differential equations. Chaos: An Interdisci- plinary Journal of Nonlinear Science, 32(7):0731 10, 07 2022.
[LH17] Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. arXiv preprint arXiv: 2017.
[Lia17] Shi-jun Liao. On the clean numerical simulation (ens) of chaotic dynamic systems. Journal of Hydrodynamics, Ser. B, 29(5):729-747, 2017.
[LK22] Yi Heng Lim and Muhammad Firmansyah Kasim. Unifying physical systems' induc- tive biases in neural ode using dynamics constraints. arXiv preprint arXiv:2208.00 026632....2.,,, 2022.
[LKA+ 20] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhat- tacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for para- metric partial differential equations. arXiv preprint arXiv:2010.08895, 2020.
[LKB22] Yuying Liu, J. Nathan Kutz, and Steven L. Brunton. Hierarchical deep learning of multiscale differential equation time-steppers. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 380(2229):0211020, 2022.
[LLSK*22] Zongyi Li, Miguel Liu-Schiaffini, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learn- ing chaotic dynamics in dissipative systems. Advances Neural Information Process- ing Systems, 35:16768-16781 2022.
[Lor63] Edward N Lorenz. Deterministic nonperiodic flow. Journal of attmopphericc sciences, 20(2):130-141) 1963.
[LP21] Kookjin Lee and Eric Parish. Parameterized neural ordinary differential equations: Applications to computational physics problems. Proceedings of the Royal Society A, 477(0253):02200122, 2021.
13

---

[LSGW 23] Remi Lam, Alvaro Baachez-Gonzaeez.. Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Oriol Vinyals, Jacklynn Stott, Alexander Pritzel, Shakir Mohamed, and Peter Battaglia. Learning skillful medium-range global weather forecasting, Science, (82(6677):
2023.
[LWCD20] Xuechen Li, Ting-Kan Kam Leonard Wong, Ricky TQ Chen, and David Duvenaud. Scalable gradients for stochastic differential equations. In International Conference on Artificial Intelligence and Statistics, pages 3870-3882. PMLR, 2020.
[MMD22] Jonas Mikhaeil, Zahra Monfared, and Daniel Durstewitz. On the difficulty of learning chaotic dynamics with rnns. Advances in Neural Information Processing Systems, 35:1 1297-11312, 2022.
[MSKF21] James Morrill, Cristopher Salvi, Patrick Kidger, and James Foster. Neural rough differential equations for long time series. In International Conference on Machine Learning, pages 7829-7838. PMLR, 2021.
[Ni21] Angxiu Ni. Approximating linear response by nonintrusive shadowing algorithms. SIAM Journal Numerical Analysis, 59()):2843-2865, 2021.
[PA Grace CY Peng, Mark Alber, Adrian Buganza Tepole, William R Cannon, Suvranu De, Savador Dura-Bernal, Krishna Garikipati, George Karniadakis, William W Lytton, Paris Perdikaris, al. Multiscale modeling meets machine learning: What learn? Archives of Computational Methods in Engineering, 28: 1017-1037, 2021.
[PHR23a] Avik Pal, Alan Edelman, and Chris Rackauckas. Continuous deep equilibrium mod- els: Training neural odes faster by integrating them to infinity. In 2023 IEEE High Performance Extreme Computing Conference (HPEC), pages 1-9. IEEE, 2023.
[PER23b] Avik Pal, Alan Edelman, and Christopher Vincent Rackauckas. Locally regularized neural differential equations: some black boxes meant to remain closed! In International Conference on Machine Learning, pages 26809-26819 PMLR, 2023.
[Pid22] Jakiw Pidstrigach. Score-based generative models detect manifolds. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances Information Processing Systems, volume 35, pages 35852-35865. Curran Associates, Inc., 2022.
Advances in Neural
[Pi106] Sergei Yu Pilyugin. Shadowing near an invariant set. Shadowing in dynamical systems, pages 1-101, 2006.
[PLH 17] Jaideep Pathak, Zhixin Lu, Brian R. Hunt, Michelle Girvan, and Edward Ott. Using machine learning to replicate chaotic attractors and calculate Yуарuпо exponents from data. Chaos, 27(12), 2017. Cited by: 364; All Open Access, Hybrid Gold Open Access.
[PMSR21] Avik Pal, Yingbo Ma, Viral Shah, and Christopher V Rackauckas. Opening the blackbox: Accelerating neural differential equations by regularizing internal solver heuristics. In International Conference on Machine Learning, pages 8325-8335.
PMLR, 2021.
[PNR 22] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57): 2021.
[2SH+22] Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay. Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Arizzadenesheli, et al. Fourcastnet: global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint 2rXiv:2202.11214 2022.
14

---

[PWF* 18] Jaideep Pathak, Alexander Wikner, Rebeckah Fussell, Sarthak Chandra, Brian R. Hunt, Michelle Girvan, and Edward Ott. Hybrid forecasting of chaotic processes: Using machine learning in conjunction with knowledge-hased model. Chaos, 28(4), 2018. Cited by: 206; All Open Access, Green Open Access.
[RBL22] Johann Rudi, Julie Bessac, and Amanda Lenzi Parameter estimation with dense and convolutional neural networks applied to the fizhugh-nagumo ode. In Mathematical and Scientific Machine Learning, pages 781-808 PMLR, 2022.
[CCD19] Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differen- tial equations for eoegularly-samplee time series, Advances in neural information processing systems, 32, 2019.
[RCVS+ + 19] Markus Reichstein, Gustau Camps-Vall Valls, Bjom Stevens, Martin Jung. Joachim Denzler, Nuno Carvalhais, and fnm Prabhat. Deep learning and process understanding for data- driven earth system science. Nature, 566(7743):195-22 954204, 2019.
[RM15] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International conference on machine learning, pages 1530-1538. PMLR, 2015.
[RM21] Alberto Racca and Luca Magri. Robust optimization and validation of echo state networks for learning chaotic dynamics. Neural Networks, 142:252-268, 2021.
[RMM+20) 20] Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint X0X0:200000043555 2020.
[Rös76] Otto Rössler. An equation for continuous chaos. Physics Letters A, 57(5):397-398, 1976.
[RPK19] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 978:686-707,2011
[Rue09] David Ruelle. A review of linear response theory for general differentiable dynamical systems. Nonlinearity, 22(4))855, 2009.
[Sau05] Tim Sauer. Computer arithmetic and sensitivity of natural measure. Journal of Difference Equations and Applications, 11(7):669-676, 2005.
[SE19] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019.
[SLST17] Tapio Schneider, Shiwei Lan, Andrew Stuart, and João Teixeira. Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations. Geophysical Research Letters, 44(24):12-396 2017.
[SSDK+20) Yang Song, Jascha Sohl-Dicksted Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:201 1.,456, 2020.
[SWP+24] Yair Schiff, Zhong Yi Wan, Jeffrey B Parker, Stephan Hoyer, Volodymyr Kuleshov, Fei Sha, and Leonardo Zépeda-Núñez. Dyslim: Dynamics stable learning by invariant measure for chaotic systems. arXiv preprint urXiv:2402.0444 2024.
[VAUK22] Pantelis R Vlachas, Georgios Arampatzis, Caroline Uhler, and Petros Koumoutsakos. Multiscale simulations of complex systems by learning their effective dynamics. Nature Machine Intelligence, 4(4):359-366, 2022.
15

---

[vdGSB + 20] Erik van der Giessen, Peter A Schultz, Nicolas Bertin, Vasily V Bulatov, Wei Cai, Gábor Csányi, Stephen M Foiles, M G D Geers, Carlos González, Markus Hütter, Woo Kyun Kim, Dennis M Kochmann, Javier L.orraa Ann E. Mattsson, Jörg Rottler, Alexander Shluger, Ryan Sills, Ingo Steinbach, Alejandro Strachan, and Ellad B Tadmor. Roadmap on multiscale materials modeling. Modelling and Simulation in Materials Science and Engineering, 28(4):44300, mar 2020.
[Wan13] Qiqi Wang. Forward and adjoint sensitivity computation of chaotic dynamical systems. Journal of Computational Physics, 235:1-13, 2013.
[Wei17] Ee Weinan. A proposal on machine learning via dynamical systems. Communications in Mathematics and Statistics, 1(5):1-11, 2017.
[WG18] Caroline L Wormell and Georg A Gottwald. On the validity of linear response theory in sigh-diimessionl deterministic dynamical systems. Journal of Statistical Physics, 172:1479-1498, 2018.
[WY21] Rui Wang and Rose Yu. Physics-guided deep learning for dynamical systems: A survey. arXiv preprint arXiv:2107.01272 2021.
[3]+23] Sungduk Yu, Walter Hannah, Liran Peng, Jerry Lin, Mohamed Aziz Bhouri, Ritwik Gupta, Bjöm Lütjens, Justus C. Will, Gunnar Behrens, Julius Busecke, Nora Loose, Charles Stern, Tom Beucler, Bryce Harrop, Benjamin Hillman, Andrea Jenney, Sa- vannah L. Ferretti, Nana Liu, Animashree Anandkumar, Noah Brenowitz, Veronika Byring, Nicholas Geneva, Pierre Gentine, Stephan Mandt, Jaideep Pathak, Akshay Subramaniam, Carl Vondrick, Rose Yu, Laure Zanna, Tian Zheng, Ryan Abernathey, Fiaz Ahmed, David Bader, Pierre Baldi, Elizabeth Barnes, Christopher Bretherton, Peter Caldwell, Wayne Chuang, Yilun Han, YU HUANG. Fernando Iglesias-Suarez Sanket Jantre, Karthik Kashinath, Marat Khairoutdinov, Thorsten Kurth, Nicholas Lutsko, Po-Lun Ma. Griffin Mooers, J. David Neelin, David Randall, Sara Shamekh, Mark Taylor, Nathan Urban, Janni Yuval, Guang Zhang, and Mike Pritchard. Climsim: A large multi-scale dataset for hybrid physics-m ml climate emulation. In Oh, Nau- mann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages Inc., 2023.
Curran Associates,
[You02] Lai-Sang Young. What are srb measures, and which dynamical systems have them? Journal of statistical physics, 108:733-754, 2002.
[You17] Lai-Sang Young. Generalizations of srb measures to nonautonomouss random, and infinite dimensional systems. Journal of Statistical Physics, 166:494-515, 2017.
[ZGUA20] Han Zhang, Xi Gao, Jacob Unterman, and Tom Arodz. Approximation capabilities of neural odes and invertible residual networks. In International Conference on Machine Learning, pages 1186-11095. PMLR, 2020.
[Zha17] Lingmei Zhang. novel 4-d butterfly hyperchaotic system. Optik, 131:215-220, 2017.
16

---

A Proofs and assumptions
Our results in section 3 make the assumption that the map F is uniformly hyperbolic. Roughly speaking, this means the uniform (on the attractor A C M) expansion and contraction of infinitesimal perturbations under the differential d/F. More precisely, in a uniformly hyperbolic system, there exists a decomposition of the tangent bundle, TM- = F. ( E", into an unstable (F") and stable (E*) subbundle such that these subbundles are d'F-innarian and moreover, for any vector field e E* ||dF & C X ||o]l t E zz. Similarly, an infinitesimal perturbation along vector field e E E" shows a uniform exponential decay in norm backward in time, ||dF'u|| S C X ||ull, L E z- Here the constants C, A are uniform over A. There are several relaxations of the uniformity in the splitting to obtain nonuniform and partial hyperbolicity; there are also several known examples of non-hy non-hyperbolic chaotic systems. Assumptions of hyperbolicity are common in dynamical systems analyses due to well-understood ergodic theoretic results, including existence of physical measures of the SRB-type [You0 You Such assumptions are also widespread in computational dynamics, wherein algorithms derived rigorously for uniformly hyperbolic systems are found to be applicable [WG18] to turbulent fluids and chaotic systems in practice [GC95] EHL00 CW22, where rigorous verification [BJS24] of these assumptions is not feasible.
We remark that we only require high-probability finite-time shadowing, and thus, our results can possibly be extended under relaxations of uniform hyperbolicity. Despite being mathematically convenient, the class of uniformly hyperbolic systems do contain examples of the pathological behaviors we address: atypicality of shadowing and large linear responses. Moreover. our numerical examples, including our primary one the Lorenz '63 system, which is only singular hyperbolic are chosen to not all be uniformly hyperbolic systems, in order for our inferences to have wider applicability.
A.1 Proof of Proposition
I
Here, we complete a proof sketch that we begin in section 3 An element V E Tn,un M with @DSnn can be identified t 02 Tego M. We define a Function Fun T Ta M -> nn M as (na())+41
u, then, fz:" + De he is an orbit of F. Writing Fun(1) = Fan(() (0) + dFF + Nun Nun(v). where N is the nonlinear part of Fun: fixed point satisfies, (Id d.Fun (0)) (0) + Nun (v). To show that exists, we follow Theorem 18. 1.3 of [KK995) and prove that (w)) = (Id-dFF dF (0))
- F(xmp + v1). If v is fixed point of Fans that is, (u)
(w) M. We have ||Tmm (v) TTn Tw)]l]] < [[( dF (0)) I||| (u) (w) S C(eo, e1) [1) -Nn (w=)l with probability > 1 To see this, note that d.F..n (0) = dF v, and F' is a hyperbolic map. by assumption, and hence Lemma 18.1.4 of [KH955 applies. Next, assuming that supzEM || [ddFF(z)]] is bounded (i.e., dF is Lipschitz), we obtain that d.N.. is Lipschitz, and hence, [() Tm (w)ll < [[(Td dF (0) Here, max, (1l (1l [w,]]] < &0.
(0)) is contraction a compact subset of M.
(v)
Nun (w)ll S C(eo: 61) Ko - will with probability > 1 which is chosen independent of n. 6o, E and 6, and K is the Lipschitz constant of dN Thus, Ta is a contraction on a do ball around 0 E Ta M when K6 < E2, for some E2 > O. Since T (0)] < C'eg, E1 Ee) and for any in 6, ball around in 0 in M, [[Tnl < ((Gg.9)6 + 62 86, when cu, 61 and dn are such that C( (1 Co + (1 - 62)5g < og, we have that Tan maps dg ball around in 0 in Tn an M to itself. Thus, the unique fixed point (from contraction mapping theorem on the ón ball) lies within the ball.
B Experimental Details
B.1
Data
In this paper, the time series data from different chaotic systems are generated by simulating the ODEs and iterated function systems below. The ODEs were numerically integrated using the fourth-order Runge- Kutta solver from the torchdiffed 3 library [CRBD18), with an absolute and relative error tolerance of 10-
17

---

We split the simulated trajectories into training and test datasets. The first 10,000 data points are used as the training data and the last 8,000 points are the data. The time step size of the simulation is included in Table 2
ID: Tent maps As first examples, choose three types of tent maps defined in [W21]. that are shown to be examples exhibiting atypical shadowing orbits: the tilted map (5), the pinched map
(6). and the plucked map 7 We use n 3 and 0.2 and 0.8 in the plucked map following
CW21
* < 1 + 8,
F(r;)) =
I+e
(5)
(2 r) I-8
2
2 1
8 <
1+ 11ot7/1+=)3-4am
F(18))
4(2-x)
2 x 1
Fs,n(x) = mi,n(n,,,, As,, (x)),,,
0<<<2
(7)
where
2x
2(1 - x))
fa(x) = min(
2
< 1,
1
[.(2>]
0.5,
(x(x)
2 f2(22a)
x2 0.5,
2
(x) 0, (2"x 12"x()
x
+2
2
2
2D: Baker's map We use the following perturbation of the classical Baker's map from (CW2
2xx Ly/r|2n
F(xx,)):ss) =
s ssin + |x/=|2m mod 2T.
3D: Lorenz '63 We conduct extensive numerical experiments in this paper with the Lorenz '63 system [Lor63], with 10, 3 8/3, and
28, which we describe in section 3D: Rössler We the parameter setting 0.2, b 0.2, and Rös76 below, which the chaotic regime.
5.7 in the Rössler system
de' dt (xx 22,,
-y-z x+ay b+z(x-c).
4D: Hyperchaos Another test case we consider the hyperchaotic system below in dimensions from [Zhal7], where the parameter values are 16, b 40, 20, d for the system to show chaotic bchavior.
+ dz yz xz by
dip' di (x)),,,w1))=
c(x 2) + xy
cly w) + xt.
127D: Kuramoto-Sivashinsk To test Neural ODE's performance in learning high dimensional chaotic system, we generate the modified Suramoto-Sivashine (KS) system's solution defined below with second order finite difference scheme used in (BW14).
18

---

Du
du
θu u
10 (u +
0x
xx²
0x1
[0, 17
E [0, 00)
u00 1) = u 1,, 1) =
0u
0
8x 0--
8x xmm.
u(x, 0) (x())
B.2 Architecture
We try two neural network architectures for learning the vector field of the systems: simple Multi- Layer Perceptron (MLP) model and ResNet (HZRS16). For the MLP model, we use GELU as the activation function, and ReL.U for ResNet model. In addition, we experiment with adding Fourier layers (LKA*20) to represent the solution operator. We use the Latent SDE code from [WCD20 for results in section
Table 2: Hyperparameter choices
Chaotic Systems
Hidden layer width
Train, Test
Neural Network
Epochs
Time step
Layers
size
[10000,
Tent map
10000
N.A.
[256]
ResNet 500
8000]
[5000,
Baker map
10000
N.A.
[512]
ResNet 100
5000]
[10000,
Lorenz '63
8000
0.01
[512]
ResNet 500
8000]
[10000,
Rössler
10000
0.01
[512]
ResNet
500
8000]
[10000,
Hyperchaos 20000
0.001
[512]
ResNet
500
8000]
Kuramoto-
[3000,
3000
0.25
[512, 256]
MLP
Sivashinsky
3000]
B.3
Hyperparameter Search
When ux(x) is true vector field and vn, (3) is learned vector field by neural network, h, given solution æ, relative error can be defined as:
relative error(z) |1r((x) vn(a))
(8)
Using grid search, hyperparameter values that yield the lowest relative error in the vector field as
per (8) were chosen. Hyperparameter search results are shown Tables and 9 Hyperparameter
values that different for each system in the Table We use the AdamW [H17 optimization algorithm implemented in the PyTorch library for all our experiments.
In addition to MLP and Resnet, for Lorenz '63, we train with FNOs LKA+ 20 and latent SDE
[LCDD2]] In FNO network, we fix the number modes to 4, with batch size 200 with Fourier Neural layers. Further details on latent SDE are discussed in section C.8
B.4 Computing
Numerical experiments were conducted using Tesla A100 GPUs with 80GB and 40GB memory capacities. All experiments were completed in under one hour, with the exception of those involving the KS system.
htnps///yyocuuuuoollaaasaaeeeeer GELU..ttmm
19

---

C Additional numerical results
In this section, we present results that described in section 2 for testing the statistical aceuracy and generalization of learned models of the Lorenz '63 system. We also present additional results for the tent maps and the KS equation, described the previous section.
C.1
Loss
0.10
Train
Train Test
0.004
Test
0.08
0.002
0.06
0.04
0.000
4000
6000
8000
10000
4000
6000
8000
10000
Epochs
Epochs
Figure 3: Training and test loss of representative Neural ODE trained with MSE and Jacobian-
matching loss
C.2 Relative Error
Figures 4 and 5 illustrate the comparison of relative error trends, as defined in Eq. 8 for vector fields learned by Neural ODE (ResNet) trained with mean squared loss and Jacobian- matching loss (as defined in Eq. (3)). for two distinct dynamical systems.
10.0
1 x
l x11
Ixl:
4
1+b:
7.5
3
5.0
2.5
1
0.0
0
5
10 time
15
20
5
10 time
15
20
Figure 4: Comparison of relative errors in the vector fields of the Lorenz '63 system produced by a Neural ODE with a ResNet, trained using mean squared loss (Left), and Jacobian-matching loss (Right) defined in and ( respectively. The vector field is evaluated on random orbit.
C.3
Tent map
In Figure 6, we show learned models representative tent map perturbations from the preceding section. We observe that the Jacobian- Jacobian-matching loss leads to reasonably accurate representations of the maps. Yet, the LE computed along Iearned orbits differ significantly at s = 0.8, as shown in Table 5 For various other perturbed tent maps, show the computed L.s alongside the true LEs in Table for most of the maps, the Jacobian-matching leads to accurate LE predictions, consistent with the results of Theorem
C.4 Lorenz '63
Here we show identical results to Figure T but with ResNet-based Neural ODE models as opposed to MLPs used in Figure]] [ We find that Jacobian meathinn-matching training leads to superior performance in
20

---

300
15
Ii x|2
Ik x12
|x|
1x|p
200
10
100
15
20
20
time
time
Figure 5: Comparison of relative errors in the vector fields of the Rössler system produced by Neural ODE with ResNet, trained using mean squared loss (Left), and Jacobian-matching loss (Right) as defined in and respectively. The vector field is evaluated on a random true orbit.
2.00
2.00
Ther
1.50
1.50
3.00
1.00
0.75
0.75
0.25
0.25
0.00
030
2.50
0.00
0.30
Figure 6: Comparison of true plucked tent map with generated with Neural ODE trained with Jacobian-matching loss defined in Equation (3) When 0.8, we observe a failure mode of training with Jacobian-matee matching loss, possibly due to atypicality of shadowing orbits observed in (C21])
Details on hyperparameter setting and other statistical results can be found in Table 2 and
terms of distributional match with μ4, when compared to MLPs. Furthermore, our overall inference about MSE models leading to atypical orbits still holds, agnostic modeling and architectural choices.
MSE
40
40
Model 0.04
20
20
True
Model
-20 20
0
100
200
300
20
100
200
300
0.00-20
Figure 7: The first columns show orbits plane obtained from RK4 integration of the Lorenz vector field ([oor63), the Neural ODE, 'MSE Res', trained with mean-square loss, and the Neural ODE, 'JAC Res", trained with columns show the probability distribution orbits generated by the true (gray), 'MSE Res' (red) and 'JAC_Res Res' (blue) models. Experimental settings are in Appendix B
loss, respectively (see section The last two
C.5 KS equation
After spatial discretization and following the scheme of [W14], we obtain 227-dimensional dynamical system representing the KS equation, which we consider to be the ground truth map F. Figure 8 shows the solutions of the KS system over physical space and time (T). Since the original model and learned models are all chaotic, we expect two solutions of even slightly different models to diverge along the time axis, even when starting with identical initial conditions. We observe that this divergence is minimal for the JAC model, while the errors the MSE model grow strikingly quickly.
21

---

700
700
800
500
600
500
500
500
400
400
400
300
2
100
100
30
100
% 50
100
50
100
X
Figure 8: Solution plot of Luramoto-SSvashinkss system when number of inner nodes is 127 and c=0.4 (see [BW14] for the parameter c). True solution (left), solution of the Neural ODE with mean squared loss (T) (center column), solution of the Neural ODE trained with Jacobian-matching loss defined in (right)
C.6 Learning Lyapunov exponents
We report I. Lyapunov exponents learned with Neural ODE trained with mean squared loss and Jacohian- matching loss in T To compute Lyapunov exponents, al [GCLP13). We find that, for most systems considered, the LEs computed by JAC models match very well with the ground truth, while the MSE models sometimes capture the leading LEs but inaccurate in the rest of them.
QR algorithm of Ginelli et
C.7 Comparison Jacobian bacobian-matching loss training with unrolling dynamics
In this section, we report experiments on the Lorenz. '63 system with an unrolled loss function:
(x(x) = h ESA
(x))2 where ve (xr) = F (x) F² Fxx).
With A 10 timesteps of unrolling time, Figure However, atypical orbits are still produced for random initializations (Figure I0).
that the attractor is reproduced well, shown in
We experiment using two types of Unn. MLP and Resnet, and varying sequence lengths, k. As shown Table we observe that the learned negative Lyapunov exponent plateaued for 40. remaining between -9 and -10 (the true value being - 14.5). with no further improvement. Also, as & increases, we observe that Neural ODEs overestimate the positive Lyapunov exponent. Overall, unrolling seems to learn more accurate representations than the MSE model but less accurate representations than the JAC models. We also observe that the unrolling time needs to fine-tuned as a hyperparameter to achieve good generalization; a small perturbation can lead to training instabilities.
To understand these results, for short times, when compared the Lyapunov time, v, are in tangent spaces along the orbit. This yields the recursive relationship, 0t 2 'dF (x+-1) when O (1,,, negligible. Thus, the unrolled loss does contain Jacobian information implicitly although does not enforce the learned trajectory to be close to the true trajectory in C' -distance. We remark, speculatively, that Theorem 2 gives a possible explanation for why the unrolling loss performs better than a one-step loss (even in some practical climate emulators, e.g., FourcastNet (PSH* 22)) at learning the physical measure. In other words, our numerical results lend support to the central thesis of this paper: adding Jacobian information improves statistical accuracy. Table 4 shows the norm difference between the reproduced invariant statisties and the true statistics of Lorenz '63.
C.8 Latent SDE: experimental details
Here we present the experimental details and results of learning the Lorenz '63 system with latent SDEs (LWCD20) that described in Section 4
As training data, we use 1024 time series of the interval j0, 6) and timestep size 0.01. With smaller sample size, empirically observe large distributional mismatch (see Figure II We also observe
22

---

Table 3: Lyapunov Spectra learned by Neural ODE models trained on the MSE 1 for multi-step prediction.
yan
k
MLP
ResNet
10
[0.89, -0.0006, -5.53]
[0.77, -0.0080, -4.48)
20
[0.89, -6.349,66.31
[0.82, 0.0128, -5.191
30
[0.906, 0.0018.-7.18]
[0.89, -0.0372,
40
[0.96, -0.0587,-9.877
[0.92, -0.0655, -9881]
50
[0.96, -0.0922..22
Table 4: True vs. learned Lorenz '63 system: comparison of statistics. (W Wasserstein-1 Distance, A [}r, 22 Aa] set of LEs. fr: empirical distribution of orbit of length I'. The subscript NN indicates quantities computed using NN models. The variable refers to the sequence length used for training.
Norm Difference
Model
k
W (isoo, (NN,5000)
|| ANNI
((r)
(00
MLP
50
1.4184
4.4824 2.9456
1.3845
ResNet
50
0.5091
0.2408
40
40
N
N
20
20
-20
-10
10
20
X
X
20
40
40
N
N
20
20
-20
-10
10
-10 0
10
20
0
20
X
X
Figure 9: Comparison of the true phase plots of Lorenz '63 with phase plots of Neural ODE trained with the unrolled loss function (9). Initial condition of the long orbits is at [-9.1164, -3.3816, 33.7482). The first row shows the orbits on the xy, xz, and the yz plane obtained from RK4 integration of the Lorenz '63 system. The second row shows the orbits on the xy, xz, and the yz plane generated from Neural ODE trained with the loss with i 50.
that increasing decreasing the time length of the trajectories can lead to unstable training. As in the supervised learning methods, the time series are simulated using a fourth-order Runge- scheme (RK4). We follow the same architecture [ as used in (LWCD20, with GRU encoder, linear decoder, and the drift and diffusion functions in the prior and posterior processes modelled by MLPs.
We present the learned Lyapunov exponents Table 6 and a comparison of the learned and true empirical measures in Figure I1 When learning a deterministic system with a latent SDE, the diffusion coefficient the learned system is small but not exactly 0, and this results in a slight difference in the Lyapunov exponents ([GHL20) compared only using the learned drift term. The latent SDE model was also tested on the stochastic Lorenz attractor [CZ221] where it reproduced the 'bimodal' distribution of the trajectories.
lorenz
23

---

25
25
N5
-25
-20
20
-20
20
-25
25
X
100
100
N
N
0
100
X
Figure 10: Comparison of the true phase plots of Lorenz '63 with phase plots of Neural ODE trained
with the unrolled loss function (9). Initial condition of the long orbits is at [-15, -15, 5]. The first row shows the orbits on the xy, and the yz plane obtained from RK4 integration of the Lorenz '63 system. The second row shows the orbits on the xy, xz, and the yz plane generated from Neural ODE trained with the loss 9 with k 50.
Table 5: Chaotic systems and the true Lyapunov Spectra, the learned Lyapunov Spectra from Neural ODEs with MSE loss (T). and the learned Lyapunov Spectra from Neural ODEs with Jacobian- matching loss (3)
Lyapunov Spectrum
Chaotic Systems
Aruee
AMSE
AJAC
Tent map (tilted, s=0.2)
[0.3188]
[0.6722]
[0.3402]
Tent map (pinched, s=0.2)
[0.6849]
[0.6840]
[0.6836]
Tent map (plucked, s=0.2)
[0.6681]
[0.6763]
[0.6348]
Tent map (tilted, s=0.8)
[0.3188]
[0.3176]
[0.3402]
Tent map (pinched, s=0.8)
[0.6215]
[0.5770]
[0.6420]
Tent map (plucked, s=0.8)
[0.3199]
[0.5819]
[0.5013]
[0.9, O.
[0.87, 0.0091,
[0888, -0.0012,
Lorenz '63
-14.52]
.4.82]
-14.54]
[0.0665, -0.0004
[0.0008.-0.022
[0.0609, -0.0004
Rössler
-5.4112]
-1.4108]
-5.3808]
[4.0039, 0.0082
[4.1393, 0.0955
[4.3789, -0.1617
Hyperchaos
19.9972, -48.0205]
15,2120, -29.9480)
19.9974, -48.0205]
[0.3036, 0.2733,
0.1652, 0.1647,
0.29944 0.2622,
0,2592, 0.2257,
0.1540, 0,1524,
0.2293, 0.1990,
0.2050, 081888,
0.1443, 0.14
0.1701, 0.1584,
Kuramoto-
0.1649, 0, 0,16,,,
0.1336, 1262
0.1320, 0.1071,
Sivashinsky
0.1288, 0.1,
0.1236, 0.1
0.0912, 0.0724,
0.0992, 0.0776,
0.1141, 0. 1091,
0.0591, 0.0442,
0.0646, 0.0492,
0.1045, 0.0971,
0.0306, 0.0157,
0.0342
0.0985 ]
0.0023, ]
C.9 Experimental details of the learned Lorenz '63 models
24

---

Table 6: Lyapunov spectra computed from the learned latent SDE model evaluated on Euler- Maruyama and RK4 solvers. Both the learned drift and the near-zero diffusion vector fields are used Fuler-Maruyama, and only the drift vector field is used RK4.
Solver
Lyapunov Spectrum
Euler- Euler-Maruyama
[0.841. 0.0125, -11.88]
RK4
-0.44, -0.508, -11.75]
Table Results of search over hyperparameters (batch size, weight decay, hidden layer depth and width) in training Neural ODEs with MLPs (with fully connected and convolution layers). We train with mean squared loss using the AdamW optimization algorithm, with two values of weight decay: 10-3 and 10-1, and an adaptive learning rate with initial value of 0.001. For each hyperparameter combination, we show the test loss and the relative error the one-timestep predictions averaged over 8000 samples; choose the hyperparameter combination that results in the least relative error. The time step of the maps (both true and NNs) are set 0.01.
Batch Size
Hidden Unit
Layer
Full
1000
2000
10-3
10
10-3
10
10-3
10
3.0577,
0.1756,
1.0858,
0.0700,
0.0300,
0.0713.
256
2.02%
2.16%
9.16%
4.92%
8,45%
6.17%
0.0967,
0.0637,
0.0262,
3.5441,
0.2689,
0.2144,
3
512
2.11%
1.84%
7.71%
6.89%
5.89%
3.98%
0.8875,
16,9632
0.0217.
0.2425,
0.0300,
0.0488.
1024
1.32%
.3.01 %
7.89%
3.12%
8.45%
11.22%
0.0232,
8.2952,
0.7407,
0.2120,
0,0700,
3.9985.
256
1.70%
1.12%
2.48%
9.03%
4.92%
10.19%
0.0936,
6.7319,
243.8366,
24.6090,
3.2698.
0.1982,
5
512
2.22%
1.29%
15.01% %
17.84%
17.54%
7.44%
158.9892,
0.0234.
0.0472,
0.0713.
0.2425,
0.0332.
1024
1.03%
2.03%
4.94%
4.02%
3.12%
11.46%
0.0442,
0.0583,
0.8257,
28.2333,
0.1761,
0.3564,
256
1.13%
1.65%
8.72%
14.39%
6.32%
4.62%
0.5062,
0.0204,
4.5527,
0.0470.
170.5588,
0.0700,
7
512
1.03%
1.41%
11.59%
12.36%
8.51%
11.94%
162.6586,
0.0502,
0.0302,
3.7282,
117.4206.
0.0502,
1024
1.87%
1.15%
7.71%
19.24%
1.20%
1.15%
25

---

Librr 46
699
Density
0.22
0.0
132
Labee
Like
Density
010
3.34
6.3
3:4
0.30
X
y
Figure 11: Comparison of empirical distributions of the x, y, and z coordinates of the true orbits of the Lorenz '63 system (black) against those of the latent SDE (red). Top row: The latent SDE model sample points. Bottom row: The latent SDE model trained with the same sample size, 10,000, as the JAC_MLP and MSE_MLP (Neural ODE) models. We generate a trajectory over the interval [0, 50] for both the Lorenz '63 model and the learned latent SDE system. Gist: empirically observe that, when training with O (100), the latent SDE model results in a worse prediction the physical distribution compared to supervised learning based on iacobian-matching loss (3).
(see section 4 is trained with 614,400
26

---

Table 8: Results of search over hyperparameters (batch size, weight decay, hidden layer depth and width) training Neural ODEs with MIPP_skip (ResNets). We train with mean squared loss using the AdamW optimization algorithm, with two values of weight decay: 10-3 and 10-1, and an adaptive learning rate with an initial value of 0.001. For each hyperparameter combination, we show the test loss and the relative error in the one-timestep predictions averaged over 8000 samples; choose the peperaarreeeer combination that results in the least relative error. The time step of the maps (both true and NNs) are set at 0.01.
Batch Size
Hidden Unit
Layer
Full
1000
2000
10-3
10
10-3
10-4
10-3
10-4
0.5759,
0.1186,
0.4406,
0.2987,
0.1015,
0.1698,
256
2.60%
3.39%
9.60%
6.64%
4.07%
5.41%
0.1550,
0.0914
0.4885,
0.1510,
0.4082,
0.3319,
3
512
1.80%
2.08%
5.57%
4.47%
2.86%
6.84%
0,4721,
0.1177.
0.0433,
0.0668,
0,1681,
0,2515,
1024
2.23%
1.76%
2.11%
4.07%
5.15%
4.89%
0.5392,
0.1348.
0. 0.1694,
0.2117,
0.1862,
0.1568,
256
2. 2.14%
2.21%
3.39%
7.68%
3.89%
6.37%
0.0437,
0,0611, 1,
0,2823,
0.1005,
0.0934,
0.0822,
5
512
1.70%
1.31 %
4.87%
4.74%
5.97%
5.23%
0.1914,
0.2093,
0.2429,
0.0563,
0.4757,
0.0377,
1024
2.44%
2.08%
8.11% %
1.94%
2.60%
4.36%
0.1701,
0.1523,
0.8257,
0.9819,
0.1255,
0.3181,
256
1.91%
1.99%
8.72%
7.22%
9.65%
6.80%
0.2088
0,0806.
0. 1857,
0.3120,
0.1538,
0.4729,
7
512
2..8%
1.87%
4.64%
5.44%
4.62%
7.19%
0.0566,
0.0237.
0.6599,
0.0326.
0.1748,
0,3441,
1024
2.59%
2.16%
8.04%
5.07%
4.07%
5.09%
27

---

Table 9: Results of search over hyperparameters (batch size, hidden layer depth and width) training Neural ODEs with MLP skip (ResNets). We train with the aacobian-matchinn loss (3) using the AdamW optimization algorithm, with weight decay of 5 x 10-4, show the loss and the relative error in the one-timestep predictions averaged over 8000 samples; we choose the hyperparameter combination that results in the least relative error. The time step of the maps (both true and NNs) set at 0.01.
initial value of 0.001. For each hyperparameter combination, and an adaptive learning rate with
Neural Architecture
Hidden Unit
Layer
MI.P
MLP __sipp
500
1000
500
1000
256
0.0535, 0.22%
0.0031, 0.26%
0.3349, 0.96%
0.3689, 1.52%
512
0.0022, 0.26%
0.0256, 0.22%
0.0780, 0.93%
0.8011, 1.83%
1024
0.3905, 0.44%
0.0237, 0.82%
0.3349, 0.96%
0.2787, 1.80%
256
0.0060, 0.38%
1.9023, 0.28%
0.4635, 0.98%
0.5906, 2.30%
512
0.0300, 0.33%
0.0612, 0.59%
0.4211, 0.76%
0.2773, 2.30%
1024
0.0885, 0.83%
0.0232, 0.75%
0.1064, 0.69 %
0.7675, 2.23%
256
0.0297, 0.20%
2.8038, 0.48%
0.1032, 1.01 %
0.5467, 2.60%
7
512
0.0437, 0.10% %
0.011, 0.45%
0.0991, 0.93%
0.0943, 0.92%
1024
0.4413, 1.3%
1.2202,0.36%
1.1442, 1.24%
0.1917, 1.43%
28

---

NeurIPS Paper Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The abstract and introduction clearly state the claims, which match with the presented experimental results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
Limitations are discussed under Related Work, Conelusion and Appendix A.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]

---

Justification:
All theorems and proofs are clearly documented in the main text along with supplementary proofs in the Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Details on experiment setting and results are thoroughly documented under Appendix B.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
30

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: To preserve anonymity, we submit main part of the code, but not the full repository which has the all the pretrained model. We aim to publish the full code with pretrained model in public later on.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please the NeurlPS code and data submission guidelines (https://nips.cc
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
Polllliicciiiioiooiiiicooooiiiii for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
We have written all the necessary details under the Appendix B.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
Justification:
We discuss all of the details the metrics used under the Appendix B.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
31

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma enror bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: Details on computing resources are written under the Appendix subsection B.4
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: ptttps://neuripstt www.urriis.cccpppeppe httpp/pp....ppp..pttttttacstttaa
Answer: [Yes]
Justification: All of the authors have reviewed the NeurlPS Code of Ethics and the research conducted in the paper conform with the Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Discussion on broader impacts are under the Conclusion.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
32

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper does not pose any such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: Citation and versions are properly cited throughout the paper.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperswithcone com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
33

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: Details given under the Appendix and the main text well. Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: This paper does not involve research with human subjects nor crowedsourcing.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve research with human subjects crowdsourcing.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
34


Reviews:
[{"review_id": "ZZpaaY2Fjq", "reviewer": "Reviewer_sRDx", "summary": "The authors focus on analyzing why MSE loss fails to capture the physical behavior of dynamical systems. Narrowing their analysis to invariant ergodic systems, they provide theoretical insights on when generalization implies statistical accuracy. They propose that for models to be statistically accurate, they must reproduce dynamical invariants, which is achieved by ensuring that the learned model closely follows the true system's orbits. Specifically, the authors provide theoretical justification for why Jacobian matching loss can better capture statistical properties than MSE loss. Empirically, they verify their analysis on the Lorenz 63 system using different architectures, including MLP, Neural ODE, and FNO.", "questions": "Q1: Figure 1 is unclear, especially when the y-axis is not aligned in the 4&5 columns, and the message shown in Figure 1 seems to contradict the authors' theoretical claims, as the probability distribution of the modeled dynamics is still mismatched when the dynamics of the attract seems to be well learned.\n\nQ2: It's confusing when the authors interchangeably used notations $h$ and $F$ for the dynamics map. Is there a particular reason for switching between these notations, and could you elaborate on why?\n\nQ3: Line 183: It's not clear why the authors state that implementing Wasserstein distance is difficult when the map $h$ is chaotic.\n\nQ4: As Lyapunov exponents are calculated using local dynamics, could you show a comparison using some long-term evaluation metrics?", "limitations": "N/A", "rating": 7, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "The authors provide a thorough theoretical analysis and propose some useful notions to characterize models' ability to reproduce the dynamics statistical measure. Overall, the paper is well-written and easy to follow.", "weaknesses": "1. Jacobian loss (Eqn.3) considered in the paper is a special case of the Sobolev norm in [1]. The authors should consider extending their analysis to the Sobolev norm.\n\n2. Recent related works have also discussed the problem of MSE and shown multiple ways to improve upon MSE loss, e.g., the Wasserstein loss used in [2,3] and the theoretical analysis provided in [4]. Although some of these works might focus more on the empirical side, my understanding is that this problem has raised more attention than the authors' claim, and some comparison to the existing work is necessary.\n\n[1] Learning Dissipative Dynamics in Chaotic Systems (https://arxiv.org/pdf/2106.06898)\n\n[2] DySLIM: Dynamics Stable Learning by Invariant Measure for Chaotic Systems (https://arxiv.org/abs/2402.04467)\n\n[3] Training neural operators to preserve invariant measures of chaotic attractors (https://arxiv.org/abs/2306.01187)\n\n[4] On the difficulty of learning chaotic dynamics with RNNs (https://arxiv.org/pdf/2110.07238)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Implications for training with Sobolev norm\n\nIt 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.$ \nHere, 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. \n\nHence, 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.\nWe 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. \n\n### Related work\n\nThese 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]. \n\n### Clarification of Figure 1\n\n\nThe 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](https://app.gemoo.com/share/image-annotation/679582047652405248?codeId=vzaQe6gZzGngO&origin=imageurlgenerator&card=679582046868070400) into one figure, revised the caption and added a Gist (please see the link). Thank you! \n\n\n### Clarifying notation of map \n\nWe 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$. \n\n### Solving OT problems with ergodic measures of chaotic maps\n\nThank 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),$  \nas \n$$ \nW^1(\\mu, \\mu_{\\rm nn}) \\geq \\lim_{t\\to\\infty} \\dfrac{1}{t} \\sum_{n\\leq t} |g(F^n(x)) -  g(F^n_{\\rm nn}(x))|, \n$$\nfor 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). \n\n### LEs are statistical measures\nLEs 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))\\|.$$ \nWe 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!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mH6BESBd37", "reviewer": "Reviewer_L3o7", "summary": "This paper extends generalization results to models trained on dynamical data, especially Neural ODEs. The paper shows and then attempts to explain why Neural ODE trained without a Jacobian matching term fail to capture physical behaviour even when they have low generalization error. Under a generalization assumption, the paper shows that the Jacobian-matching loss can lead to statistically accurate models. The paper shows a number of experiments that show the behaviour stipulated by the theoretical results.", "questions": "Related to the above. In lines 264-265 the paper states that C0 generalization is insufficient for learning shadowing orbits. Is it an obvious fact? Further explanation would be useful. Also useful would be a comparison of a standard MLP minimizing loss 2 and a neural ODE minimizing the same loss.", "limitations": "The paper does not have a section on limitations. I think a discussion of limitations would add to the value of the paper.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is technical but well written and ideas are clearly explained.\n\nThe mathematical formulation is quite clean and targets an important problem of generalization in ODE models for dynamical data.\n\nThe explanation of generalization for Neural ODEs is I think novel.", "weaknesses": "I think the main weakness is assumption 1. There should be some explanation of the difference between C1 and strong C1 generalization. Furthermore, there should at least be some justification of why this should hold and when we could not expect it to hold.\n\nAnother possible problem I see is this: The paper explains why minimizing loss 3 (with the jacobian) implies statistical accuracy. I don’t think that it properly explains why minimizing loss 2 does not imply statistical accuracy. I would understand why a standard MLP might not match derivatives. But I have a harder time understanding why a neural ODE with low generalization error isn’t able to do that. Some more intuition about this would be useful.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### When Assumption 1 is expected to hold\n\nA 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.\n\n### Why minimizing MSE loss does not yield statistical accuracy\n\nThis 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.\n\n### Neural ODEs are neural parameterizations of the vector field\n\nNeural 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. \nOur 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.\n\nIntuitively, 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. \n\n### Limitations\n\nWe 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.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "5oqZDCWPRb", "reviewer": "Reviewer_vAp3", "summary": "The paper addresses the problem that classical ERM training of dynamical systems models often fails to capture invariant measures of the observed dynamics, even when test errors are low. The authors use ergodic theory to explain this failure from a theoretical viewpoint. They further demonstrate that incorporating Jacobian information during training enables much better reconstructions of invariant properties and validate their approach on various neural flow operators trained on several common benchmark chaotic systems.", "questions": "- how exactly are the empirical measures $\\mu$ and $\\mu_{NN}$ and their $W_1$ distance computed (for Table 1 e.g.)?\n- Could the authors include comparisons of the Jacobian matching loss to the more common training method of unrolling the dynamics for at least some time $t$? \n- How can one draw practical consequences from the theoretical findings?", "limitations": "The authors do mention limitations of their approach (e.g. that the Jacobian matching loss is hard to implement in practice when ground truth is not known), but do not address how their findings might be used in real-world data settings.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- I think the authors study an interesting problem, i.e. when classical ERM for dynamical models can still lead to a reconstruction of invariant properties of the observed dynamical system, as classical ERM (i.e. here simply “one-step-ahead” predictions) comes with benefits of decreased training time and less training difficulties\n- The authors use dynamical systems and ergodic theory to accompany and explain the practical behavior of several neural flow operators (N-ODEs, RNNs, MLPs, etc.) when training them using ERM on observed dynamical systems\n- The paper is generally well written", "weaknesses": "- l. 53 - 56: Indeed there are also connections of training recurrent architectures to model chaotic DS and exploding loss gradients in gradient descent based training [1]. Moreover, recently methods for dealing with these training instabilities have been introduced and proven useful [1, 2].\n- There is also a recent paper which discuss (out-of-domain) generalization in dynamical systems reconstruction [3] which might make sense to add to related work, as it also introduces measures to assess generalization with a focus on evaluating invariant measures of the underlying dynamics and addresses the problem that classical ERM framework is not enough to assess generalization in dynamical systems reconstruction\n- While it is good to see that the Jacobian matching loss improves the reconstruction of invariant measures of the underlying dynamical system, it’s application in practice is fairly limited as Jacobians have to be estimated from data if the ground truth vector field is not known.\n- I see that the authors want to study the case where training is only perform using “one-step-ahead” predictions (i.e. no unrolling of dynamics), however, I think considering the case of unrolling dynamics as a comparison to the Jacobian matching loss would be helpful as 1) unrolling dynamics is often enough to get a good reconstruction of the invariant measures of the observed attractor and 2) it is easy to apply in practice (no need for extra knowledge of the Jacobians). I.e. I think the more interesting question is whether knowledge of the Jacobians outperforms unrolling of dynamics during training significantly\n- I think section 4, i.e. dynamic generative models, could be written a tad clearer, i.e. more explicit in terms of what exactly is now assumed to be stochastic in the specific VAE architecture and which part is deterministic. Maybe a Figure demonstrating the architecture setup in favor of moving e.g. proof of Theorem 1 into the Appx.\n- To me the paper lacks a clear connection of the theoretical findings to the practical real-world setting, i.e. how these insights influence how we should train models to reconstruct real-world dynamics\n\t\nI really like the idea and insights of this paper, but given the Weaknesses above, I do think it lacks actual practical relevance and the novelty is also very much limited to the (seemingly impractical) theoretical findings.\n\nMiscellaneous comments:\n- Figure 2 readability could be improved by removing the legend from within the plots and moving it to the outside (e.g. right, or between the plots), as the legend is the same for both plots\n- Hyperlink of Table 3 jumps to section 3\n- typo l. 192: “shdaowing”\n\n[1] [Mikhaeil et al. (NeurIPS, 2022), On the difficulty of learning chaotic dynamics with RNNs](https://proceedings.neurips.cc/paper_files/paper/2022/hash/495e55f361708bedbab5d81f92048dcd-Abstract-Conference.html)\n\n[2] [Hess et al. (ICML, 2023), Generalized teacher forcing for learning chaotic dynamics.](https://proceedings.mlr.press/v202/hess23a.html)\n\n[3] [Göring et al. (ICML, 2024), Out-of-domain generalization in dynamical systems reconstruction.](https://arxiv.org/abs/2402.18377)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Related work on invariant measures\n\nThank 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. \n\n### Practicality of Jacobian loss and implications of theory\n\nThis 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.\n\nFurthermore, 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.\n\n### Unrolling dynamics\n\nWe 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](https://ibb.co/YBKhPjT), but [atypical orbits are still produced](https://ibb.co/dr9bPr6) 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.\n\nOverall, 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.\n\nTo 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.\n\n### W1 distance\n\nWe 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.\n\n### Edits\n\nMany 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.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IZMHG6TC7X", "reviewer": "Reviewer_PvwR", "summary": "The present manuscript concerns the use of neural networks to fit physical dynamical systems of generic kind (including and focusing on chaotic maps). It is shown that adding the information of the Jacobian of the dynamical map in the supervised learning process leads to better performances.", "questions": "1) I believe that it must be carfully stated that it is assumed that the map F does not depend on time. This is an important assumption which may be relaxed.", "limitations": "The main and severe limitation is that the improvement training method needs the knowledge of the Jacobian of the ground truth dynamical system. It is unclear how this can be estimated from data contain just points in the orbits of the maps and what are realistic settings where this is known.\nI believe that this is a severe limitation of the work. This does not spoil the result in the sense that the manuscript clearly shows that the information coming from the Jacobian is essential, but on the other hand it seems to me that this information is hard to get.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "Originality: the work builds on previous ideas and problems. The original part (as claimed by the authors) concerns adding the information on the Jacobian of the dynamical map to the supervised learning process. The author show that this significantly improves the statistical properties of the points generated by the fitted dynamical systems.\nQuality: The paper is very well written, well structured and clear. Section 4 is very clear.\nSignificance: I believe that the problem tackled in this work is a relevant one and the results show an interesting perspective on it.", "weaknesses": "I believe that the most important weakness of the manuscript relies on the limitations of the results.\nThe paper clearly shows that having information on the Jacobian of a dynamical map is essential for reproducing the statistical properties of the generated orbits. However such information is typically lacking.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Results require full Jacobian but it can be approximated with available derivative information in practice\n\nYou 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.\n\n### Focus on autonomous systems\n\nIn 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!", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: All reviews agreed on the premise of the paper: existing approaches to learning dynamical systems using MSE losses fail to generalize, and the paper introduces a Jacobian-matching loss that can learn underlying properties such as orbits and invariants. The major weakness identified by two reviewers is that the Jacobian is often not possible to get from data and some other alternative signal would be needed in most practical situations. The discussion resulted in extra experiments that other practical alternatives still do not perform as well as the Jacobian-matching method, and there was agreement that the proposed method is still a valuable contribution. All reviewers agreed on the theoretical soundness of the paper.

After discussion, all the reviewers all agreed that this paper should be accepted. One reviewer was borderline, but after the discussion, the authors promised to add an additional experiment in the revision:
> “ I think the paper is of great value ...I will raise my score to 5 with the assumption that the authors will include all additional experiments in the final version (reevaluate LE spectrum calculation, making sure the true spectrum is reported, estimating Jacobians from data alone, addressing computational complexity).”

The other three reviewers did not think that these results were necessary for their accept recommendations. The AC thinks the additional results obtained in the discussion period would be easy to include in the revision and not dramatically change the paper. As such, the paper can be accepted.
Decision: Accept (poster)