# Nesting Particle Filters for Experimental Design in Dynamical Systems

Sahel Iqbal  Adrien Corenflos  Simo Särkkä  Hany Abdulsamad

###### Abstract

In this paper, we propose a novel approach to Bayesian Experimental Design
(BED) for non-exchangeable data that formulates it as risk-sensitive policy
optimization. We develop the Inside-Out SMC2 algorithm that uses a nested
sequential Monte Carlo (SMC) estimator of the expected information gain and
embeds it into a particle Markov chain Monte Carlo (pMCMC) framework to
perform gradient-based policy optimization. This is in contrast to recent
approaches that rely on biased estimators of the expected information gain
(EIG) to amortize the cost of experiments by learning a design policy in
advance. Numerical validation on a set of dynamical systems showcases the
efficacy of our method in comparison to other state-of-the-art strategies.

Bayesian experimental design, sequential Monte Carlo, non-exchangeable data,
amortized inference

##  1 Introduction

Traditionally, Bayesian inference on the parameters of a statistical model is
performed _after the fact_ by employing a posterior elicitation routine to
previously gathered data. However, in many scenarios, experimenters can
proactively _design_ experiments to acquire maximal information about the
parameter of interest. Bayesian experimental design (BED) (Lindley, 1956;
Chaloner & Verdinelli, 1995) offers a principled framework to achieve this
goal by maximizing the expected information gain (EIG) over the design space.
BED has found applications in fields as varied as physics (McMichael et al.,
2021; Caouette-Mansour et al., 2022), neuroscience (Shababo et al., 2013),
active learning (Bickford Smith et al., 2023), robotics (Schultheis et al.,
2020), and psychology (Myung et al., 2013). A recent overview of the field of
BED can be found in Rainforth et al. (2023).

In Bayesian experimental design, we are given a prior $p(\theta)$ and a
likelihood $p(x\mid\xi,\theta)$, where $\theta\in\Theta$ is the set of unknown
parameters of interest, $x\in\mathcal{X}$ is the experimental outcome, and
$\xi\in\Xi$ is a controllable design. The information gain (IG, Lindley,
1956), in the parameter $\theta$, as a function of the design $\xi$ is defined
as

$\mathcal{G}(\xi)\coloneq\mathbb{H}[p(\theta)]-\mathbb{H}[p(\theta\mid
x,\xi)],$ (1)

where $\mathbb{H}[p(.)]$ is the entropy of a random variable defined as

$\mathbb{H}[p(\cdot)]\coloneq-\mathbb{E}_{p(\cdot)}[\log p(\cdot)].$ (2)

However, given that the outcomes $x$ are themselves random variables for a
fixed design $\xi$, the quantity defined in (1) can be only evaluated in
expectation, leading to the _expected_ information gain (Lindley, 1956)

$\mathcal{I}(\xi)\coloneq\mathbb{E}_{p(x\mid\xi)}\bigl{[}\mathbb{H}[p(\theta)]-\mathbb{H}[p(\theta\mid
x,\xi)]\bigr{]},$ (3)

where $p(x\mid\xi)\coloneq\mathbb{E}_{p(\theta)}[p(x\mid\theta,\xi)]$ is the
marginal distribution of outcomes $x$ under the prior $p(\theta)$.
Consequently, the expected information gain measures the anticipated decrease
in uncertainty regarding the unknown variable $\theta$ following the
implementation of a design $\xi$. The optimal design $\xi^{*}$ is the
maximizer of the information measure $\mathcal{I}(\xi)$.

Despite its elegant nature, the BED framework involves maximising nested
expectations over intractable quantities such as the marginal likelihood of
$x$ and the posterior of $\theta$ appearing in (3). This is a computationally
intensive task (Kueck et al., 2009; Rainforth et al., 2018), even more salient
when optimizing designs for a series of experiments conducted sequentially,
where the impact of a single design needs to be accounted for across the
entire sequence of experiments. This makes the deployment of sequential BED
close to impossible on real-time systems with high-frequency data.

Huan & Marzouk (2016) addressed this limitation by introducing a parametric
policy to predict the designs as a function of the running parameter
posterior, eliminating the cost of the maximization step in each experiment.
Foster et al. (2021) extended this idea to amortize the overall cost of
sequential BED by conditioning the policy on the entire outcome-design
history, thus eliminating the cost of explicit posterior computation. In this
approach, called deep adaptive design (DAD), there is an upfront cost to learn
the policy, but experiments can be performed in real-time. While DAD is only
applicable for exchangeable models, implicit DAD (iDAD) (Ivanova et al., 2021)
generalizes the concept to accommodate non-exchangeable models that cover
dynamical systems. These methods leverage a lower bound to the expected
information gain known as the sequential prior contrastive estimator (sPCE).
This bound exhibits significant bias in low-sample regimes, leading to a
higher computational cost for accurate estimates of the EIG.

In this work, we introduce a novel amortization scheme that circumvents the
drawbacks of sPCE. Our approach is rooted in understanding sequential Bayesian
experimental design as an adaptive risk-sensitive decision-making process
(Whittle, 1990). We demonstrate that risk-sensitive decision-making can be
cast as an inference problem within an equivalent non-Markovian non-linear and
non-Gaussian state-space model (Toussaint & Storkey, 2006; Rawlik, 2013). This
insight motivates a novel nested sequential Monte Carlo (SMC) algorithm that
jointly estimates the EIG and the corresponding optimal designs. We refer to
this algorithm as Inside-Out SMC2 due to its relation to the SMC2 algorithm of
Chopin et al. (2013). Finally, we embed our technique within a particle Markov
chain Monte Carlo (pMCMC) algorithm to perform gradient-based optimization of
the amortizing policy. The code to reproduce our results is available at
_redacted, see supplementary material_.

##  2 Problem Statement

We are interested in the sequential BED problem for non-exchangeable models
and specifically dynamical systems. Accordingly, we assume a scenario of
$\,T\in\mathbb{N}$ sequential experiments to infer a parameter vector
$\theta$, starting from a prior $p(\theta)$, a Markovian likelihood
$f(x_{t+1}\mid x_{t},\xi_{t},\theta)$, and an initial distribution $p(x_{0})$,
where the subscript $t$ indexes the experiment number. We further assume that
the designs are sampled from a stochastic policy $\pi_{\phi}(\xi_{t}\mid
z_{0:t})$ parameterized by $\phi$, where we define $z_{0}\coloneq\\{x_{0}\\}$
and $z_{t}\coloneq\\{x_{t},\xi_{t-1}\\}$ for all $t\geq 1$, and denote the
outcome-design history up to time $t$ by
$z_{0:t}\coloneq\\{x_{0:t},\xi_{0:t-1}\\}$.

With the following factorization for the joint density of observations and
designs,

$\displaystyle p_{\phi}(z_{0:T}\mid\theta)$
$\displaystyle=p(z_{0})\prod_{t=1}^{T}p_{\phi}(z_{t}\mid z_{0:t-1},\theta)$
(4) $\displaystyle=\begin{aligned} p(x_{0})&\left\\{\prod_{t=1}^{T}f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)\right\\}\\\
&\times\left\\{\prod_{t=1}^{T}\pi_{\phi}(\xi_{t-1}\mid
x_{1:t-1},\xi_{1:t-2})\right\\},\end{aligned}$

the expected information gain can be written analogously to that of the single
experiment:

$\displaystyle\mathcal{I}(\phi)$
$\displaystyle\coloneq\mathbb{E}_{p_{\phi}(z_{0:T})}\Bigl{[}\mathbb{H}\bigl{[}p(\theta)\bigr{]}-\mathbb{H}\bigl{[}p(\theta\mid
z_{0:T})\bigr{]}\Bigr{]},$ (5)

where
$p_{\phi}(z_{0:T})=\mathbb{E}_{p(\theta)}\bigl{[}p_{\phi}(z_{0:T}\mid\theta)\bigr{]}$.
Note that the EIG in (5) is evaluated under the expectation of the marginal
distribution $p_{\phi}(z_{0:T})$, including the stochastic design policy. This
definition corresponds to the terminal reward framework in literature (Foster,
2021, Section 1.8), as it compares the prior entropy of $\theta$ with the
posterior entropy at the end of the experiment sequence. The resulting
experimental design objective corresponds to finding the optimal policy
parameters
$\phi^{*}\coloneq\operatorname*{arg\,max}_{\phi}\,\,\mathcal{I}(\phi)$. The
upcoming sections will present a novel interpretation of this objective in a
sequential decision-making framework that leverages the duality with inference
techniques to perform policy amortization.

##  3 Sequential Bayesian Experimental Design as Probabilistic Inference

To formulate sequential BED as an inference problem, we demonstrate a
factorization of the EIG over time steps.

######  Proposition 1.

For models specified by the joint density in (4), the expected information
gain factorizes to

$\mathcal{I}(\phi)=\mathbb{E}_{p_{\phi}(z_{0:T})}\biggl{[}\sum_{t=1}^{T}r_{t}(z_{0:t})\biggr{]},$
(6)

where $r_{t}(z_{0:t})$ is a stage reward defined as

$r_{t}(z_{0:t})=\alpha_{t}(z_{0:t})+\beta_{t}(z_{0:t}),$ (7)

with $\alpha_{t}(z_{0:t})$ and $\beta_{t}(z_{0:t})$ defined as

$\displaystyle\\!\\!\\!\alpha_{t}(z_{0:t})$ $\displaystyle=\int p(\theta\mid
z_{0:t})\log f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\,\mathrm{d}\theta,$
$\displaystyle\\!\\!\\!\beta_{t}(z_{0:t})$
$\displaystyle\\!=-\log\\!\int\\!p(\theta\mid z_{0:t-1})f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)\,\mathrm{d}\theta.$

Furthermore, for models with additive, constant noise in the dynamics, the EIG
can be written as

$\mathcal{I}(\phi)\equiv\mathbb{E}_{p_{\phi}(z_{0:T})}\left[\sum_{t=1}^{T}\beta_{t}(z_{0:t})\right],$
(8)

where ’$\equiv$’ denotes equality up to an additive constant.

###### Proof.

See Appendix A. ∎

Written in this form, the EIG resembles the expected total reward of a
discrete-time, finite-horizon non-Markovian decision-making problem (Puterman,
2014) with a stage reward $r_{t}(z_{0:t})$ that captures the information
content regarding the unknown parameters $\theta$. Moving forward, we use this
factorization of $\mathcal{I}(\phi)$ to derive a risk-sensitive objective and
a dual inference perspective, that leads to a novel amortized BED learning
scheme.

###  3.1 The Dual Inference Problem

To leverage the duality between risk-sensitive decision-making and inference,
we follow the formulation of Toussaint & Storkey (2006) and Rawlik (2013), and
introduce the potential function

$g_{t}(z_{0:t})\coloneq\exp\Big{\\{}\eta\,r_{t}(z_{0:t})\Big{\\}},$ (9)

with $\eta\in\mathbb{R}_{>0}$. If we define the potential of an entire
trajectory, $g_{1:T}(z_{0:T})$, to be the product of the potential functions
over time steps, then

$\\!\\!\\!\\!\log g_{1:T}(z_{0:T})=\sum_{t=1}^{T}\log
g_{t}(z_{0:t})\\!=\eta\left[\sum_{t=1}^{T}r_{t}(z_{0:t})\right],$ (10)

is the total reward of a trajectory scaled by $\eta$. In this context, the
potentials $g_{1:T}$ play the role of an un-normalized pseudo-likelihood
proportional to the _probability_ of the trajectory $z_{0:T}$ being optimal
(Dayan & Hinton, 1997). This perspective allows us to define a non-Markovian
state-space model characterized by the following joint density

