# Adaptive tempering schedules with approximative intermediate measures for
filtering problems

Iris Rammelmüller  Department of Statistics, University of Klagenfurt,
Universitätsstraße 65–67, 9020 Klagenfurt, Austria<EMAIL_ADDRESS>Gottfried Hastermann  Institute for Mathematics, Ilmenau University of
Technology, Weimarer Str. 25, 98693 Ilmenau, gottfried.hastermann@tu-
ilmenau.de Jana de Wiljes  Institute for Mathematics, Ilmenau University of
Technology, Weimarer Str. 25, 98693 Ilmenau<EMAIL_ADDRESS>School of Engineering Sciences, Department of Computational Engineering,
Yliopistonkatu 34 53850 Lappeenranta, Finland LUT University, Finland

(Preprint, May 2024)

###### Abstract

Data assimilation algorithms integrate prior information from numerical model
simulations with observed data. Ensemble-based filters, regarded as state-of-
the-art, are widely employed for large-scale estimation tasks in disciplines
such as geoscience and meteorology. Despite their inability to produce the
true posterior distribution for nonlinear systems, their robustness and
capacity for state tracking are noteworthy. In contrast, Particle filters
yield the correct distribution in the ensemble limit but require substantially
larger ensemble sizes than ensemble-based filters to maintain stability in
higher-dimensional spaces. It is essential to transcend traditional Gaussian
assumptions to achieve realistic quantification of uncertainties. One approach
involves the hybridisation of filters, facilitated by tempering, to harness
the complementary strengths of different filters. A new adaptive tempering
method is proposed to tune the underlying schedule, aiming to systematically
surpass the performance previously achieved. Although promising numerical
results for certain filter combinations in toy examples exist in the
literature, the tuning of hyperparameters presents a considerable challenge. A
deeper understanding of these interactions is crucial for practical
applications.

Keywords. Nonlinear filtering, adaptive tempering, approximative consistent
filters
AMS(MOS) subject classifications. 65C05, 62M20, 63G35,

##  1 Introduction

Estimation of large-scale numerical weather and climate prediction signals on
the basis of available observation and evolution models is approached with
data assimilation methods, c.f. [45, 32, 20]. These methods are typically
based on the assumption that the associated posterior distribution is a
Gaussian. This specification leads to highly robust filters, which possess
strong tracking capabilities, c.f. [17, 16]. Many of the systems analysed
exhibit characteristics closely aligned with near Gaussian distributions;
however, they also display pronounced nonlinear behaviour, especially in case
of (very) sparse and uncertain observation. Forecasting of precipitation
serves as one of the most prominent and practical examples. More generally,
nonlinear transport with little dissipative effects is the foundation of
multiple examples with strong non-Gaussian prior distributions. Depending on
the observational model the posterior might then also be substantially non-
Gaussian. Considering that numerous decisions impacting public health rely on
these predictions, it becomes imperative to refine approximations of the
underlying uncertainties in both prior and posterior distributions of each
filtering step. This refinement can be realised through non-parametric
methods, which provide empirical estimates of the distributions in question.
However, vastness of the state space often complicates the direct use of
standard techniques such as sequential Monte Carlo methods. Consequently,
substantial research has been devoted to developing methodologies that not
only yield more accurate uncertainty estimates but also remain computationally
viable, c.f. [4]. A particular line of research has focused on developing and
training cost-effective surrogate models. Various paths are taken in the
derivation of these models, including filtering on multiple levels of coarse-
grained grids [26], estimating parameters for alternative families of
functions like neural networks or systems of linear equations, c.f.[46, 25].
These approaches are often aimed at facilitating the generation of larger
ensembles of samples. To ensure that coarse-grained simulations accurately
reflect the true system, noise carrying the system’s physical information is
incorporated into the evolution equation thereby improving the efficiency and
scalability of simulations while maintaining the same accuracy in generating
samples of the posterior. While this ansatz is particularly promising for
specific large-scale systems, e.g. [13, 14], it is not trivial to adapt to
other model equations and requires a series of other tools namely jittering
(or rejuvenation/ inflation techniques) and tempering, which in general
improve the feasibility of non-parametric empirical inference approaches. One
of the key advantages of most Gaussian approximate filters opposed to standard
particle filters is that individual ensemble methods are moved towards the
observation in the equations, which can be interpreted as evolution equations
of the particles in time in a continuous setting. It is feasible to develop
analogous evolution equations for the non-parametric scenario, leading to the
introduction of filters such as the Feedback Particle Filter [49], see [42]
for an extensive and thorough generalisation. While these approaches have a
similar update structure as their Gaussian counterparts they are typically
limited by the existing computational boundaries. As the aim is to map samples
from the prior, associated with the underlying model, to a posterior
distribution conditioned on the important additional observation information
source a family of methods is focused on estimating the underlying transport
map [47, 8, 41]. These maps are calculated based on the available particles
and associated importance weights. However, in high dimensions, the samples
may not sufficiently cover the underlying state space. Therefore, it is often
necessary to either use one of the aforementioned techniques or to combine
complementing filters. The latter integration is facilitated through a
procedure known as tempering, which involves dividing the likelihood into
parts and then merging non-parametric with parametric filters, c.f. [10, 22].
This blend of filters is often termed hybrid, which should not be confused
with the usage of ’hybrid filters’ in the context of combining ensemble and
variational methods such as 3DVar or 4DVar [23]. One can view these approaches
in the context of preconditioners or, improved proposal densities. They have
been successfully numerically implemented for specific combinations of filters
for a few toy models [10, 1, 40]. Although these numerical experiments
highlight the approach’s potential, it is important to extensively examine the
impacts of the ordering of the filters, fine-tuning hyperparameters such as
the tempering schedule, and testing the underlying hypotheses across a broader
spectrum of filters to enable application in real-life data assimilation set-
ups. To explore the effects of various tempered combinations more effectively,
we designed a series of toy setups to numerically investigate underlying
systems. Specifically, we employ a Langevin model to simulate problems
characteristic of high-dimensional spaces. We also propose a novel method for
adaptively tuning the tempering schedule, which shows promising results across
all toy models and filter combinations compared to previous approaches using
the effective sample size.

##  2 Mathematical framework

Let $z_{k}\in{\mathbb{R}}^{N_{z}}$ be the state of a dynamical system and
$\Psi:{\mathbb{R}}^{N_{z}}\to{\mathbb{R}}^{N_{z}}$ be the solution operator of
a dynamical system. Given an initial condition $z_{0}\in{\mathbb{R}}^{N_{z}}$,
a dynamical system defines a solution sequence $\\{z_{k}\\}_{k\geq 0}$ via the
forward iteration

$z_{k+1}=\Psi(z_{k}).$ (1)

Let $\\{y_{k}\\}_{k\geq 0}=y\in{\mathbb{R}}^{N_{y}}$ be the noisy observations
then at each discrete time instance we observe a possibly nonlinear function
of the signal

$y_{k+1}=H(z_{k+1})+\eta_{k+1},$ (2)

where $H:{\mathbb{R}}^{N_{z}}\to{\mathbb{R}}^{N_{y}}$ is the observation
operator and $\eta_{k+1}$ is Gaussian distributed with mean zero and error
covariance matrix $R\in\mathbb{R}^{N_{y}\times N_{y}}$. Given this discrete
scenario, we would like to solve the time-discrete filtering problem, i.e. to
estimate the density $p(z_{k}|y_{1:k})$ associated with the signal $z_{k}$
conditioned on the observation up to time $k$. Employing Bayes’ theorem allows
us to leverage the following relationship:

$p(z_{k}|y_{1:k})\propto p(y_{k}|z_{k})p(z_{k}|y_{1:k-1}),$ (3)

where $p(y_{k}|z_{k})$ is the likelihood and $p(z_{k}|y_{1:k-1})$ is the prior
distribution c.f. [28]. This representation of the posterior allows one to
modify the likelihood. For instance one can introduce artificial measures

$\frac{d\mu_{r}}{d\mu_{0}}\propto p(y_{k}|z_{k})^{\tau_{r}}$ (4)

where $0=\tau_{0}<\ldots<\tau_{r}<\ldots<\tau_{T}=1$ is the schedule, which
represents the artificial time steps in between the two actual time $k-1$ and
$k$, $\mu_{T}=p(z_{k}|y_{1:k})$ and $\mu_{0}=p(z_{k}|y_{1:k-1})$. Then one can
use this split to sequentially update with the intermediate measures

$\frac{d\mu_{r}}{d\mu_{r-1}}\propto p(y_{k}|z_{k})^{\tau_{r}-\tau_{r-1}}.$ (5)

This procedure, known as tempering, is particularly beneficial when the
supports of the prior and the posterior distributions do not overlap. Although
this technique is employed broadly within the context of sequential Monte
Carlo methods, c.f. [9], it is especially crucial for filtering problems [14,
4]. Here, we consider a specialised version where different approximation
schemes are used for the intermediate updates. This approach shows very
promising results, though the benefits vary depending on the quality of the
filters used, the underlying system dynamics, and the available data [10, 22].
Therefore, for these approximate intermediate measure updates, adaptive
tempering with an appropriate combination of filters is desirable and is the
main focus of this work. One method for achieving this is introduced in
Section 5, and extensively numerically explored in Section 6. In the following
section, we will provide brief introductions to representatives of important
filtering classes such as Gaussian-approximative filters, transport-based
particle filter and McKean-Vlasov filters that we will amongst others consider
for our approximate intermediate updates.

##  3 Gaussian-approximative filters

The classical Kalman filter provides the optimum filtering solution for
Gaussian noise and linear maps $\Psi$ and $H$, e.g. [30, 28]. In general, if
either $\Psi$ or $H$ are nonlinear or the initial distribution is not
Gaussian, it is not necessarily possible to identify the (parametric) family
of distributions associated with the posterior. Empirical estimates,
linearisation, or a combination of both are usually used to enforce a Gaussian
approximation of the associated densities. A summaries of the different
variants of Kalman based filters, e.g. Ensemble Kalman Filter (EnKF), Ensemble
Square Root Filter, are given in [45, 32, 20]. The EnKF [19, 31] and its
variants use the Kalman filter update for the mean and covariance to the their
empirical approximation provided by the ensemble. In real-world applications,
a stable implementation with small numerical error propagation is desirable.
The stochastic noise of the vanilla variante of the EnKF [7] can cause
numerical challenges and is only producing correct moments in the linear case
in the ensemble limit. As an alternative the family of Ensemble Square Root
Filter was developed. These filters use rank reduction methods, and a square
root matrix to overcome the matrix oversize problem, c.f. [51] and yield
correct moments in the Gaussian case. In this article, we consider a
particular variant within the family ensemble square root filters as described
in [45], yet other members can be used interchangeably.

###  3.1 Ensemble Square Root Filter

The Ensemble Square Root Filter (ESRF) as presented in [45] utilises Monte
Carlo estimates to compute the forecast mean and anomalies as follows:

