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 σZ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _Z^2$$\end{document}, we show under assumptions that the filter and smoother are well approximated by a Gaussian with high probability when h and σZ2h\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^2_Z h$$\end{document} are sufficiently small. Based on this result we show that the maximum a posteriori (MAP) estimators are asymptotically optimal in mean square error as σZ2h\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^2_Z h$$\end{document} 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 among the most important problems for several applications, featuring contributions from mathematics, statistics, engineering and many more fields; see, for instance, [13] 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 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 [27]. The model itself has a substantial number of practical applications, including weather prediction, oceanography and oil reservoir simulation, see, for instance, [24].
In this paper we consider smoothing and filtering for partially observed deterministic dynamical systems of the general form du dt = −Au − B(u, u) + f , (1.1) where u : R + → R d is a dynamical system in R d for some d ∈ Z + , A is linear In order to ensure the existence of a solution to Eq. (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 Stuart [42] and Law et al. [27] 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 Sect. 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 where H : R d → R d o is a linear operator and (Z j ) j≥0 are i.i.d. centred random vectors taking values in R d o 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 nonparametric inference for infinite-dimensional PDE models, see [15] for an overview and references, and Giné and [21] 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, for example, [11,12,14,39,50,52]. There are also other randomization and optimization based methods that have been recently proposed in the literature, see, for example, [2,51].
A key property of these methods is that the prior is defined on the function space, and the discretisations 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 [4] 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 infinitedimensional 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, for example, [10,16,19,23,25,34,36,49]. Some other important works on similar models and/or associated filtering/smoothing algorithms include [6,22,28]. 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 to 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 et al. [38] our key assumption on the dynamics was verified in 100 trials when A, B and f were randomly 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 nonlinear setting due to the existence of local maxima for the log-likelihood. 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 [29] and Talagrand and Courtier [44] 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 [17] for some theoretical results, and Navon [35], Bannister [1] for an overview of some recent advances. The present paper offers rigorous statistical foundations for this method for the class of nonlinear systems defined by (1.1).
The structure of the paper is as follows. In Sect. 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 Sect. 3 we apply our algorithm to the Lorenz 96' model. Section 4 contains some preliminary results, and Sect. 5 contains the proofs of our main results. Finally, "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 Stuart [42] and Law et al. [27] 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 et al. [27], (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, for example, [45], or Chapter 2 of Stuart and Humphries [43]). 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,

)−w(t))−(B(v(t), v(t)−w(t)) − B(w(t)−v(t), w(t))),
and therefore by Grönwall's lemma, we have that for any t ≥ 0, for a constant G: 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 observations Y 0:k := {Y 0 , . . . , Y k }. 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 d 1 → R d 2 at point v as a k + 1-dimensional array, denoted by J k g(v) or equivalently J k v g , with elements Using ( By induction, we can show that for any i ≥ 2, and any v ∈ R d , we have From this, the following bounds follow (see Section A.1 of "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 d 1 → R d 2 be k + 1 times differentiable for some k ∈ N. Then using the one-dimensional Taylor expansion of the functions g i (a + t h) in t (where g i denotes the ith component of g), one can show that for any a, h ∈ R d 1 , we have (1.17) where h j :=(h, . . . , h) denotes the j times repetition of h and the error term whose norm can be bounded using the fact that t h) . (1.19) 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.
and using inequality (1.15), we can show that )) for t 1 + · · · + t m = t 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 "Appendix").

Lemma 1.2 For any k
Then for any k ≥ 1, T ≥ 0, we have

