# Differentiable Particle Filtering using Optimal Placement Resampling

Domonkos Csuzdi1, Olivér Törő1 and Tamás Bécsi1  1Department of Control for
Transportation and Vehicle Systems, Faculty of Transportation Engineering and
Vehicle Engineering,
Budapest University of Technology and Economics, Muegyetem rkp. 3., Budapest,
1111, Hungary
Email<EMAIL_ADDRESS>
###### Abstract

Particle filters are a frequent choice for inference tasks in nonlinear and
non-Gaussian state-space models. They can either be used for state inference
by approximating the filtering distribution or for parameter inference by
approximating the marginal data (observation) likelihood. A good proposal
distribution and a good resampling scheme are crucial to obtain low variance
estimates. However, traditional methods like multinomial resampling introduce
nondifferentiability in PF-based loss functions for parameter estimation,
prohibiting gradient-based learning tasks. This work proposes a differentiable
resampling scheme by deterministic sampling from an empirical cumulative
distribution function. We evaluate our method on parameter inference tasks and
proposal learning.

##  I Introduction

Particle filters (PFs) are numeric approximation methods based on Monte Carlo
sampling for the general nonlinear Bayesian estimation problem. Monte Carlo
methods have been around since the earliest days of computation [1], with the
Metropolis-Hastings algorithm being one of the most notable [2]. The PF,
alternatively referred to as a sequential Monte Carlo estimator, is a
recursive method introduced in [3]. To evade numerical issues and work
efficiently, a PF needs careful design [4], meaning the algorithm may include
adaptive resampling, regularization, or progressive update. Without at least a
resampling step, a PF is ineffective.

Neural networks can replace parts of a PF that are hard to design or
parameterize, e.g., the proposal distribution. To train a neural network, one
usually applies backpropagation, which is basically a gradient computation.
This means that a standard PF structure cannot be used as it almost always
involves a nondifferentiable resampling step. Recently, there has been a focus
on differentiable particle filters as a consequence.

###  I-A State-Space Models

State-space models (SSMs) are a class of sequential models, where a latent
state $x_{t}$ evolves by a Markov process, and generates observations $y_{t}$.
We are interested in SSMs, where there are no external inputs, and the
observations are also Markovian, i.e., $y_{t}$ is conditionally independent of
the other states and measurements, given the current state $x_{t}$. These
relations are shown in Figure 1.

Figure 1: Directed probabilistic graphical model of an SSM. Shaded nodes are
observed random variables.

SSMs can also be represented by a parameterized generative model, defined as

$\displaystyle x_{1}$ $\displaystyle\sim\pi_{\theta}(x_{1})\,,$ (1)
$\displaystyle x_{t}$ $\displaystyle\sim f_{\theta}(x_{t}|x_{t-1})\,,$ (2)
$\displaystyle y_{t}$ $\displaystyle\sim g_{\theta}(y_{t}|x_{t})\,,$ (3)

where $\pi_{\theta}(x_{1})$ is the initial distribution of the latent state,
$f_{\theta}(x_{t}|x_{t-1})$ is the transition model, $g_{\theta}(y_{t}|x_{t})$
is the observation model, and $\theta$ is a parameter set of interest.

There are generally two main goals when working with SSMs: performing
inference over latent states, and performing inference over model parameters
$\theta$. The former has three forms based on which probability density
function (pdf) is being computed: $p(x_{t}|y_{1:t})$ for filtering,
$p(x_{t+1}|y_{1:t})$ for prediction and $p(x_{t}|y_{1:T})$ for smoothing.
Inference over $\theta$ can be done by maximum likelihood estimation using the
log of the marginal data likelihood:

$\displaystyle\log p_{\theta}(y_{1:T})$ $\displaystyle=\log\int
p_{\theta}(x_{1:T},y_{1:T})\mathrm{d}x_{1:T}$ (4)
$\displaystyle=\sum_{t=1}^{T}\log p_{\theta}(y_{t}|y_{1:t-1})\,,$ (5)

