methods, known as Metropolis-Hastings algorithms. We will prove that the
Metropolis-Hastings kernel is invariant with respect to the desired target
distribution, and we will show _geometric ergodicity_ of the pCN Metropolis-
Hastings algorithm, meaning that the distribution $\pi_{n}$ of the $n$-th
sample approaches the invariant distribution exponentially fast in total
variation distance. The idea is illustrated in Figure 9: after an initial
number of burn-in steps the samples from the chain start to concentrate in
regions where the target distribution has greatest mass. We will not discuss
sample path ergodicity, noting simply that a general abstract theory exists to
deduce it from geometric ergodicity.

Figure 9: The Markov chain samples points from distribution $\pi_{n}$ at step
$n$, and the sampling distribution converges towards the target distribution
$\pi$ whose high density regions are represented by the dashed circles.

####  6.2 The Metropolis-Hastings Algorithm

Here we outline the Metropolis-Hastings algorithm. The algorithm has two
ingredients: a proposal distribution $q(u,v)$, which is a Markov transition
kernel; and an acceptance probability $a(u,v)$ that will be used to convert
the proposal kernel into a kernel $p_{\mbox{\tiny{\rm MH}}}(u,v)$ that is
invariant with respect to the given target $\pi.$ Given the $n^{\rm th}$
sample $u^{(n)},$ we generate $u^{(n+1)}$ by drawing $v^{\star}$ from the
distribution $q(u^{(n)},\cdot)$ and accepting the result, which means setting
$u^{(n+1)}=v^{\star}$, with probability $a(u^{(n)},v^{\star})$, or instead
setting $u^{(n+1)}=u^{(n)}$ with the remaining probability
$1-a(u^{(n)},v^{\star})$. The acceptance probability is given by

$\displaystyle
a(u,v)=\text{min}\Bigg{(}\frac{\pi(v)q(v,u)}{\pi(u)q(u,v)},\,1\Bigg{)}.$ (6.1)

Algorithm 1 makes precise the steps just described:

Algorithm 1 Metropolis-Hastings Algorithm

1: Input: Target distribution $\pi$, initial distribution $\pi_{0},$ Markov
kernel $q(u,v)$, number of samples $N.$

2: Initial Draw: Draw initial sample $u^{(0)}\sim\pi_{0}.$

3: Subsequent Samples: For $n=0,1,\dots,N-1$, perform the following steps:

  1. 1.

Sample $v^{\star}\sim q(u^{(n)},\cdot).$

  2. 2.

Calculate the acceptance probability $a_{n}:=a(u^{(n)},v^{\star}).$

  3. 3.

Update

$\displaystyle u^{(n+1)}=\begin{cases}v^{\star},&\text{w.p.}~{}~{}~{}a_{n},\\\
u^{(n)},&\text{otherwise}.\end{cases}$

4: Output Samples $u^{(1)},u^{(2)},\dots,u^{(N)}.$

The Metropolis-Hastings algorithm implicitly defines a Markov kernel
$p_{\mbox{\tiny{\rm MH}}}(u,\cdot)$ which specifies the density of the
$(n+1)^{\rm th}$ sample given that the $n^{\rm th}$ sample is located at $u.$
For $u\neq v,$ the Metropolis-Hastings kernel has the following simple
expression in terms of the proposal kernel and the acceptance probability

$p_{\mbox{\tiny{\rm MH}}}(u,v)=a(u,v)q(u,v);$ (6.2)

this expression may be deduced noting that in order to move to a new location,
the proposal needs to be accepted.

######  Remark 6.1.

  * •

In order to implement the Metropolis-Hastings algorithm one needs to be able
to sample from the proposal kernel $q(u,\cdot)$ and evaluate the acceptance
probability $a(u,v).$ Importantly, the target distribution only appears in the
acceptance probability $a(u,v),$ and only the ratio $\pi(v)/\pi(u)$ is
involved. Therefore the Metropolis-Hastings algorithm may be implemented for
target distributions that are only specified up to an unknown normalizing
constant.

  * •

If $q(u,v)=q(v,u)$ the acceptance probability simplifies to $\pi(v)/\pi(u).$
In such a case, moves to regions of higher target density are always accepted,
while moves to regions of smaller but non-zero target density are accepted
with positive probability in order to ensure exploration of the target space.
If $q$ is not symmetric, the method favors moves that are easier to be
reversed, namely moves for which $q(v,u)>q(u,v).$

  * •

The Metropolis-Hastings algorithm is extremely flexible due to the freedom in
the choice of proposal kernel $q(u,v)$. The ergodic behavior of the algorithm
is heavily dependent on the choice of proposal.

  * •

The accept-reject step may be implemented by drawing, independently from the
proposal, a uniformly distributed random variable $\theta_{n}$ in the interval
$[0,1].$ If $\theta_{n}\in[0,a_{n})$ then the proposal is accepted
($u^{(n+1)}=v^{\star}$); it is rejected ($u^{(n+1)}=u^{(n)})$ otherwise.

####  6.3 Invariance of the Target Distribution $\pi$

In this section we show that the Metropolis-Hastings kernel is invariant with
respect to the target $\pi.$ We start by introducing formal definitions of
detailed balance and invariance of a Markov kernel, and showing that the
former implies the latter. We then prove that the Metropolis-Hastings kernel
satisfies detailed balance with respect to $\pi$, and hence it is invariant.

#####  6.3.1 Detailed Balance and its Implication

The following definition makes precise the notion of detailed balance and
invariant distribution of a Markov kernel.

######  Definition 6.2.

A Markov kernel $p(u,v)$ satisfies detailed balance with respect to $\pi$ if,
for any $u,v\in\mathbb{R}^{d},$

$\pi(u)p(u,v)=\pi(v)p(v,u).$

We say that $\pi$ is an invariant distribution of the Markov kernel $p(u,v)$
if, for any $v\in\mathbb{R}^{d},$

$\int_{\mathbb{R}^{d}}\pi(u)p(u,v)\,du=\pi(v).$ (6.3)

Detailed balance of $p(u,v)$ with respect to $\pi$ implies that $\pi$ is
invariant for $p(u,v)$. To see this, note that if $p(u,v)$ satisfies detailed
balance with respect to $\pi$ then

$\int_{\mathbb{R}^{d}}\pi(u)p(u,v)\,du=\pi(v)\int_{\mathbb{R}^{d}}p(v,u)\,du=\pi(v).$

Invariance guarantees that, if the chain is distributed according to $\pi$ at
a given step, then it will also be distributed according to $\pi$ in the
following step. Detailed balance guarantees that the in/out flow between any
two states is preserved, which is a stronger condition. The concept of
detailed balance is illustrated through the flow diagram shown in Figure 10.
The flow to/from any node must be the same under detailed balance, so the
density at each node must remain constant at each step, as shown in Figure
10(b); Figure 10(a) fails to satisfy this.

#####  6.3.2 Detailed Balance and the Metropolis-Hastings Algorithm

Consider a generic Markov kernel $q(u,v)$ as shown in Figure 10(a), which does
not satisfy detailed balance, which we accept-reject according to the
Metropolis algorithm. Intuitively, this Metropolis-Hastings kernel, limits the
“flow” between nodes such that detailed balance is satisfied with respect to
the distribution $\pi$, as shown in Figure 10(b). There the red-highlighted
numbers represent where the transition probabilities of the kernel have been
reduced. The following theorem establishes the detailed balance of the
Metropolis-Hastings kernel with respect to the target $\pi$; it implies, as a
consequence, that the Metropolis-Hastings kernel is invariant with respect to
the target.

######  Theorem 6.3 (Metropolis-Hastings and Detailed Balance).

The Metropolis-Hastings kernel satisfies detailed balance with respect to the
distribution $\pi$.

###### Proof.

We need to show that, for any $u,v\in\mathbb{R}^{d},$

$\pi(u)p_{\mbox{\tiny{\rm MH}}}(u,v)=\pi(v)p_{\mbox{\tiny{\rm MH}}}(u,v).$

If $u=v$ the equality is trivial. Let us assume henceforth that $u\neq v,$ so
that by equation (6.2) we have

$p_{\mbox{\tiny{\rm MH}}}(u,v)=a(u,v)q(u,v).$

We rewrite

$\displaystyle p_{\mbox{\tiny{\rm MH}}}(u,v)$ $\displaystyle=\text{min
}\biggl{(}\frac{\pi(v)q(v,u)}{\pi(u)q(u,v)},1\biggr{)}~{}q(u,v)$
$\displaystyle=\frac{1}{\pi(u)}\times\text{min
}\Bigl{(}\pi(u)q(u,v),\pi(v)q(v,u)\Bigr{)}.$

Thus, invoking symmetry,

$\pi(u)p_{\mbox{\tiny{\rm MH}}}(u,v)=\text{min
}\bigl{(}\pi(u)q(u,v),\pi(v)q(v,u)\bigr{)}=\pi(v)p_{\mbox{\tiny{\rm
MH}}}(v,u).$

∎

The invariance of the Metropolis-Hastings $p_{\mbox{\tiny{\rm MH}}}$ with
respect to $\pi$ implies that if the initial sample is drawn from the target
($\pi_{0}=\pi)$, then all subsequent samples are also distributed according to
the target $(\pi_{n}=\pi).$

Figure 10: Toy representation of MCMC with 3 states ($u,v,w$). The numbers
associated with each term represent an example case. (a) The top chain
represents a Markov chain with Markov kernel $q$ which does not satisfy
detailed balance with respect to $\pi$. (b) The bottom chain utilizes Markov
kernel $p$, which satisfies detailed balance by modifying the transition
probabilities of $q$.

####  6.4 Convergence to the Target Distribution

In the previous section we showed that if we initialize the Metropolis-
Hastings algorithm with distribution $\pi$, all samples produced by the
algorithm will be distributed according to $\pi.$ But our original problem was
exactly that we were not able to sample from $\pi$. Our aim in this section is
to show that, for certain Metropolis-Hastings methods, the law $\pi_{n}$ of
the $n$-th sample converges to $\pi$ regardless of the initial distribution
$\pi_{0}.$ This is a strong form of ergodic behavior which does not hold in
general, as illustrated by the chain depicted in Figure 11.

Figure 11: The arrows represent transitions with probability one in a four
state Markov chain. The invariant distribution is the uniform distribution but
for $\pi_{0}=\delta_{A}$, we do not have ergodicity.

In order to understand the mechanisms behid ergodicity we will first consider
Markov chains with finite state space. We then study a specific Metropolis-
Hastings algorithm, known as the preconditioned Crank-Nicolson (pCN) method,
which applies to targets $\pi$ defined by their density with respect to a
Gaussian distribution.

#####  6.4.1 Finite State Space

We consider a Markov chain on the finite state space
$S=\mathopen{}\mathclose{{}\left\\{1,\cdots,d}\right\\}.$ We illustrate a
coupling approach to proving ergodicity and then, in the next subsection, use
it to study the pCN method on a continuous state space.

######  Theorem 6.4 (Ergodicity in Finite State Spaces).

Let $u_{n}$ be a Markov chain with state space $S$, transition kernel $p$, and
initial distribution $\pi_{0}$. Assume that

$\varepsilon\coloneqq\min_{i,j\in
S}p\mathopen{}\mathclose{{}\left(i,j}\right)>0.$ (6.4)

Then, $p$ has a unique invariant distribution $\pi$, and the following
convergence result holds:

$d_{\mbox{\tiny{\rm
TV}}}(\pi_{n},\pi)\leq\mathopen{}\mathclose{{}\left(1-\varepsilon}\right)^{n},$
(6.5)

where $\pi_{n}$ is the law of $u_{n}.$

###### Proof.

First note that the Markov kernel $p$ maps a probability distribution on $S$
into another probability distribution on $S$; it thus maps a compact, convex
set into itself. By Brouwer’s fixed point theorem it follows that $p$ has a
fixed point in this space, ensuring that an invariant distribution exists. We
will now show that for any invariant distribution $\pi$ equation (6.5) holds,
which then implies the uniqueness of the invariant distribution.

Let $\pi$ be an invariant distribution, a probability vector on $S$. Proving
convergence to equilibrium amounts to “forgetting the past”, to show that the
long time behavior of the Markov chain does not depend on the initial
distribution $\pi_{0}$ and in fact converges to $\pi$. In general, $u_{n+1}$
will be strongly dependent on $u_{n}$, but the condition given in (6.4)
implies that there is always some residual chance that the chain jumps to any
new state, at each step, independently of where it is currently located,
$u_{n}$. We will show that this residual probability of the chain to make a
“totally random” move will diminish the stochastic dependence on $u_{0}$ as
$n$ increases.

To formalize this idea, let $b_{n}$ be i.i.d. Bernoulli random variables with
$\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(b_{n}=1}\right)=\varepsilon$
and
$\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(b_{n}=0}\right)=1-\varepsilon;$
furthermore assume that the sequence $\\{b_{n}\\}$ is independent of the
randomness defining draws from $\\{p(u_{n},\cdot)\\}$. Because of the lower
bound on $p$ we may define a new Markov chain as follows:

$w_{n+1}\sim\begin{cases}&s\mathopen{}\mathclose{{}\left(w_{n},\cdot}\right),\hfill\text{
for }b_{n}=0,\\\
&r\mathopen{}\mathclose{{}\left(w_{n},\cdot}\right),\hfill\text{ for
}b_{n}=1,\end{cases}$ (6.6)

where

$s\mathopen{}\mathclose{{}\left(i,j}\right):=\frac{p\mathopen{}\mathclose{{}\left(i,j}\right)-\varepsilon
r\mathopen{}\mathclose{{}\left(i,j}\right)}{1-\varepsilon}$

and $r$ is the uniform transition kernel with equal probability of
transitioning to each state in $S:r(i,j)=d^{-1}$ for all $(i,j)\in S\times S.$
We may now compute

$\displaystyle\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(w_{n+1}=j|w_{n}=i}\right)$
$\displaystyle=\varepsilon\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(w_{n+1}=j|w_{n}=i,b_{n}=1}\right)+\mathopen{}\mathclose{{}\left(1-\varepsilon}\right)\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(w_{n+1}=j|w_{n}=i,b_{n}=0}\right)$
$\displaystyle=\varepsilon
r\mathopen{}\mathclose{{}\left(i,j}\right)+p\mathopen{}\mathclose{{}\left(i,j}\right)-\varepsilon
r\mathopen{}\mathclose{{}\left(i,j}\right)$
$\displaystyle=p\mathopen{}\mathclose{{}\left(i,j}\right).$

Thus the kernel defined by (6.6) is equivalent in law to $p$. However, by
introducing the ancillary random variables $b_{n}$, we have made explicit the
concept of “forgetting the past entirely, with a small probability” at every
step. We may now use this to complete the proof. Let $f:S\mapsto\mathbb{R}$ be
an arbitrary test function with $|f|_{\infty}\leq 1$ and
$\tau\coloneqq\min\mathopen{}\mathclose{{}\left(n\in\mathbb{N}:b_{n}=1}\right).$
Then, regardless of how $w_{n}$ is initialized,

$\displaystyle\operatorname{\mathbb{E}}\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(w_{n}}\right)}\right]$
$\displaystyle=\operatorname{\mathbb{E}}\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(w_{n}}\right)|\tau\geq
n}\right]\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(\tau\geq
n}\right)+\sum\limits_{l=0}^{n-1}\operatorname{\mathbb{E}}\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(w_{n}}\right)|\tau=l}\right]\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(\tau=l}\right)$
$\displaystyle=\underbrace{\operatorname{\mathbb{E}}\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(w_{n}}\right)|\tau\geq
n}\right]\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(\tau\geq
n}\right)}_{|\cdot|\leq\mathopen{}\mathclose{{}\left(1-\varepsilon}\right)^{n}}+\underbrace{\sum\limits_{l=0}^{n-1}\operatorname{\mathbb{E}}_{u_{0}\sim\mathfrak{u}\mathopen{}\mathclose{{}\left(\cdot}\right)}\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(w_{n-l}}\right)}\right]\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(\tau=l}\right)}_{\text{independent
of original initial distribution}},$

where $\mathfrak{u}$ denotes the uniform distribution on $S.$

Now consider two Markov chains $w_{n}$ and $w_{n}^{\prime}$ with kernel (6.6),
the first initialized from $\pi_{0}$ and the second from an invariant
distribution $\pi.$ The law of $w_{n}$ agrees with the law $\pi_{n}$ of the
original chain $u_{n},$ while for the second chain it follows from invariance
that $\pi_{n}^{\prime}=\pi.$ We will use the variational characterization of
the total variation distance established in Lemma 1.10. Employing the
preceding identity and noting that the contribution which is independent of
the initial distribution will cancel in the two different Markov chains, we
obtain

$\displaystyle d_{\mbox{\tiny{\rm
TV}}}(\pi_{n},\pi_{n}^{\prime})=\frac{1}{2}\sup\limits_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}^{\pi_{n}}\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(u}\right)}\right]-\operatorname{\mathbb{E}}^{\pi_{n}^{\prime}}\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(u}\right)}\right]}\right|\leq\mathopen{}\mathclose{{}\left(1-\varepsilon}\right)^{n}.$

Since $\pi_{n}^{\prime}=\pi$ the desired result follows. ∎

Before extending the above argument to a setting with continuous state space,
we make two remarks:

######  Remark 6.5.

The coupling proof we have just exhibited may be generalized in a number of
ways; in particular:

  * •

The distribution $r$ need not be uniform; it was only chosen so for
convenience. What is important is that $r_{i,j}$ is lower bounded,
independently of $i$, for all $j$. Adapting $r$ to the $p$ at hand, might in
some cases greatly improve the above bound – a larger $\epsilon$ might be
identified.

  * •

Convergence to equilibrium can also be shown if condition (6.4) holds with $p$
replaced by the $n$-step transition kernel $p^{n}$. Again, for some chains
this may yield faster bounds on the convergence to equilibrium.

#####  6.4.2 The pCN Method

The coupling argument used in the previous subsection for Markov chains with
finite state space may also be employed to study ergodicity of Markov chains
on a continuous state space. To illustrate this we consider a particular
Metropolis-Hastings algorithm, the pCN method, applied to a specific Bayesian
inverse problem. Before we get into the details of the inverse problem we
describe the idea behind the pCN method at a high level. The idea is this. If
the desired target distribution has the form

$\pi(u)=\frac{1}{Z}\tilde{g}(u)\tilde{\rho}(u)$ (6.7)

and if the Metropolis-Hastings proposal kernel $q$ satisfies detailed balance
with respect to $\tilde{\rho}$, then (6.1) simplifies to give

$\displaystyle
a(u,v)=\text{min}\Bigg{(}\frac{\tilde{\rho}(v)\tilde{g}(v)q(v,u)}{\tilde{\rho}(u)\
\tilde{g}(u)q(u,v)},1\Bigg{)}=\text{min}\Bigg{(}\frac{\tilde{g}(v)}{\tilde{g}(u)},1\Bigg{)}.$
(6.8)

We will apply and study this idea in the case where $\tilde{\rho}$ is a
Gaussian distribution, in which case it is straightforward to construct a
proposal kernel that satisfies detailed balance with respect to
$\tilde{\rho}$. This scenario arises naturally in Bayesian inverse problems
where the the prior is defined in terms of a Gaussian. We now formalize the
inverse problem setting that we consider by imposing certain assumptions on
the likelihood and the prior, and then relate both to the functions
$\tilde{g}$ and $\tilde{\rho}$ in Equation 6.7.

######  Assumption 6.6.

We make the following assumptions on the Bayesian inverse problem:

  * •

Bounded likelihood: there are $g^{-},g^{+}>0$ such that, for all
$u\in\mathbb{R}^{d},$

$0<g^{-}<g(u)<g^{+}.$

  * •

Truncated Gaussian prior: there is a compact set $B\subset\mathbb{R}^{d}$ with
positive Lebesgue measure such that $\rho(u)=\mathbb{1}_{B}(u)z(u),$ where
$z=\mathcal{N}\mathopen{}\mathclose{{}\left(0,\widehat{C}}\right)$.

Under Assumption 6.6 we obtain for the posterior density

$\pi\mathopen{}\mathclose{{}\left(u}\right)\propto
g(u)\mathbb{1}_{B}\mathopen{}\mathclose{{}\left(u}\right)z(u),$

