# Accelerated Inference for Partially Observed Markov Processes using
Automatic Differentiation

K. Tan1, G. Hooker1 and E. L. Ionides2
1Department of Statistics, University of Pennsylvania
2Department of Statistics, University of Michigan

(Draft compiled on )

###### Abstract

Automatic differentiation (AD) has driven recent advances in machine learning,
including deep neural networks and Hamiltonian Markov Chain Monte Carlo
methods. Partially observed nonlinear stochastic dynamical systems have proved
resistant to AD techniques because widely used particle filter algorithms
yield an estimated likelihood function that is discontinuous as a function of
the model parameters. We show how to embed two existing AD particle filter
methods in a theoretical framework that provides an extension to a new class
of algorithms. This new class permits a bias/variance tradeoff and hence a
mean squared error substantially lower than the existing algorithms. We
develop likelihood maximization algorithms suited to the Monte Carlo
properties of the AD gradient estimate. Our algorithms require only a
differentiable simulator for the latent dynamic system; by contrast, most
previous approaches to AD likelihood maximization for particle filters require
access to the system’s transition probabilities. Numerical results indicate
that a hybrid algorithm that uses AD to refine a coarse solution from an
iterated filtering algorithm show substantial improvement on current state-of-
the-art methods for a challenging scientific benchmark problem.

Many scientific models involve highly nonlinear stochastic dynamic systems
possessing significant random variation in both the process dynamics and the
measurements. Commonly, the latent system is modeled as a Markov process,
giving rise to a partially observed Markov process (POMP) model, also known as
a hidden Markov models or a state space model. POMP models arise in fields as
diverse as automated control [1], epidemiology [2, 3], ecology [4] and finance
[5, 6]. Despite their ubiquity, estimation and inference within this broad
class of models remains a challenging problem. This article concerns the use
automatic differentiation (AD) to construct improved algorithms for inference
on complex POMP models.

The particle filter, also known as sequential Monte Carlo, serves as the
foundation for various inference algorithms for POMP models. It provides an
unbiased estimate of the likelihood function [7], enabling Bayesian inference
[8, 9] and likelihood-based inference [10, 11]. Likelihood-based inference is
statistically efficient [12] and is robust to a moderate amount of Monte Carlo
error [13, 14]. An attractive feature of basic particle filter algorithms,
also known as boostrap filters, is that they do not require evaluation of the
transition density of the latent Markov process, enabling an arbitrary model
simulator to be plugged into the algorithm. This plug-and-play property is
useful for scientific applications [2]. Some plug-and-play methods depend on
the construction of low-dimensional summary statistics [15, 16], sacrificing
statistical efficiency for computational convenience. Plug-and-play methods
are often called likelihood-free [17], and it may be counter-intuitive that
likelihood-based inference is possible using likelihood-free methods. However,
iterated filtering algorithms have shown that this is practical in various
scientific investigations, such as [18, 19, 20, 21, 22, 23]. Nevertheless, the
Monte Carlo variability arising in iterated filtering applications becomes
increasingly problematic as the size of the data and the complexity of the
model increases. Algorithmic advances are needed to make plug-and-play
likelihood maximization for POMP models more numerically efficient. We develop
plug-and-play AD methods which provide a large computational advantage over
the methods used in these previous applications.

Recent advances in AD for particle filters [24, 25, 26, 27, 1] have drawn
attention to AD as a tool for inference in POMP models. However, existing
approaches are either asymptotically biased [24, 25], have high variance [28,
27], are computationally expensive [26, 29], or require access to transition
densities [28, 27, 1, 29].

Scibior and Wood [27] showed that the estimators derived by Poyiadjis et al.
[28] can be attained with standard AD software using an algorithmic procedure,
called a stop-gradient, which allows selected expressions to be evaluated but
not differentiated. Unfortunately, [27] used their estimator within an
algorithm which does not have the plug-and-play property and which has high
Monte Carlo variability, while the stop-gradient procedure is introduced
without mathematical motivation. We use [27] and [28] as our starting point,
while developing a new approach which remedies the weaknesses of these papers.

We start by presenting a new construction of the gradient estimator of [27]
and [28]. Specifically, we show how this estimator can be derived as the
direct derivative of a suitably weighted particle filter, which (for reasons
which will be discussed later) we call a Measurement Off-Parameter (MOP)
particle filter. In short, instead of directly differentiating through a basic
particle filter, MOP uses AD to differentiate through a smooth particle filter
constructed using measurement density ratios and a differentiable simulator.
Our MOP algorithm differs from previous smooth particle filters [30, 31] by
possessing the plug-and-play property. Critically, we also avoid the high
variance of [30] by generalizing the MOP particle filter to add a discount
factor, $\alpha\in[0,1]$. This discount factor enables MOP-$\alpha$ to
interpolate between the biased gradient estimator of [24], when $\alpha=0$,
and the high-variance gradient estimate from [28, 27], when $\alpha=1$
(Theorem 1). The bias-variance tradeoff induced by $\alpha$ suggests, both in
theory (Theorem 5) and in practice (Figure 1), the use of $\alpha$ values
strictly between 0 and 1.

MOP-$\alpha$ avoids the issue of asymptotic bias [26] arising from having to
drop terms arising from resampling when differentiating through the standard
particle filter as in [24] in a novel way. This particle filter is smooth by
construction, allowing us to bypass the seemingly incompatible paradox of
differentiating through discrete Monte Carlo sampling by instead using AD to
differentiate through a series of measurement density ratios and a
differentiable simulator. Yet, MOP-$\alpha$ can nevertheless provide strongly
consistent particle and log-likelihood estimates (Theorem 3), and when
$\alpha=1$, consistent score estimates (Theorem 4).

We derive a linear convergence rate for stochastic gradient descent (SGD)
using MOP-$\alpha$ in the presence of strong convexity (Theorem 2).
Critically, the estimator converges quickly within a neighborhood of the
maximum but struggles to reach this neighborhood in the presence of local
minima and saddle points. That behavior is complementary to that of iterated
filtering algorithms [10, 11] which provides a relatively fast and stable way
to identify this neighborhood. We therefore build a simple but effective
hybrid algorithm called Iterated Filtering with Automatic Differentiation
(IFAD) that warm-starts first-order or second-order gradient methods with a
preliminary solution obtained from a few rounds of iterated filtering.
Promising numerical results indicate that IFAD beats IF2 (and by the numerical
results of [11], also IF1 [10, 32] and the Liu-West filter [33]) on a
challenging problem in epidemiology, the Dhaka cholera model of [18].

These improvements also extend to Bayesian inference, as we show in Section 5
that we can use the MOP-$\alpha$ gradient estimates within a No-U-Turn Sampler
(NUTS) [34] in conjunction with a nonparametric empirical Bayes-style prior
estimated with IF2 to reduce the burn-in period of particle MCMC [8] on the
Dhaka cholera model of [18] from the $700,000$ iterations in [35] to just
$500$. To the best of our knowledge, attaining rapid mixing on this
challenging problem has not previously been achievable.

##  1 Problem Setup

Consider an unobserved Markov process $\\{X(t),t\geq t_{0}\\}$, with discrete-
time observations $Y_{1},...,Y_{N}$ realized at values
$y_{1}^{*},...,y_{N}^{*}$ at times $t_{1},...,t_{N}$. The process is
parameterized by an unknown parameter
$\theta\in\Theta\subseteq\mathbb{R}^{p}$, where the state $X(t)$ take values
in the state space ${\mathcal{X}}\subseteq\mathbb{R}^{d}$, the observations
$Y_{n}$ take values in ${\mathcal{Y}},$ and we write $X_{n}:=X(t_{n})$.

We suppose that the discete-time latent process model has a density
$f_{X_{n}|X_{n-1}}\left(x_{n}\mid x_{n-1};\theta\right)$. The existence of
this density is necessary to define the likelihood function, but plug-and-play
algorithms do not have access to evaluation of this density. Instead, they
have access to a corresponding simulator, which we write as
$\texttt{process}_{n}\left(x_{n}\mid x_{n-1};\theta\right)$. We set
$f_{Y_{n}|X_{n}}\left(y_{n}\mid x_{n};\theta\right)$ to be the measurement
density, which we suppose can be evaluated. We write $y_{n}^{*}$ for the
actual values of the observations that were observed. We call
$f_{X_{1:n}|Y_{1:n}}(x_{1:n}|y_{1:n}^{*};\theta)$ the posterior distribution
of states, and $f_{X_{n}|Y_{1:n}}(x_{n}|y_{1:n}^{*};\theta)$ the filtering
distribution at time $t_{n}$. Superscripts $x_{n,j}^{A}$ of particles denote
the ancestral trajectory of particle $j$ at time $n$, $a(\cdot)$ is the
ancestor function that maps a particle to its parent, $j\mapsto a(j)$, and
$k_{j}$ denotes the resample indices drawn for each particle $j$.

The above densities are defined on a probability space
$(\Omega,\Sigma,\mathbb{P})$ which is also assumed to enable construction of
independent replicates and all other random variables defined in our
algorithms. For the algorithmic interpretation of our theory, we identify the
random number seed with an element of the sample space $\omega\in\Omega$. The
random seed determines the sequence of pseudo-random numbers as generated by a
computer, just as the outcome $\omega\in\Omega$ generates the sequence of
random variables.

##  2 Off-Parameter Particle Filters

Our approach is as follows. Instead of directly differentiating through a
basic particle filter, we instead use AD to differentiate through a series of
measurement density ratios and a differentiable simulator. This is done
through the construction of a novel particle filtering algorithm called
Measurement Off-Parameter with discount factor $\alpha$ (MOP-$\alpha$),
defined by the pseudocode in Algorithm 1. Note that the plug-and-play property
forbids access to transition densities but permits access to the density of
the measurements conditional on the value of the latent process. The simulator
can be differentiated using AD as long as the underlying computer code for the
simulator is differentiable, and so MOP-$\alpha$ is applicable to a broad
class of POMP models.

The MOP-$\alpha$ algorithm resamples the particles according to an arbitrary
resampling rule that depends on both the target parameter value and a baseline
parameter value, and this explains the name measurement off-parameter. This is
intended as an analogy to off-policy learning in reinforcement learning.
Specifically, MOP-$\alpha$ evaluates the likelihood at some $\theta\in\Theta$,
but instead resamples the particles according to the indices generated by a
vanilla particle filter run at a baseline parameter, $\phi\in\Theta$. This
ensures that the resampling indices are invariant to $\theta$ when $\omega$
and $\phi$ are fixed, bypassing the issue of Monte Carlo resampling. We show
in Theorem 3 that MOP-$\alpha$ exactly targets the filtering distribution when
$\theta=\phi$ or $\alpha=1$.

We suppose that the simulator is a differentiable function of $\theta$ for
every fixed $\omega$, a condition that requires the latent process to be a
continuous random variable. This condition is also known as the
reparameterization trick [26]. Supposing also that the measurement model is
differentiable in $\theta$, direct differentiation of MOP-$\alpha$ is
available via AD. Theorem 1 shows that MOP-$\alpha$ is constructed so this
direct derivative obtains the score estimator of [28, 27] when $\alpha=1$ and
that of [24] when $\alpha=0$. Setting $\alpha<1$ adds bias but reduces
variance, raising opportunities for a favorable tradeoff as illustrated in
theory (Theorem 5) and in practice (Figure 1).

Algorithm 1 MOP-$\alpha$

Input: Number of particles $J$, timesteps $N$, measurement model
$f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n},\theta)$, simulator
$\texttt{process}_{n}(x_{n+1}|x_{n};\theta)$, evaluation parameter $\theta$,
behavior parameter $\phi$, seed $\omega$.

First pass: Set $\theta=\phi$ and fix $\omega$, yielding $X_{n,j}^{P,\phi}$,
$X_{n,j}^{F,\phi}$, $g^{\phi}_{n,j}$.

Second pass: Fix $\omega$, and filter at $\theta\neq\phi$:

Initialize  particles
${X}_{0,j}^{F,\theta}\sim{f}_{{X}_{0}}\left(\cdot\,;{\theta}\right)$, weights
$w^{F,\theta}_{0,j}=1$.
For $n=1,...,N$:
Accumulate discounted weights,
$w_{n,j}^{P,\theta}=\big{(}w_{n-1,j}^{F,\theta}\big{)}^{\alpha}$.
Simulate process model,
${X}_{n,j}^{P,\theta}\sim\texttt{process}_{n}\big{(}\cdot|{X}_{n-1,j}^{F,\theta};{\theta}\big{)}$.
Measurement density,
$g^{\theta}_{n,j}={f}_{{Y}_{n}|{X}_{n}}(y_{n}^{*}|{X}_{n,j}^{P,\theta}\,;{\theta})$.
Compute
$L_{n}^{B,\theta,\alpha}={\sum_{j=1}^{J}g^{\theta}_{n,j}\,w^{P,\theta}_{n,j}}\,\big{/}\,{\sum_{j=1}^{J}w^{P,\theta}_{n,j}}$.
Conditional likelihood under $\phi$,
$L_{n}^{\phi}=\frac{1}{J}\sum_{m=1}^{J}g^{\phi}_{n,m}$.
Select resampling indices $k_{1:J}$ with
$\mathbb{P}\big{(}k_{j}=m\big{)}\propto g^{\phi}_{n,m}$.
Obtain resampled particles ${X}_{n,j}^{F,\theta}={X}_{n,k_{j}}^{P,\theta}$.
Calculate resampled corrected weights
$w_{n,j}^{F,\theta}=w^{P,\theta}_{n,k_{j}}\,g^{\theta}_{n,k_{j}}\,\big{/}\,{g^{\phi}_{n,k_{j}}}$.
Compute
$L_{n}^{A,\theta,\alpha}=L_{n}^{\phi}\cdot{\sum_{j=1}^{J}w^{F,\theta}_{n,j}}\,\big{/}\,{\sum_{j=1}^{J}w^{P,\theta}_{n,j}}$.
Return: likelihood estimate
$\hat{\mathcal{L}}(\theta)=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}$ or
$\hat{\mathcal{L}}(\theta)=\prod_{n=1}^{N}L_{n}^{B,\theta,\alpha}$, filtering
distributions $\\{(X_{n,j}^{F,\theta},w^{F,\theta}_{n,j})\\}_{n,j=1}^{N,J}.$

###  2.1 Algorithm Outline

Algorithm 1 constructs two coupled sets of particles, one under
$\phi\in\Theta$, and another with the process model at $\theta\in\Theta$ but
with the resampling indices constructed from the first pass, with the baseline
parameter, $\phi$. The resampling indices, which we write as
$k_{j}\sim\text{Categorical}\big{(}g^{\phi}_{n,1},...,g^{\phi}_{n,J}\big{)}$,
are a function of $\phi$ for any value of $\theta$. For the categorical
distribution, we use systematic resampling [36, 37] which usually has superior
performance to multinomial resampling. If $\theta$ and $\phi$ coincide, one
only needs one particle filter run at $\theta=\phi$, otherwise one needs two
runs at the same seed $\omega\in\Omega$.

Algorithm 1 reweights the conditional likelihoods by a correction factor
accumulated over time to account for the resampling under $\phi$. Writing
$g^{\theta}_{n,j}={f}_{{Y}_{n}|{X}_{n}}(y_{n}^{*}|{X}_{n,j}^{P,\theta}\,;{\theta})$
for the measurement density and
$L_{n}^{\phi}=\frac{1}{J}\sum_{m=1}^{J}g^{\phi}_{n,m}$ for the conditional
likelihood estimate under $\phi$, we obtain two suitably weighted estimates of
the conditional likelihood under $\theta$,

$L_{n}^{B,\theta,\alpha}=\frac{\sum_{j=1}^{J}g^{\theta}_{n,j}\,w^{P,\theta}_{n,j}}{\sum_{j=1}^{J}w^{P,\theta}_{n,j}},\hskip
42.67912ptL_{n}^{A,\theta,\alpha}=L_{n}^{\phi}\cdot\frac{\sum_{j=1}^{J}w^{F,\theta}_{n,j}}{\sum_{j=1}^{J}w^{P,\theta}_{n,j}},$
(1)

where the weights for each particle are updated recursively to correct for the
cumulative error incurred by the off-parameter resampling:

$w_{n,j}^{P,\theta}=(w_{n-1,j}^{F,\theta})^{\alpha},\hskip
14.22636ptw^{F,\theta}_{n,j}=w^{P,\theta}_{n,k_{j}}\,g^{\theta}_{n,k_{j}}\,\big{/}g^{\phi}_{n,k_{j}},\hskip
14.22636ptw^{F,\theta}_{0,j}=1.$ (2)

The before-resampling conditional likelihood estimate
$L_{n}^{B,\theta,\alpha}$ is preferable in practice, as it has slightly lower
variance than the after-resampling estimate $L_{n}^{A,\theta,\alpha}$, but the
latter is useful in deriving properties of the MOP-$\alpha$ gradient estimate
such as that in Theorem 1.

If $\alpha=1$ in Algorithm 1, the weights, $w^{P,\theta}_{n,j}$ and
$w^{F,\theta}_{n,j}$, accumulate as $n$ increases. This leads to numerical
instability unless $N$ is small. For $\alpha<1$, the weights fom previous
timesteps are discounted, as in Equation (2). Heuristically, $\alpha$ controls
a rate of exponential decay of the memory that the filter at time $t_{n}$, and
its resulting gradient estimate, has over the ancestral trajectories prior to
$t_{n}$. With $\alpha<1$, Algorithm 1 continues to target the filtering
distribution and likelihood if $\theta=\phi$. However, we expect bias for
$\theta\neq\phi$ that shrinks as $\theta$ approaches $\phi$. This lets us
optimize a bias-variance tradeoff for the MOP-$\alpha$ score estimate.

When $\alpha=1$, MOP-$1$ maintains complete memory of each particle’s
ancestral trajectory, and Theorem 1 shows that it recovers the consistent but
high-variance gradient estimate from [28, 27]. At the other extreme, when
$\alpha=0$, MOP-$0$ considers only single-step transition dynamics, recovering
the low-variance but asymptotically biased gradient estimator of [24]. The
novel possibility of $0<\alpha<1$ is effective both in theory (Theorem 5) and
in practice (Figure 1).

Our approach has similarities with [30], who construct a deterministic local
approximation of the likelihood at the $k$-th optimization step via saving the
particles from a run at some $\theta_{k-1}$, perform a reweighting with
transition and measurement density ratios to evaluate the likelihood at any
sufficiently nearby $\theta$, and then use a $0$-th order optimizer such as
optim in R to choose a $\theta_{k}$ that maximizes the likelihood in that
neighborhood. However, due to their use of transition densities in the
reweighting, their approach is not plug-and-play. The variance of their log
likelihood estimate is also roughly $O(||\theta-\theta_{k-1}||_{2}^{2N})$,
where $N$ is the horizon. Our approach bypasses these two issues. The use of
the reparameterization trick to consider particle paths dependent on $\theta$
(and not $\theta_{k-1}$ or $\phi$) allows us to retain the plug-and-play
property. On the other hand, our use of AD allows us to only evaluate the
likelihood and score estimates when $\theta=\phi$. The use of AD bypasses the
issue of exponentially increasing variance, but nevertheless can result in
problematic polynomial variance. Our introduction of $\alpha$ avoids that
difficulty.

