2.1 Prediction

We assume under the chosen computational setting that we have access to $M$ samples $z_{0}^{i}\in \mathbb{R}^{N_{z}}$ , $i=1,\ldots ,M$ , from the filtering distribution at $t=0$ . We also assume that we know (explicitly or implicitly) the forward transition probabilities $q_{+}(z_{1}|z_{0}^{i})$ of the underlying Markovian stochastic process. This leads to the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ as given by (2.4).

Before we consider two specific examples, we introduce two concepts related to the forward transition kernel which we will need later in order to address scenarios (B) & (C) from Definition 2.4.

We first introduce the backward transition kernel $q_{-}(z_{0}|z_{1})$ , which is defined via the equation

$$\begin{eqnarray}q_{-}(z_{0}|z_{1})\,\unicode[STIX]{x1D70B}_{1}(z_{1})=q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}_{0}(z_{0}).\end{eqnarray}$$

Note that $q_{-}(z_{0}|z_{1})$ as well as $\unicode[STIX]{x1D70B}_{0}$ are not absolutely continuous with respect to the underlying Lebesgue measure, that is,

(2.15) $$\begin{eqnarray}q_{-}(z_{0}|z_{1})=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\displaystyle \frac{q_{+}(z_{1}|z_{0}^{i})}{\unicode[STIX]{x1D70B}_{1}(z_{1})}\,\unicode[STIX]{x1D6FF}(z_{0}-z_{0}^{i}).\end{eqnarray}$$

The backward transition kernel $q_{-}(z_{1}|z_{0})$ reverses the prediction process in the sense that

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{0}(z_{0})=\int q_{-}(z_{0}|z_{1})\,\unicode[STIX]{x1D70B}_{1}(z_{1})\,\text{d}z_{1}.\end{eqnarray}$$

Remark 2.7. Let us assume that the detailed balance

$$\begin{eqnarray}q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}(z_{0})=q_{+}(z_{0}|z_{1})\,\unicode[STIX]{x1D70B}(z_{1})\end{eqnarray}$$

holds for some PDF $\unicode[STIX]{x1D70B}$ and forward transition kernel $q_{+}(z_{1}|z_{0})$ . Then $\unicode[STIX]{x1D70B}_{1}=\unicode[STIX]{x1D70B}$ for $\unicode[STIX]{x1D70B}_{0}=\unicode[STIX]{x1D70B}$ (invariance of $\unicode[STIX]{x1D70B}$ ) and $q_{-}(z_{0}|z_{1})=q_{+}(z_{1}|z_{0})$ .

We next introduce a class of forward transition kernels using the concept of twisting (Guarniero et al. Reference Guarniero, Johansen and Lee2017, Heng et al. Reference Heng, Bishop, Deligiannidis and Doucet2018), which is an application of Doob’s H-transform technique (Doob Reference Doob1984).

Definition 2.8. Given a non-negative twisting function $\unicode[STIX]{x1D713}(z_{1})$ such that the modified transition kernel (1.2) is well-defined, one can define the twisted forecast PDF

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1}):=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0}^{i})=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\displaystyle \frac{\unicode[STIX]{x1D713}(z_{1})}{\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})}\,q_{+}(z_{1}|z_{0}^{i}).\end{eqnarray}$$

The PDFs $\unicode[STIX]{x1D70B}_{1}$ and $\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}$ are related by

(2.17) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D70B}_{1}(z_{1})}{\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1})}=\displaystyle \frac{\mathop{\sum }_{i=1}^{M}q_{+}(z_{1}|z_{0}^{i})}{\mathop{\sum }_{i=1}^{M}{\textstyle \frac{\unicode[STIX]{x1D713}(z_{1})}{\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})}}\,q_{+}(z_{1}|z_{0}^{i})}.\end{eqnarray}$$

Equation (2.17) gives rise to importance weights

(2.18) $$\begin{eqnarray}w^{i}\propto \displaystyle \frac{\unicode[STIX]{x1D70B}_{1}(z_{1}^{i})}{\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1}^{i})}\end{eqnarray}$$

for samples $z_{1}^{i}=Z_{1}^{i}(\unicode[STIX]{x1D714})$ drawn from the twisted forecast PDF, that is,

$$\begin{eqnarray}Z_{1}^{i}\sim q_{+}^{\unicode[STIX]{x1D713}}(\cdot \,|z_{0}^{i})\end{eqnarray}$$

and

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{1}(z)\approx \displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}w^{i}\,\unicode[STIX]{x1D6FF}(z-z_{1}^{i})\end{eqnarray}$$

in a weak sense. Here we have assumed that the normalization constant in (2.18) is chosen such that

(2.19) $$\begin{eqnarray}\mathop{\sum }_{i=1}^{M}w^{i}=M.\end{eqnarray}$$

Such twisted transition kernels will become important when looking at the filtering and smoothing as well as the Schrödinger problem later in this section.

Let us now discuss a couple of specific models which give rise to transition kernels $q_{+}(z_{1}|z_{0})$ . These models will be used throughout this paper to illustrate mathematical and algorithmic concepts.

2.1.1 Gaussian model error

Let us consider the discrete-time stochastic process

(2.20) $$\begin{eqnarray}Z_{1}=\unicode[STIX]{x1D6F9}(Z_{0})+\unicode[STIX]{x1D6FE}^{1/2}\unicode[STIX]{x1D6EF}_{0}\end{eqnarray}$$

for given map $\unicode[STIX]{x1D6F9}:\mathbb{R}^{N_{z}}\rightarrow \mathbb{R}^{N_{z}}$ , scaling factor $\unicode[STIX]{x1D6FE}>0$ , and Gaussian distributed random variable $\unicode[STIX]{x1D6EF}_{0}$ with mean zero and covariance matrix $B\in \mathbb{R}^{N_{z}\times N_{z}}$ . The associated forward transition kernel is given by

(2.21) $$\begin{eqnarray}q_{+}(z_{1}|z_{0})=\text{n}(z_{1};\unicode[STIX]{x1D6F9}(z_{0}),\unicode[STIX]{x1D6FE}B).\end{eqnarray}$$

Recall that we have introduced the shorthand $\text{n}(z;\bar{z},P)$ for the PDF of a Gaussian random variable with mean $\bar{z}$ and covariance matrix  $P$ .

Let us consider a twisting potential $\unicode[STIX]{x1D713}$ of the form

$$\begin{eqnarray}\unicode[STIX]{x1D713}(z_{1})\propto \exp \biggl(-\displaystyle \frac{1}{2}(Hz_{1}-d)^{\text{T}}R^{-1}(Hz_{1}-d)\biggr)\end{eqnarray}$$

for given $H\in \mathbb{R}^{N_{z}\times N_{d}}$ , $d\in \mathbb{R}^{N_{d}}$ , and covariance matrix $R\in \mathbb{R}^{N_{d}\times N_{d}}$ . We define

(2.22) $$\begin{eqnarray}K:=BH^{\text{T}}(HBH^{\text{T}}+\unicode[STIX]{x1D6FE}^{-1}R)^{-1}\end{eqnarray}$$

and

(2.23) $$\begin{eqnarray}\bar{B}:=B-KHB,\quad \bar{z}_{1}^{i}:=\unicode[STIX]{x1D6F9}(z_{0}^{i})-K(H\unicode[STIX]{x1D6F9}(z_{0}^{i})-d).\end{eqnarray}$$

The twisted transition kernels are given by

$$\begin{eqnarray}q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0}^{i})=\text{n}(z_{1};\bar{z}_{1}^{i},\unicode[STIX]{x1D6FE}\bar{B})\end{eqnarray}$$

and

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})\propto \exp \biggl(-\displaystyle \frac{1}{2}(H\unicode[STIX]{x1D6F9}(z_{0}^{i})-d)^{\text{T}}(R+\unicode[STIX]{x1D6FE}HBH^{\text{T}})^{-1}(H\unicode[STIX]{x1D6F9}(z_{0}^{i})-d)\biggr)\end{eqnarray}$$

for $i=1,\ldots ,M$ .

2.1.2 SDE models

Consider the (forward) SDE (Pavliotis Reference Pavliotis2014)

(2.24) $$\begin{eqnarray}\,\text{d}Z_{t}^{+}=f_{t}(Z_{t}^{+})\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{+}\end{eqnarray}$$

with initial condition $Z_{0}^{+}=z_{0}$ and $\unicode[STIX]{x1D6FE}>0$ . Here $W_{t}^{+}$ stands for standard Brownian motion in the sense that the distribution of $W_{t+\unicode[STIX]{x1D6E5}t}^{+}$ , $\unicode[STIX]{x1D6E5}t>0$ , conditioned on $w_{t}^{+}=W_{t}^{+}(\unicode[STIX]{x1D714})$ is Gaussian with mean $w_{t}^{+}$ and covariance matrix $\unicode[STIX]{x1D6E5}t\,I$ (Pavliotis Reference Pavliotis2014) and the process $Z_{t}^{+}$ is adapted to $W_{t}^{+}$ .

The resulting time- $t$ transition kernels $q_{t}^{+}(z|z_{0})$ from time zero to time $t$ , $t\in (0,1]$ , satisfy the Fokker–Planck equation (Pavliotis Reference Pavliotis2014)

$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}q_{t}^{+}(\cdot \,|z_{0})=-\unicode[STIX]{x1D6FB}_{z}\cdot (q_{t}^{+}(\cdot \,|z_{0})f_{t})+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}q_{t}^{+}(\cdot \,|z_{0})\end{eqnarray}$$

with initial condition $q_{0}^{+}(z|z_{0})=\unicode[STIX]{x1D6FF}(z-z_{0})$ , and the time-one forward transition kernel $q_{+}(z_{1}|z_{0})$ is given by

$$\begin{eqnarray}q_{+}(z_{1}|z_{0})=q_{1}^{+}(z_{1}|z_{0}).\end{eqnarray}$$

We introduce the operator ${\mathcal{L}}_{t}$ by

$$\begin{eqnarray}{\mathcal{L}}_{t}g:=\unicode[STIX]{x1D6FB}_{z}g\cdot f_{t}+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}g\end{eqnarray}$$

and its adjoint ${\mathcal{L}}_{t}^{\dagger }$ by

(2.25) $$\begin{eqnarray}{\mathcal{L}}_{t}^{\dagger }\unicode[STIX]{x1D70B}:=-\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D70B}\,f_{t})+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}\unicode[STIX]{x1D70B}\end{eqnarray}$$

(Pavliotis Reference Pavliotis2014). We call ${\mathcal{L}}_{t}^{\dagger }$ the Fokker–Planck operator and ${\mathcal{L}}_{t}$ the generator of the Markov process associated to the SDE (2.24).

Solutions (realizations) $z_{[0,1]}=Z_{[0,1]}^{+}(\unicode[STIX]{x1D714})$ of the SDE (2.24) with initial conditions drawn from $\unicode[STIX]{x1D70B}_{0}$ are continuous functions of time, that is, $z_{[0,1]}\in {\mathcal{C}}:=C([0,1],\mathbb{R}^{N_{z}})$ , and define a probability measure $\mathbb{Q}$ on ${\mathcal{C}}$ , that is,

$$\begin{eqnarray}Z_{[0,1]}^{+}\sim \mathbb{Q}.\end{eqnarray}$$

We note that the marginal distributions $\unicode[STIX]{x1D70B}_{t}$ of $\mathbb{Q}$ , given by

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{t}(z_{t})=\int q_{t}^{+}(z_{t}|z_{0})\,\unicode[STIX]{x1D70B}_{0}(z_{0})\,\text{d}z_{0},\end{eqnarray}$$

also satisfy the Fokker–Planck equation, that is,

(2.26) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70B}_{t}={\mathcal{L}}_{t}^{\dagger }\,\unicode[STIX]{x1D70B}_{t}=-\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D70B}_{t}\,f_{t})+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}\unicode[STIX]{x1D70B}_{t}\end{eqnarray}$$

for given PDF $\unicode[STIX]{x1D70B}_{0}$ at time $t=0$ .

Furthermore, we can rewrite the Fokker–Planck equation (2.26) in the form

(2.27) $$\begin{eqnarray}\unicode[STIX]{x2202}\unicode[STIX]{x1D70B}_{t}=-\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D70B}_{t}(f_{t}-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D70B}_{t}))-\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}\unicode[STIX]{x1D70B}_{t},\end{eqnarray}$$

which allows us to read off from (2.27) the backward SDE

(2.28) $$\begin{eqnarray}\displaystyle \,\text{d}Z_{t}^{-} & = & \displaystyle f_{t}(Z_{t}^{-})\,\text{d}t-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D70B}_{t}\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{-},\nonumber\\ \displaystyle & = & \displaystyle b_{t}(Z_{t}^{-})\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{-}\end{eqnarray}$$

with final condition $Z_{1}^{-}\sim \unicode[STIX]{x1D70B}_{1}$ , $W_{t}^{-}$ backward Brownian motion, and density-dependent drift term

$$\begin{eqnarray}b_{t}(z):=f_{t}(z)-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D70B}_{t}\end{eqnarray}$$

(Nelson Reference Nelson1984, Chen et al. Reference Chen, Georgiou and Pavon2014). Here backward Brownian motion is to be understood in the sense that the distribution of $W_{t-\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D70F}}^{-}$ , $\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D70F}>0$ , conditioned on $w_{t}^{-}=W_{t}^{-}(\unicode[STIX]{x1D714})$ is Gaussian with mean $w_{t}^{-}$ and covariance matrix $\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D70F}\,I$ and all other properties of Brownian motion appropriately adjusted. The process $Z_{t}^{-}$ is adapted to $W_{t}^{-}$ .

Lemma 2.9. The backward SDE (2.28) induces a corresponding backward transition kernel from time one to time $t=1-\unicode[STIX]{x1D70F}$ with $\unicode[STIX]{x1D70F}\in [0,1]$ , denoted by $q_{\unicode[STIX]{x1D70F}}^{-}(z|z_{1})$ , which satisfies the time-reversed Fokker–Planck equation

$$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}q_{\unicode[STIX]{x1D70F}}^{-}(\cdot \,|z_{1})=\unicode[STIX]{x1D6FB}_{z}\cdot (q_{\unicode[STIX]{x1D70F}}^{-}(\cdot \,|z_{1})\,b_{1-\unicode[STIX]{x1D70F}})+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}q_{\unicode[STIX]{x1D70F}}^{-}(\cdot \,|z_{1})\end{eqnarray}$$

with boundary condition $q_{0}^{-}(z|z_{1})=\unicode[STIX]{x1D6FF}(z-z_{1})$ at $\unicode[STIX]{x1D70F}=0$ (or, equivalently, at $t=1$ ). The induced backward transition kernel $q_{-}(z_{0}|z_{1})$ is then given by

$$\begin{eqnarray}q_{-}(z_{0}|z_{1})=q_{1}^{-}(z_{0}|z_{1})\end{eqnarray}$$

and satisfies (2.15).

Proof. The lemma follows from the fact that the backward SDE (2.28) implies the Fokker–Planck equation (2.27) and that we have reversed time by introducing $\unicode[STIX]{x1D70F}=1-t$ .◻

Remark 2.10. The notion of a backward SDE also arises in a different context where the driving Brownian motion is still adapted to the past, i.e.  $W_{t}^{+}$ in our notation, and a final condition is prescribed as for (2.28). See (2.54) below as well as Appendix A.4 and Carmona (Reference Carmona2016) for more details.

We note that the mean-field equation,

(2.29) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}z_{t}=f_{t}(z_{t})-\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D70B}_{t}(z_{t})=\displaystyle \frac{1}{2}(f_{t}(z_{t})+b_{t}(z_{t})),\end{eqnarray}$$

resulting from (2.26), leads to the same marginal distributions $\unicode[STIX]{x1D70B}_{t}$ as the forward and backward SDEs, respectively. It should be kept in mind, however, that the path measure generated by (2.29) is different from the path measure $\mathbb{Q}$ generated by (2.24).

Please also note that the backward SDE and the mean-field equation (2.29) become singular as $t\rightarrow 0$ for the given initial PDF (2.3). A meaningful solution can be defined via regularization of the Dirac delta function, that is,

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{0}(z)\approx \displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\text{n}(z;z_{0}^{i},\unicode[STIX]{x1D716}I),\end{eqnarray}$$

and taking the limit $\unicode[STIX]{x1D716}\rightarrow 0$ .

We will find later that it can be advantageous to modify the given SDE (2.24) by a time-dependent drift term $u_{t}(z)$ , that is,

(2.30) $$\begin{eqnarray}\,\text{d}Z_{t}^{+}=f_{t}(Z_{t}^{+})\,\text{d}t+u_{t}(Z_{t}^{+})\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{+}.\end{eqnarray}$$

In particular, such a modification leads to the time-continuous analogue of the twisted transition kernel (1.2) introduced in Section 2.1.

Lemma 2.11. Let $\unicode[STIX]{x1D713}_{t}(z)$ denote the solutions of the backward Kolmogorov equation

(2.31) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}_{t}=-{\mathcal{L}}_{t}\unicode[STIX]{x1D713}_{t}=-\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D713}_{t}\cdot f_{t}-\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}\unicode[STIX]{x1D713}_{t}\end{eqnarray}$$

for given final $\unicode[STIX]{x1D713}_{1}(z)>0$ and $t\in [0,1]$ . The controlled SDE (2.30) with

(2.32) $$\begin{eqnarray}u_{t}(z):=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D713}_{t}(z)\end{eqnarray}$$

leads to a time-one forward transition kernel $q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0})$ which satisfies

$$\begin{eqnarray}q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0})=\unicode[STIX]{x1D713}_{1}(z_{1})\,q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D713}_{0}(z_{0})^{-1},\end{eqnarray}$$

where $q_{+}(z_{1}|z_{0})$ denotes the time-one forward transition kernel of the uncontrolled forward SDE (2.24).

Proof. A proof of this lemma has, for example, been given by Dai Pra (Reference Dai Pra1991, Theorem 2.1).◻

More generally, the modified forward SDE (2.30) with $Z_{0}^{+}\sim \unicode[STIX]{x1D70B}_{0}$ generates a path measure which we denote by $\mathbb{Q}^{u}$ for given functions $u_{t}(z)$ , $t\in [0,1]$ . Realizations of this path measure are denoted by $z_{[0,1]}^{u}$ . According to Girsanov’s theorem (Pavliotis Reference Pavliotis2014), the two path measures $\mathbb{Q}$ and $\mathbb{Q}^{u}$ are absolutely continuous with respect to each other, with Radon–Nikodym derivative

(2.33) $$\begin{eqnarray}\displaystyle \frac{\text{d}\mathbb{Q}^{u}}{\text{d}\mathbb{Q}}_{|z_{[0,1]}^{u}}=\exp \biggl(\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}}\int _{0}^{1}(\Vert u_{t}\Vert ^{2}\,\text{d}t+2\unicode[STIX]{x1D6FE}^{1/2}u_{t}\cdot \text{d}W_{t}^{+})\biggr),\end{eqnarray}$$

