Learning Latent Dynamics for Partially-Observed Chaotic Systems

TL;DR

This paper introduces a neural network-based framework for identifying latent dynamical models of partially-observed chaotic systems, improving forecasting accuracy and long-term behavior prediction over existing methods.

Contribution

It presents a novel approach that jointly learns an ODE in latent space and reconstructs states, integrating ideas from Koopman theory and Takens' theorem.

Findings

Outperforms state-of-the-art methods in short-term forecasting

Accurately captures long-term asymptotic patterns

Provides a unified neural network framework for latent dynamics

Abstract

This paper addresses the data-driven identification of latent dynamical representations of partially-observed systems, i.e., dynamical systems for which some components are never observed, with an emphasis on forecasting applications, including long-term asymptotic patterns. Whereas state-of-the-art data-driven approaches rely on delay embeddings and linear decompositions of the underlying operators, we introduce a framework based on the data-driven identification of an augmented state-space model using a neural-network-based representation. For a given training dataset, it amounts to jointly learn an ODE (Ordinary Differential Equation) representation in the latent space and reconstructing latent states. Through numerical experiments, we demonstrate the relevance of the proposed framework w.r.t. state-of-the-art approaches in terms of short-term forecasting performance and long-term behaviour. We further discuss how the proposed framework relates to Koopman operator theory and Takens' embedding theorem.

Equations

$$\{\begin{array}{rl}& {\dot{\mathrm{z}}}_t=\alpha {\mathrm{z}}_t\\ & {\mathrm{z}}_{t_0}={\mathrm{z}}_0\end{array}$$

$${\mathrm{z}}_t=K{e}^{\alpha t}$$

$${{\mathrm{X}}_t}^T=[{\mathrm{x}}_t^T,{\mathrm{y}}_t^T]$$

$$\{\begin{array}{rl}& {\dot{\mathrm{X}}}_t=f_{\theta }({\mathrm{X}}_t)\\ & {\mathrm{x}}_t=G(\mathrm{X}(\mathrm{t}))\end{array}$$

$$\begin{array}[]{c}\hat{\theta}=\arg\displaystyle\min_{\theta}\displaystyle\min_{\{y_{t}\}_{t}}\sum_{t=1}^{T}\|\mathrm{x}_{t}-G\left(\Phi_{\theta,t}\left(X_{t-1})\right)\right)\|^{2}\\ \nobreak\leavevmode\hfil\\ \mbox{Subject to }\left\{\begin{array}[]{lcl}X_{t}&=&\Phi_{\theta,t}(X_{t-1})\\ \nobreak\leavevmode\\ G(X_{t})&=&\mathrm{x}_{t}\\ \end{array}\right.\end{array}$$

$${\Phi }_{\theta ,t}(X_{t-1})=X_{t-1}+{\int }_{t-1}^t{f}_{\theta }(X_u)du$$

$$\begin{gathered} \underset{\theta }{\min }\underset{\{{\mathrm{y}}_t{\}}_t}{\min }\sum _{t=1}^T\Vert {\mathrm{x}}_t-G\left({\Phi }_{\theta ,t}\left(X_{t-1})\right)\right){\Vert }^2+\lambda \Vert X_t-{\Phi }_{\theta ,t}(X_{t-1}){\Vert }^2 \\ \text{leavevmode} \end{gathered}$$

$${\hat{\mathrm{y}}}_T^n=\arg \underset{{\mathrm{y}}_T^n}{\min }\underset{\{{\mathrm{y}}_t^n{\}}_{t<T}}{\min }\sum _{t=T+1}^{T+\delta T}\Vert {\mathrm{x}}_t^n-G\left({\Phi }_{\theta ,t}\left(X_{t-1}^n)\right)\right){\Vert }^2+\lambda \Vert X_t^n-{\Phi }_{\theta ,t}(X_{t-1}^n){\Vert }^2$$

$$\{\begin{array}{rl}& {\dot{\mathrm{z}}}_t=\alpha {\mathrm{z}}_t\\ & {\mathrm{z}}_{t_0}={\mathrm{z}}_0\end{array}$$

$$\psi \left(\{{\mathrm{x}}_t{\}}_{t_0:T}\right)=\arg \underset{{\mathrm{X}}_T}{\min }\underset{\{{\mathrm{X}}_t{\}}_{t<T}}{\min }\sum _{t=1}^T\Vert {\mathrm{x}}_t-G\left({\Phi }_{\theta ,t}\left(X_{t-1})\right)\right){\Vert }^2+\lambda \Vert X_t-{\Phi }_{\theta ,t}(X_{t-1}){\Vert }^2$$

