$\displaystyle\smash{\widetilde{\mathbf{C}}}_{t}\coloneqq\sum_{s=(t-L)\vee
1}^{(t+L)\wedge T}\bm{\Sigma}_{s,t},$ (66)

for some $L\in[T]_{0}$ (note that this still requires the model assumption
(38)).

###  A.6 Relationship with other methods

The algorithms proposed above relate to existing methods as follows.

  1. 1.

Generalisation of aPCNL. For $\kappa=1$, the Particle-aPCNL algorithm (and
similarly the Particle-aPCNL+ algorithm as well as the twisted versions of
either) generalises an _auxiliary preconditioned Crank–Nicolson–Langevin
(aPCNL)_ algorithm (which was mentioned but not explicitly derived in 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 so that $\pi(\mathbf{x})\propto
M(\mathbf{x})G(\mathbf{x})$, where
$M(\mathbf{x})=\operatorname{N}(\mathbf{x};\mathbf{m},\mathbf{C})$. We also
take $\smash{\widetilde{\mathbf{C}}}\coloneqq\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 A.1 first
refreshes the auxiliary variable by sampling
$\mathbf{u}\sim\operatorname{N}(\mathbf{x}^{0}+\tfrac{\delta}{2}\mathbf{C}\nabla\log
G(\mathbf{x}^{0}),\tfrac{\delta}{2}\mathbf{C})$ and then proposes
$\mathbf{x}^{1}\sim\operatorname{N}((1-\beta)\mathbf{m}+\beta\mathbf{u},(1-\beta)\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}{APCNL}{\leavevmode
aPCNL}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})$, where

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{APCNL}{\leavevmode
aPCNL}}}(\mathbf{x}^{0},\mathbf{x}^{1};\mathbf{u})$ (67)
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}$ (68)
$\displaystyle=\frac{\pi(\mathbf{x}^{1})\operatorname{N}(\mathbf{u};\mathbf{x}^{1}+\tfrac{\delta}{2}\mathbf{C}\nabla\log
G(\mathbf{x}^{1}),\tfrac{\delta}{2}\mathbf{C})\operatorname{N}(\mathbf{x}^{0};(1-\beta)\mathbf{m}+\beta\mathbf{u},(1-\beta)\mathbf{C})}{\pi(\mathbf{x}^{0})\operatorname{N}(\mathbf{u};\mathbf{x}^{0}+\tfrac{\delta}{2}\mathbf{C}\nabla\log
G(\mathbf{x}^{0}),\tfrac{\delta}{2}\mathbf{C})\operatorname{N}(\mathbf{x}^{1};(1-\beta)\mathbf{m}+\beta\mathbf{u},(1-\beta)\mathbf{C})}.$
(69)

Otherwise, the old state
$\tilde{\mathbf{x}}\coloneqq\mathbf{x}^{0}=\mathbf{x}$ is returned as the new
state. If $N=1$ and $\kappa=1$ but $T>1$ then the aPCNL algorithm could still
be recovered as a special case of (a slightly modified version of) the twisted
Particle-aPCNL+ algorithm (and also of the twisted Particle-aPCNL algorithm if
the potential functions $G_{t}(\mathbf{x}_{t-1:t})$ do not depend on
$\mathbf{x}_{t-1}$) if the latter uses no resampling and ancestral tracing
instead of backward sampling, and if $\delta_{1}=\dotsc=\delta_{T}$ and if the
preconditioning matrices are specified via (64). However, we do not recommend
this choice of preconditioning matrix as it leads to squared computational
complexity in $T$ (also incurred by aPCNL).

  2. 2.

Generalisation of PCNL. Still taking $\kappa=1$, the Particle-PCNL algorithm
generalises the _preconditioned Crank–Nicolson–Langevin (PCNL)_ algorithm
(Cotter et al., , 2013) 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 aPCNL above. Step 1c. of Algorithm A.2 then
marginally proposes
$\mathbf{x}^{1}\sim\operatorname{N}((1-\beta)\mathbf{m}+\beta[\mathbf{x}^{0}+\tfrac{\delta}{2}\mathbf{C}\nabla\log
G(\mathbf{x}^{0})],(1-\beta^{2})\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}{PCNL}{\leavevmode
PCNL}}}(\mathbf{x}^{0},\mathbf{x}^{1})$, where

$\displaystyle\alpha_{\text{\lx@glossaries@gls@link{acronym}{PCNL}{\leavevmode
PCNL}}}(\mathbf{x}^{0},\mathbf{x}^{1})$
$\displaystyle\coloneqq\frac{1-W^{0}}{1-W^{1}}$ (70)
$\displaystyle=\frac{\pi(\mathbf{x}^{1})\operatorname{N}((1-\beta)\mathbf{m}+\beta[\mathbf{x}^{1}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{1})],(1-\beta^{2})\mathbf{C})}{\pi(\mathbf{x}^{0})\operatorname{N}((1-\beta)\mathbf{m}+\beta[\mathbf{x}^{0}+\tfrac{\delta}{2}\nabla\log
G(\mathbf{x}^{0})],(1-\beta^{2})\mathbf{C})}.$ (71)

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

##  Appendix B Twisted proposals

In this section, we detail the mutation kernel
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})$ used by the
twisted Particle-aGRAD/Particle-aGRAD+ (41) and twisted Particle-
aPCNL/Particle-aPCNL+ (61) algorithms. This mutation kernel can be thought of
as the fully-twisted particle-filter proposal for the state-space model which
is defined by the Gaussian transitions
$p(\mathbf{x}_{t}|\mathbf{x}_{t-1})=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})$
from (38) and observation densities
$p(\mathbf{u}_{t}|\mathbf{x}_{t})=\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\frac{\delta_{t}}{2}\mathbf{V}_{t})$,
where we take $\mathbf{V}_{t}\coloneqq\mathbf{I}$ in the case of the twisted
Particle-aGRAD or twisted Particle-aGRAD+ algorithm and
$\mathbf{V}_{t}\coloneqq\mathbf{C}_{t}$ in the case of the twisted Particle-
aPCNL or twisted Particle-aPCNL+ algorithm:

$\displaystyle
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})$
$\displaystyle=p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{1:T})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t}).$
(72)

##### General algorithm.

Algorithm B explains how the twisted-proposal parameters
$\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 calculated at linear complexity in $T$, independently of the total number
of particles, $N$. Notably, Algorithm B does not require $\mathbf{C}_{t}$ to
be invertible.

[twisted-proposal parameters] At the start of an iteration of the twisted
Particle-aGRAD, Particle-aGRAD+, Particle-aPCNL or Particle-aPCNL+ algorithm
(after having sampled all the auxiliary variables $\mathbf{u}_{1:T}$ upfront –
e.g., as in Algorithm D.1 from Appendix D.1 – which is possible because these
only depend on the reference path),

  1. 1.

recursively compute the moments of
$p(\mathbf{x}_{t})=\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t},\bm{\Sigma}_{t})$,
for $t=1,\dotsc,T$, as

$\displaystyle\bm{\mu}_{t}$
$\displaystyle\coloneqq\mathbf{F}_{t}\bm{\mu}_{t-1}+\mathbf{b}_{t},$ (73)
$\displaystyle\bm{\Sigma}_{t}$
$\displaystyle\coloneqq\mathbf{F}_{t}\bm{\Sigma}_{t-1}\mathbf{F}_{t}^{\mathrm{T}}+\mathbf{C}_{t},$
(74)

if $t>1$, and with initial condition $\bm{\mu}_{1}=\mathbf{b}_{1}$ and
$\bm{\Sigma}_{1}=\mathbf{C}_{1}$,

  2. 2.

recursively compute the moments of the time-reversed state transition kernels
$p(\mathbf{x}_{t}|\mathbf{x}_{t+1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}_{t}^{\leftarrow}\mathbf{x}_{t+1}+\mathbf{b}_{t}^{\leftarrow},\mathbf{C}_{t}^{\leftarrow})$,
for $t=T-1,\dotsc,1$, as

$\displaystyle\mathbf{F}_{t}^{\leftarrow}$
$\displaystyle\coloneqq\mathbf{F}_{t+1}\bm{\Sigma}_{t}\bm{\Sigma}_{t+1}^{-1},$
(75) $\displaystyle\mathbf{b}_{t}^{\leftarrow}$
$\displaystyle\coloneqq\bm{\mu}_{t}-\mathbf{F}_{t}^{\leftarrow}\bm{\mu}_{t+1},$
(76) $\displaystyle\mathbf{C}_{t}^{\leftarrow}$
$\displaystyle\coloneqq\bm{\Sigma}_{t}-\mathbf{F}_{t}^{\leftarrow}\bm{\Sigma}_{t}\mathbf{F}_{t+1}^{\mathrm{T}}.$
(77)

  3. 3.

run the Kalman filtering recursion for the time-reversed state-space model
(i.e., with observation densities
$p(\mathbf{u}_{t}|\mathbf{x}_{t})=\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\frac{\delta_{t}}{2}\mathbf{V}_{t})$,
initial distribution $p(\mathbf{x}_{t})$ and time-reversed state transitions
$p(\mathbf{x}_{t}|\mathbf{x}_{t+1})$ found in Steps 1 and 2) to compute the
moments of
$p(\mathbf{x}_{t}|\mathbf{u}_{t:T})=\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t|t}^{\leftarrow},\bm{\Sigma}_{t|t}^{\leftarrow})$,
for $t=T,\dotsc,1$,

  4. 4.

set $\mathbf{F}^{\prime}_{1}\coloneqq\bm{0}_{D\times D}$,
$\mathbf{b}^{\prime}_{1}\coloneqq\bm{\mu}_{1|1}^{\leftarrow}$ as well as
$\mathbf{C}^{\prime}_{1}\coloneqq\bm{\Sigma}_{1|1}^{\leftarrow}$, and, for
$t=2,\dotsc,T$,

$\displaystyle\mathbf{F}^{\prime}_{t}$
$\displaystyle\coloneqq\bm{\Sigma}_{t|t}^{\leftarrow}\mathbf{F}_{t-1}^{\leftarrow}(\mathbf{C}_{t-1}^{\leftarrow}+\mathbf{F}_{t-1}^{\leftarrow}\bm{\Sigma}_{t|t}^{\leftarrow}\\{\mathbf{F}_{t-1}^{\leftarrow}\\}^{\mathrm{T}})^{-1},$
(78) $\displaystyle\mathbf{b}^{\prime}_{t}$
$\displaystyle\coloneqq\bm{\mu}_{t|t}^{\leftarrow}-\mathbf{F}^{\prime}_{t}(\mathbf{F}_{t-1}^{\leftarrow}\bm{\mu}_{t|t}^{\leftarrow}+\mathbf{b}_{t-1}^{\leftarrow}),$
(79) $\displaystyle\mathbf{C}^{\prime}_{t}$
$\displaystyle\coloneqq(\mathbf{I}-\mathbf{F}^{\prime}_{t}\mathbf{F}_{t-1}^{\leftarrow})\bm{\Sigma}_{t|t}^{\leftarrow}.$
(80)

Algorithm B is justified by the decomposition

$\displaystyle
M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t:T})$
$\displaystyle\propto p(\mathbf{x}_{t-1:t}|\mathbf{u}_{t:T})$ (81)
$\displaystyle=p(\mathbf{x}_{t-1}|\mathbf{x}_{t})p(\mathbf{x}_{t}|\mathbf{u}_{t:T})$
(82)
$\displaystyle=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}_{t}^{\leftarrow}\mathbf{x}_{t+1}+\mathbf{b}_{t}^{\leftarrow},\mathbf{C}_{t}^{\leftarrow})\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t|t}^{\leftarrow},\bm{\Sigma}_{t|t}^{\leftarrow})$
(83)
$\displaystyle\propto\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t}),$
(84)

where, as described in Algorithm B,
$p(\mathbf{x}_{t}|\mathbf{u}_{t:T})=\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t|t}^{\leftarrow},\bm{\Sigma}_{t|t}^{\leftarrow})$
is the time-$t$ filter for the time-reversed state-space model with the same
observation densities
$p(\mathbf{u}_{t}|\mathbf{x}_{t})=\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t},\frac{\delta_{t}}{2}\mathbf{V}_{t})$
as before but with initial distribution
$p(\mathbf{x}_{t})=\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{T},\bm{\Sigma}_{T})$
and time-reversed state transitions
$p(\mathbf{x}_{t}|\mathbf{x}_{t+1})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}_{t}^{\leftarrow}\mathbf{x}_{t+1}+\mathbf{b}_{t}^{\leftarrow},\mathbf{C}_{t}^{\leftarrow})$.
Thus, in Algorithm B:

  * •

Step 1 calculates the marginal prior distributions of the states. To see this,
note that this step is effectively the Kalman-filter recursion without
observations. Note that if the prior dynamics are stationary with stationary
distribution $\operatorname{N}(\bm{\mu},\bm{\Sigma})$, then Step 1 can be
skipped (because then $\bm{\mu}_{t}=\bm{\mu}$ and
$\bm{\Sigma}_{t}=\bm{\Sigma}$, for any $t\in[T]$).

  * •

Step 2 computes the parameters of the time-reversed transition kernels and
follows from standard Gaussian algebra (see, e.g., Särkkä and Svensson, ,
2023, Section A.1) by noting that

$\begin{bmatrix}\mathbf{x}_{t}\\\
\mathbf{x}_{t+1}\end{bmatrix}\sim\operatorname{N}\biggl{(}\begin{bmatrix}\bm{\mu}_{t}\\\
\bm{\mu}_{t+1}\end{bmatrix},\begin{bmatrix}\bm{\Sigma}_{t}&\mathbf{F}_{t+1}\bm{\Sigma}_{t}\\\
\bm{\Sigma}_{t}\mathbf{F}_{t+1}^{\mathrm{T}}&\bm{\Sigma}_{t+1}\end{bmatrix}\biggr{)}.$
(85)

  * •

Step 4 derives (84) from (83) and corresponds to a single update step of a
Kalman filter (Särkkä and Svensson, , 2023, Chapter 6, Equation 6.21), where
$\mathbf{x}_{t-1}$ plays the rôle of an observation.

##### Alternative algorithm for invertible covariance matrices.

If $\mathbf{C}_{t}$ is invertible for all $t\in[T]$, then the twisted-proposal
parameters can be alternatively computed via Algorithm B which may be slightly
simpler to implement for some users and which may provide additional numerical
advantages in the case of explosive prior dynamics.

[twisted-proposal parameters: alternative] At the start of an iteration of the
twisted Particle-aGRAD, Particle-aGRAD+, Particle-aPCNL or Particle-aPCNL+
algorithm (after having sampled all the auxiliary variables $\mathbf{u}_{1:T}$
upfront – e.g., as in Algorithm D.1 from Appendix D.1 – which is possible
because these only depend on the reference path),

  1. 1.