provided that the Kullback–Leibler divergence $\text{KL}(\mathbb{Q}^{u}||\mathbb{Q})$ between $\mathbb{Q}^{u}$ and $\mathbb{Q}$ , given by

(2.34) $$\begin{eqnarray}\text{KL}(\mathbb{Q}^{u}||\mathbb{Q}):=\int \biggl[\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}}\int _{0}^{1}\Vert u_{t}\Vert ^{2}\,\text{d}t\biggr]\mathbb{Q}^{u}(\text{d}z_{[0,1]}^{u}),\end{eqnarray}$$

is finite. Recall that the Kullback–Leibler divergence between two path measures $\mathbb{P}\ll \mathbb{Q}$ on ${\mathcal{C}}$ is defined by

$$\begin{eqnarray}\text{KL}(\mathbb{P}||\mathbb{Q})=\int \log \displaystyle \frac{\text{d}\mathbb{P}}{\text{d}\mathbb{Q}}\,\mathbb{P}(\text{d}z_{[0,1]}).\end{eqnarray}$$

If the modified SDE (2.30) is used to make predictions, then its solutions $z_{[0,1]}^{u}$ need to be weighted according to the inverse Radon–Nikodym derivative

(2.35) $$\begin{eqnarray}\displaystyle \frac{\text{d}\mathbb{Q}}{\text{d}\mathbb{Q}^{u}}_{|z_{[0,1]}^{u}}=\exp \biggl(-\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}}\int _{0}^{1}(\Vert u_{t}\Vert ^{2}\,\text{d}t+2\unicode[STIX]{x1D6FE}^{1/2}u_{t}\cdot \,\text{d}W_{t}^{+})\biggr)\end{eqnarray}$$

in order to reproduce the desired forecast PDF $\unicode[STIX]{x1D70B}_{1}$ of the original SDE (2.24).

Remark 2.12. A heuristic derivation of equation (2.33) can be found in Section 3.1.2 below, where we discuss the numerical approximation of SDEs by the Euler–Maruyama method. Equation (2.34) follows immediately from (2.33) by noting that the expectation of Brownian motion under the path measure $\mathbb{Q}^{u}$ is zero.

2.2 Filtering and smoothing

We now incorporate the likelihood

$$\begin{eqnarray}l(z_{1})=\unicode[STIX]{x1D70B}(y_{1}|z_{1})\end{eqnarray}$$

of the data $y_{1}$ at time $t_{1}=1$ . Bayes’ theorem tells us that, given the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ at time $t_{1}$ , the posterior PDF $\widehat{\unicode[STIX]{x1D70B}}_{1}$ is given by (2.5). The distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ solves the filtering problem at time $t_{1}$ given the data $y_{1}$ . We also recall the definition of the evidence (2.7). The quantity ${\mathcal{F}}=-\log \unicode[STIX]{x1D6FD}$ is called the free energy in statistical physics (Hartmann et al. Reference Hartmann, Richter, Schütte and Zhang2017).

An appropriate transition kernel $q_{1}(\widehat{z}_{1}|z_{1})$ , satisfying (2.11), is required in order to complete the transition from $\unicode[STIX]{x1D70B}_{0}$ to $\widehat{\unicode[STIX]{x1D70B}}_{1}$ following scenario (A) from Definition 2.4. A suitable framework for finding such transition kernels is via the theory of optimal transportation (Villani Reference Villani2003). More specifically, let $\unicode[STIX]{x1D6F1}_{\text{c}}$ denote the set of all joint probability measures $\unicode[STIX]{x1D70B}(z_{1},\widehat{z}_{1})$ with marginals

$$\begin{eqnarray}\int \unicode[STIX]{x1D70B}(z_{1},\widehat{z}_{1})\,\text{d}\widehat{z}_{1}=\unicode[STIX]{x1D70B}_{1}(z_{1}),\quad \int \unicode[STIX]{x1D70B}(z_{1},\widehat{z}_{1})\,\text{d}z_{1}=\widehat{\unicode[STIX]{x1D70B}}_{1}(\widehat{z}_{1}).\end{eqnarray}$$

We seek the joint measure $\unicode[STIX]{x1D70B}^{\ast }(z_{1},\widehat{z}_{1})\in \unicode[STIX]{x1D6F1}_{\text{c}}$ which minimizes the expected Euclidean distance between the two associated random variables $Z_{1}$ and $\widehat{Z}_{1}$ , that is,

(2.36) $$\begin{eqnarray}\unicode[STIX]{x1D70B}^{\ast }=\arg \inf _{\unicode[STIX]{x1D70B}\in \unicode[STIX]{x1D6F1}_{\text{c}}}\int \int \Vert z_{1}-\widehat{z}_{1}\Vert ^{2}\,\unicode[STIX]{x1D70B}(z_{1},\widehat{z}_{1})\,\text{d}z_{1}\,\text{d}\widehat{z}_{1}.\end{eqnarray}$$

The minimizing joint measure is of the form

(2.37) $$\begin{eqnarray}\unicode[STIX]{x1D70B}^{\ast }(z_{1},\widehat{z}_{1})=\unicode[STIX]{x1D6FF}(\widehat{z}_{1}-\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D6F7}(z_{1}))\,\unicode[STIX]{x1D70B}_{1}(z_{1})\end{eqnarray}$$

with suitable convex potential $\unicode[STIX]{x1D6F7}$ under appropriate conditions on the PDFs $\unicode[STIX]{x1D70B}_{1}$ and $\widehat{\unicode[STIX]{x1D70B}}_{1}$ (Villani Reference Villani2003). These conditions are satisfied for dynamical systems with Gaussian model errors and typical SDE models. Once the potential $\unicode[STIX]{x1D6F7}$ (or an approximation) is available, samples $z_{1}^{i}$ , $i=1,\ldots ,M$ , from the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ can be converted into samples $\widehat{z}_{1}^{i}$ , $i=1,\ldots ,M$ , from the filtering distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ via the deterministic transformation

(2.38) $$\begin{eqnarray}\widehat{z}_{1}^{i}=\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D6F7}(z_{1}^{i}).\end{eqnarray}$$

We will discuss in Section 3 how to approximate the transformation (2.38) numerically. We will find that many popular data assimilation schemes, such as the ensemble Kalman filter, can be viewed as approximations to (2.38) (Reich and Cotter Reference Reich and Cotter2015).

We recall at this point that classical particle filters start from the importance weights

$$\begin{eqnarray}w^{i}\propto \displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1}^{i})}{\unicode[STIX]{x1D70B}_{1}(z_{1}^{i})}=\displaystyle \frac{l(z_{1}^{i})}{\unicode[STIX]{x1D6FD}},\end{eqnarray}$$

and obtain the desired samples $\widehat{z}^{i}$ by an appropriate resampling with replacement scheme (Doucet et al. Reference Doucet, de Freitas and Gordon2001, Arulampalam et al. Reference Arulampalam, Maskell, Gordon and Clapp2002, Douc and Cappe Reference Douc and Cappe2005) instead of applying a deterministic transformation of the form (2.38).

Remark 2.13. If one replaces the forward transition kernel $q_{+}(z_{1}|z_{0})$ with a twisted kernel (1.2), then, using (2.17), the filtering distribution (2.5) satisfies

(2.39) $$\begin{eqnarray}\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})}{\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1})}=\displaystyle \frac{l(z_{1})\mathop{\sum }_{j=1}^{M}q_{+}(z_{1}|z_{0}^{j})}{\unicode[STIX]{x1D6FD}\mathop{\sum }_{j=1}^{M}{\textstyle \frac{\unicode[STIX]{x1D713}(z_{1})}{\widehat{\unicode[STIX]{x1D713}}(z_{0}^{j})}}\,q_{+}(z_{1}|z_{0}^{j})}.\end{eqnarray}$$

Hence drawing samples $z_{1}^{i}$ , $i=1,\ldots ,M$ , from $\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}$ instead of $\unicode[STIX]{x1D70B}_{1}$ leads to modified importance weights

(2.40) $$\begin{eqnarray}w^{i}\propto \displaystyle \frac{l(z_{1}^{i})\mathop{\sum }_{j=1}^{M}q_{+}(z_{1}^{i}|z_{0}^{j})}{\unicode[STIX]{x1D6FD}\mathop{\sum }_{j=1}^{M}{\textstyle \frac{\unicode[STIX]{x1D713}(z_{1}^{i})}{\widehat{\unicode[STIX]{x1D713}}(z_{0}^{j})}}\,q_{+}(z_{1}^{i}|z_{0}^{j})}.\end{eqnarray}$$

We will demonstrate in Section 2.3 that finding a twisting potential $\unicode[STIX]{x1D713}$ such that $\widehat{\unicode[STIX]{x1D70B}}_{1}=\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}$ , leading to importance weights $w^{i}=1$ in (2.40), is equivalent to solving the Schrödinger problem (2.58)–(2.61).

The associated smoothing distribution at time $t=0$ can be defined as follows. First introduce

(2.41) $$\begin{eqnarray}\unicode[STIX]{x1D713}(z_{1}):=\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})}{\unicode[STIX]{x1D70B}_{1}(z_{1})}=\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}.\end{eqnarray}$$

Next we set

(2.42) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D713}}(z_{0}):=\int q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D713}(z_{1})\,\text{d}z_{1}=\unicode[STIX]{x1D6FD}^{-1}\int q_{+}(z_{1}|z_{0})\,l(z_{1})\,\text{d}z_{1},\end{eqnarray}$$

and introduce $\widehat{\unicode[STIX]{x1D70B}}_{0}:=\unicode[STIX]{x1D70B}_{0}\,\widehat{\unicode[STIX]{x1D713}}$ , that is,

(2.43) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D70B}}_{0}(z_{0}) & = & \displaystyle \displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})\,\unicode[STIX]{x1D6FF}(z_{0}-z_{0}^{i})\nonumber\\ \displaystyle & = & \displaystyle \displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D6FE}^{i}\,\unicode[STIX]{x1D6FF}(z_{0}-z_{0}^{i})\end{eqnarray}$$

since $\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})=\unicode[STIX]{x1D6FE}^{i}$ with $\unicode[STIX]{x1D6FE}^{i}$ defined by (2.9).

Lemma 2.14. The smoothing PDFs $\widehat{\unicode[STIX]{x1D70B}}_{0}$ and $\widehat{\unicode[STIX]{x1D70B}}_{1}$ satisfy

(2.44) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{0}(z_{0})=\int q_{-}(z_{0}|z_{1})\,\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})\,\text{d}z_{1}\end{eqnarray}$$

with the backward transition kernel defined by (2.15). Furthermore,

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})=\int \widehat{q}_{+}(z_{1}|z_{0})\,\widehat{\unicode[STIX]{x1D70B}}_{0}(z_{0})\,\text{d}z_{0}\end{eqnarray}$$

with twisted forward transition kernels

(2.45) $$\begin{eqnarray}\widehat{q}_{+}(z_{1}|z_{0}^{i}):=\unicode[STIX]{x1D713}(z_{1})\,q_{+}(z_{1}|z_{0}^{i})\,\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})^{-1}=\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}\,\unicode[STIX]{x1D6FE}^{i}}q_{+}(z_{1}|z_{0}^{i})\end{eqnarray}$$

and $\unicode[STIX]{x1D6FE}^{i}$ , $i=1,\ldots ,M$ , defined by (2.9).

Proof. We note that

$$\begin{eqnarray}q_{-}(z_{0}|z_{1})\,\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})=\displaystyle \frac{\unicode[STIX]{x1D70B}_{0}(z_{0})}{\unicode[STIX]{x1D70B}_{1}(z_{1})}q_{+}(z_{1}|z_{0})\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})=\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}_{0}(z_{0}),\end{eqnarray}$$

which implies the first equation. The second equation follows from $\widehat{\unicode[STIX]{x1D70B}}_{0}=\widehat{\unicode[STIX]{x1D713}}\,\unicode[STIX]{x1D70B}_{0}$ and

$$\begin{eqnarray}\int \widehat{q}_{+}(z_{1}|z_{0})\,\widehat{\unicode[STIX]{x1D70B}}_{0}(z_{0})\,\text{d}z_{0}=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}\,q_{+}(z_{1}|z_{0}^{i}).\end{eqnarray}$$

In other words, we have defined a twisted forward transition kernel of the form (1.2).◻

Seen from a more abstract perspective, we have provided an alternative formulation of the joint smoothing distribution

(2.46) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}(z_{0},z_{1}):=\displaystyle \frac{l(z_{1})\,q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}_{0}(z_{0})}{\unicode[STIX]{x1D6FD}}\end{eqnarray}$$

in the form of

(2.47) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D70B}}(z_{0},z_{1}) & = & \displaystyle \displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}\displaystyle \frac{\unicode[STIX]{x1D713}(z_{1})}{\unicode[STIX]{x1D713}(z_{1})}q_{+}(z_{1}|z_{0})\displaystyle \frac{\widehat{\unicode[STIX]{x1D713}}(z_{0})}{\widehat{\unicode[STIX]{x1D713}}(z_{0})}\unicode[STIX]{x1D70B}_{0}(z_{0})\nonumber\\ \displaystyle & = & \displaystyle \widehat{q}_{+}(z_{1}|z_{0})\,\widehat{\unicode[STIX]{x1D70B}}_{0}(z_{0})\end{eqnarray}$$

because of (2.41). Note that the marginal distributions of $\widehat{\unicode[STIX]{x1D70B}}$ are provided by $\widehat{\unicode[STIX]{x1D70B}}_{0}$ and $\widehat{\unicode[STIX]{x1D70B}}_{1}$ , respectively.

One can exploit these formulations computationally as follows. If one has generated $M$ equally weighted particles $\widehat{z}_{0}^{j}$ from the smoothing distribution (2.43) at time $t=0$ via resampling with replacement, then one can obtain equally weighted samples $\widehat{z}_{1}^{j}$ from the filtering distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ using the modified transition kernels (2.45). This is the idea behind the optimal proposal particle filter (Doucet et al. Reference Doucet, de Freitas and Gordon2001, Arulampalam et al. Reference Arulampalam, Maskell, Gordon and Clapp2002, Fearnhead and Künsch Reference Fearnhead and Künsch2018) and provides an implementation of scenario (B) as introduced in Definition 2.4.

Lemma 2.16. Let $\unicode[STIX]{x1D713}(z)>0$ be a twisting potential such that

$$\begin{eqnarray}l^{\unicode[STIX]{x1D713}}(z_{1}):=\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}\,\unicode[STIX]{x1D713}(z_{1})}\unicode[STIX]{x1D70B}_{0}[\widehat{\unicode[STIX]{x1D713}}]\end{eqnarray}$$

is well-defined with $\widehat{\unicode[STIX]{x1D713}}$ given by (1.3). Then the smoothing PDF (2.46) can be represented as

(2.48) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}(z_{0},z_{1})=l^{\unicode[STIX]{x1D713}}(z_{1})\,q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}(z_{0}),\end{eqnarray}$$

where the modified forward transition kernel $q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0})$ is defined by (1.2) and the modified initial PDF by

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}(z_{0}):=\displaystyle \frac{\widehat{\unicode[STIX]{x1D713}}(z_{0})\,\unicode[STIX]{x1D70B}_{0}(z_{0})}{\unicode[STIX]{x1D70B}_{0}[\widehat{\unicode[STIX]{x1D713}}]}.\end{eqnarray}$$

2.2.1 Gaussian model errors (cont.)

We return to the discrete-time process (2.20) and assume a Gaussian measurement error leading to a Gaussian likelihood

$$\begin{eqnarray}l(z_{1})\propto \exp \biggl(-\displaystyle \frac{1}{2}(Hz_{1}-y_{1})^{T}R^{-1}(Hz_{1}-y_{1})\biggr).\end{eqnarray}$$

We set $\unicode[STIX]{x1D713}_{1}=l/\unicode[STIX]{x1D6FD}$ in order to derive the optimal forward kernel for the associated smoothing/filtering problem. Following the discussion from Section 2.1.1, this leads to the modified transition kernels

$$\begin{eqnarray}\widehat{q}_{+}(z_{1}|z_{0}^{i}):=\text{n}(z_{1};\bar{z}_{1}^{i},\unicode[STIX]{x1D6FE}\bar{B})\end{eqnarray}$$

with $\bar{B}$ and $K$ defined by (2.23) and (2.22), respectively, and

$$\begin{eqnarray}\bar{z}_{1}^{i}:=\unicode[STIX]{x1D6F9}(z_{0}^{i})-K(H\unicode[STIX]{x1D6F9}(z_{0}^{i})-y_{1}).\end{eqnarray}$$

The smoothing distribution $\widehat{\unicode[STIX]{x1D70B}}_{0}$ is given by

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{0}(z_{0})=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D6FE}^{i}\,\unicode[STIX]{x1D6FF}(z-z_{0}^{i})\end{eqnarray}$$

with coefficients

$$\begin{eqnarray}\unicode[STIX]{x1D6FE}^{i}\propto \exp \biggl(-\displaystyle \frac{1}{2}(H\unicode[STIX]{x1D6F9}(z_{0}^{i})-y_{1})^{\text{T}}(R+\unicode[STIX]{x1D6FE}HBH^{\text{T}})^{-1}(H\unicode[STIX]{x1D6F9}(z_{0}^{i})-y_{1})\biggr).\end{eqnarray}$$

It is easily checked that, indeed,

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})=\int q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0})\,\widehat{\unicode[STIX]{x1D70B}}_{0}(z_{0})\,\text{d}z_{0}.\end{eqnarray}$$

The results from this subsection have been used in simplified form in Example 2.5 in order to compute (2.14). We also note that a non-optimal, i.e.  $\unicode[STIX]{x1D713}(z_{1})\not =l(z_{1})/\unicode[STIX]{x1D6FD}$ , but Gaussian choice for $\unicode[STIX]{x1D713}$ leads to a Gaussian $l^{\unicode[STIX]{x1D713}}$ and the transition kernels $q_{+}^{\unicode[STIX]{x1D713}}(z_{1}|z_{0}^{i})$ in (2.48) remain Gaussian as well. This is in contrast to the Schrödinger problem, which we discuss in the following Section 2.3 and which leads to forward transition kernels of the form (2.66) below.

2.2.2 SDE models (cont.)

The likelihood $l(z_{1})$ introduces a change of measure over path space $z_{[0,1]}\in {\mathcal{C}}$ from the forecast measure $\mathbb{Q}$ with marginals $\unicode[STIX]{x1D70B}_{t}$ to the smoothing measure $\widehat{\mathbb{P}}$ via the Radon–Nikodym derivative

(2.49) $$\begin{eqnarray}\displaystyle \frac{\text{d}\widehat{\mathbb{P}}}{\text{d}\mathbb{Q}}_{|z_{[0,1]}}=\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}.\end{eqnarray}$$

