# Differentiable Particle Filtering via Entropy-Regularized Optimal Transport

Adrien Corenflos  James Thornton  George Deligiannidis  Arnaud Doucet

###### Abstract

Particle Filtering (PF) methods are an established class of procedures for
performing inference in non-linear state-space models. Resampling is a key
ingredient of PF, necessary to obtain low variance likelihood and states
estimates. However, traditional resampling methods result in PF-based loss
functions being non-differentiable with respect to model and PF parameters. In
a variational inference context, resampling also yields high variance gradient
estimates of the PF-based evidence lower bound. By leveraging optimal
transport ideas, we introduce a principled differentiable particle filter and
provide convergence results. We demonstrate this novel method on a variety of
applications.

Monte Carlo, Particle Filtering, State-Space Models, Optimal Transport

##  1 Introduction

In this section we provide a brief introduction to state-space models (SSMs)
and PF methods. We then illustrate one of the well-known limitations of PF
(Kantas et al., 2015): resampling steps are required in order to compute low-
variance estimates, but these estimates are not differentiable w.r.t. to model
and PF parameters. This hinders end-to-end training. We discuss recent
approaches to address this problem in econometrics, statistics and machine
learning (ML), outline their limitations and our contributions.

###  1.1 State-Space Models

SSMs are an expressive class of sequential models, used in numerous scientific
domains including econometrics, ecology, ML and robotics; see e.g. (Chopin &
Papaspiliopoulos, 2020; Douc et al., 2014; Doucet & Lee, 2018; Kitagawa &
Gersch, 1996; Lindsten & Schön, 2013; Thrun et al., 2005). SSM may be
characterized by a latent $\mathcal{X}$-valued Markov process $(X_{t})_{t\geq
1}$ and $\mathcal{Y}$-valued observations $(Y_{t})_{t\geq 1}$ satisfying
$X_{1}\sim\mu_{\theta}(\cdot)$ and for $t\geq 1$

$X_{t+1}|\\{X_{t}=x\\}\sim f_{\theta}(\cdot|x),~{}~{}Y_{t}|\\{X_{t}=x\\}\sim
g_{\theta}(\cdot|x),$ (1)

where $\theta\in\Theta$ is a parameter of interest. Given observations
$(y_{t})_{t\geq 1}$ and parameter values $\theta$, one may perform state
inference at time $t$ by computing the posterior of $X_{t}$ given
$y_{1:t}\coloneqq(y_{1},...,y_{t})$ where

$\displaystyle p_{\theta}(x_{t}|y_{1:t-1})$ $\displaystyle=\int
f_{\theta}(x_{t}|x_{t-1})p_{\theta}(x_{t-1}|y_{1:t-1})\textrm{d}x_{t-1},$
$\displaystyle p_{\theta}(x_{t}|y_{1:t})$
$\displaystyle=\frac{g_{\theta}(y_{t}|x_{t})p_{\theta}(x_{t}|y_{1:t-1})}{\int
g_{\theta}(y_{t}|x_{t})p_{\theta}(x_{t}|y_{1:t-1})\textrm{d}x_{t}},$

with $p_{\theta}(x_{1}|y_{0}):=\mu_{\theta}(x_{1})$.

The log-likelihood $\ell(\theta)=\log p_{\theta}(y_{1:T})$ is then given by

$\ell(\theta)=\sum_{t=1}^{T}\log p_{\theta}(y_{t}|y_{1:t-1}),$

with $p_{\theta}(y_{1}|y_{0})\coloneqq\int
g_{\theta}(y_{1}|x_{1})\mu_{\theta}(x_{1})\textrm{d}x_{1}$ and for $t\geq 2$

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

The posteriors $p_{\theta}(x_{t}|y_{1:t})$ and log-likelihood
$p_{\theta}(y_{1:T})$ are available analytically for only a very restricted
class of SSM such as linear Gaussian models. For non-linear SSM, PF provides
approximations of such quantities.

###  1.2 Particle Filtering

PF are Monte Carlo methods entailing the propagation of $N$ weighted particles
$(w^{i}_{t},X^{i}_{t})_{i\in[N]}$, here $[N]\coloneqq\\{1,...,N\\}$, over time
to approximate the filtering distributions $p_{\theta}(x_{t}|y_{1:t})$ and
log-likelihood $\ell(\theta)$. Here $X^{i}_{t}\in\mathcal{X}$ denotes the
value of the $i^{\text{th}}$ particle at time $t$ and
$\mathbf{w}_{t}\coloneqq(w_{t}^{1},...,w_{t}^{N})$ are weights satisfying
$w^{i}_{t}\geq 0,~{}\sum_{i=1}^{N}w^{i}_{t}=1$. Unlike variational methods, PF
methods provide consistent approximations under weak assumptions as
$N\rightarrow\infty$ (Del Moral, 2004). In the general setting, particles are
sampled according to proposal distributions $q_{\phi}(x_{1}|y_{1})$ at time
$t=1$ and $q_{\phi}(x_{t}|x_{t-1},y_{t})$ at time $t\geq 2$ prior to weighting
and resampling. One often chooses $\theta=\phi$ but this is not necessarily
the case (Le et al., 2018; Maddison et al., 2017; Naesseth et al., 2018).

Algorithm 1 Standard Particle Filter

1: Sample
$X_{1}^{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}q_{\phi}(\cdot|y_{1})$
for $i\in[N]$

2: Compute
$\omega_{1}^{i}=\frac{p_{\theta}(X^{i}_{1},y_{1})}{q_{\phi}(X^{i}_{1}|y_{1})}$
for $i\in[N]$

3: $\hat{\ell}(\theta)\leftarrow\frac{1}{N}\sum_{i=1}^{N}\omega_{1}^{i}$

4: for $t=2,...,T$ do

5: Normalize weights $w^{i}_{t-1}\propto\omega^{i}_{t-1}$,
$\sum_{i=1}^{N}w^{i}_{t-1}=1$

6: Resample
$\tilde{X}^{i}_{t-1}\sim\sum_{i=1}^{N}w^{i}_{t-1}\delta_{X^{i}_{t-1}}$ for
$i\in[N]$

7: Sample $X_{t}^{i}\sim q_{\phi}(\cdot|\tilde{X}^{i}_{t-1},y_{t})$ for
$i\in[N]$

8: Compute
$\omega^{i}_{t}=\frac{p_{\theta}(X_{t}^{i},y_{t}|\tilde{X}^{i}_{t-1})}{q_{\phi}(X_{t}^{i}|\tilde{X}^{i}_{t-1},y_{t})}$
for $i\in[N]$

9: Compute
$\hat{p}_{\theta}(y_{t}|y_{1:t-1})=\frac{1}{N}\sum_{i=1}^{N}\omega^{i}_{t}$

10:
$\hat{\ell}(\theta)\leftarrow\hat{\ell}(\theta)+\log{\hat{p}}_{\theta}(y_{t}|y_{1:t-1})$

11: end for

12: Return: log-likelihood estimate
$\hat{\ell}(\theta)=\log\hat{p}_{\theta}(y_{1:T})$

A generic PF is described in Algorithm 1 where
$p_{\theta}(x_{1},y_{1})\coloneqq\mu_{\theta}(x_{1})g_{\theta}(y_{1}|x_{1})$
and $p_{\theta}(x_{t},y_{t}|x_{t-1})\coloneqq
f_{\theta}(x_{t}|x_{t-1})g_{\theta}(y_{t}|x_{t})$. Resampling is performed in
step 6 of Algorithm 1; it ensures particles with high weights are replicated
and those with low weights are discarded, allowing one to focus computational
efforts on ‘promising’ regions. The scheme used in Algorithm 1 is known as
multinomial resampling and is unbiased (as are other traditional schemes such
as stratified and systematic (Chopin & Papaspiliopoulos, 2020)), i.e.

$\mathbb{E}\left[\tfrac{1}{N}{\textstyle\sum}_{i=1}^{N}\psi(\tilde{X}^{i}_{t})\right]=\mathbb{E}\left[{\textstyle\sum}_{i=1}^{N}w_{t}^{i}\psi(X_{t}^{i})\right],$
(2)

for any $\psi:\mathcal{X}\rightarrow\mathbb{R}$. This property guarantees
$\exp\small(\hat{\ell}(\theta)\small)$ is an unbiased estimate of the
likelihood $\exp(\ell(\theta))$ for any $N$.

Henceforth, let $\mathcal{X}=\mathbb{R}^{d_{x}}$,
$\theta\in\Theta=\mathbb{R}^{d_{\theta}}$ and
$\phi\in\Phi=\mathbb{R}^{d_{\phi}}$. We assume here that
$\theta\mapsto\mu_{\theta}(x)$, $\theta\mapsto f_{\theta}(x^{\prime}|x)$ and
$\theta\mapsto g_{\theta}(y_{t}|x)$ are differentiable for all $x,x^{\prime}$
and $t\in[T]$ and $\theta\mapsto{\ell}(\theta)$ is differentiable. These
assumptions are satisfied by a large class of SSMs. We also assume that we can
use the reparameterization trick (Kingma & Welling, 2014) to sample the
particles; i.e. we have $\Gamma_{\phi}(y_{1},U)\sim
q_{\phi}(x_{1}|y_{1}),\Psi_{\phi}(y_{t},x_{t-1},U)\sim
q_{\phi}(x_{t}|x_{t-1},y_{t})$ for some mappings $\Gamma_{\phi},\Psi_{\phi}$
differentiable w.r.t. $\phi$ and $U\sim\lambda$, $\lambda$ being independent
of $\phi$.

###  1.3 Related Work and Contributions

Let $\bf{U}$ be the set of all random variables used to sample and resample
the particles. The distribution of $\bf{U}$ is $(\theta,\phi)$-independent as
we use the reparameterization trick111For example, multinomial resampling
relies on $N$ uniform random variables.. However, even if we sample and fix
$\bf{U}=\bf{u}$, resampling involves sampling from an atomic distribution and
introduces discontinuities in the particles selected when $\theta,\phi$ vary.

(a) Kalman Filter

(b) Standard PF

(c) Differentiable PF

Figure 1: Left: Log-likelihood $\ell(\theta)$ and PF estimates
$\hat{\ell}(\theta;\phi,\bf{u})$ for linear Gaussian SSM, given in Section
5.1, with $d_{\theta}=2$ $d_{x}=2$, and $T=150,N=50$. Right:
$\nabla_{\theta}\ell(\theta)$ and
$\nabla_{\theta}\hat{\ell}(\theta;\phi,\bf{u})$.

For $d_{x}=1$, Malik & Pitt (2011) make
$\theta\mapsto\hat{\ell}(\theta;\phi,\mathbf{u})$ continuous w.r.t. $\theta$
by sorting the particles and then sampling from a smooth approximation of
their cumulative distribution function. For $d_{x}>1$, Lee (2008) proposes a
smoother but only piecewise continuous estimate. DeJong et al. (2013) returns
a differentiable log-likelihood estimate $\hat{\ell}(\theta;\phi,\mathbf{u})$
by using a marginal PF (Klaas et al., 2005), where importance sampling is
performed on a collapsed state-space. However, the standard marginal PF uses
the proposal
$q_{\phi}(x_{t}):=\sum_{i=1}^{N}w_{t-1}^{i}q_{\phi}(x_{t}|X^{i}_{t-1},y_{t})$
from which one cannot generally sample smoothly for arbitrary mixture
components. As a consequence they instead suggest using a simple Gaussian
distribution for $q_{\phi}(x_{t})$, which can lead to poor estimates for
multimodal posteriors. Moreover, in contrast to standard PF, this marginal PF
cannot be applied in scenarios where the transition density can only be
sampled from (e.g. using the reparameterization trick) but not evaluated
pointwise (Murray et al., 2013), as the importance weight would be
intractable. The implicit reparameterization method of Graves (2016) may be
used to obtain low variance gradient estimates with a mixture proposal. This
method however is only compatible for component distributions with tractable
conditional CDFs, such as Gaussian distributions. An alternative unbiased
estimate of the likelihood based on dynamic programming may also be obtained
(Finke et al., 2016; Aitchison, 2019). As emphasized by Aitchison (2019), this
estimate is differentiable. This approach is again limited however to a
restricted class of proposal distributions, such as an unweighted mixture
proposal, which may perform poorly for slow-mixing time-series.

In the context of robot localization, a modified resampling scheme has been
proposed in (Karkus et al., 2018; Ma et al., 2020a, b) referred to as ‘soft-
resampling’ (SPF). SPF has parameter $\alpha\in[0,1]$ where $\alpha=1$
corresponds to regular PF resampling and $\alpha=0$ is essentially sampling
particles uniformly at random. The resulting PF-net is said to be
differentiable but computes gradients that ignore the non-differentiable
component of the resampling step. Jonschkowski et al. (2018) proposed another
PF scheme which is said to be differentiable but simply ignores the non-
differentiable resampling terms and proposes new states based on the
observation and some neural network. This approach however does not propagate
gradients through time. Finally, Zhu et al. (2020) propose a differentiable
resampling scheme based on transformers but they report that the best results
are achieved when not backpropagating through it, due to exploding gradients.
Hence no fully differentiable PF is currently available in the literature
(Kloss et al., 2020).

PF methods have also been fruitfully exploited in Variational Inference (VI)
to estimate $\theta,\phi$ (Le et al., 2018; Maddison et al., 2017; Naesseth et
al., 2018). As
$\mathbb{E}_{\mathbf{U}}[\exp\big(\hat{\ell}(\theta;\phi,\mathbf{U})\big{missing})]=\exp(\ell(\theta))$
is an unbiased estimate of $\exp(\ell(\theta))$ for any $N,\phi$ for standard
PF, then one has indeed by Jensen’s inequality

$\ell^{\mathrm{ELBO}}(\theta,\phi)\coloneqq\mathbb{E}_{\mathbf{U}}[\hat{\ell}(\theta;\phi,\mathbf{U})]\leq\ell(\theta).$
(3)

The standard ELBO corresponds to $N=1$ and many variational families for
approximating $p_{\theta}(x_{1:T}|y_{1:T})$ have been proposed in this context
(Archer et al., 2015; Krishnan et al., 2017; Rangapuram et al., 2018). The
variational family induced by a PF differs significantly as
$\ell^{\mathrm{ELBO}}(\theta,\phi)\rightarrow\ell(\theta)$ as
$N\rightarrow\infty$ and thus yields a variational approximation converging to
$p_{\theta}(x_{1:T}|y_{1:T})$. This attractive property comes at a
computational cost; i.e. the PF approach trades off fidelity to the posterior
with computational complexity. While unbiased gradient estimates of the PF-
ELBO (3) can be computed, they suffer from high variance as the resampling
steps require having to use REINFORCE gradient estimates (Williams, 1992).
Consequently, Hirt & Dellaportas (2019); Le et al. (2018); Maddison et al.
(2017); Naesseth et al. (2018) use biased gradient estimates which ignore
these terms, yet report improvements as $N$ increases over standard VI
approaches and Importance Weighted Auto-Encoders (IWAE) (Burda et al., 2016).

Finally, if one is only interested in estimating $\theta$ (and not some
distinct $\phi$), then particle techniques approximating pointwise the score
vector $\nabla_{\theta}\ell(\theta)$ are also available (Poyiadjis et al.,
2011; Kantas et al., 2015).

The contributions of this paper are four-fold.

  * •

We propose the first fully Differentiable Particle Filter (DPF) which can use
general proposal distributions. DPF provides a differentiable estimate of
$\ell(\theta)$, see Figure 1-c, and more generally differentiable estimates of
PF-based losses. Empirically, in a VI context, DPF-ELBO gradient estimates
also exhibit much smaller variance than those of PF-ELBO.

  * •

We provide quantitative convergence results on the differentiable resampling
scheme and establish consistency results for DPF.

  * •

We show that existing techniques provide inconsistent gradient estimates and
that the non-vanishing bias can be very significant, leading practically to
unreliable parameter estimates.

  * •

We demonstrate that DPF empirically outperforms recent alternatives for end-
to-end parameter estimation on a variety of applications.

Proofs of results are given in the Supplementary Material.

##  2 Resampling via Optimal Transport

###  2.1 Optimal Transport and the Wasserstein Metric

Since Optimal Transport (OT) (Peyré & Cuturi, 2019; Villani, 2008) is a core
component of our scheme, the basics are presented here. Given two probability
measures $\alpha,\beta$ on $\mathcal{X}=\mathbb{R}^{d_{x}}$ the squared
2-Wasserstein metric between these measures is given by

$\displaystyle\mathcal{W}_{2}^{2}(\alpha,\beta)$
$\displaystyle=\min_{\mathcal{P}\in\mathcal{U}(\alpha,\beta)}\mathbb{E}_{(U,V)\sim\mathcal{P}}\big{[}||U-V||^{2}\big{]},$
(4)

