Particle filters for high-dimensional geoscience applications: a review
Peter Jan van Leeuwen Hans R. Künsch Lars Nerger Roland Potthast Sebastian Reich

figures/

Peter Jan van Leeuwen, et al.Particle Filters for Applications in Geosciences





aDepartment of Meteorology and National Centre for Earth Observation, University of Reading, Reading, UK

                bDepartment of Atmospheric Sciences, Colorado State University, USA

                 cSeminar für Statistik, ETH Zurich, Swiszerland

               dAlfred Wegener Institute Helmholtz Center for Polar and Marine Research, Bremerhaven, Germany

               e Deutscher Wetterdienst, Offenbach, Germany

                fInstitut für Mathematik, Universitaet Potsdam,  Potsdam, Germany
               


Peter Jan van Leeuwen, Department of Meteorology, University of Reading, Reading RG6 6BB, UK. E-mail: p.j.vanleeuwen@reading.ac.uk







Particle filters contain the promise of fully nonlinear data assimilation.
They have been applied in numerous science areas, including the geosciences, 
but their application to high-dimensional geoscience systems
has been limited due to their inefficiency in high-dimensional systems
in standard settings. However, huge progress has been made, and this limitation
is disappearing fast due to recent developments in proposal densities, 
the use of ideas from (optimal) transportation,
the use of localisation and intelligent adaptive resampling strategies. 
Furthermore, powerful hybrids between particle filters
and ensemble Kalman filters and variational methods have been developed. 
We present a state of the art discussion of present efforts of developing particle filters
for high-dimensional nonlinear geoscience state-estimation problems with an 
emphasis on atmospheric and oceanic applications,
including many new ideas, derivations, and unifications, highlighting hidden connections,
including pseudo code,
and generating a valuable tool and guide for the community.
Initial experiments show that particle filters can be competitive with present-day methods for 
numerical weather prediction suggesting that they will become mainstream soon.





Particle filters for high-dimensional geoscience applications: A Review [2]
    Peter Jan van Leeuwena,b, Hans R. Künschc,Lars Nergerd, Roland Potthaste, Sebastian Reichf
    
==============================================================================================



[2]Please ensure that you use the most up to date class file, available from the QJRMS Home Page at 

http://onlinelibrary.wiley.com/journal/10.1002/(ISSN)1477-870X






§ INTRODUCTION




Data assimilation for geoscience applications, such as weather or ocean prediction, is a slowly maturing field.
Even the linear data-assimilation problem cannot be solved adequately because of the size of the problem.
Typically, global-scale numerical weather prediction needs estimation of over 10^9 state variables, assimilating 
over 10^7 observations every 6-12 hours. 
Existing methods like 4DVar do not provide accurate uncertainty estimates and need efficient pre-conditioners, 
while Ensemble Kalman Filters heavily
rely on somewhat ad-hoc fixes like localisation and inflation to find accurate estimates.
Hybrids of variational and ensemble Kalman filter methods are a step forward, although localisation
and inflation are still needed in realistic applications. An extra complication is localisation 
over time needed in ensemble smoothers like the Ensemble Kalman Smoother and 4DEnsVar
when the fluid flow is strong: what is local at observation time is not necessary local
at the start of the assimilation window because the observation influence is advected with the flow.
Furthermore, the recent surge of papers on accurate treatment of observation errors shows that a long way
is still ahead of us to solve even the (close to) linear data-assimilation problem.

Although these problems are formidable, another difficulty arises from the fact that the problem is typically nonlinear,
and, with increasing model resolution and more complex observation operators,  increasingly so.
Both variational and Kalman-filter-like methods have difficulty handling nonlinear problems.
Variational methods can easily fail when the cost function is multimodal, and are hampered by the assumption
that the prior probability density function (pdf) of the state is assumed to be Gaussian. 
Ensemble Kalman filters make the explicit assumption that the prior pdf and the likelihood of the observations as function of the state
are Gaussian, or, somewhat equivalently, assume that the analysis is a linear combination of prior state and 
observations.
Both methods have been shown to fail for nonlinear data-assimilation problems in low-dimensional systems, 
and both have been reported to have serious difficulties in numerical weather prediction at the convective scale
where the model resolution is only a few km. Particle filters hold the promise of fully nonlinear data assimilation 
without any assumption on prior or likelihood,
and recent text books like <cit.>, <cit.>, and <cit.> 
provide useful introductions to data-assimilation in general, and particle filters in particular.

Other fully nonlinear data-assimilation methods are Markov-Chain Monte-Carlo methods that draw directly from the posterior in a sequential way, so one sample
after the other, after a burn-in period, see e.g. <cit.>, or <cit.> for a geophysics-friendly 
introduction. The samples are correlated, often 100% when the new sample is not accepted, making them
very inefficient in high-dimensional systems. This is why we concentrate on particle filtering here.

The standard or bootstrap particle filter can be described as follows.
The starting point is an ensemble of size N of model states _i^n ∈^N_x, called particles, that represent the prior 
probability density function (pdf) p(^n), as:


    p(^n) ≈∑_i=1^N 1/Nδ(^n-_i^n)


Between observations, each of these particles is propagated forward from time n-1 to time n with the 
typically nonlinear model equations


    ^n = f(^n-1) + ^n


in which f(..) denotes the deterministic model, and ^n is a random forcing representing missing physics,
discretisation errors, etc. In this paper we assume this model noise to be additive, but one could also
consider multiplicative noise in which ^n is a function of the state of the system.
We assume that the pdf from which the ^n are drawn is known; typically a Gaussian N(0,).

At observation times the true system is observed via:


    ^n = (^n_true) + ^n


in which the observation errors ^n are random vectors representing measurement errors and possibly
representation errors. Again we assume that these errors have known characteristics, often Gaussian, so e.g. ^n ∼ N(0,).
These observations ^n ∈^N_y are assimilated by multiplying the prior pdf above with the likelihood of each possible state,
i.e. the probability density  p(^n|^n) of the observation vector given each possible model state,
following Bayes Theorem:


    p(^n|^n) = p(^n|^n)/p(^n)p(^n)


in which p(^n|^n) is the posterior pdf, the holy grail of data assimilation.
To avoid confusion, it is good to realise that the true state is not a random variable when we apply
Bayes Theorem. It is a realisation of a process, which could be random or deterministic, from which we then take 
noisy observations. Instead, Bayes Theorem is a statement of what we think the true state might be.
Since the pdf of the ^n is known and Bayes Theorem is a statement for each possible
state ^n to be the true state, p(^n|^n) is the pdf of ^n given that the true state vector would be ^n. 
In general, since for a given state ^n the observation ^n is equal to the observation error  shifted by (^n), we find (see e.g. <cit.>):


    p(^n|^n) = p_(^n-(^n))



If we insert our particle representation of the prior into this theorem we find:


    p(^n|^n) ≈∑_i=1^N w_i δ(^n-_i^n)


in which the particle weights w_i are given by:


    w_i^n = p(^n|_i^n)/Np(^n) = p(^n|_i^n)/N∫ p(^n|^n)p(^n) d^n≈p(^n|_i^n)/∑_j p(^n|_j^n)


Since all terms are known explicitly we can just calculate this as a number.
The self-normalisation in the last part of (<ref>) is consistent with the notion that for a proper representation of
a pdf the sum of the weights should be equal to one, so that the integral over the whole state
space of the particle representation of the pdf is equal to one.
Figure <ref> depicts the working of this filter.

    

Propagating the particles _i^n to the next observation time n +1 gives a weighted 
representation of the prior at time n + 1. Assimilating the observation at time n + 1 by Bayes Theorem
leads to a modification of the weights (see e.g. <cit.> or <cit.>):


    w_i^n+1 = w_i^n p(^n+1|_i^n+1)/∑_j p(^n+1|_j^n+1)


Even in low-dimensional applications, the variation of the weights
increases with the number of assimilation steps. Eventually one
particle has a much higher weight than all the others. To prevent
this, resampling can be used before propagation to obtain equally
weighted particles. This duplicates high-weight particles and abandons low-weight particles.
After resampling, some of the particles have
identical values, but if the model contains a stochastic component
and independent random forcings are used for different particles,
diversity is restored. See e.g. <cit.> or <cit.> for details.
Algorithm <ref> illustrates the steps.


A simple resampling scheme using only one draw from a uniform distribution U is presented in Algorithm <ref>.


In high-dimensional problems the weights vary enormously even at one observation time, and typically one particle 
obtains a much higher weight than all the others.
<cit.> have shown that the number of particles needed
to avoid weights collapse, in which one particle gets weight 1 and the rest weights very close to zero, 
has to grow exponentially with the dimension of the observations  for a large class of particle filters.
If the weights collapse, all particles are identical after resampling, and all diversity is lost. From this discussion it becomes clear that for particle filters to work we need to ensure that
their weights remain similar.

In this review we will discuss four basic ways to make progress on this fundamental problem
of weight degeneracy.
In the first one, we explore the so-called proposal-density freedom to 
steer particles through state space such that they obtain very similar weights, see e.g. <cit.>. 
As pointed out by e.g.  <cit.> there 
are fundamental problems when applying these techniques to the high-dimensional 
geoscience applications. We will examine the issue in detail and discuss so-called equal-weight particle filters, 
which point towards new ways to formulate and attack the degeneracy problem. 

The second approach transform the prior particles into particles from 
the posterior, either in one go, or via a more smooth transformation process, see <cit.>. 
While the one-step approaches can be shown to fail in high-dimensional settings, they do lend themselves 
very naturally to localisation. The more smooth multi-step transition variants
seem to be able to avoid the degeneracy problem without localisation, and are
an interesting new development.

The third, more straightforward from the geoscience experience, approach is to introduce localisation in particle filters.
While initial implementations were discouraging (e.g. Van Leeuwen, 2009), new formulations have shown remarkable successes,
such that localised particle filters are now tested in global operational 
numerical weather prediction systems <cit.>.