$\Gamma_{\phi}(z_{0:T})=\frac{1}{Z(\phi)}\,p(z_{0})\prod_{t=1}^{T}p_{\phi}(z_{t}\mid
z_{0:t-1})\,g_{t}(z_{0:t}),$ (11)

where $p_{\phi}(z_{t}\mid z_{0:t-1})$ are the marginal dynamics under the
running filtered posterior $p(\theta\mid z_{0:t-1})$

$\\!\\!\\!p_{\phi}(z_{t}\mid z_{0:t-1})\\!=\\!\int\\!p_{\phi}(z_{t}\mid
z_{0:t-1},\theta)\,p(\theta\mid z_{0:t-1})\,\mathrm{d}\theta,$ (12)

and $Z(\phi)$ is the normalizing constant

$Z(\phi)=\int g_{1:T}(z_{0:T})\,p_{\phi}(z_{0:T})\,\mathrm{d}z_{0:T}.$ (13)

For ease of exposition, we will henceforth refer to $Z(\phi)$ as the marginal
likelihood of _being optimal_ , even though it may not represent any
meaningful probability.

The duality principle becomes evident when we apply Jensen’s inequality to
show that the log marginal likelihood is an upper bound on the EIG scaled by
$\eta$:

$\displaystyle\log Z(\phi)$
$\displaystyle=\log\mathbb{E}_{p_{\phi}(z_{0:T})}\left[g_{1:T}(z_{0:T})\right]$
(14a) $\displaystyle\geq\mathbb{E}_{p_{\phi}(z_{0:T})}\left[\log
g_{1:T}(z_{0:T})\right]$ (14b) $\displaystyle=\eta\,\mathcal{I}(\phi).$ (14c)

Hence maximizing the marginal likelihood is equivalent to maximizing a risk-
sensitive EIG objective, which we denote as $\mathcal{I}_{\eta}(\phi)$ (Marcus
et al., 1997; Rawlik, 2013)

$\mathcal{I}_{\eta}(\phi)=\frac{1}{\eta}\log\mathbb{E}_{p_{\phi}(z_{0:T})}\left[\exp\left\\{\eta\sum_{t=1}^{T}r_{t}(z_{0:t})\right\\}\right].$
(15)

Note that $\eta$ modulates the bias-variance trade-off of this objective. This
aspect is revealed by considering a first-order expansion of the objective
around $\eta=0$

$\mathcal{I}_{\eta}(\phi)\approx\mathbb{E}\left[\sum_{t=1}^{T}r_{t}(z_{0:t})\right]-\frac{\eta}{2}\mathbb{V}\left[\sum_{t=1}^{T}r_{t}(z_{0:t})\right],$
(16)

where the expectation and variance operators are in relation to
$p_{\phi}(z_{0:T})$. Note that in the limit $\eta\rightarrow 0$, we recover
the risk-neutral EIG objective (Whittle, 1990).

In this section, we formalized the connection between a risk-sensitive
sequential BED objective and inference in an equivalent non-Markovian state-
space model. In the following, we leverage this insight to formulate a
gradient-based policy optimization technique within a particle MCMC framework.

###  3.2 Amortization as Likelihood Maximization

The duality principle demonstrated in Section 3 enables us to view amortized
BED from an inference-centric perspective and to frame policy optimization in
terms of maximum likelihood estimation (MLE) in a non-Markovian, non-linear,
and non-Gaussian state-space model, as specified by (11). Following the
literature on particle methods for MLE (Kantas et al., 2015), we employ a
stochastic gradient ascent algorithm (Robbins & Monro, 1951).

To obtain the derivative of the log marginal likelihood,
$\mathcal{S}(\phi)\coloneq\nabla_{\phi}\log Z(\phi)$, also known as the score
function, we make use of Fisher’s identity (Cappé et al., 2005),

$\displaystyle\mathcal{S}(\phi)$
$\displaystyle=\int\nabla_{\phi}\log\tilde{\Gamma}_{\phi}(z_{0:T})\,\Gamma_{\phi}(z_{0:T})\,\mathrm{d}z_{0:T}$
$\displaystyle=\int\nabla_{\phi}\log
p_{\phi}(z_{0:T})\,\Gamma_{\phi}(z_{0:T})\,\mathrm{d}z_{0:T}.$

where we define $\tilde{\Gamma}_{\phi}(z_{0:T})\coloneq
p_{\phi}(z_{0:T})\,g_{1:T}(z_{0:T})$ to be the un-normalized density from
(11). This identity provides a Monte Carlo estimate of the score
$\hat{\mathcal{S}}(\phi)$ given samples from $\Gamma_{\phi}$. It is well-known
that computing this expectation naively by first sampling from
$p_{\phi}(z_{0:T})$ and then weighting the samples by $g$ results in very high
variance estimates (see, e.g., the discussion in Doucet et al., 2009, Section
3.3). However, approximate samples from $\Gamma_{\phi}$ can still be drawn via
particle smoothing methods, and they yield consistent, albeit biased,
estimates of expectations under the smoothing distribution for a finite sample
size (Chopin & Papaspiliopoulos, 2020, Chapter 12).

Algorithm 1 Markovian score climbing for maximum likelihood estimation

0: Initial trajectory $z_{0:T}^{0}$, initial parameters $\phi_{0}$, step size
sequence $\\{\gamma_{i}\\}_{i\in\mathbb{N}}$, Markov kernel $\mathcal{K}$.

0: Local optimum $\phi^{*}$ of the marginal likelihood

1: $k\leftarrow 1$

2: while not converged do

3: Sample $z_{0:T}^{k}\sim\mathcal{K}(\cdot\mid z_{0:T}^{k-1},\phi_{k-1})$

4: Compute $\hat{\mathcal{S}}(\phi_{k-1})\leftarrow\nabla_{\phi}\log
p_{\phi}(z_{0:T})|_{\phi=\phi_{k-1}}$

5: Update
$\phi_{k}\leftarrow\phi_{k-1}+\gamma_{k}\,\hat{\mathcal{S}}(\phi_{k-1})$

6: $k\leftarrow k+1$

7: end while

8: return $\phi_{k}$

Instead, we use Markovian score climbing (MSC) (Gu & Kong, 1998; Naesseth et
al., 2020), presented in Algorithm 1. MSC uses a $\Gamma_{\phi}$-ergodic
Markov chain Monte Carlo (MCMC, see, e.g., Brooks et al., 2011, for a review
and definition) kernel, $\mathcal{K}(\cdot\mid z_{0:T},\phi)$, to compute a
Monte Carlo estimate of the score function. Contrary to simply using particle
smoother approximations within a gradient ascent procedure (see Kantas et al.,
2015, Section 5), Algorithm 1 is guaranteed to converge to a local optimum of
the marginal likelihood (Naesseth et al., 2020, Proposition 1).

In this work, we construct the MCMC kernel $\mathcal{K}$ as a variant of the
conditional sequential Monte Carlo (CSMC) kernel (Andrieu et al., 2010),
namely the Rao–Blackwellized CSMC kernel (Olsson & Westerborn, 2017; Cardoso
et al., 2023; Abdulsamad et al., 2023). In practice, CSMC (as well as its
Rao–Blackwellized modifications aforementioned) can be implemented as a simple
modification to a particle filter representation of the smoothing distribution
(Kitagawa, 1996, Section 4.1). In the next section, we describe how such a
particle filter can be implemented for (11), and, for the sake of brevity, we
defer the full description of the conditional version to Section C.4.

##  4 Inside-Out SMC2

In this section, we describe how to sample from a consistent approximation of
(11) using nested particle filters.

###  4.1 Approximating the Filtered Posterior

A _bootstrap_ particle filter samples particles from the transition density
and weights them using the potential function (Chopin & Papaspiliopoulos,
2020, Section 10.3). In our non-Markovian model, this would imply sampling
from the marginal dynamics $p_{\phi}(z_{t}\mid z_{0:t-1})$, (12), and
evaluating the potential function $g_{t}(z_{0:t})$, (9). Both of these require
computing the filtered posterior $p(\theta\mid z_{0:t-1})$. Fortunately, for
models of the form

$\theta\sim p(\theta),\quad z_{t}\sim p_{\phi}(z_{t}\mid
z_{0:t-1},\theta),\quad t\geq 1,$

the iterated batch importance sampling (IBIS) algorithm of Chopin (2002) can
be used to generate Monte Carlo samples
$\\{\theta^{m}_{t}\\}_{m=1}^{M}\eqqcolon\theta^{1:M}_{t}$ that are
approximately distributed according to $p(\theta\mid z_{0:t-1})$ at each time
step using a specialized particle filtering procedure. We summarize the method
in Algorithm 2.

Algorithm 2 Single step of IBIS

0: Any operation with superscript $m$ is to be understood as performed for all
$m=1,\dots,M$.

0: Ibis_Step$(z_{0:t},\theta^{1:M},W^{1:M})$

1: Compute $v_{t}(\theta^{m})=p_{\phi}(z_{t}\mid z_{0:t-1},\theta^{m})$.

2: Reweight: $W^{m}\propto W^{m}v_{t}(\theta^{m})$.

3: if some degeneracy criterion is fulfilled then

4: Resample: $a_{t}^{m}\sim\mathcal{M}(W^{1:M})$.

5: Move: $\tilde{\theta}^{m}\sim K_{t}(\theta^{a_{t}^{m}},\cdot)$.

6: Replace the current set of weighted particles with
$(\theta^{m},W^{m})\leftarrow(\tilde{\theta}^{m},1/M)$.

7: end if

8: return $\\{\theta^{m},W^{m}\\}_{m=1}^{M}$

8:

In Algorithm 2, $\mathcal{M}(W^{1:M})$ denotes multinomial sampling using the
normalized weights $W^{1:M}$, and $K_{t}$ is a $p(\theta\mid z_{0:t})$-ergodic
Markov chain Monte Carlo kernel. If the degeneracy criterion is met, IBIS
employs a resample-move step (lines 4–6 in Algorithm 2, see Gilks & Berzuini,
2001) to rejuvenate particles using the Markov kernel $K_{t}$. A standard
degeneracy measure, which we use in this work, is given by the effective
sample size (ESS) (Liu & Chen, 1995) of the particle representation computed
as $\textrm{ESS}=1/\sum_{m=1}^{M}(W^{m})^{2}$, corresponding roughly to the
number of equivalent independent samples one would need to compute integrals
with the same precision. The resample-move step is then triggered if the ESS
falls below a chosen fraction of the total particles $M$, taken to be $75\%$
in this work. Details on the choice of the Markov kernel $K_{t}$ are given in
Appendix C.2.

###  4.2 The Inside-Out SMC2 Algorithm

Chopin (2004) showed that, for integrable functions $\psi$,
$\sum_{m=1}^{M}W^{m}\,\psi(\theta^{m})$ is a consistent and asymptotically (as
$M\to\infty$) normal estimator of the integral $\int\psi(\theta)\,p(\theta\mid
z_{0:t})\,\mathrm{d}z_{0:t}$. Thus, a natural solution to perform inference in
(11) is to use IBIS within a standard particle filter. This idea is similar to
the SMC2 algorithm of Chopin et al. (2013) which can be seen as a particle
filter within IBIS targeting the distribution $p(\theta\mid y_{1:t})$ for a
state-space model with noisy observations $y_{1:t}$. We thus call our
algorithm Inside-Out SMC2, which we reproduce in Algorithm 3, with $N$ and $M$
denoting the numbers of samples of $z$ and $\theta$ respectively. Any
operation therein with superscripts $m$ or $n$ is to be understood as
performed for every $m=1,\dots,M$ and $n=1,\dots,N$, and we use an upper
script dot ${u}^{\bullet n}$ to denote the collection
$\\{{u}^{mn}\\}_{m=1}^{M}$.

