{'fraction_non_alphanumeric': 0.03208456973293768, 'fraction_numerical': 0.01094213649851632, 'mean_word_length': 4.414574666475398, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
DAN-An optimal Data Assimilation framework based on machine learning Recurrent Networks

Data assimilation algorithms aim at forecasting the state of a dynamical system by combining a mathematical representation of the system with noisy observations thereof. We propose a fully data driven deep learning architecture generalizing recurrent Elman networks and data assimilation algorithms which provably reaches the same prediction goals as the latter. On numerical experiments based on the well-known Lorenz system and when suitably trained using snapshots of the system trajectory (i.e. batches of state trajectories) and observations, our architecture successfully reconstructs both the analysis and the propagation of probability density functions of the system state at a given time conditioned to past observations.


Introduction.
1.1. Context. In Data assimilation (DA, [1]), the time dependent state of a system is estimated using two models that are observation model, which relates the state to physical observations, and the dynamical model, that is used to propagate the state along the time dimension. If DA was originally developed in the area of Numerical Weather Prediction [22], it is now used in many areas such as Neutronics [19] or Oceanography [23], to name only a few.
Both observational and dynamical models are described using random variables that account for observation and model errors. Hence DA algorithms are grounded on a Bayesian approach in which prior information is combined with the above statistical models to obtain the probability density of the system state, conditioned to past observations. The estimation is done in two steps: the analysis of an incoming observation at a given time step t, and the propagation to the next time step t + 1.
Data assimilation algorithms use additional assumptions to obtain closed expressions for the probability densities so that they can be handled by computers. Historically in the linear Kalman filter (KF, [15]) approach, statistical models are supposed to be Gaussian and operators linear. Hence, the propagation and analysis steps consist in updating mean and covariance matrix of the involved densities. In the Ensemble Kalman Filter (EnKF, [11]) approach, these densities are represented by a (hopefully small) set of sampling vectors, and formula are available for the analysis and the propagation steps as well.
This above presentation of KF and EnKF, intentionally hides the statistical principles that lead to the mathematical expressions implemented in algorithms (the reader interested in more details is referred to recent references, e.g. [1]). It however puts forward a functional view of DA in which the analysis and propagation steps are seen as transformations, acting on probability densities. In practice, how to choose the density representation and how to implement the related analysis and propagation steps is left to the user. He will articulate his choice according to the precision he aims for the prediction and to many other practical aspects including the computer memory he can access to or the computational time he can dispose for the prediction. For instance, large atmospheric DA systems preclude the direct manipulation of covariance matrices because of memory issues and favor ensemble approaches or diffusion operators that model the covariance matrix effect on a vector [12].
The central theoretical question we address in this paper is the following: For a given memory requirement for the density representation, is it possible to learn the analysis and propagation transformations using data consisting of a batch of system state and noisy observation trajectories ?
We will give a positive answer to this question. The bottom line of this data driven approach is to let a learning process estimate what is the best internal representation for the densities that will substitute for mean and covariance in the KF or to an ensemble of state realizations in the case of the EnKF. This representation will be stored in some computer memory H in a (hopefully not too long) vector. We propose to use machine learning to estimate the analysis and propagation transformation acting on this memory. These transformations will take the form of recurrent networks, hence the name Data Assimilation Networks (DAN) we use throughout this paper. The memory H will be connected to the system state densities with a third transformation we call the "procoder" that will also be learned. We show that this general algorithm, that is a variant of the Elman Network, outperforms state of the art weakly constrained ensemble variational techniques in the case of twin experiments with the 40 variables Lorenz system. The reason behind this good performance is that the DAN is trained to be a better approximation (in a sense that will be made precise) of the ideal Bayesian solution than standard DA alternatives.

Related work.
Combining Data Assimilation and Machine Learning is not a new idea. It has notably been explored along three directions we outline now.
In a first approach, ingredients of the observed dynamical system are learned from noisy observations and a preexisting DA algorithm. Here, the DA algorithm acts on the system state. In [6,4], a surrogate propagation operator (depending on some parameters) and the system state are conjointly estimated from noisy observations. This is done in two phases. In the first phase, the propagation operator parameters are fixed and a preexisting DA algorithm estimates the states from noisy observations. In a second phase, these states and observations are used as data in the training of the propagation operator parameters and the process is iterated.
In a second approach, a DA algorithm is emulated using machine learning techniques. The authors use the term emulated to stress the fact that the procedure computes a surrogate of an a priori given DA algorithm. In [13], an Elman neural network is trained using the outputs of an EKF. It is claimed that the resulting surrogate DA algorithm is cheaper than the original Kalman filter but fails to be more accurate. As well [7] uses an artificial neural network to emulate a localized ensemble transform Kalman filter (LETKF).
Finally, DA algorithms can be seen as optimizers to train recurrent networks. In this approach, the DA algorithm acts on the network parameters. The sequential nature of the DA algorithms is advantageously used to optimize cost functions that are sum of terms that have a backward dependency for a given ordering. In [20], an extended Kalman filter (EKF) is used to train a recurrent network. In [18], an ensemble Kalman filter (EnKF) is adapted to update the weights of a recurrent neural network (RNN) and in [9] a theoretical parallel is drawn between the Kalman filter and the FORCE learning in recurrent neural networks. In [14], the EKF is mentioned besides the truncated back propagation through time as a method to train RNNs.
The first two approaches inherit the limitations of the DA algorithm they approximate. Our objective is to challenge those DA algorithms and to use optimal ML techniques for training RNNs to outperform existing DA. The last approach is definitely closer to our work in the sense that, as we shall see, they both rely on the truncated backpropagation through time technique (TBPTT, [14]).

Outline.
In section 2, we depict the statistical assumptions on the time dependent system we address, that we term an observed dynamical system (ODS), in this paper. Using time invariant transformations, we present, three common classes of data assimilation algorithms. A general formalism is proposed in section 3 for the above DA algorithms, that is similar, though not identical to Elman networks (that is a common architecture in machine learning). In section 4, we introduce a cost function on data assimilation networks, that allows us to optimize them provided data from an ODS are given. Finally in section 5, our data assimilation networks are implemented using neural networks. Their train and test performance are evaluated on the 40 variables Lorenz95 system and compared with state of the art data assimilation techniques.
2.1. Observed Dynamical Systems. Data assimilation assumes an observed dynamical system (ODS) is given. Such an autonomous system can be described by the following equations between random variables.
where M is the propagation operator that acts on the state random variable x t at some time index t and approximately returns the state random variable x t+1 at time t + 1. These state random variables take their values in some space X of dimension n (often R n ). They are initialized at time 0 with a variable x 0 that is usually assumed Gaussian N µ b 0 , Σ b 0 with background mean µ b 0 and background covariance Σ b 0 . The propagation equation (2.1a) accounts for the fact that the propagation operator is only known up to an additive model error η t that is usually assumed Gaussian N (0 n , Q) with zero mean and covariance matrix Q. As well, H is the observation operator that acts on the state random variable x t and approximately returns the observation random variable y t at time t. This observation random variable takes its values in some space Y of dimension d (often R d ). In (2.1b), the observation equation accounts for the fact that the observation operator is only known up to an additive observation error ε t that is usually assumed Gaussian N (0 d , R) with zero mean and covariance matrix R. The family of observation and propagation errors along time is usually assumed mutually independent.
The former ODS description (2.1) is a formal relation between random variables from which they can be simulated (sampled) on a computer. However most DA algorithms aim at quantifying the uncertainty over the system state when an observation sample (realization) becomes available. Such an analysis starts by rewriting, under suitable mathematical assumptions, the ODS in terms of probability density functions (pdfs) and conditional probability density functions (cpdfs). The cpdf equivalent to the observation equation (2.1b) is The cpdf p y t |xt is a function from X to the set P Y of pdfs over Y i.e. an element of the function set X → P Y . Specifically, given a realization x t ∈ X of the state random variable x t , the observation cpdf returns the pdf over observations p y t |xt (x t ) ∈ P Y which is Gaussian with mean H (x t ) and covariance matrix R. In turn, given an element y t ∈ Y, this pdf returns the scalar p y t |xt (x t ) (y t ) but this cumbersome notation is usually avoided with the notation p y t |xt (y t |x t ). Applying these conventions to the propagation equation (2.1a) as well yields the following ODS description that is equivalent to (2.1) where the cpdf input and output spaces have been made explicit as it will improve the readability of this paper.
Using the Markovian properties equivalent to the error independence, i.e state and observations only depend on the most recent state realization, the cpdfs (2.3) appear as factors of the joint pdf over the state and observation random trajectories (x 0:T , y 0:T ) = (x 0 , . . . , x T , y 0 , . . . , y T ) ∈ X T +1 × Y T +1 until some time index T : It determines all states cpdfs given past observations which are the objects of DA algorithms. In this paper we will rather focus on the filtering problem where only the most recent state is estimated and where p xt,y 0:t is considered.
These definitions now allow us to present in a functional way three well-known DA algorithms namely the Bayesian DA algorithm (BDA), the Kalman filter (KF) and the ensemble Kalman filter (EnKF).

Sequential
Bayesian data assimilation. From the ODS equations, the random states depend explicitly on the random observations and on model and observation errors. Estimating the state amounts to quantify its uncertainty as a function of the observations. In a Bayesian framework, this is done through the computation of the state cpdf given past observations where p xt,y 0:t and p y 0:t are marginals of the joint pdf in (2.4). Highlighting input and output spaces, a rigorous definition for the state cpdf is In addition, DA algorithms take a realization of the random observation trajectory as input and recursively evaluates these conditional pdfs on it. Specifically, at cycle t, the analysis step computes the posterior pdf p xt|y 0:t (y 0:t ) ∈ P X from the prior pdf p xt|y 0:t−1 (y 0:t−1 ) ∈ P X and the current observation y t with Bayes rule The propagation step computes the next cycle prior pdf p x t+1 |y 0:t (y 0:t ) from the posterior pdf p xt|y 0:t (y 0:t ) by marginalization Using (2.7) and (2.8) recursively starting from a given initial prior p x 0 ∈ P X provides a closed expression for the filtering densities p xt|y 0:t for t ∈ {0, . . . , T }.
However, owing to the ODS structure, p xt|y t does not depend on time and will be simply denoted by p y|x . The analysis step (2.7) can be written using a function of pdfs or transformation that is invariant over time as As well, the former propagation (2.8) is related to a specific application of a more general propagation transformation that is invariant over time. Denoting the time invariant quantity p x t+1 |xt by p x|x , the propagation transformation is the function that transforms a posterior pdf in P X into the next cycle prior pdf in P X given by In the above presentation of sequential Bayesian DA, the DA cycle is the composition of the two time invariant transformations a BDA and b BDA . However, this expression is never used in practice because it involves multidimensional integral and a density parameterization to store them in a computer memory. We now present popular approximations that circumvent these issues.

The Kalman filter.
In the Kalman filter (KF, [15]), the propagation operator M is supposed affine with M as linear part and the observation operator H also affine with H as linear part. In this case, the analysis and propagation transformations preserve Gaussian pdfs that are easily characterized by their mean and covariance matrix. The analysis and propagation transformations then boil down to algebraic expressions on these pairs as we shall see in this section.
Suppose that the internal representation of a Gaussian pdf is formalized by the following injective transformation where S X is the set of mean and covariance matrix pairs over X. The KF analysis transformation is the function that transforms such a prior pair in S X and an observation in Y into the posterior pair in S X given by , and the diagram in Figure 1a commutes.
As well, the KF propagation transformation is the function that transforms a posterior pair in S X into the next cycle prior in S X given by , and the diagram in Figure 1b commmutes.    Unfortunately, operators linearity is rarely met in practice and covariance matrices may not be easy to store and manipulate in the case of large scale problems. A popular dimension reduction approach is the ensemble Kalman filter that has proven effective in several large scale applications.

The ensemble Kalman filter.
In the Ensemble Kalman Filter (EnKF, [11]), statistics (µ, Σ) ∈ S X are themselves estimated from an ensemble matrix X ∈ X m = R n×m having m columns with the empirical estimators ∈ R m×m and I m ∈ R m×m is the identity matrix. Thus, the algebra over mean and covariance matrices pairs can boils down to algebra on ensembles. The advantage is that nonlinear operators can be evaluated columnwise on ensembles and ensembles with few columns produce low-rank approximations of covariance matrices. Hence ensembles are an internal representation for the pdfs that are transformed by the following function into a Gaussian pdf The EnKF analysis transformation is the function that transforms such a prior ensemble in X m and an observation in Y into the posterior ensemble in X m given by Up to operators non-linearities and sampling errors diagram in Figure 1c commutes. As well, the EnKF propagation transformation is the function that transforms a posterior ensemble in X m into the next cycle prior ensemble in X m given by where N ∈ X m is m samples of N (0 n , Q) leading to the commutativity of the diagram in Figure 1d in the linear case, up to sampling errors.
3. Data assimilation networks. In section 2, we have seen three DA algorithms which are the Bayesian DA, the Kalman filter and the ensemble Kalman filter. In this presentation, we abstracted the analysis and the propagation steps into the three transformations a, b and c. A common pattern for DA algorithms appears while inspecting the input and output spaces of these functions. More precisely, given some set H, a Data Assimilation Network (DAN) is a triplet of transformations such that where the internal set of pdf representation H is P X in the BDA (c is the identity), S X in the KF and X m in the EnKF. A representation of a DAN is given Figure 2a. The term "procoder" is a contraction of "probability coder" as the function c transforms an internal representation into an actual pdf over X.
In machine learning, this network is similar to an Elman network (EN, [10]) which is a basic structure of recurrent network. It is made of the following pair of functions where Y is some set of inputs, X is some set of outputs and H is some set of memory. This memory makes the Elman network able to process variable length sequences of inputs to produce sequences of outputs. Indeed, given a new input y t ∈ Y in the input sequence, the function f updates the memory from h t−1 to h t = f (h t−1 , y t ) and the function g decodes the memory into the new output x t = g (h t ) in the output sequence. In a way, the memory of an Elman network is expected to gather relevant information from the past inputs to perform satisfactory predictions. In machine learning, the training process will optimally induce how to manipulate the memory from data.
The similarity between DAN and EN can be made explicit with the composition b • a ∈ H×Y → H as an analogous for f and c as an analogous for g. However the EN does not perform DA operations in all its generality. First it can not make predictions without observations, that is estimating strict future states from past observations. This is because the function f performs both the propagation and the analysis at once. In a way, the EN only produces posterior outputs and no prior outputs while the DAN produces prior or posterior outputs by applying the procoder c before or after the propagater b. It can also produces strict future predictions without observations by applying the propagater b multiple times before applying the procoder c. Second, the DAN provides a probabilistic representation of the state i.e. an element in P X instead of an element in X. Also, note that the compositions of b and c make a generalized propagation operator as it propagates in time probabilistic representations of the state rather than punctual realizations.
While the differences between DANs and ENs justify the introduction of the former, their similarities allow to adapt the training of an EN to a DAN. Specifically, we can train (e.g. by gradient descent) broader families of functions (e.g. neural networks) from data to perform DA instead of designing DA algorithms by hand from speculative hypotheses. The advantage is that we do not have to specify what are H, a, b, c, the training process will find their optimal value, provided a satisfactory cost function is given as well has a good sampling procedure of the ODS variability. We shall prove in the rest of the paper that this methodology may provide better result than existing DA algorithms that implement additional simplifications to have tractable formula. For instance in ensemble methods, the Jacobian of the trajectory is approximated by finite differences.
4. Training a data assimilation network.

Cost function definition and properties.
To train the DAN using sample from the ODS, a cost function has to be designed. Our motivation will be to quantify the information loss committed each cycle by the DAN prior and posterior approximation of the ideal BDA prior and posterior.
Let h a −1 ∈ H be some initial memory. Given a trajectory of observation y 0:T ∈ Y T +1 a DAN outputs a trajectory of prior and posterior cpdfs over X. That is any DAN (a, b, c where the element h a −1 ∈ H is confounded with the function h a −1 ∈ Y 0 → H from the singleton set Y 0 at the initialization t = 0. This recursive application is represented in Figure 2b.
The Bayesian Data assimilation from p x 0 similarly yields an element p = p b 0:T , p a 0:T ∈ P which is in this case the usual filtering cpdfs ∀t ∈ 0 : T, We want the DAN to generate a trajectory of priors and posteriors candidate cpdfs q ∈ P close enough to the ideal Bayesian one p ∈ P in the sense of information theory. The proximity from candidate pdfs to ideal pdfs is quantified by their relative entropy [8]. In our case, it is natural to sum the relative entropies between corresponding pdfs in the trajectories p and q. However this would require the intractable expression of the Bayesian cpdfs. Theorem 4.1 shows that this objective is equivalent to the minimization of the following sampled time averaged cross entropy: ∀q = q b 0:T , q a 0:T ∈ P, where the data x 0:T,0:I ∈ X (T +1)×(I+1) , y 0:T,0:I ∈ Y (T +1)×(I+1) samples the ODS joint pdf (2.4) 1 .
Theorem 4.1. Given the ODS joint pdf p x 0:T ,y 0:T (2.4), let P X be the set of pdfs over X and G X be the set of Gaussians pdfs over X, let be the sets of prior and posterior cpdf and Gaussian cpdf trajectories respectively, let p = p b 0:T , p a 0:T ∈ P be the Bayesian prior and posterior cpdfs trajectory, let g = g b 0:T , g a 0:T ∈ G be the Gaussian prior and posterior cpdf trajectory such that ∀t ∈ 0 : T, Proof. The proof consists in revealing the relative entropies [8] in the cross entropies. Then the properties of L will derive from properties of the relative entropy.
Let q = q b 0:T , q a 0:T ∈ P, the expected time averaged cross entropy writes t − ln q a t dp x 0:T ,y 0:T , (4.8) where the integration variables have been discarded in the compact notation dp x 0:T ,y 0:T for p x 0:T ,y 0:T (x 0:T , y 0:T ) dx 0:T dy 0:T . Integral and finite sum inversion followed by marginalization yields t dp x 0:T ,y 0:T − T t=0 ln q a t dp x 0:T ,y 0:T , (4.9a) ln q a t dp xt,y 0:t , (4.9b) but the Bayesian cpdfs definition (4.5a) and (4.5b) yield the factorizations p xt,y 0:t−1 = p b t p y 0:t−1 and p xt,y 0:t = p a t p y 0:t thus ln q a t dp a t dp y 0:t , (4.10) with the compact notations dp b t dp y 0:t−1 for p b t (x t |y 0:t−1 ) dx t p y 0:t−1 (y 0:t−1 ) dy 0:t−1 and dp a t dp y 0:t for p a t (x t |y 0:t ) dx t p y 0:t (y 0:t ) dy 0:t . Introducing the following dependent relative entropies (written with the integration variables) ∀t ∈ 0 : T, t dp a t dp y 0:t , (4.12a) q a t ] dp y 0:t . (4.12b) Then p ∈ arg min L because the relative entropies are positive [8]. Moreover if q ∈ arg min L then ∀t ∈ 0 : T, are quadratic functions thus their expectation only involve the mean and covariance of p b 4.2. Truncated back propagation through time. However, it will be shown in the numerical experiment subsection 5.1 that the direct optimization of L (4.3a) is impractical. The gradient backpropagation through time in repeated compositions of (4.1) generate a large computational graph, that consumes a lot of memory. This limits the time averaged cross entropy time length T and batch size I which, in turn, leads to overfitting due to limited data. A workaround is to resort to gradient descent with truncated backpropagation through time (TBPTT, [14]). TBPTT is a sequential (online) method which may be seen as a DA method in itself; as a result it consumes little memory. More precisely, let θ be some parameter to be trained in some parameter set O. A DAN a, b, c thus becomes a parametric family of DANs a θ , b θ , c θ indexed by θ. Assume at some time step t ∈ 0 : T we have a θ-constant memorỹ h a t−1 ∈ Y t → H (the tilde notation is to avoid confusion with (4.1c)) and some iterate θ t−1 ∈ O. Then redefine (4.1) with ∀θ ∈ O, Provided the θ-dependent DAN is derivable over its parameters, the next iterate and the next θ-constant memory of the truncated through time gradient descent are then given by, where α is the learning rate. Note that this new θ-constant memory is evaluated during the evaluation of L t q b t,θ t−1 , q a t,θ t−1 . The method can be readily extended to more than one time step backward. However, this would be more computationally expensive and numerical experiment subsection 5.2 shows satisfactory results for only one time step backward compared to an existing technique.
In the numerical experiments, we will implement the parametric family of DANs with neural networks and the data (i.e a batch of trajectories of states and observations) will be generated by Gaussian perturbations of Lorenz95 [17] system trajectories.

Numerical experiments.
In these experiments, data assimilation networks are trained and tested in the framework of twin experiments. That is we generate training data i.e. a batch of state and observation trajectories from some given ODS. Then a parametric family of DANs is trained on these data. After this training, new observation trajectories are generated from a new unknown (for the trained DAN) state trajectory. This testing observations together with a null initial memory vector are then given as input of the trained DAN in a test phase and its outputs are compared with the unknown state.
The experiments setting is now detailed, it essentially that of test 1 (with T = 1) in [21] in order to compare with their IEnKF-Q which is a state of the art weak constrained, ensemble, variational DA algorithm. That is the IEnKF-Q takes model error into account in its derivation, it uses ensemble representation of the pdfs, the analysis is iterated in order to find the posterior pdf mode and the model error is taking into account by construction.
The experiments use the resolvant M ∈ R n → R n of the n = 40 dimensional Lorenz95 [17] system as the propagation operator, so X = R n . This resolvant is calculated with a Runge-Kutta 4 integration over a time step of δt = 0.05. The constitution of batches starts with an initial batch of states x 1:I,−burn which samples N (3 × 1 n , I n ) with burn = 10 3 and a batch size of I = 1024 for the training and I = 1 for the test (in order to see the scores fluctuations). Then this batch is propagated to t 0 to skip any transient regime ∀i ∈ 1 : I, After this initialization, independent Gaussian propagation errors η 1:I,0:T sampling N 0 n , q 2 I n with q = 0.1 are added each subsequent propagations to get a batch of state trajectories ∀i ∈ 1 : I, ∀t ∈ 0 : T, Then independent Gaussian observation errors ε 1:I,0:T sampling N 0 n , r 2 I n with r = 1 are added to the observation operator evaluations to get a training batch of observation trajectories ∀i ∈ 1 : I, ∀t ∈ 0 : T, where H is the identity operator H = id.
Artificial Neural networks are used to construct the parameterized family of DANs. The memory set is H = R m with m = 40 × 20 in order to compare with 20 members ensemble Kalman filters. The parametric family of analyzers and propagaters, are 20 layers fully connected neural networks with LeakyReLU [3] activation functions where θ gathers all weights. We use the ReZero trick [2] to improve trainability, that is each non-terminal layer verifies output = input + αLeakyReLU (Weight × input + bias) , (5.5) where α is a trainable parameter initialized at 0. The procoders c θ are implemented as as simple linear layers without activation function (H → R n+ n(n+1) 2 ) followed by a function that transforms a vector v ∈ R n+ n(n+1) This structure ensures the covariance matrix invertibility and good conditioning. Finally, trainings are performed with the Adam optimizer [16] and a learning rate of 10 −4 on a NVIDIA GTX 1080 graphic card.

Direct optimization.
In this experiment, we found that the cost function maximum time length is about T = 25 before the graphic card runs out of memory for a batch size of I = 1024. Indeed, the evaluation tree becomes massive with the compositions of a, b and c over this time window. This illustrates a limitation of the direct optimization. Scores are displayed in Figure 3 with the following organization.
We performed 1.3 × 10 4 optimization steps. The train time averaged cross entropy L (4.3b) is plotted as a function of these optimization steps in left panel, dark blue while the train prior and posterior time averaged rmses are plotted in right panel, dark orange and dark green resp. Those time averaged rmses are defined with where the prior, posterior mean µ b i,t , µ a i,t are the mean of q b t (x i,t |y i,0:t−1 ) , q a t (x i,t |y i,0:t ) in (4.1) respectively.
Each 2 × 10 3 optimization steps (marked by a vertical gray bar), the weights θ of the DAN are saved and fixed for a test over 10 3 time steps. This test scores are plotted on the left of the vertical bar, e.g at optimization step 6 × 10 3 , the test scores are plotted from iteration 5×10 3 to iteration 6×10 3 (the tests are spaced to stress their independence and to lighten the graph). The test data are of course different from the train data and the memory is initialized to zero each test to simulate incorrect memory initialization. The test instantaneous cross entropy L t (4.3b) is plotted as a function of the time step t in the left panel, light blue while the test prior and posterior instantaneous rmses rmse b t (5.8b) and rmse a t (5.8d) are plotted in right panel, light orange and light green respectively.
It appears that the time average of the L t , which approximates the value of L on the test data, becomes much greater than the train value of L towards the end of the training (iteration 11 × 10 3 ). This implies that the DAN overfits the data. As a consequence, the right panel shows poor results. Indeed, the rmses are greater than the observation error standard deviation which is 1 in our setting. This means our DAN is ineffective because it has poorer results than the trivial DA algorithm returning observations without any treatment. The reason is that the GPU memory limits the data set to T = 25 time steps trajectories and a batch size of I = 1024 which is too small to satisfactorily sample the lorenz system. Therefore there is no point in performing more optimization steps. We need to reduce the optimization memory (and computational) cost to reach longer time intervals.

Truncated back propagation through time.
Truncated back propagation through time is a sequential (online) method; as a result it consumes little memory and the optimization is not anymore limited in time by memory constraints. Metrics are displayed in Figure 4 with the following organization.
We did 6 × 10 5 optimization steps with TBPTT. The train instantaneous cross entropy L t (4.3b) is plotted as a function of these optimization steps in the left panel (dark blue) while the train prior and posterior instantaneous rmses rmse b t , rmse a t are plotted in the right panel (dark orange, dark green) respectively. Each 2 × 10 5 optimization steps (marked by a vertical gray bar), the weights θ of the DAN are saved and kept constant for a test over 10 5 time steps. Again, the test data are of course different from the train data and the memory is initialized to zero each test. The test instantaneous cross entropy L t (4.3b) is plotted as a function of the time step t in the left panel (light blue) while the test prior and posterior instantaneous rmses rmse b t (5.8b) and rmse a t (5.8d) are plotted in the right panel (light orange, light green) respectively. Note that we do not display the training values of L as this would require a DAN application over 6×10 5 time steps for each time step. However, the time average of the instantaneous test losses L t , are close to the train losses L t . This illustrates the TBPTT good behavior that avoids both overfitting and bias. In comparison with the test instantaneous loss L t in Figure 3, the results are much better. This translates into improved rmses in the right panel. The IEnKF-Q [21] posterior instantaneous rmses (with optimal inflation) are also displayed as a reference in the right panel (red). At the beginning of each test, the DAN metrics present a peak. Specifically, its rmses attain 4 (which is known to be the rmse between lorenz95 independent trajectories) because the memory is initialized at zero to simulate poor prior information on the state; on the other hand, the IEnKF-Q culminates at around 1 because the initial ensemble as 1 as std. After the peak, those metrics quickly decrease indicating that the DAN succeeds in retrieving the state. For a precise comparison of the trained DAN and the IEnKF-Q, Table 1 displays their time averaged posterior rmses for a test over 6 × 10 5 time steps, the DAN shows the best results which is significant (for instance, it corresponds to the improvement brought by the IEnKF-Q over non iterative methods in [21], test 1).  Table 1: Time averaged posterior (filtering) rmses 6. Conclusion. The common point of most Bayesian based DA algorithms is that they are written as a succession of pdf transformations occurring in the analysis and propagation steps. We have shown that these densities can be expressed as the unique solution to an optimization problem that is convex on common examples of density spaces. Next, once an approximation space for the pdf is given, and represented by some parameterization (in some memory H), the analysis and propagation steps of DA consist in applying time-invariant transformations a and b that update the pdfs using incoming observations. Among other possible choices, if we represent a and b as neural networks, the estimation problem takes the form of the minimization of a loss using recurrent networks (DAN). We believe these theoretical connection between Bayesian DA and recurrent networks is new, and works extremely well as demonstrated on a Lorenz system in which an up-to-date ensemble technique is outperformed. This rigorous connection between Bayesian DA and Machine Learning and its numerical relevance opens of course a number of avenues for future research.
A natural question concerns scalability of the estimation. It is well known that the accuracy of the estimation of Gaussian densities in the Ensemble Kalman filter is a function of the number of samples. In our numerical study we purposely chose the same memory storage for AD and for our technique. A natural question would be to study the accuracy of both techniques in experiments where an increasingly large memory is available.
In this work we decided to have a genuine data driven approach, in which neither the time evolution, nor the observational model, nor the Bayesian formula are provided to the learning and test process once the data is sampled. This result is crucial, and its good statistical properties are thoroughly verified in our experiments. However one may argue that in some situations, reliable model information is available. This may accelerates the learning process by reducing the amount of data needed to achieve a given performance. For instance it sounds reasonable that if the observation scenario is rich enough, in some sense to be determined, state trajectories are dispensable, and the estimation should be satisfactory. One could also consider the situation in which one would like to correct a DA procedure for errors with simple signatures, such as a bias. Our estimation process should benefit from this a priori information making the estimation much easier than a full reconstruction of both a dynamical and an observation model.
Having demonstrated that a DAN can be at least as efficient as a well-established EnKF, an important question could be to better understand why such a phenomenon occurs. We could for instance think that this results from the good approximation properties of neural networks based approximation space for a and b. We could also argue that Kalman formula is essentially grounded on linear approximations and that even iterating these formulae fails to provide an accurate estimation, whereas NN are naturally nonlinear. One could for instance explore numerically this question, by introducing a series of dynamical systems with a parameterized nonlinearity. One could also compare the estimated transformations a, b and c with the corresponding quantities available in the EnKF and explore their dissimilarities.
Another question that pertains to machine learning is whether DANs are amenable to transfer learning. This question relates to the possibility of reusing an already-trained DAN on a slightly modified observed dynamical system. Transfer learning usually relies on an estimation of all or parts of the NN parameters. In standard DA application transfer learning would be useful to handle systems that are slowly varying along physical time.
Finally, practitioners may also be interested in a number of issues that are related to computational time. Of course, DANs are naturally adapted to modern platforms, as they rely on machine learning software tools, the performance of which being strategic for most most hardware designers. One very good aspect of DAN is that the propagation and analysis are computed by the inference steps of the neural architecture, without expensive iterations for handling nonlinearities. Yet remains to be shown in practical cases that computationally affordable neural architectures can be found enough to provide satisfactory results.