run the Kalman filtering recursion to compute the moments of
$p(\mathbf{x}_{t}|\mathbf{u}_{1:t-1})=\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t|t-1},\bm{\Sigma}_{t|t-1})$,
for $t=1,\dotsc,T$,

  2. 2.

run the Kalman smoothing (a.k.a. Rauch–Tung–Striebel smoothing) recursion to
compute the moments of
$p(\mathbf{x}_{t}|\mathbf{u}_{1:T})=\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t|T},\bm{\Sigma}_{t|T})$,
for $t=T,T-1,\dotsc,1$,

  3. 3.

for $t\in[T]$, set

$\displaystyle\mathbf{C}^{\prime}_{t}$
$\displaystyle\coloneqq[\mathbf{C}_{t}^{-1}+\bm{\Sigma}_{t|T}^{-1}-\bm{\Sigma}_{t|t-1}^{-1}]^{-1},$
(86) $\displaystyle\mathbf{F}^{\prime}_{t}$
$\displaystyle\coloneqq\mathbf{C}^{\prime}_{t}[\mathbf{C}_{t}^{-1}\mathbf{b}_{t}+\bm{\Sigma}_{t|T}^{-1}\bm{\mu}_{t|T}-\bm{\Sigma}_{t|t-1}^{-1}\bm{\mu}_{t|t-1}],$
(87) $\displaystyle\mathbf{b}^{\prime}_{t}$
$\displaystyle\coloneqq\mathbf{C}^{\prime}_{t}\mathbf{C}_{t}^{-1}\mathbf{F}_{t}.$
(88)

Algorithm B is justified by the decomposition

$\displaystyle p(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{u}_{1:T})$
$\displaystyle\propto p(\mathbf{x}_{t-1:t},\mathbf{u}_{1:T})$ (89)
$\displaystyle=p(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{u}_{1:T})p(\mathbf{x}_{t}|\mathbf{u}_{1:T})$
(90)
$\displaystyle\propto\frac{p(\mathbf{x}_{t}|\mathbf{x}_{t-1})}{p(\mathbf{x}_{t}|\mathbf{u}_{1:t-1})}p(\mathbf{x}_{t}|\mathbf{u}_{1:T})$
(91)
$\displaystyle=\frac{\operatorname{N}(\mathbf{x}_{t};\mathbf{F}_{t}\mathbf{x}_{t-1}+\mathbf{b}_{t},\mathbf{C}_{t})}{\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t|t-1},\bm{\Sigma}_{t|t-1})}\operatorname{N}(\mathbf{x}_{t};\bm{\mu}_{t|T},\bm{\Sigma}_{t|T})$
(92)
$\displaystyle\propto\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t}).$
(93)

##  Appendix C Integrating out the auxiliary variables

In this section, we prove a few lemmata which are used in subsequent sections.

  * •

Lemmata 1 and 2. Lemmata 1 and 2 derive the determinant and inverse of a
certain simple block matrix which appears repeatedly in the remainder of this
section and also in Appendix E.

  * •

Lemma 3. Lemma 3 will allow us to derive marginal proposal distributions of
the various algorithms, i.e., the distribution of
$\smash{\mathbf{x}_{t}^{-k_{t}}=(\mathbf{x}_{t}^{0},\dotsc,\mathbf{x}_{t}^{k_{t}-1},\mathbf{x}_{t}^{k_{t}+1},\dotsc,\mathbf{x}_{t}^{N})}$
conditional on $k_{t}$, $\mathbf{x}_{t}^{k_{t}}=\mathbf{x}_{t}$ and all the
particles and ancestor indices with time indices $s<t$, but with the auxiliary
variables $\mathbf{u}_{1:T}$ integrated out.

  * •

Lemma 4. Lemma 4 will allow us to evaluate the particle weights used in the
‘marginal’ algorithms (Particle-MALA, Particle-mGRAD and Particle-PCNL) at
linear complexity in $N$ although the weight of the $n$th particle depends on
the values of all other particles.

###  C.1 Properties of a particular block matrix

Let $\mathbf{I}_{M}$ denote the $(M\times M)$ identity matrix. When $M=D$ we
continue to leave out the subscript. Furthermore, let $\bm{1}_{M\times
N}\in\\{1\\}^{M\times N}$ and denote a matrix in which every element is $1$.
For matrices $\mathbf{A},\mathbf{B}\in\mathbb{R}^{D\times D}$, define the
block matrix

$\displaystyle\mathcal{M}_{N}(\mathbf{A},\mathbf{B})$
$\displaystyle\coloneqq\mathbf{I}_{N}\otimes\mathbf{A}+\bm{1}_{N\times
N}\otimes\mathbf{B}=\begin{bmatrix}\mathbf{A}+\mathbf{B}&\mathbf{B}&\dotsc&\mathbf{B}\\\
\mathbf{B}&\mathbf{A}+\mathbf{B}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\mathbf{B}\\\
\mathbf{B}&\dotsc&\mathbf{B}&\mathbf{A}+\mathbf{B}\end{bmatrix}\in\mathbb{R}^{(DN)\times(DN)}.$
(94)

######  Lemma 1.

For $N,D\in\mathbb{N}$, let $\mathbf{A},\mathbf{B}\in\mathbb{R}^{D\times D}$.
Then,

$\displaystyle\det(\mathcal{M}_{N}(\mathbf{A},\mathbf{B}))=\det(\mathbf{A})^{N-1}\det(\mathbf{A}+N\mathbf{B}).$
(95)

######  Proof.

Subtracting the last row of $\mathcal{M}_{N}(\mathbf{A},\mathbf{B})$ from all
other rows and then adding the sum of the first $N-1$ columns to the last
column gives the upper-triangular block matrix

$\displaystyle\begin{bmatrix}\mathbf{A}&\bm{0}&\dotsc&\bm{0}\\\
\mathbf{B}&\ddots&\ddots&\vdots\\\ \vdots&\ddots&\mathbf{A}&\bm{0}\\\
\mathbf{B}&\dotsc&\mathbf{B}&\mathbf{A}+N\mathbf{B}\end{bmatrix}.$ (96)

This proves the result. $\Box$□

######  Lemma 2.

For $N,D\in\mathbb{N}$, let $\mathbf{A},\mathbf{B}\in\mathbb{R}^{D\times D}$,
such that $\mathbf{A}$ and $(\mathbf{A}+N\mathbf{B})$ are invertible. Then

$\displaystyle\mathcal{M}_{N}(\mathbf{A},\mathbf{B})^{-1}=\mathcal{M}_{N}(\mathbf{A}^{-1},-(\mathbf{A}+N\mathbf{B})^{-1}\mathbf{B}\mathbf{A}^{-1}).$
(97)

######  Proof.

We must have
$\mathcal{M}_{N}(\mathbf{A},\mathbf{B})\mathcal{M}_{N}(\mathbf{F},\mathbf{G})=\mathbf{I}_{DN}$
and hence

$\displaystyle(\mathbf{A}+\mathbf{B})(\mathbf{F}+\mathbf{G})+(N-1)\mathbf{B}\mathbf{G}$
$\displaystyle=\mathbf{I},$ (98)
$\displaystyle\mathbf{B}(\mathbf{F}+\mathbf{G})+(\mathbf{A}+\mathbf{B})\mathbf{G}+(N-2)\mathbf{B}\mathbf{G}$
$\displaystyle=\bm{0}.$ (99)

This implies $\mathbf{F}=\mathbf{A}^{-1}$ and
$\mathbf{G}=-(\mathbf{A}+N\mathbf{B})^{-1}\mathbf{B}\mathbf{A}^{-1}$. $\Box$□

###  C.2 Conditional and marginal proposal distributions

In this section, for any tuple $(\mathbf{z}_{0},\dotsc,\mathbf{z}_{N})$ of
values in $\mathcal{X}\coloneqq\mathbb{R}^{D}$ and any $n\in[N]_{0}$, we write
$\mathbf{z}_{-n}\coloneqq(\mathbf{z}_{0},\dotsc,\mathbf{z}_{n-1},\mathbf{z}_{n+1},\dotsc,\mathbf{z}_{N})$.
Given some $N\in\mathbb{N}$, $n\in[N]_{0}$,
$\mathbf{x}_{n},\bm{\phi}_{n}\in\mathcal{X}$, we consider the following joint
distribution on $\mathcal{X}^{N+1}$:

$\displaystyle q_{-n}(\mathbf{u},\mathbf{x}_{-n}|\mathbf{x}_{n})$
$\displaystyle\coloneqq\operatorname{N}(\mathbf{u};\mathbf{x}_{n}+\bm{\phi}_{n},\mathbf{E})\prod_{\begin{subarray}{c}\mathllap{m}=\mathrlap{0}\\\
\mathllap{m}\neq\mathrlap{n}\end{subarray}}^{N}\operatorname{N}(\mathbf{x}_{m};\mathbf{v}_{m}+\mathbf{H}_{m}\mathbf{u},\mathbf{D}_{m}),$
(100)

where, for any $m\in[N]_{0}$, $\mathbf{v}_{m}\in\mathcal{X}$ and
$\mathbf{H}_{m}\in\mathbb{R}^{D\times D}$, and
$\mathbf{D}_{m},\mathbf{E}\in\mathbb{R}^{D\times D}$ are positive definite and
symmetric.

To simplify the presentation – and with some abuse of notation since we use
the same symbols for tuples and their vectorised versions – we write

$\displaystyle\mathbf{x}_{-n}$
$\displaystyle\coloneqq\begin{bmatrix}\mathbf{x}_{0}\\\ \vdots\\\
\mathbf{x}_{n-1}\\\ \mathbf{x}_{n+1}\\\ \vdots\\\
\mathbf{x}_{n}\end{bmatrix},\quad\mathbf{v}_{-n}\coloneqq\begin{bmatrix}\mathbf{v}_{0}\\\
\vdots\\\ \mathbf{v}_{n-1}\\\ \mathbf{v}_{n+1}\\\ \vdots\\\
\mathbf{v}_{n}\end{bmatrix},\quad\mathbf{H}_{-n}\coloneqq\begin{bmatrix}\mathbf{H}_{0}\\\
\vdots\\\ \mathbf{H}_{n-1}\\\ \mathbf{H}_{n+1}\\\ \vdots\\\
\mathbf{H}_{n}\end{bmatrix},\quad\mathbf{D}_{-n}\coloneqq\operatorname{diag}\left(\begin{bmatrix}\mathbf{D}_{0}\\\
\vdots\\\ \mathbf{D}_{n-1}\\\ \mathbf{D}_{n+1}\\\ \vdots\\\
\mathbf{D}_{n}\end{bmatrix}\right),$ (101)

where, in the last expression, $\operatorname{diag}$ induces a block-diagonal
matrix. With this notation, we can formulate (100) equivalently as

$\displaystyle
q_{-n}(\mathbf{u},\mathbf{x}_{-n}|\mathbf{x}_{n})=\operatorname{N}(\mathbf{u};\mathbf{x}_{n}+\bm{\phi}_{n},\mathbf{E})\operatorname{N}(\mathbf{x}_{-n};\mathbf{v}_{-n}+\mathbf{H}_{-n}\mathbf{u},\mathbf{D}_{-n}).$
(102)

######  Lemma 3.

For any $n\in[N]_{0}$ and $\mathbf{x}_{n}\in\mathcal{X}$,

  1. 1.

the marginal distribution of $\mathbf{x}_{-n}$ given $\mathbf{x}_{n}$ under
(102) is

$\displaystyle q_{-n}(\mathbf{x}_{-n}|\mathbf{x}_{n})$
$\displaystyle=\operatorname{N}(\mathbf{x}_{-n};\bm{\mu}_{-n},\bm{\Sigma}_{-n}),$
(103)

where

$\displaystyle\bm{\mu}_{-n}$
$\displaystyle\coloneqq\mathbf{v}_{-n}+\mathbf{H}_{-n}(\mathbf{x}_{n}+\bm{\phi}_{n}),$
(104) $\displaystyle\bm{\Sigma}_{-n}$
$\displaystyle\coloneqq\mathbf{D}_{-n}+\mathbf{H}_{-n}\mathbf{E}\mathbf{H}_{-n}^{\mathrm{T}};$
(105)

  2. 2.

the conditional distribution of $\mathbf{u}$ given $\mathbf{x}_{-n}$ and
$\mathbf{x}_{n}$ under (102) is

$\displaystyle\\!\\!\\!q_{-n}(\mathbf{u}|\mathbf{x}_{-n},\mathbf{x}_{n})$
$\displaystyle=\operatorname{N}(\mathbf{u};\mathbf{x}_{n}+\bm{\phi}_{n}+\mathbf{K}[\mathbf{x}_{-n}-\mathbf{v}_{-n}-\mathbf{H}_{-n}(\mathbf{x}_{n}+\bm{\phi}_{n})],(\mathbf{I}-\mathbf{K}\mathbf{H}_{-n})\mathbf{E}),\\!\\!\\!$
(106)

where
$\mathbf{K}\coloneqq\mathbf{E}\mathbf{H}_{-n}^{\mathrm{T}}(\mathbf{D}_{-n}+\mathbf{H}_{-n}\mathbf{E}\mathbf{H}_{-n}^{\mathrm{T}})^{-1}$.

######  Proof.

This follows by simple algebra (see, e.g., Särkkä and Svensson, , 2023,
Appendix A.1). $\Box$□

######  Lemma 4.

Assume now that $\mathbf{H}_{m}=\mathbf{H}$ and $\mathbf{D}_{m}=\mathbf{D}$,
for any $m\in[N]_{0}$. Then, with the notation from Lemma 3,

$\displaystyle q_{-n}(\mathbf{x}_{-n}|\mathbf{x}_{n})$ $\displaystyle\propto
H_{\bm{\phi}_{n}}(\mathbf{x}_{n},\mathbf{v}_{n},\bar{\mathbf{x}},\bar{\mathbf{v}})I((\mathbf{x}_{m}-\mathbf{v}_{m})_{m=0}^{N}),$
(107)

where

  1. 1.

$\mathbf{z}_{0:N}\mapsto I(\mathbf{z}_{0:N})$ is invariant under any
permutation of its arguments;

  2. 2.

$\bar{\mathbf{x}}\coloneqq\frac{1}{N+1}\sum_{m=0}^{N}\mathbf{x}_{n}$, and
$\bar{\mathbf{v}}\coloneqq\frac{1}{N+1}\sum_{m=0}^{N}\mathbf{v}_{n}$, and