###  2.2 MOP-$\alpha$ Encompasses the Estimators of [28, 27, 24]

[27] show that the estimate of [24] is the gradient of a vanilla particle
filter when resampling terms are dropped, and also recover the estimate of
[28] through the use of a stop-gradient operation in the AD procedure. It
turns out that both of these, when applied on the bootstrap filter, correspond
to special cases of MOP-$\alpha$.

######  Theorem 1 (MOP-$0$ and MOP-$1$ Functional Forms).

Writing $\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)$ for the gradient estimate
yielded by MOP-$\alpha$ when $\theta=\phi$ and using the after-resampling
conditional likelihood estimate so that
$\hat{\mathcal{L}}(\theta)=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}$, when
$\alpha=0$,

$\nabla_{\theta}\hat{\ell}^{0}(\theta)=\frac{1}{J}\sum_{n=1}^{N}\sum_{j=1}^{J}\nabla_{\theta}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right),\vspace*{-2.5mm}$

yielding the estimate of [24] for the bootstrap filter. When $\alpha=1$,

$\nabla_{\theta}\hat{\ell}^{1}(\theta)=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,j}^{A,F,\theta}\right),\vspace*{-2.5mm}$

yielding the estimator of [28, 27] for the bootstrap filter.

We defer the proof to Appendix A. The argument relies on a useful
decomposition of the after-resampling conditional likelihood estimate
$L_{n}^{A,\theta,\alpha}$ that yields a telescoping product in the MOP-$1$
case. Repeated applications of the log-derivative identity that
$\nabla_{x}\log(f(x))=(\nabla_{x}f(x))/f(x)$, and noting $\theta=\phi$ implies
that $w_{n,j}^{P,\theta}$ evaluates to $1$, yield the result.

This further illustrates how $\alpha$ dictates the memory of the gradient
estimate. As [27] remarks, the MOP-$0$ estimator of [24] only considers
single-step quantities, and is “memoryless” beyond a single step. This is in
contrast to the case when $\alpha=1,$ as that estimate, studied by [28], only
considers the surviving particles at time $N$ and so fully tracks dependencies
over time.

Figure 1: Illustration of the bias-variance tradeoff induced by the
discounting hyperparameter $\alpha$, on the Dhaka cholera model of [18]. We
display the MSE of score estimates for the trend in transmission, evaluated at
the MLE.

###  2.3 Summary of Theoretical Guarantees

The construction of MOP-$\alpha$ bypasses the issue of differentiating through
a Monte Carlo algorithm with discontinuous resampling by turning it into a
problem of differentiating through a simulator and a series of measurement
density ratios. We defer the theoretical analysis of MOP-$\alpha$ to Section
6, but first we highlight a few key results. MOP-$\alpha$ estimates the
likelihood and conditional distributions of latent variables, as one expects
of a particle filter (Theorem 3). When differentiated, we obtain the
estimators of [28, 27, 24] as special cases (Theorem 1). In particular,
MOP-$1$ is consistent for the score (Theorem 4), has rates for its bias and
variance under different choices of $\alpha$ (Theorem 5) that illustrate the
desirable bias-variance tradeoff observed empirically in Figure 1, and enjoys
a linear rate of convergence for gradient descent with the resulting gradient
estimate (Theorem 2).

##  3 Practical Maximum-Likelihood Estimation

If $\theta$ is evaluated at $\phi$, the particles at $\theta$ and $\phi$
coincide. One then only needs to run one particle filter at $\theta=\phi$,
setting the particles at $\phi$ to be copies of the particles at $\theta$
where gradients don’t propagate. This is done algorithmically through the
stop_gradient() function in JAX, providing a mathematical justification for
the use of the stop-gradient operation by [27].

###  3.1 Optimization

This still leaves us with the question of designing an effective procedure for
likelihood maximization with MOP-$\alpha$. We propose a simple algorithm we
call Iterated Filtering with Automatic Differentiation (IFAD) in Algorithm 2
that runs a few iterations of IF2 to warm-start an iterative (first or second-
order) method that uses the MOP-$\alpha$ gradient estimate. This leverages
IF2’s quick empirical convergence to a neighborhood of the MLE, overcoming the
tendency of gradient methods to get stuck in saddle points and local minima.
Conversely, switching to gradient ascent with MOP-$\alpha$ score estimates
lets one bypass the difficulty that IF2 has with optimizing the last few units
of log-likelihood. Combining these two methods in this way lets us enjoy the
best of both worlds.

The convergence of IF2 (and so the first stage of IFAD) to a neighborhood of
the MLE happens rapidly in practice. In the case of the Dhaka cholera model of
[18], when an aggressive geometric cooling multiplier of 0.95 and initial
random walk standard deviation of 0.02 is used, initial convergence happens
within 40 iterations. Finding the MLE itself with IF2, however, takes much
longer, as one has to use a less aggressive cooling rate to do so. For
example, [11] use 100 iterations with the Dhaka cholera model of [18], while
[38] use 200 with Model 1 in [39]. By simply requiring that the first stage of
IFAD get to a neighborhood of the MLE and not the MLE itself, we are able to
substantially reduce the number of IF2 iterations required.

###  3.2 Linear Convergence Rates

On the other hand, the second stage of IFAD enjoys a linear convergence rate
under the usual smoothness and strong convexity assumptions in the theory of
convex optimization, as we show below in Theorem 2.

######  Theorem 2 (Linear Convergence of IFAD).

Consider the second stage of IFAD (Algorithm 2) where one stops if
$\|\nabla_{\theta}\hat{\ell}^{\alpha}(\theta_{m})\|\leq(1+\sigma)\epsilon$,
where $\sigma\geq\frac{4\Gamma}{(1-\beta)}$, for some $\beta\in(0,1)$. Assume
$-\ell$ is strongly convex and smooth, $\gamma
I\preceq\nabla_{\theta}^{2}(-\ell)\preceq\Gamma I$. Choose the learning rate
$\eta$ such that $\eta\leq\frac{c(1-\beta)}{2\Gamma}$. Then, for sufficiently
large $\alpha$ and $J$ to ensure the score estimate is an
$\epsilon$-approximation and the minimum eigenvalue of $H$ is greater than
some $c>0$ for all $m$ with probability at least $1-\delta$, the second stage
of IFAD converges linearly to the MLE:

$\ell(\theta^{*})-\ell(\theta_{m+1})\leq\Big{(}1-\eta\beta\,\frac{8\gamma}{9c}\Big{)}\big{(}\ell(\theta^{*})-\ell(\theta_{m})\big{)}.$

We borrow from tools in the field of randomized numerical linear algebra [40]
to solve this problem. The proof, which is similar to that of Theorem 6 in
[40], is deferred to Appendix B. For simplicity, we only prove this for
MOP-$1$ and note that the general case follows as the required
$\epsilon$-approximation can still be obtained for sufficiently large
$\alpha$. When $\alpha=1$, the required $J$ is given in the supplementary
information in Lemmas 3 and 4. This result shows that particle estimates for
the gradient can enjoy the same linear convergence rate on well-conditioned
problems as gradient descent with access to the score. As , we note that it is
possible.

We therefore see that the second stage of IFAD converges linearly to the MLE
if (1) the log-likelihood surface is $\gamma$-strongly convex in a
neighborhood of the MLE and (2) the first stage of IFAD successfully reaches a
(high-probability) basin of attraction of the MLE. This happens fairly often
in practice, for example, when sufficient regularity conditions for local
asymptotic normality of the MLE hold. We conjecture that this applies to the
entirety of IFAD, as IF2 converges very quickly to a neighborhood of the MLE,
and will explore this in future work.

Algorithm 2 IFAD

Input: Number of particles $J$, timesteps $N$, IF2 cooling schedule
$\eta_{m}$, MOP-$\alpha$ discounting parameter $\alpha$, $\theta_{0}$, $m=0.$
Run IF2 until initial ”convergence” under cooling schedule $\eta_{m}$, or for
a fixed number or iterations, to obtain $\\{\Theta_{j},j=1,...,J\\}$, set
$\theta_{m}:=\frac{1}{J}\sum_{j=1}^{J}\Theta_{j}.$
While procedure not converged:
Run Algorithm 1 to obtain $\hat{\ell}(\theta_{m}).$
Obtain
$g(\theta_{m})=\nabla_{\theta_{m}}\big{(}-\hat{\ell}(\theta_{m})\big{)}$,
$H(\theta_{m})$ s.t. $\lambda_{\min}\big{(}H(\theta_{m})\big{)}\geq c$.
Update $\theta_{m+1}:=\theta_{m}-\eta(H(\theta_{m}))^{-1}g(\theta_{m})$,
$m:=m+1.$
Return $\hat{\theta}:=\theta_{m}.$

##  4 Application to a Cholera Transmission Model

The cholera transmission model that [18] developed for Dhaka, Bangladesh, has
been used to benchmark the performance of various POMP inference methods [11,
35, 41], and we employ it here for the same purpose. This model categorizes
individuals in a population as susceptible, $S(t)$, infected, $I(t)$, and
recovered, $R(t)$ and so is called an SIR compartmental model. In this case,
the compartment $R(t)$ is further subdivided into three recovered compartments
$R^{1}(t)$, $R^{2}(t)$, $R^{3}(t)$ denoting varying degrees of cholera
immunity. We write $P(t)$ for the total population, and $M_{n}$ for the
cholera deaths in each month. As in [18, 11], the transition dynamics follow a
series of stochastic differential equations:

$\displaystyle dS$ $\displaystyle=\big{(}k\epsilon
R^{k}+\delta(S-P)-\lambda(t)S\big{)}\,dt+dP-({\sigma SI}/{H})\,dB,$
$\displaystyle dI$
$\displaystyle=\big{(}\lambda(t)S-(m+\delta+\gamma)I\big{)}\,dt+({\sigma
SI}/{H})\,dB,$ $\displaystyle dR^{1}$ $\displaystyle=\big{(}\gamma
I-(k\epsilon+\delta)R^{1}\big{)}\,dt,\hskip 5.69054pt\dots$ $\displaystyle
dR^{k}$ $\displaystyle=\big{(}k\epsilon
R^{k-1}-(k\epsilon+\delta)R^{k}\big{)}\,dt,\vspace*{-2mm}$

with Brownian motion $B(t)$, cholera death rate $m$, recovery rate $\gamma$,
mean immunity duration $1/\epsilon$, standard deviation of the force of
infection $\sigma$, and population death rate $\delta=0.02$. The force of
infection, $\lambda_{t}$, is modeled by splines $(s_{j})_{j=1}^{6}$

$\lambda_{t}=\exp\left(\beta_{\text{trend}}(t-t_{0})+\sum_{j=1}^{6}\beta_{j}s_{j}(t)\right)\frac{I}{P}+\exp\left(\sum_{j=1}^{6}\omega_{j}s_{j}(t)\right),\vspace*{-2mm}$

where the coefficients $(\beta_{j})_{j=1}^{6}$ model seasonality in the force
of infection, $\beta_{\text{trend}}$ models the trend in the force of
infection, and the $\omega_{j}$ represent seasonality of a non-human
environmental reservoir of disease. The measurement model for observed monthly
cholera deaths is given by $Y_{n}\sim{\mathcal{N}}(M_{n},\tau^{2}M_{n}^{2})$,
where $M_{n}=\gamma\int_{t_{n-1}}^{t_{n}}I(s)\,ds$ is the modeled number of
cholera deaths in that month.

Figure 2: A compartment flow diagram for the Dhaka cholera model from [18].

###  4.1 Results

We tested IFAD against IF2 on a global search problem for the Dhaka cholera
model. We re-implemented IF2 to do so, but we compare our results with the
results of [11] (labeled ”IF2 2015”). Our re-implementation ourperforms that
of [11], likely due to a better choice of algorithmic parameters. For each
method, we performed 100 searches, initialized with 100 initial starting
parameter vectors drawn uniformly from the same wide bounding box used in
[11]. We summarize our findings below.

| Best Log-Likelihood | Rank
---|---|---
IFAD-0.97 | -3750.2 | 1
IFAD-0 | -3752.2 | 2
IFAD-1 | -3754.6 | 3
IF2 Warm Start | -3757.3 | 4
IF2 | -3758.2 | 5
IF2 2015 | -3768.6 | 6
Table 1: Maximum log-likelihood found by IF2, IFAD, and MOP alone. IFAD
performs the best among all methods. Our implementation of IF2 outperforms
that of [11], but still ultimately underperforms IFAD. IFAD manages to find
the MLE, matching the highest log-likelihood previously found in the Dhaka
cholera model implemented within the pomp package of [37].

#### IFAD Successfully Finds the MLE:

Previously, an MLE at a log-likelihood of $-3748.6$ was reported by [37]. This
MLE was obtained with much computational effort, using many global and local
IF2 searches, and with the assistance of likelihood profiling. Meanwhile, [11]
only achieve a maximum log-likelihood of $-3768.6$, while the best log-
likelihood found by [18] was only $-3793.4$. Despite being initialized for a
global search, IFAD manages to get much closer to the MLE over the $100$
searches than [11] and finds it up to $2.5$ standard deviations of Monte Carlo
error, as seen in Table 1. On this problem, the sequence of local searches,
refinement, and likelihood profiling that was previously required for finding
the MLE is not necessary with the IFAD algorithm. In other words, the
improvement in numerical efficiency of IFAD over IF2 is so large that routine
application of IFAD (consisting of a collection of Monte Carlo replications
from random starting points) generates results outside the reach of a routine
application of IF2. When interpreting Table 1, bear in mind that by Wilks’
Theorem, any difference of over $1.92$ log units has statistical relevance
when testing one parameter at the $0.05$ significance level, and therefore
potentially has scientific value.

#### IFAD Outperforms Both IF2 and MOP Alone:

While IF2 quickly approaches a neighborhood of the MLE within only 40
iterations, performing IF2 alone ultimately fails to achieve the last few log-
likelihood units, as no IF2 search comes within 7 log-likelihood units of the
MLE (as seen in Figure 3). Conversely, when we tried gradient descent with MOP
alone, we encountered many failed searches. This is a difficult, nonconvex,
and noisy problem with $18$ parameters, and the search gets stuck in local
minima and saddle points, failing to approach the MLE.

IFAD, in comparison, approaches the MLE quickly due to the IF2 warm-start (as
seen in Figure 5) and also succeeds at refining the coarse solution found by
the warm-start with MOP gradient steps to find the MLE (as seen in Figures 3
and 4). IFAD therefore successfully combines the best qualities of IFAD and
MOP, outperforming either of them alone. Monte Carlo replication is
appropriate for IFAD, or any Monte Carlo algorithm used to solve a challenging
numerical problem, but Figure 4 shows that IFAD can find higher likelihood
values, more quickly and more reliably than the previous state-of-the-art. In
particular, Figure 4 shows that IFAD with $\alpha=0.97$ can maximize the
challenging likelihood of [18] using a modest number of iterations and Monte
Carlo replications, whereas previously it was necessary to carry out an
extensive customized search or to live with an incompletely maximized
likelihood.

For the results presented here, we did not include the many common heuristics
used in the machine learning and optimization literature in the gradient
descent stage. We used constant learning rates of $0.01,0.05$, and $0.2$ for
IFAD-$1,0,$ and $0.97$ respectively, and a constant cooling rate of $0.95$ for
our IF2 implementation. While techniques such as annealing learning rates,
gradient normalization, and momentum could further improve the performance of
IFAD in other simulations we performed, we chose to report the results of the
simplest implementation to serve as a baseline for the method’s performance.
The results for this basic implementation are already sufficient to show the
high potential of the approach.

Figure 3: Scatterplots depicting the performance of IFAD against that of IF2.
Left: Paired searches from the same starting point. Controlling for initial
starting point, tuning $\alpha$ allows IFAD-$0.97$ to strictly improve on IF2,
[28], and [24], on almost every iteration. Right: Q-Q plot of ranked IFAD
searches against ranked IF2 searches. It is clear that on average, IFAD
performs best, and manages to find the MLE while no IF2 search successfully
gets within 7 log-likelihood units of it. The dotted red line shows the true
maximized log-likelihood.

Figure 4: Left: Raincloud plot depicting the performance of IFAD and IF2 where
we plot the results of the best run out of every ten runs, representing the
common procedure of running a few searches and choosing the best one. Right:
Raincloud plot of all searches. IFAD outperforms all other methods, and the
gradient steps improve on the warm-start given by running 40 IF2 iterations.
The dotted red line shows the true maximized log-likelihood. Figure 5:
Optimization progress of IFAD and IF2. The dashed orange line depicts the
median warm-start given by running 40 IF2 iterations. While running 60 more
iterations of IF2 improves upon the median warm-start, doing so ultimately
underperforms IFAD. We see that IFAD has better tail control and successfully
reaches the MLE. We use a dotted red line to display the MLE.

##  5 Application to Bayesian Inference

Particle MCMC, as introduced by [8], is arguably the most popular method for
full-information plug-and-play Bayesian inference for POMP models. The plug-
and-play property enables its use in simulation-based Bayesian inference for
complex scientific models where computing transition densities is not
feasible, such as in disease modeling. However, the particle Metropolis-
Hastings algorithm often experiences slow mixing and requires long burn-in
periods. For instance, [35] reported needing $700,000$ burn-in iterations for
effective sampling in a cholera model. This computational burden necessitated
simplifying the scientific model to reduce computational costs, by adjusting
the length of each Euler timestep to be a month instead of a day.

We employ a NUTS sampler powered by a MOP-$\alpha$, with a nonparametric
empirical prior initialized by the IF2 warm start from the previous section.
To construct the prior, we perform a kernel density estimate (KDE) on the
parameter swarm from the warm start, and use the output of that as an
empirical prior. The KDE is performed to ensure some degree of mutual
contiguity between the densities. This yields a prior that is roughly
equivalent to a parameter cloud centered at some point in the neighborhood of
the MLE, with a standard deviation of $\sqrt{(0.02\cdot 0.95^{40})^{2}\cdot
600}\approx 6.3\%$ of the MLE. Using a NUTS sampler powered by MOP-$\alpha$
with this prior significantly lowers the mixing time of the Markov chain to as
little as $500$ iterations, as seen in Figure 6.

