Data assimilation as a learning tool to infer ordinary differential equation 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 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 coping with high-dimensional models. It has recently been suggested that neural network architectures could be interpreted as dynamical systems. Recipro-cally, we show that our ODE representations are reminiscent of deep learning architectures. Furthermore, numerical analysis considerations of stability shed light on the assets and limitations of the method.


Data assimilation and model error
Data assimilation aims at estimating the state of a physical system from its observation and a numerical dynamical model for it.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 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 in matching 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 relies 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 for ensemble-variational assimilation (Amezcua et al., 2017;Sakov et al., 2018).In the weakly parametrised form, these methods should be completed by an estimation of the model error statistics (e.g.Pulido et al., 2018;Tandeo et al., 2019).Moreover, a model error's impact can be mitigated, and this is often the case in applications, 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
An alternative is to renounce physically based numerical models of the phenomenon of interest and instead to 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 a perspective may generate intense debates.
For instance, forecasting of a physical system can be done by looking up past situations and patterns using the techniques of analogues, which can be combined with present observations using data assimilation (Lguensat et al., 2017), or 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 an approach applied to the forecasting of low-order chaotic geophysical models are Park and Zhu (1994), who used a bilinear recurrent neural network and applied it to 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), and Dueben and Bauer (2018), who use a neural network on a low-order Lorenz three-scale model and on coarse two-dimensional geopotential height maps at 500 hPa.The last three contributions have to resort to local reservoir networks or convolutional layers, respectively, to cope with 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) andBrunton et al. (2016).With a view to forecasting dynamical systems, it is possible to design neural networks in order to reflect the iterative form 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 so as to mimic a fourth-order RK 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 for 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.As opposed to our proposal, described hereafter, none of the aforementioned techniques is embedded in a Bayesian framework, making them less suitable for working 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 works of Paduart et al. (2010) and 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 or process.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 N x is the state vector and x −→ φ(x) is a vector field that we shall call the flow rate.For the sake of simplicity, the dynamics in this study 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<N x .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 accommodates partial and noisy observations.Moreover, we will 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 the 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 loworder chaotic models (L63, L96, KS, and a two-scale Lorenz model) of various sizes, from a perfectly identifiable model, i.e.where the model used to generate the dataset can be retrieved exactly, to 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 as where A ∈ R N x ×N p is a matrix of real coefficients to be estimated and r : R N x −→ R N p is a map that defines regressor functions of x.R N p 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 three variables, x 0 , x 1 , and x 2 , such as L63, has 10 such regressors: Higher-order regressors, as well as regressors of different functional forms, could be included 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 .As a consequence, higher-order products of the state variables are anyhow generated by the 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 Gaussianbased methods (Hodyss, 2011(Hodyss, , 2012)).
M. Bocquet et al.: Data assimilation as a learning tool to infer dynamics

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 x N x = 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 in between A and r(x) must be changed accordingly.Nonetheless, the operation that assigns coefficients in A to the monomials in r(x) remains linear, and we write it as Let us take the example of a one-dimensional extended space as those used in Sect. 4. The domain is supposed to be periodic (circle) and the nodes are indexed by 0 ≤ n < N x .Recall that the node of index N x is associated with the extra {1}.For 0 ≤ n < N x , the stencil s n is defined as the set of 2L + 1 nodes of index n − L, n − L + 1, . ..n + L − 1, n + L, plus the extra node of index N x .The stencil s N x 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 N p = 1 + N x (2 + L) monomials.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 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 longdistance dependencies (such as teleconnections).