$\displaystyle\log
H_{\bm{\phi}}(\mathbf{x},\mathbf{v},\bar{\mathbf{x}},\bar{\mathbf{v}})$ (108)
$\displaystyle=\tfrac{1}{2}(\mathbf{x}-\mathbf{v})^{\mathrm{T}}(\mathbf{D}^{-1}+\mathbf{G})(\mathbf{x}-\mathbf{v})$
(109)
$\displaystyle\qquad-\tfrac{1}{2}N(\mathbf{x}+\bm{\phi})^{\mathrm{T}}\mathbf{H}^{\mathrm{T}}(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}(\mathbf{x}+\bm{\phi})$
(110)
$\displaystyle\qquad-(N+1)(\bar{\mathbf{x}}-\bar{\mathbf{v}})^{\mathrm{T}}[\mathbf{G}(\mathbf{x}-\mathbf{v})-(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}(\mathbf{x}+\bm{\phi})]$
(111)
$\displaystyle\qquad-(\mathbf{x}-\mathbf{v})^{\mathrm{T}}(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}(\mathbf{x}+\bm{\phi}),$
(112)

whose evaluation complexity does not depend on $N$. Here,

$\displaystyle\mathbf{G}$
$\displaystyle\coloneqq(\mathbf{D}+N\mathbf{H}\mathbf{E}\mathbf{H}^{\mathrm{T}})^{-1}\mathbf{H}\mathbf{E}\mathbf{H}^{\mathrm{T}}\mathbf{D}^{-1}$
(113)
$\displaystyle=\mathbf{D}^{-1}\mathbf{H}\mathbf{E}(\mathbf{E}+N\mathbf{E}\mathbf{H}^{\mathrm{T}}\mathbf{D}^{-1}\mathbf{H}\mathbf{E})^{-1}\mathbf{E}\mathbf{H}^{\mathrm{T}}\mathbf{D}^{-1}.$
(114)

######  Proof.

The equivalence of (113) and (114) follows from the _push-through identity_
(Henderson and Searle, , 1981). By assumption,
$\bm{\Sigma}_{-n}=\mathcal{M}_{N}(\mathbf{D},\mathbf{H}\mathbf{E}\mathbf{H}^{\mathrm{T}})$.
Thus, Lemma 2 gives

$\displaystyle\bm{\Sigma}_{-n}^{-1}=\mathcal{M}_{N}(\mathbf{D}^{-1},-\mathbf{G}).$
(115)

In particular, letting ${\otimes}$ be the Kronecker product, this implies that

$\displaystyle\bm{\Sigma}_{-n}^{-1}\mathbf{H}_{-n}$
$\displaystyle=\bm{1}_{N\times
1}\otimes[(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}],$ (116)
$\displaystyle\mathbf{H}_{-n}^{\mathrm{T}}\bm{\Sigma}_{-n}^{-1}\mathbf{H}_{-n}$
$\displaystyle=N\mathbf{H}^{\mathrm{T}}(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}.$
(117)

Therefore, defining

$\displaystyle\mathbf{x}\coloneqq\begin{bmatrix}\mathbf{x}_{0}\\\ \vdots\\\
\mathbf{x}_{N}\end{bmatrix},\quad\mathbf{v}\coloneqq\begin{bmatrix}\mathbf{v}_{0}\\\
\vdots\\\
\mathbf{v}_{N}\end{bmatrix},\quad\bm{\Sigma}\coloneqq\mathcal{M}_{N+1}(\mathbf{D}^{-1},-\mathbf{G})^{-1},$
(118)

we have

$\displaystyle q_{-n}(\mathbf{x}_{-n}|\mathbf{x}_{n})$
$\displaystyle\propto\exp\bigl{(}-\tfrac{1}{2}\bigl{[}(\mathbf{x}_{-n}-\mathbf{v}_{-n}-\mathbf{H}_{-n}(\mathbf{x}_{n}+\bm{\phi}_{n}))^{\mathrm{T}}\bm{\Sigma}_{-n}^{-1}(\mathbf{x}_{-n}-\mathbf{v}_{-n}-\mathbf{H}_{-n}(\mathbf{x}_{n}+\bm{\phi}_{n}))\bigr{]}\bigr{)}$
$\displaystyle=\exp\bigl{(}-\tfrac{1}{2}\bigl{[}(\mathbf{x}_{-n}-\mathbf{v}_{-n})^{\mathrm{T}}\bm{\Sigma}_{-n}^{-1}(\mathbf{x}_{-n}-\mathbf{v}_{-n})$
$\displaystyle\qquad\qquad\quad+(\mathbf{x}_{n}+\bm{\phi}_{n})^{\mathrm{T}}\mathbf{H}_{-n}^{\mathrm{T}}\bm{\Sigma}_{-n}^{-1}\mathbf{H}_{-n}(\mathbf{x}_{n}+\bm{\phi}_{n})$
$\displaystyle\qquad\qquad\quad-2(\mathbf{x}_{-n}-\mathbf{v}_{-n})^{\mathrm{T}}\bm{\Sigma}_{-n}^{-1}\mathbf{H}_{-n}(\mathbf{x}_{n}+\bm{\phi}_{n})\bigr{]}\bigr{)}$
$\displaystyle=\exp\bigl{(}-\tfrac{1}{2}\bigl{[}(\mathbf{x}-\mathbf{v})^{\mathrm{T}}\bm{\Sigma}^{-1}(\mathbf{x}-\mathbf{v})$
$\displaystyle\qquad\qquad\quad-(\mathbf{x}_{n}-\mathbf{v}_{n})^{\mathrm{T}}(\mathbf{D}^{-1}-\mathbf{G})(\mathbf{x}_{n}-\mathbf{v}_{n})$
$\displaystyle\qquad\qquad\quad+2(\mathbf{x}_{-n}-\mathbf{v}_{-n})^{\mathrm{T}}[\bm{1}_{N\times
1}\otimes\mathbf{G}](\mathbf{x}_{n}-\mathbf{v}_{n})$
$\displaystyle\qquad\qquad\quad+N(\mathbf{x}_{n}+\bm{\phi}_{n})^{\mathrm{T}}\mathbf{H}^{\mathrm{T}}(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}(\mathbf{x}_{n}+\bm{\phi}_{n})$
$\displaystyle\qquad\qquad\quad-2(\mathbf{x}_{-n}-\mathbf{v}_{-n})^{\mathrm{T}}\bm{\Sigma}_{-n}^{-1}\mathbf{H}_{-n}(\mathbf{x}_{n}+\bm{\phi}_{n})\bigr{]}\bigr{)}$
$\displaystyle=\exp\bigl{(}-\tfrac{1}{2}\bigl{[}(\mathbf{x}-\mathbf{v})^{\mathrm{T}}\bm{\Sigma}^{-1}(\mathbf{x}-\mathbf{v})$
$\displaystyle\qquad\qquad\quad-(\mathbf{x}_{n}-\mathbf{v}_{n})^{\mathrm{T}}(\mathbf{D}^{-1}-\mathbf{G})(\mathbf{x}_{n}-\mathbf{v}_{n})$
$\displaystyle\qquad\qquad\quad+N(\mathbf{x}_{n}+\bm{\phi}_{n})^{\mathrm{T}}\mathbf{H}^{\mathrm{T}}(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}(\mathbf{x}_{n}+\bm{\phi}_{n})$
$\displaystyle\qquad\qquad\quad+2(\mathbf{x}-\mathbf{v})^{\mathrm{T}}\left\\{\bm{1}_{(N+1)\times
1}\otimes[\mathbf{G}(\mathbf{x}_{n}-\mathbf{v}_{n})-(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}(\mathbf{x}_{n}+\bm{\phi}_{n})]\right\\}$
$\displaystyle\qquad\qquad\quad-2(\mathbf{x}_{n}-\mathbf{v}_{n})^{\mathrm{T}}[\mathbf{G}(\mathbf{x}_{n}-\mathbf{v}_{n})-(\mathbf{D}^{-1}-N\mathbf{G})\mathbf{H}(\mathbf{x}_{n}+\bm{\phi}_{n})]\bigr{]}\bigr{)}$
$\displaystyle=H_{\bm{\phi}_{n}}(\mathbf{x}_{n},\mathbf{v}_{n},\bar{\mathbf{x}},\bar{\mathbf{v}})I((\mathbf{x}_{m}-\mathbf{v}_{m})_{m=0}^{N}),$

with

$\displaystyle
I((\mathbf{x}_{m}-\mathbf{v}_{m})_{m=0}^{N})\propto\exp(-\tfrac{1}{2}(\mathbf{x}-\mathbf{v})^{\mathrm{T}}\bm{\Sigma}^{-1}(\mathbf{x}-\mathbf{v})).$
(119)

This completes the proof. $\Box$□

##  Appendix D Generic algorithms and proof of Propositions 1–13

In this section, we prove that the algorithms proposed in this work leave
$\pi_{T}$ invariant. To this end, we first prove the validity of two generic
algorithms.

  * •

Generic auxiliary algorithm. The first generic algorithm includes auxiliary
variables $\mathbf{u}_{t}$ in the space and admits the ‘auxiliary-variable’
based algorithms: Particle-aMALA, Particle-aGRAD, Particle-aPCNL as well as
their smoothing-gradient (‘+’) and twisted versions, as special cases. Its
proof extends the auxiliary-variable interpretation of the Particle-RWM
algorithm which was given in Corenflos and Särkkä,  (2023).

  * •

Generic marginal algorithm. The second generic algorithm integrates out the
auxiliary variables and admits the ‘marginal’ algorithms from the main
manuscript (Particle-MALA, Particle-mGRAD, Particle-PCNL). Its proof relies on
an argument previously given in Finke et al.,  (2016).

###  D.1 Generic auxiliary algorithm

Define an extended target 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}+\Phi_{t}(\mathbf{x}_{1:T}),\mathbf{E}_{t}(\mathbf{x}_{t-1:t})),$
(120)

where, for any $t\in[T]$, $\Phi_{t}\colon\mathcal{X}^{T}\to\mathcal{X}$ is a
function satisfying

$\displaystyle\Phi_{t}(\mathbf{x}_{1:T})=\begin{cases}\bm{\phi}_{T}(\mathbf{x}_{t-T:T}),&\text{if
$t=T$,}\\\
\bm{\phi}_{t}(\mathbf{x}_{t-1:t})+\bm{\psi}_{t}(\mathbf{x}_{t:t+1}),&\text{otherwise,}\end{cases}$
(121)

and $\mathbf{E}_{t}(\mathbf{x}_{t-1:t})\in\mathbb{R}^{D\times D}$ is some
positive-definite symmetric matrix. Additionally, let

$\displaystyle Q^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{1:T})$
$\displaystyle=M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{1:T}),$
(122)

for some mutation kernel
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})$ and some
potential function $G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{1:T})$ (both
of which may depend on some or all of $\mathbf{u}_{1},\dotsc,\mathbf{u}_{t}$)
such that

$\displaystyle\pi^{\prime}_{T}(\mathbf{x}_{1:T}|\mathbf{u}_{1:T})\propto\prod_{t=1}^{T}Q^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{1:T}).$
(123)

[generic auxiliary algorithm] Given $\mathbf{x}_{1:T}\in\mathcal{X}^{T}$,
sample

$\displaystyle\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\Phi_{t}(\mathbf{x}_{1:T}),\mathbf{E}_{t}(\mathbf{x}_{t-1:t})),$
(124)

for any $t=1,\dotsc,T$ and then

  1. 1.

for $t=1,\dotsc,T$,

    1. (a)

sample $k_{t}$ from a uniform distribution on $[N]_{0}$ and set
$\mathbf{x}_{t}^{k_{t}}\coloneqq\mathbf{x}_{t}$,

    2. (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}\\}$,

    3. (c)

sample $\smash{\mathbf{x}_{t}^{n}\sim
M^{\prime}_{t}(\,\cdot\,|\mathbf{x}_{t-1}^{a_{t-1}^{n}};\mathbf{u}_{1:T})}$
for $n\in[N]_{0}\setminus\\{k_{t}\\}$,

    4. (d)

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

    5. (e)

for $n\in[N]_{0}$, set $W_{t}^{n}\coloneqq w_{t}^{n}/\sum_{m=0}^{N}w_{t}^{m}$;

  2. 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.
$1\wedge\dfrac{1-W_{T}^{k_{T}}}{1-W_{T}^{i}}$; otherwise, set $l_{T}\coloneqq
k_{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}_{1:T})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{i},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{1:T})}{\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}_{1:T})Q^{\prime}_{t+2}((\mathbf{x}_{t}^{n},\mathbf{x}_{t+1}^{l_{t+1}},\mathbf{x}_{t+2}^{l_{t+2}});\mathbf{u}_{1:T})};$
(125)

  4. 4.

return
$\tilde{\mathbf{x}}_{1:T}\coloneqq(\mathbf{x}_{1}^{l_{1}},\dotsc,\mathbf{x}_{t}^{l_{T}})$.

######  Proposition 14 (validity of the generic auxiliary algorithm).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm D.1
induces a Markov kernel $P(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which
leaves $\pi_{T}$ invariant.

######  Proof (of Proposition 14).

The extended distribution from (120) admits $\pi_{T}(\mathbf{x}_{1:T})$ as a
marginal. Therefore, a valid MCMC update for sampling from this extended
distribution is given 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}+\Phi_{t}(\mathbf{x}_{1:T}),\mathbf{E}_{t}(\mathbf{x}_{t-1:t}))$,
for $t=1,\dotsc,T$;

  2. 2.

run a standard CSMC algorithm with backward sampling (as in Algorithm 2.1.1)
targeting $\pi^{\prime}_{T}(\mathbf{x}_{1:T}|\mathbf{u}_{1:T})$ 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
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})$,
$G^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{1:T})$ and
$Q^{\prime}_{t}(\mathbf{x}_{t-2:t};\mathbf{u}_{1:T})$, and with appropriate
adjustments (e.g., of the backward kernels) to account for the possibility
that the model may only be second-order Markov.

These to steps are equivalent to Algorithm D.1. $\Box$□

###  D.2 Generic marginal algorithm

Consider the same setting as above but now assume that for any $t\in[T]$,
$\bm{\psi}_{t}\equiv\bm{0}$, so that
$\Phi_{t}(\mathbf{x}_{1:T})=\bm{\phi}_{t}(\mathbf{x}_{t-1:t})$ as well as that
$\mathbf{E}_{t}(\mathbf{x}_{t-1:t})=\mathbf{E}_{t}$ is independent of
$\mathbf{x}_{t-1:t}$.

Furthermore, assume that
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})=M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})$
only depends on the $t$th auxiliary variable $\mathbf{u}_{t}$ and,
specifically, is a Gaussian distribution of the following form:

$\displaystyle M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{t})$
$\displaystyle\coloneqq\operatorname{N}(\mathbf{x}_{t};\mathbf{v}_{t}(\mathbf{x}_{t-1})+\mathbf{H}_{t}\mathbf{u}_{t},\mathbf{D}_{t}),$
(126)

where $\mathbf{v}_{t}(\mathbf{x})\in\mathcal{X}$ whilst
$\mathbf{H}_{t},\mathbf{D}_{t}\in\mathbb{R}^{D\times D}$ do not depend on
$\mathbf{x}_{t-1}$ and define