where $\mathcal{U}(\alpha,\beta)$ the set of distributions on
$\mathcal{X}\times\mathcal{X}$ with marginals $\alpha$ and $\beta$, and the
minimizing argument of (4) is the OT plan denoted $\mathcal{P}^{\mathrm{OT}}$.
Any element $\mathcal{P}\in\mathcal{U}(\alpha,\beta)$ allows one to
“transport” $\alpha$ to $\beta$ (and vice-versa) i.e.

$\displaystyle\beta(\textrm{d}v)=\int\mathcal{P}(\textrm{d}u,\textrm{d}v)=\int\mathcal{P}(\textrm{d}v|u)\alpha(\textrm{d}u).$

For atomic probability measures $\alpha_{N}=\sum_{i=1}^{N}a_{i}\delta_{u_{i}}$
and $\beta_{N}=\sum_{j=1}^{N}b_{j}\delta_{v_{j}}$ with weights
$\mathbf{a}=(a_{i})_{i\in[N]}$, $\mathbf{b}=(b_{j})_{j\in[N]}$, and atoms
$\mathbf{u}=(u_{i})_{i\in[N]}$, $\mathbf{v}=(v_{j})_{j\in[N]}$, one can show
that

$\displaystyle\mathcal{W}_{2}^{2}(\alpha_{N},\beta_{N})=\min_{\mathbf{P}\in\mathcal{S}(\mathbf{a},\mathbf{b})}\textstyle{\sum}_{i=1}^{N}\sum_{j=1}^{N}c_{i,j}p_{i,j},$
(5)

where any $\mathcal{P}\in\mathcal{U}(\alpha_{N},\beta_{N})$ is of the form

$\mathcal{P}(\textrm{d}u,\textrm{d}v)={\textstyle\sum}_{i,j}~{}p_{i,j}\delta_{u_{i}}(\textrm{d}u)\delta_{v_{j}}(\textrm{d}v),$

$c_{i,j}=||u_{i}-v_{j}||^{2}$, $\mathbf{P}=(p_{i,j})_{i,j\in[N]}$ and
$\mathcal{S}(\mathbf{a},\mathbf{b})=\big{\\{}\mathbf{P}\in[0,1]^{N\times
N}:\sum_{j=1}^{N}p_{i,j}=a_{i},~{}\sum_{i=1}^{N}p_{i,j}=b_{j}\big{\\}}.$ In
such cases, one has

$\mathcal{P}(\textrm{d}v|u=u_{i})={\textstyle\sum}_{j}~{}a^{-1}_{i}p_{i,j}\delta_{v_{j}}(\textrm{d}v).$
(6)

The optimization problem (5) may be solved through linear programming. It is
also possible to exploit the dual formulation

$\mathcal{W}_{2}^{2}(\alpha_{N},\beta_{N})=\max_{\mathbf{f},\mathbf{g}\in\mathcal{R}(C)}\mathbf{a}^{t}\mathbf{f}+\mathbf{b}^{t}\mathbf{g},$
(7)

where $\mathbf{f}=(f_{i})$, $\mathbf{g}=(g_{i})$, $\mathbf{C}=(c_{i,j})$ and
$\mathcal{R}(\mathbf{C})=\\{\mathbf{f},\mathbf{g}\in\mathbb{R}^{N}|f_{i}+g_{j}\leq
c_{i,j},i,j\in[N]\\}$.

###  2.2 Ensemble Transform Resampling

The use of OT for resampling in PF has been pioneered by Reich (2013). Unlike
standard resampling schemes (Chopin & Papaspiliopoulos, 2020; Doucet & Lee,
2018), it relies not only on the particle weights but also on their locations.

At time $t$, after the sampling step (Step 7 in Algorithm 1),
$\alpha^{(t)}_{N}=\tfrac{1}{N}\sum_{i=1}^{N}\delta_{X_{t}^{i}}$ is a particle
approximation of $\alpha^{(t)}\coloneqq\int
q_{\phi}(x_{t}|x_{t-1},y_{t})p_{\theta}(x_{t-1}|y_{1:t-1})\textrm{d}x_{t-1}$
and $\beta^{(t)}_{N}=\sum w_{t}^{i}\delta_{X_{t}^{i}}$ is an approximation of
$\beta^{(t)}\coloneqq p_{\theta}(x_{t}|y_{1:t})$. Under mild regularity
conditions, the OT plan minimizing $\mathcal{W}_{2}(\alpha^{(t)},\beta^{(t)})$
is of the form
$\mathcal{P}^{\mathrm{OT}}(\textrm{d}x,\textrm{d}x^{\prime})=\alpha^{(t)}(\textrm{d}x)\delta_{\mathbf{T}^{(t)}(x)}(\textrm{d}x^{\prime})$
where $\mathbf{T}^{(t)}:\mathcal{X}\rightarrow\mathcal{X}$ is a deterministic
map; i.e if $X\sim\alpha^{(t)}$ then $\mathbf{T}^{(t)}(X)\sim\beta^{(t)}$. It
is shown in (Reich, 2013) that one can one approximate this transport map with
the ‘Ensemble Transform’ (ET) denoted $\mathbf{T}^{(t)}_{N}$. This is found by
solving the OT problem (5) between $\alpha^{(t)}_{N}$ and $\beta^{(t)}_{N}$
and taking an expectation w.r.t. (6), that is

$\tilde{X}^{i}_{t}=N{\textstyle\sum}_{k=1}^{N}~{}p^{\mathrm{OT}}_{i,k}X^{k}_{t}\coloneqq\mathbf{T}^{(t)}_{N}(X^{i}_{t}),$
(8)

where we slightly abuse notation as $\mathbf{T}^{(t)}_{N}$ is a function of
$X^{1:N}_{t}$. Reich (2013) uses this update instead of using
$\tilde{X}^{i}_{t}\sim\sum_{i=1}^{N}w^{i}_{t}\delta_{X^{i}_{t}}$. This is
justified by the fact that, as $N\rightarrow\infty$,
$\mathbf{T}^{(t)}_{N}(X^{i}_{t})\rightarrow\mathbf{T}^{(t)}(X^{i}_{t})$ in
some weak sense (Reich, 2013; Myers et al., 2021). Compared to standard
resampling schemes, the ET only satisfies (2) for affine functions $\psi$.

This OT approach to resampling involves solving the linear program (4) at cost
$O(N^{3}\log N)$ (Bertsimas & Tsitsiklis, 1997). This is not only
prohibitively expensive but moreover the resulting ET is not differentiable.
To address these problems, one may instead rely on entropy-regularized OT
(Cuturi, 2013).

##  3 Differentiable Resampling via Entropy-Regularized Optimal Transport

###  3.1 Entropy-Regularized Optimal Transport

Entropy-regularized OT may be used to compute a transport matrix that is
differentiable with respect to inputs and computationally cheaper than the
non-regularized version, i.e. we consider the following regularized version of
(5) for some $\epsilon>0$ (Cuturi, 2013; Peyré & Cuturi, 2019)

$\mathcal{W}^{2}_{2,\epsilon}(\alpha_{N},\beta_{N})=\min_{\mathbf{P}\in\mathcal{S}(\mathbf{a},\mathbf{b})}\sum_{i,j=1}^{N}p_{i,j}\Big{(}c_{i,j}+\epsilon\log\frac{p_{i,j}}{a_{i}b_{j}}\Big{)}.$
(9)

The function minimized in (9) is strictly convex and hence admits a unique
minimizing argument
$\mathbf{P}^{\mathrm{OT}}_{\epsilon}=(p^{\mathrm{OT}}_{\epsilon,i,j})$.
$\mathcal{W}^{2}_{2,\epsilon}(\alpha_{N},\beta_{N})$ can also be computed
using the regularized dual; i.e.
$\mathcal{W}^{2}_{2,\epsilon}(\alpha_{N},\beta_{N})=\max_{\mathbf{f},\mathbf{g}}\mathrm{DOT}_{\epsilon}(\mathbf{f},\mathbf{g})$
with

$\displaystyle\mathrm{DOT}_{\epsilon}(\mathbf{f},\mathbf{g})\coloneqq\mathbf{a}^{t}\mathbf{f}+\mathbf{b}^{t}\mathbf{g}-\epsilon\mathbf{a}^{t}\mathbf{M}\mathbf{b}$
(10)

where $(\mathbf{M})_{i,j}\coloneqq\exp(\epsilon^{-1}(f_{i}+g_{j}-c_{i,j}))-1$
and $\mathbf{f},\mathbf{g}$ are now unconstrained. For the dual pair
$(\mathbf{f}^{*},\mathbf{g}^{*})$ maximizing (10), we have
$\nabla_{\mathbf{f},\mathbf{g}}\mathrm{DOT}_{\epsilon}(\mathbf{f},\mathbf{g})|_{(\mathbf{f^{*}},\mathbf{g^{*}})}=\mathbf{0}$.
This first-order condition leads to

$f_{i}^{*}=\mathcal{T}_{\epsilon}(\mathbf{b},\mathbf{g}^{*},\mathbf{C}_{i:}),\hskip
28.45274ptg_{i}^{*}=\mathcal{T}_{\epsilon}(\mathbf{a},\mathbf{f}^{*},\mathbf{C}_{:i}),$
(11)

where $\mathbf{C}_{i:}$ (resp. $\mathbf{C}_{:i}$) is the $i^{\text{th}}$ row
(resp. column) of $\mathbf{C}$. Here
$\mathcal{T}_{\epsilon}:\mathbb{R}^{N}\times\mathbb{R}^{N}\times\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$
denotes the mapping

$\mathcal{T}_{\epsilon}(\mathbf{a},\mathbf{f},\mathbf{C}_{:,i})=-\epsilon~{}\log\sum_{k}\exp\\{\log
a_{k}+\epsilon^{-1}\left(f_{k}-c_{k,i}\right)\big{\\}}.$

One may then recover the regularized transport matrix as

$p^{\mathrm{OT}}_{\epsilon,i,j}=a_{i}b_{j}\exp\left(\epsilon^{-1}(f^{*}_{i}+g^{*}_{j}-c_{i,j})\right).$
(12)

The dual can be maximized using the Sinkhorn algorithm introduced for OT in
the seminal paper of Cuturi (2013). Algorithm 2 presents the implementation of
Feydy et al. (2019) where the fixed point updates based on Equation 11 have
been stabilized.

Algorithm 2 Sinkhorn Algorithm

1: Function $\mathbf{Potentials}(\mathbf{a},\mathbf{b},\mathbf{u},\mathbf{v})$

2: Local variables: $\mathbf{f},\mathbf{g}\in\mathbb{R}^{N}$

3: Initialize: $\mathbf{f}=\mathbf{0}$, $\mathbf{g}=\mathbf{0}$

4: Set
$\mathbf{C}\leftarrow\mathbf{u}\mathbf{u}^{t}+\mathbf{v}\mathbf{v}^{t}-2\mathbf{u}\mathbf{v}^{t}$

5: while stopping criterion not met do

6: for $i\in[N]$ do

7:
$f_{i}\leftarrow\frac{1}{2}\left(f_{i}+\mathcal{T}_{\epsilon}(\mathbf{b},\mathbf{g},\mathbf{C}_{i:})\right)$

8:
$g_{i}\leftarrow\frac{1}{2}\left(g_{i}+\mathcal{T}_{\epsilon}(\mathbf{a},\mathbf{f},\mathbf{C}_{:i})\right)$

9: end for

10: end while

11: Return $\mathbf{f},\mathbf{g}$

The resulting dual vectors $(\mathbf{f}^{*},\mathbf{g}^{*})$ can then be
differentiated for example using automatic differentiation through the
Sinkhorn algorithm loop (Flamary et al., 2018), or more efficiently using
“gradient stitching” on the dual vectors at convergence, which we do here (see
Feydy et al. (2019) for details). The derivatives of
$\mathbf{P}^{\mathrm{OT}}_{\epsilon}$ are readily accessible by combining the
derivatives of (11) with the derivatives of (12), using automatic
differentiation at no additional cost.

###  3.2 Differentiable Ensemble Transform Resampling

We obtain a differentiable ET (DET), denoted $\mathbf{T}^{(t)}_{N,\epsilon}$,
by computing the entropy-regularized OT using Algorithm 3 for the weighted
particles ($\mathbf{X}_{t},\mathbf{w}_{t},N$) at time $t$

$\displaystyle\tilde{X}^{i}_{t}=N{\textstyle\sum}_{k=1}^{N}~{}p^{\mathrm{OT}}_{\epsilon,i,k}X^{k}_{t}\coloneqq\mathbf{T}^{(t)}_{N,\epsilon}(X^{i}_{t}).$
(13)

Algorithm 3 DET Resampling

1: Function $\mathbf{EnsembleTransform}(\mathbf{X},\mathbf{w},N$)

2: $\mathbf{f},\mathbf{g}\leftarrow$
$\mathbf{Potentials}(\mathbf{w},\frac{1}{N}\mathbf{1},\mathbf{X},\mathbf{X})$

3: for $i\in[N]$ do

4: for $j\in[N]$ do

5:
$p^{\mathrm{OT}}_{\epsilon,i,j}=\frac{w_{i}}{N}\exp\Big(\frac{f_{i}+g_{j}-c_{i,j}}{\epsilon}\Big{missing})$

6: end for

7: end for

8: Return $\tilde{\mathbf{X}}=N\mathbf{P}^{\mathrm{OT}}_{\epsilon}\mathbf{X}$

Compared to the ET, the DET is differentiable and can be computed at cost
$O(N^{2})$ as it relies on the Sinkhorn algorithm. This algorithm converges
quickly (Altschuler et al., 2017) and is particularly amenable to GPU
implementation.

The DPF proposed in this paper is similar to Algorithm 1 except that we sample
from the proposal $q_{\phi}$ using the reparameterization trick and Step 6 is
replaced by the DET. While such a differentiable approximation of the ET has
previously been suggested in ML (Cuturi & Doucet, 2014; Seguy et al., 2018),
it has never been realized before that this could be exploited to obtain a
DPF. In particular, we obtain differentiable estimates of expectations w.r.t.
the filtering distributions with respect to $\theta$ and $\phi$ and, for a
fixed “seed” $\mathbf{U}=\mathbf{u}$ 222Here $\mathbf{U}$ denotes only the set
of $\theta,\phi$-independent random variables used to generate particles as,
contrary to standard PF, DET resampling does not rely on any additional random
variable., we obtain a differentiable estimate of the log-likelihood function
$\theta\mapsto\hat{\ell}_{\epsilon}(\theta;\phi,\mathbf{u})$.

Like ET, DET only satisfies (2) for affine functions $\psi$. Unlike
$\mathbf{P}^{\mathrm{OT}}$, $\mathbf{P}^{\mathrm{OT}}_{\epsilon}$ is sensitive
to the scale of $\mathbf{X}_{t}$. To mitigate this sensitivity, one may
compute
$\delta(\mathbf{X}_{t})=\sqrt{d_{x}}\max_{k\in[d_{x}]}\text{std}_{i}(X^{i}_{t,k})$
for $\mathbf{X}_{t}\in\mathbb{R}^{N\times d_{x}}$ and rescale $\mathbf{C}$
accordingly to ensure that $\epsilon$ is approximately independent of the
scale and dimension of the problem.

##  4 Theoretical Analysis

We show here that the gradient estimates of PF-based losses ignoring gradients
terms due to resampling are not consistent and can suffer from a large non-
vanishing bias. On the contrary, we establish that DPF provides consistent and
differentiable estimates of the filtering distributions and log-likelihood
function. This is achieved by obtaining novel quantitative convergence results
for the DET.

###  4.1 Gradient Bias from Ignoring Resampling Terms

We first provide theoretical results on the asymptotic bias of the gradient
estimates computed from PF-losses, by dropping the gradient terms from
resampling, as adopted in (Hirt & Dellaportas, 2019; Jonschkowski et al.,
2018; Karkus et al., 2018; Le et al., 2018; Ma et al., 2020b; Maddison et al.,
2017; Naesseth et al., 2018). We limit ourselves here to the ELBO loss.
Similar analysis can be carried out for the non-differentiable resampling
schemes and losses considered in robotics.

######  Proposition 4.1.

Consider the PF in Algorithm 1 where $\phi$ is distinct from $\theta$ then,
under regularity conditions, the expectation of the ELBO gradient estimate
$\hat{\nabla}_{\theta}\ell^{\textup{ELBO}}(\theta,\phi)$ ignoring resampling
terms considered in (Le et al., 2018; Maddison et al., 2017; Naesseth et al.,
2018) converges as $N\rightarrow\infty$ to