The fourth approach is to abandon the idea of using pure particle filters and 
combine them with Ensemble Kalman Filters. This should not be confused with 
using Ensemble Kalman Filters in proposal densities. Several variants exist, 
such as second-order exact filters, in which only the first two moments are 
estimated, sequential versions in which first an EnKF is used and the 
posterior EnKF ensemble is used as input for the particle filter, or vice versa, 
and combinations in which localised weights are calculated and dependent on the effective ensemble
size a full particle filter, an EnKF, or a combination of both is used. 

These four variants form the basis of the following four chapters. Each chapter contains a 
critical discussion of the approximations and remaining major issues.  It should be noted that the
pseudo code provided does not give the most efficient implementation of the different
particle filters, but is rather an illustration of the computational steps involved. Efficient pseudo code for some of the more complex schemes can be found in <cit.>.
The paper is closed with
a concluding section and an outlook of what possible next steps could be.





§ PROPOSAL DENSITY PARTICLE FILTERS



Ideally we draw independent samples directly from the posterior pdf because the samples would all have equal weight
automatically. This can only be done, however, when the shape of the posterior pdf is known and when it is easy to
draw from the posterior. An example of this is a Gaussian prior combined with a linear Gaussian likelihood. Under these assumptions the
posterior is also Gaussian and the mean and covariance can be calculated directly from the prior using the Kalman update
equations. Ensemble Kalman filters make use of this result and draw directly from that pdf, which is why all posterior particles have equal
weights in an Ensemble Kalman Filter.

The standard particle filter draws particles from the prior. These then have to be modified to become particles 
of the posterior via the weighting with the likelihood. This is a general
procedure in statistics called importance sampling: one draws from
an approximation of the pdf one is interested in, and corrects for this 
via so-called importance weights. 

In the introduction we argued that drawing from the prior leads to weights that vary too much: typically, 
in high-dimensional problems with numerous independent observations one particle gets weight 1, and
all other particles have a weight very close to zero.
However, we could explore the idea of importance sampling on the transition from one time to the next.
When the numerical model is not deterministic but stochastic we have the freedom 
to change the model equations to move the particles to those parts of state space where we want them to be,
for instance closer to the observations. 

Mathematically this works as follows. Assume we have observations at time
n, so Bayes Theorem at time n is given by (<ref>).
If the model is stochastic, we can write the prior as


    p(^n) = ∫ p(^n|^n-1) p(^n-1) d^n-1


where p(^n|^n-1) is the transition density, the pdf of the state at
time n when the state at time n-1 is known. For
instance, if the model error is additive and the model equation is given by (<ref>),
it holds that 

    p(^n|^n-1) = p_(^n-f(^n-1)).


Often the model errors are assumed to be Gaussian
∼ N(0,), and we find


    p(^n|^n-1) = N(f(^n-1),).


but the method is more general than that.

Assume now that at time n-1 we have a set of weighted particles as in (<ref>),
but with weights w_i^n-1 instead of 1/N.
We can evaluate the expression (<ref>) for the prior as
a weighted mixture of transition densities


    p(^n) ≈∑_i=1^N w^n-1_i p(^n|_i^n-1)


In the following we neglect the approximation error at time n-1
and assume that (<ref>) is exact. 
This is not necessarily a good approximation, especially when the number of particles 
is small. On the other hand, it is consistent with the particle filter approximation in the first place,
and one of the few things one can do.
By Bayes formula
(<ref>), the posterior can then be written as:


    p(^n|^n) ≈∑_i=1^N w^n-1_i p(^n|^n)/p(^n)p(^n|_i^n-1)



In the standard particle filter one makes one draw from
p(^n|_i^n-1) for each i, 
and we know that this leads to ensemble collapse for high-dimensional systems.
However, now the prior particles at time n are allowed to arise from
following a different model equation. This works as follows. 
We can multiply and divide equations (<ref>)
and (<ref>) by a so-called proposal density
q(^n|^n-1,^n), leading to:


    p(^n) ≈∑_i=1^N w^n-1_i p(^n|_i^n-1)/q(^n|_i^n-1,^n) q(^n|_i^n-1,^n)
 
and

    p(^n|^n) ≈∑_i=1^N w^n-1_i p(^n|^n)/p(^n)p(^n|_i^n-1)/q(^n|_i^n-1,^n) q(^n|_i^n-1,^n)


where q(^n|_i^n-1,^n) should be non-zero whenever p(^n|_i^n-1) is.
This step is completely general. 

Now realise that drawing from p(^n|_i^n-1)
corresponds to running the original stochastic model.
We could instead draw from q(^n|_i^n-1,^n),
which would correspond to a model equation from our choosing.
Figure <ref> illustrates the basic idea.

    


For instance when the original model is given by
(<ref>), we can use 


    ^n = g(^n-1,^n) + ^n


in which g(.,.) is now the deterministic part and ^n is the
stochastic part. These can be freely chosen, and examples of these will be given below.
Note that we allowed g(..) to depend on the observations at the future
time. This means that we generate the prior particles at time n by making
one draw from q(^n|^n-1_i,^n) for each i where


    q(^n|^n-1,^n) =p_(^n-g(^n-1,^n))


In general, we draw the particles at time n from the 
alternative model q(^n|^n-1,^n) and account for this by
changing the weights of the particles. Equations (<ref>)
and (<ref>) can be written as
  
    p(^n) = ∑_i=1^N ŵ^n-1_i q(^n|_i^n-1,^n)


and

    p(^n|^n) = ∑_i=1^N ŵ^n_i q(^n|_i^n-1,^n)

where the weights are given by:


    ŵ^n-1_i  ∝ w^n-1_i p(_i^n|_i^n-1)/q(_i^n|_i^n-1,^n) .

and


    ŵ^n_i  ∝ŵ^n-1_i p(^n|_i^n)/p(^n)∝ w^n-1_i p(^n|_i^n) p(_i^n|_i^n-1)/q(_i^n|_i^n-1,^n).

Here the coefficients of
proportionality ensure that the weights sum to 1. 
In a reinterpretation of these equations, if _i^n is drawn from the
alternative model q(^n|^n-1_i,^n) we can also write

    p(^n) ≈∑_i=1^N ŵ^n-1_i δ(^n-_i^n)

and

    p(^n|^n) ≈∑_i=1^N ŵ^n_i δ(^n-_i^n).



We see that the weights now contain two factors, the likelihood weight, which also appears in the
standard particle filter, and a proposal weight. These two weights have opposing effects.
If we use a proposal density that strongly pushes the model towards the observations, the likelihood weight
will be large because the difference between observations and model states becomes smaller, but the proposal weight 
becomes smaller because the model is pushed away from where it wants to go, so p(^n|^n-1_i)
will be small. On the other hand, a weak pushing towards the observations keeps the proposal weight high,
but leads to a small likelihood weight. This suggests that there is an optimum weight related to an optimal
position _i^n for each particle as function of its position at time n-1. This will be explored
in equal-weight formulations of the particle filter.
Figure <ref> shows how typical proposal-density particle filters work. Equal-weight particle filters are discussed later.

    




 §.§ A simple relaxation scheme

To illustrate the idea of a proposal density we consider the following simple example.
We could add a relaxation or nudging term
to the original equation to steer the particles towards the observations and
make their weights more similar, as pioneered by <cit.> for geoscience 
applications. The model equation is written as:


    ^m = f(^m-1) + (^n-(^m-1))+ ^m


where we used time index m for the state vector to emphasise that there are several model time steps between observation times.
 is a relaxation matrix of our choice.
In this example, the deterministic part consists of the first two 
terms on the right-hand side of the equation, while the third term denotes the random part.
Let's assume the pdf of the random forcing is Gaussian with mean zero and covariance . Then we can immediately 
write for the proposal density


    q(^m|^m-1,^n) = N(f(^m-1)+(^n-(^m-1)),)


since the pdf of ^m is just a shift in the mean of the pdf of ^m.
For the original model, we assume that the random part is Gaussian with zero mean and covariance ,
so that


    p(^m|^m-1) = N(f(^m-1),)



The change in the model equations is compensated for in particle filters by a change in the relative weight
of each particle, and the expression for this change in weight for this case is:


    w_i^m      =      w^m-1_i p(_i^m|_i^m-1)/q(_i^m|_i^m-1,^n)
       ∝    w_i^m-1exp[ -J_p]/exp[ - J_q
                 ]

in which, for Gaussian model errors,


    J_p =1/2(_i^m-f(_i^m-1))^T^-1(_i^m-f(_i^m-1))


and 

    J_q     =    1/2(_i^m-f(_i^m-1)-(^n-H(^m-1))
                )^T ·
          ^-1(_i^m-f(_i^m-1)-(^n-(^m-1))) 
        =    1/2  (_i^m)^T  ^-1_i^m


Note that the normalisation factors of the Gaussians do not have to be calculated explicitly
if we use that the sum of the weights has to be equal to one. The scheme is depicted by Algorithm <ref>.



Simple as the scheme is, it does not solve the degeneracy problem. However, it can be used
as a simple scheme when several model time steps are used between observation times,
because the proposal is independent of the proposal at other time steps.
This scheme an easily be used in combination with other schemes that work at observation time,
to be discussed next.




 §.§ Weighted Ensemble Kalman Filter

One could also use other existing data-assimilation methods in proposal densities, like
Ensemble Kalman filters or variational methods. In the Weighted Ensemble Kalman filter
<cit.> the stochastic EnKF of <cit.> is used as follows. The Ensemble Kalman Filter 
update can be written as:


    _i^n = _i^f + (^n - _i^f - _i)


in which _i^f = f(_i^n-1)+ _i^n, the matrix K is the ensemble Kalman gain
and _i ∼ N(0,), with  the observational error covariance.
Using the expression for the forecast _i^f in the Kalman filter update equation we find:


    _i^n = f(_i^n-1) + (^n -  f(_i^n-1))+ (-)_i^n - _i