$\displaystyle
q_{t}^{-n}(\mathbf{x}_{t}^{-n},\mathbf{u}_{t}|\mathbf{x}_{t}^{n};\mathcal{H}_{t-1})\coloneqq\operatorname{N}(\mathbf{u}_{t};\mathbf{x}_{t}^{n}+\bm{\phi}_{t}(\mathbf{x}_{t-1:t}^{(n)}),\mathbf{E}_{t})\prod_{\begin{subarray}{c}\mathllap{m}=\mathrlap{0}\\\
\mathllap{m}\neq\mathrlap{n}\end{subarray}}^{N}M^{\prime}_{t}(\mathbf{x}_{t}^{m}|\mathbf{x}_{t-1}^{a_{t-1}^{m}};\mathbf{u}_{t}),$
(127)

where $\mathcal{H}_{t-1}$ is the history of the particle system up to time
$t-1$, i.e., all particles and ancestor indices with ‘time’ subscript $s\leq
t-1$. By Lemma 3 from Appendix C, we obtain a closed-form expression for

$\displaystyle
q_{t}^{-n}(\mathbf{x}_{t}^{-n}|\mathbf{x}_{t}^{n};\mathcal{H}_{t-1})\coloneqq\int_{\mathcal{X}}q_{t}^{-n}(\mathbf{x}_{t}^{-n},\mathbf{u}_{t}|\mathbf{x}_{t}^{n};\mathcal{H}_{t-1})\,\mathrm{d}\mathbf{u}_{t}.$
(128)

[generic marginal algorithm] Given $\mathbf{x}_{1:T}\in\mathcal{X}^{T}$:

  1. 1.

for $t=1,\dotsc,T$,

    1. (a)

sample $k_{t}$ from a uniform distribution on $[N]_{0}$ and set
$\mathbf{x}_{t}^{k_{t}}\coloneqq\mathbf{x}_{t}$,

    2. (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}\\}$,

    3. (c)

sample $\smash{\mathbf{x}_{t}^{-k_{t}}\sim
q_{t}^{-k_{t}}(\mathbf{x}_{t}^{-k_{t}}|\mathbf{x}_{t}^{k_{t}};\mathcal{H}_{t-1})}$
(e.g. by sampling
$\smash{\mathbf{u}_{t}\sim\operatorname{N}(\mathbf{x}_{t}+\bm{\phi}_{t}(\mathbf{x}_{t-1:t}),\mathbf{E}_{t})}$
and then $\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}\\}$),

    4. (d)

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

    5. (e)

for $n\in[N]_{0}$, set $W_{t}^{n}\coloneqq w_{t}^{n}/\sum_{m=0}^{N}w_{t}^{m}$;

  2. 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.
$1\wedge\dfrac{1-W_{T}^{k_{T}}}{1-W_{T}^{i}}$; otherwise, set $l_{T}\coloneqq
k_{T}$;

  3. 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. 4.

return
$\tilde{\mathbf{x}}_{1:T}\coloneqq(\mathbf{x}_{1}^{l_{1}},\dotsc,\mathbf{x}_{t}^{l_{T}})$.

Algorithm D.2 can be implemented in $\operatorname{\mathrm{O}}(N)$ operations
because Lemma 3 from Appendix C allows us to write the weight in Step 1d as

$\displaystyle w_{t}^{n}$ $\displaystyle\propto
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})q_{t}^{-n}(\mathbf{x}_{t}^{-n}|\mathbf{x}_{t}^{n};\mathcal{H}_{t-1})$
(129) $\displaystyle\propto
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})\smash{H_{t,\bm{\phi}_{t}(\mathbf{x}_{t-1:t}^{(n)})}(\mathbf{x}_{t}^{n},\mathbf{v}_{t}^{n},\bar{\mathbf{x}}_{t},\bar{\mathbf{v}}_{t})},$
(130)

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

$\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}}(\mathbf{D}_{t}^{-1}+\mathbf{G}_{t})(\mathbf{x}-\mathbf{v})$
$\displaystyle\qquad-\tfrac{1}{2}N(\mathbf{x}+\bm{\phi})^{\mathrm{T}}\mathbf{H}_{t}^{\mathrm{T}}(\mathbf{D}_{t}^{-1}-N\mathbf{G}_{t})\mathbf{H}_{t}(\mathbf{x}+\bm{\phi})$
$\displaystyle\qquad-(N+1)(\bar{\mathbf{x}}-\bar{\mathbf{v}})^{\mathrm{T}}[\mathbf{G}_{t}(\mathbf{x}-\mathbf{v})-(\mathbf{D}_{t}^{-1}-N\mathbf{G}_{t})\mathbf{H}_{t}(\mathbf{x}+\bm{\phi})]$
$\displaystyle\qquad-(\mathbf{x}-\mathbf{v})^{\mathrm{T}}(\mathbf{D}_{t}^{-1}-N\mathbf{G}_{t})\mathbf{H}_{t}(\mathbf{x}+\bm{\phi}),$

with
$\smash{\mathbf{G}_{t}\coloneqq(\mathbf{D}+N\mathbf{H}_{t}\mathbf{E}_{t}\mathbf{H}_{t}^{\mathrm{T}})^{-1}\mathbf{H}_{t}\mathbf{E}_{t}\mathbf{H}_{t}^{\mathrm{T}}\mathbf{D}_{t}^{-1}}$
(see (114) for an alternative expression).

######  Proposition 15 (validity of the generic marginal algorithm).

Sampling $\tilde{\mathbf{x}}_{1:T}$ given $\mathbf{x}_{1:T}$ via Algorithm D.2
induces a Markov kernel $P(\tilde{\mathbf{x}}_{1:T}|\mathbf{x}_{1:T})$ which
leaves $\pi_{T}$ invariant.

######  Proof (of Proposition 15).

We begin with a few observations.

  1. 1.

Since the unnormalised weights satisfy

$\displaystyle w_{t}^{n}\propto
Q_{t}(\mathbf{x}_{t-1:t}^{(n)})q_{t}^{-n}(\mathbf{x}_{t}^{-n}|\mathbf{x}_{t}^{n};\mathcal{H}_{t-1}),$
(131)

we have that

$\displaystyle
q_{t}^{-k_{t}}(\mathbf{x}_{t}^{-k_{t}}|\mathbf{x}_{t}^{k_{t}};\mathcal{H}_{t-1})=\frac{w_{t}^{k_{t}}}{w_{t}^{l_{t}}}\frac{Q_{t}(\mathbf{x}_{t-1:t}^{(l_{t})})}{Q_{t}(\mathbf{x}_{t-1:t}^{(k_{t})})}q_{t}^{-l_{t}}(\mathbf{x}_{t}^{-l_{t}}|\mathbf{x}_{t}^{l_{t}};\mathcal{H}_{t-1}).$
(132)

  2. 2.

For a given set of final-time weights $\\{W_{T}^{n}\\}_{n\in[N]_{0}}$, let
$R_{T}(\,\cdot\,|\,\cdot\,;\mathcal{H}_{T})$ be the
$\sum_{n=0}^{N}W_{T}^{n}\delta_{n}$-invariant Markov kernel used in Step 2 of
Algorithm D.2. That is, sampling $l_{T}\sim
R_{T}(\,\cdot\,|k_{T};\mathcal{H}_{T})$ could be the forced-move update; or,
in the more common specification of CSMC algorithms (Andrieu et al., , 2010),
i.e. without the forced-move update, we would simply have
$R_{T}(l_{T}|k_{T};\mathcal{H}_{T})=W_{T}^{l_{T}}$. In either case, it can
then be verified that

$\displaystyle
W_{T}^{k_{T}}R_{T}(l_{T}|k_{T};\mathcal{H}_{T})=W_{T}^{l_{T}}R_{T}(k_{T}|l_{T};\mathcal{H}_{T}),$
(133)

for any $k_{T},l_{T}\in[N]_{0}$.

  3. 3.

Under Algorithm D.2, we have the following identities (with probability $1$):
$\mathbf{x}_{t}=\mathbf{x}_{t}^{k_{t}}$ and
$\mathbf{x}_{t}^{\prime}=\mathbf{x}_{t}^{l_{t}}$, for $1\leq t\leq T$, as well
as $a_{t-1}^{k_{t}}=k_{t-1}$, for any $1<t\leq T$.

Putting these observations together then shows that the Algorithm D.2 targets
the following extended distribution (i.e., this is the distribution of all
random variables obtained if we first sampled $\mathbf{x}_{1:T}\sim\pi_{T}$
and then ran Algorithm D.2):

$\displaystyle\quad\frac{\pi_{T}(\mathbf{x}_{1:T})}{(N+1)^{T}}\delta_{\mathbf{x}_{1:T}}(\mathbf{x}_{1:T}^{k_{1:T}})\biggl{[}\prod_{t=1}^{T}q_{t}^{-k_{t}}(\mathbf{x}_{t}^{-k_{t}}|\mathbf{x}_{t}^{k_{t}};\mathcal{H}_{t-1})\biggr{]}\biggl{[}\prod_{t=2}^{T}\delta_{k_{t-1}}(a_{t-1}^{k_{t}})\prod_{\begin{subarray}{c}\mathllap{n}=\mathrlap{0}\\\
\mathllap{n}\neq\mathrlap{k_{t}}\end{subarray}}^{N}W_{t-1}^{a_{t-1}^{n}}\biggr{]}$
$\displaystyle\quad\times
R_{T}(l_{T}|k_{T};\mathcal{H}_{T})\biggl{[}\prod_{t=1}^{T-1}\frac{w_{t}^{l_{t}}Q_{t+1}(\mathbf{x}_{t}^{l_{t}},\mathbf{x}_{t+1}^{l_{t+1}})}{\sum_{m=0}^{N}w_{t}^{m}Q_{t+1}(\mathbf{x}_{t}^{m},\mathbf{x}_{t+1}^{l_{t+1}})}\biggr{]}\delta_{\mathbf{x}_{1:T}^{l_{1:T}}}(\tilde{\mathbf{x}}_{1:T})$
$\displaystyle=\frac{\pi_{T}(\tilde{\mathbf{x}}_{1:T})}{(N+1)^{T}}\delta_{\tilde{\mathbf{x}}_{1:T}}(\mathbf{x}_{1:T}^{l_{1:T}})\biggl{[}\prod_{t=1}^{T}q_{t}^{-l_{t}}(\mathbf{x}_{t}^{-l_{t}}|\mathbf{x}_{t}^{l_{t}};\mathcal{H}_{t-1})\biggr{]}$
$\displaystyle\qquad\times\biggl{[}\prod_{t=2}^{T}\frac{w_{t-1}^{a_{t-1}^{l_{t}}}Q_{t}(\mathbf{x}_{t-1}^{a_{t-1}^{l_{t}}},\mathbf{x}_{t}^{l_{t+1}})}{\sum_{m=0}^{N}w_{t-1}^{m}Q_{t}(\mathbf{x}_{t-1}^{m},\mathbf{x}_{t}^{l_{t}})}\prod_{\begin{subarray}{c}\mathllap{n}=\mathrlap{0}\\\
\mathllap{n}\neq\mathrlap{l_{t}}\end{subarray}}^{N}W_{t-1}^{a_{t-1}^{n}}\biggr{]}$
$\displaystyle\qquad\times
R_{T}(k_{T}|l_{T};\mathcal{H}_{T})\biggl{[}\prod_{t=1}^{T-1}\delta_{a_{t-1}^{k_{t}}}(k_{t-1})\biggr{]}\delta_{\mathbf{x}_{1:T}^{k_{1:T}}}(\mathbf{x}_{1:T}),$

where the r.h.s. is the distribution obtained if we first sampled
$\tilde{\mathbf{x}}_{1:T}\sim\pi_{T}$ and then ran Algorithm D.2 algorithm but
with ancestor sampling (Lindsten et al., , 2012) instead of backward sampling.
This is a modified version of the proof technique from Finke et al.,  (2016).
In other words, if $\mathbf{x}_{1:T}\sim\pi_{T}$ and if
$\tilde{\mathbf{x}}_{1:T}$ is sampled via Algorithm D.2, then
$\tilde{\mathbf{x}}_{1:T}\sim\pi_{T}$. This completes the proof. $\Box$□

###  D.3 Invariance of the algorithms

We can now easily verify the validity of the ‘auxiliary’ algorithms (Particle-
aMALA, Particle-aMALA+, Particle-aGRAD, Particle-aGRAD+, Particle-aPCNL,
Particle-aPCNL+, and twisted Particle-aGRAD/Particle-aGRAD+/Particle-
aPCNL/Particle-aPCNL+) by noting that these are special cases of Algorithm
D.1, and the validity of the ‘marginal’ algorithms (Particle-MALA, Particle-
mGRAD, Particle-PCNL) by noting that these are special cases of Algorithm D.2.

######  Proof (of Proposition 1).

This follows by taking
$\bm{\phi}_{t}(\mathbf{x}_{t-1:t})\coloneqq\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
Q_{t}(\mathbf{x}_{t-1:t})$, $\bm{\psi}_{t}\equiv\bm{0}$,
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})=\operatorname{N}(\mathbf{x}_{t};\mathbf{u}_{t},\tfrac{\delta_{t}}{2}\mathbf{I})$
and $\mathbf{E}_{t}\equiv\tfrac{\delta_{t}}{2}\mathbf{I}$ in Proposition 14.
$\Box$□

######  Proof (of Proposition 2).

This follows from Proposition 15 with the same setting as in Proposition 1. In
particular, in this case, $\mathbf{v}_{t}\equiv\bm{0}$,
$\mathbf{H}_{t}=\mathbf{I}$,
$\mathbf{D}_{t}=\mathbf{E}_{t}=\tfrac{\delta_{t}}{2}\mathbf{I}$. Consequently,
(D.2) then simplifies to (15), where we have used that
$\mathbf{G}_{t}=[\tfrac{\delta_{t}}{2}(N+1)]^{-1}\mathbf{I}=\mathbf{D}_{t}^{-1}/(N+1)$
and $\mathbf{D}_{t}^{-1}-N\mathbf{G}_{t}=\mathbf{G}_{t}$. $\Box$□

######  Proof (of Proposition 3).

This follows in the same way as the proof of Proposition 1 except that now
$\bm{\psi}_{t}(\mathbf{x}_{t:t+1})=\kappa\frac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
Q_{t+1}(\mathbf{x}_{t:t+1})$. $\Box$□

######  Proof (of Proposition 4).

This follows in the same way as the proof of Proposition 1 except that now
$\bm{\phi}_{t}(\mathbf{x}_{t-1:t})\coloneqq\kappa\tfrac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t}(\mathbf{x}_{t-1:t})$, and
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}^{\prime}_{t}(\mathbf{x}_{t-1},\mathbf{u}_{t}),\mathbf{C}^{\prime}_{t}(\mathbf{x}_{t-1}))$,
where $\mathbf{m}^{\prime}_{t}(\mathbf{x}_{t-1},\mathbf{u}_{t})$ and
$\mathbf{C}^{\prime}_{t}(\mathbf{x}_{t-1})$ are defined in (29) and (30).
$\Box$□

