Title: LEARNING LATENT DYNAMICS FOR PARTIALLY-OBSERVED CHAOTIC SYSTEMS

Abstract: This paper addresses the data-driven identification of latent Ordinary Differential Equation (ODE) representations of partially-observed dynamical systems, i.e. dynamical systems whose some components are never observed, with an emphasis on forecasting applications and long-term asymptotic patterns. Whereas stateof-the-art data-driven approaches rely on an explicit mapping of the observation series onto a latent space, the proposed approach only relies on the definition of an augmented space, higher-dimensional than the manifold spanned by the observed variables, where the dynamics of the observations can be fully described by an ODE. From a numerical point of view, the proposed approach exploits a neural-network ODE representation and the associated variational minimization scheme. Numerical experiments support the relevance of the proposed framework w.r.t. state-of-the-art approaches, including neural ODE schemes, both in terms of short-term forecasting errors and long-term behaviour. We further discuss how the proposed framework relates to Koopman operator theory and Takens' embedding theorem.

1. INTRODUCTION:

Learning the underlying dynamical representation of some observed variables x t ∈ R n (where t ∈ {t 0 , ..., T } is the temporal sampling time and n the dimension of the observations) is a key challenge in various scientific fields, including control theory, geoscience, fluid dynamics, economics; for applications ranging from system identification to forecasting and assimilation issues (Lai & Wei, 1982; Abarbanel & Lall, 1996; Jeong & Hussain, 1995; Koopmans, 1949) .

For fully-observed systems, i.e. when the observed variables x t relate to some underlying deterministic states s t according to a diffeomorphic mapping possibly corrupted by noise, recent advances (Brunton et al., 2016b; Fablet et al., 2018; Nguyen et al., 2019) have shown that one can identify the governing equations of the dynamics of s t from a representative dataset of observations {x ti } i . Unfortunately, when the observed variables x t only relate to some but not all the components of underlying states s t , these approaches can not apply since no ODE or, more generally, no one-toone mapping defined in the observation space can represent the time evolution of the observations. In this context, Takens's theorem states the conditions under which a delay embedding, formed by lagged versions of the observed variables, guarantees the existence of governing equations in the embedded space (Takens, 1981 ).

Takens's theorem has motivated a rich literature of machine learning schemes to identify dynamical representations of partially-observed systems using a delay embedding. This comprises both nonparametric schemes based on nearest-neighbors or analogs (Abarbanel, 1996a) as well as parametric schemes which include polynomial representations (Paduart et al., 2010) , neural network models (Wan, 1993) and Support Vector Regression (SVR) models (Mattera & Haykin, 1999; Müller et al., 1999) . For all these approaches, the identification of the appropriate delay embedding is a critical issue (Abarbanel, 1996b; c) .

From a neural network and machine learning perspective, the inference of a latent space, within a State Space Model (SSM) framework, for dynamical systems has motivated a broad literature especially for time series forecasting (Ghahramani & Roweis, 1999; Wang et al., 2006; Mirowski & LeCun, 2009; He et al., 2015; Krishnan et al., 2016) . Most of those techniques were introduced in the context of reduced order modeling (ROM) to infer low-dimensional manifolds, where the dynamics of the observations can be represented. When considering partially-observed systems, these approaches state this issue as the inference of a (non-linear) projection of an input sequence in a latent space where the observations can be modeled. This projection is usually computed in a probabilistic framework using Bayesian filtering techniques. However, recovering the attractor's dynamics using iterative predictions is still an issue for such models since the explicit modeling of latent space as a delay embedding of the observations may limit the expressiveness of the latent states, especially when considering chaotic dynamics.

In this work, we show that we do not need to rely explicitly neither on a delay embedding nor on the learning of an inference model to map the observation series to a latent space. We only assume that we may define an augmented space, higher-dimensional than the manifold spanned by the observed variables, where the dynamics of the observations can be fully described by an ODE. Using neural-network representations for the parametrization of the dynamical model, it amounts to jointly learning the governing ODE and reconstructing the augmented latent states for a given observation dataset using a variational framework. We report experiments on linear and chaotic dynamics, which illustrate the relevance of the proposed framework compared to state-of-the-art approaches. We then further discuss the key features of this framework with respect to state-of-theart dynamical systems identification tools such as Koopman operator theory (Koopman, 1931) .

2. BACKGROUND AND RELATED WORK:

This section introduces the learning of dynamical representations for partially-observed systems and links this problem to recent advances in machine learning.

Let us consider an unobserved state variable s governed by an autonomous system of m differential equations ṡt = f (s t ). Let us also assume that this system generates a flow Φ ti (s t0 ) = ti t0 f (z w )dw ∈ R m with trajectories that are asymptotic to a limit-cycle L of dimension d contained in R m . We further assume that we are provided with a measurement function H that maps our state variable s to our observations x t = H(s t ) ∈ R n .

When considering the data-driven identification of a dynamical mapping that governs some observation data, we first need to evaluate whether the dynamics in the observation space can be described using a smooth1 ODE. Another way to tackle this question is to find the conditions under which the deterministic properties of the unobserved limit-cycle L are preserved in the observation space in R n such that one can reliably perform forecasts in the observation space. The general condition under which a mapping H preserves the topological properties of the initial limit-cycle involves a differential structure. Assuming that L is a smooth compact differential manifold, the topological properties of L are preserved through a mapping H in R n if H is one-to-one and is an immersion of L in R n . Under these conditions our observation mapping is called an embedding (Sauer et al., 1991) .

The simplest example of an embedding involves an identity observation operator H. With this embedding, we have direct access to the state variable s which is governed by a deterministic ODE. This particular case has been widely studied in the literature. Parametric representations have been for decades the most popular models thanks to their simplicity and interpretability (Paduart et al., 2010) , (Brunton et al., 2016b) . Recently, these approaches have been enriched by neural network and deep learning schemes (Wiewel et al., 2018) , (Raissi et al., 2018) . In particular, the link between residual networks (Chen et al., 2018; Ouala et al., 2019) and numerical integration schemes have opened new research avenues for learning extremely accurate dynamical models even from irregularly-sampled training data. These schemes show greater interpretability and forecasting performance for the datadriven representation of systems governed by an ODE, compared with other state-of-the-art neural networks schemes, including Recurrent Neural Networks (RNN) such as LSTM (Long-Short-Term Memory).

However, for a wide range of real-world systems, we are never provided with an observation operator that forms an embedding of the unobserved dynamical system. In such situations, we do not have any guarantee on the existence of a smooth ODE that governs the temporal evolution of our observations. From this point of view, the question of finding an appropriate dynamical representa-tion of some observed data may not be this straightforward. The fact that our data may come from some unobserved governing equation may restrict the use of the above-mentioned state-of-the-art algorithms. The main difficulty lies in the ability to map observation series to a latent space that provides at least a one-to-one mapping between two successive states. From a geometrical point of view, the time delay theorem (Takens, 1981) provides a way to build a latent space that preserves the topological properties of the true (unobserved) dynamics limit-cycle. A generalization of this theorem (Sauer et al., 1991) shows that one can reconstruct topologically similar limit-cycles using any appropriate smooth composition map of the observations. The derivation of a dynamical system from such representations however encounters large disparities since no explicit relationships between the defined phase space and an ODE formulation have been clearly made. Classical state-of-the-art techniques such as polynomial representations (Brunton et al., 2016b) and K-Nearest Neighbors (KNN) (Lguensat et al., 2017) algorithms were proposed but they often fail to achieve both accurate short-term forecasting performance and long-term topologically similar reconstructed limit-cycles (see experiments for an illustration).

