sectioningparagraphhyperrefToken pgfplotsTable theorem]Algorithm
theorem]Application to state-space models

#  Particle-MALA and Particle-mGRAD:
Gradient-based MCMC methods for high-dimensional state-space models

Adrien Corenflos  Department of Electrical Engineering and Automation,
Aalto University, Finland
<EMAIL_ADDRESS>Axel Finke  Department of Mathematical Sciences,
Loughborough University, UK
<EMAIL_ADDRESS>
###### Abstract

State-of-the-art methods for Bayesian inference in state-space models are (a)
_conditional sequential Monte Carlo (CSMC)_ algorithms; (b) sophisticated
‘classical’ MCMC algorithms like _MALA ,_ or _mGRAD_ from Titsias and
Papaspiliopoulos,  (2018). The former propose $N$ particles at each time step
to exploit the model’s ‘decorrelation-over-time’ property and thus scale
favourably with the time horizon, $T$, but break down if the dimension of the
latent states, $D$, is large. The latter leverage gradient-/prior-informed
local proposals to scale favourably with $D$ but exhibit sub-optimal
scalability with $T$ due to a lack of model-structure exploitation. We
introduce methods which combine the strengths of both approaches. The first,
_Particle-MALA ,_ spreads $N$ particles locally around the current state using
gradient information, thus extending MALA to $T>1$ time steps and $N>1$
proposals. The second, _Particle-mGRAD_ , additionally incorporates
(conditionally) Gaussian prior dynamics into the proposal, thus extending the
mGRAD algorithm to $T>1$ time steps and $N>1$ proposals. We prove that
Particle-mGRAD interpolates between CSMC and Particle-MALA, resolving the
‘tuning problem’ of choosing between CSMC (superior for highly informative
prior dynamics) and Particle-MALA (superior for weakly informative prior
dynamics). We similarly extend other ‘classical’ MCMC approaches like
_auxiliary MALA_ , _aGRAD_ , and _preconditioned Crank–Nicolson–Langevin
(PCNL)_ to $T>1$ time steps and $N>1$ proposals. In experiments, for both
highly and weakly informative prior dynamics, our methods substantially
improve upon both CSMC and sophisticated ‘classical’ MCMC approaches.

##  1 Introduction

###  1.1 Feynman–Kac models

The aim of this work is to construct efficient _Markov chain Monte Carlo
(MCMC)_ updates for sampling from a continuous _joint smoothing_ distribution
$\pi_{T}(\mathbf{x}_{1:T})$ on $\mathcal{X}^{T}$, where
$\mathcal{X}\coloneqq\mathbb{R}^{D}$ and where for any $t\leq T$, we have
defined the following distributions (termed _filters_):

$\displaystyle\pi_{t}(\mathbf{x}_{1:t})\propto\prod_{s=1}^{t}Q_{s}(\mathbf{x}_{s-1:s}).$
(1)

Here, $Q_{t}(\mathbf{x}_{t-1:t})>0$ is differentiable and can be evaluated
point-wise. Throughout this work, we use the convention that quantities with
‘time’ subscripts $t\leq 0$ or $t>T$ should be ignored, so that, e.g.,
$Q_{1}(x_{0:1})\equiv Q_{1}(\mathbf{x}_{1})$ and
$Q_{T+1}(\mathbf{x}_{T:T+1})\equiv 1$. We will frequently work with some
decomposition

$\displaystyle
Q_{t}(\mathbf{x}_{t-1:t})=M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})G_{t}(\mathbf{x}_{t-1:t}),$
(2)

such that

  * •

$M_{t}(\,\cdot\,|\mathbf{x}_{t-1})$ is a density (w.r.t. a suitable version of
the Lebesgue measure) and also defines a Markov transition kernel called
_mutation kernel_ ;

  * •

$G_{t}(\mathbf{x}_{t-1:t})>0$ is called _potential function_.

We assume that these densities and potential functions are differentiable and
that they (as well as their gradients) can be evaluated point-wise. Motivated
by the following example, we will sometimes refer to
$\smash{M_{1:T}(\mathbf{x}_{1:T})\coloneqq\prod_{t=1}^{T}M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})}$
as the _prior dynamics_ of the _latent states_ $\mathbf{x}_{1:T}$ and
$\smash{G_{1:T}(\mathbf{x}_{1:T})\coloneqq\prod_{t=1}^{T}G_{t}(\mathbf{x}_{t-1:t})}$
as the _likelihood_.

######  Example 1 (state-space model).

One important special case of Feynman–Kac models are _state-space models_. A
state-space model is a bivariate Markov chain
$(\mathbf{x}_{t},\mathbf{y}_{t})_{t\geq 1}$ on $\mathcal{X}\times\mathcal{Y}$,
where $\mathcal{X}\coloneqq\mathbb{R}^{D}$ and
$\mathcal{Y}\coloneqq\mathbb{R}^{D^{\prime}}$, with initial density
$p(\mathbf{x}_{1},\mathbf{y}_{1})=f_{1}(\mathbf{x}_{1})g_{1}(\mathbf{y}_{1}|\mathbf{x}_{1})$
and transition densities
$p(\mathbf{x}_{t},\mathbf{y}_{t}|\mathbf{x}_{t-1})=f_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})g_{t}(\mathbf{y}_{t}|\mathbf{x}_{t})$
(w.r.t. a suitable version of the Lebesgue measure). State-space models assume
that only the measurements $\mathbf{y}_{1:T}$ can be observed whilst the
Markov chain $(\mathbf{x}_{t})_{t\geq 1}$ (often representing the evolution of
the phenomenon of interest) is latent. The joint smoothing distribution then
encodes our knowledge of the latent states $\mathbf{x}_{1:T}$ given the
available data $\mathbf{y}_{1:T}$:

$\displaystyle\pi_{T}(\mathbf{x}_{1:T})\coloneqq
p(\mathbf{x}_{1:T}|\mathbf{y}_{1:T})\propto\prod_{t=1}^{T}f_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})g_{t}(\mathbf{y}_{t}|\mathbf{x}_{t}).$
(3)

One possible way of casting such a state-space model as a Feynman–Kac model
(there are others) is then to take
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=f_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})$
and $G_{t}(\mathbf{x}_{t-1:t})=g_{t}(\mathbf{y}_{t}|\mathbf{x}_{t})$. In this
case, $M_{1:T}(\mathbf{x}_{1:T})=p(\mathbf{x}_{1:T})$,
$G_{1:T}(\mathbf{x}_{1:T})=p(\mathbf{y}_{1:T}|\mathbf{x}_{1:T})$, and
$\pi_{t}(\mathbf{x}_{1:t})=p(\mathbf{x}_{1:t}|\mathbf{y}_{1:t})$, for $t\leq
T$.

###  1.2 Sampling the latent states

Performing inference about the latent states $\mathbf{x}_{1:T}$ requires
calculating expectations of the form
$\operatorname{\mathbb{E}}_{\mathbf{x}_{1:T}\sim\pi_{T}}[\varphi(\mathbf{x}_{1:T})]$,
for some integrable test function $\varphi\colon\mathcal{X}^{T}\to\mathbb{R}$.
Unfortunately, such expectations do not admit closed-form expressions in most
realistic problems and must be approximated by some Monte Carlo estimate
$\smash{\frac{1}{I}\sum_{i=1}^{I}\varphi(\mathbf{x}_{1:T}[i])}$ using samples
$(\mathbf{x}_{1:T}[i])_{i=1}^{I}$ (approximately) distributed according to
$\pi_{T}$. These often come from some MCMC algorithm targeting $\pi_{T}$.

##### ‘Classical’ MCMC methods.

Unfortunately, simple MCMC approaches like the _independent
Metropolis–Hastings (IMH)_ algorithm (Hastings, , 1970) perform poorly if the
problem size: $D\times T$, is large due the difficulty of constructing
efficient _global_ (a.k.a.  _independent_) proposal distributions in high
dimensions. To circumvent this difficulty, MCMC algorithms with _local_ moves
like the _random-walk Metropolis (RWM)_ algorithm (Metropolis et al., , 1953),
propose a new state of the Markov chain near the current state. By decreasing
the proposal scale at a suitable rate with the problem size, the RWM algorithm
can circumvent this curse of dimension (Roberts et al., , 1997). Further
improved performance can be achieved by exploiting

  * •

_gradient information,_ i.e. by including gradients of the log-likelihood or
log-target density into the proposal as in the _Metropolis-adjusted Langevin
algorithm (MALA)_ from Besag,  (1994) and in the _auxiliary MALA (aMALA)_ from
Titsias and Papaspiliopoulos,  (2018); and _additionally_

  * •

_prior information,_ i.e. by explicitly incorporating the prior dependence
structure into the proposal as in the _preconditioned Crank–Nicolson–Langevin
(PCNL)_ and related algorithms (see, e.g., Cotter et al., , 2013, and
references therein) or in the _marginal gradient (mGRAD)_ and _auxiliary
gradient (aGRAD)_ 111Throughout this work, ‘aGRAD’ refers more specifically to
the ‘aGrad-z’ algorithm from Titsias and Papaspiliopoulos,  (2018). algorithms
from Titsias,  (2011); Titsias and Papaspiliopoulos,  (2018).

Figure 1(a) illustrates that ‘classical’ MCMC algorithms can scale favourably
with $D$.

However, ‘classical’ MCMC algorithms are agnostic to the ‘decorrelation-over-
time’ structure of the target distribution $\pi_{T}(\mathbf{x}_{1:T})$, i.e.,
to the fact that, for suitably regular models, the correlation of $x_{s}$ and
$\mathbf{x}_{t}$ under $\pi(\mathbf{x}_{1:T})$ decays with $\lvert t-s\rvert$.
For example, for the simple RWM algorithm and MALA, the _step size_ $\delta>0$
(i.e., proposal variance) would need to decrease at a suitable rate with $T$
($\delta\in\operatorname{\mathrm{O}}((DT)^{-1})$ and
$\delta\in\operatorname{\mathrm{O}}((DT)^{-1/3})$, respectively) even if the
model was completely independent across time steps (Roberts and Rosenthal, ,
2001). Therefore, it stands to reason that the scaling of ‘classical’ MCMC
methods like MALA, PCNL or mGRAD/aGRAD with the time horizon $T$ could be
improved by empowering them to exploit this model structure.

##### CSMC methods.

Another popular $\pi_{T}$-invariant MCMC-kernel,
$P_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode conditional
sequential Monte Carlo (CSMC)}}}$, is induced by running the _CSMC_ algorithm
proposed in the seminal work Andrieu et al.,  (2010); Whiteley,  (2010). Given
the current state $\mathbf{x}_{1:T}\in\mathcal{X}^{T}$ of the Markov chain
(then called the _reference path_) this algorithm generates
$\tilde{\mathbf{x}}_{1:T}\sim
P_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}(\,\cdot\,|\mathbf{x}_{1:T})$ as follows, where we write
$[n]\coloneqq\\{1,\dotsc,n\\}$ and $[n]_{0}\coloneqq[n]\cup\\{0\\}$:

  1. 1.

For $t=1,\dotsc,T$, sample some index $k_{t}$ from a uniform distribution on
$[N]_{0}$; set $\mathbf{x}_{t}^{k_{t}}\coloneqq\mathbf{x}_{t}$ and sample the
remaining _particles_
$\smash{\mathbf{x}_{t}^{-k_{t}}\coloneqq(\mathbf{x}_{t}^{0},\dotsc,\mathbf{x}_{t}^{k_{t}-1},\mathbf{x}_{t}^{k_{t}+1},\dotsc,\mathbf{x}_{t}^{N})}$
conditionally independently such that for $n\neq k_{t}$,

$\displaystyle\mathbf{x}_{t}^{n}\sim
M_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}}),$ (4)

for _ancestor indices_ $a_{t-1}^{n}\in[N]_{0}$ whose rôle is explained later.

  2. 2.

Return
$\smash{\tilde{\mathbf{x}}_{1:T}\coloneqq(\mathbf{x}_{1}^{l_{1}},\dotsc,\mathbf{x}_{t}^{l_{T}})}$,
for indices $l_{1},\dotsc,l_{T}\in[N]_{0}$ sampled from an appropriate
distribution.

Informally, the CSMC algorithm can be interpreted as employing $T$ separate
accept–reject steps (one at each time point) which allows it to exploit the
‘decorrelation-over-time’ property of $\pi_{T}(\mathbf{x}_{1:T})$ (akin to a
‘classical’ MCMC algorithm with blocking in the ‘time’ direction as noted by
Singh et al.,  2017). For suitably regular problems, the CSMC algorithm
therefore scales more favourably with $T$ than ‘classical’ MCMC approaches as
illustrated in Figure 1(b).

(a) Empirical scaling with $D$ for fixed time horizon $T=25$.

(b) Empirical scaling with $T$ for fixed state dimension $D=10$.

Figure 1: Toy linear-Gaussian state-space model with
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{x}_{t-1},\lambda\mathbf{I})$
and
$G_{t}(\mathbf{x}_{t-1:t})=\operatorname{N}(\mathbf{y}_{t};\mathbf{x}_{t},\mathbf{I})$,
where $\mathbf{I}$ is the ($D\times D$)-identity matrix and $\lambda=1$. For a
fair comparison, all methods use $N+1=32$ particles. The step sizes are:
$(TD)^{-1}$ for (multi-proposal) RWM, $D^{-1}$ for Particle-RWM, $(TD)^{-1/3}$
for (multi-proposal) aMALA/MALA/aGRAD, and $D^{-{1/3}}$ for the remaining
(i.e., new) methods. Panel 1(a) illustrates that as $D$ increases, some
‘classical’ MCMC algorithms (RWM, MALA, aMALA and aGRAD) are stable but the
CSMC algorithm breaks down. Conversely, Panel 1(b) illustrates that as $T$
increases, the CSMC algorithm is stable in $T$ but all ‘classical’ MCMC
algorithms (IMH, RWM, MALA, aMALA and aGRAD) break down.

Unfortunately, as shown in Finke and Thiery,  (2023), the CSMC algorithm
suffers from a curse of dimension in the state dimension $D$: as $D$
increases, it becomes increasingly likely that $\tilde{\mathbf{x}}_{1:T}$
coincides exactly with $\mathbf{x}_{1:T}$, i.e., the induced MCMC chain gets
stuck. This is unsurprising because the CSMC algorithm generalises the IMH
algorithm to $T>1$ time steps and $N>1$ proposals. Indeed, note that (4) is
again an independent (i.e. global) proposal in the sense that it does not
depend on the time-$t$ component of the current state of the Markov chain,
$\mathbf{x}_{t}$; and such proposals are known to scale poorly with dimension
(due to the difficulty of finding efficient global proposals in high
dimensions). The only potential remedy: increasing $N$ exponentially with $D$,
is prohibitively costly.

##### Existing combinations of ‘classical’ MCMC and CSMC.

To circumvent this problem, Finke and Thiery,  (2023) introduced the
_Particle-RWM_ 222Referred to as ‘random-walk CSMC’ therein. algorithm which
scatters the particles locally around the reference path (see also
Shestopaloff and Neal, , 2018; Malory, , 2021, for related approaches). That
is, conditional on the reference path $\mathbf{x}_{1:T}$, the remaining
particles $\mathbf{x}_{t}^{-k_{t}}$ are proposed from a joint distribution
under which

$\displaystyle\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{x}_{t},\delta_{t}\mathbf{I}),$
(5)

for $n\neq k_{t}$, where $\mathbf{I}$ is the $(D\times D)$-identity matrix. As
noted in Tjelmeland,  (2004), sampling from this joint proposal distribution
can be achieved by first sampling an auxiliary variable
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$
and then
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$,
for $n\neq k_{t}$. Finke and Thiery,  (2023) also showed that scaling the step
size as $\delta_{t}\in\operatorname{\mathrm{O}}(D^{-1})$ (independently of
$T$) guarantees stability in high dimensions. This is again unsurprising
because the Particle-RWM algorithm generalises the RWM algorithm with Gaussian
proposals (and proposal variance $\delta_{1}$) to $T>1$ time steps and $N>1$
proposals. Recently, Corenflos and Särkkä,  (2023) showed that the Particle-
RWM algorithm can be viewed as a Gibbs-sampling step for the auxiliary
variables $\mathbf{u}_{t}$ followed by a CSMC update which targets a modified
Feynman–Kac model which depends on $\mathbf{u}_{1:T}$, allowing for greater
flexibility in the choice proposals. Including related ‘pseudo observations’
$\mathbf{u}_{t}$ into CSMC updates had previously been suggested by Murray et
al.,  (2013); Fearnhead and Meligkotsidou,  (2016); Karppinen and Vihola,
(2021) but primarily aimed at overcoming the problem that the CSMC algorithm
mixes poorly if the initial distribution $M_{1}(\mathbf{x}_{1})$ is diffuse
(and potentially also for improving mixing in the presence of ‘static’ model
parameters).

