Year: 2019

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.

heading: 1 INTRODUCTION
text: Learning the underlying dynamical representation of some observed variables xt ∈ Rn (where t ∈ {t0, ..., 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 xt relate to some underlying deterministic states st 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 st from a representative dataset of observations {xti}i. Unfortunately, when the observed variables xt only relate to some but not all the components of underlying states st, 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; Mü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).

heading: 2 BACKGROUND AND RELATED WORK
text: 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 ṡt = f(st). Let us also assume that this system generates a flow Φti(st0) =∫ ti t0 f(zw)dw ∈ Rm with trajectories that are asymptotic to a limit-cycle L of dimension d contained in Rm. We further assume that we are provided with a measurement function H that maps our state variable s to our observations xt = H(st) ∈ Rn. 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 Rn 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 Rn if H is one-to-one and is an immersion of L in Rn. 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-
1The word smooth here stands for continuously differentiable or C1.
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.

heading: 3 LEARNING LATENT REPRESENTATIONS OF PARTIALLY-OBSERVED DYNAMICS
text: Augmented latent dynamics: Let us consider a continuous m-dimensional dynamical system st governed by an autonomous ODE ṡt = f(st) with Φt the corresponding flow Φt(st0) =∫ t t0 f(sw)dw. We assume the state st to be never fully-observed. Formally, we can define an obser-
vation function H : Rm −→ Rn such that the observation xt follows xt = H(st). An example of interest comprises the situation when one observes only the n first variables of st.
Our goal is to derive an ODE representation of xt ∈ Rn. However, the key question arising here is the extent to which the dynamics expressed in the observations space, reflect the true underlying dynamics in Rm, 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: {
ṡt = αst st0 = s0 (1)
with s ∈ C a state variable and α ∈ C a complex imaginary number. The solution of this problem is st = Ke αt (2)
with K a constant depending on s0. Let us assume now that we are only provided with the real part as direct measurements of the unobserved state i.e. H(.) = Real(.) : xt = Real(st). Proposition 1 : The flow of an ODE cannot represent the time evolution of xt.
The proof of the proposition is given in the appendix. The intuition exemplified for the 2- dimensional system given by (1) is as follows. If we assume that we are only provided the real part as direct measurements xt ∈ R of the true states st, 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 st0 and st1 correspond to two states that have the same real part but distinct imaginary parts, the associated observed states are equal xt0 = xt1 . However, the time evolution of the states st0 and st1 differ if they have different imaginary parts. The observed states xt0+δ and xt1+δ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 operatorH of a deterministic underlying dynamical system that governs st, one can define a delay embedding from an observation sequence, i.e. {xt0−aδ, xt0 , . . . , xt0+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 ut ∈ RdE define a dE-dimensional augmented latent state formed by an observed component2 xt and an unobserved latent component referred to as zt.
ut T = [xTt , z T t ] (3)
Assuming that the augmented latent state ut is governed by an ODE leads to the the following state space model : {
u̇t = fθ(ut) xt = G(ut) (4)
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 xt = G(ut). It means that operator G only returns the observed component of the augmented
2For 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.
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 xt. 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 {x0, . . . ,xT } 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 {xT+1, . . . ,xT+N}. Both issues involve the joint inference of the unobserved component zt of the augmented latent states. This formulation generalizes to reduced-order modeling (ROM) (Champion et al., 2019), where observation xt 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 {x0, . . . ,xT } 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 zt of the augmented state ut is never observed, this identification requires the joint optimization of the model parameters θ as well as of the hidden component zt. Formally, this problem is stated as the following minimization of the forecasting error on observed variables:
θ̂ = arg min θ min {zt}t T∑ t=1 ‖xt −G (Φθ,t (ut−1))‖2
Subject to
{ ut = Φθ,t(ut−1)
G(ut) = xt
(5)
with Φθ,t the one-step-ahead diffeomorphic flow associated with operator fθ such that:
Φθ,t(ut−1) = ut−1 + ∫ t t−1 fθ(uw)dw
In (5), the loss to be minimized involves the one-step-ahead forecasting error for the observed variable xt. The constraints state that the augmented latent state sequence {u0, . . . ,uT } 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 ‖xt −G (Φθ,t (ut−1))‖2 + λ‖ut − Φθ,t(ut−1)‖2 (6)
where uTt = [x T t , z T t ] and λ a weighting parameter. The term ‖ut −Φθ,t(ut−1)‖2 may be regarded as a regularization term such that the inference of the unobserved component zt of the augmented state ut 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 4th-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 {x0, . . . ,xT }, we minimize (6) jointly w.r.t. θ and unobserved variables {z0, . . . , zT }. Application to forecasting: We also apply the proposed framework to the forecasting of the observed states xt. Given a trained latent dynamical model (4), forecasting future states for xt relies on the forecasting of the entire augmented latent state ut. The latter amounts to determining an initial condition of the unobserved component zt and performing a numerical integration of the trained ODE (4).
Let us denote by xnt , t ∈ {t0, ..., T} a new series of observed states. We aim to forecast future states xnt , t ∈ {T + 1, ..., T + δT}. Following (6), we infer the unobserved component ẑT of the augmented state unT at time T from the following minimization:
ẑnT = arg min znT min {znt }t<T T∑ t=t0 ‖xnt −G ( Φθ,t ( unt−1) )∥∥2 + λ‖unt − Φθ,t(unt−1)‖2 (7)
Here, we only minimize w.r.t. latent variables {znt } 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 xnt . 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.

heading: 4 NUMERICAL EXPERIMENTS
text: 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:{
ṡt = αst st0 = s0 (8)
with α = −0.1 − 0.5j, j2 = −1 and s0 = 0.5. As α ∈ C with Real(α) < 0 and s0 6= 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 xt = Real(st). 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 ut = [xt, 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 xt = st,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 dE 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 τMI and τCorr. Regarding the dimension of the embedding representation, we used the Takens embedding condition dE = 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, dE = 3), however, it does not capture well the chaotic
3The results of the neural networks based models were averaged over 5 runs except for NbedDyn dE = 4 and dE = 5 where we only launched 3 runs. We will include other runs for the final version of the manuscript
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, dE = 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 dE = 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 : rt = M(xt) ∈ Rk with k = 15 which amounts to capture 92% of the total variance. The proposed framework is applied simply by replacing xt by the PCA components rt. Please refere to the appendix for a generalization of the proposed approache in the case of ROM.
We report forecasting performance for our model and include a comparison with analog methods (AF), Sparse regression (SR), LSTM (RNN) and a neural ODE setting (Latent-ODE) in Tab. 2 (The results of the neural networks based models were averaged over 5 runs). Regarding the proposed NbedDyn model we consider an augmented latent space with dE = 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 (dE = 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.

heading: 5 DISCUSSION
text: 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:
ψ ({xt}t0:T ) = arg minuT min{ut}t<T T∑ t=1 ‖xt −G (Φθ,t (ut−1))‖2 + λ‖ut − Φθ,t(ut−1)‖2 (9)
Given this embedding, the resulting one-step-ahead forecasting for the observed variable may written as: xT+1 = G (Φθ,t (ψ ({xt}t=t0:T ))) (10) 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 operatorG 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).

heading: A PROOF OF PROPOSITION 1
text: 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).
{ ẋt = f(xt)
xt0 = x0 (11)
This ODE generates a flow xt = Ψt(x0).
Since our observation operator is not one-to-one, we can assume the existence of some t̂, t1, t2 where Real(Φt̂(st1)) = Real(Φt̂(st2)) with Real(st1) 6= Real(st2) (Φ is the flow generated by the unobserved ODE illustrated in Eq. (1)). Projecting this equality to the observation space leads to : Ψt̂(xt1) = Ψt̂(xt2) with xt1 6= xt2 . Since the above ODE is smooth (or continuously differentiable), we can show that f is locally Lipschitz on any interval containing t0 (Sohrab, 2003) which garentees by Picard’s Existence Theorem the existance of a unique solution (Coddington & Levinson, 1955). Formally, for the times t̂, t1, t2, Ψt̂(xt1) = Ψt̂(xt2) if and only if xt1 = xt2 . This contradicts the assumption that xt1 6= xt2 and thus, there is no existence of a t̂ such that Real(Φt̂(st1)) = Real(Φt̂(st2)) with Real(st1) 6= Real(st2).

heading: B AUGMENTED REDUCED ORDER MODELING
text: 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 xt are fully characterized by a lower-dimensional state M(xt) ∈ Rk, with M a global coordinates chart, mapping each observation to its coordinatesM(xt) ∈ Rk. Typical examples for operatorM 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 :
ut T = [M(xt)T , zTt ] (12)
with the corresponding SSM : { u̇t = fθ(ut)
xt =M−1(G(ut)) (13)
Based on the above formulation, we can state the optimization problem stated in Eq. (5) as:
θ̂ = arg min θ min {zt}t T∑ t=1 ‖xt −M−1(G (Φθ,t (ut−1))) ‖2
Subject to  ut = Φθ,t(ut−1)M−1(G(ut)) = xt (14)
with the corresponding non constrained formulation :
min θ min {zt}t T∑ t=1 ‖xt −M−1(G (Φθ,t (ut−1))) ‖2 + λ‖ut − Φθ,t(ut−1)‖2 (15)
Finally, when considering forecasting applications of new observations, we can restate the equation 7 in the context of ROM as follow :
ẑnT = arg min znT min {znt }t<T T+δT∑ t=T+1 ‖xnt −M−1(G ( Φθ,t ( unt−1) )) ‖2 + λ‖unt − Φθ,t(unt−1)‖2 (16)

heading: C LIST OF SYMBOLS
text: s Unobserved state governing the underlying dynamics in Rm
x Observed component of s in Rn, usually related to some components of s
M In the case of ROM (x high dimensional),M is a projection operator from Rn to Rk with n >> k
u Augmented latent space of dimension RdE which consists on a concatenation of an observed components (M(x)) and a latent components z
z Latent component concatenated toM(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θ(.)

heading: D DIMENSIONALITY ANALYSIS OF THE NBEDDYN MODEL
text: 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 dE = 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θ.

heading: E ADDITIONAL FIGURES OF THE LORENZ 63 EXPERIMENT
text: 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 dE = 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.

heading: F ADDITIONAL FIGURES OF THE SEA LEVEL ANOMALY EXPERIMENT
text: Forecasted states of the Sea Level Anomaly are illustrated in Fig. 8 and 9. The visual analysis of the forecasted SLA states emphasize the relevance of the proposed NbedDyn model. While state of the art approaches generally overestimate the time evolution of some structures such as eddies, our model is the only one to give near perfect forecasting up to 4 days.
0.0 2.5 5.0 7.5 10.0 12.5 EOF mode
−0.25
0.00
0.25 0.50 V a lu e
Day 1 forecast Forecast True state
0.0 2.5 5.0 7.5 10.0 12.5 EOF mode
−0.25
0.00
0.25
0.50
V a lu e
Day 2 forecast Forecast True state
0.0 2.5 5.0 7.5 10.0 12.5 EOF mode
−0.25
0.00
0.25
0.50
0.75
V a lu e
Day 3 forecast Forecast True state
0.0 2.5 5.0 7.5 10.0 12.5 EOF mode
−0.25
0.00
0.25
0.50
V a lu e
Day 4 forecast Forecast True state
(a) Analog forecasting.

heading: G NEURAL NETWORKS HYPERPARAMETERS
text: 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.1 LORENZ 63 EXPERIMENTS HYPERPARAMETERS
G.2 SLA EXPERIMENTS HYPERPARAMETERS

heading: H SCOOP AND LIMITATIONS
text: 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.

heading: I CODE SAMPLE
text: # 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 as 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 seed = 0 np . random . seed ( seed ) t o r c h . m a n u a l s e e d ( seed ) # Genera te da ta c l a s s GD:
model = ’ Lorenz 63 ’ c l a s s p a r a m e t e r s :
s igma = 1 0 . 0 rho = 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 = 50 .01 # 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
# run t h e da ta g e n e r a t i o n c a t a l o g , xt , 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 . shape [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 . shape [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 . shape [ 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 . shape [ 0 ] # n e u r a l n e t params params = {} params [ ’ t r a n s i t i o n l a y e r s ’ ]=1 params [ ’ b i l i n e a r l a y e r s ’ ]=6 params [ ’ d i m h i d d e n l i n e a r ’ ] = 6 params [ ’ d i m i n p u t ’ ]=1 params [ ’ d i m l a t e n t ’ ]=5 params [ ’ d i m o u t p u t ’ ]= params [ ’ d i m i n p u t ’ ] + params [ ’ d i m l a t e n t ’ ] params [ ’ 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 # Dynamical model c l a s s FC net ( t o r c h . nn . Module ) :
def i n i t ( s e l f , params ) : super ( FC net , s e l f ) . i n i t ( ) y aug = np . random . un i fo rm ( s i z e =( d a t a s i z e , params [ ’ d i m l a t e n t ’ ] ) ) s e l f . y aug = t o r c h . nn . P a r a m e t e r ( t o r c h . from numpy ( y aug ) . 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 ( params [ ’ d i m o u t p u t ’ ]\
, params [ ’ 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 . Modu leL i s t (\ [ t o r c h . nn . L i n e a r ( params [ ’ d i m o u t p u t ’ ] , 1 )\ f o r i in range ( params [ ’ 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 . Modu leL i s t (\ [ t o r c h . nn . L i n e a r ( params [ ’ d i m o u t p u t ’ ] , 1 )\ f o r i in range ( params [ ’ 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 = params [ ’ b i l i n e a r l a y e r s ’ ]\ + params [ ’ 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 . Modu leL 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 , params [ ’ 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 ( params [ ’ d i m o u t p u t ’ ] , params [ ’ d i m o u t p u t ’ ] ) \ f o r i in range ( 1 , params [ ’ 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 ( params [ ’ d i m o u t p u t ’ ] ,\
params [ ’ d i m o u t p u t ’ ] ) def f o r w a r d ( s e l f , inp , d t ) :
i f i n p . shape [−1]<params [ ’ d i m l a t e n t ’ ]+ params [ ’ d i m i n p u t ’ ] : a u g i n p = t o r c h . c a t ( ( inp , s e l f . y aug ) , dim =1) e l s e : a u g i n p = i n p BP outp = 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 ] , \ params [ ’ b i l i n e a r l a y e r s ’ ] ) ) ) L ou tp = s e l f . l i n e a r C e l l ( a u g i n p ) f o r i in range ( ( params [ ’ b i l i n e a r l a y e r s ’ ] ) ) :
BP outp [ : , 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 outp , BP outp ) , dim =1) f o r i in range ( ( params [ ’ 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 ad = s e l f . o u t p u t L a y e r ( a u g v e c t ) re turn grad , a u g i n p
model = FC net ( params ) # compute f l o w : RK4 c l a s s INT ne t ( t o r c h . nn . Module ) :
def i n i t ( s e l f , params ) : super ( INT net , s e l f ) . i n i t ( )
# s e l f . add module ( ’ Dyn ne t ’ , FC net ( params ) ) s e l f . Dyn net = model
def f o r w a r d ( s e l f , inp , d t ) : k1 , a u g i n p = s e l f . Dyn net ( inp , d t ) i n p k 2 = a u g i n p + 0 . 5∗ params [ ’ d t i n t e g r a t i o n ’ ]∗ k1 k2 , tmp = s e l f . Dyn net ( inp k2 , d t ) i n p k 3 = a u g i n p + 0 . 5∗ params [ ’ d t i n t e g r a t i o n ’ ]∗ k2 k3 , tmp = s e l f . Dyn net ( inp k3 , d t ) i n p k 4 = a u g i n p + params [ ’ d t i n t e g r a t i o n ’ ]∗ k3 k4 , tmp = s e l f . Dyn net ( inp k4 , d t ) p r ed = a u g i n p + d t ∗ ( k1+2∗k2+2∗k3+k4 ) / 6 re turn pred , k1 , inp , a u g i n p
# I n s t a n c i a t e t h e model modelRINN = INT ne t ( params ) 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 emen twi se mean ’ ) o p t i m i z e r = t o r c h . opt im . 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 params [ ’ n t r a i n ’ ] = [ 3 0 0 0 0 0 , 1 0 0 0 0 ] f o r t in range ( params [ ’ n t r a i n ’ ] [ 0 ] ) :
pred , grad , inp , a u g i n p = modelRINN ( x , params [ ’ 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 rad [ : , : 1 ] , z ) l o s s 2 = c r i t e r i o n ( p red [ : −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 . backward ( r e t a i n g r a p h =True ) 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 in range ( params [ ’ n t r a i n ’ ] [ 1 ] ) :
pred , grad , inp , a u g i n p = modelRINN ( x , params [ ’ 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 red [ : −1 , : ] , a u g i n p [ 1 : , : ] ) l o s s 2 = c r i t e r i o n ( p red [ : −1 , 1 : ] , a u g i n p [ 1 : , 1 : ] )
l o s s =10 .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 . backward ( r e t a i n g r a p h =True ) o p t i m i z e r . s t e p ( )

Review summaries:
The paper addresses the problem of data-driven identification of latent representations for partially-observed dynamical systems. The proposed method uses an augmented state-space model where the dynamical operator is parameterized using neural networks. Given some training data, the proposed approach jointly estimates the parameters of the model as well as the corresponding unobserved components of the latent space for the training sequences. Experimental results on forecasting using real and synthetic datasets are presented.