We may also point out that the limitation of ODE-based representation in deep learning architecture has also been pointed out recently in (Dupont et al., 2019; Zhang et al., 2019) for classification issues. As ODE-derived trajectories do not intersect, it may limit the ability of neural ODE representations to reach relevant classification performance in a given feature space. To address this issue, (Dupont et al., 2019) and (Zhang et al., 2019) propose to consider an augmented state, simply by augmenting the observed state by a number of zeros to create a high-dimensional space in which an ODE representation can be identified. Such a strategy cannot apply to time series modeling as successive augmented states cannot be forced to zero for some dimensions.

Advances in the inference of latent spaces in state space models was introduced essentially, from a dynamical systems perspective, to retrieve low-dimensional manifolds, where the dynamics of the system evolve. When applied to partially-observed systems, the latent variables are typically inferred from a sequence of observations through a parametric modeling of the posterior distribution as in (He et al., 2015; Krishnan et al., 2016; Chen et al., 2018) or through marginalization with model constraints as in (Ghahramani & Roweis, 1999; Wang et al., 2006) However, such models often fail in accounting for long term patterns (as shown in the experiments). This is due to the fact that the latent space is constrained to be a non-linear projection of a sequence of observations, which limits the expressiveness of the dynamical model. Interestingly, (Mirowski & LeCun, 2009) does not involve the learning of an inference model as the reconstruction of the latent states is solved as gradient-based minimization of the dynamical prior w.r.t. an observation series. However, the dynamical prior relies on an explicit delay embedding as the dynamics of the latent state depend both on the previous latent state and on a delay embedding of the observations.

In this work, we address the identification of a latent embedding, associated with an ODE representation, for partially-observed systems. The core idea of this work is to infer an augmented latent space, governed by an ODE, which fully explain the observed time series and their dynamics. In contrast to previous work (He et al., 2015; Krishnan et al., 2016; Chen et al., 2018; Ghahramani & Roweis, 1999; Wang et al., 2006) , we do not exploit either a delay embedding or an explicit modeling of the inference model (i.e., the reconstruction of the latent states given the observed time series). As such, our scheme only involves the selection of the class of ODEs of interest. The expected benefits are as follows: (i) our model ensures the existance of a latent embedding associated with an ODE, which may not be guaranteed when considering a parametric inference model and/or a delay embedding, (ii) our model reduces the complexity of the overall scheme to the complexity of the ODE representation, (iii) our model guarantees the consistency of the reconstructed latent states w.r.t. the learnt ODE. Our experiments suggest that this may be of key importance to jointly address short-term forecasting performance and long-term chaotic patterns.

3. LEARNING LATENT REPRESENTATIONS OF PARTIALLY-OBSERVED:



DYNAMICS:

Augmented latent dynamics: Let us consider a continuous m-dimensional dynamical system s t governed by an autonomous ODE ṡt = f (s t ) with Φ t the corresponding flow Φ t (s t0 ) = t t0 f (s w )dw. We assume the state s t to be never fully-observed. Formally, we can define an obser-vation function H : R m -→ R n such that the observation x t follows x t = H(s t ). An example of interest comprises the situation when one observes only the n first variables of s t .

Our goal is to derive an ODE representation of x t ∈ R n . However, the key question arising here is the extent to which the dynamics expressed in the observations space, reflect the true underlying dynamics in R m , and consequently, the conditions on H under which the predictable deterministic dynamical behavior of the hidden states is still predictable in the observations space. To illustrate this issue, we may consider a linear dynamical system in the complex domain governed by the following linear ODE: ṡt = αs t s t0 = s 0 (1) with s ∈ C a state variable and α ∈ C a complex imaginary number. The solution of this problem is

s t = Ke αt (2)

with K a constant depending on s 0 . Let us assume now that we are only provided with the real part as direct measurements of the unobserved state i.e. H(.) = Real(.) :

x t = Real(s t ).

Proposition 1 : The flow of an ODE cannot represent the time evolution of x t .

The proof of the proposition is given in the appendix. The intuition exemplified for the 2dimensional system given by ( 1) is as follows. If we assume that we are only provided the real part as direct measurements x t ∈ R of the true states s t , no smooth autonomous ODE model in the scalar observation space can describe the trajectories of the observations as the mapping between two observations is not one-to-one. For example, assuming that s t0 and s t1 correspond to two states that have the same real part but distinct imaginary parts, the associated observed states are equal x t0 = x t1 . However, the time evolution of the states s t0 and s t1 differ if they have different imaginary parts. The observed states x t0+δ and x t1+δ are then no longer equal after any time increment δ. As a consequence, a given observation may have more than one future state and this behavior cannot be represented by a smooth ODE in the observation space. Hence, the application of an ODE mapping such as (Chen et al., 2018) and (Fablet et al., 2018) for such observations will lead to poor forecasting performance. From a naive Neural networks point of view, fitting such a model will most likely force the forecasting into an equilibrium point since we may iteratively match the same inputs to different output predictions. This simple example emphasizes the need for considering augmented latent states.

In this context, Takens's theorem guarantees the existence of an augmented space, in which a oneto-one mapping exists between successive time steps (Takens, 1981) . More precisely, for a given observation operator H of a deterministic underlying dynamical system that governs s t , one can define a delay embedding from an observation sequence, i.e. {x t0-aδ , x t0 , . . . , x t0+bδ } with δ the time delay and a, b two natural numbers. An ODE representation in a latent space can than be formulated according to (Chen et al., 2018) by sampling the latent states from a parametrized posterior distribution given the delay embedding time series. Here, we aim to identify an augmented latent space without resorting to an explicit mapping from the observed series to the latent space, such that the latent dynamics are governed by a smooth ODE and can be mapped to the observations. Formally, let us introduce the following notations. Let u t ∈ R d E define a d E -dimensional augmented latent state formed by an observed component2 x t and an unobserved latent component referred to as z t .

EQUATION

Assuming that the augmented latent state u t is governed by an ODE leads to the the following state space model :

EQUATION

where the dynamical operator f θ belongs to a family of smooth operators (in order to guarantee uniqueness (Coddington & Levinson, 1955) ) parametrized by θ. G is a projection matrix that satisfies x t = G(u t ). It means that operator G only returns the observed component of the augmented latent state. Compared with previous work, the key difference lies in the definition of the augmented latent space, where only the observed component is explicitly stated as a function of observation x t . Similarly to neural ODE frameworks (Chen et al., 2018; Fablet et al., 2018) , we consider a neural-network representation with Lipschitz nonlinearities and finite weights for ODE operator f θ .

As detailed in the next sections, we address the identification of operator f θ from a dataset of observed state series {x 0 , . . . , x T } as well as the exploitation of the identified latent dynamics for the forecasting of the time evolution of the observed states, for instance unobserved future states {x T +1 , . . . , x T +N }. Both issues involve the joint inference of the unobserved component z t of the augmented latent states. This formulation generalizes to reduced-order modeling (ROM) (Champion et al., 2019) , where observation x t fully derives from a lower-dimensional variable. We let the reader refer to the appendix for the associated generalization of the proposed framework.

Learning scheme: Given an observation time series {x 0 , . . . , x T } we aim to identify the statespace model defined by ( 4), which amounts to learning the parameters θ of the dynamical operator f θ . However, as the component z t of the augmented state u t is never observed, this identification requires the joint optimization of the model parameters θ as well as of the hidden component z t .

Formally, this problem is stated as the following minimization of the forecasting error on observed variables:

EQUATION

with Φ θ,t the one-step-ahead diffeomorphic flow associated with operator f θ such that:

