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

Authors: jayjay  park
nicole tianjiao yang
nisha  chandramoorthy
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 16:12:22

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.

Main paper: 
 introduction: 
 Learning a dynamical system from time series data is a 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 in downstream computations such as optimization and uncertainty quantification. The field of physics-guided machine learning [vdGSB + 20, FO22, LK22, KKLL21, RPK19, KKL + 21] has emerged in response, rapidly integrating neural networks into data-driven modeling and prediction workflows for a wide variety of complex dynamics, from geophysical fluid flows to phase transitions in materials (see [KMA + 21, CCC + 19] for surveys). Yet rigorous generalization analyses of neural parameterizations in these applications, wherein the underlying dynamics can exhibit chaotic behavior, have been underexplored. Figure 1: A random orbit on the x-z plane obtained from RK4 integration of the Lorenz vector field ( [Lor63] (first column), Neural ODE, 'MSE_MLP', trained with mean-squared loss (second column) and Neural ODE, 'JAC_MLP', trained with Jacobian loss (third column). The last column shows the empirical PDF of the orbit generated by the true (gray), 'MSE_MLP' (red) and 'JAC_MLP' (blue) models. Experimental settings and additional results are in Appendices B and C respectively. 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/statistical behavior accurately.Here we investigate data-driven neural parameterizations of chaotic ODEs/PDEs and maps (discretetime dynamical systems) motivated by the following observation: a neural representation that is learned well, i.e., with a small generalization error, can still produce unphysical long-term or ensemble behavior. Figure 1 (left) shows the classical Lorenz '63 chaotic attractor [Lor63] plotted using a long random trajectory/orbit (time integration of the Lorenz equation vector field starting from a random initial condition, which is indicated by a '+' sign). The second column shows an orbit from a neural network model, 'MSE_MLP', which minimizes the mean-squared 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 1 (column 2) shows that the learned 'MSE_MLP' Neural ODE [CRBD18] model produces an atypical orbit -an orbit different from almost every orbit of the true system -for the same randomly chosen initial condition. As a result, the learned empirical distribution is not close to the physical distribution -that of almost every true orbit, as shown in Figure 1 (column 4). On the other hand, the 'JAC_MLP' model, which is obtained by minimizing the mean-squared error in the vector field and its first derivative (Jacobian matrix), reproduces the Lorenz '63 attractor and the physical distribution on the attractor (Figure 1 , column 4). The 'JAC_MLP' also captures all the Lyapunov exponents (LEs) -measures of asymptotic stability to perturbations, which are 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 of 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 m samples from a time series, {(x t , x t+1 )} t , 0 ≤ t ≤ (m -1), how can we learn a model, F nn , such that, i) x t+1 ≈ F nn (x t ), for all time t, and ii) the underlying distribution of the states x t and other dynamical invariants such as Lyapunov exponents are reproduced by orbits of F nn ? It is widely accepted that learning such an F nn involves matching orbits of F nn with x t over large t during training. However, small errors propagate over orbits, by definition, in a chaotic system, leading to training instabilities. In response, a vast literature has been dedicated to developing sophisticated training models based on RNNs [PLH + 17, PWF + 18, RM21], operator learning [LLSK + 22] and regularizations [LG22, FJNO20] .We focus instead on the empirical risk minimization (ERM) for F nn that does not explicitly use the temporal correlations/dynamical structure in the data, avoiding training instabilities. Thus, we attempt to characterize when an elementary regression problem F nn 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 a model's ability to reproduce dynamical invariants.• We develop new dynamics-aware generalization bounds for minimization of errors in the C r , r = 0, 1, topology.• We leverage shadowing theory from dynamical systems to rigorously characterize failure modes in learning statistically accurate models. 