We let $\widehat{\unicode[STIX]{x1D70B}}_{t}$ denote the marginal distributions of the smoothing measure $\widehat{\mathbb{P}}$ at time  $t$ .

Lemma 2.18. Let $\unicode[STIX]{x1D713}_{t}$ denote the solution of the backward Kolmogorov equation (2.31) with final condition $\unicode[STIX]{x1D713}_{1}(z)=l(z)/\unicode[STIX]{x1D6FD}$ at $t=1$ . Then the controlled forward SDE

(2.50) $$\begin{eqnarray}\text{d}Z_{t}^{+}=(f_{t}(Z_{t}^{+})+\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D713}_{t}(Z_{t}^{+}))\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{+},\end{eqnarray}$$

with $Z_{0}^{+}\sim \widehat{\unicode[STIX]{x1D70B}}_{0}$ at time $t=0$ , implies $Z_{1}^{+}\sim \widehat{\unicode[STIX]{x1D70B}}_{1}$ at final time $t=1$ .

Proof. The lemma is a consequence of Lemma 2.11 and definition (2.45) of the smoothing kernel $\widehat{q}_{+}(z_{1}|z_{0})$ with $\unicode[STIX]{x1D713}(z_{1})=\unicode[STIX]{x1D713}_{1}(z_{1})$ and $\widehat{\unicode[STIX]{x1D713}}(z_{0})=\unicode[STIX]{x1D713}_{0}(z_{0})$ .◻

The SDE (2.50) is obviously a special case of (2.30) with control law $u_{t}$ given by (2.32). Note, however, that the initial distributions for (2.30) and (2.50) are different. We will reconcile this fact in the following subsection by considering the associated Schrödinger problem (Föllmer and Gantert Reference Föllmer and Gantert1997, Leonard Reference Leonard2014, Chen et al. Reference Chen, Georgiou and Pavon2014).

Lemma 2.19. The solution $\unicode[STIX]{x1D713}_{t}$ of the backward Kolmogorov equation (2.31) with final condition $\unicode[STIX]{x1D713}_{1}(z)=l(z)/\unicode[STIX]{x1D6FD}$ at $t=1$ satisfies

(2.51) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{t}(z)=\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{t}(z)}{\unicode[STIX]{x1D70B}_{t}(z)}\end{eqnarray}$$

and the PDFs $\widehat{\unicode[STIX]{x1D70B}}_{t}$ coincide with the marginal PDFs of the backward SDE (2.28) with final condition $Z_{1}^{-}\sim \widehat{\unicode[STIX]{x1D70B}}_{1}$ .

Proof. We first note that (2.51) holds at final time $t=1$ . Furthermore, equation (2.51) implies

$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\widehat{\unicode[STIX]{x1D70B}}_{t}=\unicode[STIX]{x1D70B}_{t}\,\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}_{t}-\unicode[STIX]{x1D713}_{t}\,\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70B}_{t}.\end{eqnarray}$$

Since $\unicode[STIX]{x1D70B}_{t}$ satisfies the Fokker–Planck equation (2.26) and $\unicode[STIX]{x1D713}_{t}$ the backward Kolmogorov equation (2.31), it follows that

(2.52) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\widehat{\unicode[STIX]{x1D70B}}_{t}=-\unicode[STIX]{x1D6FB}_{z}\cdot (\widehat{\unicode[STIX]{x1D70B}}_{t}(f-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D70B}_{t}))-\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}\widehat{\unicode[STIX]{x1D70B}}_{t},\end{eqnarray}$$

which corresponds to the Fokker–Planck equation (2.27) of the backward SDE (2.28) with final condition $Z_{1}^{-}\sim \widehat{\unicode[STIX]{x1D70B}}_{1}$ and marginal PDFs denoted by $\widehat{\unicode[STIX]{x1D70B}}_{t}$ instead of $\unicode[STIX]{x1D70B}_{t}$ .◻

Note that (2.52) is equivalent to

$$\begin{eqnarray}\unicode[STIX]{x2202}\widehat{\unicode[STIX]{x1D70B}}_{t}=-\unicode[STIX]{x1D6FB}_{z}\cdot (\widehat{\unicode[STIX]{x1D70B}}_{t}(f-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D70B}_{t}+\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \widehat{\unicode[STIX]{x1D70B}}_{t}))+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}\widehat{\unicode[STIX]{x1D70B}}_{t},\end{eqnarray}$$

which in turn is equivalent to the Fokker–Planck equation of the forward smoothing SDE (2.50) since $\unicode[STIX]{x1D713}_{t}=\widehat{\unicode[STIX]{x1D70B}}_{t}/\unicode[STIX]{x1D70B}_{t}$ .

We conclude from the previous lemma that one can either solve the backward Kolmogorov equation (2.31) with $\unicode[STIX]{x1D713}_{1}(z)=l(z)/\unicode[STIX]{x1D6FD}$ or the backward SDE (2.28) with $Z_{1}^{-}\sim \widehat{\unicode[STIX]{x1D70B}}_{1}=l\,\unicode[STIX]{x1D70B}_{1}/\unicode[STIX]{x1D6FD}$ in order to derive the desired control law $u_{t}(z)=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D713}_{t}$ in (2.50).

Remark 2.20. The notion of a backward SDE used throughout this paper is different from the notion of a backward SDE in the sense of Carmona (Reference Carmona2016), for example. More specifically, Itô’s formula

$$\begin{eqnarray}\text{d}\unicode[STIX]{x1D713}_{t}=\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}_{t}\,\,\text{d}t+\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D713}_{t}\cdot \,\text{d}Z_{t}^{+}+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}\unicode[STIX]{x1D713}_{t}\,\text{d}t\end{eqnarray}$$

and the fact that $\unicode[STIX]{x1D713}_{t}$ satisfies the backward Kolmogorov equation (2.31) imply that

(2.53) $$\begin{eqnarray}\text{d}\unicode[STIX]{x1D713}_{t}=\unicode[STIX]{x1D6FE}^{1/2}\,\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D713}_{t}\cdot \,\text{d}W_{t}^{+}\end{eqnarray}$$

along solutions of the forward SDE (2.24). In other words, the quantities $\unicode[STIX]{x1D713}_{t}$ are materially advected along solutions $Z_{t}^{+}$ of the forward SDEs in expectation or, in the language of stochastic analysis, $\unicode[STIX]{x1D713}_{t}$ is a martingale. Hence, by the martingale representation theorem, we can reformulate the problem of determining $\unicode[STIX]{x1D713}_{t}$ as follows. Find the solution $(Y_{t},V_{t})$ of the backward SDE

(2.54) $$\begin{eqnarray}\text{d}Y_{t}=V_{t}\cdot \,\text{d}W_{t}^{+}\end{eqnarray}$$

subject to the final condition $Y_{1}=l(Z_{1}^{+})/\unicode[STIX]{x1D6FD}$ at $t=1$ . Here (2.54) has to be understood as a backward SDE in the sense of Carmona (Reference Carmona2016), for example, where the solution $(Y_{t},V_{t})$ is adapted to the forward SDE (2.24), i.e. to the past $s\leq t$ , whereas the solution $Z_{t}^{-}$ of the backward SDE (2.28) is adapted to the future $s\geq t$ . The solution of (2.54) is given by $Y_{t}=\unicode[STIX]{x1D713}_{t}(Z_{t}^{+})$ and $V_{t}=\unicode[STIX]{x1D6FE}^{1/2}\,\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D713}_{t}(Z_{t}^{+})$ in agreement with (2.53) (compare Carmona Reference Carmona2016, page 42). See Appendix A.4 for the numerical treatment of (2.54).

A variational characterization of the smoothing path measure $\widehat{\mathbb{P}}$ is given by the Donsker–Varadhan principle

(2.55) $$\begin{eqnarray}\widehat{\mathbb{P}}=\arg \inf _{\mathbb{P}\ll \mathbb{Q}}\{-\mathbb{P}[\log (l)]+\text{KL}(\mathbb{P}||\mathbb{Q})\},\end{eqnarray}$$

that is, the distribution $\widehat{\mathbb{P}}$ is chosen such that the expected loss, $-\mathbb{P}[\log (l)]$ , is minimized subject to the penalty introduced by the Kullback–Leibler divergence with respect to the original path measure $\mathbb{Q}$ . Note that

(2.56) $$\begin{eqnarray}\inf _{\mathbb{P}\ll \mathbb{Q}}\{\mathbb{P}[-\log (l)]+\text{KL}(\mathbb{P}||\mathbb{Q})\}=-\log \unicode[STIX]{x1D6FD}\end{eqnarray}$$

with $\unicode[STIX]{x1D6FD}=\mathbb{Q}[l]$ . The connection between smoothing for SDEs and the Donsker–Varadhan principle has, for example, been discussed by Mitter and Newton (Reference Mitter and Newton2003). See also Hartmann et al. (Reference Hartmann, Richter, Schütte and Zhang2017) for an in-depth discussion of variational formulations and their numerical implementation in the context of rare event simulations for which it is generally assumed that $\unicode[STIX]{x1D70B}_{0}(z)=\unicode[STIX]{x1D6FF}(z-z_{0})$ in (2.3), that is, the ensemble size is $M=1$ when viewed within the context of this paper.

Remark 2.21. One can choose $\unicode[STIX]{x1D713}_{t}$ differently from the choice made in (2.51) by changing the final condition for the backward Kolmogorov equation (2.31) to any suitable $\unicode[STIX]{x1D713}_{1}$ . As already discussed for twisted discrete-time smoothing, such modifications give rise to alternative representations of the smoothing distribution $\widehat{\mathbb{P}}$ in terms of modified forward SDEs, likelihood functions and initial distributions. See Kappen and Ruiz (Reference Kappen and Ruiz2016) and Ruiz and Kappen (Reference Ruiz and Kappen2017) for an application of these ideas to importance sampling in the context of partially observed diffusion processes. More specifically, let $u_{t}$ denote the associated control law (2.32) for the forward SDE (2.30) with given initial distribution $Z_{0}^{+}\sim q_{0}$ . Then

$$\begin{eqnarray}\displaystyle \frac{\text{d}\widehat{\mathbb{P}}}{\text{d}\mathbb{Q}}_{|z_{[0,1]}^{u}}=\displaystyle \frac{\text{d}\widehat{\mathbb{P}}}{\text{d}\mathbb{Q}^{u}}_{|z_{[0,1]}^{u}}\,\displaystyle \frac{\text{d}\mathbb{Q}^{u}}{\text{d}\mathbb{Q}}_{|z_{[0,1]}^{u}},\end{eqnarray}$$

which, using (2.33) and (2.49), implies the modified Radon–Nikodym derivative

(2.57) $$\begin{eqnarray}\displaystyle \frac{\text{d}\widehat{\mathbb{P}}}{\text{d}\mathbb{Q}^{u}}_{|z_{[0,1]}^{u}}=\displaystyle \frac{l(z_{1}^{u})}{\unicode[STIX]{x1D6FD}}\displaystyle \frac{\unicode[STIX]{x1D70B}_{0}(z_{0}^{u})}{q_{0}(z_{0}^{u})}\exp \biggl(-\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}}\int _{0}^{1}(\Vert u_{t}\Vert ^{2}\,\text{d}t+2\unicode[STIX]{x1D6FE}^{1/2}u_{t}\cdot \,\text{d}W_{t}^{+})\biggr).\end{eqnarray}$$

Recall from Lemma 2.18 that the control law (2.32) with $\unicode[STIX]{x1D713}_{t}$ defined by (2.51) together with $q_{0}=\widehat{\unicode[STIX]{x1D70B}}_{0}$ leads to $\widehat{\mathbb{P}}=\mathbb{Q}^{u}$ .

2.3 Schrödinger problem

In this subsection we show that scenario (C) from Definition 2.4 leads to a certain boundary value problem first considered by Schrödinger (Reference Schrödinger1931). More specifically, we state the so-called Schrödinger problem and show how it is linked to data assimilation scenario (C).

In order to introduce the Schrödinger problem, we return to the twisting potential approach as utilized in Section 2.2, with two important modifications. These modifications are, first, that the twisting potential $\unicode[STIX]{x1D713}$ is determined implicitly and, second, that the modified transition kernel $q_{+}^{\unicode[STIX]{x1D713}}$ is applied to $\unicode[STIX]{x1D70B}_{0}$ instead of the tilted initial density $\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}$ as in (2.48). More specifically, we have the following.

Definition 2.22. We seek the pair of functions $\widehat{\unicode[STIX]{x1D713}}(z_{0})$ and $\unicode[STIX]{x1D713}(z_{1})$ which solve the boundary value problem

(2.58) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70B}_{0}(z_{0}) & = & \displaystyle \unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}(z_{0})\,\widehat{\unicode[STIX]{x1D713}}(z_{0}),\end{eqnarray}$$

(2.59) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1}) & = & \displaystyle \unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1})\,\unicode[STIX]{x1D713}(z_{1}),\end{eqnarray}$$

(2.60) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1}) & = & \displaystyle \int q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}(z_{0})\,\text{d}z_{0},\end{eqnarray}$$

(2.61) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D713}}(z_{0}) & = & \displaystyle \int q_{+}(z_{1}|z_{0})\,\unicode[STIX]{x1D713}(z_{1})\,\text{d}z_{1},\end{eqnarray}$$

for given marginal (filtering) distributions $\unicode[STIX]{x1D70B}_{0}$ and $\widehat{\unicode[STIX]{x1D70B}}_{1}$ at $t=0$ and $t=1$ , respectively. The required modified PDFs $\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}$ and $\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}$ are defined by (2.58) and (2.59), respectively. The solution $(\widehat{\unicode[STIX]{x1D713}},\unicode[STIX]{x1D713})$ of the Schrödinger system (2.58)–(2.61) leads to the modified transition kernel

(2.62) $$\begin{eqnarray}q_{+}^{\ast }(z_{1}|z_{0}):=\unicode[STIX]{x1D713}(z_{1})\,q_{+}(z_{1}|z_{0})\,\widehat{\unicode[STIX]{x1D713}}(z_{0})^{-1},\end{eqnarray}$$

which satisfies

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})=\int q_{+}^{\ast }(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}_{0}(z_{0})\,\text{d}z_{0}\end{eqnarray}$$

by construction.

The modified transition kernel $q_{+}^{\ast }(z_{1}|z_{0})$ couples the two marginal distributions $\unicode[STIX]{x1D70B}_{0}$ and $\widehat{\unicode[STIX]{x1D70B}}_{1}$ with the twisting potential $\unicode[STIX]{x1D713}$ implicitly defined. In other words, $q_{+}^{\ast }$ provides the transition kernel for going from the initial distribution (2.3) at time $t_{0}$ to the filtering distribution at time $t_{1}$ without the need for any reweighting, i.e. the desired transition kernel for scenario (C). See Leonard (Reference Leonard2014) and Chen et al. (Reference Chen, Georgiou and Pavon2014) for more mathematical details on the Schrödinger problem.

Remark 2.23. Let us compare the Schrödinger system to the twisting potential approach (2.48) for the smoothing problem from Section 2.2 in some more detail. First, note that the twisting potential approach to smoothing replaces (2.58) with

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}(z_{0})=\unicode[STIX]{x1D70B}_{0}(z_{0})\,\widehat{\unicode[STIX]{x1D713}}(z_{0})\end{eqnarray}$$

and (2.59) with

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1})=\unicode[STIX]{x1D70B}_{1}(z_{1})\,\unicode[STIX]{x1D713}(z_{1}),\end{eqnarray}$$

where $\unicode[STIX]{x1D713}$ is a given twisting potential normalized such that $\unicode[STIX]{x1D70B}_{1}[\unicode[STIX]{x1D713}]=1$ . The associated $\widehat{\unicode[STIX]{x1D713}}$ is determined by (2.61) as in the twisting approach. In both cases, the modified transition kernel is given by (2.62). Finally, (2.60) is replaced by the prediction step (2.4).

In order to solve the Schrödinger system for our given initial distribution (2.3) and the associated filter distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ , we make the ansatz

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}(z_{0})=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D6FC}^{i}\,\unicode[STIX]{x1D6FF}(z_{0}-z_{0}^{i}),\quad \mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D6FC}^{i}=M.\end{eqnarray}$$

This ansatz together with (2.58)–(2.61) immediately implies

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})=\displaystyle \frac{1}{\unicode[STIX]{x1D6FC}^{i}},\quad \unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1})=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D6FC}^{i}\,q_{+}(z_{1}|z_{0}^{i}),\end{eqnarray}$$

as well as

(2.63) $$\begin{eqnarray}\unicode[STIX]{x1D713}(z_{1})=\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})}{\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1})}=\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}\displaystyle \frac{{\textstyle \frac{1}{M}}\mathop{\sum }_{j=1}^{M}q_{+}(z_{1}|z_{0}^{j})}{{\textstyle \frac{1}{M}}\mathop{\sum }_{j=1}^{M}\unicode[STIX]{x1D6FC}^{j}\,q_{+}(z_{1}|z_{0}^{j})}.\end{eqnarray}$$

Hence we arrive at the equations

(2.64) $$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D713}}(z_{0}^{i}) & = & \displaystyle \displaystyle \frac{1}{\unicode[STIX]{x1D6FC}^{i}}\end{eqnarray}$$

(2.65) $$\begin{eqnarray}\displaystyle & = & \displaystyle \int \unicode[STIX]{x1D713}(z_{1})\,q_{+}(z_{1}|z_{0}^{i})\,\text{d}z_{1}\nonumber\\ \displaystyle & = & \displaystyle \int \displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}\displaystyle \frac{\mathop{\sum }_{j=1}^{M}q_{+}(z_{1}|z_{0}^{j})}{\mathop{\sum }_{j=1}^{M}\unicode[STIX]{x1D6FC}^{j}\,q_{+}(z_{1}|z_{0}^{j})}\,q_{+}(z_{1}|z_{0}^{i})\,\text{d}z_{1}\end{eqnarray}$$

for $i=1,\ldots ,M$ . These $M$ equations have to be solved for the $M$ unknown coefficients $\{\unicode[STIX]{x1D6FC}^{i}\}$ . In other words, the Schrödinger problem becomes finite-dimensional in the context of this paper. More specifically, we have the following result.

Lemma 2.24. The forward Schrödinger transition kernel (2.62) is given by

(2.66) $$\begin{eqnarray}\displaystyle q_{+}^{\ast }(z_{1}|z_{0}^{i}) & = & \displaystyle \displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})}{\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}(z_{1})}q_{+}(z_{1}|z_{0}^{i})\,\unicode[STIX]{x1D6FC}^{i}\nonumber\\ \displaystyle & = & \displaystyle \displaystyle \frac{\unicode[STIX]{x1D6FC}^{i}\mathop{\sum }_{j=1}^{M}\,q_{+}(z_{1}|z_{0}^{j})}{\mathop{\sum }_{j=1}^{M}\unicode[STIX]{x1D6FC}^{j}\,q_{+}(z_{1}|z_{0}^{j})}\,\displaystyle \frac{l(z_{1})}{\unicode[STIX]{x1D6FD}}\,q_{+}(z_{1}|z_{0}^{i}),\end{eqnarray}$$