where

$p_{\theta}(y_{t}|y_{1:t-1})=\int
g_{\theta}(y_{t}|x_{t})p_{\theta}(x_{t}|y_{1:t-1})\mathrm{d}x_{t}\,.$ (6)

Analytical expression for $p_{\theta}(y_{1:T})$ rarely exists, with a notable
exception of linear Gaussian SSMs. Here, the Kalman filter/smoother provides
the solution for the state inference tasks, and also for computing the
marginal data log-likelihood.

There are many approximate solutions for the state inference tasks in
nonlinear and/or non-Gaussian SSMs, with the PF being one popular method. PF
can also be used to estimate the marginal data log-likelihood, discussed in
the following subsections.

###  I-B Particle Filters

Particle filters are numeric approximations of the recursive Bayesian
estimation problem involving arbitrary density functions [3]. A pdf $p$ is
represented by $x^{i}$ particles and their weights $w^{i}$ in the form of a
Dirac sum

$\hat{p}(x;\chi,\omega)=\sum\limits_{i=1}^{N}w^{i}\delta(x-x^{i})\,,$ (7)

where $\chi\coloneqq\\{x^{i}\\}_{i=1}^{N}$ and
$\omega\coloneqq\\{w^{i}\\}_{i=1}^{N}$ are the parameters of the approximation
$\hat{p}$.

Particles are drawn from the proposal distribution $q$ and then weighted using
the measurement likelihood and the transitional density, producing the Dirac
sum representation of the posterior distribution. At timestep $t$, the
particle approximation of the posterior reads

$p(x_{t}|y_{1:t})\approx\sum_{i=1}^{N}w_{t}^{i}\delta(x_{t}-x_{t}^{i})\,,$ (8)

where the posterior weights are

$w_{t}^{i}=w_{t-1}^{i}\frac{p(y_{t}|x_{t}^{i})p(x_{t}^{i}|x_{t-1}^{i})}{q(x_{t}^{i}|x_{t-1}^{i},y_{t})}\,.$
(9)

In its most basic form, a PF uses the transitional density for the proposal
distribution and the likelihood values for the weights. A particle filter to
work in practice needs at least an adequate sampling method and a resampling
strategy to evade particle depletion [5] to concentrate particles on regions
with high probability mass. One common resampling scheme is multinomial
resampling, where particles are resampled from a categorical distribution
composed of the normalized particle weights. For other standard resampling
schemes, see [6].

Besides posterior inference, PFs can also be used as unbiased estimators of
the marginal data likelihood if the resampling scheme is unbiased (see proofs
in [7, 8]):

$\hat{p}_{\theta}(y_{1:T})=\prod_{t=1}^{T}\sum_{i=1}^{N}w_{t}^{i},$ (10)
$\mathbb{E}\left[\hat{p}_{\theta}(y_{1:T})\right]=p_{\theta}(y_{1:T})\,,$ (11)

where the expectation is taken w.r.t. the joint pdf of all the random
variables in the particle filtering process.

###  I-C Maximum Likelihood Estimation with Particle Filters

Maximum likelihood estimation (MLE) in probabilistic models can be formulated
as

$\hat{\theta}=\arg\underset{\theta}{\max}\ p_{\theta}(y_{1:T})\,.$ (12)

Unfortunately for general SSMs, $\hat{\theta}$ cannot be obtained
analytically, as the marginal data likelihood $p_{\theta}(y_{1:T})$ is
intractable, with the notable exception of linear Gaussian SSMs.

A common approach is instead the use of a variational lower bound to the
marginal data log-likelihood, serving as a surrogate objective for
optimization. One typical choice is the evidence lower bound (ELBO) [9],
defined as

$\mathcal{L}(y,\theta,\phi)=\mathbb{E}_{q_{\phi}(x|y)}\left[\log\frac{p_{\theta}(x,y)}{q_{\phi}(x|y)}\right]\leq\log
p_{\theta}(y)\,.$ (13)