which is of the form in Equation (6.7) with
$\tilde{g}(u)=g(u)\mathbb{1}_{B}(u)$ and $\tilde{\rho}(u)=z(u).$

The pCN method has proposal kernel which is invariant with respect to the
Gaussian $z=\mathcal{N}\mathopen{}\mathclose{{}\left(0,\widehat{C}}\right)$
defined by

$q(u,\cdot)\sim\mathcal{N}\Bigl{(}(1-\beta^{2})^{1/2}u,\beta^{2}\widehat{C}\Bigr{)}.$
(6.9)

Thus, given the sample $u^{(n)},$ the pCN proposes a new sample

$v^{\star}\sim\mathopen{}\mathclose{{}\left(1-\beta^{2}}\right)^{1/2}u^{(n)}+\beta\xi^{(n)},\quad\xi^{(n)}\sim\mathcal{N}\mathopen{}\mathclose{{}\left(0,\widehat{C}}\right),$

which only requires to sample a Gaussian.

######  Theorem 6.7 (Ergodicity for pCN Method).

Assume that we apply the pCN method to sample from a posterior density $\pi$
arising from Assumptions 6.6 with initial condition drawn from any density
supported on $B$. Then there exists a constant
$\varepsilon\in\mathopen{}\mathclose{{}\left(0,1}\right)$ such that

$d_{\mbox{\tiny{\rm
TV}}}(\pi_{n},\pi)\leq\mathopen{}\mathclose{{}\left(1-\varepsilon}\right)^{n},$

where $\pi_{n}$ is the law of the $n$-th sample from the pCN Metropolis-
Hastings algorithm.

Recall the notation for the covariance weighted inner-product and resulting
norm described in the introduction to these notes. Before proving the above
theorem, we will prove an auxiliary lemma.

######  Lemma 6.8.

Let the assumptions of Theorem 6.7 hold. The proposal kernel (6.9) satisfies
detailed balance with respect to the Gaussian
$z=\mathcal{N}\mathopen{}\mathclose{{}\left(0,\widehat{C}}\right)$ and we have
the following expression for the pCN accept-reject probability:

$\displaystyle
a\mathopen{}\mathclose{{}\left(u,v}\right)=\min\mathopen{}\mathclose{{}\left(\frac{g(v)}{g(u)}\mathbb{1}_{B}\mathopen{}\mathclose{{}\left(v}\right),1}\right).$

###### Proof.

We need to show that $z(v)q(v,u)$ is symmetric in $u$ and $v.$ By direct
calculation,

$\displaystyle-\log\bigl{(}z(v)q(v,u)\bigr{)}$
$\displaystyle=\frac{1}{2}|v|^{2}_{\widehat{C}}+\frac{1}{2\beta^{2}}\mathopen{}\mathclose{{}\left|u-\mathopen{}\mathclose{{}\left(1-\beta^{2}}\right)^{1/2}v}\right|_{\widehat{C}}^{2}$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{1}{2}+\frac{\mathopen{}\mathclose{{}\left(1-\beta^{2}}\right)}{2\beta^{2}}}\right)|v|^{2}_{\widehat{C}}+\frac{1}{2\beta^{2}}|u|^{2}_{\widehat{C}}-\frac{\mathopen{}\mathclose{{}\left(1-\beta^{2}}\right)^{1/2}}{\beta^{2}}\langle
u,v\rangle_{\widehat{C}}$
$\displaystyle=\frac{1}{2\beta^{2}}\mathopen{}\mathclose{{}\left(|v|^{2}_{\widehat{C}}+|u|^{2}_{\widehat{C}}}\right)-\frac{\mathopen{}\mathclose{{}\left(1-\beta^{2}}\right)^{1/2}}{\beta^{2}}\langle
u,v\rangle_{\widehat{C}}.$

For the expression for the acceptance probability note that since $u_{0}\in
B$, we have $u_{n}\in B$ for all $n$ as any proposed move out of $B$ will be
rejected. Thus $\mathbb{1}_{B}\mathopen{}\mathclose{{}\left(u}\right)$ may be
dropped from the formula for the acceptance probability. ∎

Using this lemma, we can prove ergodicity much as we did in the previous
subsection in the finite state space setting. The main idea is that,
restricted to the bounded set $B$, the probability density of the transition
kernel will be bounded away from zero by some $\varepsilon$. Splitting off a
“forgetful part” that is triggered with probability $\varepsilon$ will then
yield the result.

###### Proof of Theorem 6.7.

Note again that since $u_{0}\in B$ we have $u_{n}\in B$ for all $n\geq 1.$
Note further that since $B$ is compact and $q$ is continuous in both of its
arguments, there is $q^{-}>0$ such that, for any $u,v\in B,$

$q(u,v)\geq q^{-}.$

Let $p$ be the Markov kernel defined by the pCN Metropolis-Hastings algorithm.
It follows that, for $u,v\in B,$

$\displaystyle p(u,v)$ $\displaystyle\geq q(u,v)a(u,v)$ $\displaystyle\geq
q^{-}\frac{g^{-}}{g^{+}}=:\varepsilon\text{Leb}(B),$

where the last equation defines $\varepsilon$ and $\text{Leb}(B)$ denotes the
Lebesgue measure of $B$ (which is assumed to be positive). Analogously to the
discrete proof, we now define $b_{n}$ to be i.i.d. Bernoulli random variables
with
$\operatorname{\mathbb{P}}\mathopen{}\mathclose{{}\left(b_{n}=1}\right)=\varepsilon$,
independently of all other randomness, and consider the transition rule

$u_{n+1}\sim\begin{cases}&s\mathopen{}\mathclose{{}\left(u_{n},\cdot}\right),\hfill\text{
for }b_{n}=0,\\\ &r\mathopen{}\mathclose{{}\left(\cdot}\right),\hfill\text{
for }b_{n}=1,\end{cases}$

where $r$ denotes the uniform distribution on $B$ and, for $A\subset B$ and
$u\in B,$

$s\mathopen{}\mathclose{{}\left(u,A}\right)\coloneqq\frac{p\mathopen{}\mathclose{{}\left(u,A}\right)-\varepsilon
r\mathopen{}\mathclose{{}\left(A}\right)}{1-\varepsilon}.$

Just as in the discrete case, one can check that the resulting Markov kernel
is equal to the pCN Metropolis-Hastings kernel $p(\cdot,\cdot)$. Exponential
convergence is then concluded in exactly the same way as in the discrete case.
∎

####  6.5 Discussion and Bibliography

The book [40] is a useful basic introduction to MCMC and the book [16]
presents state of the art as of 2010. The paper [21] overviews the pCN method
and related MCMC algorithms specifically designed for inverse problems and
other sampling problems in high dimensional state spaces. The book [77]
describes the coupling method in a general setting. The book [86] contains a
wide-ranging presentation of Markov chains, and their long-time behavior,
including ergodicity, and coupling. Furthermore the book describes the general
methodology for going from convergence of expectations in (possibly weighted)
total variation metrics to sample path ergodicity and almost sure convergence
of time averages, a topic we did not cover in this chapter. The paper [85]
describes the coupling methodology in the context of stochastic differential
equations and their approximations.

###  7 Filtering and Smoothing Problems and Well-Posedness

In this chapter we introduce data assimilation problems in which the model of
interest, and the data associated with it, have a time-ordered nature. We
distinguish between the filtering problem (on-line) in which the data is
incorporated sequentially as it comes in, and the smoothing problem (off-line)
which is a specific instance of the inverse problems that have been the
subject of the preceding chapters.

####  7.1 Formulation of Filtering and Smoothing Problems

Consider the stochastic dynamics model given by

$\displaystyle v_{j+1}$
$\displaystyle=\Psi(v_{j})+\xi_{j},\>j\in\mathbb{Z}^{+},$ $\displaystyle
v_{0}$
$\displaystyle\sim\mathcal{N}(m_{0},C_{0}),\>\xi_{j}\sim\mathcal{N}(0,\Sigma)\>\text{i.i.d.},$

where we assume that $v_{0}$ is independent of the sequence $\\{\xi_{j}\\}$;
this is often written as $v_{0}\perp\\{\xi_{j}\\}.$ Now we add the data model
given by

$\displaystyle y_{j+1}$
$\displaystyle=h(v_{j+1})+\eta_{j+1},\>j\in\mathbb{Z}^{+},$
$\displaystyle\eta_{j}$
$\displaystyle\sim\mathcal{N}(0,\Gamma)\>\text{i.i.d.},$

where we assume that $\\{\eta_{j}\\}\perp v_{0}$ for all $j$ and
$\eta_{k}\perp\xi_{j}$ for all $j,k$. The following will be assumed in the
remainder of this book.

######  Assumption 7.1.

The matrices $C_{0},$ $\Sigma$ and $\Gamma$ are positive-definite. Further, we
have $\Psi\in C(\mathbb{R}^{d},\mathbb{R}^{d})$ and $h\in
C(\mathbb{R}^{d},\mathbb{R}^{k})$.

We define

$V:=\\{v_{0},\ldots,v_{J}\\},\>Y:=\\{y_{1},\ldots,y_{J}\\},\>\text{and}\>Y_{j}:=\\{y_{1},\ldots,y_{j}\\}.$

The sequence $V$ is often termed the signal and the sequence $Y$ the data.

######  Definition 7.2 (The Smoothing Problem).

The smoothing problem is to find the probability density
$\Pi(V):=\operatorname{\mathbb{P}}(V|Y)=\operatorname{\mathbb{P}}(\\{v_{0},\ldots,v_{J}\\}|\\{y_{1},\ldots,y_{J}\\})$
on $\mathbb{R}^{d(J+1)}$ for some fixed integer $J.$ We refer to $\Pi$ as the
smoothing distribution.

######  Definition 7.3 (The Filtering Problem).

The filtering problem is to find the probability densities
$\pi_{j}(v_{j}):=\operatorname{\mathbb{P}}(v_{j}|Y_{j})$ on $\mathbb{R}^{d}$
for $j=1,\dots,J.$ We refer to $\pi_{j}$ as the filtering distribution at time
$j.$

The key conceptual issue to appreciate concerning the filtering problem, in
comparison with the smoothing problem, is that interest is focused on
characterizing, or approximating, a sequence of probably distributions,
defined in an iterative fashion as the data is acquired sequentially.

######  Remark 7.4.

We note the following identity:

$\displaystyle\int\Pi(v_{0},\ldots,v_{J})dv_{0}dv_{1}\ldots
dv_{J-1}=\pi_{J}(v_{J}).$

This expresses the fact that the marginal of the smoothing distribution at
time $J$ corresponds to the filtering distribution at time $J$. Note also
that, in general, for $j<J$

$\int\Pi(v_{0},\ldots,v_{J})dv_{0}\ldots
dv_{j-1}\widehat{dv_{j}}dv_{j+1}\ldots dv_{J}\neq\pi_{j}(v_{j}),$

since the expression on the left-hand side of the equation depends on data
$Y_{J}$, whereas that on the right-hand side depends only on $Y_{j}$, and
$j<J$.

####  7.2 The Smoothing Problem

#####  7.2.1 Formula for pdf of the Smoothing Problem

The smoothing distribution can be found by combining a prior on $v$ and a
likelihood function using Bayes formula. The prior is the probability
distribution on $v$ implied by the distribution of $v_{0}$ and the stochastics
dynamics model; the likelihood function is defined by the data model. We now
derive the prior and the likelihood separately.

The prior distribution can be derived as follows:

$\displaystyle\operatorname{\mathbb{P}}(V)$
$\displaystyle=\operatorname{\mathbb{P}}(v_{J},v_{J-1},\ldots,v_{0})$
$\displaystyle=\operatorname{\mathbb{P}}(v_{J}|v_{J-1},\ldots,v_{0})\operatorname{\mathbb{P}}(v_{J-1},\ldots,v_{0})$
$\displaystyle=\operatorname{\mathbb{P}}(v_{J}|v_{J-1})\operatorname{\mathbb{P}}(v_{J-1},\ldots,v_{0}).$

The third equality comes from the Markov, or memoryless, property which
follows from the independence of the elements of the sequence $\\{\xi_{j}\\}.$
By induction, we have:

$\displaystyle\operatorname{\mathbb{P}}(V)$
$\displaystyle=\prod_{j=0}^{J-1}\operatorname{\mathbb{P}}(v_{j+1}|v_{j})\operatorname{\mathbb{P}}(v_{0})$
$\displaystyle=\frac{1}{Z_{\rho}}\exp\bigl{(}-\mathsf{R}(V)\bigr{)}$
$\displaystyle=:\rho(V),$

where $Z_{\rho}>0$ is a normalizing constant and

$\mathsf{R}(V):=\frac{1}{2}|v_{0}-m_{0}|^{2}_{C_{0}}+\frac{1}{2}\sum_{j=0}^{J-1}|v_{j+1}-\Psi(v_{j})|^{2}_{\Sigma}.$

The likelihood function, which incorporates the measurements gathered from
observing the system, depends only on the measurement model and may be derived
as follows:

$\displaystyle\operatorname{\mathbb{P}}(Y|V)$
$\displaystyle=\prod_{j=0}^{J-1}\operatorname{\mathbb{P}}(y_{j+1}|v_{0},\ldots,v_{J})$
$\displaystyle=\prod_{j=0}^{J-1}\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})$
$\displaystyle\propto\exp\bigl{(}-\mathsf{L}(V;Y)\bigr{)},$

where

$\mathsf{L}(V;Y):=\frac{1}{2}\displaystyle{\sum_{j=0}^{J-1}}|y_{j+1}-h(v_{j+1})|_{\Gamma}^{2}.$

The factorization of $\operatorname{\mathbb{P}}(Y|V)$ in terms of the product
of the $\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})$ follows from the
independence of the elements of $\\{\eta_{j}\\}$ and the fact that the
observation at time $j+1$ depends only on the state at time $j+1.$

Using Bayes Theorem 1.2 we find the smoothing distribution by combining the
likelihood and the prior

$\displaystyle\Pi(V)$
$\displaystyle\propto\operatorname{\mathbb{P}}(Y|V)\operatorname{\mathbb{P}}(V)$
$\displaystyle=\frac{1}{Z}\exp\bigl{(}-\mathsf{R}(V)-\mathsf{L}(V;Y)\bigr{)}.$

Note that $V\in\mathbb{R}^{d(J+1)}$, and $Y\in\mathbb{R}^{kJ}$.

#####  7.2.2 Well-Posedness of the Smoothing Problem

Now we study the well-posedness of the smoothing problem with respect to
perturbations in the data. To this end we consider two smoothing distributions
corresponding to different observed data sequences $Y,Y^{\prime}$:

$\displaystyle\Pi(V)$
$\displaystyle:=\operatorname{\mathbb{P}}(V|Y)=\frac{1}{Z}\exp\bigl{(}-\mathsf{R}(V)-\mathsf{L}(V;Y)\bigr{)},$
$\displaystyle\Pi^{\prime}(V)$
$\displaystyle:=\operatorname{\mathbb{P}}(V|Y^{\prime})=\frac{1}{Z^{\prime}}\exp\bigl{(}-\mathsf{R}(V)-\mathsf{L}(V;Y^{\prime})\bigr{)}.$

We make the following assumptions:

######  Assumption 7.5.

There is finite non-negative constant $R$ such that the data $Y,Y^{\prime}$
and the observation function $h$ satisfy:

  * •

$|Y|,|Y^{\prime}|\leq R;$

  * •

Letting $\varphi(V):=\Bigl{(}\sum_{j=1}^{J}(|h(v_{j})|^{2})\Bigr{)}^{1/2}$, it
holds that $\operatorname{\mathbb{E}}^{\rho}[\varphi(V)]<\infty$.

The following theorem shows well-posedness of the smoothing problem.

######  Theorem 7.6 (Well-posedness of Smoothing).

Under Assumption 7.5, there is $\kappa\in[0,\infty),$ independent of $Y$ and
$Y^{\prime}$ such that

$d_{\mbox{\tiny{\rm H}}}(\Pi,\Pi^{\prime})\leq\kappa|Y-Y^{\prime}|.$

###### Proof.

We show that the proof of Theorem 1.14, which established well-posedness for
Bayesian inverse problems under Assumption 1.12, applies in the smoothing
context as well. To do so we rewrite the problem in the same notation used in
Chapter 1, and show that Assumption 7.5 above implies Assumption 1.12. Write

$\displaystyle\Pi(V)$
$\displaystyle=\frac{1}{Z}\exp\bigl{(}-\mathsf{L}(V;Y)\bigr{)}\rho(V)=\frac{1}{Z}g(V;Y)\rho(V),$
$\displaystyle\Pi^{\prime}(V)$
$\displaystyle=\frac{1}{Z^{\prime}}\exp\bigl{(}-\mathsf{L}(V;Y^{\prime})\bigr{)}\rho(V)=\frac{1}{Z^{\prime}}g(V;Y^{\prime})\rho(V),$

where $Z,Z^{\prime}>0$ are normalization constants. Here $|Y-Y^{\prime}|$
plays the role of $\delta$ in Theorem 1.14. Since
$g(V;Y)=\exp\bigl{(}-\mathsf{L}(V;Y)\bigr{)}$ and $\mathsf{L}(V;Y)$ is
positive we have that

$\sup_{v}\Bigl{|}\sqrt{g(V;Y)}\Bigr{|}+\Bigl{|}\sqrt{g(V;Y^{\prime})}\Bigr{|}\leq
2,$

and so Assumption 1.12 (ii) is satisfied. To see that Assumption 1.12 (i) is
also satisfied, note that $e^{-x}$ is Lipschitz-1. Therefore, using the
Cauchy-Schwarz inequality and some algebraic manipulations, there is $\kappa$
independent of $Y$ and $Y^{\prime}$ such that

$\displaystyle\Bigl{|}\sqrt{g(V,Y)}-\sqrt{g(V;Y^{\prime})}\Bigr{|}$
$\displaystyle\leq\frac{1}{2}\Bigl{|}\mathsf{L}(V;Y)-\mathsf{L}(V;Y^{\prime})\Bigr{|}$
$\displaystyle=\frac{1}{2}\biggl{|}\sum_{j=0}^{J-1}\frac{1}{2}\Bigl{(}|y_{j+1}-h(v_{j+1})|_{\Gamma}^{2}-|y^{\prime}_{j+1}-h(v_{j+1})|_{\Gamma}^{2}\Bigr{)}\biggr{|}$
$\displaystyle\leq\kappa\sum_{j=0}^{J-1}|y_{j+1}-y^{\prime}_{j+1}|_{\Gamma}|y_{j+1}+y^{\prime}_{j+1}-2h(v_{j+1})|_{\Gamma}$
$\displaystyle\leq\kappa|Y-Y^{\prime}|\varphi(V),$

where $\varphi$ is defined in Assumption 7.5. This shows that under Assumption
7.5 the likelihood function of the smoothing problem satisfies Assumption 1.12
with $\delta=|Y-Y^{\prime}|$; Theorem 7.6 follows from Theorem 1.14. ∎

####  7.3 The Filtering Problem

#####  7.3.1 Formula for pdf of the Filtering Problem

Filtering consists on the iterative updating of a sequence of probability
distributions, as new data arrives. This rests on the formula

$\pi_{j+1}={\mathcal{A}}_{j}{\mathcal{P}}\pi_{j},$

where ${\mathcal{P}}$ is a linear Markov map and ${\mathcal{A}}_{j}$ is a
nonlinear likelihood map (Bayes theorem). We introduce
$\widehat{\pi}_{j+1}=\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})$ and recall that
${\pi}_{j}=\operatorname{\mathbb{P}}(v_{j}|Y_{j}).$ Then the maps
${\mathcal{P}}$ and ${\mathcal{A}}_{j}$ are defined by

The map ${\mathcal{P}}$ is sometimes termed prediction; the map
${\mathcal{A}}_{j}$ is termed analysis. Throughout this section, all integrals
are over $\mathbb{R}^{d}.$ First we derive the linear Markov map
${\mathcal{P}}$. By the Markov property of the stochastic dynamics model we
have