which we can rewrite as the sum of a deterministic and a stochastic part as:


    ^n = g(^n-1,^n) + _i^n


identifying g(^n-1) = f(_i^n-1) + (^n -  f(_i^n-1)) and _i^n = (-)_i^n - _i.
Therefore, we find for the proposal density:


    q(^n|_i^n-1,^n) = N(f(^n-1)+(^n- f(^n-1),)


with 


    = (-)(-)^T + ^T .


Strictly speaking, this is correct only if the Kalman gain is calculated using the ensemble covariance
of f(^n-1), so without the model errors ^n, otherwise the proposal is not Gaussian.
We can calculate the weights of the particles in a similar way as in the previous example. 
Algorithm <ref> shows the algorithmic steps.



The behaviour of this filter has been studied extensively in <cit.>. In
high-dimensional systems this 
filter will be degenerate, consistent with the theory of <cit.>, and
as proven in the next section.
The only way to make this work is to include localisation, not only at the EnKF level, but also
at the level of the particle filter, see e.g. <cit.>. 




 §.§ Optimal proposal density



In the class of particle filters in which the proposal density of each particle is dependent on only
that particle, an optimal proposal density can be derived, as e.g. shown in <cit.>. They defined optimality 
as the proposal density that gives a minimal variance of the weights, and
<cit.> provide an elegant proof of this optimality.
In this section we generalise this result and show that the
optimal proposal density is  
optimal even when each particle has its own proposal density which is
allowed to depend on all previous particles,
so a proposal of the form q(^n|i,_1:N^n-1,^n).

<cit.> concentrate on the case that one is interested in
an optimal  representation of p(x^n,x^n-1|y^n) in a sequential
algorithm, so in a sequential smoother. 
To this end they introduce the random variable 


    w^*(^n,^n-1) = p(^n,^n-1|^n)/q(^n,^n-1|^n)


and determine that proposal density q that minimises the variance in the weights w^*,
with the expectation taken over the density from which we draw the particles, so the proposal q.

Here we show that the optimal proposal density is also optimal for the strict filtering case,
so when we are interested in minimal variance of the weights at time n only. 
Specifically, the question is: given the set of particles at t=n-1 drawn from p(_n-1|_1:n-1),
which proposal density of the form q(^n|i,_1:N^n-1,^n) gives minimal variance of the weights at time n?

Using Bayes formula, we can write the expression 
for the weight of particle i as function of the state at time n as:


    w_i^n=w_i(_i^n)     =    p(^n|_i^n)/N p(^n) p(_i^n|_i^n-1) /  q(_i^n|i,_1:N^n-1,^n)
        =    p(^n|_i^n-1)/N p(^n) p(_i^n|_i^n-1,^n) /  q(_i^n|i,_1:N^n-1,^n)


where we assume, without loss of generality, an equally weighted ensemble at time n-1.
Note that the second equality follows from Bayes Theorem, as follows:


    p(_i^n|_i^n-1,^n)     =    p(^n|_i^n,_i^n-1)/p(^n|_i^n-1) p(_i^n|_i^n-1) 
        =    p(^n|_i^n)/p(^n|_i^n-1) p(_i^n|_i^n-1)




Consider the pair of random variables (I,^n) where Prob(I=i) =
1/N and, conditionally on I=i, ^n ∼ q(^n|i,_1:N^n-1,^n).
Furthermore, define the associated random variable


    W=w_I(^n) = p(^n|_I^n-1)/Np(^n) p(^n|_I^n-1,^n) /q(^n|I,_1:N^n-1,^n)


where 


    p(^n) = 1/N∑_j=1^N p(^n|_j^n-1)


In order to find the proposal q that minimizes the variance of W,
we use the well-known law of total variance (derived in the appendix for completeness):


    var_W(W) = var_I(E_^n|I(W)) + E_I(var_^n|I(W)).


First, we see that, under the proposal q:


    E_^n|I(W) = p(^n|_I^n-1)/Np(^n)∫ p(^n|_I^n-1,^n)
    d^n = p(^n|_I^n-1)/Np(^n)


is independent of q. Moreover, E_W(W)=E_I(E_^n|I(W))=1/N and thus
the first term in var_W(W) is


    1/N∑_i p(^n|_i^n-1)^2/N^2p(^n)^2 -
    1/N^2
    =1/N∑_i (p(^n|_i^n-1)/Np(^n) -
      1/N) ^2 ≥ 0.


For the second term
we use that var_^n|I(W)≥ 0 with equality 
if and only if W is almost surely constant in ^n,
that is if and only if


    p(^n|_i^n-1,^n) /q(^n|i,_1:N^n-1,^n) = cst(i,_1:N^n-1,^n).


in which cst(..) is this constant which can depend on other
    variables than ^n.
Because both p and q are densities (in ^n), cst=1. Combining
these results, we have a lower bound for var(W) that is determined by the variance
of p(^n|_i^n-1) over i, with equality if and only if


    q(^n|i,_1:N^n-1,^n)=p(^n|_i^n-1,^n)


Note that this is a new result as previous proofs only considered proposal densities of the form q(^n|_i^n-1,^n),
and we extended it to more general proposal densities of the form  q(^n|i,_1:N^n-1,^n).
  
This remarkable result shows that firstly the optimal proposal density, so 
p(^n|_i^n-1,^n), does indeed lead to the lowest variance in the weights 
for the class of particle filters in which the transition density is of the form q(^n|i,_1:N^n-1,^n).
Secondly, it shows that we can predict the variance in the weights without doing the 
actual experiment, for any number of particles, provided we can compute
p(^n|^n-1_i), and thirdly the weights are independent of the position
of the particles ^n. Unfortunately, this variance is zero only 
when the observations are not dependent on the state at time n-1, which
is never the case in the geosciences.

A simple case where we can compute both the optimal proposal density and
the weights p(^n|^n-1_i) is when p(^n|^n-1_i) is given
by (<ref>) and the observation operator H = is linear. By
the same argument that 
is used to derive the Kalman filter update, we find


    p(^n|^n-1_i,^n) =


    = N(f(^n-1_i) + 
    (^n - f(^n-1_i)),(-^T)),


where = ^T(^T+)^-1 is the Kalman-like gain with the background covariance , and
the weights are proportional to:


    p(^n|^n-1_i) = N( f(^n-1_i), ^T + )


This shows two things: First, in this special case, the simple relaxation
scheme of Section 
2.1 is equal to the optimal proposal when the relaxation matrix  is chosen as 
above. Second, comparing the weights of the optimal proposal
with the weights of the standard filter, they both depend on the squared distance
||^n -  f(^n-1_i)||^2, and ||^n - ^n_i||^2, respectively, but in the standard particle filter the
distance is defined w.r. to  and in the optimal proposal the distance 
it is defined is w.r.
to ^T +. Hence the weights with the optimal proposal are more similar,
but the improvement is substantial only if  is large, and the analysis
of  weight collapse by <cit.>  still applies. 

One can extend the optimal proposal density idea to more than one time step.
<cit.> show that the optimal proposal is the proposal of this form with minimal variance in the 
weights in this case too, which can also easily be seen by applying the above to


    W=w_i(^n) = p(^n|_i^m-1)/Np(^n) p(^n|_i^m-1,^n) /q(^n|_i^m-1,^n)


for m<n.

Looking back at the filters described in the previous sections we find the following.
The relaxation scheme uses a simple proposal density that is of the form q(^n|_i^n-1,^n),
so the theory holds, and that proposal will lead to degenerate results.
This is indeed the finding of <cit.>.
The Weighted Ensemble Kalman Filter has a proposal that depends on all particles at time n-1 through
 the Kalman gain , so the proposal is of the form q(^n|i,_1:N^n-1,^n).
Hence also this filter will perform worse than the optimal proposal and hence will be degenerate for high-dimensional 
systems. This was first explored in detail by <cit.>.





 §.§ Implicit Particle filter 


The Implicit Particle Filter is an indirect way to draw from the optimal proposal, even over several time steps.
Often the assumption is made that the model errors of both original model and proposal density
are Gaussian, and the observation operator  is linear. In this case, a draw from the optimal
proposal is a draw from a multivariate Gaussian, and we know how to do that.

However, when  is nonlinear, or when the proposal is used over several model time steps 
the density to draw from is not Gaussian anymore. <cit.> realised that one could
still draw from a Gaussian and then apply a transformation to that draw to find samples from the 
optimal proposal density. 
The method is explained here for one time step, but the extension to multiple time steps is straightforward.
Figure <ref> illustrates the basic idea.

    

As mentioned in Sec. <ref> on the proposal density the posterior pdf can be written as:


    p(^n|^n) = ∑_i=1^N w^n-1_i p(^n|^n)/p(^n)p(^n|_i^n-1)/q(^n|_i^n-1,^n)q(^n|_i^n-1,^n)


The scheme draws from a Gaussian proposal q() = N(0,), and we can write the transformation as 
q(^n|_i^n-1,^n) = q() _i^-1 in which _i is the Jacobian of the transformation from ^n to .
That transformation is found implicitly, hence the name of the filter, by defining


    F_i(^n) = - log[ p(^n|^n)p(^n|_i^n-1) ]


 and, after drawing _i for each particle, solving for ^n in


    F_i(^n) =1/2_i^T _i + ϕ_i


for each particle, in which ϕ_i = min_^n F_i(^n) ∝ p(^n|_i^n-1).
The weights of the particles become:


    w_i^n     =      w^n-1_i p(^n|_i^n)/p(^n)p(_i^n|_i^n-1)/q(_i^n|_i^n-1^n)
        =     w^n-1_i exp[-F_i(_i^n)]/exp[ -1/2_i^T _i ]_i 
        =     w^n-1_i exp[-F_i(_i^n)]/exp[ -F_i(^n) + ϕ_i ]_i 
        =     w^n-1_i exp[ - ϕ_i ] _i


Interestingly, while the optimal proposal density shows that the weights are only dependent on
the position of the particles at the previous time, so on _i^n-1 via ϕ_i, the
implicit map makes the weights also dependent on the positions at the current time n, so on _i^n
via the Jacobian of the transformation between  and . 
Only when the Jacobian is a constant, so when F_i is quadratic in _i, this
dependence disappears.


Solving (<ref>) is not straightforward in general.
<cit.> suggest a random map of the form


    _i^n = _i^a + λ_i(_i) ^1/2_i


in which  is a chosen covariance matrix, ideally the covariance of the posterior pdf, 
_i^a = min F_i(^n) and λ_i is a scalar. This transforms the problem into 
solving a highly nonlinear scalar equation for λ_i, which is a much simpler problem than 
finding _i^n directly. This map can be shown to be a bijection when F_i(_i^n) has only closed
contours in the high-probability regions; otherwise one would have to first choose a closed 
contour area and then perform the map. 
In general, when the optimal proposal (over several time steps if needed) is multimodal, the transformation from 
the state variable to a Gaussian is not monotonic, and the Implicit Particle Filter needs to be adapted, e.g. by
using a separate Gaussian for each mode.
The algorithm is given in Algorithm <ref>.



Of further interest is that _i^a is the same as the solution to a 4DVar problem well known in meteorology.
But it is a special 4DVar 
as the initial position of each particle is fixed and it has to be a weak-constraint 4DVar.
The latter condition is needed as a strong-constraint 4Dvar would have no possibility to move a 
particle in state space as its initial condition is fixed. 

However, also this filter will suffer from weight collapse in high-dimensional applications
as it is still a sampling scheme for the optimal proposal density. 
The following sections will discuss ways to improve on the optimal proposal.




 §.§ Equal weights by resampling at time n-1

As noted already in equation (<ref>), we can write equation
(<ref>) as

    p(^n|^n)    =   ∑_i=1^N w^n-1_i
      p(^n|^n-1_i)/p(^n) p(^n|^n-1_i,^n) 
        =    ∑_i=1^N α_i p(^n|_i^n-1,^n)


where


    α_i =  w^n-1_i  p(^n|_i^n-1)/p(^n)


This says that, assuming the pdf at the previous time can be approximated by
a set of N particles, the analysis distribution is a mixture of the
optimal proposal pdf's p(^n|_i^n-1,^n) with mixture weights
α_i. 

If we can compute the optimal proposal density and the weights α_i
in closed form, 
we can also draw samples directly from this mixture density. For this,
we first draw an index I  from the discrete distribution
with weights α_i, Prob(I=j) = α_j, followed by a draw from
the corresponding  
pdf p(^n|_I^n-1,^n). Doing this N times will lead to N
different particles with equal weights because each of them 
is an independent draw directly from the posterior. If the index I 
is equal to a value j more than once, the particle ^n-1_j 
is propagated from time n-1 to time n with independent random forcing
for each of these draws. This simple scheme provides better
samples than the optimal proposal density because all particles are different
at time n by construction.

However, this does not solve the problem of weight collapse because drawing
the index I is nothing else than resampling the particles at time n-1
with weights proportional to w^n-1_i  p(^n|_i^n-1). If
w^n-1_i=1/N, the variance of these weights is exactly equal to
the lower bound that we found in Section 2.3. The main difference is that
the collapse now happens at time n-1. The only advantage is that 
all particles will be different at time n.


If we cannot compute the optimal proposal density and the weights
α_i in closed form, we can still use the importance
sampling idea to draw from the mixture p(^n|^n) 
by drawing pairs (I,^n) consisting of an index I and a
state ^n at time n. We choose a proposal distribution 
β_i = β_i(^n) for the index and proposal distributions
q(^n|^n-1_i,^n) for the state. Then we draw the index I_i
with Prob(I_i=j)=β_j(^n) and conditionally on I_i=j we
draw ^n_i from q(^n|_j^n-1,^n). Finally, we compute 
weights w^n_i by

    w^n_i ∝w^n-1_j p(^n_i|^n-1_j)p(|^n_i)/β_j(^n) q(^n|_j^n-1,^n) if  I_i=j

The particles ^n_i with weights w^n_i provide the desired
approximation of p(^n|^n) whereas the indices I_i can be
discarded after the weights have been computed. We could produce
an evenly-weighted approximation by a further resampling step, or
take the weights w^n_i into account during the next iteration.

In this approach we can obtain equal weights w^n_i by choosing

    q(^n|_j^n-1,^n)= p(^n|_j^n-1,^n)

and 

    β_i(^n) ∝ w^n-1_i p(^n|^n-1_i).

With this choice, we draw directly from the mixture (<ref>).
As mentioned before, although the weights w^n_i are then equal to 
1/N, the algorithm contains a hidden weighting and resampling
step of particles at time n-1. It thus remains susceptible to weight
collapse in high dimensions.  

This approach of using importance sampling for the joint distribution of
(I,^n) is due to <cit.> who called it “Auxiliary Particle
Filter” (the index I is an auxiliary variable that is discarded at the
end). They discuss, in addition, approximations of the
optimal proposal density and the optimal weights α_i. One of their
suggestions is to use for the index I the proposal with weights

    β_i ∝ w^n-1_i p(^n| ^n_i)

where ^n_i is a likely value of the distribution
p(^n|^n-1_i), e.g. the mean or median or simply a draw from it.
Typically, ^n_i is found by a probing step where particles at time
n are propagated by a simplified model, e.g. by omitting
stochastic terms or with simplified subgrid-scale parameterisations or 
thermodynamics. If I_i=j and the state ^n_i at time n is proposed from
p(^n| ^n-1_j), the weights become

    w^n_i ∝p(|^n_i)/p(^n| ^n_j)

They will vary less provided ^n_i is close to ^n_j, i.e.
provided the simplified model is a good approximation to the full model
and the stochastic part of the full model is small. 





 §.§ The Equivalent-Weights Particle Filter (EWPF)

The EWPF <cit.> uses the idea 
to obtain a more evenly weighted set of particles by not sampling
from the exact posterior, but allowing for a small error. It
starts with determining the weight of each particle at the mode of
p(^n|_i^n-1,^n)
for each particle i, w_i^max∝ p(^n|_i^n-1).
Note that these weights are equal to the weights obtained in the optimal proposal density.
In the optimal proposal density case the weights do not depend on the position ^n of the particle,
but note that the proposal used here will be different.

The particles are not moved to these modes, but the weights are used to define a target weight.
This target weight w_target is chosen such that a certain fraction ρ of particles can reach that weight.
To this end we sort the weights in magnitude from high to low in an array w^*_i,  i = {1,2,...,N]} and set  w_target = w^*_N*ρ.
For instance, with 100 particles and a fraction of ρ = 0.8 we would find w_target = w^*_80.