Algorithm 3 Inside-Out SMC2

1: Sample $z_{0}^{n}\sim p(z_{0})$, $\theta_{0}^{mn}\sim p(\theta)$ and set
$W_{0,\theta}^{mn}\leftarrow 1/M$.

2: Sample $z_{1}^{n}\sim\hat{p}_{\phi}(\cdot\mid z_{0}^{n})$ and initialize
the state history $z_{0:1}^{n}\leftarrow(z_{0}^{n},z_{1}^{n})$.

3: Compute and normalize the weights
$\smash{W_{z}^{n}\propto\exp\Big{\\{}-\eta\log\hat{p}(x_{1}^{n}\mid
x_{0}^{n},\xi_{0}^{n})\Big{\\}}.}$

4: for $t\leftarrow 1,\dots,T-1$ do

5: Sample $b_{t}^{n}\sim\mathcal{M}(W_{z}^{1:N})$.

6: $\theta_{t}^{\bullet n},W^{\bullet
n}_{t,\theta},\leftarrow\textsc{Ibis\\_Step}(z_{0:t}^{b_{t}^{n}},\theta_{t-1}^{\bullet
b_{t}^{n}},W^{\bullet b_{t}^{n}}_{t-1,\theta})$

7: Sample $z_{t+1}^{n}\sim\hat{p}_{\phi}(\cdot\mid z_{0:t}^{b_{t}^{n}}),$ and
append to state history
$z_{0:t+1}^{n}\leftarrow[z_{0:t}^{b_{t}^{n}},z_{t+1}^{n}]$.

8: Compute and normalize the weights
$\smash{W_{z}^{n}\propto\exp\left\\{-\eta\log\hat{p}(x_{t+1}^{n}\mid
z_{0:t}^{b_{t}^{n}},\xi_{t}^{n})\right\\}.}$

9: end for

10: return $\\{z_{0:T}^{n},W_{z}^{n}\\}_{n=1}^{N}$.

The heuristic behind the algorithm is to use IBIS to form a mixture
distribution that approximates the marginal dynamics at each time step. More
concretely, assume we are at time step $t$, and we have a trajectory
$z_{0:t}^{n}$. Using IBIS we obtain a set of weighted particles
$\\{\theta_{t}^{\bullet n},W^{\bullet n}_{t,\theta}\\}$ that approximates the
$\theta$ posterior

$p(\theta\mid z_{0:t}^{n})\approx\hat{p}(\theta\mid
z_{0:t}^{n})\coloneq\sum_{m=1}^{M}W_{t,\theta}^{mn}\,\delta_{\theta_{t}^{mn}}(\theta),$

where $\delta$ is the Dirac delta function. We can then form an approximation
to the marginal dynamics for the state $x_{t+1}$ as

$\displaystyle\hat{p}(x_{t+1}\mid z_{0:t}^{n},\xi_{t})$ $\displaystyle=\int
f(x_{t+1}\mid x_{t}^{n},\xi_{t},\theta)\,\hat{p}(\theta\mid
z_{0:t}^{n})\,\mathrm{d}\theta$
$\displaystyle=\sum_{m=1}^{M}W_{t,\theta}^{mn}\,f(x_{t+1}\mid
x_{t}^{n},\xi_{t},\theta_{t}^{mn}),$

and consequently, for the augmented state as

$\hat{p}_{\phi}(z_{t+1}\mid z_{0:t}^{n})=\hat{p}(x_{t+1}\mid
z_{0:t}^{n},\xi_{t})\,\pi_{\phi}(\xi_{t}\mid z_{0:t}^{n}).$

If the Markovian density $f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)$ is
conditionally linear in the parameters and conjugate to the prior $p(\theta)$,
we can compute the posterior in closed form, and therefore the marginal
dynamics as well. In this case, we do not need IBIS, and this significantly
reduces the computational complexity of Algorithm 3 (see Section C.3).

Note that, for the sake of clarity, in Algorithm 3 and in its analysis in
Section 4.3, we only consider the case when the dynamics have constant noise
(corresponding to the stage reward in (8)). The modification to Algorithm 3
for the more general case is straightforward by simply replacing the weight
function, which is detailed in Section C.1.

###  4.3 Target Distribution of Inside-Out SMC2

We now show that the nested particle filter introduced in the previous section
asymptotically targets the correct distribution. Similarly to Chopin et al.
(2013, Proposition 1), Algorithm 3 is a particle filter targeting a particle
filter. Indeed, dropping the $n$ indices and the explicit dependence on
$\phi$, let $\Gamma_{t}^{M}(z_{0:t},\theta_{0:t-1}^{1:M},a_{1:t-1}^{1:M})$
denote the distribution of all stochastic variables generated by an instance
of the inner IBIS at line 6 in Algorithm 3. It is evident that

$\Gamma_{0}^{M}(z_{0})=p(z_{0}),\quad\Gamma^{M}_{0}(z_{0},\theta_{0}^{1:M})=p(z_{0})\prod_{m=1}^{M}p(\theta_{0}^{m}).$

Let us break the ratio of the distributions over successive iterations as

$\frac{\Gamma^{M}_{t+1}(z_{0:t+1},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t-1}^{1:M},a_{1:t-1}^{1:M})}=\begin{aligned}
&\frac{\Gamma^{M}_{t+1}(z_{0:t+1},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}\\\
&\times\frac{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t-1}^{1:M},a_{1:t-1}^{1:M})}.\end{aligned}$

The second fraction is the IBIS rejuvenation step:

$\frac{\Gamma_{t}^{M}(z_{0:t},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t-1}^{1:M},a_{1:t-1}^{1:M})}=\prod_{m=1}^{M}W_{t,\theta}^{a_{t}^{m}}\,K_{t}(\theta_{t-1}^{a_{t}^{m}},\theta_{t}^{m}),$
(17)

where the normalized weights (of the $\theta$ particles) $W_{t,\theta}^{m}$
are given by

$W_{t,\theta}^{m}=\frac{v_{t}^{m}}{\sum_{m=1}^{M}v_{t}^{m}},\quad
v_{t}^{m}=p_{\phi}(z_{t}\mid z_{0:t-1},\theta_{t-1}^{m}),$

and we have assumed multinomial resampling at every time step. The first
fraction is the trajectory update step:

