Optimization Based Methods for Partially Observed Chaotic Systems

In this paper we consider filtering and smoothing of partially observed chaotic dynamical systems that are discretely observed, with an additive Gaussian noise in the observation. These models are found in a wide variety of real applications and include the Lorenz 96' model. In the context of a fixed observation interval $T$, observation time step $h$ and Gaussian observation variance $\sigma_Z^2$, we show under assumptions that the filter and smoother are well approximated by a Gaussian with high probability when $h$ and $\sigma^2_Z h$ are sufficiently small. Based on this result we show that the Maximum-a-posteriori (MAP) estimators are asymptotically optimal in mean square error as $\sigma^2_Z h$ tends to $0$. Given these results, we provide a batch algorithm for the smoother and filter, based on Newton's method, to obtain the MAP. In particular, we show that if the initial point is close enough to the MAP, then Newton's method converges to it at a fast rate. We also provide a method for computing such an initial point. These results contribute to the theoretical understanding of widely used 4D-Var data assimilation method. Our approach is illustrated numerically on the Lorenz 96' model with state vector up to 1 million dimensions, with code running in the order of minutes. To our knowledge the results in this paper are the first of their type for this class of models.


Introduction
Filtering and smoothing are amongst the most important problems for several applications, featuring contributions from mathematics, statistics, engineering and many more fields; see for instance Crisan and Rozovskii [2011] and the references therein. The basic notion of such models, is the idea of an unobserved stochastic process, that is observed indirectly by data. The most typical model is perhaps where the unobserved stochastic process is a Markov chain, either in discrete time, or a diffusion process. In this paper, we are mainly concerned with the scenario when the unobserved dynamics are deterministic and moreover chaotic.
The only randomness in the unobserved system is uncertainty in the initial condition and it is this quantity that we wish to infer, on the basis of discretely and sequentially observed data; we explain the difference between filtering and smoothing in this context below. This class of problems has slowly become more important in the literature, particularly in the area of data assimilation (Law, Stuart, and Zygalakis [2015]). The model itself has a substantial number of practical applications, including weather prediction, oceanography and oil reservoir simulation, see for instance Kalnay [2003].
In this paper we consider smoothing and filtering for partially observed deterministic dynamical systems of the general form where u : R + → R d is a dynamical system in R d for some d ∈ Z + , A is linear operator in R d (i.e. A is a d × d matrix), f ∈ R d is a constant vector, and B(u, u) is a bilinear form corresponding to the nonlinearity (i.e. B is a d × d × d array). We denote the solution of equation (1.1) with initial condition u(0) := v for t ≥ 0 by v(t). The derivatives of the solution v(t) at time t = 0 will be denoted by in particular, D 0 v = v , Dv := D 1 v = −Av − B(v, v) + f (the right hand side of (1.1)), In order to ensure the existence of a solution to the equation (1.1) for every t ≥ 0, we assume that there are constants R > 0 and δ > 0 such that (1.3) We call this the trapping ball assumption. Let B R := {v ∈ R d : v ≤ R} be the ball of radius R. Using the fact that d dt v(t), v(t) = 1 2 d dt v(t) 2 , one can show that the solution to (1.1) exists for t ≥ 0 for every v ∈ B R , and satisfies that v(t) ∈ B R for t ≥ 0.
Equation (1.1) was shown in Sanz-Alonso and  and Law, Stuart, and Zygalakis [2015] to be applicable to three chaotic dynamical systems, the Lorenz 63' model, the Lorenz 96' model, and the Navier-Stokes equation on the torus; such models have many applications. We note that instead of the trapping ball assumption, these papers have considered different assumptions on A and B(v, v). As we shall explain in Section 1.1, their assumptions imply (1.3), thus the trapping ball assumption is more general.
We assume that the system is observed at time points t j = jh for j = 0, 1, . . ., with observations Y j := Hu(t j ) + Z j where H : R d → R do is a linear operator, and (Z j ) j≥0 are i.i.d. centered random vectors taking values in R do describing the noise. We assume that these vectors have distribution η that is Gaussian with i.i.d. components of variance σ 2 Z . 1 The contributions of this article are as follows. In the context of a fixed observation interval T , we show under assumptions that the filter and smoother are well approximated by a Gaussian law when σ 2 Z h is sufficiently small. Our next result, using the ideas of the first one, shows that the Maximum-a-posteriori (MAP) estimators (of the filter and smoother) are asymptotically optimal in mean square error when σ 2 Z h tends to 0. The main practical implication of these mathematical results is that we can then provide a batch algorithm for the smoother and filter, based on Newton's method, to obtain the MAP. In particular, we prove that if the initial point is close enough to the MAP, then Newton's method converges to it at a fast rate. We also provide a method for computing such an initial point, and prove error bounds for it. Our approach is illustrated numerically on the Lorenz 96' model with state vector up to 1 million dimensions. We believe that the method of this paper has a wide range of potential applications in meteorology, but we only include one example due to space considerations.
We note that in this paper, we consider finite dimensional models. There is a substantial interest in the statistics literature in recent years in non-parametric inference for infinite dimensional PDE models, see Dashti and Stuart [2017] for an overview and references, and Giné and Nickl [2016] for a comprehensive monograph on the mathematical foundations of infinite dimensional statistical models. This approach can result in MCMC algorithms that are robust with respect to the refinement of the discretisation level, see e.g. Cotter, Roberts, Stuart, and White [2013], Cotter, Dashti, and Stuart [2012], Pillai, Stuart, and Thiéry [2014], Vollmer [2015], Cui, Law, and Marzouk [2016], Yao, Hu, and Li [2016]. There are also other randomisation and optimization based methods that have been recently proposed in the 1 We believe that our results in this paper hold for non-Gaussian noise distributions as well, but proving this would be technically complex.
A key property of these methods is that the prior is defined on the function space, and the discretizations automatically define corresponding prior distributions with desirable statistical properties in a principled manner. This is related to modern Tikhonov-Phillips regularisation methods widely used in applied mathematics, see Benning and Burger [2018] for a comprehensive overview. In the context of infinite dimensional models, MAP estimators are non-trivial to define in a mathematically precise way on the infinite dimensional function space, but several definitions of MAP estimators, various weak consistency results under the small noise limit, and posterior contraction rates have been shown in recent years, see e.g. Cotter, Dashti, Robinson, and Stuart [2009], Dashti, Law, Stuart, and Voss [2013], Vollmer [2013], Helin and Burger [2015], Kekkonen, Lassas, and Siltanen [2016], Monard, Nickl, and Paternain [2017], Nickl [2017], Dunlop and Stuart [2016]. Some other important work on similar models and/or associated filtering/smoothing algorithms include Hayden, Olson, and Titi [2011], Blomker, Law, Stuart, and Zygalakis [2013], Law, Sanz-Alonso, Shukla, and Stuart [2016]. These results are very interesting from a mathematical and statistical point of view, however the intuitive meaning of some of the necessary conditions, and their algorithmic implications are difficult to grasp.
In contrast with these works, our results in this paper concern the finite dimensional setting that is the most frequently used one in the data assimilation community. By working in finite dimensions, we are able to show consistency results and convergence rates for the MAP estimators under small observation noise / high observation frequency limits under rather weak assumptions (in particular, in Section 4.4 of Paulin, Jasra, Crisan, and Beskos [2018] our key assumption on the dynamics was verified in 100 trials when A, B and f were randomly chosen chosen, and only the first component of the system was observed, and they were always satisfied). Moreover, previous work in the literature has not said anything about the computational complexity of actually finding the MAP estimators, which is a non-trivial problem in non-linear setting due to the existence of local maxima for the loglikelihood. In our paper we propose appropriate initial estimators, and show that Newton's method started from them converges to the true MAP with high probability in the small noise/high observation frequency scenario when started from this initial estimator.
It is important to mention that the MAP estimator forms the basis of the 4D-Var method introduced in Le Dimet and Talagrand [1986], Talagrand and Courtier [1987] that is widely used in weather forecasting. A key methodological innovation of this method is that the gradients of the log-likelihood are computed via the adjoint equations, so that each gradient evaluation takes a similar amount of computation effort as a single run of the model. This has allowed the application of the method on large scale models with up to d = 10 9 dimensions.
See Dimet and Shutyaev [2005] for some theoretical results, and Navon [2009], Bannister [2016] for an overview of some recent advances. The present paper offers rigorous statistical foundations for this method for the class of non-linear systems defined by (1.1).
The structure of the paper is as follows. In Section 1.1, we state some preliminary results for systems of the type (1.1). Section 2 contains our main results: Gaussian approximations, asymptotic optimality of MAP estimators, and approximation of MAP estimators via Newton's method with precision guarantees. In Section 3 we apply our algorithm to the Lorenz 96' model. Section 4 contains some preliminary results, and Section 5 contains the proofs of our main results. Finally, the Appendix contains the proofs of our preliminary results based on concentration inequalities for empirical processes.