for each particle $z_{0}^{i}$ with the coefficients $\unicode[STIX]{x1D6FC}^{j}$ , $j=1,\ldots ,M$ , defined by (2.64)–(2.65).

Proof. Because of (2.64)–(2.65), the forward transition kernels (2.66) satisfy

(2.67) $$\begin{eqnarray}\int q_{+}^{\ast }(z_{1}|z_{0}^{i})\,\text{d}z_{1}=1\end{eqnarray}$$

for all $i=1,\ldots ,M$ and

(2.68) $$\begin{eqnarray}\displaystyle \int q_{+}^{\ast }(z_{1}|z_{0})\,\unicode[STIX]{x1D70B}_{0}(z_{0})\,\text{d}z & = & \displaystyle \displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}q_{+}^{\ast }(z_{1}|z_{0}^{i})\nonumber\\ \displaystyle & = & \displaystyle \displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D6FC}^{i}q_{+}(z_{1}|z_{0}^{i})\,\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})}{{\textstyle \frac{1}{M}}\mathop{\sum }_{j=1}^{M}\unicode[STIX]{x1D6FC}^{j}\,q_{+}(z_{1}|z_{0}^{j})}\nonumber\\ \displaystyle & = & \displaystyle \widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1}),\end{eqnarray}$$

as desired. ◻

Numerical implementations will be discussed in Section 3. Note that knowledge of the normalizing constant $\unicode[STIX]{x1D6FD}$ is not required a priori for solving (2.64)–(2.65) since it appears as a common scaling factor.

We note that the coefficients $\{\unicode[STIX]{x1D6FC}^{j}\}$ together with the associated potential $\unicode[STIX]{x1D713}$ from the Schrödinger system provide the optimally twisted prediction kernel (1.2) with respect to the filtering distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ , that is, we set $\widehat{\unicode[STIX]{x1D713}}(z_{0}^{i})=1/\unicode[STIX]{x1D6FC}^{i}$ in (2.16) and define the potential $\unicode[STIX]{x1D713}$ by (2.63).

Lemma 2.25. The Schrödinger transition kernel (2.62) satisfies the following constrained variational principle. Consider the joint PDFs given by $\unicode[STIX]{x1D70B}(z_{0},z_{1}):=q_{+}(z_{1}|z_{0})\unicode[STIX]{x1D70B}_{0}(z_{0})$ and $\unicode[STIX]{x1D70B}^{\ast }(z_{0},z_{1}):=q_{+}^{\ast }(z_{1}|z_{0})\unicode[STIX]{x1D70B}_{0}(z_{0})$ . Then

(2.69) $$\begin{eqnarray}\unicode[STIX]{x1D70B}^{\ast }=\arg \inf _{\widetilde{\unicode[STIX]{x1D70B}}\in \unicode[STIX]{x1D6F1}_{\text{S}}}\text{KL}(\widetilde{\unicode[STIX]{x1D70B}}||\unicode[STIX]{x1D70B}).\end{eqnarray}$$

Here a joint PDF $\widetilde{\unicode[STIX]{x1D70B}}(z_{0},z_{1})$ is an element of $\unicode[STIX]{x1D6F1}_{\text{S}}$ if

$$\begin{eqnarray}\int \widetilde{\unicode[STIX]{x1D70B}}(z_{0},z_{1})\,\text{d}z_{1}=\unicode[STIX]{x1D70B}_{0}(z_{0}),\quad \int \widetilde{\unicode[STIX]{x1D70B}}(z_{0},z_{1})\,\text{d}z_{0}=\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1}).\end{eqnarray}$$

The constrained variational formulation (2.69) of Schrödinger’s problem should be compared to the unconstrained Donsker–Varadhan variational principle

(2.70) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}=\arg \inf \{-\widetilde{\unicode[STIX]{x1D70B}}[\log (l)]+\text{KL}(\widetilde{\unicode[STIX]{x1D70B}}||\unicode[STIX]{x1D70B})\}\end{eqnarray}$$

for the associated smoothing problem. See Remark 2.27 below.

Remark 2.26. The Schrödinger problem is closely linked to optimal transportation (Cuturi Reference Cuturi and Burges2013, Leonard Reference Leonard2014, Chen et al. Reference Chen, Georgiou and Pavon2014). For example, consider the Gaussian transition kernel (2.21) with $\unicode[STIX]{x1D6F9}(z)=z$ and $B=I$ . Then the solution (2.66) to the associated Schrödinger problem of coupling $\unicode[STIX]{x1D70B}_{0}$ and $\widehat{\unicode[STIX]{x1D70B}}_{1}$ reduces to the solution $\unicode[STIX]{x1D70B}^{\ast }$ of the associated optimal transport problem

$$\begin{eqnarray}\unicode[STIX]{x1D70B}^{\ast }=\arg \inf _{\widetilde{\unicode[STIX]{x1D70B}}\in \unicode[STIX]{x1D6F1}_{\text{S}}}\int \int \Vert z_{0}-z_{1}\Vert ^{2}\,\widetilde{\unicode[STIX]{x1D70B}}(z_{0},z_{1})\,\text{d}z_{0}\,\text{d}z_{1}\end{eqnarray}$$

in the limit $\unicode[STIX]{x1D6FE}\rightarrow 0$ .

2.3.1 SDE models (cont.)

At the SDE level, Schrödinger’s problem amounts to continuously bridging the given initial PDF $\unicode[STIX]{x1D70B}_{0}$ with the PDF $\widehat{\unicode[STIX]{x1D70B}}_{1}$ at final time using an appropriate modification of the stochastic process $Z_{[0,1]}^{+}\sim \mathbb{Q}$ defined by the forward SDE (2.24) with initial distribution $\unicode[STIX]{x1D70B}_{0}$ at $t=0$ . The desired modified stochastic process $\mathbb{P}^{\ast }$ is defined as the minimizer of

$$\begin{eqnarray}{\mathcal{L}}(\widetilde{\mathbb{P}}):=\text{KL}(\widetilde{\mathbb{P}}||\mathbb{Q})\end{eqnarray}$$

subject to the constraint that the marginal distributions $\widetilde{\unicode[STIX]{x1D70B}}_{t}$ of $\widetilde{\mathbb{P}}$ at time $t=0$ and $t=1$ satisfy $\unicode[STIX]{x1D70B}_{0}$ and $\widehat{\unicode[STIX]{x1D70B}}_{1}$ , respectively (Föllmer and Gantert Reference Föllmer and Gantert1997, Leonard Reference Leonard2014, Chen et al. Reference Chen, Georgiou and Pavon2014).

Remark 2.27. We note that the Donsker–Varadhan variational principle (2.55), characterizing the smoothing path measure $\widehat{\mathbb{P}}$ , can be replaced by

$$\begin{eqnarray}\mathbb{P}^{\ast }=\arg \inf _{\widetilde{\mathbb{ P}}\in \unicode[STIX]{x1D6F1}}\{-\widetilde{\unicode[STIX]{x1D70B}}_{1}[\log (l)]+\text{KL}(\widetilde{\mathbb{P}}||\mathbb{Q})\}\end{eqnarray}$$

with

$$\begin{eqnarray}\unicode[STIX]{x1D6F1}=\{\widetilde{\mathbb{P}}\ll \mathbb{Q}:\,\widetilde{\unicode[STIX]{x1D70B}}_{1}=\widehat{\unicode[STIX]{x1D70B}}_{1},\,\widetilde{\unicode[STIX]{x1D70B}}_{0}=\unicode[STIX]{x1D70B}_{0}\}\end{eqnarray}$$

in the context of Schrödinger’s problem. The associated

$$\begin{eqnarray}-\log \unicode[STIX]{x1D6FD}^{\ast }:=\inf _{\widetilde{\mathbb{ P}}\in \unicode[STIX]{x1D6F1}}\{-\widetilde{\unicode[STIX]{x1D70B}}_{1}[\log (l)]+\text{KL}(\widetilde{\mathbb{P}}||\mathbb{Q})\}=-\widehat{\unicode[STIX]{x1D70B}}[\log (l)]+\text{KL}(\mathbb{P}^{\ast }||\mathbb{Q})\end{eqnarray}$$

can be viewed as a generalization of (2.56) and gives rise to a generalized evidence $\unicode[STIX]{x1D6FD}^{\ast }$ , which could be used for model comparison and parameter estimation.

The Schrödinger process $\mathbb{P}^{\ast }$ corresponds to a Markovian process across the whole time domain $[0,1]$ (Leonard Reference Leonard2014, Chen et al. Reference Chen, Georgiou and Pavon2014). More specifically, consider the controlled forward SDE (2.30) with initial conditions

$$\begin{eqnarray}Z_{0}^{+}\sim \unicode[STIX]{x1D70B}_{0}\end{eqnarray}$$

and a given control law $u_{t}$ for $t\in [0,1]$ . Let $\mathbb{P}^{u}$ denote the path measure associated to this process. Then, as discussed in detail by Dai Pra (Reference Dai Pra1991), one can find time-dependent potentials $\unicode[STIX]{x1D713}_{t}$ with associated control laws (2.32) such that the marginal of the associated path measure $\mathbb{P}^{u}$ at times $t=1$ satisfies

$$\begin{eqnarray}\unicode[STIX]{x1D70B}_{1}^{u}=\widehat{\unicode[STIX]{x1D70B}}_{1}\end{eqnarray}$$

and, more generally,

$$\begin{eqnarray}\mathbb{P}^{\ast }=\mathbb{P}^{u}.\end{eqnarray}$$

We summarize this result in the following lemma.

Lemma 2.28. The Schrödinger path measure $\mathbb{P}^{\ast }$ can be generated by a controlled SDE (2.30) with control law (2.32), where the desired potential $\unicode[STIX]{x1D713}_{t}$ can be obtained as follows. Let $(\widehat{\unicode[STIX]{x1D713}},\unicode[STIX]{x1D713})$ denote the solution of the associated Schrödinger system (2.58)–(2.61), where $q_{+}(z_{1}|z_{0})$ denotes the time-one forward transition kernel of (2.24). Then $\unicode[STIX]{x1D713}_{t}$ in (2.32) is the solution of the backward Kolmogorov equation (2.31) with prescribed $\unicode[STIX]{x1D713}_{1}=\unicode[STIX]{x1D713}$ at final time $t=1$ .

Remark 2.29. As already pointed out in the context of smoothing, the desired potential $\unicode[STIX]{x1D713}_{t}$ can also be obtained by solving an appropriate backward SDE. More specifically, given the solution $(\widehat{\unicode[STIX]{x1D713}},\unicode[STIX]{x1D713})$ and the implied PDF $\widetilde{\unicode[STIX]{x1D70B}}_{0}^{+}:=\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D70B}_{0}/\widehat{\unicode[STIX]{x1D713}}$ of the Schrödinger system (2.58)–(2.61), let $\widetilde{\unicode[STIX]{x1D70B}}_{t}^{+}$ , $t\geq 0$ , denote the marginals of the forward SDE (2.24) with $Z_{0}^{+}\sim \widetilde{\unicode[STIX]{x1D70B}}_{0}^{+}$ . Furthermore, consider the backward SDE (2.28) with drift term

(2.71) $$\begin{eqnarray}b_{t}(z)=f_{t}(z)-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \widetilde{\unicode[STIX]{x1D70B}}_{t}^{+}(z),\end{eqnarray}$$

and final time condition $Z_{1}^{-}\sim \widehat{\unicode[STIX]{x1D70B}}_{1}$ . Then the choice of $\widetilde{\unicode[STIX]{x1D70B}}_{0}^{+}$ ensures that $Z_{0}^{-}\sim \unicode[STIX]{x1D70B}_{0}$ . Furthermore the desired control in (2.30) is provided by

$$\begin{eqnarray}u_{t}=\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \displaystyle \frac{\widetilde{\unicode[STIX]{x1D70B}}_{t}^{-}}{\widetilde{\unicode[STIX]{x1D70B}}_{t}^{+}},\end{eqnarray}$$

where $\widetilde{\unicode[STIX]{x1D70B}}_{t}^{-}$ denotes the marginal distributions of the backward SDE (2.28) with drift term (2.71) and $\widetilde{\unicode[STIX]{x1D70B}}_{1}^{-}=\widehat{\unicode[STIX]{x1D70B}}_{1}$ . We will return to this reformulation of the Schrödinger problem in Section 3 when considering it as the limit of a sequence of smoothing problems.

Remark 2.30. The solution to the Schrödinger problem for linear SDEs and Gaussian marginal distributions has been discussed in detail by Chen, Georgiou and Pavon (Reference Chen, Georgiou and Pavon2016b ).

2.3.2 Discrete measures

We finally discuss the Schrödinger problem in the context of finite-state Markov chains in more detail. These results will be needed in the following sections on the numerical implementation of the Schrödinger approach to sequential data assimilation.

Let us therefore consider an example which will be closely related to the discussion in Section 3. We are given a bi-stochastic matrix $Q\in \mathbb{R}^{L\times M}$ with all entries satisfying $q_{lj}>0$ and two discrete probability measures represented by vectors $p_{1}\in \mathbb{R}^{L}$ and $p_{0}\in \mathbb{R}^{M}$ , respectively. Again we assume for simplicity that all entries in $p_{1}$ and $p_{0}$ are strictly positive. We introduce the set of all bi-stochastic $L\times M$ matrices with those discrete probability measures as marginals, that is,

(2.72) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{\text{s}}:=\{P\in \mathbb{R}^{L\times M}:\,P\geq 0,\,P^{\text{T}}\unicode[STIX]{x1D7D9}_{L}=p_{0},\,P\unicode[STIX]{x1D7D9}_{M}=p_{1}\}.\end{eqnarray}$$

Solving Schrödinger’s system (2.58)–(2.61) corresponds to finding two non-negative vectors $u\in \mathbb{R}^{L}$ and $v\in \mathbb{R}^{M}$ such that

$$\begin{eqnarray}P^{\ast }:=D(u)\,Q\,D(v)^{-1}\in \unicode[STIX]{x1D6F1}_{\text{s}}.\end{eqnarray}$$

In turns out that $P^{\ast }$ is uniquely determined and minimizes the Kullback–Leibler divergence between all $P\in \unicode[STIX]{x1D6F1}_{\text{s}}$ and the reference matrix $Q$ , that is,

(2.73) $$\begin{eqnarray}P^{\ast }=\arg \min _{P\in \unicode[STIX]{x1D6F1}_{\text{s}}}\text{KL}\,(P||Q).\end{eqnarray}$$

See Peyre and Cuturi (Reference Peyre and Cuturi2018) and the following remark for more details.

Remark 2.31. If one makes the ansatz

$$\begin{eqnarray}p_{lj}=\displaystyle \frac{u_{l}q_{lj}}{v_{j}},\end{eqnarray}$$

then the minimization problem (2.73) becomes equivalent to

$$\begin{eqnarray}P^{\ast }=\arg \min _{(u,v)>0}\mathop{\sum }_{l,j}p_{lj}(\log u_{l}-\log v_{j})\end{eqnarray}$$

subject to the additional constraints

$$\begin{eqnarray}P\,\unicode[STIX]{x1D7D9}_{M}=D(u)QD(v)^{-1}\,\unicode[STIX]{x1D7D9}_{M}=p_{1},\,P^{\text{T}}\,\unicode[STIX]{x1D7D9}_{L}=D(v)^{-1}Q^{\text{T}}D(u)\,\unicode[STIX]{x1D7D9}_{L}=p_{0}.\end{eqnarray}$$

Note that these constraint determine $u>0$ and $v>0$ up to a common scaling factor. Hence (2.73) can be reduced to finding $(u,v)>0$ such that

$$\begin{eqnarray}u^{\text{T}}\unicode[STIX]{x1D7D9}_{L}=1,\quad P\unicode[STIX]{x1D7D9}_{M}=p_{1},\quad P^{\text{T}}\unicode[STIX]{x1D7D9}_{L}=p_{0}.\end{eqnarray}$$

Hence we have shown that solving the Schrödinger system is equivalent to solving the minimization problem (2.73) for discrete measures. Thus

$$\begin{eqnarray}\min _{P\in \unicode[STIX]{x1D6F1}_{\text{s}}}\text{KL}\,(P||Q)=p_{1}^{\text{T}}\log u-p_{0}^{\text{T}}\log v.\end{eqnarray}$$

Lemma 2.32. The Sinkhorn iteration (Sinkhorn Reference Sinkhorn1967)

(2.74) $$\begin{eqnarray}\displaystyle u^{k+1} & := & \displaystyle D(P^{2k}\unicode[STIX]{x1D7D9}_{M})^{-1}\,p_{1},\end{eqnarray}$$

(2.75) $$\begin{eqnarray}\displaystyle P^{2k+1} & := & \displaystyle D(u^{k+1})\,P^{2k},\end{eqnarray}$$

(2.76) $$\begin{eqnarray}\displaystyle v^{k+1} & := & \displaystyle D(p_{0})^{-1}\,(P^{2k+1})^{\text{T}}\,\unicode[STIX]{x1D7D9}_{L},\end{eqnarray}$$

(2.77) $$\begin{eqnarray}\displaystyle P^{2k+2} & := & \displaystyle P^{2k+1}\,D(v^{k+1})^{-1},\end{eqnarray}$$

$k=0,1,\ldots ,$ with initial $P^{0}=Q\in \mathbb{R}^{L\times M}$ provides an algorithm for computing $P^{\ast }$ , that is,

(2.78) $$\begin{eqnarray}\lim _{k\rightarrow \infty }P^{k}=P^{\ast }.\end{eqnarray}$$

Proof. See, for example, Peyre and Cuturi (Reference Peyre and Cuturi2018) for a proof of (2.78), which is based on the contraction property of the iteration (2.74)–(2.77) with respect to the Hilbert metric on the projective cone of positive vectors.◻

It follows from (2.78) that

$$\begin{eqnarray}\lim _{k\rightarrow \infty }u^{k}=\unicode[STIX]{x1D7D9}_{L},\quad \lim _{k\rightarrow \infty }v^{k}=\unicode[STIX]{x1D7D9}_{M}.\end{eqnarray}$$

The essential idea of the Sinkhorn iteration is to enforce

$$\begin{eqnarray}P^{2k+1}\,\unicode[STIX]{x1D7D9}_{M}=p_{1},\quad (P^{2k})^{\text{T}}\,\unicode[STIX]{x1D7D9}_{L}=p_{0}\end{eqnarray}$$

at each iteration step and that the matrix $P^{\ast }$ satisfies both constraints simultaneously in the limit $k\rightarrow \infty$ . See Cuturi (Reference Cuturi and Burges2013) for a computationally efficient and robust implementation of the Sinkhorn iteration.