The bound is tight if the variational posterior $q_{\phi}(x|y)$ is the true
posterior $p_{\theta}(x|y)$.

In time series models, the ELBO can be loose, as it does not take into
consideration the time dependency between states and measurements, such as
defined by an SSM. The concurrent works of [7, 8, 10] propose to use the
marginal data log-likelihood estimator of a PF as an objective for MLE
instead:

$\mathcal{L}(y_{1:T},\theta,\phi)=\mathbb{E}\left[\log\prod_{t=1}^{T}\sum_{i=1}^{N}w_{t}^{i}\right]\approx\frac{1}{B}\sum_{i=1}^{B}\log\prod_{t=1}^{T}\sum_{i=1}^{N}w_{t}^{i},$
(14)

where $\theta$ parameterizes the SSM, and $\phi$ parameterizes the proposal
distribution. Often, $\theta=\phi$. This bound can potentially be much tighter
than ELBO by increasing the number of particles. The ELBO used in [9] can be
extended into a family of ELBOs, where (14) is also an ELBO. To avoid
notational clutter, we will refer to (14) as ELBO from now on.

For optimization via gradient ascent, $\nabla_{\theta,\phi}\mathcal{L}(\cdot)$
is required. The analytic form of the expectation in (14) does not exist; only
its Monte Carlo estimator can be obtained. The expectation is taken w.r.t. the
joint distribution of all the random variables in the sampling and resampling
process, which may depend on the parameters ($\theta,\phi$). In most cases,
the reparameterization trick [9] can be used to move the gradient computation
inside the expectation. This requires the ability to reparameterize the
proposal distribution, and the resampling distribution. The problems arise
with the latter, as traditional resampling schemes, like multinomial
resampling, are discontinuous w.r.t. the model parameters. This means that
small changes in the model parameters can lead to abrupt changes in the
resampling outcome and thus in the Monte Carlo estimate (14), illustrated in
Figure 2. For instance, a particle with a low weight may not be resampled
initially, but slight adjustments in the model parameters could cause its
weight to increase, resulting in its occasional resampling. Ignoring this
effect leads to a high-variance Monte Carlo estimate of the gradient, and
prevents backpropagation through time.

Figure 2: Illustration of the marginal data log-likelihood estimation of a
standard PF with multinomial resampling (left), and a differentiable PF with
optimal placement resampling (right) for a 1-dimensional linear Gaussian SSM,
described in Section IV-A.

##  II Related Work and Contributions

Performing joint model and proposal learning in directed probabilistic models
using stochastic gradient ascent goes back to [9]. Their approach was further
extended by importance sampling principles in [11]. Although these methods
also work on structured time series models, the authors in [7, 8, 10]
concurrently showed that the marginal likelihood estimator of a PF provides a
tighter lower bound for such models. These works also noted the fact that the
resampling step is nondifferentiable but did not propose a solution, instead,
they used a biased gradient estimator. The nondifferentiable problem also came
up in [12, 13], where the aim was to learn motion and measurement models for
robotic localization, parameterized by neural networks and using a PF-based
supervised loss. Since then, multiple differentiable PF techniques have been
proposed [12, 14, 15, 16, 17, 18], but most notably [19], which was the first
fully differentiable particle filter. A recent overview of differentiable
particle filters is presented in [20], while [21] provides guidelines on
training them.

In this paper, we propose a differentiable particle filter structure. Our work
is similar to [22], where they also make the likelihood estimation of a PF
continuous by sampling from the cdf. Contrary to their work, we use
deterministic sampling, which satisfies a given optimality criteria. We
provide empirical results in learning system and proposal parameters on toy
examples and also on real-world data.

##  III Deterministic Sampling