###  1.3 Contributions

Recall that in the ‘classical’ MCMC setting, improved performance can often be
achieved by enhancing the proposal distribution using gradient or prior
information. Thus, in this work, we introduce a methodology which combines the
strength of CSMC methods (i.e., exploitation of the ‘decorrelation-over-time’
property of the target distribution) with the strengths of sophisticated
‘classical’ MCMC approaches (i.e., gradient-enhanced local proposals).

In the remainder of this section, we detail the contributions of this paper
(Table 1 summarises our proposed methodology).

In Section 3, we introduce the following CSMC type methods which propose
particles locally around the reference path guided by gradient information:

  * •

Particle-aMALA. In Section 3.1, we extend the Particle-RWM algorithm to
incorporate gradient information into the proposals. That is, conditional on
the reference path $\mathbf{x}_{1:T}$, the remaining particles
$\mathbf{x}_{t}^{-k_{t}}$ are proposed from a joint distribution under which

$\displaystyle\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{x}_{t}+\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:T}),\delta_{t}\mathbf{I}),$
(6)

for $n\neq k_{t}$. Sampling from this joint proposal can be achieved by first
sampling an auxiliary variable
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t}),\tfrac{\delta_{t}}{2}\mathbf{I})}$
and then
$\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})$,
for $n\neq k_{t}$. We call this method _Particle-aMALA_ because the auxiliary
variables $\mathbf{u}_{t}$ are explicitly included in the space, i.e. they
appear in the particle weights, and because the algorithm generalises a
version of _auxiliary MALA (aMALA)_ from Titsias and Papaspiliopoulos,  (2018)
to $T>1$ time steps and $N>1$ proposals.

  * •

Particle-MALA. In Section 3.2, we improve Particle-aMALA by marginalising out
the auxiliary variables $\mathbf{u}_{t}$. We call the resulting method
_Particle-MALA_ because it generalises _MALA_ (Besag, , 1994) to $T>1$ time
steps and $N>1$ proposals.

  * •

Particle-aMALA+. In Section 3.3, we improve Particle-aMALA differently by
replacing the ‘filter’ gradient
$\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t})$ in (6) with the
‘smoothing’ gradient $\nabla_{\mathbf{x}_{t}}\log\pi_{T}(\mathbf{x}_{1:T})$
which is beneficial when future observations are informative about the current
state. We call the resulting method _Particle-aMALA+_.

In Section 4, we consider the special case that the Feynman–Kac model has
conditionally Gaussian mutation kernels:
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$.
In this setting, we introduce the following methods which propose particles
locally around the reference path guided by both gradient information and
prior information:

  * •

Particle-aGRAD. In Section 4.1, we propose an algorithm which, conditional on
the reference path $\mathbf{x}_{1:T}$, proposes the remaining particles
$\mathbf{x}_{t}^{-k_{t}}$ from a joint distribution under which

$\displaystyle\mathbf{x}_{t}^{n}$
$\displaystyle\sim\operatorname{N}\bigl{(}(\mathbf{I}-\mathbf{A}_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{n}}))\mathbf{m}_{t}(\smash{\mathbf{x}_{t-1}^{a_{t-1}^{n}}})+\mathbf{A}_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{n}})[\mathbf{x}_{t}+\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t})],\mathbf{B}_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{n}})\bigr{)},$
(7)

for $n\neq k_{t}$, where
$\mathbf{A}_{t}(\mathbf{x})\coloneqq(\mathbf{C}_{t}(\mathbf{x})+\tfrac{\delta_{t}}{2}\mathbf{I})^{-1}\mathbf{C}_{t}(\mathbf{x})$
and
$\mathbf{B}_{t}(\mathbf{x})\coloneqq\tfrac{\delta_{t}}{2}\mathbf{A}_{t}(\mathbf{x})^{2}+\mathbf{A}_{t}(\mathbf{x})$.
Sampling from this joint proposal can be achieved by first sampling an
auxiliary variable
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\tfrac{\delta_{t}}{2}\mathbf{I})}$ and then
$\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{t})}$,
for $n\neq k_{t}$, where
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})=p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})$
is the fully-adapted auxiliary particle-filter proposal for the state-space
model with Gaussian transitions $\mathbf{x}_{t}|\mathbf{x}_{t-1}\sim
M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})$ and pseudo observations
$\mathbf{u}_{t}|\mathbf{x}_{t}\sim\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t};\frac{\delta_{t}}{2}\mathbf{I})$.
We call this the _Particle-aGRAD_ algorithm because the auxiliary variables
$\mathbf{u}_{t}$ again appear in the particle weights, and because it
generalises the powerful _aGRAD_ algorithm from Titsias and Papaspiliopoulos,
(2018) to $T>1$ time steps and $N>1$ proposals.

  * •

Particle-mGRAD. In Section 4.2, under the assumption that
$\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t}$ and in analogy to Section
3.2, we improve Particle-aGRAD by marginalising out the auxiliary variables
$\mathbf{u}_{t}$. We call the resulting method _Particle-mGRAD_ because it
generalises the powerful _mGRAD_ algorithm from Titsias and Papaspiliopoulos,
(2018) to $T>1$ time steps and $N>1$ proposals.

  * •

Particle-aGRAD+. In Section 4.3, in analogy to Section 3.3, we improve
Particle-aGRAD by replacing the ‘filter-potential’ gradients
$\nabla_{\mathbf{x}_{t}}\log G_{t}(\mathbf{x}_{t-1:t})$ in (7) with
‘smoothing-potential’ gradients $\nabla_{\mathbf{x}_{t}}\log
G_{1:T}(\mathbf{x}_{1:T})$ which may be beneficial if
$G_{t}(\mathbf{x}_{t-1:t})$ varies significantly in $\mathbf{x}_{t-1}$. We
call this method _Particle-aGRAD+_.

  * •

Twisted Particle-aGRAD (+). In Section 4.4, under the assumption that
$\mathbf{m}_{t}(\mathbf{x}_{t-1})=\mathbf{F}_{t}\mathbf{x}_{t-1}+\mathbf{b}_{t}$
and $\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t}$, we improve Particle-
aGRAD and Particle-aGRAD+ by instead using all future auxiliary variables
$\mathbf{u}_{t:T}$ to propose $\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{t:T})}$,
for $n\neq k_{t}$, where
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})=p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t:T})$
is the fully _twisted_ particle filter proposal for the state-space model with
Gaussian transitions and pseudo observations $\mathbf{u}_{t}$ mentioned above.
We call the resulting methods _twisted Particle-aGRAD_ and _twisted Particle-
aGRAD+_.

In Section 4.6, we prove that Particle-aGRAD and Particle-mGRAD (and their
smoothing-gradient/twisted variants) solve the ‘tuning’ problem of having to
choose between:

  1. 1.

the CSMC algorithm (which proposes particles solely based on the prior
dynamics);

  2. 2.

the Particle-aMALA, Particle-MALA or Particle-aMALA+ (which propose particles
solely locally around the reference path).

Figure 2: Empirical illustration of the ‘interpolation’ from Propositions 8
and 9 in the toy linear-Gaussian state-space model from Figure 1 (with
$D=T=10$).

This choice is not always clear: on the one hand, Choice 2 can exhibit
superior performance in high dimensions. On the other hand, if the prior
dynamics are highly informative then Choice 1 can outperform Choice 2.
Specifically, we prove that the following results hold in stationarity and
under the simplifying assumption that the model factorises over time, i.e., if
$G_{t}$, $\mathbf{m}_{t}$, $\mathbf{C}_{t}$ (and hence $\mathbf{A}_{t}$ and
$\mathbf{B}_{t}$ in (7)) do not depend on the state at time $t-1$:

  * •

Proposition 8. Particle-aGRAD and Particle-mGRAD reduce to the CSMC algorithm
as prior dynamics become _more_ informative. Informally, we then have
$\mathbf{A}_{t}\approx\bm{0}$ and $\mathbf{B}_{t}\approx\mathbf{C}_{t}$ so
that (7) reduces to (4).

  * •

Proposition 9. Particle-aGRAD and Particle-mGRAD reduce to Particle-aMALA and
Particle-MALA, respectively, as prior dynamics become _less_ informative.
Informally, we then have $\mathbf{A}_{t}\approx\mathbf{I}$ and
$\mathbf{B}_{t}\approx\delta_{t}\mathbf{I}$ so that (7) reduces to (6).

Propositions 8 and 9 are illustrated in Figure 2 for a model in which the
independence across time-steps is not verified. As a by-product, these
propositions show that the aGRAD/mGRAD algorithms from Titsias and
Papaspiliopoulos,  (2018) can be viewed as automatically interpolating between
the IMH algorithm and aMALA/MALA, depending on the ‘informativeness’ of the
prior. To our knowledge, this has not been pointed out in the literature. As
another by-product, the methodology presented in this section also addresses
the ‘tuning problem’ of having to choose whether to sample the initial latent
state $\mathbf{x}_{1}$ within the CSMC scheme (which is preferable if the
prior on the initial state is informative) or to treat it as a ‘static’
parameter to be sampled separately (which is preferable if this prior is
diffuse, see Murray et al., , 2013; Fearnhead and Meligkotsidou, , 2016;
Karppinen and Vihola, , 2021).

In Section 5, we demonstrate the performance of our methodology on a high-
dimensional multivariate stochastic volatility model, often used as a
benchmark in the particle filtering literature. The different methods proposed
in this article dramatically improve on existing CSMC and related methods and
also on ‘classical’ MCMC methods in terms of effective sample size for
different levels of prior informativeness.

All proofs (e.g., of the fact that the proposed methods leave
$\pi_{T}(\mathbf{x}_{1:T})$ invariant) are deferred to the appendix.
Additionally, in Appendix A, we introduce _Particle-PCNL_ methods which
generalise the _preconditioned Crank–Nicolson–Langevin (PCNL)_ algorithm from
Cotter et al.,  (2013) to $T>1$ time steps and $N>1$ proposals. The methods
proposed in this work and their special cases if $N=T=1$ are summarised in
Table 1. Note that for $T=1$ but $N>1$, our work implies novel _multi-
proposal_ versions of ‘classical’ MCMC kernels like MALA, aMALA, mGRAD, aGRAD
and PCNL. These may be of independent interest because they can exploit
parallel computing architectures for inference in non-dynamic models.

Importantly, and in keeping with existing CSMC methodology, the computational
cost of all our proposed algorithms is linear in both $T$ and $N$, in time and
memory.

Finally, the Python code for reproducing our experiments is publicly available
at https://github.com/AdrienCorenflos/particle_mala. It was written as a
library and can be extended to accommodate other models than the ones
considered here.

Table 1: The methods mentioned in this work (new methods are in italic).

Method |  Section |  Special case if $N=T=1$
---|---|---
CSMC † | 2.1 | IMH
Particle-RWM | 2.2 | RWM
Particle-aMALA | 3.1 | aMALA
Particle-MALA | 3.2 | MALA
Particle-aMALA+ | 3.3 | aMALA
Particle-aGRAD | 4.1 | aGRAD
Particle-mGRAD | 4.2 | mGRAD
Particle-aGRAD+ | 4.3 | aGRAD
Twisted Particle-aGRAD (+) | 4.4 | aGRAD
Particle-PCNL & more‡ |  Appendix A | PCNL

  * †

In our taxonomy, CSMC could be called ‘Particle-IMH’. However, the latter
already refers to an entirely different algorithm in Andrieu et al.,  (2010).

  * ‡

We again also describe auxiliary-variable, smoothing-gradient (‘+’) and
twisted versions.

##  2 Existing methodology

###  2.1 CSMC (particle extension of IMH)

####  2.1.1 Algorithm

Assume that we can generate independent and identically distributed (IID)
samples from the mutation kernels $M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})$. A
method for constructing a $\pi_{T}$-invariant MCMC kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ is then given by the
_CSMC_ algorithm from Andrieu et al.,  (2010) which proposes $N$ particles at
each time step to build up an efficient proposal. Algorithm 2.1.1 summarises
the scheme, where ‘w.p.’ is short for ‘with probabilitity’. We also
recursively define the $n$th surviving particle lineage at time $t$ as

$\displaystyle\mathbf{x}_{1:t}^{(n)}\coloneqq(\mathbf{x}_{1:t-1}^{(a_{t-1}^{n})},\mathbf{x}_{t}^{n}).$
(8)

In particular, therefore,
$\smash{\mathbf{x}_{t-1:t}^{(n)}=(\mathbf{x}_{t-1}^{a_{t-1}^{n}},\mathbf{x}_{t}^{n})}$.
[CSMC ] Given $\mathbf{x}_{1:T}\in\mathcal{X}^{T}$: 1. for $t=1,\dotsc,T$, (a)
sample $k_{t}$ from a uniform distribution on $[N]_{0}$ and set
$\mathbf{x}_{t}^{k_{t}}\coloneqq\mathbf{x}_{t}$, (b) if $t>1$, set
$\smash{a_{t-1}^{k_{t}}\coloneqq k_{t-1}}$ and sample $a_{t-1}^{n}=i$ w.p.
$\smash{W_{t-1}^{i}}$, for $n\in[N]_{0}\setminus\\{k_{t}\\}$, (c) sample
$\smash{\mathbf{x}_{t}^{n}\sim
M_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}})}$ for
$n\in[N]_{0}\setminus\\{k_{t}\\}$, (d) for $n\in[N]_{0}$, set
$\smash{w_{t}^{n}\propto G_{t}(\mathbf{x}_{t-1:t}^{(n)})}$. (e) for
$n\in[N]_{0}$, set $\smash{W_{t}^{n}\coloneqq
w_{t}^{n}/\sum_{m=0}^{N}w_{t}^{m}}$; 2. sample
$i\in[N]_{0}\setminus\\{k_{T}\\}$ w.p. $\dfrac{W_{T}^{i}}{1-W_{T}^{k_{T}}}$;
set $l_{T}\coloneqq i$ w.p.
$\smash{1\wedge\dfrac{1-W_{T}^{k_{T}}}{1-W_{T}^{i}}}$; otherwise, set
$l_{T}\coloneqq k_{T}$; 3. for $t=T-1,\dotsc,1$, sample $l_{t}=i\in[N]_{0}$
w.p.
$\dfrac{W_{t}^{i}Q_{t+1}(\mathbf{x}_{t}^{i},\mathbf{x}_{t+1}^{l_{t+1}})}{\sum_{n=0}^{N}W_{t}^{n}Q_{t+1}(\mathbf{x}_{t}^{n},\mathbf{x}_{t+1}^{l_{t+1}})};$
4. return
$\tilde{\mathbf{x}}_{1:T}\coloneqq(\mathbf{x}_{1}^{l_{1}},\dotsc,\mathbf{x}_{t}^{l_{T}})$.
Algorithm 2.1.1 includes two extensions to the original presentation of the
CSMC algorithm in Andrieu et al.,  (2010):

  * •

Step 2 uses the so-called _forced-move_ extension for CSMC algorithms which
was proposed in Chopin and Singh,  (2013) (see also Liu, , 1996). The
algorithm would still be valid if we instead sampled $l_{T}=i\in[N]_{0}$ with
probability $W_{T}^{i}$.

  * •

Step 3 is the _backward-sampling_ extension from Whiteley,  (2010). The
algorithm would still be valid if we instead set
$\smash{l_{t}=a_{t}^{l_{t+1}}}$ (but typically much less efficient, especially
if $T$ is large).

Importantly, sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ as
described in Algorithm 2.1.1 induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves $\pi_{T}$
invariant. For sufficiently ergodic models, this MCMC kernel can yield highly
efficient updates of the sequence of latent states, even if the time horizon
$T$ is large (Lee et al., , 2020; Karjalainen et al., , 2023).

####  2.1.2 Relationship with ‘classical’ MCMC algorithms

Interestingly, the CSMC algorithm generalises the classical _IMH_ algorithm
(Hastings, , 1970) in the sense that the former reduces to the latter if
$T=N=1$. This can be seen as follows, where we suppress the ‘time’ subscript
$t=1$ everywhere to simplify the notation. Given that the current state of the
Markov chain is $\mathbf{x}=\mathbf{x}^{0}$ (we can assume that $k=0$ without
loss of generality), Step 1c of Algorithm 2.1.1 proposes $\mathbf{x}^{1}\sim
M$. The remaining steps return $\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{1}$ as
the new state with acceptance probability
$1\wedge\alpha_{\text{\lx@glossaries@gls@link{acronym}{IMH}{\leavevmode
IMH}}}(\mathbf{x}^{0},\mathbf{x}^{1})$, where

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{IMH}{\leavevmode
IMH}}}(\mathbf{x}^{0},\mathbf{x}^{1})$
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}=\frac{G(\mathbf{x}^{1})}{G(\mathbf{x}^{0})}=\frac{\pi(\mathbf{x}^{1})M(\mathbf{x}^{0})}{\pi(\mathbf{x}^{0})M(\mathbf{x}^{1})}.$
(9)