$\displaystyle\mathbb{E}[\hat{\nabla}_{\theta}\ell^{\textup{ELBO}}(\theta,\phi)]{\rightarrow}\int\nabla_{\theta}\log
p_{\theta}(x_{1},y_{1})~{}p_{\theta}(x_{1}|y_{1})\textup{d}x_{1}$
$\displaystyle+{\textstyle\sum\limits_{t=2}^{T}}\int\nabla_{\theta}\log
p_{\theta}(x_{t},y_{t}|x_{t-1})~{}p_{\theta}(x_{t-1:t}|y_{1:t})\textup{d}x_{t-1:t}$

whereas Fisher’s identity yields

$\displaystyle\nabla_{\theta}\ell(\theta)=\int\nabla_{\theta}\log
p_{\theta}(x_{1},y_{1})~{}p_{\theta}(x_{1}|y_{1:T})\textup{d}x_{1}$ (14)
$\displaystyle+{\textstyle\sum\limits_{t=2}^{T}}\int\nabla_{\theta}\log
p_{\theta}(x_{t},y_{t}|x_{t-1})~{}p_{\theta}(x_{t-1:t}|y_{1:T})\textup{d}x_{t-1:t}.$

Hence, whereas we have
$\nabla_{\theta}\ell^{\textup{ELBO}}(\theta,\phi)\rightarrow\nabla_{\theta}\ell(\theta)$
as $N\rightarrow\infty$ under regularity assumptions, the asymptotic bias of
$\hat{\nabla}_{\theta}\ell^{\textup{ELBO}}(\theta,\phi)$ only vanishes if
$p_{\theta}(x_{t-1:t}|y_{1:t})=p_{\theta}(x_{t-1:t}|y_{1:T})$; i.e. for models
where the $X_{t}$ are independent. When $y_{t+1:T}$ do not bring significant
information about $X_{t}$ given $y_{t:T}$, as for the models considered in (Le
et al., 2018; Maddison et al., 2017; Naesseth et al., 2018), this is a
reasonable approximation which explains the good performance reported therein.
However, we show in Section 5 that this bias can also lead practically to
inaccurate parameter estimation.

###  4.2 Quantitative Bounds on the DET

Weak convergence results for the ET have been established in (Reich, 2013;
Myers et al., 2021) and the DET in (Seguy et al., 2018). We provide here the
first quantitative bound for the ET ($\epsilon=0$) and DET ($\epsilon>0$)
which holds for any $N\geq 1$ by building upon results of (Li & Nochetto,
2021) and (Weed, 2018). We use the notation
$\nu(\psi)\coloneqq\int\psi(x)\nu(\textup{d}x)$ for any measure $\nu$ and
function $\psi$.

######  Proposition 4.2.

Consider atomic probability measures
$\alpha_{N}=\sum_{i=1}^{N}a_{i}\delta_{Y^{i}}$ with $a_{i}>0$ and
$\beta_{N}=\sum_{i=1}^{N}b_{i}\delta_{X^{i}}$, with support
$\mathcal{X}\subset\mathbb{R}^{d}$. Let
$\tilde{\beta}_{N}=\sum_{i=1}^{N}a_{i}\delta_{\tilde{X}^{i}_{N,\epsilon}}$
where
$\tilde{\mathbf{X}}_{N,\epsilon}=\Delta^{-1}\mathbf{P}^{\mathrm{OT}}_{\epsilon}\mathbf{X}$
for $\Delta=\textup{diag}(a_{1},...,a_{N})$ and
$\mathbf{P}^{\mathrm{OT}}_{\epsilon}$ is the transport matrix corresponding to
the $\epsilon$-regularized OT coupling,
$\mathcal{P}^{\mathrm{OT},N}_{\epsilon}$, between $\alpha_{N}$ and
$\beta_{N}$. Let $\alpha,\beta$ be two other probability measures, also
supported on $\mathcal{X}$, such that there exists a unique
$\lambda$-Lipschitz optimal transport map $\mathbf{T}$ between them. Then for
any bounded $1$-Lipschitz function $\psi$, we have

$\left|\beta_{N}(\psi)-\tilde{\beta}_{N}(\psi)\right|\leq
2\lambda^{1/2}\mathcal{E}^{1/2}\left[\mathfrak{d}^{1/2}+\mathcal{E}\right]^{1/2}\\\
+\max\\{\lambda,1\\}\left[\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)\right],$
(15)

where $\mathfrak{d}\coloneqq\sup_{x,y\in\mathcal{X}}|x-y|$ and
$\mathcal{E}=\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)+\sqrt{2\epsilon\log
N}.$

If $\mathcal{W}_{2}(\alpha_{N},\alpha),\mathcal{W}_{2}(\beta_{N},\beta)\to 0$
and we choose $\epsilon_{N}=o(1/\log N)$ the bound given in (15) vanishes with
$N\to\infty$. This suggested dependence of $\epsilon$ on $N$ comes from the
entropic radius, see Lemma C.1 in the Supplementary and (Weed, 2018), and is
closely related to the fact that entropy-regularized OT is sensitive to the
scale of $\mathbf{X}$. Equivalently one may rescale $\mathbf{X}$ by a factor
$\log N$ when computing the cost matrix. In particular when $\alpha_{N}$ and
$\beta_{N}$ are Monte Carlo approximations of $\alpha$ and $\beta$, we expect
$\mathcal{W}_{2}(\alpha_{N},\alpha),\mathcal{W}_{2}(\beta_{N},\beta)=O(N^{-1/d})$
with high probability (Fournier & Guillin, 2015).

###  4.3 Consistency of DPF

The parameters $\theta,\phi$ are here fixed and omitted from notation. We now
establish consistency results for DPF, showing that both the resulting
particle approximations
$\tilde{\beta}^{(t)}_{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{\tilde{X}^{i}_{t}}$
of $\beta^{(t)}=p(x_{t}|y_{1:t})$ and the corresponding log-likelihood
approximation $\log\hat{p}_{N}(y_{1:T})$ of $\log p(y_{1:T})$ are consistent.
In the interest of simplicity, we limit ourselves to the scenario where the
proposal is the transition, $q=f$, so
$\omega(x_{t-1},x_{t},y_{t})=g(y_{t}|x_{t})$, known as the bootstrap PF and
study a slightly non-standard version of it proposed in (Del Moral & Guionnet,
2001); see Appendix D for details. Consistency is established under regularity
assumptions detailed in the Supplementary. Assumption B.1 is that the space
$\mathcal{X}\subset\mathbb{R}^{d}$ has a finite diameter $\mathfrak{d}$.
Assumption B.2 implies that the proposal mixes exponentially fast in the
Wasserstein sense at a rate $\kappa$, which is reasonable given compactness,
and essential for the error to not accumulate. Assumption B.3 assumes a
bounded importance weight function i.e.
$g(y_{t}|x_{t})\in[\Delta,\Delta^{-1}]$, again not unreasonable given
compactness. Assumption B.4 states that at each time step, the optimal
transport problem between $\alpha^{(t)}$ and $\beta^{(t)}$ is solved uniquely
by a deterministic, globally Lipschitz map. Uniqueness is crucial for the
quantitative stability results provided in the following proposition.

######  Proposition 4.3.

Under Assumptions B.1, B.2, B.3 and B.4, for any $\delta>0$, with probability
at least $1-2\delta$ over the sampling steps, for any bounded 1-Lipschitz
$\psi$, for any $t\in[1:T]$, the approximations of the filtering distributions
and log-likelihood computed by the bootstrap DPF satisfy

$\displaystyle|\tilde{\beta}^{(t)}_{N}(\psi)-\beta^{(t)}(\psi)|$
$\displaystyle\leq\mathfrak{G}_{\epsilon,\delta/T,N,d}^{(t)}\left(\lambda(c,C,d,T,N,\delta)\right),$
$\displaystyle\left|\log\frac{\hat{p}_{N}(y_{1:T})}{p(y_{1:T})}\right|$
$\displaystyle\leq\frac{\kappa}{\Delta}\max_{t\in[1:T]}\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]$
$\displaystyle\times\sum_{t=1}^{T}\mathfrak{G}_{\epsilon,\delta/T,N,d}^{(t)}\left(\lambda(c,C,d,T,N,\delta)\right),$

for
$\lambda(c,C,d,T,N,\delta)=\sqrt{f_{d}^{-1}\left(\frac{\log(CT/\delta)}{cN}\right)}$
where $c,C$ are finite constants independent of $T,$ and $\mathrm{Lip}[f]$ is
the Lipschitz constant of the function $f,$ and
$\mathfrak{G}^{(t)}_{N,\epsilon},f_{d}$ defined in Appendix D are two
functions such that if we set $\epsilon_{N}=o(1/\log N)$ then we have in
probability

$|\tilde{\beta}^{(t)}_{N}(\psi)-\beta^{(t)}(\psi)|\to
0,\quad\left|\log\frac{\hat{p}_{N}(y_{1:T})}{p(y_{1:T})}\right|\to 0.$

The above bounds are certainly not sharp. A glimpse into the behavior of the
above bounds in terms of $T$ can be obtained through careful consideration of
the quantities appearing in Proposition D.1 in the supplement. In particular,
for $\kappa$ small enough, it suggests that the bound on the error of the log-
likelihood estimator grows linearly with $T$ as for standard PF under mixing
assumptions. Sharper bounds are certainly possible, e.g. using a $L_{1}$
version of Theorem 3.5 in (Li & Nochetto, 2021). It would also be of interest
to weaken the assumptions, in particular, to remove the bounded space
assumption although it is very commonly made in the PF literature to obtain
quantitative bounds; see e.g. (Del Moral, 2004; Douc et al., 2014). Although
this is not made explicit in the expressions above, there is an exponential
dependence of the bounds on the state dimension $d_{x}$. This is unavoidable
however and a well-known limitation of PF methods.

Finally note that DPF provides a biased estimate of the likelihood contrary to
standard PF, so we cannot guarantee that the expectation of its logarithm,
$\ell_{\epsilon}^{\mathrm{ELBO}}(\theta,\phi)\coloneqq\mathbb{E}_{\mathbf{U}}[\hat{\ell}_{\epsilon}(\theta;\phi,\mathbf{U})]$.
is actually a valid ELBO. However in all our experiments, see e.g. Section
5.1,
$|\ell^{\mathrm{ELBO}}_{\epsilon}(\theta,\phi)-\ell^{\mathrm{ELBO}}(\theta,\phi)|$
is significantly smaller than $\ell(\theta)-\ell^{\mathrm{ELBO}}(\theta,\phi)$
so $\ell^{\mathrm{ELBO}}_{\epsilon}(\theta,\phi)<\ell(\theta)$. Hence we keep
the ELBO terminology.

##  5 Experiments

In Section 5.1, we assess the sensitivity of the DPF to the regularization
parameter $\epsilon$. All other DPF experiments presented here use the DET
Resampling detailed in Algorithm 3 with $\epsilon=0.5$, which ensures
stability of the gradient calculations while adding little bias to the
calculation of the ELBO compared to standard PF. Our method is implemented in
both PyTorch and TensorFlow, the code to replicate the experiments as well as
further experiments may be found at https://github.com/JTT94/filterflow.

###  5.1 Linear Gaussian State-Space Model

We consider here a simple two-dimensional linear Gaussian SSM for which the
exact likelihood can be computed exactly using the Kalman filter

$\displaystyle X_{t+1}|\\{X_{t}=x\\}$
$\displaystyle\sim\mathcal{N}\left(\text{diag}(\theta_{1}~{}\theta_{2})x,0.5\mathbf{I}_{2}\right),$
$\displaystyle Y_{t}|\\{X_{t}=x\\}$
$\displaystyle\sim\mathcal{N}(x,0.1\mathbf{I}_{2}).$

We simulate $T=150$ observations using
$\theta=(\theta_{1},\theta_{2})=(0.5,0.5)$, for which we evaluate the ELBO at
$\theta=(0.25,0.25)$, $\theta=(0.5,0.5)$, and $\theta=(0.75,0.75)$. More
precisely, using a standard PF with $N=25$ particles, we compute the mean and
standard deviation of
$\frac{1}{T}(\hat{\ell}(\theta;\mathbf{U})-\ell(\theta))$ over 100
realizations of $\bf{U}$. The mean is an estimate of the ELBO minus the true
log-likelihood (rescaled by $1/T$). We then perform the same calculations for
the DPF using the same number of particles and $\epsilon=0.25,0.5,0.75$. As
mentioned in Section 3.2 and Section 4.3, the DET resampling scheme is only
satisfying Equation 2 for affine functions $\psi$ so the DPF provides a biased
estimate of the likelihood. Hence we cannot guarantee that the expectation of
the corresponding log-likelihood estimate is a true ELBO. However, from Table
1, we observe that the difference between the ELBO estimates computed using PF
and DPF is negligible for the three values of $\epsilon$. The standard
deviation of the log-likelihood estimates is also similar.

Table 1: Mean & std of $\frac{1}{T}(\hat{\ell}(\theta;\mathbf{U})-\ell(\theta))$ $\theta_{1}$, $\theta_{2}$ | 0.25 | 0.5 | 0.75
---|---|---|---
PF | mean | -1.13 | -0.93 | -1.05
std | 0.20 | 0.18 | 0.17
DPF ($\epsilon=0.25$) | mean | -1.14 | -0.94 | -1.07
std | 0.20 | 0.18 | 0.19
DPF ($\epsilon=0.5$) | mean | -1.14 | -0.94 | -1.08
std | 0.20 | 0.18 | 0.18
DPF ($\epsilon=0.75$) | mean | -1.14 | -0.94 | -1.08
std | 0.20 | 0.18 | 0.18

Recall here that alternative techniques estimating the score vector
$\nabla_{\theta}\ell(\theta)$ by approximating (14) using particle smoothing
algorithms (Poyiadjis et al., 2011; Kantas et al., 2015) could also be used to
estimate $\theta$.

###  5.2 Learning the Proposal Distribution

We consider a similar example as in (Naesseth et al., 2018) where one learns
the parameters $\phi$ of the proposal using the ELBO for the following linear
Gaussian SSM:

$\displaystyle
X_{t+1}|\\{X_{t}=x\\}\sim\mathcal{N}\left(\mathbf{A}x,\mathbf{I}_{d_{x}}\right),$
(16) $\displaystyle
Y_{t}|\\{X_{t}=x\\}\sim\mathcal{N}(\mathbf{I}_{d_{y},d_{x}}x,\mathbf{I}_{d_{y}}),$
(17)

with $\mathbf{A}=(0.42^{|i-j|+1})_{1\leq i,j\leq d_{x}}$,
$\mathbf{I}_{d_{y},d_{x}}$ is a $d_{y}\times d_{x}$ matrix with $1$ on the
diagonal for the $d_{y}$ first rows and zeros elsewhere. For
$\mathbf{\phi}\in\mathbb{R}^{d_{x}+d_{y}}$, we consider

$q_{\phi}(x_{t}|x_{t-1},y_{t})=\mathcal{N}(x_{t}|\Delta_{\phi}^{-1}\left(\mathbf{A}x_{t-1}+\Gamma_{\phi}y_{t}\right),\Delta_{\phi}),$

with $\Delta_{\phi}=\text{diag}(\phi_{1},\dots,\phi_{d_{x}})$ and a
$d_{x}\times d_{y}$ matrix
$\Gamma_{\phi}=\text{diag}_{d_{x},d_{y}}(\phi_{1},\dots,\phi_{d_{x}})$ with
$\phi_{i}$ on the diagonal for $d_{x}$ first rows and zeros elsewhere. The
locally optimal proposal $p(x_{t}|x_{t-1},y_{t})\propto
g(y_{t}|x_{t})f(x_{t}|x_{t-1})$ in (Doucet & Johansen, 2009) corresponds to
$\mathbf{\phi}=\mathbf{1}$, the vector with unit entries of dimension
$d_{\phi}=d_{x}+d_{y}$.

For $d_{x}=25,d_{y}=1$, $M=100$ realizations of $T=100$ observations using
(16)-(17), we learn $\phi$ on each realization using 100 steps of stochastic
gradient ascent with learning rate $0.1$ on the $\ell^{\text{ELBO}}(\phi)$
using regular PF with biased gradients as in (Maddison et al., 2017; Le et
al., 2018; Naesseth et al., 2018) and $\ell^{\text{ELBO}}(\phi)$ with four
independent filters using DPF. We use $N=500$ for regular PF and $N=25$ for
DPF so as to match the computational complexity. While
$p(x_{t}|x_{t-1},y_{t})$ is not guaranteed to maximize the ELBO, our
experiments showed that it outperforms optimized proposals. We therefore
report the RMSE of $\mathbf{\phi}-\mathbf{1}$ and the average Effective Sample
Size (ESS) (Doucet & Johansen, 2009) as proxy performance metrics. On both
metrics, DPF outperforms regular PF. The RMSE over 100 experiments is 0.11 for
DPF vs 0.22 for regular PF while the average ESS after convergence is around
60% for DPF vs 25% for regular PF. The average time per iteration was around
15 seconds for both DPF and PF.