$\displaystyle\overline{z}^{f}$
$\displaystyle=\frac{1}{N_{ens}}\sum_{i=1}^{N_{ens}}z^{f}_{i},$ (6)
$\displaystyle A^{f}$ $\displaystyle=(z^{f}_{i}-\overline{z}^{f}).$ (7)

The rank reduction approach, which involves computing only the anomalies of
the covariance matrix instead of the full matrix, is used to prevent oversized
vectors and matrices in high-dimensional systems. The main idea of the ESRF is
to define the analysis step so that the empirical covariance of the resulting
particle ensemble exactly matches the Kalman mean and covariance in the linear
case:

$\displaystyle P^{f}-KHP^{f}=P^{a}$
$\displaystyle=\frac{1}{N_{ens}-1}A^{f}SS^{T}(A^{f})^{T}.$ (8)

Rearranging (8) leads to the following representation of the square root
matrix $S$.

$S=\left(I_{N_{ens}}+\frac{1}{N_{ens}-1}(HA^{f})^{T}R^{-1}(HA^{f})\right)^{-1/2}$
(9)

ESRF by construction avoids loss of positive definiteness of the error
covariance matrices. Note that not all matrix square roots bear the desired
relationship to the forecast ensemble: the analysis ensemble mean may not be
equal to the analysis state estimate and consequently there may be an
accompanying shortfall in the spread of the analysis ensemble as expressed by
the ensemble covariance matrix. This points to the need for a restricted
version of the notion of an ESRF, which is called an unbiased ESRF, c.f. [36].
There is a set of generic necessary and sufficient conditions for the ensemble
to remain centred on the analysis state estimate (unbiased). Hence, the
entries $t_{ij}$ of the symmetric matrix $S$ should satisfies the following
conditions:

$\sum_{i=1}^{N_{ens}}s_{ij}=\sum_{j=1}^{N_{ens}}s_{ij}=1.$ (10)

Next one calculates the weights, which can be interpreted as importance
weights, whereby $e_{i}$ denotes the i${}^{\text{th}}$ basis vector in
${\mathbb{R}}^{N_{ens}}$:

$\displaystyle\hat{w}_{i}$
$\displaystyle=\frac{1}{N_{ens}}-\frac{1}{N_{ens}-1}e_{i}^{T}S^{2}(HA^{f})^{T}R^{-1}(H\overline{z}^{f}-y^{obs}),$
(11)

with filter coefficients

$\displaystyle d_{ij}$ $\displaystyle=\hat{w}_{i}-\frac{1}{N_{ens}}+s_{ij}.$
(12)

The weights $\hat{w}_{i}$ and the coefficient $d_{ij}$ satisfy

$\sum_{i=1}^{N_{ens}}\hat{w}_{i}=1\ \text{ and }\
\sum_{i=1}^{N_{ens}}d_{ij}=1.$ (13)

Finally, we can update the ensemble members with above equations using a
general linear transformation from forecast to analysis state as described in
[45], which is given by

$z_{j}^{a}=\sum_{j=1}^{N_{ens}}z^{f}_{i}d_{ij}.$ (14)

The implementation scheme can be found in the Appendix Algorithm 4. According
to [10] the computational complexity scales like $\mathcal{O}(N_{ens}^{3})$ in
the ensemble size.

##  4 Approximative consistent filters

In the context of nonlinear evolution models, a non-parametric empirical
approach, such as a particle filter (PF), is considered state-of-the-art for
state and parameter estimation tasks, c.f. [21], if sufficiently many
particles can be generated. A key characteristic of these empirical filters is
their consistency: they converge towards the true posterior density in
nonlinear and non-Gaussian scenarios without requiring significant assumptions
or restrictions on the prior or observations. Generally, particles of the
prior are updated using individual importance weights derived from the
likelihood. In high-dimensional spaces, particles are less likely to reside in
regions of significant probability mass relative to the likelihood.
Consequently, after performing importance sampling with respect to the
likelihood, the posterior weights exhibit significant fluctuations. Many
particles end up with very low relative weights, leading to a phenomenon known
as filter collapse after a few iterations. This occurs because these low-
weight particles are effectively lost for subsequent forecasting steps,
demonstrating the filter’s vulnerability to the curse of dimensionality. We
characterise filter collapse by the effective sample size falling below a
given threshold, as discussed in [22]. The desire to mitigate these effects
has inspired the development of various families of consistent filters. While
none of the proposed filters are immune to the curse of dimensionality, they
still do not assume any specific parametric structure for the underlying
densities. However, it is important to mention that some form of approximation
may still be involved (details in the respective sections).

###  4.1 Bootstrap Particle Filter

The vanilla Bootstrap Particle Filter belongs to the family of Sequential
Monte Carlo methods and is a simulation-based methods for calculating
approximations to posterior distributions, c.f. [35, 9]. As mentioned above
this is achieved via importance weights, which are defined by:

$\displaystyle\hat{w}_{i}$
$\displaystyle\propto\exp\left(-\frac{1}{2}(Hz^{f}_{i}-y_{k}^{obs})^{\top}R^{-1}(Hz^{f}_{i}-y_{k}^{obs})\right)\cdot
w_{0},\ \ \ w_{0}=\frac{1}{N_{ens}},$ (15) $\displaystyle w$
$\displaystyle=\frac{\hat{w}_{i}}{\sum_{i=1}^{N_{ens}}\hat{w}_{i}}.$ (16)

To address the issue of weight degeneracy, a resampling procedure is employed.
This method of sequential importance resampling selectively removes particles
with small weights and duplicates particles with large weights, effectively
mitigating the degeneracy problem. Consequently, a new set of ensembles is
resampled based on the calculated weights, resulting in an updated ensemble.
However, this procedure can produce identical or near identical analysis
ensemble members in practice. As a result, the filter can collapse. To
overcome this issue, tuning tools are required c.f. [2, 13]. We present some
of the variants thereof in Section 6.1. A pseudo-code implementation of the
method is provided in the Appendix Algorithm 5.

###  4.2 Transport based Particle Filter

Transport-based filters approximate a coupling between the prior and posterior
distributions. Typically, the goal is to find couplings that optimise an
underlying optimal transport problem with respect to a given distance metric,
such as the Monge-Kantorovich problem with the Wasserstein distance. It can be
demonstrated that these optimal couplings are induced by a transport map.
Several works [18, 49, 47, 8, 41] suggest approximative discrete optimal
transport maps as surrogates of the true map. The advantage of using a map
rather than a resampling step is that particles are moved rather than
replicated, making them less dependent on the existing particles. However, it
is important to note that the map is approximated based on the existing
particles, which can potentially negatively affect the quality of the
estimated map.

####  4.2.1 Ensemble Transform Particle Filter

The Ensemble Transform Particle Filter (ETPF) relies on a linear programming
construction to approximate the optimal transport map, c.f. [44]. As indicated
above, the resulting linear transformation replaces the resampling step
typically used in a standard Sequential Importance Resampling (SIR)
filter[45]. The map is approximated as follows

$T^{\ast}=\mathop{\mathrm{argmin}}\limits_{\begin{subarray}{c}T\in\mathcal{T}\end{subarray}}\sum_{i,j=1}^{N_{ens}}t_{ij}d(z^{f}_{i},z^{f}_{j})^{2},$
(17)

where $\mathcal{T}=\\{T\in\mathbb{R}^{N_{ens}\times
N_{ens}}\colon\sum_{i=1}^{N_{ens}}t_{ij}=\frac{1}{N_{ens}}\land\sum_{j=1}^{N_{ens}}t_{ij}=\hat{w}_{i}\\}$
and $d$ is a suiteable metric on the state space i.e. on $\mathbb{R}^{N_{z}}$.
Note that the weights $\hat{w}_{i}$ are given like (15). The implementation
scheme can be found in the Appendix Algorithm 6. The ETPF is computationally
expensive since one needs to solve the linear transport problem (17) in every
assimilation step, c.f. [45]. Using the same algorithm to solve the system as
in [10], i.e. FastEMD, the computational complexity scales like
$\mathcal{O}(N_{ens}^{3}\log(N_{ens}))$ in the ensemble size. As an
alternative to the exact solution of the linear transport problem, one can
compute the corresponding Sinkhorn approximation, which is obtained by
regularising the underlying optimisation problem, c.f. [15]. As the spread is
severely underestimated, one can also compute a second-order correction
extension of the ETPF [1].

###  4.3 McKean-Vlasov Filters

The general idea of moving particles instead of resampling them is also
applied here. In particular, a modified evolution equation, i.e. in the
continuous setting, is developed, incorporating the innovation term, i.e.
describing the discrepancy between observation and projected state, weighted
by a generalised gain function. This evolution equation directly generates
samples from the posterior. To achieve this, the Fokker-Planck equation of the
modified evolution equation is equated with the Kushner-Stratonovich equation,
which describes the evolution of the posterior density in the continuous
filtering problem. As a result, the generalised gain must be a solution of a
weighted Poisson equation. Using this general approach, different families of
filters can be derived. For a derivation of the general class of McKean-Vlasov
Filters, see [42]. Our focus is on a variant of these filters, known as the
Feedback Particle Filter (FPF) [52, 48].

####  4.3.1 Feedback Particle Filter

As previously stated, McKean-Vlasov Filters seek to derive a generalised gain
that is utilised to update the prior sample with a structure similar to
Ensemble Kalman filters, resulting in samples from the posterior [52]:

$z_{i}^{a}=z_{i}^{f}-K_{i}^{(\epsilon,N_{ens})}\left(\frac{1}{2}(H(z_{i}^{f}))+H(\overline{z}^{f}))-y_{obs}\right).$
(18)

To compute the gain $K$, a weighted Poisson equation

$\displaystyle\begin{split}\nabla\cdot(p(z,k)\nabla\phi(z,k))&=-R^{-1}(H(z)-\hat{H})p(z,k),\\\
\int\phi(z,k)p(z,k)dz&=0,\end{split}$ (19)

with $p$ being the conditional distribution of the state given the observation
and $k\geq 0$, must be solved at every time step. Further, $R$ is the
covariance matrix from the data model (2). The gain function is given for
every time step $k$ by

$\displaystyle K(z,k)$ $\displaystyle=\nabla\phi(z,k),$ (20)

whereby $\phi$ is the solution of the Poisson equation. [53] pointed out that
$K$, in principle is, not uniquely determined. Yet [49] presented results on
existence, uniqueness and regularity of (19) based on a weak formulation of
the Poisson equation. Hence, the choice of the gain function in the
multivariable case requires careful consideration of the uniqueness of the
solutions of (19). A significant challenge in this context is that (19)
contains the posterior density, which is only available in the form of a
sample approximation. Unfortunately, the estimation of the gain is not
straightforward, since there is no closed-form solution of (19), and the
density $p$ is not explicitly known. Due to the limited amount of solution
strategies available for high dimensional elliptic problems, we focus on a
diffusion map [11, 12] approximation, theoretically investigated in [43]. The
underlying idea is to use diffusion map coordinates, which are obtained by
uncovering the underlying geometry spanned by the particles. This approach
transforms the problem into a lower-dimensional space, making it easier to
solve. By leveraging the geometric structure of the data, the diffusion map
approximation provides an efficient way to handle the complexities of the
weighted Poisson equation within the McKean-Vlasov Filter framework. The first
step is to construct a Markov matrix $T^{(\epsilon)}$ for $\epsilon>0$ defined
on the state space approximating the semi-group, c.f. [24], generated by
$\frac{1}{p}\nabla\cdot(p\nabla)$. This semi-group then provides the solution
of the corresponding heat equation, which is given by