Figure 6: Convergence diagnostics for NUTS with the informative nonparametric
empirical Bayes prior from applying KDE to preliminary searches using IF2,
with $4$ chains. We display the results for the trend parameter in the Dhaka
cholera model of [18]. The NUTS sampler mixes quickly, and the posterior
estimates from each chain share roughly the same posterior mode.

For completeness, we note that without the empirical Bayes-style prior from
IF2, NUTS with the MOP-$\alpha$ estimate alone tended to diverge. In contrast,
particle Metropolis-Hastings (with or without the empirical Bayes prior from
IF2) fails to effectively explore the parameter space. We display the results
of the above in Figures 7 and 8 in Appendix G. The result obtained in this
section, of a speedup in the mixing time of particle MCMC by over three orders
of magnitude, provides hope that particle Bayesian inference might finally be
practical for the complex scientific models commonly encountered in areas like
epidemiology.

##  6 Theoretical Analysis of MOP-$\alpha$

Here, we will show that MOP-$\alpha$ targets the filtering distribution and
likelihood under $\theta$, is consistent when $\alpha=1$, and characterize
rates for its bias and variance under different choices of $\alpha$. To do so,
we require the following assumptions:

  1. (A1)

Continuity of the Likelihood. $\ell(\theta)$ has more than two continuous
derivatives in a neighborhood
$\left\\{\theta:\ell(\theta)>\lambda_{1}\right\\}$ for some
$\lambda_{1}<\sup_{\varphi}\ell(\varphi)$.

  2. (A2)

Bounded Process Model. There exist $\underline{M},\bar{M}$ such that
$0<\underline{M}\leq
f_{X_{n}|X_{n-1}}(x_{n}|x_{n-1};\theta)\leq\bar{M}<\infty$.

  3. (A3)

Bounded Measurement Model. There exist $\underline{G},\bar{G}$ such that
$0<\underline{G}\leq f_{Y_{n}\mid X_{n}}\left(y_{n}^{*}\mid
x_{n};\theta\right)\leq\bar{G}<\infty$ and there exists $G^{\prime}(\theta)$
with $\|\nabla_{\theta}\log f_{Y_{n}\mid X_{n}}\left(y_{n}^{*}\mid
x_{n};\theta\right)\|_{\infty}\leq G^{\prime}(\theta)<\infty$.

  4. (A4)

Bounded Gradient Estimates. There are functions
$G(\theta),H(\theta):\Theta\to[0,\infty)$ uniformly bounded by
$G^{*},H^{*}<\infty$, so the MOP-$\alpha$ gradient and Hessian estimates at
$\theta=\phi$ are almost surely bounded by $G(\theta)$ and $H(\theta)$ for all
$\alpha$.

  5. (A5)

Differentiability of Density Ratios and Simulator. The measurement density,
$f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n};\theta)$, and simulator have more than two
continuous derivatives in $\theta$.

The fact that the likelihood estimate yielded by MOP-$\alpha$ has more than
two continuous derivatives in $\theta$ follows from the construction of the
likelihood estimate in equations 1 and 2, as well as Assumptions (A3) and
(A5).

###  6.1 MOP-$\alpha$ Targets the Filtering Distribution

We show here that MOP-$\alpha$ targets the filtering distribution under
$\theta$ and is strongly consistent for the likelihood under $\theta$ when
$\alpha=1$ or $\theta=\phi$. Employing $\alpha<1$ lead to inconsistency when
$\theta$ deviates from $\phi$, but this bias disappears as $\theta$ approaches
$\phi$.

While the result is presented here as specific to MOP-$\alpha$, we actually
prove a more general result in Appendix D. That is, we show a strong law of
large numbers for triangular arrays of particles with off-parameter
resampling, where we resample the particles according to an arbitrary
resampling rule not necessarily in proportion to the targeted distribution of
interest and employ weights that encode the cumulative discrepancy between the
resampling and the target distribution instead of equal weights.

######  Theorem 3 (MOP-$\alpha$ Targets the Filtering Distribution and
Likelihood).

When $\alpha=1$ or $\theta=\phi$, MOP-$\alpha$ targets the filtering
distribution under $\theta$ and is strongly consistent for the likelihood
under $\theta$. That is, for
$\pi_{n}(\theta)=f_{X_{n}|Y_{1:n}}(x_{n}|y_{1:n}^{*};\theta)$ and any
measurable and bounded functional $h$ and for
$\hat{\mathcal{L}}(\theta)=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}$ or
$\prod_{n=1}^{N}L_{n}^{B,\theta,\alpha}$, it holds that

$\frac{\sum_{j=1}^{J}h(x_{n,j}^{F,\theta})\,w_{n,j}^{F,\theta}}{\sum_{j=1}^{J}w_{n,j}^{F,\theta}}\stackrel{{\scriptstyle
a.s.}}{{\to}}E_{\pi_{n}(\theta)}\big{[}h(X_{n})\big{]},\hskip
14.22636pt\hat{\mathcal{L}}(\theta)\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathcal{L}(\theta).$

###### Proof.

We provide a proof sketch here, deferring most of the details, and discussion
of the after-resampling conditional likelihood estimate
$L_{n}^{A,\theta,\alpha}$, to Appendix D. When $\theta=\phi$, regardless of
the value of $\alpha$, the ratio ${g_{n,j}^{\theta}}/{g_{n,j}^{\phi}}=1,$ and
this reduces to the vanilla particle filter. When $\alpha=1$ and
$\theta\neq\phi,$ suppose inductively that
$\\{(X^{F,\theta}_{n-1,j},w^{F,\theta}_{n-1,j})\\}_{j=1}^{J}$ targets
$f_{X_{n-1}|Y_{1:n-1}}(x_{n-1}|y^{*}_{1:n-1};\theta)$. It can then be shown
that $\\{(X^{P,\theta}_{n,j},w^{P,\theta}_{n,j})\\}_{j=1}^{J}$ targets
$f_{X_{n}|Y_{1:n-1}}(x_{n}|y^{*}_{1:n-1};\theta)$, that
$\\{(X^{P,\theta}_{n,j},w^{P,\theta}_{n,j}\,g^{\theta}_{n,j})\\}_{j=1}^{J}$
targets $f_{X_{n}|Y_{1:n}}(x_{n}|y^{*}_{1:n};\theta)$, and that weighting the
particles by
$(X^{F,\theta}_{n,j},w^{F,\theta}_{n,j})=(X^{P,\theta}_{n,k_{j}},w^{P,\theta}_{n,k_{j}}\,g^{\theta}_{n,k_{j}}\big{/}g^{\phi}_{n,k_{j}})$,
resampling the $k_{j}$ with probabilities proportional to $g^{\phi}_{n,j}$,
also targets $f_{X_{n}|Y_{1:n}}(x_{n}|y^{*}_{1:n};\theta)$. If the likelihood
is estimated with the before-resampling conditional likelihoods
$\hat{\mathcal{L}}(\theta)=\prod_{n=1}^{N}L_{n}^{B,\theta,\alpha}$, the strong
consistency is a direct consequence of our earlier result that
$\\{\big{(}X^{P,\theta}_{n,j},w^{P,\theta}_{n,j}\big{)}\\}$ targets
$f_{X_{n}|Y_{1:n-1}}(x_{n}|y^{*}_{1:n-1};\theta)$. ∎

###  6.2 MOP-$1$ Is Consistent for the Score

Despite showing that the MOP-$1$ gradient estimate yields the estimate of [28,
27] when applied to the bootstrap filter, [28, 27] estimate the Fisher score
by

$\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log
f_{X_{0:N},Y_{1:N}}\left(x_{0:n,j}^{A,F,\theta},y_{1:N}^{*};\theta\right)$

and not

$\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,j}^{A,F,\theta};\theta\right).$

It is therefore not immediately apparent that these two converge to the same
thing. As such we directly show the consistency of the MOP-$1$ gradient
estimate below. We present an abbreviated proof here, postponing the full
argument to Appendix E.

######  Theorem 4 (Consistency of MOP-$1$ Gradient Estimate).

The gradient estimate of MOP-$\alpha$ when $\alpha=1$, $\theta=\phi$ is
strongly consistent for the score:
$\nabla_{\theta}\hat{\ell}_{J}^{1}(\theta)\stackrel{{\scriptstyle
a.s.}}{{\to}}\nabla_{\theta}\ell(\theta)$ as $J\to\infty$.

###### Proof.

Fix $\omega\in\Omega$, and set $\phi=\theta$, where $\theta$ is the point at
which we wish to evaluate the gradient. The sequence
$(\nabla_{\theta}\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega))_{J\in\mathbb{N}}$
is uniformly bounded over all $J$ by Assumption (A4). Again by Assumption
(A4), the second derivative of
$\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega)|_{\theta=\theta^{\prime}}$ is also
uniformly bounded over all $J$ by $H^{*}$ for almost every $\omega\in\Omega$
and every $\theta^{\prime}\in\Theta$. So
$(\nabla_{\theta}\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega))_{J\in\mathbb{N}}$
is uniformly Lipschitz, and therefore uniformly equicontinuous for almost
every $\omega\in\Omega$.

By Arzela-Ascoli, there is a uniformly convergent subsequence. But there is
only one subsequential limit, as we can treat the gradient estimate at
$\theta$ as a bounded functional of the particles by Assumption (A4), allowing
us to apply Theorem 3 to see that the sequence
$(\nabla_{\theta}\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega))_{J\in\mathbb{N}}$
converges pointwise for $\theta=\phi$ and almost every $\omega\in\Omega$. So
the whole sequence must converge uniformly to
$\lim_{J\to\infty}\nabla_{\theta}\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega).$

With uniform convergence for the derivatives established, we can swap the
limit and derivative and obtain, in conjunction with the strong consistency
$\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega)\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathcal{L}(\theta)$ in Theorem 3, that for almost every
$\omega\in\Omega$,
$\lim_{J\to\infty}\nabla_{\theta}\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega)=\nabla_{\theta}\lim_{J\to\infty}\hat{\mathcal{L}}_{J}^{1}(\theta)(\omega)=\nabla_{\theta}\mathcal{L}(\theta).$
The result then follows by the continuous mapping theorem. ∎

###  6.3 MOP-$\alpha$ Error, Bias and Variance

We now provide a result showing that MOP-$\alpha$ has a desirable bias-
variance tradeoff when $0<\alpha<1$. This is achieved because it combines
favorable properties from the low-variance but asymptotically biased estimate
of [24] and the high-variance but consistent estimate of [28]. As the bias
itself is difficult to analyze, we instead analyze the MSE and variance. The
below result applies for any $\alpha\in[0,1)$, but not $\alpha=1$, as the
gradient estimate when $\alpha=1$ has no forgetting properties.

######  Theorem 5.

When $\alpha\in(0,1)$ and $\theta=\phi$, define
$\psi(\alpha)=(\alpha^{k}+\alpha^{k+1}-\alpha)/(1-\alpha)$. There exists an
$\epsilon>0$ depending on $\bar{M},\underline{M},\bar{G},\underline{G}$ as in
[42] such that the MSE and variance of MOP-$\alpha$ are:

$\displaystyle\mathbb{E}\big{\|}\nabla_{\theta}\ell(\theta)-\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{\|}_{2}^{2}$
$\displaystyle\lesssim$ $\displaystyle\min_{k\leq
N}Np\,G^{\prime}(\theta)^{2}\left(k^{2}J^{-1}+(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor}+k+\psi(\alpha)\right),$ (3)
$\displaystyle\mathrm{Var}\big{(}\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{)}$
$\displaystyle\lesssim$ $\displaystyle\min_{k\leq
N}Np\,G^{\prime}(\theta)^{2}\left(\frac{k^{2}}{(1-\alpha)^{2}J}+\frac{\alpha^{k}}{1-\alpha}N\right).$
(4)

We defer the proof to Appendix F, but provide a brief outline here. The
variance bound can be reduced to the approximation error between MOP-$\alpha$
and a variant called MOP-$(\alpha,k)$ where the discounted weights are
truncated at lag $k$. The discount factor, $\alpha$, ensures strong mixing.
The covariance of the derivative of MOP-$(\alpha,k)$ can be controlled with
Davydov’s inequality and a standard $L^{p}$ error bound for the particle
filter. The MSE bound considers the error between the score and MOP-$(1,k)$,
and the error between MOP-$(1,k)$ and MOP-$\alpha$. The latter can be shown to
be $\tilde{O}\big{(}Np\,G^{\prime}(\theta)^{2}(k+\psi(\alpha))\big{)}$, and
the former can be controlled with a result on the the forgetting of the
particle filter from [42] and the very same $L^{p}$ error bound mentioned
above.

#### Interpretation:

Theorem 5 provides theoretical understanding of the results observed in Figure
1 that show the empirically favorable bias-variance tradeoff enjoyed by
MOP-$\alpha$. The first term in the MSE bound in Equation 3 corresponds to the
error of the particle approximation, the second mixing error, and the third
and fourth the error between the MOP-$(1,k)$ estimate and the MOP-$\alpha$
estimate. As $\alpha$ goes to $1$, choosing $k$ appropriately, the first,
third and fourth term (corresponding to the variance) increases, while the
second term (corresponding to the bias) decreases. Likewise, the first term in
the variance bound in Equation 4 corresponds to the error of the particle
approximation, while the second corresponds to mixing error. As $\alpha$ goes
to $1$, the variance increases.

#### Bias-Variance Tradeoff:

We show in Appendix F that the variance of MOP-$0$ is
$\tilde{O}\big{(}Np\,G^{\prime}(\theta)^{2}\big{/}J\big{)}$. Previously, [28]
established that the variance of MOP-$1$ is $\tilde{O}(N^{4}/J)$, ignoring
factors of $p$ and $G^{\prime}$. Combining these results with Theorem 5 shows
that MOP-$\alpha$ interpolates between MOP-$0$ and MOP-$1$, as $N\leq
Nk^{2}\leq N^{4}$, with a phase transition as soon as $\alpha<1$. The phase
transition arises as even though the particle filter has forgetting properties
[42], the resulting derivative estimate of [28] does not, and we require both
for the variance reduction.

In contrast, as $Nk^{2}\leq N^{4}$ and $\alpha,k$ can be as large as desired
to balance the impact of the last three terms, MOP-$\alpha$ achieves a lower
MSE than MOP-$1$ does. We also show that MOP-$0$ achieves a MSE of
$\tilde{O}\big{(}NpG^{\prime}(\theta)^{2}(J^{-1}+(1-\epsilon)^{\lfloor
1/c\log(J)\rfloor})\big{)}$, where the second term corresponds to uncontrolled
mixing error. Unlike MOP-$\alpha$, MOP-$0$ has no opportunity to tune $\alpha$
and $k$ to reduce said mixing error, leading to uncontrolled asymptotic bias.

##  7 Computational Efficiency

MOP-$\alpha$ and IFAD are fast algorithms, both in theory and practice. In
line with the cheap gradient principle of [43], getting a gradient estimate
from MOP-$\alpha$ takes no more than 6 times that of the runtime of the
particle filter. In our simulations, MOP-$\alpha$ took 3.75 times the runtime
of the particle filter. MOP-$\alpha$ and IFAD therefore share the same $O(NJ)$
time complexity as the particle filter, unlike the $O(NJ^{2})$ complexity of
[26] and Algorithm 2 in [28, 27].

Our implementation of the particle filter, simulator, and MOP-$\alpha$ in JAX
[44] enabled us to take advantage of just-in-time compilation and GPU
acceleration, even with a simulator written in Python. This led to a 16x
speedup (379ms vs 6.29s on a Intel i9-13900K CPU and NVIDIA RTX3090 GPU) over
the CPU-only implementation of the particle filter (with a simulator written
in C++) in the pomp package of [37].

##  8 Discussion

If the simulator is discontinuous in $\theta$, MOP-$\alpha$ no longer applies.
We are working on a variation of MOP-$\alpha$ applicable to this case.
Specifically, differentiable transition density ratios can be used instead of
a differentiable simulator. The off-parameter treatment of the measurement
model is extended to an off-parameter simulation for the dynamic process. In
that case, one requires access to the transition densities, or at least their
ratios. The discounting parameter, $\alpha$, can be incorporated as for
MOP-$\alpha$. This algorithm may be better suited to large discrete state
spaces than either MOP-$\alpha$ or the Baum-Welch algorithm. Additionally, our
gradient estimate can be used for variational inference to approximate the
posterior distribution over latent states [24] and parameters.

Discounting the weights by some $\alpha\in[0,1]$ is not the only way to
interpolate between the estimators of [24] and [28]. As the proof of Theorem 5
implies, we can also truncate the weights at a fixed lag, corresponding to the
MOP-$(1,k)$ estimate mentioned. The analysis is similar, with comparable rates
but with a slightly less convenient implementation.

Practical likelihood-based data analysis involves many likelihood
optimizations [18, 19, 20, 21, 22, 23]. In particular, profile likelihood
confidence intervals require a sequence of optimizations. Additionally, a
careful scientist should consider many model variations, to see whether the
conclusions of the study are sensitive to alternative model specifications.
Coding model variations is relatively simple when using plug-and-play
inference methodology. However, assessing the scientific potential of these
variations requires likelihood optimization. If many rounds of laborious
searching are required to attain the maximum (as arises in previous methods)
that has practical consequences for the rate at which scientists can evaluate
hypotheses. We have demonstrated, for the first time, a plug-and-play maximum
likelihood approach which can directly maximize the likelihood for the complex
benchmark model of [18], requiring some Monte Carlo replication but not
prolonged sequences of explorations. We anticipate that this will promote
scientific advances in all domains where complex POMP models arise.

## Acknowledgments

This research was supported by National Science Foundation grant DMS-1761603.
We thank Nicolas Chopin for helpful communications regarding the strong law of
large numbers for off-parameter resampled particle filters. We also thank
Arnaud Doucet for helpful discussion regarding the manuscript.