Homogeneity
Furthermore, a symmetry hypothesis could optionally be used by assuming translational invariance of the ODEs, called homogeneity 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 N a .
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 (1) 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 N a = 18 coefficients of A are 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 twodimensional turbulence satisfy both constraints and yet generate non-uniform 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, a 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, t k+1 − t k = N k c h, 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).We define intermediate state vectors in between [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 an RK scheme applied to x ≡ x k,l , with 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 (N RK = 1, 2, and 4, respectively).Provided the dynamics are autonomous, a general RK scheme reads as 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 N RK −1 , followed by the sum Eq. (12a) to get f A (x).
In the following, h will be absorbed into 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 to be Gaussian with 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 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 , supposed to be white in time and uncorrelated from the observation errors.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 of identifying a model or building 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,  which would require finding 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 with this joint pdf is by definition J (A, x 0:K ) = − ln(p(A, x 0:K |y 0:K )).Because Eq. ( 14) is Markovian and given the Gaussian form of both model and observational errors, the cost function reads as 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 and directly 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 notably similar to a traditional least-square function used in machine and deep learning regression.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.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 p(y 0 , A) = p(y 0 |A)p(A) is neglected in Eq. ( 18), then the maximum a posteriori should not depend on a global rescaling of Q 1:K .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 analytically derive the gradient of Eq. ( 17) with respect to both A and x 0:K .As for x 0:K , we have where assuming x 0 is independent of 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 adjoints ∇ x k H k T of the tangent linear operators of the observation operators 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.

Discussion of the theoretical points
In this section, we discuss 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 factually 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 becoming accessory.In both cases, but particularly in the latter one, 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 with any part of the surrogate model.For instance, assume that the reference model is partially identifiable, which means that part of the reference model could 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 N x ×N p out of which the evaluation of J (A, x 0:K ) is hazardous.A failure to evaluate the cost function also depends on the number of integration steps since the instabilities are likely to increase exponentially with N k c .This says 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 says that we should strike an empirical compromise between the composition numbers N k c and the easiness in evaluating J (A, x 0:K ).On 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 sensitivity 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 in variables is the one implemented for all numerical applications described in Sect.4, together with the normalisation.From these experiments, we learned that these tricks often turned critical in the first iterates of the optimisation when the estimate of A progressively migrated to the A-stability domain.After a few iterations, however, 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.Weinan, 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) and 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 the resolvent to exhibit more realistic long-range correlations.Even when using a reasonably small stencil, long-range 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 a 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 is 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 is 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, addition of extra parameters, adaptive batch normalisation, and 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 correlations 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.

Model 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 models.
2. The L96 model as defined by ODEs defined over a periodic domain of variables indexed by n = 0, . .., N x − 1 where N x = 40: where , and F = 8.This model is an idealised representation of a onedimensional latitude band of the Earth's atmosphere.Its Lyapunov time is 0.60.It is integrated using the RK4 scheme and with δt r = 0.05.
3. The KS model, as defined by the PDE: 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 stands for the integer division of m by 10.We use the original values for the parameters: c = 10 for the timescale 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 timescale 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 ODE 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 ODE 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 ODE model A obtained from the minimisation.With a view to comparing 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 truth (reference model) is δt r over the time window [t 0 , t K ].This parameter only matters for the reference model integration since only the training time steps t k+1 − t k and the output of the model y 0:K (which may include knowledge of the observation operator) are known to the observer.
The integration time step of the surrogate model within the 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 and 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 and that 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: perfectly identifiable models
In the first couple of experiments, we consider a densely observed2 reference model with noiseless observations.In this case, H k is the identity operator, i.e. each grid point value is observed, and r k ≡ 0 so that a uniform rescaling of the Q k , all chosen to be I x , is irrelevant, assuming p(y 0 , A) = p(y 0 |A)p(A) can be neglected, which is hypothesised here and is generally true for large K.Moreover, we use the same numerical scheme with 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, meaning that it belongs to the set of all 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 pre- cision 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 Lyapunov times.
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 (see Sect. 2.2.1 and 2.2.2).The stencil has a width of 5 (i.e. the local grid points with two points on its left and two 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 , advection terms, dissipation) are retrieved with a precision of a 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 Lyapunov times, and reaches a saturation for a lead time of 25 Lyapunov times.