###  5.3 Variational Recurrent Neural Network (VRNN)

A VRNN is an SSM introduced by (Chung et al., 2015) to improve upon LSTMs
(Long Short Term Memory networks) with the addition of a stochastic component
to the hidden state, this extends variational auto-encoders to a sequential
setting. Indeed let latent state be $X_{t}=(R_{t},Z_{t})$ where $R_{t}$ is an
RNN state and $Z_{t}$ a latent Gaussian variable, here $Y_{t}$ is a vector of
binary observations. The VRNN is detailed as follows. $\text{RNN}_{\theta}$
denotes the forward call of an LSTM cell which at time $t$ emits the next RNN
state $R_{t+1}$ and output $O_{t+1}$. $E_{\theta}$, $h_{\theta}$,
$\mu_{\theta}$, $\sigma_{\theta}$ are fully connected neural networks;
detailed fully in the Supplementary Material. This model is trained on the
polyphonic music benchmark datasets (Boulanger-Lewandowski et al., 2012),
whereby $Y_{t}$ represents which notes are active. The observation sequences
are capped to length $150$ for each dataset, with each observation of
dimension $88$. We chose latent states $Z_{t}$ and $R_{t}$ to be of dimension
$d_{z}=8$ and $d_{r}=16$ respectively so $d_{x}=24$. We use
$q_{\phi}(x_{t}|x_{t-1},y_{t})=f_{\theta}(x_{t}|x_{t-1})$.

$\displaystyle(R_{t+1},O_{t+1})$
$\displaystyle=\text{RNN}_{\theta}(R_{t},Y_{1:t},E_{\theta}(Z_{t})),$
$\displaystyle Z_{t+1}$
$\displaystyle\sim\mathcal{N}(\mu_{\theta}(O_{t+1}),\sigma_{\theta}(O_{t+1})),$
$\displaystyle\hat{p}_{t+1}$
$\displaystyle=h_{\theta}(E_{\theta}(Z_{t+1}),O_{t+1}),$ $\displaystyle
Y_{t}|X_{t}$ $\displaystyle\sim\text{Ber}(\hat{p}_{t}).$

Table 2: ELBO $\pm$ Standard Deviation evaluated using Test Data. | MuseData | JSB | Nottingham
---|---|---|---
DPF | $-\mathbf{7.59}_{\pm 0.01}$ | $-\mathbf{7.67}_{\pm 0.08}$ | $\mathbf{-3.79}_{\pm 0.02}$
PF | $-7.60_{\pm 0.06}$ | $-7.92_{\pm 0.13}$ | $-3.81_{\pm 0.02}$
SPF | $-7.73_{\pm 0.14}$ | $-8.17_{\pm 0.07}$ | $-3.91_{\pm 0.05}$

The VRNN model is trained by maximizing
$\ell_{\epsilon}^{\text{ELBO}}(\theta)$ using DPF. We compare this to the same
model trained by maximizing $\ell^{\text{ELBO}}(\theta)$ computed with regular
PF (Maddison et al., 2017) and also trained with ‘soft-resampling’ (SPF)
introduced by (Karkus et al., 2018) and described in Section 1.3, SPF is used
here with parameter $\alpha=0.1$. Unlike regular resampling, SPF partially
preserves a gradient through the resampling step, however SPF still involves a
non-differentiable operation, again resulting in a biased gradient. SPF also
produces higher variance estimates as the resampled approximation is not
uniformly weighted, essentially interpolating between PF and IWAE. Each of the
methods are performed with $N=32$ particles. Although DET is computationally
more expensive than the other resampling schemes, the computational times of
DPF, PF, and SPF are very similar due to most of the complexity coming from
neural network operations. The learned models are then evaluated on test data
using multinomial resampling for comparable ELBO results. Due to the fact that
our observation model is $\text{Ber}(\hat{p}_{t})$, this recovers the negative
log-predictive cross-entropy.

Table 2 illustrates the benefit of using DPF over regular PF and SPF for the
JSB dataset. Although DPF remains competitive compared to other heuristic
approaches, the difference is relatively minor for the other datasets. We
speculate that the performance of the heuristic methods is likely due to low
predictive uncertainty for the next observation given the previous one.

###  5.4 Robot Localization

Consider the setting of a robot/agent in a maze (Jonschkowski et al., 2018;
Karkus et al., 2018). Given the agent’s initial state, $S_{1}$, and inputs
$a_{t}$, one would like to infer the location of the agent at any specific
time given observations $O_{t}$. Let the latent state be denoted
$S_{t}=(X^{(1)}_{t},X^{(2)}_{t},\gamma_{t})$ where $(X^{(1)}_{t},X^{(2)}_{t})$
are location coordinates and $\gamma_{t}$ the robot’s orientation. In our
setting observations $O_{t}$ are images, which are encoded to extract useful
features using a neural network $E_{\theta}$, where $Y_{t}=E_{\theta}(O_{t})$.
This problem requires learning the relationship between the robot’s location,
orientation and the observations. Given actions
$a_{t}=(v^{(1)}_{t},v^{(2)}_{t},\omega_{t})$, we have

$\displaystyle S_{t+1}$
$\displaystyle=F_{\theta}(S_{t},a_{t})+\nu_{t},\quad\nu_{t}\overset{\textup{i.i.d.}}{\sim}\mathcal{N}(\mathbf{0},\Sigma_{F}),$
$\displaystyle Y_{t}$
$\displaystyle=G_{\theta}(S_{t})+\epsilon_{t},\quad\epsilon_{t}\overset{\textup{i.i.d.}}{\sim}\mathcal{N}(\mathbf{0},\sigma_{G}^{2}\mathbb{I}_{e_{d}}),$

where
$\Sigma_{F}=\textup{diag}(\sigma_{x}^{2},\sigma_{x}^{2},\sigma_{\theta}^{2})$
and the relationship between state $S_{t}$ and image encoding $Y_{t}$ may be
parameterized by another neural network $G_{\theta}$. We consider here a
simple linear model of the dynamics

$F(S_{t},a_{t})=\begin{bmatrix}X^{(1)}_{t}+v^{(1)}_{t}\cos(\gamma_{t})+v^{(2)}_{t}\sin(\gamma_{t})\\\
X^{(2)}_{t}+v^{(1)}_{t}\sin(\gamma_{t})-v^{(2)}_{t}\cos(\gamma_{t})\\\
\gamma_{t}+\omega_{t}\end{bmatrix}.$

$D_{\theta}$ denotes a decoder neural network, mapping the encoding back to
the original image. $E_{\theta}$, $G_{\theta}$ and $D_{\theta}$ are trained
using a loss function consisting of the PF-estimated log-likelihood
$\hat{\mathcal{L}}_{\mathrm{PF}}$; PF-based mean squared error (MSE),
$\hat{\mathcal{L}}_{\mathrm{MSE}}$; and auto-encoder loss,
$\hat{\mathcal{L}}_{\mathrm{AE}}$, given per-batch as in (Wen et al., 2020):

$\displaystyle\hat{\mathcal{L}}_{\mathrm{MSE}}$
$\displaystyle\coloneqq\frac{1}{T}\sum_{t=1}^{T}||X^{*}_{t}-\sum_{i=1}^{N}w^{i}_{t}X^{i}_{t}||^{2},~{}~{}\hat{\mathcal{L}}_{\mathrm{PF}}:=-\frac{1}{T}\hat{\ell}(\theta),$
$\displaystyle\hat{\mathcal{L}}_{\mathrm{AE}}$
$\displaystyle\coloneqq\sum_{t=1}^{T}||D_{\theta}(E_{\theta}(O_{t}))-O_{t}||^{2},$

where $X^{\star}_{t}$ are the true states available from training data and
$\sum_{i=1}^{N}w^{i}_{t}X^{i}_{t}$ are the PF estimates of
$\mathbb{E}[X_{t}|y_{1:t}]$. The auto-encoder / reconstruction loss
$\hat{\mathcal{L}}_{\mathrm{AE}}$ ensures the encoder is informative and
prevents the case whereby networks $G_{\theta}$, $E_{\theta}$ map to a
constant. The PF-based loss terms $\hat{\mathcal{L}}_{\mathrm{MSE}}$ and
$\hat{\mathcal{L}}_{\mathrm{PF}}$ are not differentiable w.r.t. $\theta$ under
traditional resampling schemes.

Figure 2: Left: Particles $(X^{(1),i}_{t},X^{(2),i}_{t})$ (green), PF estimate
of $\mathbb{E}[X_{t}|y_{1:t}]$ (blue), true state $X_{t}^{*}$ (red). Right:
Observation, $O_{t}$. Table 3: MSE and $\pm$ Standard Deviation evaluated on
Test Data: Lower is better
| Maze 1 | Maze 2 | Maze 3
---|---|---|---
DPF | $\mathbf{3.55}_{\pm\mathbf{0.20}}$ | $\mathbf{4.65}_{\pm 0.50}$ | $\mathbf{4.44}_{\pm\mathbf{0.26}}$
PF | $10.71_{\pm 0.45}$ | $11.86_{\pm 0.57}$ | $12.88_{\pm 0.65}$
SPF | $9.14_{\pm 0.39}$ | $10.12_{\pm\mathbf{0.40}}$ | $11.42_{\pm 0.37}$

We use the setup from (Jonschkowski et al., 2018) with data from DeepMind Lab
(Beattie et al., 2016). This consists of $3$ maze layouts of varying sizes. We
have access to ‘true’ trajectories of length $1,000$ steps for each maze. Each
step has an associated state, action and observation image, as described
above. The visual observation $O_{t}$ consists of $32\times 32$ RGB pixel
images, compressed to $24\times 24$, as shown in Figure 2. Random, noisy
subsets of fixed length are sampled at each training iteration. To illustrate
the benefits of our proposed method, we select the random subsets to be of
length $50$ as opposed to length $20$ as chosen in (Jonschkowski et al.,
2018). Training details in terms of learning rates, number of training steps
and neural network architectures for $E_{\theta}$, $G_{\theta}$ and
$D_{\theta}$ are given in the Appendices.

We compare our method, DPF, to regular PF used in (Maddison et al., 2017) and
Soft PF (SPF) used in (Karkus et al., 2018; Ma et al., 2020a, b), whereby the
soft resampling is used with $\alpha=0.1$. As most of the computational
complexity arises from neural network operations, DPF is of similar overall
computational cost to SPF and PF. As shown in Table 3 and Figure 3, DPF
significantly outperforms previously considered PF methods in this experiment.
The observation model becomes increasingly important for longer sequences due
to resampling and weighting operations. Indeed, as shown in Figure 4, the
error is small for both models at the start of the sequence, however the error
at later stages in the sequence is visibly smaller for the model trained using
DPF.

(a) Maze 1

(b) Maze 2

(c) Maze 3

Figure 3: MSE of PF (red), SPF (green) and DPF (blue) estimates, evaluated on
test data during training.

(a) Standard PF

(b) Differentiable PF

Figure 4: Illustrative Example: PF estimate of path compared to true path
(black) on a single $50$-step trajectory from test data.

##  6 Discussion

This paper introduces the first principled, fully differentiable PF (DPF)
which can use general proposal distributions. It provides a differentiable
estimate of the log-likelihood function and more generally differentiable
estimates of PF-based losses. This permits parameter inference in state-space
models and proposal distributions, using end-to-end gradient based
optimization. This also allows the use of PF routines in general
differentiable programming pipelines, in particular as a differentiable
sampling method for inference in probabilistic programming languages (Dillon
et al., 2017; Ge et al., 2018; van de Meent et al., 2018).

For a given number of particles $N$, existing PF methods ignoring resampling
gradient terms have computational complexity $O(N)$. Training with these
resampling schemes however is unreliable and performance cannot be improved by
increasing $N$ as gradient estimates are inconsistent and the limiting bias
can be significant. DPF has complexity $O(N^{2})$ during training. However,
this cost is dwarfed when training large neural networks. Additionally, once
the model is trained, standard PF may be ran at complexity $O(N)$. The
benefits of DPF are confirmed by our experimental results where it was shown
to outperform existing techniques, even when an equivalent computational
budget was used. Moreover, recent techniques have been proposed to speed up
the Sinkhorn algorithm (Altschuler et al., 2019; Scetbon & Cuturi, 2020) at
the core of DPF and could potentially be used here to reduce its complexity.

Regularization parameter $\epsilon$ was not fine-tuned in our experiments. In
future work, it would be interesting to obtain sharper quantitative bounds on
DPF to propose principled guidelines on choosing $\epsilon,$ further improving
its performance. Finally, we have focused on the use of the differentiable
ensemble transform to obtain a differentiable resampling scheme. However,
alternative OT approaches could also be proposed such as a differentiable
version of the second order ET presented in (Acevedo et al., 2017), techniques
based on point cloud optimization (Cuturi & Doucet, 2014; Peyré & Cuturi,
2019) relying on the Sinkhorn divergence (Genevay et al., 2018) or the sliced-
Wasserstein metric. Alternative non-entropic regularizations, such as the
recently proposed Gaussian smoothed OT (Goldfeld & Greenewald, 2020), could
also lead to DPFs of interest.

## Acknowledgments

Adrien Corenflos was supported by the Academy of Finland (projects 321900 and
321891). Arnaud Doucet is supported by the EPSRC CoSInES (COmputational
Statistical INference for Engineering and Security) grant EP/R034710/1, James
Thornton by the OxWaSP CDT through grant EP/L016710/1. Computing resources
were provided through the Google Cloud Research Credits Programme. The authors
thank Valentin De Bortoli and Yuyang Shi for their comments and Laurence
Aitchison for bringing attention to reference (Aitchison, 2019).