$${\mathrm{x}}_{T+1}=G\left({\Phi }_{\theta ,t}\left(\psi \left(\{{\mathrm{x}}_t{\}}_{t=t_0:T}\right)\right)\right)$$

Full text

1 Introduction

Learning the underlying dynamical representation of some observed variables ${\mathrm{x}}_{t}\in\mathbb{R}^{n}$, where $t\in\left\{t_{0},...,T\right\}$ is the temporal sampling and $n$ the dimension of our observation space 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 [1, 2, 3, 4].

For fully-observed systems, i.e. when observed variables ${\mathrm{x}}_{t}$ relate to some underlying deterministic states ${\mathrm{z}}_{t}$ according to a diffeomorphic mapping possibly corrupted by some noise, recent advances [5, 6, 7] have shown that one can identify the governing equations of the dynamics of ${\mathrm{z}}$ from a representative dataset of observations $\{{\mathrm{x}}_{t_{i}}\}_{i}$. When the observed variables ${\mathrm{x}}_{t}$ only relate to some but not all the components of underlying states ${\mathrm{z}}_{t}$, these approaches do not apply since no ODE or, more generally, no one-to-one 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 [8].

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 non-parametric schemes based on nearest-neighbors or analogs [9] as well as parametric schemes which include polynomial representations [10], neural network models [11], Support Vector Regression (SVR) models [12]. For all these approaches, the identification of the appropriate delay embedding is a critical issue [13, 14].

In this work, we do not rely explicitly on a delay embedding. We rather aim to identify an augmented space of higher dimension than the manifold spanned by the observed variables, where the dynamics of the observations can be fully described by an ODE. Numerically, the joint identification of the augmented space and the governing ODE for a given observation dataset exploits neural-network representations. We report experiments for linear and chaotic dynamics, which illustrate the relevance of the proposed framework compared to state-of-the-art approaches, and then further discuss the key features of this framework with respect to state-of-the-art dynamical systems identification tools such as Koopman operator theory [15].

2 Background and Related Work

Let us consider an unobserved state variable $\boldsymbol{\mathrm{z}}$ governed by an autonomous system of $s$ differential equations $\dot{{\mathrm{z}}}_{t}=f({\mathrm{z}}_{t})$. Let us also assume that this system generates a flow $\Phi_{t_{i}}({\mathrm{z}}_{0})=\int_{t_{0}}^{t_{i}}f({\mathrm{z}}_{u})du\in\mathbb{R}^{s}$ with trajectories that are characterized by an asymptotic limit-cycle $L$ of dimension $d$ contained in $\mathbb{R}^{s}$. We further assume that we are provided with a measurement function $\mathcal{H}$ that maps our state variable $\mathrm{z}$ to our observation space $\mathrm{x}_{t}=\mathcal{H}(\mathrm{z}_{t})\in\mathbb{R}^{n}$.

When considering the data-driven identification of a dynamical mapping that governs some observation data, we first need to evaluate if the dynamics expressed in our observation space are deterministic and thus, can be described using a smooth111The word smooth here stands for continuously differentiable or $\mathcal{C}^{1}$. 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 $\mathbb{R}^{n}$ so one can reliably perform forecasts in the observation space.

The general condition under which a mapping $\mathcal{H}$ will preserve 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 projection $\mathcal{H}$ in $\mathbb{R}^{n}$ if $\mathcal{H}$ is one-to-one and an immersion of $L$ in $\mathbb{R}^{n}$. Under these conditions our observation mapping is called an embedding [16].

The simplest example of an embedding is the case where our observation operator is an identity matrix. With such embedding we have direct access to the state variable $z$ which is governed by a deterministic ODE. This particular case has been widely studied in the literature, parametric representations based on augmented state formulations have been for decades the most popular models due to their simplicity and interpretability [10], [5]. In the last years, neural network and deep learning have enriched the state-of-the-art of parametric representations [17], [18]. In particular, the link between residual networks [19], [20] and numerical integration schemes have opened new research avenues for learning extremely accurate dynamical models even from irregularly-sampled training data.

However, for a wide range of real-world systems, we are never provided with an observation operator that forms an embedding of the latent state variables. The key issue is that we do not have any guarantees on the existence of a smooth ODE governing the temporal evolution of our observations in $\mathbb{R}^{n}$.