Remark 2.33. One can introduce a similar iteration for the Schrödinger system (2.58)–(2.61). For example, pick $\widehat{\unicode[STIX]{x1D713}}(z_{0})=1$ initially. Then (2.58) implies $\unicode[STIX]{x1D70B}_{0}^{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D70B}_{0}$ and (2.60) $\unicode[STIX]{x1D70B}_{1}^{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D70B}_{1}$ . Hence $\unicode[STIX]{x1D713}=l/\unicode[STIX]{x1D6FD}$ in the first iteration. The second iteration starts with $\widehat{\unicode[STIX]{x1D713}}$ determined by (2.61) with $\unicode[STIX]{x1D713}=l/\unicode[STIX]{x1D6FD}$ . We again cycle through (2.58), (2.60) and (2.59) in order to find the next approximation to $\unicode[STIX]{x1D713}$ . The third iteration takes now this $\unicode[STIX]{x1D713}$ and computes the associated $\widehat{\unicode[STIX]{x1D713}}$ from (2.61) etc. A numerical implementation of this procedure requires the approximation of two integrals which essentially leads back to a Sinkhorn type algorithm in the weights of an appropriate quadrature rule.

3.1 Prediction

Generating samples from the forecast distributions $q_{+}(\cdot \,|z_{0}^{i})$ is in most cases straightforward. The computational expenses can, however, vary dramatically, and this impacts on the choice of algorithms for sequential data assimilation. We demonstrate in this subsection how samples from the prediction PDF $\unicode[STIX]{x1D70B}_{1}$ can be used to construct an associated finite-state Markov chain that transforms $\unicode[STIX]{x1D70B}_{0}$ into an empirical approximation of $\unicode[STIX]{x1D70B}_{1}$ .

Definition 3.1. Let us assume that we have $L\geq M$ independent samples $z_{1}^{l}$ from the $M$ forecast distributions $q_{+}(\cdot \,|z_{0}^{j})$ , $j=1,\ldots ,M$ . We introduce the $L\times M$ matrix $Q$ with entries

(3.1) $$\begin{eqnarray}q_{lj}:=q_{+}(z_{1}^{l}|z_{0}^{j}).\end{eqnarray}$$

We now consider the associated bi-stochastic matrix $P^{\ast }\in \mathbb{R}^{L\times M}$ , as defined by (2.73), with the two probability vectors in (2.72) given by $p_{1}=\unicode[STIX]{x1D7D9}_{L}/L\in \mathbb{R}^{L}$ and $p_{0}=\unicode[STIX]{x1D7D9}_{M}/M\in \mathbb{R}^{M}$ , respectively. The finite-state Markov chain

(3.2) $$\begin{eqnarray}Q_{+}:=MP^{\ast }\end{eqnarray}$$

provides a sample-based approximation to the forward transition kernel $q_{+}(z_{1}|z_{0})$ .

More precisely, the $i$ th column of $Q_{+}$ provides an empirical approximation to $q_{+}(\cdot |z_{0}^{i})$ and

$$\begin{eqnarray}Q_{+}\,p_{0}=p_{1}=\displaystyle \frac{1}{L}\unicode[STIX]{x1D7D9}_{L},\end{eqnarray}$$

which is in agreement with the fact that the $z_{1}^{l}$ are equally weighted samples from the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ .

Remark 3.2. Because of the simple relation between a bi-stochastic matrix $P\in \unicode[STIX]{x1D6F1}_{\text{s}}$ with $p_{0}$ in (2.72) given by $p_{0}=\unicode[STIX]{x1D7D9}_{M}/M$ and its associated finite-state Markov chain $Q_{+}=MP$ , one can reformulate the minimization problem (2.73) in those cases directly in terms of Markov chains $Q_{+}\in \unicode[STIX]{x1D6F1}_{\text{M}}$ with the definition of $\unicode[STIX]{x1D6F1}_{\text{s}}$ adjusted to

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{\text{M}}:=\bigg\{Q\in \mathbb{R}^{L\times M}:\,Q\geq 0,\,Q^{\text{T}}\unicode[STIX]{x1D7D9}_{L}=\unicode[STIX]{x1D7D9}_{M},\,\displaystyle \frac{1}{M}Q\unicode[STIX]{x1D7D9}_{M}=p_{1}\bigg\}.\end{eqnarray}$$

Remark 3.3. The associated backward transition kernel $Q_{-}\in \mathbb{R}^{M\times L}$ satisfies

$$\begin{eqnarray}Q_{-}\,D(p_{1})=(Q_{+}\,D(p_{0}))^{\text{T}}\end{eqnarray}$$

and is hence given by

$$\begin{eqnarray}Q_{-}=(Q_{+}\,D(p_{0}))^{\text{T}}\,D(p_{1})^{-1}=\displaystyle \frac{L}{M}Q_{+}^{\text{T}}.\end{eqnarray}$$

Thus

$$\begin{eqnarray}Q_{-}\,p_{1}=D(p_{0})\,Q_{+}^{\text{T}}\,\unicode[STIX]{x1D7D9}_{L}=D(p_{0})\,\unicode[STIX]{x1D7D9}_{M}=p_{0},\end{eqnarray}$$

as desired.

Definition 3.4. We can extend the concept of twisting to discrete Markov chains such as (3.2). A twisting potential $\unicode[STIX]{x1D713}$ gives rise to a vector $u\in \mathbb{R}^{L}$ with normalized entries

$$\begin{eqnarray}u_{l}=\displaystyle \frac{\unicode[STIX]{x1D713}(z_{1}^{l})}{\mathop{\sum }_{k=1}^{L}\unicode[STIX]{x1D713}(z_{1}^{k})},\end{eqnarray}$$

$l=1,\ldots ,L$ . The twisted finite-state Markov kernel is now defined by

(3.4) $$\begin{eqnarray}Q_{+}^{\unicode[STIX]{x1D713}}:=D(u)\,Q_{+}\,D(v)^{-1},\quad v:=(D(u)\,Q_{+})^{\text{T}}\,\unicode[STIX]{x1D7D9}_{L}\in \mathbb{R}^{M},\end{eqnarray}$$

and thus $\unicode[STIX]{x1D7D9}_{L}^{\text{T}}\,Q_{+}^{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D7D9}_{M}^{\text{T}}$ , as required for a Markov kernel. The twisted forecast probability is given by

$$\begin{eqnarray}p_{1}^{\unicode[STIX]{x1D713}}:=Q_{+}^{\unicode[STIX]{x1D713}}\,p_{0}\end{eqnarray}$$

with $p_{0}=\unicode[STIX]{x1D7D9}_{M}/M$ . Furthermore, if we set $p_{0}=v$ then $p_{1}^{\unicode[STIX]{x1D713}}=u$ .

3.1.1 Gaussian model errors (cont.)

The proposal density is given by (2.21) and it is easy to produce $K>1$ samples from each of the $M$ proposals $q_{+}(\cdot \,|z_{0}^{j})$ . Hence we can make the total sample size $L=K\,M$ as large as desired. In order to produce $M$ samples, $\widetilde{z}_{1}^{j}$ , from a twisted finite-state Markov chain (3.4), we draw a single realization from each of the $M$ associated discrete random variables $\widetilde{Z}_{1}^{j}$ , $j=1,\ldots ,M$ , with probabilities

$$\begin{eqnarray}\mathbb{P}[\widetilde{Z}_{1}^{j}(\unicode[STIX]{x1D714})=z_{1}^{l}]=(Q_{+}^{\unicode[STIX]{x1D713}})_{lj}.\end{eqnarray}$$

We will provide more details when discussing the Schrödinger problem in the context of Gaussian model errors in Section 3.4.1.

3.1.2 SDE models (cont.)

The Euler–Maruyama method (Kloeden and Platen Reference Kloeden and Platen1992)

(3.5) $$\begin{eqnarray}Z_{n+1}^{+}=Z_{n}^{+}+f_{t_{n}}(Z_{n}^{+})\,\unicode[STIX]{x1D6E5}t+(\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6E5}t)^{1/2}\,\unicode[STIX]{x1D6EF}_{n},\quad \unicode[STIX]{x1D6EF}_{n}\sim \text{N}(0,I),\end{eqnarray}$$

$n=0,\ldots ,N-1$ , will be used for the numerical approximation of (2.24) with step-size $\unicode[STIX]{x1D6E5}t:=1/N$ , $t_{n}=n\unicode[STIX]{x1D6E5}t$ . In other words, we replace $Z_{t_{n}}^{+}$ with its numerical approximation $Z_{n}^{+}$ . A numerical approximation (realization) of the whole solution path $z_{[0,1]}$ will be denoted by $z_{0:N}=Z_{0:N}^{+}(\unicode[STIX]{x1D714})$ and can be computed recursively due to the Markov property of the Euler–Maruyama scheme. The marginal PDFs of $Z_{n}^{+}$ are denoted by $\unicode[STIX]{x1D70B}_{n}$ .

For any finite number of time-steps $N$ , we can define a joint PDF $\unicode[STIX]{x1D70B}_{0:N}$ on ${\mathcal{U}}_{N}=\mathbb{R}^{N_{z}\times (N+1)}$ via

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D70B}_{0:N}(z_{0:N})\propto \exp \biggl(-\displaystyle \frac{1}{2\unicode[STIX]{x1D6E5}t}\mathop{\sum }_{n=0}^{N-1}\Vert \unicode[STIX]{x1D702}_{n}\Vert ^{2}\biggr)\unicode[STIX]{x1D70B}_{0}(z_{0})\end{eqnarray}$$

with

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{n}:=\unicode[STIX]{x1D6FE}^{-1/2}(z_{n+1}-z_{n}-f_{t_{n}}(z_{n})\unicode[STIX]{x1D6E5}t)\end{eqnarray}$$

and $\unicode[STIX]{x1D702}_{n}=\unicode[STIX]{x1D6E5}t^{1/2}\unicode[STIX]{x1D6EF}_{n}(\unicode[STIX]{x1D714})$ . Note that the joint PDF $\unicode[STIX]{x1D70B}_{0:N}(z_{0:N})$ can also be expressed in terms of $z_{0}$ and $\unicode[STIX]{x1D702}_{0:N-1}$ .

The numerical approximation of SDEs provides an example for which the increase in computational cost for producing $L>M$ samples from the PDF $\unicode[STIX]{x1D70B}_{0:N}$ versus $L=M$ is non-trivial, in general.

We now extend Definition 3.1 to the case of temporally discretized SDEs in the form of (3.5).

Definition 3.5. Let us assume that we have $L=M$ independent numerical solutions $z_{0:N}^{i}$ of (3.5). We introduce an $M\times M$ matrix $Q_{n}$ for each $n=1,\ldots ,N$ with entries

$$\begin{eqnarray}q_{lj}=q_{+}(z_{n}^{l}|z_{n-1}^{j}):=\text{n}(z_{n}^{l};z_{n-1}^{j}+\unicode[STIX]{x1D6E5}tf(z_{n-1}^{j}),\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6E5}t\,I).\end{eqnarray}$$

With each $Q_{n}$ we associate a finite-state Markov chain $Q_{n}^{+}$ as defined by (3.2) for general transition densities $q_{+}$ in Definition 3.1. An approximation of the Markov transition from time $t_{0}=0$ to $t_{1}=1$ is now provided by

(3.8) $$\begin{eqnarray}Q_{+}:=\mathop{\prod }_{n=1}^{N}Q_{n}^{+}.\end{eqnarray}$$

Remark 3.6. The approximation (3.2) can be related to the diffusion map approximation of the infinitesimal generator of Brownian dynamics

(3.9) $$\begin{eqnarray}\text{d}Z_{t}^{+}=-\unicode[STIX]{x1D6FB}_{z}U(Z_{t}^{+})\,\text{d}t+\sqrt{2}\,\,\text{d}W_{t}^{+}\end{eqnarray}$$

with potential $U(z)=-\log \unicode[STIX]{x1D70B}^{\ast }(z)$ in the following sense. First note that $\unicode[STIX]{x1D70B}^{\ast }$ is invariant under the associated Fokker–Planck equation (2.26) with (time-independent) operator ${\mathcal{L}}^{\dagger }$ written in the form

$$\begin{eqnarray}{\mathcal{L}}^{\dagger }\unicode[STIX]{x1D70B}=\unicode[STIX]{x1D6FB}_{z}\cdot \biggl(\unicode[STIX]{x1D70B}^{\ast }\unicode[STIX]{x1D6FB}_{z}\displaystyle \frac{\unicode[STIX]{x1D70B}}{\unicode[STIX]{x1D70B}^{\ast }}\biggr).\end{eqnarray}$$

Let $z^{i}$ , $i=1,\ldots ,M$ , denote $M$ samples from the invariant PDF $\unicode[STIX]{x1D70B}^{\ast }$ and define the symmetric matrix $Q\in \mathbb{R}^{M\times M}$ with entries

$$\begin{eqnarray}q_{lj}=\text{n}(z^{l};z^{j},2\unicode[STIX]{x1D6E5}t\,I).\end{eqnarray}$$

Then the associated (symmetric) matrix (3.2), as introduced in Definition 3.1, provides a discrete approximation to the evolution of a probability vector $p_{0}\propto \unicode[STIX]{x1D70B}_{0}/\unicode[STIX]{x1D70B}^{\ast }$ over a time-interval $\unicode[STIX]{x1D6E5}t$ and, hence, to the semigroup operator $\text{e}^{\unicode[STIX]{x1D6E5}t{\mathcal{L}}}$ with the infinitesimal generator ${\mathcal{L}}$ given by

(3.10) $$\begin{eqnarray}{\mathcal{L}}g=\displaystyle \frac{1}{\unicode[STIX]{x1D70B}^{\ast }}\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D70B}^{\ast }\unicode[STIX]{x1D6FB}_{z}g).\end{eqnarray}$$

We formally obtain

(3.11) $$\begin{eqnarray}{\mathcal{L}}\approx \displaystyle \frac{Q_{+}-I}{\unicode[STIX]{x1D6E5}t}\end{eqnarray}$$

for $\unicode[STIX]{x1D6E5}t$ sufficiently small. The symmetry of $Q_{+}$ reflects the fact that ${\mathcal{L}}$ is self-adjoint with respect to the weighted inner product

$$\begin{eqnarray}\langle f,g\rangle _{\unicode[STIX]{x1D70B}^{\ast }}=\int f(z)\,g(z)\,\unicode[STIX]{x1D70B}^{\ast }(z)\,\text{d}z.\end{eqnarray}$$

See Harlim (Reference Harlim2018) for a discussion of alternative diffusion map approximations to the infinitesimal generator ${\mathcal{L}}$ and Appendix A.1 for an application to the feedback particle filter formulation of continuous-time data assimilation.

We also consider the discretization

(3.12) $$\begin{eqnarray}Z_{n+1}^{+}=Z_{n}^{+}+(f_{t_{n}}(Z_{n}^{+})+u_{t_{n}}(Z_{n}^{+}))\unicode[STIX]{x1D6E5}t+(\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6E5}t)^{1/2}\,\unicode[STIX]{x1D6EF}_{n},\end{eqnarray}$$

$n=0,\ldots ,N-1$ , of a controlled SDE (2.30) with associated PDF $\unicode[STIX]{x1D70B}_{0:N}^{u}$ defined by

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D70B}_{0:N}^{u}(z_{0:N}^{u})\propto \exp \biggl(-\displaystyle \frac{1}{2\unicode[STIX]{x1D6E5}t}\mathop{\sum }_{n=0}^{N-1}\Vert \unicode[STIX]{x1D702}_{n}^{u}\Vert ^{2}\biggr)\unicode[STIX]{x1D70B}_{0}(z_{0}),\end{eqnarray}$$

where

$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{n}^{u} & := & \displaystyle \unicode[STIX]{x1D6FE}^{-1/2}\{z_{n+1}^{u}-z_{n}^{u}-(f_{t_{n}}(z_{n}^{u})+u_{t_{n}}(z_{n}^{u}))\unicode[STIX]{x1D6E5}t\}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D702}_{n}-\displaystyle \frac{\unicode[STIX]{x1D6E5}t}{\unicode[STIX]{x1D6FE}^{1/2}}u_{t_{n}}(z_{n}^{u}).\nonumber\end{eqnarray}$$

Here $z_{0:N}^{u}$ denotes a realization of the discretization (3.12) with control laws $u_{t_{n}}$ . We find that

$$\begin{eqnarray}\displaystyle \displaystyle \frac{1}{2\unicode[STIX]{x1D6E5}t}\Vert \unicode[STIX]{x1D702}_{n}^{u}\Vert ^{2} & = & \displaystyle \displaystyle \frac{1}{2\unicode[STIX]{x1D6E5}t}\Vert \unicode[STIX]{x1D702}_{n}\Vert ^{2}-\displaystyle \frac{1}{\unicode[STIX]{x1D6FE}^{1/2}}u_{t_{n}}(z_{n}^{u})^{\text{T}}\unicode[STIX]{x1D702}_{n}+\displaystyle \frac{\unicode[STIX]{x1D6E5}t}{2\unicode[STIX]{x1D6FE}}\Vert u_{t_{n}}(z_{n}^{u})\Vert ^{2}\nonumber\\ \displaystyle & = & \displaystyle \displaystyle \frac{1}{2\unicode[STIX]{x1D6E5}t}\Vert \unicode[STIX]{x1D702}_{n}\Vert ^{2}-\displaystyle \frac{1}{\unicode[STIX]{x1D6FE}^{1/2}}u_{t_{n}}(z_{n}^{u})^{\text{T}}\unicode[STIX]{x1D702}_{n}^{u}-\displaystyle \frac{\unicode[STIX]{x1D6E5}t}{2\unicode[STIX]{x1D6FE}}\Vert u_{t_{n}}(z_{n}^{u})\Vert ^{2},\nonumber\end{eqnarray}$$

and hence

(3.14) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D70B}_{0:N}^{u}(z_{0:N}^{u})}{\unicode[STIX]{x1D70B}_{0:N}(z_{0:N}^{u})}=\exp \biggl(\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}}\mathop{\sum }_{n=0}^{N-1}(\Vert u_{t_{n}}(z_{n}^{u})\Vert ^{2}\unicode[STIX]{x1D6E5}t+2\unicode[STIX]{x1D6FE}^{1/2}u_{t_{n}}(z_{n}^{u})^{\text{T}}\unicode[STIX]{x1D702}_{n}^{u})\biggr),\end{eqnarray}$$

which provides a discrete version of (2.33) since $\unicode[STIX]{x1D702}_{n}^{u}=\unicode[STIX]{x1D6E5}t^{1/2}\,\unicode[STIX]{x1D6EF}_{n}(\unicode[STIX]{x1D714})$ are increments of Brownian motion over time intervals of length $\unicode[STIX]{x1D6E5}t$ .

Remark 3.7. Instead of discretizing the forward SDE (2.24) in order to produce samples from the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ , one can also start from the mean-field formulation (2.29) and its time discretization, for example,