######  Proof (of Proposition 5).

This follows from Proposition 15 with the same setting as in Proposition 4. In
particular, in this case,
$\mathbf{H}_{t}=\mathbf{A}_{t}\coloneqq(\mathbf{C}_{t}+\tfrac{2}{\delta_{t}}\mathbf{I})^{-1}\mathbf{C}_{t}$,
$\mathbf{D}_{t}=\tfrac{\delta_{t}}{2}\mathbf{A}_{t}$ and
$\mathbf{E}_{t}=\tfrac{\delta_{t}}{2}\mathbf{I}$. Consequently, (D.2) then
simplifies to (4.2), where we have used that $\mathbf{A}_{t}$ is symmetric,
that
$\mathbf{H}_{t}^{\mathrm{T}}\mathbf{D}_{t}^{-1}=\mathbf{D}_{t}^{-1}\mathbf{H}_{t}=\tfrac{2}{\delta_{t}}\mathbf{I}$
and hence

$\displaystyle\mathbf{D}_{t}^{-1}-N\mathbf{G}_{t}$
$\displaystyle=\mathbf{A}_{t}^{-1}\mathbf{G}_{t}=\mathbf{G}_{t}\mathbf{A}_{t}^{-1}.$
(134)

This completes the proof. $\Box$□

######  Proof (of Proposition 6).

This follows in the same way as the proof of Proposition 4 except that now
$\bm{\psi}_{t}(\mathbf{x}_{t:t+1})=\kappa\frac{\delta_{t}}{2}\nabla_{\mathbf{x}_{t}}\log
G_{t+1}(\mathbf{x}_{t:t+1})$. $\Box$□

######  Proof (of Proposition 7).

This follows in the same way as the proof of Propositions 4 and 6,
respectively, but with
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t})$.
$\Box$□

######  Proof (of Proposition 10).

This follows in the same way as the proof of Proposition 1 except that now
$\bm{\phi}_{t}(\mathbf{x}_{t-1:t})\coloneqq\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})$,
$\mathbf{E}_{t}(\mathbf{x}_{t-1:t})\coloneqq\tfrac{\delta_{t}}{2}\mathbf{C}_{t}(\mathbf{x}_{t-1})$
and
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})=\operatorname{N}(\mathbf{x}_{t};\mathbf{m}^{\prime}_{t}(\mathbf{x}_{t-1},\mathbf{u}_{t}),\mathbf{C}^{\prime}_{t}(\mathbf{x}_{t-1}))$,
where $\mathbf{m}^{\prime}_{t}(\mathbf{x}_{t-1},\mathbf{u}_{t})$ and
$\mathbf{C}^{\prime}_{t}(\mathbf{x}_{t-1})$ are defined in (52) and (53).
$\Box$□

######  Proof (of Proposition 11).

This follows from Proposition 15 with the same setting as in Proposition 10.
In particular, in this case, $\mathbf{H}_{t}=\beta_{t}\mathbf{I}$,
$\mathbf{D}_{t}=(1-\beta_{t})\mathbf{C}_{t}$ and
$\mathbf{E}_{t}=\tfrac{\delta_{t}}{2}\mathbf{C}_{t}$. Consequently, (D.2) then
simplifies to (A.2), where we have used that
$\mathbf{E}\mathbf{H}_{t}^{\mathrm{T}}\mathbf{D}_{t}^{-1}=\mathbf{D}_{t}^{-1}\mathbf{H}_{t}\mathbf{E}=\mathbf{I}$
and hence

$\displaystyle\mathbf{D}_{t}^{-1}+\mathbf{G}_{t}$
$\displaystyle=(\beta_{t}^{-1}+N+1)\mathbf{G}_{t},$ (135)
$\displaystyle\mathbf{D}_{t}^{-1}-N\mathbf{G}_{t}$
$\displaystyle=\beta_{t}^{-1}\mathbf{G}_{t}.$ (136)

This completes the proof. $\Box$□

######  Proof (of Proposition 12).

This follows in the same way as the proof of Proposition 10 except that now
$\bm{\psi}_{t}(\mathbf{x}_{t:t+1})=\kappa\frac{\delta_{t}}{2}\smash{\widetilde{\mathbf{C}}}_{t}(\mathbf{x}_{t})\nabla_{\mathbf{x}_{t}}\log
G_{t+1}(\mathbf{x}_{t:t+1})$. $\Box$□

######  Proof (of Proposition 13).

This follows in the same way as the proof of Propositions 10 and 12,
respectively, but with
$M^{\prime}_{t}(\mathbf{x}_{t}|\mathbf{x}_{t-1};\mathbf{u}_{1:T})=\operatorname{N}(\mathbf{x}_{t};\mathbf{F}^{\prime}_{t}\mathbf{x}_{t-1}+\mathbf{b}^{\prime}_{t},\mathbf{C}^{\prime}_{t})$.
$\Box$□

##  Appendix E Proof of Propositions 8 and 9

###  E.1 Preliminaries

For some given $N\in\mathbb{N}$, let $\Psi^{n}$ denote either the _Boltzmann
selection function_ (with the convention $h^{0}\coloneqq 0$):

$\displaystyle\Psi^{n}(h^{1:N})$
$\displaystyle\coloneqq\frac{\exp(h^{n})}{1+\sum_{m=0}^{N}\exp(h^{m})},$ (137)

or the _Rosenbluth–Teller selection function_ :

$\displaystyle\Psi^{n}(h^{1:N})$
$\displaystyle\coloneqq\begin{dcases}\frac{\exp(h^{n})}{1+\sum_{m=1}^{N}\exp(h^{m})-1\wedge\exp(h^{n})},&\text{if
$n>0$,}\\\ 1-\sum_{l=1}^{N}\Psi^{l}(h^{1:N}),&\text{if $n=0$.}\end{dcases}$
(138)

In either case, $\Psi^{n}$ is Lipschitz continuous with constant denoted
$[\Psi^{n}]_{\textnormal{{lip}}}$.

###  E.2 Marginal MCMC kernels in the special case: $T=1$

For the moment, we assume that $T=1$. To simplify the notation, we drop the
‘time’ subscripts $t=1$. With this convention, for some bounded and
differentiable $G:\mathbb{R}^{D}\to(0,\infty)$, define

$\displaystyle\pi(\mathbf{x})$
$\displaystyle\propto\operatorname{N}(\mathbf{x};\mathbf{m},\mathbf{C})G(\mathbf{x}).$
(139)

The $\pi$-invariant Markov kernels induced by the (non-auxiliary variable
based) algorithms discussed in this work can then be written as

$\displaystyle P_{a}(\mathrm{d}\tilde{\mathbf{x}}|\mathbf{x})$
$\displaystyle=\sum_{l=0}^{N}\int_{\mathcal{X}^{N+2}}\delta_{\mathbf{x}}(\,\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}|\mathbf{x}^{0})\Psi^{l}(\\{h_{a}^{n}(\mathbf{x}^{0:N})\\}_{n=1}^{N})\delta_{\mathbf{x}^{l}}(\mathrm{d}\tilde{\mathbf{x}}),$
(140)

where have appealed to symmetry to always place the reference ‘path’ in
position $0$, and with

$\displaystyle h_{a}^{n}(\mathbf{x}^{0:N})$ $\displaystyle\coloneqq\log
q_{a}^{-n}(\mathbf{x}^{-n}|\mathbf{x}^{n})-\log
q_{a}^{-0}(\mathbf{x}^{-0}|\mathbf{x}^{0}),$ (141) $\displaystyle
q_{a}^{-n}(\mathbf{x}^{-n}|\mathbf{x}^{n})$
$\displaystyle=\operatorname{N}(\mathbf{x}^{-n};\mathbf{m}_{a}(\mathbf{x}^{n}),\mathbf{C}_{a}),$
(142)

where $\mathbf{m}_{a}(\mathbf{x}^{n})\in\mathbb{R}^{ND}$ is a suitable mean
vector (which may depend on $\mathbf{x}^{n}\in\mathbb{R}^{D}$),
$\mathbf{C}_{a}\in\mathbb{R}^{(ND)\times(ND)}$ a suitable variance, and where
we again slightly abuse notation to let $\mathbf{x}^{-n}$ represent both the
tuple
$(\mathbf{x}^{0},\dotsc,\mathbf{x}^{n-1},\mathbf{x}^{n+1},\dotsc,\mathbf{x}^{N})$
and its vectorised form

$\displaystyle\mathbf{x}^{-n}\coloneqq\operatorname{vec}(\mathbf{x}^{-n})=\begin{bmatrix}\mathbf{x}^{0}\\\
\vdots\\\ \mathbf{x}^{n-1}\\\ \mathbf{x}^{n+1}\\\ \vdots\\\
\mathbf{x}^{N}\end{bmatrix}\in\mathbb{R}^{ND}.$ (143)

Additionally, ‘$a$’ is a placeholder for ‘CSMC’, ‘Particle-MALA’, or
‘Particle-mGRAD’. Specifically, by the developments from Section C (Lemma 3
and its proof), and recalling that the block matrix operator $\mathcal{M}_{N}$
was defined in (94),

$\displaystyle\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}}}(\mathbf{x}^{n})$ $\displaystyle=\bm{1}_{N\times
1}\otimes[\mathbf{m}+\mathbf{A}(\mathbf{x}^{n}+\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})],$
(144)
$\displaystyle\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}}}$
$\displaystyle=\tfrac{\delta}{2}\mathcal{M}_{N}(\mathbf{A},\mathbf{A}^{2})=\tfrac{\delta}{2}[\mathbf{I}_{N}\otimes\mathbf{A}+\bm{1}_{N\times
N}\otimes\mathbf{A}^{2}],$ (145)
$\displaystyle\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}(\mathbf{x}^{n})$ $\displaystyle=\bm{1}_{N\times 1}\otimes\mathbf{m},$
(146)
$\displaystyle\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}$ $\displaystyle=\mathcal{M}_{N}(\mathbf{C},\bm{0}_{D\times
D})=\mathbf{I}_{N}\otimes\mathbf{C},$ (147)
$\displaystyle\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-MALA}}}(\mathbf{x}^{n})$ $\displaystyle=\bm{1}_{N\times
1}\otimes[\mathbf{x}^{n}+\bm{\phi}(\mathbf{x}^{n})+\bm{\varphi}(\mathbf{x}^{n})],$
(148)
$\displaystyle\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-MALA}}}$
$\displaystyle=\tfrac{\delta}{2}\mathcal{M}_{N}(\mathbf{I},\mathbf{I})=\tfrac{\delta}{2}[\mathbf{I}_{ND}+\bm{1}_{N\times
N}\otimes\mathbf{I}],$ (149)

where $\bm{\phi}(\mathbf{x})\coloneqq\kappa\tfrac{\delta}{2}\nabla\log
G(\mathbf{x})$ and
$\bm{\varphi}(\mathbf{x})\coloneqq\kappa\tfrac{\delta}{2}\nabla\log
M(\mathbf{x})$ and with
$\mathbf{A}\coloneqq(\tfrac{\delta}{2}\mathbf{C}^{-1}+\mathbf{I})^{-1}=\mathbf{C}(\mathbf{C}+\tfrac{\delta}{2}\mathbf{I})^{-1}=(\mathbf{C}+\tfrac{\delta}{2}\mathbf{I})^{-1}\mathbf{C}$.

Key to our proofs will be the following bound which follows from the triangle
inequality and a telescoping-sum decomposition (here: $a$ and $b$ are again
placeholders which take values in
$\\{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}}\\}$):

$\displaystyle\lVert
P_{a}(\,\cdot\,|\mathbf{x})-P_{b}(\,\cdot\,|\mathbf{x})\rVert_{\textnormal{{tv}}}$
(150) $\displaystyle\leq\lVert
q_{a}^{-0}(\,\cdot\,|\mathbf{x})-q_{b}^{-0}(\,\cdot\,|\mathbf{x})\rVert_{\textnormal{{tv}}}$
(151)
$\displaystyle\quad+\sum_{l=0}^{N}\int_{\mathcal{X}^{N+1}}\delta_{\mathbf{x}}(\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}|\mathbf{x}^{0})\lvert\Psi^{l}(\\{h_{a}^{n}(\mathbf{x}^{0:N})\\}_{n=1}^{N})-\Psi^{l}(\\{h_{b}^{n}(\mathbf{x}^{0:N})\\}_{n=1}^{N})\rvert$
(152)
$\displaystyle\leq\sqrt{\operatorname{KL}(q_{a}^{-0}(\,\cdot\,|\mathbf{x})\|q_{b}^{-0}(\,\cdot\,|\mathbf{x}))}$
(153)
$\displaystyle\quad+\sum_{l=0}^{N}[\Psi^{l}]_{\textnormal{{lip}}}\int_{\mathcal{X}^{N+1}}\delta_{\mathbf{x}}(\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}|\mathbf{x}^{0})\sum_{n=1}^{N}\lvert
h_{a}^{n}(\mathbf{x}^{0:N})-h_{b}^{n}(\mathbf{x}^{0:N})\rvert$ (154)
$\displaystyle\leq
C\Bigl{[}\sqrt{D_{a,b}^{0}(\mathbf{x})}+\sum_{n=0}^{N}D_{a,b}^{n}(\mathbf{x})\Bigr{]}.$
(155)

Here, the penultimate line follows from Pinsker’s inequality and the Lipschitz
continuity of the selection function; $C\geq 0$ is some constant which may
depend on these Lipschitz constants and $N$ and $D$; for the last inequality,
we have defined

$\displaystyle D_{a,b}^{n}(\mathbf{x})$
$\displaystyle\coloneqq\int_{\mathcal{X}^{N+1}}\delta_{\mathbf{x}}(\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}|\mathbf{x}^{0})\lvert\log
q_{a}^{-n}(\mathbf{x}^{-n}|\mathbf{x}^{n})-\log
q_{b}^{-n}(\mathbf{x}^{-n}|\mathbf{x}^{n})\rvert.$ (156)

###  E.3 Proofs of Part 1

######  Proof (of Part 1 of Proposition 8).

By Assumption A1, the model factorises over time and so do the CSMC and
Particle-mGRAD algorithms. Hence, without loss of generality, we prove the
result in the case that $T=1$ (and we drop the ‘time’ subscript $t=1$
hereafter). Throughout the proof, we will also make repeated use of the fact
that the eigenvalues of $\mathbf{A}_{k}$ are given by
$(2\lambda_{k,d})/(2\lambda_{k,d}+\delta)$, for $d\in[D]$.

For $\varepsilon\geq 1$ the result is trivially true but meaningless. Fix
$\varepsilon\in(0,1)$.