From this point of view, the question of finding an appropriate dynamical representation of some observed data may not be this straightforward. The fact that our data may come from some higher-dimensional unobserved governing equation may restrict the above-mentioned state-of-the-art algorithms. The main difficulty is how to go from a low-dimensional observation space to higher-dimensional spaces that provides at least a one-to-one mapping between two successive observations. From a geometrical point of view, the time delay theorem [8] provides a way to build higher dimensional spaces that preserve the topological properties of the true (unobserved) dynamics limit-cycle from a lower-dimensional phase space. A generalization of this theorem [16] shows that we can reconstruct topologically similar limit-cycles using any appropriate smooth composition map of the observations. The derivation of a dynamical system from such representations on the other hand encountered 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 [5] and K-Nearest Neighbors (KNN) algorithms were proposed but they often fail to achieve both accurate short-term forecasting performance and topologically similar reconstructed limit-cycles.

In this work, we address the identification of a latent embedding that fits an ODE representation of the underlying dynamics of our observations. As detailed in the next Section, the core idea of this work is to replace the augmented observable vector formulated in [8] by augmented states that are outputs of an optimization problem with respect to an ODE formulation.

3 Learning latent representations of partially-observed dynamics

Augmented latent dynamics: Let us assume a continuous $s$-dimensional dynamical system $\mathrm{z}_{t}$ governed by an autonomous ODE $\dot{{\mathrm{z}}}_{t}=f({\mathrm{z}}_{t})$ with $\Phi_{t}$ the corresponding flow $\Phi_{t}(z_{t_{0}})=\int_{t_{0}}^{t}f(z_{u})du$. In many applications, one cannot fully access to state $\mathrm{z}$ and the observed variables, also called observables, only relate to some components of state $\mathrm{z}$. Formally, we can define an observation function $\mathcal{H}:\mathbb{R}^{s}\longrightarrow\mathbb{R}^{n}$ such that observed state $\mathrm{x}_{t}$ follows $\mathrm{x}_{t}=\mathcal{H}(\mathrm{z}_{t})$

The key question arising here is the extent to which the dynamics expressed in the observation space $\mathbb{R}^{n}$ reflect the true underlying dynamics in $\mathbb{R}^{s}$, and consequently, the conditions on $\mathcal{H}$ under which the predictable deterministic dynamical behavior of the hidden states is still predictable in the observation space. To illustrate this issue, we may consider a linear dynamical system in the complex domain governed by the following linear ODE:

[TABLE]

with ${\mathrm{z}}\in\mathbb{C}$ a state variable and $\alpha\in\mathbb{C}$ a complex imaginary number. The solution of this problem is

[TABLE]

with $K$ a constant depending on ${\mathrm{z}}_{0}$. Assuming that we are only provided the real part as direct measurements ${\mathrm{x}}_{t}\in\mathbb{R}$ of the true states ${\mathrm{z}}_{t}$, no smooth autonomous ODE model in the scalar observation space can describe the trajectories of the observed states as the mapping between two successive observed states is not one-to-one. For example, assuming that $\mathrm{z}_{t_{0}}$ and $\mathrm{z}_{t_{1}}$ correspond to two states that have the same real part but distinct imaginary parts, the associated observed states are equal $\mathrm{x}_{t_{0}}=\mathrm{x}_{t_{1}}$. However, the time evolution of the states $\mathrm{z}_{t_{0}}$ and $\mathrm{z}_{t_{1}}$ differ if they have different imaginary parts, such that the observed states $\mathrm{x}_{t_{0}+\delta}$ and $\mathrm{x}_{t_{1}+\delta}$ after any time increment $\delta$ are no more equal. As a consequence, a given observation may have more than one future state and this behavior can not be represented by a smooth ODE in the observation space.

For a given observation operator $\mathcal{H}$ of a deterministic underlying dynamical system governing $\mathrm{z}_{t}$, Takens’s theorem guarantees the existence of an augmented space, defined as a delay embedding of observed state series, in which a one-to-one mapping exists between successive time steps of the observation series [8]. Rather than exploring such delay embeddings, we aim to identify an augmented latent space, where latent dynamics are governed by a smooth ODE and can be mapped to observed states. Let us define $\mathrm{X}_{t}\in\mathbb{R}^{d_{E}}$ the $d_{E}$-dimensional augmented latent state with $d_{E}>n$ as follows:

[TABLE]

with $\mathrm{y}_{t}\in\mathbb{R}^{d_{E}-n}$ the unobserved component of latent state $\mathrm{X}_{t}$. We aim to define the augmented latent space according to the following state space model:

[TABLE]