Main Results
In this section, we present our main results. We start by introducing our assumptions. In Sect. 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 Sect. 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 Sect. 2.3 we propose estimators to use as initial point x 0 that satisfy this criteria when σ 2 Z h and h are 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 v ∈ B R , for every h ≤ h max (u, T ) (or equivalently, every As we shall see in Proposition 2.1, this assumption follows from the following assumption on the derivatives (introduced in Paulin et al. [38]).

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 H D i u k refers to coordinate k of the vector H D i u ∈ R d o and ∇ denotes the gradient of the function in u.

Proposition 2.1 Assumption 2.2 implies Assumption 2.1.
The proof is given in Section A.1 of "Appendix". Assumption 2.2 was verified for the Lorenz 63' and 96' models in Paulin et al. [38] (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 Sect. 2.1 we are going to state Gaussian approximation results, then in Sect. 2.2 we state various results about the MAP estimators, and in Sect. 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 where JΦ t i and J 2 Φ t i denote the first and second Jacobian of Φ t i , respectively, and If A k is positive definite, then we define the centre of the Gaussian approximation of the smoother as and define the Gaussian approximation of the smoother as 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 first 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 (2.14) We are also going to use the constant C A defined as W (u, T ), and C (2) W (u, T ) independent of σ Z , h and ε such that for any 0 < ε ≤ 1, σ Z > 0, and W (u, T ), and D (2) W (u, T ) independent of σ Z , h and ε such that for any 0 < ε ≤ 1, σ Z > 0, and 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, 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∈B R μ sm (v|Y 0:k ). (2.23) In case there are multiple maxima, we choose any of them. For the filter, we will use the push-forward MAP estimatorû Based on the Gaussian approximation results, we prove the following two theorems about these estimators.
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.
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 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.29) 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)].
satisfy that P for every i ∈ N, x i is well defined and

Remark 2.2
The bound (2.31) means that the 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 H D 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. (2.32) For j max ∈ N, we define M ( j max |k) ∈ R ( j max +1)×(k+1) as a matrix with rows v (0|k) , . . ., v ( j max |k) . We denote by I j max +1 the identity matrix of dimension j max + 1, and by e (l| j max ) a column vector in R j max +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 H D l u. The fact that the matrix M ( j max |k) (M ( j max |k) ) is invertible follows from the fact that v (0|k) , . . . , v ( j max |k) are linearly independent (since the matrix with rows v (0|k) , . . . , v (k|k) is the so-called Vandermonde matrix whose determinant is nonzero). From (2.33), it follows that the norm of c (l| j max |k) can be expressed as . (2.35) To lighten the notation, for j max ≥ l andk ≥ j max , we will denote . (2.36) 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 ≥ 2 j max + 3. Then for any 0 < ε ≤ 1, (2.37) The following lemma shows that ask → ∞, the constant C (l| j max |k) M tends to a limit.
The proofs of the above two results are included in Sect. 5.4. Based on these results, we choosek ∈ {2 j 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 several k, 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 (2.39) Based on this, we havê } 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 H D l u aŝ The following theorem bounds the error of this estimator.

Theorem 2.7
Suppose that u ∈ B R , and T = kh. Then for any l ∈ N, there exist some positive constants h (l) max , s (l) max (T ) and S (l) max (T ) such that for any choice of the parameter J The following theorem proposes a way of estimating u from estimates for the derivatives H D i u 0≤i≤ j . This will be used as our initial estimator for Newton's method.
In particular, under Assumption 2.2 and for j as determined therein, function F defined as satisfies conditions (2.43) and (2.44).
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 optimization theory (see [26] for a theoretical overview), and several toolboxes are available (see [40,47]). Besides (2.46), other problem-specific choices of F satisfying conditions (2.43) and (2.44) can also be used, as we explain in Sect. 3.2 for the Lorenz 96' model.

Optimization 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.43) and (2.44) 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.46).
Step 1: We compute the estimatorsΦ (0) , . . . ,Φ ( j) based on (2.42), and set the initial point as Step 2: We compute the iterates x i for i ≥ 1 based on (2.30) recursively until x i − x i−1 becomes smaller than Δ min , and returnû = x n for n := min i∈Z The following algorithm returns an online estimator of (u(t i )) i≥0 given Y 0:i at time t i .

