Data assimilation as a deep learning tool to infer ODE representations of dynamical models

Recent progress in machine learning has shown how to forecast and, to some extent, learn the dynamics of a model from its output, resorting in particular to neural networks and deep learning techniques. We will show how the same goal can be directly achieved using data assimilation techniques without leveraging on machine learning software libraries, with a view to high-dimensional models. The dynamics of a model are learned from its observation and an ordinary differential equation (ODE) representation of this model is inferred using a recursive nonlinear regression. Because the method is embedded in a 5 Bayesian data assimilation framework, it can learn from partial and noisy observations of a state trajectory of the physical model. Moreover, a space-wise local representation of the ODE system is introduced and is key to cope with high-dimensional models. It has recently been suggested that neural network architectures could be interpreted as dynamical systems. Reciprocally, we show that our ODE representations are reminiscent of deep learning architectures. Furthermore, numerical analysis considera10 tions on stability shed light on the assets and limitations of the method. The method is illustrated on several chaotic discrete and continuous models of various dimensions, with or without noisy observations, with the goal to identify or improve the model dynamics, build a surrogate or reduced model, or produce forecasts from mere observations of the physical model. Copyright statement. 15


Data assimilation and model error
Data assimilation aims at estimating the state of a physical system from observations and a numerical dynamical model of that system.It has been successfully applied to numerical weather and ocean prediction, where it often consisted in estimating the initial conditions for the state trajectory of chaotic geofluids (Kalnay, 2002;Asch et al., 2016;Carrassi et al., 2018).
This objective is mostly impeded by the deficiencies of the numerical model (discretrisation, approximate physical schemes, unresolved scales and their uncertain parametrisations, e.g., Harlim, 2017) and the difficulty to match numerical representations of the system with the observations (representation error, Janjić et al., 2018).As a result, the quality of numerical weather predictions based on a methodologically sound data assimilation method crucially depends on both the sensitivity to the initial condition due to the chaotic unstable dynamics and on model error (Magnusson and Källén, 2013).
Model errors can take many forms and accounting for them depends on the chosen data assimilation scheme.A first class of solutions rely on parametrising model error by, for instance, transforming the problem into a physical parameter estimation problem (e.g., Bocquet, 2012;Aster et al., 2013).Other solutions are based on a weakly parametrised form of model error, for instance when it is assumed to be additive noise.Such techniques have been developed for variational data assimilation (e.g., Trémolet, 2006;Carrassi and Vannitsem, 2010), for ensemble Kalman filters and smoothers (e.g., Ruiz et al., 2013;Raanes et al., 2015) and ensemble-variational assimilation (Amezcua et al., 2017;Sakov et al., 2018).In the weakly parametrised form, these techniques should be completed by an estimation of the model error statistics (e.g., Pulido et al., 2018;Tandeo et al., 2018).Moreover, model error's impact can be mitigated by multiplicative and additive inflation (e.g., Whitaker and Hamill, 2012;Grudzien et al., 2018;Raanes et al., 2019) and by physically-driven stochastic perturbations of model simulations in ensemble approaches (e.g., Buizza et al., 1999), or by stochastic subgrid parametrisations (e.g., Resseguier et al., 2017).This account is very far from exhaustive as this is a vast, multiform and very active subject of research.
These approaches essentially seek to correct, calibrate or improve an existing model using observations.Hence, they all primarily make use of data assimilation techniques.

Data-driven forecast of a physical system
Another approach is to renounce to physically-based numerical models of the phenomenon of interest and instead only use observations of that system.Given the huge required datasets, this may seem a far-reaching goal for operational weather and ocean forecasting systems, but recent progress in data-driven methods and convincing applications to geophysical problems of small to intermediate complexity are strong incentives to investigate this bolder approach.Eventually, the perspective of putting numerical models away has a strong practical appeal, even though such perspective may generate intense debates.
For instance, forecasting of a physical system can be done by looking up at past situations and patterns using the techniques of analogues, which can be combined with present observations using data assimilation (Lguensat et al., 2017).It can rely on a representation of the physical system based on diffusion maps that look for a spectral representation of the dataset (see chapter 6 of Harlim, 2018).An original data-driven stochastic modelling approach has been developed by Kondrashov et al. (2015).The method, recently extended to deal with multi-scale datasets (Kondrashov and Chrekroun, 2017), has been applied to successfully estimate reduced models of geophysical phenomena (see e.g., Kondrashov et al., 2018, and references therein).
A fourth route relies on neural networks and deep learning to represent the hidden model and make forecasts from this representation.Examples of such approach applied to the forecasting of low-order chaotic geophysical models are: Park and Zhu (1994) who use a bilinear recurrent neural network and applied it on the three-variable Lorenz model (Lorenz, 1963, hereafter L63), Pathak et al. (2017Pathak et al. ( , 2018) ) who use reservoir network techniques on the L63 model and on the Kuramoto-Sivashinski model (Kuramoto and Tsuzuki, 1976;Sivashinsky, 1977, hereafter KS); Dueben and Bauer (2018) who use a neural network on a low-order Lorenz three-scale model, and on coarse 2D geopotential height maps at 500 hPa.The last two contributions have to resort to local reservoir networks or convolutional layers, respectively, to account for the dimensionality of the models.However, all these representations are not mechanistic and the neural network becomes a surrogate for the hidden model.This marks a key distinction with respect to our approach where the dynamics to be determined are explicitly formulated, as will be clarified later.

Learning the dynamics of a model from its output
Data-driven techniques that seek to represent the model in a more explicit manner, and therefore with a greater interpretability, may use specific classes of nonlinear regression as advocated by Paduart et al. (2010); Brunton et al. (2016).With a view to forecast dynamical systems, it is possible to design neural networks in order to reflect the iterative scheme of a Runge-Kutta (RK) integration scheme.Wang and Lin (1998) proposed and achieved such a goal using classical activation functions, which may however blur the interpretation of the underlying dynamics.Fablet et al. (2018) went further and used a bilinear residual neural network structured to mimic a fourth-order Runge-Kutta scheme (RK4) and noise-free data.Using the Keras tool with the TensorFlow backend, their approach proved to be a very effective tool for the L63 model and to a lesser extent to the 40-variable Lorenz model (Lorenz and Emanuel, 1998, hereafter L96).In particular, they retrieved the parameters of the L63 equations to a high precision.Long et al. (2018) sought the operators of the partial differential equations (PDEs) of a physical system by identifying differentiations with convolution operators of a feed-forward neural network.They successfully applied their method to advection-diffusion problems.Note that, as opposed to our proposal, none of the aforementioned techniques are embedded in a Bayesian framework, making them less suitable to work with noisy and partial data.