## References

  * [1] A. Singh, O. Makhlouf, M. Igl, J. Messias, A. Doucet, and S. Whiteson, “Particle-based score estimation for state space model learning in autonomous driving,” arXiv, vol. 2212.06968, 2022.
  * [2] D. He, E. L. Ionides, and A. A. King, “Plug-and-play inference for disease dynamics: Measles in large and small towns as a case study,” Journal of the Royal Society Interface, vol. 7, pp. 271–283, 2010.
  * [3] T. Stocks, T. Britton, and M. Höhle, “Model selection and parameter estimation for dynamic epidemic models via iterated filtering: Application to rotavirus in Germany,” Biostatistics, vol. 21, pp. 400–416, 09 2018.
  * [4] J. Knape and P. de Valpine, “Fitting complex population models by combining particle filters with markov chain monte carlo,” Ecology, vol. 93, no. 2, pp. 256–263, 2012.
  * [5] J. Kim and D. S. Stoffer, “Fitting stochastic volatility models in the presence of irregular sampling via particle methods and the EM algorithm,” Journal of Time Series Analysis, vol. 29, no. 5, pp. 811–833, 2008.
  * [6] C. Bretó, “On idiosyncratic stochasticity of financial leverage effects,” Statistics and Probability Letters, vol. 91, pp. 20–26, 2014.
  * [7] P. Del Moral, Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications.  New York: Springer, 2004.
  * [8] C. Andrieu, A. Doucet, and R. Holenstein, “Particle Markov chain Monte Carlo methods,” Journal of the Royal Statistical Society, Series B, vol. 72, pp. 269–342, 2010.
  * [9] N. Chopin, P. E. Jacob, and O. Papaspiliopoulos, “SMC2: An efficient algorithm for sequential analysis of state space models,” Journal of the Royal Statistical Society, Series B, vol. 75, no. 3, pp. 397–426, 2013.
  * [10] E. L. Ionides, C. Bretó, and A. A. King, “Inference for nonlinear dynamical systems,” Proceedings of the National Academy of Sciences of the USA, vol. 103, pp. 18438–18443, 2006.
  * [11] E. L. Ionides, D. Nguyen, Y. Atchadé, S. Stoev, and A. A. King, “Inference for dynamic and latent variable models via iterated, perturbed Bayes maps,” Proceedings of the National Academy of Sciences of the USA, vol. 112, no. 3, pp. 719––724, 2015.
  * [12] Y. Pawitan, In all likelihood: statistical modelling and inference using likelihood.  Oxford University Press, 2001.
  * [13] E. L. Ionides, C. Breto, J. Park, R. A. Smith, and A. A. King, “Monte Carlo profile confidence intervals for dynamic systems,” Journal of the Royal Society Interface, vol. 14, pp. 1–10, 2017.
  * [14] N. Ning, E. L. Ionides, and Y. Ritov, “Scalable Monte Carlo inference and rescaled local asymptotic normality,” Bernoulli, vol. 27, pp. 2532–2555, 2021.
  * [15] S. N. Wood, “Statistical inference for noisy nonlinear ecological dynamic systems,” Nature, vol. 466, pp. 1102–1104, 2010.
  * [16] T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. P. Stumpf, “Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems,” Journal of the Royal Society Interface, vol. 6, pp. 187–202, 2009.
  * [17] J. Owen, D. J. Wilkinson, and C. S. Gillespie, “Likelihood free inference for Markov processes: A comparison,” Statistical Applications in Genetics and Molecular Biology, vol. 14, no. 2, pp. 189–209, 2015.
  * [18] A. A. King, E. L. Ionides, M. Pascual, and M. J. Bouma, “Inapparent infections and cholera dynamics,” Nature, vol. 454, pp. 877–880, 2008.
  * [19] I. M. Blake, R. Martin, A. Goel, N. Khetsuriani, J. Everts, C. Wolff, S. Wassilak, R. B. Aylward, and N. C. Grassly, “The role of older children and adults in wild poliovirus transmission,” Proceedings of the National Academy of Sciences of the USA, vol. 111, no. 29, pp. 10604–10609, 2014.
  * [20] M. Pons-Salort and N. C. Grassly, “Serotype-specific immunity explains the incidence of diseases caused by human enteroviruses,” Science, vol. 361, no. 6404, pp. 800–803, 2018.
  * [21] R. Subramanian, Q. He, and M. Pascual, “Quantifying asymptomatic infection and transmission of COVID-19 in New York City using observed cases, serology, and testing capacity,” Proceedings of the National Academy of Sciences of the USA, vol. 118, no. 9, 2021.
  * [22] S. J. Fox, M. Lachmann, M. Tec, R. Pasco, S. Woody, Z. Du, X. Wang, T. A. Ingle, E. Javan, M. Dahan, K. Gaither, M. E. Escott, S. I. Adler, S. C. Johnston, J. G. Scott, and L. A. Meyers, “Real-time pandemic surveillance using hospital admissions and mobility data,” Proceedings of the National Academy of Sciences of the USA, vol. 119, no. 7, p. e2111870119, 2022.
  * [23] J. M. Drake, A. Handel, É. Marty, E. B. O’Dea, T. O’Sullivan, G. Righi, and A. T. Tredennick, “A data-driven semi-parametric model of SARS-CoV-2 transmission in the United States,” PLOS Computational Biology, vol. 19, no. 11, p. e1011610, 2023.
  * [24] C. Naesseth, S. Linderman, R. Ranganath, and D. Blei, “Variational sequential Monte Carlo,” in Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics (A. Storkey and F. Perez-Cruz, eds.), vol. 84 of Proceedings of Machine Learning Research, pp. 968–977, PMLR, 09–11 Apr 2018.
  * [25] R. Jonschkowski, D. Rastogi, and O. Brock, “Differentiable particle filters: End-to-end learning with algorithmic priors,” arXiv, vol. 1805.11122, 2018.
  * [26] A. Corenflos, J. Thornton, G. Deligiannidis, and A. Doucet, “Differentiable particle filtering via entropy-regularized optimal transport,” in Proceedings of the 38th International Conference on Machine Learning (M. Meila and T. Zhang, eds.), vol. 139 of Proceedings of Machine Learning Research, pp. 2100–2111, PMLR, 18–24 Jul 2021.
  * [27] A. Ścibior and F. Wood, “Differentiable particle filtering without modifying the forward pass,” arXiv, vol. 2106.10314, 2021.
  * [28] G. Poyiadjis, A. Doucet, and S. S. Singh, “Particle approximations of the score and observed information matrix in state space models with application to parameter estimation,” Biometrika, vol. 98, no. 1, pp. 65–80, 2011.
  * [29] X. Chen and Y. Li, “Normalising flow-based differentiable particle filters,” arXiv, vol. 2403.01499, 2024.
  * [30] A. Svensson, F. Lindsten, and T. B. Schön, “Learning nonlinear state-space models using smooth particle-filter-based likelihood approximations,” IFAC–PapersOnLine, vol. 51, no. 15, pp. 652–657, 2018.
  * [31] S. Malik and M. K. Pitt, “Particle filters for continuous likelihood evaluation and maximisation,” Journal of Econometrics, vol. 165, no. 2, pp. 190–209, 2011.
  * [32] E. L. Ionides, A. Bhadra, Y. Atchadé, and A. A. King, “Iterated filtering,” Annals of Statistics, vol. 39, pp. 1776–1802, 2011.
  * [33] J. Liu and M. West, “Combining parameter and state estimation in simulation-based filtering,” in Sequential Monte Carlo Methods in Practice (A. Doucet, N. de Freitas, and N. Gordon, eds.), pp. 197–224, New York: Springer, 2001.
  * [34] M. D. Homan and A. Gelman, “The No-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo,” J. Mach. Learn. Res., vol. 15, p. 1593–1623, jan 2014.
  * [35] M. Fasiolo, N. Pya, and S. N. Wood, “A comparison of inferential methods for highly nonlinear state space models in ecology and epidemiology,” Statistical Science, vol. 31, pp. 96–118, 02 2016.
  * [36] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear, non-Gaussian Bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, pp. 174 – 188, 2002.
  * [37] A. A. King, D. Nguyen, and E. L. Ionides, “Statistical inference for partially observed Markov processes via the R package pomp,” Journal of Statistical Software, vol. 69, pp. 1–43, 2016.
  * [38] J. Wheeler, A. Rosengart, J. Zhuxun, K. H. E. Tan, N. Treutle, and E. L. Ionides, “Informing policy via dynamic models: Cholera in Haiti,” PLOS Computational Biology, vol. 20, p. e1012032, 2024.
  * [39] E. C. Lee, D. L. Chao, J. C. Lemaitre, L. Matrajt, D. Pasetto, J. Perez-Saez, F. Finger, A. Rinaldo, J. D. Sugimoto, M. E. Halloran, I. M. Longini, R. Ternier, K. Vissieres, A. S. Azman, J. Lessler, and L. C. Ivers, “Achieving coordinated national immunity and cholera elimination in Haiti through vaccination: A modelling study,” The Lancet Global Health, vol. 8, no. 8, pp. e1081–e1089, 2020.
  * [40] F. Roosta-Khorasani and M. W. Mahoney, “Sub-sampled Newton methods I: Globally convergent algorithms,” arXiv, vol. 1601.04737, 2016.
  * [41] N. Wycoff, J. W. Smith, A. S. Booth, and R. B. Gramacy, “Voronoi candidates for Bayesian optimization,” 2024.
  * [42] J. Karjalainen, A. Lee, S. S. Singh, and M. Vihola, “On the forgetting of particle filters,” arXiv, vol. 2309.08517, 2023.
  * [43] S. Kakade and J. D. Lee, “Provably correct automatic subdifferentiation for qualified programs,” arXiv, vol. 1809.08530, 2019.
  * [44] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, “JAX: composable transformations of Python+NumPy programs,” 2018.
  * [45] P. Del Moral and A. Doucet, “On a class of genealogical and interacting Metropolis models,” in Séminaire de Probabilités XXXVII (J. Azéma, M. Émery, M. Ledoux, and M. Yor, eds.), pp. 415–446, Berlin: Springer, 2003.
  * [46] P. Del Moral and E. Rio, “Concentration inequalities for mean field particle models,” The Annals of Applied Probability, vol. 21, jun 2011.
  * [47] N. Chopin and O. Papaspiliopoulos, An Introduction to Sequential Monte Carlo.  Springer, 2020.
  * [48] N. Chopin, “Central limit theorem for sequential Monte Carlo methods and its application to Bayesian inference,” Annals of Statistics, vol. 32, pp. 2385–2411, 2004.

##  Appendix A MOP-$\alpha$ Functional Forms

Theorem 1 follows immediately as a consequence of the following results,
Lemmas 1 and 2.

######  Lemma 1.

Write $\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)$ for the gradient estimate
yielded by MOP-$\alpha$ when $\theta=\phi$. Consider the case where we use the
after-resampling conditional likelihood estimate so that
$\hat{\mathcal{L}}(\theta)=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}$. When
$\alpha=1$,

$\nabla_{\theta}\hat{\ell}^{1}(\theta)=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,j}^{A,F,\theta}\right),$ (5)

yielding the estimator of [28, 27] with the bootstrap filter.

###### Proof.

Consider the case of MOP-$\alpha$ when $\alpha=1$ and $\theta=\phi$. We then
have a nice telescoping product property for the after-resampling likelihood
estimate:

$\hat{\mathcal{L}}^{1}(\theta):=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}=\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{\sum_{j=1}^{J}w_{n,j}^{F,\theta}}{\sum_{j=1}^{J}w_{n,j}^{P,\theta}}=\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{\sum_{j=1}^{J}w_{n,j}^{F,\theta}}{\sum_{j=1}^{J}w_{n-1,j}^{F,\theta}}=\left(\frac{1}{J}\sum_{j=1}^{J}w_{N,j}^{F,\theta}\right)\prod_{n=1}^{N}L_{n}^{\phi},$
(6)

where the third equality follows from the choice of $\alpha=1$, and the fourth
equality is the resulting telescoping property. The log-derivative identity
lets us decompose the score estimate as

$\nabla_{\theta}\hat{\ell}^{1}(\theta)=\frac{\nabla_{\theta}\hat{\mathcal{L}}^{1}(\theta)}{\hat{\mathcal{L}}^{1}(\theta)}=\frac{\nabla_{\theta}\left(\frac{1}{J}\sum_{j=1}^{J}w_{N,j}^{F,\theta}\right)\prod_{n=1}^{N}L_{n}^{\phi}}{\prod_{n=1}^{N}L_{n}^{\phi}}=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{N,j}^{F,\theta}.$
(7)

From (7), we see that the derivative of the log-likelihood estimate is

$\nabla_{\theta}\hat{\ell}^{1}(\theta):=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{N,j}^{F,\theta}.$
(8)

We proceed to decompose (8). First, observe that as $\alpha=1$,

$w_{n,j}^{P,\theta}=w_{n-1,j}^{F,\theta}\frac{g_{n,j}^{\theta}}{g_{n,j}^{\phi}}=\prod_{i=1}^{n}\frac{g_{i,j}^{A,P,\theta}}{g_{i,j}^{A,P,\phi}},$
(9)

where we use the $(\cdot)^{A}$ superscript to denote the ancestral trajectory
of the $j$-th prediction or filtering particle at timestep $n$. Note that this
quantity is the cumulative product of measurement density ratios over the
ancestral trajectory of the $j$-th prediction particle at timestep $n$. We
then use the log-derivative identity again, yielding the following expression
for the gradient of the log-weights as the sum of the log measurement
densities over the ancestral trajectory:

$\displaystyle\frac{\nabla_{\theta}w_{n,j}^{P,\theta}}{w_{n,j}^{P,\theta}}=\nabla_{\theta}\log
w_{n,j}^{P,\theta}$ $\displaystyle=$
$\displaystyle\nabla_{\theta}\log\left(\prod_{i=1}^{n}\frac{g_{i,j}^{A,P,\theta}}{g_{i,j}^{A,P,\phi}}\right)$
(10) $\displaystyle=$ $\displaystyle\nabla_{\theta}\sum_{i=1}^{n}\left(\log
g_{i,j}^{A,P,\theta}-\log g_{i,j}^{A,P,\phi}\right)$ (11) $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\nabla_{\theta}\log g_{i,j}^{A,P,\theta}.$ (12)

This is equal to the gradient of the logarithm of the conditional density of
the observed measurements given the ancestral trajectory of the $j$-th
prediction particle up to timestep $n$:

$\displaystyle\nabla_{\theta}\sum_{n=1}^{N}\log g_{n,j}^{A,\theta}$
$\displaystyle=$
$\displaystyle\nabla_{\theta}\log\left(\prod_{n=1}^{N}g_{n,j}^{A,P,\theta}\right)$
(13) $\displaystyle=$
$\displaystyle\nabla_{\theta}\log\left(\prod_{n=1}^{N}f_{Y_{n}|X_{n}}\left(y_{n}^{*}|x_{n,j}^{A,P,\theta}\right)\right)$
(14) $\displaystyle=$ $\displaystyle\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,j}^{A,P,\theta}\right).$ (15)

Multiplying both sides of the expression by $w_{N,j}^{P,\theta}$ yields an
expression for the gradient of the weights at timestep $N$:

$\nabla_{\theta}w_{N,j}^{P,\theta}=w_{N,j}^{P,\theta}\sum_{n=1}^{N}\nabla_{\theta}\log
g_{n,j}^{A,P,\theta}=w_{N,j}^{P,\theta}\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,j}^{A,P,\theta}\right).$ (16)

Substituting the above identity into the log-likelihood decomposition obtained
earlier in Equation 7 yields

$\nabla_{\theta}\hat{\ell}^{1}(\theta):=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{N,j}^{F,\theta}=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{N,k_{j}}^{P,\theta}=\frac{1}{J}\sum_{j=1}^{J}w_{N,k_{j}}^{P,\theta}\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,k_{j}}^{A,P,\theta}\right).$ (17)

Finally, observing that $\theta=\phi$ implies $w_{N,j}^{F,\theta}=1$, we
obtain

$\nabla_{\theta}\hat{\ell}^{1}(\theta):=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,j}^{A,F,\theta}\right).$ (18)

This yields the gradient estimators of [28, 27] when applied to the bootstrap
filter. ∎

Note that the variance of the MOP-$\alpha$ log-likelihood estimate scales
poorly with $N$ the moment $\theta\neq\phi$. This can be seen by observing
that

$\nabla_{\theta}\hat{\ell}^{1}(\theta):=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{N,j}^{F,\theta}=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{N,k_{j}}^{P,\theta}=\frac{1}{J}\sum_{j=1}^{J}w_{N,k_{j}}^{P,\theta}\nabla_{\theta}\log
f_{Y_{1:N}|X_{1:N}}\left(y_{1:N}^{*}|x_{1:n,k_{j}}^{A,P,\theta}\right).$ (19)

When $\theta\neq\phi$, we see that $w_{N,k_{j}}^{P,\theta}=O(c^{N})$. When
$\theta=\phi$, this is a special case of the [28] estimator, which has
$O(N^{4})$ variance by a property of functionals of the particle filter [45].

######  Lemma 2.

Write $\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)$ for the gradient estimate
yielded by MOP-$\alpha$ when $\theta=\phi$. Consider the case where we use the
after-resampling conditional likelihood estimate so that
$\hat{\mathcal{L}}(\theta)=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}$. When
$\alpha=0$,

$\nabla_{\theta}\hat{\ell}^{0}(\theta)=\frac{1}{J}\sum_{n=1}^{N}\sum_{j=1}^{J}\nabla_{\theta}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right),$
(20)

yielding the estimate of [24] when applied to the bootstrap filter.

###### Proof.

First, write

$s_{n,j}=\frac{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\theta})}{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\phi})}$

as shorthand for the measurement density ratios. Observe that, when
$\alpha=0,$, the likelihood estimate becomes

$\displaystyle\hat{\mathcal{L}}^{0}(\theta):=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}$
$\displaystyle=$
$\displaystyle\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{\sum_{j=1}^{J}w_{n,j}^{F,\theta}}{\sum_{j=1}^{J}w_{n,j}^{P,\theta}}$
(21) $\displaystyle=$
$\displaystyle\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{1}{J}\sum_{j=1}^{J}s_{n,j}$
(22) $\displaystyle=$
$\displaystyle\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{1}{J}\sum_{j=1}^{J}\frac{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\theta})}{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\phi})}.$
(23)

We lose the nice telescoping property observed in the MOP-$1$ case, but this
expression still yields something useful. This is because its gradient when
$\theta=\phi$ is therefore

$\displaystyle\nabla_{\theta}\hat{\ell}^{0}(\theta)$
$\displaystyle:=\sum_{n=1}^{N}\nabla_{\theta}\log\left(L_{n}^{\phi}\frac{1}{J}\sum_{j=1}^{J}s_{n,j}\right)$
(24)
$\displaystyle=\sum_{n=1}^{N}\frac{\nabla_{\theta}\left(L_{n}^{\phi}\frac{1}{J}\sum_{j=1}^{J}s_{n,j}\right)}{\left(L_{n}^{\phi}\frac{1}{J}\sum_{j=1}^{J}s_{n,j}\right)}$
(25)
$\displaystyle=\sum_{n=1}^{N}\frac{\sum_{j=1}^{J}\nabla_{\theta}s_{n,j}}{\sum_{j=1}^{J}s_{n,j}}$
(26)
$\displaystyle=\sum_{n=1}^{N}\frac{1}{J}\sum_{j=1}^{J}\frac{\nabla_{\theta}f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)}{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\phi};\phi)}$
(27)
$\displaystyle=\frac{1}{J}\sum_{n=1}^{N}\sum_{j=1}^{J}\nabla_{\theta}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right),$
(28)