generalization of parameterized ergodic dynamics: 
 In this section, we motivate, through illustrative examples, the need for a dynamics-aware 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 a discrete-time dynamical system, F, is a function on a closed and bounded (compact) set, M ⊂ R d , 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 ∈ C 1 (M ) is a differentiable function on M. We denote by F t , t ∈ Z + , the iterates of the dynamics, or the compositions of F with itself t times, i.e., F t = F • F t-1 . An orbit or trajectory of F t starting at an initial state x ∈ M is the sequence {F t (x)} t∈Z + . In practice, F may be a numerical ODE solver that approximates the continuous-time solutions, φ t , t ∈ R + , of the ODE (written in dynamical systems notation2 ):dφ t (x)/dt = v(φ t (x)), where v : M → T M is the true vector field describing the governing equations. Fixing some τ ∈ R + , F := φ τ . We say a map F is chaotic if there exists a subbundle of T M , called the unstable subbundle, where infinitesimal perturbations grow exponentially under the linearization (Jacobian map), dF t , of the dynamics. The true map F generates a deterministic, autonomous system (see Appendix A).Physical measures. A probability measure µ :M → R + is a physical measure [You02] for the dynamics F if it is a) F -invariant, b) ergodic and c) observable through F . A measure µ is observable if time-averages along any orbit starting from a randomly chosen initial point converge to constants that are expectations (phase space average) with respect to µ. That is, for anyf ∈ C(M ), (1/T ) t≤T f (F t (x))T →∞ ----→ Ef (x), for any initial point x chosen Lebesgue a.e. on a set U ⊆ M. The orbits starting almost everywhere on the basin of attraction, U, asymptotically enter a set, Λ, called the attractor. In dissipative chaotic systems, the attractor, Λ, which is the compact support of the physical measure µ, has Lebesgue measure 0. Consequently, µ may not have a probability density, that is, µ is not absolutely continuous, or is singular, with respect to Lebesgue measure.Neural ODE. Introduced in [CRBD18] , a Neural ODE, denoted here by, v θ : M → R d , with parameters, θ, is a vector field represented by a neural network. The vector field can be time integrated to obtain solutions φ t θ : R d → R d , t ∈ R + to the ODE, dφ t θ (x)/dt = v θ (φ t θ (x)). As noted in the introduction, suppose we have n distinct pairs S = {(x i , F (x i ))} i∈[n] , 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, ℓ, that can be chosen to be a square loss, e.g., ℓ(x, φ τ θ ) = ∥φ τ θ (x) -F (x)∥ 2 . That is, we solve for θ that minimizes the training loss, (1/n) x∈S ℓ(x, θ). 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 → φ τ (x), as F nn , or neural representation of the target map, F. Data and optimization. Suppose we use an m-length orbit as our training data, i.e., x i+1 = F (x i ), then, the data are not, strictly speaking, independent. However, a feature of chaotic systems is an exponential decay of correlations, and so training data of the form, {(x i , F ω (x i ))} starting from an initial condition x 0 ∼ µ, can be thought of as iid, for a large enough ω ∈ N. In that case, the learned map F nn represents the function F ω . During optimization, infinitesimal linear perturbations, will have to be evolved for time ω. Since adjoint solutions blow up exponentially, a longer ω will lead to training instabilities. To avoid numerical difficulties in training and focus on issues surrounding generalization, we choose τ := δt, a small time step, to define the target map F and learn a 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,EQUATIONwhere F nn := φ δt is a neural network representing the map. For a map h : M → M, we define the training and generalization errors in the usual way:EQUATIONFigure 2 : Learned and true LEs computed over 30,000 time steps using the QR algorithm of Ginelli et al [GCLP13] 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 ( 1) and (3) respectively.where the expectation is over the data distribution, which may be µ or any initial probability distribution 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 dφ t (x)/dt = v(φ t (x)), where the vector field v is given by v 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 µ 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 if 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 if the training data include transient dynamics off the attractor.(x) = [σ(y -x), x(ρ -z) -y,Lyapunov exponents, roughly speaking, measure the asymptotic exponential growth/decay of infinitesimal 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 to µ. 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 a center direction, or the vector field v itself) and one negative LE. The ground truth map F, which is a time δt approximation of the flow φ t 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 [GCLP13] . The Neural ODE model, marked 'MSE', produces a reasonably close approximation of the ground truth value (≈ 0.9) for the positive LE but obtains an incorrect approximation of the stable LE (true value ≈ -14.5). We note also that standard deviations in the LE values computed over different orbits is quite large, indicating some orbits with atypical behaviors. We show the ℓ 2 error in the computed LEs by the Neural ODE model in Table 3 , and full details in Table 5 . Jacobian-matching. We now consider Neural ODE models trained with the following loss function.EQUATIONwhereF nn : M → M is a Neural ODE map, dF nn (x) : T x M → T x M its (d × d) Jacobian matrix atx, and λ > 0 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:EQUATIONWe 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 ℓ, 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 1 , a Neural ODE trained with the Jacobianmatching loss produces an attractor (column 3) that is visually similar to the Lorenz attractor. The empirical distributions match closely with the ground truth, as shown in column 5 of Figure 1 and in Table 3 , where the model is marked with a 'JAC', while one trained with the loss ℓ in (1) is indicated with an 'MSE'. The LEs and statistical averages of the state are all accurately represented, although the model is not explicitly designed to learn temporal patterns in the data. The MSE models learn vector fields with comparable accuracy with the Jacobian-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 is simply stated as an ERM for regression, the generalization errors R(F nn ) and R λ (F nn ) cannot determine whether F nn 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-accurate 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 dynamics-aware generalization bounds for Neural ODEs.Our ultimate goal is to minimize a statistical loss function rather than (1) or (3) since this measures how accurately a model h can reproduce ergodic averages associated with F. For instance,ℓ stat (x, h) = sup f ∈Lip 1 |E[f (x)] -lim T →∞ (1/T ) t≤T f (h t (x))|. This is the Wasserstein (W 1 ) distance (expressed in its dual form) between the ergodic measure µ and the ergodic measure associated with the orbit of h, a learned model, starting at x, limT →∞ Unif{x, h(x), h 2 (x), • • • , h T (x)}.As noted, ERMs for this loss function do not have a straightforward implementation when the map h sought is chaotic. Hence, we seek conditions under which solving ERMs defined with losses (1) and (3) can still minimize ℓ stat . To understand when solving regression problems for F nn lead to statistically accurate physical representations, we first assume that a notion of "shadowing" applies to F.Definition 1. We say that the shadowing property applies to a map F if for any δ > 0, there exists an ϵ = O(δ) so that for any map G with ∥G-F ∥ 1 := sup x∈M (∥G(x)-F (x)∥+∥dG(x)-dF (x)∥) ≤ ϵ, there exists a map τ : M → M close to the identity such that ∥G t (τ (x)) -F t (x)∥ ≤ δ for all t.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 uniform-in-time shadowing is a classical result for a mathematically ideal class of chaotic systems, called uniformly hyperbolic systems (see Katok and Hasselblatt [KKH95] Ch 18; Appendix A). For a textbook presentation and extension of shadowing for dynamical systems with some hyperbolicity, see [Pil06] .For a uniformly hyperbolic F (see Appendix A), we now assume that a neural representation F nn of F trained with n samples generalizes well in terms of C 1 -distance. That is, an ERM solution for the loss 3 ('JAC' models in sections 1 and 2) generalizes so that R λ (F nn ) is small.Definition 2 (C 1 generalization). Given δ > 0, there exist E 0 , E 1 > 0 and a function (δ, E 0 , E 1 ) → τ (δ, E 0 , E 1 ) ∈ N such that E x∼µ ∥F (x) -F nn (x)∥ ≤ E 0 and E x∼µ ∥dF (x) -dF nn (x)∥ ≤ E 1 for all m ≥ τ with probability ≥ 1 -δ over the randomness of the training data from µ m .We now make an optimistic assumption on a learned model, F nn , that satisfies the above definition of C 1 generalization. Using Hoeffding's inequality, we know that for any δ 0 > 0, with probability at least 1 -δ 0 over the randomness of x, ∥F (x) -F nn (x)∥ ≤ E 0 + (sup x∈M ∥F (x) - F nn (x)∥) log(2/δ 0 ) and ∥dF (x) -dF nn (x)∥ ≤ E 1 + (sup x∈M ∥dF (x) -dF nn (x)∥) log(2/δ 0 ). Given δ > 0, let ϵ 0 := 2E 0 and ϵ 1 := 2E 1 . Fixing δ 0 > 0, suppose that the trained model F nn is such that (sup x∈M ∥F (x) -F nn (x)∥) < ϵ 0 /(2 log(2/δ 0 )) and (sup x∈M ∥dF (x) -dF nn (x)∥) ≤ ϵ 1 /(2 log(2/δ 0 )). Taking a union bound, with probability > 1 -(δ + δ 0 ), ∥F (x) -F nn (x)∥ ≤ ϵ 0 and ∥dF (x) -dF nn (x)∥ ≤ ϵ 1 .We enhance this inequality to obtain a stronger assumption on F nn . Assumption 1 (C 1 strong generalization). Given δ > 0, there exist ϵ 0 , ϵ 1 > 0 and a function(ϵ 0 , ϵ 1 , m) → n(ϵ 0 , ϵ 1 , m) ∈ N such that ∥F (F t nn (x)) -F t+1 nn (x)∥ ≤ ϵ 0 and ∥dF (F t nn (x)) - dF nn (F t nn (x))∥ ≤ ϵ 1 for all t ≤ n, m ≥ τ (δ, ϵ), and n → ∞ as m → ∞, with probability ≥ 1 -δ 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 a spin-off time to achieve a state on the attractor. Given a tuple, (ϵ 0 , ϵ 1 ), an orbit with initial condition x that satisfies, ∥F (F t nn (x)) -F t+1 nn (x)∥ ≤ ϵ 0 and ∥dF (F t nn (x)) -dF nn (F t nn (x))∥ ≤ ϵ 1 for all t ≤ m will be referred to as an (ϵ 0 , ϵ 1 ) orbit. That is, at each time, the neural representation F nn is (ϵ 0 , ϵ 1 )-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 (ϵ 0 , ϵ 1 ) orbit. Proposition 1 (Shadowing). Let F nn be an approximation of F that satisfies the C 1 strong generalization (Assumption 1). Given any δ > 0, there exist ϵ 0 , ϵ 1 , n such that every (ϵ 0 , ϵ 1 ) orbit is δ-shadowed by an orbit of F. That is, there is a true orbit, say, {F t (x)}, corresponding to every orbit,{G t (x ′ )} such that ∥F t (x) -G t (x ′ )∥ ≤ δ, for all t ≤ n.See section A.1 for the proof. Let {x nn t } t≤n be an n-length orbit of F nn , i.e., x nn t+1 = F nn (x nn t ). We use T n,nn M to denote the direct sum ⊕ n i=1 T x nn i M of tangent spaces along the orbit, x nn t . The proof follows Theorem 18.1.3 of [KKH95] to apply contraction mapping on a compact ball in T n,nn M . The above result defines, for each (ϵ 0 , ϵ 1 )-orbit, {x nn t } t , a shadowing orbit, x sh := {F (x nn t + v t )} t , where v = ⊕ t v t is the fixed point of the contraction map in the proof (section A.1). Let µ sh n (x nn 0 ) be the empirical measure defined on the shadowing orbit corresponding to an (ϵ 0 , ϵ 1 )-orbit, {x nn t } t , i.e.,µ sh n (x nn 0 ) = Unif{x sh 0 , • • • , x sh t , • • • , x sh n-1 }.A shadowing orbit is indeed an orbit of the true map F, but, unexpectedly, it may be atypical for the physical measure, µ. That is, for an atypical shadowing orbit, the time average, (1/n) t≤n f (x sh t ) does not converge to the expected value E x∼µ f (x), as n → ∞. This means that the Wasserstein distance, W 1 (µ sh n (x nn 0 ), µ), does not converge to 0 as n → ∞.Given an initial condition, x nn 0 , of an (ϵ 0 , ϵ 1 )-orbit, the corresponding shadowing orbit may be typical with some probability (over the distribution of x nn 0 ), and this probability of finding typical shadowing orbits is a property of the true dynamics, F. When µ, the physical measure of F, is highly sensitive to perturbations of F, (see [Rue09, Bal14] 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 1 generalization (see Definition 2). A neural model F nn which generalizes well in the C 1 sense (Definition 2) can be thought of as a C 1 -smooth perturbation of the true dynamics F. If smooth perturbations of F can cause a large change in µ, then even when the training size m → ∞, and the orbit length n → ∞, W 1 (µ sh n (x nn 0 ), µ), 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 us a characterization of statistically accurate learning.Theorem 1 (Statistically accurate learning). Let F nn be a model of F that satisfies C 1 strong generalization. In addition, let F nn and F be such that for any δ > 0, there exists an ϵ 2 > 0 so that lim n→∞ W 1 (µ sh n (x), µ) ≤ ϵ 2 with probability (over the randomness of x) ≥ 1-δ. Then, for any δ > 0, there exists an ϵ > 0 such thatlim n→∞ W 1 (Unif{x, F nn (x), • • • , F t nn (x), • • • , F n nn (x)}, µ) ≤ ϵ with probability ≥ 1 -δ. Proof : Let µ nn n (x) := Unif{x, F nn (x), • • • , F t nn (x), • • • , F n-1 nn (x)} be the empirical measure of an n-length orbit of F nn starting at x. Given a δ > 0, we choose (ϵ 0 , ϵ 1 ) so that C 1 strong generalization is satisfied with probability ≥ 1 -δ/2. Thus, an initial condition x is such that {F t nn (x)} t≤n is an (ϵ 0 , ϵ 1 ) orbit, where the tuple (ϵ 0 , ϵ 1 ) is as defined in Proposition 1, with probability ≥ 1 -δ/2. Applying Proposition 1, we have, for any 1-Lipschitz function f :M → R, (1/n) t≤n |f (F t nn (x)) -f (F t (x))| ≤ (1/n) t≤n ∥F t nn (x)) -F t (x)∥ ≤ δ/2. Taking a supre- mum over f , W 1 (µ sh n (x), µ nn n (x)) ≤ δ/2, with probability ≥ 1 -δ/2.By assumption, there exists some ϵ 2 > 0 such that lim n→∞ W 1 (µ, µ sh n (x)) < ϵ 2 with probability 1 -δ/2. Hence,lim n→∞ W 1 (µ, µ nn n (x)) ≤ lim n→∞ W 1 (µ sh n (x), µ nn n (x)) + W 1 (µ, µ sh n (x)) ≤ ϵ 2 + δ/2 := ϵ,with probability > 1 -δ, using triangle inequality and taking union bound. This result explains why training to minimize Jacobian-matching loss (3) can lead to statistically accurate models, even though, long-time temporal patterns in the data are not learned explicitly by regression for the one-time map F. Since C 0 generalization (i.e., small errors in(2)) is insufficient for learning shadowing orbits, and thus for Proposition 1 and Theorem 1 to hold, the models trained on MSE loss 1 are not expected to learn ergodic/statistical averages with respect to µ.When shadowing orbits are atypical with high probability, we observe numerically that C 1 generalization, i.e., training with Jacobian-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 learned neural representations with Jacobian-matching do have good C 1 generalization (Figure 6 ), but do not exhibit good statistical accuracy and learn incorrect Lyapunov exponents (see Table 5 , Plucked Tent map).Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? • The method for calculating the error bars should be explained (closed form formula, call to 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. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 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 out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
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 Jacobian-matching information can yield samples 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, we can use several popular generative models for our target physical measure here, including scorebased methods [SE19, SSDK + 20], Variational Autoencoders (VAE) [RM15] or normalizing flows [RM15, 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 to faithfully represent both the underlying dynamics as well as the physical measure.In a dynamic generative model, we approximate F t with a stochastic map, that can written as,F ls := f θ • Φ t ϕ • g ϕ ,where the subscript ls stands for "latent SDE". Here, the function g ϕ : R d → R d l , with learnable parameters, ϕ, is a (possibly stochastic) embedding from the data to latent space (R d l ), such that, g ϕ♯ µ = q ϕ,0 . Recall the pushforward notation (♯), i.e., if x ∼ µ, then, g ϕ (x) ∼ q ϕ,0 . The dynamical system Φ t ϕ : R d l → R d l acts on the latent space, and defines a sequence of pushforward distributions Φ t ϕ♯ q ϕ,0 = q ϕ,t . A special case of this setup is a "latent ODE", where Φ t ϕ is a 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, Φ t ϕ is a solution map of a Neural SDE:dΦ t ϕ (z) = w ϕ (t, Φ t ϕ (z))dt + σ ϕ (t, Φ t ϕ (z)) • dW t .Here the drift term, w ϕ , and the diffusion term, σ ϕ , are represented as neural networks. The decoder f θ : R d l → R d is a deterministic map that defines the conditional, f θ♯ q ϕ,t = p θ (•|Z t ). 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, z → q ϕ,0 (z|X 1:m ), of Z 0 given X 1:m := {F t (x)} t≤m is modeled as a Gaussian distribution whose learnable parameters are also denoted by ϕ. Similarly, the conditional p θ (•|Z t ) is again modeled as a Gaussian with parameters θ. The parameters, ϕ and θ respectively, of the encoder and the decoder are trained by maximizing the following dynamic version of the evidence lower bound (ELBO):ℓ ls (X 1:m , (ϕ, θ)) := m t=1 E zt∼q ϕ,t (•|X1:m) [-log p θ (x t |z t )] + KL(q ϕ,0 (•|X 1:m )∥p Z0 )],where the prior p Z0 follows a standard Gaussian distribution in the latent dimension. For alternatives to the ELBO objective above, such as the Wasserstein-GAN objective, we refer the reader to [Kid22] .We now evaluate both the learned dynamics, F ls , and the learned generative model, p θ , which approximates µ. In Figure 11 , we present the empirical distribution of a 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-approximated. First, even though the system is deterministic, the learned Φ t ϕ with minimum generalization error (Eℓ stat ) encountered in the hyperparameter search (Appendix C.8) is not, i.e, the diffusion term σ ϕ is not zero. The learned stochastic map, F ls , produces an incorrect stable LE for the Lorenz system (≈ -11.8), while the leading unstable LE matches reasonably well.Since the map F ls , or its underlying vector field, on the latent space is not unique, we may obtain maps that do not preserve dynamical structure or invariants, even if the generated samples from p θ approximately capture µ. As noted in section 2, the physical measure µ 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 [BCV13] ). 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 high-dimensional chaotic systems. Even in the Lorenz '63 system, we require O(10 6 ) samples for training reasonably accurate model in Figure 11 , while the Jacobian-matching (Figure 1 column 5) produces smaller Wasserstein distances with fewer samples (10 4 ), and a simpler regression problem as opposed to variational optimization above. 
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 R λ 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 1, 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 R λ 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://github.com/ni-sha-c/stacNODE. 
related work: 
 Neural ODEs and generalization. Introduced as continuous-time analogues [Wei17, HR17, CRBD18] of Residual neural networks (ResNets) [HZRS16] , Neural ODEs [CRBD18, HR17] offer a vector field parameterization that can be time-integrated with an ODE solver. Augmented Neural ODEs [DDT19, RCD19, DA22] demonstrate improved expressivity for complex dynamics, and some universal approximation results appear in [ZGUA20, Kid22] . Neural ODEs [CRBD18] as normalizing flows [RM15] 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, FJNO20, PER23a, PER23b, PMSR21, GBD + 20, MSKF21]. Our focus is different: we characterize when an elementary Jacobian-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 physicsinformed machine learning[VAUK22, LKB22, PABT + 21, HZB + 21, vdGSB + 20, WY21, 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 [HKUT20, PABT + 21, HZB + 21, RBL22, VAUK22, vdGSB + 20, YHP + 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 + 19, SLST17, LSGW + 23] and both faster training and inference making them suitable for optimization [BB21, RBL22, LP21, JD21] and inverse problems [HVT23, AEOV23] . Since the fundamental regression problem we study underlies both hybrid and data-driven methods, we provide insight into dynamics-aware 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, HMBD23, JLOW24] in chaotic systems when using recurrent architectures have been proposed recently. The failure of generalization notions based on mean-squared error have also been noted in [SWP + 24, JLOW24] 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 to 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 a future work.Ergodic theory and shadowing. Hyperbolic dynamics and ergodic theory (see e.g. the textbook [KKH95] ) lay the foundation for understanding the long-time/statistical physical behavior [You02] of complex systems. The scientific computing community has leveraged ergodic theory and shadowing [Ano67, Bow75, Pil06] for rigorous computations that use high-fidelity numerical simulations of chaotic systems [GL24, Wan13, Ni21] and to analyze the correctness of numerical simulations [HYG87, CW21, Lia17, Sau05, 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 + 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 [LLSK + 22]. 
conclusion: 
 Our dynamics-aware generalization (Theorem 1) 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-practice 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 not imply statistical accuracy and preservation of dynamical invariants and properties such as Lyapunov exponents, which is crucial for trustworthy ML for science. Thus, we reinforce the need to move toward a context-aware theory of generalization, organically unifying complex dynamics with learning theory.Limitations: A limitation of the theoretical results about C 1 generalization is that we need to assume the typicality of shadowing. Empirically, the Jacobian can 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 a future work. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?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 or 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). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
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 is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness 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 a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an 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 a 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 as 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 or 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 to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe 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)? 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 a 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 a best faith effort. 12. Licenses for existing assets Question: 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? Answer: [Yes] Justification: Citation and versions are properly cited throughout the paper. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a 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 of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. 
 
 Appendix: Our results in section 3 make the assumption that the map F is uniformly hyperbolic. Roughly speaking, this means the uniform (on the attractor Λ ⊆ M ) expansion and contraction of infinitesimal perturbations under the differential dF. More precisely, in a uniformly hyperbolic system, there exists a decomposition of the tangent bundle, T M = E u ⊕ E s , into an unstable (E u ) and stable (E s ) subbundle such that these subbundles are dF -invariant, and moreover, for any vector field v ∈ E s , ∥dF t v∥ ≤ C λ t ∥v∥, t ∈ Z + . Similarly, an infinitesimal perturbation along vector field v ∈ E u shows a uniform exponential decay in norm backward in time, ∥dF t v∥ ≤ C λ |t| ∥v∥, t ∈ Z -. Here the constants C, λ are uniform over Λ. There are several relaxations of the uniformity in the splitting to obtain nonuniform and partial hyperbolicity; there are also several known examples of 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 [You02, You17] . 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, EHL04, CW22], where a 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.Here, we complete a proof sketch that we begin in section 3. An element v ∈ T n,nn M with v = ⊕ t≤n v t can be identified as t → v t ∈ T x nn t M. We define a functionF nn : T n,nn M → T n,nn M as F nn (v) t+1 = x nn t+1 -F (x nn t + v t ). If v is a fixed point of F nn , that is, F nn (v) = v, then, {x nn t + v t } t≤n is an orbit of F. Writing F nn (v) = F nn (0) + dF nn (0)v + N nn (v), where N is the nonlinear part of F nn , a fixed point v satisfies, (Id -dF nn (0))v = F nn (0) + N nn (v). To show that v exists, we follow Theorem 18.1.3 of [KKH95] and prove that T nn (w) = (Id-dF nn (0)) -1 (N nn (w)+ F nn (0)) is a contraction on a compact subset of T n,nn M . We have ∥T nn (v) -T nn (w)∥ ≤ ∥(Id -dF nn (0)) -1 ∥∥N nn (v) -N nn (w)∥ ≤ C(ϵ 0 , ϵ 1 )∥N nn (v) -N nn (w)∥ with probability > 1 -δ. To see this, note that dF nn (0)(v) = dF v, and F is a hyperbolic map, by assumption, and hence Lemma 18.1.4 of [KKH95] applies. Next, assuming that sup x∈M ∥d 2 F (x)∥ is bounded (i.e., dF is Lipschitz), we obtain that dN nn is Lipschitz, and hence, ∥T nn (v) -T nn (w)∥ ≤ ∥(Id -dF nn (0)) -1 ∥∥N nn (v) -N nn (w)∥ ≤ C(ϵ 0 , ϵ 1 )Kδ 0 ∥v -w∥ with probability > 1 -δ. Here, max t {∥v t ∥, ∥w t ∥} ≤ δ 0 , which is chosen independent of n, ϵ 0 , ϵ 1 and δ, and K is the Lipschitz constant of dN nn . Thus, T nn is a contraction on a δ 0 ball around 0 ∈ T n,nn M when C(ϵ)Kδ 0 < 1 -ϵ 2 , for some ϵ 2 > 0. Since ∥T nn (0)∥ ≤ C(ϵ 0 , ϵ 1 )ϵ 0 , and for any v in a δ 0 ball around in 0 in T n,nn M, ∥T nn (v)∥ ≤ C(ϵ 0 , ϵ 1 )ϵ 0 + (1 -ϵ 2 )δ 0 , when ϵ 0 , ϵ 1 and δ 0 are such that C(ϵ 0 , ϵ 1 )ϵ 0 + (1 -ϵ 2 )δ 0 < δ 0 , we have that T nn maps a δ 0 ball around in 0 in T n,nn M to itself. Thus, the unique fixed point (from contraction mapping theorem on the δ 0 ball) lies within the ball.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 torchdiffeq3 library [CRBD18] , with an absolute and relative error tolerance of 10 -8 . 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 test data. The time step size of the simulation is included in Table 2 .1D: Tent maps As our first examples, we choose three types of tent maps defined in [CW21] , 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 s = 0.2 and s = 0.8 in the plucked map following [CW21] .F (x; s) = 2 1+s x, x < 1 + s, 2 1-s (2 -x),x ≥ 1 + s;(5)EQUATIONwheref s (x) = min( 2x 1 -s , 2 - 2(1 -x) 1 + s ), x < 1, o s (x) = fs(2x) 2 , x < 0.5, 2 -fs(2-2x) 2 , x ≥ 0.5, λ s,n (x) = o s (2 n x -⌊2 n x⌋) 2 n + 2 ⌊2 n x⌋ 2 n .2D: Baker's map We use the following perturbation of the classical Baker's map from [CW22] . 3D: Lorenz '63 We conduct extensive numerical experiments in this paper with the Lorenz '63 system [Lor63] , with σ = 10, β = 8/3, and ρ = 28, which we describe in section 2.3D: Rössler We use the parameter setting a = 0.2, b = 0.2, and c = 5.7 in the Rössler system [Rös76] below, which is in the chaotic regime.dφ t dt ([x, y, z] ⊤ ) = -y-z x+ay b+z(x-c).4D: Hyperchaos Another test case we consider is the hyperchaotic system below in 4 dimensions from [Zha17] , where the parameter values are set as a = 16, b = 40, c = 20, d = 8 for the system to show chaotic behavior.dφ t dt ([x, y, z, w] ⊤ ) =    ax + dz -yz xz -by c(x -z) + xy c(y -w) + xz.127D: Kuramoto-Sivashinsky To test Neural ODE's performance in learning high dimensional chaotic system, we generate the modified Kuramoto-Sivashinsky (KS) system's solution defined below with a second order finite difference scheme used in [BW14] .∂u ∂t = -(u + c) ∂u ∂x - ∂ 2 u ∂x 2 - ∂ 4 u ∂x 4 x ∈ [0, L], t ∈ [0, ∞) u(0, t) = u(L, t) = 0 ∂u ∂x x=0 = ∂u ∂x x=L = 0 u(x, 0) = u 0 (x)We try two neural network architectures for learning the vector field of the systems: a simple Multi-Layer Perceptron (MLP) model and a ResNet [HZRS16] . For the MLP model, we use GELU4 as the activation function, and ReLU 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 [LWCD20] for results in section 4. When v F (x) is a true vector field and v h (x) is a learned vector field by a neural network, h, given solution x, relative error can be defined as:relative error(x) = ∥v F (x) -v h (x)∥ ∥v F (x)∥ (8)Using a grid search, hyperparameter values that yield the lowest relative error in the vector field as per (8) were chosen. Hyperparameter search results are shown in Tables 7, 8, and 9. Hyperparameter values that are different for each system are in the Table 2 . We use the AdamW [LH17] 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 [LWCD20] . In the FNO network, we fix the number of modes to 4, with a batch size of 200 with 4 Fourier Neural layers. Further details on latent SDE are discussed in section C.8.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.In this section, we present results that are described in section 2 for testing the statistical accuracy and generalization of learned models of the Lorenz '63 system. We also present additional results for the tent maps and the KS equation, described in the previous section. Figures 4 and 5 illustrate the comparison of relative error trends, as defined in Eq. ( 8), for vector fields learned by a Neural ODE (ResNet) trained with mean squared loss and Jacobian-matching loss (as defined in Eq. ( 3)), for two distinct dynamical systems.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) as defined in ( 1) and (3) respectively. The vector field is evaluated on a random true orbit.In Figure 6 , we show learned models of representative tent map perturbations from the preceding section. We observe that the Jacobian-matching loss leads to reasonably accurate representations of the maps. Yet, the LE computed along learned orbits differ significantly at s = 0.8, as shown in Table 5 . For various other perturbed tent maps, we show the computed LEs alongside the true LEs in Table 5 ; for most of the maps, the Jacobian-matching leads to accurate LE predictions, consistent with the results of Theorem 1.Here we show identical results to Figure 1 but with ResNet-based Neural ODE models as opposed to MLPs used in Figure 1 . We find that Jacobian-matching training leads to superior performance in Figure 5 : Comparison of relative errors in the vector fields of the Rössler system produced by a Neural ODE with a ResNet, trained using mean squared loss (Left), and Jacobian-matching loss (Right) as defined in ( 1) and (3) respectively. The vector field is evaluated on a random true orbit. Details on hyperparameter setting and other statistical results can be found in Table 2 and 5 .terms of distributional match with µ, when compared to MLPs. Furthermore, our overall inference about MSE models leading to atypical orbits still holds, agnostic to modeling and architectural choices. After spatial discretization and following the scheme of [BW14] , we obtain a 127-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 x and time (T). Since the original model and the 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 in the MSE model grow strikingly quickly.Figure 8 : Solution plot of Kuramoto-Sivashinksy 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 (1) (center column), solution of the Neural ODE trained with Jacobian-matching loss defined in (3) (right)We report Lyapunov exponents learned with Neural ODE trained with mean squared loss and Jacobianmatching loss in Table 5 . To compute Lyapunov exponents, we use a QR algorithm of Ginelli et 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 are inaccurate in the rest of them.In this section, we report experiments on the Lorenz '63 system with an unrolled loss function:EQUATIONWith k = 10 timesteps of unrolling time, we see that the attractor is reproduced well, as shown in Figure 9 . However, atypical orbits are still produced for random initializations (Figure 10 ).We experiment using two types of v nn , MLP and Resnet, and varying sequence lengths, k. As shown in Table 3 , we observe that the learned negative Lyapunov exponent plateaued for k ≥ 40, remaining between -9 and -10 (the true value being ∼ -14.5), with no further improvement. Also, as k 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 be 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 to the Lyapunov time, v t are in tangent spaces along the orbit. This yields the recursive relationship, v t ≈ dF (x t-1 )v t-1 , when O(∥v t ∥ 2 ) is negligible. 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. 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 statistics and the true statistics of Lorenz '63.Here we present the experimental details and results of learning the Lorenz '63 system with latent SDEs [LWCD20] that are described in Section 4.As training data, we use 1024 time series of the interval [0, 6] and timestep size 0.01. With a smaller sample size, we empirically observe a large distributional mismatch (see Figure 11 ). We also observe that increasing or 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-Kutta scheme (RK4). We follow the same architecture5 as used in [LWCD20] , with a GRU encoder, a linear decoder, and the drift and diffusion functions in the prior and posterior processes modelled by MLPs.We present the learned Lyapunov exponents in Table 6 , and a comparison of the learned and true empirical measures in Figure 11 . When learning a deterministic system with a latent SDE, the diffusion coefficient of the learned system is small but not exactly 0, and this results in a slight difference in the Lyapunov exponents ([GHL20]) compared to only using the learned drift term. The latent SDE model was also tested on the stochastic Lorenz attractor [CZH21] , where it reproduced the 'bimodal' distribution of the trajectories. 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 in 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 on 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 a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.We use this notation rather than the customary dx(t)/dt = v(x(t)), x(0) = x0, since in dynamical systems, we are interested in all possible orbits, as opposed to a particular orbit/path, for which defining the flow, φ t , is necessary.https://github.com/rtqichen/torchdiffeqhttps://pytorch.org/docs/stable/generated/torch.nn.GELU.htmlhttps://github.com/google-research/torchsde/blob/master/examples/latent_sde_ lorenz.py 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