$\frac{\Gamma^{M}_{t+1}(z_{0:t+1},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}\propto\frac{1}{M}\sum_{m=1}^{M}p_{\phi}(z_{t+1}\mid
z_{0:t},\theta_{t}^{m})\\\
\times\exp\left\\{-\eta\log\frac{1}{M}\sum_{m=1}^{M}f(x_{t+1}\mid
x_{t},\xi_{t},\theta_{t}^{m})\right\\}.$ (18)

The following proposition, akin to the law of large numbers, proven in
Appendix B ensures that Algorithm 3 asymptotically targets the correct
distribution.

######  Proposition 2 (Consistency of the target distribution).

Let
$\tilde{\Gamma}^{M}_{t}(z_{0:t},\theta_{t})=\mathbb{E}\left[\Gamma_{t}^{M}(z_{0:t},\theta_{t}^{1:M},a_{1:t}^{1:M})\right]$
be the joint expected empirical distribution111Note that
$\Gamma_{t}^{M}(z_{0:t},\theta_{t}^{1:M},a_{1:t}^{1:M})$ is a random measure
and hence this expectation is a measure. over $(z_{0:t},\theta_{t})$ taken by
integrating over $a_{1:t}^{1:M}$ and $\theta_{0:t}^{1:M}$, then under
technical conditions listed in Appendix B, as $M\to\infty$, empirical
expectations under $\tilde{\Gamma}^{M}_{t}$ converge to expectations under
$\Gamma(z_{0:t},\theta_{t})$, that is, for any bounded test function
$\psi(z_{0:t},\theta_{t})$, we have

$\mathbb{E}_{\tilde{\Gamma}^{M}_{t}}\left[\psi(z_{0:t},\theta_{t})\right]\to\mathbb{E}_{\Gamma(z_{0:t},\theta_{t})}\left[\psi(z_{0:t},\theta_{t})\right]$
(19)

almost surely.

##  5 Related Work

Foster et al. (2021) optimizes the following sequential Prior Contrastive
Estimation (sPCE) lower bound

$\mathcal{L}_{T}^{\textrm{sPCE}}(\phi,L)=\mathbb{E}_{p(\theta_{0},z_{0:T}\mid\phi)p(\theta_{1:L})}\bigl{[}g_{L}(\theta_{0:L},z_{0:T}\bigr{]},$
(20)

where

$g_{L}(\theta_{0:L},z_{0:T})=\log\frac{p(z_{0:T}\mid\theta_{0},\phi)}{\frac{1}{L+1}\sum_{\ell=0}^{L}p(z_{0:T}\mid\theta_{\ell},\phi)}.$
(21)

A naive nested Monte Carlo estimator would exclude $\theta_{0}$ from the
denominator, but this estimator would have a large variance, particularly for
a large number of experiments $T$. By including $\theta_{0}$ in the estimate
for $p(z_{0:T}\mid\phi)$, the authors show that $g_{L}$ is upper bounded by
$\log(L+1)$, resulting in a low variance estimator amenable to optimization
(Poole et al., 2019). However, this implies that, to not be artificially
restricted by this upper bound of the EIG estimator, one needs to take an
exponentially large value of $L$. The InfoNCE bound (van den Oord et al.,
2019) used by iDAD (Ivanova et al., 2021) exhibits the same $\log(L+1)$ bound.
The proposed solution is to introduce an amortized parametric proposal
density, $q(\theta;z_{0:T})$, from which the regularizing samples
$\theta_{1:L}$ are drawn instead of the prior (Foster et al., 2021, Appendix
B). By explicitly leveraging the Markovian structure of our model, our Inside-
Out SMC2 algorithm does not suffer from this particular drawback and only
requires a small number of $\theta$ particles to estimate the EIG with
reasonable precision.

Drovandi et al. (2013), Drovandi et al. (2014), and Moffat et al. (2020) have
previously proposed the use of IBIS for sequential experimental design.
However, they are only concerned with exchangeable models, and IBIS is used
solely as a way to track the posterior over the parameters given the past
experiments. We do not use IBIS on its own, but rather as one element of a
broader, general particle smoothing algorithm for dynamical models with
marginalized dynamics. Furthermore, their design optimization is myopic, only
optimizing for the next experiment, in contrast to our approach that optimizes
across the entire horizon of experiments.

Our work also shares similarities in formulation with Blau et al. (2022) which
formulates the sequential BED problem as a hidden parameter Markov decision
process (Doshi-Velez & Konidaris, 2016). Unlike our work, they adopt a
reinforcement learning approach that applies only to exchangeable models and
they optimize for the sPCE bound.

In the context of SMC, a related approach to ours is given by Wigren et al.
(2019) who too propagate the filtering posterior for the parameter at hand to
perform parameter-marginalized inference. However because they are not
interested in experimental design, and in contrast to our Inside-Out SMC2,
they only need to compute the _marginal_ posterior distribution, whereas we
reuse the _pathwise_ posterior distribution to compute our potential function
$g$. Nonetheless, we believe the two approaches are related and could be
combined in future works.

##  6 Empirical Evaluation

For the empirical evaluation of our method, we consider the setting of
designing 50 sequential experiments to identify the parameters of input
(design) dependent dynamical systems described by stochastic differential
equations discretized using the Euler–Maruyama method (see, e.g., Särkkä &
Svensson, 2023, Section 4.3). Hence all scenarios involve likelihoods that are
non-exchangeable Gaussian densities.

We construct the stochastic policy $\pi_{\phi}$ using a mean function
$m_{\phi}$ parameterized by a long short-term memory architecture (Hochreiter
& Schmidhuber, 1997), with an additional learnable parameter $\Sigma_{\phi}$
for the variance. Then the stochastic policy $\pi_{\phi}$ is constructed as
the law of the random variable $\xi_{t}=a\cdot\tanh(s_{t})+b$, where
$s_{t}\sim\mathcal{N}(m_{\phi}(z_{0:t}),\Sigma_{\phi})$. Here $(a,b)$ are
scale and shift parameters that reflect the design constraints. Exhaustive
details of the network architecture and hyperparameters are given in Section
D.

During evaluation, we compare to a baseline random policy and iDAD (Ivanova et
al., 2021), as it can accommodate non-exchangeable sequential experiments.
Since we have access to the conditional transition densities in our
experiments, this eliminates the need to use a dedicated critic network and we
train iDAD on the sPCE lower bound, facilitating a fair comparison.

As evaluation metrics, we report both the sPCE bound and a Monte Carlo
estimate based on Algorithm 3. The latter can be retrieved by turning off the
resampling step for the $z$ particles (equivalent to setting $b_{t}^{n}=n$ in
line 5), which leads to a routine that generates trajectories from the
marginal distribution $p_{\phi}(z_{0:T})$. Assuming a particle approximation
$\sum_{n=1}^{N}W_{z}^{n}\,\delta_{z_{0:T}^{n}}(\cdot)\approx p_{\phi}(\cdot)$,
we can form a Monte Carlo estimate of the EIG as follows

$\mathcal{I}(\phi)\approx\sum_{n=1}^{N}W_{z}^{n}\biggl{[}\sum_{t=1}^{T}\hat{r}_{t}(z_{0:t}^{n})\biggr{]}.$
(22)

Here $\hat{r}_{t}$ denotes the log potential functions for $\eta=1$. By
accumulating this quantity over experiments, we can use Algorithm 3 as an
estimator of the EIG.

$10$$20$$30$$40$$50$$0$$2$$4$$6$$8$$10$$12$Number of ExperimentsInformation
GainIO-SMC2IO-SMC2 (Exact)Random Figure 1: Information gain of different
policies computed in closed-form on the conditionally-linear stochastic
pendulum. We report the mean and standard deviation over 512 realizations.
$0$$2$$4$$6$$8$$10$$12$$14$$16$$0$$2$$4$$6$$8$$10$$12$Training EpochsEIG
EstimateIO-SMC2IO-SMC2 (Exact) Figure 2: Training progression of the IO-SMC2
policy and its exact variant on the conditionally-linear stochastic pendulum.
At every epoch, we evaluate the EIG estimate using the mean policy. We report
the mean and standard deviation over 25 training runs.

###  6.1 Stochastic Pendulum

Our first experiment deals with two different representations of the
stochastic dynamics of a compound pendulum. In this experiment, the aim is to
infer the mass and length of the pendulum. The outcome of an experiment is an
observation of the state of the pendulum characterized by its angular position
and velocity, while the design is the torque applied as input to the system.

####  6.1.1 Conditionally Linear formulation

First, we consider a conditionally linear formulation of the dynamics of the
pendulum, see Appendix D.2.1. Here, we assume a Gaussian prior over the
parameters. This setting allows us to compare IO-SMC2 with exact posterior
computation against the approximate posteriors delivered by IBIS. Details on
exact posterior inference are found in Section C.3.

In Figure 1, we see that policies trained using both algorithms outperform the
random baseline. In addition, we also observe that training policies using
exact marginal computations result in slightly better policies, and this is
confirmed in Figure 2. This result is in line with expectations since the
inner particle filter adds variance to the score function estimate. Thus,
where applicable, we recommend using the exact version of our algorithm.

$0$$10$$20$$30$$40$$50$$-5$$0$$5$Number of Experiments$q$$\dot{q}$$\xi$ Figure
3: A sample experiment trajectory generated by the amortizing policy during
deployment on the non-linear stochastic pendulum environment. $q$ is the angle
of the pendulum from the vertical, $\dot{q}$ is the angular velocity and $\xi$
is the design.

####  6.1.2 Non-linear formulation

We now consider the standard, non-linear version of the stochastic pendulum.
We place a log-normal prior that ensures these parameters can take only
positive values. Hence the exact version of our algorithm is no longer
applicable and we have to use IO-SMC2. EIG estimates and lower bounds for
different policies in this environment are given in Table 1. An example of a
trajectory generated by the mean of the trained policy is given in Figure 3.
The policy has learnt to swing the pendulum to achieve greater and greater
angular velocities, and consequently explore the sample space by alternating
between maximum positive and negative designs.

Table 1: EIG estimates and sPCE lower bounds on the non-linear pendulum
experiment for the considered methods. The EIG estimate was computed with
$N=256$ and $M=64$ samples. The lower bound was evaluated using $L=5\times
10^{4}$ regularizing samples and $16$ outer trajectories. We report the mean
$\pm$ standard deviation for 25 evaluations.
Policy | EIG Estimate | sPCE
---|---|---
Random | $3.75\pm 0.30$ | $3.66\pm 0.32$
iDAD | - | $2.85\pm 0.31$
Ours | $\bf{4.75\pm 0.20}$ | $\bf{4.57\pm 0.34}$

###  6.2 Stochastic Cart-Pole

For our final experiment, we consider a cart-pole system with additive noise
in the acceleration of the cart. The unknown parameters are the mass and
length of the pole. The outcome of every experiment is an observation of the
cart-pole state consisting of the position and velocity of the cart and the
angular position and angular velocity of the pole. The design is the force
applied to the cart at discrete intervals. We place a log-normal prior for the
parameters. Further details on the experiment setup are given in Appendix D.3.

For this experiment, we only compare our algorithm to the random policy
baseline. Due to the very high final value of the EIG, iDAD requires at least
$\exp{(12)}\approx 160000$ regularizing samples for a sufficiently low-bias
sPCE estimate, which was computationally unfeasible. Figure 4 shows how the
EIG estimate of the mean policy improves over training epochs. Finally, a
sample trajectory generated by the trained policy mean is depicted in Figure
5.

Table 2: EIG estimates the stochastic cart-pole experiment for the considered
methods. The EIG estimate was computed with $N=256$ and $M=128$ samples. We
report the mean $\pm$ standard deviation for 25 evaluations.
Policy | EIG Estimate
---|---
Random | $7.10\pm 0.36$
Ours | $\bf{11.65\pm 0.50}$
$2$$4$$6$$8$$10$$5$$10$Training EpochsEIG EstimateIO-SMC2 Figure 4: Training
evolution of the policy for the stochastic cart-pole experiment. We estimate
the EIG after each training epoch, and we report the mean and standard
deviation over 25 seeds. $0$$10$$20$$30$$40$$50$$-4$$-2$$0$$2$$4$Number of
Experiments$q$$\dot{q}$$\xi$ Figure 5: A sample experiment trajectory
generated by the amortizing policy during deployment on the stochastic cart-
pole environment. Here, $q$ is the angle of the pole from the vertical,
$\dot{q}$ is the angular velocity and $\xi$ is the action.

##  7 Discussion and Limitations

We have introduced a novel method of amortized sequential Bayesian
experimental design, taking inspiration from the control-as-inference
framework. We formulated the problem of sequential BED as that of maximum
likelihood estimation in a non-Markovian time series model. To perform
inference in this model, we constructed a novel particle smoothing algorithm,
which we call the Inside-Out SMC2 as it forms the converse structure of the
SMC2 of Chopin et al. (2013). While we have used this algorithm within the
context of BED, we believe it may find uses in other settings, where one
wishes to obtain pathwise smoothing of filtering trajectories under parameter-
marginalized models.

The primary limitation of our approach is the requirement to evaluate the
conditional transition densities in closed form. Thus our method cannot be
applied to address sequential experimental design problems arising in
dynamical models with intractable densities, for example, for choosing optimal
measurement times in a compartmental epidemic model that is modeled as a
Markov jump process (Whitehouse et al., 2023). We also do not address the
question of how to choose an optimal tempering parameter $\eta$.

Nonetheless, our experimental evaluation shows that, despite these
limitations, our approach shows promise as a generic and efficient way to
learn experimental design policies: it is easy to train, and the policies it
learns outperform the main alternative (iDAD, Ivanova et al., 2021) at a
fraction of the cost. Indeed, with the hyperparameter settings as outlined in
Appendix D.2.2, our algorithm converges in the order of minutes on a CPU,
while iDAD requires multiple hours on a GPU (with their implementation). We
are hopeful that this opens a door to a new approach to experimental design in
non-exchangeable data.

## Impact Statement

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none of
which we feel must be specifically highlighted here.

## References

  * Abdulsamad et al. (2023) Abdulsamad, H., Iqbal, S., Corenflos, A., and Särkkä, S.  Risk-sensitive stochastic optimal control as Rao-Blackwellized Markovian score climbing.  _arXiv preprint arXiv:2312.14000_ , 2023.
  * Andrieu et al. (2010) Andrieu, C., Doucet, A., and Holenstein, R.  Particle Markov chain Monte Carlo methods.  _Journal of the Royal Statistical Society: Series B_ , 72(3):269–342, 2010.
  * Belousov et al. (2019) Belousov, B., Abdulsamad, H., Schultheis, M., and Peters, J.  Belief space model predictive control for approximately optimal system identification.  In _4th Multidisciplinary Conference on Reinforcement Learning and Decision Making (RLDM)_ , 2019.
  * Bickford Smith et al. (2023) Bickford Smith, F., Kirsch, A., Farquhar, S., Gal, Y., Foster, A., and Rainforth, T.  Prediction-oriented Bayesian active learning.  In _Proceedings of the 26th International Conference on Artificial Intelligence and Statistics_ , 2023.
  * Blau et al. (2022) Blau, T., Bonilla, E. V., Chades, I., and Dezfouli, A.  Optimizing sequential experimental design with deep reinforcement learning.  In _International Conference on Machine Learning_ , 2022.
  * Brooks et al. (2011) Brooks, S., Gelman, A., Jones, G., and Meng, X.-L.  _Handbook of Markov chain Monte Carlo_.  CRC press, 2011.
  * Caouette-Mansour et al. (2022) Caouette-Mansour, M., Solyom, A., Ruffolo, B., McMichael, R. D., Sankey, J., and Childress, L.  Robust spin relaxometry with fast adaptive Bayesian estimation.  _Physical Review Applied_ , 17:064031, 2022.
  * Cappé et al. (2005) Cappé, O., Moulines, E., and Rydén, T.  _Inference in Hidden Markov Models_.  Springer Series in Statistics. Springer, 2005.
  * Cardoso et al. (2023) Cardoso, G., Janati El Idrissi, Y., Le Corff, S., Moulines, E., and Olsson, J.  State and parameter learning with PaRIS particle Gibbs.  In _International Conference on Machine Learning_ , volume 202, pp. 3625–3675. PMLR, 2023.
  * Chaloner & Verdinelli (1995) Chaloner, K. and Verdinelli, I.  Bayesian experimental design: A review.  _Statistical Science_ , 10(3):273–304, 1995.
  * Chopin (2002) Chopin, N.  A sequential particle filter method for static models.  _Biometrika_ , 89(3):539–552, 2002.
  * Chopin (2004) Chopin, N.  Central limit theorem for sequential Monte Carlo methods and its application to Bayesian inference.  _The Annals of Statistics_ , 32(6):2385–2411, 2004.
  * Chopin & Papaspiliopoulos (2020) Chopin, N. and Papaspiliopoulos, O.  _An introduction to sequential Monte Carlo_ , volume 4.  Springer, 2020.
  * Chopin et al. (2013) Chopin, N., Jacob, P. E., and Papaspiliopoulos, O.  SMC2: An efficient algorithm for sequential analysis of state space models.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 75(3):397–426, 2013.
  * Dayan & Hinton (1997) Dayan, P. and Hinton, G. E.  Using expectation-maximization for reinforcement learning.  _Neural Computation_ , 9(2):271–278, 1997.
  * Doshi-Velez & Konidaris (2016) Doshi-Velez, F. and Konidaris, G.  Hidden parameter Markov decision processes: A semiparametric regression approach for discovering latent task parametrizations.  In _International Joint Conference on Artificial Intelligence_ , 2016.
  * Doucet et al. (2009) Doucet, A., Johansen, A. M., et al.  A tutorial on particle filtering and smoothing: Fifteen years later.  _Handbook of nonlinear filtering_ , 12(656-704):3, 2009.
  * Drovandi et al. (2013) Drovandi, C. C., McGree, J. M., and Pettitt, A. N.  Sequential Monte Carlo for Bayesian sequentially designed experiments for discrete data.  _Computational Statistics & Data Analysis_, 57(1):320–335, 2013.
  * Drovandi et al. (2014) Drovandi, C. C., McGree, J. M., and Pettitt, A. N.  A sequential Monte Carlo algorithm to incorporate model uncertainty in Bayesian sequential design.  _Journal of Computational and Graphical Statistics_ , 23(1):3–24, 2014.
  * Foster et al. (2021) Foster, A., Ivanova, D. R., Malik, I., and Rainforth, T.  Deep adaptive design: Amortizing sequential Bayesian experimental design.  In _International Conference on Machine Learning_ , 2021.
  * Foster (2021) Foster, A. E.  _Variational, Monte Carlo and Policy-Based Approaches to Bayesian Experimental Design_.  PhD thesis, University of Oxford, 2021.
  * Gelman et al. (1997) Gelman, A., Gilks, W. R., and Roberts, G. O.  Weak convergence and optimal scaling of random walk Metropolis algorithms.  _The Annals of Applied Probability_ , 7(1):110 – 120, 1997.
  * Gilks & Berzuini (2001) Gilks, W. R. and Berzuini, C.  Following a moving target-Monte Carlo inference for dynamic Bayesian models.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 63(1):127–146, 2001.
  * Gu & Kong (1998) Gu, M. G. and Kong, F. H.  A stochastic approximation algorithm with Markov chain Monte-Carlo method for incomplete data estimation problems.  _Proceedings of the National Academy of Sciences_ , 95(13):7270–7274, 1998.
  * Hastings (1970) Hastings, W. K.  Monte Carlo sampling methods using Markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long Short-Term Memory.  _Neural Computation_ , 9(8):1735–1780, 1997.
  * Huan & Marzouk (2016) Huan, X. and Marzouk, Y. M.  Sequential Bayesian optimal experimental design via approximate dynamic programming.  _arXiv preprint arXiv:1604.08320_ , 2016.
  * Ivanova et al. (2021) Ivanova, D. R., Foster, A., Kleinegesse, S., Gutmann, M., and Rainforth, T.  Implicit deep adaptive design: Policy–based experimental design without likelihoods.  In _Advances in Neural Information Processing Systems_ , volume 34, pp. 25785–25798. Curran Associates, Inc., 2021.
  * Kantas et al. (2015) Kantas, N., Doucet, A., Singh, S. S., Maciejowski, J., and Chopin, N.  On particle methods for parameter estimation in state-space models.  _Statistical Science_ , 30(3), 2015.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kitagawa (1996) Kitagawa, G.  Monte Carlo filter and smoother for non-Gaussian nonlinear state space models.  _Journal of Computational and Graphical Statistics_ , 5(1):1–25, 1996.
  * Kueck et al. (2009) Kueck, H., Hoffman, M., Doucet, A., and de Freitas, N.  Inference and learning for active sensing, experimental design and control.  In _Pattern Recognition and Image Analysis_. Springer Berlin Heidelberg, 2009.
  * Lindley (1956) Lindley, D. V.  On a measure of the information provided by an experiment.  _The Annals of Mathematical Statistics_ , 27(4):986–1005, 1956.
  * Liu & Chen (1995) Liu, J. S. and Chen, R.  Blind deconvolution via sequential imputations.  _Journal of the American Statistical Association_ , 90(430):567–576, 1995.
  * Marcus et al. (1997) Marcus, S. I., Fernández-Gaucherand, E., Hernández-Hernandez, D., Coraluppi, S., and Fard, P.  Risk sensitive Markov decision processes.  In _Systems and Control in the Twenty-First Century_ , pp. 263–279, Boston, MA, 1997. Birkhäuser Boston.
  * McMichael et al. (2021) McMichael, R. D., Dushenko, S., and Blakley, S. M.  Sequential Bayesian experiment design for adaptive Ramsey sequence measurements.  _Journal of Applied Physics_ , 130(14):144401, 2021.
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E.  Equation of State Calculations by Fast Computing Machines.  _The Journal of Chemical Physics_ , 21(6):1087–1092, 1953.
  * Moffat et al. (2020) Moffat, H., Hainy, M., Papanikolaou, N. E., and Drovandi, C.  Sequential experimental design for predator–prey functional response experiments.  _Journal of the Royal Society, Interface_ , 17(166):20200156, 2020.
  * Myung et al. (2013) Myung, J. I., Cavagnaro, D. R., and Pitt, M. A.  A tutorial on adaptive design optimization.  _Journal of Mathematical Psychology_ , 57(3):53–67, 2013.
  * Naesseth et al. (2020) Naesseth, C., Lindsten, F., and Blei, D.  Markovian score climbing: Variational inference with KL(P||Q).  In _Advances in Neural Information Processing Systems_ , volume 33, pp. 15499–15510. Curran Associates, Inc., 2020.
  * Olsson & Westerborn (2017) Olsson, J. and Westerborn, J.  Efficient particle-based online smoothing in general hidden Markov models: The PaRIS algorithm.  _Bernoulli_ , 23(3):1951–1996, 2017.
  * Poole et al. (2019) Poole, B., Ozair, S., Van Den Oord, A., Alemi, A., and Tucker, G.  On variational bounds of mutual information.  In _International Conference on Machine Learning_ , 2019.
  * Puterman (2014) Puterman, M. L.  _Markov Decision Processes: Discrete Stochastic Dynamic Programming_.  John Wiley & Sons, 2014.
  * Rainforth et al. (2018) Rainforth, T., Cornish, R., Yang, H., Warrington, A., and Wood, F.  On nesting Monte Carlo estimators.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80, pp. 4267–4276, 2018.
  * Rainforth et al. (2023) Rainforth, T., Foster, A., Ivanova, D. R., and Smith, F. B.  Modern Bayesian experimental design.  _arXiv preprint arXiv:2302.14545_ , 2023.
  * Rawlik (2013) Rawlik, K. C.  _On Probabilistic Inference Approaches to Stochastic Optimal Control_.  PhD thesis, University of Edinburgh, 2013.
  * Robbins & Monro (1951) Robbins, H. and Monro, S.  A stochastic approximation method.  _The Annals of Mathematical Statistics_ , 22(3):400–407, 1951.
  * Särkkä & Solin (2019) Särkkä, S. and Solin, A.  _Applied Stochastic Differential Equations_.  Cambridge University Press, 2019.
  * Särkkä & Svensson (2023) Särkkä, S. and Svensson, L.  _Bayesian Filtering and Smoothing_.  Cambridge University Press, 2023.
  * Schultheis et al. (2020) Schultheis, M., Belousov, B., Abdulsamad, H., and Peters, J.  Receding horizon curiosity.  In _Conference on robot learning_ , pp. 1278–1288. PMLR, 2020.
  * Shababo et al. (2013) Shababo, B., Paige, B., Pakman, A., and Paninski, L.  Bayesian inference and online experimental design for mapping neural microcircuits.  In _Proceedings of the 26th International Conference on Neural Information Processing Systems_ , pp. 1304–1312. Curran Associates Inc., 2013.
  * Tedrake (2023) Tedrake, R.  _Underactuated Robotics_.  2023\.  URL https://underactuated.csail.mit.edu.  Course Notes for MIT 6.832, Accessed 21.12.2023.
  * Toussaint & Storkey (2006) Toussaint, M. and Storkey, A.  Probabilistic inference for solving discrete and continuous state Markov decision processes.  In _International Conference on Machine Learning_ , pp. 945–952, 2006.
  * van den Oord et al. (2019) van den Oord, A., Li, Y., and Vinyals, O.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_ , 2019.
  * Whitehouse et al. (2023) Whitehouse, M., Whiteley, N., and Rimella, L.  Consistent and fast inference in compartmental models of epidemics using Poisson approximate likelihoods.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 85(4):1173–1203, 2023.
  * Whittle (1990) Whittle, P.  A risk-sensitive maximum principle.  _Systems & Control Letters_, 15(3):183–192, 1990.
  * Wigren et al. (2019) Wigren, A., Risuleo, R. S., Murray, L., and Lindsten, F.  Parameter elimination in particle Gibbs sampling.  _Advances in Neural Information Processing Systems_ , 32, 2019.

##  Appendix A Proof of Proposition 1

###### Proof.

We start with the definition of the expected information gain in the terminal
reward framework,

$\displaystyle\mathcal{I}(\phi)$
$\displaystyle\coloneq\mathbb{E}_{p(z_{0:T}\mid\phi)}\Bigl{[}\mathbb{H}\bigl{[}p(\theta)\bigr{]}-\mathbb{H}\bigl{[}p(\theta\mid
z_{0:T})\bigr{]}\Bigr{]}$
$\displaystyle=\mathbb{E}_{p(z_{0:T},\theta\mid\phi)}\left[\log\frac{p(\theta\mid
z_{0:T})}{p(\theta)}\right].$

With repeated applications of Bayes’ rule, we can write this in an equivalent
form as

$\mathcal{I}(\phi)=\mathbb{E}_{p(z_{0:T},\theta\mid\phi)}\bigl{[}\log
p(z_{0:T}\mid\theta,\phi)\bigr{]}-\mathbb{E}_{p(z_{0:T}\mid\phi)}\bigl{[}\log
p(z_{0:T}\mid\phi)\bigr{]}.$ (23)

Let us look at the first term in (23). From (4), we know that the conditional
trajectory likelihood is

$p(z_{0:T}\mid\theta,\phi)=p(x_{0})\prod_{t=1}^{T}f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)\,\pi_{\phi}(\xi_{t-1}\mid z_{0:t-1}).$

We can then evaluate

$\displaystyle T_{1}$
$\displaystyle\coloneq\mathbb{E}_{p(z_{0:T},\theta\mid\phi)}\bigl{[}\log
p(z_{0:T}\mid\theta,\phi)\bigr{]}$
$\displaystyle=\mathbb{E}_{p(z_{0:T},\theta\mid\phi)}\left[\log
p(x_{0})+\sum_{t=1}^{T}\log f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)\right]+\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\log\pi_{\phi}(\xi_{t-1}\mid
z_{0:t-1})\right]$
$\displaystyle=-\mathbb{H}[p(x_{0})]+\mathbb{E}_{p(z_{0:T},\theta\mid\phi)}\left[\sum_{t=1}^{T}\log
f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)\right]+\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\log\pi_{\phi}(\xi_{t-1}\mid
z_{0:t-1})\right].$ (24)