Algorithm 2 Optimization based filtering
Step 1: This algorithm can be modified to run Step 2 only at every K step for some K ∈ Z + (i.e. for i = k + l K for l ∈ N) and propagate forward the estimate of the previous time we ran Step 2 at the intermediate time points. This increases the execution speed at the cost of the loss of some precision (depending on the choice of K ).
Based on our results in the previous sections, we can see that if the assumptions of the results hold and Δ min is chosen sufficiently small then the estimation errors are of O(σ Z √ h) 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 [31]. 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 [32,33].) As shown on page 16 of Sanz-Alonso and Stuart [42], this system can be written in the form (1.1), and the bilinear form B(u, u) satisfies the energy-conserving property (i.e. 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 nonobserved 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 = 1,000,002) 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. [38], 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 (optimization based smoother) to each of these cases. Note that Algorithm 2 (optimization based filter) applied to this setting yields very similar results.
In Fig. 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. [38]. 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 Fig. 2, we present results for a d = 1,000,002-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 Fig. 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 requirements on the observations noise σ Z and observation time step h for the applicability of our method depend heavily on the parameters of the model (such as the dimension d) and on the observation matrix H. In the first simulation ( Fig. 1), relatively large noise (σ Z = 10 −3 ) and large time step (h = 10 −2 ) were possible. For the second simulation ( Fig. 1), due to the high dimensionality of the model (d = 1,000,002), and the sparse approximation used in the solver, we had to choose smaller time step (h ≤ 10 −6 ) and smaller observation noise (σ Z ≤ 10 −6 ). Finally, in the third simulation (Fig. 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 Law et al. [28] using the 3D-Var 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 simulations. First, in Sect. 3.1, we bound some constants in our theoretical results for the Lorenz 96' model. In Sect. 3.2, we explain the choice of the function F in our initial estimator (see Theorem 2.8) in the two observation scenarios. In Sect. 3.3, we adapt the Taylor expansion method for numerically solving ODEs to our setting. Finally, based on these preliminary results, we give the technical details of the simulations in Sect. 3.4.

Some Bounds for Lorenz 96' Model
In this section, we will bound some of the constants in Sect. 1 For B, given any u, v ∈ R d such that u , v ≤ 1, by the arithmetic mean-rootmean-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 definitions (1.16) and (1.6), we have

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.8 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 :

Notice that
so for m = 0, we have In general, for m ≥ 1, by differentiating (3.3) m times, we obtain that (3.5) 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−2 j , u d+1−2 j . 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, . . . , H D 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. [38] 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 "Appendix", we state a simple modification of the initial estimator of Theorem 2.8 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 system (1.1) and to implement Newton's method as described in (2.30).

Lemma 3.1 For any
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 be the projection of v on B R .
Proof Using (3.11) and the fact that the projection P B R decreases distances we know that for every j ∈ {0, 1, By inequality (1.6), it follows 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 Fig. 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. The running time of Algorithm 1 was approximately 1 s. The RMSEs were numerically approximated from 20 parallel runs. The parameters J (0) max and J (1) max of the initial estimator were chosen as 1.
In the case of Fig. 1 (d = 1,000,002, half of the coordinates observed), we could not use the same simulation technique as previously (finite difference approximation of the gradient and Hessian of g sm ) because of the huge computational and memory requirements. Instead, we have computed the Newton's method iterations described in (2.30) based on preconditioned conjugate gradient solver, with the gradient and the product of the Hessian with a vector evaluated based on adjoint methods as described by Eqs. (3.5)-(3.7) and Section 3.2.1 of Paulin et al. [37] (see also Le Dimet et al. [30]). This means that the Hessians were approximated using products of Jacobian matrices that were stored in sparse format due to the local dependency of Eq. (3.1). This efficient storage has allowed us to run Algorithm 1 in approximately 20-40 min in the simulations. We made 4 parallel runs to estimate the MSEs. The parameters J (0) max and J (1) max of the initial estimator were chosen as 2. Finally, in the case of Fig. 3 (d = 60, first 3 coordinates observed), we used the same method as in the first example (finite difference approximation of the gradient and Hessian of g sm ). The running time of Algorithm 1 was approximately 1 hour (in part due to the need of using arbitrary precision arithmetics with hundreds of digits of precision). The MSEs were numerically approximated from 2 parallel runs. The parameters J (0) max , . . . , J (19) max of the initial estimator were chosen as 24.

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)) If Assumption 2.1 holds for u, then for any The following lemma is useful for controlling the total variation distance of two distributions that are only known up to normalising constants.
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.13)) satisfies that for any y ∈ R d , 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 , (2) f,g (E). (4.5) It is easy to check that this distribution has marginals f and g, so by (2.13), and the fact that 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 "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.8) 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.7) 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. Proof Let M d denote the space of d × d complex matrices. By the chain rule, we can write the derivatives of the log-determinant as functions of the derivatives of the determinant, which were shown to exist in Bhatia and Jain [5]. Following their notation, we define the kth derivative of det at a point M ∈ M d as a map D k det M from (M d ) k to C with value They have defined the norm of kth derivative of the determinant as From Theorem 4 of Bhatia and Jain [5], it follows that for any k ≥ 1, the norm of the kth derivative can be bounded as 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 (1) TV (u, T ) and C (2) 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.15). 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 Eqs. (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 , (5.13) then based on the assumption that σ Z √ h ≤ C TV (u, T, ε) −1 , (5.2) and (5.3), it follows that In the following six lemmas, we bound the three terms in inequalities (5.16) and (5.17).