$\displaystyle\widehat{\pi}_{j+1}(v_{j+1})$
$\displaystyle=\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})$ (7.4)
$\displaystyle=\int\operatorname{\mathbb{P}}(v_{j+1}|Y_{j},v_{j})\operatorname{\mathbb{P}}(v_{j}|Y_{j})dv_{j}$
$\displaystyle=\int\operatorname{\mathbb{P}}(v_{j+1}|v_{j})\operatorname{\mathbb{P}}(v_{j}|Y_{j})dv_{j}$
$\displaystyle=\int\operatorname{\mathbb{P}}(v_{j+1}|v_{j})\pi_{j}(v_{j})dv_{j}$
$\displaystyle=\frac{1}{(2\pi)^{d/2}({\rm
det}\Sigma)^{1/2}}\int\exp\Bigl{(}-\frac{1}{2}|v_{j+1}-\Psi(v_{j})|^{2}_{\Sigma}\Bigr{)}\pi_{j}(v_{j})dv_{j}.$

This defines the linear integral operator ${\mathcal{P}}$.

Now we derive the nonlinear likelihood map ${\mathcal{A}}_{j}$ by using Bayes
theorem:

$\pi_{j+1}={\mathcal{A}}_{j}\widehat{\pi}_{j+1}.$

We note that

$\displaystyle\pi_{j+1}(v_{j+1})$
$\displaystyle=\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})$ (7.5)
$\displaystyle=\operatorname{\mathbb{P}}(v_{j+1}|Y_{j},y_{j+1})$
$\displaystyle=\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1},Y_{j})\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})}{\operatorname{\mathbb{P}}(y_{j+1}|Y_{j})}$
$\displaystyle=\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})}{\operatorname{\mathbb{P}}(y_{j+1}|Y_{j})}$
$\displaystyle=\frac{\exp(-\frac{1}{2}|y_{j+1}-h(v_{j+1})|^{2}_{\Gamma})\widehat{\pi}_{j+1}(v_{j+1})}{\int\exp(-\frac{1}{2}|y_{j+1}-h(v_{j+1})|^{2}_{\Gamma})\widehat{\pi}_{j+1}(v_{j+1})dv_{j+1}}.$

This defines the nonlinear map ${\mathcal{A}}_{j+1}$ through multiplication by
the likelihood, and then normalization to a probability measure.

#####  7.3.2 Well-Posedness of the Filtering Problem

Now we establish the well-posedness of the filtering problem. We let

$\pi_{J}=\operatorname{\mathbb{P}}(v_{J}|Y),\quad\quad\pi_{J}^{\prime}=\operatorname{\mathbb{P}}(v_{J}|Y^{\prime})$

be two filtering distributions arising from observed data $Y=Y_{J}$ and
$Y^{\prime}=Y_{J}^{\prime}.$ As noted in Remark 7.4, the filtering
distribution at time $J$ is the $J$-th marginal of the smoothing distribution;
using this observation, the well-posedness of the filtering problem is a
direct consequence of the well-posedness of the smoothing problem.

######  Corollary 7.7 (Well-posedness of Filtering).

Under Assumption 7.5, there exists $\kappa=\kappa(R)$, such that
$d_{\mbox{\tiny{\rm TV}}}(\pi_{J},\pi_{J}^{\prime})\leq\kappa|Y-Y^{\prime}|$.

###### Proof.

Let $\Pi,\Pi^{\prime}$ be the posterior distributions from the smoothing
problems $\Pi=\operatorname{\mathbb{P}}(V|Y)$ and
$\Pi^{\prime}=\operatorname{\mathbb{P}}(V|Y^{\prime})$. We note that there
exists $\kappa$ such that $d_{\mbox{\tiny{\rm
TV}}}(\Pi,\Pi^{\prime})\leq\kappa|Y-Y^{\prime}|$ by Theorem 7.6 and by the
fact that the Hellinger metric bounds the total variation metric (Lemma 1.9).
Let $f:\mathbb{R}^{d}\to\mathbb{R}$ and
$F:\mathbb{R}^{d(J+1)}\mapsto\mathbb{R}.$ Then

$\displaystyle d_{\mbox{\tiny{\rm TV}}}(\pi_{J},\pi_{J}^{\prime})$
$\displaystyle=\frac{1}{2}\sup_{|f|_{\infty}\leq
1}\Bigl{|}\operatorname{\mathbb{E}}^{\pi_{J}}[f(v_{J})]-\operatorname{\mathbb{E}}^{\pi_{J^{\prime}}}[f(v_{J})]\Bigr{|}$
$\displaystyle=\frac{1}{2}\sup_{|f|_{\infty}\leq
1}\Bigl{|}\operatorname{\mathbb{E}}^{\Pi}[f(v_{J})]-\operatorname{\mathbb{E}}^{\Pi^{\prime}}[f(v_{J})]\Bigr{|}$
$\displaystyle\leq\frac{1}{2}\sup_{|F|_{\infty}\leq
1}\Bigl{|}\operatorname{\mathbb{E}}^{\Pi}[F(V)]-\operatorname{\mathbb{E}}^{\Pi^{\prime}}[F(V)]\Bigr{|}$
$\displaystyle=d_{\mbox{\tiny{\rm TV}}}(\Pi,\Pi^{\prime})$
$\displaystyle\leq\kappa|Y-Y^{\prime}|.$

Here the first inequality follows from the fact that $\\{|f|\leq 1\\}$ can be
viewed as a subset of $\\{|F|\leq 1\\}$. ∎

####  7.4 Discussion and Bibliography

The book [73] gives a mathematical introduction to data assimilation; for
further information on the smoothing problem as presented here, see section
2.3 in that book; for further information on the filtering problem as
presented here, see section 2.4. The books [1, 99] give alternative
foundational presentations of the subject. The books [66, 91, 82, 17] study
data assimilation in the context of weather forecasting, oil reservoir
simulation, turbulence modeling and geophysical sciences, respectively.

###  8 The Linear-Gaussian Setting

Recall the stochastic dynamics and data models introduced in the previous
chapter:

$\displaystyle\begin{split}v_{j+1}&=\Psi(v_{j})+\xi_{j},\quad\quad\,\,\,\,\;\xi_{j}\sim\mathcal{N}(0,\Sigma)\text{
i.i.d.},\\\
y_{j+1}&=h(v_{j+1})+\eta_{j},\quad\quad\,\eta_{j}\sim\mathcal{N}(0,\Gamma)\text{
i.i.d.},\end{split}$ (8.1)

with $v_{0}\sim\mathcal{N}(m_{0},C_{0}),$ $C_{0},\Sigma$ and $\Gamma$ positive
definite and $v_{0}\perp\\{\xi_{j}\\}\perp\\{\eta_{j}\\}.$ Here we study the
filtering and smoothing problems under the assumption that both the state-
transition operator $\Psi(\cdot)$ and the observation operator $h(\cdot)$ are
linear. Throughout, we will assume the following:

######  Assumption 8.1.

The stochastic dynamics and the data models defined by equation (8.1) hold.
Moreover,

1

Linear dynamics: $v_{j+1}=Mv_{j}+\xi_{j}$ for some ${M}\in\mathbb{R}^{d\times
d}$;

2

Linear obervation: $y_{j+1}=Hv_{j+1}+\eta_{j+1}$ for some
$H\in\mathbb{R}^{k\times d}.$

We will be mostly concerned with the case where $d>k.$ Under the _linear-
Gaussian_ assumption, the filtering and smoothing distributions are Gaussian
and therefore are fully characterized by their mean and covariance. We
consider first the _Kalman filter_ , which gives explicit formulae for the
iterative update of the mean and covariance of the filtering distribution, and
then the _Kalman smoother_ , which characterizes the smoothing distribution.
While the Kalman filter and the Kalman smoother only characterize the
filtering and smoothing distributions in the linear-Gaussian setting, their
importance extends beyond this setting, as will be demonstrated in the next
two chapters.

####  8.1 Kalman Filter

The filtering problem is to estimate the state at time $j$ given the data from
the past up to the present time $j$. That is, we want to determine the pdf
$\pi_{j}=\operatorname{\mathbb{P}}(v_{j}|Y_{j})$, where
$Y_{j}:=\\{y_{1},\ldots,y_{j}\\}.$ We define
$\widehat{\pi}_{j+1}=\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})$ and recall the
evolution

$\pi_{j+1}={\mathcal{A}}_{j}{\mathcal{P}}\pi_{j},\quad\pi_{0}=\mathcal{N}(m_{0},C_{0}),$

which can be decomposed in terms of the prediction and analysis steps
(LABEL:eq:pna). Note that ${\mathcal{P}}$ does not depend on $j$ because the
same Markov process defined by the state dynamics governs each prediction
step, whereas ${\mathcal{A}}_{j}$ depends on $j$ because at each step $j$ the
likelihood sees different data. The linear dynamics assumption implies that
applying the operator ${\mathcal{P}}$ to a Gaussian distribution gives again a
Gaussian, and the linear observation assumption implies that applying the
operator ${\mathcal{A}}_{j}$ to a Gaussian gives again a Gaussian. Thefore we
have the following:

######  Theorem 8.2.

(Gaussianity of Filtering Distributions) Under Assumption 8.1, $\pi_{0}$,
$\\{\pi_{j+1}\\}_{j\in\mathbb{Z}^{+}}$ and
$\\{\widehat{\pi}_{j+1}\\}_{j\in\mathbb{Z}^{+}}$ are all Gaussian
distributions.

As a consequence, the filtering distributions can be entirely characterized by
their mean and covariance. We write

$\displaystyle\widehat{\pi}_{j+1}$
$\displaystyle=\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})=\mathcal{N}(\widehat{m}_{j+1},\widehat{C}_{j+1}),$
(prediction) $\displaystyle\pi_{j+1}$
$\displaystyle=\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})=\mathcal{N}(m_{j+1},C_{j+1}),$
(analysis)

and aim to find update formulae for these means and covariances. The Kalman
filter algorithm achieves this.

######  Theorem 8.3 (Kalman filter).

Under Assumption 8.1 for all $j\in\mathbb{Z}^{+}$, $C_{j}$ is positive
definite and

$\displaystyle\widehat{m}_{j+1}$ $\displaystyle=Mm_{j},$
$\displaystyle\widehat{C}_{j+1}$ $\displaystyle=MC_{j}M^{T}+\Sigma,$
$\displaystyle C_{j+1}^{-1}$
$\displaystyle=(MC_{j}M^{T}+\Sigma)^{-1}+H^{T}\Gamma^{-1}H,$ $\displaystyle
C_{j+1}^{-1}m_{j+1}$
$\displaystyle=(MC_{j}M^{T}+\Sigma)^{-1}Mm_{j}+H^{T}\Gamma^{-1}y_{j+1}.$

###### Proof.

The proof proceeds by breaking the Kalman filter step above into the
prediction and the analysis steps.

Prediction: The mean and variance of the prediction step may be calculated as
follows. The mean is given by:

$\displaystyle\widehat{m}_{j+1}$
$\displaystyle=\operatorname{\mathbb{E}}\big{[}v_{j+1}|Y_{j}\big{]}$
$\displaystyle=\operatorname{\mathbb{E}}\big{[}Mv_{j}+\xi_{j}|Y_{j}\big{]}$
$\displaystyle=M\operatorname{\mathbb{E}}\big{[}v_{j}|Y_{j}\big{]}+\operatorname{\mathbb{E}}\big{[}\xi_{j}|Y_{j}\big{]}$
$\displaystyle=Mm_{j},$

where we used that $\xi_{j}$ and $Y_{j}$ are independent. The covariance is
given by:

$\displaystyle\widehat{C}_{j+1}=$
$\displaystyle\operatorname{\mathbb{E}}\big{[}(v_{j+1}-\widehat{m}_{j+1})\otimes(v_{j+1}-\widehat{m}_{j+1})|Y_{j}\big{]}$
$\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\big{[}M(v_{j}-m_{j})\otimes
M(v_{j}-m_{j})|Y_{j}\big{]}+\operatorname{\mathbb{E}}\big{[}\xi_{j}\otimes\xi_{j}|Y_{j}\big{]}$
$\displaystyle+\operatorname{\mathbb{E}}\big{[}\xi_{j}\otimes
M(v_{j}-m_{j})|Y_{j}\big{]}+\operatorname{\mathbb{E}}\big{[}M(v_{j}-m_{j})\otimes\xi_{j}|Y_{j}\big{]}$
$\displaystyle=$ $\displaystyle
M\operatorname{\mathbb{E}}\big{[}(v_{j}-m_{j})\otimes(v_{j}-m_{j})|Y_{j}\big{]}M^{T}+\Sigma$
$\displaystyle=$ $\displaystyle MC_{j}M^{T}+\Sigma,$

where we used that $\xi_{j}$ and $v_{j}$ are independent. Thus in the linear
Gaussian setting the prediction operator ${\mathcal{P}}$ from
$\pi_{j}=\mathcal{N}(m_{j},C_{j})$ to
$\widehat{\pi}_{j+1}=\mathcal{N}(\widehat{m}_{j+1},\widehat{C}_{j+1})$ is
given by

$\displaystyle\widehat{m}_{j+1}$ $\displaystyle=Mm_{j},$
$\displaystyle\widehat{C}_{j+1}$ $\displaystyle=MC_{j}M^{T}+\Sigma.$

Analysis: The analysis step may be derived as follows, using Bayes Theorem:

$\displaystyle\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})$
$\displaystyle=\operatorname{\mathbb{P}}(v_{j+1}|y_{j+1},Y_{j})$
$\displaystyle\propto\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1},Y_{j})\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})$
$\displaystyle=\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|Y_{j}).$

This gives

$\displaystyle\begin{split}\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})&\propto\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|v_{j+1}-m_{j+1}|_{C_{j+1}}^{2}}\right)\\\
&\propto\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|y_{j+1}-Hv_{j+1}|_{\Gamma}^{2}}\right)\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|v_{j+1}-\widehat{m}_{j+1}|_{\widehat{C}_{j+1}}^{2}}\right)\\\
&=\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|y_{j+1}-Hv_{j+1}|_{\Gamma}^{2}-\frac{1}{2}|v_{j+1}-\widehat{m}_{j+1}|_{\widehat{C}_{j+1}}^{2}}\right).\end{split}$
(8.2)

Taking logartihms and matching quadratic and linear terms in $v_{j+1}$ from
either side of this identity gives the update operator ${\mathcal{A}}_{j}$
from $\widehat{\pi}_{j+1}=\mathcal{N}(\widehat{m}_{j+1},\widehat{C}_{j+1})$ to
$\pi_{j+1}=\mathcal{N}(m_{j+1},C_{j+1})$:

$\displaystyle C_{j+1}^{-1}$
$\displaystyle=\widehat{C}_{j+1}^{-1}+H^{T}\Gamma^{-1}H,$ $\displaystyle
C_{j+1}^{-1}m_{j+1}$
$\displaystyle=\widehat{C}_{j+1}^{-1}\widehat{m}_{j+1}+H^{T}\Gamma^{-1}y_{j+1}.$

Combining the prediction operator ${\mathcal{P}}$ and update operator
${\mathcal{A}}_{j}$ yields the desired update formulae.

Positive-definiteness of $C_{j}$: It remains to show that $C_{j}>0$ for all
$j\in\mathbb{Z}^{+}$. We will use induction. By assumption the result holds
true for $j=0$. Assume that it is true for $C_{j}$. For the prediction
operator ${\mathcal{P}}$ we have, for $u\neq 0$,

$\displaystyle\langle u,\widehat{C}_{j+1}u\rangle$ $\displaystyle=\langle
u,MC_{j}M^{T}u\rangle+\langle u,\Sigma u\rangle$ $\displaystyle=\langle
M^{T}u,C_{j}M^{T}u\rangle+\langle u,\Sigma u\rangle$ $\displaystyle\geq\langle
u,\Sigma u\rangle$ $\displaystyle>0,$

where we used that $C_{j}>0$ and $\Sigma>0.$ Therefore
$\widehat{C}_{j+1},\widehat{C}_{j+1}^{-1}>0$. Then for the update operator
${\mathcal{A}}_{j}$:

$\displaystyle\langle u,C_{j+1}^{-1}u\rangle$ $\displaystyle=\langle
u,\widehat{C}_{j+1}^{-1}u\rangle+\langle u,H^{T}\Gamma^{-1}Hu\rangle$
$\displaystyle=\langle u,\widehat{C}_{j+1}^{-1}u\rangle+\langle
Hu,\Gamma^{-1}Hu\rangle$ $\displaystyle\geq\langle
u,\widehat{C}_{j+1}^{-1}u\rangle$ $\displaystyle>0,$

where we used that $\Gamma>0.$ Therefore, $C_{j+1},C_{j+1}^{-1}>0$, which
concludes the proof. ∎

######  Remark 8.4.

The previous proof reveals two interesting facts about the structure of the
Kalman filter updates. The first is that the covariance update does not
involve the observed data; this can be thought of as a consequence of the fact
that the posterior covariance in the linear-Gaussian setting for inverse
problems does not depend on the observed data, as noted in Chapter 2. The
second is that the update formulae for the covariance are linear in the
prediction step, but nonlinear in the analysis step; specifically the analysis
step is linear in the precisions (inverse covariances).

#####  8.1.1 Kalman Filter: Alternative Formulation

We now rewrite the Kalman filter in an alternative form. This formulation is
written in terms of the covariances directly, and does not involve the
precisions. Furthermore, the formulation in Theorem 8.3 involves a matrix
inversion in the state space while the one that we will present here requires
only inversion the data space, namely $S_{j+1}^{-1}$ in what follows. In many
applications the observation space dimension is much smaller than the state
space dimension ($k\ll d$), and then the formulation given in this section is
much cheaper to compute than the one given in Theorem 8.3.

######  Corollary 8.5 (Standard Form of Kalman Filter).

Under Assumption 8.3, the Kalman update formulae may be written as

$\displaystyle\begin{split}m_{j+1}&=\widehat{m}_{j+1}+K_{j+1}d_{j+1},\\\
C_{j+1}&=(I-K_{j+1}H)\widehat{C}_{j+1},\end{split}$ (8.3)

where

$\displaystyle d_{j+1}$ $\displaystyle=y_{j+1}-H\widehat{m}_{j+1},$
$\displaystyle S_{j+1}$ $\displaystyle=H\widehat{C}_{j+1}H^{T}+\Gamma,$
$\displaystyle K_{j+1}$ $\displaystyle=\widehat{C}_{j+1}H^{T}S_{j+1}^{-1},$

and where $\widehat{m}_{j+1},\widehat{C}_{j+1}$ are defined as in the proof of
Theorem 8.3.

######  Remark 8.6.

The vector $d_{j+1}$ is known as the innovation and the matrix $K_{j+1}$ as
the Kalman gain. Note that $d_{j+1}$ measures the mismatch of the predicted
state from the given data.

The corollary may be proved by an application of the following:

######  Lemma 8.7.

(Woodbury Matrix Identity) Let $A\in\mathbb{R}^{p\times p}$,
$U\in\mathbb{R}^{p\times q}$, $B\in\mathbb{R}^{q\times q}$,
$V\in\mathbb{R}^{q\times p}$. If $A,B>0$, then $A+UBV$ is invertible and

$(A+UBV)^{-1}=A^{-1}-A^{-1}U(B^{-1}+VA^{-1}U)^{-1}VA^{-1}.$

Combining the form of $d_{j+1}$ and $\widehat{m}_{j+1}$ shows that the update
formula for the Kalman mean can be written as

$m_{j+1}=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1},\quad\widehat{m}_{j+1}=Mm_{j}.$
(8.4)

This update formula has the very natural interpretation that the mean update
is formed as a linear combination of the evolution of the noise-free dynamics
and of the data. Equations (8.3) and (8.4) show that the Kalman gain $K_{j+1}$
determines the weight given to the new observation $y_{j+1}$ in the state
estimation. The update formula (8.4) may be also derived from an optimization
perspective, the topic of the next subsection.

#####  8.1.2 Optimization Perspective: Mean of Kalman Filter

Since $\pi_{j+1}$ is Gaussian, its mean agrees with its mode. Thus, formulae
(8.2) implies that