$\frac{d}{dt}u=\frac{1}{p}\nabla\cdot(p\nabla)u.$ (21)

Next we consider solutions to the initial value problem posed by (21) and
$u(0)=\phi$ and denote them by $\phi^{\epsilon}=T^{(\epsilon)}\phi$. This, in
turn, is used to formulate the following fixed-point problem derived
(formally) from the fundamental theorem of calculus and (19)

$\phi^{(\epsilon)}=T^{(\epsilon)}\phi^{(\epsilon)}+R^{-1}\int_{0}^{(\epsilon)}T^{(s)}(H-\hat{H})ds,$
(22)

with

$\displaystyle T^{(\epsilon)}\phi(z_{i})$ $\displaystyle:=\frac{\int
k^{(\epsilon)}(z_{i},z_{j})\phi(z_{j})d\mu(z_{j})}{\int
k^{(\epsilon)}(z_{i},z_{j})d\mu(z_{j})},$ (23) $\displaystyle
k^{(\epsilon)}(z_{i},z_{j})$
$\displaystyle:=\frac{g^{(\epsilon)}(z_{i},z_{j})}{\sqrt{\int
g^{(\epsilon)}(z_{i},z_{j})d\mu(z_{j})}\sqrt{\int
g^{(\epsilon)}(z_{i},z_{j})d\mu(z_{i})}},$ (24) $\displaystyle
g^{(\epsilon)}(z_{i},z_{j})$
$\displaystyle=\exp\left(-\frac{|z_{i}-z_{j}|^{2}}{4\epsilon}\right),$ (25)

whereby $z_{i}$ and $z_{j}$ are ensembles of the state vector, $\mu$ is an
absolute continuous probability measure with respect to the posterior $p$, and
$g$ is the Gaussian kernel.

For a numerical implementation one uses
$\epsilon(H-\hat{H})\sim\int_{0}^{\epsilon}T^{(s)}(H-\hat{H})ds$. The particle
approximation is obtained by Monte Carlo integration and empirical
approximation of $T^{(\epsilon)}\phi$. To the knowledge of the authors, there
is a gap in the literature when implementing the FPF in time-discrete
scenarios. For the time continuous case, $z_{i}$ are already samples
consistent with the filtering distribution, i.e. the posterior. However, in
the time discrete case, only the prior samples $z_{i}^{f}$ are available. In
general, they are not consistent with $\mu$, for every $\Delta t_{obs}>0$.
Therefore, we propose to use an empirical approximation, to our understanding
different from [52, 53, 48, 50]. To this end, we apply Monte Carlo
integration, importance sampling, and (3)

$\displaystyle T^{(\epsilon)}\phi(z_{i})$ $\displaystyle=\frac{\int
k^{(\epsilon)}(z_{i},z_{j})\phi(z_{j})d\mu(z_{j})}{\int
k^{(\epsilon)}(z_{i},z_{j})d\mu(z_{j})}$ (26)
$\displaystyle=\frac{\sum_{j}k^{(\epsilon)}(z_{i},z_{j})\Phi^{(\epsilon,N_{ens})}_{j}e^{-\frac{1}{2}(H(z_{j})-y_{obs})^{T}R^{-1}(H(z_{j})-y_{obs})}}{\sum_{j}k^{(\epsilon)}(z_{i},z_{j})e^{-\frac{1}{2}(H(z_{j})-y_{obs})^{T}R^{-1}(H(z_{j})-y_{obs})}}$

and subsequently state the following empirical approximation of the Markov
transition matrix

$T^{(\epsilon,N_{ens})}_{i,j}=\frac{k^{(\epsilon)}(z_{i},z_{j})e^{-\frac{1}{2}(H(z_{j})-y_{obs})^{T}R^{-1}(H(z_{j})-y_{obs})}}{\sum_{j}k^{(\epsilon)}(z_{i},z_{j})e^{-\frac{1}{2}(H(z_{j})-y_{obs})^{T}R^{-1}(H(z_{j})-y_{obs})}}.$
(27)

Note that $\Phi^{(\epsilon,N_{ens})}_{j}=\phi(z_{j})$ for $j=1,...,N_{ens}$. A
pseudo implementation scheme is given in Appendix Algorithm 7. It is important
to state, that despite the consistency theory, provided in [50], the numerical
resolution of (22) is unstable for fixed particle number $N_{ens}$ and
$\epsilon\rightarrow 0$. This is not a contradiction to other well-posedness
statements in [50], as the corresponding limits in $\epsilon$ and $N_{ens}$ do
not commute. When considering (22) and (25), it comes at no surprise that
depending on the dimension of the state space, we must scale the bandwith of
the kernel accordingly, otherwise $T^{(\epsilon)}$ coincides with the identity
in every reasonable finite floating point precision and (22) becomes ill-
posed. For the actual implementation in our time discrete scenario, we adapt
this estimation To avoid another tuning parameter, we scale the bandwidth/time
parameter with the maximum distance between particles. This obviously leads to
higher approximation errors of the integral term, but considering the arguably
low number of particles we do not expect this to be a relevant source of
approximation errors. However, one could improve the latter by using better
quadrature formulas to approximate the integral.

##  5 Adaptive tempering

Here, we propose an adaptive approximative tempering method, which is a
variant of the traditional tempering introduced in Section 2. Updates of the
intermediate measures (5) are performed using members of different families of
filters, including approximative Gaussian filters. This approach introduces an
error, as the true posterior cannot be fully preserved within the tempering
procedure. However, the advantage of this approximative ansatz lies in its
enhanced robustness in scenarios where individual representatives of
approximative consistent filter families encounter difficulties. For our
tempering schedule, we restrict ourselves to $T=2$, which corresponds to
$\tau_{0}=0$, $\tau_{1}=\alpha$, and $\tau_{2}=1$. While this method can be
extended to include more tempering steps (either through traditional
consistent tempering or approximative updates), the increased number of
possible combinations significantly raises the tuning overhead, which is
already non-trivial in the simplified case considered here. Using the notation
in (5) for our simplified schedule yields the following tempered likelihood:

$\displaystyle\begin{split}p({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}y_{k}}|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}z_{k}})&\propto\exp\left(\frac{\alpha}{2}(Hz^{f}_{i,k}-y_{k}^{obs})^{\top}R^{-1}(Hz^{f}_{i,k}-y_{k}^{obs})\right)\\\
&\times\exp\left(\frac{1-\alpha}{2}(Hz^{f}_{i,k}-y_{k}^{obs})^{\top}R^{-1}(Hz^{f}_{i,k}-y_{k}^{obs})\right).\end{split}$
(28)

This or similar splitting approaches have been explored in the literature for
specific filters [10, 22]. This strategy is often referred to as a hybrid
approach, aiming to combine filters with complementary properties. Typically,
this involves selecting one filter from the family of Gaussian approximative
filters and one from the family of (approximative) consistent filters. We
introduce key representatives from the important classes of consistent filters
and focus on the family of Ensemble Square Root Filters within the Gaussian
approximative filter family. While it is possible to consider different
filters from these two classes, the key lies in the ability to adjust the
likelihood in the filter under consideration. Since the likelihood is
incorporated differently in the filter family examined in this paper, we will
provide a brief overview of how the split affects the updates.

###  5.1 Tempered ESRF

For the ESRF, the square root matrix (9) and the weights (11) change into:

$\displaystyle S(\alpha)$
$\displaystyle=\left(I_{N_{ens}}+\frac{1-\alpha}{N_{ens}-1}(HA^{f})^{T}R^{-1}(HA^{f})\right)^{-1/2}$
(29) $\displaystyle\hat{w}_{i}(\alpha)$
$\displaystyle=\frac{1}{N_{ens}}-\frac{1-\alpha}{N_{ens}-1}e_{i}^{T}S^{2}(HA^{f})^{T}R^{-1}(H\overline{z}^{f}-y^{obs}).$
(30)

###  5.2 Tempered Bootstrap and ETPF

For the Bootstrap and ETPF the weights in (15) change into:

$\displaystyle\hat{w}_{i,k}(\alpha)$
$\displaystyle\propto\exp\left(-\frac{\alpha}{2}(Hz^{f}_{i,k}-y_{k}^{obs})^{\top}R^{-1}(Hz^{f}_{i,k}-y_{k}^{obs})\right)\cdot
w_{0}.$ (31)

Additionally, the transport map (17) from the ETPF changes as follows:

$T^{\ast}(\alpha)=\mathop{\mathrm{argmin}}\limits_{\begin{subarray}{c}T\in\mathcal{T}\end{subarray}}\sum_{i,j=1}^{N_{ens}}t_{ij}d(z^{f}_{i},z^{f}_{j})^{2},$
(32)

with the constraints

$\sum_{i=1}^{N_{ens}}t_{ij}=\frac{1}{N_{ens}},\ \ \
\sum_{j=1}^{N_{ens}}t_{ij}=\hat{w}_{i}(\alpha).$ (33)

###  5.3 Tempered FPF

For the FPF we insert the bridging parameter in (27) as follows:

$T^{(\epsilon,N_{ens})}_{i,j}=\frac{k^{(\epsilon)}(z_{i},z_{j})e^{-\frac{\alpha}{2}(H(z_{j})-y_{obs})^{T}R^{-1}(H(z_{j})-y_{obs})}}{\sum_{j}k^{(\epsilon)}(z_{i},z_{j})e^{-\frac{\alpha}{2}(H(z_{j})-y_{obs})^{T}R^{-1}(H(z_{j})-y_{obs})}}.$
(34)

###  5.4 Tuning algorithm

At this point, it remains unclear in which filtering cycles tempering is
beneficial and how to choose the bridging parameter $\alpha$. One brute-force
approach to selecting the value for $\alpha$ is to evaluate the bridging
parameter in every assimilation cycle, as initially discussed in [22].
Following this approach, [40] provides a detailed study on the choice of the
bridging parameter for the tempered combination of the local ensemble
transform Kalman filter with the nonlinear ensemble transform filter.
Traditionally, the effective sample size is considered state-of-the-art for
adjusting the tempering schedule. While these techniques can easily be
combined with our adaptive tempering approach, we limit ourselves to fixed
values of $\alpha\in\\{0.2,1\\}$. This choice simplifies the evaluation of the
different combinations already being considered and aligns with the intuition
that small updates with beneficial properties are sufficient to adjust the
main filter being used. Additionally, this choice is consistent with the
values found to be optimal by numerical experiments for a complementary filter
combination considered in [10]. In previous works on the type of tempered
filters considered in here, tempering was fixed to occur in every data
assimilation cycle. Our main focus is to adaptively decide in each filtering
step whether a tempering step is required. Intuitively, the effective sample
size offers a measure of the number of particles actively contributing to the
update procedure. This approach is not only intuitively justified but also
provides valuable statistical insights for tempering and the decision process
of whether to resample. We aim to establish a even more incisive and
informative criterion for determining whether an approximative tempering step
is necessary. The novel criterion we propose relies on general descriptive
statistics of the particles mapped to the observational space. To benchmark
against the state-of-the-art approach, we compare our numerical results with
the fixed tempering method employed in earlier works on hybrid filtering (for
details see Section 6).