Otherwise, the old state
$\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{0}=\mathbf{x}$ is returned as the new
state.

####  2.1.3 Breakdown in high dimensions

Unfortunately, as shown in Finke and Thiery,  (2023), Algorithm 2.1.1 suffers
from a curse of dimension if $D$ is large (unless the number of proposed
particles, $N$, grows exponentially in $D$ but that is prohibitive). This is
not surprising since the IMH algorithm is known to break down in high
dimensions (due to the difficulty of finding an efficient global proposal
distribution $M$ in high dimensions).

###  2.2 Particle-RWM

####  2.2.1 Algorithm

To circumvent the curse of dimension, Finke and Thiery,  (2023) (see also
Shestopaloff and Neal, , 2018; Malory, , 2021, for related methods) developed
the _particle random-walk Metropolis (Particle-RWM)_ algorithm which scatters
the proposed particles locally around the reference path using Gaussian
perturbations as outlined in Algorithm 2.2.1.

[Particle-RWM ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c) and the weight calculation (Step 1d) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$,
and
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$,
for $n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

for $n\in[N]_{0}$, set $\smash{w_{t}^{n}\propto
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})}$.

Notably, Step 1c. marginally samples
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{x}_{t},\delta_{t}\mathbf{I})}$,
for $n\neq k_{t}$.

####  2.2.2 Interpretation as a CSMC update on an extended space

Corenflos and Särkkä,  (2023) showed that Algorithm 2.2.1 can be derived by
including the auxiliary variables $\mathbf{u}_{t}$ into the space and thus
considering the extended distribution

$\displaystyle\pi^{\prime}_{T}(\mathbf{x}_{1:T},\mathbf{u}_{1:T})$
$\displaystyle\coloneqq\pi_{T}(\mathbf{x}_{1:T})\prod_{t=1}^{T}\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\tfrac{\delta_{t}}{2}\mathbf{I}),$
(10)

which admits $\pi_{T}(\mathbf{x}_{1:T})$ as a marginal and which can be
targeted by alternating the following two steps. Given
$\mathbf{x}_{1:T}\in\mathcal{X}^{T}$,

  1. 1.

sample
$\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})$,
for $t=1,\dotsc,T$;

  2. 2.

run the CSMC algorithm (Algorithm 2.1.1) but with
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})$, $G_{t}(\mathbf{x}_{t-1:t})$, and
$Q_{t}(\mathbf{x}_{t-1:t})$ replaced by
$\smash{M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})\coloneqq\operatorname{N}(\mathbf{x}_{t};\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$,
$\smash{G^{\prime}_{t}(\mathbf{x}_{t-1:t})\coloneqq
Q_{t}(\mathbf{x}_{t-1:t})}$ and
$\smash{Q^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})}\coloneqq\smash{M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})G^{\prime}_{t}(\mathbf{x}_{t-1:t})}$.

In particular, this shows that sampling $\tilde{\mathbf{x}}_{1:T}$ given
$\mathbf{x}_{1:T}$ via Algorithm 2.2.1 induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLERWM}{\leavevmode Particle-
RWM}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves $\pi_{T}$
invariant.

####  2.2.3 Relationship with ‘classical’ MCMC algorithms

The Particle-RWM algorithm generalises the classical (Gaussian) _RWM_
algorithm of Metropolis et al.,  (1953) in the sense that the former reduces
to the latter if $T=N=1$. This can be seen as follows, where we again suppress
the ‘time’ subscript $t=1$ everywhere to simplify the notation. Given that the
current state of the Markov chain is $\mathbf{x}=\mathbf{x}^{0}$ (we can again
assume that $k=0$ without loss of generality), Step 1c. of Algorithm 2.2.1
proposes
$\mathbf{x}^{1}\sim\operatorname{N}(\mathbf{x}^{0},\delta\mathbf{I})$. The
remaining steps return $\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{1}$ as the new
state with acceptance probability
$1\wedge\alpha_{\text{\lx@glossaries@gls@link{acronym}{RWM}{\leavevmode
RWM}}}(\mathbf{x}^{0},\mathbf{x}^{1})$, where

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{RWM}{\leavevmode
RWM}}}(\mathbf{x}^{0},\mathbf{x}^{1})$
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}=\frac{\pi(\mathbf{x}^{1})}{\pi(\mathbf{x}^{0})}.$
(11)

Otherwise, the old state
$\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{0}=\mathbf{x}$ is returned as the new
state.

####  2.2.4 Stability in high dimensions

Finke and Thiery,  (2023) proved that the Particle-RWM algorithm circumvents
the curse of dimensionality if the proposal variance is scaled as
$\delta_{t}\in\operatorname{\mathrm{O}}(D^{-1})$ (see also Malory, , 2021, for
a proof for non-Gaussian exchangeable proposals but in the case where the
model factorises over time). However, from the literature on classical MCMC
algorithms, it is well known that faster convergence rates can be achieved by
incorporating gradient information into the proposal (Roberts and Rosenthal, ,
1998). Thus, in the next section, we extend the Particle-RWM to allow for
gradient-informed proposals.

##  3 Particle extensions of MALA and aMALA

###  3.1 Particle-aMALA

We now propose _Particle-aMALA_ , a method which extends the Particle-RWM
algorithm from Finke and Thiery,  (2023) by allowing for the use of gradient
information in the proposal. For the moment, gradients are taken w.r.t. the
filtering densities and we employ an indicator $\kappa\in\\{0,1\\}$ to permit
switching off the use of gradient information.

We now write

$\displaystyle M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})$
$\displaystyle\coloneqq\operatorname{N}(\mathbf{x}_{t};\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I}),$
(12) $\displaystyle G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$
$\displaystyle\coloneqq
Q_{t}(\mathbf{x}_{t-1:t})\frac{\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t}),\tfrac{\delta_{t}}{2}\mathbf{I})}{\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})},$
(13)

as well as $Q^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$,
where we note that

$\displaystyle\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t})=\nabla_{\mathbf{x}_{t}}\log
Q_{t}(\mathbf{x}_{t-1:t}).$ (14)

A single iteration of the Particle-aMALA is then as follows.

[Particle-aMALA ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c) and the weight calculation (Step 1d) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t}),\tfrac{\delta_{t}}{2}\mathbf{I})}$,
and
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$,
for $n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

for $n\in[N]_{0}$, set $\smash{w_{t}^{n}\propto
G^{\prime}_{t}(\mathbf{x}_{t-1:t}^{(n)};\mathbf{u}_{t}),}$

and also replace $Q_{t+1}(\,\cdot\,)$ in the backward kernel in Step 3 by
$\smash{Q^{\prime}_{t+1}(\,\cdot\,;\mathbf{u}_{t+1})}$.

Step 1c. marginally samples
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t}),\delta_{t}\mathbf{I})}$,
for $n\neq k_{t}$. This follows from Lemma 3 in Appendix C.

######  Proposition 1 (validity of Particle-aMALA).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm 3.1
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-aMALA}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves
$\pi_{T}$ invariant.

###  3.2 Particle-MALA

In this section, we analytically integrate out the auxiliary variables
$\mathbf{u}_{t}$ appearing in the weights of the Particle-aMALA. A single
iteration of the resulting methodology – which we term the _Particle-MALA_ –
is as follows, where we write

$\displaystyle\log
H_{t,\bm{\phi}}(\mathbf{x},\bar{\mathbf{x}})\coloneqq\frac{1}{\delta_{t}}\bigl{[}2\bm{\phi}^{\mathrm{T}}(\bar{\mathbf{x}}-\mathbf{x})-\tfrac{N}{N+1}\bm{\phi}^{\mathrm{T}}\bm{\phi}\bigr{]}.$
(15)

[Particle-MALA ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c) and the weight calculation (Step 1d) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t}),\tfrac{\delta_{t}}{2}\mathbf{I})}$,
and
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$,
for $n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

set
$\smash{\bar{\mathbf{x}}_{t}\coloneqq\tfrac{1}{N+1}\sum_{n=0}^{N}\mathbf{x}_{t}^{n}}$
and, for $n\in[N]_{0}$,

$\displaystyle w_{t}^{n}\propto
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})H_{t,\kappa\frac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}^{n}}\log
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})}(\mathbf{x}_{t}^{n},\bar{\mathbf{x}}_{t}).$
(16)

Step 1d. pre-computes $\smash{\bar{\mathbf{x}}_{t}}$ to ensure that the
algorithm can still be implemented in $\operatorname{\mathrm{O}}(N)$
operations even though the weight of the $n$th particle now depends on the
values of all $N+1$ particles. However, note that the auxiliary variables
$\mathbf{u}_{t}$ no longer appear in the weights.

######  Remark 1 (Particle-aMALA ‘exactly approximates’ Particle-MALA).

Note that the Particle-aMALA differs from the Particle-MALA only in the
definition of the weights (and the backward-sampling weights). This allows us
to interpret the former as a ‘noisy’ version of the latter. Indeed, write the
unnormalised weight of the $n$th particle at time-$t$ in the Particle-aMALA as
$w_{t}^{n}(\mathbf{u}_{t})$, whilst $w_{t}^{n}$ denotes the corresponding
weight under the Particle-MALA (which does not depend on the auxiliary
variable $\mathbf{u}_{t}$). Then we have

$\displaystyle\frac{w_{t}^{n}(\mathbf{u}_{t})}{w_{t}^{k_{t}}(\mathbf{u}_{t})}$
$\displaystyle=\frac{w_{t}^{n}}{w_{t}^{k_{\mathrlap{t}}}}\times\frac{q_{t}^{-n}(\mathbf{u}_{t}|\mathbf{x}_{t}^{-n},\mathbf{x}_{t}^{n};\mathcal{H}_{t-1})}{q_{t}^{-k_{t}}(\mathbf{u}_{t}|\mathbf{x}_{t}^{-k_{t}},\mathbf{x}_{t}^{k_{t}};\mathcal{H}_{t-1})},$
(17)

where
$\smash{q_{t}^{-n}(\mathbf{u}_{t}|\mathbf{x}_{t}^{-n},\mathbf{x}_{t}^{n};\mathcal{H}_{t-1})=\operatorname{N}(\mathbf{u}_{t};\bar{\mathbf{x}}_{t}+\kappa\tfrac{\delta_{t}}{2(N+1)}\nabla_{\mathbf{x}_{t}^{n}}\log
Q_{t}(\mathbf{x}_{t-1:t}^{(n)}),\tfrac{\delta_{t}}{2(N+1)}\mathbf{I})}$ is the
conditional distribution of $\mathbf{u}_{t}$ under the joint distribution of
all random variables generated by Algorithm 3.1 up to (and including) time $t$
assuming the reference particle at time $t$ is placed in position $n$ (and
$\mathcal{H}_{t-1}$ denotes the history of the particle system, i.e. all
particles and ancestor indices up to time $t-1$). This conditional
distribution follows from Lemma 3 in Appendix C. In particular, we therefore
have

$\displaystyle\operatorname{\mathbb{E}}\biggl{[}\frac{w_{t}^{n}(\mathbf{u}_{t})}{w_{t}^{k_{t}}(\mathbf{u}_{t})}\biggr{]}=\frac{w_{t}^{n}}{w_{t}^{k_{\mathrlap{t}}}},$
(18)

where the expectation is taken w.r.t.
$\smash{q_{t}^{-k_{t}}(\mathbf{u}_{t}|\mathbf{x}_{t}^{-k_{t}},\mathbf{x}_{t}^{k_{t}};\mathcal{H}_{t-1})}$.
Interestingly, for the Particle-RWM algorithm (recovered by setting
$\kappa=0$), the ‘auxiliary’ and ‘marginal’ variants are statistically
equivalent.

######  Proposition 2 (validity of Particle-MALA).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm 3.2
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode Particle-
MALA}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves $\pi_{T}$
invariant.

###  3.3 Particle-aMALA+

In this section, we extend the Particle-aMALA in a different manner: we now
modify the algorithm so that the proposal distributions incorporate gradients
w.r.t. the joint smoothing distribution $\pi_{T}$ rather than w.r.t. the
filters, $\pi_{t}$. This can be beneficial if there is a significant
discrepancy between the marginal distribution of $\mathbf{x}_{t}$ under the
former and the latter as is typically the case if $D$ is large. Indeed, this
discrepancy is likely the reason for the decay in performance of Particle-
aMALA and Particle-MALA for very large $D$ visible in Figure 1(a).

For $M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})$ and
$G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$ still defined as in the
Particle-aMALA algorithm (i.e., as in (12) and (13)), we now write

$\displaystyle G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})$
$\displaystyle\coloneqq
G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})\frac{\operatorname{N}(\mathbf{u}_{t-1};\mathbf{x}_{t-1}+\kappa\tfrac{\delta_{t-1}}{2}\nabla_{\mathbf{x}_{t-1}}\log\pi_{T}(\mathbf{x}_{1:T}),\tfrac{\delta_{t-1}}{2}\mathbf{I})}{\operatorname{N}(\mathbf{u}_{t-1};\mathbf{x}_{t-1}+\kappa\tfrac{\delta_{t-1}}{2}\nabla_{\mathbf{x}_{t-1}}\log\pi_{t-1}(\mathbf{x}_{1:t-1}),\tfrac{\delta_{t-1}}{2}\mathbf{I})},\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$

as well as $Q^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})$,
where we note that

$\displaystyle\nabla_{\mathbf{x}_{t}}\log\pi_{T}(\mathbf{x}_{1:T})=\nabla_{\mathbf{x}_{t}}[\log
Q_{t}(\mathbf{x}_{t-1:t})+\log Q_{t+1}(\mathbf{x}_{t:t+1})].$ (19)

A single iteration of the resulting ‘smoothing-gradient’ methodology – which
we term the _Particle-aMALA+_ – is then as follows.

[Particle-aMALA+ ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c), the weight calculation (Step 1d), and backward sampling (Step 3) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{T}(\mathbf{x}_{1:T}),\tfrac{\delta_{t}}{2}\mathbf{I})}$,
and
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})}$,
for $n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

for $n\in[N]_{0}$, set $\smash{w_{t}^{n}\propto
G^{\prime}_{t}(\mathbf{x}_{t-2:t}^{(n)};\mathbf{u}_{t-1:t}),}$

  3. 3.

for $t=T-1,\dotsc,1$, sample $l_{t}=i\in[N]_{0}$ w.p.

$\displaystyle\dfrac{W_{t}^{i}Q^{\prime}_{t+1}((\mathbf{x}_{t-1:t}^{(i)},\mathbf{x}_{t+1}^{l_{t+1}});\mathbf{u}_{t:t+1})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{i},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{t+1:t+2})}{\sum_{n=0}^{N}W_{t}^{n}Q^{\prime}_{t+1}((\mathbf{x}_{t-1:t}^{(n)},\mathbf{x}_{t+1}^{l_{t+1}});\mathbf{u}_{t:t+1})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{n},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{t+1:t+2})}.$
(20)

In Step 3., we recall the convention that any quantity with ‘time’ index $t>T$
should be ignored, so that $Q^{\prime}_{T+1}\equiv 1$. Some comments about
Algorithm 3.3 are in order.

  * •

Step 1c. marginally samples
$\smash{\mathbf{x}_{t}^{n}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{T}(\mathbf{x}_{1:T}),\delta_{t}\mathbf{I})}$,
for $n\neq k_{t}$. This is in contrast to the Particle-aMALA and Particle-
MALA, whose (marginal) proposal distribution is centred around
$\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log\pi_{t}(\mathbf{x}_{1:t})$.

  * •

Steps 1d. and 3. are similar to the weight-calculation and backward-sampling
steps in the previous algorithms. The only difference here is that the model
is now no longer (first-order) Markov in the sense that the (incremental)
weights at time $t$ now also depend on the state at time $t-2$.

######  Proposition 3 (validity of Particle-aMALA+).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm 3.3
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode Particle-
MALA}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves $\pi_{T}$
invariant.

###  3.4 Relationship with other methods

We end this section by relating the proposed algorithms to existing
methodologies.

  1. 1.

Generalisation of Particle-RWM and RWM. If $\kappa=0$, then the algorithms
introduced in this section (Particle-aMALA, Particle-MALA and Particle-aMALA+)
do not make use of any gradient information and reduce to the Particle-RWM
algorithm. In particular, if $T=N=1$, they thus reduce to the RWM algorithm.

  2. 2.