where we use the log-derivative trick in the second equality, observe that
$\sum_{j=1}^{J}s_{n,j}=J$ when $\theta=\phi$ in the fourth equality, and use
the log-derivative trick again while noting that $\theta=\phi$ in the fifth
equality. This yields the desired result. ∎

##  Appendix B Optimization Convergence Analysis

The analysis in this section roughly follows the analysis in [40], except with
the caveat that none of the matrix concentration bounds they use apply here as
the particles are dependent. We instead use the concentration inequality from
[46] to bound the gradient and Hessian estimates. In this section, we fix
$\omega\in\Omega$ only within each filtering iteration, evaluate Algorithm 1
at $\theta=\phi$, and analyze Algorithm 2 post-iterated filtering.

The convergence analysis in Theorem 2 is limited to the case where $-\ell$ is
$\gamma$-strongly convex. Though it is true that in a neighborhood of the
optimum local asymptotic normality holds and the log-likelihood is strongly
convex in this neighborhood, in practice likelihood surfaces for POMPs are
often highly nonconvex globally. The convergence to an optimum, local or
global, must therefore be sensitive to initialization.

###  B.1 Bounding the Gradient

######  Lemma 3 (Concentration of Measure for Gradient Estimate).

Consider the gradient estimate obtained by MOP-1, which we know by Theorem 4
is consistent for the score, where $\theta=\phi$. For
$\|\nabla_{\theta}\hat{\ell}^{1}(\theta)-\nabla_{\theta}\ell(\theta)\|_{2}$ to
be bounded by $\epsilon$ with probability $1-\delta$, we require

$\displaystyle
J>\max\left\\{2G(\theta)\frac{r_{N}\sqrt{p}}{\epsilon}\left(1+h^{-1}\left(\log\left(\frac{2p}{\delta}\right)\right)\right),8G(\theta)^{2}\beta_{N}^{2}\frac{p\log(2p/\delta)}{\epsilon^{2}}\right\\},$
(29)

where $NG^{\prime}(\theta)\leq G(\theta)$ are defined in Assumptions (A3) and
(A4), $h(t)=\frac{1}{2}(t-\log(1+t))$, and $\beta_{N}$ and $r_{N}$ are two
additional finite model-specific constants that do not depend on $J$, but do
depend on $N$ and $p$, as defined in [46]. Equivalently, with probability at
least $1-\delta$, it holds that

$\displaystyle\|\nabla_{\theta}\hat{\ell}^{1}(\theta)-\nabla_{\theta}\ell(\theta)\|_{2}\leq
G(\theta)\left(\frac{r_{N}\sqrt{p}}{J}(1+h^{-1}(\log(2p/\delta)))+\sqrt{\frac{2p\log(2p/\delta)}{J}}\beta_{N}\right).$
(30)

#### Remark:

According to [46], under some regularity conditions, $r_{N}$ and $\beta_{N}$
are linear in the trajectory length $N$. This corresponds to the finding by
[28] that the variance of the estimate is at least quadratic in the trajectory
length, and their remark that the result of [45] establishes that the $L_{p}$
error is bounded by $O(N^{2}J^{-1/2})$ (equivalently, the variance is bounded
by $O(N^{4}J^{-1})$) after accounting for the sum over timesteps. The MOP-$1$
variance upper bound is therefore in fact $O(N^{4})$, in contrast to the
MOP-$\alpha$, where $\alpha<1$, upper bound of $O(N)$.

###### Proof.

We will seek to use the concentration inequality of [46] to bound the
deviation of the gradient estimate from the gradient of the negative log-
likelihood in the sup norm with a union bound. Fix $\theta=\phi$. From the
decomposition in the proof of Lemma 1, as $w_{N,j}^{F,\theta}=1$ when
$\theta=\phi$, we have that

$\nabla_{\theta}\hat{\ell}(\theta):=\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}w_{N,j}^{F,\theta}=\frac{1}{J}\sum_{j=1}^{J}\sum_{n=1}^{N}w_{N,k_{j}}^{P,\theta}\nabla_{\theta}\log
g_{n,k_{j}}^{A,\theta}=\frac{1}{J}\sum_{j=1}^{J}\sum_{n=1}^{N}\nabla_{\theta}\log
g_{n,k_{j}}^{A,\theta}.$ (31)

Define
$\varphi_{n}^{i}(x_{n,j}^{F,\theta}):=\frac{\partial}{\partial\theta_{i}}\log
g_{n,k_{j}}^{A,\theta}$, which is a functional of the filtering particles
$x_{n,j}^{F,\theta}=x_{n,k_{j}}^{P,\theta}$. These are bounded measurable
functionals bounded by $G^{\prime}(\theta)$ by Assumption (A3). Therefore,
these have bounded oscillation, satisfying the requirement that
$\text{osc}\left(\frac{\partial}{\partial\theta_{i}}\varphi_{i}(x_{n,j}^{P,\theta})\right)\leq
G^{\prime}(\theta)$. Note that [46] in fact assume $\text{osc}(f)\leq 1$, so
we simply scale their bound accordingly.

Now we apply the Hoeffding-type concentration inequality from Del Moral and
Rio [46] and a union bound over each $\varphi_{n}^{i}(x_{n,j}^{F,\theta})$,
totaling $N$ timesteps and $p$ parameters, to find that

$\displaystyle\max_{n=1,...,N}\left\lVert\frac{1}{J}\sum_{j=1}^{J}\nabla_{\theta}\log
g_{n,k_{j}}^{A,\theta}-\nabla_{\theta}\ell_{n}(\theta)\right\rVert_{\infty}\leq
G^{\prime}(\theta)\left(\frac{r_{N}}{J}(1+h^{-1}(t))+\sqrt{\frac{2t}{J}}\beta_{N}\right)$
(32)

with probability at least $1-2Np\exp(-t)$. Although the above concentration
inequality only considers the error from the expectation under the filtering
distribution, we invoke the consistency of MOP-$1$ shown in Theorem 4 to
establish that the expectation under the filtering distribution is in fact the
score. It therefore holds that with the same probability, that when summing
over $N$, as $NG^{\prime}(\theta)\leq G(\theta)$,

$\displaystyle\left\lVert\frac{1}{J}\sum_{j=1}^{J}\sum_{n=1}^{N}\nabla_{\theta}\log
g_{n,k_{j}}^{A,\theta}-\nabla_{\theta}\ell(\theta)\right\rVert_{\infty}$
$\displaystyle\leq
G^{\prime}(\theta)N\left(\frac{r_{N}}{J}(1+h^{-1}(t))+\sqrt{\frac{2t}{J}}\beta_{N}\right)$
(33) $\displaystyle\leq
G(\theta)\left(\frac{r_{N}}{J}(1+h^{-1}(t))+\sqrt{\frac{2t}{J}}\beta_{N}\right).$
(34)

We split the $\delta$ failure probability among these $2Np$ terms, to find
$\delta\leq 2Np\exp(-t)$, and therefore, $t\leq\log(2Np/\delta)$, where
$h(t)=\frac{1}{2}(t-\log(1+t))$. The two additional model-specific parameters
are $\beta_{t}$ and $r_{t}$, which do not depend on $J$. The analogous bound
for the 2-norm follows from scaling the right-hand side by $\sqrt{p}$, to
require

$\displaystyle\|\nabla_{\theta}\hat{\ell}(\theta)-\nabla_{\theta}\ell(\theta)\|_{2}\leq
G(\theta)\left(\frac{r_{N}\sqrt{p}}{J}(1+h^{-1}(\log(2p/\delta)))+\sqrt{\frac{2p\log(2p/\delta)}{J}}\beta_{N}\right).$
(35)

We therefore need

$\displaystyle
J>\max\left\\{2G(\theta)\frac{r_{N}\sqrt{p}}{\epsilon}\left(1+h^{-1}\left(\log\left(\frac{2p}{\delta}\right)\right)\right),8G(\theta)^{2}\beta_{N}^{2}\frac{p\log(2p/\delta)}{\epsilon^{2}}\right\\}.$
(36)

∎

###  B.2 Bounding Hessian Estimates

Should one choose to use a second-order method involving a particle Hessian
estimate, we provide a guarantee for its positive-definiteness below.

######  Lemma 4 (Minimum Eigenvalue Bound for Hessian Estimate).

Assume that the Hessian of the negative log-likelihood
$H=\sum_{j=1}^{J}\mathbb{E}H_{j}$ has a minimum eigenvalue $0<\gamma<1$, and
that $\mathbb{E}\lambda_{\min}(H_{j})=\gamma^{\prime}>0$. If

$J>\max\left\\{\frac{2r_{t}(1+h^{-1}(t))+2c}{\gamma^{\prime}},\frac{2(2t\beta_{t}^{2}+c)^{2}}{\gamma^{\prime
2}}\right\\}\geq\frac{r_{t}(1+h^{-1}(t))}{\gamma^{\prime}}+\sqrt{2tJ}\beta_{t}/\gamma^{\prime}+c/\gamma^{\prime}$
(37)

then $\hat{H}(\theta)$ is invertible and positive definite with minimum
eigenvalue greater than or equal to
$c\in(0,\sum_{j=1}^{J}\mathbb{E}\lambda_{\min}(H_{j}))$, with probability at
least $1-\exp(-t)$.

###### Proof.

Write $\hat{H}(\theta)=\hat{H}=\sum_{j=1}^{J}H_{j}$ for the estimate of the
negative of the Hessian, where $H_{j}$ is an element of the outer sum over the
$J$ particles.

As the negative log-likelihood is convex, we want to bound the minimum
eigenvalue of $\hat{H}(\theta)$ from below with high probability, so that all
the eigenvalues of $\hat{H}(\theta)$ are positive with high probability. This
ensures that the estimated Hessian is invertible and positive-definite.

It is known that the minimum eigenvalue of a symmetric matrix is concave.
Therefore, it suffices to show that the first inequality in the below
expression

$0<\sum_{j=1}^{J}\lambda_{\min}(H_{j})\leq\lambda_{\min}\left(\sum_{j=1}^{J}H_{j}\right)=\lambda_{\min}(\hat{H})$
(38)

holds with high probability. We apply the particle Hoeffding concentration
inequality from [46] to find that

$\displaystyle\frac{1}{J}\sum_{j=1}^{J}\lambda_{\min}(H_{j})-\mathbb{E}_{\tilde{\pi}_{t}}\lambda_{\min}(H_{j})$
$\displaystyle=\frac{1}{J}\sum_{j=1}^{J}\lambda_{\min}(H_{j})-\gamma^{\prime}\geq-\frac{r_{t}}{J}\big{(}1+h^{-1}(t)\big{)}-\sqrt{\frac{2t}{J}}\beta_{t}$
(39) $\displaystyle\sum_{j=1}^{J}\lambda_{\min}(H_{j})$ $\displaystyle\geq-
r_{t}\big{(}1+h^{-1}(t)\big{)}-\sqrt{2tJ}\beta_{t}+J\gamma^{\prime},$ (40)

with probability at least $1-\exp(-t)$. Here, $h(t)=\frac{1}{2}(t-\log(1+t))$.
The two additional model-specific parameters are $\beta_{t}$ and $r_{t}$,
which do not depend on $J$.

We additionally require, for
$c\in\big{(}0,\sum_{j=1}^{J}\mathbb{E}\lambda_{\min}(H_{j})\big{)}$,

$\displaystyle\sum_{j=1}^{J}\lambda_{\min}(H_{j})\geq-
r_{t}\big{(}1+h^{-1}(t)\big{)}-\sqrt{2tJ}\beta_{t}+J\gamma^{\prime}\geq c,$
(41) $\displaystyle J\gamma^{\prime}\geq
c+r_{t}\big{(}1+h^{-1}(t)\big{)}+\sqrt{2tJ}\beta_{t}.$ (42)

It is therefore sufficient to have

$J>\max\left\\{\frac{2r_{t}\big{(}1+h^{-1}(t)\big{)}+2c}{\gamma^{\prime}},\frac{2(2t\beta_{t}^{2}+c)^{2}}{\gamma^{\prime
2}}\right\\}\geq\frac{r_{t}\big{(}1+h^{-1}(t)\big{)}}{\gamma^{\prime}}+\sqrt{2tJ}\beta_{t}/\gamma^{\prime}+c/\gamma^{\prime}$
(43)

for $\hat{H}(\theta)$ to be invertible and positive definite with minimum
eigenvalue greater than or equal to $c$ with probability at least
$1-\exp(-t)$.

∎

###  B.3 Convergence Analysis of Theorem 2

###### Proof.

In this analysis, we largely follow the proof of Theorem 6 in [40]. Define
$\theta_{\eta}=\theta_{m}+\eta p_{m}$, where
$p_{m}=-(H(\theta_{m}))^{-1}g(\theta_{m})$. As in Roosta-Khorasani and Mahoney
[40], we want to show there is some iteration-independent $\tilde{\eta}>0$
such that the Armijo condition

$f(\theta_{m}+\eta p_{m})\leq
f(\theta_{m})+\eta\beta\,p_{m}^{T}\,g(\theta_{m}),$ (44)

holds for any $0<\eta<\tilde{\eta}$ and some $\beta\in(0,1)$. By an argument
found in the beginning of the proof of Theorem 6 in [40], we have that
choosing $J$ such that
$\|\nabla_{\theta}\hat{\ell}(\theta_{m})-\nabla_{\theta}\ell(\theta_{m})\|\leq\epsilon$
and $\lambda_{\min}(H(\theta_{m}))\geq c>0$ for each $m$, yields

$\displaystyle f(\theta_{\eta})-f(\theta_{m})\leq\eta
p_{m}^{T}\,g(\theta_{m})+\epsilon\eta\|p_{m}\|+\eta^{2}\Gamma\|p_{m}\|^{2}/2,$
(45)

with probability $1-\delta/2$. From now on, we assume that we are on the
success event of this high-probability statement. Consequently, we have

$p_{m}^{T}g(\theta_{m})=-p_{m}^{T}\,H(\theta_{m})\,p_{m}\geq-c\|p_{m}\|^{2},$
(46)

and we can obtain a decrease in the objective. Substituting this into the
previous expression,

$\displaystyle f(\theta_{\eta})-f(\theta_{m})\leq-\eta
p_{m}^{T}H(\theta_{m})\,p_{m}+\epsilon\eta\|p_{m}\|+\eta^{2}\Gamma\|p_{m}\|^{2}/2,$
(47)

the Armijo condition becomes

$\displaystyle-\eta
p_{m}^{T}H(\theta_{m})p_{m}+\epsilon\eta\|p_{m}\|+\eta^{2}\Gamma\|p_{m}\|^{2}/2$
$\displaystyle\leq\eta\beta p_{m}^{T}g(\theta_{m})=-\eta\beta
p_{m}^{T}H(\theta_{m})p_{m}$ (48)
$\displaystyle\epsilon\|p_{m}\|+\eta\Gamma\|p_{m}\|^{2}/2$
$\displaystyle\leq(1-\beta)p_{m}^{T}H(\theta_{m})p_{m}$ (49)
$\displaystyle\epsilon+\eta\Gamma\|p_{m}\|/2$ $\displaystyle\leq
c(1-\beta)\|p_{m}\|.$ (50)

This holds and guarantees an iteration-independent lower bound if

$\eta\leq\frac{c(1-\beta)}{\Gamma},\;\;\epsilon\leq\frac{c(1-\beta)}{2\Gamma}\|g(\theta_{m})\|\leq\frac{c(1-\beta)}{2}\|p_{m}\|,$
(51)

which is given by our choice of $\eta$. Now, first note that

$\|g(\theta_{m})\|-\|\nabla_{\theta}f(\theta_{m})\|\leq\|g(\theta_{m})-\nabla_{\theta}f(\theta_{m})\|\leq\epsilon\implies\|\nabla_{\theta}f(\theta_{m})\|\geq\|g(\theta_{m})\|-\epsilon$
(52)

and

$\|\nabla_{\theta}f(\theta_{m})\|-\|g(\theta_{m})\|\leq\|\nabla_{\theta}f(\theta_{m})-g(\theta_{m})\|\leq\epsilon\implies\|g(\theta_{m})\|\geq\|\nabla_{\theta}f(\theta_{m})\|-\epsilon.$
(53)

There are now two cases. If the algorithm terminates and
$\|g(\theta_{m})\|\leq\sigma\epsilon$, we can derive

$\|\nabla_{\theta}f(\theta_{m})\|\leq\|g(\theta_{m})\|+\epsilon=\sigma\epsilon+\epsilon=(\sigma+1)\epsilon.$
(54)

If the algorithm does not terminate, then $\|g(\theta_{m})\|>\sigma\epsilon$.
Notice that

$\displaystyle\epsilon\geq\|g(\theta_{m})-\nabla_{\theta}f(\theta_{m})\|$
$\displaystyle\geq$
$\displaystyle\|g(\theta_{m})\|-\|\nabla_{\theta}f(\theta_{m})\|$ (55)
$\displaystyle\|\nabla_{\theta}f(\theta_{m})\|+\epsilon$ $\displaystyle\geq$
$\displaystyle\|g(\theta_{m})\|\geq\sigma\epsilon$ (56)
$\displaystyle\|\nabla_{\theta}f(\theta_{m})\|$ $\displaystyle\geq$
$\displaystyle\sigma\epsilon-\epsilon=(\sigma-1)\epsilon$ (57)
$\displaystyle\frac{\|\nabla_{\theta}f(\theta_{m})\|}{\sigma-1}$
$\displaystyle\geq$ $\displaystyle\epsilon,$ (58)

and now