####  5.4.1 Tempering criterion I

In a general particle filter setting the effective sample size (ESS) is used
as an effective criterion to resample and is given as

$ESS=\frac{1}{\sum_{i=1}^{N_{ens}}w_{i}^{2}}.$ (35)

This approach simplifies the comparison of filters that require an update of
the particles at each step. Rather than using the ESS for resampling, we use
it to decide if tempering, which entails two complementary filtering updates,
should be performed. Specifically, the default is to perform updates only with
one filter, unless the ESS drops below a certain threshold $\theta$.

$ESS<\theta N_{ens}.$ (36)

The threshold is set at $\theta=0.5$ times the ensemble size in this case,
which is a commonly suggested value in the resampling context. Naturally, the
hyperparameter $\theta$ can also be adjusted to more extreme values and this
is certainly a topic of interest for future research.

####  5.4.2 Tempering criterion II

Since the decision-making process should be based on the discrepancy between
observed particle positions and the observational data, we propose using a
criterion grounded in descriptive statistics for this property. Specifically,
we apply an individual filter of our choice if the current observation falls
within the interquartile range (IQR) of the observed particles. If the
observation is outside the IQR, we employ a tempered combination approach: one
update is performed with the currently considered filter, and the other update
is performed with a filter from a different family. The main aim is to use
consistent filters to enhance the uncertainty quantification of the estimate.
These filters typically require periodic adjustment through robust alternative
filters. However, we also consider Gaussian approximative filters as primary
filters combined with various other filters to investigate the associated
benefits. An illustration of a box plot is given in Figure 1, whereby the
observation is outside the IQR, which illustrates the case when tempering is
applied. Box plots and the underlying descriptive statistics offer an
effective graphical representation of data concentration, as illustrated in
Figure 1. They also accurately depict the extent to which extreme values
deviate from the central tendency of the data. A box plot represents five key
statistics: the minimum (lower whisker), the first quartile (Q1), the median
(Q2), the third quartile (Q3), and the maximum (upper whisker). Further,
outliers are marked as points as depicted in Figure 1. The position of the
quantiles are determined via:

$Q_{i}=\frac{i\cdot(N_{ens}+1)}{4}\ \text{ for }i=1,2,3.$ (37)

An outlier is a data point that differs numerically from the rest of the data.
They are determined by a quantile factor $1.5$ times the interquartile range
(IQR), which is given as

$IQR(\\{\xi_{i}\\}_{i\in\\{1,\dots,N_{ens}\\}})=[\xi_{\sigma(Q_{3})},\xi_{\sigma(Q_{1})}].$
(38)

where $\xi_{i}\in\mathbb{R}$ is an ensemble of real numbers and $\sigma$ is
the permutation such that $\xi_{\sigma(i)}$ is monotonously ordered for
$i\in\\{1,\dots,N_{ens}\\}$. The IQR interval relates to fifty per cent of the
score. The usual value of the quartile factor, associated with detecting minor
outliers of the sample set, is $1.5$. Hence, we differ the length of the IQR
with the quantile factor $\\{0,1.5\\}$. Note that the lower and upper whisker
excludes the outliers. Twenty-five per cent of scores fall below the $Q1$.
Contrary, Twenty-five percent of data are above $Q3$. Meaning that seventy-
five per cent of the scores fall below the third quartile. Figure 1 depicts a
line that divides the box into two portions to indicate the median, which
represents the midpoint of the data. Half of the scores exceed or equal this
value, while the other half fall below it. The box plot shape indicates
whether a statistical dataset is normally distributed or skewed. When the
median is in the centre of the box and the whiskers are roughly equal on both
sides, the distribution is symmetric. When the median approaches $Q1$ and the
whisker on the lower end of the box is shorter, the distribution is positively
skewed, i.e. skewed right. When the median is closer to $Q3$ and the whisker
is shorter on the top end of the box, the distribution is negatively skewed,
i.e. left skewed. Furthermore, the box plot is useful as it shows the
dispersion of a data set. The data is more evenly distributed as the IQR
length increases. The smaller the IQR, the lower the spread of the data. One
can determine how similar other data values are. To apply the criterion to an
ensemble of particles in $z_{i}\in\mathbb{R}^{N_{z}}$ we first observe the
states and subsequently project onto each individual component of the
observation i.e. we compute

$IQR(\\{e_{j}^{T}Hz_{i}\\}_{i\in\\{1,\dots,N_{ens}\\}})$ (39)

for each component $j\in\\{1,\dots,N_{y}\\}$ and with $e_{j}$ being the
canonical basis vector in $\mathbb{R}^{N_{y}}$. Subsequently, the criterion
then reads as in Algorithm 1.

Algorithm 1 IQR criterion for ensembles in $\mathbb{R}^{N_{ens}\times N_{y}}$

if $e_{j}^{T}y_{obs}\in
IQR(\\{e_{j}^{T}Hz_{i}\\}_{i\in\\{1,\dots,N_{ens}\\}})\quad\forall
j\in\\{1,\dots,N_{y}\\}$ then

$\alpha=1$

else

$\alpha=0.2$

end if

MedianQ1Q3LowerWhiskerUpperWhiskerIQR1.5$\times$IQROutliers Figure 1:
Illustration of a typical box plot. The red dot denotes the observation.

##  6 Numerical Examples

We demonstrate the numerical behaviour of the proposed combination of the
aforementioned filters and the related complementary filtering techniques. Our
choice of experiments includes benchmarks based on the classical models
introduced in [37, 38] and the reference setup provided in [10]. Additionally,
we investigate the behaviour of existing and the proposed tempering criteria
using the Langevin-type stochastic differential equation originating from a
double well potential. We consider this scenario as prototypical for many
situations in high-dimensional state spaces, whenever the (particle) filter is
far from the asymptotic regime i.e. the number of particles is small in
comparison to the state space and the regularity of the filtering
distribution. Finally, we consider the one-dimensional shallow water equations
to provide a descriptive illustration and physical interpretation of the
proposed tempering criteria using a spatio-temporal model. Additionally, to
the baseline of each individual filter, we test three different tempering
approaches, c.f. Section 5. In the first one, we use the tempered combination
in each data assimilation step, denoted by ETPF-ESRF. The second one includes
the ESS as criteria to temper, i.e. ESS-ETPF-ESRF Section 5.4.1. Meaning that
the tempered version is only used in case the $ESS<0.5N_{ens}$. Last but not
least, we test the IQR approach, denoted by IQR-ETPF-ESRF, c.f. Section 5.4.2.
Meaning that we only use the tempered filter combination in case the
observation is not within the IQR. Before providing the first example, we
would like to point out the crucial additional strategies used to ensure the
filters are both robust and computationally practical.

###  6.1 Add-on filter techniques

Due to computational constraints that limit ensemble size, a critical
consideration is mitigating the effects of small ensembles. Specifically,
filters are susceptible to sampling errors, which diminish accuracy and lead
to underestimation of error covariances. This can result in filter divergence,
where the filter fails to track the reference solution effectively. Techniques
such as localization or inflation/rejuvenation can address these issues.

####  6.1.1 Inflation

Inflation, for instance, increases forecast uncertainty to enhance the
robustness of ensemble filter implementations against the effects of finite
ensemble sizes. Given a forecast ensemble $z_{i,k}^{f}$, we first compute its
empirical mean $\overline{z_{k}}^{f}$ and the ensemble anomalies $A_{k}^{f}$.
An inflation ensemble is defined according to [45] by

$z^{f}=\overline{z}^{f}+\gamma A^{f}$ (40)

whereby $\gamma>1$ is the inflation factor. Inflation is used to avoid filter
divergence for small ensemble sizes. However, inflation does not directly
account for spatial structures. [13] presented two further tuning tools
related to inflation, namely jittering and nudging. Jittering improves the
diversity of the ensemble by computing new ensembles, and nudging corrects the
solution term of the evolution model to keep ensembles closer to the true
state. An optimal nudging procedure in particle filters is given in [34].

####  6.1.2 Rejuvenation

A general sequential importance resampling filter can produce identical or
near identical ensemble members. This issue can be overcome according to [10]
by using particle rejuvenation in the following way:

$z_{j}^{a}=z_{j}^{a}+\sum_{i=1}^{N_{ens}}(z_{i}^{f}-\overline{z}^{f})\frac{\tau\xi_{ij}}{\sqrt{N_{ens}-1}},$
(41)

where $\tau>0$ denotes the rejuvenation parameter and $\xi_{ij}$ are iid
Gaussian random variables with mean zero and standard deviation one.
Rejuvenation is employed to prevent the occurrence of identical or nearly
identical ensemble members, while also increasing the ensemble’s spread.

####  6.1.3 Localisation

Difficulties arising when filtering spatio-temporal processes are twofold.
Firstly the ensemble size is much smaller than the state space and secondly,
spurious correlations are introduced as not all observations are physically
relevant for each point in the domain. Both issues can at least partly
mitigated by localisation, a technique, which determines the analysis state by
conducting ”local analyses” at each model grid point, considering only
observations within a specified nearby region, i.e. Figure 2. The seminal work
[27] provides a standard implementation for this technique applied to Ensemble
Kalman filter and is commonly used in numerical weather prediction and climate
analysis, as it realistically aligns with the underlying dynamics. The general
idea can be implemented in multiple ways. A prevalent technique involves
modifying the ensemble’s covariance matrix to mitigate spurious long-range
correlations. Known as B-localization, this method is widely acknowledged yet
challenging to apply uniformly across various filtering techniques, c.f. [45].
Conversely, R-localization incorporates only observations within a predefined
vicinity of the grid point under consideration. Similar to B-localization, the
selection of the effective radius for local regions in R-localization is
critical, as it must accurately reflect the spatial extent of significant
dynamical correlations manifested by the ensemble. Considering the differences
in how the R matrix and individual updates are adjusted across different
filtering families, we will specify the localisation schemes employed for the
families of filters addressed in this study. According to [45], the
localisation matrix $\tilde{C_{g}}\in\mathbb{R}^{N_{y}\times N_{y}}$ with
localisation radius $r_{loc}$ is given by:

$(\tilde{C_{g}})_{ll}=\rho\left(\frac{|x_{g}-\tilde{x_{l}}|_{L}}{r_{loc}}\right),$
(42)