Figure 3: Illustration of optimal placement resampling. The goal is to perform
resampling, i.e., create an unweighted particle set from a weighted particle
set (7) by moving particles to positions with high probability mass. The true
pdf (black) and true cdf (dashed black) do not change during resampling, only
their representation by particles. The two rows show two consecutive
timesteps. Left: weighted particle set before resampling, with stem length
proportional to particle weight. Middle: resampled particle set by optimal
placement resampling. Right: empirical cdfs. In the second row, the resampled
particles in the previous timestep are weighted according to (9), and the
resampling is performed again. Due to the particles’ finite representational
power, their histograms slightly differ before and after resampling.  Figure
4: Illustration of empirical pdf (left) and cdf (right) based on the optimal
particle locations. Black line: analytical pdf and cdf. Blue line: empirical
pdf and cdf. Blue triangles: particle locations. Red line: traditional step-
wise empirical cdf, which we replaced by the blue line approximation given in
(23).

Resampling in PFs prevents particle degeneration by eliminating low-weight
particles and moving particles to highly-weighted particle positions. This
procedure essentially replaces the weighted Dirac mixture representation of
the posterior with an unweighted Dirac mixture representation. It is important
to mention that the underlying posterior pdf does not change, only its
representation by the particles. Standard resampling schemes have two main
properties: they are stochastic, and they do not introduce new particle
positions.

Schrempf et al. [23] proposed a deterministic alternative for approximating a
pdf by unweighted particles. Their approach builds on the integral quadratic
distance

$D=\int_{-\infty}^{\infty}\left(F(x)-G(x)\right)^{2}\mathop{}\\!\mathrm{d}{x}\,,$
(15)

where $F$ and $G$ are the cumulative density functions of the distributions
involved. If a distribution is approximated by particles in the form (7), the
corresponding empirical cdf is

$\displaystyle\hat{F}(x;\chi,\omega)$
$\displaystyle=\int_{-\infty}^{x}\hat{f}(s;\chi,\omega)\mathop{}\\!\mathrm{d}{s}=\sum_{i=1}^{N}w^{i}H(x-x^{i})\,,$
(16)

where $H$ is the Heaviside function defined as