Φ θ,t (u t-1 ) = u t-1 + t t-1 f θ (u w )dw

In (5), the loss to be minimized involves the one-step-ahead forecasting error for the observed variable x t . The constraints state that the augmented latent state sequence {u 0 , . . . , u T } defined by (3) should be a solution of the ODE (4). Here, we numerically minimize the equivalent formulation:

min θ min {zt}t T t=1 x t -G (Φ θ,t (u t-1 )) 2 + λ u t -Φ θ,t (u t-1 ) 2 (6)

where u T t = [x T t , z T t ] and λ a weighting parameter. The term u t -Φ θ,t (u t-1 ) 2 may be regarded as a regularization term such that the inference of the unobserved component z t of the augmented state u t is not solved independently for each time step.

Using a neural-network parametrization for the ODE operator f θ , the corresponding forecasting flow Φ θ,t is also stated as a neural network based on a numerical integration scheme formulation (typically a 4 th -order Runge-Kutta scheme). This architecture, very much similar to a ResNet (He et al., 2015) , allows very accurate identification of ODE models (Fablet et al., 2018; Ouala et al., 2019) . Hence, for a given observation sequence {x 0 , . . . , x T }, we minimize (6) jointly w.r.t. θ and unobserved variables {z 0 , . . . , z T }.

Application to forecasting: We also apply the proposed framework to the forecasting of the observed states x t . Given a trained latent dynamical model (4), forecasting future states for x t relies on the forecasting of the entire augmented latent state u t . The latter amounts to determining an initial condition of the unobserved component z t and performing a numerical integration of the trained ODE (4).

Let us denote by x n t , t ∈ {t 0 , ..., T } a new series of observed states. We aim to forecast future states x n t , t ∈ {T + 1, ..., T + δT }. Following (6), we infer the unobserved component z T of the augmented state u n T at time T from the following minimization:

ẑn T = arg min z n T min {z n t } t<T T t=t0 x n t -G Φ θ,t u n t-1 ) 2 + λ u n t -Φ θ,t (u n t-1 ) 2 (7)

Here, we only minimize w.r.t. latent variables {z n t } given the trained forecasting operator Φ θ,t . This minimization relates to a variational assimilation issue with partially-observed states and known dynamical and observation operators (Lynch & Huang, 2010) . Similarly to the learning step, we benefit from the neural-network parameterization of operator Φ θ,t and from the associated automatic differentiation tool to compute the solution of the above minimization using a gradient descent.

We may consider different initialization strategies for this minimization problem. Besides a simple random initialization, we may benefit from the information gained on the manifold spanned by the unobserved components during the training stage. The basic idea comes to assume that the training dataset is likely to comprise state trajectories which are similar to the new ones. As the training step embeds the inference of the whole latent state sequence, we may pick as initialization for minimization (7) the inferred augmented latent state in the training dataset which leads to the observed state trajectory that is the most similar (in the sense of the L2 norm) to the new observed sequence x n t . The interest of this initialization scheme is two-fold: (i) speeding-up the convergence of minimization (7) as we expect to be closer to the minimum; (ii) considering an initial condition which is in the basin of attraction of the reconstructed limit-cycle. The latter may be critical as we cannot guarantee that the learnt model does not involve other limit-cycles than the ones truly revealed by the training dataset, which may lead to a convergence to a local and poorly relevant minimum.

4. NUMERICAL EXPERIMENTS:

In this section, we report numerical experiments to illustrate the key features of proposed framework. We consider three case-studies: a linear ODE case-study; a chaotic system, namely Lorenz-63 dynamics, and real upper ocean data.

Application to a linear ODE: In order to illustrate the key principles of the proposed framework, we consider the following linear ODE in the complex domain: ṡt = αs t s t0 = s 0 (8) with α = -0.1 -0.5j, j 2 = -1 and s 0 = 0.5. As α ∈ C with Real(α) < 0 and s 0 = 0, the solution of this ODE is an ellipse in the complex plane (Fig. 1 ).

As observations, we consider the real part of the underlying state, i.e. the observation function H : C -→ R is given by x t = Real(s t ). This is a typical example, where the mapping between two successive observations is not a one-to-one mapping since all the states that have the same real part lead to the same observation. As explained in section 3, one cannot identify an autonomous ODE model that will reproduce the dynamical behavior of the observations in the observations space.

We apply the proposed framework to this toy example. We consider a 2-dimensional augmented state

u t = [x t , z 1 t ].

As neural-network parametrization for operator f θ , we consider a neural network with a single linear fully-connected layer. We use an observation series of 10000 time steps as training data. As illustrated in Fig. 1 , given the same initial condition over the observable state, the inferred latent state dynamics, though different from the true ones, depicts a similar spiral pattern. This result is in agreement with the geometrical reconstruction techniques (Takens, 1981) of the latent dynamics up to a diffeomorphic mapping. Overall, our model learns a dynamical behavior similar to the true model represented by an elliptic transient and an equilibrium point limit-set. Furthermore, the projection of the augmented latent space and the true solution of Eq. ( 8) in the real axis illustrate the relevance of the proposed framework in forecasting the observations dynamics (mean square error < 1E -6).

Lorenz-63 dynamics: Lorenz-63 dynamical system is a 3-dimensional model that involves, under some specific parametrizations (Lorenz, 1963) , chaotic dynamics with a strange attractor. We simulate chaotic Lorenz-63 state sequences with the same model parameters as proposed in (Lorenz, 1963) using the LOSDA ODE solver (Hindmarsh, 1983) with an integration step of 0.01. We assume that only the first Lorenz-63 variable is observed x t = s t,1 . We apply the proposed framework to this experimental setting using a training sequence of 4000 time-steps.