(3.15) $$\begin{eqnarray}z_{n+1}^{i}=z_{n}^{i}+(f_{t_{n}}(z_{n}^{i})+u_{t_{n}}(z_{n}^{i}))\unicode[STIX]{x1D6E5}t\end{eqnarray}$$

for $i=1,\ldots ,M$ and

$$\begin{eqnarray}u_{t_{n}}(z)=-\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6FB}_{z}\log \widetilde{\unicode[STIX]{x1D70B}}_{n}(z).\end{eqnarray}$$

Here $\widetilde{\unicode[STIX]{x1D70B}}_{n}$ stands for an approximation to the marginal PDF $\unicode[STIX]{x1D70B}_{t_{n}}$ based on the available samples $z_{n}^{i}$ , $i=1,\ldots ,M$ . A simple approximation is obtained by the Gaussian PDF

$$\begin{eqnarray}\widetilde{\unicode[STIX]{x1D70B}}_{n}(z)=\text{n}(z;\bar{z}_{n},P_{n}^{zz})\end{eqnarray}$$

with empirical mean

$$\begin{eqnarray}\bar{z}_{n}=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}z_{n}^{i}\end{eqnarray}$$

and empirical covariance matrix

$$\begin{eqnarray}P_{n}^{zz}=\displaystyle \frac{1}{M-1}\mathop{\sum }_{i=1}^{M}z_{n}^{i}(z_{n}^{i}-\bar{z}_{n})^{\text{T}}.\end{eqnarray}$$

The system (3.15) becomes

$$\begin{eqnarray}z_{n+1}^{i}=z_{n}^{i}+(f_{t_{n}}(z_{n}^{i})+\unicode[STIX]{x1D6FE}(P_{n}^{zz})^{-1}(z_{n}^{i}-\bar{z}_{n}))\unicode[STIX]{x1D6E5}t,\end{eqnarray}$$

$i=1,\ldots ,M$ , and provides an example of an interacting particle approximation. Similar mean-field formulations can be found for the backward SDE (2.28).

3.2 Filtering

Let us assume that we are given $M$ samples, $z_{1}^{i}$ , from the forecast PDF using forward transition kernels $q_{+}(\cdot \,|z_{0}^{i})$ , $i=1,\ldots ,M$ . The likelihood function $l(z)$ leads to importance weights

(3.16) $$\begin{eqnarray}w^{i}\propto l(z_{1}^{i}).\end{eqnarray}$$

We also normalize these importance weights such that (2.19) holds.

Remark 3.8. The model evidence $\unicode[STIX]{x1D6FD}$ can be estimated from the samples, $z_{1}^{i}$ , and the likelihood $l(z)$ as follows:

$$\begin{eqnarray}\widetilde{\unicode[STIX]{x1D6FD}}:=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}l(z_{1}^{i}).\end{eqnarray}$$

If the likelihood is of the form

$$\begin{eqnarray}l(z)\propto \exp \biggl(-\displaystyle \frac{1}{2}(y_{1}-h(z))^{\text{T}}R^{-1} (y_{1}-h(z)\biggr)\end{eqnarray}$$

and the prior distribution in $y=h(z)$ can be approximated as being Gaussian with covariance

$$\begin{eqnarray}P^{hh}:=\displaystyle \frac{1}{M-1}\mathop{\sum }_{i=1}^{M}h(z_{1}^{i})(h(z_{1}^{i})-\bar{h})^{\text{T}},\quad \bar{h}:=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}h(z_{1}^{i}),\end{eqnarray}$$

then the evidence can be approximated by

$$\begin{eqnarray}\widetilde{\unicode[STIX]{x1D6FD}}\approx \displaystyle \frac{1}{(2\unicode[STIX]{x1D70B})^{N_{y}/2}|P^{yy}|^{1/2}}\exp \biggl(-\displaystyle \frac{1}{2}(y_{1}-\bar{h})^{\text{T}}(P^{yy})^{-1}(y_{1}-\bar{h})\biggr)\end{eqnarray}$$

with

$$\begin{eqnarray}P^{yy}:=R+P^{hh}.\end{eqnarray}$$

Such an approximation has been used, for example, in Carrassi, Bocquet, Hannart and Ghil (Reference Carrassi, Bocquet, Hannart and Ghil2017). See also Reich and Cotter (Reference Reich and Cotter2015) for more details on how to compute and use model evidence in the context of sequential data assimilation.

Sequential data assimilation requires that we produce $M$ equally weighted samples $\widehat{z}_{1}^{j}\sim \widehat{\unicode[STIX]{x1D70B}}_{1}$ from the $M$ weighted samples $z_{1}^{i}\sim \unicode[STIX]{x1D70B}_{1}$ with weights $w^{i}$ . This is a standard problem in Monte Carlo integration and there are many ways to tackle this problem, among which are multinomial, residual, systematic and stratified resampling (Douc and Cappe Reference Douc and Cappe2005). Here we focus on those resampling methods which are based on a discrete Markov chain $P\in \mathbb{R}^{M\times M}$ with the property that

(3.17) $$\begin{eqnarray}w=\displaystyle \frac{1}{M}P\,\unicode[STIX]{x1D7D9}_{M},\quad w=\biggl(\displaystyle \frac{w^{1}}{M},\ldots ,\displaystyle \frac{w^{M}}{M}\biggr)^{\text{T}}.\end{eqnarray}$$

The Markov property of $P$ implies that $P^{\text{T}}\unicode[STIX]{x1D7D9}_{M}=\unicode[STIX]{x1D7D9}_{M}$ . We now consider the set of all Markov chains $\unicode[STIX]{x1D6F1}_{\text{M}}$ , as defined by (3.3), with $p_{1}=w$ . Any Markov chain $P\in \unicode[STIX]{x1D6F1}_{\text{M}}$ can now be used for resampling, but we seek the Markov chain $P^{\ast }\in \unicode[STIX]{x1D6F1}_{\text{M}}$ which minimizes the expected distance between the samples, that is,

(3.18) $$\begin{eqnarray}P^{\ast }=\arg \min _{P\in \unicode[STIX]{x1D6F1}_{\text{M}}}\mathop{\sum }_{i,j=1}^{M}p_{ij}\Vert z_{1}^{i}-z_{1}^{j}\Vert ^{2}.\end{eqnarray}$$

Note that (3.18) is a special case of the optimal transport problem (2.36) with the involved probability measures being discrete measures. Resampling can now be performed according to

(3.19) $$\begin{eqnarray}\mathbb{P}[\widehat{Z}_{1}^{j}(\unicode[STIX]{x1D714})=z_{1}^{i}]=p_{ij}^{\ast }\end{eqnarray}$$

for $j=1,\ldots ,M$ .

Since it is known that (3.18) converges to (2.36) as $M\rightarrow \infty$ (McCann Reference McCann1995) and since (2.36) leads to a transformation (2.38), the resampling step (3.19) has been replaced by

(3.20) $$\begin{eqnarray}\widehat{z}_{1}^{j}=\mathop{\sum }_{i=1}^{M}z_{1}^{i}\,p_{ij}^{\ast }\end{eqnarray}$$

in the so-called ensemble transform particle filter (ETPF) (Reich Reference Reich2013, Reich and Cotter Reference Reich and Cotter2015). In other words, the ETPF replaces resampling with probabilities $p_{ij}^{\ast }$ by its mean (3.20) for each $j=1,\ldots ,M$ . The ETPF leads to a biased approximation to the resampling step which is consistent in the limit $M\rightarrow \infty$ .

The general formulation (3.20) with the coefficients $p_{ij}^{\ast }$ chosen appropriatelyFootnote ² leads to a large class of so-called ensemble transform particle filters (Reich and Cotter Reference Reich and Cotter2015). Ensemble transform particle filters generally result in biased and inconsistent but robust estimates, which have found applications to high-dimensional state space models (Evensen Reference Evensen2006, Vetra-Carvalho et al. Reference Vetra-Carvalho, van Leeuwen, Nerger, Barth, Altaf, Brasseur, Kirchgessner and Beckers2018) for which traditional particle filters fail due to the ‘curse of dimensionality’ (Bengtsson, Bickel and Li Reference Bengtsson, Bickel and Li2008). More specifically, the class of ensemble transform particle filters includes the popular ensemble Kalman filters (Evensen Reference Evensen2006, Reich and Cotter Reference Reich and Cotter2015, Vetra-Carvalho et al. Reference Vetra-Carvalho, van Leeuwen, Nerger, Barth, Altaf, Brasseur, Kirchgessner and Beckers2018, Carrassi, Bocquet, Bertino and Evensen Reference Carrassi, Bocquet, Bertino and Evensen2018) and so-called second-order accurate particle filters with coefficients $p_{ij}^{\ast }$ in (3.20) chosen such that the weighted ensemble mean

$$\begin{eqnarray}\bar{z}_{1}:=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}w^{i}z_{1}^{i}\end{eqnarray}$$

and the weighted ensemble covariance matrix

$$\begin{eqnarray}\widetilde{P}^{zz}:=\displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}w^{i}(z_{1}^{i}-\bar{z}_{1})(z_{1}^{i}-\bar{z}_{1})^{\text{T}}\end{eqnarray}$$

are exactly reproduced by the transformed and equally weighted particles $\widehat{z}_{1}^{j}$ , $j=1,\ldots ,M$ , defined by (3.20), that is,

$$\begin{eqnarray}\displaystyle \frac{1}{M}\mathop{\sum }_{j=1}^{M}\widehat{z}_{1}^{j}=\bar{z}_{1},\quad \displaystyle \frac{1}{M-1}\mathop{\sum }_{j=1}^{M}(\widehat{z}_{1}^{i}-\bar{z}_{1})(\widehat{z}_{1}^{i}-\bar{z}_{1})^{\text{T}}=\widetilde{P}^{zz}.\end{eqnarray}$$

See the survey paper by Vetra-Carvalho et al. (Reference Vetra-Carvalho, van Leeuwen, Nerger, Barth, Altaf, Brasseur, Kirchgessner and Beckers2018) and the paper by Acevedo et al. (Reference Acevedo, de Wiljes and Reich2017) for more details. A summary of the ensemble Kalman filter can be found in Appendix A.3.

In addition, hybrid methods (Frei and Künsch Reference Frei and Künsch2013, Chustagulprom, Reich and Reinhardt Reference Chustagulprom, Reich and Reinhardt2016), which bridge between classical particle filters and the ensemble Kalman filter, have recently been successfully applied to atmospheric fluid dynamics (Robert, Leuenberger and Künsch Reference Robert, Leuenberger and Künsch2018).

Remark 3.9. Another approach for transforming samples, $z_{1}^{i}$ , from the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ into samples, $\widehat{z}_{1}^{i}$ , from the filtering PDF $\widehat{\unicode[STIX]{x1D70B}}_{1}$ is provided through the mean-field interpretation

(3.21) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}s}\breve{Z}_{s}=-\unicode[STIX]{x1D6FB}_{z}\log \displaystyle \frac{\breve{\unicode[STIX]{x1D70B}}_{s}(\breve{Z}_{s})}{\widehat{\unicode[STIX]{x1D70B}}_{1}(\breve{Z}_{s})}\end{eqnarray}$$

of the Fokker–Planck equation (2.26) for a random variable $\breve{Z}_{s}$ with law $\breve{\unicode[STIX]{x1D70B}}_{s}$ , drift term $f_{s}(z)=\unicode[STIX]{x1D6FB}_{z}\log \widehat{\unicode[STIX]{x1D70B}}_{1}$ and $\unicode[STIX]{x1D6FE}=2$ , that is,

$$\begin{eqnarray}\unicode[STIX]{x2202}_{s}\breve{\unicode[STIX]{x1D70B}}_{s}=\unicode[STIX]{x1D6FB}_{z}\cdot \biggl(\breve{\unicode[STIX]{x1D70B}}_{s}\unicode[STIX]{x1D6FB}_{z}\log \displaystyle \frac{\breve{\unicode[STIX]{x1D70B}}_{s}}{\widehat{\unicode[STIX]{x1D70B}}_{1}}\biggr)\end{eqnarray}$$

in artificial time $s\geq 0$ . It holds under fairly general assumptions that

$$\begin{eqnarray}\lim _{s\rightarrow \infty }\breve{\unicode[STIX]{x1D70B}}_{s}=\widehat{\unicode[STIX]{x1D70B}}_{1}\end{eqnarray}$$

(Pavliotis Reference Pavliotis2014) and one can set $\breve{\unicode[STIX]{x1D70B}}_{0}=\unicode[STIX]{x1D70B}_{1}$ . The more common approach would be to solve Brownian dynamics

$$\begin{eqnarray}\,\text{d}\breve{Z}_{s}=\unicode[STIX]{x1D6FB}_{z}\log \widehat{\unicode[STIX]{x1D70B}}_{1}(Z_{s})\,\text{d}t+\sqrt{2}\,\text{d}W_{s}^{+}\end{eqnarray}$$

for each sample $z_{1}^{i}$ from $\unicode[STIX]{x1D70B}_{1}$ , i.e.  $\breve{Z}_{0}(\unicode[STIX]{x1D714})=z_{1}^{i}$ , $i=1,\ldots ,M$ , at initial time and

$$\begin{eqnarray}\widehat{z}_{1}^{i}=\lim _{s\rightarrow \infty }\breve{Z}_{s}(\unicode[STIX]{x1D714}).\end{eqnarray}$$

In other words, formulation (3.21) replaces stochastic Brownian dynamics with a deterministic interacting particle system. See Appendix A.1 and Remark 4.3 for further details.

3.3 Smoothing

Recall that the joint smoothing distribution $\widehat{\unicode[STIX]{x1D70B}}(z_{0},z_{1})$ can be represented in the form (2.47) with modified transition kernel (2.10) and smoothing distribution (2.8) at time $t_{0}$ with weights $\unicode[STIX]{x1D6FE}^{i}$ determined by (2.9).

Let us assume that it is possible to sample from $\widehat{q}_{+}(z_{1}|z_{0}^{i})$ and that the weights $\unicode[STIX]{x1D6FE}^{i}$ are available. Then we can utilize (2.47) in sequential data assimilation as follows. We first resample the $z_{0}^{i}$ at time $t_{0}$ using a discrete Markov chain $P\in \mathbb{R}^{M\times M}$ satisfying

(3.22) $$\begin{eqnarray}\widehat{p}_{0}=\displaystyle \frac{1}{M}P\,\unicode[STIX]{x1D7D9}_{M},\quad \widehat{p}_{0}:=\unicode[STIX]{x1D6FE},\end{eqnarray}$$

with $\unicode[STIX]{x1D6FE}$ defined in (2.12). Again optimal transportation can be used to identify a suitable $P$ . More explicitly, we now consider the set of all Markov chains $\unicode[STIX]{x1D6F1}_{\text{M}}$ , as defined by (3.3), with $p_{1}=\unicode[STIX]{x1D6FE}$ . Then the Markov chain $P^{\ast }$ arising from the associated optimal transport problem (3.18) can be used for resampling, that is,

$$\begin{eqnarray}\mathbb{P}[\widetilde{Z}_{0}^{j}(\unicode[STIX]{x1D714})=z_{0}^{i}]=p_{ij}^{\ast }.\end{eqnarray}$$

Once equally weighted samples $\widehat{z}_{0}^{i}$ , $i=1,\ldots ,M$ , from $\widehat{\unicode[STIX]{x1D70B}}_{0}$ have been determined, the desired samples $\widehat{z}_{1}^{i}$ from $\widehat{\unicode[STIX]{x1D70B}}_{1}$ are simply given by

$$\begin{eqnarray}\widehat{z}_{1}^{i}:=\widehat{Z}_{1}^{i}(\unicode[STIX]{x1D714}),\quad \widehat{Z}_{1}^{i}\sim \widehat{q}_{+}(\cdot \,|\widehat{z}_{0}^{i}),\end{eqnarray}$$

for $i=1,\ldots ,M$ .

The required transition kernels (2.10) are explicitly available for state space models with Gaussian model errors and Gaussian likelihood functions. In many other cases, these kernels are not explicitly available or are difficult to draw from. In such cases, one can resort to sample-based transition kernels.

For example, consider the twisted discrete Markov kernel (3.4) with twisting potential $\unicode[STIX]{x1D713}(z)=l(z)$ . The associated vector $v$ from (3.4) then gives rise to a probability vector $\widehat{p}_{0}=c\,v\in \mathbb{R}^{M}$ with $c>0$ an appropriate scaling factor, and

(3.23) $$\begin{eqnarray}\widehat{p}_{1}:=Q_{+}^{\unicode[STIX]{x1D713}}\widehat{p}_{0}\end{eqnarray}$$

approximates the filtering distribution at time $t_{1}$ . The Markov transition matrix $Q_{+}^{\unicode[STIX]{x1D713}}\in \mathbb{R}^{L\times M}$ together with $\widehat{p}_{0}$ provides an approximation to the smoothing kernel $\widehat{q}_{+}(z_{1}|z_{0})$ and $\widehat{\unicode[STIX]{x1D70B}}_{0}$ , respectively.

The approximations $Q_{+}^{\unicode[STIX]{x1D713}}\in \mathbb{R}^{L\times M}$ and $\widehat{p}_{0}\in \mathbb{R}^{M}$ can be used to first generate equally weighted samples $\widehat{z}_{0}^{i}\in \{z_{0}^{1},\ldots ,z_{0}^{M}\}$ with distribution $\widehat{p}_{0}$ via, for example, resampling with replacement. If $\widehat{z}_{0}^{i}=z_{0}^{k}$ for an index $k=k(i)\in \{1,\ldots ,M\}$ , then

$$\begin{eqnarray}\mathbb{P}[\widehat{Z}_{1}^{i}(\unicode[STIX]{x1D714})=z_{1}^{l}]=(Q_{+}^{\unicode[STIX]{x1D713}})_{lk}\end{eqnarray}$$

for each $i=1,\ldots ,M$ . The $\widehat{z}_{1}^{i}$ are equally weighted samples from the discrete filtering distribution $\widehat{p}_{1}$ , which is an approximation to the continuous filtering PDF $\widehat{\unicode[STIX]{x1D70B}}_{1}$ .

Remark 3.10. One has to take computational complexity and robustness into account when deciding whether to utilize methods from Section 3.2 or this subsection to advance $M$ samples $z_{0}^{i}$ from the prior distribution $\unicode[STIX]{x1D70B}_{0}$ into $M$ samples $\widehat{z}_{1}^{i}$ from the posterior distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ . While the methods from Section 3.2 are easier to implement, the methods of this subsection benefit from the fact that

$$\begin{eqnarray}M>\displaystyle \frac{1}{\Vert \unicode[STIX]{x1D6FE}\Vert ^{2}}\geq \displaystyle \frac{1}{\Vert w\Vert ^{2}}\geq 1,\end{eqnarray}$$