where the dynamical operator $f_{\theta}$ belongs to a family of smooth operators (in order to guarantee uniqueness [21]) parametrized by a parameter vector $\theta$. We typically consider a neutral-network representation with Lipschitz nonlinearities and finite weights. As detailed in the next sections, we address the identification of the operator $f_{\theta}$ and associated latent space $\mathrm{X}$ from a dataset of observed state series $\{\mathrm{x}_{0},\ldots,\mathrm{x}_{T}\}$ as well as the exploitation of identified latent dynamics for the forecasting of the time evolution of the observed states, for instance unobserved future states $\{\mathrm{x}_{T+1},\ldots,\mathrm{x}_{T+N}\}$.

Learning scheme: Given an observation time series $\{\mathrm{x}_{0},\ldots,\mathrm{x}_{T}\}$, we aim to identify the state-space model defined by (4), which amounts to learning the parameters $\theta$ of the dynamical operator $f_{\theta}$. However, as the component $\mathrm{y}_{t}$ of the latent state ${\mathrm{X}_{t}}$ is never observed, this identification requires the joint optimization of the model parameters $\theta$ as well as the hidden component $\mathrm{y}_{t}$. Formally, this problem is then stated as the following minimization of the forecasting error on observed variables:

[TABLE]

with $\Phi_{\theta,t}$ the one-step-ahead diffeomorphic mapping associated with operator $f_{\theta}$ such that:

[TABLE]

In (5), the loss to be minimized involves the one-step-ahead forecasting error for the observed variable $\mathrm{x}_{t}$. The constraints state that the augmented state $\{X_{t}\}_{t}$ is composed of observed component $G(X_{t})$ and should be a solution of ODE (4). Here, we numerically minimize the equivalent formulation:

[TABLE]

where $X_{t}^{T}=[\mathrm{x}_{t}^{T},\mathrm{y}_{t}^{T}]$ and $\lambda$ a weighting parameter. The term $\|X_{t}-\Phi_{\theta,t}(X_{t-1})\|^{2}$ may be regarded as a regularization term such that the inference of the unobserved component $\mathrm{y}_{t}$ of the augmented state $X_{t-1}$ is not solved independently for each time step.