with the current grid point $x_{g}=g\Delta x$ and the observation grid point
$\tilde{x_{l}}=l\Delta x$. $\rho$ is known as filtering function and can have
linear or higher polynomial order, e.g. [45]. Figure 2 shows an example of the
localisation technique described above. The green point is the current grid
point. If an observation is within the radius of the current grid point, it is
taken into account when calculating the localisation matrix. All other
observations are not of interest.

Figure 2:  Illustration of localization, Blue are the grid points, green the
currently updated point, orange the available observation and red the selected
obs

####  6.1.4 LESRF

The localised Ensemble Square Root Filter (LESRF) inserts the localisation
matrix (42) directly into the square root matrix $S$, c.f. [45]. In addition,
each grid point can now be treated separately. Parallel processing can now
speed up the algorithm. Using inflation as well as localisation as a tuning
tool, the algorithm of the ESRF, Appendix Algorithm 4, changes into:

Algorithm 2 Localised Ensemble Square Root Filter (LESRF) (LETKF, [27]) see
also the presentation in [45]

1:ensemble forecast $z^{f}_{i}\in\mathbb{R}^{N_{x}},\quad\forall i\in\\{1\dots
N_{ens}\\}$ and observation $y^{obs}\in\mathbb{R}^{N_{y}},\gamma=1.05$

2:$\overline{z}^{f}=\frac{1}{N_{ens}}\sum_{i=1}^{N_{ens}}z^{f}_{i}$

3:$A^{f}=(z^{f}_{i}-\overline{z}^{f})$

4:inflation: $z_{i}^{f}=\overline{z}^{f}+\gamma A^{f}$

5:compute step 1 and 2 again

6:for all $x_{g}\in\Omega_{N_{z}}$ do

7:
$S[x_{g}]=\left(I_{N_{ens}}+\frac{1}{N_{ens}-1}(HA^{f})^{T}\tilde{C_{g}}R^{-1}(HA^{f})\right)^{-1/2}$

8:
$\hat{w}_{i}[x_{g}]=\frac{1}{N_{ens}}-\frac{1}{N_{ens}-1}e_{i}^{T}T^{2}[x_{g}](HA^{f})^{T}\tilde{C_{g}}R^{-1}(H\overline{z}^{f}-y^{obs})$

9: $d_{ij}[x_{g}]=\hat{w}_{i}[x_{g}]-\frac{1}{N_{ens}}+t_{ij}[x_{g}]$

10: $z_{j}^{a}[x_{g}]=\sum_{i=1}^{N_{ens}}z^{f}_{i}[x_{g}]d_{ij}[x_{g}]$

11:end for

12: return $z^{a}_{j}$ for every $i\in\\{1\dots N_{ens}\\}$

####  6.1.5 LETPF

For the localised Ensemble Transform Particle Filter (LETPF), a localised
distance between two ensemble members at each grid point $x_{g}$ is defined as
described in [10].

$\tilde{d}[x_{g}](z_{i_{1}}^{f},z_{i_{2}}^{f})^{2}=\sum_{l=1}^{N_{d}}\rho\left(\frac{|x_{g}-\tilde{x_{l}}|_{L}}{r_{loc}}\right)d_{2}\left(z_{i_{1}}^{f}(\tilde{x_{l}}),z_{i_{2}}^{f}(\tilde{x_{l}})\right)^{2}$
(43)

This distance is a weighted average of the differences between the two
ensemble members over nearby gird points [45]. Using rejuvenation as well as
localisation as tuning tool, the algorithm of the ETPF, Appendix Algorithm 6,
changes into:

Algorithm 3 Localised Ensemble Transform Particle Filter (LETPF), see also the
presentation in [45]

1:ensemble forecast $z^{f}_{i}\in\mathbb{R}^{N_{x}},\quad\forall i\in\\{1\dots
N_{ens}\\}$ and observation $y^{obs}\in\mathbb{R}^{N_{y}},\tau=0.2$

2:$\overline{z_{k}}^{f}=\frac{1}{N_{ens}}\sum_{i=1}^{N_{ens}}z^{f}_{i}$

3:$A^{f}=(z^{f}_{i}-\overline{z}^{f})$

4:for all $x_{g}\in\Omega_{N_{z}}$ do

5:
$w_{i}[x_{g}]\propto\exp\left(-\frac{1}{2}(Hz^{f}_{i}-y^{obs})^{\top}\tilde{C_{g}}R^{-1}(Hz^{f}_{i}-y^{obs})\right)$

6: $\hat{w}_{i}[x_{g}]=\frac{w_{i}[x_{g}]}{\sum_{i=1}^{N_{ens}}w_{i}[x_{g}]}$

7: $T^{\ast}[x_{g}]\leftarrow$ solution of (17) using the localised metric
i.e. $d=\tilde{d}[x_{g}]$.

8: $rejuvenation=\frac{\tau}{N_{ens}-1}A^{f}\xi,\ \xi\sim\mathcal{N}(0,1)$

9:
$z^{a}_{i}[x_{g}]=\sum_{i=1}^{N_{ens}}z^{f}_{i}[x_{g}]N_{ens}t_{ij}^{\ast}[x_{g}]+rejuvenation$

10:end for

11: return $z^{a}_{i}$ for every $i\in\\{1\dots N_{ens}\\}$

###  6.2 Langevin Dynamics

As the first prototypical scenario, we consider the evolution of the dampened
Langevin equation:

$\displaystyle\begin{split}q_{t}&=M^{-1}p_{t}\,dt,\\\
p_{t}&=\nabla\phi(q_{t})\,dt-\gamma p_{t}+\sigma
M^{\frac{1}{2}}dW_{t}.\end{split}$ (44)

The deterministic forces are given by the gradient of the double well
potential $\phi(z)=(z-5)^{2}(z+5)^{2}$, and the stochastic forces by the
standard normal Wiener process $W_{t}$. We obtain the time discrete numerical
approximation by the BAOAB splitting strategy as presented in [33]. As
stepsize we choose $\Delta t=0.01$ and for the physical parameters we choose
$M=1$, $\gamma=10$, $\sigma=\sqrt{5000}$. Further, we observe all variables in
observation intervals of $\Delta t_{obs}=0.8$ with observation error
covariance $R=diag(0.5,0.5)$. For our purpose, a total of $50$ assimilation
steps are performed. We include ensemble inflation and particle rejuvenation,
see Sections 6.1.1 and 6.1.2. We set the inflation parameter $\gamma=1.05$,
the rejuvenation parameter $\tau=0.2$. In this setup, we compare ESRF, ETPF,
Bootstrap, FPF and their combination controlled by the tempering criteria
discussed in Section 5.
To understand the intricacies of the filter behaviour, we consider a concrete
scenario: Consider a system of particles with mean exit time from a potential
well comparable to the length of the observation interval. To investigate the
performance of the filters and their tempered versions under a change to
another potential well we initialise an ensemble of $N_{ens}=35$ particles.
Figure 3 provide a visual representation of the comparative performance
between the ESRF, ETPF, and their tempered combination in the described
scenario. Remarkably, the data indicates a notable advantage of the ESRF over
the ETPF, particularly evident in situations where particles undergo tail
switching dynamics, i.e. Figure 3. Further, in Figure 3 we observe the
tempered filter (ETPF-ESRF) to be competitive with the ESRF more likely than
the single ETPF. Here, we apply the three criteria for the tempering approach
as described in Section 5. The threshold for the ESS criterion is set to $0.5$
and we choose a factor of $1.5$ for the IQR criterion. All numerical results
are collected in Table 1

Figure 3: Visualisation of the prior and posterior distributions as well as
the truth and observation over two time steps (from left to right),
highlighting the notable tail switching dynamics between time steps.

Filter | RMSE
---|---
ESRF | 0.41761
ETPF | 1.99017
Bootstrap | 1.25118
FPF | 5.03897
ETPF-ESRF | 0.46150
Bootstrap-ESRF | 0.41130
FPF-ESRF | 0.41151
ESS-ETPF-ESRF | 0.46150
ESS-Bootstrap-ESRF | 0.41130
ESS-FPF-ESRF | 0.41151
IQR-ETPF-ESRF | 0.46150
IQR-Bootstrap-ESRF | 0.41130
IQR-FPF-ESRF | 0.41151

Table 1: Root Mean Squared Error (RMSE) for all filter combinations in the
Langevin experiment.

###  6.3 Lorenz 63

We use the chaotic Lorenz 63 system with the standard parameter setting,
whereby $\sigma=10,\rho=28,$ and $\beta=8/3$. We observe the first component
of the three-dimensional system in observation intervals of $\Delta
t_{obs}=0.12$ with observation error variance $R=8$. Our experiment covers
$50000$ assimilation steps, whereby the first $500$ are neglected for the
evaluation, to ensure the particles to be sufficiently close to the attractor.
We include ensemble inflation and particle rejuvenation, see Sections 6.1.1
and 6.1.2. Note that $\gamma=1.05$, $\tau=0.2$ as in the previous example.
This data assimilation setting has already been used in [1]. As in the
previous example, we test several combinations of filters and present the
results in Table 2. Remarkably, the approach using the IQR performs best when
comparing the RMSE of the single ETPF with the tempered version and beats even
the single ESRF. In general, the IQR approach holds better results then the
other two approaches. Figure 4 depicts first the box plots of the prior
observations for the single ETPF. As visible in assimilation cycle $5-13$, the
particles do not cover the observations. The second graph in Figure 4 shows
the tempered filter using the IQR criteria with quantile factor $1.5$.

Filter | RMSE
---|---
ESRF | 2.10011
ETPF | 3.55604
Bootstrap | 5.94889
FPF | 4.94511
ETPF-ESRF | 2.06520
Bootstrap-ESRF | 2.08702
FPF-ESRF | 2.22518
ESS-ETPF-ESRF | 1.75024
ESS-Bootstrap-ESRF | 3.12627
ESS-FPF-ESRF | 2.22518
IQR-ETPF-ESRF | 1.64179
IQR-Bootstrap-ESRF | 2.01076
IQR-FPF-ESRF | 2.22518

Table 2: Root Mean Squared Error (RMSE) for all filter combinations in the
Lorenz 63 model.

Figure 4: The box plots include the observations and the ensemble state at the
first state of the Lorenz 63 over a range of 50 assimilating cycles. The upper
and lower diagram depicts the statistics computed by the plain ETPF and the
tempered (IQR-criterion) version respectively.

###  6.4 Lorenz 96