$\displaystyle m_{j+1}$ $\displaystyle={\rm argmax}_{v}\pi_{j+1}(v)$
$\displaystyle={\rm argmin}_{v}{\mathsf{J}}(v),$

where

${\mathsf{J}}(v):=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{m}_{j+1}|_{\widehat{C}_{j+1}}^{2}.$

In other words, $m_{j+1}$ is chosen to fit both the observed data $y_{j+1}$
and the predictions $\widehat{m}_{j+1}$ as well as possible. The covariances
$\Gamma$ and $\widehat{C}_{j+1}$ determine the relative weighting between the
two quadratic terms. The solution of the minimization problem is given by
(8.4), as may be verified by direct differentiation of ${\mathsf{J}}.$

An alternative derivation which is helpful in more sophisticated contexts is
to cast the problem in terms of constrained minimization. Write
$v^{\prime}=v-\widehat{m}_{j+1}$, $y^{\prime}=y_{j+1}-H\widehat{m}_{j+1}$ and
$C^{\prime}=\widehat{C}_{j+1}.$ Then minimization of ${\mathsf{J}}$ may be
reformulated as

$\displaystyle m_{j+1}=\widehat{m}_{j+1}+{\rm
argmin}_{v^{\prime}}\mathopen{}\mathclose{{}\left(\frac{1}{2}|y^{\prime}-Hv^{\prime}|_{\Gamma}^{2}+\frac{1}{2}\langle
v^{\prime},b\rangle}\right),$

where the minimization is now subject to the constraint
$C^{\prime}b=v^{\prime}.$ Using Lagrange multipliers we write

${\mathsf{I}}(v^{\prime})=\frac{1}{2}|y^{\prime}-Hv^{\prime}|_{\Gamma}^{2}+\frac{1}{2}\langle
v^{\prime},b\rangle+\langle\lambda,C^{\prime}b-v^{\prime}\rangle;$ (8.5)

computing the dervative and setting to zero gives

$\displaystyle-H^{T}\Gamma^{-1}(y^{\prime}-Hv^{\prime})+\frac{1}{2}b-\lambda$
$\displaystyle=0,$ $\displaystyle\frac{1}{2}v^{\prime}+C^{\prime}\lambda$
$\displaystyle=0,$ $\displaystyle v^{\prime}-C^{\prime}b$ $\displaystyle=0.$

The last two equations imply that $C^{\prime}(2\lambda+b)=0.$ Thus we set
$\lambda=-\frac{1}{2}b$ and drop the second equation, replacing the first by

$-H^{T}\Gamma^{-1}(y^{\prime}-HC^{\prime}b)+b=0.$

Solving for $b$ gives

$\displaystyle v$ $\displaystyle=\widehat{m}_{j+1}+v^{\prime}$
$\displaystyle=\widehat{m}_{j+1}+C^{\prime}b$
$\displaystyle=\widehat{m}_{j+1}+C^{\prime}(H^{T}\Gamma^{-1}HC^{\prime}+I)^{-1}H^{T}\Gamma^{-1}y^{\prime}$
$\displaystyle=\widehat{m}_{j+1}+C^{\prime}(H^{T}\Gamma^{-1}HC^{\prime}+I)^{-1}H^{T}\Gamma^{-1}(y_{j+1}-H\widehat{m}_{j+1})$
$\displaystyle=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1},$

where we have defined

$K_{j+1}=C^{\prime}(H^{T}\Gamma^{-1}HC^{\prime}+I)^{-1}H^{T}\Gamma^{-1}.$

It remains to show that $K_{j+1}$ agrees with the definition given in
Corollary 8.5. To see this we note that if we choose $S$ to be any matrix
satisfying $K_{j+1}=C^{\prime}H^{T}S^{-1}$ then

$H^{T}S^{-1}=(H^{T}\Gamma^{-1}HC^{\prime}+I)^{-1}H^{T}\Gamma^{-1}$

so that

$(H^{T}\Gamma^{-1}HC^{\prime}+I)H^{T}=H^{T}\Gamma^{-1}S.$

Thus

$H^{T}\Gamma^{-1}HC^{\prime}H^{T}+H^{T}=H^{T}\Gamma^{-1}S$

which may be achieved by choosing any $S$ so that

$\Gamma^{-1}(HC^{\prime}H^{T}+\Gamma)=\Gamma^{-1}S$

and multiplication by $\Gamma$ gives the desired formula for $S$.

#####  8.1.3 Optimality of Kalman Filter

The following theorem states that the Kalman filter gives the best estimator
of the mean in an online setting. In the following $\operatorname{\mathbb{E}}$
denotes expectation with respect to all randomness present in the problem
statement, through the initial condition, the noisy dynamical evolution, and
the noisy data. Furthermore $\operatorname{\mathbb{E}}[\cdot|Y_{j}]$ denotes
conditional expectation, given the data $Y_{j}$ upto time $j$.

######  Theorem 8.8 (Optimality of Kalman Filter).

Let $\\{m_{j}\\}$ be the sequence computed using the Kalman filter, and
$\\{z_{j}\\}$ be any sequence in $\mathbb{R}^{d}$ such that $z_{j}$ is $Y_{j}$
measurable.111For practical purposes this means $z_{j}$ is a fixed non-random
function of given observed $Y_{j}.$ Then, for all $j\in\mathbb{N},$

$\operatorname{\mathbb{E}}\Big{[}|v_{j}-m_{j}|^{2}\mid
Y_{j}\Big{]}\leq\operatorname{\mathbb{E}}\Big{[}|v_{j}-z_{j}|^{2}\mid
Y_{j}\Big{]}.$

###### Proof.

Note that $m_{j}$ and $z_{j}$ are fixed and non-random, given $Y_{j}$. Thus we
have:

$\displaystyle\operatorname{\mathbb{E}}\Big{[}|v_{j}-z_{j}|^{2}\mid
Y_{j}\Big{]}=$
$\displaystyle\operatorname{\mathbb{E}}\Big{[}|v_{j}-m_{j}+m_{j}-z_{j}|^{2}\mid
Y_{j}\Big{]}$ $\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\Big{[}|v_{j}-m_{j}|^{2}+2\Big{\langle}v_{j}-m_{j},m_{j}-z_{j}\Big{\rangle}+|m_{j}-z_{j}|^{2}\mid
Y_{j}\Big{]}$ $\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\Big{[}|v_{j}-m_{j}|^{2}\mid
Y_{j}\Big{]}+2\Big{\langle}\operatorname{\mathbb{E}}\Big{[}v_{j}-m_{j}\mid
Y_{j}\Big{]},m_{j}-z_{j}\Big{\rangle}+|m_{j}-z_{j}|^{2}$ $\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\Big{[}|v_{j}-m_{j}|^{2}\mid
Y_{j}\Big{]}+2\Big{\langle}\operatorname{\mathbb{E}}\Big{[}v_{j}\mid
Y_{j}\Big{]}-m_{j},m_{j}-z_{j}\Big{\rangle}+|m_{j}-z_{j}|^{2}$
$\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\Big{[}|v_{j}-m_{j}|^{2}\mid
Y_{j}\Big{]}+0+|m_{j}-z_{j}|^{2}$ $\displaystyle\geq$
$\displaystyle\operatorname{\mathbb{E}}\Big{[}|v_{j}-m_{j}|^{2}\mid
Y_{j}\Big{]}.$

The fifth step follows since $m_{j}=\operatorname{\mathbb{E}}\big{[}v_{j}\mid
Y_{j}\big{]}.$ ∎

####  8.2 Kalman Smoother

We next discuss the Kalman smoother, which refers to the smoothing problem in
the linear-Gaussian setting. As with the Kalman filter, it is possible to
solve the problem explicitly because the smoothing distribution is itself a
Gaussian. The explicit formulae computed help to build intuition about the
smoothing distribution more generally. We recall remark 7.4, which implies
that the filtering distribution at time $j=J$ determines the marginal of the
Kalman smoother on its last coordinate. However, the filtering distributions
do not determine the Kalman smoother in its entirety.

#####  8.2.1 Defining Linear System

Let $V=\\{v_{0},\ldots,v_{J}\\}$ and $Y=\\{y_{1},\ldots,y_{J}\\}$. Using Bayes
theorem and the fact that $\\{\xi_{j}\\}$, $\\{\eta_{j}\\}$ are mutually
independent i.i.d. sequences, independent of $v_{0}$, we have

$\operatorname{\mathbb{P}}(V|Y)\propto\operatorname{\mathbb{P}}(Y|V)\operatorname{\mathbb{P}}(V)=\prod\limits_{j=1}^{J}\operatorname{\mathbb{P}}(y_{j}|v_{j})\times\prod\limits_{j=0}^{J-1}\operatorname{\mathbb{P}}(v_{j+1}|v_{j})\times\operatorname{\mathbb{P}}(v_{0}).$

Noting that

$v_{j+1}|v_{j}\sim\mathcal{N}(Mv_{j},\Sigma),\quad
y_{j}|v_{j}\sim\mathcal{N}(Hv_{j},\Gamma)$

the Kalman smoothing distribution can be expressed as

$\operatorname{\mathbb{P}}(V|Y)\propto\exp\bigl{(}-{\mathsf{J}}(V)\bigr{)},$
(8.7)

where

${\mathsf{J}}(V):=\frac{1}{2}|v_{0}-m_{0}|^{2}_{C_{0}}+\frac{1}{2}\sum\limits_{j=0}^{J-1}|v_{j+1}-Mv_{j}|^{2}_{\Sigma}+\frac{1}{2}\sum\limits_{j=0}^{J-1}|y_{j+1}-Hv_{j+1}|^{2}_{\Gamma}.$
(8.8)

######  Theorem 8.9 (Characterization of the Kalman Smoother).

Suppose that Assumption 8.1 holds. Then $\operatorname{\mathbb{P}}(V|Y)$ is
Gaussian with a block tridiagonal precision matrix $\Omega>0$ and mean $m$
solving $\Omega m=r$, where

$\displaystyle\Omega=$
$\displaystyle\begin{bmatrix}\Omega_{0,0}&\Omega_{0,1}&&&&\\\
\Omega_{1,0}&\Omega_{1,1}&...&&0&\\\ 0&...&...&&&\\\ &0&...&...&...&\\\
&&&...&\Omega_{J-1,J-1}&\Omega_{J-1,J}\\\
&&&&\Omega_{J,J-1}&\Omega_{J,J}\end{bmatrix}$ (8.9)

with

$\displaystyle\Omega_{0,0}$ $\displaystyle=C_{0}^{-1}+M^{T}\Sigma^{-1}M,$
$\displaystyle\Omega_{j,j}$
$\displaystyle=\Sigma^{-1}+M^{T}\Sigma^{-1}M+H^{T}\Gamma^{-1}H,\quad 1\leq
j\leq J-1,$ $\displaystyle\Omega_{J,J}$
$\displaystyle=\Sigma^{-1}+H^{T}\Gamma^{-1}H,$ $\displaystyle\Omega_{j,j+1}$
$\displaystyle=-M^{T}\Sigma^{-1},\quad 1\leq j\leq J-1,$ $\displaystyle r_{0}$
$\displaystyle=C_{0}^{-1}m_{0},$ $\displaystyle r_{j}$
$\displaystyle=H^{T}\Gamma^{-1}y_{j},\quad 1\leq j\leq J.$

###### Proof.

We may write ${\mathsf{J}}(V)=\frac{1}{2}|\Omega^{1/2}(V-m)|^{2}+q$ with $q$
independent of $V$, by definition. Note that $\Omega$ is then the Hessian of
${\mathsf{J}}(V)$, and differentiating in equation (8.8) we obtain that

$\displaystyle\Omega_{0,0}$
$\displaystyle=\partial_{v_{0}}^{2}{\mathsf{J}}(V)=C_{0}^{-1}+M^{T}\Sigma^{-1}M,$
$\displaystyle\Omega_{j,j}$
$\displaystyle=\partial_{v_{j}}^{2}{\mathsf{J}}(V)=\Sigma^{-1}+M^{T}\Sigma^{-1}M+H^{T}\Gamma^{-1}H,$
$\displaystyle\Omega_{J,J}$
$\displaystyle=\partial_{v_{J}}^{2}{\mathsf{J}}(V)=\Sigma^{-1}+H^{T}\Gamma^{-1}H,$
$\displaystyle\Omega_{j-1,j}$
$\displaystyle=\Omega_{j,j-1}=\partial_{v_{j},v_{j-1}}^{2}{\mathsf{J}}(V)=-M^{T}\Sigma^{-1}.$

Otherwise, for all other values of indices $\\{k,l\\}$, $\Omega_{k,l}=0.$ This
proves that the matrix $\Omega$ has a block tridiagonal structure.

Now we focus on finding $m$. We have that
$\nabla_{V}{\mathsf{J}}(V)=\Omega(V-m)$, so that
$-\nabla_{V}{\mathsf{J}}(V)|_{V=0}=\Omega m$. Thus, we find $r$ as

$\displaystyle r_{0}=$
$\displaystyle-\nabla_{v_{0}}{\mathsf{J}}(V)|_{V=0}=-(-C_{0}^{-1}m_{0})=C_{0}^{-1}m_{0},$
$\displaystyle r_{j}=$
$\displaystyle-\nabla_{v_{j}}{\mathsf{J}}(V)|_{V=0}=-(-H^{T}\Gamma^{-1}y_{j})=H^{T}\Gamma^{-1}y_{j}.$

We have shown that $\Omega$ is symmetric and that $\Omega\geq 0$; to prove
that $\Omega$ is a precision matrix, we need to show that $\Omega>0.$ Take,
for the sake of argument, $Y=0$ and $m_{0}=0$ in equation (8.8), so that every
term in the expansion of ${\mathsf{J}}(V)$ involves $V$. It is evident that in
such case ${\mathsf{J}}(V)=V^{T}\Omega V$. Suppose that $V^{T}\Omega V=0$ for
some nonzero $V$. Then by positive-definiteness of $C_{0},\Sigma$, and
$\Gamma$, it must be that $v_{0}=0$ and $v_{j+1}=Mv_{j}$ for $j=0,1,\dots,J$.
Thus we must have $V=0$. This proves that $\Omega$ is positive-definite. ∎

######  Remark 8.10.

We note once again that since the smoothing distribution in the linear-
Gaussian setting is itself Gaussian, its mean agrees with its mode. Therefore,
the posterior mean found above is the unique minimizer of ${\mathsf{J}}(V),$
that is, the MAP estimator.

#####  8.2.2 Solution of the Linear System

We may also obtain $m$ via Gaussian elimination, using the block tridiagonal
structure of $\Omega$, as follows. First we form the matrix sequence
$\\{\Omega_{j}\\}$:

$\displaystyle\Omega_{0}$ $\displaystyle=\Omega_{0,0},$
$\displaystyle\Omega_{j+1}$
$\displaystyle=\Omega_{j+1,j+1}-\Sigma^{-1}M^{T}\Omega_{j}^{-1}M^{T}\Sigma^{-1},\quad
j=0,\dots,J-1;$ (8.10)

and we form the vector sequence $\\{z_{j}\\}$:

$\displaystyle z_{0}$ $\displaystyle=C_{0}^{-1}m_{0},$ $\displaystyle z_{j+1}$
$\displaystyle=H^{T}\Gamma^{-1}y_{j+1}-\Sigma^{-1}M^{T}\Omega_{j}^{-1}z_{j}.$

We may then read off $m_{J}$ from the equation $\Omega_{J}m_{J}=z_{J}$ and
finally we perform back-substitution to obtain

$\Omega_{j}m_{j}=z_{j}-\Omega_{j,j+1}m_{j+1},\quad j=J-1,\ldots,1.$

Note that $m_{J}$ found this way coincides with the mean of the Kalman filter
at $j=J.$

######  Proposition 8.11.

The matrices $\\{\Omega_{j}\\}$ in (8.10) are positive definite.

###### Proof.

The proof of this theorem relies on the following two lemmas:

######  Lemma 8.12.

If

$X:=\begin{bmatrix}X_{1}&\times&\times&\times\\\ \times&X_{2}&\times&\times\\\
\times&\times&...&\times\\\ \times&\times&\times&X_{d}\end{bmatrix}$

is positive-definite symmetric then $X_{i}$ is positive-definite symmetric for
all $i\in\\{1,\dots,d\\}.$

######  Lemma 8.13.

Let $B$ be a block lower (or upper) triangular matrix with identity on the
diagonal. Then, $B$ is an invertible matrix.

Using Lemma 8.12, we deduce that $\Omega_{0}=\Omega_{0,0}$ is positive-
definite symmetric. Consider the matrix $B\in\mathbb{R}^{d(J+1)\times d(J+1)}$
defined as:

$B=\begin{bmatrix}I&0&&0\\\ -\Omega_{1,0}\Omega_{0}^{-1}&I&...&...\\\
0&&...&0\\\ &&0&I\end{bmatrix}.$

We compute

$B\Omega B^{T}=\begin{bmatrix}\Omega_{0}&0&&&0\\\
0&\Omega_{1}&\Omega_{1,2}&...&0\\\ &\Omega_{2,1}&\Omega_{2,2}&...&\\\ &&&&\\\
0&&&\Omega_{J-1,J-1}&\Omega_{J-1,J}\\\ 0&&&\Omega_{J,J-1}&\Omega_{J,J}\\\
\end{bmatrix}.$

By Lemma 8.12, the matrix

$\tilde{\Omega}=\begin{bmatrix}\Omega_{1}&\Omega_{1,2}&...&0\\\
\Omega_{2,1}&\Omega_{2,2}&...&\\\ &&&\\\ &&\Omega_{J-1,J-1}&\Omega_{J-1,J}\\\
&&\Omega_{J,J-1}&\Omega_{J,J}\\\ \end{bmatrix}$

is positive definite, and so is $\Omega_{0}.$

Lemma 8.12 and the positive definiteness of $\tilde{\Omega}$ imply that
$\Omega_{1}$ is positive definite. Therefore, by Lemma 8.13 the matrix

$B_{2}=\begin{bmatrix}I&0&&0\\\ -\Omega_{2,1}\Omega_{1}^{-1}&I&...&...\\\
&&...&0\\\ 0&&&I\\\ \end{bmatrix}$

is invertible. Thus we have

$B_{2}\tilde{\Omega}B_{2}^{T}=\begin{bmatrix}\Omega_{1}&0&&&0\\\
0&\Omega_{2}&\Omega_{2,3}&...&\\\ &\Omega_{3,2}&\Omega_{3,3}&&0\\\ &&&&\\\
&&&\Omega_{J-1,J-1}&\Omega_{J-1,J}\\\ 0&&&\Omega_{J,J-1}&\Omega_{J,J}\\\
\end{bmatrix},$

giving the positive definiteness of $\Omega_{2}$. Iterating the argument shows
that all the $\Omega_{j}$ are positive definite. ∎

####  8.3 Discussion and Bibliography

The original paper of Kalman, which is arguably the first systematic
presentation of a methodology to combine models with data, is [64]. The
continuous time analogue of that work may be found in [65]. The book [53]
overviews the subject in the context of time-series analysis and economics.
The proof of Corollary 8.5 may be found in [73]. The paper [101] contains an
application of the optimality property of the Kalman filter (which applies
beyond the linear Gaussian setting to the mean of the filtering distribution
in quite general settings).

The subject of Kalman smoothing is overviewed in the text [5]; see also [96].
A link between the standard implementation of the smoother and Gauss-Newton
methods for MAP estimation is made in [9]. Fur further details on the Kalman
smoother, in both discrete and continuous time, see [73] and [51].

###  9 Optimization for Filtering and Smoothing: 3DVAR and 4DVAR

In the previous chapter we showed how the mean of the Kalman filter could be
derived through an optimization principle, once the predictive covariance is
known. In this chapter we discuss two optimization based approaches to
filtering and smoothing, namely the 3DVAR and 4DVAR methodologies. We
emphasize that the methods we present in this chapter do not provide
approximations of the filtering and smoothing probability distributions; they
simply provide estimates of the signal, given data, in the filtering (on-line)
and smoothing (off-line) data scenarios.

####  9.1 The Setting