$\displaystyle\|\nabla_{\theta}f(\theta_{m})\|-\epsilon$ $\displaystyle\geq$
$\displaystyle\|\nabla_{\theta}f(\theta_{m})\|-\frac{\|\nabla_{\theta}f(\theta_{m})\|}{\sigma-1}$
(59) $\displaystyle=$
$\displaystyle\left(1-\frac{1}{\sigma-1}\right)\|\nabla_{\theta}f(\theta_{m}\|$
(60) $\displaystyle=$
$\displaystyle\frac{\sigma-2}{\sigma-1}\|\nabla_{\theta}f(\theta_{m})\|$ (61)
$\displaystyle\geq$
$\displaystyle\frac{2}{3}\|\nabla_{\theta}f(\theta_{m})\|.$ (62)

Since $\|A^{-1}\|=1/\sigma_{\min}(A)$,

$\displaystyle p_{m}^{T}H(\theta_{m})p_{m}$
$\displaystyle=\big{(}-(H(\theta_{m}))^{-1}g(\theta_{m})\big{)}^{T}H(\theta_{m})\big{(}-(H(\theta_{m}))^{-1}g(\theta_{m})\big{)}$
(63) $\displaystyle=g(\theta_{m})^{T}(H(\theta_{m}))^{-1}g(\theta_{m})$ (64)
$\displaystyle\geq\frac{1}{c}\|g(\theta_{m})\|^{2}$ (65)
$\displaystyle\geq\frac{1}{c}\big{(}\|\nabla_{\theta}f(\theta_{m})\|-\epsilon\big{)}^{2}$
(66) $\displaystyle\geq\frac{4}{9c}\|\nabla_{\theta}f(\theta_{m})\|^{2}.$ (67)

From the assumption that $f$ is $\gamma$-strongly convex, $\gamma
I\preceq\nabla_{\theta}^{2}-\ell\preceq\Gamma I$, by an implication of
$\gamma$-strong convexity we have

$\displaystyle
f(\theta_{m})-f(\theta^{*})\leq\frac{1}{2\gamma}\big{\|}\nabla_{\theta}f(\theta_{m})\big{\|}^{2},$
(68)

and we put together:

$f(\theta_{m})-f(\theta^{*})\leq\frac{1}{2\gamma}\big{\|}\nabla_{\theta}f(\theta_{m})\big{\|}^{2}\leq\frac{9c}{4}\frac{1}{2\gamma}\,p_{m}^{T}H(\theta_{m})\,p_{m},$
(69)
$\frac{8\gamma}{9c}\big{(}f(\theta_{m})-f(\theta^{*})\big{)}\leq\frac{4}{9c}\,\big{\|}\nabla_{\theta}f(\theta_{m})\big{\|}^{2}\leq
p_{m}^{T}H(\theta_{m})\,p_{m},$ (70)
$f(\theta_{m})-f(\theta^{*})\leq\frac{9c}{8\gamma}\,p_{m}^{T}H(\theta_{m})\,p_{m},$
(71)
$-\frac{8\gamma}{9c}\big{(}f(\theta_{m})-f(\theta^{*})\big{)}\geq-\frac{4}{9c}\,\big{\|}\nabla_{\theta}f(\theta_{m})\big{\|}^{2}\geq-
p_{m}^{T}H(\theta_{m})\,p_{m}.$ (72)

From earlier, as the Armijo condition is fulfilled with our choice of $\eta$
and $\epsilon$,

$\displaystyle f(\theta_{m+1})-f(\theta_{m})$
$\displaystyle\leq-\eta\,p_{m}^{T}H(\theta_{m})p_{m}+\epsilon\eta\|p_{m}\|+\eta^{2}\Gamma\|p_{m}\|^{2}/2$
(73) $\displaystyle\leq-\eta\beta\,p_{m}^{T}H(\theta_{m})\,p_{m}$ (74)
$\displaystyle\leq-\eta\beta\,\frac{8\gamma}{9c}\,\big{(}f(\theta_{m})-f(\theta^{*})\big{)}.$
(75)

Therefore,

$\displaystyle f(\theta_{m+1})-f(\theta^{*})$
$\displaystyle=f(\theta_{m+1})-f(\theta_{m})+f(\theta_{m})-f(\theta^{*})$ (76)
$\displaystyle\leq
f(\theta_{m})-f(\theta^{*})-\eta\beta\,\frac{8\gamma}{9c}\,\big{(}f(\theta_{m})-f(\theta^{*})\big{)}$
(77)
$\displaystyle=\Big{(}1-\eta\beta\frac{8\gamma}{9c}\Big{)}\big{(}f(\theta_{m})-f(\theta^{*})\big{)}.$
(78)

∎

##  Appendix C Feynman-Kac Models and Monte Carlo Approximations

In this section, we introduce the Feynman-Kac convention of [7] that has since
become commonplace [42] for the analysis of the particle filter. The
mathematical formalization and notation introduced here will be adopted in the
remainder of the analysis, in order to prove Theorems 3, 4, and 5. Let
$(\eta_{n})_{n=1}^{N},(\pi_{n})_{n=1}^{N},(\rho_{n})_{n=1}^{N}$ be sequences
of probability measures on the state space ${\mathcal{X}}$. This is the
sequence of prediction distributions $f_{X_{n}|Y_{1:n-1}}$, filtering
distributions $f_{X_{n}|Y_{1:n}}$, and posterior distributions
$f_{X_{1:n}|Y_{1:n}}$ that we seek to approximate with the particle filter.
For any measurable bounded functional $h$, we adopt the following functional-
analytic notation, borrowed from [7, 47, 42]. We choose our specific choice of
notation and definitions to be in line with that of [42].

#### Markov kernels and the process model:

A Markov kernel $M$ with source ${\mathcal{X}}_{1}$ and target
${\mathcal{X}}_{2}$ is a map
$M:{\mathcal{X}}_{1}\times\mathcal{B}({\mathcal{X}}_{2})\to[0,1]$ such that
for every set $A\in\mathcal{B}({\mathcal{X}}_{2})$ and every point
$x\in{\mathcal{X}}_{1}$, the map $x\mapsto M(x,A)$ is a measurable function of
$x$, and the map $A\mapsto M(x,A)$ is a probability measure on
${\mathcal{X}}_{2}$. The quantity $M(x,A)$ can be thought of as the
probability of transitioning to the set $A$ given that we are at the point
$x$. If this yields a density, this then corresponds to the process density
$f_{X_{2}|X_{1}}$ conditional on $x$ and integrated over $A$.

#### Markov kernels and measures:

For any measure $\eta$, any Markov kernel $M$ on ${\mathcal{X}}$, any point
$x\in{\mathcal{X}}$ and any measurable subset $A\subseteq{\mathcal{X}}$, let

$\displaystyle\eta(h)$ $\displaystyle=\int h\,d\eta=\int h(x)\eta(dx),$ (79)
$\displaystyle(\eta M)(A)$ $\displaystyle=\int\eta(dx)M(x,A),$ (80)
$\displaystyle(Mh)(x)$ $\displaystyle=\int M(x,dy)h(y).$ (81)

#### Compositions of Markov kernels:

The composition of a Markov kernel $M_{1}$ with another Markov kernel $M_{2}$
is another Markov kernel, given by

$(M_{1}M_{2})(x,A)=\int M_{1}(x,dy)M_{2}(y,A).$ (82)

#### Total variation distance:

The total variation distance between two measures $\mu$ and $\nu$ on
${\mathcal{X}}$ is

$\|\mu-\nu\|_{\mathrm{TV}}=\sup_{\|h\|_{\infty}\leq
1/2}|\mu(\phi)-\nu(\phi)|=\sup_{\operatorname{osc}(h)\leq 1}|\mu(h)-\nu(h)|.$
(83)

#### Dobrushin contraction:

The Dobrushin contraction coefficient $\beta_{\text{TV}}$ of a Markov kernel
$M$ is given by

$\beta_{\mathrm{TV}}(M)=\sup_{x,y\in{\mathcal{X}}}\big{\|}M(x,\cdot)-M(y,\cdot)\big{\|}_{\mathrm{TV}}=\sup_{\mu,\nu\in\mathcal{P},\mu\neq\nu}\frac{\|\mu
M-\nu M\|_{\mathrm{TV}}}{\|\mu-\nu\|_{\mathrm{TV}}}.$ (84)

#### Potential functions and the measurement model:

A potential function $G:{\mathcal{X}}\to[0,\infty)$ is a non-negative function
of an element of the state space $x\in{\mathcal{X}}$. In our case, this
corresponds to the measurement model, and in our previous notation is written
as $g_{n,j}=f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F})=G_{n}(x_{n,j}^{F})$, where
in a slight abuse of notation we suppress the dependence on $\theta$ for
notational simplicity. Note that
$G_{n}(\cdot)=f_{Y_{n}|X_{n}}(y_{n}^{*}|\;\cdot\;)$ is the conditional density
of the observed measurement at time $n$, where we condition on the filtering
particle $x_{n,j}^{F}$ as an element of the state space.

#### Feynman-Kac models:

A Feynman-Kac model on ${\mathcal{X}}$ is a tuple
$(\pi_{0},(M_{n})_{n=1}^{N},(G_{n})_{n=1}^{N})$ of an initial probability
measure on the state space $\pi_{0}$, a sequence of transition kernels
$(M_{n})_{n=1}^{N}$, and a sequence of potential functions
$(G_{n})_{n=1}^{N}$. In the notation used in the main text, this corresponds
to the starting distribution $f_{X_{0}}$, the sequence of transition densities
$f_{X_{n}|X_{n-1}}$, and the measurement densities $f_{Y_{n}|X_{n}}$. This
induces a set of mappings from the set of probability measures on
${\mathcal{X}}$ to itself,
$\mathcal{P}({\mathcal{X}})\to\mathcal{P}({\mathcal{X}})$, as follows:

  * •

The update from the prediction to the filtering distributions is given by

$\pi_{n}(dx)=\Psi_{n}(\eta_{n})(dx)=\frac{G_{n}(x)\cdot\eta_{n}(dx)}{\eta_{n}(G_{n})}.$
(85)

  * •

The map from the prediction distribution at timestep $n$ to timestep $n+1$ is
given by

$\Phi_{n+1}(\eta_{n})=\Psi_{n}(\eta_{n})M_{n+1}.$ (86)

  * •

The composition of maps between prediction distributions yields the map from
the prediction distribution at time $k$ to the prediction distribution at time
$n$ where $k\leq n$,

$\Phi_{k,n}=\Phi_{n}\circ...\circ\Phi_{k+1}.$ (87)

#### The particle filter:

The particle filter then yields a Monte Carlo approximation to the above
Feynman-Kac model, via a sequence of mixture Dirac measures. When one
resamples at every timestep, the prediction measure at timestep $n$ is then
given by

$\eta_{n}^{J}=\frac{1}{J}\sum_{j=1}^{J}\delta_{x_{n,j}^{P}},$ (88)

and the filtering measure at timestep $n$ is given by

$\pi_{n}^{J}=\frac{\sum_{j=1}^{J}g_{n,j}\delta_{x_{n,j}^{P}}}{\sum_{j=1}^{J}g_{n,j}}\approx\frac{1}{J}\sum_{j=1}^{J}\delta_{x_{n,j}^{F}}.$
(89)

In a slight abuse of notation, we will identify
$x_{n,1:J}^{P}\equiv\eta_{n}^{J}$, and $x_{n,1:J}^{F}\equiv\pi_{n}^{J}$. As in
[42], one can view this as an inhomogenous Markov process evolving on
${\mathcal{X}}^{J}$. The corresponding Markov transition kernel is then

$\textbf{M}_{n}(x_{n-1,1:J}^{P},\cdot)=\left(\Phi_{n}\left(\eta_{n}^{J}\right)\right)^{\otimes
J}=\left(\Phi_{n}\left(\frac{1}{J}\sum_{j=1}^{J}\delta_{x_{n,j}^{P}}\right)\right)^{\otimes
J},$ (90)

and the composition of Markov kernels on particles from timestep $n$ to
timestep $n+k$ is written

$\textbf{M}_{n,n+k}=\textbf{M}_{n+k}\circ...\circ\textbf{M}_{n}.$ (91)

One may wonder why [42] require this process to evolve on ${\mathcal{X}}^{J}$.
This is because at every timestep $n$, we in fact draw
$X_{n,j}^{P}|\\{X_{n-1,1:J}^{P}=x_{n-1,1:J}^{P}\\}\sim\textbf{M}_{n}(x_{n-1,1:J}^{P},\cdot)=\eta_{0}^{\otimes
J}\textbf{M}_{0,n}$ for $j=1,...,J$.

#### Forgetting of the particle filter:

The above formalization yields a result from [42] on the forgetting of the
particle filter that we require for our analysis of the bias, variance, and
error of MOP-$\alpha$. That is, [42] show that

$\beta_{\mathrm{TV}}\left(\mathbf{M}_{n,n+k}\right)\leq(1-\epsilon)^{\lfloor
k/(O(\log J))\rfloor},$ (92)

for some $\epsilon$ dependent on $\bar{G},\underline{G},\bar{M},\underline{M}$
in Assumptions (A3) and (A2). As a result, the mixing time of the particle
filter is only on the order of $O(\log(J))$ timesteps.

Equipped with the above formalisms and results, we are now in a position to
provide guarantees on the performance of MOP-$\alpha$ itself.

##  Appendix D A Strong Law of Large Numbers for Triangular Arrays of
Particles With Off-Parameter Resampling

Here, we prove a more general result, from which it will be clear that
MOP-$\alpha$ targets the filtering distribution. We will prove a strong law of
large numbers for triangular arrays of particles with off-parameter
resampling, meaning that we resample the particles according to an arbitrary
resampling rule that is not necessarily in proportion to the target
distribution of interest. Using weights that encode the cumulative discrepancy
between the resampling distribution and the target distribution (instead of
resampling to equal weights, as in the basic particle filter) provides a
sufficient correction to ensure almost sure convergence. We now introduce the
precise definition of an off-parameter resampled particle filter.

######  Definition 1 (Off-Parameter Resampled Particle Filters).

An off-parameter resampled particle filter is a Monte Carlo approximation to a
Feynman-Kac model
$\left(\pi_{0},\left(M_{n}\right)_{n=1}^{N},\left(G_{n}\right)_{n=1}^{N}\right)$,
where we inductively define given some $x_{n-1,j}^{F}$, $w_{n-1,j}^{F}$
comprising the filtering measure approximation $\pi_{n}^{J}$:

  1. 1.

The prediction particles at timestep $n$, $x_{n,j}^{P}$, are drawn from
$X_{n,j}^{P}\sim\pi_{n}^{J}M_{n}$, and the prediction weights are given by
$w_{n,j}^{P}=w_{n-1,j}^{F}$.

  2. 2.

The prediction measure at timestep $n$ is given by

$\eta_{n}^{J}=\frac{\sum_{j=1}^{J}w_{n,j}^{P}\,\delta_{x_{n,j}^{P}}}{\sum_{j=1}^{J}w_{n,j}^{P}}.$
(93)

  3. 3.

The particles are resampled at every timestep $n$, yielding indices $k_{j}$
according to some arbitrary probabilities $(p_{n,j})_{j=1}^{J}$.

  4. 4.

the filtering measure at timestep $n$ is given by

$\pi_{n}^{J}=\frac{\sum_{j=1}^{J}\delta_{x_{n,j}^{P}}w_{n,j}^{P}\,g_{n,j}}{\sum_{j=1}^{J}w_{n,j}^{P}\,g_{n,j}}\approx\frac{\sum_{j=1}^{J}\delta_{x_{n,k_{j}}^{P}}w_{n,k_{j}}^{P}\,g_{n,k_{j}}\big{/}p_{n,k_{j}}}{\sum_{j=1}^{J}w_{n,k_{j}}^{P}\,g_{n,k_{j}}\big{/}p_{n,k_{j}}}=\frac{\sum_{j=1}^{J}w_{n,j}^{F}\,\delta_{x_{n,j}^{F}}}{\sum_{j=1}^{J}w_{n,j}^{F}}.$
(94)

  5. 5.

The prediction weights at the next timestep are given by
$w_{n+1,j}^{P}=w_{n,j}^{F}=w_{n,k_{j}}^{P}\,g_{n,k_{j}}\big{/}p_{n,k_{j}}$.

To prove that the weight correction is sufficient for almost sure convergence,
we first introduce what it means for a triangular array of particles to target
a given target distribution.

######  Definition 2 (Targeting).

A pair of random vectors $(X,W)$ drawn from some measure $g$ targets another
measure $\pi$ if for any measurable and bounded functional $h$,

$E_{g}\big{[}h(X)\cdot W\big{]}=E_{\pi}\big{[}h(X)\big{]}.$ (95)

A set of particles targeting $\pi$ is a triangular array of pairs of random
vectors $(X^{J}_{j},W^{J}_{j}),j=1,2,\ldots,J$ such that for any measurable
and bounded functional $h$,

$\frac{\sum_{j=1}^{J}h(X^{J}_{j})\,W^{J}_{j}}{\sum_{j=1}^{J}W^{J}_{j}}\stackrel{{\scriptstyle
a.s.}}{{\to}}E_{\pi}(h(X))$ (96)

as $J\to\infty$.

[48] asserted without proof that common particle filter algorithms targets the
filtering distribution in this sense, while [47] proved a related result
assuming bounded densities. We follow a similar approach to [47], based on
showing strong laws of large numbers for triangular arrays, noting that
triangular array strong laws do not hold without an additional regularity
condition such as boundedness in general.

In order to prove the consistency of our variation on the particle filter, we
now present three helper lemmas. The first follows from standard importance
sampling arguments, the second from integrating out the marginal, and the
third from Bayes’ theorem. We state Lemma 5 assuming multinomial resampling,
which is convenient for the proof though other resampling strategies may be
preferable in practice.

######  Lemma 5 (Change of Weight Measure).

Suppose that $\\{(\tilde{X}_{j}^{J},U_{j}^{J}),j=1,\dots,J\\}$ targets
$f_{X}$. Now, let $\\{(Y_{j}^{J},V_{j}^{J}),j=1,\dots,J\\}$ be a multinomial
sample with indices $k_{j}$ drawn from $\\{(\tilde{X}_{j}^{J},U_{j}^{J})\\}$
where $(\tilde{X}_{j}^{J},U_{j}^{J})$ is represented, on average, proportional
to $\pi^{J}_{j}J$ times. Write

$(Y_{j}^{J},V_{j}^{J})=\big{(}\tilde{X}^{J}_{k_{j}},U^{J}_{k_{j}}\big{/}\pi^{J}_{k_{j}}\big{)}.$
(97)

If the importance sampling weights $U_{j}/\pi_{j}$ are bounded, then
$\\{(Y^{J}_{j},V^{J}_{j}),j=1,\dots,J\\}$ targets $f_{X}$.

###### Proof.

Note that as the $Y_{j}^{J}$ are a subsample from $X_{j}^{J}$, $h$ can be a
function of $Y$ as well as it is one for $X$. We then expand

$\frac{\sum_{j}h(Y_{j}^{J})\,V_{j}^{J}}{\sum_{j}V_{j}^{J}}=\frac{\sum_{j}h(\tilde{X}_{k_{j}}^{J})\,\frac{U_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}}{\sum_{j}\frac{U_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}}.$
(98)

By hypothesis,

$\frac{\sum_{j}h(\tilde{X}_{j}^{J})\,U_{j}^{J}}{\sum_{j}U_{j}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{X}}\big{[}f(X)\big{]}.$ (99)

We want to show

$\frac{\sum_{j}h(X_{k_{j}}^{J})\,\frac{U_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}}{\sum_{j}\frac{U_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}}-\frac{\sum_{j}h(\tilde{X}_{j}^{J})\,U_{j}^{J}}{\sum_{j}U_{j}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}0.$ (100)

For this, it is sufficient to show that

$\sum_{j}h(\tilde{X}_{k_{j}}^{J})\frac{U_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}-\sum_{j}h(\tilde{X}_{j}^{J})\frac{U_{j}^{J}}{\pi_{j}^{J}}\pi_{j}^{J}\stackrel{{\scriptstyle
a.s.}}{{\to}}0$ (101)