For the second term in the above equation, we get

$\displaystyle\mathbb{E}_{p(z_{0:T},\theta\mid\phi)}\left[\sum_{t=1}^{T}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\right]$
$\displaystyle=\mathbb{E}_{p(\theta)p(z_{0:T}\mid\theta,\phi)}\left[\sum_{t=1}^{T}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\right]$
$\displaystyle=\mathbb{E}_{p(\theta)}\left[\sum_{t=1}^{T}\mathbb{E}_{p(z_{0:T}\mid\theta,\phi)}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\right]$
$\displaystyle=\mathbb{E}_{p(\theta)}\left[\sum_{t=1}^{T}\mathbb{E}_{p(z_{0:t}\mid\theta,\phi)}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\right]$
$\displaystyle=\sum_{t=1}^{T}\mathbb{E}_{p(z_{0:t},\theta\mid\phi)}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)$
$\displaystyle=\sum_{t=1}^{T}\mathbb{E}_{p(z_{0:t}\mid\phi)\,p(\theta\mid
z_{0:t})}\log f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)$
$\displaystyle=\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\mathbb{E}_{p(\theta\mid
z_{0:t})}\log f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\right]$
$\displaystyle\eqqcolon\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\alpha_{t}(z_{0:t})\right]$

where we have made use of the fact that the conditional dynamics is Markovian
in multiple places. We are then left with