$\displaystyle
F_{k}\coloneqq\bigl{\\{}\mathbf{x}\in\mathbb{R}^{D}\,\big{|}\,\lVert\mathbf{x}-\mathbf{m}\rVert_{2}\leq\lambda_{k}^{(1-\varepsilon)/2}\bigr{\\}},$
(157)

denote a ball of radius $\smash{\lambda_{k}^{(1-\varepsilon)/2}}$ around
$\mathbf{m}$, for any $k\geq 1$. We then have $\pi_{k}(F_{k})=(1+H_{k})^{-1}$,
where, letting $F_{k}^{\mathrm{c}}\coloneqq\mathcal{X}\setminus F_{k}$:

$\displaystyle H_{k}$
$\displaystyle\coloneqq\frac{\int_{F_{k}^{\mathrm{c}}}G(\mathbf{x})\operatorname{N}(\mathrm{d}\mathbf{x};\mathbf{m},\mathbf{C}_{k})}{\int_{F_{k}}G(\mathbf{x})\operatorname{N}(\mathrm{d}\mathbf{x};\mathbf{m},\mathbf{C}_{k})}$
(158)
$\displaystyle\leq\frac{\sup_{x\in\mathcal{X}}G(\mathbf{x})}{\inf_{\mathbf{x}\in
F_{k}}G(\mathbf{x})}\frac{\int_{F_{k}^{\mathrm{c}}}\operatorname{N}(\mathrm{d}\mathbf{x};\mathbf{m},\mathbf{C}_{k})}{\int_{F_{k}}\operatorname{N}(\mathrm{d}\mathbf{x};\mathbf{m},\mathbf{C}_{k})}$
(159)
$\displaystyle\leq\frac{\sup_{x\in\mathcal{X}}G(\mathbf{x})}{\inf_{\mathbf{x}\in
F_{k}}G(\mathbf{x})}\frac{\int_{\mathcal{X}}\operatorname{N}(\mathrm{d}\mathbf{x};\bm{0},\mathbf{I})\operatorname{\mathbb{I}}\\{\lVert\mathbf{x}\rVert_{2}>\lambda_{k}^{-\varepsilon/2}\\}}{\int_{\mathcal{X}}\operatorname{N}(\mathrm{d}\mathbf{x};\bm{0},\mathbf{I})\operatorname{\mathbb{I}}\\{\lVert\mathbf{x}\rVert_{2}\leq\lambda_{k}^{-\varepsilon/2}\\}}$
(160) $\displaystyle\to 0,$ (161)

as $k\to\infty$, where we have used that $G$ is bounded and that
$(\inf_{\mathbf{x}\in F_{k}}G(\mathbf{x}))_{k\geq 1}$ is an increasing
sequence in $(0,\infty)$ (since $(F_{k})_{k\geq 1}$ is decreasing and $F_{1}$
is compact).

By the decomposition from (155), all that remains is to control the terms

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}D_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}^{n}(\mathbf{x}^{0}),$ (162)

for arbitrary $n\in[N]_{0}$.

Firstly, by Lemma 1 from Appendix C, letting
$\lambda(\mathbf{C}_{k})=\\{\lambda_{k,1},\dotsc,\lambda_{k,D}\\}$ denote the
eigenvalues of $\mathbf{C}_{k}$ and noting that $\mathbf{A}_{k}$ is
simultaneously diagonalisable with $\mathbf{A}_{k}^{2}$:

$\displaystyle\lvert\log(\det(\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},k}))-\log(\det(\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}))\rvert$ (163)
$\displaystyle=\Bigl{\lvert}\sum_{d=1}^{D}N\log\lambda_{k,d}-(N-1)\log\Bigl{(}\frac{\delta\lambda_{k,d}}{2\lambda_{k,d}+\delta}\Bigr{)}-\log\left(\frac{\delta\lambda_{k,d}}{2\lambda_{k,d}+\delta}+\frac{2\delta
N\lambda_{k,d}^{2}}{(2\lambda_{k,d}+\delta)^{2}}\right)\Bigr{\rvert}$ (164)
$\displaystyle=\Bigl{\lvert}\sum_{d=1}^{D}N\log\Bigl{(}\frac{2\lambda_{k,d}+\delta}{\delta}\Bigr{)}+\log\Bigl{(}\frac{\delta\lambda_{k,d}}{2\lambda_{k,d}+\delta}\Bigr{)}-\log\left(\frac{\delta\lambda_{k,d}}{2\lambda_{k,d}+\delta}+\frac{2\delta
N\lambda_{k,d}^{2}}{(2\lambda_{k,d}+\delta)^{2}}\right)\Bigr{\rvert}$ (165)
$\displaystyle=\Bigl{\lvert}\sum_{d=1}^{D}N\log\Bigl{(}\frac{2\lambda_{k,d}+\delta}{\delta}\Bigr{)}+\log\Bigl{(}\frac{2\lambda_{k,d}+\delta}{2\lambda_{k,d}(N+1)+\delta}\Bigr{)}\Bigr{\rvert}\in\operatorname{\mathrm{O}}(\lambda_{k}).$
(166)

Secondly, by Lemma 2 from Appendix C,

$\displaystyle\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-
mGRAD}},k}^{-1}=\frac{2}{\delta}\mathcal{M}_{N}(\mathbf{A}_{k}^{-1},-(\mathbf{I}+N\mathbf{A}_{k})^{-1}),$
(167)

and with the conventions that the sum symbol $\sum_{i}$ is shorthand for
$\sum_{i\in[N]_{0}\setminus\\{n\\}}$, that $\sum_{j}$ is shorthand for
$\sum_{j\in[N]_{0}\setminus\\{n\\}}$, that $\sum_{i\neq j}$ is shorthand for
$\sum_{j\in[N]_{0}\setminus\\{n,i\\}}$, and again writing
$\bm{\phi}(\mathbf{x})=\kappa\tfrac{\delta}{2}\nabla\log G(\mathbf{x})$ we
obtain:

$\displaystyle\bigl{\lvert}(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}(\mathbf{x}^{n}))^{\mathrm{T}}\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},k}^{-1}(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}(\mathbf{x}^{n}))$ (168)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!-(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-
mGRAD}},k}(\mathbf{x}^{n}))^{\mathrm{T}}\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-
mGRAD}},k}^{-1}(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}(\mathbf{x}^{n}))\bigr{\rvert}$ (169)
$\displaystyle=\Bigl{\lvert}\sum_{i}(\mathbf{x}^{i}-\mathbf{m})^{\mathrm{T}}\mathbf{C}_{k}^{-1}(\mathbf{x}^{i}-\mathbf{m})$
(170)
$\displaystyle\quad-\frac{2}{\delta}\sum_{i}(\mathbf{x}^{i}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})^{\mathrm{T}}\mathbf{A}_{k}^{-1}(\mathbf{x}^{i}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})$
(171)
$\displaystyle\quad+\frac{2}{\delta}\sum_{i}\sum_{j}(\mathbf{x}^{i}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})^{\mathrm{T}}(\mathbf{I}+N\mathbf{A}_{k})^{-1}(\mathbf{x}^{j}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})\Bigr{\lvert}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(172)
$\displaystyle=\Bigl{\lvert}\frac{2}{\delta}\sum_{i}\sum_{j}(\mathbf{x}^{i}-\mathbf{m})^{\mathrm{T}}(\mathbf{I}+N\mathbf{A}_{k})^{-1}(\mathbf{x}^{j}-\mathbf{m})$
(173)
$\displaystyle\quad+\frac{2}{\delta}\sum_{i}(\mathbf{x}^{i}-\mathbf{m})^{\mathrm{T}}\Bigl{(}\frac{\delta}{2}\mathbf{C}_{k}^{-1}-\mathbf{A}_{k}^{-1}\Bigr{)}(\mathbf{x}^{i}-\mathbf{m})$
(174)
$\displaystyle\quad+\frac{2}{\delta}\sum_{i}(\mathbf{x}^{i}-\mathbf{m})^{\mathrm{T}}\Bigl{[}\mathbf{I}+\Bigl{(}\frac{4(N-1)}{\delta}-1\Bigr{)}(\mathbf{I}+N\mathbf{A}_{k})^{-1}\mathbf{A}_{k}\Bigr{]}(\mathbf{x}^{n}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})$
(175)
$\displaystyle\quad+N(\mathbf{x}^{n}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})^{\mathrm{T}}\Bigl{[}\mathbf{A}_{k}+\Bigl{(}\frac{2(N-1)}{\delta}-1\Bigr{)}\mathbf{A}_{k}(\mathbf{I}+N\mathbf{A}_{k})^{-1}\mathbf{A}_{k}\Bigr{]}(\mathbf{x}^{n}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m})\Bigr{\lvert}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(176)
$\displaystyle\leq\frac{2}{\delta}\sum_{i}\sum_{j}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\lVert\mathbf{x}^{j}-\mathbf{m}\rVert_{2}\lVert(\mathbf{I}+N\mathbf{A}_{k})^{-1}\rVert_{2,2}$
(177)
$\displaystyle\quad+\frac{2}{\delta}\sum_{i}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\Bigl{\lVert}\frac{\delta}{2}\mathbf{C}_{k}^{-1}-\mathbf{A}_{k}^{-1}\Bigr{\rVert}_{2,2}$
(178)
$\displaystyle\quad+\frac{2}{\delta}\sum_{i}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\lVert\mathbf{x}^{n}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m}\rVert_{2}\Bigl{\lVert}\mathbf{I}+\Bigl{(}\frac{4(N-1)}{\delta}-1\Bigr{)}(\mathbf{I}+N\mathbf{A}_{k})^{-1}\mathbf{A}_{k}\Bigr{\rVert}_{2,2}\\!\\!\\!\\!\\!\\!\\!\\!$
(179)
$\displaystyle\quad+N\lVert\mathbf{x}^{n}-\bm{\phi}(\mathbf{x}^{n})-\mathbf{m}\rVert_{2}^{2}\Bigl{\lVert}\mathbf{A}_{k}+\Bigl{(}\frac{2(N-1)}{\delta}-1\Bigr{)}\mathbf{A}_{k}(\mathbf{I}+N\mathbf{A}_{k})^{-1}\mathbf{A}_{k}\Bigr{\rVert}_{2,2}$
(180) $\displaystyle\leq
C\Bigl{[}\sum_{i}\sum_{j}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\lVert\mathbf{x}^{j}-\mathbf{m}\rVert_{2}$
(181)
$\displaystyle\qquad\quad+(1+\lVert\mathbf{x}^{n}-\mathbf{m}\rVert_{2})\sum_{i}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}$
(182)
$\displaystyle\qquad\quad+\lambda_{k}(1+\lVert\mathbf{x}^{n}-\mathbf{m}\rVert_{2})^{2}\Bigr{]},$
(183)

for some constant $C\geq 0$ which only depends on $N$, $\delta$ and
$\mathbf{m}$. Here, we have used that all the matrices inside the operator
norms are simultaneously diagonalisable with $\mathbf{C}_{k}$ (so that the
operator norms can be bounded above by some function of $\lambda_{k}$):

$\displaystyle\lVert(\mathbf{I}+N\mathbf{A}_{k})^{-1}\rVert_{2,2}$
$\displaystyle\leq 1,$ (184)
$\displaystyle\Bigl{\lVert}\frac{\delta}{2}\mathbf{C}_{k}^{-1}-\mathbf{A}_{k}^{-1}\Bigr{\rVert}_{2,2}$
$\displaystyle=1,$ (185)
$\displaystyle\Bigl{\lVert}\mathbf{A}_{k}+\Bigl{(}\frac{2(N-1)}{\delta}-1\Bigr{)}\mathbf{A}_{k}(\mathbf{I}+N\mathbf{A}_{k})^{-1}\mathbf{A}_{k}\Bigr{\rVert}_{2,2}$
$\displaystyle\leq C^{\prime}\frac{2\lambda_{k}}{2\lambda_{k}+\delta}\leq
C^{\prime}\frac{2}{\delta}\lambda_{k},$ (186)
$\displaystyle\Bigl{\lVert}\mathbf{I}+\Bigl{(}\frac{4(N-1)}{\delta}-1\Bigr{)}(\mathbf{I}+N\mathbf{A}_{k})^{-1}\mathbf{A}_{k}\Bigr{\rVert}_{2,2}$
$\displaystyle\leq C^{\prime\prime},$ (187)

for other constants $C^{\prime},C^{\prime\prime}\geq 0$.

Furthermore, by definition of $(F_{k})_{k\geq 1}$,

$\displaystyle\sup_{\mathbf{x}\in
F_{k}}\lVert\mathbf{x}-\mathbf{m}\rVert_{2}\in\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)/2}).$
(188)

Consequently, for $i,j\in[N]_{0}$:

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}\int_{\mathcal{X}^{N}}\operatorname{N}(\mathrm{d}\mathbf{x}^{-0};\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}},\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},k})\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\lVert\mathbf{x}^{j}-\mathbf{m}\rVert_{2}$
(189)
$\displaystyle\in\begin{cases}\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)}),&\text{if
$i=j=0$,}\\\
\operatorname{\mathrm{O}}(\lambda_{k}^{(2-\varepsilon)/2}),&\text{if either
$i=0$ or $j=0$,}\\\ \operatorname{\mathrm{O}}(\lambda_{k}),&\text{if neither
$i=0$ nor $j=0$,}\\\ \end{cases}$ (190)

as $\lambda_{k}\to 0$, and where the last two cases follow from the
Cauchy–Schwarz inequality. Similarly, for $i\in[N]_{0}$,

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}\int_{\mathcal{X}^{N}}\operatorname{N}(\mathrm{d}\mathbf{x}^{-0};\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}},\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},k})\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}$
$\displaystyle\in\begin{cases}\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)/2}),&\text{if
$i=0$,}\\\ \operatorname{\mathrm{O}}(\lambda_{k}^{1/2}),&\text{if $i\neq
0$,}\end{cases}$ (191)

as $\lambda_{k}\to 0$. Combining the bounds from (166)–(191) then shows that

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}D_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-
mGRAD}},k}^{n}(\mathbf{x}^{0})\in\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)/2}),$
(192)

for any $n\in[N]_{0}$. Plugging these bounds into (155) completes the proof.
$\Box$□

######  Proof (of Part 1 of Proposition 9).

By Assumption A1, the model factorises over time and so do the Particle-MALA
and Particle-mGRAD algorithms. Hence, without loss of generality, we again
only prove the result in the case that $T=1$ (and we again drop the ‘time’
subscript $t=1$ hereafter).

For $\varepsilon\geq 1$ the result is trivially true but meaningless. Fix
$\varepsilon\in(0,1)$. Since $G$ is integrable (by Assumption A3) and since
$\pi_{k}$ is invariant to scaling of $G$ by a positive constant factor, we
assume that $\int_{\mathcal{X}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}$ = 1,
without loss of generality, so that $G$ can be viewed as a density (and we
will also use the symbol $G$ to denote the corresponding distribution). Let
$\mathbf{m}_{G}$ and $\mathbf{C}_{G}$ be mean and variance of $G$ (which exist
by Assumption A3) and define