Generalisation of aMALA. For $\kappa=1$, the Particle-aMALA (and similarly the
Particle-aMALA+) algorithm generalise the _auxiliary MALA (aMALA)_ from
Titsias and Papaspiliopoulos,  (2018) in the sense that the former reduces to
the latter if $T=N=1$. This can be seen as follows, where we again suppress
the ‘time’ subscript $t=1$ everywhere. Given that the current state of the
Markov chain is $\mathbf{x}=\mathbf{x}^{0}$ (we can assume that $k=0$ without
loss of generality), Step 1c. of Algorithm 3.1 first refreshes the auxiliary
variable by sampling
$\mathbf{u}\sim\operatorname{N}(\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log\pi(\mathbf{x}^{0}),\tfrac{\delta}{2}\mathbf{I})$
and then proposes
$\mathbf{x}^{1}\sim\operatorname{N}(\mathbf{u},\tfrac{\delta}{2}\mathbf{I})$.
The remaining steps return $\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{1}$ as the
new state with acceptance probability
$1\wedge\alpha_{\text{\lx@glossaries@gls@link{acronym}{AMALA}{\leavevmode
aMALA}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})$, where

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{AMALA}{\leavevmode
aMALA}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})$
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}=\frac{\pi(\mathbf{x}^{1})\operatorname{N}(\mathbf{u};\mathbf{x}^{1}+\tfrac{\delta}{2}\nabla\log\pi(\mathbf{x}^{1}),\tfrac{\delta}{2}\mathbf{I})\operatorname{N}(\mathbf{x}^{0};\mathbf{u},\tfrac{\delta}{2}\mathbf{I})}{\pi(\mathbf{x}^{0})\operatorname{N}(\mathbf{u};\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log\pi(\mathbf{x}^{0}),\tfrac{\delta}{2}\mathbf{I})\operatorname{N}(\mathbf{x}^{1};\mathbf{u},\tfrac{\delta}{2}\mathbf{I})}.$
(21)

Otherwise, the old state
$\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{0}=\mathbf{x}$ is returned as the new
state. This induces the same Markov chain on $\mathcal{X}$ as the aMALA from
Titsias and Papaspiliopoulos,  (2018) (the only difference relates to a re-
centring of the auxiliary variables $\mathbf{u}$ previously discussed in
Corenflos and Särkkä,  (2023) but this does not change the law of the Markov
chain on the marginal space which does not include the auxiliary variable).

  3. 3.

Generalisation of MALA. Still taking $\kappa=1$, the Particle-MALA generalises
the _Metropolis-adjusted Langevin algorithm (MALA)_ (Besag, , 1994) in the
sense that the former reduces to the latter if $T=N=1$. This can be seen as
follows, where use the same notational conventions as in the case of aMALA
above. Step 1c. of Algorithm 3.2 then marginally proposes
$\mathbf{x}^{1}\sim\operatorname{N}(\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log\pi(\mathbf{x}^{0}),\delta\mathbf{I})$.
The remaining steps return $\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{1}$ as the
new state with acceptance probability
$1\wedge\alpha_{\text{\lx@glossaries@gls@link{acronym}{MALA}{\leavevmode
MALA}}}(\mathbf{x}^{0},\mathbf{x}^{1})$, where

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{MALA}{\leavevmode
MALA}}}(\mathbf{x}^{0},\mathbf{x}^{1})$
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}=\frac{\pi(\mathbf{x}^{1})\operatorname{N}(\mathbf{x}^{0};\mathbf{x}^{1}+\tfrac{\delta}{2}\nabla\log\pi(\mathbf{x}^{1}),\delta\mathbf{I})}{\pi(\mathbf{x}^{0})\operatorname{N}(\mathbf{x}^{1};\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log\pi(\mathbf{x}^{0}),\delta\mathbf{I})}.$
(22)

Otherwise, the old state
$\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{0}=\mathbf{x}$ is returned as the new
state.

In particular, Remark 1 shows that we can view the aMALA as a ‘noisy’ version
of MALA (as already mentioned in Titsias and Papaspiliopoulos, , 2018)
because, dropping the time subscript again, by Lemma 3:

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{AMALA}{\leavevmode
aMALA}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})$
$\displaystyle=\alpha_{\text{\lx@glossaries@gls@link{acronym}{MALA}{\leavevmode
MALA}}}(\mathbf{x}^{0},\mathbf{x}^{1})\frac{\operatorname{N}(\mathbf{u};\bar{\mathbf{x}}+\tfrac{\delta}{4}\nabla\log\pi(\mathbf{x}^{1}),\tfrac{\delta}{4}\mathbf{I})}{\operatorname{N}(\mathbf{u};\bar{\mathbf{x}}+\tfrac{\delta}{4}\nabla\log\pi(\mathbf{x}^{0}),\tfrac{\delta}{4}\mathbf{I})},$
(23)

where $\bar{\mathbf{x}}=(\mathbf{x}^{0}+\mathbf{x}^{1})/2$, and hence

$\displaystyle\operatorname{\mathbb{E}}[\alpha_{\text{\lx@glossaries@gls@link{acronym}{AMALA}{\leavevmode
aMALA}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})]=\alpha_{\text{\lx@glossaries@gls@link{acronym}{MALA}{\leavevmode
MALA}}}(\mathbf{x}^{0},\mathbf{x}^{1}),$ (24)

where the expectation is w.r.t. the conditional distribution of $\mathbf{u}$
under the joint distribution of the random variables sampled in Step 1c. of
the Particle-aMALA, i.e. w.r.t.
$\operatorname{N}(\bar{\mathbf{x}}+\tfrac{\delta}{4}\nabla\log\pi(\mathbf{x}^{0}),\tfrac{\delta}{4}\mathbf{I})$.
In other words, this algorithm is the same as MALA except that the acceptance
ratio is ‘randomised’ in the sense that it is multiplied by a non-negative
random variable whose expectation is $1$. Other examples of such algorithms
can be found in Ceperley and Dewing,  (1999); Nicholls et al.,  (2012); see
also Finke,  (2015, Section 3.3.3) for a discussion as well as Andrieu and
Vihola,  (2016, page 2669) for a simple argument showing that the asymptotic
variance of aMALA cannot be smaller than that of MALA.

##  4 Particle extensions of mGRAD and aGRAD

###  4.1 Particle-aGRAD

The gradient-informed algorithms (Particle-MALA, etc) developed in Section 3
can be expected to improve upon the Particle-RWM algorithm in the same way
that aMALA/MALA improve upon the RWM algorithm. However, they may underperform
compared to the CSMC algorithm when the prior dynamics of the latent states
are highly informative in the same way that MALA can underperform relative to
the IMH algorithm (with prior as proposal) if the prior is highly informative.
Additionally, note that the algorithms from Section 3 employ proposals that
are _separable_ in the sense that, given the reference path, the marginal
proposal distribution of $\mathbf{x}_{t}^{n}$ does not depend on the ancestor
particle $\smash{\mathbf{x}_{t-1}^{a_{t-1}^{n}}}$ (that is, separability
implies that the weight-calculation and resampling steps could be postponed
until _after_ all particles have been proposed); such separable proposals can
be expected to perform poorly if the latent states are highly correlated
across time.

In this section, we further incorporate (conditionally) Gaussian prior
dynamics into the particle proposals and thus interpolate between the CSMC
algorithm and the gradient-informed algorithms of Section 3. Our construction
generalises the aGRAD and mGRAD algorithms of Titsias and Papaspiliopoulos,
(2018). In particular, the algorithms introduced in this section do not imply
separable proposals, i.e., the proposal kernel for particle
$\mathbf{x}_{t}^{n}$ will generally depend on its ancestor particle
$\smash{\mathbf{x}_{t-1}^{a_{t-1}^{n}}}$.

Specifically, in this section, we consider the special case of the generic
Feynman–Kac model from (1) in which we can find a decomposition
$Q_{t}(\mathbf{x}_{t-1:t})=M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})G_{t}(\mathbf{x}_{t-1:t})$,
such that

$\displaystyle
M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1})),$
(25)

is a Gaussian transition density whose mean $\mathbf{m}_{t}(\mathbf{x}_{t-1})$
and non-singular covariance matrix $\mathbf{C}_{t}(\mathbf{x}_{t-1})$ may
depend on the previous state $\mathbf{x}_{t-1}$, for $t>1$; and that
$M_{1}(\mathbf{x}_{1})=\operatorname{N}(\mathbf{x}_{1};\mathbf{m}_{1},\mathbf{C}_{1})$.

######  Example 2 (state-space model, continued).

The methods proposed in this section immediately apply with
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})\coloneqq
f_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})$ if the state-space model has
conditionally Gaussian dynamics, i.e. if
$f_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$,
by taking $G_{t}(\mathbf{x}_{t-1:t})=g_{t}(\mathbf{y}_{t}|\mathbf{x}_{t})$.
However, they may often still apply to state-space models with non-Gaussian
dynamics via a change of measure, i.e., by taking
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})\coloneqq\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$
and
$G_{t}(\mathbf{x}_{t-1:t})=f_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})g_{t}(\mathbf{y}_{t}|\mathbf{x}_{t})/\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$,
or through a suitable transformation.

The first method proposed in this section is termed _Particle-aGRAD_.
Conditional on the auxiliary variables $\mathbf{u}_{1:T}$, it can be viewed as
a CSMC algorithm whose proposal kernels are those of the fully-adapted
auxiliary particle filter for the state-space model defined by the Gaussian
transitions
$p(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$
from (25) and ‘pseudo observations’ $\mathbf{u}_{t}$ with
$p(\mathbf{u}_{t}|\mathbf{x}_{t})=\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\frac{\delta_{t}}{2}\mathbf{I})$.
We now write

$\displaystyle M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})$
$\displaystyle\coloneqq p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})$
(26)
$\displaystyle\propto\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})$
(27)
$\displaystyle\propto\operatorname{N}(\mathbf{x}_{t};\mathbf{m}^{\prime}_{t}(\mathbf{x}_{t-1},\mathbf{u}_{t}),\mathbf{C}^{\prime}_{t}(\mathbf{x}_{t-1})),\
$ (28) with $\displaystyle\mathbf{m}^{\prime}_{t}(\mathbf{x},\mathbf{u})$
$\displaystyle\coloneqq\mathbf{m}_{t}(\mathbf{x})+\mathbf{A}_{t}(\mathbf{x})[\mathbf{u}-\mathbf{m}_{t}(\mathbf{x})],$
(29) $\displaystyle\mathbf{C}^{\prime}_{t}(\mathbf{x})$
$\displaystyle\coloneqq(\mathbf{I}-\mathbf{A}_{t}(\mathbf{x}))\mathbf{C}_{t}(\mathbf{x})=\tfrac{\delta_{t}}{2}\mathbf{A}_{t}(\mathbf{x}),$
(30) $\displaystyle\mathbf{A}_{t}(\mathbf{x})$
$\displaystyle\coloneqq(\mathbf{C}_{t}(\mathbf{x})+\tfrac{\delta_{t}}{2}\mathbf{I})^{-1}\mathbf{C}_{t}(\mathbf{x}),$
(31) as well as $\displaystyle
G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$ $\displaystyle\coloneqq
Q_{t}(\mathbf{x}_{t-1:t})\frac{\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\frac{\delta_{t}}{2}\mathbf{I})}{M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})},$
(32)

and $Q^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$.
A single iteration of the Particle-aGRAD algorithm is as follows.

[Particle-aGRAD ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c) and the weight calculation (Step 1d) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\tfrac{\delta_{t}}{2}\mathbf{I})}$, and
$\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{t})}$, for
$n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

for $n\in[N]_{0}$, set $\smash{w_{t}^{n}\propto
G^{\prime}_{t}(\mathbf{x}_{t-1:t}^{(n)};\mathbf{u}_{t}),}$

and also replace $Q_{t+1}(\,\cdot\,)$ in the backward kernel in Step 3 by
$\smash{Q^{\prime}_{t+1}(\,\cdot\,;\mathbf{u}_{t})}$.

######  Proposition 4 (validity of Particle-aGRAD).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm 4.1
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves
$\pi_{T}$ invariant.

###  4.2 Particle-mGRAD

In this section, we analytically integrate out the auxiliary variables
$\mathbf{u}_{t}$ which appeared in the weights of the Particle-aGRAD
algorithm. Here we consider the case when the covariance matrices appearing in
the conditionally Gaussian mutation kernel (25) do not depend on the previous
state, i.e.,

$\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t},$ (33)

which then also implies that
$\mathbf{A}_{t}(\mathbf{x}_{t-1})=\mathbf{A}_{t}$. A single iteration of the
resulting methodology – which we term the _Particle-mGRAD_ algorithm – is as
follows, where we write

$\displaystyle\log
H_{t,\bm{\phi}}(\mathbf{x},\mathbf{v},\bar{\mathbf{x}},\bar{\mathbf{v}})$
$\displaystyle=\tfrac{1}{2}(\mathbf{x}-\mathbf{v})^{\mathrm{T}}((\tfrac{\delta_{t}}{2}\mathbf{A}_{t})^{-1}+\mathbf{G}_{t})(\mathbf{x}-\mathbf{v})$
$\displaystyle\qquad-[\tfrac{1}{2}N(\mathbf{x}+\bm{\phi})^{\mathrm{T}}\mathbf{A}_{t}+(\mathbf{x}-\mathbf{v})^{\mathrm{T}}]\mathbf{G}_{t}(\mathbf{x}+\bm{\phi})$
$\displaystyle\qquad+(N+1)(\bar{\mathbf{x}}-\bar{\mathbf{v}})^{\mathrm{T}}\mathbf{G}_{t}(\mathbf{v}+\bm{\phi}),$

for
$\mathbf{G}_{t}\coloneqq\tfrac{2}{\delta_{t}}(\mathbf{I}+N\mathbf{A}_{t})^{-1}$.

[Particle-mGRAD ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c) and the weight calculation (Step 1d) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\tfrac{\delta_{t}}{2}\mathbf{I})}$ and
$\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{t})}$, for
$n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

set
$\smash{\bar{\mathbf{x}}_{t}\coloneqq\tfrac{1}{N+1}\sum_{n=0}^{N}\mathbf{x}_{t}^{n}}$,
$\mathbf{v}_{t}^{n}\coloneqq(\mathbf{I}-\mathbf{A}_{t})\mathbf{m}_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{n}})$,
$\smash{\bar{\mathbf{v}}_{t}\coloneqq\tfrac{1}{N+1}\sum_{n=0}^{N}\mathbf{v}_{t}^{n}}$,
and, for $n\in[N]_{0}$,

$\displaystyle w_{t}^{n}\propto
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})H_{t,\kappa\frac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}^{n}}\log
G_{t}(\mathbf{x}_{t-1:t}^{(n)})}(\mathbf{x}_{t}^{n},\mathbf{v}_{t}^{n},\bar{\mathbf{x}}_{t},\bar{\mathbf{v}}_{t}).$
(34)

######  Remark 2 (Particle-aGRAD ‘exactly approximates’ Particle-mGRAD).

In analogue to the relationship between Particle-aMALA and Particle-MALA
discussed in Remark 1, Particle-aGRAD is a noisy version of Particle-mGRAD.
That is, letting $w_{t}^{n}(\mathbf{u}_{t})$ and $w_{t}^{n}$ be the
unnormalised weights under Particle-aGRAD and Particle-mGRAD, respectively, we
have

$\displaystyle\operatorname{\mathbb{E}}\biggl{[}\frac{w_{t}^{n}(\mathbf{u}_{t})}{w_{t}^{k_{t}}(\mathbf{u}_{t})}\biggr{]}=\frac{w_{t}^{n}}{w_{t}^{k_{\mathrlap{t}}}},$
(35)

where the expectation is taken with respect to the conditional distribution of
$\mathbf{u}_{t}$ under the joint distribution of all random variables
generated by Algorithm 4.1 up to (and including) time $t$.

######  Proposition 5 (validity of Particle-mGRAD).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm 4.2
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves
$\pi_{T}$ invariant.

###  4.3 Particle-aGRAD+

While the algorithm of Section 4.1 incorporates information from the smoothing
distribution by merit of not modifying the latent dynamics, it may happen that
the potential $G_{t}(\mathbf{x}_{t-1:t})$ strongly depends on
$\mathbf{x}_{t-1}$. In this case, considering the ‘myopic’ gradient
information $\nabla_{\mathbf{x}_{t}}\log G_{t}(\mathbf{x}_{t-1:t})$ may not
suffice to improve the mixing of the algorithm and information from
$\mathbf{x}_{t-1}$ may then be beneficial. Similarly to Section 3.3, in this
section, we extend the Particle-aGRAD algorithm to incorporate gradients
w.r.t. the ‘smoothing potential’
$G_{1:T}(\mathbf{x}_{1:T})=\prod_{t=1}^{T}G_{t}(\mathbf{x}_{t-1:t})$ rather
than w.r.t. the ‘filtering potential’
$\prod_{s=1}^{t}G_{s}(\mathbf{x}_{s-1:s})$.

For $M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})$ and
$G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$ still defined as in the
Particle-aGRAD algorithm (i.e., as in (28) and (32)), we now write