Goal and outline
From this point on, the physical system under scrutiny will be called the reference model.It will be assumed to be known only from observations.We follow a data-driven approach inspired by the work of Paduart et al. (2010); Fablet et al. (2018) in the sense that we will consider an observed physical reference model, which might be generated by a hidden mathematical model.This work is focused on either one or a combination of the following goals: (i) to build a surrogate model for the dynamics, (ii) to produce forecasts that emulate those of the reference model, and (iii) to identify the underlying dynamics of the reference model given by a mathematical model.The reference model could be totally unknown or only partially specified.
To achieve these goals, we introduce a surrogate model defined by a set of ordinary differential equations (ODEs): where x ∈ R Nx is the state vector, and x → φ(x) is a vector field that we shall call flow rate.For the sake of simplicity, the dynamics in this paper are supposed to be autonomous, i.e. do not explicitly depend on time.Our technique seeks a fit φ given observations of the reference model.This is a rather general representation since, for instance, PDEs can be discretised into ODEs.We will restrict ourselves to the case where φ is at most quadratic in {x n } 0≤n<Nx .The numerical integration of Eq. ( 1) could be based on any RK scheme, but should additionally rely on the composition of such integration steps.As a result, quite general resolvents of Eq. ( 1) can be built (the resolvent is the model, i.e. the flow rate, integrated in time over a finite time interval).
Importantly, we will not require any machine learning software tool since the adjoint of the model resolvent can be derived without a lot of effort.As opposed to the contributions mentioned in the previous subsections, we embed the technique in a data assimilation framework.From a data assimilation standpoint, the technique can be seen as meant to deal with model error (with or without some prior on the model) and it naturally allows to use partial and noisy observations.We will also build representations of the dynamics that are either invariant by spatial translation (homogeneous) and/or local (i.e. the flow rate of a variable x n only depends on neighbouring variables whose perimeter is defined by a stencil).These properties make our technique scalable and thus potentially applicable to high-dimensional systems.
In Sect.2, we present model identification as a Bayesian data assimilation problem.We first choose an ODE representation of the dynamics, introduce a nonlinear regressor basis, and define integration schemes we will work with.We describe the local and homogeneous representations as physically-based simplifications of the most general case, and we derive the gradient of the problem's cost function based on these representations.We then introduce the Bayesian problem and the resulting cost function used for joint supervised learning of the optimal representation and estimation of the state trajectory.The latter is the standard goal of data assimilation while the former is that of machine learning.Our approach blends them together using the formalism of data assimilation.
In Sect.3, we discuss several theoretical issues: the prior of the model, the convergence of the training step, the connection with numerical analysis of integration schemes, the connection with deep learning architectures, and, finally, the pros and cons of our approach.
In Sect.4, we illustrate the method with several low-order chaotic models (L63, L96, KS and a two-scale Lorenz model) of various sizes, from a a perfectly identifiable model, i.e.where the model used to generate the dataset can be retrieved exactly, to the design of a reduced-order model where the model used to generate the dataset cannot be retrieved exactly, using full or partial, noiseless or noisy observations.Conclusions are given in Sect. 5.
2 Model identification as a data assimilation problem

Ordinary differential equation representation
Our surrogate model is chosen to be represented by an ODE system as described by Eq. (1).We additionally assume that the flow rate can be written where A ∈ R Nx×Np is a matrix of real coefficients to be estimated and r : R Nx → R Np is a map that defines regressor functions of x.R Np is the latent space of the regressors in which the flow rate is linear.
In the absence of any peculiar symmetry, we choose this map to list all the monomials up to second-order built on x, i.e the constant, linear and bilinear monomials.Let us call D = {0, 1, . . ., N x − 1} the set of all variable indices, and P the set of all pairs of variable indices.We introduce the augmented state vector extend D to D = D ∪ {N x }, and define P as the distinct pairs of variable indices in D.
As a result, the regressors are compactly defined by where the scalars in the bracket are the entries of the vector r(x).We count regressors, i.e. the cardinal of P. For instance, a model with 3 variables, x 0 , x 1 , x 3 , such as L63, has 10 such regressors: Higher-order regressors could be included, or of different functional form as in Brunton et al. (2016).However, it is important to keep in mind that we do not seek an expansion of the resolvent of the reference model but of the flow rate φ A .Furthermore, higher-order products of the state variables are later generated by integration schemes and their composition.It is worth mentioning that nonlinear regressions are not widespread in geophysical data assimilation.We are nonetheless aware of at least one noticeable exception that extends traditional Gaussian-based methods (Hodyss, 2011(Hodyss, , 2012)).

Local and homogeneous representations
At least two useful simplifications for the ODEs could be exploited if the state x is assumed to be the discretisation of a spatial field.