$\displaystyle F_{k}$
$\displaystyle\coloneqq\bigl{\\{}\mathbf{x}\in\mathcal{X}\,\big{|}\,\lVert\mathbf{x}-\mathbf{m}\rVert_{2}\vee\sqrt{(\mathbf{x}-\mathbf{m}_{G})^{\mathrm{T}}\mathbf{C}_{G}^{-1}(\mathbf{x}-\mathbf{m}_{G})}<\lambda_{k}^{\varepsilon/2}\bigr{\\}}.$
(193)

We then have $\pi_{k}(F_{k})=(1+H_{k})^{-1}$, where, letting $\mathbf{Y}\sim
G$ and letting $F_{k}^{\mathrm{c}}\coloneqq\mathcal{X}\setminus F_{k}$:

$\displaystyle H_{k}$
$\displaystyle\coloneqq\frac{\int_{F_{k}^{\mathrm{c}}}G(\mathbf{x})\operatorname{N}(\mathrm{d}\mathbf{x};\mathbf{m},\mathbf{C}_{k})}{\int_{F_{k}}G(\mathbf{x})\operatorname{N}(\mathrm{d}\mathbf{x};\mathbf{m},\mathbf{C}_{k})}$
(194)
$\displaystyle\leq\frac{\int_{F_{k}^{\mathrm{c}}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}}{\inf_{\mathbf{x}\in
F_{k}}\exp(-\tfrac{1}{2}\lVert\mathbf{x}-\mathbf{m}\rVert_{2}^{2}\lambda_{k}^{-1})\int_{F_{k}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}}$
(195)
$\displaystyle\leq\frac{\int_{F_{k}^{\mathrm{c}}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}}{\int_{F_{k}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}}\exp(\tfrac{1}{2}\lambda_{k}^{\varepsilon-1})$
(196) $\displaystyle=\operatorname{\mathbb{P}}(\mathbf{Y}\in
F_{k}^{\mathrm{c}})\frac{\exp(\tfrac{1}{2}\lambda_{k}^{\varepsilon-1})}{\int_{F_{k}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}}$
(197)
$\displaystyle\leq\operatorname{\mathbb{P}}\bigl{(}\sqrt{(\mathbf{Y}-\mathbf{m}_{G})^{\mathrm{T}}\mathbf{C}_{G}^{-1}(\mathbf{Y}-\mathbf{m}_{G})}\geq\lambda_{k}^{\varepsilon/2}\bigr{)}\frac{\exp(\tfrac{1}{2}\lambda_{k}^{\varepsilon-1})}{\int_{F_{k}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}}$
(198)
$\displaystyle\leq\frac{D}{\lambda_{k}^{\varepsilon}}\frac{\exp(\tfrac{1}{2}\lambda_{k}^{\varepsilon-1})}{\int_{F_{k}}G(\mathbf{x})\,\mathrm{d}\mathbf{x}}$
(199) $\displaystyle\to 0.$ (200)

The penultimate line follows from the (multidimensional) Chebyshev’s
inequality and the last line uses that $F_{k}\to\mathcal{X}$ as $k\to\infty$.

By the decomposition from (155), all that remains is to control the terms

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}D_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}^{n}(\mathbf{x}^{0}),$ (201)

for arbitrary $n\in[N]_{0}$

Firstly, by Lemma 1 from Appendix C, letting
$\lambda(\mathbf{C}_{k})=\\{\lambda_{k,1},\dotsc,\lambda_{k,D}\\}$ denote the
eigenvalues of $\mathbf{C}_{k}$ and noting that $\mathbf{A}_{k}$ is
simultaneously diagonalisable with $\mathbf{A}_{k}^{2}$:

$\displaystyle\lvert\log(\det(\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}}}))-\log(\det(\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}))\rvert$ (202)
$\displaystyle=\Bigl{\lvert}\sum_{d=1}^{D}N\log\Bigl{(}\frac{2\lambda_{k,d}+\delta}{2\lambda_{k,d}}\Bigr{)}+\log\Bigl{(}\frac{2\lambda_{k,d}+\delta}{2\lambda_{k,d}+\delta/(N+1)}\Bigr{)}\Bigr{\rvert}\in\operatorname{\mathrm{O}}(\lambda_{k}^{-1}).$
(203)

Secondly, by Lemma 2 from Appendix C, and again with the conventions that
$\sum_{i}$ is shorthand for $\sum_{i\in[N]_{0}\setminus\\{n\\}}$, that
$\sum_{j}$ is shorthand for $\sum_{j\in[N]_{0}\setminus\\{n\\}}$, that
$\sum_{i\neq j}$ is shorthand for $\sum_{j\in[N]_{0}\setminus\\{n,i\\}}$, and
writing $\bm{\phi}(\mathbf{x})\coloneqq\kappa\tfrac{\delta}{2}\nabla\log
G(\mathbf{x})$ as well as
$\bm{\varphi}_{k}(\mathbf{x})\coloneqq\kappa\tfrac{\delta}{2}\nabla\log
M_{k}(\mathbf{x})=\kappa\tfrac{\delta}{2}\mathbf{C}_{k}^{-1}(\mathbf{m}-\mathbf{x})$,
so that
$\bm{\phi}(\mathbf{x})+\bm{\varphi}_{k}(\mathbf{x})=\kappa\tfrac{\delta}{2}\nabla\log\pi_{k}(\mathbf{x})$:

$\displaystyle\bigl{\lvert}(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},k}(\mathbf{x}^{n}))^{\mathrm{T}}\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}}}^{-1}(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-MALA}},k}(\mathbf{x}^{n})$ (204)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!-(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-
mGRAD}},k}(\mathbf{x}^{n}))^{\mathrm{T}}\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-
mGRAD}},k}^{-1}(\mathbf{x}^{-n}-\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}(\mathbf{x}^{n}))\bigr{\rvert}$ (205) $\displaystyle\leq
C\lambda_{k}^{-1}\Bigl{[}\sum_{i}\sum_{j}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\lVert\mathbf{x}^{j}-\mathbf{m}\rVert_{2}+(1+\lVert\mathbf{x}^{n}-\mathbf{m}\rVert_{2})\sum_{i=0}^{N}\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\Bigr{]},$
(206)

for some constant $C\geq 0$ which only depends on $N$, $\delta$ and
$\mathbf{m}$. Here, we have followed the same steps as for (183) and used that
all the matrices inside the operator norms are simultaneously diagonalisable
with $\mathbf{C}_{k}$ (so that the operator norms can be bounded above by some
function of $\lambda_{k}^{-1}$).

Furthermore, by definition of $F_{k}$, we have

$\displaystyle\sup_{\mathbf{x}\in
F_{k}}\lVert\mathbf{x}-\mathbf{m}\rVert_{2}\in\operatorname{\mathrm{O}}(\lambda_{k}^{\varepsilon/2}),$
(207)

as $\lambda_{k}\to\infty$. Consequently, for $i,j\in[N]_{0}$, by the
Cauchy–Schwarz inequality:

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}\int_{\mathcal{X}^{N}}\operatorname{N}(\mathrm{d}\mathbf{x}^{-0};\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},k},\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}}})\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}\lVert\mathbf{x}^{j}-\mathbf{m}\rVert_{2}$
(208)
$\displaystyle\in\begin{cases}\operatorname{\mathrm{O}}(\lambda_{k}^{\varepsilon}),&\text{if
$i=j=0$,}\\\ \operatorname{\mathrm{O}}(\lambda_{k}^{\varepsilon/2}),&\text{if
either $i=0$ or $j=0$,}\\\ \operatorname{\mathrm{O}}(1),&\text{if neither
$i=0$ nor $j=0$,}\\\ \end{cases}$ (209)

as $\lambda_{k}\to\infty$. Similarly, for $i\in[N]_{0}$,

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}\int_{\mathcal{X}^{N}}\operatorname{N}(\mathrm{d}\mathbf{x}^{-0};\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},k},\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-MALA}}})\lVert\mathbf{x}^{i}-\mathbf{m}\rVert_{2}$
$\displaystyle\in\begin{cases}\operatorname{\mathrm{O}}(\lambda_{k}^{\varepsilon/2}),&\text{if
$i=0$,}\\\ \operatorname{\mathrm{O}}(1),&\text{if $i\neq 0$,}\end{cases}$
(210)

as $\lambda_{k}\to\infty$. Combining the bounds from (203)–(210) then shows
that

$\displaystyle\sup_{\mathbf{x}^{0}\in
F_{k}}D_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-
mGRAD}},k}^{n}(\mathbf{x}^{0})\in\operatorname{\mathrm{O}}(\lambda_{k}^{-(1-\varepsilon)/2}),$
(211)

for any $n\in[N]_{0}$. Plugging these bounds into (155) completes the proof.
$\Box$□

###  E.4 Auxiliary MCMC kernels in the special case: $T=1$

The $\pi$-invariant Markov kernels induced by the auxiliary-variable based
algorithms discussed in this work can then be written as

$\displaystyle P_{a}(\mathrm{d}\tilde{\mathbf{x}}|\mathbf{x})$
$\displaystyle=\sum_{l=0}^{N}\int_{\mathcal{X}^{N+3}}\delta_{\mathbf{x}}(\,\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}\times\mathrm{d}\mathbf{u}|\mathbf{x}^{-0},\mathbf{x}^{0})\Psi^{l}(\\{h_{a}^{n}(\mathbf{x}^{0:N},\mathbf{u})\\}_{n=1}^{N})\delta_{\mathbf{x}^{l}}(\mathrm{d}\tilde{\mathbf{x}}),$
(212)

where have appealed to symmetry to always place the reference ‘path’ in
position $0$, where ‘$a$’ is now a placeholder for ‘Particle-aGRAD’,
‘Particle-aMALA’, or ‘CSMC’ and with

$\displaystyle q_{a}^{-n}(\mathbf{x}^{-n},\mathbf{u}|\mathbf{x}^{n})$
$\displaystyle=q_{a}^{-n}(\mathbf{x}^{-n}|\mathbf{x}^{n})q_{a}^{-n}(\mathbf{u}|\mathbf{x}^{-n},\mathbf{x}^{n}),$
(213) $\displaystyle h_{a}^{n}(\mathbf{x}^{0:N},\mathbf{u})$
$\displaystyle\coloneqq\log q_{a}^{-n}(\mathbf{x}^{-n}|\mathbf{x}^{n})-\log
q_{a}^{-0}(\mathbf{x}^{-0}|\mathbf{x}^{0})$ (214)
$\displaystyle\quad+\operatorname{\mathbb{I}}\\{a\neq\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}\\}[\log q_{a}^{-n}(\mathbf{u}|\mathbf{x}^{-n},\mathbf{x}^{n})-\log
q_{a}^{-0}(\mathbf{u}|\mathbf{x}^{-0},\mathbf{x}^{0})],$ (215) $\displaystyle
q_{a}^{-n}(\mathbf{x}^{-n}|\mathbf{x}^{n})$
$\displaystyle=\operatorname{N}(\mathbf{x}^{-n};\mathbf{m}_{a}(\mathbf{x}^{n}),\mathbf{C}_{a}),$
(216) $\displaystyle q_{a}^{-n}(\mathbf{u}|\mathbf{x}^{-n},\mathbf{x}^{n})$
$\displaystyle=\operatorname{N}(\mathbf{u};\bm{\nu}_{a}(\mathbf{x}^{-n},\mathbf{x}^{n}),\bm{\Upsilon}_{a}),$
(217)

where again $\mathbf{m}_{a}(\mathbf{x}^{n})\in\mathbb{R}^{ND}$ and
$\bm{\nu}_{a}(\mathbf{x}^{-n},\mathbf{x}^{n})\in\mathbb{R}^{D}$ are suitable
mean vector, and $\mathbf{C}_{a}\in\mathbb{R}^{(ND)\times(ND)}$,
$\bm{\Upsilon}_{a}\in\mathbb{R}^{D\times D}$ are suitable covariance variance
matrices, and we again write
$\mathbf{x}^{-n}\coloneqq\operatorname{vec}(\mathbf{x}^{-n})$. Specifically,

$\displaystyle\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}(\mathbf{x}^{n})$
$\displaystyle=\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}}}(\mathbf{x}^{n})$ (218)
$\displaystyle\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}$
$\displaystyle=\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}}},$ (219)
$\displaystyle\bm{\nu}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}(\mathbf{x}^{-n},\mathbf{x}^{n})$
$\displaystyle=(\mathbf{I}+N\mathbf{A})^{-1}[(N+1)\bar{\mathbf{x}}+\bm{\phi}(\mathbf{x}^{n})+N(\mathbf{A}-\mathbf{I})\mathbf{m}],$
(220)
$\displaystyle\bm{\Upsilon}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}$
$\displaystyle=\tfrac{\delta}{2}(\mathbf{I}+N\mathbf{A})^{-1},$ (221)
$\displaystyle\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-aMALA}}}(\mathbf{x}^{n})$
$\displaystyle=\mathbf{m}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-aMALA}}}(\mathbf{x}^{n}),$ (222)
$\displaystyle\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-MALA}}}$
$\displaystyle=\mathbf{C}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-MALA}}},$ (223)
$\displaystyle\bm{\nu}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-aMALA}}}(\mathbf{x}^{-n},\mathbf{x}^{n})$
$\displaystyle=\bar{\mathbf{x}}+\tfrac{1}{N+1}\bm{\phi}(\mathbf{x}^{n}),$
(224)
$\displaystyle\bm{\Upsilon}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-aMALA}}}$ $\displaystyle=\tfrac{\delta}{2(N+1)}\mathbf{I},$ (225)

by Part 2 of Lemma 3 and Lemma 2 from Appendix C. Of course, the standard CSMC
algorithm does not make use of the auxiliary variable $\mathbf{u}$, so we
extend the space to include $\mathbf{u}$ with

$\displaystyle\bm{\nu}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}(\mathbf{x}^{-n},\mathbf{x}^{n})$
$\displaystyle=\bm{\nu}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}(\mathbf{x}^{-n},\mathbf{x}^{n}),$ (226)
$\displaystyle\bm{\Upsilon}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}$
$\displaystyle=\bm{\Upsilon}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}.$ (227)

Key to our proofs will be the following bound which follows by the triangle
inequality and a telescoping-sum decomposition (here: $a$ is again a
placeholder which takes values in
$\\{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-MALA}}\\}$ whilst we will always set
$b=\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode Particle-
aGRAD}}$ and $q_{a}^{-m}(\,\cdot\,|\mathbf{x})$; and, unless otherwise stated,
$q_{b}^{-m}(\,\cdot\,|\mathbf{x})$ denote the joint distributions on the space
that includes the auxiliary variable $\mathbf{u}$):