The next step is to find a position in state space for each particle that can reach this weight
such that its weight is exactly equal to the target weight. This means we solve for ^n in


    w_i(^n) = w_target


for each particle i that can reach this weight. There are many solutions
of this equation, but we choose the one which is on the line through
^a_i and f(^n-1_i) and is closest to f(^n-1_i).
Denote this position as _i^*. Note that this is purely deterministic move, so a stochastic part still has to be added.
The final position of these particles is then determined by adding a very small random perturbation 
from a chosen density, so


    _i^n = _i^* + _i^n


This stochastic move  ensures that the proposal has full support and is not
a delta function centred at _i^*. The density of _i should on
the one hand have most of its mass concentrated around 0 in order 
not to change the weights of the particles too much, and on the other hand
it should be relatively constant since we divide by the value of the
proposal density. Both requirements cannot be fulfilled exactly, but
we can take some error in the sampling into account and choose 
a narrow uniform distribution. The scheme is depicted in Algorithm
<ref> for the special case that of Gaussian model errors and a linear observation operator.
If these conditions do not hold, one will typically need iterations 
to solve for a_i and b_i.




It is common knowledge, see e.g. <cit.>, that the proposal should be wider or at least as wide as the
target, while the width of the stochastic part of the proposal is chosen very small here.
The reason that we can do this is that the position of the centres of these proposal densities
are typically further away from the observations than e.g. in the optimal proposal
because the target weight forces particles away from their optimal positions, so away from the observations.
This means that the deterministic moves of the particles ensure a large spread in the full proposal.

A formal way to avoid such an error has been described by 
<cit.>. They choose the proposal to be a mixture of a uniform density and a Gaussian which is also used in Alg. <ref>. 
Both have small variance, and the mixture coefficient of the uniform density is chosen to be much larger 
than that of the Gaussian. This means that drawing from the Gaussian and also drawing from its tails becomes
highly unlikely. In practice, since we always work with small ensemble sizes the chance of filter degeneracy
by drawing from the Gaussian, and then drawing from the tail of the Gaussian is indeed highly unlikely.

Finally, the full weights for the new particles are calculated and the whole ensemble is resampled, including
those particles that were unable to reach the target weight. Because of the target-weight construction 
the weights of the particles are very similar, and filter degeneracy is avoided.
This filter has been used in a reduced-gravity ocean model by <cit.>,
and in the same system studied for the gravity-wave production by the scheme in <cit.>.
It has also been applied in a climate model by <cit.>.

To analyse the scheme further, we can again look at the variance of the weights.
For this it is important to note that this scheme does not see the weight of a particle as a 
function of the state  and particle index I, but rather the state as function of the weight W
and index I, so (W,I). 
Specifically, W|I has values in two ranges. For the particles with I=i
that can reach the target weight we find w|I = w_target + ϵ_i in which ϵ_i is a small perturbation 
from the target weight due to the small stochastic move discussed above.
For those particles that cannot reach the target weight their weights are very close to zero.
So we find:


    E_I[W] ≈ρ (w_target+ϵ̅)+(1-ρ)0 = ρ (w_target+ϵ̅)