$\displaystyle G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})$
$\displaystyle\coloneqq
G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})\frac{\operatorname{N}(\mathbf{u}_{t-1};\mathbf{x}_{t-1}+\kappa\tfrac{\delta_{t-1}}{2}\nabla_{\mathbf{x}_{t-1}}\log
G_{1:T}(\mathbf{x}_{1:T}),\tfrac{\delta_{t-1}}{2}\mathbf{I})}{\operatorname{N}(\mathbf{u}_{t-1};\mathbf{x}_{t-1}+\kappa\tfrac{\delta_{t-1}}{2}\nabla_{\mathbf{x}_{t-1}}\log
G_{t-1}(\mathbf{x}_{t-2:t-1}),\tfrac{\delta_{t-1}}{2}\mathbf{I})},\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$

as well as $Q^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})$,
where we note that

$\displaystyle\nabla_{\mathbf{x}_{t}}\log
G_{1:T}(\mathbf{x}_{1:T})=\nabla_{\mathbf{x}_{t}}\log[G_{t}(\mathbf{x}_{t-1:t})+\log
G_{t+1}(\mathbf{x}_{t:t+1})].$ (36)

A single iteration of the resulting ‘smoothing-gradient’ methodology – which
we term the _Particle-aGRAD+_ algorithm – is as follows.

[Particle-aGRAD+ ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c), the weight calculation (Step 1d), and backward sampling (Step 3) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{1:T}(\mathbf{x}_{1:T}),\tfrac{\delta_{t}}{2}\mathbf{I})}$, and
$\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})}$, for
$n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

for $n\in[N]_{0}$, set $\smash{w_{t}^{n}\propto
G^{\prime}_{t}(\mathbf{x}_{t-2:t}^{(n)};\mathbf{u}_{t-1:t}),}$

  3. 3.

for $t=T-1,\dotsc,1$, sample $l_{t}=i\in[N]_{0}$ w.p.

$\displaystyle\dfrac{W_{t}^{i}Q^{\prime}_{t+1}((\mathbf{x}_{t-1:t}^{(i)},\mathbf{x}_{t+1}^{l_{t+1}});\mathbf{u}_{t:t+1})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{i},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{t+1:t+2})}{\sum_{n=0}^{N}W_{t}^{n}Q^{\prime}_{t+1}((\mathbf{x}_{t-1:t}^{(n)},\mathbf{x}_{t+1}^{l_{t+1}});\mathbf{u}_{t:t+1})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{n},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{t+1:t+2})}.$
(37)

Note that if $G_{t}(\mathbf{x}_{t-1:t})=G_{t}(\mathbf{x}_{t})$ does not depend
on $\mathbf{x}_{t-1}$, then the Particle-aGRAD+ algorithm coincides with the
Particle-aGRAD algorithm. However, when $G_{t}(\mathbf{x}_{t-1:t})$ varies
highly in $\mathbf{x}_{t-1}$, their behaviours may differ substantially.

######  Proposition 6 (validity of Particle-aGRAD+).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm 4.3
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRADPLUS}{\leavevmode
Particle-aGRAD+}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves
$\pi_{T}$ invariant.

###  4.4 Twisted Particle-aGRAD(+)

Recall that, conditionally on the auxiliary variables $\mathbf{u}_{1:T}$, the
Particle-aGRAD algorithm could be viewed as a CSMC algorithm whose proposal
kernels
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})=p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})$
are those of the fully-adapted auxiliary particle filter for the state-space
model which is defined by the Gaussian transitions
$p(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$
from (25) and observation densities
$p(\mathbf{u}_{t}|\mathbf{x}_{t})=\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\frac{\delta_{t}}{2}\mathbf{I})$.

In this section (and in this section only), we make the more restrictive
assumption that the transition kernel from (25) is not only Gaussian but also
affine, i.e.,

$\displaystyle\mathbf{m}_{t}(\mathbf{x}_{t-1})=\mathbf{F}_{t}\mathbf{x}_{t-1}+\mathbf{b}_{t},\quad\text{and}\quad\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t},$
(38)

for some $\mathbf{F}_{t}\in\mathbb{R}^{D\times D}$,
$\mathbf{b}_{t}\in\mathbb{R}^{D}$, and some covariance matrix
$\mathbf{C}_{t}\in\mathbb{R}^{D\times D}$. Under (38), we can then go one step
further and implement the fully _twisted_ particle filter (Whiteley and Lee, ,
2014; Guarniero et al., , 2017; Heng et al., , 2020) proposal which conditions
on _all_ future pseudo observations $\mathbf{u}_{t:T}$. That is, we now write

$\displaystyle
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})$
$\displaystyle\coloneqq p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t:T})$
(39)
$\displaystyle\propto\int_{\mathcal{X}^{T-t}}\biggl{[}\prod_{s=t}^{T}\operatorname{N}(\mathbf{x}_{s};\mathbf{F}_{s}\mathbf{x}_{s-1}+\mathbf{b}_{s},\mathbf{C}_{s})\operatorname{N}(\mathbf{u}_{s};\mathbf{x}_{s},\tfrac{\delta_{s}}{2}\mathbf{I})\biggr{]}\,\mathrm{d}\mathbf{x}_{t+1:T}$
(40)
$\displaystyle\propto\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t}),$
(41) $\displaystyle G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t:T})$
$\displaystyle\coloneqq
Q_{t}(\mathbf{x}_{t-1:t})\frac{\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\frac{\delta_{t}}{2}\mathbf{I})}{M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})},$
(42)

as well as $Q^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t:T})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t:T})$.
Here, $\mathbf{b}^{\prime}_{t}\in\mathbb{R}^{D}$ and
$\mathbf{F}^{\prime}_{t},\mathbf{C}^{\prime}_{t}\in\mathbb{R}^{D\times D}$ can
be obtained via Kalman-filtering recursions as explained in Appendix B.

A single iteration of the resulting methodology – which we term the _twisted
Particle-aGRAD_ algorithm – is then exactly as the Particle-aGRAD (Algorithm
4.1), except that $M^{\prime}_{t}(\,\cdot\,|\,\cdot\,;\mathbf{u}_{t})$,
$G^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t})$ and
$Q^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t})$ from Section 4.1 are replaced by
$M^{\prime}_{t}(\,\cdot\,|\,\cdot\,;\mathbf{u}_{t:T})$,
$G^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t:T})$ and
$Q^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t:T})$ from this section. When the
potential functions $G_{t}(\mathbf{x}_{t-1:t})$ vary in $\mathbf{x}_{t-1}$,
then we can further construct a _twisted Particle-aGRAD+_ algorithm by
replacing $M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})$ in
Algorithm 4.3 and in the denominator of
$G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})$ by
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t:T})$.

######  Proposition 7 (validity of the twisted Particle-aGRAD/Particle-
aGRAD+).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via the twisted
Particle-aGRAD or twisted Particle-aGRAD+ algorithm induces a Markov kernel
which leaves $\pi_{T}$ invariant.

###  4.5 Relationship with other methods

The algorithms proposed above relate to existing methods as follows.

  1. 1.

Generalisation of aGRAD. For $\kappa=1$, the Particle-aGRAD algorithm (and
similarly the Particle-aGRAD+ algorithm as well as the twisted versions of
either) generalises the _auxiliary gradient (aGRAD)_ algorithm from Titsias
and Papaspiliopoulos,  (2018, called ‘aGrad-z’ therein) in the sense that the
former reduces to the latter if $T=N=1$. This can be seen as follows, where we
again suppress the ‘time’ subscript $t=1$ everywhere so that
$\pi(\mathbf{x})\propto M(\mathbf{x})G(\mathbf{x})$, where
$M(\mathbf{x})=\operatorname{N}(\mathbf{x};\mathbf{m},\mathbf{C})$. Given that
the current state of the Markov chain is $\mathbf{x}=\mathbf{x}^{0}$ (we can
assume that $k=0$ without loss of generality), Step 1c. of Algorithm 4.1 first
refreshes the auxiliary variable by sampling
$\mathbf{u}\sim\operatorname{N}(\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{0}),\tfrac{\delta}{2}\mathbf{I})$ and then proposes
$\mathbf{x}^{1}\sim\operatorname{N}((\mathbf{I}-\mathbf{A})\mathbf{m}+\mathbf{A}\mathbf{u},\tfrac{\delta}{2}\mathbf{A})$,
for $\mathbf{A}=(\mathbf{C}+\tfrac{\delta}{2}\mathbf{I})^{-1}\mathbf{C}$. The
remaining steps return $\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{1}$ as the new
state with acceptance probability
$1\wedge\alpha_{\text{\lx@glossaries@gls@link{acronym}{AGRAD}{\leavevmode
aGRAD}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})$, where

$\displaystyle\\!\\!\\!\\!\alpha_{\text{\lx@glossaries@gls@link{acronym}{AGRAD}{\leavevmode
aGRAD}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})$
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}$ (43)
$\displaystyle=\frac{\pi(\mathbf{x}^{1})\operatorname{N}(\mathbf{u};\mathbf{x}^{1}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{1}),\tfrac{\delta}{2}\mathbf{I})\operatorname{N}(\mathbf{x}^{0};(\mathbf{I}-\mathbf{A})\mathbf{m}+\mathbf{A}\mathbf{u},\tfrac{\delta}{2}\mathbf{A})}{\pi(\mathbf{x}^{0})\operatorname{N}(\mathbf{u};\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{0}),\tfrac{\delta}{2}\mathbf{I})\operatorname{N}(\mathbf{x}^{1};(\mathbf{I}-\mathbf{A})\mathbf{m}+\mathbf{A}\mathbf{u},\tfrac{\delta}{2}\mathbf{A})}.\\!\\!\\!\\!$
(44)

Otherwise, the old state
$\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{0}=\mathbf{x}$ is returned as the new
state.

  2. 2.

Generalisation of mGRAD. Still taking $\kappa=1$, the Particle-mGRAD algorithm
generalises the _marginal gradient (mGRAD)_ algorithm from Titsias and
Papaspiliopoulos,  (2018) in the sense that the former reduces to the latter
if $T=N=1$. This can be seen as follows, where we use the same notational
conventions as in the case of aGRAD above. Step 1c. of Algorithm 4.2 then
marginally proposes
$\mathbf{x}^{1}\sim\operatorname{N}((\mathbf{I}-\mathbf{A})\mathbf{m}+\mathbf{A}[\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{0})],\mathbf{B})$, where
$\mathbf{B}\coloneqq\tfrac{\delta}{2}\mathbf{A}^{2}+\mathbf{A}$. The remaining
steps return $\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{1}$ as the new state with
acceptance probability
$1\wedge\alpha_{\text{\lx@glossaries@gls@link{acronym}{MGRAD}{\leavevmode
mGRAD}}}(\mathbf{x}^{0},\mathbf{x}^{1})$, where

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{MGRAD}{\leavevmode
mGRAD}}}(\mathbf{x}^{0},\mathbf{x}^{1})$
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}=\frac{\pi(\mathbf{x}^{1})\operatorname{N}((\mathbf{I}-\mathbf{A})\mathbf{m}+\mathbf{A}[\mathbf{x}^{1}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{1})],\mathbf{B})}{\pi(\mathbf{x}^{0})\operatorname{N}((\mathbf{I}-\mathbf{A})\mathbf{m}+\mathbf{A}[\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{0})],\mathbf{B})}.$ (45)

Otherwise, the old state
$\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{0}=\mathbf{x}$ is returned as the new
state. In particular, by Remark 2, in analogue to Section 3.4, we can again
interpret aGRAD as a version of mGRAD with ‘randomised’ acceptance ratio.

  3. 3.

Generalisation of a ‘preconditioned’ Particle-RWM algorithm. If $\kappa=0$,
then the Particle-aGRAD and Particle-aGRAD+ algorithms reduce to a method
recently proposed in Corenflos and Särkkä,  (2023, Section 4.3), which can be
seen as a ‘preconditioned’ version of the Particle-RWM algorithm.

###  4.6 Interpolation between CSMC and Particle-MALA/Particle-aMALA

The Particle-MALA (and related methods) proposed in Section 3 may be
outperformed by the CSMC algorithm in the case when the prior dynamics are
highly informative – in the same way that MALA may be outperformed by the IMH
algorithm (with prior as proposal) if the prior dominates the posterior. For
instance, in the extreme case that all the potential functions are constant,
the CSMC algorithm proposes $N$ trajectories (in addition to the reference
path) that are IID samples from $\pi_{T}$ (assuming an adaptive or low-
variance conditional resampling scheme is used) while the $N$ trajectories
proposed by Particle-MALA are still highly correlated with the reference path.

Put differently, the user is faced with the ‘tuning problem’ of having to
decide between the CSMC algorithm on the one hand and the Particle-MALA (and
related methods) on the other hand. In this section, we show that the
Particle-mGRAD algorithm resolves this tuning problem in the sense that it can
be viewed as interpolating between CSMC and Particle-MALA. Specifically,
Proposition 8 shows that Particle-mGRAD reduces to the CSMC algorithm if the
prior dynamics are highly informative. Conversely, Proposition 9 shows that
Particle-mGRAD reduces to the Particle-MALA if the prior dynamics are
uninformative. The same results hold for the auxiliary-variable versions:
Particle-aMALA and Particle-aGRAD.

We make the following assumptions (assumed to hold for all $t\in[T]$):

  1. A1

For any $\mathbf{x}_{t}\in\mathcal{X}$,
$\mathbf{m}_{t}(\mathbf{x}_{t-1})=\mathbf{m}_{t}$,
$\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t}$ and
$G_{t}(\mathbf{x}_{t-1:t})=G_{t}(\mathbf{x}_{t})$ are constant in
$\mathbf{x}_{t-1}$, with $G_{t}$ uniformly bounded on $\mathcal{X}$ and
$\mathbf{C}_{t}$ invertible.

  2. A2

There exist $C_{0},C_{1}\geq 0$ such that $\kappa\lVert\nabla\log
G_{t}(\mathbf{x}_{t})\rVert_{2}\leq
C_{0}+C_{1}\lVert\mathbf{x}_{t}\rVert_{2}$.

  3. A3

$\max_{d\in[D]}\int_{\mathcal{X}}x_{t,d}^{2}G_{t}(\mathbf{x}_{t})\,\mathrm{d}\mathbf{x}_{t}<\infty$,
where $x_{t,d}$ is the $d$th component of $\mathbf{x}_{t}$.

Whenever $T>1$, Assumption A1 is strong because it requires the Feynman–Kac
model to factorise over time. However, we expect that it could be relaxed at
the cost of greatly complicating the arguments. Indeed, note that the model
used in Figure 2 does not satisfy this assumption. Assumption A2 is rather
mild, e.g. it holds in a state-space model with Gaussian measurement errors.

In the following, for each $t\in[T]$, we will consider a sequence of prior
covariance matrices $(\mathbf{C}_{t,k})_{k\geq 1}$. We will therefore add the
subscript $k$ to any quantity which depends on $\mathbf{C}_{t,k}$. We also let
$\lambda(\mathbf{A})$ denote the set of eigenvalues of some matrix
$\mathbf{A}$. The following propositions are proved in Appendix E.

######  Proposition 8.

For some $D,T,N\geq 1$, assume A1–A2, and assume that there exists a sequence
$(\lambda_{k})_{k\geq 1}$ in $(0,\infty)$ with
$\max\\{\lambda(\mathbf{C}_{1,k}),\dotsc,\lambda(\mathbf{C}_{T,k})\\}\leq\lambda_{k}\to
0$ as $k\to\infty$. Then for any $\varepsilon>0$, there exists a sequence
$(F_{T,k})_{k\geq 1}$ of subsets of $\mathcal{X}^{T}$ with
$\lim_{k\to\infty}\pi_{T,k}(F_{T,k})=1$ such that

  1. 1.

$\smash{\sup_{\mathbf{x}_{1:T}\in F_{T,k}}\lVert
P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode Particle-
mGRAD}},k}(\,\cdot\,|\mathbf{x}_{1:T})-P_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},k}(\,\cdot\,|\mathbf{x}_{1:T})\rVert_{\textnormal{{tv}}}\in\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)/4})}$;

  2. 2.

$\smash{\sup_{\mathbf{x}_{1:T}\in F_{T,k}}\lVert
P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode Particle-
aGRAD}},k}(\,\cdot\,|\mathbf{x}_{1:T})-P_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},k}(\,\cdot\,|\mathbf{x}_{1:T})\rVert_{\textnormal{{tv}}}\in\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)/4})}$.

######  Proposition 9.

For some $D,T,N\geq 1$, assume A1–A3, and assume that there exists a sequence
$(\lambda_{k})_{k\geq 1}$ in $(0,\infty)$ with
$\min\\{\lambda(\mathbf{C}_{1,k}),\dotsc,\lambda(\mathbf{C}_{T,k})\\}\geq\lambda_{k}\to\infty$
as $k\to\infty$. Then for any $\varepsilon>0$, there exists a sequence
$(F_{T,k})_{k\geq 1}$ of subsets of $\mathcal{X}^{T}$ with
$\lim_{k\to\infty}\pi_{T,k}(F_{T,k})=1$ such that

  1. 1.