For benchmarking purposes, we perform a quantitative comparison with state-of-the-art approaches using delay embedding representations (Takens, 1981) . The parameters of the delay embedding representation, namely the lag τ and the dimension d E of the augmented space were computed using state-of-the-art techniques. Specifically, the lag parameter was computed using both the mutual information and correlation techniques (Abarbanel, 1996b) , respectively denoted as τ M I and τ Corr . Regarding the dimension of the embedding representation, we used the Takens embedding condition d E = 2d + 1 with d the dimension of the hidden limit-cycle. The delay embedding dimension was also computed using the False Nearest Neighbors (FNN) method (Abarbanel, 1996c) . We also tested arbitrary parameters for the delay embedding dimension. Given the delay embedding representation, we tested two state-of-the-art data-driven representations of the dynamics. The Analog Forecasting technique (AF) which is based on the nearest neighbours algorithm (Lguensat et al., 2017) and the Sparse Regression (SR) method on a second order polynomial representation of the delay embedding states. Regarding deep learning models, we compare our method to a stacked Bidirectional LSTM (RNN) and to the Latent-ODE model as proposed in (Chen et al., 2018) . Finally, the proposed framework, referred to as Neural embedding for Dynamical Systems (NbedDyn) was tested for different dimensions of the augmented state space, namely from 3 to 6 (please refer to the apendix for details on the neural networks based architectures)3 . Fig. 2 illustrates the learning process for the latent space from the initialization to the last training epoch. We also report the analysis of short-term forecasting performance as well as the long-term asymptotic behavior characterized by the largest Lyapunov exponent of the benchmarked models in Tab 1. The proposed model leads to significant improvements in terms of short term forecasting performance with respect to the other approaches. Surprisingly, The Latent-ODE and RNN models lead to the poorest performance both in terms of forecasting error and asymptotic behavior. This is mainly due, in the Latent-ODE case, to the fact that the latent space is seen as a non linear projection of the observed variables through the optimization of the ELBO loss (Krishnan et al., 2016) . By contrast, our latent embedding formulation optimizes the latent states to forecast the observed variables which explicitly constrain the latent space to be an embedding of the true underlying dynamics. The RNN model in the other hand converges to a periodic solution (please refer to the appendix for forecasting figures) with still a poor short term forcasting performances. Overall, this results on deep learning models suggest that one should use such tools with care to guarantee satisfying the specifications of the underlying system.The SR model seems to lead to better short term forecast (using ad hoc parameters (τ = 6, d E = 3), however, it does not capture well the chaotic Under review as a conference paper at ICLR 2020 Figure 2 : Evolution of the learnt latent space: starting from a random initialization of the augmented states z i , the latent space is optimized according to the minimization of Eq. ( 6) to form a limit-cycle similar to the true Lorenz 63 attractor. We depict 3-dimensional projections of the learnt latent space for the proposed model with different embedding dimensions from

Model t 0 + h t 0 + 4h λ 1 AF τ M I =16 d E (F N N ) = 3 5.6E -3 1.3E -2 0.85 τ M I =16 d E (T akens) = 6 9.9E -3 2.4E -2 N aN τ Corr =27 d E (F N N ) = 3 8.9E -3 2.3E -2 12.35 τ Corr =27 d E (T akens) = 6 8.5E -3 1.9E -2 N aN τ = 6 d E = 3 8.0E -4 9.0E -4 0.87 τ = 10 d E = 3 2.1E -3 4.9E -3 0.60 SR τ M I =16 d E (F N N ) = 3 7.8E -2 2.5E -1 0.12 τ M I =16 d E (T akens) = 6 4.5E -2 1.7E -1 N aN τ Corr =27 d E (FNN) = 3 1.4E -1 4.6E -1 N aN τ Corr =27 d E (Takens) = 6 2.1E -1 8.4E -1 N aN τ = 6 d E = 3 7.6E -3 7.4E -3 N aN τ = 10 d E = 3 2.5E -2 5.7E -2 0.2535 Latent-ODE 6.9E -2 ± 2.9E -2 1.5E -1 ± 3E -2 N aN RNN 6.9E -2 ± 4.6E -2 1.5E -1 ± 1.1E -1 -6.79 ± 0.0 NbedDyn d E = 3 3.2E -4 ± 1.3E -4 1.7E -3 ± 7.5E -4 0.81 ± 0.09 d E = 4 1.3E -4 ± 5.2E -5 7.3E -4 ± 2.2E -4 0.82 ± 0.06 d E = 5 3.8E -4 ± 7.4E -4 2.0E -3 ± 3.4E -4 0.80 ± 0.02 d E = 6 3.7E -4 ± 2.8E -4 2.0E -3 ± 1.7E -3 0.92 ± 0.02 d E = 6 (Best) 9.1E -5 4.7E -4 0.92

d E = 3 to d E = 6.

patterns, which are associated to a positive largest Lyapunov exponent. This may suggest the combination of the SR model and a delay embedding may require additional investigation as a good geometrical reconstruction of the phase space as stated in Takens' theorem does not guarantee the existence of a parametric ODE model based on the corresponding delay embedding variables. Better performance is reported using an analog forecasting approach. The performance however greatly varies depending on the considered definition of the delay embedding. Using ad hoc parameters (τ = 6, d E = 3), one may retrieve the expected long-term chaotic behavior (λ 1 = 0.87) with a relatively low short-term forecasting error (8.0e -4 for a one-step-ahead forecast). When considering the proposed model, we report for all parametrizations, augmented space dimensions from 3 to 6, performance at least in the same range as the best analog forecasting setting. Besides, when increasing the dimension of the augmented space, we significantly decrease short-term forecasting errors (< 1.e -4 for a one-step-ahead forecast when considering the best fit for d E = 6, i.e. one order of magnitude compared to the best benchmark model) while keeping an appropriate chaotic long-term pattern (λ 1 = 0.92).

Modeling Sea Level Anomaly (SLA): The data driven identification of dynamical representations of real data is an extremely difficult task especially when the underlying processes involve non stable behaviors such as chaotic attractors. This is mainly due to the fact that we do not have any exact knowledge of the closed form of the equations governing the temporal evolution of our variables. Furthermore, the measured quantity may depend on other unobserved variables which makes the exploitation of data-driven techniques highly challenging.

In this context, we report an application to SLA (Sea Level Anomaly) dynamics, which relate to upper ocean dynamics and are monitored by satellite altimeters (Calmant et al., 2008) . Sea surface dynamics are chaotic and clearly involve latent processes, typically subsurface and atmospheric processes. The dataset used in our experiments is a SLA time series obtained using the WMOP product (Juza et al., 2016) . The spatial resolution of our data is a 0.05 • and the temporal resolution h = 1 day. We use the data from January 2009 to December 2014 as training data and we tested our approach on the last month of the year 2014. The considered region is located on south Mallorca (2.5 • E to 4.25 • E, 37.25 • N to 39.5 • N ). Finally, and in order to identify a ROM, our observations are mapped to a low dimensional space using a projection matrix defined offline using a PCA as follow : Regarding the proposed NbedDyn model we consider an augmented latent space with d E = 60. Our model clearly outperforms the three benchmarked schemes with a very significant gain for the forecasting performance at one day (relative gain greater than 90 %) and two days (relative gain greater than 90 %). For a 4-day-ahead forecasting, our model still outperforms the other ones though the gain is lower (relative gain of 40%). Finally, and in order to illustrate the influence of adding extra dimensions to define an augmented latent space on real data, we also tested the proposed NbedDyn model directly on the PCA space (d E = k = 15) this model is referred to as NbedDynZERO and the influence of the latent components is clear from the results in Tab. 2 witch allows a relative gain up to 90 % with respect to the same model directly on the PCA space. We let the reader refer to the Supplementary Material for a more detailed analysis of these experiments, including visual comparisons of the forecasts.

r t = M(x t ) ∈ R k with k =

5. DISCUSSION:

In this work, we address the data-driven identification of latent dynamics for systems which are only partially observed, i.e. when some components of the system of interest are never observed. The reported forecasting performance for Lorenz-63 dynamics is in line with the forecasting performance of state-of-the-art learning-based approaches for a noise-free and fully-observed setting. This is of key interest for real-world applications, where observing systems most often monitor only some components of the underlying dynamics. As a typical example, the SLA forecasting experiment clearly motivates the proposed framework in the context of ocean dynamics for which neither in situ nor satellite observing systems can provide direct observations for all state variables (e.g., subsurface velocities, fine-scale sea surface currents).

We may also further discuss how the proposed framework relates to state-of-the-art dynamical system theory approaches. Most of these approaches rely on delay embedding, as Takens's theorem states the existence of a delay embedding in which the topological properties of the hidden dynamical system are equivalent to those of the true systems up to a diffeomorphic mapping. Hence, stateof-the-art approaches typically combine the selection of a delay embedding representation within classic regression models to represent the one-step-ahead mapping in the considered embedding.

Here, we consider latent dynamics governed by an unknown ODE (4) but we do not explicitly state the latent space. This is however implicit in our forecasting framework. By construction, the considered forecasting model relies on the integration of the learnt ODE ( 4) from an initial condition given as the solution of minimization ( 7). Let us consider the following embedding ψ such that:

ψ ({x t } t0:T ) = arg min u T min {ut} t<T T t=1 x t -G (Φ θ,t (u t-1 )) 2 + λ u t -Φ θ,t (u t-1 ) 2 (9)

Given this embedding, the resulting one-step-ahead forecasting for the observed variable may written as:

EQUATION

Hence, ψ defines a delay embedding representation implicitly stated through minimization ( 7). In this embedding, the dynamics of the observed system x is governed by the composition of observation operator G and forecasting operator Φ θ,t . Regarding the literature on Koopman operator theory, most approaches rely on the explicit identification of eigenfunctions and eigenvalues of the Koopman operator (Koopman, 1931; Brunton et al., 2016a; Tu et al., 2014) . Our framework relates to the identification of the infinitesimal generator f θ of the one-parameter subgroup defined by Koopman operator through the ODE representation (4). By construction, the Koopman operator associated with the identified operator f θ is also diagonalizable, such that the identification of infinitesimal generator f θ provides an implicit decomposition of the Koopman operator of the underlying and unknown dynamical system onto the eigenbasis of the learnt latent dynamics governed by ODE (4).

Future work will further explore methodological aspects, especially the application to highdimensional and stochastic systems. In the considered framework, when considering ROM of the SLA experiment, the projection from the observations space to the low dimensional space was carried using simple PCA projection. Although for the geociences community, using PCA to reduce the dimensionality is motivated by the Galerkin derivation of reduced order models from complex high dimensional governing partial differential equations (Holmes et al., 2012) , using auto-encoders have shown promising results in discovering optimal coordinates when trained jointly with a dynamical system. The combination of the proposed framework with the variational setting considered in the Latent-ODE model (Chen et al., 2018) also appears as an interesting direction for future work. The extension to stochastic systems through the identification of a Stochastic ODE is also of key interest, for instance for future applications of the proposed framework to geophysical random flows, especially to the simulation and forecasting of ocean-atmosphere dynamics in which stochastic components naturally arise (Chapron et al., 2018) .

APPENDIX A PROOF OF PROPOSITION 1:

This proposition can be easily extended to any observation function that doesn't form an embedding of the initial unobserved ODE. However, for the sake of simplicity, we will consider the example given in Eq. ( 1).

Lets suppose a a smooth ODE in the observation space that governs the time evolution of x from Eq. ( 1).

ẋt = f (x t ) x t0 = x 0 (11)

This ODE generates a flow x t = Ψ t (x 0 ).

Since our observation operator is not one-to-one, we can assume the existence of some t, t 1 , t 2 where Real(Φ t(s t1 )) = Real(Φ t(s t2 )) with Real(s t1 ) = Real(s t2 ) (Φ is the flow generated by the unobserved ODE illustrated in Eq. ( 1)). Projecting this equality to the observation space leads to :

Ψ t(x t1 ) = Ψ t(x t2 ) with x t1 = x t2 .

Since the above ODE is smooth (or continuously differentiable), we can show that f is locally Lipschitz on any interval containing t 0 (Sohrab, 2003) which garentees by Picard's Existence Theorem the existance of a unique solution (Coddington & Levinson, 1955) . Formally, for the times t, t 1 , t 2 , Ψ t(x t1 ) = Ψ t(x t2 ) if and only if x t1 = x t2 . This contradicts the assumption that x t1 = x t2 and thus, there is no existence of a t such that Real(Φ t(s t1 )) = Real(Φ t(s t2 )) with Real(s t1 ) = Real(s t2 ).

B AUGMENTED REDUCED ORDER MODELING:

With a view to providing a general formulation which also accounts for reduced-order modeling (ROM) (Champion et al., 2019) , we assume that the observations x t are fully characterized by a lower-dimensional state M(x t ) ∈ R k , with M a global coordinates chart, mapping each observation to its coordinates M(x t ) ∈ R k . Typical examples for operator M includes the identity as in the equations of section 3. In ROM, M may include PCA projections (Benner et al., 2005) for Galerkin derivation of ROM for complex high dimensional partial differential equations (Rowley et al., 2004) . Formally, the augmented state space is formulated as :

EQUATION

with the corresponding SSM :

EQUATION

Based on the above formulation, we can state the optimization problem stated in Eq. ( 5) as:

EQUATION

with the corresponding non constrained formulation :

EQUATION

Finally, when considering forecasting applications of new observations, we can restate the equation 7 in the context of ROM as follow :

ẑn T = arg min z n T min {z n t } t<T T +δT t=T +1 x n t -M -1 (G Φ θ,t u n t-1 ) 2 + λ u n t -Φ θ,t (u n t-1 ) 2 (16)

C LIST OF SYMBOLS s Unobserved state governing the underlying dynamics in R m

x Observed component of s in R n , usually related to some components of s

M:

In the case of ROM (x high dimensional), M is a projection operator from R n to R k with n >> k u Augmented latent space of dimension R d E which consists on a concatenation of an observed components (M(x)) and a latent components z z Latent component concatenated to M(x) to form the augmented state u f (.) Unknown dynamical system governing the unobserved state s Φ t (.) Flow map of the Unknown dynamical system f (.) f θ (.) Dynamical system governing the augmented state u Φ t,θ (.) Flow map of the dynamical system f θ (.)

D DIMENSIONALITY ANALYSIS OF THE NBEDDYN MODEL:

One of the Key parameters of the proposed model is the dimension of the latent space. Despite the fact that it is extremely challenging to get a prior idea of the dimension of the model in the case of real data experiments, one can analyze the spawned manifold of the learnt latent states to get an idea of the true dimension of the underlying model (true here stands for a sufficient dimension of the latent space). The idea here is to compute the modulus of the eigenvalues of the Jacobian matrix for each input of the training data. An eigenvalue does not influence the temporal evolution of the latent state if it has a modulus that tend to zero. The number of non-zero eigenvalues can then be seen as a sufficient dimension of the latent space.

Regarding the identification of an ODE model governing the first state variable of the Lorenz 63 model, Fig. 3 illustrates the eigenvalues of the Jacobian matrix and their modulus for a dimension of the latent space d E = 6. Interestingly, only 3 eigenvalues have non-zero modulus and are effectively influencing the underlying dynamics. This result shows that one can use a 3 dimensional latent-space as a sufficient dimension to identify an ODE model governing the first state of the Lorenz 63 system which is the same dimension as the true Lorenz 63 model. The analysis of the eigenvalues of the Sea Level Anomaly model in the other hand are not as straightforward as in the case of the Lorenz model since we do not have any idea on the analytical form of the underlying dynamical model. Fig. 4 illustrates that using a 60 dimensional latent space for the NbedDyn model, only 50 eigenvalues have non-zero modulus and thus, are effectively influencing the underlying dynamics. The conclusion in this case is that the observed SLA data evolve in a 50 dimensional latent space parametrised by the dynamical model f θ .

E ADDITIONAL FIGURES OF THE LORENZ 63 EXPERIMENT:

We illustrate the forecasting performance of the tested models for the Lorenz-63 experiment through an example of forecasted trajectories in Fig. 5 . Our model with d E = 6 leads to a trajectory similar to the true one up to 7 Lyapunov times, when the best alternative approach diverge from the true trajectory beyond 4 Lyapunov times. An other interesting experiment is to find the initial condition for new observation data. This issue is addressed as presented in section 3 as follow. Given a new noisy and partial observation sequence (Fig. 6 ), we first look for a potential initial condition in the inferred training latent state sequence. This initial condition is then optimized using the cost function described by equation ( 7) to minimize the forecasting error of the new observation sequence. 

G NEURAL NETWORKS HYPERPARAMETERS:

The neural-network parametrization for the dynamical model operator f θ is a bilinear architecture as presented in (Fablet et al., 2018) . The idea of this model is to embed second order polynomial repre-sentations of the input states by combining fully-connected layers and an element-wise product. This bilinear layer are then concatenated to classical fully connected layers to perform regression. This architecture incorporates physical knowledge on the dynamics since physical dynamical systems often involve bilinear non-linearities. Among classic physical models, one may cite for instance advection-diffusion dynamics or shallow water equation.

This model is than integrated numerically to compute the flow-map Φ θ,t using a fourth order Runge Kutta scheme. In practice, the integration scheme is implemented as a ResNet with four layers sharing the dynamical model.

G. 

H SCOOP AND LIMITATIONS:

Constraining limit cycles The proposed augmented ODE formulation does not suppose any prior knowledge on the underlying dynamics responsible for the temporal evolution of the observations. This can lead in some cases (especially when working on chaotic dynamics) to output a dynamical representation that has several attracting regions in addition to the one leading to the observations limit cycle. This can lead to inappropriate results when trying to find an initial condition that forecasts a given observation sequence. The Idea of using the manifold spanned by the augmented training data allows to bypass this issue but we believe that adding additional constraints (energy preserving constraints, known symmetries in the models ...etc.) can significantly improve the quality of the learnt dynamical models. 

I CODE SAMPLE:



:

The word smooth here stands for continuously differentiable or C 1 .

For simplicity, the observed component is considered in the equations as the direct observations. However, in general, when considering Reduced Order Models (ROMs), the observed component can be seen as a low dimensional projection of the observations. Please refer to the appendix for a similar derivation of the equations in the case of reduced order modeling.

The results of the neural networks based models were averaged over 5 runs except for NbedDyn dE =

and dE =

where we only launched 3 runs. We will include other runs for the final version of the manuscript

Figures and Tables:

Figure 1: Illustration for a 2-dimensional linear ODE: Forecasted augmented latent space with respect to the true states given (left) the same initial condition as the training sequence and (right) a different initial condition. Given the initial condition we illustrate both the prediction (forecast up to the end of the training time) of the trained model and the extrapolation (forecast beyond the training time) performances with respect to the true trajectory .The projection of the solutions in the real plane illustrates the forecasting of the observations.

Figure 3: Analysis of the eigenvalues of the NbedDyn model Jacobian matrix.: Lorenz-63 casestudy with d E = 6

Figure 4: Analysis of the eigenvalues of the NbedDyn model Jacobian matrix.: Sea Level Anomaly case-study with d E = 60

Figure 5: Generated time series of the proposed models. : Given the same initial condition, we generated a time series of 1000 time steps.

Figure6: Forecasted Lorenz state sequence given noisy and partial observations: Given noisy and partial observations, our model optimizes equation (7) to infer an initial condition that minimize the forecasting of the observations.

Figure7: Generated time series of the proposed models for the forecasting of the SLA dynamics. : Given the same initial condition, we generated a time series of 150 days.

Figure 9: Forecasted SLA states of the proposed models.

i m p o r t l i b s from g e n e r a t e d a t a import g e n e r a t e d a t a import numpy a s np import t o r c h from t o r c h . a u t o g r a d import V a r i a b l e s e e d = 0 np . random . s e e d ( s e e d ) t o r c h . m a n u a l s e e d ( s e e d ) # G e n e r a t e d a t a c l a s s GD: model = ' L o r e n z 6 3 ' c l a s s p a r a m e t e r s : s i g m a = 1 0 . 0 r h o = 2 8 . 0 b e t a = 8 . 0 / 3 d t i n t e g r a t i o n = 0 . 0 1 # i n t e g r a t i o n t i m e n b l o o p d a t a = 5 0 . 0 1 # s i z e o f t h e c a t a l o g t e s t s a m p l e s = 1000 # r u n t h e d a t a g e n e r a t i o n c a t a l o g , x t , yo = g e n e r a t e d a t a (GD) X t e s t = c a t a l o g . a n a l o g s [ c a t a l o g . a n a l o g s . s h a p e [0] -GD. t e s t s a m p l e s : , : 1 ] X t r a i n = c a t a l o g . a n a l o g s [ : c a t a l o g . a n a l o g s . s h a p e [0] -GD. t e s t s a m p l e s , : 1 ] G r a d t = np . g r a d i e n t ( X t r a i n [ : , 0 ] ) . r e s h a p e ( X t r a i n . s h a p e [ 0 ] , 1 ) \ /GD. d t i n t e g r a t i o n x = V a r i a b l e ( t o r c h . from numpy ( X t r a i n ) . f l o a t ( ) ) z = V a r i a b l e ( t o r c h . from numpy ( G r a d t ) . f l o a t ( ) ) d a t a s i z e = X t r a i n . s h a p e [ 0 ] # n e u r a l n e t params p a r a m s = {} p a r a m s [ ' t r a n s i t i o n l a y e r s ' ] = 1 p a r a m s [ ' b i l i n e a r l a y e r s ' ] = 6 p a r a m s [ ' d i m h i d d e n l i n e a r ' ] = 6 p a r a m s [ ' d i m i n p u t ' ] = 1 p a r a m s [ ' d i m l a t e n t ' ] = 5 p a r a m s [ ' d i m o u t p u t ' ] = p a r a m s [ ' d i m i n p u t ' ] + p a r a m s [ ' d i m l a t e n t ' ] p a r a m s [ ' d t i n t e g r a t i o n ' ] = GD. d t i n t e g r a t i o n # D y n a m i c a l model c l a s s F C n e t ( t o r c h . nn . Module ) : d e f i n i t ( s e l f , p a r a m s ) : s u p e r ( FC ne t , s e l f ) .i n i t ( ) y a u g = np . random . u n i f o r m ( s i z e = ( d a t a s i z e , p a r a m s [ ' d i m l a t e n t ' ] ) ) s e l f . y a u g = t o r c h . nn . P a r a m e t e r ( t o r c h . from numpy ( y a u g ) . f l o a t ( ) ) s e l f . l i n e a r C e l l = t o r c h . nn . L i n e a r ( p a r a m s [ ' d i m o u t p u t ' ] \ , p a r a m s [ ' d i m h i d d e n l i n e a r ' ] ) s e l f . B l i n e a r C e l l 1 = \ t o r c h . nn . M o d u l e L i s t ( \ [ t o r c h . nn . L i n e a r ( p a r a m s [ ' d i m o u t p u t ' ] , 1 ) \ f o r i i n range ( p a r a m s [ ' b i l i n e a r l a y e r s ' ] ) ] ) s e l f . B l i n e a r C e l l 2 = \ t o r c h . nn . M o d u l e L i s t ( \ [ t o r c h . nn . L i n e a r ( p a r a m s [ ' d i m o u t p u t ' ] , 1 ) \ f o r i i n range ( p a r a m s [ ' b i l i n e a r l a y e r s ' ] ) ] ) a u g m e n t e d s i z e = p a r a m s [ ' b i l i n e a r l a y e r s ' ] \ + p a r a m s [ ' d i m h i d d e n l i n e a r ' ] s e l f . t r a n s L a y e r s = \ t o r c h . nn . M o d u l e L i s t ( \ [ t o r c h . nn . L i n e a r ( a u g m e n t e d s i z e , p a r a m s [ ' d i m o u t p u t ' ] ) ] ) s e l f . t r a n s L a y e r s . e x t e n d ( \ [ t o r c h . nn . L i n e a r ( p a r a m s [ ' d i m o u t p u t ' ] , p a r a m s [ ' d i m o u t p u t ' ] ) \ f o r i i n range ( 1 , p a r a m s [ ' t r a n s i t i o n l a y e r s ' ] ) ] ) s e l f . o u t p u t L a y e r = t o r c h . nn . L i n e a r ( p a r a m s [ ' d i m o u t p u t ' ] , \ p a r a m s [ ' d i m o u t p u t ' ] ) d e f f o r w a r d ( s e l f , i n p , d t ) : i f i n p . s h a p e [-1]< p a r a m s [ ' d i m l a t e n t ' ] + p a r a m s [ ' d i m i n p u t ' ] : a u g i n p = t o r c h . c a t ( ( i n p , s e l f . y a u g ) , dim = 1 ) e l s e : a u g i n p = i n p B P o u t p = V a r i a b l e ( t o r c h . z e r o s ( ( a u g i n p . s i z e ( ) [ 0 ] , \ p a r a m s [ ' b i l i n e a r l a y e r s ' ] ) ) ) L o u t p = s e l f . l i n e a r C e l l ( a u g i n p ) f o r i i n range ( ( p a r a m s [ ' b i l i n e a r l a y e r s ' ] ) ) : B P o u t p [ : , i ] = s e l f . B l i n e a r C e l l 1 [ i ] ( a u g i n p ) [ : , 0 ] * \ s e l f . B l i n e a r C e l l 2 [ i ] ( a u g i n p ) [ : , 0 ] a u g v e c t = t o r c h . c a t ( ( L o u t p , B P o u t p ) , dim = 1 ) f o r i i n range ( ( p a r a m s [ ' t r a n s i t i o n l a y e r s ' ] ) ) : a u g v e c t = ( s e l f . t r a n s L a y e r s i ] ( a u g v e c t ) ) g r a d = s e l f . o u t p u t L a y e r ( a u g v e c t ) r e t u r n g r a d , a u g i n p model = F C n e t ( p a r a m s ) # c o m p u t e f l o w : RK4 c l a s s I N T n e t ( t o r c h . nn . Module ) : f . a d d m o d u l e ( ' D y n n e t ' , F C n e t ( params ) ) s e l f . D y n n e t = model d e f f o r w a r d ( s e l f , i n p , d t ) : k1 , a u g i n p = s e l f . D y n n e t ( i n p , d t ) i n p k 2 = a u g i n p + 0 . 5 * p a r a m s [ ' d t i n t e g r a t i o n ' ] * k1 k2 , tmp = s e l f . D y n n e t ( i n p k 2 , d t ) i n p k 3 = a u g i n p + 0 . 5 * p a r a m s [ ' d t i n t e g r a t i o n ' ] * k2 k3 , tmp = s e l f . D y n n e t ( i n p k 3 , d t ) i n p k 4 = a u g i n p + p a r a m s [ ' d t i n t e g r a t i o n ' ] * k3 k4 , tmp = s e l f . D y n n e t ( i n p k 4 , d t ) p r e d = a u g i n p + d t * ( k1 +2 * k2 +2 * k3+k4 ) / 6 r e t u r n p r e d , k1 , i n p , a u g i n p # I n s t a n c i a t e t h e model modelRINN = I N T n e t ( p a r a m s ) c r i t e r i o n = t o r c h . nn . MSELoss ( r e d u c t i o n = ' e l e m e n t w i s e m e a n ' ) o p t i m i z e r = t o r c h . o p t i m . Adam ( model . p a r a m e t e r s ( ) ) # P r e t r a i n i n g : f i t t h e g r a d i e n t p a r a m s [ ' n t r a i n ' ] = [ 3 0 0 0 0 0 , 1 0 0 0 0 ] f o r t i n range ( p a r a m s [ ' n t r a i n ' ] [ 0 ] ) : p r e d , g r a d , i n p , a u g i n p = modelRINN ( x , p a r a m s [ ' d t i n t e g r a t i o n ' ] ) l o s s 1 = c r i t e r i o n ( g r a d [ : , : 1 ] , z ) l o s s 2 = c r i t e r i o n ( p r e d [ : -1 , : ] , a u g i n p [ 1 : , : ] ) l o s s = 0 . 9 * l o s s 1 + 0 . 1 * l o s s 2 p r i n t ( t , l o s s ) o p t i m i z e r . z e r o g r a d ( ) l o s s . b a c k w a r d ( r e t a i n g r a p h = T r u e ) o p t i m i z e r . s t e p ( ) # t r a i n i n g f o r t i n range ( p a r a m s [ ' n t r a i n ' ] [ 1 ] ) : p r e d , g r a d , i n p , a u g i n p = modelRINN ( x , p a r a m s [ ' d t i n t e g r a t i o n ' ] ) l o s s 1 = c r i t e r i o n ( p r e d [ : -1 , : ] , a u g i n p [ 1 : , : ] ) l o s s 2 = c r i t e r i o n ( p r e d [ : -1 , 1 : ] , a u g i n p [ 1 : , 1 : ] ) l o s s = 1 0 . 0 * l o s s 1 + 1 . 0 * l o s s 2 p r i n t ( t , l o s s ) o p t i m i z e r . z e r o g r a d ( ) l o s s . b a c k w a r d ( r e t a i n g r a p h = T r u e ) o p t i m i z e r . s t e p ( )

Forecasting performance on the test set of data-driven models for Lorenz-63 dynamics where only the first variable is observed: first two columns : mean RMSE for different forecasting time steps, third column : largest Lyapunov exponent of a predicted series of length of 10000 timesteps (The true largest Lyapunov exponent of the Lorenz 63 model is 0.91(Sprott, 2003)).

Table content: <table/>

15 which amounts to capture 92% of the total variance. The proposed framework is applied simply by replacing x t by the PCA components r t . Please refere to the appendix for a generalization of the proposed approache in the case of ROM.

Table content: <table><tr><td>Model</td><td/><td>t0 + h</td><td>t0 + 2h</td><td>t0 + 4h</td></tr><tr><td>AF</td><td>RMSE Corr</td><td>0.036 98.93%</td><td>0.049 96.97%</td><td>0.067 93.99%</td></tr><tr><td>SR</td><td>RMSE Corr</td><td>0.014 99.42%</td><td>0.021 97.63%</td><td>xx 0.037 90.91%</td></tr><tr><td>Latent-ODE RNN NbedDynZERO NbedDyn</td><td>RMSE Corr RMSE Corr RMSE Corr RMSE Corr</td><td colspan="3">0.030 ± 0.05 98.20% ± 0.39% 0.026 ± 0.003 98.36% ± 0.40% 0.016 ± 0.0 99.44% ± 0.0% 0.002 ± 0.0003 99.99% ± 0.0017% % 99.91% ± 0.01% 99.01% ± 0.04% 0.031 ± 0.031 0.040 ± 0.040 97.39% ± 0.36% 93.42% ± 0.55% 0.038 ± 0.007 0.053 ± 0.016 95.29% ± 1.73% 74.97% ± 5.75% 0.023 ± 0.0 0.038 ± 0.0 97.71% ± 0.0% 91.18% ± 0.0% 0.006 ± 0.001 0.020 ± 0.004</td></tr></table>

SLA Forecasting performance on the test set of data-driven models: RMSE and correlation coefficients for different forecasting time steps.

Table content: <table><tr><td>We report forecasting performance</td></tr><tr><td>for our model and include a com-</td></tr><tr><td>parison with analog methods (AF),</td></tr><tr><td>Sparse regression (SR), LSTM</td></tr><tr><td>(RNN) and a neural ODE setting</td></tr><tr><td>(Latent-ODE) in Tab. 2 (The re-</td></tr><tr><td>sults of the neural networks based</td></tr><tr><td>models were averaged over 5 runs).</td></tr></table>

RNN parameters in the Lorenz 63 Experiment.

Table content: <table><tr><td colspan="2">1 LORENZ 63 EXPERIMENTS HYPERPARAMETERS</td><td/><td/></tr><tr><td>Parameter</td><td>Value</td><td>Parameter</td><td>Value</td></tr><tr><td>Number of LSTM layers Hidden size Sequence length Learning rate Optimizer Training data</td><td>10 10 30 0.001 Adam 4000</td><td>Latent dimension Hidden size RNN hidden size Learning rate Optimizer Training data</td><td>4 15 100 0.01 Adam 4000</td></tr></table>

Latent-ODE parameters in the Lorenz 63 Experiment, please refer to (Chen et al., 2018) for more details.

Table content: <table><tr><td>Parameter</td><td>Value</td></tr><tr><td>Augmented Latent dimension</td><td>6</td></tr><tr><td>Number of bilinear layers</td><td>6</td></tr><tr><td>Number of linear layers</td><td>6</td></tr><tr><td>Integration scheme</td><td>Runge Kutta 4</td></tr><tr><td>Learning rate</td><td>0.001</td></tr><tr><td>Optimizer</td><td>Adam</td></tr><tr><td>Training data</td><td>4000</td></tr></table>

NbedDyn parameters in the Lorenz 63 Experiment, please refer to(Fablet et al., 2018) for more details.

Table content: <table><tr><td colspan="2">G.2 SLA EXPERIMENTS HYPERPARAMETERS</td><td/><td/></tr><tr><td>Parameter</td><td>Value</td><td>Parameter</td><td>Value</td></tr><tr><td>Number of LSTM layers Hidden size Sequence length Learning rate Optimizer Training data</td><td>5 20 40 0.001 Adam 2000</td><td>Latent dimension Hidden size RNN hidden size Learning rate Optimizer Training data</td><td>60 70 200 0.01 Adam 2000</td></tr><tr><td colspan="2">Table 6: RNN parameters in the SLA Exper-</td><td/><td/></tr><tr><td>iment.</td><td/><td/><td/></tr></table>

Latent-ODE parameters in the SLA Experiment, please refer to(Chen et al., 2018) for more details.

Table content: <table><tr><td>Parameter</td><td>Value</td></tr><tr><td>Augmented Latent dimension</td><td>60</td></tr><tr><td>Number of bilinear layers</td><td>60</td></tr><tr><td>Number of linear layers</td><td>60</td></tr><tr><td>Integration scheme</td><td>Runge Kutta 4</td></tr><tr><td>Learning rate</td><td>0.001</td></tr><tr><td>Optimizer</td><td>Adam</td></tr><tr><td>Training data</td><td>2000</td></tr></table>

NbedDyn parameters in the SLA Experiment, please refer to(Fablet et al., 2018) for more details.

Table content: <table/>

Peer review (iclr2020):
Update: I raised the score from 1 to 3 to acknowledge the authors* consideration for the 2000-2010 literature on learning dynamical systems from partial observations. Unfortunately, the writing is still confusing, some of the claims in the introduction and rebuttal are inexact ([5] does not embed the observations and does work with partially observed environments), and the method lacks originality compared to existing work. Other work that relies on ODE integration and in finding high-dimensional state variables has recently been published and tested on more ambitious datasets, e.g. Ayed et al, *Learning Dynamical Systems from Partial Observations*, arXiv 2019. *** TL;DR: relatively well written (if sometimes confusing) paper that reinvents the inference of latent variables in nonlinear dynamical systems that has been published in the 2000s, and that misses an important chunk of literature (and experiments on dynamical systems such as Lorenz-63) from that time. This paper proposes an approach for learning dynamical systems from partial observations x, by using an augmented state variable z that follows dynamics that can be described by an ordinary differential equation (ODE) with dynamics f. The authors motivate their work by the problem of dynamical system identification when only partial observations are available. The authors claim that was to date primarily addressed using time-delay embedding, following Takens* theorem. The authors introduce s-dimensional unknown state variables z, dynamical function f (for the ODE on z), flow phi on z, limit cycles on z, observation function H: z -> x that goes from z to n-dimensional observations x, low k-dimensional manifold r (with a map M: x -> r), and state augmentation variable y. The reconstructed state variable u is the concatenation of r and y. One key ingredient of the method is to infer the optimal value of state augmentation variable y during learning (see equations 5 and 6) and inference for forecasting (7); this is not well explained in the abstract and introduction. I would note that the problem of state space modeling (SSM) and dynamical system identification has been well studied, and the notation and reformulation in this paper is somewhat confusing for those who are used to the notation in SSMs (specifically, expressing the observation approximation as M^{-1}(G(phi(u_{t-1}))). Learning a state-space model involves both learning parameters and inferring the latent states representation (or, in graphical models, the distribution of these latent states) given the parametric models. One approach has been to formulate the state-space model learning by maximum likelihood learning of the model parameters so that the generative model fits observed data x, and this would involve factoring out the distribution of latent states z; the algorithm would rely on Expectation Maximisation, and could involve variational approximations or sampling. While the state space models were hampered by their linearity, several papers in 2000s showed how it is possible to learn nonlinear dynamical models, e.g. [4], [5], [6] and [7] to cite earlier ones. Equations (5) and (6) are similar to the standard equations for a dynamical system expressed in continuous time, with the only difference that the optimisation is with respect to y only, rather than w.r.t. z or u (why not 	ilde z or hat z?). The paper mentions various initialisation strategies for y (last paragraph of section 3). Why not predict from the past of the observations, like is done in many other similar work? The literature review mixes older and newer references. For example, on page 1, I would note that the Takens* theorem has been applied in conjuction with Support Vector Regression as early as 1999 [1][2], and with neural networks in 1993 [3]. Most importantly, the ideas of this paper have already been published in [4] (with architecture constraints on the neural network state-space model), in [5] (with any nonlinear neural network state-space model), in [6] (using Restricted Boltzmann Machines) and in [7] (using Gaussian Process latent variable models). The model is illustrated with experiments on a 2D linear attractor, on the Lorenz-63 attractor. Given the results published in [1] and [2] using SVR on 1D observations of that attractor, and in [5] using a recurrent neural network, I am unconvinced by these results. It seems in particular that the number of training points (around 4000) limits the performance of RNN / LSTM models. The application to Sea Level Anomaly is interesting. Minor comments: *Unfortunately, When* (page 1) There is a missing -1 after M in equation (5) and (10) In equation (7), should not the sum go from t=0 to T, as x_t is unknown for t>T? What is prediction and what is extrapolation on Figure 1? The caption of Fig 1 contains (left) The figures seem squeezed with the captions / titles un-aesthetically wide. Labels on Figure 5 in the appendix seem mixed, and red should be the ground truth [1] Mattera & Haykin (1999) *Support vector machines for dynamic reconstruction of a chaotic system* [2] Muller, Smola, Ratsch, Scholkopf, Kohlmorgen & Vapnik (1999) *Using support vector machines for time-series prediction* [3] Wan (1994) *Time series prediction by using a connectionist network with internal delay lines* [4] Ghahramani, and Roweis (1999) *Learning nonlinear dynamical systems using an EM algorithm* [5] Mirowski & LeCun (2009) *Dynamic Factor Graphs for Time Series Modeling* [6] Taylor, Hinton & Roweis (2006) *Modeling human motion using binary latent variables* [7] Wang, Fleet & Hertzmann (2006) *Gaussian process dynamical models*