Inferring the dynamics from dense and noiseless observations: non-identifiable models
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 K = 50.
As opposed to the reference model, in these nonidentifiable 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 the 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 is 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 be.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 ODE 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 either might 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 , given that Eq. ( 17) is now used instead of Eq. ( 18).For the sake of simplicity, we have chosen them to be of the scalar forms r k ≡ σ 2 y I y and Q k ≡ qI x .In these synthetic experiments, σ y is supposed to be known, while q is not.We only have a qualitative view of the potential mismatch between the reference and the surrogate model.Moreover, a Gaussian additive noise might not even be the best statistical model for such error.Nonetheless holding to the above Gaussian assumptions for model 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).The use of advanced methods of that sort to estimate model error will be considered in future works, while, in the following, we have chosen values of q that yield nearto-optimal skill scores (typically q ∈ [10 −4 , 1]σ 2 y ).Moreover, note that we have chosen the relatively small K = 50.While we expect increasing K to be beneficial, especially with noisy observations, it would force us to address issues related 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 aspect 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.This is confirmed by Fig. 8, where the precisions in identifying the model, measured by either the spectral norm Using the same setup, we have also reduced the number of observations.The observations of grid point values 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 when 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 the 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 shows the same but in the significantly noisy case where σ y = 1 and with 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 were encountered when scaling up the method.The large-scale variable x of the reference model is noiselessly and fully observed over a short training window (K = 50, which corresponds to about 0.35 Lyapunov time), i.e. all slow variable values are observed, 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 δt a = δt f = t/N c .
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 sub-scale perturbations that are not properly represented by the surrogate model.Reciprocally there are long-lasting stripes of low error not yet impacted by sub-scale 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 in 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 composition 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 highdimensional 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 models yielding a perfect reconstruction close to machine precision (L63 and L96 models).It has also been applied to densely or par-    tially observed, identifiable or non-identifiable models with or without noise in the observations (L96 and KS models).
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 multiscale 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 built on an ensemblebased 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 use of advanced weak constraint variational optimisation techniques.
In this paper, only autonomous dynamics have been considered.We could at least partially extend the method to non-autonomous 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 schemes 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.
All of these results, Eqs.(B10), (B14), (B20), and (B24), allow us to efficiently compute the gradients of the cost function J (A, x 0:K ) with respect to both A and x 0:K .Note, however, that they have been derived under the simplifying assumption that φ A is given by Eq. ( 2) with a traditional matrix multiplication between A and r(x), but not by the compact Eq. ( 8).When relying on homogeneity and/or locality, the calculation of the gradient with respect to A follows the principle described above but requires further adaptations, which can be derived using Eq.(A2), with the asset of strongly reducing the computational burden.

Appendix C: Adjoint differentiation with PyTorch and TensorFlow
As an alternative to the explicit computation of the gradients of Eq. ( 17) and the associated adjoint models, we have used PyTorch and TensorFlow as automatic differentiation tools.Only the cost function code needs to be implemented.We made a few tests of the experiments of Sect. 4 that showed that the fastest code is a C++ implementation using the explicit gradients, followed by an implementation using Python/Numpy/Numba using the explicit gradients, followed by a much slower implementation with TensorFlow (graph execution), followed by an implementation with TensorFlow (eager execution) and finally with PyTorch.Our experiments made use of a multi-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 used in the way we built our cost function with relatively few but significant parameters, and (ii) the time spent in the development of each approach scales with their speed.

Figure 1 .
Figure 1.Representation of the data assimilation system as a hidden Markov chain model and of the model resolvents F k A and f A .

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 δt r , δt a , and δt f are arbitrary in this schematic.

Figure 3 .
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 4 .
Figure 4. Density plot of the L96 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 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 deviations σ y .

Figure 8 .
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 9 .
Figure 9. L96 is the reference model, which is fully observed without noise: plot of the N a = 18 coefficients of the surrogate model as a function of the minimisation iteration number.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.

Figure 10 .
Figure10.L96 is the reference model, which is fully observed with observation error standard deviation σ y = 1: plot of the N a = 18 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.

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).Panel (d) shows a zoom of the difference between times 4 and 5.