## References

  * Acevedo et al. (2017) Acevedo, W., de Wiljes, J., and Reich, S.  Second-order accurate ensemble transform particle filters.  _SIAM Journal on Scientific Computing_ , 39(5):A1834–A1850, 2017.
  * Aitchison (2019) Aitchison, L.  Tensor Monte Carlo: particle methods for the GPU era.  _Advances in Neural Information Processing Systems_ , 2019.
  * Altschuler et al. (2017) Altschuler, J., Niles-Weed, J., and Rigollet, P.  Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration.  In _Advances in Neural Information Processing Systems_ , pp. 1964–1974, 2017.
  * Altschuler et al. (2019) Altschuler, J., Bach, F., Rudi, A., and Niles-Weed, J.  Massively scalable Sinkhorn distances via the Nyström method.  In _Advances in Neural Information Processing Systems_ , pp. 4429–4439, 2019.
  * Archer et al. (2015) Archer, E., Park, I. M., Buesing, L., Cunningham, J., and Paninski, L.  Black box variational inference for state space models.  _arXiv preprint arXiv:1511.07367_ , 2015.
  * Beattie et al. (2016) Beattie, C., Leibo, J. Z., Teplyashin, D., Ward, T., Wainwright, M., Küttler, H., Lefrancq, A., Green, S., Valdés, V., Sadik, A., Schrittwieser, J., Anderson, K., York, S., Cant, M., Cain, A., Bolton, A., Gaffney, S., King, H., Hassabis, D., Legg, S., and Petersen, S.  DeepMind Lab, 2016.
  * Bertsimas & Tsitsiklis (1997) Bertsimas, D. and Tsitsiklis, J. N.  _Introduction to Linear Optimization_.  Athena Scientific Belmont, MA, 1997.
  * Boulanger-Lewandowski et al. (2012) Boulanger-Lewandowski, N., Bengio, Y., and Vincent, P.  Modeling temporal dependencies in high-dimensional sequences: Application to polyphonic music generation and transcription.  In _International Conference on Machine Learning_ , pp. 1881–1888, 2012.
  * Burda et al. (2016) Burda, Y., Grosse, R. B., and Salakhutdinov, R.  Importance weighted autoencoders.  In _International Conference on Learning Representations_ , 2016.
  * Chopin & Papaspiliopoulos (2020) Chopin, N. and Papaspiliopoulos, O.  _An Introduction to Sequential Monte Carlo_.  Springer, 2020.
  * Chung et al. (2015) Chung, J., Kastner, K., Dinh, L., Goel, K., Courville, A. C., and Bengio, Y.  A recurrent latent variable model for sequential data.  In _Advances in Neural Information Processing Systems_ , pp. 2980–2988, 2015.
  * Cuturi (2013) Cuturi, M.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _Advances in Neural Information Processing Systems_ , pp. 2292–2300, 2013.
  * Cuturi & Doucet (2014) Cuturi, M. and Doucet, A.  Fast computation of Wasserstein barycenters.  In _International Conference on Machine Learning_ , pp. 685–693, 2014.
  * DeJong et al. (2013) DeJong, D. N., Liesenfeld, R., Moura, G. V., Richard, J.-F., and Dharmarajan, H.  Efficient likelihood evaluation of state-space representations.  _Review of Economic Studies_ , 80(2):538–567, 2013.
  * Del Moral (2004) Del Moral, P.  _Feynman-Kac Formulae_.  Springer, 2004.
  * Del Moral & Guionnet (2001) Del Moral, P. and Guionnet, A.  On the stability of interacting processes with applications to filtering and genetic algorithms.  In _Annales de l’Institut Henri Poincaré (B) Probability and Statistics_ , volume 37, pp. 155–194, 2001.
  * Dillon et al. (2017) Dillon, J. V., Langmore, I., Tran, D., Brevdo, E., Vasudevan, S., Moore, D., Patton, B., Alemi, A., Hoffman, M., and Saurous, R. A.  Tensorflow distributions.  _arXiv preprint arXiv:1711.10604_ , 2017.
  * Douc et al. (2014) Douc, R., Moulines, E., and Stoffer, D.  _Nonlinear Time Series: Theory, Methods and Applications with R Examples_.  CRC press, 2014.
  * Doucet & Johansen (2009) Doucet, A. and Johansen, A. M.  A tutorial on particle filtering and smoothing: Fifteen years later.  _Handbook of Nonlinear Filtering_ , 12:656–704, 2009.
  * Doucet & Lee (2018) Doucet, A. and Lee, A.  Sequential Monte Carlo methods.  _Handbook of Graphical Models_ , pp. 165–189, 2018.
  * Feydy et al. (2019) Feydy, J., Séjourné, T., Vialard, F.-X., Amari, S.-I., Trouvé, A., and Peyré, G.  Interpolating between optimal transport and MMD using Sinkhorn divergences.  In _International Conference on Artificial Intelligence and Statistics_ , 2019.
  * Finke et al. (2016) Finke, A., Doucet, A., and Johansen, A. M.  On embedded hidden Markov models and particle Markov chain Monte Carlo methods.  _arXiv preprint arXiv:1610.08962_ , 2016.
  * Flamary et al. (2018) Flamary, R., Cuturi, M., Courty, N., and Rakotomamonjy, A.  Wasserstein discriminant analysis.  _Machine Learning_ , 107(12):1923–1945, 2018\.
  * Fournier & Guillin (2015) Fournier, N. and Guillin, A.  On the rate of convergence in Wasserstein distance of the empirical measure.  _Probability Theory and Related Fields_ , 162(3):707–738, 2015.
  * Ge et al. (2018) Ge, H., Xu, K. X., and Ghahramani, Z.  Turing: A language for flexible probabilistic inference.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1682–1690, 2018.
  * Genevay et al. (2018) Genevay, A., Peyré, G., and Cuturi, M.  Learning generative models with Sinkhorn divergences.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1608–1617, 2018.
  * Goldfeld & Greenewald (2020) Goldfeld, Z. and Greenewald, K.  Gaussian-smooth optimal transport: Metric structure and statistical efficiency.  _arXiv preprint arXiv 2001.09206_ , 2020.
  * Graves (2016) Graves, A.  Stochastic backpropagation through mixture density distributions.  _arXiv preprint arXiv:1607.05690_ , 2016.
  * Hirt & Dellaportas (2019) Hirt, M. and Dellaportas, P.  Scalable Bayesian learning for state space models using variational inference with SMC samplers.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 76–86, 2019.
  * Jonschkowski et al. (2018) Jonschkowski, R., Rastogi, D., and Brock, O.  Differentiable particle filters: End-to-end learning with algorithmic priors.  In _Proceedings of Robotics: Science and Systems_ , 2018.
  * Kantas et al. (2015) Kantas, N., Doucet, A., Singh, S. S., Maciejowski, J., and Chopin, N.  On particle methods for parameter estimation in state-space models.  _Statistical Science_ , 30(3):328–351, 2015.
  * Karkus et al. (2018) Karkus, P., Hsu, D., and Lee, W. S.  Particle filter networks with application to visual localization.  In _Conference on Robot Learning_ , 2018.
  * Kingma & Welling (2014) Kingma, D. P. and Welling, M.  Auto-encoding variational Bayes.  In _International Conference on Learning Representations_ , 2014.
  * Kitagawa & Gersch (1996) Kitagawa, G. and Gersch, W.  _Smoothness Priors Analysis of Time Series_ , volume 116.  Springer Science & Business Media, 1996.
  * Klaas et al. (2005) Klaas, M., De Freitas, N., and Doucet, A.  Toward practical ${N}^{2}$ Monte Carlo: the marginal particle filter.  _Uncertainty in Artificial Intelligence_ , 2005.
  * Kloss et al. (2020) Kloss, A., Martius, G., and Bohg, J.  How to train your differentiable filter.  _arXiv preprint arXiv:2012.14313_ , 2020.
  * Krishnan et al. (2017) Krishnan, R. G., Shalit, U., and Sontag, D.  Structured inference networks for nonlinear state space models.  In _AAAI Conference on Artificial Intelligence_ , pp. 2101–2109, 2017.
  * Le et al. (2018) Le, T. A., Igl, M., Rainforth, T., Jin, T., and Wood, F.  Auto-encoding sequential Monte Carlo.  In _International Conference on Learning Representations_ , 2018.
  * Lee (2008) Lee, A.  Towards smooth particle filters for likelihood estimation with multivariate latent variables.  Master’s thesis, University of British Columbia, 2008.
  * Li & Nochetto (2021) Li, W. and Nochetto, R. H.  Quantitative stability and error estimates for optimal transport plans.  _IMA Journal of Numerical Analysis_ , 2021.
  * Lindsten & Schön (2013) Lindsten, F. and Schön, T. B.  Backward simulation methods for Monte Carlo statistical inference.  _Foundations and Trends® in Machine Learning_ , 6(1):1–143, 2013.
  * Ma et al. (2020a) Ma, X., Karkus, P., Hsu, D., and Lee, W. S.  Particle filter recurrent neural networks.  In _AAAI Conference on Artificial Intelligence_ , 2020a.
  * Ma et al. (2020b) Ma, X., Karkus, P., Ye, N., Hsu, D., and Lee, W. S.  Discriminative particle filter reinforcement learning for complex partial observations.  In _International Conference on Learning Representations_ , 2020b.
  * Maddison et al. (2017) Maddison, C. J., Lawson, D., Tucker, G., Heess, N., Norouzi, M., Mnih, A., Doucet, A., and Teh, Y. W.  Filtering variational objectives.  In _Advances in Neural Information Processing Systems_ , 2017.
  * Malik & Pitt (2011) Malik, S. and Pitt, M. K.  Particle filters for continuous likelihood evaluation and maximisation.  _Journal of Econometrics_ , 165(2):190–209, 2011\.
  * Murray et al. (2013) Murray, L. M., Jones, E. M., and Parslow, J.  On disturbance state-space models and the particle marginal Metropolis–Hastings sampler.  _SIAM/ASA Journal on Uncertainty Quantification_ , 1(1):494–521, 2013.
  * Myers et al. (2021) Myers, A., Thiery, A. H., Wang, K., and Bui-Thanh, T.  Sequential ensemble transform for Bayesian inverse problems.  _Journal of Computational Physics_ , 427:110055, 2021.
  * Naesseth et al. (2018) Naesseth, C. A., Linderman, S. W., Ranganath, R., and Blei, D. M.  Variational sequential Monte Carlo.  In _International Conference on Artificial Intelligence and Statistics_ , 2018.
  * Peyré & Cuturi (2019) Peyré, G. and Cuturi, M.  Computational optimal transport.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * Poyiadjis et al. (2011) Poyiadjis, G., Doucet, A., and Singh, S. S.  Particle approximations of the score and observed information matrix in state space models with application to parameter estimation.  _Biometrika_ , 98(1):65–80, 2011.
  * Rangapuram et al. (2018) Rangapuram, S. S., Seeger, M. W., Gasthaus, J., Stella, L., Wang, Y., and Januschowski, T.  Deep state space models for time series forecasting.  In _Advances in Neural Information Processing Systems_ , pp. 7785–7794, 2018.
  * Reich (2013) Reich, S.  A nonparametric ensemble transform method for Bayesian inference.  _SIAM Journal on Scientific Computing_ , 35(4):A2013–A2024, 2013.
  * Scetbon & Cuturi (2020) Scetbon, M. and Cuturi, M.  Linear time Sinkhorn divergences using positive features.  In _Advances in Neural Information Processing Systems_ , 2020.
  * Seguy et al. (2018) Seguy, V., Damodaran, B. B., Flamary, R., Courty, N., Rolet, A., and Blondel, M.  Large-scale optimal transport and mapping estimation.  In _International Conference on Learning Representations_ , 2018.
  * Thrun et al. (2005) Thrun, S., Burgard, W., and Fox, D.  _Probabilistic Robotics_.  MIT Press, 2005.
  * van de Meent et al. (2018) van de Meent, J.-W., Paige, B., Hongseok, Y., and Wood, F.  An introduction to probabilistic programming.  _arXiv preprint arXiv:1809.10756_ , 2018.
  * Villani (2008) Villani, C.  _Optimal Transport: Old and New_ , volume 338.  Springer Science & Business Media, 2008.
  * Weed (2018) Weed, J.  An explicit analysis of the entropic penalty in linear programming.  In _Proceedings of the 31st Conference On Learning Theory_ , 2018\.
  * Wen et al. (2020) Wen, H., Chen, X., Papagiannis, G., Hu, C., and Li, Y.  End-to-end semi-supervised learning for differentiable particle filters.  _arXiv preprint arXiv:2011.05748_ , 2020.
  * Williams (1992) Williams, R. J.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine Learning_ , 8(3-4):229–256, 1992.
  * Zhu et al. (2020) Zhu, M., Murphy, K., and Jonschkowski, R.  Towards differentiable resampling.  _arXiv preprint arXiv:2004.11938_ , 2020.

##  Appendix A Proof of Proposition 4.1

A particle filter with multinomial resampling is defined by the following
joint distribution

$\displaystyle\overline{q}_{\theta,\phi}(x_{1:T}^{1:N},a_{1:T-1}^{1:N})$
$\displaystyle=\prod_{i=1}^{N}q_{\phi}\left(x_{1}^{i}\right)\prod_{t=2}^{T}\prod_{i=1}^{N}w_{t-1}^{a_{t-1}^{i}}q_{\phi}\left(x_{t}^{i}|x{}_{t-1}^{a_{t-1}^{i}}\right)$

where $a_{t-1}^{i}\in\\{1,...,N\\}$ is the ancestral index of particle
$x_{t}^{i}$ and

$\displaystyle\omega_{\theta,\phi}(x_{1},y_{1})$
$\displaystyle=\frac{p_{\theta}(x_{1},y_{1})}{q_{\phi}(x_{1})},\quad\omega_{\theta,\phi}(x_{t-1},x_{t},y_{t})=\frac{p_{\theta}(x_{t},y_{t}|x_{t-1})}{q_{\phi}\left(x_{t}|x_{t-1}\right)}.$

Finally, we have
$w_{t}^{i}\propto\omega_{\theta,\phi}(x_{t-1}^{a_{t-1}^{i}},x_{t}^{i},y_{t}),~{}\sum_{i=1}^{N}w_{t}^{i}=1$.
We do not emphasize notationally that the weights $w_{t-1}^{a_{t-1}^{i}}$ are
$\theta,\phi$ and observations dependent.

The ELBO is given by

$\displaystyle\ell^{\text{ELBO}}(\theta,\phi)=$
$\displaystyle\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\log\widehat{p}_{\theta}(y_{1:T})\right]=\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\log\left(\frac{1}{N}\sum_{i=1}^{N}\omega_{\theta,\phi}(X_{1}^{i},y_{1})\right)+\sum_{t=2}^{T}\log\left(\frac{1}{N}\sum_{i=1}^{N}\omega_{\theta,\phi}(X_{t-1}^{A_{t-1}^{i}},X_{t}^{i},y_{t})\right)\right].$

We now compute $\nabla_{\theta}\ell^{\text{ELBO}}(\theta,\phi)$. We assume
from now on that the regularity conditions allowing us to swap the expectation
and differentiation operators are satisfied as in (Maddison et al., 2017; Le
et al., 2018; Naesseth et al., 2018). We can split the gradient using the
product rule and apply the log-derivative trick:

$\displaystyle\nabla_{\theta}\ell^{\text{ELBO}}(\theta,\phi)=$
$\displaystyle\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\nabla_{\theta}\log\widehat{p}_{\theta}(y_{1:T})\right]+\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\log\widehat{p}_{\theta}(y_{1:T})\nabla_{\theta}\log\overline{q}_{\theta,\phi}(X_{1:T}^{1:N},A_{1:T-1}^{1:N})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\nabla_{\theta}\log\left(\frac{1}{N}\sum_{i=1}^{N}\omega_{\theta,\phi}(X_{1}^{i},y_{1})\right)+\sum_{t=2}^{T}\nabla_{\theta}\log\left(\frac{1}{N}\sum_{i=1}^{N}\omega_{\theta,\phi}(X_{t-1}^{A_{t-1}^{i}},X_{t}^{i},y_{t})\right)\right]$
(18)
$\displaystyle+\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\log\widehat{p}_{\theta}(y_{1:T})\left\\{\sum_{t=2}^{T}\sum_{i=1}^{N}\nabla_{\theta}\log
w_{t-1}^{A_{t-1}^{i}}\right\\}\right]$ (19)

For the first part of the ELBO gradient (18), we have

$\displaystyle\nabla_{\theta}\log\left(\frac{1}{N}\sum_{i=1}^{N}\omega_{\theta,\phi}(X_{1}^{i}),y_{1}\right)$
$\displaystyle=\sum_{i=1}^{N}w_{1}^{i}\nabla_{\theta}\log
w_{\theta,\phi}(X_{1}^{i},y_{1})=\sum_{i=1}^{N}w_{1}^{i}\nabla_{\theta}\log
p_{\theta}(X_{1}^{i},y_{1})$

and

$\displaystyle\nabla_{\theta}\log\left(\frac{1}{N}\sum_{i=1}^{N}\omega_{\theta,\phi}(X_{t-1}^{A_{t-1}^{i}},X_{t}^{i},y_{t})\right)$
$\displaystyle=\sum_{i=1}^{N}w_{t}^{i}\nabla_{\theta}\log\omega_{\theta,\phi}(X_{t-1}^{A_{t-1}^{i}},X_{t}^{i},y_{t})=\sum_{i=1}^{N}w_{t}^{i}\nabla_{\theta}\log
p_{\theta}(X_{t}^{i},y_{t}|X_{t-1}^{A_{t-1}^{i}}).$

This gives

$\displaystyle\nabla_{\theta}\ell^{\text{ELBO}}(\theta,\phi)$
$\displaystyle=\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\sum_{i=1}^{N}w_{1}^{i}\nabla_{\theta}\log
p_{\theta}(X_{1}^{i},y_{1})+\sum_{t=2}^{T}\sum_{i=1}^{N}w_{t}^{i}\nabla_{\theta}\log
p_{\theta}(X_{t}^{i},y_{t}|X_{t-1}^{A_{t-1}^{i}})\right]$ (20)
$\displaystyle+\mathbb{E}_{\overline{q}_{\theta,\phi}}\left[\log\widehat{p}_{\theta}(y_{1:T})\left\\{\sum_{t=2}^{T}\sum_{i=1}^{N}\nabla_{\theta}\log
w_{t-1}^{A_{t-1}^{i}}\right\\}\right].$ (21)

When we ignore the gradient terms due to resampling corresponding to (21) as
proposed in (Naesseth et al., 2018; Le et al., 2018; Maddison et al., 2017;
Hirt & Dellaportas, 2019), we only use an unbiased estimate of the first term
(20), i.e.

$\hat{\nabla}_{\theta}\ell^{\textup{ELBO}}(\theta,\phi)\coloneqq\sum_{i=1}^{N}w_{1}^{i}\nabla_{\theta}\log
p_{\theta}(X_{1}^{i},y_{1})+\sum_{t=2}^{T}\sum_{i=1}^{N}w_{t}^{i}\nabla_{\theta}\log
p_{\theta}(X_{t}^{i},y_{t}|X_{t-1}^{A_{t-1}^{i}}),\quad\text{where~{}}(X^{1:N}_{1:T},A^{1:N}_{1:T-1})\sim\overline{q}_{\theta,\phi}(\cdot).$
(22)

Now we assume that the mild assumptions ensuring almost sure convergence of
the PF estimates are satisfied (see e.g. (Del Moral, 2004)). Under these
assumptions, the estimator (22) converges almost surely as
$N\rightarrow\infty$ towards