2. 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.

[1] Learning Dissipative Dynamics in Chaotic Systems (https://arxiv.org/pdf/2106.06898)

[2] DySLIM: Dynamics Stable Learning by Invariant Measure for Chaotic Systems (https://arxiv.org/abs/2402.04467)

[3] Training neural operators to preserve invariant measures of chaotic attractors (https://arxiv.org/abs/2306.01187)

[4] On the difficulty of learning chaotic dynamics with RNNs (https://arxiv.org/pdf/2110.07238)
questions: 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.

Q2: 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?

Q3: Line 183: It's not clear why the authors state that implementing Wasserstein distance is difficult when the map $h$ is chaotic.

Q4: As Lyapunov exponents are calculated using local dynamics, could you show a comparison using some long-term evaluation metrics?
limitations: limitations: N/A
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

Another 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.
questions: 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: limitations: The paper does not have a section on limitations. I think a discussion of limitations would add to the value of the paper.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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].
- 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
- 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.
- 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
- 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.
- 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
	
I 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.

Miscellaneous comments:
- 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
- Hyperlink of Table 3 jumps to section 3
- typo l. 192: “shdaowing”

[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)

[2] [Hess et al. (ICML, 2023), Generalized teacher forcing for learning chaotic dynamics.](https://proceedings.mlr.press/v202/hess23a.html)

[3] [Göring et al. (ICML, 2024), Out-of-domain generalization in dynamical systems reconstruction.](https://arxiv.org/abs/2402.18377)
questions: questions: - how exactly are the empirical measures $\mu$ and $\mu_{NN}$ and their $W_1$ distance computed (for Table 1 e.g.)?
- 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$? 
- How can one draw practical consequences from the theoretical findings?
limitations: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: I believe that the most important weakness of the manuscript relies on the limitations of the results.
The 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.
questions: 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: 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.
I 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.

Comments: title: Thank you for the comments and discussion on learning chaotic dynamics
Thank you and request for response!
Thank you!
Response to your rebuttal
Thank you!
Re: Jacobian and LEs
Re: Jacobian and LEs
Request to read rebuttals
Jacobian and LEs
Re: Rebuttal
Request for response
Thank you for the feedback
comment: Dear Reviewers and ACs,

We have truly grateful for all the time you have spent in giving us a thorough review. The subsequent discussions (including the references) have been very fruitful and led to substantial improvements in the revision. We are excited about both directions our work contributes to and is influenced by: practical methods for learning physical measures (ensemble/long-term behavior) of chaotic dynamics using supervised learning and introducing dynamical systems/ergodic theory to analyze and understand generalization.  Overall, the intersection of nonlinear dynamics with machine learning theory and practice is growing and we are very happy for the opportunity to be able to improve our work at this intersection through these discussions with you.
Dear Reviewers,

Thank you for reading our rebuttals, seeking further clarification and re-evaluating our paper! 
We are extremely grateful for all the time you have spent on improving our work.
If you have not done so already, could you please read our rebuttal, and re-assess our work in light of the clarifications in the rebuttal?
Sorry for the multiple reminders -- we are almost at the end of the discussion period.

Thank you again!
We really appreciate your time and for re-assessing our work based on our rebuttal! 
Thank you!
Thank you for your response! Your response regarding the Sobolev norm was valuable, and your further clarifications cleared up my concerns. I appreciate your discussion regarding the related works, which helped me better assess the status of your work. I will increase my score to accept.
Thank you very much -- we especially appreciate your questions!
The purpose of this paper, as you indeed correctly point out, is dynamics-aware generalization bounds. We show that learning the Jacobian leads to statistical accuracy and *explain why* by deriving a generalization bound via shadowing theory. 
Computing the Jacobian is actually not impractical compared to more sophisticated training approaches using recurrent architectures or generative modeling of the SRB measure (physical measure in the paper). 
We appreciate your time and effort! Please let us know if we can provide any further clarification!
Thank you for reading through our comment and for your response! We are so grateful to you for responding promptly and thoughtfully!
You are absolutely right. Since LEs are statistical quantities, and for a singular hyperbolic system like the Lorenz system, the Central limit theorem holds, we do see a 1/sqrt{N} convergence of ergodic averages and LEs computed along trajectories of N length. The "true" label does refer to long averaging times conducted with the ground truth system -- i.e., the true vector field. It is indeed a fair comparison to maintain an equally long time averaging but compute the LEs with the Neural ODE as a vector field instead. These are what are referred to as "JAC" or "MSE" values -- the LEs computed with these learned models. A 5% difference in LEs can happen even in between runs with different initial conditions and this should, for instance, not be considered as an indication of lack of ergodicity. Overall, we want to clarify that the LEs are expectations of random variables and therefore cannot be estimated like deterministic quantities. It would not be meaningful to report them (in a chaotic system) with many digits of precision, e.g.
Please let us know if any further clarification is needed!
Thanks for the clarification - makes sense that the runtime complexity increases similar to PINN training! 

Regarding the LEs; First of all, thanks for providing the code. I think this needs more careful evaluation. If you call something "true" LEs, you should indeed compute the right quantity - for the true spectrum, the statistical error or noise should be negligible. I.e. you need to let the trajectories converge into their limit set and draw a long enough trajectory to pin down the correct spectrum (you may plot time steps/integration time vs error in the spectrum to verify you draw long enough orbits to converge). But if the reported "true" LE is indeed ~0.905 instead of ~0.85, then this changes the interpretation of Table 2 (is it just that the models need to be evolved for longer or are they then indeed underestimating the max. LE even with Jac. loss?).

To be clear, I really appreciate the effort the authors put into the paper and rebuttal and I think the paper is of great value, but to me it seems it could really benefit from another revision and submission at a later point. 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).
Thank you for the clarifications.
From the other reviews, I see there are some concerns regarding practicality of obtaining the Jacobian. However, I think the explanation of generalization also has value. I maintain my accept recommendation.
Dear Reviewers,

Since the discussion period is set to close soon, we request you to please read our rebuttals and ask for clarification. If we have answered all your questions satisfactorily, we request you to please re-evaluate our paper.
Thank you very much for reading through the rebuttal! We really appreciate your interest and excellent questions!

> Does runtime/complexity still only scale with problem dimensionality squared?

This is a very keen observation -- thank you. We meant that computing the Jacobian (using Finite Difference or AD) has quadratic complexity, but you are absolutely right that training time increases also because of differentiating the model at the training points. During training the Jacobian (with respect to the input state) is differentiated with respect to the parameters of the neural network. The complexity with respect to the dimension of the problem still remains quadratic but the training complexity increases just like a physics-informed neural training; the factor of increase therefore also depends on the complexity of the network. Inference time is the same as a vanilla neural ODE however since only the training is modified. Practically we found that KS system (128-dimensional system) was still training in about 10 hours on a single RTX30xx GPU. 

> estimation of the Jacobians from time series data, where no ground truth knowledge is available

This is exactly the scenario we had in mind for the rebuttal. We can estimate the Jacobian from timeseries data using finite difference as long as the data sample the system frequently (there are classical results known for the impossibility of system identification for longer observation frequency and very noisy observations of chaotic systems). As you correctly point out, there is an error made by finite difference (which is on order of epsilon^2, where epsilon is the step size norm). However, even though our theoretical results only use exact Jacobians, the results with unrolling dynamics, e.g., where only indirect Jacobian information is available, suggest that even inexact Jacobians can help us learn the dynamical invariants better. 
For instance, even Lyapunov vectors and exponents can be estimated with high accuracy using just finite difference estimates of the Jacobian (this is easy to try by replacing the AD with inexact Jacobians in the code snippet below).

But as you suggest, we should test this on real world high-dimensional examples (even ones where the attractor is only known through delay embedding), which we defer to a future work. Our present paper emphasizes on understanding why Jacobian information leads to statistical accuracy. Our main contribution is using shadowing theory to provide an explanation and suggest that there is hope of learning chaotic systems with elementary supervised learning methods, and that these can still learn the invariant measure.

About the Lyapunov exponents, small differences (statistical error/noise) may arise depending on the integration time (len(traj_gpu) in the code below), time integration method (RK4 in the code below), and subsequently how we estimate Jacobians. Here is the code snippet we wrote for computing 3 LEs:

```
def lyap_exps(dyn_sys_info, traj, iters):
    model, dim, time_step = dyn_sys_info
    LE = torch.zeros(dim).to(device)
    traj_gpu = traj.to(device)
    f = lambda x: rk4(x, model, time_step)
    Jac = torch.vmap(torch.func.jacrev(f))(traj_gpu)
    Q = torch.rand(dim,dim).to(device)
    eye_cuda = torch.eye(dim).to(device)
    for i in range(iters):
        if i > 0 and i % 1000 == 0:
            print("Iteration: ", i, ", LE[0]: ", LE[0].detach().cpu().numpy()/i/time_step)
        Q = torch.matmul(Jac[i], Q)
        Q, R = torch.linalg.qr(Q)
        LE += torch.log(abs(torch.diag(R)))
    return LE/iters/time_step
```

Thank you very much for your questions and please feel free to let us know if anything is unclear!

Thank you again for your time and help!
> Related work on invariant measures

Happy to hear that the authors add the references to the related work and thanks for the additional clarification and connection to their own work!

> Practicality of Jacobian loss and implications of theory

Indeed, the complexity of computing the *model Jacobian* using AD/Finite differences is one part, however, did the authors also consider the added complexity by including this Jacobian in the loss function, which means one has to a backward-over-forward AD problem (-> second order differentiation)? Does runtime/complexity still only scale with problem dimensionality squared?

My main concern with the practicality is the estimation of the *Jacobians from time series data*, where **no** ground truth knowledge is available (I appreciate the geophysics examples and I see that in these cases the Jacobian matching loss might improve things). I really miss any experiments on trying to make this loss more applicable in real-world scenarios (e.g. estimating Jacobians from time delay embeddings of real-world data -> do crude approximations already help learning dynamical invariants?).

> Unrolling dynamics

Thanks for conducting the additional experiments! That zero/negative LEs are still underesimated when unrolling dynamics is indeed a problem I encountered in practice, too. Hence it is great to see that knowledge of the Jacobians can eliminate this problem. However, I checked the paper again and found that the authors report a ground-truth max. LE for the Lorenz63 system with standard settings of ~0.86 (cf. Table 2). However, to my knowledge the Lyapunov spectrum is given by ~[0.905, 0, -14.57]. I checked this with minimal code using Julia (see below). Do the authors mind to share how they estimated the ground-truth Lyapunov spectrum or whether they used specific literature values?
```
using DynamicalSystems, ChaosTools

# Lorenz63
ds = Systems.lorenz(ρ=28.0, σ=10.0, β=8 / 3)

# Lyapunov spectrum
lyap = lyapunovspectrum(ds, 100_000, Ttr=1000) # ~ [0.905, 6e-6, -14.57]
```
> W1 distance

Thanks for the clarification!

> Edits

I see, thanks!
Dear Reviewers,

We request you to please find time to read our rebuttal and help us identify any concerns that we have not yet addressed.
We really appreciate your time and several insightful comments, which have considerably enhanced our paper!
Dear Reviewers,

We are extremely grateful to you for spending the time to understand our paper and raise many interesting questions. Due to lack of space, we have not quoted from your review in our responses, but we have answered all your comments to the best of our knowledge under easily identifiable headings. We outline the revisions/edits that we have already made to our paper, and as soon as the option becomes available, we will upload the revised version. For new experiments that we have conducted, we have shared **anonymous links to the figures** since we cannot yet share the revised paper. We request you to visit these links to see the plots for the new experiments on unrolled loss functions. We have indicated in our response wherever there are similarities in the questions/comments by more than one Reviewer. We request you to therefore read the *relevant sections of the rebuttal to other reviewers* as well. 
We thank you once again for your time and valuable feedback that has improved our paper substantially!