3DVAR borrows from the Kalman filter optimization principle, but substitutes a
fixed given covariance for the predictive covariance. Here “VAR” refers to
variational, and encodes the concept of optimization. The 3D and 4D,
respectively, refer to three Euclidean spatial dimensions and to three
Euclidean spatial dimensions plus a time dimension; this nomenclature reflects
the historical derivation of these problems in the geophysical sciences, but
the specific structure of fields over three dimensional Euclidean space plays
no role in the generalized form of the methods described here. The key
distinction is that 3DVAR solves a sequence of optimization problems at each
point in time (hence is an on-line filtering method); in contrast 4DVAR solves
an optimization problem which involves data distributed over time (and is an
off-line smoothing method).

Throughout we consider the set-up in which the dynamics model is nonlinear,
but the observation operator is linear, commonly occurring in applications. We
thus have a discrete-time dynamical system with noisy state transitions and
noisy observations given by

$\displaystyle\text{Dynamics Model:}\quad v_{j+1}$
$\displaystyle=\Psi(v_{j})+\xi_{j},\quad j\in\mathbb{Z}^{+}.$
$\displaystyle\text{Data Model:}\quad y_{j+1}$
$\displaystyle=Hv_{j+1}+\eta_{j+1},\quad j\in\mathbb{Z}^{+},\text{ for
some}\,\,H\in\mathbb{R}^{k\times d}.$ $\displaystyle\text{Probabilistic
Structure:}\quad v_{0}$
$\displaystyle\sim\mathcal{N}(m_{0},C_{0}),\quad\xi_{j}\sim\mathcal{N}(0,\Sigma),\quad\eta_{j}\sim\mathcal{N}(0,\Gamma).$
$\displaystyle\text{Probabilistic Structure:}\quad v_{0}$
$\displaystyle\perp\\{\xi_{j}\\}\perp\\{\eta_{j}\\}\text{ independent.}$

####  9.2 3DVAR

We introduce 3DVAR by analogy with the update formula (8.4) for the Kalman
filter, and its derivation through optimization from section 8.1.2. The
primary differences between 3DVAR and the Kalman filter mean update is that
$\Psi(\cdot)$ can be nonlinear for 3DVAR, and that for 3DVAR we have no closed
update formula for the covariances. To deal with this second issue 3DVAR uses
a fixed predicted covariance, independent of time $j$, and pre-specified. The
resulting minimization problem, and its solution, is described in Table 1,
making the analogy with the Kalman filter. Note that the minimization itself
is of a quadratic functional, and so may be solved by means of linear algebra.
The constraint formulation used for the Kalman filter, in section 8.1.2, may
also be applied and used to derive the mean update formula.

Kalman Filter | 3DVAR
---|---
$m_{j+1}=\arg\min_{v}{\mathsf{J}}(v)$ | $m_{j+1}=\arg\min_{v}{\mathsf{J}}(v)$
${\mathsf{J}}(v)=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{m}_{j+1}|_{\widehat{C}_{j+1}}^{2}$ | ${\mathsf{J}}(v)=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{m}_{j+1}|_{\widehat{C}}^{2}$
$\widehat{m}_{j+1}=Mm_{j}$ | $\widehat{m}_{j+1}=\Psi(m_{j})$
$m_{j+1}=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1}$ | $m_{j+1}=(I-KH)\widehat{m}_{j+1}+Ky_{j+1}$
Table 1: Comparison of Kalman Filter and 3DVAR update formulae.

The Kalman gain $K$ for 3DVAR is fixed, because the predicted covariance
$\widehat{C}$ is fixed. By analogy with Corollary 8.5 we have the following
formulae for the 3DVAR gain matrix $K$, and the update formula for the
estimator $m_{j}$:

$\displaystyle S$ $\displaystyle=H\widehat{C}H^{T}+\Gamma,$ $\displaystyle K$
$\displaystyle=\widehat{C}H^{T}S^{-1},$ $\displaystyle m_{j+1}$
$\displaystyle=(I-KH)\Psi({m}_{j})+Ky_{j+1}.$

The method also delivers an implied analysis covariance $C=(I-KH)\widehat{C}.$
Note that the resulting algorithm which maps $m_{j}$ to $m_{j+1}$ may be
specified directly in terms of the gain $K$, without need to introduce
$\widehat{C},C$ and $S$. In the remainder of this section we simply view $K$
as fixed and given. In this setting we show that the 3DVAR algorithm produces
accurate state estimation under vanishing noise assumptions in the
dynamics/data model. The governing assumptions concerning the dynamics/data
model are encapsulated in:

######  Assumption 9.1.

Consider the dynamics/data model under the assumptions that $\xi_{j}\equiv
0,\Gamma=\gamma^{2}\Gamma_{0},|\Gamma_{0}|=1$ and assume that the data
$y_{j+1}$ used in the 3DVAR algorithm is found from observing a true signal
$v_{j}^{\dagger}$ given by

$\displaystyle\emph{Dynamics Model:}\quad v_{j+1}^{\dagger}$
$\displaystyle=\Psi(v_{j}^{\dagger}),\quad j\in\mathbb{Z}^{+}.$
$\displaystyle\emph{Data Model:}\quad y_{j+1}$
$\displaystyle=Hv_{j+1}^{\dagger}+\gamma\eta^{\dagger}_{j+1,0},\quad
j\in\mathbb{Z}^{+}.$

With this assumption of noise-free dynamics ($\xi_{j}\equiv 0$) we deduce that
the 3DVAR filter produces output which, asymptotically, has an error of the
same size as the observational noise error $\gamma.$ The key additional
assumption in the theorem that allows this deduction is a relationship between
the Kalman gain $K$ and the derivative $D\Psi(\cdot)$ of the dynamics model.
Encoded in the assumption are two ingredients: that the observation operator
$H$ is rich enough in principle to learn enough components of the system to
synchronize the whole system; and that $K$ is designed cleverly enough to
effect this synchronization. The proof of the theorem is simply using these
two ingredients and then controlling the small stochastic perturbations,
arising from Assumption 9.1.

######  Theorem 9.2 (Accuracy of 3DVAR).

Let Assumption 9.1 hold with
$\eta^{\dagger}_{j,0}\sim\mathcal{N}(0,\Gamma_{0}).$ an i.i.d. sequence.
Assume that, for the gain matrix $K$ appearing in the 3DVAR method, there
exists a norm $\lVert\cdot\rVert$ on $\mathbb{R}^{d}$ and constant
$\lambda\in(0,1)$ such that, for all $v\in\mathbb{R}^{d},$

$\lVert(I-KH)D\Psi(v)\rVert\leq\lambda.$

Then there is constant $c>0$ such that the 3DVAR algorithm satisfies the
following large-time asymptotic error bound:

$\limsup_{j\to\infty}\operatorname{\mathbb{E}}[\lVert
m_{j}-v^{\dagger}_{j}\rVert]\leqslant\frac{c\gamma}{1-\lambda},$

where the expectation is taken with respect to the sequence
$\\{\eta^{\dagger}_{j,0}\\}.$

###### Proof.

We have

$\displaystyle v^{\dagger}_{j+1}$ $\displaystyle=\Psi(v^{\dagger}_{j}),$
$\displaystyle m_{j+1}$ $\displaystyle=(I-KH)\Psi(m_{j})+Ky_{j+1},$

and hence that

$\displaystyle v^{\dagger}_{j+1}$
$\displaystyle=(I-KH)\Psi(v^{\dagger}_{j})+KH\Psi(v^{\dagger}_{j}),$
$\displaystyle m_{j+1}$
$\displaystyle=(I-KH)\Psi(m_{j})+KH\Psi(v^{\dagger}_{j})+\gamma
K\eta^{\dagger}_{j+1,0}.$

Define $e_{j}=m_{j}-v_{j}^{\dagger}.$ By subtracting the evolution equation
for $v_{j}^{\dagger}$ from that for $m_{j}$ we obtain, using the mean value
theorem,

$\displaystyle e_{j+1}$ $\displaystyle=m_{j+1}-v^{\dagger}_{j+1}$
$\displaystyle=(I-KH)(\Psi(m_{j})-\Psi(v^{\dagger}_{j}))+\gamma
K\eta^{\dagger}_{j+1,0}$
$\displaystyle=\mathopen{}\mathclose{{}\left((I-KH)\int_{0}^{1}{D}\Psi\mathopen{}\mathclose{{}\left(sm_{j}+(1-s)v_{j}^{\dagger}}\right)ds}\right)e_{j}+\gamma
K\eta_{j+1,0}^{\dagger}.$

As a result, by the triangle inequality,

$\displaystyle\lVert e_{j+1}\rVert$
$\displaystyle\leq\mathopen{}\mathclose{{}\left\lVert\mathopen{}\mathclose{{}\left(\int_{0}^{1}(I-KH)D\Psi\mathopen{}\mathclose{{}\left(sm_{j}+(1-s)v_{j}^{\dagger}}\right)ds}\right)e_{j}}\right\rVert+\lVert\gamma
K\eta_{j+1,0}^{\dagger}\rVert$
$\displaystyle\leq\mathopen{}\mathclose{{}\left(\int_{0}^{1}\mathopen{}\mathclose{{}\left\lVert(I-KH)D\Psi\mathopen{}\mathclose{{}\left(sm_{j}+(1-s)v_{j}^{\dagger}}\right)}\right\rVert
ds}\right)\lVert e_{j}\rVert+\lVert\gamma K\eta_{j+1,0}^{\dagger}\rVert$
$\displaystyle\leq\lambda\lVert e_{j}\rVert+\gamma\lVert
K\eta_{j+1,0}^{\dagger}\rVert.$

Taking expectations on both sides, we obtain, for
$c:=\operatorname{\mathbb{E}}[\lVert K\eta_{0,j+1}^{\dagger}\rVert]>0$,

$\begin{split}\operatorname{\mathbb{E}}[\lVert
e_{j+1}\rVert]&\leq\lambda\operatorname{\mathbb{E}}[\lVert
e_{j}\rVert]+\gamma\operatorname{\mathbb{E}}[\lVert
K\eta_{j+1,0}^{\dagger}\rVert]\\\ &\leq\lambda\operatorname{\mathbb{E}}[\lVert
e_{j}\rVert]+\gamma c.\end{split}$ (9.1)

Using the discrete Gronwall inequality we have that:

$\begin{split}\operatorname{\mathbb{E}}[\lVert
e_{j}\rVert]&\leq\lambda^{j}\operatorname{\mathbb{E}}[\lVert
e_{0}\rVert]+\sum_{i=0}^{j-1}c\lambda^{i}\gamma\\\
&\leq\lambda^{j}\operatorname{\mathbb{E}}[\lVert
e_{0}\rVert]+c\gamma\frac{1-\lambda^{j}}{1-\lambda},\end{split}$ (9.2)

where $e_{0}=m_{0}-v_{0}$. As $\lambda<1$, the desired statement follows. ∎

####  9.3 4DVAR

Recall that 3DVAR differs from 4DVAR because, whilst also based on an
optimization principle, 4DVAR is applied in a distributed fashion over all
data in the time interval $j=1,\ldots,J$; in contrast 3DVAR is applied
sequentially from time $j-1$ to time $j,$ for $j=1,\ldots,J$. We consider two
forms of the methodology: weak constraint 4DVAR (w4DVAR), in which the fact
that the dynamics model contains randomness is allowed for in the
optimization; and 4DVAR (sometimes known as strong constraint 4DVAR) which can
be derived from w4DVAR in the limit of $\Sigma\to 0$ (no randomness in the
dynamics).

The objective function minimized in w4DVAR is

${\mathsf{J}}(V)=\frac{1}{2}|v_{0}-m_{0}|^{2}_{C_{0}}+\frac{1}{2}\sum^{J-1}_{j=0}|v_{j+1}-\Psi(v_{j})|^{2}_{\Sigma}+\frac{1}{2}\sum^{J-1}_{j=0}|y_{j+1}-H{v_{j+1}}|^{2}_{\Gamma},$
(9.3)

where $V=\\{v_{j}\\}^{J}_{j=0}\in\mathbb{R}^{d(J+1)}$,
$Y=\\{y_{j}\\}^{J}_{j=1}\in\mathbb{R}^{kJ}$, $v_{j}\in\mathbb{R}^{d}$,
$y_{j}\in\mathbb{R}^{k}$, $H$ is the observation operator, $\Sigma$ is the
random dynamical system covariance, $\Gamma$ is the data noise covariance and
$m_{0}$ and $C_{0}$ are the mean and covariance of the initial state. The
three terms in the objective function enforce, in turn, the initial condition
$v_{0}$, the dynamics model and the data model. Note that, because $\Psi$ is
nonlinear, the objective is not quadratic and cannot be optimized in closed
form. In contrast, each step of 3DVAR required solution of a quadratic
optimization, tractable in closed form.

######  Theorem 9.3 (Minimizer Exists for w4DVAR).

Assume that $\Psi$ is bounded and continuous. Then ${\mathsf{J}}$ has a
minimizer, which is a MAP estimator for the smoothing problem.

###### Proof.

Recall Theorem 3.5, which shows that the MAP estimator based on the smoothing
distribution
$\operatorname{\mathbb{P}}(V|Y)\propto\exp\bigl{(}-{\mathsf{J}}(V)\bigr{)}$ is
attained provided that ${\mathsf{J}}$ is guaranteed to be non-negative,
continuous and satisfy ${\mathsf{J}}(V)\to\infty$ as $|V|\to\infty.$ Now, the
objective ${\mathsf{J}}$ defined by equation (9.3) is clearly non-negative,
and it is continuous since $\Psi$ is assumed to be continuous. It remains to
show that ${\mathsf{J}}(V)\to\infty$ as $|V|\to\infty.$ Let $R$ be a bound for
$\Psi,$ so that $|\Psi(v_{j})|_{\Sigma}\leq R$ for all
$v_{j}\in\mathbb{R}^{d}.$ Then, since

${\mathsf{J}}(V)\geq\frac{1}{2}|v_{0}|_{C_{0}}^{2}-2|v_{0}|_{C_{0}}|m_{0}|_{C_{0}}+\frac{1}{2}\sum_{j=0}^{J-1}\Bigl{(}|v_{j+1}|^{2}_{\Sigma}-2R|v_{j+1}|_{\Sigma}\Bigr{)}$

it follows that ${\mathsf{J}}(V)\to\infty$ as $|V|\to\infty$ and the proof is
complete. ∎

We now consider the vanishing dynamical noise limit of w4DVAR. This is to
minimize

${\mathsf{J}}_{0}(V)=\frac{1}{2}|v_{0}-m_{0}|^{2}_{C_{0}}+\frac{1}{2}\sum^{J-1}_{j=0}|y_{j+1}-H{v_{j+1}}|^{2}_{\Gamma}$

subject to the hard constraint that

$v_{j+1}=\Psi(v_{j}),\quad j=0,\ldots,J-1.$

This is 4DVAR. Note that by using the constraint 4DVAR can be written as a
minimization over $v_{0}$, rather than over the entire sequence
$\\{v_{j}\\}_{j=0}^{J}$ as is required in w4DVAR.

We let ${\mathsf{J}}_{\sigma}$ denote the objective function ${\mathsf{J}}$
from w4DVAR in the case where $\Sigma\to\sigma^{2}\Sigma_{0}$. Roughly
speaking, the following result shows that minimizers of
${\mathsf{J}}_{\sigma}$ converge as $\sigma\to 0$ to points in
${\mathbb{R}}^{k(J+1)}$ which satisfy the hard constraint associated with
4DVAR.

######  Theorem 9.4 (Small Signal Noise Limit of w4DVAR).

Suppose that $\Psi$ is bounded and continuous and let $V^{\sigma}$ be a
minimizer of ${\mathsf{J}}_{\sigma}.$ Then as $\sigma\to 0$ there is a
convergent subsequence of $V^{\sigma}$ with limit $V^{*}$ satisfying
$v^{*}_{j+1}=\Psi(v^{*}_{j})$.

###### Proof.

Throughout this proof $c$ is a constant which may change from instance to
instance, but is independent of $\sigma.$ Consider $V\in\mathbb{R}^{d(J+1)}$
defined by $v_{0}=m_{0}$ and $v_{j+1}=\Psi(v_{j})$. Then $V$ is bounded, as
$\Psi(\cdot)$ is bounded, and the bound is independent of $\sigma$.
Furthermore

${\mathsf{J}}_{\sigma}(V)=\frac{1}{2}\sum^{J-1}_{j=0}|y_{j+1}-Hv_{j+1}|^{2}_{\Gamma}\leq
c,$

where $c$ is independent of $\sigma$. It follows that

${\mathsf{J}}_{\sigma}(V^{\sigma})\leq{\mathsf{J}}_{\sigma}(V)\leq c.$

Thus

$\displaystyle\frac{1}{2}|v^{\sigma}_{j+1}-\Psi(v^{\sigma}_{j})|^{2}_{\Sigma_{0}}$
$\displaystyle=\frac{\sigma^{2}}{2}|v^{\sigma}_{j+1}-\Psi(v^{\sigma}_{j})|^{2}_{\Sigma}=\sigma^{2}{\mathsf{J}}_{\sigma}(V^{\sigma})\leq\sigma^{2}c,$
$\displaystyle\frac{1}{2}|v^{\sigma}_{0}-m_{0}|^{2}_{C_{0}}$
$\displaystyle\leq{\mathsf{J}}_{\sigma}(V^{\sigma})\leq c.$

Since $\Psi$ is bounded these bounds imply that $|V^{\sigma}|\leq c$.
Therefore there is a limit $V^{*}:V^{\sigma}\to V^{*}$ along a subsequence. By
continuity

$0\leq\frac{1}{2}|v^{*}_{j+1}-\Psi(v^{*}_{j})|^{2}_{\Sigma_{0}}\leftarrow\frac{1}{2}|v^{\sigma}_{j+1}-\Psi(v^{\sigma}_{j})|^{2}_{\Sigma_{0}}\leq\sigma^{2}c$

Letting $\sigma\to 0$ we obtain that $v^{*}_{j+1}=\Psi(v^{*}_{j}).$ ∎

####  9.4 Discussion and Bibliography

The 3DVAR and 4DVAR methodologies, in the context of weather forecasting, are
discussed in [78] and [38] respectively. The accuracy analysis presented here
is similar to that which first appeared in the papers [14], [87] and was
developed further in [72, 74]. It arises from considering stochastic
perturbations of the seminal work of Titi and co-workers, exemplified by the
paper [54]. For an overview of variational methods, and their links to
problems in physics and mechanics, see the book [1], and the references
therein; see also the paper [15].

###  10 The Extended and Ensemble Kalman Filters

In this chapter we describe the Extended Kalman Filter (ExKF)222The extended
Kalman filter is often termed the EKF in the literature, a terminology
introduced before the existence of the EnKF; we find it useful to write ExKF
to uneqivocally distinguish it from the EnKF. and the Ensemble Kalman Filter
(EnKF). The ExKF approximates the predictive covariance by linearization,
while the EnKF approximates it by the empirical covariance of a collection of
particles. The status of the two methods in relation to the true filtering
distribution is as follows. The ExKF is a provably accurate approximation of
the filtering distribution in situations in which small noise is present in
both signal and data, and the filtering distribution is well-approximated by a
Gaussian. The EnKF is also in principle a good approximation of the filtering
distribution in this situation, if a large number of particles is used.
However the EnKF is typically deployed for problems where the use of a
sufficiently large number of particles is impractical; it is then better
viewed as an online optimizer, in the spirit of 3DVAR, but using multiple
particles to better estimate the covariances appearing in the quadratic
objective function which is minimized to find particle updates.

####  10.1 The Setting

Throughout this chapter we consider the setting in which 3DVAR was introduced
and may be applied: the dynamics model is nonlinear, but the observation
operator is linear. For purposes of exposition we summarize it again here:

$\displaystyle v_{j+1}$ $\displaystyle=\Psi(v_{j})+\xi_{j},$
$\displaystyle\xi_{j}\sim\mathcal{N}(0,\Sigma)\text{ i.i.d.},$ $\displaystyle
y_{j+1}$ $\displaystyle=Hv_{j+1}+\eta_{j+1},$
$\displaystyle\eta_{j}\sim\mathcal{N}(0,\Gamma)\text{ i.i.d.},$