$\displaystyle\int$ $\displaystyle\nabla_{\theta}\log
p_{\theta}(x_{1},y_{1})p_{\theta}(x_{1}|y_{1})\mathrm{d}x_{1}+\sum_{t=2}^{T}\int\nabla_{\theta}\log
p_{\theta}(x_{t},y_{t}|x_{t-1})p_{\theta}(x_{t-1:t}|y_{1:t-1})\mathrm{d}x_{t-1:t}.$
(23)

Under an additional uniform integrability condition on
$\hat{\nabla}_{\theta}\ell^{\textup{ELBO}}(\theta,\phi)$, we thus have that
$\mathbb{E}_{\overline{q}_{\theta,\phi}}[\hat{\nabla}_{\theta}\ell^{\textup{ELBO}}(\theta,\phi)]$
converges towards (23). We recall that the true score is given by Fisher’s
identity and satisfies

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

This concludes the proof of Proposition 4.1.

##  Appendix B Notation and Assumptions

###  B.1 Filtering Notation

Recall $\mathcal{X}=\mathbb{R}^{d_{x}}$, denote the Borel sets of
$\mathcal{X}$ by $\mathcal{B}(\mathcal{X})$ and $\mathcal{P}(\mathcal{X})$ the
set of Borel probability measures on $(\mathcal{X},\mathcal{B}(\mathcal{X}))$.
In an abuse of notation, we shall use the same notation for a probability
measure and its density w.r.t. Lebesgue measure; i.e.
$\nu(\mathrm{d}x)=\nu(x)\mathrm{d}x$. We also use the standard notation
$\nu(\psi)=\int\psi(x)\nu(x)\mathrm{d}x$ for any test function $\psi$. In the
interest of notational clarity, we will remove subscript $\theta,\phi$ where
unnecessary in further workings.

We denote $\\{\alpha^{(t)}\\}_{t\geq 0}$ the predictive distributions where
$\alpha^{(t)}(x_{t})=p(x_{t}|y_{1:{t-1}})$ for $t>1$ and
$\alpha^{(1)}(x_{1})=\mu(x_{1})$ while $\\{\beta^{(t)}\\}_{t\geq 1}$ denotes
the filtering distributions; i.e. $\beta^{(t)}(x_{t})=p(x_{t}|y_{1:t})$ for
$t\geq 1$.

Using this notation, we have

$\displaystyle\alpha^{(t)}(\psi)=\int\psi(x_{t})f(x_{t}|x_{t-1})\beta^{(t-1)}(x_{t-1})\mathrm{d}x_{t-1}\mathrm{d}x_{t}\coloneqq\beta^{(t-1)}f(\psi),$
(24)
$\displaystyle\beta^{(t)}(\psi)=\frac{\alpha^{(t)}(g(y_{t}|\cdot)\psi)}{\alpha^{(t)}(g(y_{t}|\cdot))}=\frac{\beta^{(t-1)}(f(g(y_{t}|\cdot)\psi))}{\beta^{(t-1)}(f(g(y_{t}|\cdot)))}.$
(25)

More generally, for a proposal distribution $q(x_{t}|x_{t-1},y_{t})\neq
f(x_{t}|x_{t-1})$ with parameter $\phi\neq\theta$, the following recursion
holds

$\displaystyle\beta^{(t)}(\psi)=\frac{\beta^{(t-1)}(q(\omega_{t}\,\psi))}{\beta^{(t-1)}(q(\omega_{t}))}$
(26)

$\displaystyle\omega_{t}(x_{t-1},x_{t})\coloneqq\omega(x_{t-1},x_{t},y_{t})=\frac{g(y_{t}|x_{t})f(x_{t}|x_{t-1})}{q(x_{t}|x_{t-1},y_{t})}.$
(27)

To simplify the presentation, we will present the analysis in the scenario
where $\phi=\theta$ and $q(x_{t}|x_{t-1},y_{t})=f(x_{t}|x_{t-1})$ so we will
analyze (24) for which $\omega_{t}(x_{t-1},x_{t})=g(y_{t}|x_{t})$. In this
case, the particle approximations of $\mu$ is denoted $\mu_{N}$ and for $t>1$,
$\alpha^{(t)}$ and $\beta^{(t)}$ are given by the random measures

$\alpha^{(t)}_{N}(\psi)=\frac{1}{N}\sum_{i=1}^{N}\psi(X^{i}_{t}),\quad\beta_{N}^{(t)}(\psi)=\sum_{i=1}^{N}w^{i}_{t}\psi(X_{t}^{i}),\quad\tilde{\beta}_{N}^{(t)}(\psi)=\frac{1}{N}\sum_{i=1}^{N}\psi(\tilde{X}_{t}^{i}),$
(28)

where $w^{i}_{t}\propto g(y_{t}|X_{t}^{i})$ with $\sum_{i=1}^{N}w^{i}_{t}=1$
and particles are drawn from $X^{i}_{t}\sim f(\cdot|\tilde{X}^{i}_{t-1})$.

Here $\beta_{N}^{(t)}$ denotes the weighted particle approximation of
$\beta^{(t)}$ while $\tilde{\beta}_{N}^{(t)}$ is the uniformly weighted
approximation obtained after the DET transformation described in Section 3.2.

###  B.2 Optimal Transport Notation

Recall from Section 2.1, $\mathcal{P}_{t}^{\text{OT}}$ denotes a transport
between $\alpha^{(t)}$ and $\beta^{(t)}$ with accompanying map
$\mathbf{T}^{(t)}$. $\mathcal{P}_{t}^{\text{OT},N}$ denotes an optimal
transport between particle approximations $\alpha_{N}^{(t)}$ and
$\beta_{N}^{(t)}$ with corresponding transport matrix,
$\mathbf{P}^{\text{OT}}$ with $i,j$ entry $p^{\text{OT}}_{i,j}$. To simplify
notation, we remove script $t$ when not needed.

Similarly from Section 3.1, $\mathcal{P}^{\text{OT,N}}_{\epsilon}$ denotes the
regularized transport between $\alpha^{(t)}_{N}$ and $\beta^{(t)}_{N}$ with
accompanying matrix $\mathbf{P}^{\text{OT}}_{\epsilon}$ with $i,j$ entry
$p^{\text{OT}}_{\epsilon,i,j}$. Recall
$\tilde{\beta}^{(t)}_{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{\tilde{X}^{i}}$ is
the uniformly weighted particle approximation for $\beta^{(t)}$ under the DET,
i.e. $\tilde{X}^{i}=\mathbf{T}^{(t)}_{N,\epsilon}(X^{i})=\int
y\mathcal{P}^{\mathrm{OT},N}_{\epsilon}(\mathrm{d}y|x^{i})$. Note that
$\tilde{X}^{i}_{N,\epsilon}$ will be used where necessary to avoid ambiguity
when comparing to other resampling schemes.

Recall also for $p>0$:

$\displaystyle\mathcal{W}_{p}^{p}(\alpha,\beta)$
$\displaystyle=\min_{\mathcal{P}\in\mathcal{U}(\alpha,\beta)}\mathbb{E}_{(U,V)\sim\mathcal{P}}\big{[}||U-V||^{p}\big{]}$
(29)

where $\mathcal{U}(\alpha,\beta)$ is the collection of couplings with
marginals $\alpha$ and $\beta$.

###  B.3 Assumptions

Our results will rely on the following four assumptions.

######  Assumption B.1.

$\mathcal{X}\subset\mathbb{R}^{d}$ is a compact subset with diameter

$\mathfrak{d}\coloneqq\sup_{x,y\in\mathcal{X}}|x-y|.$

######  Assumption B.2.

There exists $\kappa\in(0,1)$ such that for any two probability measures
$\pi,\rho$ on $\mathcal{X}$

$\mathcal{W}_{k}(\pi f,\rho f)\leq\kappa\mathcal{W}_{k}(\pi,\rho),\qquad
k=1,2.$

######  Assumption B.3.

The weight function $\omega^{(t)}:\mathcal{X}\to[\Delta,\Delta^{-1}]$ is
1-Lipschitz for all $t$.

######  Assumption B.4.

There exists a $\lambda>0$, such that for all $t\geq 0$ the unique optimal
transport plan between $\alpha^{(t)}$ and $\beta^{(t)}$ is given by a
deterministic, $\lambda$-Lipschitz map $\mathbf{T}^{(t)}$.

##  Appendix C Auxiliary Results and Proof of Proposition 4.2

We start by establishing a couple of key auxiliary results which will be then
used subsequently to establish Proposition 4.2.

###  C.1 Auxiliary Results

As per section 2.1, let $\mathcal{S}(\alpha_{N},\beta_{N})$ denote the
collection of coupling matrices between
$\alpha_{N}=\sum_{i=1}^{N}a_{i}\delta_{Y^{i}}$ with $a_{i}>0$ and
$\beta_{N}=\sum_{i=1}^{N}b_{i}\delta_{X^{i}}$. We also denote entropy by $H$
where $H(\mathbf{P})=\sum_{i,j}p_{i,j}\log(1/p_{i,j})$ for
$\mathbf{P}=(p_{i,j})_{i,j}\in\mathcal{S}(\alpha_{N},\beta_{N})$.

######  Lemma C.1.

The entropic radius, $R_{H}$, of simplex $\mathcal{U}(\alpha_{N},\beta_{N})$
may be bounded above as follows

$R_{H}\coloneqq\max_{\mathbf{P}_{1},\mathbf{P}_{2}\in\mathcal{S}(\alpha_{N},\beta_{N})}H(\mathbf{P}_{1})-H(\mathbf{P}_{2})\leq
2\log(N)$

###### Proof.

Notice that $-H(\mathbf{P})$ is convex, so that $H(\mathbf{P})$ is concave.

$\displaystyle\sum_{i,j}p_{i,j}\log\left(\frac{1}{p_{i,j}}\right)$
$\displaystyle=N^{2}\sum_{i,j}\frac{1}{N^{2}}p_{i,j}\log\left(\frac{1}{p_{i,j}}\right)$
$\displaystyle\leq
N^{2}H\left(\frac{1}{N^{2}}\sum_{i,j}p_{i,j}\right)=N^{2}H(1/N^{2})=N^{2}\frac{1}{N^{2}}\log(N^{2})=2\log(N).$

In addition since $p_{i,j}\leq 1$ for all $i,j$, we have that
$H(\mathbf{P})\geq 0$ and therefore we can bound

$R_{H}=\max_{P_{1},P_{2}\in\mathcal{P}}H(\mathbf{P}_{1})-H(\mathbf{P}_{2})\leq\max_{\mathbf{P}_{1}\in\mathcal{S}(\alpha_{N},\beta_{N})}H(\mathbf{P}_{1})\leq
2\log(N).$

∎

######  Lemma C.2.

Let $\mathcal{X}\subset\mathbb{R}^{d}$ be compact with diameter
$\mathfrak{d}>0$. Suppose we are given two probability measures $\alpha,\beta$
on $\mathcal{X}$ with a unique deterministic, $\lambda$-Lipschitz optimal
transport map $\mathbf{T}$ while
$\alpha_{N}=\sum_{i=1}^{N}a_{i}\delta_{Y^{i}}$ with $a_{i}>0$ and
$\beta_{N}=\sum_{i=1}^{N}b_{i}\delta_{X^{i}}$. We write
$\mathcal{P}^{\mathrm{OT},N}$, resp. $\mathcal{P}^{\mathrm{OT},N}_{\epsilon}$,
for an optimal coupling between $\alpha_{N}$ and $\beta_{N}$, resp. the
$\epsilon$-regularized optimal transport plan, between $\alpha_{N}$ and
$\beta_{N}$. Then

$\displaystyle\left[\int||y-\mathbf{T}(x)||^{2}\mathcal{P}^{\mathrm{OT},N}_{\epsilon}(\mathrm{d}x,\mathrm{d}y)\right]^{\frac{1}{2}}\leq
2\lambda^{1/2}\mathcal{E}^{1/2}\left[\mathfrak{d}^{1/2}+\mathcal{E}\right]^{1/2}+\max\\{\lambda,1\\}\left[\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)\right],$

where

$\mathcal{E}\coloneqq\mathcal{E}(N,\epsilon,\alpha,\beta)\coloneqq\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)+\sqrt{2\epsilon\log(N)}.$

###### Proof.

From Corollary 3.8 from (Li & Nochetto, 2021)

$\displaystyle\left[\int\|\mathbf{T}(x)-y\|^{2}\mathcal{P}^{\mathrm{OT},N}_{\epsilon}(\mathrm{d}x,\mathrm{d}y)\right]^{1/2}$
$\displaystyle\leq
2\lambda^{1/2}\sqrt{\tilde{e}_{N,\epsilon}}\left[\mathcal{W}_{2}(\alpha,\beta)+\tilde{e}_{N,\epsilon}\right]^{1/2}+\lambda\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta),$

where $\lambda$ is the Lipschitz constant of the optimal transport map
$\mathbf{T}$ sending $\alpha$ to $\beta$, and

$\tilde{e}_{N,\epsilon}\coloneqq\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)+\left[\int\|x-y\|^{2}\mathcal{P}^{\mathrm{OT},N}_{\epsilon}(\mathrm{d}x,\mathrm{d}y)\right]^{1/2}-\mathcal{W}_{2}(\alpha_{N},\beta_{N}).$
(30)

From Proposition 4 of (Weed, 2018),

$\sum_{i,j=1,\dots,N}p^{\mathrm{OT}}_{\epsilon,i,j}|Y_{i}-X_{j}|^{2}-\mathcal{W}^{2}_{2}(\alpha_{N},\beta_{N})\leq\epsilon
R_{H},$

where $R_{H}$ is the entropic radius as defined in Lemma C.1.

By Lemma C.1 we therefore have that

$\int\|x-y\|^{2}\mathcal{P}^{\mathrm{OT},N}_{\epsilon}(\mathrm{d}x,\mathrm{d}y)-\mathcal{W}^{2}_{2}(\alpha_{N},\beta_{N})\leq
2\epsilon\log(N).$

Since $x\mapsto\sqrt{x}$ is sub-additive, for $r,s>0$ we have that
$\sqrt{r}-\sqrt{s}\leq\sqrt{r-s}$, whence

$\left[\int\|x-y\|^{2}\mathcal{P}^{\mathrm{OT},N}_{\epsilon}(\mathrm{d}x,\mathrm{d}y)\right]^{1/2}-\mathcal{W}_{2}(\alpha_{N},\beta_{N})\leq\sqrt{2\epsilon\log
N}.$

We thus have

$\tilde{e}_{N,\epsilon}\leq\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)+\sqrt{2\epsilon\log(N)}.$

In addition, by Assumption B.1 we have that
$\mathcal{W}_{2}(\alpha,\beta)\leq\mathfrak{d}^{1/2}$ and the result follows.
∎

###  C.2 Proof of Proposition 4.2

###### Proof of Proposition 4.2.

By definition, we have
$\tilde{\beta}_{N}(\mathrm{d}\tilde{x})=\int\alpha_{N}(\mathrm{d}x)\delta_{\mathbf{T}_{N,\epsilon}(x)}(\mathrm{d}\tilde{x})$
with
$\mathbf{T}_{N,\epsilon}(x)\coloneqq\int\tilde{x}\mathcal{P}^{\mathrm{OT,N}}_{\epsilon}(\mathrm{d}\tilde{x}|x)$
while, as $\mathcal{P}^{\mathrm{OT},N}_{\epsilon}$ belongs to
$\mathcal{U}(\alpha_{N},\beta_{N})$, we also have
$\beta_{N}(\mathrm{d}\tilde{x})=\int\alpha_{N}(\mathrm{d}x)\mathcal{P}^{\mathrm{OT,N}}_{\epsilon}(\mathrm{d}\tilde{x}|x)$.
We then have for any 1-Lipschitz function

$\displaystyle\left|\beta_{N}(\psi)-\tilde{\beta}_{N}(\psi)\right|$
$\displaystyle=\left|\int\left[\int(\psi(\tilde{x})-\psi(\mathbf{T}_{N,\epsilon}(x)))\mathcal{P}^{\mathrm{OT,N}}_{\epsilon}(\mathrm{d}\tilde{x}|x)\right]\alpha_{N}(\mathrm{d}x)\right|$
$\displaystyle\leq\iint\left|\psi(\tilde{x})-\psi(\mathbf{T}_{N,\epsilon}(x))\right|\alpha_{N}(\mathrm{d}x)\mathcal{P}^{\mathrm{OT,N}}_{\epsilon}(\mathrm{d}\tilde{x}|x)$
$\displaystyle\leq\iint||\tilde{x}-\mathbf{T}_{N,\epsilon}(x)||\mathcal{P}^{\mathrm{OT,N}}_{\epsilon}(\mathrm{d}x,\mathrm{d}\tilde{x})$
$\displaystyle\leq\left(\iint||\tilde{x}-\mathbf{T}_{N,\epsilon}(x)||^{2}\mathcal{P}^{\mathrm{OT,N}}_{\epsilon}(\mathrm{d}x,\mathrm{d}\tilde{x})\right)^{\frac{1}{2}}$
$\displaystyle\leq\left(\iint||\tilde{x}-\mathbf{T}(x)||^{2}\mathcal{P}^{\mathrm{OT,N}}_{\epsilon}(\mathrm{d}x,\mathrm{d}\tilde{x})\right)^{\frac{1}{2}},$