Using a neural-network parametrization for the ODE operator $f_{\theta}$, the corresponding forecasting operator $\Phi_{\theta,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 [22], allows very accurate identification of ODE models [6, 20]. Hence, for a given observed state series $\{x_{0},\ldots,x_{T}\}$, we minimize (6) jointly w.r.t. $\theta$ and unobserved variables $\{\mathrm{y}_{0},\ldots,\mathrm{y}_{T}\}$. In the experiments reported in Section 4, we consider bilinear architectures [6]. The proposed framework applies to any neural-network architecture.

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

Let us denote by $\mathrm{x}^{n}_{t},$ $t\in\left\{t_{0},...,T\right\}$ a new series of observed states. We aim to forecast future states $\mathrm{x}^{n}_{t},$ $t\in\left\{T+1,...,T+\delta T\right\}$. Following (6), we infer the unobserved component $\widehat{\mathrm{y}}_{T}$ of latent state $X^{n}_{T}$ at time $T$ from the minimization of:

[TABLE]

Here, we only minimize w.r.t. latent variables $\{\mathrm{y}^{n}_{t}\}$ given the trained forecasting operator $\Phi_{\theta,t}$. This minimization relates to a variational assimilation issue with partially-observed states and known dynamical and observation operators [23]. Similarly to the learning step, we benefit from the neural-network parameterization of operator $\Phi_{\theta,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 one. 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 $\mathrm{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

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:

[TABLE]

with $\alpha=-0.1-0.5j$, $j^{2}=-1$ and $\mathrm{z}_{0}=0.5$. As $\alpha\in\mathbb{C}$ with $Real(\alpha)<0$ and $\mathrm{z}_{0}\neq 0$, the solution of this ODE is an ellipse in the complex plane (Fig. 1).

We assume that the observed state is the real part of the underlying state, i.e. the observation function $\mathcal{H}:\mathbb{C}\longrightarrow\mathbb{R}$ is given by $\mathrm{x}_{t}=Real(\mathrm{z}_{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 $\mathrm{X}_{t}=[\mathrm{x}_{t},\mathrm{y}^{1}_{t}]$. As neural-network parametrization for operator $f_{\theta}$, 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, 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 [8] 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.

We also illustrate the forecasting performance of the proposed framework in Fig. 1. Given the same initial condition over the observable state, we compare the forecasting of the observed state by the true equation and the proposed model. For the proposed model, the initial condition over the unobservable component results from minimization (7). The resulting forecast of the observation sequence is similar to the true observation sequence (mean square error < 1E-5).

Lorenz-63 dynamics: Lorenz-63 dynamical system is a 3-dimensional model that involves, under some specific parametrizations [24], chaotic dynamics with a strange attractor. We simulate chaotic Lorenz-63 state sequences with the same model parameters as proposed in [24] using the LOSDA ODE solver [25] with an integration step of 0.01. We assume that only the first Lorenz-63 variable is observed $\mathrm{x}_{t}=\mathrm{z}_{t,1}$. We apply the proposed framework to this experimental setting using a training sequence of $10000$ time-steps.

For benchmarking purposes, we perform a quantitative comparison with state-of-the-art approaches using delay embedding representations [8]. The parameters of the delay embedding representation, namely the lag $\tau$ 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 [13], respectively denoted as $\tau_{MI}$ and $\tau_{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 [14]. 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 [26] and the Sparse Regression (SR) method on a second order polynomial representation of the delay embedding states. Regarding deep learning models, we evaluate the Latent-ODE model as proposed in [19] for the modeling of time series. The architecture of the Latent-ODE model is the following: The recognition network is an RNN with 100 hidden units. The dimension of the latent space is 4 and the dynamical model is a fully-connected network with bilinear layers as proposed in [6] (the same as used in our model). The decoder is a fully connected network with one layer and 15 units.

Regarding the proposed framework, referred to as Neural embedding for Dynamical Systems (NbedDyn), we tested the model for different dimensions of the augmented state space, namely from 3 to 6. The neural-network parametrization for operator $f_{\theta}$ is a fully-connected architecture with bilinear layers as proposed in [6] 222The code will be available on GitHub soon..

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 model lead to the poorest performance both in terms of forecasting error and asymptotic behavior. This is mainly due 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 [28]. 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 SR model seems to lead to better short term forecast, however, it does not capture well the chaotic 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 ($\tau_{MI}$ = 6, $d_{E}=3)$, one may retrieve the expected long-term chaotic behavior ($\lambda_{1}=0.84$) with a relatively low short-term forecasting error (1.2e-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-5 for a one-step-ahead forecast when considering $d_{E}=6$, i.e. one order of magnitude compared to the second best model) while keeping an appropriate chaotic long-term pattern ($\lambda_{1}=0.87$).

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 [29]. Sea surface dynamics are chaotic dynamics 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 [30]. The spatial resolution of our data is a $0.05\degree$ 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\degree E$ to $4.25\degree E$, $37.25\degree N$ to $39.5\degree N$). The spatial scales of the data were decomposed through a PCA decomposition of dimension $N_{E}=15$, which amounts to capture $92\%$ of the total variance.

We report forecasting performance for our model and include a comparison with analog methods and a neural ODE setting in Tab. 2. For our model, we consider the following parameterization: a latent space with $d_{E}=60$ and a fully connected bilinear architecture with 60 bilinear neurons and 4 fully-connected layers for the dynamical model $f_{\theta}$. 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%). We let the reader refer to the Supplementary Material for a more detailed analysis of these experiments, including visual comparisons of the forecasts.

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 systems. 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, state-of-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 $\psi$ such that:

[TABLE]

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

[TABLE]

Hence, $\psi$ defines a delay embedding representation implicitly stated through minimization (7). In this embedding, the dynamics of the observed system $\mathrm{x}$ is governed by the composition of observation operator $G$ and forecasting operator $\Phi_{\theta,t}$. Regarding the literature on Koopman operator theory, most approaches rely on the explicit identification of eigenfunctions and eigenvalues of the Koopman operator [15, 31, 32]. Our framework relates to the identification of the infinitesimal generator $f_{\theta}$ 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_{\hat{\theta}}$ is also diagonalizable, such that the identification of infinitesimal generator $f_{\hat{\theta}}$ 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 high-dimensional and stochastic systems. In the considered framework, operator $G$ is stated as an identity operator on the observed component of state $X_{t}$. For high-dimensional systems, operator $G$ may be stated as an encoder to keep the latent space low-dimensional and make minimization (6) tractable. The combination of the proposed framework with the variational setting considered in the Latent-ODE model [19] 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 [33].

Appendix

Appendix A 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_{\theta}$.

Appendix B 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 $9$ Lyapunov times, when the best alternative approach diverge from the true trajectory beyond $5$ 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.

Appendix C Additional figures of the Sea level Anomaly experiment

Forecasted states of the Sea Level Anomaly and their gradients are illustrated in Fig. 7, 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.