with $v_{0}\sim\mathcal{N}(m_{0},C_{0})$ independent of the independent i.i.d.
sequences $\\{\xi_{j}\\}$ and $\\{\eta_{j}\\}$. Throughout this chapter we
assume that $v_{j}\in\mathbb{R}^{d},\,y_{j}\in\mathbb{R}^{k}.$

####  10.2 The Extended Kalman Filter

This method is derived by applying the Kalman methodology, using linearization
to propagate the covariance $C_{j}$ to the predictive covariance
$\widehat{C}_{j+1}.$ Table 2 summarizes the idea, after which we calculate the
formulae required in full detail.

Kalman Filter | ExKF
---|---
$m_{j+1}=\arg\min_{v}{\mathsf{J}}(v)$ | $m_{j+1}=\arg\min_{v}{\mathsf{J}}(v)$
${\mathsf{J}}(v)=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{m}_{j+1}|_{\widehat{C}_{j+1}}^{2}$ | ${\mathsf{J}}(v)=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{m}_{j+1}|_{\widehat{C}_{j+1}}^{2}$
$\widehat{m}_{j+1}=Mm_{j}$ | $\widehat{m}_{j+1}=\Psi(m_{j})$
$\widehat{C}_{j+1}$ update exact |  $\widehat{C}_{j+1}$ update by linearization
$m_{j+1}=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1}$ | $m_{j+1}=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1}$
Table 2: Comparison of Kalman Filter and ExKF update formulae

We first recall the Kalman filter update formulae and their derivation. We
have

$\displaystyle\widehat{v}_{j+1}=Mv_{j}+\xi_{j},\quad
v_{j}\sim\mathcal{N}(m_{j},C_{j}),\quad\xi_{j}\sim\mathcal{N}(0,\Sigma).$
(10.1)

From this we deduce, by taking expectations, that

$\widehat{m}_{j+1}=\operatorname{\mathbb{E}}[\widehat{v}_{j+1}\mid
Y_{j}]=\operatorname{\mathbb{E}}[Mv_{j}+\xi_{j}\mid
Y_{j}]=\operatorname{\mathbb{E}}[Mv_{j}\mid
Y_{j}]+\operatorname{\mathbb{E}}[\xi_{j}\mid Y_{j}]=Mm_{j}.$ (10.2)

The covariance update is derived as follows.

$\displaystyle\widehat{C}_{j+1}$
$\displaystyle=\operatorname{\mathbb{E}}\Bigl{[}(\widehat{v}_{j+1}-\widehat{m}_{j+1})\otimes(\widehat{v}_{j+1}-\widehat{m}_{j+1})\mid
Y_{j}\Bigr{]}$
$\displaystyle=\operatorname{\mathbb{E}}\Bigl{[}(M(v_{j}-m_{j})+\xi_{j})\otimes(M(v_{j}-m_{j})+\xi_{j})\mid
Y_{j}\Bigr{]}$
$\displaystyle=\operatorname{\mathbb{E}}\Bigl{[}(M(v_{j}-m_{j}))\otimes(M(v_{j}-m_{j}))\mid
Y_{j}\Bigr{]}+\operatorname{\mathbb{E}}\Bigl{[}\xi_{j}\otimes\xi_{j}\mid
Y_{j}\Bigr{]}$ (10.3)
$\displaystyle\quad+\operatorname{\mathbb{E}}\Bigl{[}(M(v_{j}-m_{j}))\otimes\xi_{j}\mid
Y_{j}\Bigr{]}+\operatorname{\mathbb{E}}\Bigl{[}\xi_{j}\otimes(M(v_{j}-m_{j}))\mid
Y_{j}\Bigr{]}$
$\displaystyle=M\operatorname{\mathbb{E}}\Bigl{[}(v_{j}-m_{j})\otimes(v_{j}-m_{j})\mid
Y_{j}\Bigr{]}M^{T}+\Sigma$ $\displaystyle=MC_{j}M^{T}+\Sigma.$

For the ExKF, the prediction map $\Psi$ is no longer linear. But since
$\xi_{j}$ is independent of $Y_{j}$ and $v_{j}$, we obtain

$\widehat{m}_{j+1}=\operatorname{\mathbb{E}}\Bigl{[}\Psi(v_{j})+\xi_{j}\mid
Y_{j}\Bigr{]}=\operatorname{\mathbb{E}}\Bigl{[}\Psi(v_{j})\mid
Y_{j}\Bigr{]}+\operatorname{\mathbb{E}}\Bigl{[}\xi_{j}\mid
Y_{j}\Bigr{]}=\operatorname{\mathbb{E}}\Bigl{[}\Psi(v_{j})\mid Y_{j}\Bigr{]}.$

If we assume that the fluctuations of $v_{j}$ around its mean $m_{j}$
(conditional on data) are small then a reasonable approximation is to take
$\Psi(v_{j})\approx\Psi(m_{j})$ so that

$\widehat{m}_{j+1}=\Psi(m_{j}).$ (10.4)

For the predictive covariance we use linearization; we have

$\displaystyle\widehat{C}_{j+1}$
$\displaystyle=\operatorname{\mathbb{E}}\Bigl{[}(\widehat{v}_{j+1}-\widehat{m}_{j+1})\otimes(\widehat{v}_{j+1}-\widehat{m}_{j+1})\mid
Y_{j}\Bigr{]}$
$\displaystyle=\operatorname{\mathbb{E}}\Bigl{[}(\Psi(v_{j})-\Psi(m_{j})+\xi_{j})\otimes(\Psi(v_{j})-\Psi(m_{j})+\xi_{j})\mid
Y_{j}\Bigr{]}$
$\displaystyle=\operatorname{\mathbb{E}}\Bigl{[}(\Psi(v_{j})-\Psi(m_{j}))\otimes(\Psi(v_{j})-\Psi(m_{j}))\mid
Y_{j}\Bigr{]}+\Sigma$ $\displaystyle\approx
D\Psi(m_{j})\operatorname{\mathbb{E}}\Bigl{[}(v_{j}-m_{j})\otimes(v_{j}-m_{j})\mid
Y_{j}\Bigr{]}D\Psi(m_{j})^{T}+\Sigma,$

and so, again assuming that fluctuations of $v_{j}$ around its mean $m_{j}$
(conditional on data) are small, we invoke the approximation

$\widehat{C}_{j+1}=D\Psi(m_{j})C_{j}D\Psi(m_{j})^{T}+\Sigma.$ (10.5)

To be self-consistent $\Sigma$ itself should be small.

The analysis step is the same as for the Kalman filter:

$\displaystyle\begin{split}S_{j+1}&=H\widehat{C}_{j+1}H^{T}+\Gamma,\\\
K_{j+1}&=\widehat{C}_{j+1}H^{T}S_{j+1}^{-1},\\\
m_{j+1}&=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1},\\\
C_{j+1}&=(I-K_{j+1}H)\widehat{C}_{j+1}.\\\ \end{split}$ (10.6)

Thus the overall ExKF comprises equations (10.4), (10.5) and (10.6). Unlike
the Kalman filter, for the extended Kalman filter the maps
$C_{j}\mapsto\widehat{C}_{j+1}\mapsto C_{j+1}$ depend on the observed data,
through the dependence of the predictive covariance on the filter mean. To be
self-consistent with the “small fluctuations around the mean” assumptions made
in the derivation of the ExKF, $\Sigma$ and $\Gamma$ should both be small.

The analysis step can also be defined by

$\displaystyle C_{j+1}^{-1}$
$\displaystyle=\widehat{C}_{j+1}^{-1}+H^{T}\Gamma^{-1}H,$ $\displaystyle
m_{j+1}$ $\displaystyle=\arg\min_{v}{\mathsf{J}}(v),$

where

${\mathsf{J}}(v)=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{m}_{j+1}|^{2}_{\widehat{C}_{j+1}}$
(10.7)

and $\widehat{m}_{j+1},\widehat{C}_{j+1}$ are calculated as above in the
prediction steps (10.4), (10.5). The constraint formulation of the
minimization problem, derived for the Kalman filter in section 8.1.2, may also
be used to derive the update formulae above.

####  10.3 Ensemble Kalman Filter

When the dynamical system is in high dimension, evaluation and storage of the
predicted covariance, and in particular the Jacobian required for the update
formula (10.5), becomes computationally inefficient and expensive for the
ExKF. The EnKF was developed to overcome this issue. The basic idea is to
maintain an ensemble of particles, and to use their empirical covariance
within a Kalman-type update. The method is summarized in Table 3. It may be
thought of as an ensemble 3DVAR technique in which a collection of particles
are generated similarly to 3DVAR, but interact through an ensemble estimate of
their covariance.

Kalman Filter | EnKF
---|---
$m_{j+1}=\arg\min_{v}{\mathsf{J}}(v)$ | $m_{j+1}=\arg\min_{v}{\mathsf{J}}(v)$
${\mathsf{J}}(v)=\frac{1}{2}|y_{j+1}-Hv|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{m}_{j+1}|_{\widehat{C}_{j+1}}^{2}$ | ${\mathsf{J}}_{n}(v)=\frac{1}{2}|y_{j+1}^{(n)}-Hm|_{\Gamma}^{2}+\frac{1}{2}|v-\widehat{v}_{j+1}^{(n)}|_{\widehat{C}_{j+1}}^{2}$
$\widehat{m}_{j+1}=Mm_{j}$ | $\widehat{v}_{j+1}^{(n)}=\Psi(v_{j}^{(n)})+\xi_{j}^{(n)}$
$\widehat{C}_{j+1}$ update exact |  $\widehat{C}_{j+1}$ update by ensemble estimate
$m_{j+1}=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1}$ | $m_{j+1}=(I-K_{j+1}H)\widehat{m}_{j+1}+K_{j+1}y_{j+1}$
Table 3: Comparison of Kalman Filter and EnKF update formulae

####  10.4 Derivation of EnKF

In the basic form which we present here, the EnKF is applied when $\Psi$ is
nonlinear, while the observation operator $H$ is linear. The $N$ particles
used at step $j$ are denoted $\\{v_{j}^{(n)}\\}^{N}_{n=1}.$ They are all given
equal weight so it is possible, in principle, to make an approximation to the
filtering distribution of the form

$\pi_{j}^{N}(v_{j})\approx\frac{1}{N}\sum_{n=1}^{N}\delta(v_{j}-v_{j}^{(n)}).$

However the EnKF is typically used with a relatively small number $N$ of
particles and may be far from approximating the desired distribution in
$\mathbb{R}^{d}.$ It is then better understood as a sequential optimization
method, similar in spirit to 3DVAR, as described above; this is our
perspective.

The state of all the particles at time $j+1$ are predicted to give
$\\{\widehat{v}_{j+1}^{(n)}\\}^{N}_{n=1}$ using the dynamical model. The
resulting empirical covariance is then used to define the objective function
(10.7c) which is minimized in order to perform the analysis step and obtain
$\\{{v}_{j+1}^{(n)}\\}^{N}_{n=1}.$ The updates are denoted schematically by

$\\{v_{j}^{(n)}\\}^{N}_{n=1}\xmapsto{ⓟ}\\{\widehat{v}_{j+1}^{(n)}\\}^{N}_{n=1}\xmapsto{ⓐ}\\{v_{j+1}^{(n)}\\}^{N}_{n=1}.$

We now detail these two steps.

ⓟ prediction:

$\displaystyle\widehat{v}_{j+1}^{(n)}$
$\displaystyle=\Psi(v_{j}^{(n)})+\xi^{(n)}_{j},\quad n=1,\ldots,N,$ (10.14a)
$\displaystyle\widehat{m}_{j+1}$
$\displaystyle=\frac{1}{N}\sum^{N}_{n=1}\widehat{v}_{j+1}^{(n)},$ (10.14b)
$\displaystyle\widehat{C}_{j+1}$
$\displaystyle=\frac{1}{N}\sum^{N}_{n=1}\bigl{(}\widehat{v}^{(n)}_{j+1}-\widehat{m}_{j+1}\bigr{)}\otimes\bigl{(}\widehat{v}^{(n)}_{j+1}-\widehat{m}_{j+1}\bigr{)}.$
(10.14c)

Here we have

$\xi_{j}^{(n)}\sim\mathcal{N}(0,\Sigma),\quad\text{i.i.d.}$

ⓐ analysis:

$\displaystyle S_{j+1}$ $\displaystyle=H\widehat{C}_{j+1}H^{T}+\Gamma,$
(10.15a) $\displaystyle K_{j+1}$
$\displaystyle=\widehat{C}_{j+1}H^{T}S_{j+1}^{-1},$ (10.15b) $\displaystyle
y_{j+1}^{(n)}$ $\displaystyle=y_{j+1}+s\eta_{j+1}^{(n)},\quad n=1,\ldots,N,$
(10.15c) $\displaystyle v_{j+1}^{(n)}$
$\displaystyle=(I-K_{j+1}H)\widehat{v}_{j+1}^{(n)}+K_{j+1}y_{j+1}^{(n)},\quad
n=1,\ldots,N.$ (10.15d)

Here we take

$\eta_{j}^{(n)}\sim\mathcal{N}(0,\Gamma),\quad\text{ i.i.d.}$

The constant $s$ takes value $0$ or $1$. When $s=1$ the $y_{j+1}^{(n)}$ are
referred to as perturbed observations. The analysis step may be written as

$v_{j+1}^{(n)}=\arg\min_{v}{\mathsf{J}}_{n}(v),$ (10.16)

where

${\mathsf{J}}_{n}(v):=\frac{1}{2}|y_{j+1}^{(n)}-Hv|^{2}_{\Gamma}+\frac{1}{2}|v-\widehat{v}_{j+1}^{(n)}|^{2}_{\widehat{C}_{j+1}}$
(10.17)

and the predicted mean and covariance are given by step ⓟ. Note that
$\widehat{C}_{j+1}$ is typically not invertible as it is a rank $N$ matrix and
$N$ is usually less than the dimension $d$ of the space on which
$\widehat{C}_{j+1}$ acts; this is since the typical use of ensemble methods is
for high dimensional state space estimation, with a small ensemble size. The
minimizing solution can be found by regularizing $\widehat{C}_{j+1}$ by
additions of $\epsilon I$ for $\epsilon>0$, deriving the update equations as
above for ⓐ, and then letting $\epsilon\to 0.$ Alternatively, the constraint
formulation of the minimization problem, derived for the Kalman filter in
section 8.1.2, may also be used to derive the update formulae above.

####  10.5 Subspace Property of EnKF

We now give another way to think of, and exploit in algorithms, the low rank
property of $\widehat{C}_{j+1}.$ Note that ${\mathsf{J}}_{n}(v)$ is undefined
unless

$v-\widehat{v}_{j+1}^{(n)}=\widehat{C}_{j+1}a$

for some $a\in\mathbb{R}^{d}.$ From the structure of $\widehat{C}_{j+1}$ given
in (10.14c) it follows that

$v=\widehat{v}_{j+1}^{(n)}+\frac{1}{N}\sum_{n=1}^{N}b_{n}\bigl{(}\widehat{v}^{(n)}_{j+1}-\widehat{m}_{j+1}\bigr{)}$
(10.18)

for some unknown parameters $\\{b_{n}\\}_{n=1}^{N},$ to be determined. (Note
that the vector $\\{b_{n}\\}_{n=1}^{N}$ depends on enesemble index $n$; we
have supressed this dependence for notational convenience.) This form for $v$
can be substituted into (10.17) to obtain a functional ${\mathsf{I}}_{n}(b)$
to be minimized over $b\in\mathbb{R}^{N}.$ We re-emphasize that $N$ will
typically be much smaller than $d$, the state-space dimension. Once $b$ is
determined it may be substituted back into (10.18) to obtain the solution to
the minimization problem.

To dig a little deeper into this calculation we define

$e^{(n)}=\widehat{v}^{(n)}_{j+1}-\widehat{m}_{j+1}$

and note that then

$\widehat{C}_{j+1}=\frac{1}{N}\sum_{n=1}^{N}e^{(n)}\otimes e^{(n)}.$

Since

$\widehat{C}_{j+1}a=\frac{1}{N}\sum_{n=1}^{N}b_{n}e^{(n)}$

we deduce that this is solved by taking

$b_{n}=\langle e^{(n)},a\rangle.$

Now note that

$\frac{1}{2}|v-\widehat{v}_{j+1}^{(n)}|^{2}_{\widehat{C}_{j+1}}=\frac{1}{2}\langle
a,\widehat{C}_{j+1}a\rangle=\frac{1}{2N}\sum_{n=1}^{N}b_{n}^{2}.$

We define

${\mathsf{I}}_{n}(b):=\frac{1}{2}\Bigl{|}y_{j+1}^{(n)}-H\widehat{v}_{j+1}^{(n)}-\frac{1}{N}\sum_{n=1}^{N}b_{n}H(\widehat{v}_{j+1}^{(n)}-\widehat{m}_{j+1})\Bigr{|}^{2}_{\Gamma}+\frac{1}{2N}\sum_{n=1}^{N}b_{n}^{2}.$
(10.19)

We have shown:

######  Theorem 10.1 (Implementation of EnKF in $N$ Dimensional Subspace).

Given the prediction ⓟ defined by (10.14a), the Kalman update formulae ⓐ may
be found by minimizing ${\mathsf{I}}_{n}(b)$ with respect to $b$ and
substituting into (10.18).

####  10.6 Discussion and Bibliography

The development and theory of the extended Kalman filter is documented in the
text [61]. A methodology for analyzing evolving probability distributions with
small variance, and establishing the validity of the Gaussian approximation,
is described in [102]. The use of the ExKF for weather forecasting was
proposed in [45]. However the dimension of the state space in most geophysical
applications renders the extended Kalman filter impractical. The ensemble
Kalman filter provided an innovation with far reaching consequences in
geophysical applications, because it allowed for the use of partial empirical
correlation information, without the computation of the full covariance. An
overview of ensemble Kalman methods may be found in the book [37], including a
historical perspective on the subject, originating from papers of Evensen and
Van Leeuwen in the mid 1990s [33, 35]; a similar idea was also developed by
Houtkamer within the Canadian meteorological service, around the same time;
[34, 56].

The presentation of the ensemble Kalman filter as a smart optimization tool is
also developed in [73], but the derivation of the update equations in a space
whose dimension is that of the ensemble is not described there. The analysis
of ensemble methods is difficult and theory is only just starting to emerge.
In the linear case the method converges in the large ensemble limit to the
Kalman filter [49], but in the nonlinear case the limit does not reproduce the
filtering distribution [32]. In any case the primary advantage of ensemble
methods is that they can provide good state estimation when the number of
particles is not large; this subject is discussed in [50, 69, 111, 112].

###  11 Particle Filter

This chapter is devoted to the particle filter, a method that approximates the
filtering distribution by a sum of Dirac measures. Particle filters provably
converge to the filtering distribution as the number of Dirac measures
approaches infinity. We focus on the bootstrap particle filter, also known as
sequential importance resampling; it is linked to the material on Monte Carlo
and importance sampling described in chapter 5. We note that the Kalman filter
completely characterizes the filtering distribution, but only in the linear
Gaussian setting; the Kalman-based methods introduced in the two previous
chapters outside the linear-Gaussian assumption are built by approximating the
predictive distribution via a Gaussian ansatz, and then using the Kalman
formulae for the analysis step. Similarly, particle filters approximate the
predictive distribution by a particle approximation measure and then solve
exactly the analysis step. In this light both Kalman-based methods (with
linear observations) and the bootstrap particle filter do exact application of
Bayes formula with wrong priors. However, particle filters have the potential
of recovering an accurate approximation to the filtering distribution in
nonlinear, non-Gaussian settings provided that the number of particles is
large enough; their caveat is that they tend to struggle in high dimensional
problems.

####  11.1 The Setting

Let us return to the setting in which we introduced filtering and smoothing,
with nonlinear stochastic dynamical system and nonlinear observation operator,
namely the model

$\displaystyle v_{j+1}$ $\displaystyle=\Psi(v_{j})+\xi_{j},$
$\displaystyle\xi_{j}\sim\mathcal{N}(0,\Sigma)\text{ i.i.d.},$ $\displaystyle
y_{j+1}$ $\displaystyle=h(v_{j+1})+\eta_{j},$
$\displaystyle\eta_{j}\sim\mathcal{N}(0,\Gamma)\text{ i.i.d.},$