$\smash{\sup_{\mathbf{x}_{1:T}\in F_{T,k}}\lVert
P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode Particle-
mGRAD}},k}(\,\cdot\,|\mathbf{x}_{1:T})-P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},k}(\,\cdot\,|x)\rVert_{\textnormal{{tv}}}\in\operatorname{\mathrm{O}}(\lambda_{k}^{-(1-\varepsilon)/4})}$;

  2. 2.

$\smash{\sup_{\mathbf{x}_{1:T}\in F_{T,k}}\lVert
P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode Particle-
aGRAD}},k}(\,\cdot\,|\mathbf{x}_{1:T})-P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-
aMALA}},k}(\,\cdot\,|x)\rVert_{\textnormal{{tv}}}\in\operatorname{\mathrm{O}}(\lambda_{k}^{-(1-\varepsilon)/4})}$.

As per Sections 2.1.2, 3.4 and 4.5, taking $T=N=1$ in Propositions 8 and 9
immediately imply that the aGRAD/mGRAD algorithm can be viewed as
automatically interpolating between the IMH algorithm with prior as proposal
(if the prior is highly informative) and aMALA/MALA (if the prior is highly
diffuse). To our knowledge, this interpretation has not been pointed out in
the literature. It provides new intuition for the noteworthy performance of
aGRAD/mGRAD in Titsias and Papaspiliopoulos,  (2018).

###  4.7 Complexity

An iteration of Particle-aGRAD or Particle-aGRAD+ requires computing $T(N+1)$
gain matrices
$\smash{\mathbf{A}_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{n}})\in\mathbb{R}^{D\times
D}}$; and all of these, in general, have a cubic cost in the latent-state
dimension $D$. While this may be reasonable for small enough systems and will
be helpful for informative likelihoods, the computational quickly outweighs
the statistical benefits of the method. However, when the dynamics have
additive noise (33), $\mathbf{A}_{t}$ does not depend on $\mathbf{x}_{t-1}$.
In this case, only $T$ gain matrices are needed and these can be pre-computed,
only paying the cubic cost in the dimension upfront rather than at each
iteration.

The same applies for the Particle-mGRAD algorithm for which we always require
(33) to hold (the auxiliary variables could still be integrated out if (33) is
relaxed, but only at the cost of a cubic computational complexity in the
number of particles).

However, as for the Particle-RWM algorithm and Particle-MALA-type methods, we
need to calibrate the step-size parameters $\delta_{t}$ which changes the gain
matrices (so that pre-computation is not possible during the calibration
stage). Thankfully, because $\mathbf{A}_{t}$ and $\mathbf{C}_{t}$ have the
same eigenvectors no matter what $\delta_{t}$ is, it is possible to use
similar spectral methods as in Titsias and Papaspiliopoulos,  (2018) to reduce
the complexity of changing $\delta_{t}$ to quadratic.

At first sight, the complexity of the twisted Particle-aGRAD seems quadratic
in $T$ as the proposal kernel
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t:T})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t})$
requires processing $T-t$ auxiliary variables for each time $t$. However, in
Appendix B, we show how $\mathbf{F}^{\prime}_{t}$, $\mathbf{b}^{\prime}_{t}$
and $\mathbf{C}^{\prime}_{t}$ can all be pre-computed based on standard Kalman
filter recursions (Kalman, , 1960), preserving the linear cost in $T$ and $N$.

##  5 Experimental validation and comparison

###  5.1 Multivariate stochastic volatility model

In this section, we illustrate the efficiency of our methods on a multivariate
stochastic volatility model often used as a benchmark for high-dimensional
sequential Monte Carlo methodology (see, e.g., Guarniero et al., , 2017). This
model is a state-space model with a non-linear observation equation:

$\displaystyle
g_{t}(\mathbf{y}_{t}|\mathbf{x}_{t})=\operatorname{N}(\mathbf{y}_{t};\bm{0},\operatorname{diag}(\exp\mathbf{x}_{t})),$
(46)

where $\exp$ is applied element-wise and where $\bm{0}$ is a $D$-dimensional
vector of zeros. The prior on the latent variables is defined through auto-
regressive Gaussian dynamics, i.e. for $t>1$:

$\displaystyle
f_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t})$
(47)

where $\mathbf{m}_{t}(\mathbf{x}_{t-1})\coloneqq\varphi\mathbf{x}_{t-1}$ and
$\mathbf{C}_{t}\in\smash{\mathbb{R}^{D\times D}}$ has diagonal entries $\tau$
and off-diagonal entries $\tau\rho$. The initial distribution
$f_{1}(\mathbf{x}_{1})=\operatorname{N}(\mathbf{x}_{1};\mathbf{m}_{1},\mathbf{C}_{1})$
is the stationary distribution under the dynamics (47), i.e.,
$\mathbf{m}_{1}\coloneqq\bm{0}$ and
$\mathbf{C}_{1}\coloneqq\mathbf{C}_{t}/(1-\varphi^{2})$. Here,
$\varphi\in(-1,1)$ is some autocorrelation coefficient, $\rho\in(-1,1)$ is
some intra-asset correlation coefficient and $\tau>0$.

Throughout our experiments, we take $\varphi=0.9$, $\rho=0.25$, and
$\tau\in\\{0.1,0.5,1,2\\}$. The eigenvalues of $\mathbf{C}_{t}$ are then
proportional to $\tau$, i.e., a small value of $\tau$ corresponds to highly
informative prior dynamics (as in Proposition 8) while a large value of $\tau$
corresponds to weakly informative prior dynamics (as in Proposition 9). To
make our observations robust to the choice of data set, for each $\tau$, we
simulated $M=5$ independent sets of $T=128$ observations from the multivariate
stochastic volatility model with $D=30$, i.e., each state $\mathbf{x}_{t}$
takes values in $\mathcal{X}=\mathbb{R}^{30}$. To make results more easily
comparable, experiments for different values of $\tau$ use the same random
number generator seed.

###  5.2 Simulation study setup

In addition to the methods proposed in Sections 3 and 4 – potentially without
the use of gradient information by taking $\kappa=0$ – we consider the
following benchmark methods:

  1. 1.

CSMC. The CSMC algorithm with bootstrap proposals (Algorithm 2.1.1).

  2. 2.

Particle-RWM. The Particle-RWM algorithm (Algorithm 2.2.1) from Finke and
Thiery,  (2023) (the special case of Particle-aMALA/Particle-MALA/Particle-
aMALA+ if $\kappa=0$).

  3. 3.

MALA and aMALA. The $N$-proposal MALA and aMALA which correspond to the
Particle-aMALA and Particle-MALA proposed in this work with a single time step
(applied to the path-space representation of the Feynman–Kac model, i.e. with
a single $(D\times T)$-dimensional state).

  4. 4.

aGRAD. The $N$-proposal aGRAD algorithm, which corresponds to the Particle-
aGRAD proposed in this work with a single time step (again on the path space).
We note that we implemented aGRAD using the auxiliary Kalman perspective of
Corenflos and Särkkä,  (2023), making the method complexity scale linearly
with $T$ rather than quadratically with $T$ as in the original version of
Titsias and Papaspiliopoulos,  (2018). We do not compare to mGRAD because
computing its particle weights (and hence acceptance ratio) has quadratic
complexity in $T$.

All algorithms use $N+1=32$ particles, and those employing resampling use the
conditional ’killing’ resampling method (Karppinen et al., , 2023), more
stable than multinomial resampling, especially with highly informative priors.
In each of $M=5$ independent experiments, algorithms start from the same
trajectory generated by a bootstrap particle filter using $32$ particles. The
samplers run for $10\,000$ steps to calibrate step-size parameters
$\delta_{t}$, detailed below (note that calibration stabilises much faster).
For CSMC, which requires no calibration, the initial $10\,000$ steps are
discarded as warm-up. After calibration, $J=4$ independent chains start at the
final calibration sample, running for $K=$50\,000$$ iterations, with the first
$5000$ discarded as burn-in to decorrelate the chains. Reported statistics are
based on these $J$ independent chains.

The step-size parameters $\delta_{t}$ are calibrated for a
$75\text{\,}\mathrm{\char 37\relax}\text{/}$ acceptance rate, as explained in
Appendix F. This slightly exceeds recommendations by, e.g., Roberts and
Rosenthal,  (2001); Titsias and Papaspiliopoulos,  (2018). This is because we
use multiple proposals and the optimal acceptance rate is expected to increase
accordingly. Here, ’acceptance rate at time $t$’ refers to the relative
frequency of with which the state $\mathbf{x}_{t}$ is updated. Figure 7 in
Appendix G.1 shows stable acceptance rates around $75\text{\,}\mathrm{\char
37\relax}\text{/}$ for all methods except CSMC across all time steps. Figure 6
in Appendix G.1 displays calibrated $\delta_{t}$ values.

Experiments ran on a shared computational cluster with identical
configurations (32 GB RAM, four processor cores, on shared machines with
$$2$\times$64$$-core AMD EPYC 7713 CPUs, clock speed 2.0 GHz). Nonetheless,
cluster idiosyncrasies may be present, potentially impacting slower methods
like Particle-aMALA+ and Particle-mGRAD.

###  5.3 Breakdown of CSMC, aMALA and MALA

Our results indicate that CSMC, aMALA, and MALA failed to explore the right
regions of the space for all of our chosen levels of informativeness of the
latent dynamics ($\tau\in\\{0.1,0.5,1,2\\}$). Specifically, Figures 8 and 9 in
Appendix G.2 show that both the estimated marginal posterior means and also
the energy traces of CSMC, aMALA and MALA differ substantially from those of
all the other algorithms. Here, ‘energy trace’ refers to
$\log\pi_{T}(\mathbf{x}_{1:T})+\mathrm{const}$ computed on the sampled
trajectories throughout the sampling procedure. Since CSMC, aMALA and MALA
thus do not produce reliable approximations of the distribution of interest,
we omit these methods from our discussions in the sequel.

In the remainder of this section, we compare the remaining algorithms in terms
of the _effective sample size (ESS)_ computed using the method of Vehtari et
al.,  (2021) with $J=4$ independent chains. We also compare the algorithms in
terms of ESS per second (ESS/s). The latter corresponds to the time it would
take to obtain a ‘perfect’ sample using the Markov chain. In the main
manuscript, we only show results for the median ESS and averaged over all $T$
time steps. Appendix G.3 shows detailed results for the minimum and maximum
ESS and ESS/s (which are qualitatively similar to the median case) separately
for each time step $t=1,\dotsc,T$.

###  5.4 Benefits of exploiting gradient information

Figure 3 compares the median ESS (‘unnormalised’) and median ESS/s (‘per
second’) of Particle-aMALA, Particle-MALA and Particle-aMALA+, i.e., for those
methods which do not make any Gaussian assumption about the prior dynamics.
Recall that these differ from the baseline: the Particle-RWM algorithm, only
in the use of gradient information. Thus, the left panel in Figure 3
illustrates the benefits (in terms of ESS) of exploiting gradient information.
Notably:

  * •

the improvement of Particle-MALA over Particle-aMALA is marginal at best.
Possibly, the difference between both algorithms decreases with $N$ but this
calls for further investigation;

  * •

the ‘smoothing-gradient’ variant Particle-aMALA+ dominates all other
alternatives for all values of $\tau$, with up to three times the performance
of Particle-RWM and twice that of the ‘filter-gradient’ variants Particle-
aMALA and Particle-MALA;

  * •

the performance of all shown methods improves as $\tau$ increases: this is
because the posterior distribution then decorrelates in time, and, therefore,
the fact that they all use proposals which are separable (in the sense
discussed in Section 4) stops being penalising.

The right panel in Figure 3 shows that the use of gradient information is
still beneficial even when accounting for the cost of gradient calculation.
However, the relative performance of the gradient-based methods is now less
clear: whilst Particle-aMALA+ has the highest sampling efficiency, it incurs
additional overheads due to computing twice as many gradients as Particle-
aMALA and Particle-MALA and due to dealing with non-Markovian potentials.

Figure 3: Performance of those proposed methods which do not require
(conditionally or unconditionally) Gaussian prior dynamics compared with the
existing Particle-RWM algorithm as a baseline.

###  5.5 Benefits of exploiting Gaussian prior dynamics

In this section, we demonstrate that exploiting the latent (conditionally)
Gaussian dynamics of the model (as done by Particle-aGRAD, Particle-mGRAD and
twisted Particle-aGRAD) can improve the sampling efficiency.

First, in Figure 4, we illustrate the performance of those methods which
require (at most) _conditionally_ Gaussian prior dynamics as in (25), i.e., of
Particle-aGRAD and Particle-mGRAD (note that the later also requires
$\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t}$ (33) to retain linear
computational complexity in $N$). In terms of ESS, these methods improve upon
the ‘filter-gradient’ methods Particle-aMALA and Particle-MALA but they are
still dominated by the ‘smoothing-gradient’ method Particle-aMALA+. However,
the picture is less clear when accounting for computation time.

Figure 4: Performance of the proposed methods which require only
_conditionally_ Gaussian prior dynamics (25), i.e.,
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$.
The Particle-mGRAD algorithm (with any $\kappa\in\\{0,1\\}$) also requires
that $\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t}$ is constant to avoid
superlinear computational complexity in $N$. Results that were already shown
in the previous figure are greyed out.

Second, in Figure 5, we illustrate the performance of the twisted Particle-
aGRAD which requires _unconditionally_ Gaussian prior dynamics as in (38). As
a baseline, we use the aGRAD algorithm from Titsias and Papaspiliopoulos,
(2018) as it makes the same assumption. The twisted Particle-aGRAD strongly
outperforms this baseline and also all the other algorithms. Furthermore, the
dominance of the twisted Particle-aGRAD algorithm does not disappear when
accounting for the computation time. This is because, in contrast to Particle-
aMALA+, its modified model is still Markovian and because it only requires the
computation of a single gradient per particle and time step.

Figure 5: Performance of the proposed methods which require _unconditionally_
Gaussian prior dynamics (38), i.e.,
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}_{t}\mathbf{x}_{t-1}+\mathbf{b}_{t},\mathbf{C}_{t})$,
compared with aGRAD (which also requires (38)) as baseline. Results that were
already shown in the previous two figures are greyed out. The abnormally large
computation time of the twisted Particle-aGRAD for $\kappa=\tau=1$ was likely
caused by some computational-cluster idiosyncrasies.

##  6 Conclusion

###  6.1 Summary

We have proposed a methodology for Bayesian inference about the latent states
in high-dimensional state-space models and beyond. Our methodology combines
the CSMC algorithm (Andrieu et al., , 2010) with sophisticated ‘classical’
MCMC algorithms like MALA (Besag, , 1994), aMALA (Titsias and
Papaspiliopoulos, , 2018), aGRAD/mGRAD (Titsias, , 2011; Titsias and
Papaspiliopoulos, , 2018) or PCNL (Cotter et al., , 2013) to retain _the best
of both worlds:_

  * •

from the CSMC algorithm, our methods retain the ability to exploit the model’s
‘decorrelation-over-time’ structure which permits favourable scaling with the
number of time steps, $T$;

  * •

from ‘classical’ MCMC algorithms, our methods retain the ability to use
gradient-informed, local proposals which permits favourable scaling with the
dimension of the states, $D$.

Most of our proposed algorithms (except the ‘marginal’ ones) leverage an
auxiliary-variable perspective recently proposed in Corenflos and Särkkä,
(2023). We name our algorithms Particle-aMALA, Particle-MALA, Particle-aGRAD,
Particle-mGRAD and Particle-PCNL. This is motivated by the fact that if
$T=N=1$ (where $N\in\mathbb{N}$ is the number of particles), they reduce to
the ‘classical’ MCMC algorithms: aMALA, MALA, PCNL, aGRAD and mGRAD,
respectively. Furthermore, if $T=1$ but $N>1$, our methods constitute novel
multi-proposal versions of such ‘classical’ MCMC algorithms which may
themselves be of interest with a view to exploiting parallelisation.

The generalisation of such ‘classical’ MCMC algorithms to $T>1$ time steps is,
however, not unique. And so we have presented additional variants named
Particle-aMALA+ and Particle-aGRAD+ and _twisted_ Particle-aGRAD/Particle-
aGRAD+. These can be viewed as ‘lookahead’ methods because their proposals
employ ‘smoothing’ rather than ‘filter’ gradients or utilise information
contained in future auxiliary variables. Notably, if $N=1$ but $T>1$, then the
IMH, RWM, aMALA and aGRAD algorithm can still be recovered as a special case
of slightly modified versions of the CSMC, Particle-RWM, Particle-aMALA+, and
twisted Particle-aGRAD+ algorithms (and also of the twisted Particle-aGRAD
algorithm if $G_{t}(\mathbf{x}_{t-1:t})$ is constant in $\mathbf{x}_{t-1}$).
Specifically, this modification would entail that the latter use no resampling
(i.e., they instead set $\smash{a_{t}^{n}=n}$ for all $n\in[N]_{0}$ and all
$t\in[T-1]$), use _ancestral tracing_ instead of backward sampling (i.e., they
instead set $\smash{l_{t}=a_{t}^{l_{t+1}}}$ for all $t\in[T-1]$) and use
$\delta_{1}=\dotsc=\delta_{T}$.