The discrete dynamical system introduced in the seminal work by [38] is a
widely adopted model to provide spatio-temporal benchmarks for data
assimilation methods. The model originates from an unstable discretization of
the dampened Burgers’ equation with constant forcing. In contrast to a stable
discretization, the evolution in time becomes chaotic. We choose $120$ grid
points, the forcing to be $8$ and observe every second grid point in
observation intervals of $\Delta t_{obs}=0.11$. The observation error
covariance is $R=diag(8,\dots,8)$ and the localisation radius is equal to two
grid points. This experiment cover 5000 assimilation steps, whereby the first
50 are neglected. Further, $\gamma=1.05$, $\tau=0.2$, and $\alpha=0.2$.
Inflation, rejuvenation, and localisation can also be applied in the tempered
filters according to Sections 6.1.1, 6.1.2 and 6.1.3. Comparing the localised
filters, the RMSE varies more strongly. The minimum RMSE is obtained with the
localised ETPF (LETPF) as visible in Table 3. In general the tempered versions
show a slightly better RMSE than the single filters alone. Focusing on the
ETPF in Table 3, marked in blue, also shows that the different choices of
tempering lead to similar values of the RMSE. This holds true for all other
combinations as well and suggests the following heuristic. If the RMSE of the
single filters differs too much or too less, then their tempered version is
unable to substantially reduce the RMSE. Therefore, we refrain from
combinations of the LETPF and LESRF, as they both perform equally well
individually. We emphasise the substantial improvement for the FPF, marked in
orange in Table 3, especially in contrast to the other particle filters and
their tempered versions. Finally, we observe again the increased spread in the
ensemble in Figure 5 for the plain ETPF and the tempered version using the
LESRF and the IQR criteria with quantile factor zero.

Filter | RMSE
---|---
ESRF | 2.48515
LESRF | 1.08808
ETPF | 2.96073
LETPF | 1.05996
Bootstrap | 2.96868
FPF | 2.15463
ETPF-ESRF | 2.86639
Bootstrap-ESRF | 2.88486
FPF-ESRF | 1.98065
ESS-ETPF-ESRF | 2.86633
ESS-Bootstrap-ESRF | 2.91596
ESS-FPF-ESRF | 1.98306
IQR-ETPF-ESRF | 2.86639
IQR-Bootstrap-ESRF | 2.88486
IQR-FPF-ESRF | 1.98065
ETPF-LESRF | 2.51371
Bootstrap-LESRF | 2.60774
FPF-LESRF | 1.37967
ESS-ETPF-LESRF | 2.53425
ESS-Bootstrap-LESRF | 2.62288
ESS-FPF-LESRF | 1.37967
IQR-ETPF-LESRF | 2.51371
IQR-Bootstrap-LESRF | 2.60774
IQR-FPF-LESRF | 1.37967

Table 3: Root Mean Squared Error (RMSE) for all filter combinations in the
Lorenz 96 model.

Figure 5: The box plots include the observations and the ensemble state at the
first grid point of the Lorenz 96 discretization over a range of 50
assimilating cycles. The upper and lower diagram depicts the statistics
computed by the plain ETPF and the tempered (IQR-criterion) version
respectively.

###  6.5 Shallow water

The shallow water system is a set of nonlinear hyperbolic partial differential
equations that describe the fluid flow below a pressure surface. The equations
are derived from depth-integrating the Navier-Stokes equations. Hence,
conservation of mass implies that the vertical velocity scale of the fluid is
small compared to the horizontal velocity scale. The characteristic speed of
the shallow water waves depends on the water height above ground and which in
turn leads to nonlinear effects e.g. wave breaking.

Let $\Omega=[0,L]\subset\mathbb{R}$ then shallow water system with gravity $g$
and orography $z$ reads

$\displaystyle\begin{split}h_{t}+\nabla\cdot hu&=0,\\\
(hu)_{t}+\nabla\cdot(hu\otimes u)+gh\nabla h&=-gh\nabla z,\end{split}$ (45)

where $h$ denotes the height over ground i.e. the total height is given by
$H=h+z$ and $u$ denotes the two-dimensional velocity vector. On the
boundaries, we assume a wall at $x=L$ and homogenous Neumann boundary
conditions at $x=0$. A solution to the initial value problem, satisfying
$u=\vec{0}$ and $\nabla(z+h)=\vec{0}$ for all times is called the lake at rest
steady state.

Despite being reasonably popular for the modelling of nonlinear water waves
and internal gravity or Rossby waves in the atmosphere [39], the shallow water
equations, a-priori, are a questionable choice as a model for data
assimilation due to the lack of any dissipative mechanism as e.g. in the
Navier-Stokes case (c.f. [29]). Nevertheless, the model is sufficient to serve
as a high-dimensional example for strong nonlinear effects. In contrast to the
Lorenz 96 model, the spatial solutions are stable and allow for a meaningful
physical interpretation.

$x$$H,z$$0$$x_{w}$$x_{r}$$x_{d}$$x_{b}$$L$$z_{d}$$H_{0}$$H_{w}$$z_{b}$

Figure 6: Initial condition of the used shallow water system.

We compute a numerical approximation of solutions to the initial valued
problem of the shallow water equations by a second-order centred finite volume
scheme presented in [6] using hydrostatic reconstruction developed in [3, 5].
This enables us to sample from the prior distribution efficiently and
robustly, even in the given scenario with wet-dry interfaces, i.e. when the
non-negative $h$ vanishes in some parts of the domain. Furthermore, the finite
volume method preserves the mass of the fluid exactly. It is important to
denote that although the actual filtering distribution is non-negative, some
filter modifications as rejuvenation might set samples of the posterior to
negative values. In this case, we cut off the values at zero. We choose the
initial condition of the reference solution according to Figure 6. The
discretization is spatially uniform on $[0,L-x_{b})$ with $49$ grid cells and
contains a single grid cell covering $[L-x_{b},L]$. The specific values are
$x_{r}=10$, $x_{d}=20$, $x_{b}=24$, $L=28$, $H_{0}=1.5$, $H_{w}=1.8$,
$z_{b}=0.5$, $z_{d}=2.2$. For the reference solution $x_{w}=3$ and the
ensemble height is distributed according to $x_{w}\approx\mathcal{N}(6,1)$,
i.e. the initial condition is misspecified on the $99.87\%$ percentile. The
velocity of the initial ensemble is generated by $v=\sqrt{g(H_{w}-H_{0})}$,
where $g=9.81$. The observational window is given by $\Delta t_{obs}=2.5$ and
we observe the fluid height at the centre of every 8th grid cell. We run 3
assimilation cylces, but skip the first for collecting the statistics. Note
that $\gamma=1.05$, $\tau=0.2$, and $\alpha=0.2$ as usual. We evaluate the
RMSE in the large cell covering the basin behind the dam and average over the
two assimilation cycles. The collected results are listed in Table 4.
Furthermore, we depict the ensemble statistics for the rightmost point of
observation in Figure 7. Although neither the RMSE nor the ensemble statistics
change drastically, we observe a similar trend as in previous examples.
Interestingly enough, we can spatially locate the differences between the
plain and the tempered version of the filter, where the particle filters, by
construction of the orography, have reduced spread in the distribution.

Filter | RMSE
---|---
ESRF | 0.038228
ETPF | 0.049417
Bootstrap | 0.056183
FPF | 0.071165
LESRF | 0.048350
LETPF | 0.056134
ETPF-ESRF | 0.040833
Bootstrap-ESRF | 0.030875
FPF-ESRF | 0.053365
ESS-ETPF-ESRF | 0.051588
ESS-Bootstrap-ESRF | 0.056429
ESS-FPF-ESRF | 0.053365
IQR-ETPF-ESRF | 0.040833
IQR-Bootstrap-ESRF | 0.030875
IQR-FPF-ESRF | 0.053365

Table 4: Root Mean Squared Error (RMSE) for all filter combinations in the
shallow water experiment. Figure 7: The box plots include the observations and
the ensemble state at the height $x=23.8$ of the shallow water system
comparing the ETPF and the IQR version respectively.

##  7 Discussion

The focus of this study is to better the performance of particle filters,
especially those described in Sections 4.1, 4.2.1 and 4.3.1, by the means of
tempering Section 5. The novel approach presented in Section 5.4.2, receives
the greatest attention and differs from the approaches developed in [22, 10,
40]. Specifically, [40] posed two choices of the bridging parameter for the
tempering approach, First, depending on the effective sample size as
recommended in [22], second, on skewness and kurtosis of the empirical
distribution in observational space. This should quantify the non-Gaussianity
of the ensemble distribution and as pointed out in, this shows promising
results, but requires further fine-tuning. This work is not primarily
concerned with the bridging parameter, but rather with the question whether to
temper or not. Targeted to the challenging situations for particle type
filters, we proposed a novel tempering criterion, based on descriptive
statistics, which is straightforward to compute and has a clear
interpretation, i.e. if the observation is within the interquartile range of
the observed particles. This criterion is complementary to a threshold on the
effective sample size, as the latter does not take the observation into
account. As particle filters in high-dimensional systems usually suffer from
the lack of overlapping support of the likelihood and the empirical prior
distribution, we propose to temper the likelihood by a Gaussian-approximative
filter to push the particles possibly towards the observation. This
prototypical situation is illustrated in our experiment by applying several
particle filters to a filtering scenario based on the double well potential
c.f. Figure 3 and Table 1. As Figure 3 visualises, the ETPF struggles clearly
to switch the modes of the double well. Furthermore, it is visible that the
tempered approach shifts the ensembles closer to the observation. Table 1
points out significant improvements for all three particle type filters,
highlighted in blue for the ETPF, green for the Bootstrap and orange for the
FPF. Note that there is no indication which of the tempering approaches serves
as best one. Further, we test the tempering approaches on different particle
filters, tuning them with rejuvenation and localisation as needed, and
provided tests on the defacto standard benchmarks based on both Lorenz models.
In the case of the Lorenz 63 model, we observe improvements for the time-
averaged root mean square error, for all tempered versions, but find the IQR
criterion to perform best as listed in Table 2. IQR-ETPF-ESRF demonstrated the
lowest RMSE ($1.64179$), surpassing even the performance of the single ESRF
($2.10011$). The IQR criteria also outperform the other two tempering
approaches. This finding counts also for the other two particle type filters
as highlighted in Table 2. In the Lorenz 96 experiments, the improvements are
clearly visible, but not as pronounced as in the previously mentioned case.
Considering Table 3, the LETPF exhibited the lowest RMSE ($1.05996$),
demonstrating that localisation approaches can increase filter accuracy. As
previously mentioned, we focus on the particle type filters (without
localisation) and their tempered variations. Table 3 shows very small
improvements for the ETPF and Bootstrap. Comparing the RMSE of the IQR-FPF-
LESRF ($1.37967$), which outperforms the FPF ($2.15463$), exhibited a
substantial improvement in performance, suggesting that the incorporation of
IQR alongside FPF and LESRF can yield significant enhancements. Further,
Figure 5 illustrates the advantages of using IQR, highlighting its
effectiveness in enhancing filtering accuracy. Finally, we provide a test case
based on the shallow water equations with orography. Here, we investigate
essentially two cycles of a constructed filtering scenario and would like to
highlight, that this scenario provides a situation, where the criterion based
on the ESS provides worse results than the constant tempering or the IQR based
criterion, i.e. Table 4. This points out the robustness of the IQR criterion
across all experiments in contrast to the ESS, which does not target the
specific difficulties under investigation. Furthermore, Figure 7 illustrates
again the shift towards the observation, localised at the specific subspace
(gridpoint).

##  8 Conclusion