with $v_{0}\sim\pi_{0}:=\mathcal{N}(m_{0},C_{0})$ independent of the i.i.d.
sequences $\\{\xi_{j}\\}$ and $\\{\eta_{j}\\}$. Here $\Psi(\cdot)$ drives the
dynamical system and $h(\cdot)$ is the observation operator. Recall that we
denote by $Y_{j}=\\{y_{1},\ldots,y_{j}\\}$ all the data up to time $j$ and by
$\pi_{j}$ the pdf of $v_{j}|Y_{j}$, that is,
$\pi_{j}=\operatorname{\mathbb{P}}(v_{j}|Y_{j})$. The filtering problem is to
determine $\pi_{j+1}$ from $\pi_{j}$. We may do so in two steps: first, we run
forward the Markov chain generated by the stochastic dynamical system
(prediction), and second, we incorporate the data by an application of Bayes
Theorem (analysis).

For the prediction step, we define the operator ${\mathcal{P}}$ acting on a
pdf $\pi$ as an application of a Markov kernel defined by

$({\mathcal{P}}\pi)(v)=\int_{\mathbb{R}^{d}}p(u,v)\pi(u)du,$ (11.1)

where $p(u,v)$ is the associated pdf of the stochastic dynamics, so that

$p(u,v)=\frac{1}{\sqrt{(2\pi)^{d}{\rm
det}\Sigma}}\text{exp}\mathopen{}\mathclose{{}\left(-\frac{1}{2}|v-\Psi(u)|_{\Sigma}^{2}}\right).$

Thus we obtain

$\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})=\widehat{\pi}_{j+1}={\mathcal{P}}\pi_{j}.$

We then define the analysis operator ${\mathcal{A}}_{j}$ acting on a pdf $\pi$
to correspond to an application of Bayes theorem, namely

$({\mathcal{A}}_{j}\pi)(v)=\frac{\text{exp}(-\frac{1}{2}|y_{j+1}-h(v)|_{\Gamma}^{2})\pi(v)}{\int_{\mathbb{R}^{d}}\text{exp}(-\frac{1}{2}|y_{j+1}-h(v)|_{\Gamma}^{2})\pi(v)dv}$
(11.2)

and so finally we obtain

$\pi_{j+1}={\mathcal{A}}_{j}\widehat{\pi}_{j+1}={\mathcal{A}}_{j}{\mathcal{P}}\pi_{j}.$

We now describe a way to numerically approximate, and update, the pdfs
$\pi_{j}.$

####  11.2 The Bootstrap Particle Filter

The bootstrap particle filter (BPF) can be thought of as performing sequential
importance resampling. Let $S^{N}$ be an operator acting on a pdf $\pi$ by
producing an $N-$samples Dirac approximation of $\pi$, that is

$(S^{N}\pi)(u)=\sum_{n=1}^{N}w_{n}\delta(u-u^{(n)}),$

where $u^{(1)},\ldots,u^{(N)}$ are i.i.d samples from $\pi$ that are weighted
uniformly i.e. $w_{n}=\frac{1}{N}$. Note that $S^{N}\pi=\pi^{N}_{MC}$, as
introduced in Chapter 5, equation (5.3). We will use the operator $S^{N}$ to
approximate the measure produced by the Markov kernel step ${\mathcal{P}}$
within the overall filtering map ${\mathcal{A}}_{j}{\mathcal{P}}$. Note that
$S^{N}$ is a random map taking pdfs into pdfs if we interpret weighted sums of
Diracs as a pdf.

Let $\pi_{0}^{N}=\pi_{0}=\mathcal{N}(m_{0},C_{0})$ and let $\pi^{N}_{j}$
denote a particle approximation of the pdf $\pi_{j}$ that we will determine in
what follows. We define

$\widehat{\pi}^{N}_{j+1}=S^{N}{\mathcal{P}}\pi^{N}_{j};$

this is an approximation of $\widehat{\pi}_{j+1}$ from the previous section.
We then apply the operator ${\mathcal{A}}_{j}$ to act on
$\widehat{\pi}^{N}_{j+1}$ by appropriately reconfiguring the weights $w_{j}$
according to the data.

To understand this reconfiguration of the weights we use the fact that, if

$({\mathcal{A}}\pi)(v)=\frac{g(v)\pi(v)}{\int_{\mathbb{R}^{d}}g(v)\pi(v)dv}$

and if

$\pi(v)=\frac{1}{N}\sum_{n=1}^{N}\delta(v-v^{(n)}),$

then

$({\mathcal{A}}\pi)(v)=\sum_{n=1}^{N}w^{(n)}\delta(v-v^{(n)}),$

where

$\bar{w}^{(n)}=g(v^{(n)})$

and the ${w}^{(n)}$ are found from the $\bar{w}^{(n)}$ by renormalizing them
to sum to one. We use this calculation concerning the application of Bayes
formula to sums of Diracs within the following desired approximation of the
filtering update formula:

$\pi_{j+1}\approx\pi_{j+1}^{N}={\mathcal{A}}_{j}\widehat{\pi}_{j+1}^{N}={\mathcal{A}}_{j}S^{N}{\mathcal{P}}\pi_{j}^{N}.$

The steps for the method are summarized in Algorithm 2.

Algorithm 2 Bootstrap Particle Filter

1: Input: Initial distribution $\pi_{0}^{N}=\pi_{0}$, observations $Y_{J}$,
number of particles $N.$

2: Particle Generation: For $j=0,1,\dots,J-1$, perform

  1. 1.

Draw $v_{j}^{(n)}\sim\pi_{j}^{N}$ for $n=1,\ldots,N$ i.i.d.

  2. 2.

Set $\widehat{v}^{(n)}_{j+1}=\Psi(v_{j}^{(n)})+\xi_{j}^{(n)}$ with
$\xi_{j}^{(n)}$ i.i.d. $\mathcal{N}(0,\Sigma).$

  3. 3.

Set
$\bar{w}_{j+1}^{(n)}=\text{exp}(-\frac{1}{2}|y_{j+1}-h(\widehat{v}^{(n)}_{j+1})|_{\Gamma}^{2}).$

  4. 4.

Set $w_{j+1}^{(n)}=\bar{w}_{j+1}^{(n)}/\sum_{n=1}^{N}\bar{w}_{j+1}^{(n)}.$

  5. 5.

Set
$\pi_{j+1}^{N}(u)=\sum_{n=1}^{N}w_{j+1}^{(n)}\delta(u-\widehat{v}^{(n)}_{j+1}).$

3: Output: pdf $\pi_{J}^{N}$ that approximates the distribution
$\operatorname{\mathbb{P}}(v_{J}|Y_{J}).$

####  11.3 Bootstrap Particle Filter Convergence

We will now show that under certain conditions, the BPF converges to the true
particle filter distribution in the limit $N\rightarrow\infty$. The proof is
similar to that of the Lax-Equivalence Theorem from the numerical
approximation of evolution equations, part of which is the statement that
consistency and stability together imply convergence. For the BPF consistency
refers to a Monte Carlo error estimate, similar to that derived in the chapter
on importance sampling, and stability manifests in bounds on the Lipschitz
constants for the operators ${\mathcal{P}}$ and ${\mathcal{A}}_{j}.$

Our first step is to define what we mean by convergence, that is, we need a
metric on probability measures. Notice that the operators ${\mathcal{P}}$ and
${\mathcal{A}}_{j}$ are deterministic, but the operator $S^{N}$ is random
since it requires sampling. As a consequence the approximate pdfs $\pi_{j}$
are also random. Thus, in fact, we need a metric on random probability
measures. To this end, for random probability measures $\pi$ and
$\pi^{\prime}$, we define

$d(\pi,\pi^{\prime})=\sup_{|f|_{\infty}\leq
1}\Bigl{(}\operatorname{\mathbb{E}}\Bigl{[}\bigl{(}\pi(f)-\pi^{\prime}(f)\bigr{)}^{2}\Bigr{]}\Bigr{)}^{1/2},$

where the expectation is taken over the random variable, in our case, the
randomness from sampling with $S^{N}$, and the supremum is taken over all
functions $f:\mathbb{R}^{d}\to[-1,1].$ Here we have used the notation defined
in equation (5.2). The following lemma is straightforward to prove, and
provides some useful intuition about the metric.

######  Lemma 11.1.

$d(\cdot,\cdot)$ as defined above does indeed define a metric on random
probability measures. Furthermore, when $\pi,\pi^{\prime}$ are deterministic,
then we have $d(\pi,\pi^{\prime})=2d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})$.

We now prove three lemmas which together will enable us to prove convergence
of the BPF. The first shows consistency; the second and third show stability
estimates for ${\mathcal{P}}$ and ${\mathcal{A}}_{j}$ respectively.

######  Lemma 11.2.

Let ${\mathscr{P}}$ be the set of probability densities on $\mathbb{R}^{d}$
then

$\sup_{\pi\in{\mathscr{P}}}d(\pi,S^{N}\pi)\leq\frac{1}{\sqrt{N}}.$

###### Proof.

Note that $S^{N}\pi$ agrees with $\pi_{MC}^{N}$ as defined in Chapter 5.
Therefore Theorem 5.1 shows that for any $f$ with $|f|_{\infty}\leq 1$ and for
any probability density $\pi$ it holds that

$\operatorname{\mathbb{E}}\mathopen{}\mathclose{{}\left[\bigl{(}S^{N}\pi(f)-\pi(f)\bigr{)}^{2}}\right]\leq\frac{1}{N}.$

Taking the square root and the supremum over $f$ and $\pi$ gives the desired
result. ∎

Now we prove a stability bound for the operator ${\mathcal{P}}$ defined in
equation (11.1).

######  Lemma 11.3.

$d({\mathcal{P}}\pi,{\mathcal{P}}\pi^{\prime})\leq d(\pi,\pi^{\prime}).$

###### Proof.

For $|f|_{\infty}\leq 1$ define a function $q$ on $\mathbb{R}^{d}$ by

$q(v^{\prime})=\int_{\mathbb{R}^{d}}p(v^{\prime},v)f(v)dv,$

where, recall, $p$ denotes the transition pdf associated to the stochastic
dynamics. Note that

$|q(v^{\prime})|\leq\int_{\mathbb{R}^{d}}p(v^{\prime},v)dv=1,$

and so

$\displaystyle\pi(q)$
$\displaystyle=\int_{\mathbb{R}^{d}}q(v^{\prime})\pi(v^{\prime})dv^{\prime}=\int_{\mathbb{R}^{d}}\Bigg{[}\int_{\mathbb{R}^{d}}p(v^{\prime},v)f(v)dv\Bigg{]}\pi(v^{\prime})dv^{\prime}$
$\displaystyle=\int_{\mathbb{R}^{d}}\Bigg{[}\int_{\mathbb{R}^{d}}p(v^{\prime},v)\pi(v^{\prime})dv^{\prime}\Bigg{]}f(v)dv$
$\displaystyle=\int_{\mathbb{R}^{d}}({\mathcal{P}}\pi)(v)f(v)dv$

by exchanging the order of integration. Consequently, we have

$\pi(q)=({\mathcal{P}}\pi)(f).$

Finally,

$\displaystyle d({\mathcal{P}}\pi,{\mathcal{P}}\pi^{\prime})$
$\displaystyle=\sup_{|f|_{\infty}\leq
1}\Bigl{(}\operatorname{\mathbb{E}}\Bigl{[}\bigl{(}({\mathcal{P}}\pi)(f)-({\mathcal{P}}\pi^{\prime})(f)\bigr{)}^{2}\Bigr{]}\Bigr{)}^{1/2}$
$\displaystyle\leq\sup_{|q|_{\infty}\leq
1}\Bigr{(}\operatorname{\mathbb{E}}\Bigl{[}\bigl{(}\pi(q)-\pi^{\prime}(q)\bigr{)}^{2}\Bigr{]}\Bigr{)}^{1/2}$
$\displaystyle=d(\pi,\pi^{\prime}).$

∎

To prove the next lemma and the main convergence theorem of the BPF below, we
will make the following assumption which encodes the idea of a bound on the
observation operator:

######  Assumption 11.4.

There exists $\kappa\in(0,1)$ such that, for all $v\in\mathbb{R}^{d},$ and
$j\in\\{1,\ldots,J\\},$

$\kappa\leq g_{j}(v)\leq\kappa^{-1}.$

It may initially appear strange to use the same constant $\kappa$ in the upper
and lower bounds, but recall that $g$ is undefined up to a multiplicative
constant. Consequently, given any upper and lower bounds, $g$ can be scaled to
achieve the bound as stated. Relatedly it is $\kappa^{-2}$ which appears in
the stability constant in the next lemma; if $g$ is not scaled to produce the
same constant $\kappa$ in the upper and lower bounds in Assumption 11.4, then
it is the ratio of the upper and lower bounds which would appear in the
stability bound.

######  Lemma 11.5.

Let Assumption 11.4 hold. Then

$d({\mathcal{A}}_{j}\pi,{\mathcal{A}}_{j}\pi^{\prime})\leq\frac{2}{\kappa^{2}}d(\pi,\pi^{\prime}).$

###### Proof.

$\displaystyle({\mathcal{A}}\pi)(f)-({\mathcal{A}}\pi^{\prime})(f)$
$\displaystyle=\frac{\pi(fg)}{\pi(g)}-\frac{\pi^{\prime}(fg)}{\pi^{\prime}(g)}$
$\displaystyle=\frac{\pi(fg)}{\pi(g)}-\frac{\pi^{\prime}(fg)}{\pi(g)}+\frac{\pi^{\prime}(fg)}{\pi(g)}-\frac{\pi^{\prime}(fg)}{\pi^{\prime}(g)}$
$\displaystyle=\frac{1}{\kappa}\bigg{(}\frac{\pi(\kappa
fg)-\pi^{\prime}(\kappa
fg)}{\pi(g)}+\frac{\pi^{\prime}(fg)}{\pi^{\prime}(g)}\frac{\pi^{\prime}(\kappa
g)-\pi(\kappa g)}{\pi(g)}\bigg{)}.$

Applying Bayes Theorem we obtain

$\Bigl{|}\frac{\pi^{\prime}(fg)}{\pi^{\prime}(g)}\Bigr{|}=|({\mathcal{A}}\pi^{\prime})(f)|\leq
1.$

Therefore,

$\bigl{|}({\mathcal{A}}\pi)(f)-({\mathcal{A}}\pi^{\prime})(f)\bigr{|}\leq\frac{1}{\kappa^{2}}\Bigl{(}\bigl{|}\pi(\kappa
fg)-\pi^{\prime}(\kappa fg)\bigr{|}+\bigl{|}\pi^{\prime}(\kappa g)-\pi(\kappa
g)\bigr{|}\Bigr{)}.$

It follows that

$\operatorname{\mathbb{E}}\Bigl{[}\bigl{(}({\mathcal{A}}\pi)(f)-({\mathcal{A}}\pi^{\prime})(f)\bigr{)}^{2}\Bigr{]}\leq\frac{2}{\kappa^{2}}\Bigl{(}\operatorname{\mathbb{E}}\Bigl{[}(\pi(\kappa
fg)-\pi^{\prime}(\kappa
fg))^{2}\Bigr{]}+\operatorname{\mathbb{E}}\Bigl{[}\bigl{(}\pi^{\prime}(\kappa
g)-\pi(\kappa g)\bigr{)}^{2}\Bigr{]}\Bigr{)}.$

Since $|\kappa g|\leq 1$ we find that

$\sup_{|f|_{\infty}\leq
1}\operatorname{\mathbb{E}}\Bigl{[}\bigl{(}({\mathcal{A}}\pi)(f)-({\mathcal{A}}\pi^{\prime})(f)\bigr{)}^{2}\Bigr{]}\leq\frac{4}{\kappa^{4}}\sup_{|f|_{\infty}\leq
1}\operatorname{\mathbb{E}}\Bigl{[}\bigl{(}\pi(f)-\pi^{\prime}(f)\bigr{)}^{2}\Bigr{]},$

and hence

$d({\mathcal{A}}_{j}\pi,{\mathcal{A}}_{j}\pi^{\prime})\leq\frac{2}{\kappa^{2}}d(\pi,\pi^{\prime}).$

∎

######  Theorem 11.6 (Convergence of the BPF).

Let Assumption 11.4 hold. Then there exists a $c=c(J,\kappa)$ such that, for
all $j=1,\ldots,J,$

$d(\pi_{j},\pi_{j}^{N})\leq\frac{c}{\sqrt{N}}.$

###### Proof.

Let $e_{j}=d(\pi_{j},\pi_{j}^{N})$, then

$\displaystyle e_{j+1}$
$\displaystyle=d(\pi_{j+1},\pi_{j+1}^{N})=d({\mathcal{A}}_{j}{\mathcal{P}}\pi_{j},{\mathcal{A}}_{j}S^{N}{\mathcal{P}}\pi_{j}^{N})$
$\displaystyle\leq
d({\mathcal{A}}_{j}{\mathcal{P}}\pi_{j},{\mathcal{A}}_{j}{\mathcal{P}}\pi_{j}^{N})+d({\mathcal{A}}_{j}{\mathcal{P}}\pi_{j}^{N},{\mathcal{A}}_{j}S^{N}{\mathcal{P}}\pi_{j}^{N})$

by the triangle inequality. Applying the stability bound for
${\mathcal{A}}_{j}$, we have

$e_{j+1}\leq\frac{2}{\kappa^{2}}\Bigl{[}d({\mathcal{P}}\pi_{j},{\mathcal{P}}\pi_{j}^{N})+d(\widehat{\pi}_{j}^{N},S^{N}\widehat{\pi}_{j}^{N})\Bigr{]},$

where $\widehat{\pi}_{j}^{N}={\mathcal{P}}\pi_{j}^{N}$. By the stability bound
for ${\mathcal{P}}$,

$d({\mathcal{P}}\pi_{j},{\mathcal{P}}\pi_{j}^{N})\leq d(\pi_{j},\pi_{j}^{N})$

and by the consistency bound for $S^{N}$

$d(\widehat{\pi}_{j}^{N},S^{N}\widehat{\pi}_{j}^{N})\leq\frac{1}{\sqrt{N}}.$

Therefore,

$\displaystyle e_{j+1}$
$\displaystyle\leq\frac{2}{\kappa^{2}}\Bigl{(}d(\pi_{j},\pi_{j}^{N})+\frac{1}{\sqrt{N}}\Bigr{)}$
$\displaystyle\leq\frac{2}{\kappa^{2}}\Bigl{(}e_{j}+\frac{1}{\sqrt{N}}\Bigr{)}.$

We let $\lambda=2/\kappa^{2}$ and note that $\lambda\geq 2$ since
$\kappa\in(0,1]$. Then the discrete Gronwall inequality gives

$e_{j}\leq\lambda^{j}e_{0}+\frac{\lambda}{\sqrt{N}}\frac{1-\lambda^{j}}{1-\lambda}.$

Recall that $\pi_{0}^{N}=\pi_{0}$ hence $e_{0}=0$. Thus letting

$c=\frac{\lambda(1-\lambda^{J})}{1-\lambda}$

completes the proof since $\lambda(1-\lambda^{j})/(1-\lambda)$ is increasing
in $j$. ∎

####  11.4 The Bootstrap Particle Filter as a Random Dynamical System

A nice interpretation of the BPF is to view it as a random dynamical system
for a set of interacting particles $\\{v_{j}^{(n)}\\}_{n=1}^{N}.$ To this end,
a measure

$\bar{\pi}_{j}^{N}(u)=\frac{1}{N}\sum_{n=1}^{N}\delta(u-v_{j}^{(n)})\approx\pi_{j}^{N}(u)\approx\pi_{j}(u)$

with equally weighted particles may be naturally defined after the resampling
step from $\pi_{j}^{N}.$ It can then be seen that the BPF updates the particle
positions

$\\{v_{j}^{(n)}\\}_{n=1}^{N}\mapsto\\{v_{j+1}^{(n)}\\}_{n=1}^{N}$

via the random map