$H(x)=\left\\{\begin{array}[]{cl}0\,,&x<0\\\ \nicefrac{{1}}{{2}}\,,&x=0\\\
1\,,&x>0\end{array}\right.\,.$ (17)

The approach proposed in [23] is the following. Let the particle approximation
use equal weights, i.e., $w^{i}=1/N$. Inserting an exact cdf $F$ and its
particle approximation $\hat{F}$ from (16) into the distance function (15) and
minimizing it in terms of the $x^{i}$ particle positions one finds that

$F(x^{i})=\frac{2i-1}{2N}\,\,(i=1\dots N)\,,$ (18)

which is an implicit relation that gives the optimal $x^{i}$ positions of the
particles to represent the exact distribution $F$ while minimizing the
integral squared distance (15) between $F$ and $\hat{F}$.

To compute the $x^{i}$ positions, the cdf must be inverted, which is not
generally feasible. For the uniform distribution on the interval [0,1], the
solution is

$x^{i}=\frac{2i-1}{2N}\,.$ (19)

For an arbitrary pdf $f(x)$, Schrempf et al. introduced a homotopy progression
from an initial pdf $g(x)$ to $f(x)$, which generates a flow of the particles
to occupy the optimal positions.

###  III-A Optimal Placement Resampling

Instead of using a flow, we propose to place the particles directly to the
optimal positions according to (18). To this end, we create a cdf
approximation constructed from particle weights and positions using easily
invertible functions.

Consider a weighted particle representation of a distribution with particle
positions $x^{i}$ and weights $w^{i}$. The pdf will be approximated, instead
of a Dirac sum, with a weighted Heaviside sum in the form

$f(x)\\!=\\!\left\\{\\!\\!\begin{array}[]{ll}\\!\dfrac{w^{1}}{2}\exp(x-x^{1})\,,&x\leq
x^{1}\\\\[5.0pt]
\\!\varpi^{1}+\sum_{i=2}^{N-1}(\varpi^{i}\\!-\\!\varpi^{i-1})H(x\\!-\\!x^{i}),\\!\\!\\!&x^{1}\\!<\\!x\\!<\\!x^{N}\\\\[1.99997pt]
\\!\dfrac{w^{N}}{2}\exp(-x+x^{N}),&x^{N}\leq x\end{array}\right.$ (20)

where

$\varpi^{i}=\frac{1}{2}\frac{w^{i}+w^{i-1}}{x^{i}-x^{i-1}}\,,$ (21)

which is also the steepness of the ramp function in the cdf.

The ramp function is defined as

$R(x)=\int_{-\infty}^{x}H(s)\mathop{}\\!\mathrm{d}{s}\,,$ (22)

which allows writing

$F(x)=\left\\{\begin{array}[]{ll}\dfrac{w^{1}}{2}\exp(x-x^{1})\,,&x\leq
x^{1}\\\\[8.99994pt] \begin{aligned} \varpi^{1}&R(x-x^{1})+\\\
&\sum_{i=2}^{N-1}(\varpi^{i}\\!-\\!\varpi^{i-1})R(x-x^{i})\,,\end{aligned}&x^{1}<x<x^{N}\\\\[8.99994pt]
1-\frac{w^{N}}{2}\exp(-x+x^{N})\,,&x^{N}\leq x\end{array}\right.$ (23)

The exponential leading and trailing parts in the piecewise definitions (20)
and (23) ensure that these functions can be evaluated at any point and insert
a minimal amount of arbitrary probability mass to regions that are not covered
by particles (Fig. 4).

To be able to use the empirical cdf for computing the optimal particle
positions in (18), the inverse cdf is needed, which consists of linear and
logarithmic terms:

$F^{-1}(w)=\left\\{\begin{array}[]{ll}x^{1}+\log(\frac{2w}{w^{1}})\,,&w\leq
w^{1}\\\\[5.0pt]
\sum_{i=2}^{N-1}R(w-w^{i})/\varpi^{i}\,,&w^{1}<w<x^{N}\\\\[5.0pt]
x^{N}+\log(\frac{w^{N}}{2-2w})\,,&w^{N}\leq w\end{array}\right.$ (24)

We will use the ideas presented above. In the PF update step, no stochastic
resampling is performed after the weighted particle set is obtained. Instead,
each particle will move deterministically to the optimal position according to
(18) using the inverse cdf constructed as in (24). Note that the particles
have to be sorted in an increasing order. We will call this scheme optimal
placement resampling (OPR). OPR obtains optimally distributed samples from the
weighted empirical posterior distribution, thus not having duplicate
particles, maintaining diversity in the particle set. However, high-
probability regions in the posterior still concentrate more particles. Figure
3. illustrates the process.

##  IV Empirical Study

###  IV-A Linear Gaussian State-Space Model

For our first study, we examine a one-dimensional linear Gaussian state-space
model (LGSSM), defined as

$\displaystyle x_{t}$ $\displaystyle=\alpha x_{t-1}+v_{t}$ (25) $\displaystyle
y_{t}$ $\displaystyle=\gamma x_{t}+e_{t}\,,$ (26)

with $v_{t}\sim\mathcal{N}(0,\sigma_{x}^{2})$,
$e_{t}\sim\mathcal{N}(0,\sigma_{y}^{2})$ and
$x_{1}\sim\mathcal{N}(0,\sigma_{x}^{2})$. An analytic expression of
$p_{\theta}(y_{1:T})$ exists for this simple model, and the true marginal data
log-likelihood can be computed using the Kalman filter. We generated a
measurement dataset of $T=100$ steps, using
$\theta^{*}=[\alpha^{*},\gamma^{*}]=[0.5,1],\sigma_{x}^{2}=0.3,\sigma_{y}^{2}=0.1$.
We then tried to learn the maximum likelihood estimation of the parameters
$\theta=[\alpha,\gamma]$ using stochastic gradient ascent on the approximation
of the ELBO from a PF with multinomial resampling (PF-MR) and optimal
placement resampling (PF-OPR). The particle filters were run with $N=50$
particles and in $B=50$ batches to evaluate the ELBO (14). The Adam optimizer
[24] was used with a learning rate of 0.01 for 200 epochs, and
$\theta_{0}=[1,1.5]$ as initial parameters. Similarly to the result in [19],
we obtain that, in this simplistic case, the nondifferentiability of PF-MR
does not cause problems in learning, and the results compared to PF-OPR are
similar. The ELBO from both PF methods has a relative error of 1.5$\%$ to the
true marginal data log-likelihood.

###  IV-B Proposal Distribution Learning

In particle filters, the proposal distribution is a key choice. Based on the
experiments in [7], we learn a time-varying proposal distribution in the form
of

$r_{\lambda}(x_{t}|x_{t-1})=\mathcal{N}(\mu_{t}+\beta_{t}\alpha
x_{t-1},\sigma_{t}^{2})$ (27)

for the LGSSM defined in (25)-(26), with
$\lambda=\\{\mu_{t},\beta_{t},\sigma^{2}_{t}\\}_{t=1}^{T}$. As the parameters
are time-varying in this case, efficient learning requires the ability to
backpropagate through time.

Using $\alpha=0.42,\gamma=1,\sigma^{2}_{x}=1,\sigma^{2}_{y}=0.1$, we performed
the maximum likelihood learning of the proposal parameters using the ELBO
estimate from PF-MR and PF-OPR. The Adam optimizer was chosen with a learning
rate of 0.1, with initial parameters $\mu_{t}=0,\beta_{t}=1,\log\sigma_{t}=0$
for $t=1...T$. The particle filters were run in $B=50$ batches with $N=100$
particles. Figure 5 shows the ELBO estimate of PF-MR in orange and PF-OPR in
blue, along with the maximal component of the gradient of the ELBO w.r.t
$\lambda$. It is clearly apparent that multinomial resampling performs worse
for this task, which is likely due to its inability to perform backpropagation
through time as its resampling scheme is nondifferentiable. The execution time
per epoch was 83.4 ms for PF-MR and 113.7 ms for PF-OPR. We note that the most
time-consuming process in PF-OPR is sorting particles required for
constructing the cdf. However, the OPR algorithm still runs in
$\mathcal{O}(N)$ time complexity.

Figure 5: Proposal distribution learning for an LGSSM using PF-MR (orange) and
PF-OPR (blue). Higher ELBO is better.

###  IV-C Stochastic Volatility Model

The stochastic volatility model [25] is a common choice for modeling financial
processes, like price changes in the stock market. The model is defined as

$\displaystyle x_{t}$ $\displaystyle=\mu+\phi(x_{t-1}-\mu)+v_{t}$ (28)
$\displaystyle y_{t}$ $\displaystyle=\exp\left(\frac{x_{t}}{2}\right)e_{t}\,,$
(29)

with $v_{t}\sim\mathcal{N}(0,\sigma_{x}^{2})$,
$e_{t}\sim\mathcal{N}(0,\sigma_{y}^{2})$,
$x_{1}\sim\mathcal{N}\left(\mu,\frac{\sigma_{x}^{2}}{1-\phi^{2}}\right)$. The
parameter set of interest is $\theta=[\mu,\phi,\sigma_{x},\sigma_{y}]$. The
observations are the log-returns, which are defined as

$y_{t}=100\log\left(\frac{s_{t}}{s_{t-1}}\right)\,,$ (30)

where $s_{t}$ is the price at $t$. Our observation dataset is constructed from
daily EUR/HUF exchange rates of the European Central Bank from 01.01.2017 to
01.01.2023. As the analytic expression for $p_{\theta}(y_{1:T})$ is not
tractable, we used the approximates produced by PF-MR and PF-OPR. The training
results with $N=50$ particles, $B=50$ batches, and using the Adam optimizer
with a learning rate of $0.01$ is shown in Figure 6. The ELBO for PF-OPR is
-634.9, and the ELBO for PF-MR is -640.0, meaning that PF-OPR is capable of
providing a tighter and therefore better ELBO estimate.

Figure 6: The estimated ELBO and its maximum gradient component w.r.t.
$\theta$ for a parameter inference task in a stochastic volatility model.
Blue: PF-OPR, orange: PF-MR. Higher ELBO is better.

##  V Conclusions and Future Work

We introduced optimal placement resampling for particle filters as a way of
solving the nondifferentiability problem of standard resampling algorithms.
Our method uses optimally placed samples from a hand-crafted empirical
cumulative distribution function based on the prior weighted particle set. We
empirically showed that optimal placement resampling outperforms the
nondifferentiable multinomial resampling in model and proposal learning tasks.

In the current form, the optimal placement resampling works only in one
dimension due to the fact that it relies on the cdf. The cdf is not a uniquely
defined function since the relation in $F(x)=P(X\leq x)$ is arbitrary. It does
not cause a problem in one dimension because the cdfs $P(X\leq x)$ and
$P(X\geq x)$ lead to the same results. Starting from dimension two, this
invariance no longer holds, and alternatives to the conventional cdf have to
be used [26]. A future task is to design easily computable optimal placement
strategies in multiple dimensions using either an alternative cdf or another
placement strategy.

## Acknowledgment

This research was supported by the European Union within the framework of the
National Laboratory for Autonomous Systems. (RRF-2.3.1-21-2022-00002).

DC was supported by ÚNKP-23-3-I-BME-146 New National Excellence Program of the
Ministry for Culture and Innovation.

The authors would also like to thank Adrien Corenflos and James Thornton for
making their code in [19] publicly available.

## References

  * [1] P. M. Mazhdrakov, D. Benov, and P. N. Valkanov, _The Monte Carlo Method_. ACMO Academic Press, Aug. 2018.
  * [2] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of State Calculations by Fast Computing Machines,” _The Journal of Chemical Physics_ , vol. 21, no. 6, pp. 1087–1092, Jun. 1953.
  * [3] N. Gordon, D. Salmond, and A. Smith, “Novel approach to nonlinear/non-Gaussian Bayesian state estimation,” _IEE Proceedings F Radar and Signal Processing_ , vol. 140, no. 2, p. 107, 1993.
  * [4] F. Daum and J. Huang, “Curse of dimensionality and particle filters,” in _2003 IEEE Aerospace Conference Proceedings (Cat. No.03TH8652)_ , vol. 4, Mar. 2003, pp. 4_1979–4_1993, iSSN: 1095-323X.
  * [5] M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” _IEEE Transactions on Signal Processing_ , vol. 50, no. 2, pp. 174–188, Feb. 2002.
  * [6] T. Li, M. Bolic, and P. M. Djuric, “Resampling Methods for Particle Filtering: Classification, implementation, and strategies,” _IEEE Signal Processing Magazine_ , vol. 32, no. 3, pp. 70–86, May 2015.
  * [7] C. Naesseth, S. Linderman, R. Ranganath, and D. Blei, “Variational Sequential Monte Carlo,” in _Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics_. PMLR, Mar. 2018, pp. 968–977, iSSN: 2640-3498.
  * [8] C. J. Maddison, J. Lawson, G. Tucker, N. Heess, M. Norouzi, A. Mnih, A. Doucet, and Y. Teh, “Filtering Variational Objectives,” in _Advances in Neural Information Processing Systems_ , vol. 30. Curran Associates, Inc., 2017.
  * [9] D. P. Kingma and M. Welling, “Auto-Encoding Variational Bayes,” _CoRR_ , Dec. 2013.
  * [10] T. Le, M. Igl, T. Rainforth, T. Jin, and F. Wood, “Auto-encoding sequential Monte Carlo,” _International Conference on Learning Representations (ICLR)_ , 2018.
  * [11] Y. Burda, R. Grosse, and R. Salakhutdinov, “Importance Weighted Autoencoders,” _International Conference on Learning Representations_ , 2015.
  * [12] P. Karkus, D. Hsu, and W. S. Lee, “Particle Filter Networks with Application to Visual Localization,” in _Proceedings of The 2nd Conference on Robot Learning_. PMLR, Oct. 2018, pp. 169–178, iSSN: 2640-3498.
  * [13] R. Jonschkowski, D. Rastogi, and O. Brock, “Differentiable Particle Filters: End-to-End Learning with Algorithmic Priors,” in _Robotics: Science and Systems XIV_. Robotics: Science and Systems Foundation, Jun. 2018.
  * [14] A. Ścibior and F. Wood, “Differentiable Particle Filtering without Modifying the Forward Pass,” Oct. 2021, arXiv:2106.10314 [cs, stat].
  * [15] A. Younis and E. Sudderth, “Differentiable and Stable Long-Range Tracking of Multiple Posterior Modes,” _Advances in Neural Information Processing Systems_ , vol. 36, Dec. 2023.
  * [16] A. Singh, O. Makhlouf, M. Igl, J. Messias, A. Doucet, and S. Whiteson, “Particle-Based Score Estimation for State Space Model Learning in Autonomous Driving,” in _Proceedings of The 6th Conference on Robot Learning_. PMLR, Mar. 2023, pp. 1168–1177, iSSN: 2640-3498.
  * [17] M. Zhu, K. Murphy, and R. Jonschkowski, “Towards Differentiable Resampling,” Apr. 2020, arXiv:2004.11938 [cs, stat].
  * [18] C. Rosato, L. Devlin, V. Beraud, P. Horridge, T. B. Schön, and S. Maskell, “Efficient Learning of the Parameters of Non-Linear Models Using Differentiable Resampling in Particle Filters,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 3676–3692, 2022.
  * [19] A. Corenflos, J. Thornton, G. Deligiannidis, and A. Doucet, “Differentiable Particle Filtering via Entropy-Regularized Optimal Transport,” in _Proceedings of the 38th International Conference on Machine Learning_. PMLR, Jul. 2021, pp. 2100–2111, iSSN: 2640-3498.
  * [20] X. Chen and Y. Li, “An overview of differentiable particle filters for data-adaptive sequential Bayesian inference,” _Foundations of Data Science_ , Dec. 2023.
  * [21] A. Kloss, G. Martius, and J. Bohg, “How to train your differentiable filter,” _Autonomous Robots_ , vol. 45, no. 4, pp. 561–578, May 2021.
  * [22] S. Malik and M. K. Pitt, “Particle filters for continuous likelihood evaluation and maximisation,” _Journal of Econometrics_ , vol. 165, no. 2, pp. 190–209, Dec. 2011.
  * [23] O. C. Schrempf, D. Brunn, and U. D. Hanebeck, “Density Approximation Based on Dirac Mixtures with Regard to Nonlinear Estimation and Filtering,” in _Proceedings of the 45th IEEE Conference on Decision and Control_. San Diego, CA, USA: IEEE, 2006, pp. 1709–1714.
  * [24] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,” Jan. 2017, arXiv:1412.6980 [cs].
  * [25] J. Hull and A. White, “The Pricing of Options on Assets with Stochastic Volatilities,” _The Journal of Finance_ , vol. 42, no. 2, pp. 281–300, 1987.
  * [26] U. D. Hanebeck and V. Klumpp, “Localized Cumulative Distributions and a multivariate generalization of the Cramér-von Mises distance,” in _2008 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems_ , Aug. 2008, pp. 33–39.