$\displaystyle\lVert
P_{a}(\,\cdot\,|\mathbf{x})-P_{b}(\,\cdot\,|\mathbf{x})\rVert_{\textnormal{{tv}}}$
$\displaystyle\leq\lVert
q_{a}^{-0}(\,\cdot\,|\mathbf{x})-q_{b}^{-0}(\,\cdot\,|\mathbf{x})\rVert_{\textnormal{{tv}}}$
$\displaystyle\quad+\sum_{l=0}^{N}\int_{\mathcal{X}^{N+2}}\delta_{\mathbf{x}}(\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}\times\mathrm{d}\mathbf{u}|\mathbf{x}^{0})$
$\displaystyle\qquad\qquad\times\lvert\Psi^{l}(\\{h_{a}^{n}(\mathbf{x}^{0:N},\mathbf{u})\\}_{n=1}^{N})-\Psi^{l}(\\{h_{b}^{n}(\mathbf{x}^{0:N},\mathbf{u})\\}_{n=1}^{N})\rvert$
$\displaystyle\leq\sqrt{\operatorname{KL}(q_{a}^{-0}(\,\cdot\,|\mathbf{x})\|q_{b}^{-0}(\,\cdot\,|\mathbf{x}))}$
$\displaystyle\quad+\sum_{l=0}^{N}[\Psi^{l}]_{\textnormal{{lip}}}\int_{\mathcal{X}^{N+2}}\delta_{\mathbf{x}}(\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}\times\mathrm{d}\mathbf{u}|\mathbf{x}^{0})\sum_{n=1}^{N}\lvert
h_{a}^{n}(\mathbf{x}^{0:N},\mathbf{u})-h_{b}^{n}(\mathbf{x}^{0:N},\mathbf{u})\rvert$
$\displaystyle\leq
C\Bigl{[}\smash{\sqrt{D_{a,b}^{0}(\mathbf{x})+E_{a,b}^{0}(\mathbf{x})}}+\sum_{n=0}^{N}D_{a,b}^{n}(\mathbf{x})+\widetilde{E}_{a,b}^{n}(\mathbf{x})\Bigr{]}.$

Here, the penultimate line follows from Pinsker’s inequality and the Lipschitz
continuity of the selection function; $C\geq 0$ is some constant which may
depend on these Lipschitz constants and on $N$ and $D$;
$D_{a,b}^{n}(\mathbf{x})$ is defined exactly as in the marginal case (156).
Furthermore, we have defined

$\displaystyle E_{a,b}^{n}(\mathbf{x})$
$\displaystyle\coloneqq\int_{\mathcal{X}^{N+1}}\delta_{\mathbf{x}}(\mathrm{d}\mathbf{x}^{0})q_{a}^{-0}(\mathrm{d}\mathbf{x}^{-0}\times\mathrm{d}\mathbf{u}|\mathbf{x}^{0})\lvert\log
q_{a}^{-n}(\mathbf{u}|\mathbf{x}^{-n},\mathbf{x}^{n})-\log
q_{b}^{-n}(\mathbf{u}|\mathbf{x}^{-n},\mathbf{x}^{n})\rvert.$ (228)

Finally, if $a\neq\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}$ and $b\neq\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}$, we we have defined

$\displaystyle\widetilde{E}_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-
aMALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}^{n}(\mathbf{x})$ $\displaystyle\coloneqq
E_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode Particle-
aMALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}^{n}(\mathbf{x}),$ (229)

whilst

$\displaystyle\widetilde{E}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}^{n}(\mathbf{x})$ (230)
$\displaystyle\coloneqq\int_{\mathcal{X}^{N+2}}\delta_{\mathbf{x}}(\mathrm{d}\mathbf{x}^{0})q_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}}}^{-0}(\mathrm{d}\mathbf{x}^{-0}\times\mathrm{d}\mathbf{u}|\mathbf{x}^{0})$
(231) $\displaystyle\qquad\times\lvert\log
q_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode Particle-
aGRAD}}}^{-n}(\mathbf{u}|\mathbf{x}^{-n},\mathbf{x}^{n})-\log
q_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode Particle-
aGRAD}}}^{-0}(\mathbf{u}|\mathbf{x}^{-0},\mathbf{x}^{0})\rvert.$ (232)

###  E.5 Proofs of Part 2

######  Proof (of Part 2 of Proposition 8).

Assume the same setting as in the proof of Part 1 of Proposition 8 with the
same definition of $F_{k}$.

We proceed by controlling the terms in (E.4). Note that
$D_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-
aGRAD}},k}^{n}=D_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}^{n}$. Hence, by the arguments from the proof of Part 1 of
Proposition 8,

$\displaystyle\sup_{\mathbf{x}\in
F_{k}}D_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-
aGRAD}},k}^{n}(\mathbf{x})\in\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)/2}).$
(233)

Additionally, due to (226)–(227),
$E_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-aGRAD}}}^{n}(\mathbf{x})=0$. Finally, using similar arguments as in
the proofs for the ‘marginal’ algorithm, we can verify that

$\displaystyle\sup_{\mathbf{x}\in
F_{k}}\widetilde{E}_{\text{\lx@glossaries@gls@link{acronym}{CSMC}{\leavevmode
CSMC}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-
aGRAD}},k}^{n}(\mathbf{x})\in\operatorname{\mathrm{O}}(\lambda_{k}^{(1-\varepsilon)/2}).$
(234)

This completes the proof. $\Box$□

######  Proof (of Part 2 of Proposition 9).

Assume the same setting as in the proof ofPart 1 of Proposition 9 with the
same definition of $F_{k}$.

We proceed by controlling the terms in (E.4). Note that
$D_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-
aMALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-
aGRAD}},k}^{n}=D_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEMALA}{\leavevmode
Particle-
MALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEMGRAD}{\leavevmode
Particle-mGRAD}},k}^{n}$. Hence, by the arguments from the proof Part 1 of
Proposition 9,

$\displaystyle\sup_{\mathbf{x}\in
F_{k}}D_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-
aMALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-
aGRAD}},k}^{n}(\mathbf{x})\in\operatorname{\mathrm{O}}(\lambda_{k}^{-(1-\varepsilon)/2}).$
(235)

Additionally, using similar arguments as in the proofs for the ‘marginal’
algorithm, we can verify that

$\displaystyle\sup_{\mathbf{x}\in
F_{k}}E_{\text{\lx@glossaries@gls@link{acronym}{PARTICLEAMALA}{\leavevmode
Particle-
aMALA}},\text{\lx@glossaries@gls@link{acronym}{PARTICLEAGRAD}{\leavevmode
Particle-
aGRAD}},k}^{n}(\mathbf{x})\in\operatorname{\mathrm{O}}(\lambda_{k}^{-(1-\varepsilon)}).$
(236)

This completes the proof. $\Box$□

##  Appendix F Step-size adaptation

All our algorithms involve the calibration of several step sizes $\delta_{t}$,
one for each time step. To calibrate these, we implement a routine that
recursively increases or decreases $\delta_{t}$ if the running average of the
acceptance rate $\alpha_{t}$ (i.e., the relative frequency with which
$\mathbf{x}_{t}$ is updated) is respectively above or below a pre-specified
target acceptance rate (in our experiments, we picked this to be
$\alpha^{*}=$75\text{\,}\mathrm{\char 37\relax}\text{/}$$). The only exception
to this lies in the twisted algorithms of Section 4.4 which we calibrate using
a single step-size $\delta$ (so that $\delta=\delta_{1}=\dotsc=\delta_{T}$),
and for which the target relates to the overall acceptance rate averaged
across time steps. The reason for this difference stems from the fact that the
twisting causes the acceptance rate at time $s$ additionally depend on
_future_ auxiliary variables $\mathbf{u}_{t}$, and therefore the future step-
size parameters $\delta_{t}$ (for $t>s$), thereby making the behaviour per
time-step harder to control. In practice, our calibration of the twisted
Particle-aGRAD is therefore more similar to that of aGRAD than that of our
other algorithms. The adaptation procedure is summarised in the following
algorithm.

[step-size adaptation]

  1. 1.

Initialise the trajectory $\mathbf{x}_{1:T}[0]$, the initial step sizes
$\delta_{t}[0]$ (for $t\in[T]$), the initial learning rate
$\rho[0]=\frac{1}{2}$.

  2. 2.

Initialise the history of accepted time steps
$A\coloneqq(A_{w,t})\in\\{0,1\\}^{W\times T}$, with $0$ everywhere.

  3. 3.

For $k=1,\dotsc,K$,

    1. (a)

sample $\mathbf{x}_{1:T}[k]\sim P(\,\cdot\,|\mathbf{x}_{1:T}[k-1])$, where $P$
denotes the Markov kernel induced by one of the algorithms discussed in this
work with step sizes $\delta_{1:T}$ set equal to $\delta_{1:T}[k-1]$,

    2. (b)

roll the array $A$ by one: set $A_{2:\min\\{W,k\\},t}\coloneqq
A_{1:\min\\{W-1,k-1\\},t}$, and
$A_{1,t}=\operatorname{\mathbb{I}}\\{\mathbf{x}_{t}[k]=\mathbf{x}_{t}[k-1]\\}$,
for $t\in[T]$,

    3. (c)

compute $\alpha_{t}=\frac{1}{\min\\{W,k\\}}\sum_{w=1}^{\min\\{W,k\\}}A_{w,t}$,
for $t\in[T]$,

    4. (d)

if $\lvert\alpha_{t}-\alpha^{*}\rvert<\sigma$ then keep
$\delta_{t}[k]=\delta_{t}[k-1]$ unchanged;
otherwise, set

$\displaystyle\delta_{t}[k]\coloneqq\delta_{t}[k-1]+\max\\{k^{\gamma}\rho,\rho_{\min}\\}(\alpha_{t}-\alpha)/\alpha^{*}.$
(237)

In our experiments, we took $\sigma=$5\text{\,}\mathrm{\char
37\relax}\text{/}$$, $K=$10\,000$$, $\delta_{t}[0]=10^{-2}$, $W=$100$$,
$\rho=\frac{1}{2}$, $\rho_{\min}=10^{-3}$, $\gamma=-\frac{1}{2}$.

##  Appendix G Additional experimental results

In this section, we provide additional simulation results for the multivariate
stochastic volatility model experiments from Section 5.

###  G.1 Calibrated step sizes and acceptance rates

Recall that the step sizes $\delta_{t}$ were calibrated to achieve an
acceptance rate of $75\text{\,}\mathrm{\char 37\relax}\text{/}$. Here, the
‘acceptance rate’ at time $t$ refers to the relative frequency with which the
state $\mathbf{x}_{t}$ is updated. The calibrated step sizes are shown in
Figure 6; and the corresponding acceptance rates are shown in Figure 7.

The results are averaged over the four chains and five simulated data sets. We
do not report CSMC as it does not require calibration. All methods
consistently resulted in acceptance rates close to the target
$75\text{\,}\mathrm{\char 37\relax}\text{/}$. Only the twisted Particle-aGRAD
algorithm showed more instability as the informativeness of the prior
decreased: this is because, contrary to the methods, only a single step-size
is used for all time steps, so calibrating for the informativeness of
individual observations is not feasible. This seems to hint to the fact that
the twisted Particle-aGRAD, _under our proposed calibration,_ is less robust
than alternatives to heterogeneous levels of informativeness.

Figure 6: Adaptation of the step-size parameters $\delta_{t}$, averaged across
all four chains and all five simulated data sets (per value of $\tau$) in the
multivariate stochastic volatility model. Figure 7: Acceptance rates (i.e.,
relative frequencies with which states are updated), averaged across all four
chains and all five simulated data sets (per value of $\tau$) in the
multivariate stochastic volatility model.

###  G.2 Breakdown of CSMC, aMALA and MALA

In this section, we illustrate the breakdown of CSMC, aMALA and MALA.

Firstly, Figure 8 illustrates that the estimates of the marginal posterior
means of $x_{t,15}$ (the $15$th component of the state at time $t$) produced
by CSMC, aMALA and MALA differ substantially from those produced by all the
other algorithms. We emphasise that the $15$th component was arbitrarily
chosen as an example and is representative of the other components.

Figure 8: Estimated posterior mean of $x_{t,15}$ minus the estimated posterior
mean of $x_{t,15}$ under the aGRAD algorithm, averaged across all four chains
for each of the five simulated data sets (per value of $\tau$) in the
multivariate stochastic volatility model. The figure shows that the estimated
posterior means of CSMC, aMALA and MALA differ substantially from those of all
the other algorithms.

Secondly, Figure 9 illustrates that the energy traces of CSMC, aMALA and MALA
differ substantially from those of all the other algorithms. Here, the energy
is defined as $\log\pi_{T}(\mathbf{x}_{1:T}[i])$, where $\mathbf{x}_{1:T}[i]$
is the sample from the $i$th iteration after burn-in. Such energy traces serve
as a visual illustration of both stationarity and mixing speed: if the energy
trace of a sampler differs too much from the others, or is not consistent
across the independent Markov chains we used, the sampler is unlikely to
perform correctly.

Figure 9: Energy (i.e., $\log\pi_{T}(\mathbf{x}_{1:T}[i])+\mathrm{const}$,
where $\mathbf{x}_{1:T}[i]$ is the sample from the $i$th iteration after burn-
in) minus the energy under the aGRAD algorithm, averaged across all four
chains for each of the five simulated data sets (per value of $\tau$) in the
multivariate stochastic volatility model. The figure shows that the energy
traces of CSMC, aMALA and MALA differ substantially from those of all the
other algorithms.

###  G.3 Effective sample sizes

In this section, in Figures 10–12 report the minimum, median and maximum ESS
and ESS per second (averaged across all four chains and all five simulated
data sets) individually for each time step $t=1,\dotsc,T$.

Figure 10: Minimum ESS and ESS per second averaged across all four chains and
all five simulated data sets (per value of $\tau$) in the multivariate
stochastic volatility model. Figure 11: Medium ESS and ESS per second averaged
across all four chains and all five simulated data sets (per value of $\tau$)
in the multivariate stochastic volatility model. Figure 12: Maximum ESS and
ESS per second averaged across all four chains and all five simulated data
sets (per value of $\tau$) in the multivariate stochastic volatility model.

###  G.4 Autocorrelation

Figure 13 shows the autocorrelation (corrected using Vehtari et al., , 2021)
of the energy from Figure 9. This serves as a visual confirmation of the
statistical performance of the different algorithms considered under several
prior dispersion regimes: as expected, the twisted Particle-aGRAD dominates
all other alternatives, while Particle-aMALA+ dominates other alternatives,
including aGRAD as soon as the prior variance is large enough, followed by
Particle-aGRAD/Particle-aGRAD, and then by Particle-aMALA/Particle-MALA, with
Particle-RWM being the least efficient.

Figure 13: Autocorrelation of the energy from Figure 9 in the multivariate
stochastic volatility model.