$T_{1}=-\mathbb{H}[p(x_{0})]+\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\alpha_{t}(z_{0:t})\right]+\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\log\pi_{\phi}(\xi_{t-1}\mid
z_{0:t-1})\right].$ (25)

Let’s now look at the second term in (23). For that, we first need the
following factorization of the marginal trajectory distribution

$\displaystyle\log p(z_{0:T}\mid\phi)$
$\displaystyle=\log\left[p(z_{0})\prod_{t=1}^{T}p(z_{t}\mid
z_{0:t-1},\phi)\right]$
$\displaystyle=\log\left[p(x_{0})\prod_{t=1}^{T}p(x_{t},\xi_{t-1}\mid
z_{0:t-1},\phi)\right]$ $\displaystyle=\log
p(x_{0})+\sum_{t=1}^{T}\Bigl{[}\log p(x_{t}\mid
z_{0:t-1},\xi_{t-1})+\log\pi_{\phi}(\xi_{t-1}\mid z_{0:t-1})\Bigr{]},$

where

$\displaystyle p(x_{t}\mid z_{0:t-1},\xi_{t-1})$
$\displaystyle=\int_{\Theta}p(x_{t},\theta\mid z_{0:t-1},\xi_{t-1})\,d\theta$
$\displaystyle=\int_{\Theta}p(x_{t}\mid
z_{0:t-1},\xi_{t-1},\theta)\,p(\theta\mid z_{0:t-1},\xi_{t-1})\,d\theta$
$\displaystyle=\int_{\Theta}f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)\,p(\theta\mid z_{0:t-1})\,d\theta.$

In the last line, we have used the fact that $\theta$ is conditionally
independent of $\xi_{t-1}$ given $z_{0:t-1}$ (since the policy is independent
of $\theta$). We can now compute the second term of the EIG as

$\displaystyle T_{2}$
$\displaystyle\coloneq\mathbb{E}_{p(z_{0:T}\mid\phi)}\bigl{[}\log
p(z_{0:T}\mid\phi)\bigr{]}$
$\displaystyle=\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\log
p(x_{0})+\sum_{t=1}^{T}\Bigl{[}\log p(x_{t}\mid
z_{0:t-1},\xi_{t-1})+\log\pi_{\phi}(\xi_{t-1}\mid z_{0:t-1})\Bigr{]}\right]$
$\displaystyle=-\mathbb{H}[p(x_{0})]+\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\log
p(x_{t}\mid
z_{0:t-1},\xi_{t-1})\right]+\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\log\pi_{\phi}(\xi_{t-1}\mid
z_{0:t-1})\right].$ (26)

The full expression for $\mathcal{I}(\phi)$ is hence

$\displaystyle\mathcal{I}(\phi)$ $\displaystyle=T_{1}-T_{2}$
$\displaystyle=\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\alpha_{t}(z_{0:t})\right]-\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\log
p(x_{t}\mid z_{0:t-1},\xi_{t-1})\right]$
$\displaystyle=\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\Bigl{\\{}\alpha_{t}(z_{0:t})+\beta_{t}(z_{0:t})\Bigr{\\}}\right]$
(27)
$\displaystyle=\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}r_{t}(z_{0:t})\right]$
(28)

where we have defined $\beta_{t}(z_{0:t})=-\log p(x_{t}\mid
z_{0:t-1},\xi_{t-1})$ and
$r_{t}(z_{0:t})=\alpha_{t}(z_{0:t})+\beta_{t}(z_{0:t})$. This concludes the
proof for the first part of Proposition 1.

Let us now assume that our model has additive, constant noise in the dynamics
(noise that is independent of the state, design and $\theta$ parameters).
Under this assumption, the entropy $\mathbb{H}[f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)]$ is a constant. Let us now go back to the second
term in (24),

$\displaystyle\mathbb{E}_{p(z_{0:T},\theta\mid\phi)}\left[\sum_{t=1}^{T}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\right]$
$\displaystyle=\sum_{t=1}^{T}\mathbb{E}_{p(\theta)\,p(z_{0:t}\mid\theta,\phi)}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)$
$\displaystyle=\sum_{t=1}^{T}\mathbb{E}_{p(\theta)\,p(x_{0:t},\xi_{0:t-1}\mid\theta,\phi)}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)$
$\displaystyle=\sum_{t=1}^{T}\mathbb{E}_{p(\theta)\,p(x_{0:t-1},\xi_{0:t-1}\mid\theta,\phi)}\int_{\mathcal{X}}\log
f(x_{t}\mid x_{t-1},\xi_{t-1},\theta)\,f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)\,\mathrm{d}x_{t}$
$\displaystyle=\sum_{t=1}^{T}\mathbb{E}_{p(\theta)\,p(x_{0:t-1},\xi_{0:t-1}\mid\theta,\phi)}\bigl{[}-\mathbb{H}[f(x_{t}\mid
x_{t-1},\xi_{t-1},\theta)]\bigr{]}$ $\displaystyle=\textrm{constant}.$

Using ’$\equiv$’ to denote equality up to an additive constant, we now have

$T_{1}\equiv-\mathbb{H}[p(x_{0})]+\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\log\pi_{\phi}(\xi_{t-1}\mid
z_{0:t-1})\right],$

and hence,

$\displaystyle\mathcal{I}(\phi)$ $\displaystyle=T_{1}-T_{2}$
$\displaystyle\equiv\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}-\log
p(x_{t}\mid z_{0:t-1},\xi_{t-1})\right]$ (29)
$\displaystyle=\mathbb{E}_{p(z_{0:T}\mid\phi)}\left[\sum_{t=1}^{T}\beta_{t}(z_{0:t})\right]$
(30)

as required. ∎

##  Appendix B Proof of Proposition 2

We prove this result by induction over $t$. We note that this result is not
directly implied by the existing classical sequential Monte Carlo theory due
to the dependency of the likelihood term on the filtering distribution over
$\theta$.

We will make use of the following assumptions.

######  Assumption 3.

For all $z_{0:t+1},t\geq 1$ there exists $\alpha>0$ such that, for all
$\theta$, $0<f(x_{t+1}\mid x_{t},\xi_{t},\theta)<\alpha$.

######  Assumption 4.

For all $z_{0:t},t\geq 1$ there exists an integrable function
$\beta_{t+1}(z_{t+1})$ such that, for all $\theta$, $0<p(z_{t+1}\mid
z_{0:t},\theta)<\beta_{t+1}(z_{t+1})$.

###### Proof.

For simplicity, we assume that resampling happens at each step of IBIS in
Algorithm 2. The result is clear for $t=0$ under the law of large numbers. Now
assume it is true for a given $t$, then, applying Proposition 2 to the test
function $\theta\mapsto f(x_{t+1}\mid x_{t},\xi_{t},\theta)$ we have, thanks
to Assumption 3

$\frac{1}{M}\sum_{m=1}^{M}f(x_{t+1}\mid
x_{t},\xi_{t},\theta_{t}^{m})\to\mathbb{E}_{p(\theta_{t}\mid
z_{0:t})}\left[f(x_{t+1}\mid x_{t},\xi_{t},\theta_{t})\right]=p(x_{t+1}\mid
z_{0:t},\xi_{t})$ (31)

so that, with probability $1$,

$\exp\left\\{-\eta\log\frac{1}{M}\sum_{m=1}^{M}f(x_{t+1}\mid
x_{t},\xi_{t},\theta_{t}^{m})\right\\}\to\exp\left\\{-\eta\log p(x_{t+1}\mid
z_{0:t},\xi_{t})\right\\}.$ (32)

Applying the induction hypothesis to $\theta\mapsto p(z_{t+1}\mid
z_{0:t},\theta)$, we have with probability $1$,

$\frac{1}{M}\sum_{m=1}^{M}p(z_{t+1}\mid z_{0:t},\theta_{t}^{m})\to
p(z_{t+1}\mid z_{0:t}).$ (33)

As a consequence,

$\left[\frac{1}{M}\sum_{m=1}^{M}p(z_{t+1}\mid
z_{0:t},\theta_{t}^{m})\right]\,\exp\left\\{-\eta\log\frac{1}{M}\sum_{m=1}^{M}f(x_{t+1}\mid
x_{t},\xi_{t},\theta_{t}^{m})\right\\}\to p(z_{t+1}\mid
z_{0:t})\exp\left\\{-\eta\log p(x_{t+1}\mid z_{0:t},\xi_{t})\right\\}$ (34)

almost surely. Similarly, using the positivity of $p(z_{t+1}\mid
z_{0:t},\theta_{t}^{m})$ and by Lebesgue’s dominated convergence theorem, the
normalizing constant of the right-hand side of (18)

$\int\left[\frac{1}{M}\sum_{m=1}^{M}p(z_{t+1}\mid
z_{0:t},\theta_{t}^{m})\right]\,\exp\left\\{-\eta\log\frac{1}{M}\sum_{m=1}^{M}f(x_{t+1}\mid
x_{t},\xi_{t},\theta_{t}^{m})\right\\}\,\mathrm{d}z_{t+1}$ (35)

converges to

$\int p(z_{t+1}\mid z_{0:t})\,\exp\left\\{-\eta\log p(x_{t+1}\mid
z_{0:t},\xi_{t})\right\\}\mathrm{d}z_{t+1}$ (36)

and we have