in which ϵ̅=E_I[ϵ].
If H is linear and the errors in the observations and the model equations are Gaussian we find
ϵ̅=0, but if any of these three conditions does not hold this is not necessarily so.
However, we do know that by construction |ϵ̅|<<1.
Since the sum of the weights should be equal to 1, we find that w_target≈1/(Nρ),
and hence E_I[W]=1/N, as expected.
Furthermore


    var_I(W)     =    ρ∑_i=1^ρ N(w_target+ϵ_i)^2  - (ρ w_target)^2  
       ≈   1/N^21-ρ/ρ

 
This expression shows that the variance in the weights ranges between
0 for ρ=1, so when all particles are kept, to (N-1)/N^2 ≈ 1/N for ρ=1/N,
so when one particle is kept. 
We can compare this with the optimal proposal when the number of independent observations is large.
In that case one particle will have a weight very close to one, and the rest will have weights very
close to zero. The variance in the weights is then (N-1)/N^2 ≈ 1/N, indeed equal to
the ρ=1/N case in the EWPF scheme, as expected. The EWPF can, however, reduce that
variance, even to zero, depending on the choice of the tuning parameter ρ. 

When this tuning parameter is chosen close to one, the target weight will be low, and hence
particles will be moved further away from the mode of the optimal proposal density.
In practise this means that the particles are pushed further away from each other, leading to
a wider posterior pdf. A small value for the fraction will have the opposite effect.  
Since we do not know a-priori what the width of the posterior should be, this is a clear drawback of this method.
We will come back to this later.




 §.§ The Implicit Equal-Weights Particle Filter


In the Implicit Equal Weights Particle Filter (IEWPF) we set the target weight equal to the minimum of the optimal
proposal weights for all particles. Then, the position of each particle is set to the mode of the optimal proposal 
density plus a scaled random perturbation. The scale factor is chosen such that the weight of each particle is equal to the target weight.
Note that in the standard setting no resampling is needed, but see <cit.> for other possibilities.

The implicit part of the scheme follows from drawing samples implicitly from a standard Gaussian distributed proposal 
density q() instead of the original q(^n|^n-1,^n), following the same procedure as in the Implicit Particle Filter.
We define a relation


    _i^n = _i^a + α_i^1/2^1/2_i^n


where _i^a is the mode of p(^n|_i^n-1,^n),  is a measure of the width of that pdf, _i^n ∈^N_x is a standard Gaussian-distributed random
vector, and α_i is a scalar. 

The IEWPF scheme is different from the Implicit Particle Filter in that it chooses the α_i such that all particles get the same weight w_target, so the scalar α_i is determined for each particle from:

    w_i(α_i)=   p(_i^n|_i^n-1,^n)p(^n|_i^n-1)/Np(^n)q(^n|i,_1:N^n-1,^n)   = w_target

This target weight is equal to the lowest weight over all particles in an optimal proposal.
This ensures that the filter is not degenerate in systems with arbitrary dimensions and an arbitrary number of independent observations.
The resulting equation for each α_i is nonlinear and complex because it will contain the Jacobian of the transformation
from ^n to ^n, similar to the Implicit Particle Filter. The Jacobian will contain
the derivative of α_i to _i, which is the main source of the complexity in this scheme.
Algorithm <ref> depicts the scheme for the case of a linear observation operator.
A nonlinear observation operator will lead to more complicated equations for the α's.



The scheme is similar to the optimal proposal density using the Implicit Particle Filter by first determining the mode
of the proposal and then adding a random vector. The difference is that in the IEWPF the size of the vector is
determined such that the each particle reaches the target weight. It turns out that this construction excludes part
of state space for all but one particle. For each particle the excluded part is different, so the ensemble samples the 
whole space, but the individual particles do not.
Details of the method can be found in <cit.>.

Analysing the scheme in more detail, the proposal density used in this scheme is of one dimension lower than that
of the state itself. The direction of the random vector in state space is determined by the proposal density, but 
the size of the random vector is then determined deterministically, dependent on that direction. So the proposal
density misses one degree of freedom for all but one particle, the particle with the lowest weight that has α_i=1.
Although missing one degree of freedom in a very high dimensional system might seem acceptable it does lead
to a bias. 
Figure <ref> shows how the implicit equal-weights particle filter works.

    





 §.§ Discussion


We first note that the optimal proposal is only optimal in a very limited sense, as has been known a long time
with the invention of the auxiliary particle filter.
We have seen that it is not difficult to generate particle filters
that even have zero variance in the weights.
In the optimal proposal setting one forces Prob(I=i)=1/N, while the simple choice Prob(I=i)∝
p(y^n|x_i^n-1) leads to an equal-weight particle filter. 
Furthermore, schemes have been introduced that consider the state as function of 
the state at the previous time and the weight the state at the current time should obtain, so instead of 
working with W(,I) we choose (W,I), which 
opens up a whole new range of efficient particle filters in high dimensional systems. 

The EWPF and the IEWPF are by construction particle filters that are not degenerate in high-dimensional systems 
 and that do not rely on localisation.
However, it is easy to see that both filters are biased, or inconsistent. In the limit of an infinite number of particles
the target-weight constructions will prevent the schemes to converge to the full posterior pdf.
The schemes are only of interest when the ensemble size is limited. 
As long as the bias from the target-weight construction is
smaller than the Monte-Carlo error this bias is of no direct consequence. 
It will be clear that the number of possible methods
that have this property is huge, and much more research is needed to explore the best possibilities.






§ TRANSPORTATION PARTICLE FILTERS


In resampling particle filters the prior particles are first weighted to represent the posterior and then transformed to 
unweighted particles simply by duplicating high-weight particles and abandoning low-weight particles. 
In transformation particle filters one tries to find a transformation that moves particles from the prior
to particles of the posterior in a deterministic manner. A related approach, which uses random transformation steps, 
is based on tempering the likelihood, which we also discuss in this section.




 §.§ One-step transportation

In one-step transportation one tries to transform samples from the prior into samples from the posterior in 
one transformation step. An example is the Ensemble Transform Particle Filter <cit.>, in which the unweighted particles are linear combinations
of the weighted particles, so one writes:


    ^a = ^f


in which the matrix ^f = (_1^f,⋯,_N^f) and similar for ^a, and in which  is a transformation 
matrix. The only conditions on  are that d_ij≥ 0,  ∑_i d_ij=1 and ∑_j d_ij=w_i N. These three conditions
leave a lot of freedom for all N^2 elements of D, and a useful way to determine them is to ensure minimal
overall movement in state space of the particles from prior to posterior. This leads to an optimal transportation 
problem and is typically solved by minimizing a cost function that penalises movement of particles. 

We can see immediately that this method will not work when the weights are degenerate as the solution will be 
degenerate and all particles have no other choice than move to the prior particle with weight (close to) one. 
However, the strength of this filter is that it allows for localisation in a very natural way by making the weights,
and hence the matrix , space dependent. The method will be discussed in more detail in Section <ref> on localisation.
Here we provide the basic algorithm in Algorithm <ref>.




The ETPF provides a direct map from prior to posterior particles without explicitly constructing a transformation map. 
An alternative approach has been suggested in <cit.>, where an approximate transportation map  is constructed 
such that  belongs to certain family of maps and  is chosen such that the Kullbeck-Leibler 
divergence between the pdf generated by  and the posterior pdf is
minimized. See <cit.> for an 
efficient implementation in the context of filtering and smoothing for low-dimensional systems.





 §.§ Tempering of the likelihood

Instead of trying to transform the particles from the prior to particles from the posterior in one step
one can also make this a smoother transition. In tempering (<cit.>, see also <cit.> 
and  <cit.>) one factorises the likelihood as follows:


    p(|) = p(|)^γ_1... p(|)^γ_m


with 0<γ_i<1 and ensuring that the sum of the γ's is equal to 1.
Then the weighting of the particle filter is first done with the first factor, so


    p_1(|) = p(|)^γ_1/p()^γ_1p()


The reason for this is that the likelihood is much less peaked, and hence the degeneracy
can be avoided when γ_1 is small enough. 
Figure <ref> illustrates the basic idea.

    


The particles are resampled, and now the weighting is performed using the second factor,
followed by resampling, etc. In this way the scheme slowly moves all particles
towards the high-probability regions of the posterior.
Of course, after resampling several particles will be identical, 
so one needs to jitter the particles, so perturb them slightly, to regain diversity.

This jittering should be a move of the particles that preserves the posterior
pdf. It could be implemented as a Markov-Chain Monte-Carlo method with
the posterior as the target density, e.g. exploring resample-move strategies, see
e.g. <cit.>. A problem is, however, that in sequential filtering we only have a representation
of the posterior density in terms of the present particles, and this representation is 
very poor due to the small number of particles. Possible avenues are to fit a pdf of a certain
shape to the present particles, e.g. a Gaussian mixture model, and use that as target
density. 

A problem in the geosciences is that this posterior fit needs to preserve the delicate balances 
between the model variables that are present in each particle, and an extra complication is that these
balances can even be nonlinear. Also the transition kernel of the Markov Chain
should somehow preserve these balances. An example of its use in the geosciences 
is the Multiple Data Assimilations (MDA) method of <cit.>, in which
the intermediate pdf's are assumed to be Gaussian. See also <cit.> for 
a comparison of this method to other iterative implementations of the Ensemble Kalman Filter/Smoother.

If, however, one allows for model error in the model equations, the following scheme proposed by <cit.>
does not have this problem. In that case the prior at observation time can be written as (see equation
<ref>):


    p(^n) ≈1/N∑_i=1^N p(^n|_i^n-1)


in which we assume equal-weight particles at time n-1 for ease of presentation.
In this case the MCMC method that has the posterior as invariant density is easy to find
as the transition densities defined above, followed by an accept/reject step. 

When several model time steps are performed between observation times one can also perform 
tempering in the time domain, as explored in <cit.>
and <cit.> in the Guided Particle Filter. The idea is to assimilate the observations ahead of time,
with using as likelihood p(^*|^m)^γ), in which ^* is taken equal to the value ^n, and γ<<1.
Here m<n is the present time of the model. This is then followed by a resampling step.
The procedure can be followed over
several time instances during the forward integration of the particles, increasing γ_i each time.
At the observation time γ=1 is used.
This will force the particles towards the observations and does not need extra jittering because each particle
will see a different model noise realisation  in the model integration after the resampling steps.