Preliminaries
Some notations and basic properties of systems of the form (1.1) are now detailed below.
The one parameter solution semigroup will be denoted by Ψ t , thus for a starting point v = (v 1 , . . . , v d ) ∈ R d , the solution of (1.1) will be denoted by Ψ t (v), or equivalently, v(t). Sanz-Alonso and  and Law, Stuart, and Zygalakis [2015] have assumed that the nonlinearity is energy conserving, i.e. B(v, v), v = 0 for every v ∈ R d . They also assume that the linear operator A is positive definite, i.e. there is a λ A > 0 such that Av, v ≥ λ A v, v for every v ∈ R d . As explained on page 50 of Law, Stuart, and Zygalakis [2015], (1.1) together with these assumptions above implies that for every v ∈ R d , From (1.4) one can show that B R is an absorbing set for any thus all paths enter into this set, and they cannot escape from it once they have reached it.
This in turn implies the existence of a global attractor (see e.g. Temam [1997], or Chapter 2 of Stuart and Humphries [1996]). Moreover, the trapping ball assumption (1.3) holds.
For t ≥ 0, let v(t) and w(t) denote the solutions of (1.1) started from some points v, w ∈ R d . Based on (1.1), we have that for any two points v, w ∈ B R , any t ≥ 0, and therefore by Grönwall's lemma, we have that for any t ≥ 0, is a one-to-one mapping, which has an inverse that we denote as The main quantities of interest of this paper are the smoothing and filtering distributions corresponding to the conditional distribution of u(t 0 ) and u(t k ), respectively, given the The densities of these distributions will be denoted by µ sm (v|Y 0:k ) and µ fi (v|Y 0:k ). To make our notation more concise, we define the observed part of the dynamics as for any t ∈ R and v ∈ B R . Using these notations, the densities of the smoothing and filtering distributions can be expressed as where det stands for determinant, and Z sm k , Z fi k are normalising constants independent of v. Since the determinant of the inverse of a matrix is the inverse of its determinant, we have the equivalent formulation (1.10) We assume a prior q on the initial condition that is absolutely continuous with respect to the Lebesgue measure, and zero outside the ball B R (where the value of R is determined by the trapping ball assumption (1.3)).
For k ≥ 1, we define the kth Jacobian of a function g : R d1 → R d2 at point v as a k + 1 dimensional array, denoted by J k g(v) or equivalently J k v g , with elements We define the norm of this kth Jacobian as Using ( (1.12) By induction, we can show that for any i ≥ 2, and any v ∈ R d , we have (1.13) From this, the following bounds follow (see Section A.1 of the Appendix for a proof).
In some of our arguments we are going to use the multivariate Taylor expansion for vector valued functions. Let g : R d1 → R d2 be k + 1 times differentiable for some k ∈ N.
Then using the one dimensional Taylor expansion of the functions g i (a + th) in t (where g i denotes the ith component of g), one can show that for any a, h ∈ R d1 , we have where h j := (h, . . . , h) denotes the j times repetition of h, and the error term R k+1 (a, h) (1.18) whose norm can be bounded using the fact that In order to be able to use such multivariate Taylor expansions in our setting, the existence and finiteness of J k Ψ t (v) can be shown rigorously in the following way. Firstly, for 0 < t < and using the inequality (1.15), we can show that and establish the existence of the partial derivatives by the chain rule.
After establishing the existence of the partial derivatives J k Ψ t (v), we are going to bound their norm in the following lemma (proven in Section A.1 of the Appendix).

Main results
In this section, we present our main results. We start by introducing our assumptions. In Section 2.1 we show that the smoother and the filter can be well approximated by Gaussian distributions when σ 2 Z h is sufficiently small. This is followed by Section 2.2 where based on the Gaussian approximation result we show that the Maximum-a-posteriori (MAP) estimators are asymptotically optimal in mean square error in the σ 2 Z h → 0 limit. We also show that Newton's method can be used for calculating the MAP estimators if the initial point x 0 can be chosen sufficiently close to the true starting position u. Finally, in Section 2.3 we propose estimators to use as initial point x 0 that satisfy this criteria when σ 2 Z h and h is sufficiently small.
We start with an assumption that will be used in these results.
Assumption 2.1. Let T > 0 be fixed, and suppose that T = kh, where k ∈ N. Suppose that u < R, and that there exist constants h max (u, T ) > 0 and c(u, T ) > 0 such that for every As we shall see in Proposition 2.1, this assumption follows from the following assumption on the derivatives (introduced in Paulin, Jasra, Crisan, and Beskos [2018]).
Assumption 2.2. Suppose that u < R, and there is an index j ∈ N such that the system of equations in v defined as has a unique solution v := u in B R , and where HD i u k refers to coordinate k of the vector HD i u ∈ R do , and ∇ denotes the gradient of the function in u.
The proof is given in Section A.1 of the Appendix. Assumption 2.2 was verified for the Lorenz 63' and 96' models in Paulin, Jasra, Crisan, and Beskos [2018] (for certain choices of the observation matrix H), thus Assumption 2.1 is also valid for these models.
We denote int(B R ) := {v ∈ R d : v < R} the interior of B R . In most of our results, we will make the following assumption about the prior q.
Assumption 2.3. The prior distribution q is assumed to be absolutely continuous with respect to the Lebesgue measure on R d , supported on B R . We assume that v → q(v) is strictly positive and continuous on B R , and that it is 3 times continuously differentiable at every interior point of B R . Let We assume that these are finite.
After some simple algebra, the smoothing distribution for an initial point u and the filtering distribution for the current position u(T ) can be expressed as where C sm k is a normalising constant independent of v (but depending on (Z j ) j≥0 ). In the following sections, we will present our main results for the smoother and the filter.
First, in Section 2.1 we are going to state Gaussian approximation results, then in Section 2.2 we state various results about the MAP estimators, and in Section 2.3 we propose an initial estimator for u based on the observations Y 0 , . . . , Y k , to be used as a starting point for Newton's method.

Gaussian approximation
We define the matrix A k ∈ R d×d and vector B k ∈ R d as (2.6) where JΦ ti and J 2 Φ ti denotes the first and second Jacobian of Φ ti , respectively, and If A k is positive definite, then we define the center of the Gaussian approximation of the smoother as and define the Gaussian approximation of the smoother as (2.9) If A k is not positive definite, then we define the Gaussian approximation of the smoother µ sm G (·|Y 0:k ) to be the d-dimensional standard normal distribution (an arbitrary choice), and u G := 0. If A k is positive definite, and u G ∈ B R , then we define the Gaussian approximation of the filter as Alternatively, if A k is not positive definite, or u G / ∈ B R , then we define the Gaussian approximation of the smoother µ fi G (·|Y 0:k ) to be the d-dimensional standard normal distribution. In order to compare the closeness between the target distributions and their Gaussian approximation, we are going to use two types of distance between distributions. The total variation distance of two distributions µ 1 , µ 2 on R d that are absolutely continuous with respect to the Lebesgue measure is defined as where µ 1 (x), and µ 2 (x) denote the densities of the distributions.
The Wasserstein distance (also called 1st Wasserstein distance) of two distributions µ 1 , µ 2 on R d (with respect to the Euclidean distance) is defined as where Γ(µ 1 , µ 2 ) is the set of all measures on R d × R d with marginals µ 1 and µ 2 .
The following two theorems bound the total variation and Wasserstein distances between the smoother, the filter, and their Gaussian approximations. In some of our bounds, the quantity T + h max (u, T ) appears. For brevity, we denote this as We are also going to use the constant C A defined as (2.14) Theorem 2.1 (Gaussian approximation of the smoother). Suppose that Assumptions 2.1 and 2.3 hold for the initial point u and the prior q. Then there are constants C TV (u, T ) + C  (2.20) Note that the Gaussian approximations µ sm G and µ fi G as defined above are not directly computable based on the observations Y 0:k , since they involve the true initial position u in both their mean and covariance matrix. However, we believe that with some additional straightforward calculations one could show that results similar to Theorems 2.1 and 2.2 also hold for the Laplace approximations of the smoothing and filtering distributions (i.e. when the mean and covariance of the normal approximation is replaced by the MAP and the inverse Hessian of the log-likelihood at the MAP, respectively), which are directly computable based on the observations Y 0:k .

MAP estimators
Let u sm be the mean of the smoothing distribution, and u fi be the mean of the filtering distribution, andû sm MAP be the maximum-a-posteriori of the smoothing distribution, i.e. u sm MAP := arg max v∈BR µ sm (v|Y 0:k ). (2.21) In case there are multiple maxima, we choose any of them. For the filter, we will use the push-forward MAP estimatorû fi := Ψ T (û sm MAP ). (2.22) Based on the Gaussian approximation results, we prove the following two theorems about these estimators. Suppose that Assumptions 2.1 and 2.3 hold for the initial point u and the prior q. Then there where the expectations are taken with respect to the random observations, and C sm (u, T ), C sm (u, T ), and C sm MAP (u, T ) are finite positive constants independent of σ Z and h.
Theorem 2.4 (Comparison of mean square error of MAP and posterior mean for filter).
Suppose that Assumptions 2.1 and 2.3 hold for the initial point u and the prior q. Then there is a constant S fi max (u, T ) > 0 independent of σ Z and h such that for 0 < h < h max (u, T ), where the expectations are taken with respect to the random observations, and C fi (u, T ), C fi (u, T ), and C fi MAP (u, T ) are finite positive constants independent of σ Z and h.
Remark 2.1. Theorems 2.3 and 2.4 in particular imply that when T is fixed, and σ Z √ h tends to 0, the ratio between the mean square errors of the posterior mean and MAP estimators conditioned on the initial position u tends to 1. Since the mean of the posterior distributions, u sm (or u fi for the filter), is the estimator U (Y 0:k ) that minimises E( U (Y 0:k ) − u 2 ) (or E( U (Y 0:k ) − u(T ) 2 ) for the filter), our results imply that the mean square error of the MAP estimators is close to optimal when σ Z √ h is sufficiently small.
Next we propose a method to compute the MAP estimators. Let g sm : (2.27) Then −g sm (v) is the log-likelihood of the smoother, except that it does not contain the normalising constant term.
The following theorem shows that Newton's method can be used to computeû sm MAP to arbitrary precision if it is initiated from a starting point x 0 that is sufficiently close to the initial position u. The proof is based on the concavity properties of the log-likelihood near u. Based on this, an approximation for the push-forward MAP estimatorû fi can be then computed by moving forward the approximation ofû sm MAP by time T according to the dynamics Ψ T (this will not increase the error by more than a factor of exp(GT ) according to (1.6)).
Theorem 2.5 (Convergence of Newton's method to the MAP). Suppose that Assumptions 2.1 and 2.3 hold for the initial point u and the prior q. Then for every 0 < ε ≤ 1, there exist finite constants S sm max (u, T, ε), N sm (u, T ) and D sm max (u, T ) ∈ (0, N sm (u, T )] (defined in (5.62) and (5.63)) such that the following holds. If σ Z √ h ≤ S sm max (u, T, ε), and the initial point x 0 ∈ B R satisfies that x 0 − u < D sm max (u, T ), then the iterates of Newton's method defined recursively as (2.29) Remark 2.2. The bound (2.29) means that number of digits of precision essentially doubles in each iteration. In other words, only a few iterations are needed to approximate the MAP estimator with high precision if x 0 is sufficiently close to u.

Initial estimator
First we are going to estimate the derivatives HD l u for l ∈ N based on observations Y 0:k .
For technical reasons, the estimators will depend on Y 0:k for some 0 ≤k ≤ k (which will be chosen depending on l). For any j ∈ N, , with the convention that 0 0 := 1.
We denote by I jmax+1 the identity matrix of dimension j max +1, and by e (l|jmax) a column vector in R jmax+1 whose every component is zero except the l + 1th one which is 1.
For any l ∈ N, j max ≥ l,k ≥ j max , we define the vector is an estimator of HD l u. The fact that the matrix M (jmax|k) (M (jmax|k) ) ′ is invertible follows from the fact that v (0|k) , . . . , v (jmax|k) are linearly independent (since the matrix with rows v (0|k) , . . . , v (k|k) is a so-called Vandermonde matrix whose determinant is nonzero). From (2.31), it follows that the norm of c (l|jmax|k) can be expressed as . (2.33) To lighten the notation, for j max ≥ l andk ≥ j max , we will denote . (2.34) The next proposition gives an error bound for this estimator, which we will use for choosing the valuesk and j max given l.
Proposition 2.2. Suppose that j max ≥ l andk ≥ 2j max + 3. Then for any 0 < ε ≤ 1, The following lemma shows that ask → ∞, the constant C (l|jmax|k) M tends to a limit. . (2.36) The proofs of the above two results are included in Section 5.4. Based on these results, we choosek ∈ {2j max +3, . . . , k} such that the function g(l, j max ,k) is minimised. We denote this choice ofk byk opt (l, j max ) (if g(l, j max ,k) takes the same value for severalk, then we choose the smallest of them). By taking the derivative of g(l, j max ,k) ink, it is easy to see that it has a single minimum among positive real numbers achieved at Based on this, we havê Finally, based on the definition ofk opt (l, j max ), we choose j opt max (l) as is a parameter to be tuned. We choose the smallest possible j max where the minimum is taken. Based on these notations, we define our estimator for HD l u asΦ (l) :=Φ (l|j opt max (l)) (Y 0:kopt(l,j opt max ) ). (2.40) The following theorem bounds the error of this estimator.
Theorem 2.6. Suppose that u ∈ B R , and T = kh. Then for any l ∈ N, there exist some positive constants h The following theorem proposes a way of estimating u from estimates for the derivatives HD i u 0≤i≤j . This will be used as our initial estimator for Newton's method.
Theorem 2.7. Suppose that for some j ∈ N there is a function F : (2.43) In particular, under Assumption 2.2 and for j as determined therein, function F defined as satisfies conditions (2.41) and (2.42).
Thus, the initial estimator can simply be chosen as with the above two theorems implying that the estimate gets close to u for decreasing σ Z √ h.
Solving polynomial sum of squares minimisation problems of this type is a well-studied problem in optimisation theory (see Lasserre [2010] for a theoretical overview), and several toolboxes are available (see Prajna, Papachristodoulou, and Parrilo [2002], Tütüncü, Toh, and Todd [2003]). Besides (2.44), other problem-specific choices of F satisfying conditions (2.41) and (2.42) can also be used, as we explain in Section 3.2 for the Lorenz 96' model.

Optimisation based smoothing and filtering
The following algorithm provides an estimator of u given Y 0:k . We assume that either there is a problem-specific F satisfying conditions (2.41) and (2.42) for some j ∈ N, or we suppose that Assumption 2.2 is satisfied for the true initial point u, and use F as defined in (2.44).
Step 1: We compute the estimatorsΦ (0) , . . . ,Φ (j) based on (2.40), and set the initial point Step 2: We compute the iterates x i for i ≥ 1 based on (2.28) recursively until x i − x i−1 becomes smaller than ∆ min , and returnû = x n for n := min i∈Z+ x i − x i−1 < ∆ min .
The following algorithm returns an online estimator of (u(t i )) i≥0 given Y 0:i at time t i .
Step 1: For i < k, return the estimate u(t i ) = 0.
Step 2: This algorithm can be modified to run Step 2 only at every K step for some K ∈ Z + (i.e. for i = k + lK for l ∈ N), and propagate forward the estimate of the previous time we with high probability for both algorithms (in the second algorithm, for i ≥ k).

Application to the Lorenz 96' model
The Lorenz 96' model is a d dimensional chaotic dynamical system that was introduced in Lorenz [1996]. In its original form, it is written as where the indices are understood modulo d, and f is the so-called forcing constant. In this paper we are going to fix this as f = 8 (this is a commonly used value that is experimentally known to cause chaotic behaviour, see Majda and Harlim [2012] and Majda, Harlim, and Gershgorin [2010]). As shown on page 16 of Sanz-Alonso and , this system can be written in the form (1.1), and the bilinear form B(u, u) satisfies the energy conserving We consider 2 observation scenarios for this model. In the first scenario, we assume that d is divisible by 6, and choose H such that coordinates 1, 2, 3, 7, 8, 9, . . ., d − 5, d − 4, d − 3 are observed directly, i.e. each observed batch of 3 is followed by a non-observed batch of 3. In this case, the computational speed is fast, and we are able to obtain simulation results for high dimensions. We consider first a small dimensional case (d = 12) to show the dependence of the MSE of the MAP estimator on the parameter k (the amount of observations), when the parameters σ Z and h are fixed. After this, we consider a high dimensional case (d = 1000002), and look at the dependence of the MSE of the MAP estimator on σ Z √ h.
In the second scenario, we choose H such that we observe the first 3 coordinates directly.
We present some simulation results for d = 60 dimensions for this scenario.
In Paulin et al. [2018], we have shown that in the second scenario, the system satisfies Assumption 2.2 for Lebesgue-almost every initial point u ∈ R d . A simple modification of that argument shows that Assumption 2.2 holds for Lebesgue-almost every initial point u ∈ R d in the first observation scenario too.
In each case, we have set the initial point as u = d+1 2d , d+2 2d , . . . , 1 (we have tried different randomly chosen initial points and obtained similar results). Figures 1, 2, and 3 show the simulation results when applying Algorithm 1 (optimisation based smoother) to each of these cases. Note that Algorithm 2 (optimisation based filter) applied to this setting yields very similar results.
In Figure 1, we can see that the MAP estimators RMSE (root mean square error) does not seem to decrease significantly after a certain amount of observations. This is consistent with the non-concentration of the smoother due to the existence of leaf sets, described in Paulin et al. [2018]. Moreover, by increasing k above 100, we have observed that the Newton's method often failed to improve significantly over the initial estimator, and the RMSE of the estimator became of order 10 −1 , significantly worse than for smaller values of k. We believe that this is due to the fact that as we increase k, while keeping σ Z and h fixed, the normal approximation of the smoother breaks down, and the smoother becomes more and more multimodal. Due to this, we are unable to find the true MAP when starting from the initial estimator, and settle down at another mode. To conclude, for optimal performance, it is important to tune the parameter k of the algorithm.
In Figure 2, we present results for a d = 1000002 dimensional Lorenz 96' system, with half of the coordinates observed. The observation time is T = 10 −5 . The circles correspond to data points with h = 10 −6 (so k = 10), while the triangles correspond to data points with h = 2 · 10 −7 (so k = 50). The plots show that the method works as expected for this high dimensional system, and that the RMSE of the estimator is proportional to σ Z √ h.
Finally, in Figure 3, we present results for a d = 60 dimensional system with the first 3 coordinates observed. The observation time is T = 10 −3 . The circles correspond to data points with h = 5 · 10 −5 , while the triangles correspond to data points with h = 2.5 · 10 −5 .
We can see that Algorithm 1 is able to handle a system which has only a small fraction of its coordinates observed. Note that the calculations are done with numbers having 360 decimal digits of precision. Such high precision is necessary because the interaction between the 3 observed coordinates of the system and some of the non-observed coordinates is weak.
The smaller observation noise (σ Z ≤ 10 −6 ). Finally, in the third simulation ( Figure 3), since only a very small fraction of the d = 60 coordinates is observed, the observation noise has to be very small in order for us to be able to recover the unobserved coordinates (σ Z ≤ 10 −120 ).
In all of the above cases, the error of the MAP estimator is several orders of magnitude less than the error of the initial estimator. When comparing these results with the simulation results of  using the 3DVAR and Extended Kalman Filter methods for the Lorenz 96' model, it seems that our method improves upon them, since it allows for larger dimensions, and smaller fraction of coordinates observed.
In the following sections, we describe the theoretical and technical details of these sim-

Some bounds for Lorenz 96' model
In this section, we will bound some of the constants in Section 1.1 for the Lorenz 96' model.
For B, given any u, v ∈ R d such that u , v ≤ 1, by the arithmetic-mean-root mean square inequality, and the inequality (ab) 2 ≤ a 2 +b 2 2 , we have shows that this bound is sharp, and B = 2. For simplicity, we have chosen the prior q as the uniform distribution on B R . Based on these, and the definitions (1.16) and (1.6), we have 3.2 Choice of the function F in the initial estimator In this section, we will construct a computationally simple function F satisfying the conditions of Theorem 2.7 for the two observation scenarios.
First, we look at the second scenario, when only the first 3 coordinates are observed.
We are going to show that for j = d−3

3
, it is possible to construct a function F : u, and F Hu, . . . , HD j v = u, thus satisfies the conditions of Theorem 2.7. Notice that In general, for m ≥ 1, by differentiating (3.3) m times, we obtain that thus for any m ≥ 1, Thus for any m ∈ N, we have a recursion for the mth derivative of u i based on the first m+1 derivatives of u i−1 and the first m derivatives of u i−2 and u i−3 . Based on this recursion, and the knowledge of the first j derivatives of u 1 , u 2 and u 3 , we can compute the first j − 1 derivatives of u 4 , then the first j − 2 derivatives of u 5 , etc. and finally the zeroth derivative of u 3+j (i.e. u 3+j itself).
In the other direction, therefore for m = 0, we have By differentiating (3.7) m times, we obtain that Thus for any m ∈ N, we have a recursion allowing us to compute the first m derivatives of u i based on the first m + 1 derivatives of u i+2 and the first m derivatives of u i+1 and u i+3 (with indices considered modulo d). This means that given the first j derivatives of u 1 , u 2 and u 3 , we can compute the first j − 1 derivatives of u d and u d−1 , then the first j − 2 derivatives of u d−2 and u d−3 , etc. and finally the zeroth derivatives of u d+2−2j , u d+1−2j .
Based on the choice j := d−3

3
, these recursions together define a function F for this case. From the recursion formulas, and the boundedness of u it follows that F is Lipschitz in a neighbourhood of Hu, . . . , HD j u as long as none of the coordinates of u is 0 (thus for Lebesgue-almost every u ∈ B R ).
We note that in the above argument, F might be not defined if some of the components of u are 0 (and the proof of Assumption 2.2 in Paulin et al. [2018] also requires that none of the components are 0). Moreover, due to the above formulas, some numerical instability might arise when some of the components of u are very small in absolute value. In Section A.2 of the Appendix, we state a simple modification of the initial estimator of Theorem 2.7 based on the above F that is applicable even when some of the components of u are zero.

Numerical solution of chaotic ODEs based on Taylor expansion
Let v ∈ B R , and i max ∈ N. The following lemma provides some simple bounds that allow us to approximate the quantities v(t) = Ψ t (v) for sufficiently small values of t by the sum of the first i max terms in their Taylor expansion. These bounds will be used to simulate the system (1.1), and to implement Newton's method as described in (2.28).
where D 0 = v, and for i ≥ 1, we have the recursion Proof. The bounds on the error in the Taylor expansion follow from (1.14). The recursion equation is just (1.13).
The following proposition shows a simple way of estimating v(t) for larger values t (the bound (3.11) is not useful for t ≥ 1 C der ). For v ∈ R d , we let be the projection of v on B R .
Then the error is bounded as Proof. Using (3.11) and the fact that the projection P BR decreases distances we know that and using the triangle inequality, we have so the claim follows.

Simulation details
The algorithms were implemented in Julia, and ran on a computer with a 2.5Ghz Intel Core i5 CPU. In all cases, the convergence of Newton's method up to the required precision (chosen to be much smaller than the RMSE) occurred in typically 3-8 steps.
In the case of Figure 1 (d = 12, half of the coordinates observed) the observation time T is much larger than 1 C der , so we have used the method of Proposition 3.1 to simulate from the system. The gradient and Hessian of the function g sm were approximated numerically based on finite difference formulas (requiring O(d 2 ) simulations from the ODE). We have noticed that in this case, the Hessian has elements with significantly large absolute value even far away from the diagonal. and Section 3.2.1 of Paulin, Jasra, Beskos, and Crisan [2017] (see also Le Dimet, Navon, and Daescu [2002]). This means that the Hessians were approximated using products of Jacobian matrices that were stored in sparse format due to the local dependency of the equations (3.1). This efficient storage has allowed us to run Algorithm 1 in approximately 20-40 minutes in the simulations. We made 4 parallel runs to estimate the MSEs. The

Preliminary results
The proof of our main theorems are based on several preliminary results. Let These quantities are related to the log-likelihoods of the smoothing distribution and its Gaussian approximation as Proposition 4.2 (A bound on the difference between l sm (v) and l sm G (v)). Suppose that Assumption 2.1 holds for u, then for any Proposition 4.3 (A bound on the difference between ∇l sm (v) and ∇l sm G (v)). Suppose that Assumption 2.1 holds for u, then for any 0 < ε ≤ 1, The following lemma is useful for controlling the total variation distance of two distributions that are only known up to normalising constants.
Lemma 4.1. Let f and g be two probability distributions which have densities on R d with respect to the Lebesgue measure. Then their total variation distance satisfies that for any Proof. We are going to use the following characterisation of the total variation distance, Based on this, the result is trivial for c = 1. If c > 1, then we have and the c < 1 case is similar.
The following lemma is useful for controlling the Wasserstein distance of two distributions.
Lemma 4.2. Let f and g be two probability distributions which have densities on R d with respect to the Lebesgue measure, also denoted by f (x) and g(x). Then their Wasserstein distance (defined as in (2.12)) satisfies that for any y ∈ R d , . Suppose first that γ = 0 and 1 − γ = 0. Let ν f,g denote the distribution of a random vector (X, X) on R d × R d for a random variable X with distribution with density m(x)/γ (the two components are equal). Let ν (2) f,g be a distribution on R d × R d with density ν (2) f,g (x 1 , x 2 ) :=f (x1) 1−γ ·ĝ (x2) 1−γ (the two components are independent).
We define the optimal coupling of f and g as a probability distribution ν f,g on R d × R d as a mixture of ν (1) f,g and ν (2) f,g , that is, for any Borel-measurable E ∈ R d × R d , It is easy to check that this distribution has marginals f and g, so by (2.12), and the fact and the result follows from the fact thatf (x) +ĝ(x) = |f (x) − g(x)|. Finally, if γ = 0 then we can set ν f,g as ν ( 2) f,g , and the same argument works, while if γ = 1, then both sides of the claim are zero (since f (x) = g(x) Lebesgue-almost surely).
The following two lemmas show concentration bounds for the norm of B k and the smallest eigenvalue of A k . The proofs are included in Section A.1 of the Appendix (they are based on matrix concentration inequalities).
Lemma 4.3. Suppose that Assumption 2.1 holds for u, then the random vector B k defined in (2.7) satisfies that for any t ≥ 0, (4.6) Thus for any 0 < ε ≤ 1, we have (4.7) Lemma 4.4. Suppose that Assumption 2.1 holds for u, then the random matrix A k defined in (2.6) satisfies that for any t ≥ 0, thus for any 0 < ε ≤ 1, we have The following proposition bounds the derivatives of the log-determinant of the Jacobian. det M + t 1 X 1 + . . . + t k X k .
They have defined the norm of kth derivative of the determinant as From Theorem 4 of Bhatia and Jain [2009], it follows that for any k ≥ 1, the norm of the kth derivative can be bounded as (4.10) Based on the chain rule, the norm first derivative of the log-determinant can be bounded as

Gaussian approximation
In the following two subsections, we are going to prove our Gaussian approximation results for the smoother, and the filter, respectively.

Gaussian approximation for the smoother
In this section, we are going to describe the proof of Theorem 2.1. First, we will show that the result can be obtained by bounding 3 separate terms, then bounds these in 3 lemmas, and finally combine them.
By choosing the constants C TV (u, T ) and C TV (u, T ) sufficiently large, we can assume that C TV (u, T, ε) satisfies the following bounds, Based on the assumption that σ Z √ h ≤ 1 2 C TV (u, T, ε) −1 , (5.1), and Lemma 4.4, we have where C A is defined as in (2.14). The event λ min (A k ) > c(u,T ) 2h implies in particular that A k is positive definite. From Lemma 4.3, we know that From Proposition 4.1, we know that Finally, from Proposition 4.2, it follows that In the rest of the proof, we are going to assume that all four of the events in the equations (5.6),(5.7),(5.8), and (5.9) hold. From the above bounds, we know that this happens with probability at least 1 − ε.
Let Z be a d dimensional standard normal random vector, then from the definition of µ sm G , it follows that when conditioned on B k and A k , has distribution µ sm G (·|Y 0:k ). This fact will be used in the proof several times. Since the normalising constant of the smoother, C sm k of (2.4) is not known, it is not easy to bound the total variation distance of the two distributions directly. Lemma 4.1 allows us to deal with this problem by rescaling the smoothing distribution suitably. We define the rescaled smoothing distribution (which is not a probability distribution in general) as which is of similar form as the Gaussian approximation then based on the assumption that σ Z √ h ≤ C TV (u, T, ε) −1 , (5.2) and (5.3), it follows that (5.14) Let B ρ := {v ∈ R d : v − u ≤ ρ(h, σ Z )}, and denote by B c ρ its complement in R d . Then by Lemma 4.1, we have By Lemma 4.2, we can bound the Wasserstein distance of µ sm (·|Y 0:k ) and µ sm G (·|Y 0:k ) as In the following six lemmas, we bound the three terms in inequalities (5.16) and (5.17).
Lemma 5.2. Using the notations and assumptions of this section, we have Proof. Note that if Z is a d dimensional standard normal random vector, then by Theorem 4.1.1 of Tropp [2015], we have for any t ≥ 0. (5.20) Since the random variable W defined in (5.10) is distributed as µ sm G (·|Y 0:k ) when conditioned on A k , B k , we have

Based on the assumption that σ
and the result follows by (5.20).
Lemma 5.3. Using the notations and assumptions of this section, we havẽ Proof. Let q max := sup v∈BR q(v). By our assumption that the event in (5.8) holds, we have where in the last step we have used the fact that det(A k ) 1/2 ≤ A k d/2 . Based on the assumption that σ Z √ h ≤ C TV (u, T, ε) −1 , and (5.5), we have therefore by (5.20), we havẽ and the claim of the lemma follows by (5.20).
From inequality (5.16) and Lemmas 5.1, 5.2 and 5.3, it follows that under the assumptions of this section, we can set C (1) TV (u, T ) and C (2) TV (u, T ) sufficiently large such that we have (5.21) Now we bound the three terms needed for the Wasserstein distance.
Lemma 5.4. Using the notations and assumptions of this section, we have for some finite positive constants C * 1 (u, T ), C * 2 (u, T ).
Proof. Note that for any v ∈ B ρ , we have From (5.21), and the fact thatμ sm (·|Y 0:k ) is a rescaled version ofμ sm (·|Y 0:k ), it follows that where Z is a d dimensional standard normal random vector. The claimed result now follows Lemma 5.5. Using the notations and assumptions of this section, we have Proof. Using (5.22), Lemma 5.3, and the fact that σ Z √ h ≤ 1 2 C TV (u, T, ε) −1 , we have , and the result follows from v∈B c ρ v − u µ sm (v|Y 0:k )dv ≤ 2R · µ sm (B c ρ |Y 0:k ).
Lemma 5.6. Using the notations and assumptions of this section, we have Proof. Let W be defined as in (5.10), and Z be d-dimensional standard normal. Using the fact that for a non-negative valued random variable X, we have , and the claim of the lemma follows (we have used Lemma 5.2 in the last step).
From inequality (5.17) and Lemmas 5.1, 5.2 and 5.3, it follows that under the assumptions of this section, for some appropriate choice of constants C (1) W (u, T ) and C (2) Proof of Theorem 2.1. The claim of the theorem follows from inequalities (5.21) and (5.23), and the fact that the assumption that all four of the events in the equations (5.6),(5.7),(5.8), and (5.9) hold happens with probability at least 1 − ε.

Gaussian approximation for the filter
In this section, we are going to describe the proof of Theorem 2.2. We start by some notation. We define the restriction of µ sm G (·|Y 0:k ) to B R , denoted by µ sm G|BR (·|Y 0:k ) as This is a probability distribution which is supported on B R . We denote its push-forward map by Ψ T as η fi G (·|Y 0:k ), i.e. if a random vector X is distributed as µ sm G|BR (·|Y 0:k ), then η fi G (·|Y 0:k ) denotes the distribution of Ψ T (X). The proof uses a coupling argument stated in the next two lemmas that allows us to deduce the results based on the Gaussian approximation of the smoother (Theorem 2.1).
Lemma 5.7 (Coupling argument for total variation distance bound). The total variation distance of the filtering distribution and its Gaussian approximation can be bounded as follows, Proof. First, notice that by Proposition 3(f) of Roberts and Rosenthal [2004], we have By Proposition 3(g) of Roberts and Rosenthal [2004], there is a coupling (X 1 , X 2 ) of random vectors such that X 1 ∼ µ sm (·|Y 0:k ), X 2 ∼ µ sm G|BR (·|Y 0:k ), and P(X 1 = X 2 |Y 0:k ) = d TV (µ sm (·|Y 0:k ), µ sm G|BR (·|Y 0:k )). Given this coupling, we look at the coupling of the transformed random variables (Ψ T (X 1 ), Ψ T (X 2 )). This obviously satisfies that P(Ψ T (X 1 ) = The statement of the lemma now follows by the triangle inequality.
Moreover, by the definitions, for a given Y 0:k , it is easy to see we can couple random variablesX 2 andX 2 ∼ µ sm G|BR (·|Y 0:k ) such that they are the same with probability at least 1 − µ sm G (B c R ). Since the maximum distance between two points in B R is at most 2R, it follows that By the triangle inequality, we obtain that d W (µ sm (·|Y 0:k ), µ sm G|BR (·|Y 0:k ))) ≤ d W (µ sm (·|Y 0:k ), µ sm G (·|Y 0:k )) + 2Rµ sm G (B c R ), and by (1.6), it follows that The claim of the lemma now follows by the triangle inequality.
As we can see, the above results still require us to bound the total variation and Wasserstein distances between the distributions η fi G (·|Y 0:k ) and µ fi G (·|Y 0:k ). Let (5.26) then the density of µ fi G (·|Y 0:k ) can be written as Since the normalising constant is not known for the case of η fi G (·|Y 0:k ), we define a rescaled versionη fi G (·|Y 0:k ) with densitỹ (5.28) The following lemma bounds the difference between the logarithms ofη fi G (v|Y 0:k ) and µ fi G (v|Y 0:k ).
Proof of Lemma 5.9. By (5.27) and (5.28), we have The absolute value of the first difference can be bounded by Propositon 4.4 as For any two vectors x, y ∈ R d , we have so the second difference can be bounded as Using (1.6), we have so by (1.6), it follows that We obtain the claim of the lemma by combining the stated bounds.
Now we are ready to prove our Gaussian approximation result for the filter.
Proof of Theorem 2.2. We suppose that D (1) TV (u, T ) and D TV (u, T ), thus D TV (u, T, ε) ≥ C TV (u, T, ε). We also assume that D Based on these assumptions on D TV (u, T, ε), and the assumption that σ Z √ h ≤ 1 2 D TV (u, T, ε) −1 , it follows that the probability that all the four events in the equations (5.6),(5.7),(5.8), and (5.9) hold is at least 1 − ε. We are going to assume that this is the case for the rest of the proof. We define (5.31) Based on (5.29), (5.30), and the assumption that σ Z √ h ≤ 1 2 D TV (u, T, ε) −1 , it follows that (5.32) The surface of a ball of radius R − u centered at u is contained in B R , and it will be transformed by Ψ T to a closed continuous manifold whose points are at least (R − u ) exp(−GT ) away from Ψ T (u) (based on (1.6)). This implies that the ball of radius By Lemma 4.1, we have (5.33) By Lemma 5.9, (5.32), and the fact that | exp(x) − 1| ≤ 2|x| for x ∈ [−1, 1], it follows that Therefore the first term of (5.33) can be bounded as This in turn can be bounded as in Lemma 5.1. The termsη fi G (B c ρ ′ |Y 0:k ) and µ fi G (B c ρ ′ |Y 0:k ) can be bounded in a similar way as in Lemmas 5.2 and 5.3. Therefore by Lemma 5.7 we obtain that under the assumptions of this section, there are some finite constants D For the Wasserstein distance bound, the proof is based on Lemma 5.8. Note that by the proof of Theorem 2.1, under the assumptions on this section, we have Therefore we only need to bound the last two terms of (5.25). The fact that µ sm G (B c R |Y 0:k )) = o(σ 2 Z h) can be shown similarly to the proof of Lemma 5.3. Finally, the last term can be bounded by applying Lemma 4.2 for y := Ψ T (u G ). This implies that These terms can be bounded in a similar way as in Lemmas 5.4, 5.5 and 5.6, and the claim of the theorem follows.

Comparison of mean square error of MAP and posterior mean
In the following two subsections, we are going to prove our results concerning the mean square error of the MAP estimator for the smoother, and the filter, respectively.

Comparison of MAP and posterior mean for the smoother
In this section, we are going to prove Theorem 2.3. First, we introduce some notation. Let The proof is based on several lemmas which are described as follows. All of them implicitly assume that the assumptions of Theorem 2.3 hold.
Lemma 5.10 (A bound on u sm − u G ). There are some finite constants D 5 (u, T ) and Proof. This is a direct consequence of the Wasserstein distance bound of Theorem 2.1, since Lemma 5.11 (A bound on û sm MAP − u ). For any 0 < ε ≤ 1, we have Proof. From Proposition 4.1, and (4.3), it follows that for any 0 < ε ≤ 1, Sinceû sm MAP is the maximizer of log µ sm (v|Y 0:k ) on B R , our claim follows.

Proof. By choosing S
(1) MAP and S (2) MAP sufficiently large, we can assume that (5.37) From Lemma 5.11, we know that Using (5.37), it follows that if the above event happens, then û sm MAP < R, and thus Using the fact that ∇l sm , and Proposition 4.3, it follows that (5.40) Moreover, by Lemma 4.4, we know that for any 0 < ε ≤ 1, we have By combining the four equations (5.38), (5.39), (5.40) and (5.41), it follows that with probability at least 1 − ε, we have and the claim of the lemma follows by rearrangement.
If this event happens, then in particular, we have By the definition of B k in (2.7), it follows that conditioned on u, B k has d-dimensional multivariate normal distribution with covariance matrix Σ B k := σ 2 Z k i=0 JΦ ti (u) ′ JΦ ti (u). This means that if Z is a d dimensional standard normal random vector, then (Σ B k ) 1/2 · Z has the same distribution as B k (conditioned on u). By Assumption 2.1, we have It is not difficult to show that for any d ≥ 1, (i) −1) ) = (1 + 2λ) −d/2 · e λd , and the claim follows by applying ) · e −λt for t = 3 4 d and λ = 1). Therefore, we have and thus by (5.43) and (5.44), it follows that By choosing D 10 (u, T ) sufficiently small, we have that for σ Z √ h ≤ D 10 (u, T ), and the result follows.
Lemma 5.14 (A bound on the difference of E u sm − u 2 u and E û sm MAP − u 2 u ). There are some finite constants D 11 (u, T ) and Proof. We define the event E k as (5.46) Under this event, we have, in particular A k ≻ 0 and u G < R. Let E c k denote the complement of E k . Then the difference in the variances can be bounded as for some finite constant C E (u, T ) independent of h and σ Z .
The term E 1 E k u G − u 2 u can be bounded as then by Lemma 5.10, it follows that for σ Z √ h < 1 2 (C (1) (5.50) By writing and using the fact that 1 E k u sm − u G 2 < 4R 2 , one can show that for σ Z √ h < S(u, T ), for some constants S(u, T ) > 0, C(u, T ) < ∞, that are independent of σ Z and h.
Finally, for bounding the term E 1 E k û sm MAP − u G 2 u , we define By Lemma 5.12, it follows that for The result now follows by (5.47) and the bounds (5.48), (5.49), (5.51) and (5.53).
Proof of Theorem 2.3. The lower bound on E( u sm −u 2 |u) σ 2 Z h follows by Lemma 5.13. Let E k be defined as in (5.46), then we have follows by (5.48), (5.49) and (5.51). Finally, the bound on E û sm MAP − u 2 |u − E u sm − u 2 |u follows directly from Lemma 5.14.

Comparison of push-forward MAP and posterior mean for the filter
The main idea of proof is similar to the proof of Theorem 2.3. We are going to use the following Lemmas (variants of Lemmas 5.10-5.14).
Lemma 5.15 (A bound on u fi − Ψ T (u G ) ). There are some finite constants D ′ 5 (u, T ) and D ′ 6 (u, T ) such that for any 0 < ε ≤ 1, for σ Z √ h ≤ 1 2 · D TV (u, T, ε) −1 , we have P c(u, T ) 2h Proof. This is a direct consequence of the Wasserstein distance bound of Theorem 2.2.
Lemma 5.17 (A lower bound on E u fi − u(T ) 2 u ). There are positive constants D ′ 9 (u, T ) and Proof. The proof is similar to the proof of Lemma 5.13. By applying Lemma 5.15 for ε = 0.1, If this event happens, then by (1.6), we have The rest of the argument is the same as in the proof of Lemma 5.13, so it is omitted.
Lemma 5.18 (A bound on the difference of E u fi − u 2 u and E û fi − u 2 u ). There are some finite constants D ′ 11 (u, T ) and Proof. We define the event E k as in (5.46). Under this event, A k ≻ 0 and u G < R, so µ fi G (·|Y 0:k ) is defined according to (2.10). The proof of the claim of the lemma follows the same lines as the proof of Lemma 5.14. In particular, we obtain from (5.49) and (1.6) that Based on Lemma 5.15, we obtain that for σ Z √ h < S ′ (u, T ), (5.59) for some constants S ′ (u, T ) > 0, C ′ (u, T ) < ∞, that are independent of σ Z and h. We omit the details.
Proof of Theorem 2.4. The lower bound on follows by Lemma 5.17. Similarly to the case of the smoother, we have which can be further bounded by (5.48), (5.58) and (5.59) to yield the upper bound on directly from Lemma 5.18.

Convergence of Newton's method to the MAP
In this section, we will prove Theorem 2.5. The following proposition shows a classical bound on the convergence of Newton's method (this is a reformulation of Theorem 5.3 of Bubeck [2015] to our setting). For v ∈ R d , r > 0, we denote the ball of radius r centered at v by B(v, r) Proposition 5.1. Suppose that Ω ⊂ R d is an open set, and g : Ω → R is a 3 times continuously differentiable function satisfying that 1. g has a local minimum at a point x * ∈ Ω, 2. there exists a radius r * > 0 and constants C H > 0, L H < ∞ such that B(x * , r * ) ⊂ Ω, ∇ 2 g(x) C H · I d for every x ∈ B(x * , r * ), and ∇ 2 g(x) is L H -Lipschitz on B(x * , r * ).
Suppose that the starting point x 0 ∈ Ω satisfies that x 0 − x * < min r * , 2 CH LH . Then the iterates of Newton's method defined recursively for every i ∈ N as always stay in B(x * , r * ) (thus they are well defined), and satisfy that Proof. We will show that x i ∈ B(x * , r * ) by induction in i. This is true for i = 0. Assuming that x i ∈ B(x * , r * ), we can write the gradient ∇g(x i ) as By the L H -Lipschitz property of ∇ 2 g(x) on B(x * , r * ), we have By combining this with the fact that (∇ 2 g( The following proposition gives a lower bound on the Hessian near u. The proof is included in Section A.1 of the Appendix. Then for any 0 < ε ≤ 1, σ Z > 0, 0 < h ≤ h max (u, T ), we have The following proposition bounds the Lipschitz coefficient of the Hessian. The proof is included in Section A.1 of the Appendix.
Proposition 5.3. Suppose that Assumptions 2.1 and 2.3 hold for the initial point u and the prior q. Then for any 0 < ε ≤ 1, σ Z > 0, 0 < h ≤ h max (u, T ), we have , Then by the assumption that σ Z √ h ≤ S sm max (u, T, ε), using Propositions 5.2, 5.3 and inequality (5.35) we know that with probability at least 1 − ε, all three of the following events hold at the same time,

Initial estimator
In this section, we will prove our results about the initial estimator that we have proposed in Section 2.3.
Proof of Proposition 2.2. By Taylor series expansion ofΦ (l|jmax) (Y 0:k ) with remainder term of order j max + 1, we obtain that where by (1.14) the remainder term R l,jmax+1 can be bounded using the Cauchy-Schwarz inequality as Due to the particular choice of the coefficients of c (l|jmax|k) , we can see that all the terms up to order j max disappear, and we are left with the remainder term that can be bounded as using the assumption thatk ≥ 2j max + 3. The concentration bound now follows directly from this bias bound, (5.20) and the fact that the estimatorΦ (l|jmax) (Y 0:k ) has Gaussian distribution with covariance matrix c (l|jmax|k) · σ Z · I do .
Proof of Lemma 2.1. Let P denote the space of finite degree polynomials with real coefficients on [0, 1]. For a, b ∈ P, we let a, b P = 1 x=0 a(x)b(x)dx. Then the elements of the matrix K jmax can be written as If K jmax would not be invertible, then its rows would be linearly dependent, that is, there would exist a non-zero vector α ∈ R jmax+1 such that However, this is not possible, since by the fundamental theorem of algebra, can have at most j max roots, so it cannot be zero Lebesgue almost everywhere in [0, 1].
Therefore K jmax is invertible. The result now follows from the continuity of the matrix inverse and the fact that for any 1 ≤ i, j ≤ j max + 1, Proof of Theorem 2.6. Let Based on Lemma 2.1, this is finite. With the choice j max = l, by (2.37), we obtain that and thus by the definition (2.38), we have By substituting this into (2.35), and applying some algebra, we obtain that g(l, l,k opt (l, l)) = C 0 H C l+1 and the claim of the theorem now follows by substituting this into Proposition 2.2.

Conclusion
In this paper, we have proven consistency and asymptotic efficiency results for MAP estimators for smoothing and filtering a class of partially observed non-linear dynamical systems.
We have also shown that the smoothing and filtering distributions are approximately Gaussian in the low observation noise / high observation frequency regime when the length of the assimilation window is fixed. These results contribute to the statistical understanding of the widely used 4D-Var data assimilation method (Le Talagrand [1986], Talagrand andCourtier [1987]). The precise size of the observation noise σ Z and assimilation step h under which the Gaussian approximation approximately holds, and the MAP estimator is close to the posterior mean is strongly dependent on the model parameters and the size of the assimilation window. However, we have found in simulations on Figure 2 that even for relatively large values of σ Z and h, for large dimensions, and not very short assimilation windows, these approximations seem to be working reasonably well. Besides theoretical importance, the Gaussian approximation of the smoother can be also used to construct the prior (background) distributions for the subsequent intervals in a flow-dependent way, as we have shown in Paulin et al. [2017] for the non-linear shallow-water equations, even for realistic values of σ Z and h. These flow-dependent prior distributions can considerably improve filtering accuracy. Going beyond the approximately Gaussian case (for example when σ Z , h, and T are large, or the system is highly non-linear) in a computationally efficient way is a challenging problem for future research (see Bocquet, Pires, and Wu [2010] for some examples where this situation arises).
Indeed, this was shown for k = 1 above, and for k ≥ 2, from (A.2), it follows that Proof of Lemma 4.3. Let (JΦ ti (u)) ·,j denote the jth column of the Jacobian, and letZ j i := Z j i /σ Z (which is a standard normal random variable). Then we can write B k as This is of the same form as in equation (4.1.2) of Theorem 4.1.1 of Tropp [2015]. Since JΦ ti (u) ≤ M 1 (T ), one can see that we also have (JΦ ti (u)) ·,j ≤ M 1 (T ), and thus the variance statistics v(Z) of Theorem 4.1.1 can be bounded as v(Z) ≤ σ 2 Z M 1 (T ) 2 (k + 1)d o . The result now follows from equation (4.1.6) of Tropp [2015] (with d 1 = 1 and d 2 = d).
Proof of Lemma 4.4. First, note that from Assumption 2.1 and looking at the Taylor expansion near u it follows that the first term in the definition of A k satisfies that λ min k i=0 JΦ ti (u) ′ JΦ ti (u) ≥ c(u, T ) h .
We can rewrite the second term in the definition (2.6) as where HΦ j ti (u) denotes the Hessian of the function Φ j ti at point u. This is of the same form as in equation (4.1.2) of Theorem 4.1.1 of Tropp [2015]. Since J 2 Φ ti (u) ≤ M 2 (T ), one can see that we also have HΦ j ti (u) ≤ M 2 (T ), and thus the variance statistics v(Z) of Theorem 4.1.1 can be bounded as v(Z) ≤ σ 2 Z M 2 (T ) 2 (k + 1)d o , and thus for any t ≥ 0, we have Proof of Proposition 2.1. The proof is quite similar to the proof of Proposition 4.1 of Paulin, Jasra, Crisan, and Beskos [2018]. With a slight modification of that argument, we will show that for any δ ∈ [0,h), for some constant c ′ (u,h) > 0 independent of δ, that is monotone increasing inh for 0 <h ≤h max . This result allows us to decouple the summation in (2.1) into sets of size j + 1 as follows. Let h 0 := min T j+1 ,h max , and seth := ⌊h 0 /h⌋ · h. Then for h ≤ h 0 , we have ⌊h 0 /h⌋ > h0 2h andh > h 0 /2, so using (A.6), we obtain that Thus Assumption 2.1 holds with h min (u, T ) := h 0 /2 and c(u, T ) := c ′ (u, h 0 /2) · h 0 /2.
To complete the proof, we will now show (A.6). Using inequality (1.14), we can see that the Taylor expansion difference between the approximation (A.7) and the derivative explicitly aŝ Proof. For d ≥ 3, this follows from Theorem 1 of Dumer [2007]. For d = 1, the result follows from the fact that 1 δ + 1 intervals suffice, and 1 δ + 1 ≤ 2 δ . For d = 2, we know that the circles of radius δ contains a square of edge length √ 2δ, and in order to cover a square of edge length 2 containing the unit ball, it suffices to use ⌈2/( √ 2δ)⌉ 2 ≤ √ 2 δ + 1 2 < 6 δ 2 squares of edge length √ 2δ, thus the result follows.
By the definition of W ′ l , this implies that E(W l ) ≤ 11(l + 1) √ ld + 1 · L √ k + 1. Moreover, By its definition, it is easy to see that for every 0 ≤ i ≤ k, ϕ i is L-Lipschitz with respect to the d 2 distance for L := M 2 (T ) + M 3 (T )R, and that ϕ i (v, s 1 , s 2 ) ≤ M for M := M 2 (T ) for every (v, s 1 , s 2 ) ∈ T 2 . The claim of the proposition now follows from Lemma A.2.
Proof of Proposition 5.3. By (4.1) and (4.3), we have Based on the assumption on q, and (1.19), the deterministic terms can be bounded as For the random terms, we first define ϕ i : T 3 → R (T 3 was defined in Lemma A.2) for 0 ≤ i ≤ k as Based on this, we let W 3 := sup (v,s1,s2,s3)∈T 3 k i=0 ϕ i (v, s 1 , s 2 , s 3 ), Z i , then one can see that the random terms can be bounded as By its definition, it is easy to see that for every 0 ≤ i ≤ k, ϕ i is L-Lipschitz with respect to the d 3 distance for L := M 3 (T )+ M 4 (T )R, and that ϕ i (v, s 1 , s 2 , s 3 ) ≤ M for M := M 3 (T ) for every (v, s 1 , s 2 , s 3 ) ∈ T 3 . The claim of the proposition now follows from Lemma A.2.
A.2 Initial estimator when some of the components are zero In Section 3.2, we have proposed a function F that allows us to express the un-observed coordinates of u from the observed coordinates and their derivatives (in the two observation scenarios described in Section 3). By substituting appropriate estimators of the derivatives, we obtained an initial estimator based on Theorem 2.7. Unfortunately, this function F was not defined when some of coordinates of u are 0. In this section we propose a modified version of this estimator that overcomes this difficulty.
We start by a lemma allowing us to run the ODE (1.1) backwards in time (for a while).
Lemma A.3. Suppose that v ∈ B R , and the trapping ball assumption (1.3) holds. Then for any 0 ≤ t < C −1 der , the series is convergent, well defined, and satisfies that Ψ t (Ψ −t (v)) = v and that for any i max ∈ N, Proof. The result follows from the bounds (1.14), and the definition of equation (1.1).
Based on this lemma, given the observations Y 0:k , we propose the following initial estimator. First, select some intermediate indices 0 = i 1 < i 2 < . . . < i m < k satisfying that i m · h < C −1 der . For each index i r , we compute the derivative estimatesΦ (l) of D l (u(t ir )), and then use the function F described in Section 3.2 to obtain initial estimators u(t ir ) of u(t ir ) for 0 ≤ r ≤ m. After this, we project these estimators to B R (see (3.13)), and run them backwards by t ir time units via the approximation (A.27), and project them back to B R , that is, for some sufficiently large i max ∈ N, we let The final initial estimatorû is then chosen as the one among (û r ) 0≤r≤m that has the largest the a-posteriori probability µ sm (û r |Y 0:k ) (see (2.4)). Based on (2.7), and some algebra, one can show that this estimator will satisfy the conditions required for the convergence of Newton's method (Theorem 2.5) if σ Z √ h and h are sufficiently small, and i max is sufficiently large, as long as at least one of the vectors (u(t ir )) 0≤r≤m has no zero coefficients. Moreover, by a continuity argument, it is possible to show that Assumption 2.2 holds as long as there is a t ∈ [0, T ] such that none of the coefficients of u(t) are 0.