$\displaystyle\widehat{v}_{j+1}^{(n)}$
$\displaystyle=\Psi(v_{j}^{(n)})+\xi_{j}^{(n)},$
$\displaystyle\xi_{j}^{(n)}\sim\mathcal{N}(0,\Sigma)\text{ i.i.d.,}$
$\displaystyle v_{j+1}^{(n)}$
$\displaystyle=\sum_{m=1}^{N}\mathbb{1}_{I_{j+1}^{(m)}}\mathopen{}\mathclose{{}\left(r_{j+1}^{(n)}}\right)\widehat{v}_{j+1}^{(m)},$
$\displaystyle r_{j+1}^{(n)}\sim\text{Uniform}(0,1)\text{ i.i.d.}$

Here the supports $I_{j}^{(m)}$ of the indicator functions have widths given
by the weights appearing in ${\pi}_{j}^{N}(u).$ Specifically we have

$\displaystyle I_{j+1}^{(m)}$
$\displaystyle=\mathopen{}\mathclose{{}\left[\alpha_{j+1}^{(m-1)},\alpha_{j+1}^{(m)}}\right),$
$\displaystyle\alpha_{j+1}^{(m+1)}=\alpha_{j+1}^{(m)}+w_{j+1}^{(m)},\quad\alpha_{j+1}^{(0)}=0.$

Note that, by construction, $\alpha_{j}^{(N)}=1.$

Thus the underlying dynamical system on particles comprises $N$ particles
governed by two steps: (i) the underlying stochastic dynamics model, in which
the particles do not interact; (ii) a resampling of the resulting collection
of particles, to reflect the different weights associated with them, in which
the particles do then interact. The interaction is driven by the weights which
see all the particle positions, and measure their goodness of fit to the data.

####  11.5 Discussion and Bibliography

Particle filters are overviewed from an algorithmic viewpoint in [29, 28], and
from a more mathematical perspective in [25]. The convergence of particle
filters was addressed in [22]; the clean proof presented here originates in
[97] and may also be found in [73]. For problems in which the dynamics evolve
in relatively low dimensional spaces they have been an enormously successful.
Generalizing them so that they work for the high dimensional problems that
arise, for example, in geophysical applications, provides a major challenge
[75].

###  12 Optimal Particle Filter

This chapter is devoted to the optimal particle filter. Like the bootstrap
filter from the previous chapter, the optimal particle filter approximates the
filtering distribution by a sum of Dirac measures. The setting will initially
be the same (nonlinear stochastic dynamics and nonlinear observations), but we
will see that the optimal particle filter cannot be implemented, in general,
in the fully nonlinear case. For this reason, we will specify to the case of
linear observation operator, where the optimal particle filter can be
implemented in a straightforward fashion and may be characterized as a set of
interacting 3DVAR filters. We conclude this chapter by summarizing and
comparing all of the filtering methods introduced in this and preceding
chapters, highlighting the setting in which they may be practically applied.

####  12.1 Introduction

The bootstrap particle filter from the previous chapter is based on
approximating the two components of a specific factorization of the filtering
step. The factorization is

$\displaystyle\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})$
$\displaystyle={\mathcal{P}}\operatorname{\mathbb{P}}(v_{j}|Y_{j}),$
$\displaystyle\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})$
$\displaystyle={\mathcal{A}}_{j}\operatorname{\mathbb{P}}(v_{j+1}|Y_{j}).$

This gives the factorization

$\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})={\mathcal{A}}_{j}{\mathcal{P}}\operatorname{\mathbb{P}}(v_{j}|Y_{j}),$

which is the basis for the bootstrap particle filter. It is natural to ask if
there are other factorizations of the filtering update and whether they might
lead to improved particle filters. In this lecture we derive the Optimal
Particle Filter (OPF) which does just that. We demonstrate a connection with
3DVAR, and we discuss the sense in which the OPF has desirable properties in
comparison with the BPF. As usual we denote by $\pi_{j}$ the filtering
distribution $\operatorname{\mathbb{P}}(v_{j}|Y_{j}).$

####  12.2 The Bootstrap and Optimal Particle Filters Compared

The fundamental filtering problem that we are interested in is determination
of $\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})$ from
$\operatorname{\mathbb{P}}(v_{j}|Y_{j}).$ In the BPF, we are approximating the
following manipulation:

$\displaystyle\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})$
$\displaystyle=\operatorname{\mathbb{P}}(v_{j+1}|y_{j+1},Y_{j})$
$\displaystyle=\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1},Y_{j})\operatorname{\mathbb{P}}(v_{j+1}|Y_{j})}{\operatorname{\mathbb{P}}(y_{j+1}|Y_{j})}$
$\displaystyle=\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1},Y_{j})}{\operatorname{\mathbb{P}}(y_{j+1}|Y_{j})}\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j},Y_{j})\operatorname{\mathbb{P}}(v_{j}|Y_{j})dv_{j}$
$\displaystyle=\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})}{\operatorname{\mathbb{P}}(y_{j+1}|Y_{j})}\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j})\operatorname{\mathbb{P}}(v_{j}|Y_{j})dv_{j}$
$\displaystyle={\mathcal{A}}_{j}{\mathcal{P}}\operatorname{\mathbb{P}}(v_{j}|Y_{j}).$

The Markov kernel ${\mathcal{P}}^{BPF}$ acts on arbitrary density $\pi$ by

${\mathcal{P}}^{BPF}\pi(v_{j+1})=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j})\pi(v_{j})dv_{j},$

and ${\mathcal{A}}_{j}$ acts on an arbitrary density $\pi$ by application of
Bayes theorem, taking into account the likelihood of the data

${\mathcal{A}}_{j}\pi(v_{j+1})=\frac{1}{Z}\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\pi(v_{j+1}),$

with $Z$ normalization to a probability density. The above manipulations are
summarized by the relationship

$\pi_{j+1}={\mathcal{A}}_{j}{\mathcal{P}}\pi_{j}.$ (12.1)

Note that in this factorization we apply a Markov kernel and then Bayes
theorem. In contrast, to derive the OPF we perform the following manipulation:

$\displaystyle\operatorname{\mathbb{P}}(v_{j+1}|Y_{j+1})$
$\displaystyle=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1},v_{j}|Y_{j+1})\,dv_{j}$
$\displaystyle=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j},Y_{j+1})\operatorname{\mathbb{P}}(v_{j}|Y_{j+1})dv_{j}\,$
$\displaystyle=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1},Y_{j})\operatorname{\mathbb{P}}(v_{j}|y_{j+1},Y_{j})\,dv_{j}$
$\displaystyle=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1})\operatorname{\mathbb{P}}(v_{j}|y_{j+1},Y_{j})\,dv_{j}$
$\displaystyle=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1})\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j},Y_{j})}{\operatorname{\mathbb{P}}(y_{j+1}|Y_{j})}\operatorname{\mathbb{P}}(v_{j}|Y_{j})\,dv_{j}$
$\displaystyle=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1})\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j})}{\operatorname{\mathbb{P}}(y_{j+1}|Y_{j})}\operatorname{\mathbb{P}}(v_{j}|Y_{j})\,dv_{j}$
$\displaystyle={\mathcal{P}}^{OPF}_{j}{\mathcal{A}}_{j}^{OPF}\operatorname{\mathbb{P}}(v_{j}|Y_{j}),$

with Markov kernel for particle update

${\mathcal{P}}_{j}^{OPF}\pi(v_{j+1})=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1})\pi(v_{j})dv_{j}$

and application of Bayes theorem to include the likelihood

${\mathcal{A}}_{j}^{OPF}\pi(v_{j})=\frac{1}{Z}\operatorname{\mathbb{P}}(y_{j+1}|v_{j})\pi(v_{j}).$

Thus we have

$\pi_{j+1}={\mathcal{P}}_{j}^{OPF}{\mathcal{A}}_{j}^{OPF}\pi_{j}.$ (12.2)

Note that in the factorization given by OPF we apply Bayes theorem and then a
Markov kernel, the opposite order to the BPF. Moreover the propagation
mechanism is different–it sees the data through the Markov kernel
${\mathcal{P}}_{j}^{OPF}$– and hence the weighting of the particles is also
different: the OPF weights are proportional to the likelihood
$\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})$ and the BPF weights are
proportional to $\operatorname{\mathbb{P}}(y_{j+1}|v_{j})$ which may be, in
general, not available in closed form. We will see that the particle updates
use a 3DVAR procedure. In the BPF, the evolution of the particles and the
observation of the data are kept separate from each other – the Markov kernel
${\mathcal{P}}$ depends only on the dynamics and not the observed data and is
thus independent of $j$.

####  12.3 Implementation of Optimal Particle Filter: Linear Observation
Setting

In general it is not possible to implement the OPF in the fully nonlinear
setting because of two computational bottlenecks:

  * •

There may not be a closed formula for evaluating the likelihood
$\operatorname{\mathbb{P}}(y_{j+1}|v_{j}),$ making unfeasible the computation
of the particle weights.

  * •

It may not be possible to sample from the Markov kernel
$\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1}),$ making unfeasible the
propagation of particles.

However, when the observation function $h(\cdot)$ is linear, i.e.
$h(\cdot)=H\cdot$ for some $H\in R^{k\times d}$ both bottlenecks are overcome.
We thus consider the following setting, which arises in many applications:

$\displaystyle v_{j+1}$ $\displaystyle=\Psi(v_{j})+\xi_{j},$
$\displaystyle\xi_{j}\sim\mathcal{N}(0,\Sigma)\text{ i.i.d.},$ $\displaystyle
y_{j+1}$ $\displaystyle=Hv_{j+1}+\eta_{j+1},$
$\displaystyle\eta_{j+1}\sim\mathcal{N}(0,\Gamma)\text{ i.i.d.},$

with $v_{0}\sim\mathcal{N}(m_{0},C_{0})$ and
$v_{0},\\{\xi_{j}\\},\\{\eta_{j}\\}$ independent. We will now show that in
this linear observation both bottlenecks are overcome. First, note that
combining the dynamics and data models we may write

$y_{j+1}=H\Psi(v_{j})+H\xi_{j}+\eta_{j+1},$

which shows that the conditional distribution for $y_{j+1}$ given $v_{j}$ is

$\operatorname{\mathbb{P}}(y_{j+1}|v_{j})=\mathcal{N}(H\Psi(v_{j}),S),$

where $S=H\Sigma H^{T}+\Gamma.$ We will use this formula to compute the
weights, thus overcoming the first computational bottleneck.

We now show that under the linear observation assumption
${\mathcal{P}}^{OPF}_{j}$ is a Gaussian kernel, and hence the second
computational bottleneck is overcome too. We have

$\displaystyle\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1})$
$\displaystyle\propto\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1},v_{j})\operatorname{\mathbb{P}}(v_{j+1}|v_{j})$
$\displaystyle=\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|v_{j})$
$\displaystyle\propto\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|y_{j+1}-Hv_{j+1}|^{2}_{\Gamma}-\frac{1}{2}|v_{j+1}-\Psi(v_{j})|^{2}_{\Sigma}}\right)$
$\displaystyle=\exp\bigl{(}-{\mathsf{J}}_{\textsc{opt}}(v_{j+1})\bigr{)}.$

This is a Gaussian distribution for $v_{j+1}$ as
${\mathsf{J}}_{\textsc{opt}}(v_{j+1})$ is quadratic with respect to
$v_{j+1}$.333${\mathsf{J}}_{\textsc{opt}}$ is identical to ${\mathsf{J}}$ on
the right-hand side of Table 1, with $\widehat{C}$ replaced by $\Sigma.$
Consequently, we can compute the mean $m_{j+1}$ and covariance $C$ (which,
note, is independent of $j$) of this Gaussian by matching the mean and
quadratic terms in the relevant quadratic forms:

$\displaystyle C^{-1}$ $\displaystyle=H^{T}\Gamma^{-1}H+\Sigma^{-1},$
$\displaystyle C^{-1}m_{j+1}$
$\displaystyle=\Sigma^{-1}\Psi(v_{j})+H^{T}\Gamma^{-1}y_{j+1}.$

Then
$\operatorname{\mathbb{P}}(v_{j+1}|y_{j+1},v_{j})=\mathcal{N}(m_{j+1},C).$
This is hence a special case of 3DVAR in which the analysis covariance is
fixed at $C$; note that when we derived 3DVAR we fixed the predictive
covariance $\widehat{C}$ which, here, is fixed at $\Sigma.$ As with the Kalman
filter, and with 3DVAR, it is possible to implement the prediction step
through the following mean and covariance formulae which avoid inversion in
state space, and require inversion only in data space:

$\displaystyle m_{j+1}$ $\displaystyle=(I-KH)\Psi(v_{j})+Ky_{j+1},$
$\displaystyle C$ $\displaystyle=(I-KH)\Sigma,$ $\displaystyle K$
$\displaystyle=\Sigma H^{T}S^{-1},$ $\displaystyle S$ $\displaystyle=H\Sigma
H^{T}+\Gamma.$

Furthermore, as for 3DVAR, the inversion of $S$ need only be performed once in
a pre-processing step before the algorithm is run. Since the expression for
$\operatorname{\mathbb{P}}(v_{j+1}|v_{j},y_{j+1})$ is Gaussian we now have the
ability to sample directly from ${\mathcal{P}}_{j}^{OPF}$. The OPF is thus
given by the following update algorithm for approximations
$\pi_{j}^{N}\approx\pi_{j}$ in which we generalize the notational conventions
used in the previous chapter to formulate particle filters as random dynamical
systems:

Algorithm 3 Algorithm for the Optimal Particle Filter with linear observation
map

1: Input: Initial distribution $\operatorname{\mathbb{P}}(v_{0})=\pi_{0}$,
observations $Y_{J}$, number of particles $N.$

2: Initial Sampling: Draw $N$ particles $v_{0}^{(n)}\sim\pi_{0}$ so that
$\pi_{0}^{N}=S^{N}\pi_{0}.$

3: Subsequent Sampling For $j=0,1,\dots,J-1$, perform:

  1. 1.

Set
$\widehat{v}_{j+1}^{(n)}=(I-KH)\Psi(v_{j}^{(n)})+Ky_{j+1}+\zeta_{j+1}^{(n)}$
with $\zeta_{j+1}^{(n)}$ i.i.d. $\mathcal{N}(0,C).$

  2. 2.

Set
$\bar{w}_{j+1}^{(n)}=\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|y_{j+1}-H\Psi(v_{j}^{(n)})|_{S}^{2}}\right).$

  3. 3.

Set $w_{j+1}^{(n)}=\bar{w}_{j+1}^{(n)}/\sum_{n=1}^{N}\bar{w}_{j+1}^{(n)}.$

  4. 4.

Set
$v_{j+1}^{(n)}=\sum_{m=1}^{N}{\mathbb{1}}_{I_{j+1}^{(m)}}(r_{j+1}^{(n)})\widehat{v}_{j+1}^{(m)}.$

  5. 5.

Set
$\pi_{j+1}^{N}(v_{j+1})=\frac{1}{N}\sum_{n=1}^{N}\delta\mathopen{}\mathclose{{}\left(v_{j+1}-{v}_{j+1}^{(n)}}\right).$

4: Output: $N$ particles $v_{J}^{(1)},v_{J}^{(2)},\dots,v_{J}^{(N)}.$

It would be desirable to interpret this algorithm as an approximation of the
filter update (12.2) in the form

$\pi_{j+1}^{N}=S^{N}{\mathcal{P}}_{j}^{OPF}{\mathcal{A}}_{j}^{OPF}\pi_{j}^{N},\quad\pi_{j}^{0}=S^{N}\pi_{0}.$

However the order in which the resampling and the particle propagation occurs
implies that this is not possible. The following slight modification of the
OPF, however, may indeed be thought of as an approximation of this form; we
simply reorder the resampling and the propagation. We refer to the resulting
algorithm as the Gaussianized Optimal Particle filter (GOPF). We may write the
resulting algorithm as follows:

Algorithm 4 Algorithm for the Gaussianized Optimal Particle Filter

1: Input: Initial distribution $\operatorname{\mathbb{P}}(v_{0})=\pi_{0}$,
observations $Y_{J}$, number of particles $N.$

2: Initial Sampling: Draw $N$ particles $v_{0}^{(n)}\sim\pi_{0}$ so that
$\pi_{0}^{N}=S^{N}\pi_{0}.$

3: Subsequent Sampling For $j=0,1,\dots,J-1$, perform:

  1. 1.

Set
$\bar{w}_{j+1}^{(n)}=\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|y_{j+1}-H\Psi(v_{j}^{(n)})|_{S}^{2}}\right).$

  2. 2.

Set $w_{j+1}^{(n)}=\bar{w}_{j+1}^{(n)}/\sum_{n=1}^{N}\bar{w}_{j+1}^{(n)}.$

  3. 3.

Set
$\widehat{v}_{j}^{(n)}=\sum_{m=1}^{N}{\mathbb{1}}_{I_{j+1}^{(m)}}(r_{j+1}^{(n)}){v}_{j}^{(m)}.$

  4. 4.

Set
${v}_{j+1}^{(n)}=(I-KH)\Psi(\widehat{v}_{j}^{(n)})+Ky_{j+1}+\zeta_{j+1}^{(n)}$
with $\zeta_{j+1}^{(n)}$ i.i.d $\mathcal{N}(0,C).$

  5. 5.

Set
$\pi_{j+1}^{N}(v_{j+1})=\frac{1}{N}\sum_{n=1}^{N}\delta\mathopen{}\mathclose{{}\left(v_{j+1}-{v}_{j+1}^{(n)}}\right).$

4: Output: $N$ particles $v_{J}^{(1)},v_{J}^{(2)},\dots,v_{J}^{(N)}.$

####  12.4 “Optimality” of the Optimal Particle Filter

Particle filter methods rely on approximating the distribution for the model
by a swarm of point Dirac functions; it is clear that the distribution will
not be well approximated by only a small number of particles in most cases.
Consequently, a performance requirement for particle filter methods is that
they do not lead to degeneracy of the particles. Resampling leads to
degeneracy if a small number of particles have all the weights. Conversely,
non-degeneracy may be promoted by ensuring that the weights $w_{j}^{(n)}$ are
similar in magnitude, so that a small number of particles are not overly
favoured during the resampling step. This condition can be formulated as a
requirement that the variance of the weights be minimized; doing this results
in the OPF.

To understand this perspective we consider an arbitrary particle update kernel
of the form $\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})$ and we study the resulting
particle filter without resampling. It is then the case that the particle
weights are updated according to the formula

$\bar{w}^{(n)}_{j+1}=\bar{w}^{(n)}_{j}\frac{\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|v^{(n)}_{j})}{\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})}.$

######  Theorem 12.1 (Meaning of Optimality).

The choice of $\operatorname{\mathbb{P}}(v_{j+1}|v_{j}^{(n)},y_{j+1})$ as the
particle update kernel $\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})$ results in the
minimal variance of the weight $w^{(n)}_{j+1}$ with respect to all possible
choices of the particle update kernel $\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1}).$

###### Proof.

We calculate the variance of the unnormalized weights (treated as random
variables) $\bar{w}^{(n)}_{j+1}$ with respect to the transition density
$\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})$ and obtain

$\displaystyle{\rm Var}_{\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})}$
$\displaystyle[\bar{w}^{(n)}_{j+1}]=\int_{\mathbb{R}^{d}}\mathopen{}\mathclose{{}\left(\bar{w}^{(n)}_{j+1}}\right)^{2}\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})dv_{j+1}$
$\displaystyle-\mathopen{}\mathclose{{}\left[\int_{\mathbb{R}^{d}}\bar{w}^{(n)}_{j+1}\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})dv_{j+1}}\right]^{2}$
$\displaystyle=\mathopen{}\mathclose{{}\left(\bar{w}^{(n)}_{j}}\right)^{2}\int_{\mathbb{R}^{d}}\frac{\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|v^{(n)}_{j})}\right)^{2}}{\pi(v_{j+1}|v_{j}^{(n)},Y_{j+1})}dv_{j+1}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad-\mathopen{}\mathclose{{}\left(\bar{w}^{(n)}_{j}}\right)^{2}\mathopen{}\mathclose{{}\left[\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(y_{j+1}|v_{j+1})\operatorname{\mathbb{P}}(v_{j+1}|v^{(n)}_{j})dv_{j+1}}\right]^{2}$