We have further proved that the Particle-aGRAD/Particle-mGRAD algorithms have
the desirable property that they naturally recover (a) the CSMC algorithm if
the prior dynamics are highly informative (i.e., if the target posterior
distribution is dominated by the prior); (b) the Particle-aMALA/Particle-MALA
if the prior dynamics are completely uninformative (i.e., if the target
posterior distribution is dominated by the likelihood). This property
independently helps explain the impressive performance of aGRAD and mGRAD
reported in Titsias and Papaspiliopoulos,  (2018).

Our methods have enabled Bayesian inference in a multivariate stochastic
volatility model with $D=30$ assets and $T=128$ observations ($3840$ unknowns
in total) in which neither CSMC nor aMALA/MALA gave reliable estimates. In
particular, in this application, our twisted Particle-aGRAD algorithm strongly
outperformed the existing sophisticated aGRAD algorithm – even when accounting
for computation time.

###  6.2 Limitations

The main limitations of our methods are the same as in all gradient-based
‘classical’ MCMC algorithms. First, they require continuously differentiable
target densities (more precisely, the densities $Q_{t}(\mathbf{x}_{t-1:t})$
need to be computable and differentiable pointwise). This requirement is
slightly softened for the methods of Section 4 where only the likelihood
$G_{t}(\mathbf{x}_{t-1:t})$ is required to be differentiable, at the cost of
needing (at least conditionally) Gaussian prior dynamics
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})$. The favourable scaling with the
dimension $D$ also typically requires target densities to be sufficiently
smooth (see, e.g., Vogrinc and Kendall, , 2021, for counterexamples). Second,
while it improves mixing properties, locality in MCMC is often detrimental
when exploring multi-modal posteriors. This is inherited by our methods which,
too, explore the space by local moves.

###  6.3 Extensions

Our work opens up multiple avenues for further research.

  * •

The algorithms proposed in this work can be extended to more general graphical
models, i.e., they can be combined with suitably ‘conditional’ versions of the
divide-&-conquer sequential Monte Carlo algorithm from Lindsten et al.,
(2017). For instance, for a particular graphical model, such a ‘conditional’
scheme was recently described in Corenflos et al.,  (2022, Section 3).

  * •

Particle-aMALA can be incorporated straightforwardly into the methodology from
Corenflos et al.,  (2022) to reduce the computation time per MCMC update from
$\operatorname{\mathrm{O}}(T)$ to $\operatorname{\mathrm{O}}(\log_{2}T)$ (for
some fixed dimension $D$) on parallel architectures. While less directly
obvious (because of the non-Markovianity of the auxiliary target), the
smoothing-gradient version Particle-aMALA+ is likely parallelisable, too, by
simply extending the framework to compute weight functions over three time
steps rather than two. It is however less clear that Particle-MALA is
parallelisable, as the marginalisation has to be done across two time steps
rather than one as presented in Section 3.2.

  * •

All our algorithms can be straightforwardly extended to use other resampling
schemes than conditional multinomial resampling, e.g., conditional systematic
resampling. In fact, in our experiments, we used the conditional killing
resampling which is stable under low-informative likelihoods (Karppinen et
al., , 2023), a regime that may happen in our case when $\delta_{t}$ takes
very small values at calibration time.

  * •

In this work, we have left aside the question of choosing $\delta_{t}$ and
have elected to take it to correspond to a $75\text{\,}\mathrm{\char
37\relax}\text{/}$ acceptance rate throughout. It is however clear that its
optimal value (and the optimal value of the acceptance rate) depends on the
number of proposals $N$ and on the dimension $D$. An optimal-scaling analysis
(see Roberts and Rosenthal, , 2001, and references therein) of the methods
proposed in this work is therefore needed. An optimal-scaling analysis for a
related algorithm without backward sampling and without gradient or prior-
informed proposals can be found in Malory,  (2021).

  * •

In Section 4 (in which we propose Particle-aGRAD and variations thereof), we
have assumed that the covariance matrices $\mathbf{C}_{t}(\mathbf{x}_{t-1})$
(or $\mathbf{C}_{t})$ are non-singular. However, it is worth noting that
proposal kernels used by the methods from Section 4 remain valid if the
covariance matrices are singular, in the sense that they are _still_
absolutely continuous w.r.t. the true dynamics. However, the use of backward
sampling is no longer possible for such degenerate dynamics. Instead, one must
resort to ancestral tracing, i.e., taking $l_{t}\coloneqq a_{t}^{l_{t+1}}$,
for $t=T-1,\dotsc,1$. However, in this case, $N$ needs to grow with $T$ at a
suitable rate which depends on the stability properties of the model, but at
least linearly (Andrieu et al., , 2018; Lindsten et al., , 2015). An
alternative is to fix $N$ but decrease the step sizes $\delta_{t}$ with $t$
(which would automatically occur when using adaptation based on acceptance
rates as considered in work), as considered in Malory,  (2021) for a related
method.

  * •

Our proposed algorithms consider solely first-order gradient information. A
natural extension would therefore be to incorporate second-order expansions or
preconditioned and adaptive versions of the Particle-MALA variants. Another
obvious direction of study is to extend our methodology to other MCMC kernels,
such as the recently proposed Barker’s robust proposal (Livingstone and
Zanella, , 2022), or non-reversible discrete-time kernels such as the discrete
bouncy particle sampler Sherlock and Thiery,  (2022). Other natural extensions
would consist of adapting the methodology to non-continuous spaces, e.g.,
using methods from Zanella,  (2020); Rhodes and Gutmann,  (2022), or
constrained spaces.

### Author contributions

A.C. and A.F. jointly developed the methodology, writing was primarily done by
A.F., A.C. implemented and conducted the experiments, after which both A.C.
and A.F. edited and reviewed the final manuscript.

## References

  * Andrieu et al., (2010) Andrieu, C., Doucet, A., and Holenstein, R. (2010).  Particle Markov chain Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342.  With discussion.
  * Andrieu et al., (2018) Andrieu, C., Lee, A., and Vihola, M. (2018).  Uniform ergodicity of the iterated conditional SMC and geometric ergodicity of particle Gibbs samplers.  Bernoulli, 24(2):842–872.
  * Andrieu and Vihola, (2016) Andrieu, C. and Vihola, M. (2016).  Establishing some order amongst exact approximations of MCMCs.  Annals of Applied Probability, 26(5):2661–2696.
  * Besag, (1994) Besag, J. E. (1994).  Contribution to the discussion on ‘Representations of knowledge in complex systems’ by Grenander, U and Miller, M. I..  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 56(4):549–581.
  * Ceperley and Dewing, (1999) Ceperley, D. M. and Dewing, M. (1999).  The penalty method for random walks with uncertain energies.  The Journal of Chemical Physics, 110(20):9812–9820.
  * Chopin and Singh, (2013) Chopin, N. and Singh, S. S. (2013).  On particle Gibbs sampling.  arXiv e-prints, arXiv:1304.1887v1.
  * Corenflos et al., (2022) Corenflos, A., Chopin, N., and Särkkä, S. (2022).  De-sequentialized Monte Carlo: A parallel-in-time particle smoother.  Journal of Machine Learning Research, 23(283):1–39.
  * Corenflos and Särkkä, (2023) Corenflos, A. and Särkkä, S. (2023).  Auxiliary MCMC and particle Gibbs samplers for parallelisable inference in latent dynamical systems.  arXiv preprint arXiv:2303.00301.
  * Cotter et al., (2013) Cotter, S. L., Roberts, G. O., Stuart, A. M., and White, D. (2013).  MCMC methods for functions: Modifying old algorithms to make them faster.  Statistical Science, 28(3):424–446.
  * Fearnhead and Meligkotsidou, (2016) Fearnhead, P. and Meligkotsidou, L. (2016).  Augmentation schemes for particle MCMC.  Statistics and Computing, 26:1293–1306.
  * Finke, (2015) Finke, A. (2015).  On Extended State-Space Constructions for Monte Carlo Methods.  PhD thesis, Department of Statistics, University of Warwick, UK.
  * Finke et al., (2016) Finke, A., Doucet, A., and Johansen, A. M. (2016).  On embedded hidden Markov models and particle Markov chain Monte Carlo methods.  arXiv e-prints, arXiv:1610.08962.
  * Finke and Thiery, (2023) Finke, A. and Thiery, A. H. (2023).  Conditional sequential Monte Carlo in high dimensions.  The Annals of Statistics, 51(2):437–463.
  * Guarniero et al., (2017) Guarniero, P., Johansen, A. M., and Lee, A. (2017).  The iterated auxiliary particle filter.  Journal of the American Statistical Association, 112(520):1636–1647.
  * Hastings, (1970) Hastings, W. K. (1970).  Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, 57(1):97–109.
  * Henderson and Searle, (1981) Henderson, H. V. and Searle, S. R. (1981).  On deriving the inverse of a sum of matrices.  SIAM Review, 23(1):53–60.
  * Heng et al., (2020) Heng, J., Bishop, A. N., Deligiannidis, G., and Doucet, A. (2020).  Controlled sequential Monte Carlo.  The Annals of Statistics, 48(5):2904 – 2929.
  * Kalman, (1960) Kalman, R. E. (1960).  A new approach to linear filtering and prediction problems.  Journal of Basic Engineering, 82:35–45.
  * Karjalainen et al., (2023) Karjalainen, J., Lee, A., Singh, S. S., and Vihola, M. (2023).  Mixing time of the conditional backward sampling particle filter.  arXiv e-prints, arXiv:2312.17572.
  * Karppinen et al., (2023) Karppinen, S., Singh, S. S., and Vihola, M. (2023).  Conditional particle filters with bridge backward sampling.  Journal of Computational and Graphical Statistics, 0(0):1–15.
  * Karppinen and Vihola, (2021) Karppinen, S. and Vihola, M. (2021).  Conditional particle filters with diffuse initial distributions.  Statistics and Computing, 31:1–14.
  * Lee et al., (2020) Lee, A., Singh, S. S., and Vihola, M. (2020).  Coupled conditional backward sampling particle filter.  Annals of Statistics, 48(5):3066–3089.
  * Lindsten et al., (2015) Lindsten, F., Douc, R., and Moulines, E. (2015).  Uniform ergodicity of the particle Gibbs sampler.  Scandinavian Journal of Statistics, 42(3):775–797.
  * Lindsten et al., (2017) Lindsten, F., Johansen, A. M., Naesseth, C. A., Kirkpatrick, B., Schön, T. B., Aston, J. A., and Bouchard-Côté, A. (2017).  Divide-and-conquer with sequential Monte Carlo.  Journal of Computational and Graphical Statistics, 26(2):445–458.
  * Lindsten et al., (2012) Lindsten, F., Jordan, M. I., and Schön, T. B. (2012).  Ancestor sampling for particle Gibbs.  In Proceedings of the 2012 Conference on Neural Information Processing Systems, Lake Tahoe, NV.
  * Liu, (1996) Liu, J. S. (1996).  Peskun’s theorem and a modified discrete-state Gibbs sampler.  Biometrika, 83(3):681–682.
  * Livingstone and Zanella, (2022) Livingstone, S. and Zanella, G. (2022).  The Barker proposal: Combining robustness and efficiency in gradient-based MCMC.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 84(2):496–523.
  * Malory, (2021) Malory, S. (2021).  Bayesian inference for stochastic processes.  PhD thesis, Lancaster University.
  * Metropolis et al., (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. (1953).  Equation of state calculations by fast computing machines.  Journal of Chemical Physics, 21(6):1087–1092.
  * Murray et al., (2013) Murray, L. M., Jones, E. M., and Parslow, J. (2013).  On disturbance state-space models and the particle marginal Metropolis–Hastings sampler.  SIAM/ASA Journal on Uncertainty Quantification, 1(1):494–521.
  * Nicholls et al., (2012) Nicholls, G. K., Fox, C., and Muir Watt, A. (2012).  Coupled MCMC with a randomized acceptance probability.  arXiv e-prints, arXiv:1205.6857.
  * Rhodes and Gutmann, (2022) Rhodes, B. and Gutmann, M. (2022).  Enhanced gradient-based MCMC in discrete spaces.  arXiv e-prints, arXiv:2208.00040.
  * Roberts et al., (1997) Roberts, G. O., Gelman, A., and Gilks, W. R. (1997).  Weak convergence and optimal scaling of random walk Metropolis algorithms.  The Annals of Applied Probability, 7(1):110–120.
  * Roberts and Rosenthal, (1998) Roberts, G. O. and Rosenthal, J. S. (1998).  Optimal scaling of discrete approximations to Langevin diffusions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268.
  * Roberts and Rosenthal, (2001) Roberts, G. O. and Rosenthal, J. S. (2001).  Optimal scaling for various Metropolis–Hastings algorithms.  Statistical Science, 16(4):351–367.
  * Särkkä and Svensson, (2023) Särkkä, S. and Svensson, L. (2023).  Bayesian filtering and smoothing, volume 17.  Cambridge University Press.
  * Sherlock and Thiery, (2022) Sherlock, C. and Thiery, A. H. (2022).  A discrete bouncy particle sampler.  Biometrika, 109(2):335–349.
  * Shestopaloff and Neal, (2018) Shestopaloff, A. Y. and Neal, R. M. (2018).  Sampling latent states for high-dimensional non-linear state space models with the embedded HMM method.  Bayesian Analysis, 13(3):797–822.
  * Singh et al., (2017) Singh, S. S., Lindsten, F., and Moulines, E. (2017).  Blocking strategies and stability of particle Gibbs samplers.  Biometrika, 104(4):953–969.
  * Titsias, (2011) Titsias, M. K. (2011).  Contribution to the discussion on ‘Riemann manifold Langevin and Hamiltonian Monte Carlo methods’ by Girolami, M., and Calderhead, b.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 73(2):123–214.
  * Titsias and Papaspiliopoulos, (2018) Titsias, M. K. and Papaspiliopoulos, O. (2018).  Auxiliary gradient-based sampling algorithms.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(4):749–767.
  * Tjelmeland, (2004) Tjelmeland, H. (2004).  Using all Metropolis–Hastings proposals to estimate mean values.  preprint 4/2004, Norwegian University of Science and Technology, Trondheim, Norway.
  * Vehtari et al., (2021) Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2021).  Rank-normalization, folding, and localization: An improved $\smash{\widehat{R}}$ for assessing convergence of MCMC (with discussion).  Bayesian Analysis, 16(2):667–718.
  * Vogrinc and Kendall, (2021) Vogrinc, J. and Kendall, W. S. (2021).  Counterexamples for optimal scaling of Metropolis–Hastings chains with rough target densities.  The Annals of Applied Probability, 31(2):972–1019.
  * Whiteley, (2010) Whiteley, N. (2010).  Contribution to the discussion on ‘Particle Markov chain Monte Carlo methods’ by Andrieu, C., Doucet, A., and Holenstein, R.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):306–307.
  * Whiteley and Lee, (2014) Whiteley, N. and Lee, A. (2014).  Twisted particle filters.  The Annals of Statistics, 42(1):115–141.
  * Zanella, (2020) Zanella, G. (2020).  Informed proposals for local MCMC in discrete spaces.  Journal of the American Statistical Association, 115(530):852–865.

##  Appendix A Particle extensions of PCN(L)

###  A.1 Particle-aPCNL

In this section, we extend the _preconditioned Crank–Nicolson–Langevin (PCNL)_
algorithm (and also the _preconditioned Crank–Nicolson (PCN)_ algorithm
recovered by setting $\kappa=0$) (Cotter et al., , 2013) to $T>1$ time steps
and $N>1$. As a by-product, we derive an ‘auxiliary-variable’ version of PCNL
which was mentioned, but not explicitly stated, in Titsias and
Papaspiliopoulos,  (2018). Throughout this section, we assume the prior
dynamics are conditionally Gaussian,
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$
as in (25).

Throughout this section, we use the parametrisation of the PCNL algorithm from
Titsias and Papaspiliopoulos,  (2018)333In the parametrisation from Cotter et
al.,  (2013), we would have $\delta_{t}\in[0,2]$
$\beta_{t}=\frac{2-\delta_{t}}{2+\delta_{t}}$, i.e., we set

$\displaystyle\beta_{t}\coloneqq\frac{2}{2+\delta_{t}}\in(0,1),$ (48)