since an application of this result with $h(x)=1$ provides almost sure
convergence of the denominator in (100). Write
$g(\tilde{X}_{j}^{J})=h(\tilde{X}_{j}^{J})\frac{U_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}$.
We therefore need to show that

$\sum_{j}Z_{j}^{J}:=\sum_{j}\left(g(\tilde{X}_{k_{j}}^{J})-g(\tilde{X}_{j}^{J})\,\pi_{j}^{J}\right)\stackrel{{\scriptstyle
a.s.}}{{\to}}0.$ (102)

Because the functional $h$ and importance sampling weights
$u_{k_{j}}^{J}/\pi_{k_{j}}^{J}$ are bounded, we have that
$\mathbb{E}\big{[}(Z_{j}^{J})^{4}\big{]}<\infty$. We can then follow the
argument of [47] from this point on, where noting that the $Z_{j}^{J}$ are
conditionally independent given the $\tilde{X}_{j}^{J}$ and $U_{j}^{J}$,

$\displaystyle\mathbb{E}\left[\left(\sum_{j}Z_{j}^{J}\right)^{4}\Bigg{|}(\tilde{X}_{j}^{J},U_{j}^{J})_{j=1}^{J}\right]$
$\displaystyle=$ $\displaystyle
J\,\mathbb{E}\left[(Z_{1}^{J})^{4}|(\tilde{X}_{j}^{J},U_{j}^{J})_{j=1}^{J}\right]$
(103) $\displaystyle\hskip 28.45274pt+\hskip
5.69054pt3J(J-1)\left(\mathbb{E}\big{[}(Z_{j}^{J})^{2}\big{|}(\tilde{X}_{j}^{J},U_{j}^{J})_{j=1}^{J}\big{]}\right)$
$\displaystyle\leq$ $\displaystyle CJ^{2},$ (104)

for some $C>0$. Taking expectations on both sides yields

$\mathbb{E}\left[\left(\sum_{j}Z_{j}^{J}\right)^{4}\right]\leq CJ^{2}$ (105)

by the tower property. Now by Markov’s inequality,

$\mathbb{P}\left(\left|\frac{1}{J}\sum_{j=1}^{J}Z_{j}^{J}\right|>\epsilon\right)\leq\frac{1}{\epsilon^{4}J^{4}}\mathbb{E}\left[\left(\sum_{j=1}^{J}Z_{j}^{J}\right)^{4}\right]\leq\frac{C}{\epsilon^{4}J^{2}},$
(106)

and as these terms are summable we can apply Borel-Cantelli to conclude that
these deviations happen only finitely often for every $\epsilon>0,$ giving us
the almost-sure convergence for

$\sum_{j}h(X_{k_{j}}^{J})\,\frac{u_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}-\sum_{j}h(X_{j}^{J})\,\frac{u_{j}^{J}}{\pi_{j}^{J}}\,\pi_{j}^{J}\stackrel{{\scriptstyle
a.s.}}{{\to}}0.$ (107)

Similarly, we also have that

$\sum_{j}\frac{u_{j}^{J}}{\pi_{j}^{J}}\pi_{j}^{J}-\sum_{j}\frac{u_{k_{j}}^{J}}{\pi_{k_{j}}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}0,$ (108)

and the result is proved.

∎

Remark: Note that Lemma 5 permits $\pi_{1:J}$ to depend on
$\\{(X_{j},u_{j})\\}$ as long as the resampling is carried out independently
of $\\{(X_{j},u_{j})\\}$, conditional on $\pi_{1:J}$.

######  Lemma 6 (Particle Marginals).

Suppose that $\\{(\tilde{X}_{j}^{J},U_{j}^{J}),j=1,\dots,J\\}$ targets
$f_{X}$. Also suppose that $\tilde{Z}_{j}^{J}\sim
f_{Z|X}(\cdot|\tilde{X}_{j}^{J})$ where $f_{Z|X}$ is a conditional probability
density function corresponding to a joint density $f_{X,Z}$ with marginal
densities $f_{X}$ and $f_{Z}$. Then, if the $U_{j}^{J}$ are bounded,
$\\{(\tilde{Z}_{j}^{J},U_{j}^{J})\\}$ targets $f_{Z}$.

###### Proof.

We want to show that, for any measurable bounded $h$,

$\frac{\sum_{j}h(\tilde{Z}_{j}^{J})\,U_{j}^{J}}{\sum_{j}U_{j}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{Z}}[h(Z)]=\mathbb{E}_{f_{X}}\big{[}\mathbb{E}_{f_{Z|X}}[h(Z)|X]\big{]}.$
(109)

By assumption, for any measurable and bounded functional $g$ with domain
${\mathcal{X}}$,

$\frac{\sum_{j}g(\tilde{X}_{j}^{J})\,U_{j}^{J}}{\sum_{j}U_{j}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{X}}[g(X)].$ (110)

Let $\bar{U}_{j}^{J}=\frac{J\,U_{j}^{J}}{\sum_{j}U_{j}^{J}}$. Examine the
numerator and denominator of the quantity

$\frac{J^{-1}\sum_{j}h(\tilde{Z}_{j}^{J})\,U_{j}^{J}}{J^{-1}\sum_{j}U_{j}^{J}}=\frac{J^{-1}\sum_{j}h(\tilde{Z}_{j}^{J})\,\bar{U}_{j}^{J}}{J^{-1}\sum_{j}\bar{U}_{j}^{J}}.$
(111)

The denominator converges to $1$ almost surely. The numerator, on the other
hand, is

$\frac{1}{J}\sum_{j}h(\tilde{Z}_{j}^{J})\,\bar{U}_{j}^{J},$ (112)

and by the same fourth moment argument to the above lemma, it converges almost
surely to the limit of its expectation,

$\displaystyle\lim_{J\to\infty}\mathbb{E}\left[\frac{1}{J}\sum_{j}h(\tilde{Z}_{j}^{J})\,\bar{U}_{j}^{J}\right]$
$\displaystyle=$
$\displaystyle\lim_{J\to\infty}\mathbb{E}\left[\frac{1}{J}\sum_{j}\mathbb{E}\Big{[}h(\tilde{Z}_{j}^{J})\,\bar{U}_{j}^{J}\Big{|}\tilde{X}_{j}^{J},\bar{U}_{j}^{J}\Big{]}\right]$
(113) $\displaystyle=$
$\displaystyle\lim_{J\to\infty}\mathbb{E}\left[\frac{1}{J}\sum_{j}\mathbb{E}\Big{[}h(Z)\Big{|}X=\tilde{X}_{j}^{J}\Big{]}\,\bar{U}_{j}^{J}\right].$
(114)

Applying (110) with $g(x)=\mathbb{E}\big{[}h(Z)|X=x\big{]}$, the average on
the right hand side converges almost surely to
$\mathbb{E}\big{\\{}\mathbb{E}[h(Z)|X]\big{\\}}=\mathbb{E}[h(Z)].$ It remains
to swap the limit and expectations. We can do so with the bounded convergence
theorem, and therefore obtain

$\frac{1}{J}\sum_{j}h(\tilde{Z}_{j}^{J})\,\bar{U}_{j}^{J}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{Z}}[h(Z)].$ (115)

∎

######  Lemma 7 (Particle Posteriors).

Suppose that $\\{(X_{j}^{J},U_{j}^{J}),j=1,\dots,J\\}$ targets $f_{X}$. Also
suppose that $(X^{\prime J}_{j},U^{\prime
J}_{j})=\big{(}X_{j}^{J},U_{j}^{J}\,f_{Z|X}(z^{*}|X_{j}^{J})\big{)}$. Then, if
$U_{j}^{J}\,f_{Z|X}(z^{*}|X_{j}^{J})$ and
$U_{j}^{J}\,f_{Z|X}(z^{*}|X_{j}^{J})\big{/}f_{Z}(z^{*})$ are bounded,
$\\{(X^{\prime J}_{j},U^{\prime J}_{j})\\}$ targets $f_{X|Z}(\cdot|z^{*})$.

###### Proof.

Again, we want to show that

$\frac{\sum_{j}h(X_{j}^{J})\cdot U_{j}^{J}\cdot
f_{Z|X}(z^{*}|X_{j}^{J})}{\sum_{j}U_{j}^{J}\cdot
f_{Z|X}(z^{*}|X_{j}^{J})}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{X|Z}}\big{[}h(X)\big{|}z^{*}\big{]}.$ (116)

We already have that for any measurable bounded $g$,

$\frac{\sum_{j}g(X_{j}^{J})\,U_{j}^{J}}{\sum_{j}U_{j}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{X}}\big{[}g(X)\big{]}.$ (117)

Consider the following:

$\frac{J^{-1}\sum_{j}h(X_{j}^{J})\,f_{Z|X}(z^{*}|X_{j}^{J})\,{U}_{j}^{J}}{J^{-1}\sum_{j}{U}_{j}^{J}}\times\left(\frac{J^{-1}\sum_{j}f_{Z|X}(z^{*}|X_{j}^{J})\,{U}_{j}^{J}}{J^{-1}\sum_{j}{U}_{j}^{J}}\right)^{-1}.$
(118)

We will apply Equation (117) to the numerator and the denominator in the ratio
above individually. The numerator converges to

$\frac{J^{-1}\sum_{j}h(X_{j}^{J})\,f_{Z|X}(z^{*}|X_{j}^{J})\,{U}_{j}^{J}}{J^{-1}\sum_{j}{u}_{j}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{X}}\big{[}h(X)\,f_{Z|X}(z^{*}|X)\big{]},$ (119)

while the reciprocal of the denominator converges to

$\frac{J^{-1}\sum_{j}f_{Z|X}(z^{*}|X_{j}^{J})\,{U}_{j}^{J}}{J^{-1}\sum_{j}{U}_{j}^{J}}\stackrel{{\scriptstyle
a.s.}}{{\to}}\mathbb{E}_{f_{X}}\big{[}f_{Z|X}(z^{*}|X)\big{]}=f_{Z}(z^{*}).$
(120)

Now we take advantage of the identities

$\displaystyle\frac{\mathbb{E}_{f_{X}}\big{[}h(X)\,f_{Z|X}(z^{*}|X)\big{]}}{f_{Z}(z^{*})}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{f_{X}}\left[\frac{h(X)\,f_{Z|X}(z^{*}|X)}{f_{Z}(z^{*})}\right]$
(121) $\displaystyle=$
$\displaystyle\mathbb{E}_{f_{X}}\left[h(X)\frac{f_{X|Z}(X|z^{*})}{f_{X}(X)}\right]=\mathbb{E}_{f_{X|Z}}\big{[}h(X)|z^{*}\big{]},$
(122)

to give the desired result.

∎

######  Theorem 6 (Off-Parameter Resampled Particle Filters Target the
Filtering Distribution).

The off-parameter resampled particle filter as outlined in Definition 1
targets the filtering distribution.

###### Proof.

Recursively applying Lemmas 5, 6, and 7, we obtain that the off-parameter
resampled particle filter targets the posterior. Specifically, suppose
inductively that $\big{\\{}\big{(}X^{F}_{n-1,j},w^{F}_{n-1,j}\big{)}\big{\\}}$
targets $\pi_{n-1}$. Then, Lemma 6 tells us that
$\big{\\{}\big{(}X^{P}_{n,j},w^{P}_{n,j}\big{)}\big{\\}}$ targets $\eta_{n}$.
Lemma 7 tells us that
$\big{\\{}\big{(}X^{P}_{n,j},w^{P}_{n,j}g^{\theta}_{n,j}\big{)}\big{\\}}$
therefore targets $\pi_{n}$. Lemma 5 guarantees that the resampling rule,
given by

$\big{(}X^{F}_{n,j},w^{F}_{n,j}\big{)}=\big{(}X^{P}_{n,k_{j}},w^{P}_{n,k_{j}}\,g_{n,k_{j}}\big{/}p_{n,k_{j}}\big{)},$
(123)

with resampling probabilities proportional to $p_{n,j}$, therefore also
targets $\pi_{n}$. ∎

######  Proposition 1 (MOP-1 Targets the Filtering Distribution).

When $\alpha=1$ or $\phi=\theta$, MOP-$\alpha$ targets the filtering
distribution.

###### Proof.

When $\theta=\phi$, regardless of the value of $\alpha$, the ratio
$\frac{g_{n,j}^{\theta}}{g_{n,j}^{\phi}}=1,$ and this reduces to the vanilla
particle filter estimate.

When $\alpha=1$, and $\theta\neq\phi,$ the proof is identical to that of 6.
Recursively applying Lemmas 5, 6, and 7, we obtain that the MOP-1 filter
targets the posterior. Specifically, suppose inductively that
$\big{\\{}\big{(}X^{F,\theta}_{n-1,j},w^{F,\theta}_{n-1,j}\big{)}\big{\\}}$ is
properly weighted for $f_{X_{n-1}|Y_{1:n-1}}(x_{n-1}|y^{*}_{1:n-1};\theta)$.
Then, Lemma 6 tells us that
$\big{\\{}\big{(}X^{P,\theta}_{n,j},w^{P,\theta}_{n,j}\big{)}\big{\\}}$
targets $f_{X_{n}|Y_{1:n-1}}(x_{n}|y^{*}_{1:n-1};\theta)$. Lemma 7 tells us
that
$\big{\\{}\big{(}X^{P,\theta}_{n,j},w^{P,\theta}_{n,j}\,g^{\theta}_{n,j}\big{)}\big{\\}}$
therefore targets $f_{X_{n}|Y_{1:n}}(x_{n}|y^{*}_{1:n};\theta)$. Lemma 5
guarantees that the resampling rule, given by

$\big{(}X^{F,\theta}_{n,j},w^{F,\theta}_{n,j}\big{)}=\big{(}X^{P,\theta}_{n,k_{j}},w^{P,\theta}_{n,k_{j}}\,g^{\theta}_{n,k_{j}}\big{/}g^{\phi}_{n,k_{j}}\big{)},$
(124)

with resampling weights proportional to $g^{\phi}_{n,j}$, therefore also
targets $f_{X_{n}|Y_{1:n}}(x_{n}|y^{*}_{1:n};\theta)$. ∎

This has addressed filtering, but not quite yet the likelihood evaluation. For
this we use the following lemma.

######  Lemma 8 (Likelihood Proper Weighting).

$f_{Y_{n}|Y_{1:n-1}}(y_{n}^{*}|y_{1_{n}-1}^{*};\theta)$ is consistently
estimated by either the before-resampling estimate,

$L_{n}^{B,\theta}=\frac{\sum_{j=1}^{J}g^{\theta}_{n,j}\,w^{P,\theta}_{n,j}}{\sum_{j=1}^{J}w^{P,\theta}_{n,j}},$
(125)

or by the after-resampling estimate,

$L_{n}^{A,\theta}=L^{\phi}_{n}\frac{\sum_{j=1}^{J}w^{F,\theta}_{n,j}}{\sum_{j=1}^{J}w^{P,\theta}_{n,j}}.$
(126)

where $L^{\phi}_{n}$ is as defined in the various algorithms.

Here, (125) is a direct consequence of our earlier result that
$\\{\big{(}X^{P,\theta}_{n,j},w^{P,\theta}_{n,j}\big{)}\\}$ targets
$f_{X_{n}|Y_{1:n-1}}(x_{n}|y^{*}_{1:n-1};\theta)$. To see (126), we write the
numerator of (5) as

$L^{\phi}_{n}\sum_{j=1}^{J}\left[\frac{g^{\theta}_{n,j}}{g^{\phi}_{n,j}}\,w^{P,\theta}_{n,j}\right]\frac{g^{\phi}_{n,j}}{L_{n}^{\phi}}=L^{\phi}_{n}\sum_{j=1}^{J}w_{n,j}^{FC,\theta}\,\frac{g^{\phi}_{n,j}}{L_{n}^{\phi}}$
(127)

Using Lemma 5, we resample according to probabilities
$\frac{g^{\phi}_{n,j}}{L_{n}^{\phi}}$ to see this is properly estimated by

$L^{\phi}_{n}\sum_{j=1}^{J}w^{F,\theta}_{n,j},$ (128)

from which we obtain (126). Using Lemma 8, we obtain a likelihood estimate,

$L^{A,\theta}=\prod_{n=1}^{N}\left(L^{\phi}_{n}\,\frac{\sum_{j=1}^{J}w^{F,\theta}_{n,j}}{\sum_{j=1}^{J}w^{P,\theta}_{n,j}}\right).$
(129)

Since $w^{F,\theta}_{n,j}=w^{P,\theta}_{n+1,j}$, this is a telescoping
product. The remaining terms are $\sum_{j=1}^{J}w^{P,\theta}_{0,j}=J$ on the
denominator and $\sum_{j=1}^{J}w^{F,\theta}_{N,j}$ on the numerator. This
derives the MOP likelihood estimates.

$L^{B,\theta}$ should generally be preferred in practice, since there is no
reason to include the extra variability from resampling when calculating the
conditional log likelihood, but it lacks the nice telescoping product that
lets us derive exact expressions for the gradient in Theorem 1.

##  Appendix E Consistency of Off-Parameter Resampled Gradient Estimates

We now provide a result showing that under sufficient regularity conditions,
one can interchange the order of differentiation and expectation of
functionals of off-parameter resampled particle estimates, showing that if an
off-parameter resampled particle estimate is consistent for some estimand, its
derivative is consistent for the derivative of the estimand as well.

One may wonder why this may be necessary, given that we have already shown
strong consistency for measurable bounded functionals in Theorem 6. If we
require the gradient to be bounded as well, then the gradient is also a
bounded functional. The answer is that the strong consistency is for
expectations under the filtering distribution, so Theorem 6 only establishes
strong consistency of the gradient of the estimate to its expectation under
the filtering distribution,
$\nabla_{\theta}\eta_{n}^{J}(h_{\theta})\stackrel{{\scriptstyle
a.s.}}{{\to}}\eta_{n}(\nabla_{\theta}h_{\theta})$, which is not a-priori the
gradient of the estimand $\nabla_{\theta}\eta_{n}(h_{\theta})$. We require an
interchange of the derivative and expectation, which we show is possible when
the particle estimates have two continuous uniformly bounded derivatives over
all $J$ below.

######  Theorem 7 (Off-Parameter Particle Filters Yield Strongly Consistent
Estimates of Derivatives of Functionals).