in general, where the importance weights $\unicode[STIX]{x1D6FE}\in \mathbb{R}^{M}$ and $w\in \mathbb{R}^{M}$ are defined in (2.12) and (3.17), respectively. In other words, the methods from this subsection lead to larger effective sample sizes (Liu Reference Liu2001, Agapiou et al. Reference Agapiou, Papaspipliopoulos, Sanz-Alonso and Stuart2017).

3.3.1 SDE models (cont.)

After discretization in time, smoothing leads to a change from the forecast PDF (3.6) to

$$\begin{eqnarray}\displaystyle \widehat{\unicode[STIX]{x1D70B}}_{0:N}(z_{0:N}) & := & \displaystyle \displaystyle \frac{l(z_{N})\unicode[STIX]{x1D70B}_{0:N}(z_{0:N})}{\unicode[STIX]{x1D70B}_{0:N}[l]}\nonumber\\ \displaystyle & \propto & \displaystyle \exp \biggl(-\displaystyle \frac{1}{2\unicode[STIX]{x1D6E5}t}\mathop{\sum }_{n=0}^{N-1}\Vert \unicode[STIX]{x1D709}_{n}\Vert ^{2}\biggr)\unicode[STIX]{x1D70B}_{0}(z_{0})\,l(z_{N})\nonumber\end{eqnarray}$$

with $\unicode[STIX]{x1D709}_{n}$ given by (3.7), or, alternatively,

$$\begin{eqnarray}\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{0:N}}{\unicode[STIX]{x1D70B}_{0:N}}(z_{0:N})=\displaystyle \frac{l(z_{N})}{\unicode[STIX]{x1D70B}_{0:N}[l]}.\end{eqnarray}$$

Remark 3.12. Efficient MCMC methods for sampling high-dimensional smoothing distributions can be found in Beskos et al. (Reference Beskos, Girolami, Lan, Farrell and Stuart2017) and Beskos, Pinski, Sanz-Serna and Stuart (Reference Beskos, Pinski, Sanz-Serna and Stuart2011). Improved sampling can also be achieved by using regularized Störmer–Verlet time-stepping methods (Reich and Hundertmark Reference Reich and Hundertmark2011) in a hybrid Monte Carlo method (Liu Reference Liu2001). See Appendix A.2 for more details.

3.4 Schrödinger problem

Recall that the Schrödinger system (2.58)–(2.61) reduces in our context to solving equations (2.64)–(2.65) for the unknown coefficients $\unicode[STIX]{x1D6FC}^{i}$ , $i=1,\ldots ,M$ . In order to make this problem tractable we need to replace the required expectation values with respect to $q_{+}(z_{1}|z_{0}^{j})$ by Monte Carlo approximations. More specifically, let us assume that we have $L\geq M$ samples $z_{1}^{l}$ from the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ . The associated $L\times M$ matrix $Q$ with entries (3.1) provides a discrete approximation to the underlying Markov process defined by $q_{+}(z_{1}|z_{0})$ and initial PDF (2.3).

The importance weights in the associated approximation to the filtering distribution

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{1}(z)=\displaystyle \frac{1}{L}\mathop{\sum }_{l=1}^{L}w^{l}\,\unicode[STIX]{x1D6FF}(z-z_{1}^{l})\end{eqnarray}$$

are given by (3.16) with the weights normalized such that

(3.24) $$\begin{eqnarray}\mathop{\sum }_{l=1}^{L}w^{l}=L.\end{eqnarray}$$

Finding the coefficients $\{\unicode[STIX]{x1D6FC}^{i}\}$ in (2.64)–(2.65) can now be reformulated as finding two vectors $u\in \mathbb{R}^{L}$ and $v\in \mathbb{R}^{M}$ such that

(3.25) $$\begin{eqnarray}P^{\ast }:=D(u)QD(v)^{-1}\end{eqnarray}$$

satisfies $P^{\ast }\in \unicode[STIX]{x1D6F1}_{\text{M}}$ with $p_{1}=w$ in (3.3), that is, more explicitly

(3.26) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{\text{M}}=\bigg\{P\in \mathbb{R}^{L\times M}:\,p_{lj}\geq 0,\,\mathop{\sum }_{l=1}^{L}p_{lj}=1,\,\displaystyle \frac{1}{M}\mathop{\sum }_{j=1}^{M}p_{lj}=\displaystyle \frac{w^{l}}{L}\bigg\}.\end{eqnarray}$$

We note that (3.26) are discrete approximations to (2.67) and (2.68), respectively. The scaling factor $\widehat{\unicode[STIX]{x1D713}}$ in (2.58) is approximated by the vector $v$ up to a normalization constant, while the vector $u$ provides an approximation to $\unicode[STIX]{x1D713}$ in (2.59). Finally, the desired approximations to the Schrödinger transition kernels $q_{+}^{\ast }(z_{1}|z_{0}^{i})$ , $i=1,\ldots ,M$ , are provided by the columns of $P^{\ast }$ , that is,

$$\begin{eqnarray}\mathbb{P}[\widehat{Z}_{1}^{i}(\unicode[STIX]{x1D714})=z_{1}^{l}]=p_{li}^{\ast }\end{eqnarray}$$

characterizes the desired equally weighted samples $\widehat{z}_{1}^{i}$ , $i=1,\ldots ,M$ , from the filtering distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ . See the following subsection for more details.

The required vectors $u$ and $v$ can be computed using the iterative Sinkhorn algorithm (2.74)–(2.76) (Cuturi Reference Cuturi and Burges2013, Peyre and Cuturi Reference Peyre and Cuturi2018).

Remark 3.14. The approximation (3.25) can be extended to an approximation of the Schrödinger forward transition kernels (2.66) in the following sense. We use $\unicode[STIX]{x1D6FC}^{i}=1/v_{i}$ in (2.66) and note that the resulting approximation satisfies (2.68) while (2.67) now longer holds exactly. However, since the entries $u_{l}$ of the vector $u$ appearing in (3.25) satisfy

$$\begin{eqnarray}u_{l}=\displaystyle \frac{w^{l}}{L}\displaystyle \frac{1}{\mathop{\sum }_{j=1}^{M}q_{+}(z_{1}^{l}|z_{0}^{j})/v_{j}},\end{eqnarray}$$

it follows that

$$\begin{eqnarray}\displaystyle \int q_{+}^{\ast }(z_{1}|z_{0}^{i})\,\text{d}z_{1} & {\approx} & \displaystyle \displaystyle \frac{1}{L}\mathop{\sum }_{l=1}^{L}\displaystyle \frac{l(z_{1}^{l})}{\unicode[STIX]{x1D6FD}}\displaystyle \frac{q_{+}(z_{1}^{l}|z_{0}^{i})/v_{i}}{\mathop{\sum }_{j=1}^{M}q_{+}(z_{1}^{l}|z_{0}^{j})/v_{j}}\nonumber\\ \displaystyle & {\approx} & \displaystyle \displaystyle \frac{1}{L}\mathop{\sum }_{l=1}^{L}w^{l}\displaystyle \frac{q_{+}(z_{1}^{l}|z_{0}^{i})/v_{i}}{\mathop{\sum }_{j=1}^{M}q_{+}(z_{1}^{l}|z_{0}^{j})/v_{j}}=\mathop{\sum }_{l=1}^{L}p_{li}^{\ast }=1.\nonumber\end{eqnarray}$$

Furthermore, one can use such continuous approximations in combination with Monte Carlo sampling methods which do not require normalized target PDFs.

3.4.1 Gaussian model errors (cont.)

One can easily generate $L$ , $L\geq M$ , i.i.d. samples $z_{1}^{l}$ from the forecast PDF (2.21), that is,

$$\begin{eqnarray}Z_{1}^{l}\sim \displaystyle \frac{1}{M}\mathop{\sum }_{j=1}^{M}\text{n}(\cdot \,;\unicode[STIX]{x1D6F9}(z_{0}^{j}),\unicode[STIX]{x1D6FE}B),\end{eqnarray}$$

and with the filtering distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ characterized through the importance weights (3.16).

We define the distance matrix $D\in \mathbb{R}^{L\times M}$ with entries

$$\begin{eqnarray}d_{lj}:=\displaystyle \frac{1}{2}\Vert z_{1}^{l}-\unicode[STIX]{x1D6F9}(z_{0}^{j})\Vert _{B}^{2},\quad \Vert z\Vert _{B}^{2}:=z^{\text{T}}B^{-1}z,\end{eqnarray}$$

and the matrix $Q\in \mathbb{R}^{L\times M}$ with entries

$$\begin{eqnarray}q_{lj}:=\text{e}^{-d_{lj}/\unicode[STIX]{x1D6FE}}.\end{eqnarray}$$

The Markov chain $P^{\ast }\in \mathbb{R}^{L\times M}$ is now given by

$$\begin{eqnarray}P^{\ast }=\arg \min _{P\in \unicode[STIX]{x1D6F1}_{\text{M}}}\text{KL}(P||Q)\end{eqnarray}$$

with the set $\unicode[STIX]{x1D6F1}_{\text{M}}$ defined by (3.26).

Once $P^{\ast }$ has been computed, the desired Schrödinger transitions from $\unicode[STIX]{x1D70B}_{0}$ to $\widehat{\unicode[STIX]{x1D70B}}_{1}$ can be represented as follows. The Schrödinger transition kernels $q_{+}^{\ast }(z_{1}|z_{0}^{i})$ are approximated for each $z_{0}^{i}$ by

(3.27) $$\begin{eqnarray}\tilde{q}_{+}^{\ast }(z_{1}|z_{0}^{i}):=\mathop{\sum }_{l=1}^{L}p_{li}^{\ast }\,\unicode[STIX]{x1D6FF}(z_{1}-z_{1}^{l}),\quad i=1,\ldots ,M.\end{eqnarray}$$

The empirical measure in (3.27) converges weakly to the desired $q_{+}^{\ast }(z_{1}|z_{0}^{i})$ as $L\rightarrow \infty$ and

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1})\approx \displaystyle \frac{1}{M}\mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D6FF}(z-\widehat{z}_{1}^{i}),\end{eqnarray}$$

with

(3.28) $$\begin{eqnarray}\widehat{z}_{1}^{i}=\widehat{Z}_{1}^{i}(\unicode[STIX]{x1D714}),\quad \widehat{Z}_{1}^{i}\sim \tilde{q}_{+}^{\ast }(\cdot \,|z_{0}^{i}),\end{eqnarray}$$

provides the desired approximation of $\widehat{\unicode[STIX]{x1D70B}}_{1}$ by $M$ equally weighted particles $\widehat{z}_{1}^{i}$ , $i=1,\ldots ,M$ .

We remark that (3.27) has been used to produce the Schrödinger transition kernels for Example 2.5 and Figure 2.3(b) in particular. More specifically, we have $M=11$ and used $L=11\,000$ . Since the particles $z_{1}^{l}\in \mathbb{R}$ , $l=1,\ldots ,L$ , are distributed according to the forecast PDF $\unicode[STIX]{x1D70B}_{1}$ , a function representation of $\tilde{q}_{+}^{\ast }(z_{1}|z_{0}^{i})$ over all of $\mathbb{R}$ is provided by interpolating $p_{li}^{\ast }$ onto $\mathbb{R}$ and multiplication of this interpolated function by $\unicode[STIX]{x1D70B}_{1}(z)$ .

For $\unicode[STIX]{x1D6FE}\ll 1$ , the measure in (3.27) can also be approximated by a Gaussian measure with mean

$$\begin{eqnarray}\bar{z}_{1}^{i}:=\mathop{\sum }_{l=1}^{L}z_{1}^{l}p_{li}^{\ast }\end{eqnarray}$$

and covariance matrix $\unicode[STIX]{x1D6FE}B$ , that is, we replace (3.28) with

$$\begin{eqnarray}\widehat{Z}_{1}^{i}\sim \text{N}(\bar{z}_{1}^{i},\unicode[STIX]{x1D6FE}B)\end{eqnarray}$$

for $i=1,\ldots ,M$ .

3.4.2 SDE (cont.)

One can also apply (3.25) in order to approximate the Schrödinger problem associated with SDE models. We typically use $L=M$ in this case and utilize (3.8) in place of $Q$ in (3.25). The set $\unicode[STIX]{x1D6F1}_{\text{M}}$ is still given by (3.26).

Figure 3.1. Histograms produced from $M=200$ Monte Carlo samples of the initial PDF $\unicode[STIX]{x1D70B}_{0}$ , the forecast PDF $\unicode[STIX]{x1D70B}_{2}$ at time $t=2$ , the filtering distribution $\widehat{\unicode[STIX]{x1D70B}}_{2}$ at time $t=2$ , and the smoothing PDF $\widehat{\unicode[STIX]{x1D70B}}_{0}$ at time $t=0$ for a Brownian particle moving in a double well potential.

Example 3.15. We consider scalar-valued motion of a Brownian particle in a bimodal potential, that is,

(3.29) $$\begin{eqnarray}\,\text{d}Z_{t}^{+}=Z_{t}^{+}\,\text{d}t-(Z_{t}^{+})^{3}\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{+}\end{eqnarray}$$

with $\unicode[STIX]{x1D6FE}=0.5$ and initial distribution $Z_{0}\sim \text{N}(-1,0.3)$ . At time $t=2$ we measure the location $y=1$ with measurement error variance $R=0.2$ . We simulate the dynamics using $M=200$ particles and a time-step of $\unicode[STIX]{x1D6E5}t=0.01$ in the Euler–Maruyama discretization (3.5). One can find histograms produced from the Monte Carlo samples in Figure 3.1. The samples from the filtering and smoothing distributions are obtained by resampling with replacement from the weighted distributions with weights given by (3.16). Next we compute (3.8) from the $M=200$ Monte Carlo samples of (3.29). Eleven out of the 200 transition kernels from $\unicode[STIX]{x1D70B}_{0}$ to $\unicode[STIX]{x1D70B}_{2}$ (prediction problem) and $\unicode[STIX]{x1D70B}_{0}$ to $\widehat{\unicode[STIX]{x1D70B}}_{2}$ (Schrödinger problem) are displayed in Figure 3.2.

Figure 3.2. (a) Approximations of typical transition kernels from time $\unicode[STIX]{x1D70B}_{0}$ to $\unicode[STIX]{x1D70B}_{2}$ under the Brownian dynamics model (3.29). (b) Approximations of typical Schrödinger transition kernels from $\unicode[STIX]{x1D70B}_{0}$ to $\widehat{\unicode[STIX]{x1D70B}}_{2}$ . All approximations were computed using the Sinkhorn algorithm and by linear interpolation between the $M=200$ data points.

The Sinkhorn approach requires relatively large sample sizes $M$ in order to lead to useful approximations. Alternatively we may assume that there is an approximative control term $u_{t}^{(0)}$ with associated forward SDE

(3.30) $$\begin{eqnarray}\,\text{d}Z_{t}^{+}=f_{t}(Z_{t})\,\text{d}t+u_{t}^{(0)}(Z_{t}^{+})\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{+},\quad t\in [0,1],\end{eqnarray}$$

and $Z_{0}^{+}\sim \unicode[STIX]{x1D70B}_{0}$ . We denote the associated path measure by $\mathbb{Q}^{(0)}$ . Girsanov’s theorem implies that the Radon–Nikodym derivative of $\mathbb{Q}$ with respect to $\mathbb{Q}^{(0)}$ is given by (compare (2.33))

$$\begin{eqnarray}\displaystyle \frac{\text{d}\mathbb{Q}}{\text{d}\mathbb{Q}^{(0)}}_{|z_{[0,1]}^{(0)}}=\exp (-V^{(0)}),\end{eqnarray}$$

where $V^{(0)}$ is defined via the stochastic integral

$$\begin{eqnarray}V^{(0)}:=\displaystyle \frac{1}{2\unicode[STIX]{x1D6FE}}\int _{0}^{1}(\Vert u_{t}^{(0)}\Vert ^{2}\,\text{d}t+2\unicode[STIX]{x1D6FE}^{1/2}u_{t}^{(0)}\cdot \,\text{d}W_{t}^{+})\end{eqnarray}$$

along solution paths $z_{[0,1]}^{(0)}$ of (3.30). Because of

$$\begin{eqnarray}\displaystyle \frac{\text{d}\widehat{\mathbb{P}}}{\text{d}\mathbb{Q}^{(0)}}_{|z_{[0,1]}^{(0)}}=\displaystyle \frac{\text{d}\widehat{\mathbb{P}}}{\text{d}\mathbb{Q}}_{|z_{[0,1]}^{(0)}}\,\displaystyle \frac{\text{d}\mathbb{Q}}{\text{d}\mathbb{Q}^{(0)}}_{|z_{[0,1]}^{(0)}}\propto l(z_{1}^{(0)})\,\exp (-V^{(0)}),\end{eqnarray}$$

we can now use (3.30) to importance-sample from the filtering PDF $\widehat{\unicode[STIX]{x1D70B}}_{1}$ . The control $u_{t}^{(0)}$ should be chosen such that the variance in the modified likelihood function

(3.31) $$\begin{eqnarray}l^{(0)}(z_{[0,1]}^{(0)}):=l(z_{1}^{(0)})\,\exp (-V^{(0)})\end{eqnarray}$$

is reduced compared to the uncontrolled case $u_{t}^{(0)}\equiv 0$ . In particular, the filter distribution $\widehat{\unicode[STIX]{x1D70B}}_{1}$ at time $t=1$ satisfies

$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D70B}}_{1}(z_{1}^{(0)})\propto l^{(0)}(z_{[0,1]}^{(0)})\,\unicode[STIX]{x1D70B}_{1}^{(0)}(z_{1}^{(0)}),\end{eqnarray}$$

where $\unicode[STIX]{x1D70B}_{t}^{(0)}$ , $t\in (0,1]$ , denote the marginal PDFs generated by (3.30).

We now describe an iterative algorithm for the associated Schrödinger problem in the spirit of the Sinkhorn iteration from Section 2.3.2.

Lemma 3.16. The desired optimal control law can be computed iteratively,

(3.32) $$\begin{eqnarray}u_{t}^{(k+1)}=u_{t}^{(k)}+\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}_{z}\log \unicode[STIX]{x1D713}_{t}^{(k)},\quad k=0,1,\ldots ,\end{eqnarray}$$

for given $u_{t}^{(k)}$ and the potential $\unicode[STIX]{x1D713}_{t}^{(k)}$ obtained as the solutions to the backward Kolmogorov equation

(3.33) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}_{t}^{(k)}=-{\mathcal{L}}_{t}^{(k)}\unicode[STIX]{x1D713}_{t}^{(k)},\quad {\mathcal{L}}_{t}^{(k)}g:=\unicode[STIX]{x1D6FB}_{z}g\cdot (f_{t}+u_{t}^{(k)})+\displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\unicode[STIX]{x1D6E5}_{z}g,\end{eqnarray}$$

with final time condition