Lemma 5.1 Using the notations and assumptions of this section, we have
Proof Note that by (5.14), we know that B ρ ⊂ B R , and Now using (5.14), we can see that sup v∈B ρ log q(v) (1) q ≤ 1 2 , and using (5.9), we have Using the fact that |1 − exp(x)| ≤ 2|x| for −1 ≤ x ≤ 1, and the bounds | log(q(v)/q(u))| ≤ C (1) q v − u and (5.9), we can see that for every v ∈ B ρ , Let Z denote a d-dimensional standard normal random vector, then it is easy to see that Since we have assumed that the events in (5.6) and (5.7) hold, we know that thus, the result follows.

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 [46], 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∈B R 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 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 for By (5.18) and (5.14), it follows thatμ sm (v|Y 0:k ) μ sm G (v|Y 0:k ) ≤ 3 for every v ∈ B ρ . By using these and bounding μ sm (v|Y 0:k ) μ sm Based on (5.19), we have where Z is a d-dimensional standard normal random vector. The claimed result now follows using the fact that

Lemma 5.5 Using the notations and assumptions of this section, we have
Proof Using (5.22), Lemma 5.3, and the fact that and the result follows from v∈B c

Lemma 5.6 Using the notations and assumptions of this section, we have
for some finite positive constants C * 3 (u, T ), C * 4 (u, T ). Proof Let W be defined as in (5.10) and Z be d-dimensional standard normal. Using the fact that for a nonnegative-valued random variable X , we have

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|B R (·|Y 0:k ) as for any Borel-measurable S ⊂ R d . (5.24) 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|B R (·|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 [41], we have By Proposition 3(g) of Roberts and Rosenthal [41], there is a coupling (X 1 , X 2 ) of random vectors such that X 1 ∼ μ sm (·|Y 0:k ), X 2 ∼ μ sm G|B R (·|Y 0:k ), and P(X 1 = X 2 |Y 0:k ) = d TV (μ sm (·|Y 0:k ), μ sm G|B R (·|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 ( The statement of the lemma now follows by the triangle inequality.

Lemma 5.8 (Coupling argument for Wasserstein distance bound) The Wasserstein distance of the filtering distribution and its Gaussian approximation satisfies that
Proof By Theorem 4.1 of Villani [48], there exists a coupling of random variables (X 1 , X 2 ) (called the optimal coupling) such that X 1 ∼ μ sm (·|Y 0:k ), X 2 ∼ μ sm G (·|Y 0:k ), and . ] denote the projection of X 2 on the ball B R . Then using the fact that X 1 ∈ B R , it is easy to see that X 2 − X 1 ≤ X 2 − X 1 , and therefore, where L(X 2 |Y 0:k ) denotes the distribution ofX 2 conditioned on Y 0:k . 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|B R (·|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 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 26) then the density of μ fi G (·|Y 0:k ) can be written as (5.27) Since the normalising constant is not known for the case of η fi G (·|Y 0:k ), we define a rescaled versionη fi The following lemma bounds the difference between the logarithms ofη fi G (v|Y 0:k ) and μ fi G (v|Y 0:k ).

Lemma 5.9 For any
Proof of Lemma 5.9 By (5.27) and (5.28), we have The absolute value of the first difference can be bounded by Proposition 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 (2) TV (u, T ) ≥ C (2) TV (u, T ), and thus D TV (u, T, ε) ≥ C TV (u, T, ε). We also assume that D (1) 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 Eqs. (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 Based on (5.29), (5.30), and the assumption that The surface of a ball of radius R − u centred 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 (R − u ) exp(−GT ) centred at Ψ T (u) is contained in Ψ T (B R ), and thus by (5.32), B ρ ⊂ Ψ T (B R ). By Lemma 4.1, we have By Lemma 5.9, (5.32), and the fact that | exp(x) − 1| ≤ 2|x| for x ∈ [−1, 1], it follows that for v ∈ B ρ , we have 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 (1) TV (u, T ) and D (2) TV (u, T ) such that 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 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 notations. Let u G := u − A −1 k B k denote the centre of the Gaussian approximation μ sm G (defined when A k is positive definite). Let ρ(h, σ Z ) be as in (5.13), 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.
Proof This is a direct consequence of the Wasserstein distance bound of Theorem 2.1, since Proof From Proposition 4.1, and (4.3), it follows that for any 0 < ε ≤ 1, Sinceû sm MAP is the maximiser of log μ sm (v|Y 0:k ) on B R , our claim follows.
MAP , D 7 (u, T ) and D 8 (u, T ) such that for any 0 < ε ≤ 1, for MAP and S (2) MAP sufficiently large, we can assume that From Lemma 5.11, we know that Using (5.37), it follows that if the above event happens, then û sm MAP < R, and thus ∇ log(μ sm (û sm MAP |Y 0:k ))) = ∇l sm (û sm MAP ) − 2σ 2 Z ∇ log q(û sm MAP ) = 0. (5.39) Using the fact that ∇l sm , and Proposition 4.3, it follows that Moreover, by Lemma 4.4, we know that for any 0 < ε ≤ 1, we have By combining the four Eqs. (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.
Proof By applying Lemma 5.10 for ε = 0.1, we obtain that for If this event happens, then in particular, we have By the definition of B k in (2.8), it follows that conditioned on u, B k has d-dimensional multivariate normal distribution with covariance matrix 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, P Z ≥ √ d 2 ≥ 1 4 . (Indeed, if (Z (i) ) 1≤i≤d are i.i.d. standard normal random variables, then for any λ > 0, and the claim follows by applying Markov's inequality P( ) · 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.
Proof We define the event E k as 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 where in the last step we have used the Cauchy-Schwarz inequality. The above terms can be further bounded as follows. By Lemmas 4.4 and 4.3 it follows that 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 For bounding the term E 1 E k u sm − u G 2 u , we define t min := D 5 (u, T )σ 2 Z h, and and then by Lemma 5.10, it follows that for By writing and using the fact that 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 t min :=D 7 (u, T )σ 2 Z h, and which implies that for σ Z √ h < S M (u, T ), 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).
Proof This is a direct consequence of the Wasserstein distance bound of Theorem 2.2.
MAP , D 7 (u, T ) and D 8 (u, T ) such that for any 0 < ε ≤ 1, for Proof The proof is similar to the proof of Lemma 5.13. By applying Lemma 5.15 for (5.56) 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.
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.11). 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 ), 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 follows 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 [9] to our setting.) For v ∈ R d , r > 0, we denote the ball of radius r centred at v by B(v, r ) Suppose that the starting point x 0 ∈ Ω satisfies that x 0 − x * < min r * , 2 C H L H . 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 By combining this with the fact that ( 2C H x i − x * 2 for every i ∈ N, and by rearrangement, it follows that for every i ∈ N, hence the result.
The following proposition gives a lower bound on the Hessian near u. The proof is included in Section A.1 of "Appendix".
Then for any (5.61) The following proposition bounds the Lipschitz coefficient of the Hessian. The proof is included in Section A.1 of "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
Now we are ready to prove Theorem 2.5. , c(u, T )

Proof of Theorem 2.5 Let
, c(u, T ) min r H (u, T ), c(u,T ) 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, (5.63)

Initial Estimator
In this section, we will prove our results about the initial estimator that we have proposed in Sect. 2.3.
Proof of Proposition 2.2 By Taylor series expansion ofΦ (l| j max ) (Y 0:k ) with remainder term of order j max + 1, we obtain that where by (1.14) the remainder term R l, j max +1 can be bounded using the Cauchy-Schwarz inequality as Due to the particular choice of the coefficients of c (l| j max |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 ≥ 2 j max + 3. The concentration bound now follows directly from this bias bound, (5.20) and the fact that the estimatorΦ (l| j max ) (Y 0:k ) has Gaussian distribution with covariance matrix c (l| j max |k) · σ Z · I d o .
Proof of Lemma 2.6 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 j max can be written as If K j max would not be invertible, then its rows would be linearly dependent, that is, there would exist a nonzero vector α ∈ R j max +1 such that However, this is not possible, since by the fundamental theorem of algebra, the polynomial α i x i−1 can have at most j max roots, so it cannot be zero Lebesgue almost everywhere in [0, 1]. Therefore, K j max 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 By substituting this into (2.37), and applying some algebra, we obtain that g(l, l,k opt (l, l)) and the claim of the theorem now follows by substituting this into Proposition 2.2.

Conclusion
In this paper, we have proved the consistency and asymptotic efficiency results for MAP estimators for smoothing and filtering a class of partially observed nonlinear 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 [29,44]. 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 shown in Fig. 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. [37] for the nonlinear 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 nonlinear) in a computationally efficient way is a challenging problem for future research (see Bocquet et al. [7] for some examples where this situation arises). so (4.8) follows by bounds (A.4) and (A.5), and (4.9) follows by rearrangement.

Proof of Proposition 2.1
The proof is quite similar to the proof of Proposition 4.1 of Paulin et al. [38]. 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 : To complete the proof, we will now show (A.6). Using inequality (1.14), we can see that the Taylor expansion is valid for times 0 ≤ t < C −1 der . Based on this expansion, assuming that ih + δ < C −1 der , H D i v can be approximated by a finite difference formula depending on the values of Φ δ (v), Φh +δ (v), . . ., Φ ih+δ , with error of O(h). This finite difference formula will be denoted asΦ The coefficients a (i,δ) l are explicitly defined in Fornberg [20], and they only depend on i and the ratio δ/h. Based on the definition of these coefficients on page 700 of Fornberg [20] By definitions (A.7) and (A.8), it follows that (A.12) and thus (A.6) follows by rearrangement.
The following lemma bounds the number of balls of radius δ required to cover a d-dimensional unit ball. It will be used in the proofs of Proof For d ≥ 3, this follows from Theorem 1 of Dumer [18]. 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 δ contain 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.
This means that Dudley's entropy integral expectation bound (Corollary 13.2 of Boucheron et al. [8]) is applicable here. To apply that result, we first need to upperbound the packing number N (δ, T l ), which is the maximum number of points that can be selected in T l such that all of them are further away from each other than δ in and (A. 16) follows. The proof of (A.17) is similar.
Proof of Proposition 4.1 Using Assumption 2.1, we know that thus, it suffices to lower-bound the terms 2 k for (r, s) ∈ T 1 , r > 0. (T 1 was defined as in Lemma A.2.) We continuously extend it to r = 0 as Based on Lemma A.2, the lower bound of the random part can be obtained based on the upper bound on the quantity W 1 := sup (r,s)∈T Proof of Proposition 4.2 From the definitions, we have The first term in the right-hand side of the above inequality can be upper-bounded as where we have used the multivariate Taylor's expansion bound (1.19). Similarly to Lemma A.2, we define W 1 := sup (r,s)∈T 1 k i=0 ϕ i (r, s), Z i , and W 1 := sup (r,s)∈T 1 k i=0 −ϕ i (r, s), Z i , then the second term in (A.21) can be bounded as Based on (1.19), the partial derivatives of ϕ i (r, s) satisfy that  The claim of the proposition now follows by applying Lemma A.2 to W 1 and W 1 separately (for ε/2 instead of ε), and then using the union bound.

Proof of Proposition 4.3 From the definitions, we have
Using the multivariate Taylor's expansion bound (1.19), the first two in the right-hand side of the above inequality can be upper-bounded as We extend this continuously to r = 0 as We define W 2 as in Lemma A.2 as W 2 := sup (r,s 1 ,s 2 )∈T 2 k i=0 ϕ i (r, s 1 , s 2 ), Z i , and then the last term of (A.23) can be bounded as (A.25) By (1.19), for any (r, s 1 , s 2 ) ∈ T 2 , the partial derivatives of ϕ i satisfy that is convergent, is well defined and satisfies that Ψ t (Ψ −t (v)) = v and that for any i max ∈ N, Proof The result follows from bounds (1.14), and the definition of Eq. (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 estimateŝ Φ (l) of D l (u(t i r )) and then use the function F described in Sect. 3.2 to obtain initial estimators u(t i r ) of u(t i r ) for 0 ≤ r ≤ m. After this, we project these estimators to B R (see (3.13)), and run them backwards by t i r time units via 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 a posteriori probability μ sm (û r |Y 0:k ) (see (2.4)). Based on (2.8), 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 i r ) 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.