Optimization Based Methods for Partially Observed Chaotic Systems

TL;DR

This paper develops a theoretical framework for filtering and smoothing in partially observed chaotic systems, demonstrating asymptotic optimality of MAP estimators and providing efficient algorithms, with applications to high-dimensional Lorenz 96' models.

Contribution

It introduces a Gaussian approximation for filters and smoothers in chaotic systems, proving MAP estimator optimality and offering a fast Newton-based algorithm for high-dimensional models.

Findings

Gaussian approximation is valid for small h and noise variance

MAP estimators are asymptotically optimal as noise diminishes

Efficient Newton-based algorithm converges rapidly in high dimensions

Abstract

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.

Equations

$$\frac{d\mathbf{u}}{dt}=-\mathbf{A}\mathbf{u}-\mathbf{B}(\mathbf{u},\mathbf{u})+\mathbf{f},$$

$${\mathbf{D}}^i\mathbf{v}:={\frac{d^i\mathbf{v}(t)}{d{t}^i}|}_{t=0}\text{ for }i\in \mathbb{N},$$

$$⟨\mathbf{D}\mathbf{v},\mathbf{v}⟩\le 0\text{ for every }\mathbf{v}\in {\mathbb{R}}^d\text{ with }\Vert \mathbf{v}\Vert \in [R,R+\delta ].$$

$${\mathbf{Y}}_j:=\mathbf{H}\mathbf{u}(t_j)+{\mathbf{Z}}_j$$

$${\frac{1}{2}\frac{d}{dt}\Vert \mathbf{v}(t){\Vert }^2|}_{t=0}\le \frac{1}{2{\lambda }_{\mathbf{A}}}\Vert \mathbf{f}{\Vert }^2-\frac{{\lambda }_{\mathbf{A}}}{2}\Vert \mathbf{v}{\Vert }^2.$$

$$R\ge \frac{\Vert \mathbf{f}\Vert }{{\lambda }_{\mathbf{A}}},$$

$$\frac{d}{dt}(\mathbf{v}(t)-\mathbf{w}(t))=-\mathbf{A}(\mathbf{v}(t)-\mathbf{w}(t))-(\mathbf{B}(\mathbf{v}(t),\mathbf{v}(t)-\mathbf{w}(t))-\mathbf{B}(\mathbf{w}(t)-\mathbf{v}(t),\mathbf{w}(t))),$$

$$\exp (-Gt)\Vert \mathbf{v}-\mathbf{w}\Vert \le \Vert \mathbf{v}(t)-\mathbf{w}(t)\Vert \le \exp (Gt)\Vert \mathbf{v}-\mathbf{w}\Vert ,$$

$$\Vert \mathbf{A}\Vert :=\underset{\mathbf{v}\in {\mathbb{R}}^d:\Vert \mathbf{v}\Vert =1}{\sup }\Vert \mathbf{A}\mathbf{v}\Vert \text{ and }\Vert \mathbf{B}\Vert :=\underset{\mathbf{v},\mathbf{w}\in {\mathbb{R}}^d:\Vert \mathbf{v}\Vert =1,\Vert \mathbf{w}\Vert =1}{\sup }\Vert \mathbf{B}(\mathbf{v},\mathbf{w})\Vert .$$

$${\Phi }_t(\mathbf{v}):=\mathbf{H}{\Psi }_t(\mathbf{v}),$$

$${\mu }^{\mathrm{sm}}(\mathbf{v}|{\mathbf{Y}}_{0:k})=\left[\prod _{i=0}^k\eta \left({\mathbf{Y}}_i-{\Phi }_{t_i}(\mathbf{v})\right)\right]\frac{q(\mathbf{v})}{Z_k^{\mathrm{sm}}}\text{ for }\mathbf{v}\in {\mathcal{B}}_R,\text{ and }0\text{ for }\mathbf{v}\notin {\mathcal{B}}_R$$

$${\mu }^{\mathrm{fi}}(\mathbf{v}|{\mathbf{Y}}_{0:k})=\left[\prod _{i=0}^k\eta \left({\mathbf{Y}}_i-{\Phi }_{t_i-t_k}(\mathbf{v})\right)\right]\left|\det (\mathbf{J}{\Psi }_{-t_k}(\mathbf{v}))\right|\frac{q({\Psi }_{-t_k}(\mathbf{v}))}{Z_k^{\mathrm{fi}}}\text{ for }\mathbf{v}\in {\Psi }_{t_k}({\mathcal{B}}_R),$$

$$\text{and }0\text{ for }\mathbf{v}\notin {\Psi }_{t_k}({\mathcal{B}}_R),$$

$$\det (\mathbf{J}{\Psi }_{-t_k}(\mathbf{v}))={\left(\det ({\mathbf{J}}_{{\Psi }_{-t_k}(\mathbf{v})}{\Psi }_{t_k})\right)}^{-1}.$$

$$({\mathbf{J}}^{k}g(\mathbf{v}){)}_{i_1,\dots ,i_{k+1}}:=\frac{{\partial }^k}{\partial v_{i_1}\dots \partial v_{i_k}}{g}_{i_{k+1}}(\mathbf{v}),1\le i_1,\dots ,i_k\le d_1,1\le i_{k+1}\le d_2.$$

$$\Vert {\mathbf{J}}^{k}g(\mathbf{v})\Vert :=\underset{{\mathbf{v}}^{(1)}\in {\mathbb{R}}^{d_1},\dots ,{\mathbf{v}}^{(k)}\in {\mathbb{R}}^{d_1},{\mathbf{v}}^{(k+1)}\in {\mathbb{R}}^{d_2}:\Vert {\mathbf{v}}^{(j)}\Vert \le 1,1\le j\le k+1}{\sup }({\mathbf{J}}^{k}g(\mathbf{v}))[{\mathbf{v}}^{(1)},\dots ,{\mathbf{v}}^{(k+1)}],$$

$$\mathbf{M}[{\mathbf{v}}^{(1)},\dots ,{\mathbf{v}}^{(k+1)}]:=\sum _{1\le i_1\le d_1,\dots ,1\le i_{k+1}\le d_{k+1}}{M}_{i_1,\dots ,i_{k+1}}\cdot v_{i_1}^{(1)}\cdot \dots \cdot v_{i_{k+1}}^{(k+1)}.$$

$$\underset{\mathbf{v}\in {\mathcal{B}}_R,t\ge 0}{\sup }\Vert \frac{d\mathbf{v}(t)}{dt}\Vert \le v_{\max }:=\Vert \mathbf{A}\Vert R+\Vert \mathbf{B}\Vert R^2+\Vert \mathbf{f}\Vert ,$$

$$\underset{\mathbf{v}\in {\mathcal{B}}_R,t\ge 0}{\sup }\Vert {\mathbf{J}}_{\mathbf{v}(t)}\left(\frac{d\mathbf{v}(t)}{dt}\right)\Vert \le a_{\max }:=\Vert \mathbf{A}\Vert +2\Vert \mathbf{B}\Vert R.$$

$${\mathbf{D}}^i\mathbf{v}=-\mathbf{A}\cdot {\mathbf{D}}^{i-1}\mathbf{v}-\sum _{j=0}^{i-1}\left(\genfrac{}{}{0ex}{}{i-1}{j}\right)\mathbf{B}\left({\mathbf{D}}^j\mathbf{v},{\mathbf{D}}^{i-1-j}\mathbf{v}\right).$$

$$\Vert {\mathbf{D}}^i\mathbf{v}\Vert \le C_0{\left(C_{\mathrm{der}}\right)}^i\cdot i!,$$

$$\Vert {\mathbf{J}}_{\mathbf{v}}^k\left({\mathbf{D}}^i\mathbf{v}\right)\Vert \le {\left(C_{\mathbf{J}}^{(k)}\right)}^i\cdot i!,\text{ where}$$

$$C_0:=R+\frac{\Vert \mathbf{f}\Vert }{\Vert \mathbf{A}\Vert },C_{\mathrm{der}}:=\Vert \mathbf{A}\Vert +\Vert \mathbf{B}\Vert R+\frac{\Vert \mathbf{B}\Vert }{\Vert \mathbf{A}\Vert }\Vert \mathbf{f}\Vert ,C_{\mathbf{J}}^{(k)}:=2^k(C_{\mathrm{der}}+\Vert \mathbf{B}\Vert ),k\ge 1.$$

$$g(\mathbf{a}+\mathbf{h})=g(\mathbf{a})+\sum _{1\le j\le k}\frac{1}{j!}\cdot \left({\mathbf{J}}^{j}g(\mathbf{a})[{\mathbf{h}}^j,\cdot ]\right)+{\mathbf{R}}_{k+1}(\mathbf{a},\mathbf{h}),$$

$${\mathbf{R}}_{k+1}(\mathbf{a},\mathbf{h}):=\frac{k+1}{(k+1)!}\cdot {\int }_{t=0}^1(1-t{)}^k{\mathbf{J}}^{k+1}g(\mathbf{a}+t\mathbf{h})[{\mathbf{h}}^{k+1},\cdot ]dt,$$

$$\Vert {\mathbf{R}}_{k+1}(\mathbf{a},\mathbf{h})\Vert \le \frac{\Vert \mathbf{h}{\Vert }^{k+1}}{(k+1)!}\cdot \underset{0\le t\le 1}{\sup }\Vert {\mathbf{J}}^{k+1}g(\mathbf{a}+t\mathbf{h})\Vert .$$

$${\mathbf{J}}^k{\Psi }_t(\mathbf{v})={\mathbf{J}}_{\mathbf{v}}^k\left(\sum _{i=0}^{\infty }{\mathbf{D}}^i\mathbf{v}\cdot \frac{t^i}{i!}\right),$$

$${\mathbf{J}}^k{\Psi }_t(\mathbf{v})=\sum _{i=0}^{\infty }{\mathbf{J}}_{\mathbf{v}}^k\left({\mathbf{D}}^i\mathbf{v}\right)\cdot \frac{t^i}{i!}$$

$$D_{\mathbf{J}}^{(k)}:=2^k\left(\Vert \mathbf{A}\Vert +\Vert \mathbf{B}\Vert +2\Vert \mathbf{B}\Vert R\right).$$

$$M_k(T)$$

$${\hat{M}}_k(T)$$

$$\sum _{i=0}^k\Vert {\Phi }_{t_i}(\mathbf{v})-{\Phi }_{t_i}(\mathbf{u}){\Vert }^2\ge \frac{c(\mathbf{u},T)}{h}\Vert \mathbf{u}-\mathbf{v}{\Vert }^2.$$

$$\mathbf{H}{\mathbf{D}}^i\mathbf{u}=\mathbf{H}{\mathbf{D}}^i\mathbf{v}\text{ for every }0\le i\le j$$

Full text

Key words: Filtering; Smoothing; Chaotic Dynamical Systems; Gaussian Approximation; Newton’s Method; Concentration inequalities; 4D-Var.

**Optimization Based Methods for Partially Observed Chaotic Systems

**

BY DANIEL PAULIN1, AJAY JASRA1, DAN CRISAN2 & ALEXANDROS BESKOS3

1Department of Statistics & Applied Probability, National University of Singapore, Singapore, 117546, SG. E-Mail: paulindani@gmail.com, staja@nus.edu.sg

2Department of Mathematics, Imperial College London, London, SW7 2AZ, UK.

E-Mail: d.crisan@ic.ac.uk

3Department of Statistical Science, University College London, London, WC1E 6BT, UK. E-Mail: a.beskos@ucl.ac.uk