(3.34) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{1}^{(k)}(z):=\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{1}(z)}{\unicode[STIX]{x1D70B}_{1}^{(k)}(z)}.\end{eqnarray}$$

Here $\unicode[STIX]{x1D70B}_{1}^{(k)}$ denotes the time-one marginal of the path measure $\mathbb{Q}^{(k)}$ induced by (2.30) with control term $u_{t}=u_{t}^{(k)}$ and initial PDF $\unicode[STIX]{x1D70B}_{0}$ . The recursion (3.32) is stopped whenever the final time condition (3.34) is sufficiently close to a constant function.

Proof. The extension of the Sinkhorn algorithm to continuous PDFs and its convergence has been discussed by Chen, Georgiou and Pavon (Reference Chen, Georgiou and Pavon2016a ).◻

Remark 3.17. Note that $\unicode[STIX]{x1D713}_{t}^{(k)}$ needs to be determined up to a constant of proportionality only since the associated control law is determined from $\unicode[STIX]{x1D713}_{t}^{(k)}$ by (3.32). One can also replace (3.33) with any other method for solving the smoothing problem associated to the SDE (2.30) with $Z_{0}^{+}\sim \unicode[STIX]{x1D70B}_{0}$ , control law $u_{t}=u_{t}^{(k)}$ , and likelihood function $l(z)=\unicode[STIX]{x1D713}_{1}^{(k)}(z)$ . See Appendix A.4 for a forward–backward SDE formulation in particular.

We need to restrict the class of possible control laws $u_{t}^{(k)}$ in order to obtain a computationally feasible implementations in practice. For example, a simple class of control laws is provided by linear controls of the form

$$\begin{eqnarray}u_{t}^{(k)}(z)=-B_{t}^{(k)}(z-m_{t}^{(k)})\end{eqnarray}$$

with appropriately chosen symmetric positive definite matrices $B_{t}^{(k)}$ and vectors $m_{t}^{(k)}$ . Such approximations can, for example, be obtained from the smoother extensions of ensemble transform methods mentioned earlier. See also the recent work by Kappen and Ruiz (Reference Kappen and Ruiz2016) and Ruiz and Kappen (Reference Ruiz and Kappen2017) on numerical methods for the SDE smoothing problem.

4.1 Smooth data

We start from a forward SDE model (2.24) with associated path measure $\mathbb{Q}$ over the space of continuous functions ${\mathcal{C}}$ . However, contrary to the previous sections, the likelihood $l$ is defined along a whole solution path $z_{[0,1]}$ as follows:

$$\begin{eqnarray}\displaystyle \frac{\text{d}\widehat{\mathbb{P}}}{\text{d}\mathbb{Q}}_{|z_{[0,1]}}\propto l(z_{[0,1]}),\quad l(z_{[0,1]}):=\exp \biggl(-\int _{0}^{1}V_{t}(z_{t})\,\text{d}t\biggr)\end{eqnarray}$$

with the assumption that $\mathbb{Q}[l]<\infty$ and $V_{t}(z)\geq 0$ . A specific example of a suitable $V_{t}$ is provided by

(4.1) $$\begin{eqnarray}V_{t}(z)=\displaystyle \frac{1}{2}\Vert h(z)-y_{t}\Vert ^{2},\end{eqnarray}$$

where the data function $y_{t}\in \mathbb{R}$ , $t\in [0,1]$ , is a smooth function of time and $h(z)$ is a forward operator connecting the model states to the observations/data. The associated estimation problem has, for example, been addressed by Mortensen (Reference Mortensen1968) and Fleming (Reference Fleming1997) from an optimal control perspective and has led to what is called the minimum energy estimator. Recall that the filtering PDF $\widehat{\unicode[STIX]{x1D70B}}_{1}$ is the marginal PDF of $\widehat{\mathbb{P}}$ at time $t=1$ .

The associated time-continuous smoothing/filtering problems are based on the time-dependent path measures $\widehat{\mathbb{P}}_{t}$ defined by

$$\begin{eqnarray}\displaystyle \frac{\text{d}\widehat{\mathbb{P}}_{t}}{\text{d}\mathbb{Q}}_{|z_{[0,1]}}\propto l(z_{[0,t]}),\quad l(z_{[0,t]}):=\exp \biggl(-\int _{0}^{t}V_{s}(z_{s})\,\text{d}s\biggr)\end{eqnarray}$$

for $t\in (0,1]$ . We let $\widehat{\unicode[STIX]{x1D70B}}_{t}$ denote the marginal PDF of $\widehat{\mathbb{P}}_{t}$ at time $t$ . Note that $\widehat{\unicode[STIX]{x1D70B}}_{t}$ is the filtering PDF, i.e. the marginal PDF at time $t$ conditioned on all the data available until time $t$ . Also note that $\widehat{\unicode[STIX]{x1D70B}}_{t}$ is different from the marginal (smoothing) PDF of $\widehat{\mathbb{P}}$ at time $t$ .

We now state a modified Fokker–Planck equation which describes the time evolution of the filtering PDFs $\widehat{\unicode[STIX]{x1D70B}}_{t}$ .

Lemma 4.1. The marginal distributions $\widehat{\unicode[STIX]{x1D70B}}_{t}$ of $\widehat{\mathbb{P}}_{t}$ satisfy the modified Fokker–Planck equation

(4.2) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\widehat{\unicode[STIX]{x1D70B}}_{t}={\mathcal{L}}_{t}^{\dagger }\widehat{\unicode[STIX]{x1D70B}}_{t}-\widehat{\unicode[STIX]{x1D70B}}_{t}(V_{t}-\widehat{\unicode[STIX]{x1D70B}}_{t}[V_{t}])\end{eqnarray}$$

with ${\mathcal{L}}_{t}^{\dagger }$ defined by (2.25).

Proof. This can be seen by setting $\unicode[STIX]{x1D6FE}=0$ and $f_{t}\equiv 0$ in (2.24) for simplicity and by considering the incremental change of measure induced by the likelihood, that is,

$$\begin{eqnarray}\displaystyle \frac{\widehat{\unicode[STIX]{x1D70B}}_{t+\unicode[STIX]{x1D6FF}t}}{\widehat{\unicode[STIX]{x1D70B}}_{t}}\propto \text{e}^{-V_{t}\unicode[STIX]{x1D6FF}t}\approx 1-V_{t}\,\unicode[STIX]{x1D6FF}t,\end{eqnarray}$$

and taking the limit $\unicode[STIX]{x1D6FF}t\rightarrow 0$ under the constraint that $\widehat{\unicode[STIX]{x1D70B}}_{t}[1]=1$ is preserved.◻

We now derive a mean-field interpretation of (4.2) and rewrite (4.2) in the form

(4.3) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\widehat{\unicode[STIX]{x1D70B}}_{t}={\mathcal{L}}_{t}^{\dagger }\widehat{\unicode[STIX]{x1D70B}}_{t}+\unicode[STIX]{x1D6FB}_{z}\cdot (\widehat{\unicode[STIX]{x1D70B}}_{t}\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D719}_{t}),\end{eqnarray}$$

where the potential $\unicode[STIX]{x1D719}_{t}:\mathbb{R}^{N_{z}}\rightarrow \mathbb{R}$ satisfies the elliptic PDE

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{z}\cdot (\widehat{\unicode[STIX]{x1D70B}}_{t}\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D719}_{t})=-\widehat{\unicode[STIX]{x1D70B}}_{t}(V_{t}-\widehat{\unicode[STIX]{x1D70B}}_{t}[V_{t}]).\end{eqnarray}$$

With (4.3) in place, we formally obtain the mean-field equation

(4.5) $$\begin{eqnarray}\text{d}Z_{t}^{+}=\{f_{t}(Z_{t}^{+})-\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D719}_{t}(Z_{t}^{+})\}\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{+},\end{eqnarray}$$

and the marginal distributions $\unicode[STIX]{x1D70B}_{t}^{u}$ of this controlled SDE agree with the marginals $\widehat{\unicode[STIX]{x1D70B}}_{t}$ of the path measures $\widehat{\mathbb{P}}_{t}$ at times $t\in (0,1]$ .

The control $u_{t}$ is not uniquely determined. For example, one can replace (4.4) with

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D70B}_{t}M_{t}\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D719}_{t})=-\unicode[STIX]{x1D70B}_{t}(V_{t}-\unicode[STIX]{x1D70B}_{t}[V_{t}]),\end{eqnarray}$$

where $M_{t}$ is a symmetric positive definite matrix. More specifically, let us assume that $\unicode[STIX]{x1D70B}_{t}$ is Gaussian with mean $\bar{z}_{t}$ and covariance matrix $P_{t}^{zz}$ and that $h(z)$ is linear, i.e.  $h(z)=Hz$ . Then (4.6) can be solved analytically for $M_{t}=P_{t}^{zz}$ with

$$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D719}_{t}(z)=\displaystyle \frac{1}{2}H^{\text{T}}(Hz+H\bar{z}_{t}-2y_{t}).\end{eqnarray}$$

The resulting mean-field equation becomes

(4.7) $$\begin{eqnarray}\text{d}Z_{t}^{+}=\bigg\{f_{t}(Z_{t}^{+})-\displaystyle \frac{1}{2}P_{t}^{zz}H^{\text{T}}(HZ_{t}^{+}+H\bar{z}_{t}-2y_{t})\bigg\}\,\text{d}t+\unicode[STIX]{x1D6FE}^{1/2}\,\text{d}W_{t}^{+},\end{eqnarray}$$

which gives rise to the ensemble Kalman–Bucy filter upon Monte Carlo discretization (Bergemann and Reich Reference Bergemann and Reich2012). See Section 5 below and Appendix A.3 for further details.

Remark 4.3. The approach described in this subsection can also be applied to standard Bayesian inference without model dynamics. More specifically, let us assume that we have samples $z_{0}^{i}$ , $i=1,\ldots ,M$ , from a prior distribution $\unicode[STIX]{x1D70B}_{0}$ which we would like to transform into samples from a posterior distribution

$$\begin{eqnarray}\unicode[STIX]{x1D70B}^{\ast }(z):=\displaystyle \frac{l(z)\,\unicode[STIX]{x1D70B}_{0}(z)}{\unicode[STIX]{x1D70B}_{0}[l]}\end{eqnarray}$$

with likelihood $l(z)=\unicode[STIX]{x1D70B}(y|z)$ . One can introduce a homotopy connecting $\unicode[STIX]{x1D70B}_{0}$ with $\unicode[STIX]{x1D70B}^{\ast }$ , for example, via

(4.8) $$\begin{eqnarray}\breve{\unicode[STIX]{x1D70B}}_{s}(z):=\displaystyle \frac{l(z)^{s}\,\unicode[STIX]{x1D70B}_{0}(z)}{\unicode[STIX]{x1D70B}_{0}[l^{s}]}\end{eqnarray}$$

with $s\in [0,1]$ . We find that

(4.9) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\breve{\unicode[STIX]{x1D70B}}_{s}}{\unicode[STIX]{x2202}s}=\breve{\unicode[STIX]{x1D70B}}_{s}(\log l-\breve{\unicode[STIX]{x1D70B}}_{s}[\log l]).\end{eqnarray}$$

We now seek a differential equation

(4.10) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}s}\breve{Z}_{s}=u_{s}(\breve{Z}_{s})\end{eqnarray}$$

with $\breve{Z}_{0}\sim \unicode[STIX]{x1D70B}_{0}$ such that its marginal distributions $\breve{\unicode[STIX]{x1D70B}}_{s}$ satisfy (4.9) and, in particular, $\breve{Z}_{1}\sim \unicode[STIX]{x1D70B}^{\ast }$ . This condition together with Liouville’s equation for the time evolution of marginal densities under a differential equation (4.10) leads to

(4.11) $$\begin{eqnarray}-\unicode[STIX]{x1D6FB}_{z}\cdot (\breve{\unicode[STIX]{x1D70B}}_{s}\,u_{s})=\breve{\unicode[STIX]{x1D70B}}_{s}(\log l-\breve{\unicode[STIX]{x1D70B}}_{s}[\log l]).\end{eqnarray}$$

In order to define $u_{s}$ in (4.10) uniquely, we make the ansatz

(4.12) $$\begin{eqnarray}u_{s}(z)=-\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D719}_{s}(z)\end{eqnarray}$$

which leads to the elliptic PDE

(4.13) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{z}\cdot (\breve{\unicode[STIX]{x1D70B}}_{s}\,\unicode[STIX]{x1D6FB}_{z}\unicode[STIX]{x1D719}_{s})=\breve{\unicode[STIX]{x1D70B}}_{s}(\log l-\breve{\unicode[STIX]{x1D70B}}_{s}[\log l])\end{eqnarray}$$

in the potential $\unicode[STIX]{x1D719}_{s}$ . The desired samples from $\unicode[STIX]{x1D70B}^{\ast }$ are now obtained as the time-one solutions of (4.10) with ‘control law’ (4.12) satisfying (4.13) and initial conditions $z_{0}^{i}$ , $i=1,\ldots ,M$ . There are many modifications of this basic procedure (Daum and Huang Reference Daum and Huang2011, Reich Reference Reich2011, El Moselhy and Marzouk Reference El Moselhy and Marzouk2012), some of them leading to explicit expressions for (4.10) such as Gaussian PDFs (Bergemann and Reich Reference Bergemann and Reich2010) and Gaussian mixture PDFs (Reich Reference Reich2012). We finally mention that the limit $s\rightarrow \infty$ in (4.8) leads, formally, to the PDF $\breve{\unicode[STIX]{x1D70B}}_{\infty }=\unicode[STIX]{x1D6FF}(z-z_{\text{ML}})$ , where $z_{\text{ML}}$ denotes the minimizer of $V(z)=-\log \unicode[STIX]{x1D70B}(y|z)$ , i.e. the maximum likelihood estimator, which we assume here to be unique, for example, $V$ is convex. In other words these homotopy methods can be used to solve optimization problems via derivative-free mean-field equations and their interacting particle approximations. See, for example, Zhang, Taghvaei and Mehta (Reference Zhang, Taghvaei and Mehta2019) and Schillings and Stuart (Reference Schillings and Stuart2017) as well as Appendices A.1 and A.3 for more details.

4.2 Random data

We now replace (4.1) with an observation model of the form

$$\begin{eqnarray}\,\text{d}Y_{t}=h(Z_{t}^{+})\,\text{d}t+\,\text{d}V_{t}^{+},\end{eqnarray}$$

where we set $Y_{t}\in \mathbb{R}$ for simplicity and $V_{t}^{+}$ denotes standard Brownian motion. The forward operator $h:\mathbb{R}^{N_{z}}\rightarrow \mathbb{R}$ is also assumed to be known. The marginal PDFs $\widehat{\unicode[STIX]{x1D70B}}_{t}$ for $Z_{t}$ conditioned on all observations $y_{s}$ with $s\in [0,t]$ satisfy the Kushner–Stratonovitch equation (Jazwinski Reference Jazwinski1970)

(4.14) $$\begin{eqnarray}\text{d}\widehat{\unicode[STIX]{x1D70B}}_{t}={\mathcal{L}}_{t}^{\dagger }\widehat{\unicode[STIX]{x1D70B}}_{t}\,\text{d}t+(h-\widehat{\unicode[STIX]{x1D70B}}_{t}[h])(\text{d}Y_{t}-\widehat{\unicode[STIX]{x1D70B}}_{t}[h]\,\text{d}t)\end{eqnarray}$$

with ${\mathcal{L}}^{\dagger }$ defined by (2.25). The following observation is important for the subsequent discussion.

Remark 4.4. Consider state-dependent diffusion

(4.15) $$\begin{eqnarray}\text{d}Z_{t}^{+}=\unicode[STIX]{x1D6FE}_{t}(Z_{t}^{+})\circ \,\text{d}U_{t}^{+},\end{eqnarray}$$

in its Stratonovitch interpretation (Pavliotis Reference Pavliotis2014), where $U_{t}^{+}$ is scalar-valued Brownian motion and $\unicode[STIX]{x1D6FE}_{t}(z)\in \mathbb{R}^{N_{z}\times 1}$ . Here the Stratonovitch interpretation is to be applied to the implicit time-dependence of $\unicode[STIX]{x1D6FE}_{t}(z)$ through $Z_{t}^{+}$ only, that is, the explicit time-dependence of $\unicode[STIX]{x1D6FE}_{t}$ remains to be Itô-interpreted. The associated Fokker–Planck equation for the marginal PDFs $\unicode[STIX]{x1D70B}_{t}$ takes the form

(4.16) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70B}_{t}=\displaystyle \frac{1}{2}\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D6FE}_{t}\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D70B}_{t}\unicode[STIX]{x1D6FE}_{t})),\end{eqnarray}$$

and expectation values $\bar{g}=\unicode[STIX]{x1D70B}_{t}[g]$ evolve in time according to

(4.17) $$\begin{eqnarray}\unicode[STIX]{x1D70B}_{t}[g]=\unicode[STIX]{x1D70B}_{0}[g]+\int _{0}^{t}\unicode[STIX]{x1D70B}_{s}[{\mathcal{A}}_{t}g]\,\text{d}s\end{eqnarray}$$

with operator ${\mathcal{A}}_{t}$ defined by

$$\begin{eqnarray}{\mathcal{A}}_{t}g=\displaystyle \frac{1}{2}\unicode[STIX]{x1D6FE}_{t}^{\text{T}}\unicode[STIX]{x1D6FB}_{z}(\unicode[STIX]{x1D6FE}_{t}^{\text{T}}\unicode[STIX]{x1D6FB}_{z}g).\end{eqnarray}$$

Now consider the mean-field equation

(4.18) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\widetilde{Z}_{t}=-\displaystyle \frac{1}{2}\unicode[STIX]{x1D6FE}_{t}(\widetilde{Z}_{t})\,J_{t},\quad J_{t}:=\widetilde{\unicode[STIX]{x1D70B}}_{t}^{-1}\unicode[STIX]{x1D6FB}_{z}\cdot (\widetilde{\unicode[STIX]{x1D70B}}_{t}\unicode[STIX]{x1D6FE}_{t}),\end{eqnarray}$$

with $\widetilde{\unicode[STIX]{x1D70B}}_{t}$ the law of $\widetilde{Z}_{t}$ . The associated Liouville equation is

$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\widetilde{\unicode[STIX]{x1D70B}}_{t}=\displaystyle \frac{1}{2}\unicode[STIX]{x1D6FB}_{z}\cdot (\widetilde{\unicode[STIX]{x1D70B}}_{t}\unicode[STIX]{x1D6FE}_{t}J_{t})=\displaystyle \frac{1}{2}\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D6FE}_{t}\unicode[STIX]{x1D6FB}_{z}\cdot (\unicode[STIX]{x1D6FE}_{t}\widetilde{\unicode[STIX]{x1D70B}}_{t})).\end{eqnarray}$$

In other words, the marginal PDFs and the associated expectation values evolve identically under (4.15) and (4.18), respectively.