Locality
First, we use a locality assumption based on the physical locality of the system: all multivariate monomials in the ODEs have variables x n that belong to a stencil, i.e. a local arrangement of grid points around a given node.This can significantly reduce the number of bilinear monomials in r(x).We assume that s n is the stencil around node n, the pattern being the same for all nodes, except the last one.For the node N x corresponding to the extra variable xNx = 1, we assume that its stencil consists of all the N x + 1 nodes.We then define P s ⊂ P as the sub-set of all pairs (n, m) of variables for which m ∈ s n .The set of required monomials can therefore be reduced to Under these conditions, A becomes sparse.Indeed, for each node n, we assume that ẋn , the time derivative of x n , is impacted only by linear terms x m such that m ∈ s n and quadratic terms x m x l such that m ∈ s n , l ∈ s n and m ∈ s l .However, to keep a dense matrix, we choose to compactly redefine and shrink A by eliminating all a priori zero entries due to the locality assumption.The number of columns of A is then significantly reduced from N p to N a .As a consequence of this redefinition of A, the matrix multiplication inbetween A and r(x) must be changed accordingly.Nonetheless, the operation that assigns coefficients in A to the monomials in r(x) remain linear and we write it as Let us take the example of a one-dimensional extended space as used in Sect. 4. The domain is supposed to be periodic (circle) and the nodes are indexed by 0 ≤ n < N x .The node of index N x is associated to the extra {1}.For 0 ≤ n < N x , the stencil s n is defined as the set of The stencil s Nx consists of all the nodes, i.e.D. We assume 2L + 1 ≤ N x .In that case r(x) as defined by Eq. ( 7) has For instance, there are 161 such regressors for a 40−variable model defined on a circular domain, such as L96, with L = 2: [1, 39 , x 39 x 0 , x 39 x 1 . (9) The row [A] n of the dense A contains the following coefficients for each 0 ≤ n < N x .First there are 2L + 2 regressors built with {1} (the constant and linear regressors).Second, we consider the square monomials x 2 m with m ∈ s n , i.e. x 2 m n−L≤m≤n+L whose number is 2L + 1.Then we consider those separated by one space step: {x m x m+1 } n−L≤m≤n+L−1 whose number is 2L, followed by those separated by two space steps whose number is 2L − 1, and so on until a separation of L is reached.Quadratic monomials of greater separation are discarded since they do not belong to a common stencil as per the above definition reflecting the locality assumption.Hence there is a total of N a = 2L+2 l=L+1 l = 3 2 (L + 1)(L + 2) coefficients per grid-cell.
In appendix A, we show in the one-dimensional space case how to compute the reduced form of the product between A and r(x), assuming locality.This type of technical parametrisation is required for a parsimonious representation of the control variables, i.e. the coefficients of A, and is key for a successful implementation with high-dimensional models.
Note that this locality assumption is hardly restrictive.Indeed, owing to the absence of long-range instantaneous interactions (which are precluded in geophysical fluids), farther distance correlations between state variables can be generated by small stencils in the definition of φ A through time integrations.This would not prevent potential specific long-distance dependencies (such as teleconnections).

Homogeneity
Furthermore, a symmetry hypothesis could optionally be used by assuming translational invariance of the ODEs, called the homogeneity assumption in the following.Because our control parameters, i.e. the coefficients of A parametrise the flow rate, the symmetry simply translates into the rows [A] n of the dense A being the same for all n.Hence A simply becomes a vector in R Na .
Let us enumerate its coefficients in the case of the L96 model with L = 2 and assuming both locality and homogeneity.The coefficients are partitioned into: A (0) for the bias, A l for the linear sector, and A (2) l,m for the bilinear sector.In the linear sector, l = −2, . . ., 2 is the relative position with respect to the current grid point.In the bilinear sector, l, m are the relative positions with respect to the current grid point of the two variables in the product.Proceeding in the same way we counted them, the Note that, while both constraints, locality and homogeneity, apply to the ODEs, they do not apply to the states per se.For instance ODEs for discretised homogeneous two-dimensional turbulence satisfy both constraints and yet hardly generate trivial flows.
For realistic geofluids, the forcing fields (solar irradiance, bathymetry, boundary conditions, friction, etc.) are heterogeneous, so that the homogeneity assumption should be dropped.Nonetheless, the fluid dynamics part of the model would remain homogeneous.As a result, an hybrid approach could be enforced.

Integration scheme and cycling
The reference model will be observed at time steps t k , indexed by integer 0 ≤ k ≤ K. Hence, we need to be able to express the resolvent of the surrogate model from t k to t k+1 .We assume that t k+1 − t k is a multiple of the integration time step of the surrogate model: , where h is the integration time step and N k c is the number of integrations.The time steps t k+1 − t k can be uneven, which is reflected in the dependence of N k c on k.Hence, the resolvent of the surrogate model from t k to t k+1 can be written as i.e. the integration of Eq. ( 1) from t k to t k+1 using the representation Eq. ( 2).
x 0 Representation of the data assimilation system as a hidden Markov chain model, and of the model resolvents F k A and f A .
We define intermediate state vectors inbetween [t k , t k+1 ]: x k,l is the state vector defined at time Figure 1 is a schematic of the composition of the integration steps, along with the state vectors x k and x k,l .
The operator f A is meant to be an explicit numerical integration scheme.In the following, we shall consider a Runge-Kutta (RK) scheme applied to x ≡ x k,l .It has N RK steps.This number of steps coincides with the accuracy of the schemes that we will consider: first-order for the Euler scheme, second-order for RK2, and fourth-order for RK4.Provided the dynamics are autonomous, a general RK scheme reads where the coefficients β i and α i,j entirely specify the scheme and h = (t k+1 − t k )/N k c .Note that α i,j are zero for j ≥ i, so that Eq. ( 12b) can be computed iteratively from k 0 to k NRK−1 , followed by the sum Eq. (12a) to get f A (x).
In the following, h will be absorbed in the definition of A and hence φ A , so that we can take h = 1 without loss of generality.

Bayesian analysis
We consider a sequence of observation vectors y k ∈ R N k y of the physical system at t k indexed by 0 ≤ k ≤ K.The system state is observed through where H k is the observation operator at time t k .The observation error k will be assumed Gaussian of zero mean and covariance matrix R k , It is also assumed to be white in time.The flow rate φ A is given by the approximation Eq. ( 2), so that the resolvent F A of the surrogate model should also be considered as an approximation of the reference model's resolvent.Hence, we generalise Eq. ( 11) to where η k are unbiased Gaussian errors of covariance matrices Q k , also supposed to be white in time.Note that, in all generality, the state space of the surrogate model does not have to match that of the reference model.We will nonetheless take them to coincide here merely for simplicity.
With the goal to identify a model, or build a surrogate of the reference one, we are interested in estimating the probability density function (pdf) p(A|y 0:K ) where y 0:K stands for all observations in the window [t 0 , t K ].To obtain a tractable expression for this conditional likelihood, we need to marginalise over the state variables x 0:K within the window: An approximate maximum a posteriori for A could be obtained by using the Laplace approximation of this integral, and hence getting the maximum of Nonetheless, maximising Eq. ( 16) rigorously yields the maximum a posteriori of the joint variables A, x 0:K .The cost function associated to this joint pdf is by definition J (A, x 0:K ) = − ln(p(A, x 0:K |y 0:K )).Under the Markovian assumption of Eq. ( 14) and the Gaussian form of both model and observational errors, the cost function reads up to a constant depending on Q 1:K and R 0:K only.The vector norm z P is defined as √ z T Pz.This is the cost function of a weak constraint 4D-Var (see Sect. 2.4.3.2 of Asch et al., 2016) with A and x 0:K as control variables.
In the case where the reference model is fully observed, i.e.H k ≡ I x and in the absence of observation noise, i.e.R k ≡ 0, we have x 0:K ≡ y 0:K and the cost function simplifies to where y 0:K is the fully and perfectly observed state trajectory of the reference model.This is similar to a traditional leastsquare function used in machine and deep learning regression.Reciprocally, when the aforementioned hypotheses of noiseless and complete observations do not hold, Eq. ( 17) can be seen as a natural data assimilation extension of Eq. ( 18).Note that Eq. ( 18) only depends on the sequence Q 1:K .If, in addition, the dependence on the prior p(y 0 , A) is neglected in Eq. ( 18), then the maximum a posteriori should not depend on a global rescaling of Q 1:K .This connection between machine learning and data assimilation cost functions had been previously put forward by Hsieh and Tang (1998) and Abarbanel et al. (2018) although in a different form.
The data assimilation system is represented in Fig. 1 as a hidden Markov chain model.This Bayesian view highlights the choice that must be made for R 0:K and/or Q 1:K and provides an interpretation in terms of errors.Furthermore, one could implement an objective estimation of these error statistics as in Pulido et al. (2018).

Gradients and adjoint of the representation
To efficiently minimise the cost function Eq. ( 17) with a gradient-based optimisation tool, we need to derive the gradient of Eq. ( 17) with respect to both A and x 0:K .As for x 0:K , we have: where As for A, we have: assuming x 0 is independent from A. Hence, a key technical aspect of the problem is to compute the tangent linear and adjoint operators required by these gradients.In this paper, we assume that the adjoint of the tangent linear operators of the observation T are known, for instance if the latter are linear as in Sect. 4.
The computations of the gradients and the required adjoints are developed in appendix B. These technicalities can be skipped since they are not required to understand the method.Nonetheless, they are critical to its numerical efficiency.
3 Discussion of the theoretical points In this section, we discuss about the prior pdf p(A), the optimisation of the cost function J (A, x 0:K ), and the connections with deep learning techniques.

Prior information on the reference model
The goal is either to reconstruct an ODE for the reference model, characterised by the coefficients A through φ A , or to build a surrogate model of it.The estimation of x 0:K is then accessory, even though practically critical to the estimation of A. The alternative would have been to consider the estimation of x 0:K as the primary problem, under model error of a prescribed ODE form, the estimation of A being accessory.In both cases, but mostly in this latter case, one may benefit from an informative prior pdf p(A).
The prior pdf p(A) can be used to encode any prior knowledge on the reference model, such as pieces of it that would be known.Indeed, p(A) can formally quantity the uncertainty associated to any part of the surrogate model.For instance, assume that the reference model is partially identifiable, which means that part of the reference model can be represented by an up to bilinear flow-rate of the form Eq. ( 2) and Eq. ( 4).Moreover, assume that there is one such part of the reference model which is known, i.e. that elements of A are actually known while others need to be estimated.Then, the known coefficients can formally be encoded in p(A) with Dirac factors.In practice it could be implemented as a constrained optimisation problem, for instance using an augmented Lagrangian, in order to avoid significantly altering the gradients with respect to A. More generally, assigning a non-trivial prior likelihood, such as a Gaussian one for A, is certainly appealing, but may not be practical.

Numerical optimisation: issues and solutions
The success of the optimisation of J (A, x 0:K ) depends above all on the ability to evaluate it robustly.In particular, it depends on the stability of the numerical integration scheme x = f A (x).In this paper, we chose to rely on one-step explicit schemes which are much simpler to describe and efficient to integrate (a family to which the RK schemes of any N RK belong).These schemes are 0-stable, which means that the finite time error growth goes to zero as the integration step goes to zero.But, as a major drawback, they have a limited absolute (or A-)stability domain (see e.g., chapters 5 and 6 of Gautschi, 2012).For a given state trajectory, there exists a stability domain D s ∈ R Nx×Np out of which the evaluation of J (A, x 0:K ) is hazardous.A failure to evaluate the cost function depends on the number of integration steps.Instabilities are likely to increase exponentially with N k c .This tells that, in the absence of a strong prior p(A), it is safer to start with A = 0 likely to lie close to D s .Alternatively, if stability constraints are known about A, they could be encoded in p(A).It also tells that we should strike an empirical compromise between the composition numbers N k c and the easiness to evaluate J (A, x 0:K ).In the one hand, the larger N k c the more the iterates of A in the optimisation must be kept confined in D s .On the other hand, the longer N k c the broader the class of achievable resolvents, and hence the ability to build a good surrogate.Moreover, the higher the stability of the integration scheme, the larger D s , and the easier the optimisation in spite of an increase in its numerical cost.
As for the sensitiviy on K, the longer the time window, the more observations are available to constrain the problem.
However, the longer K the higher the chances of having a significant instability: the chances of a successful integration typically decrease exponentially with the length K.
This stability issue can be somehow alleviated by normalising the observations y k by their mean and variance in order to avoid excessively large value ranges of the regressors.This will not change the fundamental stability of the schemes, yet may delay the occurrences of instabilities due to the nonlinear terms.
Moreover, instabilities can significantly be mitigated by replacing the monomials with smoothed or truncated ones: One can for instance choose ζ(x) = λ tanh(x/λ), in order to cut-off too large values of |x| and hence delay the growth of instabilities.The parameter λ > 0 is roughly chosen as the typical maximum amplitude of |x| as approximately inferred from the observations.If tanh is deemed to be numerically too costly, one can choose instead This latter change of variables is the one implemented for all numerical applications described in Sect.4, together with the normalisation.From these experiments, we learnt that these tricks often turned critical in the first iterates of the optimisation as the estimate of A progressively migrates to the A-stability domain.After a few iterations, the integrations are stabilised and the nonlinear regime of the truncations in Eq. ( 21) is not tapped into anymore.

Connection and analogies with deep learning architectures
It has recently been advocated that residual deep learning architectures of neural networks can roughly be interpreted as dynamical systems (e.g., E, 2017; Chang et al., 2018).Each layer of the network contributes marginally to the output, so that there exists an asymptotic continuum limit representation of the neural network.Furthermore, as mentioned in the introduction, Wang and Lin (1998); Fablet et al. (2018) have shown that the architecture of the network can follow that of an integration scheme.
By contrast, we have started here from a pure dynamical system standpoint, and proposed to use data assimilation techniques.
In order to explore complex model resolvents, applied to each interval [t k , t k+1 ] between observations, we need the composition of N c integration steps.In particular, this allows for the resolvent to exhibit more realistic long-range correlations.Even when using a reasonably small stencil, long-correlations will arise as a result of the integration steps.Nonetheless the stencil might not be too small so as to model discretised higher-order differential operators.As noted by Abarbanel et al. (2018), each application of f A could be seen as a layer of the neural network.Moreover, within such layer, there are sublayers corresponding to the steps of the integration scheme.The larger N c is the deeper this network is and the richer the class of resolvents to optimise on.
Following this analogy, the analysis step where J (A, x 0:K ) is optimised can be called the training phase.Backpropagation in the network, as coined in machine learning (Goodfellow et al., 2016), corresponds to the computation of the gradient of the network with respect to A and of the model adjoint derived in Sect. 2. This a shortcut for the use of machine learning software libraries such as TensorFlow or PyTorch (see appendix C for a brief discussion).
Because of our complete control of the backpropagation, we hope for a gain in efficiency.However, our method does not have the flexibility of deep learning through established tools.For instance, adding extra parameters, adaptive batch normalisation, dropouts are not granted in our approach without further considerations.
Convolutional layers play the role of localisation in neural architecture.In our approach this role is played by the locality assumption and its stencil prescription.Recall that a tight stencil does not prevent longer-range correlation that are built up through the integration scheme and their composition.This is similar to stacking several convolutional layers to account for multiple scales from the reference model which the neural network is meant to learn from.
Finally, we note that, as opposed to most practical deep learning strategies with a huge amount of weights to estimate, we have reduced the number of control variables (i.e.A) as much as possible.

Models setup and forecast skill
In this section, we shall consider four low-order chaotic models defined on a physical one-dimensional space, except for L63 which is 0-dimensional.They will serve as reference model: 1. the L63 model as defined by the ODEs: with the canonical values (σ, ρ, β) = (10, 28, 8/3).Its Lyapunov time1 is about 1.10.This model is introduced for benchmarking with e.g., Fablet et al. (2018).It is integrated relying on an RK4 scheme with δt r = 0.01 as the integration time-step.
2. the L96 model as defined by ODEs defined over a periodic domain of variables indexed by n = 0, . . ., N x − 1 where where , and F = 8.This model is an idealised representation of a one-dimensional latitude band of the Earth atmosphere.Its Lyapunov time is 0.60.It is integrated with an RK4 scheme and δt r = 0.05.
3. the KS model, as defined by the partial differential equations: over the periodic domain α ∈ [0, 32π] on which we apply a spectral decomposition with N x = 128 modes.The Lyapunov time of our KS model is 10.2 time units.This model is of interest because, even though it has dynamical properties comparable to that of L96, it is much steeper so that much more stringent numerical integration schemes are required to efficiently integrate it.It is defined by a PDE, not an ODE system.It is integrated using the EDTRK4 scheme (Kassam and Trefethen, 2005) and δt r = 0.05.
4. the two-scale Lorenz model (L05III, Lorenz, 2005) is given by the two-scale ODEs: for n = 0, . . ., N x − 1 with N x = 36 and m = 0, . . ., N u − 1 with N u = 360.The indices apply periodically over their domain; m/10 is the integer division of m by 10.We use the original values for the parameters: c = 10 for the time-scale ratio, b = 10 for the space-scale ratio, h = 1 for the coupling, and F = 10 for the forcing.When uncoupled (h = 0), the Lyapunov time of the slow variables x sector of the model Eq. ( 25) is 0.72, which will be the key time scale when focusing on the slow variables (see e.g., Carlu et al., 2019).
This model is of interest because the variable u is meant to represent unresolved scales and hence model error when only considering the slow variables x.For this reason, it has been used in data assimilation papers focusing on estimating model error (e.g., Mitchell and Carrassi, 2015;Pulido et al., 2018).It is integrated with an RK4 scheme and δt r = 0.005 since it is steeper than the L96 model.
The numerical experiments consist of three main steps.First, the truth is generated; i.e. a trajectory of the reference model is computed.The reference model equations are supposed to be unknown, but the trajectory is observed through Eq. ( 13) to generate the observation vector sequence y 0:K .
Next, estimators of the ODEs model and state trajectory x 0:K are learned by minimising the cost function J (A, x 0:K ).
We choose to minimise it using an implementation of the quasi-Newton BFGS algorithm (Byrd et al., 1995), which critically relies on the gradients obtained in Sect. 2. The default choices for the initial ODEs model are A = 0 and x 0:K defined as the space-wise linear interpolation of y 0:K .Note that the minimisation could converge to a local minimum, which may or may not yield satisfactory estimates.
Finally, we can make forecasts using the tentative optimal ODEs model A obtained from the minimisation.With a view to compare it to the reference model used to generate the data, we will consider a set of forecasts with (approximately) independent initial conditions.Both the reference model and the surrogate one will be forecasted from these initial conditions.
The departure from their trajectories, as measured by a root mean square error (RMSE) over the observed variables, will be computed for several forecast lead times.The RMSE is then averaged over all the initial conditions.We will also display the state trajectories of the reference and surrogate models starting from one of the initial conditions.The integration time step of the surrogate model within this training time window [t 0 , t K ] is δt a .It is assumed to be an integer divisor of the training time step ∆t = t k+1 −t k supposed to be constant, i.e. ∆t/δt a is a constant integer, the number of compositions N c , and that is why the index k on N k c has been dropped.The integration time step of the surrogate model within the forecast time window [T, T + T f ] is denoted δt f .Note that δt a and δt f can be distinct.They are critical to the stability of the training and the forecast step, respectively.
The three steps of the numerical experiments are depicted in Fig. 2. Except when explicitly mentioned, the prior p(A) is disregarded, which means that no explicit regularisation on A is introduced.

Inferring the dynamics from dense and noiseless observations of perfectly identifiable models
In the first couple of experiments, we consider a densely observed2 reference model with noiseless observations.In this case, R k ≡ 0 and a uniform rescaling of the Q k , all chosen to be I x , is irrelevant, assuming the prior p(y 0 , A) can be neglected which is hypothesised here and is generally true for large K.Moreover, we use the same numerical scheme and the same integration time step to generate the reference model trajectory as the one used by the surrogate model.In principle, we should be able to retrieve the reference model, since the reference is identifiable among all the possible surrogate models.
Let us first experiment with the L63 model, using an RK4 integration scheme, with ∆t = 0.01 and K = 10 4 (this corresponds to about 91 Lyapunov times).We have N x = 3 and N p = 10.We choose δt a = δt f = 0.01.A convergence to the highest possible precision is achieved after about 120 iterations.The cost function value reaches 0 to machine precision at A .The estimated A is given by A a = A /δt a , because, as mentioned above, the optimised A matrix absorbs the time step.The accuracy of A a is measured by the uniform norm A a − A r ∞ , i.e. the absolute values of the entries of the difference A a −A r , where A r is the matrix of the flow rate of L63 (including the zero coefficients).We obtain A a − A r ∞ = 8.46 × 10 −18 .To compute the RMSE as a function of the forecast lead time, we average over N e = 10 3 runs (each one starting from a difference initial condition).The RMSE (not shown) starts significantly diverging from 0 after 16 Lyapunov time units and reaches a saturation for a lead time of 23.
A similar experiment is carried out with the L96 model, using an RK4 integration scheme, with ∆t = 0.05 and K = 50 (this corresponds to about 4.2 Lyapunov times).We choose here to implement the locality and homogeneity assumptions.The stencil has a width of 5 (i.e. the local grid points with 2 points on its left and 2 points on its right).We have N x = 40, N p = 161 and N a = 18.We choose δt a = δt f = 0.05.Through the minimisation, the main coefficients of the L96 model (forcing F , advections terms, dissipation) are retrieved with a precision of a least 8.88 × 10 −15 .
To compute the RMSE as a function of the forecast lead time, we average over N e = 10 3 runs.The RMSE starts significantly diverging from 0 after 12 time units and reaches a saturation for a lead time of 25.

Inferring the dynamics from dense and noiseless observations of a non-identifiable model
In this second couple of experiments, we consider again a densely observed reference model with noiseless observations.The reference model trajectory is generated by the L96 model (N x = 40) integrated with the RK4 scheme, with ∆t = 0.05 and As opposed to the reference model, in these non-identifiable model experiments, the surrogate model is based on the RK2 scheme, with N c compositions.We choose to implement the locality and homogeneity assumptions, with a stencil of width 5.
We have N p = 161 and N a = 18.We choose δt a = δt f = ∆t/N c .In all cases, the convergence is reached within a few dozens of iterations.The error on the coefficients of A (i.e.A a − A r ∞ ) is about 4 × 10 −2 but with the dominant contribution from F .The RMSE as a function of the forecast lead time is computed for N c = 1, 2, 3, 4, 5 and shown in Fig. 3.The error is reduced as N c is increased.But the improvement saturates at about N c = 5.
Figure 4 shows the trajectories of the reference and surrogate models starting from the same initial condition, as well as their difference, as a function of the forecast lead time.Their divergence becomes significant after 4 Lyapunov times and saturates after 8 Lyapunov times.
Next, the reference model trajectory is generated by the KS model (N x = 128) integrated with the ETDRK4 scheme, with ∆t = 0.05 and K = 50 (this corresponds to about 0.25 Lyapunov time).We choose to implement the locality and homogeneity assumptions, with a stencil of width 9.The surrogate model is based on the RK4 scheme, with N c = 2 compositions.Note that in this experiment, the reference and surrogate models and their integration schemes significantly differ.We have N p = 769 and N a = 45.We choose δt a = ∆t/N c and δt f = 10 −3 .The forecast time step δt f is somehow smaller than δt a because the  KS equations are stiff and so will the surrogate model.This emphasises once again that we have learned about the intrinsic flow rate of the reference model, and not a resolvent thereof.Alternatively, we could use a more robust integration scheme than RK4 such as ETDRK4 for the forecast.
Figure 5 shows the trajectories of the reference and surrogate models starting from the same initial condition, as well as their difference, as a function of the forecast lead time, for a stencil of 9. Their divergence becomes significant after 4 Lyapunov times and saturates after 8 Lyapunov times.
To check whether the PDE of the KS model could be retrieved in spite of the differences in the method of integrations and representations, we have computed a Taylor expansion of all monomials in the surrogate ODEs flow rate up to order 4 so as to obtain an approximately PDE equivalent.The coefficients of this PDE (up to order 4 in the expansion) are displayed in Fig. 6 and compared to the coefficients of the reference model's PDE.The match is good and the terms −x∂ α x, −∂ 2 α x and −∂ 4 α x are correctly identified as the dominant ones.Nonetheless, there are three non-negligible coefficients for higher-order terms that may have been generated by the Taylor expansion, or may originate from a degeneracy among the higher-order operators, or may simply be identified with a shortcoming of our specific ODE representation.

Inferring the dynamics from partial and noisy observations
We come back to the L96 model which is densely observed but with noisy observations that are generated using an independently identically distributed normal noise.The surrogate model is based on an RK4 scheme, N c = 1 and a stencil of length 5, which makes the reference model identifiable.In this case, the outcome theoretically depends on the choice for R k and Q k , as Eq. ( 17) is now used instead of Eq. ( 18).For the sake of simplicity, we have chosen them to be of the scalar form R k ≡ σ 2 y I y and Q k ≡ qI x .In this synthetic experiments, σ y is supposed to be known.However, q is not.We only have a qualitative view on the potential mismatch between the reference and the surrogate model.A Gaussian additive noise might not even be the best statistical model for such error.The optimal value of q could be determined using an empirical Bayes approach based on, for instance, the expectation-maximisation technique in order to determine the maximum a posteriori of the conditional density of q (see e.g., Dreano et al., 2017;Pulido et al., 2018).However, this would make us deviate too much from the objective of this work.In the following we have chosen values of q that yield near to optimal skill scores (typically q ∈ [10 −4 , 1]σ 2 y ).Moreover, we have chosen the relatively small K = 50.It must be beneficial to increase K, especially with noisy observations, but this would force us to address issues relative to weak constraint 4D-Var optimisation for long time windows, a topic which is also beyond the scope of this paper.Preliminary results on this topic are nonetheless discussed later in this section.
Figure 7 shows the forecast skill of the surrogate model as a function of the forecast lead time, and for increasing noise in the observations.Even though, in this configuration, the model is identifiable, the reference value A 0 for A may not correspond to a minimum of the cost function.The cost function might have several local minima.As a consequence, there is no guarantee, starting from a non-trivial initial value for A, that the model will be identified.Indeed, as seen in Fig. 7, the forecast skill degrades significantly as the observation error standard deviation is increased.Using the same setup, we have also reduced the number of observations.The observations are regularly spaced and shifted by one grid cell at each observation time step.The initial A in the optimisation remains 0 while the initial state x 0:K is taken as a cubic spline interpolation of the observations over the whole surrogate model grid.
If the observations are noiseless, the reference model is easily retrieved to a high precision down to a density of 1 site over 4. If the observations are noisy, the performance slowly degrades down as the density is decreased down to about 1 site over 4, below which the minimisation, trapped in a deceiving local minimum, yields an improper surrogate model.
We would like to point out that in the case of noiseless observations, the performance depends little on the length of the training window, beyond a relatively short length, typically K = 50.However, in presence of noisy observations, the overall performance improves with longer K, as expected since the information content of the observations linearly increases with the length of the window.
Figure 9 displays the values of the coefficients in A with respect to the minimisation iteration index for the noiseless and fully observed case.As expected, 4 coefficients converge to the value specified by the exact L96 flow rate while the 14 other coefficients collapse to 0. Figure 10 does the same but in the significantly noisy case where σ y = 1 and a significantly longer window K = 5000 (about 417 Lyapunov times).

Inferring reduced dynamics of a multiscale model
In this experiment, we consider the L05III model.With the locality and the homogeneity assumptions, the scalability is typically linear with the size of the system, and we actually consider the 10-fold model where N x = 360 and N u = 3600 to demonstrate that no issues was encountered when scaling up the method.The large scale variable x of the reference model is noiselessly coefficients of the surrogate model as a function of the minimisation iteration number.Note that the index axis is in logarithmic scale.The coefficient of the forcing (F ) is in green, the coefficients of the convective terms are in cyan, and the dampening coefficient is in magenta.
and fully observed over a short training window (K = 50, which corresponds to about 0.35 Lyapunov time) whereas the short scale variable u is not observed.The surrogate model is based on the RK4 scheme and N c = 2 compositions.We choose to implement the locality and homogeneity assumptions, with a stencil of width 5. We have N p = 161 and N a = 18.We choose Figure 11 shows the trajectories of the reference and surrogate models starting from the same initial condition, as well as their difference, as a function of time.
The emergence of error, i.e. the divergence from the reference, appears as long darker stripes on the density plot of the difference (close to zero difference values appear as white or light colour).We argue that these stripes result from the emergence of subscale perturbations that are not properly represented by the surrogate model.Reciprocally there are long lasting stripes of low error not yet impacted by subscale perturbations.As expected and similarly to the L96 model, the perturbations are transported eastward as shown by the upward tilt of the stripes in Fig. 11.Clearly, in this case, a flow rate of the form Eq. ( 2) could be insufficient.Adding a stochastic parametrisation with parameters additionally inferred might offer a solution, as in Pulido et al. (2018).Because of this mixed performance, the RMSE slowly degrades (compared to the other experiments reported so far) with the increase of the forecast lead time (not shown).

Conclusions
We have proposed to infer the dynamics of a reference model from its observation using Bayesian data assimilation, which is a new and original scope for data assimilation.Over a given training time window, the control variables are the state trajectory and the coefficients of an ODE representation for the surrogate model.We have chosen the surrogate model to be the com- position of an explicit integration scheme (Runge-Kutta typically) applied to this ODE representation.Time-invariance, space homogeneity and locality of the dynamics can be enforced making the method suitable for high-dimensional systems.The cost function of the data assimilation problem is minimised using the adjoint of the surrogate resolvent which is explicitly derived.
Analogies between the surrogate resolvent and a deep neural network have been discussed, as well as the impact of stability issues of the reference and surrogate dynamics.
The method has been applied to densely noiseless observed systems with identifiable reference model yielding a perfect reconstruction close to machine precision (L63 and L96 models).It has also been applied to densely or partially observed, identifiable or non-identifiable model with or without noise in the observations (L96 and KS model).For moderate noise and sufficiently dense observation, the method is successful in the sense that the forecast is accurate beyond several Lyapunov times.The method has also been used as a way to infer a reduced model for a multi-scale observed system (L05-III model).
The reduced model was successful in emulating slow dynamics of the reference model but could not properly account for the impact of the fast unresolved scale dynamics on the slow ones.A subgrid parametrisation would be required or would have to be inferred.
Two potential obstacles have been left aside on purpose but should later be addressed.First, the model error statistics have not been estimated.This could be achieved using for instance an empirical Bayesian analysis based on a ensemble-based stochastic expectation maximisation technique.This is an especially interesting problem since the potential discrepancy between the reference and the surrogate dynamics is in general non-trivial.Second, we have used relatively short training time windows.Numerical efficient training on longer windows will likely require using advanced weak constraint variational optimisation techniques.
In this paper, only autonomous dynamics have been considered.We could at least partially extent the method to nonautonomous systems by keeping a static part for the pure dynamics and consider time-dependent forcing fields.We have not numerically explored non-homogeneous dynamics but we have shown how to learn from them using non-homogeneous local representations.
A promising yet challenging path would be to consider implicit or semi-implicit scheme following for instance the idea in Chen et al. (2018).This idea is known in geophysical data assimilation as the continuous adjoint (see e.g., Bocquet, 2012).
This would considerably strengthen the stability of the training and forecast steps at the cost of more intricate mathematical developments.
If observations keep coming after the training time window, then one can perform data assimilation using the ODE surrogate model of the reference model.This data assimilation scheme could only focus on state estimation, or it could continue to update the ODE surrogate model for the forecast.
the coefficients of the RK scheme as defined in Eq. ( 12a), with the state space identity matrix, and where κ is the vector of size N RK N x defined after Eq. (B2).Looking first at the gradient with respect to the state variable, and using Eq.(B3), we have which yields the adjoint operator Let us consider an arbitrary vector d ∈ R Nx ; we have To avoid computing G −T explicitly, let us define the vector z ∈ R NxNRK such that Because G T is upper triangular with diagonal entries of value 1, z is the solution of a linear system easily solvable iteratively, which stands as an adjoint/dual to the RK iterative construction.Hence, we finally obtain a formula and algorithm to evaluate which is key to computing Eq. ( 19).Indeed Eq. ( 19) now reads where z k is the iterative solution of the system Second, let us look at the gradient of J (A, x 0:K ) with respect to A. From Eq. (B4) and Eq.(B3), and now considering variations with respect to A, we obtain: which yields, using z as defined by Eq. (B8), For each 0 ≤ i < N RK , let us introduce r i = r(x i ) ∈ R Np , and let us denote z i ∈ R Nx the subvector of z for the i-th block of the Runge-Kutta scheme.Then, we have for 0 ≤ n < N x and 0 ≤ p < N p : This is key to efficiently computing Eq. ( 20), which now reads where z k is the solution of The index k of G k indicates that the operators defined in the entries of G are evaluated at x k .

B3 Composition of integration steps
We now consider a resolvent which is the composition of N k c ≥ 2 integration steps over [t k , t k+1 ]: where x is an alias to x k .Let us first look at the gradient with respect to the state variable.Within the scope of this section, we define x 0 ≡ x and for 1 ≤ l ≤ N k c : x l ≡ f A (x l−1 ).Hence, x = x Nc .We also define ∇ x l f A to be the tangent linear of f A at x l .By Leibniz rule, we obtain (B15) We can now apply Eq. (B9) to each individual integration step and obtain for any d ∈ R Nx by PyTorch.Our experiments made use of a muti-core CPU and/or a GPU.This purely qualitative ranking is not surprising since (i) PyTorch and TensorFlow excel in tensor algebra operations which are not massively use in the way we built our cost function with relatively few but significant parameters (ii) the time spent in the development of each method scales with their speed.
Author contributions.MB first developed the theory, implemented, run and interpreted the numerical experiments, and wrote the original version of the manuscript.All authors have discussed the theory, the interpretation of the results and edited the manuscript.The four authors approved the manuscript for publication.
Competing interests.The authors declare that they have no conflict of interest.

Figure 2 .
Figure 2. Schematic of the three steps of the experiments, with the associated time steps (see main text).The beginning of the forecast window may or may not coincide with the end of the training window.The lengths of the segments δtr, δta, and δt f are arbitrary in this schematic.

Figure 3 .Figure 4 .
Figure 3. Average RMSE of the surrogate model (L96 with an RK2 structure) compared to the reference model (L96 with an RK4 integration scheme) as a function of the forecast lead time (in Lyapunov time unit) for an increasing number of compositions.

Figure 5 .
Figure 5. Density plot of the KS reference and surrogate model trajectories, as well as their difference trajectory, as a function of the forecast lead time (in Lyapunov time unit).The observations are noiseless and dense; the model is not identifiable.

Figure 6 .
Figure 6.Coefficients of the surrogate PDE model (blue) resulting from the expansion of the surrogate ODEs, and compared to the reference PDE's coefficients (orange).

Figure 7 .
Figure 7. Average RMSE of the surrogate model (L96 with an RK4 structure) compared to the reference model (L96 with an RK4 integration scheme) as a function of the forecast lead time (in Lyapunov time unit) for a range of observation error standard deviation σy.

Figure 8 .Figure 9 .
Figure8.Gap between the surrogate (L96 with an RK4 structure) and the (identifiable) reference dynamics (L96 with an RK4 integration scheme) as a function of the observation error standard deviation σy.Note the use of logarithmic scales.

Figure 10 .
Figure 10.L96 is the reference model, which is fully observed with observation error standard deviation σy = 1: plot of the Na = 18

Figure 11 .
Figure 11.Density plot of the L05III reference and surrogate model trajectories, as well as their difference trajectory, as a function of the forecast lead time (in Lyapunov time unit).The right panel shows a zoom of the difference between times 4 and 5.

(
∇ x l f A ) T d = d + (∇ x l ϕ) T z l , (B16)where z l is the solution ofG T l z l = b T d. (B17) Hence, to compute ∇ x F k A T • d, we define xNc = d and for N c − 1 ≥ l ≥ 0: xl = (∇ x l f A ) T xl+1.This finally reads:xl = xl+1 + (∇ x l ϕ l ) T z l ,(B18)for N c − 1 ≥ l ≥ 0, where z l is the solution ofG T l z l = b T xl+1 .(B19)