1 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 et al. (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

[TABLE]

where $\bm{u}:\mathbb{R}^{+}\to\mathbb{R}^{d}$ is a dynamical system in $\mathbb{R}^{d}$ for some $d\in\mathbb{Z}_{+}$, $\bm{A}$ is linear operator in $\mathbb{R}^{d}$ (i.e. $\bm{A}$ is a $d\times d$ matrix), $\bm{f}\in\mathbb{R}^{d}$ is a constant vector, and $\bm{B}(\bm{u},\bm{u})$ is a bilinear form corresponding to the nonlinearity (i.e. $\bm{B}$ is a $d\times d\times d$ array). We denote the solution of equation (1.1) with initial condition $\bm{u}(0)\mathrel{\mathop{:}}=\bm{v}$ for $t\geq 0$ by $\bm{v}(t)$. The derivatives of the solution $\bm{v}(t)$ at time $t=0$ will be denoted by

[TABLE]

in particular, $\bm{D}^{0}\bm{v}=\bm{v}$ , $\bm{D}\bm{v}\mathrel{\mathop{:}}=\bm{D}^{1}\bm{v}=-\bm{A}\bm{v}-\bm{B}(\bm{v},\bm{v})+\bm{f}$ (the right hand side of (1.1)), and $\bm{D}^{2}\bm{v}=-\bm{A}\bm{D}^{1}\bm{v}-\bm{B}(\bm{D}^{1}\bm{v},\bm{v})-\bm{B}(\bm{v},\bm{D}^{1}\bm{v})$.

In order to ensure the existence of a solution to the equation (1.1) for every $t\geq 0$, we assume that there are constants $R>0$ and $\delta>0$ such that

[TABLE]

We call this the trapping ball assumption. Let $\mathcal{B}_{R}\mathrel{\mathop{:}}=\{\bm{v}\in\mathbb{R}^{d}:\|\bm{v}\|\leq R\}$ be the ball of radius $R$. Using the fact that $\left<\frac{d}{dt}\bm{v}(t),\bm{v}(t)\right>=\frac{1}{2}\frac{d}{dt}\|\bm{v}(t)\|^{2}$, one can show that the solution to (1.1) exists for $t\geq 0$ for every $\bm{v}\in\mathcal{B}_{R}$, and satisfies that $\bm{v}(t)\in\mathcal{B}_{R}$ for $t\geq 0$.

Equation (1.1) was shown in Sanz-Alonso and Stuart (2015) 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 $\bm{A}$ and $\bm{B}(\bm{v},\bm{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,\ldots$, with observations

[TABLE]

where $\bm{H}:\mathbb{R}^{d}\to\mathbb{R}^{d_{o}}$ is a linear operator, and $(\bm{Z}_{j})_{j\geq 0}$ are i.i.d. centered random vectors taking values in $\mathbb{R}^{d_{o}}$ describing the noise. We assume that these vectors have distribution $\eta$ that is Gaussian with i.i.d. components of variance $\sigma_{Z}^{2}$. 111We believe that our results in this paper hold for non-Gaussian noise distributions as well, but proving this would be technically complex.

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 $\sigma^{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 $\sigma^{2}_{Z}h$ tends to [math]. 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 et al. (2013, 2012); Pillai et al. (2014); Vollmer (2015); Cui et al. (2016); Yao et al. (2016). There are also other randomisation and optimization based methods that have been recently proposed in the literature, see e.g. Bardsley et al. (2014); Wang et al. (2018).

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 et al. (2009); Dashti et al. (2013); Vollmer (2013); Helin and Burger (2015); Kekkonen et al. (2016); Monard et al. (2017); Nickl (2017); Dunlop and Stuart (2016). Some other important work on similar models and/or associated filtering/smoothing algorithms include Hayden et al. (2011); Blomker et al. (2013); Law et al. (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 et al. (2018) our key assumption on the dynamics was verified in 100 trials when $\bm{A}$, $\bm{B}$ and $\bm{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 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 (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.

1.1 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 $\Psi_{t}$, thus for a starting point $\bm{v}=(v_{1},\ldots,v_{d})\in\mathbb{R}^{d}$, the solution of (1.1) will be denoted by $\Psi_{t}(\bm{v})$, or equivalently, $\bm{v}(t)$. Sanz-Alonso and Stuart (2015) and Law, Stuart, and Zygalakis (2015) have assumed that the nonlinearity is energy conserving, i.e. $\left<\bm{B}(\bm{v},\bm{v}),\bm{v}\right>=0$ for every $\bm{v}\in\mathbb{R}^{d}$. They also assume that the linear operator $\bm{A}$ is positive definite, i.e. there is a $\lambda_{\bm{A}}>0$ such that $\left<\bm{A}\bm{v},\bm{v}\right>\geq\lambda_{\bm{A}}\left<\bm{v},\bm{v}\right>$ for every $\bm{v}\in\mathbb{R}^{d}$. As explained on page 50 of Law, Stuart, and Zygalakis (2015), (1.1) together with these assumptions above implies that for every $\bm{v}\in\mathbb{R}^{d}$,

[TABLE]

From (1.4) one can show that $\mathcal{B}_{R}$ is an absorbing set for any

[TABLE]

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\geq 0$, let $\bm{v}(t)$ and $\bm{w}(t)$ denote the solutions of (1.1) started from some points $\bm{v},\bm{w}\in\mathbb{R}^{d}$. Based on (1.1), we have that for any two points $\bm{v},\bm{w}\in\mathcal{B}_{R}$, any $t\geq 0$,

[TABLE]

and therefore by Grönwall’s lemma, we have that for any $t\geq 0$,

[TABLE]

for a constant $G\mathrel{\mathop{:}}=\|\bm{A}\|+2\|\bm{B}\|R$, where

[TABLE]

For $t\geq 0$, let $\Psi_{t}(\mathcal{B}_{R})\mathrel{\mathop{:}}=\{\Psi_{t}(\bm{v}):\bm{v}\in\mathcal{B}_{R}\}$, then by (1.6), it follows that $\Psi_{t}:\mathcal{B}_{R}\to\Psi_{t}(\mathcal{B}_{R})$ is a one-to-one mapping, which has an inverse that we denote as $\Psi_{-t}:\Psi_{t}(\mathcal{B}_{R})\to\mathcal{B}_{R}$.

The main quantities of interest of this paper are the smoothing and filtering distributions corresponding to the conditional distribution of $\bm{u}(t_{0})$ and $\bm{u}(t_{k})$, respectively, given the observations $\bm{Y}_{0:k}\mathrel{\mathop{:}}=\left\{\bm{Y}_{0},\ldots,\bm{Y}_{k}\right\}$. The densities of these distributions will be denoted by $\mu^{\mathrm{sm}}(\bm{v}|\bm{Y}_{0:k})$ and $\mu^{\mathrm{fi}}(\bm{v}|\bm{Y}_{0:k})$. To make our notation more concise, we define the observed part of the dynamics as

[TABLE]

for any $t\in\mathbb{R}$ and $\bm{v}\in\mathcal{B}_{R}$. Using these notations, the densities of the smoothing and filtering distributions can be expressed as

[TABLE]

where $\det$ stands for determinant, and $Z_{k}^{\mathrm{sm}},Z_{k}^{\mathrm{fi}}$ are normalising constants independent of $\bm{v}$. Since the determinant of the inverse of a matrix is the inverse of its determinant, we have the equivalent formulation

[TABLE]

We assume a prior $q$ on the initial condition that is absolutely continuous with respect to the Lebesgue measure, and zero outside the ball $\mathcal{B}_{R}$ (where the value of $R$ is determined by the trapping ball assumption (1.3)).

For $k\geq 1$, we define the $k$th Jacobian of a function $g:\mathbb{R}^{d_{1}}\to\mathbb{R}^{d_{2}}$ at point $\bm{v}$ as a $k+1$ dimensional array, denoted by $\bm{J}^{k}g(\bm{v})$ or equivalently $\bm{J}^{k}_{\bm{v}}g$ , with elements

[TABLE]

We define the norm of this $k$th Jacobian as

[TABLE]

where for a $k+1$ dimensional $d_{1}\times\ldots\times d_{k+1}$ sized array $\bm{M}$, we denote

[TABLE]

Using (1.1) and (1.3), we have that

[TABLE]

By induction, we can show that for any $i\geq 2$, and any $\bm{v}\in\mathbb{R}^{d}$, we have

[TABLE]

From this, the following bounds follow (see Section A.1 of the Appendix for a proof).

Lemma 1.1**.**

For any $i\geq 0$, $k\geq 1$, $\bm{v}\in\mathcal{B}_{R}$, we have

[TABLE]

In some of our arguments we are going to use the multivariate Taylor expansion for vector valued functions. Let $g:\mathbb{R}^{d_{1}}\to\mathbb{R}^{d_{2}}$ be $k+1$ times differentiable for some $k\in\mathbb{N}$. Then using the one dimensional Taylor expansion of the functions $g_{i}(\bm{a}+t\bm{h})$ in $t$ (where $g_{i}$ denotes the $i$th component of $g$), one can show that for any $\bm{a},\bm{h}\in\mathbb{R}^{d_{1}}$, we have

[TABLE]

where $\bm{h}^{j}\mathrel{\mathop{:}}=(\bm{h},\ldots,\bm{h})$ denotes the $j$ times repetition of $\bm{h}$, and the error term $\bm{R}_{k+1}(\bm{a},\bm{h})$ is of the form

[TABLE]

whose norm can be bounded using the fact that $\int_{t=0}^{1}(1-t)^{k}dt=\frac{1}{k+1}$ as

[TABLE]

In order to be able to use such multivariate Taylor expansions in our setting, the existence and finiteness of $\bm{J}^{k}\Psi_{t}(\bm{v})$ can be shown rigorously in the following way. Firstly, for $0<t<\left(C_{\bm{J}}^{(k)}\right)^{-1}$, one has

[TABLE]

and using the inequality (1.15), we can show that

[TABLE]

is convergent and finite. For $t\geq\left(C_{\bm{J}}^{(k)}\right)^{-1}$, we can express $\Psi_{t}(\bm{v})$ as a composition $\Psi_{t_{1}}(\ldots(\Psi_{t_{m}}(\bm{v})))$ for $t_{1}+\cdots+t_{m}=t$, and establish the existence of the partial derivatives by the chain rule.

After establishing the existence of the partial derivatives $\bm{J}^{k}\Psi_{t}(\bm{v})$, we are going to bound their norm in the following lemma (proven in Section A.1 of the Appendix).

Lemma 1.2**.**

For any $k\geq 1$, let

[TABLE]

Then for any $k\geq 1$, $T\geq 0$, we have

[TABLE]

2 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 $\sigma_{Z}^{2}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 $\sigma_{Z}^{2}h\to 0$ limit. We also show that Newton’s method can be used for calculating the MAP estimators if the initial point $\bm{x}_{0}$ can be chosen sufficiently close to the true starting position $\bm{u}$. Finally, in Section 2.3 we propose estimators to use as initial point $\bm{x}_{0}$ that satisfy this criteria when $\sigma_{Z}^{2}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\in\mathbb{N}$. Suppose that $\|\bm{u}\|<R$, and that there exist constants $h_{\max}(\bm{u},T)>0$ and $c(\bm{u},T)>0$ such that for every $\bm{v}\in\mathcal{B}_{R}$, for every $h\leq h_{\max}(\bm{u},T)$ (or equivalently, every $k\geq T/h_{\max}(\bm{u},T)$), we have

[TABLE]

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 $\|\bm{u}\|<R$, and there is an index $j\in\mathbb{N}$ such that the system of equations in $\bm{v}$ defined as

[TABLE]

has a unique solution $\bm{v}:=\bm{u}$ in $\mathcal{B}_{R}$, and

[TABLE]

where $\left(\bm{H}\bm{D}^{i}\bm{u}\right)_{k}$ refers to coordinate $k$ of the vector $\bm{H}\bm{D}^{i}\bm{u}\in\mathbb{R}^{d_{o}}$, and $\nabla$ denotes the gradient of the function in $\bm{u}$.

Proposition 2.1**.**

Assumption 2.2 implies Assumption 2.1.

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 $\bm{H}$), thus Assumption 2.1 is also valid for these models.

We denote $\mathrm{int}(\mathcal{B}_{R})\mathrel{\mathop{:}}=\{\bm{v}\in\mathbb{R}^{d}:\|\bm{v}\|<R\}$ the interior of $\mathcal{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 $\mathbb{R}^{d}$, supported on $\mathcal{B}_{R}$. We assume that $v\to q(\bm{v})$ is strictly positive and continuous on $\mathcal{B}_{R}$, and that it is 3 times continuously differentiable at every interior point of $\mathcal{B}_{R}$. Let

[TABLE]

We assume that these are finite.

After some simple algebra, the smoothing distribution for an initial point $\bm{u}$ and the filtering distribution for the current position $\bm{u}(T)$ can be expressed as

[TABLE]

where $C_{k}^{\mathrm{sm}}$ is a normalising constant independent of $\bm{v}$ (but depending on $(\bm{Z}_{j})_{j\geq 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 $\bm{u}$ based on the observations $\bm{Y}_{0},\ldots,\bm{Y}_{k}$, to be used as a starting point for Newton’s method.

2.1 Gaussian approximation

We define the matrix $\bm{A}_{k}\in\mathbb{R}^{d\times d}$ and vector $\bm{B}_{k}\in\mathbb{R}^{d}$ as

[TABLE]

where $\bm{J}\Phi_{t_{i}}$ and $\bm{J}^{2}\Phi_{t_{i}}$ denotes the first and second Jacobian of $\Phi_{t_{i}}$, respectively, and $\bm{J}^{2}\Phi_{t_{i}}(\bm{u})[\cdot,\cdot,\bm{Z}_{i}]$ denotes the $d\times d$ matrix with elements

[TABLE]

If $\bm{A}_{k}$ is positive definite, then we define the center of the Gaussian approximation of the smoother as

[TABLE]

and define the Gaussian approximation of the smoother as

[TABLE]

If $\bm{A}_{k}$ is not positive definite, then we define the Gaussian approximation of the smoother $\mu^{\mathrm{sm}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ to be the $d$-dimensional standard normal distribution (an arbitrary choice), and $\bm{u}^{\mathcal{G}}:=\bm{0}$. If $\bm{A}_{k}$ is positive definite, and $\bm{u}^{\mathcal{G}}\in\mathcal{B}_{R}$, then we define the Gaussian approximation of the filter as

[TABLE]

Alternatively, if $\bm{A}_{k}$ is not positive definite, or $\bm{u}^{\mathcal{G}}\notin\mathcal{B}_{R}$, then we define the Gaussian approximation of the smoother $\mu^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{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 $\mu_{1},\mu_{2}$ on $\mathbb{R}^{d}$ that are absolutely continuous with respect to the Lebesgue measure is defined as

[TABLE]

where $\mu_{1}(x)$, and $\mu_{2}(x)$ denote the densities of the distributions.

The Wasserstein distance (also called 1st Wasserstein distance) of two distributions $\mu_{1},\mu_{2}$ on $\mathbb{R}^{d}$ (with respect to the Euclidean distance) is defined as

[TABLE]

where $\Gamma(\mu_{1},\mu_{2})$ is the set of all measures on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $\mu_{1}$ and $\mu_{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}(\bm{u},T)$ appears. For brevity, we denote this as

[TABLE]

We are also going to use the constant $C_{\|\bm{A}\|}$ defined as

[TABLE]

Theorem 2.1** (Gaussian approximation of the smoother).**

Suppose that Assumptions 2.1 and 2.3 hold for the initial point $\bm{u}$ and the prior $q$. Then there are constants $C^{(1)}_{\mathrm{TV}}(\bm{u},T)$, $C^{(2)}_{\mathrm{TV}}(\bm{u},T)$, $C^{(1)}_{\mathrm{W}}(\bm{u},T)$, and $C^{(2)}_{\mathrm{W}}(\bm{u},T)$ independent of $\sigma_{Z}$, $h$ and $\varepsilon$ such that for any $0<\varepsilon\leq 1$, $\sigma_{Z}>0$, and $0\leq h\leq h_{\max}(\bm{u},T)$ satisfying that $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, we have

[TABLE]

[TABLE]

Theorem 2.2** (Gaussian approximation of the filter).**

Suppose that Assumptions 2.1 and 2.3 hold for the initial point $\bm{u}$ and the prior $q$. Then there are constants $D^{(1)}_{\mathrm{TV}}(\bm{u},T)$, $D^{(2)}_{\mathrm{TV}}(\bm{u},T)$, $D^{(1)}_{\mathrm{W}}(\bm{u},T)$, and $D^{(2)}_{\mathrm{W}}(\bm{u},T)$ independent of $\sigma_{Z}$, $h$ and $\varepsilon$ such that for any $0<\varepsilon\leq 1$, $\sigma_{Z}>0$, and $0\leq h\leq h_{\max}(\bm{u},T)$ satisfying that $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}D_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, we have

[TABLE]

[TABLE]

Note that the Gaussian approximations $\mu^{\mathrm{sm}}_{\mathcal{G}}$ and $\mu^{\mathrm{fi}}_{\mathcal{G}}$ as defined above are not directly computable based on the observations $\bm{Y}_{0:k}$, since they involve the true initial position $\bm{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 $\bm{Y}_{0:k}$.

2.2 MAP estimators

Let $\overline{\bm{u}}^{\mathrm{sm}}$ be the mean of the smoothing distribution, and $\overline{\bm{u}}^{\mathrm{fi}}$ be the mean of the filtering distribution, and $\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}$ be the maximum-a-posteriori of the smoothing distribution, i.e.

[TABLE]

In case there are multiple maxima, we choose any of them. For the filter, we will use the push-forward MAP estimator

[TABLE]

Based on the Gaussian approximation results, we prove the following two theorems about these estimators.

Theorem 2.3** (Comparison of mean square error of MAP and posterior mean for smoother).**

Suppose that Assumptions 2.1 and 2.3 hold for the initial point $\bm{u}$ and the prior $q$. Then there is a constant $S^{\mathrm{sm}}_{\max}(\bm{u},T)>0$ independent of $\sigma_{Z}$ and $h$ such that for $0<h<h_{\max}(\bm{u},T)$, $\sigma_{Z}\sqrt{h}\leq S^{\mathrm{sm}}_{\max}(\bm{u},T)$, we have that

[TABLE]

where the expectations are taken with respect to the random observations, and $\underline{C}^{\mathrm{sm}}(\bm{u},T)$, $\overline{C}^{\mathrm{sm}}(\bm{u},T)$, and $C^{\mathrm{sm}}_{\mathrm{MAP}}(\bm{u},T)$ are finite positive constants independent of $\sigma_{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 $\bm{u}$ and the prior $q$. Then there is a constant $S^{\mathrm{fi}}_{\max}(\bm{u},T)>0$ independent of $\sigma_{Z}$ and $h$ such that for $0<h<h_{\max}(\bm{u},T)$, $\sigma_{Z}\sqrt{h}\leq S^{\mathrm{fi}}_{\max}(\bm{u},T)$, we have that

[TABLE]

where the expectations are taken with respect to the random observations, and $\underline{C}^{\mathrm{fi}}(\bm{u},T)$, $\overline{C}^{\mathrm{fi}}(\bm{u},T)$, and $C^{\mathrm{fi}}_{\mathrm{MAP}}(\bm{u},T)$ are finite positive constants independent of $\sigma_{Z}$ and $h$.

Remark 2.1**.**

Theorems 2.3 and 2.4 in particular imply that when $T$ is fixed, and $\sigma_{Z}\sqrt{h}$ tends to [math], the ratio between the mean square errors of the posterior mean and MAP estimators conditioned on the initial position $\bm{u}$ tends to 1. Since the mean of the posterior distributions, $\overline{\bm{u}}^{\mathrm{sm}}$ (or $\overline{\bm{u}}^{\mathrm{fi}}$ for the filter), is the estimator $U(\bm{Y}_{0:k})$ that minimises $\mathbb{E}(\|U(\bm{Y}_{0:k})-\bm{u}\|^{2})$ (or $\mathbb{E}(\|U(\bm{Y}_{0:k})-\bm{u}(T)\|^{2})$ for the filter), our results imply that the mean square error of the MAP estimators is close to optimal when $\sigma_{Z}\sqrt{h}$ is sufficiently small.

Next we propose a method to compute the MAP estimators. Let $g^{\mathrm{sm}}:\mathcal{B}_{R}\to\mathbb{R}$ be

[TABLE]

Then $-g^{\mathrm{sm}}(\bm{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 $\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}$ to arbitrary precision if it is initiated from a starting point $\bm{x}_{0}$ that is sufficiently close to the initial position $\bm{u}$. The proof is based on the concavity properties of the log-likelihood near $\bm{u}$. Based on this, an approximation for the push-forward MAP estimator $\hat{\bm{u}}^{\mathrm{fi}}$ can be then computed by moving forward the approximation of $\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}$ by time $T$ according to the dynamics $\Psi_{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 $\bm{u}$ and the prior $q$. Then for every $0<\varepsilon\leq 1$, there exist finite constants $S_{\max}^{\mathrm{sm}}(\bm{u},T,\varepsilon)$, $N^{\mathrm{sm}}(\bm{u},T)$ and $D_{\max}^{\mathrm{sm}}(\bm{u},T)\in(0,N^{\mathrm{sm}}(\bm{u},T)]$ (defined in (5.62) and (5.63)) such that the following holds. If $\sigma_{Z}\sqrt{h}\leq S_{\max}^{\mathrm{sm}}(\bm{u},T,\varepsilon)$, and the initial point $\bm{x}_{0}\in\mathcal{B}_{R}$ satisfies that $\|\bm{x}_{0}-\bm{u}\|<D_{\max}^{\mathrm{sm}}(\bm{u},T)$, then the iterates of Newton’s method defined recursively as

[TABLE]

satisfy that

[TABLE]

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 $\bm{x}_{0}$ is sufficiently close to $\bm{u}$.

2.3 Initial estimator

First we are going to estimate the derivatives $\bm{H}\bm{D}^{l}\bm{u}$ for $l\in\mathbb{N}$ based on observations $\bm{Y}_{0:k}$. For technical reasons, the estimators will depend on $\bm{Y}_{0:\hat{k}}$ for some $0\leq\hat{k}\leq k$ (which will be chosen depending on $l$). For any $j\in\mathbb{N}$, we define $\bm{v}^{(j|\hat{k})}\in\mathbb{R}^{\hat{k}+1}$ as

[TABLE]

For $j_{\max}\in\mathbb{N}$, we define $\bm{M}^{(j_{\max}|\hat{k})}\in\mathbb{R}^{(j_{\max}+1)\times(\hat{k}+1)}$ as a matrix with rows $\bm{v}^{(0|\hat{k})},\ldots,\bm{v}^{(j_{\max}|\hat{k})}$.

We denote by $\bm{I}_{j_{\max}+1}$ the identity matrix of dimension $j_{\max}+1$, and by $\bm{e}^{(l|j_{\max})}$ a column vector in $\mathbb{R}^{j_{\max}+1}$ whose every component is zero except the $l+1$th one which is 1.

For any $l\in\mathbb{N}$, $j_{\max}\geq l$, $\hat{k}\geq j_{\max}$, we define the vector

[TABLE]

is an estimator of $\bm{H}\bm{D}^{l}\bm{u}$. The fact that the matrix $\bm{M}^{(j_{\max}|\hat{k})}(\bm{M}^{(j_{\max}|\hat{k})})^{\prime}$ is invertible follows from the fact that $\bm{v}^{(0|\hat{k})},\ldots,\bm{v}^{(j_{\max}|\hat{k})}$ are linearly independent (since the matrix with rows $\bm{v}^{(0|\hat{k})},\ldots,\bm{v}^{(\hat{k}|\hat{k})}$ is a so-called Vandermonde matrix whose determinant is non-zero). From (2.31), it follows that the norm of $\bm{c}^{(l|j_{\max}|\hat{k})}$ can be expressed as

[TABLE]

To lighten the notation, for $j_{\max}\geq l$ and $\hat{k}\geq j_{\max}$, we will denote

[TABLE]

The next proposition gives an error bound for this estimator, which we will use for choosing the values $\hat{k}$ and $j_{\max}$ given $l$.

Proposition 2.2**.**

Suppose that $j_{\max}\geq l$ and $\hat{k}\geq 2j_{\max}+3$. Then for any $0<\varepsilon\leq 1$,

[TABLE]

where

[TABLE]

The following lemma shows that as $\hat{k}\to\infty$, the constant $C_{\bm{M}}^{(l|j_{\max}|\hat{k})}$ tends to a limit.

Lemma 2.1**.**

Let $\bm{K}^{(j_{\max})}\in\mathbb{R}^{j_{\max}+1\times j_{\max}+1}$ be a matrix with elements $\bm{K}^{(j_{\max})}_{i,j}:=\frac{1}{i+j-1}$ for $1\leq i,j\leq j_{\max}+1$. Then for any $l\in\mathbb{N}$, $j_{\max}\geq l$, the matrix $\bm{K}^{(j_{\max})}$ is invertible, and

[TABLE]

The proofs of the above two results are included in Section 5.4. Based on these results, we choose $\hat{k}\in\{2j_{\max}+3,\ldots,k\}$ such that the function $g(l,j_{\max},\hat{k})$ is minimised. We denote this choice of $\hat{k}$ by $\hat{k}_{\mathrm{opt}}(l,j_{\max})$ (if $g(l,j_{\max},\hat{k})$ takes the same value for several $\hat{k}$, then we choose the smallest of them). By taking the derivative of $g(l,j_{\max},\hat{k})$ in $\hat{k}$, it is easy to see that it has a single minimum among positive real numbers achieved at

[TABLE]

Based on this, we have

[TABLE]

Finally, based on the definition of $\hat{k}_{\mathrm{opt}}(l,j_{\max})$, we choose $j_{\max}^{\mathrm{opt}}(l)$ as

[TABLE]

where $J_{\max}^{(l)}\in\{l,l+1,\ldots,\lfloor(k-3)/2\rfloor\}$ 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 $\bm{H}\bm{D}^{l}\bm{u}$ as

[TABLE]

The following theorem bounds the error of this estimator.

Theorem 2.6**.**

Suppose that $\bm{u}\in\mathcal{B}_{R}$, and $T=kh$. Then for any $l\in\mathbb{N}$, there exist some positive constants $h_{\max}^{(l)}$, $s_{\max}^{(l)}(T)$ and $S_{\max}^{(l)}(T)$ such that for any choice of the parameter $J_{\max}^{(l)}\in\{l,l+1,\ldots,\lfloor(k-3)/2\rfloor\}$, any $\varepsilon>0$, $0<s\leq s_{\max}^{(l)}(T)$, $0<h\leq\frac{h_{\max}^{(l)}\cdot s}{\sqrt{1+\frac{\log\left(1/\varepsilon\right)}{\log(d_{o}+1)}}}$, $0\leq\sigma_{Z}\sqrt{h}\leq S_{\max}^{(l)}(T)\cdot\left(\frac{s}{\sqrt{1+\frac{\log\left(1/\varepsilon\right)}{\log(d_{o}+1)}}}\right)^{l+3/2}$,

[TABLE]

The following theorem proposes a way of estimating $\bm{u}$ from estimates for the derivatives $\left(\bm{H}\bm{D}^{i}\bm{u}\right)_{0\leq i\leq j}$. This will be used as our initial estimator for Newton’s method.

Theorem 2.7**.**

Suppose that for some $j\in\mathbb{N}$ there is a function $F:(\mathbb{R}^{d_{o}})^{j+1}\to\mathbb{R}^{d}$ independent of $\bm{u}$ such that

[TABLE]

for $\sum_{i=0}^{j}\left\|\bm{H}\bm{D}^{i}\bm{u}-\bm{x}^{(i)}\right\|^{2}\leq D_{F}(\bm{u})$, for some positive constants $C_{F}(\bm{u})$, $D_{F}(\bm{u})$. Then the estimator $F(\hat{\Phi}^{(0)},\ldots,\hat{\Phi}^{(j)})$ satisfies that if $\sum_{i=0}^{j}\left\|\bm{H}\bm{D}^{i}\bm{u}-\hat{\Phi}^{(i)}\right\|^{2}\leq D_{F}(\bm{u})$, then

[TABLE]

In particular, under Assumption 2.2 and for $j$ as determined therein, function $F$ defined as

[TABLE]

satisfies conditions (2.41) and (2.42).

Thus, the initial estimator can simply be chosen as

[TABLE]

with the above two theorems implying that the estimate gets close to $u$ for decreasing $\sigma_{Z}\sqrt{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 et al. (2002), Tütüncü et al. (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.

2.4 Optimisation based smoothing and filtering

The following algorithm provides an estimator of $\bm{u}$ given $\bm{Y}_{0:k}$. We assume that either there is a problem-specific $F$ satisfying conditions (2.41) and (2.42) for some $j\in\mathbb{N}$, or we suppose that Assumption 2.2 is satisfied for the true initial point $\bm{u}$, and use $F$ as defined in (2.44).

The following algorithm returns an online estimator of $\left(\bm{u}(t_{i})\right)_{i\geq 0}$ given $\bm{Y}_{0:i}$ at time $t_{i}$.

This algorithm can be modified to run Step 2 only at every $K$ step for some $K\in\mathbb{Z}_{+}$ (i.e. for $i=k+lK$ for $l\in\mathbb{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 $\Delta_{\min}$ is chosen sufficiently small then the estimation errors are of $O(\sigma_{Z}\sqrt{h})$ with high probability for both algorithms (in the second algorithm, for $i\geq k$).

3 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

[TABLE]

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 et al. (2010)). As shown on page 16 of Sanz-Alonso and Stuart (2015), this system can be written in the form (1.1), and the bilinear form $\bm{B}(\bm{u},\bm{u})$ satisfies the energy conserving property (i.e. $\left<\bm{B}(\bm{v},\bm{v}),\bm{v}\right>=0$ for any $\bm{v}\in\mathbb{R}^{d}$).

We consider 2 observation scenarios for this model. In the first scenario, we assume that $d$ is divisible by 6, and choose $\bm{H}$ such that coordinates $1,2,3$, $7,8,9$, $\ldots$, $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 $\sigma_{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 $\sigma_{Z}\sqrt{h}$.

In the second scenario, we choose $\bm{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 $\bm{u}\in\mathbb{R}^{d}$. A simple modification of that argument shows that Assumption 2.2 holds for Lebesgue-almost every initial point $\bm{u}\in\mathbb{R}^{d}$ in the first observation scenario too.

In each case, we have set the initial point as $\bm{u}=\left(\frac{d+1}{2d},\frac{d+2}{2d},\ldots,1\right)$ (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 $\sigma_{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\cdot 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 $\sigma_{Z}\sqrt{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\cdot 10^{-5}$, while the triangles correspond to data points with $h=2.5\cdot 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 $\sigma_{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 $\bm{H}$. In the first simulation (Figure 1), relatively large noise ($\sigma_{Z}=10^{-3}$) and large time step ($h=10^{-2}$) were possible. For the second simulation (Figure 1), due to the high dimensionality of the model ($d=1000002$), and the sparse approximation used in the solver, we had to choose smaller time step ($h\leq 10^{-6}$) and smaller observation noise ($\sigma_{Z}\leq 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 ($\sigma_{Z}\leq 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. (2016) 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 simulations. First, in Section 3.1, we bound some constants in our theoretical results for the Lorenz 96’ model. In Section 3.2, we explain the choice of the function $F$ in our initial estimator (see Theorem 2.7) in the two observation scenarios. In Section 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 Section 3.4.

3.1 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. Since $\bm{A}$ is the identity matrix, we have $\|\bm{A}\|=1$ and satisfies that $\left<\bm{v},\bm{A}\bm{v}\right>\geq\lambda_{\bm{A}}\bm{v}$ for every $\bm{v}\in\mathbb{R}^{d}$ for $\lambda_{\bm{A}}=1$. The condition that $\left<\bm{B}(\bm{v},\bm{v}),\bm{v}\right>=0$ for every $\bm{v}\in\mathbb{R}^{d}$ was verified for the Lorenz 96’ model, see Property 3.1 of Law et al. (2016). For our choice $f=8$, we have $\|\bm{f}\|=8\sqrt{d}$. Thus based on (1.5), the trapping ball assumption (1.3) is satisfied for the choice $R:=\frac{\|\bm{f}\|}{\lambda_{\bm{A}}}=8\sqrt{d}$.

For $\bm{B}$, given any $\bm{u},\bm{v}\in\mathbb{R}^{d}$ such that $\|\bm{u}\|,\|\bm{v}\|\leq 1$, by the arithmetic-mean-root mean square inequality, and the inequality $(ab)^{2}\leq\frac{a^{2}+b^{2}}{2}$, we have

[TABLE]

thus $\|\bm{B}\|\leq 2$. If $d$ is divisible by 2, then the choice $u_{i}=v_{i}=\frac{(-1)^{i}}{\sqrt{d}}$ shows that this bound is sharp, and $\|\bm{B}\|=2$. For simplicity, we have chosen the prior $q$ as the uniform distribution on $\mathcal{B}_{R}$. Based on these, and the definitions (1.16) and (1.6), we have

[TABLE]

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=\left\lceil\frac{d-3}{3}\right\rceil$, it is possible to construct a function $F:(\mathbb{R}^{d_{o}})^{j+1}\to\mathbb{R}^{d}$ such that $F$ is computationally simple, Lipschitz in a neighbourhood of $\bm{u}$, and $F\left(\bm{H}\bm{u},\ldots,\bm{H}\bm{D}^{j}\bm{v}\right)=\bm{u}$, thus satisfies the conditions of Theorem 2.7. Notice that

[TABLE]

so for $m=0$, we have

[TABLE]

In general, for $m\geq 1$, by differentiating (3.3) $m$ times, we obtain that

[TABLE]

thus for any $m\geq 1$,

[TABLE]

Thus for any $m\in\mathbb{N}$, we have a recursion for the $m$th 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,

[TABLE]

therefore for $m=0$, we have

[TABLE]

By differentiating (3.7) $m$ times, we obtain that

[TABLE]

thus

[TABLE]

Thus for any $m\in\mathbb{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:=\left\lceil\frac{d-3}{3}\right\rceil$, these recursions together define a function $F$ for this case. From the recursion formulas, and the boundedness of $\bm{u}$ it follows that $F$ is Lipschitz in a neighbourhood of $\left(\bm{H}\bm{u},\ldots,\bm{H}\bm{D}^{j}\bm{u}\right)$ as long as none of the coordinates of $\bm{u}$ is 0 (thus for Lebesgue-almost every $\bm{u}\in\mathcal{B}_{R}$).

Now we look at the first observation scenario, i.e. suppose that $d$ is divisible by 6, and we observe coordinates $(6i+1,6i+2,6i+3)_{0\leq i\leq d/6-1}$. In this case, we choose $j=1$, and define $F\left(\bm{H}\bm{u},\bm{H}\bm{D}\bm{u}\right)$ based on formulas (3.4) and (3.8), so that we can express $u_{6i+4},u_{6i-1},u_{6i-2}$ based on $u_{6i+1},u_{6i+2},u_{6i+3}$ and $\bm{D}u_{6i+1},\bm{D}u_{6i+2},\bm{D}u_{6i+3}$ for $0\leq i\leq d/6-1$ (with indices counted modulo $d$). Based on the equations (3.4) and (3.8), we can see that $F$ defined as above satisfies the conditions of Theorem 2.7 as long as none of the coordinates of $u$ is 0 (thus for Lebesgue-almost every $\bm{u}\in\mathcal{B}_{R}$). In both scenarios, $F$ is computationally simple.

We note that in the above argument, $F$ might be not defined if some of the components of $\bm{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 $\bm{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 $\bm{u}$ are zero.

3.3 Numerical solution of chaotic ODEs based on Taylor expansion

Let $\bm{v}\in\mathcal{B}_{R}$, and $i_{\max}\in\mathbb{N}$. The following lemma provides some simple bounds that allow us to approximate the quantities $\bm{v}(t)=\Psi_{t}(\bm{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).

Lemma 3.1**.**

For any $\bm{v}\in\mathcal{B}_{R}$, we have

[TABLE]

where $\bm{D}^{0}=\bm{v}$, and for $i\geq 1$, we have the recursion

[TABLE]

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 $\bm{v}(t)$ for larger values $t$ (the bound (3.11) is not useful for $t\geq\frac{1}{C_{\mathrm{der}}}$). For $\bm{v}\in\mathbb{R}^{d}$, we let

[TABLE]

be the projection of $\bm{v}$ on $\mathcal{B}_{R}$.

Proposition 3.1**.**

Let $\bm{v}\in\mathcal{B}_{R}$, and $\Delta<\frac{1}{C_{\mathrm{der}}}$. Let $\widehat{\bm{v}}(0)=\bm{v}$, $\widehat{\bm{v}}(\Delta):=P_{\mathcal{B}_{R}}\left(\sum_{i=0}^{i_{\max}}\frac{\Delta^{i}}{i!}\bm{D}^{i}\bm{v}\right)$, and similarly, given $\widehat{\bm{v}}(j\Delta)$, define $\widehat{\bm{v}}((j+1)\Delta):=P_{\mathcal{B}_{R}}\left(\sum_{i=0}^{i_{\max}}\frac{\Delta^{i}}{i!}\bm{D}^{i}(\widehat{\bm{v}}(j\Delta))\right).$ Finally, let $\delta:=t-\lfloor\frac{t}{\Delta}\rfloor\Delta$ and $\widehat{\bm{v}}(t):=P_{\mathcal{B}_{R}}\left(\sum_{i=0}^{i_{\max}}\frac{\delta^{i}}{i!}\bm{D}^{i}(\widehat{\bm{v}}(\lfloor\frac{t}{\Delta}\rfloor\Delta))\right).$ Then the error is bounded as

[TABLE]

Proof.

Using (3.11) and the fact that the projection $P_{\mathcal{B}_{R}}$ decreases distances we know that for every $j\in\left\{0,1,\ldots,\lceil t/\Delta\rceil-1\right\}$, $\|\widehat{\bm{v}}((j+1)\Delta)-\Psi_{\Delta}(\widehat{\bm{v}}(j\Delta))\|\leq C_{0}\left(C_{\mathrm{der}}\Delta\right)^{i_{\max}+1}$. By inequality (1.6), it follows that

[TABLE]

and using the triangle inequality, we have

[TABLE]

so the claim follows. ∎

3.4 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 $\frac{1}{C_{\mathrm{der}}}$, so we have used the method of Proposition 3.1 to simulate from the system. The gradient and Hessian of the function $g^{\mathrm{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 second. The RMSEs were numerically approximated from 20 parallel runs. The parameters $J_{\max}^{(0)}$ and $J_{\max}^{(1)}$ of the initial estimator were chosen as 1.

In the case of Figure 1 ($d=1000002$, 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^{\mathrm{sm}}$) because of the huge computational and memory requirements. Instead, we have computed the Newton’s method iterations described in (2.28) based on preconditioned conjugate gradient solver, with the gradient and the product of the Hessian with a vector were evaluated based on adjoint methods as described by equations (3.5)-(3.7) and Section 3.2.1 of Paulin et al. (2017) (see also Le Dimet et al. (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 parameters $J_{\max}^{(0)}$ and $J_{\max}^{(1)}$ of the initial estimator were chosen as 2.

Finally, in the case of Figure 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^{\mathrm{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_{\max}^{(0)},\ldots,J_{\max}^{(19)}$ of the initial estimator were chosen as 24.

4 Preliminary results

The proof of our main theorems are based on several preliminary results. Let

[TABLE]

These quantities are related to the log-likelihoods of the smoothing distribution and its Gaussian approximation as

[TABLE]

where $\succ$ denotes the positive definite order (i.e. $\bm{A}\succ\bm{B}$ if and only if $\bm{A}-\bm{B}$ is positive definite). Similarly, $\succeq$ denotes the positive semidefinite order.

The next three propositions show various bounds on the log-likelihood related quantity $l^{\mathrm{sm}}(\bm{v})$. Their proof is included in Section A.1 of the Appendix.

Proposition 4.1** (A lower bound on the tails of $l^{\mathrm{sm}}(\bm{v})$).**

Suppose that Assumption 2.1 holds for $\bm{u}$, then for any $0<\varepsilon\leq 1$, for every $\sigma_{Z}>0$, $h\leq h_{\max}(\bm{u},T)$, we have

[TABLE]

where

[TABLE]

Proposition 4.2** (A bound on the difference between $l^{\mathrm{sm}}(\bm{v})$ and $l^{\mathrm{sm}}_{\mathcal{G}}(\bm{v})$).**

Suppose that Assumption 2.1 holds for $\bm{u}$, then for any $0<\varepsilon\leq 1$, $\sigma_{Z}>0$, $0<h\leq h_{\max}(\bm{u},T)$, we have

[TABLE]

where

[TABLE]

Proposition 4.3** (A bound on the difference between $\nabla l^{\mathrm{sm}}(\bm{v})$ and $\nabla l^{\mathrm{sm}}_{\mathcal{G}}(\bm{v})$).**

Suppose that Assumption 2.1 holds for $\bm{u}$, then for any $0<\varepsilon\leq 1$, $\sigma_{Z}>0$, $0<h\leq h_{\max}(\bm{u},T)$, we have

[TABLE]

where

[TABLE]

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 $\mathbb{R}^{d}$ with respect to the Lebesgue measure. Then their total variation distance satisfies that for any constant $c>0$,

[TABLE]

Proof.

We are going to use the following characterisation of the total variation distance,

[TABLE]

Based on this, the result is trivial for $c=1$. If $c>1$, then we have

[TABLE]

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 $\mathbb{R}^{d}$ with respect to the Lebesgue measure, also denoted by $f(\bm{x})$ and $g(\bm{x})$. Then their Wasserstein distance (defined as in (2.12)) satisfies that for any $\bm{y}\in\mathbb{R}^{d}$,

[TABLE]

Proof.

Let $m(\bm{x}):=\min(f(\bm{x}),g(\bm{x}))$, $\gamma:=\int_{\bm{x}\in\mathbb{R}^{d}}m(\bm{x})d\bm{x}$, $\hat{f}(\bm{x}):=f(\bm{x})-m(\bm{x})$, $\hat{g}(\bm{x}):=g(\bm{x})-m(\bm{x})$. Suppose first that $\gamma\neq 0$ and $1-\gamma\neq 0$. Let $\nu^{(1)}_{f,g}$ denote the distribution of a random vector $(\bm{X},\bm{X})$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ for a random variable $\bm{X}$ with distribution with density $m(\bm{x})/\gamma$ (the two components are equal). Let $\nu^{(2)}_{f,g}$ be a distribution on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with density $\nu^{(2)}_{f,g}(\bm{x}_{1},\bm{x}_{2}):=\frac{\hat{f}(\bm{x}_{1})}{1-\gamma}\cdot\frac{\hat{g}(\bm{x}_{2})}{1-\gamma}$ (the two components are independent).

We define the optimal coupling of $f$ and $g$ as a probability distribution $\nu_{f,g}$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ as a mixture of $\nu^{(1)}_{f,g}$ and $\nu^{(2)}_{f,g}$, that is, for any Borel-measurable $E\in\mathbb{R}^{d}\times\mathbb{R}^{d}$,

[TABLE]

It is easy to check that this distribution has marginals $f$ and $g$, so by (2.12), and the fact that $\int_{\bm{x}_{1}\in\mathbb{R}^{d}}\hat{f}(\bm{x}_{1})d\bm{x}_{1}=\int_{\bm{x}_{2}\in\mathbb{R}^{d}}\hat{g}(\bm{x}_{2})d\bm{x}_{2}=1-\gamma$, we have

[TABLE]

and the result follows from the fact that $\hat{f}(\bm{x})+\hat{g}(\bm{x})=|f(\bm{x})-g(\bm{x})|$. Finally, if $\gamma=0$ then we can set $\nu_{f,g}$ as $\nu^{(2)}_{f,g}$, and the same argument works, while if $\gamma=1$, then both sides of the claim are zero (since $f(\bm{x})=g(\bm{x})$ Lebesgue-almost surely). ∎

The following two lemmas show concentration bounds for the norm of $\bm{B}_{k}$ and the smallest eigenvalue of $\bm{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 $\bm{u}$, then the random vector $\bm{B}_{k}$ defined in (2.7) satisfies that for any $t\geq 0$,

[TABLE]

Thus for any $0<\varepsilon\leq 1$, we have

[TABLE]

Lemma 4.4**.**

Suppose that Assumption 2.1 holds for $\bm{u}$, then the random matrix $\bm{A}_{k}$ defined in (2.6) satisfies that for any $t\geq 0$,

[TABLE]

thus for any $0<\varepsilon\leq 1$, we have

[TABLE]

The following proposition bounds the derivatives of the log-determinant of the Jacobian.

Proposition 4.4**.**

The function $\log\det\bm{J}\Psi_{T}:\mathcal{B}_{R}\to\mathbb{R}$ is continuously differentiable on $\mathrm{int}(\mathcal{B}_{R})$ and its derivative can be bounded as

[TABLE]

Proof.

Let $\mathbb{M}^{d}$ denote the space of $d\times 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 (2009). Following their notation, we define the $k$th derivative of $\det$ at a point $\bm{M}\in\mathbb{M}^{d}$ as a map $D^{k}\det\bm{M}$ from $(\mathbb{M}^{d})^{k}$ to $\mathbb{C}$ with value

[TABLE]

They have defined the norm of $k$th derivative of the determinant as

[TABLE]

From Theorem 4 of Bhatia and Jain (2009), it follows that for any $k\geq 1$, the norm of the $k$th derivative can be bounded as

[TABLE]

Based on the chain rule, the norm first derivative of the log-determinant can be bounded as

[TABLE]

and the result follows by (1.21) and (4.10) for $k=1$ and $\bm{M}=\bm{J}\Psi_{T}(\bm{v})$. ∎

5 Proof of the main results

5.1 Gaussian approximation

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

5.1.1 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)}_{\mathrm{TV}}(\bm{u},T)$ and $C^{(2)}_{\mathrm{TV}}(\bm{u},T)$ sufficiently large, we can assume that $C_{\mathrm{TV}}(\bm{u},T,\varepsilon)$ satisfies the following bounds,

[TABLE]

Based on the assumption that $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, (5.1), and Lemma 4.4, we have

[TABLE]

where $C_{\|\bm{A}\|}$ is defined as in (2.14). The event $\lambda_{\min}(\bm{A}_{k})>\frac{c(\bm{u},T)}{2h}$ implies in particular that $\bm{A}_{k}$ is positive definite. From Lemma 4.3, we know that

[TABLE]

From Proposition 4.1, we know that

[TABLE]

Finally, from Proposition 4.2, it follows that

[TABLE]

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-\varepsilon$.

Let $\bm{Z}$ be a $d$ dimensional standard normal random vector, then from the definition of $\mu^{\mathrm{sm}}_{\mathcal{G}}$, it follows that when conditioned on $\bm{B}_{k}$ and $\bm{A}_{k}$,

[TABLE]

has distribution $\mu^{\mathrm{sm}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$. This fact will be used in the proof several times.

Since the normalising constant of the smoother, $C_{k}^{\mathrm{sm}}$ 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

[TABLE]

which is of similar form as the Gaussian approximation

[TABLE]

Let

[TABLE]

then based on the assumption that $\sigma_{Z}\sqrt{h}\leq C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, (5.2) and (5.3), it follows that

[TABLE]

Let $B_{\rho}:=\{\bm{v}\in\mathbb{R}^{d}:\|\bm{v}-\bm{u}\|\leq\rho(h,\sigma_{Z})\}$, and denote by $B_{\rho}^{c}$ its complement in $\mathbb{R}^{d}$. Then by Lemma 4.1, we have

[TABLE]

By Lemma 4.2, we can bound the Wasserstein distance of $\mu^{\mathrm{sm}}(\cdot|\bm{Y}_{0:k})$ and $\mu^{\mathrm{sm}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ as

[TABLE]

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

[TABLE]

Proof.

Note that by (5.14), we know that $B_{\rho}\subset\mathcal{B}_{R}$, and

[TABLE]

Now using (5.14), we can see that $\sup_{\bm{v}\in B_{\rho}}\left|\log\left(\frac{q(\bm{v})}{q(\bm{u})}\right)\right|\leq C_{q}^{(1)}\cdot\frac{1}{2C_{q}^{(1)}}\leq\frac{1}{2}$, and using (5.9), we have $\left|\frac{(l^{\mathrm{sm}}(\bm{v})-l^{\mathrm{sm}}_{\mathcal{G}}(\bm{v}))}{2\sigma_{Z}^{2}}\right|\leq\frac{1}{2}$. Using the fact that $|1-\exp(x)|\leq 2|x|$ for $-1\leq x\leq 1$, and the bounds $|\log(q(\bm{v})/q(\bm{u}))|\leq C_{q}^{(1)}\|\bm{v}-\bm{u}\|$ and (5.9), we can see that for every $\bm{v}\in B_{\rho}$,

[TABLE]

therefore

[TABLE]

Let $\bm{Z}$ denote a $d$-dimensional standard normal random vector, then it is easy to see that $\mathbb{E}(\|\bm{Z}\|)\leq(\mathbb{E}(\|\bm{Z}\|^{2}))^{1/2}\leq d^{1/2}$, and $\mathbb{E}(\|\bm{Z}\|^{3})\leq(\mathbb{E}(\|\bm{Z}\|^{4}))^{3/4}\leq(3d^{2})^{3/4}\leq 3d^{3/2}$. Since we have assumed that the events in (5.6) and (5.7) hold, we know that

[TABLE]

Finally, it is not difficult to show that for any $a,b\geq 0$, $(a+b)^{3}\leq 4(a^{3}+b^{3})$. Therefore

[TABLE]

thus the result follows. ∎

Lemma 5.2**.**

Using the notations and assumptions of this section, we have

[TABLE]

Proof.

Note that if $\bm{Z}$ is a $d$ dimensional standard normal random vector, then by Theorem 4.1.1 of Tropp (2015), we have

[TABLE]

Since the random variable $\bm{W}$ defined in (5.10) is distributed as $\mu^{\mathrm{sm}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ when conditioned on $\bm{A}_{k},\bm{B}_{k}$, we have

[TABLE]

Based on the assumption that $\sigma_{Z}\sqrt{h}\leq C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, and (5.4), we have $\frac{2C_{\bm{B}}\left(\frac{\varepsilon}{4}\right)}{c(\bm{u},T)}\cdot\sigma_{Z}\sqrt{h}\leq\frac{\rho(h,\sigma_{Z})}{2}$, and

[TABLE]

and the result follows by (5.20). ∎

Lemma 5.3**.**

Using the notations and assumptions of this section, we have

[TABLE]

Proof.

Let $q_{\max}:=\sup_{\bm{v}\in\mathcal{B}_{R}}q(\bm{v})$. By our assumption that the event in (5.8) holds, we have

[TABLE]

where in the last step we have used the fact that $\det(\bm{A}_{k})^{1/2}\leq\|\bm{A}_{k}\|^{d/2}$. Based on the assumption that $\sigma_{Z}\sqrt{h}\leq C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, and (5.5), we have

[TABLE]

therefore by (5.20), we have

[TABLE]

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_{\mathrm{TV}}^{(1)}(\bm{u},T)$ and $C_{\mathrm{TV}}^{(2)}(\bm{u},T)$ sufficiently large such that we have

[TABLE]

Now we bound the three terms needed for the Wasserstein distance.

Lemma 5.4**.**

Using the notations and assumptions of this section, we have

[TABLE]

for some finite positive constants $C_{1}^{*}(\bm{u},T)$, $C_{2}^{*}(\bm{u},T)$.

Proof.

Note that for any $\bm{v}\in\mathcal{B}_{\rho}$, we have

[TABLE]

From (5.21), and the fact that $\tilde{\mu}^{\mathrm{sm}}(\cdot|\bm{Y}_{0:k})$ is a rescaled version of $\tilde{\mu}^{\mathrm{sm}}(\cdot|\bm{Y}_{0:k})$, it follows that for $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, we have

[TABLE]

By (5.18) and (5.14), it follows that $\frac{\tilde{\mu}^{\mathrm{sm}}(\bm{v}|\bm{Y}_{0:k})}{\mu^{\mathrm{sm}}_{\mathcal{G}}(\bm{v}|\bm{Y}_{0:k})}\leq 3$ for every $\bm{v}\in\mathcal{B}_{\rho}$. By using these and bounding $\left|\frac{\tilde{\mu}^{\mathrm{sm}}(\bm{v}|\bm{Y}_{0:k})}{\mu^{\mathrm{sm}}_{\mathcal{G}}(\bm{v}|\bm{Y}_{0:k})}-1\right|$ via (5.18), we obtain that

[TABLE]

Based on (5.19), we have

[TABLE]

where $\bm{Z}$ is a $d$ dimensional standard normal random vector. The claimed result now follows using the fact that $\mathbb{E}(\|\bm{Z}\|)\leq\sqrt{d}$, $\mathbb{E}(\|\bm{Z}\|^{2})\leq d$, and $\mathbb{E}(\|\bm{Z}\|^{4})\leq 3d^{2}$. ∎

Lemma 5.5**.**

Using the notations and assumptions of this section, we have

[TABLE]

Proof.

Using (5.22), Lemma 5.3, and the fact that $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, we have

[TABLE]

and the result follows from $\int_{\bm{v}\in B_{\rho}^{c}}\|\bm{v}-\bm{u}\|\mu^{\mathrm{sm}}(\bm{v}|\bm{Y}_{0:k})d\bm{v}\leq 2R\cdot\mu^{\mathrm{sm}}(B_{\rho}^{c}|\bm{Y}_{0:k})$. ∎

Lemma 5.6**.**

Using the notations and assumptions of this section, we have

[TABLE]

for some finite positive constants $C_{3}^{*}(\bm{u},T)$, $C_{4}^{*}(\bm{u},T)$.

Proof.

Let $\bm{W}$ be defined as in (5.10), and $\bm{Z}$ be $d$-dimensional standard normal. Using the fact that for a non-negative valued random variable $X$, we have $\mathbb{E}(X)=\int_{t=0}^{\infty}\mathbb{P}(X\geq t)dt$, it follows that

[TABLE]

[TABLE]

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_{\mathrm{W}}^{(1)}(\bm{u},T)$ and $C_{\mathrm{W}}^{(2)}(\bm{u},T)$, we have

[TABLE]

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-\varepsilon$. ∎

5.1.2 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 $\mu^{\mathrm{sm}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ to $\mathcal{B}_{R}$, denoted by $\mu^{\mathrm{sm}}_{\mathcal{G}|\mathcal{B}_{R}}(\cdot|\bm{Y}_{0:k})$ as

[TABLE]

This is a probability distribution which is supported on $\mathcal{B}_{R}$. We denote its push-forward map by $\Psi_{T}$ as $\eta^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$, i.e. if a random vector $\bm{X}$ is distributed as $\mu^{\mathrm{sm}}_{\mathcal{G}|\mathcal{B}_{R}}(\cdot|\bm{Y}_{0:k})$, then $\eta^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ denotes the distribution of $\Psi_{T}(\bm{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,

[TABLE]

Proof.

First, notice that by Proposition 3(f) of Roberts and Rosenthal (2004), we have

[TABLE]

By Proposition 3(g) of Roberts and Rosenthal (2004), there is a coupling $(\bm{X}_{1},\bm{X}_{2})$ of random vectors such that $\bm{X}_{1}\sim\mu^{\mathrm{sm}}(\cdot|\bm{Y}_{0:k})$, $\bm{X}_{2}\sim\mu^{\mathrm{sm}}_{\mathcal{G}|\mathcal{B}_{R}}(\cdot|\bm{Y}_{0:k})$, and $\mathbb{P}(\bm{X}_{1}\neq\bm{X}_{2}|\bm{Y}_{0:k})=d_{\mathrm{TV}}(\mu^{\mathrm{sm}}(\cdot|\bm{Y}_{0:k}),\mu^{\mathrm{sm}}_{\mathcal{G}|\mathcal{B}_{R}}(\cdot|\bm{Y}_{0:k}))$. Given this coupling, we look at the coupling of the transformed random variables $(\Psi_{T}(\bm{X}_{1}),\Psi_{T}(\bm{X}_{2}))$. This obviously satisfies that $\mathbb{P}(\Psi_{T}(\bm{X}_{1})\neq\Psi_{T}(\bm{X}_{2})|\bm{Y}_{0:k})\leq d_{\mathrm{TV}}(\mu^{\mathrm{sm}}(\cdot|\bm{Y}_{0:k}),\mu^{\mathrm{sm}}_{\mathcal{G}|\mathcal{B}_{R}}(\cdot|\bm{Y}_{0:k}))$. Moreover, we have $\Psi_{T}(\bm{X}_{1})\sim\mu^{\mathrm{fi}}(\cdot|\bm{Y}_{0:k})$ and $\Psi_{T}(\bm{X}_{2})\sim\eta^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$, thus

[TABLE]

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 can be bounded as follows,

[TABLE]

Proof.

By Theorem 4.1 of Villani (2009), there exists a coupling of random variables $(\bm{X}_{1},\bm{X}_{2})$ (called the optimal coupling) such that $\bm{X}_{1}\sim\mu^{\mathrm{sm}}(\cdot|\bm{Y}_{0:k})$, $\bm{X}_{2}\sim\mu^{\mathrm{sm}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$, and

[TABLE]

Let $\hat{\bm{X}}_{2}:=\bm{X}_{2}\cdot 1_{[\|\bm{X}_{2}\|\leq R]}+\frac{\bm{X}_{2}}{\|\bm{X}_{2}\|}\cdot R\cdot 1_{[\|\bm{X}_{2}\|>R]}$ denote the projection of $\bm{X}_{2}$ on the ball $\mathcal{B}_{R}$. Then using the fact that $\bm{X}_{1}\in\mathcal{B}_{R}$, it is easy to see that $\|\hat{\bm{X}}_{2}-\bm{X}_{1}\|\leq\|\bm{X}_{2}-\bm{X}_{1}\|$, and therefore

[TABLE]

where $\mathcal{L}(\hat{\bm{X}}_{2}|\bm{Y}_{0:k})$ denotes the distribution of $\hat{\bm{X}}_{2}$ conditioned on $\bm{Y}_{0:k}$.

Moreover, by the definitions, for a given $\bm{Y}_{0:k}$, it is easy to see we can couple random variables $\hat{\bm{X}}_{2}$ and $\tilde{\bm{X}}_{2}\sim\mu^{\mathrm{sm}}_{\mathcal{G}|\mathcal{B}_{R}}(\cdot|\bm{Y}_{0:k})$ such that they are the same with probability at least $1-\mu^{\mathrm{sm}}_{\mathcal{G}}(\mathcal{B}_{R}^{c})$. Since the maximum distance between two points in $\mathcal{B}_{R}$ is at most $2R$, it follows that

[TABLE]

By the triangle inequality, we obtain that

[TABLE]

and by (1.6), it follows that

[TABLE]

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 $\eta^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ and $\mu^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$. Let

[TABLE]

then the density of $\mu^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ can be written as

[TABLE]

Since the normalising constant is not known for the case of $\eta^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$, we define a rescaled version $\tilde{\eta}^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{Y}_{0:k})$ with density

[TABLE]

The following lemma bounds the difference between the logarithms of $\tilde{\eta}^{\mathrm{fi}}_{\mathcal{G}}(\bm{v}|\bm{Y}_{0:k})$ and $\mu^{\mathrm{fi}}_{\mathcal{G}}(\bm{v}|\bm{Y}_{0:k})$.

Lemma 5.9**.**

For any $\bm{v}\in\Psi_{T}(\mathcal{B}_{R})$, we have

[TABLE]

Proof of Lemma 5.9.

By (5.27) and (5.28), we have

[TABLE]

The absolute value of the first difference can be bounded by Propositon 4.4 as

[TABLE]

For any two vectors $\bm{x},\bm{y}\in\mathbb{R}^{d}$, we have

[TABLE]

so the second difference can be bounded as

[TABLE]

Using (1.6), we have

[TABLE]

By (1.19), we have

[TABLE]

so by (1.6), it follows that

[TABLE]

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)}_{\mathrm{TV}}(\bm{u},T)\geq C^{(1)}_{\mathrm{TV}}(\bm{u},T)$ and $D^{(2)}_{\mathrm{TV}}(\bm{u},T)\geq C^{(2)}_{\mathrm{TV}}(\bm{u},T)$, thus $D_{\mathrm{TV}}(\bm{u},T,\varepsilon)\geq C_{\mathrm{TV}}(\bm{u},T,\varepsilon)$. We also assume that $D^{(1)}_{\mathrm{TV}}(\bm{u},T)$ satisfies that

[TABLE]

Based on these assumptions on $D_{\mathrm{TV}}(\bm{u},T,\varepsilon)$, and the assumption that $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}D_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-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-\varepsilon$. We are going to assume that this is the case for the rest of the proof. We define

[TABLE]

Based on (5.29), (5.30), and the assumption that $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}D_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, it follows that

[TABLE]

The surface of a ball of radius $R-\|\bm{u}\|$ centered at $\bm{u}$ is contained in $\mathcal{B}_{R}$, and it will be transformed by $\Psi_{T}$ to a closed continuous manifold whose points are at least $(R-\|\bm{u}\|)\exp(-GT)$ away from $\Psi_{T}(\bm{u})$ (based on (1.6)). This implies that the ball of radius $(R-\|\bm{u}\|)\exp(-GT)$ centered at $\Psi_{T}(\bm{u})$ is contained in $\Psi_{T}(\mathcal{B}_{R})$, and thus by (5.32), $\mathcal{B}_{\rho^{\prime}}\subset\Psi_{T}(\mathcal{B}_{R})$.

By Lemma 4.1, we have

[TABLE]

By Lemma 5.9, (5.32), and the fact that $|\exp(x)-1|\leq 2|x|$ for $x\in[-1,1]$, it follows that for $\bm{v}\in\mathcal{B}_{\rho^{\prime}}$, we have

[TABLE]

Therefore the first term of (5.33) can be bounded as

[TABLE]

This in turn can be bounded as in Lemma 5.1. The terms $\tilde{\eta}^{\mathrm{fi}}_{\mathcal{G}}(\mathcal{B}_{\rho^{\prime}}^{c}|\bm{Y}_{0:k})$ and $\mu^{\mathrm{fi}}_{\mathcal{G}}(\mathcal{B}_{\rho^{\prime}}^{c}|\bm{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)}_{\mathrm{TV}}(\bm{u},T)$ and $D^{(2)}_{\mathrm{TV}}(\bm{u},T)$ such that

[TABLE]

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

[TABLE]

Therefore we only need to bound the last two terms of (5.25). The fact that $\mu^{\mathrm{sm}}_{\mathcal{G}}(\mathcal{B}_{R}^{c}|\bm{Y}_{0:k}))=o(\sigma_{Z}^{2}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 $\bm{y}:=\Psi_{T}(\bm{u}^{\mathcal{G}})$. This implies that

[TABLE]

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. ∎

5.2 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.

5.2.1 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 $\bm{u}^{\mathcal{G}}:=\bm{u}-\bm{A}_{k}^{-1}\bm{B}_{k}$ denote the center of the Gaussian approximation $\mu^{\mathrm{sm}}_{\mathcal{G}}$ (defined when $\bm{A}_{k}$ is positive definite). Let $\rho(h,\sigma_{Z})$ be as in (5.13), $B_{\rho}:=\{\bm{v}\in\mathbb{R}^{d}:\|\bm{v}-\bm{u}\|\leq\rho(h,\sigma_{Z})\}$, and $B_{\rho}^{c}$ be the complement of $B_{\rho}$. 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 $\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}^{\mathcal{G}}\|$).**

There are some finite constants $D_{5}(\bm{u},T)$ and $D_{6}(\bm{u},T)$ such that for any $0<\varepsilon\leq 1$, for $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}\cdot C_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, we have

[TABLE]

Proof.

This is a direct consequence of the Wasserstein distance bound of Theorem 2.1, since

[TABLE]

Lemma 5.11** (A bound on $\|\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}-\bm{u}\|$).**

For any $0<\varepsilon\leq 1$, we have

[TABLE]

Proof.

From Proposition 4.1, and (4.3), it follows that for any $0<\varepsilon\leq 1$,

[TABLE]

Since $\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}$ is the maximizer of $\log\mu^{\mathrm{sm}}(\bm{v}|\bm{Y}_{0:k})$ on $\mathcal{B}_{R}$, our claim follows. ∎

Lemma 5.12** (A bound on $\|\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}-\bm{u}^{\mathcal{G}}\|$).**

There are finite constants $S_{\mathrm{MAP}}^{(1)}>0$, $S_{\mathrm{MAP}}^{(2)}$, $D_{7}(\bm{u},T)$ and $D_{8}(\bm{u},T)$ such that for any $0<\varepsilon\leq 1$, for $\sigma_{Z}\sqrt{h}<\left(S_{\mathrm{MAP}}^{(1)}+S_{\mathrm{MAP}}^{(2)}\left(\log\left(\frac{1}{\varepsilon}\right)\right)^{1/2}\right)^{-1}$, we have

[TABLE]

Proof.

By choosing $S_{\mathrm{MAP}}^{(1)}$ and $S_{\mathrm{MAP}}^{(2)}$ sufficiently large, we can assume that

[TABLE]

From Lemma 5.11, we know that

[TABLE]

Using (5.37), it follows that if the above event happens, then $\|\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}\|<R$, and thus

[TABLE]

Using the fact that $\nabla l^{\mathrm{sm}}_{\mathcal{G}}(\bm{v})=2\bm{A}_{k}(\bm{v}-\bm{u}^{\mathcal{G}})$, and Proposition 4.3, it follows that

[TABLE]

Moreover, by Lemma 4.4, we know that for any $0<\varepsilon\leq 1$, we have

[TABLE]

By combining the four equations (5.38), (5.39), (5.40) and (5.41), it follows that with probability at least $1-\varepsilon$, we have

[TABLE]

and the claim of the lemma follows by rearrangement. ∎

Lemma 5.13** (A lower bound on $\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}\|^{2}\right|\bm{u}\right)$).**

There are positive constants $D_{9}(\bm{u},T)$ and $D_{10}(\bm{u},T)$ such that for $\sigma_{Z}\sqrt{h}\leq D_{10}(\bm{u},T)$, we have

[TABLE]

Proof.

By applying Lemma 5.10 for $\varepsilon=0.1$, we obtain that for $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}\cdot C_{\mathrm{TV}}(\bm{u},T,0.1)^{-1}$, we have

[TABLE]

If this event happens, then in particular, we have

[TABLE]

By the definition of $\bm{B}_{k}$ in (2.7), it follows that conditioned on $\bm{u}$, $\bm{B}_{k}$ has $d$-dimensional multivariate normal distribution with covariance matrix $\Sigma_{\bm{B}_{k}}:=\sigma_{Z}^{2}\sum_{i=0}^{k}\bm{J}\Phi_{t_{i}}(\bm{u})^{\prime}\bm{J}\Phi_{t_{i}}(\bm{u})$. This means that if $\bm{Z}$ is a $d$ dimensional standard normal random vector, then $(\Sigma_{\bm{B}_{k}})^{1/2}\cdot\bm{Z}$ has the same distribution as $\bm{B}_{k}$ (conditioned on $\bm{u}$). By Assumption 2.1, we have $\lambda_{\min}\left(\Sigma_{\bm{B}_{k}}\right)\geq\sigma_{Z}^{2}\cdot\frac{c(\bm{u},T)}{h}$, and thus $\|(\Sigma_{\bm{B}_{k}})^{1/2}\cdot\bm{Z}\|\geq\sigma_{Z}\cdot\sqrt{\frac{c(\bm{u},T)}{h}}\cdot\|\bm{Z}\|$.

It is not difficult to show that for any $d\geq 1$, $\mathbb{P}\left(\|\bm{Z}\|\geq\frac{\sqrt{d}}{2}\right)\geq\frac{1}{4}$ (indeed, if $(Z_{(i)})_{1\leq i\leq d}$ are i.i.d. standard normal random variables, then for any $\lambda>0$, $\mathbb{E}(e^{-\lambda Z_{1}^{2}})=\sqrt{\frac{1}{1+2\lambda}}$, so $\mathbb{E}(e^{-\lambda(\|\bm{Z}\|^{2}-d)})=\mathbb{E}(e^{-\lambda\sum_{i=1}^{d}(Z_{(i)}^{2}-1)})=(1+2\lambda)^{-d/2}\cdot e^{\lambda d}$, and the claim follows by applying Markov’s inequality $\mathbb{P}(\|\bm{Z}\|^{2}-d\leq-t)\leq\mathbb{E}(e^{-\lambda(\|\bm{Z}\|^{2}-d)})\cdot e^{-\lambda t}$ for $t=\frac{3}{4}d$ and $\lambda=1$). Therefore, we have

[TABLE]

and thus by (5.43) and (5.44), it follows that

[TABLE]

By choosing $D_{10}(\bm{u},T)$ sufficiently small, we have that for $\sigma_{Z}\sqrt{h}\leq D_{10}(\bm{u},T)$,

[TABLE]

and the result follows. ∎

Lemma 5.14** (A bound on the difference of $\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}\|^{2}\right|\bm{u}\right)$ and $\mathbb{E}\left(\left.\|\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}-\bm{u}\|^{2}\right|\bm{u}\right)$).**

There are some finite constants $D_{11}(\bm{u},T)$ and $D_{12}(\bm{u},T)>0$ such that for $\sigma_{Z}\sqrt{h}\leq D_{12}(\bm{u},T)$, we have

[TABLE]

Proof.

We define the event $E_{k}$ as

[TABLE]

Under this event, we have, in particular $\bm{A}_{k}\succ\bm{0}$ and $\|\bm{u}^{\mathcal{G}}\|<R$. Let $E_{k}^{c}$ denote the complement of $E_{k}$. Then the difference in the variances can be bounded as

[TABLE]

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

[TABLE]

for some finite constant $C_{E}(\bm{u},T)$ independent of $h$ and $\sigma_{Z}$.

The term $\mathbb{E}\left(\left.1_{E_{k}}\|\bm{u}^{\mathcal{G}}-\bm{u}\|^{2}\right|\bm{u}\right)$ can be bounded as

[TABLE]

For bounding the term $\mathbb{E}\left(\left.1_{E_{k}}\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}^{\mathcal{G}}\|^{2}\right|\bm{u}\right)$, we define

[TABLE]

then by Lemma 5.10, it follows that for $\sigma_{Z}\sqrt{h}<\frac{1}{2}(C_{\mathrm{TV}}^{(1)})^{-1}$, for $t\in[t_{\min},t_{\max}]$, we have

[TABLE]

By writing

[TABLE]

and using the fact that $1_{E_{k}}\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}^{\mathcal{G}}\|^{2}<4R^{2}$, one can show that for $\sigma_{Z}\sqrt{h}<\overline{S}(\bm{u},T)$,

[TABLE]

for some constants $\overline{S}(\bm{u},T)>0$, $\overline{C}(\bm{u},T)<\infty$, that are independent of $\sigma_{Z}$ and $h$.

Finally, for bounding the term $\mathbb{E}\left(\left.1_{E_{k}}\|\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}-\bm{u}^{\mathcal{G}}\|^{2}\right|\bm{u}\right)$, we define

[TABLE]

By Lemma 5.12, it follows that for $\sigma_{Z}\sqrt{h}<\left(S_{\mathrm{MAP}}^{(1)}\right)^{-1}$, for $t\in[t_{\min}^{\prime},t_{\max}^{\prime}]$, we have

[TABLE]

which implies that for $\sigma_{Z}\sqrt{h}<S_{\mathrm{M}}(\bm{u},T)$,

[TABLE]

for some constants $S_{\mathrm{M}}(\bm{u},T)>0$, $C_{\mathrm{M}}(\bm{u},T)>0$.

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 $\frac{\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}\|^{2}\right|\bm{u}\right)}{\sigma_{Z}^{2}h}$ follows by Lemma 5.13. Let $E_{k}$ be defined as in (5.46), then we have

[TABLE]

so the upper bound on $\frac{\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}\|^{2}\right|\bm{u}\right)}{\sigma_{Z}^{2}h}$ follows by (5.48), (5.49) and (5.51). Finally, the bound on $\left|\mathbb{E}\left[\|\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}-\bm{u}\|^{2}|\bm{u}\right]-\mathbb{E}\left[\|\overline{\bm{u}}^{\mathrm{sm}}-\bm{u}\|^{2}|\bm{u}\right]\right|$ follows directly from Lemma 5.14. ∎

5.2.2 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 $\|\overline{\bm{u}}^{\mathrm{fi}}-\Psi_{T}(\bm{u}^{\mathcal{G}})\|$).**

There are some finite constants $D_{5}^{\prime}(\bm{u},T)$ and $D_{6}^{\prime}(\bm{u},T)$ such that for any $0<\varepsilon\leq 1$, for $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}\cdot D_{\mathrm{TV}}(\bm{u},T,\varepsilon)^{-1}$, we have

[TABLE]

Proof.

This is a direct consequence of the Wasserstein distance bound of Theorem 2.2. ∎

Lemma 5.16** (A bound on $\|\hat{\bm{u}}^{\mathrm{fi}}-\Psi_{T}(\bm{u}^{\mathcal{G}})\|$).**

There are finite constants $S_{\mathrm{MAP}}^{(1^{\prime})}>0$, $S_{\mathrm{MAP}}^{(2^{\prime})}$, $D_{7}^{\prime}(\bm{u},T)$ and $D_{8}^{\prime}(\bm{u},T)$ such that for any $0<\varepsilon\leq 1$, for $\sigma_{Z}\sqrt{h}<\left(S_{\mathrm{MAP}}^{(1^{\prime})}+S_{\mathrm{MAP}}^{(2^{\prime})}\left(\log\left(\frac{1}{\varepsilon}\right)\right)^{1/2}\right)^{-1}$, we have

[TABLE]

Proof.

The result follows by Lemma 5.12, Theorem 2.2, and (1.6). ∎

Lemma 5.17** (A lower bound on $\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{fi}}-\bm{u}(T)\|^{2}\right|\bm{u}\right)$).**

There are positive constants $D_{9}^{\prime}(\bm{u},T)$ and $D_{10}^{\prime}(\bm{u},T)$ such that for $\sigma_{Z}\sqrt{h}\leq D_{10}^{\prime}(\bm{u},T)$, we have

[TABLE]

Proof.

The proof is similar to the proof of Lemma 5.13. By applying Lemma 5.15 for $\varepsilon=0.1$, we obtain that for $\sigma_{Z}\sqrt{h}\leq\frac{1}{2}\cdot D_{\mathrm{TV}}(\bm{u},T,0.1)^{-1}$, we have

[TABLE]

If this event happens, then by (1.6), we have

[TABLE]

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 $\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{fi}}-\bm{u}\|^{2}\right|\bm{u}\right)$ and $\mathbb{E}\left(\left.\|\hat{\bm{u}}^{\mathrm{fi}}-\bm{u}\|^{2}\right|\bm{u}\right)$).**

There are some finite constants $D_{11}^{\prime}(\bm{u},T)$ and $D_{12}^{\prime}(\bm{u},T)>0$ such that for $\sigma_{Z}\sqrt{h}\leq D_{12}^{\prime}(\bm{u},T)$, we have

[TABLE]

Proof.

We define the event $E_{k}$ as in (5.46). Under this event, $\bm{A}_{k}\succ\bm{0}$ and $\|\bm{u}^{\mathcal{G}}\|<R$, so $\mu^{\mathrm{fi}}_{\mathcal{G}}(\cdot|\bm{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

[TABLE]

Based on Lemma 5.15, we obtain that for $\sigma_{Z}\sqrt{h}<\overline{S}^{\prime}(\bm{u},T)$,

[TABLE]

for some constants $\overline{S}^{\prime}(\bm{u},T)>0$, $\overline{C}^{\prime}(\bm{u},T)<\infty$, that are independent of $\sigma_{Z}$ and $h$. We omit the details. ∎

Proof of Theorem 2.4.

The lower bound on $\frac{\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{fi}}-\bm{u}(T)\|^{2}\right|\bm{u}\right)}{\sigma_{Z}^{2}h}$ follows by Lemma 5.17. Similarly to the case of the smoother, we have

[TABLE]

which can be further bounded by (5.48), (5.58) and (5.59) to yield the upper bound on $\frac{\mathbb{E}\left(\left.\|\overline{\bm{u}}^{\mathrm{fi}}-\bm{u}(T)\|^{2}\right|\bm{u}\right)}{\sigma_{Z}^{2}h}$. Finally, the bound on $\left|\mathbb{E}\left[\|\hat{\bm{u}}^{\mathrm{fi}}-\bm{u}(T)\|^{2}|\bm{u}\right]-\mathbb{E}\left[\|\overline{\bm{u}}^{\mathrm{fi}}-\bm{u}(T)\|^{2}|\bm{u}\right]\right|$ follows directly from Lemma 5.18. ∎

5.3 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 $\bm{v}\in R^{d}$, $r>0$, we denote the ball of radius $r$ centered at $\bm{v}$ by $B(\bm{v},r):=\{\bm{x}\in\mathbb{R}^{d}:\|\bm{v}-\bm{x}\|\leq r^{*}\}$)

Proposition 5.1**.**

Suppose that $\Omega\subset\mathbb{R}^{d}$ is an open set, and $g:\Omega\to\mathbb{R}$ is a 3 times continuously differentiable function satisfying that

$g$* has a local minimum at a point $\bm{x}^{*}\in\Omega$,* 2. 2.

there exists a radius $r^{*}>0$ and constants $C_{H}>0,L_{H}<\infty$ such that $B(\bm{x}^{*},r^{*})\subset\Omega$, $\nabla^{2}g(\bm{x})\succeq C_{H}\cdot\bm{I}_{d}$ for every $\bm{x}\in B(\bm{x}^{*},r^{*})$, and $\nabla^{2}g(\bm{x})$ is $L_{H}$-Lipschitz on $B(\bm{x}^{*},r^{*})$.

Suppose that the starting point $x_{0}\in\Omega$ satisfies that $\|\bm{x}_{0}-\bm{x}^{*}\|<\min\left(r^{*},2\frac{C_{H}}{L_{H}}\right)$. Then the iterates of Newton’s method defined recursively for every $i\in\mathbb{N}$ as

[TABLE]

always stay in $B(\bm{x}^{*},r^{*})$ (thus they are well defined), and satisfy that

[TABLE]

Proof.

We will show that $\bm{x}_{i}\in B(\bm{x}^{*},r^{*})$ by induction in $i$. This is true for $i=0$. Assuming that $\bm{x}_{i}\in B(\bm{x}^{*},r^{*})$, we can write the gradient $\nabla g(\bm{x}_{i})$ as

[TABLE]

By the $L_{H}$-Lipschitz property of $\nabla^{2}g(\bm{x})$ on $B(\bm{x}^{*},r^{*})$, we have

[TABLE]

By combining this with the fact that $\|(\nabla^{2}g(\bm{x}_{i}))^{-1}\|\leq\frac{1}{C_{H}}$, we obtain that $\|\bm{x}_{i+1}-\bm{x}^{*}\|\leq\frac{L_{H}}{2C_{H}}\|\bm{x}_{i}-\bm{x}^{*}\|^{2}$ for every $i\in\mathbb{N}$, and by rearrangement, it follows that $\log\left(\frac{L_{H}}{2C_{H}}\|\bm{x}_{i+1}-\bm{x}^{*}\|\right)\leq 2\log\left(\frac{L_{H}}{2C_{H}}\|\bm{x}_{i}-\bm{x}^{*}\|\right)$ for every $i\in\mathbb{N}$, hence the result. ∎

The following proposition gives a lower bound on the Hessian near $\bm{u}$. The proof is included in Section A.1 of the Appendix.

Proposition 5.2**.**

Suppose that Assumptions 2.1 and 2.3 hold for the initial point $\bm{u}$ and the prior $q$. Let

[TABLE]

Then for any $0<\varepsilon\leq 1$, $\sigma_{Z}>0$, $0<h\leq h_{\max}(\bm{u},T)$, we have

[TABLE]

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 $\bm{u}$ and the prior $q$. Then for any $0<\varepsilon\leq 1$, $\sigma_{Z}>0$, $0<h\leq h_{\max}(\bm{u},T)$, we have

[TABLE]

[TABLE]

Now we are ready to prove Theorem 2.5.

Proof of Theorem 2.5.

Let

[TABLE]

Then by the assumption that $\sigma_{Z}\sqrt{h}\leq S_{\max}^{\mathrm{sm}}(\bm{u},T,\varepsilon)$, using Propositions 5.2, 5.3 and inequality (5.35) we know that with probability at least $1-\varepsilon$, all three of the following events hold at the same time,

$\nabla^{2}g^{\mathrm{sm}}(\bm{v})\succeq C_{g}\cdot\bm{I}_{d}$ for every $\bm{v}\in B(\bm{u},r_{H}(\bm{u},T))$ for $C_{g}:=\frac{1}{2}\frac{c(\bm{u},T)}{\sigma_{Z}^{2}h}$, 2. 2.

$\nabla^{2}g^{\mathrm{sm}}$ is $L$-Lipschitz in $\mathcal{B}_{R}$ for $L:=\frac{2C_{7}(\bm{u},T)}{\sigma_{Z}^{2}h}$, 3. 3.

$\|\bm{u}-\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}\|\leq\min\left(\frac{r_{H}(\bm{u},T)}{4},\frac{1}{4}\frac{c(\bm{u},T)}{C_{7}(\bm{u},T)}\right)$.

If these events hold, then the conditions of Proposition 5.1 are satisfied for the function $g$ with $\bm{x}^{*}=\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}$, $r^{*}:=\frac{3}{4}r_{H}(\bm{u},T)$, $C_{H}:=\frac{1}{2}\frac{c(\bm{u},T)}{\sigma_{Z}^{2}h}$ and $L_{H}:=\frac{2C_{7}(\bm{u},T)}{\sigma_{Z}^{2}h}$. Therefore, (5.60) holds if $\|\bm{x}_{0}-\hat{\bm{u}}^{\mathrm{sm}}_{\mathrm{MAP}}\|\leq\min\left(\frac{3}{4}r_{H}(\bm{u},T),\frac{1}{2}\frac{c(\bm{u},T)}{C_{7}(\bm{u},T)}\right)$. By the triangle inequality, and the third event above, this is satisfied if $\|\bm{x}_{0}-\bm{u}\|\leq\min\left(\frac{1}{2}r_{H}(\bm{u},T),\frac{1}{4}\frac{c(\bm{u},T)}{C_{7}(\bm{u},T)}\right)$. Thus the claim of the theorem follows for

[TABLE]

5.4 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 $\hat{\Phi}^{(l|j_{\max})}(\bm{Y}_{0:\hat{k}})$ with remainder term of order $j_{\max}+1$, we obtain that

[TABLE]

where by (1.14) the remainder term $R_{l,j_{\max}+1}$ can be bounded using the Cauchy-Schwarz inequality as

[TABLE]

Due to the particular choice of the coefficients of $\bm{c}^{(l|j_{\max}|\hat{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

[TABLE]

using the assumption that $\hat{k}\geq 2j_{\max}+3$. The concentration bound now follows directly from this bias bound, (5.20) and the fact that the estimator $\hat{\Phi}^{(l|j_{\max})}(\bm{Y}_{0:\hat{k}})$ has Gaussian distribution with covariance matrix $\left\|\bm{c}^{(l|j_{\max}|\hat{k})}\right\|\cdot\sigma_{Z}\cdot\bm{I}_{d_{o}}$. ∎

Proof of Lemma 2.1.

Let $\mathcal{P}$ denote the space of finite degree polynomials with real coefficients on $[0,1]$. For $a,b\in\mathcal{P}$, we let $\left<a,b\right>_{\mathcal{P}}=\int_{x=0}^{1}a(x)b(x)dx$. Then the elements of the matrix $\bm{K}^{j_{\max}}$ can be written as

[TABLE]

If $\bm{K}^{j_{\max}}$ would not be invertible, then its rows would be linearly dependent, that is, there would exist a non-zero vector $\bm{\alpha}\in\mathbb{R}^{j_{\max}+1}$ such that $\sum_{i=1}^{j_{\max}+1}K^{j_{\max}}_{i,j}=0$ for every $1\leq j\leq j_{\max}+1$. This would imply that $\left<\sum_{i=1}^{j_{\max}+1}\alpha_{i}x^{i-1},x^{j-1}\right>_{\mathcal{P}}=0$ for every $1\leq j\leq j_{\max}+1$, and thus

[TABLE]

However, this is not possible, since by the fundamental theorem of algebra, $\sum_{i=1}^{j_{\max}+1}\alpha_{i}x^{i-1}$ can have at most $j_{\max}$ roots, so it cannot be zero Lebesgue almost everywhere in $[0,1]$. Therefore $\bm{K}^{j_{\max}}$ is invertible. The result now follows from the continuity of the matrix inverse and the fact that for any $1\leq i,j\leq j_{\max}+1$,

[TABLE]

Proof of Theorem 2.6.

Let

[TABLE]

Based on Lemma 2.1, this is finite. With the choice $j_{\max}=l$, by (2.37), we obtain that

[TABLE]

By choosing $s_{\max}^{(l)}$ sufficiently small, we can ensure that for $\sigma_{Z}\sqrt{h}\leq s_{\max}^{(l)}$, $h\hat{k}_{\mathrm{min}}(l,l)\leq T$, and thus by the definition (2.38), we have

[TABLE]

By substituting this into (2.35), and applying some algebra, we obtain that

[TABLE]

and the claim of the theorem now follows by substituting this into Proposition 2.2. ∎

6 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 Dimet and Talagrand (1986); Talagrand and Courtier (1987)). The precise size of the observation noise $\sigma_{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 $\sigma_{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 $\sigma_{Z}$ and $h$. These flow-dependent prior distributions can considerably improve filtering accuracy. Going beyond the approximately Gaussian case (for example when $\sigma_{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 et al. (2010) for some examples where this situation arises).

Acknowledgements

DP & AJ were supported by an AcRF tier 2 grant: R-155-000-161-112. AJ is affiliated with the Risk Management Institute, the Center for Quantitative Finance and the OR & Analytics cluster at NUS. DC was partially supported by the EPSRC grant: EP/N023781/1. We thank the anonymous referees for their valuable comments and suggestions that have substantially improved the paper.

Appendix A Appendix

A.1 Proof of preliminary results

Proof of Lemma 1.1.

To prove (1.14), it suffices to first verify (1.14) and for $i=0$ and $i=1$, and then use induction and the recursion formula (1.13) for $i\geq 2$. For (1.15), by taking the $k$th derivative of (1.13), we obtain that

[TABLE]

For $k=1$, (1.15) can be verified by first checking it for $i=0$ and $i=1$, and then using mathematical induction and (1.13) for $i\geq 2$. Suppose that (1.15) holds for $k=1\ldots,k^{\prime}-1$, then by mathematical induction and (A.1), we only need to show that

[TABLE]

which is straightforward to check since $\sum_{l=0}^{k^{\prime}}{k^{\prime}\choose l}=2^{k^{\prime}}$ and ${i-1\choose j}j!(i-1-j)!=(i-1)!$. ∎

Proof of Lemma 1.2.

Notice that the derivative $\frac{d}{dt}(\bm{J}^{k}\Psi_{t}(\bm{v}))$ can be rewritten by exchanging the order of derivation (which can be justified by the Taylor series expansion and the bounds (1.15)) as

[TABLE]

For $k=1$, the above equation implies that

[TABLE]

thus using the fact that $\|\bm{J}\Psi_{0}(\bm{v})\|=\|\bm{J}_{\bm{v}}(\bm{v})\|=1$, by Grönwall’s lemma, we have

[TABLE]

Now we are going to show that

[TABLE]

Indeed, this was shown for $k=1$ above, and for $k\geq 2$, from (A.2), it follows that

[TABLE]

thus (A.3) can be proven by mathematical induction and Grönwall’s lemma. This implies in particular our first claim, (1.21). Our second claim, (1.22) follows by the fact that $\Phi_{t}(\bm{v})=\bm{H}\Psi_{t}(\bm{v})$ is a linear transformation. ∎

Proof of Lemma 4.3.

Let $(\bm{J}\Phi_{t_{i}}(\bm{u}))_{\cdot,j}$ denote the $j$th column of the Jacobian, and let $\tilde{\bm{Z}}_{i}^{j}:=\bm{Z}_{i}^{j}/\sigma_{Z}$ (which is a standard normal random variable). Then we can write $\bm{B}_{k}$ as

[TABLE]

This is of the same form as in equation (4.1.2) of Theorem 4.1.1 of Tropp [2015]. Since $\|\bm{J}\Phi_{t_{i}}(\bm{u})\|\leq\widehat{M}_{1}(T)$, one can see that we also have $\|(\bm{J}\Phi_{t_{i}}(\bm{u}))_{\cdot,j}\|\leq\widehat{M}_{1}(T)$, and thus the variance statistics $v(\bm{Z})$ of Theorem 4.1.1 can be bounded as $v(\bm{Z})\leq\sigma_{Z}^{2}\widehat{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 $\bm{u}$ it follows that the first term in the definition of $\bm{A}_{k}$ satisfies that

[TABLE]

We can rewrite the second term in the definition (2.6) as

[TABLE]

where $H\Phi_{t_{i}}^{j}(\bm{u})$ denotes the Hessian of the function $\Phi_{t_{i}}^{j}$ at point $\bm{u}$. This is of the same form as in equation (4.1.2) of Theorem 4.1.1 of Tropp [2015]. Since $\|\bm{J}^{2}\Phi_{t_{i}}(\bm{u})\|\leq\widehat{M}_{2}(T)$, one can see that we also have $\|H\Phi_{t_{i}}^{j}(\bm{u})\|\leq\widehat{M}_{2}(T)$, and thus the variance statistics $v(\bm{Z})$ of Theorem 4.1.1 can be bounded as $v(\bm{Z})\leq\sigma_{Z}^{2}\widehat{M}_{2}(T)^{2}(k+1)d_{o}$, and thus for any $t\geq 0$, we have

[TABLE]

By the Bauer-Fike theorem (see Bauer and Fike [1960]), we have

[TABLE]

and

[TABLE]

so (4.8) follows by the bounds (A.4) and (A.5), and (4.4) follows by rearrangement. ∎

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 $\delta\in[0,\tilde{h})$,

[TABLE]

for some constant $c^{\prime}(\bm{u},\tilde{h})>0$ independent of $\delta$, that is monotone increasing in $\tilde{h}$ for $0<\tilde{h}\leq\tilde{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\left(\frac{T}{j+1},\tilde{h}_{\max}\right)$, and set $\tilde{h}:=\lfloor h_{0}/h\rfloor\cdot h$. Then for $h\leq h_{0}$, we have $\lfloor h_{0}/h\rfloor>\frac{h_{0}}{2h}$ and $\tilde{h}>h_{0}/2$, so using (A.6), we obtain that

[TABLE]

Thus Assumption 2.1 holds with $h_{\min}(\bm{u},T):=h_{0}/2$ and $c(\bm{u},T):=c^{\prime}(\bm{u},h_{0}/2)\cdot h_{0}/2$.

To complete the proof, we will now show (A.6). Using inequality (1.14), we can see that the Taylor expansion

[TABLE]

is valid for times $0\leq t<C_{\mathrm{der}}^{-1}$. Based on this expansion, assuming that $i\tilde{h}+\delta<C_{\mathrm{der}}^{-1}$, $\bm{H}\bm{D}^{i}\bm{v}$ can be approximated by a finite difference formula depending on the values of $\Phi_{\delta}(\bm{v}),\Phi_{\tilde{h}+\delta}(\bm{v}),\ldots,\Phi_{i\tilde{h}+\delta}$, with error of $O(\tilde{h})$. This finite difference formula will be denoted as

[TABLE]

The coefficients $a_{l}^{(i,\delta)}$ are explicitly defined in Fornberg [1988], and they only depend on $i$ and the ratio $\delta/\tilde{h}$. Based on the definition of these coefficients on page 700 of Fornberg [1988], we can see that

[TABLE]

i.e. they can be bounded by a finite constant independently of $\delta$. By Taylor’s expansion of the terms $\Phi_{l\tilde{h}+\delta}(\bm{v})$ around time point 0, for $l\tilde{h}+\delta<C_{\mathrm{der}}^{-1}$, we have

[TABLE]

with $b_{m}^{(i,\delta)}:=\frac{1}{m!}\cdot\sum_{l=0}^{i}a_{l}^{(i,\delta)}(l+\delta/\tilde{h})^{m}$. Due to the particular choice of the constants $a_{l}^{(i,\delta)}$, we have $b_{m}^{(i,\delta)}=0$ for $0\leq m<i$ and $b_{m}^{(i,\delta)}=1$ for $m=i$. Based on this, we can write the difference between the approximation (A.7) and the derivative explicitly as

[TABLE]

Let us denote $\tilde{\Phi}^{(i,\delta)}(\bm{v},\tilde{h}):=\sum_{m=i+1}^{\infty}\tilde{h}^{m-i-1}\cdot b_{m}^{(i,\delta)}\cdot\bm{H}\bm{D}^{m}\bm{v}$. Using inequality (1.15), and the bound $|b_{m}^{(i,\delta)}|\leq\frac{\overline{a}\cdot(i+1)^{m+1}}{m!}$, we have that for $0\leq i\leq j$, $\tilde{h}\leq\frac{1}{2(j+1)C_{\bm{J}}^{(1)}}$,

[TABLE]

Denote $C_{\mathrm{Lip}}:=2\|\bm{H}\|\overline{a}\cdot\max_{0\leq i\leq j}\left((i+1)\cdot\left((i+1)C_{\bm{J}}^{(1)}\right)^{i+1}\right)$, then we know that for every $0\leq i\leq j$, $\tilde{h}<\frac{1}{2(j+1)C_{\bm{J}}^{(1)}}$, the functions $\tilde{\Phi}^{(i)}(\bm{v},\tilde{h})$ are $C_{\mathrm{Lip}}$ - Lipschitz in $\bm{v}$ with respect to the $\|\cdot\|$ norm, and thus for every $\bm{v}\in\mathcal{B}_{R}$,

[TABLE]

By Assumption 2.2, and the boundedness of $\mathcal{B}_{R}$, it follows that there is constant $C_{\bm{D}}(\bm{u},T)>0$ such that for every $\bm{v}\in\mathcal{B}_{R}$,

[TABLE]

From equations (A.9), (A.10), and the boundedness of $\mathcal{B}_{R}$, it follows that there is a constant $C_{j}(\bm{u},T,\tilde{h})>0$ that is non-decreasing in $\tilde{h}$ such that for every $\bm{v}\in\mathcal{B}_{R}$,

[TABLE]

By the definitions (A.7) and (A.8), it follows that

[TABLE]

and thus (A.6) follows by rearrangement. ∎

The following lemma bounds the number of balls of radius $\delta$ required to cover a $d$-dimensional unit ball. It will be used in the proofs of Propositions 4.1, 4.2 and 4.3.

Lemma A.1**.**

For any $d\geq 1$, $0<\delta\leq 1$, a $d$-dimensional unit ball can be fully covered by the union of $c(d)\cdot(\frac{1}{\delta})^{d}$ balls of radius $\delta$, where

[TABLE]

Proof.

For $d\geq 3$, this follows from Theorem 1 of Dumer [2007]. For $d=1$, the result follows from the fact that $\frac{1}{\delta}+1$ intervals suffice, and $\frac{1}{\delta}+1\leq\frac{2}{\delta}$. For $d=2$, we know that the circles of radius $\delta$ contains a square of edge length $\sqrt{2}\delta$, and in order to cover a square of edge length 2 containing the unit ball, it suffices to use $\lceil 2/(\sqrt{2}\delta)\rceil^{2}\leq\left(\frac{\sqrt{2}}{\delta}+1\right)^{2}<\frac{6}{\delta^{2}}$ squares of edge length $\sqrt{2}\delta$, thus the result follows. ∎

Our next lemma shows some concentration inequalities that will be used in the proof of our propositions. It is a reformulation of Corollary 13.2 and Theorem 5.8 of Boucheron et al. [2013] to our setting.

Lemma A.2**.**

For every $l\in\mathbb{N}$, define the sets

[TABLE]

For any two elements $(r,\bm{s}_{1},\ldots,\bm{s}_{l}),(r,\bm{s}_{1}^{\prime},\ldots,\bm{s}_{l}^{\prime})\in\mathcal{T}_{l}$, we define the distance

[TABLE]

Similarly, for any two elements $(\bm{v},\bm{s}_{1},\ldots,\bm{s}_{l}),(\bm{v}^{\prime},\bm{s}_{1}^{\prime},\ldots,\bm{s}_{l}^{\prime})\in\overline{\mathcal{T}}_{l}$, we define

[TABLE]

Suppose that $\bm{Z}_{0},\ldots,\bm{Z}_{k}$ are i.i.d. $d_{o}$ dimensional standard normal random vectors, and $\varphi_{0},\ldots,\varphi_{k}:\mathcal{T}_{l}\to\mathbb{R}^{d_{o}}$ are functions that are $L$-Lipschitz with respect to the distance $d_{l}$ on $\mathcal{T}_{l}$, and satisfy that $\|\varphi_{i}(r,\bm{s}_{1},\ldots,\bm{s}_{l})\|\leq M$ for any $0\leq i\leq k$, $(r,\bm{s}_{1},\ldots,\bm{s}_{l})\in\mathcal{T}_{l}$ (the constants $M$ and $L$ can depend on $l$). Then $W_{l}:=\sup_{(r,\bm{s}_{1},\ldots,\bm{s}_{l})\in\mathcal{T}_{l}}\sum_{i=0}^{k}\left<\varphi_{i}(r,\bm{s}_{1},\ldots,\bm{s}_{l}),\bm{Z}_{i}\right>$ satisfies that for any $0<\varepsilon\leq 1$,

[TABLE]

Similarly, if $\overline{\varphi}_{0},\ldots,\overline{\varphi}_{k}:\overline{\mathcal{T}}_{l}\to\mathbb{R}^{d_{o}}$ are $L$-Lipschitz with respect to $\overline{d}_{l}$, and satisfy that $\|\overline{\varphi}_{i}(\bm{v},\bm{s}_{1},\ldots,\bm{s}_{l})\|\leq M$ for any $0\leq i\leq k$, $(\bm{v},\bm{s}_{1},\ldots,\bm{s}_{l})\in\overline{\mathcal{T}}_{l}$, then the quantity $\overline{W}_{l}:=\sup_{(\bm{v},\bm{s}_{1},\ldots,\bm{s}_{l})\in\overline{\mathcal{T}}_{l}}\sum_{i=0}^{k}\left<\overline{\varphi}_{i}(\bm{v},\bm{s}_{1},\ldots,\bm{s}_{l}),\bm{Z}_{i}\right>$ satisfies that for any $0<\varepsilon\leq 1$,

[TABLE]

Proof.

For $(r,\bm{s}_{1},\ldots,\bm{s}_{l})\in\mathcal{T}_{l}$, let us denote

[TABLE]

Then we have $W_{l}^{\prime}(0,\bm{0},\ldots,\bm{0})=0$, and by the $L$-Lipschitz assumption on $\varphi_{i}$, we can see that $W_{l}^{\prime}(r,\bm{s}_{1},\ldots,\bm{s}_{l})-W_{l}^{\prime}(0,\bm{0},\ldots,\bm{0})$ is a one dimensional Gaussian random variable whose variance is bounded by $\left(d_{l}((r,\bm{s}_{1},\ldots,\bm{s}_{l}),(0,\bm{0}))\right)^{2}$. Therefore its moment generating function can be bounded as

[TABLE]

This means that Dudley’s entropy integral expectation bound (Corollary 13.2 of Boucheron, Lugosi, and Massart [2013]) is applicable here. To apply that result, we first need to upper bound the packing number $N(\delta,\mathcal{T}_{l})$, which is the maximum number of points that can be selected in $\mathcal{T}_{l}$ such that all of them are further away from each other than $\delta$ in $d_{l}$ distance. It is easy to show that $N(\delta,\mathcal{T}_{l})\leq N^{\prime}(\delta/2,\mathcal{T}_{l})$, where $N^{\prime}(\delta/2,\mathcal{T}_{l})$ is the number of spheres of radius $\frac{\delta}{2}$ in $d_{l}$ distance needed to cover $\mathcal{T}_{l}$.

Since $\mathcal{T}_{l}\subset[0,2R]\times\mathcal{B}_{1}^{l}$, it follows that $N^{\prime}(\delta/2,\mathcal{T}_{l})\leq N^{\prime}(\delta/2,[0,2R]\times\mathcal{B}_{1}^{l})$. Moreover, due to the product nature of the space $[0,2R]\times\mathcal{B}_{1}^{l}$ and the definition of the distance $d_{l}$, if we first cover $[0,2R]$ with intervals of length $R\delta/(l+1)$ (in 1 dimensional Euclidean distance), and then cover each sphere $\mathcal{B}_{1}$ in $\mathcal{B}_{1}^{l}$ with spheres of radius $\delta/(2l+2)$ (in $d$ dimensional Euclidean distance), then the product of any such interval and spheres will be contained in a sphere of $d_{l}$-radius less than or equal to $\delta/2$, and the union of all such spheres will cover $[0,2R]\times\mathcal{B}_{1}^{l}$. Therefore using Lemma A.1, we obtain that for $0<\delta\leq 2l+2$,

[TABLE]

and using the fact that $\log(c(d))\leq d$ for any $d\in\mathbb{Z}_{+}$, we have

[TABLE]

Using this, and the fact that the maximum $d_{l}$ distance of any point in $\mathcal{T}_{l}$ from $(0,\bm{0},\ldots,\bm{0})$ is bounded by $l+1$, by Corollary 13.2 of Boucheron, Lugosi, and Massart [2013], we obtain that

[TABLE]

By the definition of $W_{l}^{\prime}$, this implies that $\mathbb{E}(W_{l})\leq 11(l+1)\sqrt{ld+1}\cdot L\sqrt{k+1}$. Moreover, it is easy to check that the conditions of Theorem 5.8 of Boucheron, Lugosi, and Massart [2013] hold for $W_{l}$, and thus for any $t\geq 0$,

[TABLE]

with

[TABLE]

and (A.16) follows. The proof of (A.17) is similar. ∎

Proof of Proposition 4.1.

Using Assumption 2.1, we know that

[TABLE]

thus it suffices to lower bound the terms $2\sum_{i=0}^{k}\left<\Phi_{t_{i}}(\bm{v})-\Phi_{t_{i}}(\bm{u}),\bm{Z}_{i}\right>$.

Let $\varphi_{i}(r,\bm{s}):=\frac{\Phi_{t_{i}}(\bm{u})-\Phi_{t_{i}}(\bm{u}+r\bm{s})}{r}$ for $(r,\bm{s})\in\mathcal{T}_{1},r>0$ ($\mathcal{T}_{1}$ was defined as in Lemma A.2). We continuously extend it to $r=0$ as

[TABLE]

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,\bm{s})\in\mathcal{T}}\sum_{i=0}^{k}\left<\varphi_{i}(r,\bm{s}),\bm{Z}_{i}\right>$, since

[TABLE]

Now we are going to obtain bounds on the constants $L$ and $M$ of Lemma A.2. We have

[TABLE]

thus the $\varphi_{i}(r,\bm{s})$ is $L$-Lipschitz with respect to the $d_{1}$ distance for $L:=\widehat{M}_{2}(T)R+\widehat{M}_{1}(T)$. Moreover, from the definition of $\varphi_{i}(r,\bm{s})$, by (1.19), it follows that $\|\varphi_{i}(r,\bm{s})\|\leq M$ for $M:=\widehat{M}_{1}(T)$. The claim of the proposition now follows by Lemma A.2. ∎

Proof of Proposition 4.2.

From the definitions, we have

[TABLE]

thus

[TABLE]

The first term in the right hand side of the above inequality can be upper bounded as

[TABLE]

where we have used the multivariate Taylor’s expansion bound (1.19).

For the second term in (A.21), for $(r,\bm{s})\in\mathcal{T}_{1}$ (defined as in Lemma A.2), $r>0$, let

[TABLE]

For $r=0$, this can be continuously extended as

[TABLE]

Similarly to Lemma A.2, we define $W_{1}:=\sup_{(r,\bm{s})\in\mathcal{T}_{1}}\sum_{i=0}^{k}\left<\varphi_{i}(r,\bm{s}),\bm{Z}_{i}\right>$, and

$W_{1}^{\prime}:=\sup_{(r,\bm{s})\in\mathcal{T}_{1}}\sum_{i=0}^{k}\left<-\varphi_{i}(r,\bm{s}),\bm{Z}_{i}\right>$, then the second term in (A.21) can be bounded as

[TABLE]

Based on (1.19), the partial derivatives of $\varphi_{i}(r,\bm{s})$ satisfy that

[TABLE]

therefore $\varphi_{i}$ is $L$-Lipschitz with respect to the distance $d_{1}$ for $L:=\frac{1}{2}\widehat{M}_{3}(T)+\frac{1}{2}\widehat{M}_{4}(T)R$, and we have $\|\varphi_{i}(r,\bm{s})\|\leq M$ for $M:=\frac{1}{6}\widehat{M}_{3}(T)$. The claim of the proposition now follows by applying Lemma A.2 to $W_{1}$ and $W_{1}^{\prime}$ separately (for $\varepsilon/2$ instead of $\varepsilon$), and then using the union bound. ∎

Proof of Proposition 4.3.

From the definitions, we have

[TABLE]

thus

[TABLE]

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

[TABLE]

For the last term in (A.23), for $(r,\bm{s}_{1},\bm{s}_{2})\in\mathcal{T}_{2}$ (defined as in Lemma A.2), $r>0$, we let

[TABLE]

We extend this continuously to $r=0$ as

[TABLE]

We define $W_{2}$ as in Lemma A.2 as $W_{2}:=\sup_{(r,\bm{s}_{1},\bm{s}_{2})\in\mathcal{T}_{2}}\sum_{i=0}^{k}\left<\varphi_{i}(r,\bm{s}_{1},\bm{s}_{2}),\bm{Z}_{i}\right>$, then the last term of (A.23) can be bounded as

[TABLE]

By (1.19), for any $(r,\bm{s}_{1},\bm{s}_{2})\in\mathcal{T}_{2}$, the partial derivatives of $\varphi_{i}$ satisfy that

[TABLE]

Based on these bounds, we can see that $\varphi_{i}$ is $L$-Lipschitz with respect to the $d_{2}$ distance for $L:=2R\widehat{M}_{4}(T)+\widehat{M}_{3}(T)$, and it satisfies that $\|\varphi_{i}(r,\bm{s}_{1},\bm{s}_{2})\|\leq M$ for $M:=\frac{1}{2}\widehat{M}_{3}(T)$. The claim of the proposition now follows by Lemma A.2. ∎

Proof of Proposition 5.2.

By (4.1) and (4.3), we have

[TABLE]

We first study the deterministic terms. Notice that

[TABLE]

By Assumption 2.1, it follows that $\sum_{i=0}^{k}\bm{J}\Phi_{t_{i}}(\bm{u})^{\prime}\cdot\bm{J}\Phi_{t_{i}}(\bm{u})\succeq\frac{c(\bm{u},T)}{h}\cdot\bm{I}_{d}$, thus for every $\bm{v}\in B(\bm{u},r_{H}(\bm{u},T))$, we have

[TABLE]

For the random terms, we first define $\overline{\varphi}_{i}:\overline{\mathcal{T}}_{2}\to\mathbb{R}$ ($\overline{\mathcal{T}}_{2}$ was defined in Lemma A.2) for $0\leq i\leq k$ as

[TABLE]

Based on these, we let $\overline{W}_{2}:=\sup_{(\bm{v},\bm{s}_{1},\bm{s}_{2})\in\overline{\mathcal{T}}_{2}}\sum_{i=0}^{k}\left<\overline{\varphi}_{i}(\bm{v},\bm{s}_{1},\bm{s}_{2}),\bm{Z}_{i}\right>$, then the random terms can be bounded as

[TABLE]

By its definition, it is easy to see that for every $0\leq i\leq k$, $\overline{\varphi}_{i}$ is $L$-Lipschitz with respect to the $\overline{d}_{2}$ distance for $L:=\widehat{M}_{2}(T)+\widehat{M}_{3}(T)R$, and that $\|\overline{\varphi}_{i}(\bm{v},\bm{s}_{1},\bm{s}_{2})\|\leq M$ for $M:=\widehat{M}_{2}(T)$ for every $(\bm{v},\bm{s}_{1},\bm{s}_{2})\in\overline{\mathcal{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

[TABLE]

Based on the assumption on $q$, and (1.19), the deterministic terms can be bounded as

[TABLE]

For the random terms, we first define $\overline{\varphi}_{i}:\overline{\mathcal{T}}_{3}\to\mathbb{R}$ ($\overline{\mathcal{T}}_{3}$ was defined in Lemma A.2) for $0\leq i\leq k$ as

[TABLE]

Based on this, we let $\overline{W}_{3}:=\sup_{(\bm{v},\bm{s}_{1},\bm{s}_{2},\bm{s}_{3})\in\overline{\mathcal{T}}_{3}}\sum_{i=0}^{k}\left<\overline{\varphi}_{i}(\bm{v},\bm{s}_{1},\bm{s}_{2},\bm{s}_{3}),\bm{Z}_{i}\right>$, then one can see that the random terms can be bounded as

[TABLE]

By its definition, it is easy to see that for every $0\leq i\leq k$, $\overline{\varphi}_{i}$ is $L$-Lipschitz with respect to the $\overline{d}_{3}$ distance for $L:=\widehat{M}_{3}(T)+\widehat{M}_{4}(T)R$, and that $\|\overline{\varphi}_{i}(\bm{v},\bm{s}_{1},\bm{s}_{2},\bm{s}_{3})\|\leq M$ for $M:=\widehat{M}_{3}(T)$ for every $(\bm{v},\bm{s}_{1},\bm{s}_{2},\bm{s}_{3})\in\overline{\mathcal{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 $\bm{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 $\bm{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 $\bm{v}\in\mathcal{B}_{R}$, and the trapping ball assumption (1.3) holds. Then for any $0\leq t<C_{\mathrm{der}}^{-1}$, the series

[TABLE]

is convergent, well defined, and satisfies that $\Psi_{t}(\Psi_{-t}(\bm{v}))=\bm{v}$ and that for any $i_{\max}\in\mathbb{N}$,

[TABLE]

Proof.

The result follows from the bounds (1.14), and the definition of equation (1.1). ∎

Based on this lemma, given the observations $\bm{Y}_{0:k}$, we propose the following initial estimator. First, select some intermediate indices $0=i_{1}<i_{2}<\ldots<i_{m}<k$ satisfying that $i_{m}\cdot h<C_{\mathrm{der}}^{-1}$. For each index $i_{r}$, we compute the derivative estimates $\hat{\Phi}^{(l)}$ of $\bm{D}^{l}(\bm{u}(t_{i_{r}}))$, and then use the function $F$ described in Section 3.2 to obtain initial estimators $\widehat{\bm{u}(t_{i_{r}})}$ of $\bm{u}(t_{i_{r}})$ for $0\leq r\leq m$. After this, we project these estimators to $\mathcal{B}_{R}$ (see (3.13)), and run them backwards by $t_{i_{r}}$ time units via the approximation (A.27), and project them back to $\mathcal{B}_{R}$, that is, for some sufficiently large $i_{\max}\in\mathbb{N}$, we let

[TABLE]

The final initial estimator $\hat{\bm{u}}$ is then chosen as the one among $(\hat{\bm{u}}^{r})_{0\leq r\leq m}$ that has the largest the a-posteriori probability $\mu^{\mathrm{sm}}(\hat{\bm{u}}^{r}|\bm{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 $\sigma_{Z}\sqrt{h}$ and $h$ are sufficiently small, and $i_{\max}$ is sufficiently large, as long as at least one of the vectors $\left(\bm{u}(t_{i_{r}})\right)_{0\leq r\leq 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\in[0,T]$ such that none of the coefficients of $\bm{u}(t)$ are 0.