Tempering significantly enhances the performance of filtering methods,
particularly for approximate consistent filters. Additionally, it can help
stabilise numerically unstable filters through occasional robust updates,
using a combination of filters rather than a single fixed choice. However,
determining the optimal timing for employing these tools during filtering
procedures remains unclear. We propose a novel approach to adaptively decide
on the use of an approximate tempering step, using the interquartile range
(IQR). To numerically validate this tempering criterion across various
filters, we consider three representative filters from key families of
consistent filters: standard sequential Monte Carlo resampling, Transport-map-
induced, and McKean-Vlasov filters. Each of these filters is paired with an
established representative from the class of Gaussian approximate filters.
Overall, the new adaptive tempering approach of the considered consistent
filters paired with adaptive tempering steps resulted in lower or equal RMSE
compared to the performance of the individual filter, a fixed tempering cycle,
and an additional effective sample size criterion. This effect is evident
across three different types of toy models, each with a specific focus. The
promising results obtained in this study suggest further investigation into
more realistic filtering scenarios based on atmospheric fluid dynamics.
Another important aspect will be a careful examination of the numerical
performance of the Feedback Particle Filter in the non-asymptotic regime with
small ensemble sizes.

## Acknowledgments

This research was funded in whole or in part by the Austrian Science Fund
(FWF) [10.55776/DOC78]. For open access purposes, the author has applied a CC
BY public copyright license to any author-accepted manuscript version arising
from this submission.
The research of Gottfried Hastermann, Iris Rammelmüller and Jana de Wiljes has
been partially funded by the Deutsche Forschungsgemeinschaft (DFG)- Project-ID
318763901 - SFB1294.
The program code is available as ancillary file from the arXiv page of this
paper.

## References

  * [1] Walter Acevedo, Jana De Wiljes, and Sebastian Reich.  Second-order Accurate Ensemble Transform Particle Filters.  SIAM Journal on Scientific Computing, 39(5):A1834–A1850, January 2017.  doi:10.1137/16M1095184.
  * [2] Ö. D. Akyildiz and J. Míguez.  Nudging the particle filter.  Statistics and Computing, 30:305–330, 2020.  doi:10.1007/s11222-019-09884-y.
  * [3] Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoît Perthame.  A Fast and Stable Well-Balanced Scheme with Hydrostatic Reconstruction for Shallow Water Flows.  SIAM Journal on Scientific Computing, 25(6):2050–2065, January 2004.  doi:10.1137/S1064827503431090.
  * [4] Alexandros Beskos, Dan Crisan, Ajay Jasra, Kengo Kamatani, and Yan Zhou.  A stable particle filter for a class of high-dimensional state-space models.  Advances in Applied Probability, 49(1):24–48, 2017.  doi:10.1017/apr.2016.77.
  * [5] N. Botta, R. Klein, S. Langenberg, and S. Lützenkirchen.  Well Balanced Finite Volume Methods for Nearly Hydrostatic Flows.  Journal of Computational Physics, 196(2):539–565, May 2004.  doi:10.1016/j.jcp.2003.11.008.
  * [6] François Bouchut.  Chapter 4 Efficient Numerical Finite Volume Schemes for Shallow Water Models.  In Edited Series on Advances in Nonlinear Science and Complexity, volume 2, pages 189–256. Elsevier, 2007.  doi:10.1016/S1574-6909(06)02004-1.
  * [7] Gerrit Burgers, Peter Jan van Leeuwen, and Geir Evensen.  Analysis Scheme in the Ensemble Kalman Filter.  Monthly Weather Review, 126(6):1719 – 1724, 1998.  URL: https://journals.ametsoc.org/view/journals/mwre/126/6/1520-0493_1998_126_1719_asitek_2.0.co_2.xml, doi:10.1175/1520-0493(1998)126<1719:ASITEK>2.0.CO;2.
  * [8] Y. Cheng and S. Reich.  A McKean optimal transportation perspective on Feynman-Kac formulae with application to data assimilation.  arXiv preprint arXiv:1311.6300, 2013.  doi:10.48550/arXiv.1311.6300.
  * [9] N. Chopin and O. Papaspiliopoulos.  An introduction to sequential Monte Carlo.  Berlin: Springer, Berlin, 2020.  doi:10.1007/978-3-030-47845-2.
  * [10] Nawinda Chustagulprom, Sebastian Reich, and Maria Reinhardt.  A Hybrid Ensemble Transform Particle Filter for Nonlinear and Spatially Extended Dynamical Systems.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):592–608, January 2016.  doi:10.1137/15M1040967.
  * [11] R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner, and S. W. Zucker.  Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps.  Proceedings of the National Academy of Sciences, 102(21):7426–7431, May 2005.  doi:10.1073/pnas.0500334102.
  * [12] Ronald R. Coifman and Stéphane Lafon.  Diffusion Maps.  Applied and Computational Harmonic Analysis, 21(1):5–30, July 2006.  doi:10.1016/j.acha.2006.04.006.
  * [13] Colin Cotter, Dan Crisan, Darryl Holm, Wei Pan, and Igor Shevchenko.  Data Assimilation for a Quasi-Geostrophic Model with Circulation-Preserving Stochastic Transport Noise.  Journal of Statistical Physics, 179(5-6):1186–1221, June 2020.  doi:10.1007/s10955-020-02524-0.
  * [14] Colin Cotter, Dan Crisan, Darryl D Holm, Wei Pan, and Igor Shevchenko.  A Particle Filter for Stochastic Advection by Lie Trancnort: A Case Study for the Damped and Forced Incomprebbible Two-Dimemiunal Euler Equation.  SIAM/ASA Journal on Uncertainty Quantification, 8(4):1446–1492, 2020.  doi:10.1137/19M1277606.
  * [15] M. Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  Advances in neural information processing systems, 26, 2013.
  * [16] J. de Wiljes and X.T. Tong.  Analysis of a localised nonlinear Ensemble Kalman Bucy Filter with complete and accurate observations.  Nonlinearity, 33(9):4752–4782, 2020.  doi:10.1088/1361-6544/ab8d14.
  * [17] Jana de Wiljes, Sebastian Reich, and Wilhelm Stannat.  Long-Time Stability and Accuracy of the Ensemble Kalman–Bucy Filter for Fully Observed Processes and Small Measurement Noise.  SIAM Journal on Applied Dynamical Systems, 17(2):1152–1181, January 2018.  doi:10.1137/17M1119056.
  * [18] Tarek A. El Moselhy and Youssef M. Marzouk.  Bayesian Inference with Optimal Maps.  Journal of Computational Physics, 231(23):7815–7850, October 2012.  doi:10.1016/j.jcp.2012.07.022.
  * [19] G. Evensen.  Sequential data assimilation with a nonlinear quasi‐geostrophic model using Monte Carlo methods to forecast error statistics.  Journal of Geophysical Research: Oceans, 99(C5):10143–10162, 1994.  doi:10.1029/94JC00572.
  * [20] Geir Evensen, Femke C. Vossepoel, and Peter Jan Van Leeuwen.  Data Assimilation Fundamentals: A Unified Formulation of the State and Parameter Estimation Problem.  Springer Textbooks in Earth Sciences, Geography and Environment. Springer International Publishing, Cham, 2022.  doi:10.1007/978-3-030-96709-3.
  * [21] P. Fearnhead and H. R. Künsch.  Particle filters and data assimilation.  Annual Review of Statistics and Its Application, 5:421–449, 2018.  doi:10.1146/annurev-statistics-031017-100232.
  * [22] M. Frei and H. R. Künsch.  Bridging the Ensemble Kalman and Particle Filters.  Biometrika, 100(4):781–800, December 2013.  doi:10.1093/biomet/ast020.
  * [23] Mohamad El Gharamti.  Hybrid Ensemble–Variational Filter: A Spatially and Temporally Varying Adaptive Algorithm to Estimate Relative Weighting.  Mobthly weather review, 149(1):65–76, 2021.  doi:10.1175/MWR-D-20-0101.1.
  * [24] A. Grigor’yan.  Heat kernels on weighted manifolds and applications.  Cont. Math, 398(2006):93–191, 2006.  doi:10.1090/conm/398/07486.
  * [25] Truong-Vinh Hoang, Sebastian Krumscheid, Hermann G. Matthies, and Raúl Tempone.  Machine learning-based conditional mean filter: A generalization of the ensemble Kalman filter for nonlinear data assimilation.  Foundations of Data Science, 5(1):56–80, 2023.  doi:10.3934/fods.2022016.
  * [26] Hakon Hoel, Gaukhar Shaimerdenova, and Raúl Tempone.  Multilevel Ensemble Kalman Filtering based on a sample average of independent EnKF estimators.  Foundations of Data Science, 2(4):351–390, 2020.  doi:10.3934/fods.2020017.
  * [27] Brian R. Hunt, Eric J. Kostelich, and Istvan Szunyogh.  Efficient Data Assimilation for Spatiotemporal Chaos: A Local Ensemble Transform Kalman Filter.  Physica D: Nonlinear Phenomena, 230(1-2):112–126, June 2007.  doi:10.1016/j.physd.2006.11.008.
  * [28] A. H. Jazwinski.  Stochastic processes and filtering theory.  Courier Corporation, 2022.
  * [29] Don A Jones and Edriss S Titi.  On the Number of Determining Nodes for the 2D Navier-Stokes Equations.  Journal of Mathematical Analysis and Applications, 168(1):72–88, July 1992.  doi:10.1016/0022-247X(92)90190-O.
  * [30] R. E. Kalman.  A new approach to linear filtering and prediction problems.  Journal of Basic Engineering, pages 35–45, 1960.  doi:10.1115/1.3662552.
  * [31] Matthias Katzfuss, Jonathan R. Stroud, and Christopher K. Wikle.  Understanding the Ensemble Kalman Filter.  The American Statistician, 70(4):350–357, October 2016.  doi:10.1080/00031305.2016.1141709.
  * [32] Kody Law, Andrew Stuart, and Konstantinos Zygalakis.  Data Assimilation: A Mathematical Introduction, volume 62 of Texts in Applied Mathematics.  Springer International Publishing, Cham, 2015.  doi:10.1007/978-3-319-20325-6.
  * [33] B. Leimkuhler and C. Matthews.  Rational Construction of Stochastic Numerical Methods for Molecular Sampling.  Applied Mathematics Research eXpress, page abs010, June 2012.  doi:10.1093/amrx/abs010.
  * [34] N. Lingala, N. Perkowski, H.C. Yeong, N. Sri Namachchivaya, and Z. Rapti.  Optimal Nudging in Particle Filters.  Probabilistic Engineering Mechanics, 37:160–169, July 2014.  doi:10.1016/j.probengmech.2013.08.007.
  * [35] J. S. Liu and J. S. Liu.  Monte Carlo strategies in Scientific Computing.  New York: springer, 2001.  doi:10.1007/978-0-387-76371-2.
  * [36] David M. Livings, Sarah L. Dance, and Nancy K. Nichols.  Unbiased Ensemble Square Root Filters.  Physica D: Nonlinear Phenomena, 237(8):1021–1028, June 2008.  doi:10.1016/j.physd.2008.01.005.
  * [37] E. N. Lorenz.  Deterministic nonperiodic flow.  Deterministic nonperiodic flow. Journal of atmospheric sciences, 20(2):130–141, 1963.  doi:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2.
  * [38] E. N. Lorenz.  Predictability: A problem partly solved.  In Proc. Seminar on predictability, Vol. 1(No. 1), 1996.  doi:10.1017/CBO9780511617652.004.
  * [39] Andrew Majda.  Introduction to PDEs and Waves for the Atmosphere and Ocean, volume 9 of Courant Lecture Notes.  American Mathematical Society, Providence, Rhode Island, January 2003.  doi:10.1090/cln/009.
  * [40] L. Nerger.  Data assimilation for nonlinear systems with a hybrid nonlinear Kalman ensemble transform filter.  Quarterly Journal of the Royal Meteorological Society, 148(743):620–640, December 2021.  doi:10.1002/qj.4221.
  * [41] G. Papamakarios, E. Nalisnick, D. Rezende, S. Mohamed, and B. Lakshminarayanan.  Normalizing Flows for Probabilistic Modeling and Inference.  Journal of Machine Learning Research, 20:1–64, 2021.  doi:10.5555/3546258.3546315.
  * [42] Sahani Pathiraja, Sebastian Reich, and Wilhelm Stannat.  McKean-Vlasov SDEs in nonlinear filtering.  SIAM J. Control Optim., 59(6):4188–4215, 2021.  doi:10.1137/20M1355197.
  * [43] Sahani Pathiraja and Wilhelm Stannat.  Analysis of the feedback particle filter with diffusion map based approximation of the gain.  Foundations of Data Science, 3(3):615, 2021.  doi:10.3934/fods.2021023.
  * [44] Sebastian Reich.  A Nonparametric Ensemble Transform Method for Bayesian Inference.  SIAM Journal on Scientific Computing, 35(4):A2013–A2024, January 2013.  doi:10.1137/130907367.
  * [45] Sebastian Reich and Colin Cotter.  Probabilistic Forecasting and Bayesian Data Assimilation.  Cambridge University Press, Cambridge, 2015.  doi:10.1017/CBO9781107706804.
  * [46] Francesco A. B. Silva, Cecilia Pagliantini, Martin Grepl, and Karen Veroy.  A reduced basis ensemble Kalman method.  GEM - International Journal on Geomathematics, 14(1), 2023.  doi:10.1007/s13137-023-00235-8.
  * [47] Alessio Spantini, Ricardo Baptista, and Youssef Marzouk.  Coupling Techniques for Nonlinear Ensemble Filtering.  SIAM Review, 64(4):921–953, 2022.  arXiv:https://doi.org/10.1137/20M1312204, doi:10.1137/20M1312204.
  * [48] Amirhossein Taghvaei, Jana De Wiljes, Prashant G. Mehta, and Sebastian Reich.  Kalman Filter and Its Modern Extensions for the Continuous-Time Nonlinear Filtering Problem.  Journal of Dynamic Systems, Measurement, and Control, 140(3):030904, March 2018.  doi:10.1115/1.4037780.
  * [49] Amirhossein Taghvaei and Prashant G. Mehta.  Gain Function Approximation in the Feedback Particle Filter.  In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 5446–5452, Las Vegas, NV, USA, December 2016. IEEE.  doi:10.1109/CDC.2016.7799105.
  * [50] Amirhossein Taghvaei, Prashant G. Mehta, and Sean P. Meyn.  Diffusion Map-based Algorithm for Gain Function Approximation in the Feedback Particle Filter.  SIAM/ASA Journal on Uncertainty Quantification, 8(3):1090–1117, January 2020.  doi:10.1137/19M124513X.
  * [51] M. K. Tippett, J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker.  Ensemble Square Root Filters.  Monthly Weather Review, 131(7):1485–1490, 2003.  doi:10.1175/1520-0493(2003)131<1485:ESRF>2.0.CO;2.
  * [52] T. Yang, P. G. Mehta, and S. P. Meyn.  Feedback Particle Filter with Mean-field Coupling.  50th IEEE conference on decision and control and European control conference, pages 7909–7916, 2011.  doi:10.1109/CDC.2011.6160950.
  * [53] T. Yang, P. G. Mehta, and S. P. Meyn.  Feedback Particle Filter.  IEEE transactions on Automatic control, 58(10):2465–2480, 2013.  doi:10.1109/TAC.2013.2258825.