Of course one has to compensate for the fact that the transition density has been changed,
and the way to do that is to realise that we have used importance sampling.
Instead of sampling from p(^m|_i^m-1), we sample from a pdf 
q(^m|_i^m-1,^n) ∝ p(^m|_i^m-1)p(^n|^m)^γ,
in which ^* is equal to ^n taken at time m, and with larger observation uncertainty related to γ. 
This means that we have to compensate for the weights
created by this sampling, so we need to introduce particle weights w_i^m = p(_i^m|_i^m-1)/q(_i^m|_i^m-1,^*) ∝ 1/p(^*|_i^m)^γ at each model time step we use this scheme.

The scheme generates extra weights during the model integration, but corrects for them 
at each new time when we resample, ensuring much better positioned particles at the
actual observation time n. It has been used in a reduced-gravity primitive equation model in <cit.>,
but not in high-dimensional settings.




 §.§ Particle flow filters

There is a recent surge in methods that dynamically move the particles in state space
from equal-weight particles representing the prior, p(), to equal-weight particles representing the posterior,
p( |). In other words, one seeks a differential equation 


    d/ds = _s()


in artificial time s≥ 0 with the flow map defining the desired transformation. 
If the initial conditions of the differential equation (<ref>) are chosen from
a pdf p_0(), then the solutions follow a distribution characterized by the Liouville equation


    ∂_s p_s = - ∇_· (p_s _s) .


with initial condition p_0() = p() and final condition p_s_final()=p(|).

Two classes of particle flow filters arise. In the first we start from the tempering approach,
such that s_final=1.
We now take the limit of more and more tempering steps by choosing γ_i=1/n = Δ s
with lim_n →∞, so lim_γ_i → 0, or lim_Δ s →  0, see
<cit.>. This leads to:


    lim_Δ s → 0 p_s+Δ s()     =     p_s() (p(|)/p())^Δ s
        =     p_s() exp[ Δ s (log p(|)-log p() ) ] 
       ≈    p_s() [ 1 - Δ s log p(|) - Δ s log p() ]


Hence we find:


    ∂_s p_s()  = -∇_· (p_s _s) = p_s() (log p(|) - c_s)


with c_s = ∫  p_s() log p(|)d.  Explicit expression for f_s are available for certain pdfs such as 
Gaussians and Gaussian mixtures <cit.>. These particle flow filters can be viewed as a continuous limit of 
the tempering methods described in the previous subsection, avoiding the need for resampling and jittering. 
Note that the elliptic partial differential equation (<ref>) does not determine f_s uniquely. Optimal choices
in the sense of minimizing the L_2(p_s)–norm of f_s lead to the theory of optimal transportation, 
see <cit.> and <cit.>.

Figure <ref> shows the basic idea behind particle flow filters.
    

Alternatively, one can explore ideas from Markov-Chain Monte Carlo (MCMC). 
One MCMC method that generates samples from the posterior is the Langevin Monte-Carlo sampling,
in which a sequence of samples is generated by


    x^j+1 = x^j - Δ s ∇_log p(|) + √(2Δ s)β^j


in which β^j a random forcing term drawn from N(0,).
One can show that in the limit of j →∞ these samples will be samples from the posterior.
The corresponding Fokker-Planck equation for this stochastic PDE reads
(see, for example, <cit.>):


    ∂_s p_s    = ∇_·(p_s ∇_ (-log p(|)) ) +  ∇_·∇_ p_s 
       = -∇_·(p_s {∇_.log p(|) - ∇_log p_s })


This equation corresponds to the deterministic PDE (<ref>)
in which _s() is given by:


    _s() := ∇_log p(|) - ∇_log p_s() = -∇_logp_s()/p( | )




Many other choices are possible that use


    lim_s→∞ p_s = p( | )


in (<ref>). An alternative approach, called Stein variational descent, 
has recently been proposed by <cit.>. Stein variational descent can be viewed as a numerical approximation
to a particle flow (<ref>) with vector field

    _s() := p_s ( ∇_log p(|by) - ∇_log p_s() )

<cit.>. We come back to this method below.

In general, to use any of these methods we need to be able to evaluate p_s(_i), which is typically unknown
as we only know the particle representation of p_s().
One way to solve this issue is to explore kernel embedding.
A numerical implementation of the two formulations (<ref>) and (<ref>) 
can be based on a reproducing-kernel Hilbert space (RKHS) F with reproducing kernel K(.,.), typically taken as a Gaussian. 
In the sequel, we will therefore assume that the kernel is symmetric K(,) = K(,).
The inner product ⟨ g,f⟩_ F in F satisfies the reproducing property

    g() = ⟨ K(,·),g ⟩_ F .

A computational approximation to (<ref>) can now be obtained as follows <cit.>.
One approximates the pdf p_s by


    p_s() = 1/N∑_j=1^N K(_j,) ,


the vector field _s by


    _s() = ∑_j=1^N K(_j,) _s^j/p_s() ,


and the N particles _j move under the differential equations


    d/ds_j= _s^j .


Since the drift term (<ref>) gives 
rise to a gradient flow in the space of pdfs with respect to the 
Kullback–Leibler divergence KL =  KL(p_s||p(·|)) 
between p_s and the posterior pdf <cit.>, it is natural to introduce 
the following particle approximation of the Kullback–Leibler divergence:


    V({_l}) := ⟨ p_s,logp_s/p(· |)⟩_ F .


in the RKHS F and to set


    _s^j := - N ∇__j V({_l})


in (<ref>), which leads to a gradient flow in the particles {_l} minimising V.
Details on the numerical implementation of this approach can be found in <cit.>.

The above formulation restricts the pdf p_s, and hence the prior and the posterior,
to be of the form (<ref>).
Alternatively, one can embed the vector field of the flow in an appropriate reproducing kernel Hilbert space
and not the density itself.
With that we can derive a practical implementation of the Stein variational formulation (<ref>) as follows.
First, note that the change in KL due to the flow field _s can easily be found as:


    dKL     =    lim_ϵ→ 0 KL(p_s+ϵ)- KL(p_s)/ϵ
        =      - ∫ p_s() [ _s()^T∇_log p(|) + ∇_·_s()] d . 
         =    ⟨∇ KL, _s ⟩_ F.


where ∇ KL is the gradient of KL, the maximal functional derivative of KL at every state vector  in the RKHS. Note that F here is different from the Hilbert space used earlier. 
Maximising this change in KL as function of the flow field _s is not trivial in general.
However, with the reproducing kernel property of _s we have


    _s() = ⟨K (·, ) , _s (·) ⟩


in which K is a vector-valued kernel, typically taken as K =  K. Using this in  (<ref>), 
the gradient of the KL divergence is found as


    ∇ KL () = -
    ∫ p_s() [ K(,)∇_log p( | )+ ∇_ K(,)] d .


The important point is that this gradient is independent from _s.
One now chooses _s along this direction, which gives the steepest descent, as

    _s() = - ϵ∇ KL ()

Finally, one replaces the integral in (<ref>) by its empirical approximation, to obtain

    _s(_j) = ϵ1/N∑_l=1^N [K(_l,_j)∇_log p(_l | ) + ∇_ K(_l,_j) ]

 for the dynamics (<ref>) of the N particles _j.
 
The intuition behind Stein variational descent is that the first term in (<ref>) pulls the particles
 towards the mode of the posterior, while the second term acts as a repulsive force that allows
 for particle diversity. <cit.> derived this formulation for a steady-state problem, and <cit.>
 have extended the method to sequential particle filters.
 The scheme is given in Algorithm <ref>.
 
 
 
The free parameter of these methods is the reproducing kernel K(.,.),
which needs to be chosen such that the particles sample the posterior and
that physical (and potentially other) balances are retained. 
One also needs to select a proper time stepping scheme, typically chosen as a forward Euler scheme with variable 
time step ϵ, which can now be viewed as the step length in a gradient descent optimisation algorithm.
 



 §.§ Discussion

Viewing particle filters as a transportation problem from equal-weight particles of the prior 
to equal-weight particles of the posterior has led to an interesting set of filters. 
None of them have been implemented yet in high-dimensional settings, but some of them
are ready to do so. 
The strong involvement of the machine learning community in problems 
of this kind also suggests rapid progress here.
Finally we mention that the equal-weight particle filters from section 2 can be viewed as
one-step transportation filters that explore the proposal density freedom, and in fact
transform equal-weight prior particles at time n-1 to equal-weight posterior particles at observation time n.







§ LOCALISATION IN PARTICLE FILTERS


Localisation is a standard technique in Ensemble Kalman filtering to increase the rank of the
ensemble perturbation matrix, allowing for more observations to be assimilated, and to suppress
spurious correlations where real correlations are very small, but ensemble correlations are 
larger because of sampling noise. Localisation limits the influence of each observation to a localisation area that is much smaller than the full model domain.
This idea can easily be incorporated when calculating the particle weights locally, as pioneered by 
<cit.>, and <cit.>, and used in a high-dimensional parameter estimation problem
in <cit.>. The difficulty, as we shall see, lies in the resampling step:
how does one generate 'smooth' global particles from locally resampled particles.
Smooth is not well defined here, but it is related to the particles having realistic
physical relations (balances) between the model variables. For example, if geostrophic balance 
is dominant, the resampling procedure should not generate particles that are completely
out of geostrophic balance as that would lead to spurious adjustment processes via spurious
gravity waves. Up to now localisation is mainly used in connection with the standard particle
filter, while more advanced proposals, apart from the optimal proposal, have not been explored.
<cit.> provide an excellent review of localisation in particle filtering,
treating a subset of the methods presented here, but including interesting extensions of 
the methods they describe.