where $\delta_{t}>0$ is again the step size at time $t$. Note that this
implies that $\tfrac{1-\beta_{t}}{\beta_{t}}=\tfrac{\delta_{t}}{2}$.

The first method proposed in this section is termed _Particle-aPCNL_.
Conditional on the auxiliary variables $\mathbf{u}_{1:T}$, it can be viewed as
a CSMC algorithm whose proposal kernels are those of the fully-adapted
auxiliary particle filter for the state-space model defined by the Gaussian
transitions
$p(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}_{t}(\mathbf{x}_{t-1}),\mathbf{C}_{t}(\mathbf{x}_{t-1}))$
from (25) and ‘pseudo observations’ $\mathbf{u}_{t}$ with
$p(\mathbf{u}_{t}|\mathbf{x}_{t})=\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\frac{\delta_{t}}{2}\mathbf{C}_{t}(\mathbf{x}_{t-1}))$.
We now write

$\displaystyle M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})$
$\displaystyle\coloneqq p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})$
(49) $\displaystyle\propto
M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\tfrac{\delta_{t}}{2}\mathbf{C}_{t}(\mathbf{x}_{t-1}))$
(50)
$\displaystyle\propto\operatorname{N}(\mathbf{x}_{t};\mathbf{m}^{\prime}_{t}(\mathbf{x}_{t-1},\mathbf{u}_{t}),\mathbf{C}^{\prime}_{t}(\mathbf{x}_{t-1})),$
(51) with
$\displaystyle\mathbf{m}^{\prime}_{t}(\mathbf{x}_{t-1},\mathbf{u}_{t})$
$\displaystyle\coloneqq\beta_{t}\mathbf{u}_{t}+(1-\beta_{t})\mathbf{m}_{t}(\mathbf{x}_{t-1}),$
(52) $\displaystyle\mathbf{C}^{\prime}_{t}(\mathbf{x}_{t-1})$
$\displaystyle\coloneqq(1-\beta_{t})\mathbf{C}_{t}(\mathbf{x}_{t-1}),$ (53)

as well as

$\displaystyle G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$
$\displaystyle\coloneqq
Q_{t}(\mathbf{x}_{t-1:t})\frac{\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\tfrac{\delta_{t}}{2}\mathbf{C}_{t}(\mathbf{x}_{t-1}))}{M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})},$
(54)

and $Q^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t})$.
Here,
$\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\in\mathbb{R}^{D\times
D}$ is some preconditioning matrix whose choice is discussed in Section A.5
below.

A single iteration of the Particle-aPCNL algorithm is then as follows.

[Particle-aPCNL ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c) and the weight calculation (Step 1d) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\tfrac{\delta_{t}}{2}\mathbf{C}_{t}(\mathbf{x}_{t-1}))}$,
and $\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{t})}$, for
$n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

for $n\in[N]_{0}$, set $\smash{w_{t}^{n}\propto
G^{\prime}_{t}(\mathbf{x}_{t-1:t}^{(n)};\mathbf{u}_{t}),}$

and also replace $Q_{t+1}(\,\cdot\,)$ in the backward kernel in Step 3 by
$\smash{Q^{\prime}_{t+1}(\,\cdot\,;\mathbf{u}_{t})}$.

######  Proposition 10 (validity of Particle-aPCNL).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm A.1
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAPCNL}{\leavevmode
Particle-aPCNL}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves
$\pi_{T}$ invariant.

###  A.2 Particle-PCNL

In this section, in analogy to the Particle-MALA and Particle-mGRAD algorithms
from the main manuscript, we analytically integrate out the auxiliary
variables $\mathbf{u}_{t}$ which appeared in the weights of the Particle-aPCNL
algorithm. As in the case of the Particle-mGRAD algorithm, we assume that the
covariance matrices appearing in the conditionally Gaussian mutation kernel
(25) do not depend on the previous state, i.e.,
$\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t}$ (33).

A single iteration of the resulting methodology – which we term the _Particle-
PCNL_ algorithm – is as follows, where we write

$\displaystyle\log
H_{t,\bm{\phi}}(\mathbf{x},\mathbf{v},\bar{\mathbf{x}},\bar{\mathbf{v}})$
$\displaystyle=\tfrac{1}{2}(\beta_{t}^{-1}+N+1)(\mathbf{x}-\mathbf{v})^{\mathrm{T}}\mathbf{G}_{t}(\mathbf{x}-\mathbf{v})$
$\displaystyle\qquad-\tfrac{1}{2}N\beta_{t}(\mathbf{x}+\bm{\phi})^{\mathrm{T}}\mathbf{G}_{t}(\mathbf{x}+\bm{\phi})$
$\displaystyle\qquad+(N+1)(\bar{\mathbf{x}}-\bar{\mathbf{v}})^{\mathrm{T}}\mathbf{G}_{t}(\mathbf{v}+\bm{\phi})$
$\displaystyle\qquad-(\mathbf{x}-\mathbf{v})^{\mathrm{T}}\mathbf{G}_{t}(\mathbf{x}+\bm{\phi}),$

for

$\displaystyle\mathbf{G}_{t}\coloneqq\frac{\beta_{t}}{(1-\beta_{t})(1+N\beta_{t})}\mathbf{C}_{t}^{-1}=\frac{2(\delta_{t}+2)}{\delta_{t}(\delta_{t}+2+N)}\mathbf{C}_{t}^{-1}.$
(55)

[Particle-PCNL ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c) and the weight calculation (Step 1d) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\tfrac{\delta_{t}}{2}\mathbf{C}_{t})}$, and
$\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{t})}$, for
$n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

set
$\smash{\bar{\mathbf{x}}_{t}\coloneqq\tfrac{1}{N+1}\sum_{n=0}^{N}\mathbf{x}_{t}^{n}}$,
$\mathbf{v}_{t}^{n}\coloneqq(1-\beta_{t})\mathbf{m}_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{n}})$,
$\smash{\bar{\mathbf{v}}_{t}\coloneqq\tfrac{1}{N+1}\sum_{n=0}^{N}\mathbf{v}_{t}^{n}}$,
and, for $n\in[N]_{0}$,

$\displaystyle w_{t}^{n}\propto
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})H_{t,\kappa\frac{\delta_{t}}{2}\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{n}})\nabla_{\mathbf{x}_{t}^{n}}\log
G_{t}(\mathbf{x}_{t-1:t}^{(n)})}(\mathbf{x}_{t}^{n},\mathbf{v}_{t}^{n},\bar{\mathbf{x}}_{t},\bar{\mathbf{v}}_{t}).$
(56)

In the same way as outlined in Remarks 1 and 2, the Particle-aPCNL algorithm
can be viewed as an ‘exact approximation’ of the Particle-PCNL algorithm.

######  Proposition 11 (validity of Particle-PCNL).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm A.2
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEPCNL}{\leavevmode Particle-
PCNL}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves $\pi_{T}$
invariant.

###  A.3 Particle-aPCNL+

In this section, similar to the Particle-aMALA+ and Particle-aGRAD+ algorithms
from the main manuscript, we extend the Particle-aPCNL algorithm to
incorporate gradients w.r.t. the ‘smoothing’ potential
$G_{1:T}(\mathbf{x}_{1:T})=\prod_{t=1}^{T}G_{t}(\mathbf{x}_{t-1:t})$ rather
than w.r.t. the ‘filtering’ potential
$\prod_{s=1}^{t}G_{s}(\mathbf{x}_{s-1:s})$.

For $M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})$ and
$G^{\prime}_{t}(\mathbf{x}_{t-1:t},\mathbf{u}_{t})$ still defined as in the
Particle-aPCNL algorithm (i.e., as in (51) and (54)), we now write

$\displaystyle G^{\prime}_{t}(\mathbf{x}_{t-2:t},\mathbf{u}_{1:T})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\coloneqq
G^{\prime}_{t}(\mathbf{x}_{t-1:t},\mathbf{u}_{t})\frac{\operatorname{N}(\mathbf{u}_{t-1};\mathbf{x}_{t-1}+\kappa\tfrac{\delta_{t-1}}{2}\smash{\widetilde{\mathbf{C}}}_{t-1}(\mathbf{x}_{t-2})\nabla_{\mathbf{x}_{t-1}}\log
G_{1:T}(\mathbf{x}_{1:T}),\tfrac{\delta_{t-1}}{2}\mathbf{C}_{t-1}(\mathbf{x}_{t-2}))}{\operatorname{N}(\mathbf{u}_{t-1};\mathbf{x}_{t-1}+\kappa\tfrac{\delta_{t-1}}{2}\smash{\widetilde{\mathbf{C}}}_{t-1}(\mathbf{x}_{t-2})\nabla_{\mathbf{x}_{t-1}}\log
G_{t-1}(\mathbf{x}_{t-2:t-1}),\tfrac{\delta_{t-1}}{2}\mathbf{C}_{t-1}(\mathbf{x}_{t-2}))},$

as well as $Q^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})$,
where we note that

$\displaystyle\nabla_{\mathbf{x}_{t}}\log
G_{1:T}(\mathbf{x}_{1:T})=\nabla_{\mathbf{x}_{t}}[\log
G_{t}(\mathbf{x}_{t-1:t})+\log G_{t+1}(\mathbf{x}_{t:t+1})].$ (57)

A single iteration of the resulting methodology – which we term the _Particle-
aPCNL+_ algorithm – is as follows.

[Particle-aPCNL+ ] Implement Algorithm 2.1.1 but replace the particle proposal
(Step 1c), the weight calculation (Step 1d), and backward sampling (Step 3) by

  1. 1c.

sample
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\nabla_{\mathbf{x}_{t}}\log
G_{T}(\mathbf{x}_{1:T}),\tfrac{\delta_{t}}{2}\mathbf{C}_{t}(\mathbf{x}_{t-1}))}$,
and $\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{t})}$, for
$n\in[N]_{0}\setminus\\{k_{t}\\}$,

  2. 1d.

for $n\in[N]_{0}$, set $\smash{w_{t}^{n}\propto
G^{\prime}_{t}(\mathbf{x}_{t-2:t}^{(n)};\mathbf{u}_{t-1:t}),}$

  3. 3.

for $t=T-1,\dotsc,1$, sample $l_{t}=i\in[N]_{0}$ w.p.

$\displaystyle\dfrac{W_{t}^{i}Q^{\prime}_{t+1}((\mathbf{x}_{t-1:t}^{(i)},\mathbf{x}_{t+1}^{l_{t+1}});\mathbf{u}_{t:t+1})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{i},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{t+1:t+2})}{\sum_{n=0}^{N}W_{t}^{n}Q^{\prime}_{t+1}((\mathbf{x}_{t-1:t}^{(n)},\mathbf{x}_{t+1}^{l_{t+1}});\mathbf{u}_{t:t+1})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{n},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{t+1:t+2})}.$
(58)

Note that if $G_{t}(\mathbf{x}_{t-1:t})=G_{t}(\mathbf{x}_{t})$ does not depend
on $\mathbf{x}_{t-1}$, then the Particle-aPCNL+ algorithm coincides with the
Particle-aPCNL algorithm. However, when $G_{t}(\mathbf{x}_{t-1:t})$ varies
highly in $\mathbf{x}_{t-1}$, their behaviours may differ substantially.

######  Proposition 12 (validity of Particle-aPCNL+).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm A.3
induces a Markov kernel
$P_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAPCNLPLUS}{\leavevmode
Particle-aPCNL+}}}(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which leaves
$\pi_{T}$ invariant.

###  A.4 Twisted Particle-aPCNL(+)

In analogue to the twisted Particle-aGRAD and twisted Particle-aGRAD+
algorithms, we can again construct ‘twisted’ versions of the Particle-aPCNL
and Particle-aPCNL+ algorithms, under the assumption that
$M_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}_{t}\mathbf{x}_{t-1}+\mathbf{b}_{t},\mathbf{C}_{t})$,
i.e., (38).

We start with the twisted Particle-aPCNL algorithm. We now write

$\displaystyle
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})$
$\displaystyle\coloneqq p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t:T})$
(59)
$\displaystyle\propto\int_{\mathcal{X}^{T-t}}\biggl{[}\prod_{s=t}^{T}\operatorname{N}(\mathbf{x}_{s};\mathbf{F}_{s}\mathbf{x}_{s-1}+\mathbf{b}_{s},\mathbf{C}_{s})\operatorname{N}(\mathbf{u}_{s};\mathbf{x}_{s},\tfrac{\delta_{s}}{2}\mathbf{C}_{t})\biggr{]}\,\mathrm{d}\mathbf{x}_{t+1:T}$
(60)
$\displaystyle\propto\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t}),$
(61) $\displaystyle G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t:T})$
$\displaystyle\coloneqq
Q_{t}(\mathbf{x}_{t-1:t})\frac{\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t}+\kappa\tfrac{\delta_{t}}{2}\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t}),\tfrac{\delta_{t}}{2}\mathbf{C}_{t})}{M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t:T})},$
(62)

as well as $Q^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t:T})\coloneqq
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})G^{\prime}_{t}(\mathbf{x}_{t-1:t};\mathbf{u}_{t:T})$.
Here, $\mathbf{b}^{\prime}_{t}\in\mathbb{R}^{D}$ and
$\mathbf{F}^{\prime}_{t},\mathbf{C}^{\prime}_{t}\in\mathbb{R}^{D\times D}$ can
again be obtained via the Kalman-filtering recursions given in Appendix B.

A single iteration of the resulting methodology – which we term the _twisted
Particle-aPCNL_ algorithm – is then exactly as the Particle-aPCNL algorithm
(Algorithm A.1), except that
$M^{\prime}_{t}(\,\cdot\,|\,\cdot\,;\mathbf{u}_{t})$,
$G^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t})$ and
$Q^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t})$ from Section A.1 are replaced by
$M^{\prime}_{t}(\,\cdot\,|\,\cdot\,;\mathbf{u}_{t:T})$,
$G^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t:T})$ and
$Q^{\prime}_{t}(\,\cdot\,;\mathbf{u}_{t:T})$ from this section. When the
potential functions $G_{t}(\mathbf{x}_{t-1:t})$ vary in $\mathbf{x}_{t-1}$,
then we can further construct a _twisted Particle-aGRAD+_ algorithm by
replacing $M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t})$ in
Algorithm A.3 and in the denominator of
$G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{t-1:t})$ by
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{t:T})$.

######  Proposition 13 (validity of the twisted Particle-aPCNL/Particle-
aPCNL+).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via the twisted
Particle-aPCNL or twisted Particle-aPCNL+ algorithm induces a Markov kernel
which leaves $\pi_{T}$ invariant.

###  A.5 Choice of preconditioning matrix

There is some degree of freedom in choosing the preconditioning matrices
$\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\in\mathbb{R}^{D\times
D}$ in the algorithms presented above.

  1. 1.

A simple option which does not require further assumptions is to take

$\displaystyle\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})\coloneqq\mathbf{C}_{t}(\mathbf{x}_{t-1}).$
(63)

  2. 2.

If we make the stronger model assumption that
$\mathbf{m}_{t}(\mathbf{x}_{t-1})=\mathbf{F}_{t}\mathbf{x}_{t-1}+\mathbf{b}_{t}$
and $\mathbf{C}_{t}(\mathbf{x}_{t-1})=\mathbf{C}_{t}$ (38) (which is assumed
to hold for the twisted versions of the Particle-aPCNL and Particle-aPCNL+
algorithms anyway), then we could alternatively set

$\displaystyle\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t-1})=\smash{\widetilde{\mathbf{C}}}_{t}\coloneqq\sum_{s=1}^{T}\bm{\Sigma}_{s,t},$
(64)

where $\bm{\Sigma}_{s,t}\in\mathbb{R}^{D\times D}$ is the block $(s,t)$ in the
covariance matrix $\bm{\Sigma}\in\mathbb{R}^{TD\times TD}$ of the prior
dynamics $M_{1:T}(\mathbf{x}_{1:T})$, i.e.,

$\displaystyle\bm{\Sigma}_{s,t}=\begin{cases}\mathbf{F}_{s}\bm{\Sigma}_{s-1,t},&\text{if
$s>t$,}\\\ \bm{\Sigma}_{t,s}^{\mathrm{T}},&\text{if $s<t$,}\\\
\bm{\Sigma}_{t},&\text{if $s=t$,}\end{cases}$ (65)

where $\bm{\Sigma}_{t}$ can be found via the recursion from Step 1 of
Algorithm B from Appendix B. As discussed below, this specification has the
potentially useful implication that the algorithm reduces to an ‘auxiliary-
variable’ version of the PCNL algorithm on the path space in the absence of
resampling and backward sampling. Unfortunately, evaluating the
preconditioning matrices $\smash{\widetilde{\mathbf{C}}}_{t}$ is likely to
incur a quadratic computational complexity in $T$ which we prefer to avoid.

  3. 3.

A compromise (which retains linear computational complexity in $T$) may be to
truncate the above sum by setting