Let $h_{\theta}:{\mathcal{X}}\to\mathbb{R}$ be a measurable bounded functional
of particles, where $\eta_{n}^{J}(h_{\theta})$ has two continuous derivatives
uniformly bounded over all $J$ by $H^{*}$ for almost every $\omega\in\Omega$
and $\theta\in\Theta$. If it holds that
$\eta_{n}^{J}(h_{\theta})\stackrel{{\scriptstyle
a.s.}}{{\to}}\eta(h_{\theta})=h^{*}_{\theta}$, then we also have that
$\nabla_{\theta}\eta_{n}^{J}(h_{\theta})\stackrel{{\scriptstyle
a.s.}}{{\to}}\eta_{n}(\nabla_{\theta}h_{\theta})=\nabla_{\theta}\eta_{n}(h_{\theta})=\nabla_{\theta}h^{*}_{\theta}$.

###### Proof.

Fix $\omega\in\Omega$. The sequence
$(\nabla_{\theta}\eta_{n}^{J}(h_{\theta})(\omega))_{J\in\mathbb{N}}$ is
uniformly bounded over all $J$, by assumption. The sequence is also uniformly
equicontinuous. To see this, by assumption, the second derivative of
$\eta_{n}^{J}(h_{\theta})(\omega)|_{\theta=\theta^{\prime}}$ is also uniformly
bounded over all $J$ for almost every $\omega\in\Omega$ and every
$\theta^{\prime}\in\Theta$. A set of functions with derivatives bounded by the
same constant is uniformly Lipschitz, and therefore uniformly equicontinuous.
So the sequence $(\eta_{n}^{J}(h_{\theta})\omega))_{J\in\mathbb{N}}$ is
uniformly equicontinuous over $\theta$ for almost every $\omega\in\Omega$.
Explicitly, for almost every $\omega\in\Omega$ and every $\epsilon>0$, there
exists some $\delta(\omega)>0$ such that for every
$||\theta-\theta^{\prime}||_{\infty}<\delta$ and every $J\in\mathbb{N}$ we
have that

$\big{|}|\eta_{n}^{J}(h_{\theta})(\omega)-\eta_{n}^{J}(h_{\theta^{\prime}})(\omega)\big{|}|_{\infty}<\epsilon.$
(130)

Then, by Arzela-Ascoli, there exists a uniformly convergent subsequence. We
claim that there is only one subsequential limit. When the gradient is
bounded, we can treat the gradient as a bounded functional. So by Theorem 6
the sequence $(\eta_{n}^{J}(h_{\theta})(\omega))_{J\in\mathbb{N}}$ converges
pointwise for $\theta=\phi$ and almost every $\omega\in\Omega$, and there is
therefore only one subsequential limit. The sequence therefore converges
uniformly to its limit $\lim_{J\to\infty}\eta_{n}^{J}(h_{\theta})(\omega).$
Therefore, with uniform convergence for the derivatives established, we can
swap the limit and derivative, and obtain that for almost every
$\omega\in\Omega$,

$\lim_{J\to\infty}\eta_{n}^{J}(h_{\theta})(\omega)=\nabla_{\theta}\lim_{J\to\infty}\eta_{n}^{J}(h_{\theta})(\omega).$
(131)

Again from Theorem 6, we know that

$\eta_{n}^{J}(h_{\theta})(\omega)\to\eta_{n}(h_{\theta})=h^{*}_{\theta}$ (132)

for almost every $\omega\in\Omega$. We then have that for almost every
$\omega\in\Omega$,

$\lim_{J\to\infty}\eta_{n}^{J}(h_{\theta})(\omega)=\nabla_{\theta}\eta_{n}(h_{\theta})=\nabla_{\theta}h^{*}_{\theta},$
(133)

as we wanted. ∎

The proof of Theorem 4 is now merely a corollary, where we apply
$\eta_{n}^{J}(h_{\theta})=\hat{\mathcal{L}}^{1}_{J}(\theta)$,
$\eta_{n}(h_{\theta})=\mathcal{L}(\theta)=h_{\theta}^{*}$, and then use the
continuous mapping theorem.

##  Appendix F Bias-Variance Analysis

In this section, we prove various rates on the bias, variance, and MSE of
MOP-$\alpha$. First, we note the following relation between the Dobrushin
contraction coefficient and the alpha-mixing coefficients in our context
below.

######  Lemma 9.

Setting $X$ to be the particle collection at time $m$, and $Y$ at time $n$, we
have that the alpha mixing coefficients,

$\int\big{|}f_{XY}(x,y)-f_{X}(x)\,f_{Y}(y)\big{|}\,dx\,dy<\alpha,$ (134)

are bounded by the Dobrushin coefficient, i.e we have

$\int\big{|}f_{Y|X}(y|x_{1})-f_{Y|X}(y|x_{2})\big{|}\,dy<\alpha$ (135)

for all $x_{1}$, $x_{2}$.

###### Proof.

We rewrite the alpha-mixing assertion as

$\int{\int\big{|}f_{Y|X}(y|x)-f_{Y}(y)\big{|}\,dy}\,f_{X}(x)\,dx.$ (136)

We claim that the Dobrushin coefficient implies

$\int\big{|}f_{Y|X}(y|x_{1})-f_{Y}(y)\big{|}\,dy<\alpha$ (137)

for all $x_{1}$. This is shown as follows:

$\displaystyle\int\big{|}f_{Y|X}(y|x)-f_{Y}(y)\big{|}\,dy$
$\displaystyle=\int\left|\int\big{[}f_{Y|X}(y|x_{1})-f_{Y|X}(y|x)\big{]}f_{X}(x)\,dx\right|dy$
(138)
$\displaystyle<\int\int\big{|}f_{Y|X}(y|x_{1})-f_{Y|X}(y|x)\big{|}\,dy\,f_{X}(x)\,dx$
(139) $\displaystyle<\int\alpha f_{X}(x)\,dx$ (140) $\displaystyle=\alpha$
(141)

We then have the desired result. ∎

###  F.1 Warm Up: MOP-$0$ Variance Bound

Note that we made Assumptions (A2) and (A3) to leverage the results from [42]
on the forgetting of the particle filter. This is required to show error
bounds on the gradient estimates that we provide, namely that the error of the
MOP-$1$ estimator, corresponding to the estimator of [28], is $O(N^{4})$, and
that the variance of MOP-$\alpha$ for any $\alpha<1$ is $O(N)$.

We can decompose the MOP-$\alpha$ estimator as follows. When $\theta=\phi$,

$\hat{\mathcal{L}}(\theta):=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}=\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{\sum_{j=1}^{J}w_{n,j}^{F,\theta}}{\sum_{j=1}^{J}w_{n,j}^{P,\theta}}=\prod_{n=1}^{N}L_{n}^{A,\theta,\alpha}=\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{\sum_{j=1}^{J}w_{n,j}^{F,\theta}}{\sum_{j=1}^{J}(w_{n-1,j}^{F,\theta})^{\alpha}}.$
(142)

Now, to illustrate the proof strategy for the general case in an easier
context, we first analyze the special case when $\alpha=0.$ We now prove the
variance bound for this case, presented below. To our knowledge, this result
for the variance of the gradient estimate of [24] is novel.

######  Theorem 8 (Variance of MOP-$0$ Gradient Estimate).

The variance of the gradient estimate from MOP-$0$, i.e. the algorithm of
[24], is:

$\mathrm{Var}(\nabla_{\theta}\hat{\ell}^{0}(\theta))\lesssim\frac{Np\,G^{\prime}(\theta)^{2}}{J}.$
(143)

###### Proof.

When $\alpha=0$, we have:

$\displaystyle\hat{\mathcal{L}}^{0}(\theta):=\prod_{n=1}^{N}L_{n}^{A,\theta,0}$
$\displaystyle=\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{\sum_{j=1}^{J}w_{n,j}^{F,\theta,0}}{\sum_{j=1}^{J}w_{n,j}^{P,\theta,0}}$
(144)
$\displaystyle=\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{1}{J}\sum_{j=1}^{J}s_{n,j}=\prod_{n=1}^{N}L_{n}^{\phi}\cdot\frac{1}{J}\sum_{j=1}^{J}\frac{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\theta})}{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\phi})}.$
(145)

Similarly to the proof of Lemma 2, we define

$s_{n,j}=\frac{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\theta})}{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{P,\phi})},$
(146)

and from the exact same arguments conclude that its gradient when
$\theta=\phi$ is therefore

$\displaystyle\nabla_{\theta}\hat{\ell}^{0}(\theta)$
$\displaystyle:=\sum_{n=1}^{N}\nabla_{\theta}\log\left(L_{n}^{\phi}\frac{1}{J}\sum_{j=1}^{J}s_{n,j}\right)$
(147)
$\displaystyle=\sum_{n=1}^{N}\frac{\nabla_{\theta}\left(L_{n}^{\phi}\frac{1}{J}\sum_{j=1}^{J}s_{n,j}\right)}{\left(L_{n}^{\phi}\frac{1}{J}\sum_{j=1}^{J}s_{n,j}\right)}$
(148)
$\displaystyle=\sum_{n=1}^{N}\frac{\sum_{j=1}^{J}\nabla_{\theta}s_{n,j}}{\sum_{j=1}^{J}s_{n,j}}$
(149)
$\displaystyle=\sum_{n=1}^{N}\frac{1}{J}\sum_{j=1}^{J}\frac{\nabla_{\theta}f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)}{f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\phi};\phi)}$
(150)
$\displaystyle=\frac{1}{J}\sum_{n=1}^{N}\sum_{j=1}^{J}\nabla_{\theta}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right),$
(151)

where we use the derivative of the logarithm, observe that
$\sum_{j=1}^{J}s_{n,j}=J$ when $\theta=\phi$, and use the derivative of the
logarithm where $\theta=\phi$ again. We do this to establish that the gradient
of the log-likelihood estimate is given by the sum of terms over all $N$ and
$J$:

$\nabla_{\theta}\hat{\ell}^{0}(\theta)=\frac{1}{J}\sum_{n=1}^{N}\sum_{j=1}^{J}\nabla_{\theta}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right).$
(152)

Therefore, for a given $\theta_{i}$ in the parameter vector $\theta$,

$\displaystyle\mathrm{Var}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}^{0}(\theta)\right)$
$\displaystyle=\frac{1}{J^{2}}\mathrm{Var}\left(\sum_{n=1}^{N}\sum_{j=1}^{J}\frac{\partial}{\partial_{\theta_{i}}}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right)\right)$
(153)
$\displaystyle=\frac{1}{J^{2}}\sum_{n=1}^{N}\mathrm{Var}\left(\sum_{j=1}^{J}\frac{\partial}{\partial_{\theta_{i}}}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right)\right)$
(154)
$\displaystyle+2\sum_{m<n}\mathrm{Cov}\left(\frac{1}{J}\sum_{j=1}^{J}\frac{\partial}{\partial_{\theta_{i}}}\log\left(f_{Y_{m}|X_{m}}(y_{m}^{*}|x_{m,j}^{F,\theta};\theta)\right),\frac{1}{J}\sum_{j=1}^{J}\frac{\partial}{\partial_{\theta_{i}}}\log\left(f_{Y_{n}|X_{n}}(y_{n}^{*}|x_{n,j}^{F,\theta};\theta)\right)\right)$
(155)
$\displaystyle=\sum_{n=1}^{N}\mathrm{Var}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)+2\sum_{m<n}\mathrm{Cov}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{m}^{0}(\theta),\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right).$
(156)

Here, we use Assumptions (A2) and (A3) that ensure strong mixing. We know from
Theorem 3 of [42] that when $\textbf{M}_{n,n+k}$ is the $k$-step Markov
operator from timestep $n$ and $\beta_{\text{TV}}(M)=\sup_{x,y\in
E}\|M(x,\cdot)-M(y,\cdot)\|_{\mathrm{TV}}=\sup_{\mu,\nu\in\mathcal{P},\mu\neq\nu}\frac{\|\mu
M-\nu M\|_{\mathrm{TV}}}{\|\mu-\nu\|_{\mathrm{TV}}}$ is the Dobrushin
contraction coefficient of a Markov operator,

$\beta_{\text{TV}}(\textbf{M}_{n,n+k})\leq(1-\epsilon)^{\lfloor
k/(c\log(J))\rfloor},$ (157)

i.e. the mixing time of the particle filter is $O(\log(J))$, where $\epsilon$
and $c$ depend on $\bar{M},\underline{M},\bar{G},\underline{G}$ in (A2) and
(A3).

By Lemma 9, the particle filter itself is strong mixing, with $\alpha$-mixing
coefficients $a(k)\leq(1-\epsilon)^{\lfloor k/(c\log(J))\rfloor}$. Therefore,
functions of particles are strongly mixing as well, with $\alpha$-mixing
coefficients bounded by the original (to see this, observe that the
$\sigma$-algebra of the functionals is contained within the original
$\sigma$-algebra). Therefore, by Davydov’s inequality, noting that
$\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\leq
G^{\prime}(\theta)$ by Assumption (A3), and without loss of generality
labeling $m$ and $n$ such that
$\mathbb{E}\big{[}(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{m}^{0}(\theta))^{4}\big{]}^{1/4}\leq\mathbb{E}\big{[}(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta))^{4}\big{]}^{1/4}$,
we see that

$\displaystyle\mathrm{Cov}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{m}^{0}(\theta),\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)$
$\displaystyle\leq
a(n-m)^{1/2}\,\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{m}^{0}(\theta)\right)^{4}\right]^{1/4}\,\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{4}\right]^{1/4}$
(158) $\displaystyle\leq
a(n-m)^{1/2}\,\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{4}\right]^{1/2}.$
(159)

To bound this, we use the fact that

$\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{4}\right]=\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{4}-\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{2}\right]^{2}\right]+\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{2}\right]^{2}$
(160)

alongside Lemma 2 of [42], which shows that

$\displaystyle\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{4}-\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{2}\right]^{2}\right]$
$\displaystyle\lesssim$ $\displaystyle\frac{G^{\prime}(\theta)^{4}}{J^{2}},$
(161)
$\displaystyle\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)-\mathbb{E}\left[\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right]\right)^{2}\right]$
$\displaystyle\lesssim$ $\displaystyle\frac{G^{\prime}(\theta)^{2}}{J}.$ (162)

It follows that

$\mathbb{E}\left[\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)^{4}\right]^{1/2}\lesssim\sqrt{\frac{G^{\prime}(\theta)^{4}}{J^{2}}+\left(\frac{G^{\prime}(\theta)^{2}}{J}\right)^{2}}=\frac{G^{\prime}(\theta)^{2}}{J},$
(163)

and we conclude that

$\mathrm{Cov}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{m}^{0}(\theta),\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)\leq(1-\epsilon)^{\frac{1}{2}\lfloor\frac{|n-m|}{c\log(J)}\rfloor}\frac{G^{\prime}(\theta)^{2}}{J}.$
(164)

Putting it all together, we see that

$\displaystyle\mathrm{Var}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}^{\alpha}(\theta)\right)$
$\displaystyle=\sum_{n=1}^{N}\mathrm{Var}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)+2\sum_{m<n}\mathrm{Cov}\left(\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{m}^{0}(\theta),\frac{\partial}{\partial_{\theta_{i}}}\hat{\ell}_{n}^{0}(\theta)\right)$
(165)
$\displaystyle\leq\frac{NG^{\prime}(\theta)^{2}}{J}+2\sum_{n=1}^{N-1}(N-n)(1-\epsilon)^{\frac{1}{2}\lfloor\frac{n}{c\log(J)}\rfloor}\frac{G^{\prime}(\theta)^{2}}{J}$
(166)
$\displaystyle\leq\frac{G^{\prime}(\theta)^{2}}{J}\left(N+2\sum_{n=1}^{N-1}(N-n)(1-\epsilon)^{\frac{1}{2}\lfloor\frac{n}{c\log(J)}\rfloor}\right)$
(167) $\displaystyle\lesssim\frac{G^{\prime}(\theta)^{2}N}{J}.$ (168)

It then follows that as $\theta\in\Theta\subseteq\mathbb{R}^{p}$,

$\mathrm{Var}\big{(}\nabla_{\theta}\hat{\ell}^{\alpha}(\theta)\big{)}\lesssim\frac{Np\,G^{\prime}(\theta)^{2}}{J}.$
(169)

∎

#### Remark:

Note that the factor of $N$ that pops up here is due to the use of the
unnormalized gradient. If one divides the gradient estimate by $\sqrt{N}$
before usage, the variance does not depend on the horizon. If one divides by
$N$, the error is $O(1/\sqrt{NJ})$.

The variance here of $\frac{Np\,G^{\prime}(\theta)^{2}}{J}$ can be thought of
as a lower bound for the variance of MOP-$\alpha$. We will later see that the
variance of MOP-$\alpha$ contains this term, and an additional term that
corresponds to the memory of the MOP-$\alpha$ gradient estimate.

###  F.2 MOP-$\alpha$ Variance Bound

We are now in a position to tackle the MOP-$\alpha$ variance bound. Here, we
analyze the case for $\alpha\in(0,1)$. The proof strategy is as follows.
Instead of analyzing the variance of the MOP-$\alpha$ gradient estimate
proper, we will analyze the variance of a modified estimator where we truncate
the weights at $k$ timesteps, so this estimator only looks at most $k$
timesteps back while accumulating the importance weights. We write

$\hat{s}_{k}^{\alpha}(\theta):=\sum_{n=1}^{N}\log\left(\frac{\sum_{j=1}^{J}\prod_{i=n-k}^{n}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}{\sum_{j=1}^{J}\prod_{i=n-k}^{n-1}\left(\frac{g_{i,j}^{A,F,\theta}}{g_{i,j}^{A,F,\phi}}\right)^{\alpha^{(n-i)}}}\right)$
(170)

for the log-likelihood estimate from the $k$-truncated estimator with discount
factor $\alpha$. It is the sum of components

$\hat{s}_{n,k}^{\alpha}(\theta):=\frac{1}{J}\sum_{j=1}^{J}\sum_{i=n-k}^{n}\left(\alpha^{n-i}\log\left(g_{i,j}^{A,F,\theta}\right)-\alpha^{n-i+1}\nabla_{\theta}\log\left(g_{i-1,j}^{A,F,\theta}\right)\right)$
(171)

over timesteps $N$.

This enables us to establish strong mixing for the truncated estimator, in
order to get a bound on the variance that is $O(NJ^{-1})$. We emphasize that
this truncated estimator is only a theoretical construct – we do not actually
use the truncated estimator. The truncated estimator, as we note in the
discussion, possesses similar theoretical guarantees as the MOP-$\alpha$
estimator, but would require the user to specify the number of timesteps to
truncate at instead of a discount factor.

Coupled with a bound ensuring that the $k$-truncated estimator provides an
estimate that is close to that of MOP-$\alpha$ proper, we get a bound on the
variance of MOP-$\alpha$ that comprises of the variance of the $k$-truncated
estimator plus the error, for any $k\leq N$. It then holds that the final
variance bound on MOP-$\alpha$ is given by the minimum over all $k$, and is
never larger than $O(\psi(\alpha)+NJ^{-1})$ for some function $\psi$
increasing in $\alpha$.