The formal way localisation can be introduced in particle filtering is as follows.
Let us denote the state at grid point k as ^k. Hence in contrast to other sections a superscript here denotes 
not the time index, but the grid point. Note that in geoscience applications
each grid point typically has several model variables, so ^k is a vector in general.
Physically it makes sense to assume
that the posterior of the state at this grid point depends only on a subset
of the observations. Let us denote that subset as ^[k]. We can then write:


    p(^k|) ≈ p(^k | ^[k])


In turn, these observations do not depend on the whole state vector but
only on part of it, denoted by ^(k):


    p(^[k]|) = p(^[k]|^(k))

Introduce the notation ^(k)∖ k to denote all those grid points in that part of the state vector
excluding grid point k. Then we can rewrite the above as an integral over the 
joint pdf:


    p(^k | ^[k])= ∫ p(^(k)|^[k])  d^(k)∖ k

Exploring Bayes Theorem we find:


    p(^(k)|^[k])    =   p(^[k]|^(k))/p(^[k])
    p(^(k)) 
       ≈   1/N∑_i^N p(^[k]|_i^(k))/p(^[k])δ(^(k) - _i^(k)) 
        =    ∑_i^N w_i^(k)δ(^(k) - _i^(k))


Taken together, this shows that

    p(^k| ^[k]) ≈∑_i^N w_i^(k)δ(^k - _i^k)

The weights w^k_i thus depend only on the local observations ^[k] and
the local prior particles ^(k)_i, so that the variance
of the weights will be much smaller.
Figure <ref> illustrates how this local weighting could look for two different particles.

   

The approximation (<ref>) is not unrealistic:
a temperature observation in New York is not expected to change our pdf of
the temperature in London at the moment of the observation. There will be,
of course, an effect at later times, but that is not relevant here.
The same assumption underlies the use of localisation in Ensemble Kalman Filters,
and in variational methods when the background error covariance is
constructed. 

However, mathematically it does not follow from the assumption
that under the prior the values of the state at grid points separated by
more than a certain distance are independent. There can be an indirect flow
of information from observations far apart over observations between 
neigboring grid points. In Ensemble Kalman Filters, the Kalman gain is
generally a dense matrix even if ^b ^T + is sparse, because
its inverse (^b ^T +  )^-1 can be dense. On the other hand,
if ^b ^T + is diagonally dominant, then often its inverse is too.

Repeating the localisation procedure for all grid points, we obtain all marginals
of the posterior pdf. However, because the weights w^(k)_i change
from one grid point to the next, it is non-trivial to obtain a consistent
posterior for pairs of state values (^k,^ℓ) (and similarly for
triplets etc.). This can easily be seen using Fig. 5: we would like to retain the
left particle in the central upper half of the domain, and abandon elsewhere.
That would mean that where ever it is abandoned we need to replace it
with another particle, perhaps partly with the particle in the right part of the
figure. At the boundary between particles a discontinuity will exist,
which will lead to unphysical behaviour when this new particle is propagated
forward in time.

This means that to obtain global
particles that can be forwarded with the model equations  
one would need to somehow smoothly glue different particles together. This is a
major problem 
and has hampered localisation in particle filtering since the early 2000's. 
However, recently clever smoothing schemes have been constructed that seem
to work well in high-dimensional geophysical applications. We will report on those
below.

Another issue is that the localisation area cannot be too large to avoid filter collapse. 
As a rule of thumb, when there are more than say 10 independent observations inside 
a local area, the particle filter will still tend to be degenerate for the number of O(10-1000) particles one can typically
afford. This means that when
the observation density is high the localisation areas have to become unphysically small, 
or observations have to be discarded. This issue might be solved using tempering techniques
as discussed earlier, but is often avoided by artificially enforcing a minimal weight of the particles,
or by changing the observations, for instance by projecting them on a lower dimensional
space favoured by the prior.

Setting a minimal weight or projecting observations to a lower dimensional space favoured by the prior
has as consequence that not all information will be extracted from the
observations, as observations that are very different from the existing particles will be largely ignored. 
This is not directly equivalent to the standard quality control measures used by operational
weather forecasting centres, in which observations that are a few standard deviations away from the forecast are ignored.
The issue here is that a distance of less then one standard deviation for a few observations can 
already lead to weight collapse, and artificially setting minimum values for the weights avoids that.






 §.§ Localisation based on resampling


Several localisation schemes have been proposed and discussed in the review by <cit.> and those will not be repeated here.
The most obvious thing to do is to weight and resample locally, and somehow glue the resampled particles together
via averaging at the edges between resampled local particles <cit.>.
In the following, several schemes in this category are discussed.




  §.§.§ The Localized Particle Filter

Recently, <cit.> used this idea with more extensive averaging, and their scheme runs as follows.
First, for each grid point j the observations close to that grid point are found and the weight of each particle i
is calculated based on the likelihood of only those observations:


    w_i,j= p(_j | _i,j)/∑_k=1^N p(_j | _k,j)


in which _j denotes the set of observations within the localisation area. Note the change of notation 
from the previous section, related to the explicit use of the particle index in all the following.
This is followed by resampling via Stochastic Universal Resampling to provide ensemble members ^a_i,j
with i=1,...,N for each grid point j.