where the final inequality follows from the fact that for any random vector
$V$ the mapping $v\mapsto\mathbb{E}[\|V-v\|^{2}]$ is minimized at
$v=\mathbb{E}[V]$. The stated result is then obtained using Lemma C.2. ∎

##  Appendix D Proof of Proposition 4.3

For technical reasons, we analyse here a slightly modified PF algorithm where

$\alpha_{N}^{(t)}=\frac{1}{N}\sum_{j=1}^{N}\delta_{X^{j}_{t}},\qquad
X^{j}_{t}\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\tilde{\beta}_{N}^{(t-1)}f=\frac{1}{N}\sum_{j=1}^{N}f\left(\cdot\Big{|}\tilde{X}^{j}_{t-1}\right).$
(31)

instead of the standard version where one has

$\alpha_{N}^{(t)}=\frac{1}{N}\sum_{j=1}^{N}\delta_{X^{j}_{t}},\qquad
X^{j}_{t}\sim f\left(\cdot\Big{|}\tilde{X}^{j}_{t-1}\right).$

This slightly modified version of the bootstrap PF was analyzed for example in
(Del Moral & Guionnet, 2001). The analysis does capture the additional error
arising from the use of DET instead of resampling. Similar results should hold
for the standard PF algorithm. The main technical reason for analysing this
modified algorithm is our reliance on Theorem 2 of (Fournier & Guillin, 2015);
analysing the standard PF algorithm requires a version of (Fournier & Guillin,
2015) for stratified sampling and will be done in future work.

######  Proposition D.1.

Suppose that Assumptions B.1, B.2 and B.3 hold. Suppose also that given
$\tilde{\beta}_{N}^{(t-1)}$, $\alpha_{N}^{(t)}$ is defined through (31).
Define the functions

$\displaystyle\mathcal{F}(x)$ $\displaystyle\coloneqq
x+\sqrt{\mathfrak{d}K_{1}(\Delta,\mathfrak{d})x}$ $\displaystyle f_{d}(x)$
$\displaystyle\coloneqq\begin{cases}x,&d<4\\\ \frac{x}{\log(2+1/x)},&d=4\\\
x^{d/2},&d>4.\end{cases}$
$\displaystyle\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)$
$\displaystyle\coloneqq\mathcal{F}\left(\kappa
x+\sqrt{f_{d}^{-1}\left(\frac{\log(C/\delta)}{cN}\right)}\right),$
$\displaystyle\frac{1}{\mathfrak{d}}\mathfrak{G}^{2}_{\epsilon,\delta,N,d}(x)$
$\displaystyle\coloneqq
2\lambda^{1/2}\left[\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)+\sqrt{2\epsilon\log
N}\right]^{1/2}\left[\mathfrak{d}^{1/2}+\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)+\sqrt{2\epsilon\log
N}\right]^{1/2}$
$\displaystyle\qquad+\lambda\kappa\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)+\max\\{\lambda,1\\}\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right).$
(32)

Then for any $\epsilon,\delta>0$ we have with probability at least $1-\delta$,
over the sampling step in (31), that

$\mathcal{W}_{2}\left(\tilde{\beta}_{N}^{(t)},\beta^{(t)}\right)\leq\mathfrak{G}_{\epsilon,\delta,N,d}\left[\mathcal{W}_{2}\left(\tilde{\beta}_{N}^{(t-1)},\beta^{(t-1)}\right)\right]$
(33)

In particular if $\mathcal{W}_{2}(\tilde{\beta}_{N}^{(t-1)},\beta^{(t-1)})\to
0$ and $\epsilon_{N}=o(1/\log(N))$ as $N\to\infty$ we have that

$\mathcal{W}_{2}\big{(}\tilde{\beta}_{N}^{(t)},\beta^{(t)}\big{)}\to 0,$

in probability.

###### Proof of Proposition D.1.

To keep notation concise we write for $N\geq 1$

$\alpha_{N}\coloneqq\alpha_{N}^{(t)},\quad\alpha_{N}^{\prime}\coloneqq\tilde{\beta}_{N}^{(t-1)}f,\quad\beta_{N}\coloneqq\beta_{N}^{(t)},\quad\tilde{\beta}_{N}\coloneqq\tilde{\beta}_{N}^{(t-1)}.$

Controlling $\mathcal{W}_{1}(\beta_{N},\beta)$. Let $\psi$ be 1-Lipschitz.
Without loss of generality we may assume that $\psi(0)=0$ since otherwise we
can remove a constant.

$\displaystyle\left|\beta_{N}(\psi)-\beta(\psi)\right|$
$\displaystyle=\left|\frac{\alpha_{N}(\omega\psi)}{\alpha_{N}(\omega)}-\frac{\alpha(\omega\psi)}{\alpha(\omega)}\right|$
$\displaystyle\leq\left|\frac{\alpha_{N}(\omega\psi)}{\alpha_{N}(\omega)}-\frac{\alpha(\omega\psi)}{\alpha_{N}(\omega)}\right|+\left|\frac{\alpha(\omega\psi)}{\alpha_{N}(\omega)}-\frac{\alpha(\omega\psi)}{\alpha(\omega)}\right|$
$\displaystyle\leq\Delta^{-1}\left|\alpha_{N}(\omega\psi)-\alpha(\omega\psi)\right|+\Delta^{-2}\alpha(\omega\psi)|\alpha_{N}(\omega)-\alpha(\omega)|.$

At this stage notice that

$|(\omega\psi)^{\prime}|\leq|\omega^{\prime}\psi|+|\omega\psi^{\prime}|\leq\|\psi\|_{\infty}+\|\omega\|_{\infty}.$

Notice that

$|\psi(x)|=|\psi(x)-\psi(0)|\leq|x-0|\leq\mathfrak{d}.$

Therefore we have that

$|(\omega\psi)^{\prime}|\leq\mathfrak{d}+\Delta^{-1},$

and thus $\omega\psi$ is $(\mathfrak{d}+\Delta^{-1})$-Lipschitz. It follows
that

$\displaystyle\left|\beta_{N}(\psi)-\beta(\psi)\right|$
$\displaystyle\leq\Delta^{-1}\left|\alpha_{N}(\omega\psi)-\alpha(\omega\psi)\right|+\Delta^{-2}\alpha(\omega\psi)|\alpha_{N}(\omega)-\alpha(\omega)|$
$\displaystyle\leq\Delta^{-1}(\mathfrak{d}+\Delta^{-1})\mathcal{W}_{1}(\alpha_{N},\alpha)+\Delta^{-3}\mathfrak{d}\mathcal{W}_{1}(\alpha_{N},\alpha)$
$\displaystyle=:K_{1}(\Delta,\mathfrak{d})\mathcal{W}_{1}(\alpha_{N},\alpha).$

Therefore we have that

$\mathcal{W}_{1}(\beta_{N},\beta)\leq
K_{1}(\Delta,\mathfrak{d})\mathcal{W}_{1}(\alpha_{N},\alpha).$ (34)

Notice that using the compactness of the state space we easily get also that

$\mathcal{W}_{2}(\beta_{N},\beta)\leq\sqrt{\mathfrak{d}\mathcal{W}_{1}(\beta_{N},\beta)}\leq\sqrt{\mathfrak{d}K_{1}(\Delta,\mathfrak{d})\mathcal{W}_{1}(\alpha_{N},\alpha)}\leq\sqrt{\mathfrak{d}K_{1}(\Delta,\mathfrak{d})\mathcal{W}_{2}(\alpha_{N},\alpha)},$
(35)

since clearly $\mathcal{W}_{1}(\rho,\sigma)\leq\mathcal{W}_{2}(\rho,\sigma)$
for any two probability measures $\rho,\sigma$.

Controlling $\mathcal{W}_{1}(\tilde{\beta}_{N,\epsilon},{\beta})$. Again
supposing $\psi$ is 1-Lipschitz, and $\psi(0)=0$, consider

$\displaystyle\left|\tilde{\beta}_{N}(\psi)-\tilde{\beta}(\psi)\right|$
$\displaystyle=\left|\int\psi(\mathbf{T}_{N,\epsilon}(x))\alpha_{N}(\mathrm{d}x)-\int\psi(\mathbf{T}(x))\alpha(\mathrm{d}x)\right|$
$\displaystyle\leq\left|\int\psi(\mathbf{T}_{N,\epsilon}(x))\alpha_{N}(\mathrm{d}x)-\int\psi(\mathbf{T}(x))\alpha_{N}(\mathrm{d}x)\right|$
$\displaystyle\quad+\left|\int\psi(\mathbf{T}(x))\alpha_{N}(\mathrm{d}x)-\int\psi(\mathbf{T}(x))\alpha(\mathrm{d}x)\right|$

For the second term, using the fact that $\mathbf{T}$ and $\psi$ are
$\lambda$\- and 1-Lipschitz respectively, we have that $\psi\circ\mathbf{T}$
is $\lambda$-Lipschitz and therefore

$\displaystyle\left|\int\psi(\mathbf{T}(x))\alpha_{N}(\mathrm{d}x)-\int\psi(\mathbf{T}(x))\alpha(\mathrm{d}x)\right|$
$\displaystyle\leq\lambda\mathcal{W}_{1}(\alpha_{N},\alpha)\leq\lambda\mathcal{W}_{2}(\alpha_{N},\alpha),$

where we used Assumption B.2 for that last inequality For the first term
recall that using Cauchy-Schwarz and Jensen we get

$\displaystyle\left|\int\psi(\mathbf{T}_{N,\epsilon}(x))\alpha_{N}(\mathrm{d}x)-\int\psi(\mathbf{T}(x))\alpha_{N}(\mathrm{d}x)\right|$
$\displaystyle\leq\int\left|\mathbf{T}_{N,\epsilon}(x)-\mathbf{T}(x)\right|\alpha_{N}(\mathrm{d}x)$
$\displaystyle\leq\int\left|\int
y\mathcal{P}_{N,\epsilon}(x,\mathrm{d}y)-\mathbf{T}(x)\right|\alpha_{N}(\mathrm{d}x)$
$\displaystyle\leq\iint\left|y-\mathbf{T}(x)\right|\alpha_{N}(\mathrm{d}x)\mathcal{P}_{N,\epsilon}(x,\mathrm{d}y)$
$\displaystyle\leq\left[\iint\left|y-\mathbf{T}(x)\right|^{2}\alpha_{N}(\mathrm{d}x)\mathcal{P}_{N,\epsilon}(x,\mathrm{d}y)\right]^{1/2}.$

Here we can directly apply Lemma C.2 to obtain

$\displaystyle\left[\iint\left|y-\mathbf{T}(x)\right|^{2}\alpha_{N}(\mathrm{d}x)\mathcal{P}_{N,\epsilon}(x,\mathrm{d}y)\right]^{1/2}$
$\displaystyle\leq
2\lambda^{1/2}\mathcal{E}^{1/2}\left[\mathfrak{d}^{1/2}+\mathcal{E}\right]^{1/2}+\max\\{\lambda,1\\}\left[\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)\right],$

where

$\mathcal{E}\coloneqq\mathcal{E}(n,\epsilon,\alpha,\beta)\coloneqq\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)+\sqrt{2\epsilon\log(N)}.$

From (35) we have that

$\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)\leq\mathcal{W}_{2}(\alpha_{N},\alpha)+\sqrt{\mathfrak{d}K_{1}(\Delta,\mathfrak{d})\mathcal{W}_{2}(\alpha_{N},\alpha)}.$

Next we want to bound $\mathcal{W}_{2}(\alpha_{N},\alpha)$. Notice first that

$\mathcal{W}_{2}(\alpha_{N},\alpha)\leq\mathcal{W}_{2}(\alpha_{N},\alpha_{N}^{\prime})+\mathcal{W}_{2}(\alpha_{N}^{\prime},\alpha)\leq\mathcal{W}_{2}(\alpha_{N},\alpha_{N}^{\prime})+\kappa\mathcal{W}_{2}\left(\tilde{\beta}_{N}^{(t-1)},\beta^{(t-1)}\right),$

by Assumption B.2.

To control the other term we use (Fournier & Guillin, 2015) to obtain a high
probability bound on $\mathcal{W}_{2}(\alpha_{N},\alpha_{N}^{\prime})$. In
particular, using Theorem 2 from (Fournier & Guillin, 2015), with
$\alpha=\infty$ since we are in a compact domain, that for some positive
constants $C,c$ we have

$\mathbb{P}\left[\mathcal{W}_{2}^{2}(\alpha_{N},\alpha_{N}^{\prime})\geq
x\right]\leq C\exp\left[-cNf^{2}_{d}(x)\right],$ (36)

where

$f_{d}(x)\coloneqq\begin{cases}x,&d<4\\\ \frac{x}{\log(2+1/x)},&d=4\\\
x^{d/2},&d>4.\end{cases}$ (37)

In particular, for any $\delta>0$, with probability at least $1-\delta$ over
the sampling step in $F_{N}$ we have that

$\mathcal{W}_{2}(\alpha_{N},\alpha_{N}^{\prime})\leq\sqrt{f_{d}^{-1}\left(\frac{\log(C/\delta)}{cN}\right)}.$
(38)

Assuming that $d\geq 4$ the rate then is of order $N^{-1/d}$ as expected.

Therefore with probability at least $1-\delta$ over the sampling step we have
that

$\mathcal{W}_{2}(\alpha_{N},\alpha)+\mathcal{W}_{2}(\beta_{N},\beta)\leq\mathcal{F}_{N,\epsilon,\delta,d}\left(\mathcal{W}_{2}\left(\tilde{\beta}_{N}^{(t-1)},\beta^{(t-1)}\right)\right),$

where

$\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)=\mathcal{F}\left(\kappa
x+\sqrt{f_{d}^{-1}\left(\frac{\log(C/\delta)}{cN}\right)}\right),\qquad\mathcal{F}(x)\coloneqq
x+\sqrt{\mathfrak{d}K_{1}(\Delta,\mathfrak{d})x}$ (39)

Thus overall we have with probability at least $1-\delta$ over the sample

$\displaystyle\mathcal{W}_{2}(\tilde{\beta}_{N,\epsilon},\tilde{\beta})$
$\displaystyle\leq\sqrt{\mathfrak{d}\mathcal{W}_{1}(\tilde{\beta}_{N,\epsilon},\tilde{\beta})}\leq\mathfrak{G}_{\epsilon,\delta,N,d}\left(\mathcal{W}_{2}\left(\tilde{\beta}_{N}^{(t-1)},\beta^{(t-1)}\right)\right),$

where

$\displaystyle\frac{1}{\mathfrak{d}}\mathfrak{G}^{2}_{\epsilon,\delta,N,d}(x)$
$\displaystyle\coloneqq
2\lambda^{1/2}\left[\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)+\sqrt{2\epsilon\log
N}\right]^{1/2}\left[\mathfrak{d}^{1/2}+\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)+\sqrt{2\epsilon\log
N}\right]^{1/2}$
$\displaystyle\qquad+\lambda\kappa\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right)+\max\\{\lambda,1\\}\mathcal{F}_{N,\epsilon,\delta,d}\left(x\right).$

In particular notice that if we set $\epsilon_{N}=o(1/\log N)$ and
$x_{N}=o(1)$ we have

$\mathfrak{G}_{\epsilon_{N},\delta,N,d}(x_{N})\to 0.$

Therefore, notice that if $\epsilon_{N}=o(1/\log N)$ and
$\mathcal{W}_{2}(\mu_{N},\mu)\to 0$, then for any $x>0$ we have that

$\mathbb{P}\left[\mathcal{W}_{2}\big{(}\tilde{\beta}_{N,\epsilon},\tilde{\beta}\big{)}\geq
x\right]\leq\mathbb{P}[\mathcal{W}_{2}(\alpha_{N}^{\prime},\alpha_{N})\geq
x^{\prime}],$

for some $x^{\prime}$ that does not depend on $N$, where the probability is
over the sampling step. The convergence in probability follows. ∎

######  Proposition D.2.