##  9 Appendix

Algorithm 4 Ensemble Square Root Filter (ESRF), see also the presentation in
[45] and Section 3.1

1:ensemble forecast $z^{f}_{i}\in\mathbb{R}^{N_{x}},\quad\forall
i\in\\{1,\dots,N_{ens}\\}$ and observation $y^{obs}\in\mathbb{R}^{N_{y}}$,
$\gamma=1.05$

2:$\overline{z}^{f}=\frac{1}{N_{ens}}\sum_{i=1}^{N_{ens}}z^{f}_{i}$

3:$A^{f}=(z^{f}_{i}-\overline{z}^{f})$

4:inflation: $z_{i}^{f}=\overline{z}^{f}+\gamma A^{f}$

5:compute step 1 and 2 again

6:$S=\left(I_{N_{ens}}+\frac{1}{N_{ens}-1}(HA^{f})^{T}R^{-1}(HA^{f})\right)^{-1/2}$

7:$\hat{w}_{i}=\frac{1}{N_{ens}}-\frac{1}{N_{ens}-1}e_{i}^{T}S^{2}(HA^{f})^{T}R^{-1}(H\overline{z}^{f}-y^{obs})$

8:$d_{ij}=\hat{w}_{i}-\frac{1}{N_{ens}}+s_{ij}$

9:$z_{j}^{a}=\sum_{i=1}^{N_{ens}}z^{f}_{i}d_{ij}$

10: return $z^{a}_{j}$ for every $j\in\\{1,\dots,N_{ens}\\}$

Algorithm 5 Bootstrap Particle Filter, see Section 4.1

1:initial ensemble $z^{f}_{i}\in\mathbb{R}^{N_{x}},\quad\forall
i\in\\{1,\dots,N_{ens}\\}$ and observation
$y^{obs}\in\mathbb{R}^{N_{y}},\tau=0.2$

2:$\overline{z}^{f}=\frac{1}{N_{ens}}\sum_{i=1}^{N_{ens}}z^{f}_{i}$

3:$A^{f}=(z^{f}_{i}-\overline{z}^{f})$

4:$\hat{w}_{i}\propto\exp\left(-\frac{1}{2}(Hz^{f}_{i}-y^{obs})^{\top}R^{-1}(Hz^{f}_{i}-y^{obs})\right)\cdot
w_{0},\ \ \ w_{0}=\frac{1}{N_{ens}}$

5:$w=\frac{\hat{w}_{i}}{\sum_{i=1}^{N_{ens}}\hat{w}_{i}}$

6:$rejuvenation=\frac{\tau}{N_{ens}-1}A^{f}\xi,\ \xi\sim\mathcal{N}(0,1)$

7:$z^{a}_{i}=rng.choice(z^{f}_{i},N_{ens},p=w,axis=1)+rejuvenation$

8: return $z^{a}_{i}$ for every $i\in\\{1,\dots,N_{ens}\\}$

Algorithm 6 Ensemble Transform Particle Filter (ETPF), see also the
presentation in [45] and Section 4.2.1

1:ensemble forecast $z^{f}_{i}\in\mathbb{R}^{N_{x}},\quad\forall
i\in\\{1,\dots,N_{ens}\\}$ and observation
$y^{obs}\in\mathbb{R}^{N_{y}},\tau=0.2$

2:$\overline{z}^{f}=\frac{1}{N_{ens}}\sum_{i=1}^{N_{ens}}z^{f}_{i}$

3:$A^{f}=(z^{f}_{i}-\overline{z}^{f})$

4:$w_{i}\propto\exp\left(-\frac{1}{2}(Hz^{f}_{i}-y^{obs})^{\top}R^{-1}(Hz^{f}_{i}-y^{obs})\right)$

5:$\hat{w}_{i}=\frac{w_{i}}{\sum_{i=1}^{N_{ens}}w_{i}}$

6:$T^{\ast}\leftarrow$ solution of (17) using the Euclidean metric i.e.
$d=d_{2}$.

7:$rejuvenation=\frac{\tau}{N_{ens}-1}A^{f}\xi,\ \xi\sim\mathcal{N}(0,1)$

8:$z^{a}_{i}=\sum_{i=1}^{N_{ens}}z^{f}_{i}N_{ens}t_{ij}^{\ast}+rejuvenation$

9: return $z^{a}_{i}$ for every $i\in\\{1,\dots,N_{ens}\\}$

Algorithm 7 Feedback Particle Filter (FPF), see also the presentation in [50]
and Section 4.3.1

1:ensemble forecast $z^{f}_{i}\in\mathbb{R}^{N_{x}},\quad\forall
i\in\\{1,\dots,N_{ens}\\}$ and observation $y^{obs}\in\mathbb{R}^{N_{y}}$,
$\epsilon=0.01$, and observation operator $H(\cdot)$, $\phi_{0}$,
$\mathcal{T}$ iteration, $\tau=0.2$

2:$g_{ij}=\exp{\frac{-|z^{f}_{i}-z^{f}_{j}|^{2}}{4\epsilon}}$

3:$l_{i}=\exp\left(-\frac{1}{2}(Hz^{f}_{i}-y^{obs})^{\top}R^{-1}(Hz^{f}_{i}-y^{obs})\right)$

4:$k_{ij}=\frac{g_{ij}}{\sqrt{\sum_{l}g_{il}}\sqrt{\sum_{l}g_{lj}}}\ l_{j}$

5:$d_{i}=\sum_{j}k_{ij}$

6:$T_{ij}=\frac{k_{ij}}{d_{i}}$

7:$\pi_{i}=\frac{d_{i}}{\sum_{j}d_{j}}$

8:$\hat{H}=\sum_{i=1}^{N_{ens}}\pi_{j}H(z^{f}_{i})$

9:$\phi_{t}=\textnormal{solution of
}(\mathop{id}-T+\mathbf{1})\phi_{t}=\varepsilon R^{-1}(H(z^{f})-\hat{H})$

10:$r_{i}=\phi_{\mathcal{T}}+\epsilon H(z^{f})$

11:$s_{ij}=\frac{1}{2\epsilon}T_{ij}(r_{j}-\sum_{k=1}^{N_{ens}}T_{ik}r_{k})$

12:$K_{i}=\sum_{j}s_{ij}z^{f}_{j}$

13:$A^{f}=(z^{f}_{i}-\overline{z}^{f})$

14:$rejuvenation=\frac{\tau}{N_{ens}-1}A^{f}\xi,\ \xi\sim N(0,1)$

15:$z^{a}=z^{f}-K_{i}(\frac{1}{2}(H(z^{f})-H(\overline{z_{k}}^{f}))-y^{obs})+rejuvenation$

16: return $z^{a}_{i}$ for every $i\in\\{1,\dots,N_{ens}\\}$