<cit.> extended this methodology by updating blocks of grid points locally, and introduce
a smoothing operator in the weights (similar to <cit.>), as:


    w_i,j= ∑_k=1^N_jG(d_j,k/h) (p(_k | _i,k)/∑_m=1^N ∑_k=1^N_jG(d_j,k/h) p(_k | _m,k)


in which G(..) is a distance weighting function, e.g. a Gaussian or an approximation of that, 
d_j,k is the distance between grid points j and k, for each observation _k 
at grid point k in the neighbourhood of grid point j. The parameter h is a distance radius, 
another tuning parameter. This formulation can be used for each grid
point j, but also for each block of grid points j. They note that G can also be a Gaussian 
of a Gaussian, such that it works directly on -log p(_k|_i,k).

As mentioned before, the issue is that two neighbouring grid points can have different sets of particles, and
smoothing is needed to ensure that the posterior ensemble consists of smooth particles.
This smoothing is performed by <cit.> for each grid point j for each particle i by averaging over the Np neighbouring points
within the localisation area around grid point j:


    x^a_i,j= 1/2x^a_i,j + 1/2Np∑_k=1^N_p x^a_i,j_k


in which j_k for k=1,...,N_p denotes the grid point index for those points in the localisation area around grid point j.
The resampling via Stochastic Universal Resampling is done such that the weights are sorted before the resampling, so that
high-weight particles are joined up to reduce spurious gradients.

<cit.> also suggest to smooth this operation, as follows:


    x^a_i,j= α x^a_i,j + (1-α)  ∑_k=1^N_p G(d_j,j_k/h) x^a_i,j_k


with α a tuning parameter. Note that choosing α=1/2 and
G(d_j,j_k/h) = 1/N_p we recover the scheme by <cit.>.

While these schemes have been shown to solve the degeneracy problem in intermediate dimensional systems 
with fixed balances, like the barotropic vorticity model, it is unclear 
how they will perform in complex systems such as
the atmosphere in which fronts can easily be smoothed out, and nonlinear balances broken, see e.g. the discussion
in <cit.>.  




  §.§.§ The Local Particle Filter

A different scheme that involves a very careful process of 
ensuring smooth posterior particles and retaining nonlinear relations has recently been proposed by <cit.>. 
An important difference with the state-space localisation methods discussed above is that
observations are assimilated sequentially to avoid the discontinuity issues of the state-space localisation.
This makes the algorithm non-parallel, so slower than the state-space localisation methods,
but <cit.> demonstrate that a lower root-mean square error (RMSE) can be achieved. 

The scheme proceeds as follows.
First, adapted weights are calculated for the first element y_1 of the observation vector, as


    w̃_i = α p(y_1 | _i)+ 1-α


These weights are then normalised by their sum W̃.
Then the ensemble is resampled according to these normalised weights to form particles _k_i.

The scalar α is an important parameter is this scheme, with α=1 leading to standard weighting, 
and α=0 leading to all weights being equal to 1 (before normalisation).
Its importance lies in the fact that the weights are always larger than 1-α, so even
a value close to 1, say α=0.99, leads to a minimum weight of 0.01 that might seem 
small, but it means that particles that are more
then 1.7 observational standard deviations away from the observations have their weights cut off to a value close to  1-α.
This limits the influence the observation can have on the ensemble. 
Furthermore, the influence of α does depend on the size of the observational error, which
is perhaps not what one would like.
It is included to avoid loosing any particle.

Now the following is done for each grid point j. 
For each member i a weight is calculated as


    ω̃_i = αρ(1,j,r) p(y_1 |_i)  + 1- αρ(1,j,r)


in which ρ(..) is the localisation function with localisation radius r.
These weights are normalised with their sum over the particles, so a normalised weight ω_i for this grid point is obtained.
Note, again, the role played by α.
Then the posterior mean for this observation at this grid point is calculated as


    _j=∑_i=1^N ω_i _i,j


in which _i,j is the state at grid point j of particle i.
Next a number of scalars are calculated that ensure smooth posterior fields <cit.>
as detailed in Algorithm <ref>.

The final estimate becomes:

    ^a_i,j  = _j+r_1j(_k_i,j-_j) + r_2j(_i,j-_j)

where k_i is the index of the i's sampled particle.

This procedure is followed for each grid point so that at the end an updated set of particles is obtained
that have incorporated the first observation. 
As a next step the whole process is repeated for the next observation, with the small change that
ω̃_i is multiplied by ω̃_i from the previous observation, until all observations have been 
assimilated. In this way, the full weight of all observations is accumulated in the algorithm.
Now the importance of α comes to full light: without α the ensemble would collapse
because the ω̃'s would be degenerate when observations are accumulated.

The final estimate shows that each particle at grid point j is the posterior mean at that point
plus a contribution from the deviation of the posterior resampled particle from that mean
and a contribution from the deviation of the prior particle from that mean. 
So each particle is a mixture of posterior and prior particles,
and departures from the prior are suppressed.
When α=1, so for a full particle filter, we find for grid points at the observation location,
for which ρ(1,j,r)=1, that c_j=0, so r_2j=0, and r_1j≈ 1, so
indeed the scheme gives back the full particle filter. 
The basic elements of the scheme are depicted in Algorithm <ref>.



At grid points between observations it can be shown that the particles have the correct first and 
second order moments, but higher-order moments are not conserved.
(<cit.> generate a scheme that is quite similar, but they ensure 
correct first and second moment by exploring the localised covariances between
observed and unobserved grid points directly in a regression step.)
To remedy this a probabilistic correction is applied at each grid point, as follows.
The prior particles are dressed by Gaussians with width 1 and weighted by 
the likelihood weights to generate the correct posterior pdf.
The posterior particles are dressed in the same way, each with weight 1/N.
Then, the cumulative density functions (cdf's) for the two densities are calculated using a trapezoidal rule integration.
A cubic spline is used to find the prior cdf values at each prior particle i, denoted by cdf_i.
Then a cubic spline is fitted to the other cdf, and the posterior particle i is found as the inverse 
of its cdf at value cdf_i. See <cit.> for details. The result of this procedure is that
higher-order moments are brought back into the ensemble between observation points.

This scheme, although rather complicated, is one of the two local particle filter scheme that has been
applied to a high-dimensional geophysical system based on primitive equations in <cit.>.
The other is the Localised Adaptive Particle Filter discussed below.
(<cit.> applied a local particle filter to a high-dimensional quasi-geostrophic system, 
but that system is quite robust to sharp gradients as it does not allow for gravity waves.)






  §.§.§ The Localised Adaptive Particle Filter



The localized adaptive particle filter (LAPF) is based on the localized
version of the ensemble 
transform (<ref>) following the LETKF described in 
<cit.>, see also <cit.>, with localization in observation space,
and resampling in the spirit of Gaussian Mixture filters <cit.>. 
Localization is carried out around each grid point, and a transform matrix  is calculated for each
localization box. We note that, as for the LETKF, the weights given by (<ref>)
depend continuously on the box location and the observations. 

In a first step, the observations are projected into the space spanned by the 
prior particles. As mentioned above, this will reduce the information extracted 
from the observations, but is perhaps less ad-hoc than setting a lower bound
on the weights, as for instance used in the LPF. 
The LAPF carries out local resampling using universal resampling (see e.g. <cit.>).

In a second step, a careful 
adaptive sampling is carried out in ensemble space around each of 
the N temporary particles. 
This scheme runs as follows:

(a) Resampling is carried out based on a 
(radial) basis function centered at each particle. A simple case would be a
Gaussian mixture, where the covariance of each of the centered Gaussians is taken 
as a scaled version c of the local dynamical ensemble covariance . 

(b) The scaling
factor c is individually calculated for each box based on the local observation
minus background error statistics. For details we refer to <cit.>.
By this, the LAPF guarantees to obtain
a spread of the analysis ensemble which is consistent with the local dynamical
observation minus background (o-b) statistics and the observation error covariance . Further standard
tools from the LETKF literature to control ensemble spread can be employed if needed.

(c) To obtain sufficient smoothness of the fields in physical
space, the LAPF uses N global random draws to generate the resampling 
vectors around each particle in the space of ensemble coefficients. 
In combination with the fact that the LAPF draws in each box around each 
particle only – in a globally uniform way modulated by the ensemble 
covariance  and the factor c only –, 
consistency and balance of the fields is achieved with sufficient
precision. The scheme is depicted in Algorithm <ref>.




The LAPF is the first particle filter that has been implemented and tested in an
operational numerical weather prediction context, and we provide a short description of 
the procedure.
The method has been implemented in the data assimilation 
system DACE (Data Assimilation Coding Environment) of Deutscher Wetterdienst 
(DWD) <cit.>. The DACE environment includes a Local 
Ensemble Transform Kalman Filter (LETKF) based on <cit.> both for 
the global ICON model system and the convection permitting COSMO model system 
of DWD, see <cit.>), both of which are run operationally at 
DWD[since January 20, 2016 for the global ICON model with 40km global 
ensemble resolution including a 20km resolved two-way nest over Europe; and 
since March 21, 2017 for the COSMO model with 2.8km resolution over central 
Europe] and build a basis, framework and reference for the LAPF particle 
filter implementation. 

The ensemble data assimilation system is equipped with a variety of tools to control the spread of 
the ensemble, such as multiplicative inflation and additive 
inflation, relaxation to prior spread (RTPS), relaxation to prior 
perturbations (RTPP) and stochastic schemes to add spread to soil moisture 
and sea surface temperature (SST) when needed (details are described in <cit.>).  

Tests with the LAPF for the global ICON model 
with 40 particles of 40km global resolution have been successfully 
and stably run over a duration of one month.
Extensive tests on how many
particles form the basis for resampling in each localization box have been
carried out, the numbers vary strongly over the globe and all heights of the
atmosphere, ranging from 1 to N, with relatively flat distribution. 
Diagnostics and tuning of the system 
is under development and discussed in <cit.>.
Results show that
the quality of the LAPF does not yet reach the scores of the operational 
global LETKF-EnVAR system, but the system runs stably and forecast scores
are about 10-15% behind the current operational system.






 §.§ The Local Ensemble Transform Particle Filter


This filter uses a classic sequential importance resampling particle filter 
from a set of forecast particles _i^ f, which can be obtained employing either the standard
or the optimal proposals (or any other) and their associated importance weights w_i^ f. 
The particles are then resampled in a statistically consistent manner, which can be characterized
by an N× N stochastic transition matrix D with the following properties: (i) all entries
d_ij of D are non-negative and


    ∑_i=1^N d_ij = 1 ,   1/N∑_j=1^N d_ij = w_i^ f .


Let us denote the set of all such matrices by 𝒟. Then any D∈𝒟 leads to a resampling scheme by randomly drawing an element 
j^∗∈{1,…,N} according to the probability vector p_j 
= (p_1j,…,p_Nj) ∈ℝ^N for each j=1,…,N. 
The jth forecast particle _j^ f is then replaced by _j^∗^ f and the new particles 
_j^n = x_j^∗^ f, j=1,…,N, provide an equally weighted set of particles from
the posterior distribution. Note that multinomial resampling corresponds to the simple choice


    d_ij = w_i^ f.


The ensemble transform particle filter (ETPF) <cit.> 
is based on the particular choice D∈ D that minimizes 
the expected squared Euclidian distance between forecast particles, i.e.,


    D = min_ D∈ D∑_i,j=1^N d_ij_i^ f - _j^ f^2 .


It has been shown under appropriate conditions that the variance of a resampling step 
based on D vanishes as N→∞ <cit.>. 
This fact is utilized by the ETPF and one defines

    _j^n = ∑_i=1^N _i^ fd_ij

even for finite particles numbers. Of course, by its very construction, the ETPF underestimates 
the posterior covariance. However, there are corrections available that lead to second-order
accurate implementations <cit.>. See Section <ref> 
for more details.

Following previously introduced notations, localization can now be implemented into the
ETPF as follows. For each grid point k, we extract the values of the forecast particle 
_i^ f  at that grid point and denote them by _i^k. 
Using the observations local to this grid point, we calculate localized importance weights w_i^k for _i^k. 
Then (<ref>) gives rise to a localized transformation matrix 


    D^k = min_ D∈ D^k∑_i,j=1^N d_ij_i^k - _j^k^2


at grid point k with the set 𝒟^k defined by


    𝒟^k ={ D∈ℝ_+^N× N:  ∑_i=1^N d_ij = 1,
    ∑_j=1^N d_ij =w_i^k N }.


Note that the transport cost (distance) t_ij = _i^k - _j^k^2 can be replaced by any
other localized cost function. See <cit.> for more details.  
The transport problem (<ref>) at each grid point can be computationally
expensive.  Less expensive approximations, such as the Sinkhorn approximation, and their
implementation into the localized ETPF (LETPF) are discussed in <cit.>.
<cit.> have extended this algorithm to block weighting, similar to their extension 
of the Local Particle Filter.

The latter authors also defined a local transform particle filter in state space. 
This involves a transformation, at each grid point, from prior to posterior 
particles by a transformation, which essentially becomes an anamorphosis step. 
The prior and posterior probability densities need to be known as continuous densities,
and <cit.> use kernel density estimation with the particles as basis.
The interesting suggestion is that since the transformation is deterministic and 
expected to be smooth over the space coordinates, no specific smoothing is needed
after the transformation.
We refer to their paper for details on this methodology. 






 §.§ Space-Time Particle Filters


The idea to run a particle filter over the spatial domain was introduced by <cit.>,
and the first algorithm, the Location Bootstrap Filter,  was published by <cit.>.
The Space-Time Particle Filter by <cit.>  improves on this algorithm
by removing the jitter step, as explained below.
In the following we assume observations at every grid point, but the algorithms can easily be 
adapted to other observation networks.

The Location Particle Filter of <cit.> runs as follows. 
The grid points are ordered 1,...,L, such that points l and l+1 are neighbouring grid points for each l ∈1,...,L.
In each grid point l we have a sample _i,l for i ∈1,...,N, and l denotes 
the grid point number.
We start the spatial particle filter at location l=1 by calculating the weight p(_1 | _i,1) (where the time index is suppressed) 
for each prior particle i, and perform resampling using these weights over the whole spatial domain.
This means that the resampled particles are now samples of p(^1:L | ^1). 
A small amount of jitter is added to avoid identical particles. The choice of this jitter density is again not clear for geophysical applications, 
more research is needed on this issue.

Then, the algorithm moves to the next grid point, calculates the weights p(_2 | _i,2), and resamples the 
full state particles using this weight, generating samples from p(_1:L | _1,_2).
Again some jitter is needed to avoid ensemble collapse, and the algorithm moves to the next grid point, until
all grid points are treated this way. Algorithm <ref> describes the computational steps.