Let $\mu_{N}=\tfrac{1}{N}\sum_{i=1}^{N}\delta_{X_{1}^{i}}$ where
$X^{i}_{1}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\mu\coloneqq
q(\cdot|y_{1})$ for $i\in[N]$ and suppose that for $t\geq 1$,
$\alpha_{N}^{(t)}$ is defined through (31). Under Assumptions B.1, B.2, B.3
and B.4, for any $\delta>0$, with probability at least $1-2\delta$ over the
sampling steps, for any bounded 1-Lipschitz $\psi$, for any $t\in[1:T]$, the
approximations of the filtering distributions and log-likelihood computed by
DPF satisfy

$\displaystyle|\tilde{\beta}^{(t)}_{N}(\psi)-\beta^{(t)}(\psi)|$
$\displaystyle\leq\mathfrak{G}_{\epsilon,\delta/T,N,d}^{(t)}\left(\sqrt{f_{d}^{-1}\left(\frac{\log(CT/\delta)}{cN}\right)}\right)$
(40) $\displaystyle\left|\log\frac{\hat{p}_{N}(y_{1:T})}{p(y_{1:T})}\right|$
$\displaystyle\leq\frac{\kappa}{\Delta}\max_{t\in[1:T]}\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]\sum_{t=1}^{T}\mathfrak{G}_{\epsilon,\delta/T,N,d}^{(t)}\left(\sqrt{f_{d}^{-1}\left(\frac{\log(CT/\delta)}{cN}\right)}\right)$
(41)

where $C$ is a finite constant independent of $T$,
$\mathfrak{G}_{\epsilon,\delta/T,N,d},f_{d}$ are defined in (32), and
$\mathrm{Lip}[f]$ is the Lipschitz constant of the function $f$.
$\mathfrak{G}^{(t)}_{\epsilon,\delta/T,N,d}$ denotes the $t$-repeated
composition of function $\mathfrak{G}_{\epsilon,\delta/T,N,d}$. In particular,
if we set $\epsilon_{N}=o(1/\log N)$

$\left|\log\frac{\hat{p}_{N}(y_{1:T})}{p(y_{1:T})}\right|\to 0,$

in probability.

###### Proof of Proposition D.2.

Following the proof of Proposition D.1, we define
${\alpha_{N}^{(t)}}^{\prime}=\tilde{\beta}_{N}^{(t-1)}f$ and for $t\in[1:T]$,
the events

$A_{t}\coloneqq\mathcal{W}_{2}\left(\alpha_{N}^{(t)},{\alpha_{N}^{(t)}}^{\prime}\right)\leq\sqrt{f_{d}^{-1}\left(\frac{\log(CT/\delta)}{cN}\right)}.$

We know from Theorem 2 in (Fournier & Guillin, 2015) that
$\mathbb{P}(A_{t})\geq 1-\delta/T$, where the probability is over the sampling
step. In particular we have that

$\mathbb{P}\left[\bigcap_{t=1}^{T}A_{t}\right]=1-\mathbb{P}\left[\bigcup_{t=1}^{T}A_{t}^{\mathtt{C}}\right]\geq
1-\sum_{t=1}^{N}\mathbb{P}\left[A_{t}^{\mathtt{C}}\right]\geq
1-T\frac{\delta}{T}=1-\delta.$

Notice that on the event $\cap_{t=1}^{T}A_{t}$, iterating the bound (33) we
have

$\displaystyle\mathcal{W}_{2}\left(\tilde{\beta}_{N}^{(t)},\beta^{(t)}\right)$
$\displaystyle\leq\mathfrak{G}_{\epsilon,\delta/T,N,d}^{(t)}\left(\mathcal{W}_{2}\left(\mu_{N},\mu\right)\right),$

with probability at least $1-\delta$. Again by Theorem 2 in (Fournier &
Guillin, 2015) we have that with probability at least $1-\delta$

$\mathcal{W}_{2}(\mu_{N},\mu)\leq\sqrt{f_{d}^{-1}\left(\frac{\log(CT/\delta)}{cN}\right)}.$

Therefore with probability at least $1-2\delta$ we have

$\mathfrak{G}_{\epsilon,\delta/T,N,d}^{(t)}\left(\sqrt{f_{d}^{-1}\left(\frac{\log(CT/\delta)}{cN}\right)}\right).$

It remains to prove (41). Note that
$|\log(x)-\log(y)|\leq\frac{|x-y|}{\min\\{x,y\\}}$ for any $x,y>0$ so

$\displaystyle\big{|}\log\hat{p}(y_{1:T})-\log p(y_{1:T})\big{|}$
$\displaystyle\leq\sum_{t=1}^{T}\big{|}\log\hat{p}(y_{t}|y_{1:t-1})-\log
p(y_{t}|y_{1:t-1})\big{|}$
$\displaystyle\leq\sum_{t=1}^{T}\big{|}\frac{\hat{p}(y_{t}|y_{1:t-1})-p(y_{t}|y_{1:t-1})}{\min(\hat{p}(y_{t}|y_{1:t-1}),p(y_{t}|y_{1:t-1}))}\big{|}$
$\displaystyle\leq\Delta^{-1}\sum_{t=1}^{T}|\hat{p}(y_{t}|y_{1:t-1})-p(y_{t}|y_{1:t-1})\big{|}$
(42)

where $\Delta$ is defined in Assumption B.3.

The term in line (42) may be written as follows

$\displaystyle\hat{p}(y_{t}|y_{1:t-1})-p(y_{t}|y_{1:t-1})$ $\displaystyle=$
$\displaystyle\iint
g(y_{t}|x_{t})f(\mathrm{d}x_{t}|\tilde{x}_{t-1})\tilde{\beta}^{(t-1)}_{N}(\mathrm{d}\tilde{x}_{t-1})-\iint
g(y_{t}|x_{t})f(\mathrm{d}x_{t}|\tilde{x}_{t-1})\tilde{\beta}^{(t-1)}(\mathrm{d}\tilde{x}_{t-1})$
$\displaystyle=$ $\displaystyle\tilde{\beta}^{(t-1)}_{N}(h)-\beta^{(t-1)}(h)$

for $\Delta^{2}\leq h(x):=\int
g(y_{t}|x^{\prime})f(x^{\prime}|x)\mathrm{d}x^{\prime}\leq\Delta^{-2}$. At
this point notice also that

$\displaystyle h(x)-h(x^{\prime})$ $\displaystyle=\int
f(\mathrm{d}w|x)g(y_{t}\mid w)-\int f(\mathrm{d}w|x^{\prime})g(y_{t}\mid w)$
$\displaystyle=\int\delta_{x}(\mathrm{d}z)\int f(\mathrm{d}w|z)g(y_{t}\mid
w)-\int\delta_{x^{\prime}}(\mathrm{d}z)\int f(\mathrm{d}w|z)g(y_{t}\mid w)$
$\displaystyle=[\delta_{x}f][g(y_{t}\mid\cdot)]-[\delta_{x^{\prime}}f][g(y_{t}\mid\cdot)]$
$\displaystyle\leq\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]\mathcal{W}_{1}\left(\delta_{x}f,\delta_{x^{\prime}}f\right)\leq\kappa\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]\mathcal{W}_{1}\left(\delta_{x},\delta_{x^{\prime}}\right)=\kappa\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]|x-x^{\prime}|,$

by Assumption B.2. It follows therefore that $h$ is Lipschitz and therefore
that

$\displaystyle\hat{p}(y_{t}|y_{1:t-1})-p(y_{t}|y_{1:t-1})$
$\displaystyle=\tilde{\beta}^{(t-1)}_{N}(h)-\beta^{(t-1)}(h)\leq\kappa\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]\mathcal{W}_{1}\left(\beta^{(t-1)}_{N},\beta^{(t-1)}\right).$

Combining (40) and (42), and using the fact that
$\mathcal{W}_{1}\leq\mathcal{W}_{2}$, we thus get

$\displaystyle\big{|}\log\hat{p}(y_{1:T})-\log p(y_{1:T})\big{|}$
$\displaystyle\leq\Delta^{-1}\kappa\sum_{t=1}^{T}\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]\mathcal{W}_{1}\left(\tilde{\beta}_{N}^{(t-1)},\beta^{(t)}\right)$
$\displaystyle\leq\Delta^{-1}\kappa\max_{t\in[1:T]}\mathrm{Lip}\left[g(y_{t}\mid\cdot)\right]\sum_{t=1}^{T}\mathfrak{G}_{\epsilon,\delta/T,N,d}^{(t)}\left(\sqrt{f_{d}^{-1}\left(\frac{\log(CT/\delta)}{cN}\right)}\right),$

where the last inequality holds with probability at least $1-\delta$ over the
sampling steps.

The convergence in probability follows from the corresponding statement of
Proposition D.1. ∎

##  Appendix E Additional Experiments and Details

###  E.1 Linear Gaussian model

We first consider the following 2-dimensional linear Gaussian SSM for which
exact inference can be carried out using Kalman techniques:

$X_{t}|\\{X_{t-1}=x\\}\sim\mathcal{N}\left(\text{diag}(\theta_{1}~{}\theta_{2})x,0.5\mathbf{I}_{2}\right),\quad
Y_{t}|\\{X_{t}=x\\}\sim\mathcal{N}(x,0.1\cdot\mathbf{I}_{2}).$ (43)

We simulate $T=150$ observations using
$\theta=(\theta_{1},\theta_{2})=(0.5,0.5)$. As a result, we expect in these
scenarios that the filtering distribution $p_{\theta}(x_{t}|y_{1:t})$ is not
too distinct from the smoothing distribution $p_{\theta}(x_{t}|y_{1:T})$ as
the latent process is mixing quickly. From Proposition 4.1, this is thus a
favourable scenario for methods ignoring resampling terms in the gradient as
the bias should not be very large. Figure 1, displayed earlier, shows
$\ell(\theta)$ obtained by Kalman and $\hat{\ell}(\theta;\bf{u})$ computed
regular PF and DPF for the same number $N=25$ of particles using
$q_{\phi}(x_{t}|x_{t-1},y_{t})=f_{\theta}(x_{t}|x_{t-1})$. The corresponding
gradient vector fields are given in Figure 1, where the gradient is computed
using the biased gradient from (Maddison et al., 2017; Naesseth et al., 2018;
Le et al., 2018) for regular PF.

We now compare the performance of the estimators $\hat{\theta}_{\text{SMLE}}$
(for DPF) and $\hat{\theta}_{\text{ELBO}}$ (for both regular PF and DPF)
learned using gradient with learning rate $10^{-4}$ on 100 steps, using $N=25$
for DPF and $N=500$ for regular PF, to $\hat{\theta}_{\text{MLE}}$ computed
using Kalman derivatives. We simulate $M=50$ realizations of $T=150$
observations using $\theta=(\theta_{1},\theta_{2})=(0.5,0.5)$. The ELBO
stochastic gradient estimates are computed using biased gradient estimates of
$\ell_{\text{ELBO}}(\theta)$ ignoring the contributions of resampling steps as
in (Maddison et al., 2017; Naesseth et al., 2018; Le et al., 2018) (we recall
that unbiased estimates suffer from very high variance) and unbiased gradients
of $\ell^{\text{ELBO}}(\theta)$ using DPF. We average $B$ parallel PFs to
reduce the variance of these gradients of the ELBO and also $B$ PFs (with
fixed random seeds) to compute the gradient of
$\hat{\ell}_{\text{SMLE}(\theta;\mathbf{u}_{1:B})}:=\tfrac{1}{B}\sum_{b=1}^{B}\hat{\ell}(\theta;\mathbf{u}_{b})$.
The results are given in Table 4. For this example,
$\hat{\theta}^{\text{DPF}}_{\text{ELBO}}$ maximizing
$\ell^{\text{ELBO}}_{\textrm{DPF}}(\theta)$ outperforms
$\hat{\theta}^{\textrm{PF}}_{\text{ELBO}}$ and $\hat{\theta}_{\text{SMLE}}$.
However, as $B$ increases, $\hat{\theta}_{\text{SMLE}}$ gets closer to
$\hat{\theta}^{\text{DPF}}_{\text{ELBO}}$ which is to be expected as
$\hat{\ell}_{\text{SMLE}}(\theta;\mathbf{u}_{1:B})\longrightarrow\ell^{\text{ELBO}}(\theta)$.
In Table 4, the Root Mean Square Error (RMSE) is defined as
$\sqrt{\tfrac{1}{M}\sum_{i=1}^{2}\sum_{k=1}^{M}(\hat{\theta}^{k}_{i}-\hat{\theta}^{k}_{\text{MLE},i})^{2}}$.

Table 4: $10^{3}\times$ RMSE444The Root Mean Square Error (RMSE) is defined as $\sqrt{\tfrac{1}{M}\sum_{i=1}^{2}\sum_{k=1}^{M}(\hat{\theta}^{k}_{i}-\hat{\theta}^{k}_{\text{MLE},i})^{2}}$. over 50 datasets - lower is better $B$ | $\hat{\theta}^{\text{PF}}_{\text{ELBO}}$ | $\hat{\theta}^{\text{DPF}}_{\text{ELBO}}$ | $\hat{\theta}_{\text{SMLE}}$
---|---|---|---
1 | 1.94 | 1.30 | 7.94
4 | 2.40 | 1.35 | 3.28
10 | 2.80 | 1.37 | 2.18

###  E.2 Variational Recurrent Neural Network

$N=32$ particles were used for training, with a regularization parameter of
$\epsilon=0.5$. The ELBO (scaled by sequence length) was used as the training
objective to maximise for each resampling/ DET procedure. The ELBO evaluated
on test data using $N=500$ particles and multinomial resampling. Resampling /
DET operations were carried out when effective sample (ESS) size fell below
$N/2$. Learning rate $0.001$ was used with the Adam optimizer.

Recall the state-space model is given by

$\displaystyle(R_{t},O_{t})$
$\displaystyle=\text{RNN}_{\theta}(R_{t-1},Y_{1:t-1},E_{\theta}(Z_{t-1})),$
$\displaystyle Z_{t}$
$\displaystyle\sim\mathcal{N}(\mu_{\theta}(O_{t}),\sigma_{\theta}(O_{t})),$
$\displaystyle\hat{p}_{t}$
$\displaystyle=h_{\theta}(E_{\theta}(Z_{t}),O_{t}),$ $\displaystyle
Y_{t}|X_{t}$ $\displaystyle\sim\text{Ber}(\hat{p}_{t}).$

Network architectures and data preprocessing steps were based loosely on
(Maddison et al., 2017). Given the low volume of data and sparsity of the
observations, relatively small neural networks were considered to prevent
overfitting, larger neural networks are considered in the more complex
robotics experiments. $R_{t}$ is of dimension $d_{r}=16$, $Z_{t}$ is of
dimension $d_{z}=8$. $E_{\theta}$ is a single layer fully connected network
with hidden layer of width $16$, output of dimension $16$ and RELU activation.

$\mu_{\theta}$ and $\sigma_{\theta}$ are both fully connected neural networks
with two hidden layers, each of $16$ units and RELU activation, the activation
function is not applied to the final output of $\mu_{\theta}$ but the softplus
is applied to the output of $\sigma_{\theta}$, which is the diagonal entries
of the covariance matrix of the normal distribution that is used to sample
$Z_{t}$.

$h_{\theta}$ is a single layer fully connected network with two hidden layers,
each of width $16$ and RELU activation. The final output is not put through
the RELU and is instead used as the logits for the Bernoulli distribution of
observations.

###  E.3 Robot Localization

Similar to the VRNN example, $N=32$ particles were used for training, with a
regularization parameter of $\epsilon=0.5$ and resampling / DET operations
were carried out when ESS size fell below $N/2$. Learning rate $0.001$ was
used with the Adam optimizer.

Network architectures and data preprocessing were based loosely on
(Jonschkowski et al., 2018). There are $3$ neural networks being considered:

  * •

Encoder $E_{\theta}$ maps RBG $24\times 24$ pixel images, hence dimension
$3\times 24\times 24$, to encoding of size $d_{E}=128$. This network consists
of a convolutional network (CNN) of kernel size $3$ and a single layer fully
connected network of hidden width $128$ and RELU activation.

  * •

Decoder $D_{\theta}$ maps encoding back to original image. This consists of a
fully connected neural network with three hidden layers of width $128$ and
RELU activation function. This is followed by a transposed convolution network
with matching specification to the CNN in the encoder, to return an output
with the same dimension as observation images, $3\times 24\times 24$.

  * •

Network $G_{\theta}$ maps the state
$S_{t}=(X^{(1)}_{t},X^{(2)}_{t},\gamma_{t})$ to encoding of dimension $128$.
First angle $\gamma_{t}$ was converted to $\sin(\gamma_{t}),\cos(\gamma_{t})$.
Then the augmented state
$(X^{(1)}_{t},X^{(2)}_{t},\sin(\gamma_{t}),\cos(\gamma_{t}))$ was passed to a
$3$ layer fully connected network with hidden layers of dimensions $16,32,64$
and RELU activation function, with final output of dimension $128$.