$\frac{\Gamma^{M}_{t+1}(z_{0:t+1},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}\to\frac{p(z_{t+1}\mid
z_{0:t})\exp\left\\{-\eta\log p(x_{t+1}\mid z_{0:t},\xi_{t})\right\\}}{\int
p(z_{t+1}\mid z_{0:t})\exp\left\\{-\eta\log p(x_{t+1}\mid
z_{0:t},\xi_{t})\right\\}\mathrm{d}z_{t+1}}.$ (37)

The recursion is then obtained by noticing that the IBIS step (17) corresponds
to a particle filter update targeting $p(\theta\mid z_{0:t+1})$, so that,
under Assumption 4, we can follow Proposition 11.4 in Chopin &
Papaspiliopoulos (2020) to obtain that, for any bounded test function $\psi$,
writing

$\Gamma^{M,\mathrm{IBIS}}_{t}(\theta^{1:M}_{t})=\frac{\Gamma_{t}^{M}(z_{0:t},\theta_{0:t}^{1:M},a_{1:t}^{1:M})}{\Gamma^{M}_{t}(z_{0:t},\theta_{0:t-1}^{1:M},a_{1:t-1}^{1:M})},$

we have

$\mathbb{E}_{\Gamma^{M,\mathrm{IBIS}}_{t}}\left[\psi(\theta_{t})\right]\to\mathbb{E}_{p(\theta_{t}\mid
z_{0:t})}\left[\psi(\theta_{t})\right]$ (38)

almost surely. Putting it all together, and noticing that
$\Gamma(z_{0:t})p(\theta_{t}\mid z_{0:t})=\Gamma(z_{0:t},\theta_{t})$, we
obtain the result. ∎

##  Appendix C Details on the algorithms

###  C.1 Reweight step in IO-SMC2 for the general potential function

Algorithm 4 Reweight function corresponding to (7)

0: Any operation with superscript $m$ is to be understood as performed for all
$m=1,\dots,M$.

0: Reweight($t$)

1: $v_{t}^{mn}=f(x_{t}^{n}\mid x_{t-1}^{n},\xi_{t-1}^{n},\theta_{t-1}^{mn})$.

2: $W_{t,\theta}^{mn}\propto W_{t-1,\theta}^{mn}\,v_{t}^{mn}$.

3:
$r_{t}^{n}=\sum_{m=1}^{M}W_{t,\theta}^{mn}\,v_{t}^{mn}-\sum_{m=1}^{M}W_{t-1,\theta}^{mn}\,v_{t}^{mn}$.

4: return $g_{t}^{n}=\exp\\{\eta\,r_{t}^{n}\\}$.

4:

Let us reproduce the general expression for the stage reward at time $t+1$
from (7).

$\displaystyle r_{t}(z_{0:t+1})$ $\displaystyle=\mathbb{E}_{p(\theta\mid
z_{0:t+1})}\bigl{[}\log f(x_{t+1}\mid x_{t},\xi_{t},\theta)\bigr{]}-\log
p(x_{t+1}\mid z_{0:t},\xi_{t})$ $\displaystyle=\mathbb{E}_{p(\theta\mid
z_{0:t+1})}\bigl{[}\log f(x_{t+1}\mid
x_{t},\xi_{t},\theta)\bigr{]}-\log\mathbb{E}_{p(\theta\mid
z_{0:t})}\left[f(x_{t+1}\mid x_{t},\xi_{t},\theta)\right].$ (39)

We see that we have expectations with respect to the filtering posteriors of
$\theta$ at times $t$ and $t+1$. For a trajectory $z_{0:t+1}^{n}$, at
iteration $t$, after the IBIS step at line 6 of Algorithm 3, we have a
particle representation
$\sum_{m=1}^{M}W_{t,\theta}^{mn}\,\delta_{\theta_{t}^{mn}}(\theta)\approx
p(\theta\mid z_{0:t}^{n})$, and the second term of the reward function can be
estimated as

$\log\mathbb{E}_{p(\theta\mid z_{0:t}^{n})}\left[f(x_{t+1}^{n}\mid
x_{t}^{n},\xi_{t}^{n},\theta)\right]\approx\log\sum_{m=1}^{M}W_{t,\theta}^{mn}\,f(x_{t+1}^{n}\mid
x_{t}^{n},\xi_{t}^{n},\theta_{t}^{mn}).$

Now to compute the first term, we need to approximate $p(\theta\mid
z_{0:t+1}^{n})$. For this we perform the reweighting step of IBIS (line 2 from
Algorithm 2) to get updated weights

$\displaystyle W_{t+1,\theta}^{mn}$ $\displaystyle\propto
W_{t,\theta}^{mn}\,p_{\phi}(z_{t+1}^{n}\mid z_{0:t}^{n},\theta_{t}^{mn})$
$\displaystyle=W_{t,\theta}^{mn}\,f(x_{t+1}^{n}\mid
x_{t}^{n},\xi_{t}^{n},\theta_{t}^{mn})\,\pi_{\phi}(\xi_{t}^{n}\mid
z_{0:t}^{n})$ $\displaystyle\propto W_{t,\theta}^{mn}\,f(x_{t+1}^{n}\mid
x_{t}^{n},\xi_{t}^{n},\theta_{t}^{mn}).$

The distribution
$\sum_{m=1}^{M}W_{t+1,\theta}^{mn}\,\delta_{\theta_{t}^{mn}}(\theta)$ now
approximates the posterior $p(\theta\mid z_{0:t+1}^{n})$ as required. With
this, we can estimate the first term of (39) as

$\mathbb{E}_{p(\theta\mid z_{0:t+1}^{n})}\bigl{[}\log f(x_{t+1}^{n}\mid
x_{t}^{n},\xi_{t}^{n},\theta)\bigr{]}\approx\sum_{m=1}^{M}W_{t+1,\theta}^{mn}\,\log
f(x_{t+1}^{n}\mid x_{t}^{n},\xi_{t}^{n},\theta_{t}^{mn}).$

The entire reweighting procedure is outlined in Algorithm 4. We choose not to
perform the resample-move step before computing the above estimate so as not
to introduce additional variance.

###  C.2 Choice of the Markov kernel for IBIS

For the Markov kernel $K_{t}$ in IBIS (Algorithm 2), we follow the choice in
Chopin et al. (2013) and use a Metropolis-Hastings kernel (Metropolis et al.,
1953; Hastings, 1970). If the prior is Gaussian, we use a Gaussian random walk
proposal

$\tilde{\theta}^{m}\mid\theta^{m}\sim\mathcal{N}(\theta^{m},cI),$ (40)

where $c\in\mathbb{R}_{>0}$ is a constant that can be tuned to achieve a
desired acceptance ratio. In our experiments, we choose a value of $c$ that
gets the average acceptance ratio reasonably close to the optimal value of
0.234 for random walk Metropolis (Gelman et al., 1997). For log-normal priors,
we use a similar random walk proposal

$\tilde{\theta}^{m}\mid\theta^{m}\sim\mathrm{LogNormal}(\theta^{m},cI).$ (41)

We also perform multiple move steps per IBIS step to get a richer
representation of samples.

###  C.3 Inside-Out SMC2 for conditionally linear models with conjugate prior
and likelihood

Algorithm 5 Inside-Out SMC2 (Exact)

1: Sample $z_{0}^{n}\sim p(\cdot)$.

2: Sample $z_{1}^{n}\sim p(\cdot\mid z_{0}^{n})$ and initialize the state
history $z_{0:1}^{n}\leftarrow(z_{0}^{n},z_{1}^{n})$.

3: Compute and normalize the weights $W_{z}^{n}\propto g_{1}(z_{0:1}^{n})$.

4: for $t\leftarrow 1,\dots,T-1$ do

5: Sample $b_{t}^{n}\sim\mathcal{M}(W_{z}^{1:N})$.

6: Sample $z_{t+1}^{n}\sim p(\cdot\mid z_{0:t}^{b_{t}^{n}})$ and append to
state history $z_{0:t+1}^{n}\leftarrow(z_{0:t}^{b_{t}^{n}},z_{t+1}^{n})$.

7: Compute and normalize the weights $W_{z}^{n}\propto g_{t}(z_{0:t}^{n})$.

8: end for

9: return $\\{z_{0:T}^{n},W_{z}^{n}\\}_{n=1}^{N}$.

Let $\theta\in\mathbb{R}^{d_{\theta}}$, $x_{t}\in\mathbb{R}^{d_{x}}$ and
$\xi_{t}\in\mathbb{R}^{d_{\xi}}$. Let us consider conditionally linear,
Gaussian transition dynamics for $x$,

$f(x_{t+1}\mid x_{t},\xi_{t},\theta)=\mathcal{N}(x_{t+1}\mid
H(x_{t},\xi_{t})\,\theta,\Sigma(x_{t},\xi_{t})),$ (42)

where $H$ is a map $\mathbb{R}^{d_{x}+d_{\xi}}\mapsto\mathbb{R}^{d_{x}\times
d_{\theta}}$ and $\Sigma:\mathbb{R}^{d_{x}+d_{\xi}}\to\mathbb{R}^{d_{x}\times
d_{x}}$ maps to positive definite matrices. Let us also assume that at time
$t$, the filtered posterior of $\theta$ is also Gaussian with mean $m_{t}$ and
covariance matrix $P_{t}$,

$p(\theta\mid z_{0:t})=\mathcal{N}(\theta\mid m_{t},P_{t}).$ (43)

Then, using basic identities of the multivariate normal distribution, the
marginal density is

$\displaystyle p(x_{t+1}\mid z_{0:t},\xi_{t})$ $\displaystyle=\int
f(x_{t+1}\mid x_{t},\xi_{t},\theta)\,p(\theta\mid z_{0:t})\,\mathrm{d}\theta$
(44a) $\displaystyle=\mathcal{N}(x_{t+1}\mid Hm_{t},HP_{t}H^{T}+\Sigma),$
(44b)

where the functional dependence of $H$ and $\Sigma$ on $(x_{t},\xi_{t})$ has
been hidden for conciseness. Furthermore, upon observing the next augmented
state $z_{t+1}=(x_{t+1},\xi_{t})$, the filtered posterior of $\theta$ can be
updated using Bayes’ rule:

$\displaystyle p(\theta\mid z_{0:t+1})$
$\displaystyle=\frac{p_{\phi}(x_{t+1},\xi_{t}\mid
z_{0:t},\theta)\,p(\theta\mid z_{0:t})}{p_{\phi}(x_{t+1},\xi_{t}\mid
z_{0:t})}$ $\displaystyle=\frac{f(x_{t+1}\mid
x_{t},\xi_{t},\theta)\,\pi_{\phi}(\xi_{t}\mid z_{0:t})\,p(\theta\mid
z_{0:t})}{p(x_{t+1}\mid z_{0:t},\xi_{t})\,\pi_{\phi}(\xi_{t}\mid z_{0:t})}$
$\displaystyle=\frac{f(x_{t+1}\mid x_{t},\xi_{t},\theta)\,p(\theta\mid
z_{0:t})}{p(x_{t+1}\mid z_{0:t},\xi_{t})}$
$\displaystyle=\mathcal{N}(\theta\mid m_{t+1},P_{t+1}),$

where

$m_{t+1}=m_{t}+G(x_{t+1}-Hm_{t}),\quad P_{t+1}=P_{t}-GHP_{t},\quad
G=P_{t}H^{T}(HP_{t}H^{T}+\Sigma)^{-1}.$ (45)

Thus we see that for a conditionally linear model with Gaussian priors and
likelihoods, we can compute the marginal density and the $\theta$ posterior in
closed form. The same holds for any conjugate prior-likelihood pair. As a
result, we can compute the stage reward in (7) and hence the potential
function in closed form, and the resulting version of the IO-SMC2 algorithm
that does not use the inner particle filter is given in Algorithm 5.

###  C.4 Conditional SMC

In Section 4.3, we saw that IO-SMC2 is a nested particle filter that targets
the distribution $\Gamma_{T}^{M}$. In this section, we construct a conditional
version of our algorithm that keeps $\Gamma_{T}^{M}$ invariant. The basic idea
behind CSMC is that given a _reference trajectory_ from the target
distribution, at each time step in the forward pass, we sample $N-1$ samples
conditionally on the reference particle surviving the resampling step (Andrieu
et al., 2010). The CSMC kernel for $\Gamma_{t}^{M}$ is outlined in Algorithm
6, where the potential function estimates $g_{t}^{n}$ are computed as detailed
in Algorithm 4.

Algorithm 6 Conditional Inside-Out SMC2 kernel

0: Reference trajectory
$\\{z_{0:T},\\{W_{t,\theta}^{\bullet},\theta_{t}^{\bullet}\\}_{t=0}^{T-1}\\}$.

0: New reference trajectory
$\\{z_{0:T}^{*},\\{W_{t,\theta}^{\bullet*},\theta_{t}^{\bullet*}\\}_{t=0}^{T-1}\\}$.

1: Set $z_{0}^{1}\leftarrow z_{0}$, $\theta_{0}^{\bullet
1}\leftarrow\theta_{0}^{\bullet}$ and $W_{0,\theta}^{\bullet 1}\leftarrow
W_{0,\theta}^{\bullet}$.

2: for $n=2,\dots,N$ do

3: Sample $z_{0}^{n}\sim p(z_{0})$, $\theta_{0}^{\bullet n}\sim p(\theta)$ and
set $W_{0,\theta}^{\bullet n}\leftarrow 1/M$.

4: end for

5: Set $z_{0:1}^{1}\leftarrow z_{0:1}$.

6: for $n=2,\dots,N$ do

7: Sample $z_{1}^{n}\sim\hat{p}_{\phi}(\cdot\mid z_{0}^{n})$ and set
$z_{0:1}^{n}\leftarrow(z_{0}^{n},z_{1}^{n})$.

8: end for

9: Compute and normalize the weights $W_{z}^{n}\propto g_{1}^{n}$ for all
$n=1,\dots,N$.

10: for $t\leftarrow 1,\dots,T-1$ do

11: Set $z_{t+1}^{1}\leftarrow z_{t+1}$, $\theta_{t}^{\bullet
1}\leftarrow\theta_{t}^{\bullet}$, $W_{t,\theta}^{\bullet 1}\leftarrow
W_{t,\theta}^{\bullet}$, and $z_{0:t+1}^{1}\leftarrow z_{0:t+1}$.

12: for $n=2,\dots,N$ do

13: Sample $b_{t}^{n}\sim\mathcal{M}(W_{z}^{1:N})$.

14: $\theta_{t}^{\bullet n},W^{\bullet
n}_{t,\theta},\leftarrow\textsc{Ibis\\_Step}(z_{0:t}^{b_{t}^{n}},\theta_{t-1}^{\bullet
b_{t}^{n}},W^{\bullet b_{t}^{n}}_{t-1,\theta})$

15: Sample $z_{t+1}^{n}\sim\hat{p}_{\phi}(\cdot\mid z_{0:t}^{b_{t}^{n}}),$ and
append to state history
$z_{0:t+1}^{n}\leftarrow[z_{0:t}^{b_{t}^{n}},z_{t+1}^{n}]$.

16: end for

17: Compute and normalize the weights $W_{z}^{n}\propto g_{t}^{n}$ for all
$n=1,\dots,N$.

18: end for

19: Sample an index $j\sim\mathcal{M}(W_{z}^{1:N})$.

20: return $\\{z_{0:T}^{j},\\{W_{t,\theta}^{\bullet j},\theta_{t}^{\bullet
j}\\}_{t=0}^{T-1}\\}$.

While Algorithm 6 may look more complicated than “classical” CSMC algorithms
(Andrieu et al., 2010), its complexity may be abstracted away by remembering
that a ‘particle’ object is, in our case, an instance of the approximate inner
distribution $\Gamma_{t}^{M}$, which is associated with its weights, particles
and ancestors, noting that the ancestors do not appear in the computation of
Algorithm 4 and are therefore omitted.

##  Appendix D Experiment Details

###  D.1 Network Architectures and Hyperparameters

We use the same network architecture for all our experiments. The architecture
of our design policy network is similar to that in iDAD (Ivanova et al.,
2021), with an encoder network transforming the augmented states into a
64-dimensional representation stacked representation
$\\{R(z_{s})\\}_{s=0}^{t}$ before passing it to an LSTM cell with two hidden
layers. The encoder networks for the augmented states are given in Table 3 and
Table 4. The policy network is trained using the Adam optimizer (Kingma & Ba,
2014) with a constant learning rate that is specified for each experiment. The
policies were trained on a single 9th Gen Intel Core i9 processor.

An additional detail that deserves to be mentioned is that we add a slew rate
penalty on the designs to the weight function in IO-SMC2. This penalizes
extreme changes in the design from one experiment to the next, resulting in
the algorithm recovering smoother design trajectories. We found that this was
required for the policy network to learn.

For iDAD, we used the same policy architecture, but we used an Nvidia GeForce
3090 GPU for training. We use the code implementation made available by
Ivanova et al. (2021).

Layer | Description | Size | Activation
---|---|---|---
Input | Augmented state $z$ | Dim($z$) | -
Hidden layer 1 | Dense | 256 | ReLU
Hidden layer 2 | Dense | 256 | ReLU
Output | Dense | 64 | -
Table 3: The encoder architecture. Layer | Description | Size | Activation
---|---|---|---
Input |  $\\{R(z_{s})\\}_{s=0}^{t}$ $z$ | $64\cdot(t+1)$ | -
Hidden layer 1 | LSTM | 64 | -
Hidden layer 2 | LSTM | 64 | -
Hidden layer 3 | Dense | 256 | ReLU
Hidden layer 4 | Dense | 256 | ReLU
Output | Dense | Dim($\xi$) | -
Table 4: The LSTM architecture.

###  D.2 Stochastic Pendulum

We consider two different dynamics for the Compound pendulum. The following
problem specifications are the same for both settings. Let
$x_{t}=(q_{t},\dot{q}_{t})^{T}$ denote the state of the pendulum, with $q_{t}$
being the angle from the vertical and $\dot{q}_{t}$ the angular velocity at
time $t$. $m$ and $l$ are the mass and length of the pendulum respectively,
$g=9.81$ is the acceleration due to gravity and $d$ is the drag coefficient.
The design, $\xi_{t}$, in this experiment, is the torque applied to the
pendulum. The scale and shift parameters to the policy are $a=2.5$ and
$b=0.0$. We discretize the respective SDEs that describe the dynamics of each
system using Euler-Maruyama with a step size $dt=0.05$ and consider a horizon
of $T=50$ time steps. The initial state is fixed at $(0,0)$.

####  D.2.1 Conditionally Linear formulation

The unknown parameters are
$\theta=\displaystyle\left(\frac{3g}{2l},\frac{3d}{ml^{2}},\frac{3}{ml^{2}}\right)$.
The dynamics is described by the following Ito SDE (Särkkä & Solin, 2019,
Chapter 3)

$dx_{t}=h(x_{t},\xi_{t})^{T}\,\theta\,\mathrm{d}t+L\,\mathrm{d}\beta$ (46)

where $h(x_{t},\xi_{t})=(-\sin(q),-\dot{q},\xi_{t})^{T}$ (Belousov et al.,
2019), $L=(0,0.1)^{T}$ and $\beta$ denotes Brownian motion.

We use a Gaussian prior with mean $(10.0,0.0,5.0)$ and identity covariance
matrix. Other hyperparameters used for this experiment are given in Table 5.

Table 5: Hyperparameters for the conditionally linear pendulum experiment. Hyperparameter | IO-SMC2 | IO-SMC2 (Exact)
---|---|---
N | 256 | 256
M | 64 | -
Tempering | 0.5 | 0.5
Slew rate penalty | 0.1 | 0.1
IBIS moves | 1 | -
Learning rate | $10^{-3}$ | $10^{-3}$
Training iterations | 15 | 15

####  D.2.2 Non-linear formulation

The unknown parameters are $\theta=(m,l)$. The dynamics is described by the
SDE

$dx_{t}=h(x_{t},\xi_{t},\theta)\,\mathrm{d}t+L\,\mathrm{d}\beta,$

where $h(x_{t},\xi_{t},\theta)=(\dot{q}_{t},\ddot{q}_{t})^{T}$ with

$\ddot{q}_{t}=-\frac{3g}{2l}\sin{q_{t}}+\frac{(\xi_{t}-d\,\dot{q}_{t})}{ml^{2}},$

with $d=10^{-3}$ and $L=(0,0.1)^{T}$. The prior is a log-normal distribution

$p(\theta)=\textrm{LogNormal}\left(\begin{pmatrix}0\\\
0\end{pmatrix},\begin{pmatrix}0.25&0\\\ 0&0.25\end{pmatrix}\right).$

Other hyperparameters used for this experiment for our algorithm are given in
Table 6, and for iDAD in Table 7.

Table 6: Hyperparameters for the non-linear pendulum experiment for IO-SMC2. Hyperparameter | IO-SMC2
---|---
N | 256
M | 64
Tempering | 0.5
Slew rate penalty | 0.1
IBIS moves | 1
Learning rate | $10^{-3}$
Training iterations | 15
Table 7: Hyperparameters for the non-linear pendulum experiment for iDAD. Hyperparameter | iDAD
---|---
Batch size | 512
Number of contrastive samples | 1023
Number of gradient steps | 50000
Initial learning rate (LR) | $5\times 10^{-4}$
LR scheduler | Reduce on plateau
LR annealing factor | $0.96$
LR annealing frequency (if no improvement) | 2000

###  D.3 Stochastic Cart-Pole

We describe the state of a cart-pole system using the vector
$x_{t}=(s_{t},q_{t},\dot{s}_{t},\dot{q}_{t})^{T}$, where $s_{t}$ and
$\dot{s}_{t}$ are the position and velocity of the cart, and $q_{t}$ and
$\dot{q}_{t}$ are the position and velocity of the pole respectively (Tedrake,
2023). The unknowns are the mass $m_{p}$ and length $l$ of the pendulum,
$\theta=(m_{p},l)^{T}$. Other constants appearing in the equations are the
acceleration due to gravity, $g=9.81$, and the mass of the cart, $m_{c}=2.0$.
The design in this experiment is the force applied to the cart, with the scale
and shift parameters being $a=5.0$ and $b=0.0$. We consider the following SDE

$dx_{t}=h(x_{t},\xi_{t},\theta)\,\mathrm{d}t+L\,\mathrm{d}\beta,$

where
$h(x_{t},\xi_{t},\theta)=(\dot{s}_{t},\dot{q}_{t},\ddot{s}_{t},\ddot{q}_{t})^{T}$,
with

$\displaystyle\ddot{s}_{t}$
$\displaystyle=\frac{1}{m_{c}+m_{p}\sin^{2}q_{t}}\bigl{[}\xi_{t}+m_{p}\sin
q_{t}(l\dot{q}_{t}^{2}+g\cos{q_{t}}\bigr{]}$ $\displaystyle\ddot{q}_{t}$
$\displaystyle=\frac{1}{l(m_{c}+m_{p}\sin^{2}q_{t})}\bigl{[}-\xi_{t}\cos{q_{t}}-m_{p}l\dot{q}_{t}^{2}\cos{q_{t}}\sin{q_{t}}-(m_{c}+m_{p})g\sin{q_{t}}\bigr{]}.$

and $L=(0,0,0.1,0)^{T}$. We discretize the SDE with a step size $dt=0.05$ and
consider a horizon of $T=50$ time steps. The initial state is fixed at
$(0,0,0,0)$. The prior for $\theta$ is

$p(\theta)=\textrm{LogNormal}\left(\begin{pmatrix}0.2\\\
0.2\end{pmatrix},\begin{pmatrix}0.04&0\\\ 0&0.04\end{pmatrix}\right).$

Other hyperparameters used for this experiment are given in Table 8.

Table 8: Hyperparameters for the stochastic cart-pole experiment for IO-SMC2. Hyperparameter | IO-SMC2
---|---
N | 256
M | 128
Tempering | 0.25
Slew rate penalty | 0.1
IBIS moves | 7
Learning rate | $10^{-3}$
Training iterations | 10