# Supplementary Material

In this document, we provide proofs for the theoretical results in our paper
and additional details on the implementation and numerical experiments.

##  1 Particle EIG contrastive bound optimization

The algorithm below summarizes the procedure developed for design
optimization. It is the result of the combination of Algorithms 1 and 2 in the
main body. It results from the replacement in the former of the current
posterior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ by its particle
approximation resulting from the adaptive tempering in Algorithm 2.

1 Set $T$ iterations, ${\boldsymbol{\xi}}_{0}$, stepsizes
$(\gamma_{t})_{t=1:T}$

2 Run tempered SMC Algorithm 2 (main body) at step $k$ with multinomial
resampling, and set $L+1$ independent particle approximations
$P_{k,\ell}^{N}$, for $\ell=0:L$ (see main paper Section 5 for details)

3 while  _$t\leq T$_ do

4 Sample ${\boldsymbol{\theta}}^{i}_{\ell,t}\\!\sim P_{k,\ell}^{N},\mbox{for
$\ell\\!=\\!0\\!:\\!L,i\\!=\\!1\\!:\\!N_{t}$ }$

5 Sample $u^{i}_{t}\sim p(u),\quad\mbox{ for $i=1:N_{t}$}$

6 Set
$\\!\nabla_{\\!t+1}\\!\\!=\\!\\!\frac{1}{N_{t}}\\!\\!\sum\limits_{i=1}^{N_{t}}\\!\\!\nabla_{\boldsymbol{\xi}}\\!F({\boldsymbol{\xi}},\\!{\boldsymbol{\theta}}^{i}_{0t},\cdot,\\!{\boldsymbol{\theta}}^{i}_{Lt},\\!T^{{\boldsymbol{\xi}}}_{{\boldsymbol{\theta}}^{i}_{0t}}\\!\\!(u^{i}_{t})\\!)|_{{\boldsymbol{\xi}}={\boldsymbol{\xi}}_{t}}\\!\\!$

7 Update
${\boldsymbol{\xi}}_{t+1}={\boldsymbol{\xi}}_{t}+\gamma_{t}{\nabla}_{t+1}$

8 end while

return ${\boldsymbol{\xi}}_{k+1,N}^{*}\\!=\\!{\boldsymbol{\xi}}_{T}$ or a
Polyak averaging value

Algorithm 1 Particle EIG contrastive bound stochastic optimization at step
$k+1$

##  2 Differentiation under the integral sign

We briefly recall standard conditions under which it is possible to exchange
differentiation and expectation operators. In our setting this implies
conditions that are also useful for other results in Section 3.3. Indeed to
carry out a stochastic gradient algorithm (Algorithm 1 above and in the
paper), we assume that all quantities are well defined. In particular,

$\displaystyle\nabla_{\\!{\boldsymbol{\xi}}}\\!I_{\\!PCE}({\boldsymbol{\xi}})\\!=\\!\mathbb{E}_{p(u)\prod_{\ell=0}^{L}p({\boldsymbol{\theta}}_{\ell})}\\!\left[\nabla_{\\!{\boldsymbol{\xi}}}F({\boldsymbol{\xi}},\\!{\boldsymbol{\theta}}_{0},\cdot,\\!{\boldsymbol{\theta}}_{L},\\!T^{{\boldsymbol{\xi}}}_{{\boldsymbol{\theta}}_{0}}\\!(U))\right].$
(1)

Denoting $p_{u}\otimes p_{\theta}^{\otimes L+1}$ the product probability
distribution $p(u)\prod_{\ell=0}^{L}p({\boldsymbol{\theta}}_{\ell})$ on ${\cal
U}\times{\boldsymbol{\Theta}}^{L+1}$, sufficient conditions for (1) are that
the function $F$ is differentiable in ${\boldsymbol{\xi}}$ for $p_{u}\otimes
p_{\theta}^{\otimes L+1}$-almost all
$(u,{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})$ and its
gradient satisfies
$\left|\nabla_{\\!{\boldsymbol{\xi}}}F({\boldsymbol{\xi}},\\!{\boldsymbol{\theta}}_{0},\cdot,\\!{\boldsymbol{\theta}}_{L},\\!T^{{\boldsymbol{\xi}}}_{{\boldsymbol{\theta}}_{0}}\\!(u))\right|\leq
H({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L},u)$ where $H$ has
a bounded expectation. In particular this implies the continuity properties
$(ii)$ and $(iii)$ required in Lemma 4 and 5 below.

In practice, it can also be useful to notice that the gradient required within
the expectation (1) can be expressed using only the log-likelihood gradient.
To simplify, in the expression below we denote
$T^{{\boldsymbol{\xi}}}_{{\boldsymbol{\theta}}_{0}}(u)$ by $y$. Using the
definition of $F$, it comes,

$\displaystyle\nabla_{\\!{\boldsymbol{\xi}}}F({\boldsymbol{\xi}},\\!{\boldsymbol{\theta}}_{0},\cdot,\\!{\boldsymbol{\theta}}_{L},y)$
$\displaystyle=\nabla_{\\!{\boldsymbol{\xi}}}\log
p(y|{\boldsymbol{\xi}},{\boldsymbol{\theta}}_{0})-\sum_{\ell=0}^{L+1}w_{\ell}({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L},{\boldsymbol{\xi}},y)\;\nabla_{\\!{\boldsymbol{\xi}}}\log
p(y|{\boldsymbol{\xi}},{\boldsymbol{\theta}}_{\ell})$ (2)

where the
$w_{\ell}({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L},{\boldsymbol{\xi}},y)=\frac{p(y|{\boldsymbol{\xi}},{\boldsymbol{\theta}}_{\ell})}{\sum\limits_{\ell^{\prime}=0}^{L}p(y|{\boldsymbol{\xi}},{\boldsymbol{\theta}}_{\ell^{\prime}})}$
sum to 1 and act as weights. Thus we only need to compute for all
${\boldsymbol{\theta}}_{0},\theta_{\ell}$, the gradient
$\nabla_{\\!{\boldsymbol{\xi}}}\log
p(T^{{\boldsymbol{\xi}}}_{{\boldsymbol{\theta}}_{0}}(u)|{\boldsymbol{\xi}},{\boldsymbol{\theta}}_{\ell})$.

##  3 Proofs of main and intermediate results

We first specify the steps that lead to Proposition 1 in the main body. The
overall goal is to show that products of particle approximations have similar
convergence properties than single particle approximations.

###  3.1 Product form estimators

Product form estimators have been used in various places in the literature but
only studied in a more formal and general way in some recent work [13]. Their
advantages are clearly highlighted in [13] but the theoretical results therein
do not cover the use of product form estimators within SMC samplers. In
contrast, [17, 15, 12, 1] consider SMC settings that may look similar to ours
but differ in key aspects. Indeed, [17] proposed a so-called block particle
filter consisting of a product of local smaller particle filters to
approximate particle filters in high dimensional models. The last line of
their Algorithm 2 clearly shows that the block particle filter is equivalent
to running several independent particle filters on smaller dimensional spaces.
In contrast, we run a single SMC from which we extract at each step a product
form estimate of the target distribution. This allows to compute expectations
and other inferential quantities in a more efficient way, using the advantages
of product form estimators, notably their lower variance, as explained in
[13]. The Divide-and-Conquer SMC (DaC-SMC) methodology introduced in [15] and
further studied in [12] targets a different objective, which is to generalize
the sequential nature of standard SMC to a more general tree-structured
execution flow of the sampler. Hierarchical decompositions of variables of
interest into subsets are considered and product form estimators are used to
recombined the splitted subsets of particles. Typical use cases include sets
of variables exhibiting a graph dependence structure such as Markov random
fields. In contrast, in our work, the execution flow remains sequential,
corresponding to a tree that is a simple line, and we do not address complex
parameter structures. At last, the Tensor Monte Carlo technique proposed in
[1] shares similarity with DaC-SMC but focuses on computational aspects and
does not cover as many theoretical properties. The work in [1] proposes to
exploit dependence structures in variables and conditional independence
properties to address the issue of computing the very large sums induced by
product form estimators. This is not an issue we encounter as, in our work,
the computation of very large sums is avoided via stochastic approximation.

###  3.2 Convergence of product form particle approximations

The result is intuitive and generalizes standard SMC proofs but as mentioned
in the previous section, we believe it is not covered in previous work. It
requires additional care as illustrated in Lemma 1 below. We follow the
presentation in Chapter 11 of [5] which is mainly based on the presentation in
[6]. As mentioned in [5] it is not the most general presentation but it has
the advantage of using standard tools while presenting the general idea. For
simplicity, we only present the L2 convergence result, which is enough to
prove Proposition 1.

We use $\tau$ as the time or step index and $\mu_{\tau}$, $M_{\tau}$ for
$,\mu_{\lambda_{\tau}},M_{\lambda_{\tau}}$ as in Algorithm 2, to denote the
step and successive targets and Markov kernels, but this extends to $k$ and
$p_{k}$, as we repeat the same process sequentially. Recall that for a SMC
procedure on ${\boldsymbol{\Theta}}$, with $G_{\tau}$ denoting the so-called
generator involved in the weights normalization, for a function $\phi$ on
${\boldsymbol{\Theta}}$, we have

$\displaystyle\mu_{\tau}[\phi]=\mathbb{E}_{\mu_{\tau}}[\phi]$
$\displaystyle=\mu_{\tau-1}M_{\tau-1}[\phi\;G_{\tau}]=\int
G_{\tau}({\boldsymbol{\theta}})\phi({\boldsymbol{\theta}})\left\\{\int
M_{\tau-1}({\boldsymbol{\theta}}^{\prime},{\boldsymbol{\theta}})\;\mu_{\tau-1}({\boldsymbol{\theta}}^{\prime})\;d{\boldsymbol{\theta}}^{\prime}\right\\}d{\boldsymbol{\theta}}.$
(3)

As a particular case, we have

$\displaystyle\mu_{\tau}({\boldsymbol{\theta}}_{\ell})=\mu_{\tau}[\delta_{{\boldsymbol{\theta}}_{\ell}}]=G_{\tau}({\boldsymbol{\theta}}_{\ell})\left\\{\int
M_{\tau-1}({\boldsymbol{\theta}}^{\prime},{\boldsymbol{\theta}}_{\ell})\;\mu_{\tau-1}({\boldsymbol{\theta}}^{\prime})\;d{\boldsymbol{\theta}}^{\prime}\right\\}\;d{\boldsymbol{\theta}}$
(4) $\displaystyle\mu_{\tau-1}M_{\tau-1}[G_{\tau}]=\mu_{\tau}[1]=1$ (5)

We generalize these definitions and results on the product space
${\boldsymbol{\Theta}}^{L+1}$ with $\phi$ being now a function of $L+1$
variables. For a function $\phi$ on ${\boldsymbol{\Theta}}^{L+1}$, we use
notation $M_{\tau}\phi$ for the function

$\displaystyle
M_{\tau}\phi:({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})\rightarrow\int\cdot\int\phi({\boldsymbol{\theta}}^{\prime}_{0},\cdot,{\boldsymbol{\theta}}^{\prime}_{L})\left(\prod_{\ell=0}^{L}M_{\tau}({\boldsymbol{\theta}}_{\ell},{\boldsymbol{\theta}}^{\prime}_{\ell})\;d{\boldsymbol{\theta}}^{\prime}_{\ell}\right).$
(6)

Similarly, we denote $\mu_{\tau}^{\otimes
L+1}[\phi]=\mathbb{E}_{\mu_{\tau}^{\otimes
L+1}}\\!\left[\phi\right]=\int\cdot\int\phi({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})\left(\prod_{\ell=0}^{L}\mu_{\tau}({\boldsymbol{\theta}}_{\ell})\;d{\boldsymbol{\theta}}_{\ell}\right).$
Defining,

$\displaystyle(\mu_{\tau-1}M_{\tau-1})^{\otimes L+1}[\phi]$
$\displaystyle=\int\cdot\int\phi({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})\prod\limits_{\ell=0}^{L}\left\\{\int
M_{\tau-1}({\boldsymbol{\theta}}^{\prime}_{\ell},{\boldsymbol{\theta}}_{\ell})\;\mu_{\tau-1}({\boldsymbol{\theta}}^{\prime}_{\ell})\;d{\boldsymbol{\theta}}^{\prime}_{\ell}\right\\}d{\boldsymbol{\theta}}_{\ell}.$

we can deduce from (3) that

$\displaystyle\mu_{\tau}^{\otimes L+1}[\phi]=(\mu_{\tau-1}M_{\tau-1})^{\otimes
L+1}\\!\left[\phi\prod_{\ell=0}^{L}G_{\tau}({\boldsymbol{\theta}}_{\ell})\right],$
(7) $\displaystyle(\mu_{\tau-1}M_{\tau-1})^{\otimes
L+1}\\!\left[\prod_{\ell=0}^{L}G_{\tau}({\boldsymbol{\theta}}_{\ell})\right]=\mu_{\tau}^{\otimes
L+1}[1]=1.$ (8)

The strategy is to establish Proposition 1 by decomposing the error at step
$\tau$ as a sum of contributions from each step (sampling, reweighting,
resampling, etc.) of the current and previous steps. We first need a general
result on variances of sum of random variables.

######  Lemma 1.

Let $\\{X_{\ell}^{i}\\}_{\ell=0:L,i=1:N}$ a collection of i.i.d. random
variables on ${\cal X}$. For any function $\phi$ on ${\cal X}^{L+1}$ so that
the variances exist we have, for all $(i_{0},\cdot,i_{L})\in[1:N]^{L+1}$,

$Var\left(\phi(X_{0}^{i_{0}},\cdot,X_{L}^{i_{L}})\right)=Var\left(\phi(X_{0}^{1},\cdot,X_{L}^{1})\right)$

and we denote this common variance by $Var(\phi_{1})$. We can then bound the
following variance

$\displaystyle
Var\left(\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi(X_{0}^{i_{0}},\cdot,X_{L}^{i_{L}})\right)\leq
N^{L+1}\left(N^{L+1}-(N-1)^{L+1}\right)Var(\phi_{1})$

#### Proof.

Let us denote $I=[1:N]$, $J=I^{L+1}$, $j=(i_{0},\cdot,i_{L})$ for $i_{\ell}\in
I$, and use the short notation
$\phi_{j}=\phi(X_{0}^{i_{0}},\cdot,X_{L}^{i_{L}})$. Then

$\displaystyle
Var\left(\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi(X_{0}^{i_{0}},\cdot,X_{L}^{i_{L}})\right)$
$\displaystyle=Var(\sum\limits_{j\in J}\phi_{j})=\sum_{j\in
J}Var(\phi_{j})+\sum\limits_{j\in
J}\sum_{j^{\prime}\not=j}Cov(\phi_{j^{\prime}},\phi_{j})$
$\displaystyle=N^{L+1}Var(\phi_{1})+\sum\limits_{j\in
J}\sum_{j^{\prime}\not=j}Cov(\phi_{j^{\prime}},\phi_{j})\;.$

The second covariance term is a sum of $N^{2(L+1)}-N^{L+1}$ pairwise
covariances among which a number are zero due to the independence of the
$\\{X_{\ell}^{i}\\}_{\ell=0:L,i=1:N}$. To count them, we notice that for each
$j=(i_{0},\cdot,i_{L})\in J$, $Cov(\phi_{j},\phi_{j}^{\prime})=0$ for every
$j^{\prime}\in J$ such that $j^{\prime}=(i^{\prime}_{0},\cdot,i^{\prime}_{L})$
with $i^{\prime}_{\ell}\in I\backslash\\{i_{\ell}\\}$. There are thus
$(N-1)^{L+1}$ such $j^{\prime}$. It follows that there are altogether
$N^{L+1}(N-1)^{L+1}$ null covariance terms and
$N^{L+1}(N^{L+1}-(N-1)^{L+1}-1)$ non null ones. For each non null covariance,
the Cauchy-Schwartz inequality implies that
$Cov(\phi_{j^{\prime}},\phi_{j})\leq\sqrt{Var(\phi_{j^{\prime}})Var(\phi_{j})}=Var(\phi_{1})$,
which concludes the proof. $\square$

We then recall some notation and definitions specific to our product form
particle approximation. The tempered SMC procedure in Algorithm 2 produces
$M=N(L+1)$ particles denoted by ${\boldsymbol{\theta}}_{\tau}^{1:M}$ at each
step indexed by $\tau$. When partitioning these particles into $L+1$ disjoint
subset we will denote by ${\boldsymbol{\theta}}_{\ell,\tau}^{1:N}$ the $N$
particles in each subset with $\ell=0:L+1$. We denote by $\zeta_{k,L}^{N}$ the
collection of random variables
$\zeta_{k,L}^{N}=\\{W_{k,\ell}^{i},{\boldsymbol{\theta}}^{i}_{k,\ell}\\}_{i=1:N,\ell=0:L}$,
and use the short notation $p_{k}^{\otimes L+1}[\phi]$ for $p_{k}^{\otimes
L+1}[\phi]=\mathbb{E}_{\prod_{\ell=0}^{L}p({\boldsymbol{\theta}}_{\ell}|{\mathb{D}}_{k})}[\phi({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})]$.
We write ${\cal F}_{\tau-1}$ for the $\sigma$-algebra ${\cal
F}_{\tau-1}=\sigma\\!\left({\boldsymbol{\theta}}^{1:M}_{\tau-1}\right)$. Let
${\cal C}_{b}({\boldsymbol{\Theta}}^{L+1})$ denote the set of functions
$\phi:{\boldsymbol{\Theta}}^{L+1}\rightarrow\mathbb{R}$ that are measurable
and bounded, and let $||\phi||_{\infty}$ denote the supremum norm
$||\phi||_{\infty}=\sup_{{\boldsymbol{\theta}}\in{\boldsymbol{\Theta}}^{L+1}}|\phi({\boldsymbol{\theta}})|$.
The following lemma bounds the Monte Carlo error at step $\tau$.

######  Lemma 2 (Monte Carlo error).

Using Algorithm 2 with a multinomial resampling scheme, for all functions
$\phi\in{\cal C}_{b}({\boldsymbol{\Theta}}^{L+1})$,

$\displaystyle\mathbb{E}_{\zeta_{k,L}^{N}}\\!\left[\left\\{\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\\!\\!\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\\!\left(\prod_{\ell=0}^{L}W_{\ell,\tau-1}^{i_{\ell}}\right)\\!M_{\tau-1}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau-1},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau-1})\right\\}^{2}\right]\\!\\!\leq\\!c_{N,L}\|\phi\|_{\infty}^{2}$
(9)

with $c_{N,L}=1-\left(1-\frac{1}{N}\right)^{L+1}$ and where the expectation is
taken over all the realizations of the random tempered SMC method, or
equivalently on $\zeta_{k,L}^{N}$.

#### Proof.

Multinomial resampling preserves conditional independence so that the
particles ${\boldsymbol{\theta}}_{\ell,\tau}^{1:N}$ are i.i.d. conditionally
on ${\cal F}_{\tau-1}$ and their conditional distribution is

${\boldsymbol{\theta}}^{i}_{\ell,\tau}|{\cal
F}_{\tau-1}\sim\sum_{i=1}^{N}W^{i}_{\ell,\tau-1}M_{\tau-1}({\boldsymbol{\theta}}_{\ell,\tau-1}^{i},\cdot).$

Thus for all $i_{0},\cdot,i_{L}$,

$\displaystyle\mathbb{E}_{\zeta_{k,L}^{N}}\\!\left[\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})\left|\right.{\cal
F}_{\tau-1}\right]$
$\displaystyle=\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod_{\ell=0}^{L}W_{\ell,\tau-1}^{i_{\ell}}\right)M_{\tau-1}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau-1},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau-1})$
$\displaystyle=\mathbb{E}_{\zeta_{k,L}^{N}}\\!\left[\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})\left|\right.{\cal
F}_{\tau-1}\right].$

It follows that

$\displaystyle\mathbb{E}_{\zeta_{\tau,L}^{N}}\\!\left[\\!\left\\{\frac{1}{N^{L+1}}\\!\\!\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})\\!-\\!\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\\!\\!\left(\prod_{\ell=0}^{L}W_{\ell,\tau-1}^{i_{\ell}}\right)M_{\tau-1}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau-1},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau-1})\right\\}^{2}\\!\left|\right.{\cal
F}_{\tau-1}\\!\right]$
$\displaystyle=Var\\!\left[\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})\left|\right.{\cal
F}_{\tau-1}\right]\;=\;\frac{1}{N^{2(L+1)}}Var\\!\left[\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})\left|\right.{\cal
F}_{\tau-1}\right]$ $\displaystyle\leq
c_{N,L}\;Var\\!\left[\phi({\boldsymbol{\theta}}^{1}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{1}_{L,\tau})\left|\right.{\cal
F}_{\tau-1}\right]$ $\displaystyle\leq
c_{N,L}\;\mathbb{E}_{\zeta_{\tau,L}^{N}}\\!\left[\phi({\boldsymbol{\theta}}^{1}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{1}_{L,\tau})^{2}\left|\right.{\cal
F}_{\tau-1}\right]\leq c_{N,L}\;\|\phi\|^{2}_{\infty}$

The line before last above results from Lemma 1. Using the tower property, we
get (9). $\square$

The constant $c_{N,L}$ is equivalent to $\frac{L+1}{N}$ and tends to 0 when
$N$ tends to $\infty$. The next Lemma involves the so-called Generator
$G_{k,\tau}(\theta)=p(\mathbf{y}_{k}|\theta,{\boldsymbol{\xi}})^{\gamma}$
which is involved in the normalization of the weights. Note that $\gamma$
usually also depends on $\tau$ and $k$ as it is found adaptively by solving
equation (11) in the paper. However all we need is that the generator and then
the likelihood is upper bounded.

######  Lemma 3 (Weights normalization error).

If $G_{k,\tau}$ is upper bounded, then for all functions $\phi\in{\cal
C}_{b}({\boldsymbol{\Theta}}^{L+1})$,

$\displaystyle\mathbb{E}_{\zeta_{\tau,L}^{N}}\\!\left[\left\\{\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau}^{i_{\ell}}\right)\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\\!\left(\prod_{\ell=0}^{L}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell,\tau}^{i_{\ell}})\right)\\!\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})\right\\}^{2}\right]$
$\displaystyle\leq\|\phi\|_{\infty}^{2}\;\mathbb{E}_{\zeta_{k,L}^{N}}\\!\left[\left\\{\prod\limits_{\ell=0}^{L}\left(\frac{1}{N}\sum_{i=1}^{N}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell,\tau}^{i})\right)-1\right\\}^{2}\right]$

#### Proof.

By definition

$W_{\ell,\tau}^{i_{\ell}}=\frac{{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell,\tau}^{i_{\ell}})}{\sum_{i=1}^{N}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell,\tau}^{i})}$

Thus the term in the square in the left hand side is

$\left(\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau}^{i_{\ell}}\right)\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})\right)\left(1-\frac{1}{N^{L+1}}\prod\limits_{\ell=0}^{L}(\sum_{i=1}^{N}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell,\tau}^{i}))\right)$

and the first factor is bounded by $\|\phi\|_{\infty}$ since
$\sum\limits_{i_{0}=1}^{N}\cdot\sum\limits_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau}^{i_{\ell}}\right)=1$.
$\square$

The two Lemmas can now be combined to establish the proposition below.

######  Proposition 1 (L2 convergence).

Using Algorithm 2 with multinomial resampling and assuming that potential
functions ${G}_{k,\tau}$ are upper bounded, there exists constants
$c_{\tau},c^{\prime}_{\tau}$ so that, for all functions $\phi\in{\cal
C}_{b}({\boldsymbol{\Theta}}^{L+1})$,

$\displaystyle\mathbb{E}_{\zeta_{k,L}^{N}}\\!\left[\left\\{\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\\!\\!(\mu_{\tau-1}M_{\tau-1})^{\otimes
L+1}[\phi]\right\\}^{2}\right]\\!\\!\leq\\!c_{\tau}\;c_{N,L}\;\|\phi\|_{\infty}^{2}$
(10)
$\displaystyle\mathbb{E}_{\zeta_{k,L}^{N}}\\!\left[\left\\{\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau}^{i_{\ell}}\right)\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\\!\\!\mu_{\tau}^{\otimes
L+1}[\phi]\right\\}^{2}\right]\\!\\!\leq\\!c^{\prime}_{\tau}\;c_{N,L}\;\|\phi\|_{\infty}^{2}$
(11)

with $c_{N,L}=1-(1-\frac{1}{N})^{L+1}$.

#### Proof

The proof works by induction on $\tau$. Assume (10) holds at step $\tau$, we
first show that (11) also holds at step $\tau$. The left hand side in (11) can
be decomposed into

$\displaystyle\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau}^{i_{\ell}}\right)\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\mu_{\tau}^{\otimes
L+1}[\phi]=$
$\displaystyle\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau}^{i_{\ell}}\right)\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\\!\left(\prod_{\ell=0}^{L}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell,\tau}^{i_{\ell}})\right)\\!\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})$
$\displaystyle+$
$\displaystyle\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\\!\left(\prod_{\ell=0}^{L}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell,\tau}^{i_{\ell}})\right)\\!\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\mu_{\tau}^{\otimes
L+1}[\phi]\;.$

For all $\phi$, when
$\phi^{\prime}({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})=\left(\prod_{\ell=0}^{L}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell})\right)\\!\phi({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})$,
we have $(\mu_{\tau-1}M_{\tau-1})^{\otimes
L+1}[\phi^{\prime}]=\mu_{\tau}^{\otimes L+1}[\phi]$ by (7). The second term
can then be bounded by $c_{\tau}c_{N,L}\|\phi^{\prime}\|_{\infty}^{2}$,
applying (10) with $\phi^{\prime}$. For the first term, we apply Lemma 3 and
then (10) with
$\phi({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})=\prod_{\ell=0}^{L}{G}_{k,\tau}({\boldsymbol{\theta}}_{\ell})$
since then $(\mu_{\tau-1}M_{\tau-1})^{\otimes L+1}[\phi]=\mu_{\tau}^{\otimes
L+1}[1]=1$ by (8). Finally we get (11) with
$c^{\prime}_{\tau}=4c_{\tau}\|G\|^{2(L+1)}_{\infty}$ where we have used that
$\mathbb{E}[(X+Y)^{2}]\leq 2(\mathbb{E}[X^{2}]+\mathbb{E}[Y^{2}])$. We then
show that (11) at step $\tau-1$ implies (10) at step $\tau$. Similarly, we
have,

$\displaystyle\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-(\mu_{\tau-1}M_{\tau-1})^{\otimes
L+1}[\phi]=$
$\displaystyle\frac{1}{N^{L+1}}\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau})-\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau-1}^{i_{\ell}}\right)M_{\tau-1}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau-1},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau-1})$
$\displaystyle+$
$\displaystyle\sum_{i_{0}=1}^{N}\cdot\sum_{i_{L}=1}^{N}\left(\prod\limits_{\ell=0}^{L}W_{\ell,\tau-1}^{i_{\ell}}\right)M_{\tau-1}\phi({\boldsymbol{\theta}}^{i_{0}}_{0,\tau-1},\cdot,{\boldsymbol{\theta}}^{i_{L}}_{L,\tau-1})-(\mu_{\tau-1}M_{\tau-1})^{\otimes
L+1}[\phi]\;.$

The first term can be bounded using Lemma 2 and for the second term by
applying (11) at $\tau-1$ to function $M_{\tau-1}\phi$ defined in (6). Then
(10) holds with $c_{\tau}=2(1+c^{\prime}_{\tau-1})$. $\square$

Proposition 1 in the paper is just the result above but reduced to the
presentation of (11).

###  3.3 Consistency of the sequential design estimators

We specify the steps that lead to the proof of Proposition 2 in the main body.
We first recall some notation and definitions. The sequential design values
produced by Algorithm 1 above can be seen as realizations of random estimators
$\\{{\boldsymbol{\xi}}_{k+1,N}^{*}\\}_{N\geq 1}$ targeting points of maximum
of random criterion functions
$\\{I_{PCE}^{k+1,N}({\boldsymbol{\xi}})\\}_{N\geq 1}$,

$\displaystyle I^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}})$
$\displaystyle=\mathbb{E}_{p(u)\prod_{\ell=0}^{L}{p}^{N}_{k,\ell}({\boldsymbol{\theta}}_{\ell})}\\!\left[F({\boldsymbol{\xi}},\\!{\boldsymbol{\theta}}_{0},\cdot,\\!{\boldsymbol{\theta}}_{L},\\!T^{{\boldsymbol{\xi}}}_{{\boldsymbol{\theta}}_{0}}\\!(U))\right]\;.$

A natural question is to study the limiting distributions of these random
quantities when the number of particles $N$ tends to infinity. The main
property is the following Theorem 1, which results from a general result from
M-estimator theory [19, 20]. The two following lemmas are then useful to
provide simpler sufficient conditions to satisfy the Theorem’s assumptions and
to establish our main result in Proposition 2 in the main body, also recalled
below.

######  Theorem 1 (Theorem 5.7 in [19]).

Assume

  * (i) For all $\epsilon>0$,

$\displaystyle\lim_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(\sup_{{\boldsymbol{\xi}}\in{\cal
E}}|{I}^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}})-I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})|\geq\epsilon\right)=0.$
(12)

  * (ii) For all $\epsilon>0$,

$\displaystyle\sup_{\|{\boldsymbol{\xi}}-{\boldsymbol{\xi}}^{*}_{k+1}\|\geq\epsilon}I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})<I^{k+1}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1})\;.$
(13)

  * (iii) There exists a sequence of positive random variables $\\{\rho_{N}\\}_{N\geq 1}$ and a sequence of random variables $\\{{\boldsymbol{\xi}}^{*}_{k+1,N}\\}_{N\geq 1}$ in ${\cal E}$ that satisfy

$\forall\epsilon>0,\quad\lim_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(\rho_{N}\geq\epsilon\right)=0$
$\lim\inf_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left({I}^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1,N})\geq
I^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1})-\rho_{N}\right)=1.$

Then the sequence of estimators $\\{{\boldsymbol{\xi}}^{*}_{k+1,N}\\}_{N\geq
1}$ is consistent, i.e. for all $\epsilon>0$,

$\lim_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(\|{\boldsymbol{\xi}}^{*}_{k+1,N}-{\boldsymbol{\xi}}^{*}_{k+1}\|\geq\epsilon\right)=0.$

#### Proof.

The proof is a special case of Theorem 5.7 in [19]. We reproduce it using our
notation. In all the following proofs, to simplify, we drop the $k+1$
notation, so that $I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})$ and
$I^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}})$ are now simply denoted by
$I_{\\!PCE}({\boldsymbol{\xi}})$ and $I^{N}_{\\!PCE}({\boldsymbol{\xi}})$ and
their respective maximizers by $\\{{\boldsymbol{\xi}}^{*}\\}$ and
$\\{{\boldsymbol{\xi}}^{*}_{N}\\}$. Since ${\boldsymbol{\xi}}^{*}$ is the
maximizer of $I_{\\!PCE}({\boldsymbol{\xi}})$, it comes that for all $N\geq
1$,

$\displaystyle 0$ $\displaystyle\leq
I_{\\!PCE}({\boldsymbol{\xi}}^{*})-I_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})$
$\displaystyle=(I_{\\!PCE}({\boldsymbol{\xi}}^{*})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*}))+(I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N}))+(I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})-I_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N}))\;.$

The first and third terms in the sum are bounded by
$\sup\limits_{{\boldsymbol{\xi}}\in{\cal
E}}|I_{\\!PCE}({\boldsymbol{\xi}})-I^{N}_{\\!PCE}({\boldsymbol{\xi}})|$ while
the second term can be bounded by
$\rho_{N}+(I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})-\rho_{N})\;\delta_{\\{I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})>\rho_{N}\\}}$
where $\delta_{{\cal A}}$ is the indicator function which is 1 if ${\cal A}$
is satisfied and 0 otherwise. It follows that

$\displaystyle 0$ $\displaystyle\leq
I_{\\!PCE}({\boldsymbol{\xi}}^{*})-I_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})$
$\displaystyle\leq 3\max\\!\left(\\!2\sup\limits_{{\boldsymbol{\xi}}\in{\cal
E}}|I_{\\!PCE}({\boldsymbol{\xi}})-I^{N}_{\\!PCE}({\boldsymbol{\xi}})|,\rho_{N},\left(I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})-\rho_{N}\right)\\!\delta_{\\{I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})>\rho_{N}\\}}\\!\\!\right).$

Assumptions $(i)$ and $(iii)$ imply that all three terms in the max tend to 0
in probability so that for all $\eta>0$,

$\displaystyle\lim_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(I_{\\!PCE}({\boldsymbol{\xi}}^{*})-I_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})\geq\eta\right)=0.$
(14)

Then, for all $\epsilon>0$, using $(ii)$, for all ${\boldsymbol{\xi}}$
satisfying $\|{\boldsymbol{\xi}}-{\boldsymbol{\xi}}^{*}\|\geq\epsilon$, there
exists $\eta>0$ so that $I_{\\!PCE}({\boldsymbol{\xi}})\leq
I_{\\!PCE}({\boldsymbol{\xi}}^{*})-\eta$. This implies that
$\\{\|{\boldsymbol{\xi}}^{*}_{N}-{\boldsymbol{\xi}}^{*}\|\geq\epsilon\\}\subset\\{I_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})\leq
I_{\\!PCE}({\boldsymbol{\xi}}^{*})-\eta\\}$ and

$p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(\|{\boldsymbol{\xi}}^{*}_{N}-{\boldsymbol{\xi}}^{*}\|\geq\epsilon\right)\leq
p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(I_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})\leq
I_{\\!PCE}({\boldsymbol{\xi}}^{*})-\eta\right)=p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(I_{\\!PCE}({\boldsymbol{\xi}}^{*})-I_{\\!PCE}({\boldsymbol{\xi}}^{*}_{N})\geq\eta\right).$

The limit in $N$ of this last term tends to 0 by (14), which concludes the
proof. $\square$

In general, the uniform convergence in $(i)$ is the most difficult assumption
to check but in our setting, when ${\cal E}$ is compact, it is easily derived
from previous assumptions and the pointwise convergence in probability ($(iv)$
below), which can be derived from the L2 convergence in Proposition 1 (see
comments in the main body). Recall the following definition,

$f_{PCE}({\boldsymbol{\xi}},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})=\mathbb{E}_{p(u)}\left[F({\boldsymbol{\xi}},\\!{\boldsymbol{\theta}}_{0},\cdot,\\!{\boldsymbol{\theta}}_{L},\\!T^{{\boldsymbol{\xi}}}_{{\boldsymbol{\theta}}_{0}}\\!(U))\right]$

and the shortcut notation $p_{k}^{\otimes
L+1}=\prod\limits_{\ell=0}^{L}p({\boldsymbol{\theta}}_{\ell}|{\mathb{D}}_{k})$.

######  Lemma 4.

Assume

  * (i) ${\cal E}\in\mathbb{R}^{d}$ is a compact set.

  * (ii) $I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})$ is a continuous function in ${\boldsymbol{\xi}}$.

  * (iii) $f_{PCE}({\boldsymbol{\xi}},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})$ is a continuous function in ${\boldsymbol{\xi}}$ for $p_{k}^{\otimes L+1}$-almost all $({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})$.

  * (iv) Convergence in probability pointwise: For all ${\boldsymbol{\xi}}\in{\cal E}$ and all $\epsilon>0$,

$\lim_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(|{I}^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}})-I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})|\geq\epsilon\right)=0.$

Then the uniform convergence in Theorem 1 $(i)$ is satisfied, that is, for all
$\epsilon>0$,

$\displaystyle\lim_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(\sup_{{\boldsymbol{\xi}}\in{\cal
E}}|{I}^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}})-I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})|\geq\epsilon\right)=0.$

#### Proof.

Continuous functions on a compact set are uniformly continuous. It follows
from $(ii)$ and $(iii)$ that for all $\epsilon>0$, there exists $\eta>0$ so
that for all ${\boldsymbol{\xi}}^{\prime}\in{\cal E}$,

$\sup_{\|{\boldsymbol{\xi}}-{\boldsymbol{\xi}}^{\prime}\|\leq\eta}|I_{PCE}({\boldsymbol{\xi}})-I_{PCE}({\boldsymbol{\xi}}^{\prime})|\leq\epsilon$
(15)

and for $p_{k}^{\otimes L+1}$-almost all
$({\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})$,

$\sup_{\|{\boldsymbol{\xi}}-{\boldsymbol{\xi}}^{\prime}\|\leq\eta}|f_{PCE}({\boldsymbol{\xi}},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})-f_{PCE}({\boldsymbol{\xi}}^{\prime},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})|\leq\epsilon.$
(16)

Let ${\cal B}({\boldsymbol{\xi}},\eta)$ be a ball centered at
${\boldsymbol{\xi}}$ with radius $\eta$. As ${\cal E}$ is compact, for all
$\eta>0$, it is possible to extract, from the cover set
$\bigcup\limits_{{\boldsymbol{\xi}}\in{\cal E}}{\cal
B}({\boldsymbol{\xi}},\eta)$, a finite subcover $\bigcup\limits_{b=1:B}{\cal
B}({\boldsymbol{\xi}}^{(b)},\eta)$ so that ${\cal
E}\subset\bigcup\limits_{b=1:B}{\cal B}({\boldsymbol{\xi}}^{(b)},\eta)$ and

$\displaystyle\sup_{{\boldsymbol{\xi}}\in{\cal
E}}|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}})-I_{\\!PCE}({\boldsymbol{\xi}})|$
$\displaystyle=$
$\displaystyle\max\limits_{b=1:B}\sup\limits_{{\boldsymbol{\xi}}\in{\cal
B}({\boldsymbol{\xi}}^{(b)},\eta)}|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}})-I_{\\!PCE}({\boldsymbol{\xi}})|\;.$
(17)

For all $b=1:B$, and all ${\boldsymbol{\xi}}\in{\cal
B}({\boldsymbol{\xi}}^{(b)},\eta)$, we also have that

$|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}})-I_{\\!PCE}({\boldsymbol{\xi}})|\leq|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{(b)})|+|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}}^{(b)})-I_{\\!PCE}({\boldsymbol{\xi}}^{(b)})|+|{I}_{\\!PCE}({\boldsymbol{\xi}}^{(b)})-I_{\\!PCE}({\boldsymbol{\xi}})|$

which implies

$\displaystyle\sup\limits_{{\boldsymbol{\xi}}\in{\cal
B}({\boldsymbol{\xi}}^{(b)},\eta)}|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}})-I_{\\!PCE}({\boldsymbol{\xi}})|$
$\displaystyle\leq$ $\displaystyle\sup\limits_{{\boldsymbol{\xi}}\in{\cal
B}({\boldsymbol{\xi}}^{(b)},\eta)}\left\\{|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{(b)})|\right\\}$
$\displaystyle+$
$\displaystyle|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}}^{(b)})-I_{\\!PCE}({\boldsymbol{\xi}}^{(b)})|$
$\displaystyle+$ $\displaystyle\sup\limits_{{\boldsymbol{\xi}}\in{\cal
B}({\boldsymbol{\xi}}^{(b)},\eta)}\left\\{|{I}_{\\!PCE}({\boldsymbol{\xi}}^{(b)})-I_{\\!PCE}({\boldsymbol{\xi}})|\right\\}$

For the first term in the right-hand side,

$\sup\limits_{{\boldsymbol{\xi}}\in{\cal
B}({\boldsymbol{\xi}}^{(b)},\eta)}\\!\\!\left\\{|{I}^{N}_{\\!PCE}({\boldsymbol{\xi}})-I^{N}_{\\!PCE}({\boldsymbol{\xi}}^{(b)})|\right\\}\leq\mathbb{E}_{p_{k}^{N\otimes
L+1}}\\!\\!\left[\sup\limits_{\|{\boldsymbol{\xi}}-{\boldsymbol{\xi}}^{(b)}\|\leq\eta}|f_{PCE}({\boldsymbol{\xi}},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})-f_{PCE}({\boldsymbol{\xi}}^{(b)},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})|\right]$

By Proposition 1 (11), the right-hand side above tends in L2-norm and then in
probability when $N$ tends to $\infty$ to

$\mathbb{E}_{p_{k}^{\otimes
L+1}}\\!\\!\left[\sup\limits_{\|{\boldsymbol{\xi}}-{\boldsymbol{\xi}}^{(b)}\|\leq\eta}|f_{PCE}({\boldsymbol{\xi}},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})-f_{PCE}({\boldsymbol{\xi}}^{(b)},{\boldsymbol{\theta}}_{0},\cdot,{\boldsymbol{\theta}}_{L})|\right],$

which is smaller than $\epsilon$ by (16). The second term in the right-hand
side of (3.3) tends in probability to 0 by $(iv)$ for all
${\boldsymbol{\xi}}^{(b)}$. Finally, using (15), the third term in (3.3) is
smaller than $\epsilon$. To conclude, (17) implies the uniform convergence
(12). $\square$

The following result gives simpler sufficient conditions for Assumption $(ii)$
in Theorem 1 to hold.

######  Lemma 5.

Assume

  * (i) ${\cal E}\in\mathbb{R}^{d}$ is a compact set.

  * (ii) $I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})$ is a continuous function in ${\boldsymbol{\xi}}$ .

  * (iii) For all ${\boldsymbol{\xi}}\not={\boldsymbol{\xi}}^{*}_{k+1}$, $I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})<I^{k+1}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1})$ .

Then Assumption $(ii)$ of Theorem 1 is satisfied.

#### Proof.

A continuous function reaches its maximum on a compact set. It follows that
$I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})$ reaches its maximum on the compact
subset ${\cal E}\backslash{\cal B}({\boldsymbol{\xi}}^{*}_{k+1},\epsilon)$.
Let ${\boldsymbol{\xi}}_{\epsilon}$ denote a value at which this maximum is
reached. Then $(iii)$ implies that
$I^{k+1}_{\\!PCE}({\boldsymbol{\xi}}_{\epsilon})<I^{k+1}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1})$,
which leads to (13). $\square$

Using Lemmas 4 and 5, we can now replace in Theorem 1 Assumptions $(i)$ and
$(ii)$ by simpler conditions. It follows the Proposition 2 presented in the
paper and recalled below.

######  Proposition 2.

Assume

######  A 1.

${\cal E}\in\mathbb{R}^{d}$ is a compact set.

######  A 2.

For all ${\boldsymbol{\xi}}\not={\boldsymbol{\xi}}^{*}_{k+1}$,
$I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})<I^{k+1}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1})$

######  A 3.

There exists a sequence of positive random variables $\\{\rho_{N}\\}_{N\geq
1}$ and a sequence of random variables
$\\{{\boldsymbol{\xi}}^{*}_{k+1,N}\\}_{N\geq 1}$ in ${\cal E}$ that satisfy

$\forall\epsilon>0,\quad\lim\limits_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(\rho_{N}\geq\epsilon\right)=0$
$\lim\inf\limits_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left({I}^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1,N})\geq
I^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}}^{*}_{k+1})-\rho_{N}\right)=1.$

Then the sequence of estimators $\\{{\boldsymbol{\xi}}^{*}_{k+1,N}\\}_{N\geq
1}$ is consistent, i.e. for all $\epsilon>0$,

$\lim\limits_{N\rightarrow\infty}p_{\zeta_{k,L}^{N}}\\!\\!\\!\left(\|{\boldsymbol{\xi}}^{*}_{k+1,N}-{\boldsymbol{\xi}}^{*}_{k+1}\|\geq\epsilon\right)=0.$

#### Proof.

With ${\cal E}$ compact, we can use Lemma 4. The continuity of
$I^{k+1}_{\\!PCE}({\boldsymbol{\xi}})$ and $f_{PCE}$ has been already assumed
earlier as specified in Section 2 and Lemma 4 $(iv)$ is a consequence of
Proposition 1. It follows the uniform convergence property $(i)$ in Theorem 1.
Then (A1-2) and Lemma 5 imply $(ii)$ in Theorem 1. With (A3), Theorem 1 leads
to the result. Note that if we assume that ${\boldsymbol{\xi}}^{*}_{k+1,N}$ is
the exact maximizer of ${I}^{k+1,N}_{\\!PCE}({\boldsymbol{\xi}})$ then (A3) is
trivially satified with $\rho_{N}=0$.

##  4 Numerical experiments

###  4.1 Sequential prior contrastive estimation (SPCE) criterion

We specify the SPCE introduced by [9] and used in our experiments and those of
[3] to assess the design sequence quality in our comparison. For a number $K$
of experiments,
${\mathb{D}}_{K}=\\{(\mathbf{y}_{1},{\boldsymbol{\xi}}_{1}),\cdot,(\mathbf{y}_{K},{\boldsymbol{\xi}}_{K})\\}$
and $L$ contrastive variables, SPCE is defined as

$\displaystyle SPCE({\boldsymbol{\xi}}_{1},\cdot,{\boldsymbol{\xi}}_{K})$
$\displaystyle=\mathbb{E}_{\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\xi}}_{k},{\boldsymbol{\theta}}_{0})\;\prod\limits_{\ell=0}^{L}p({\boldsymbol{\theta}}_{\ell})}\left[\log\frac{\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}}_{0},{\boldsymbol{\xi}}_{k})}{\frac{1}{L+1}\sum\limits_{\ell=0}^{L}\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}}_{\ell},{\boldsymbol{\xi}}_{k})}\right]\;.$
(19)

SPCE is a lower bound of the total EIG which is the expected information
gained from the entire sequence of design parameters
${\boldsymbol{\xi}}_{1},\ldots,{\boldsymbol{\xi}}_{K}$ and it becomes tight
when $L$ tends to $\infty$. In addition, SPCE has the advantage to use only
samples from the prior $p({\boldsymbol{\theta}})$ and not from the successive
posterior distributions. It makes it a fair criterion to compare methods on
design sequences only. Considering a true parameter value denoted by
${\boldsymbol{\theta}}^{*}$, given a sequence of design values
$\\{{\boldsymbol{\xi}}_{k}\\}_{k=1:K}$, observations
$\\{\mathbf{y}_{k}\\}_{k=1:K}$ are simulated using
$p(\mathbf{y}|{\boldsymbol{\theta}}^{*},{\boldsymbol{\xi}}_{k})$ respectively.
Therefore, for a given ${\mathb{D}}_{k}$, the corresponding SPCE is estimated
numerically by sampling
${\boldsymbol{\theta}}_{1},\cdot,{\boldsymbol{\theta}}_{L}$ from the prior,

$\displaystyle SPCE({\mathb{D}}_{K})$
$\displaystyle=\frac{1}{N}\sum\limits_{i=1}^{N}\left\\{\log\frac{\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}}^{*},{\boldsymbol{\xi}}_{k})}{\frac{1}{L+1}\left(\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}}^{*},{\boldsymbol{\xi}}_{k})+\sum\limits_{\ell=1}^{L}\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}}^{i}_{\ell},{\boldsymbol{\xi}}_{k})\right)}\right\\}\;.$

As shown in [9] (Appendix A), SPCE increases with $L$ to reach the total EIG
$I({\boldsymbol{\xi}}_{1},\ldots,{\boldsymbol{\xi}}_{K})$ when
$L\rightarrow\infty$ at a rate ${\cal O}(L^{-1})$ of convergence. More
specifically, it is shown in [9] that

$\displaystyle 0\leq
I({\boldsymbol{\xi}}_{1},\ldots,{\boldsymbol{\xi}}_{K})-SPCE({\boldsymbol{\xi}}_{1},\ldots,{\boldsymbol{\xi}}_{K})\leq\frac{C}{L+1}$
(20)

where
$C=\mathbb{E}_{p({\mathb{D}}_{K})p({\boldsymbol{\theta}}|{\mathb{D}}_{K})}\left[\frac{p({\mathb{D}}_{K}|{\boldsymbol{\theta}})}{p({\mathb{D}}_{K})}\right]-1$
with the notation
$p({\mathb{D}}_{K}|{\boldsymbol{\theta}})=\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}},{\boldsymbol{\xi}}_{k})$.

It is also shown in [9] that for a given $L$, SPCE is bounded by $\log(L+1)$
while the upper bound SNMC below is potentially unbounded. As in [3], if we
use $L=10^{7}$ to compute SPCE and SNMC, the bound is $\log(L+1)=16.12$ for
SPCE. In practice this does not impact the numerical methods comparison as the
intervals [SPCE, SNMC] containing the total EIG remain clearly distinct.

###  4.2 Sequential nested Monte Carlo (SNMC) criterion

Similarly, an upper bound on the total EIG, with similar features, has also
been introduced by [9]. Its expression is very similar to that of SPCE,
varying only through the sum in the denominator,

$\displaystyle SNMC({\boldsymbol{\xi}}_{1},\cdot,{\boldsymbol{\xi}}_{K})$
$\displaystyle=\mathbb{E}_{\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\xi}}_{k},{\boldsymbol{\theta}}_{0})\;\prod\limits_{\ell=0}^{L}p({\boldsymbol{\theta}}_{\ell})}\left[\log\frac{\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}}_{0},{\boldsymbol{\xi}}_{k})}{\frac{1}{L}\sum\limits_{\ell=1}^{L}\prod\limits_{k=1}^{K}p(\mathbf{y}_{k}|{\boldsymbol{\theta}}_{\ell},{\boldsymbol{\xi}}_{k})}\right]\;.$

###  4.3 Implementation details

For VPCE [10] and RL-BOED [3], we use the code available at github.com/csiro-
mlai/RL-BOED, using the settings recommended therein to reproduce the results
in the respective papers. From the obtained sequences of observations and
design values, we compute SPCE and SNMC as explained above and retrieve the
same results as in their respective papers.

Our code is implemented in Jax [4] and provided as a supplementary zip file.
Several packages are used through the repository. Namely, we used Optax [2] to
run Gradient Descents, the Sequential Monte Carlo part was heavily inspired
and built using Kernels from BlackJax [14] and we used OTT [7] to compute
Wasserstein distances.

Table 1 summarizes the main features and running times of the compared
methods. The RL-BOED method has the advantage to be both non-myopic and
amortized in the sense that a policy is learnt upfront and then used
straightforwardly at each new experiment. It follows a much longer training
time, which does not exist for the other methods. Note that in comparison the
deployment times of all methods are neglible (see Table 1 in the paper). In
contrast RL-BOED does not provide approximations for the posterior
distributions.

Method | Posterior | Amortized | Non-myopic | Training Time
---|---|---|---|---
PASOA | ✓ | ✗ | ✗ | —
SMC | ✓ | ✗ | ✗ | —
RL-BOED [3] | ✗ | ✓ | ✓ |  | CES: $\sim 20$h
---
Sources: $\sim 10$h
VPCE [10] | ✓ | ✗ | ✗ | —
Table 1: Main features and training times of the compared methods: the second
column indicates whether a method also provides approximation of posterior
distribution, the third if it is amortized and the fourth if it is non-myopic.
The last column shows training times for the amortized method RL-BOED and a
sequence of $K$ experiments run on a single Nvidia V100 GPU, for the source
finding and CES examples.

###  4.4 Hardware details

Our method can be run on a local machine and was tested on a Apple M1 Pro 16Gb
chip. However, for a faster running time, each experiment was finally produced
by running our method on a single Nvidia V100 GPU. One other advantage of
tempering and of our PASOA method is that by reducing the number of needed
particles for an accurate procedure, it lowers the hardware requirements for
this method as it becomes feasible to run it on CPUs.

###  4.5 Source location example

For the 2D location finding example used in [9, 3] and tested in the paper,
with 2 sources, $K=30$ successive design optimisations, and 100 repetitions of
the experiment, the number of gradient steps was set to 5000 and the ESS for
the SMC procedure to 0.9. Figure 3 in the paper shows the SPCE, SNMC and the
L2 Wasserstein distances between weighted particles and the true source
locations, providing three quantitative assessment and comparison of methods.
As an additional, visual assessment of the quality of the posterior
approximation provided by our method, Figure 1 in the paper and Figure 1 below
illustrate the evolution of the particles over the design steps, starting from
a sample following the prior to a sample concentrating around the true source
locations. In particular, the $k=0$ step shows particles simulated according
to the prior. In most use cases, plain SMC already gives better results than
other reference methods. Figure 1 and Table 1 in the paper show that tempering
allows to reduce the number of particles. In Figure 1 below, the source
locations, indicated by red crosses in the plots, are chosen in a part of the
space not well covered by the prior to illustrate the robustness of our
approach to a potential prior misspecification. After some iterations PASOA is
able to explore the parameter space to finally concentrate the posterior on
the true source locations. In contrast, SMC may miss some of the sources when
they are outside the prior mass or when there are too many of them. We suspect
that prior misspecification is a typical very common feature that jeopardizes
SMC performance while impacting much less PASOA. This is actually the same
problem encountered with IS that tempering aims at solving. Similarly, Figure
3 shows, in terms of SPCE, SNMC and Wasserstein distance, that SMC is more
robust with tempering. Figure 2 indicates the number of tempering steps taken
on average. The median (over 100 rollouts) number of tempering steps varies
from 14 to 2 and is globally decreasing, being under 5 after 15 experiments.
The number of tempering steps reduces when the posterior concentrates and when
adding new observations becomes less informative.

Figure 1: PASOA evolution of particles (in purple) over some selected steps
$k$. Particles correspond initially to a sample from the prior
$p({\boldsymbol{\theta}})$ and progressively evolve to a sample of particles
located around the initially unknown true source positions indicated by red
crosses. Green crosses indicate the optimal measurement locations
${\boldsymbol{\xi}}_{k}^{*}$ obtained at each step $k$.  Figure 2: Source
location example: median (over 100 rollouts) number of tempering steps with
respect to the number of particles.

|  |
---|---|---

Figure 3: Source location example. Prior misspecification: median (over 100
rollouts) (a) SPCE, (b) SNMC and (c) Wasserstein distances for SMC (blue) and
PASOA (red).

###  4.6 Constant Elasticity of Substitution example

This other model, used in [3, 10], comes from behavioral economics. In this
model, an agent compares two baskets of goods and gives a rating $y$ on a
sliding scale from 0 to 1. The goal is to design the two baskets of goods so
as to infer the agent’s utility function, which depends on some unknown
parameters. The designs are 6-dimensional vectors
${\boldsymbol{\xi}}=({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2})$
corresponding to the two baskets with 3 values each
${\boldsymbol{\xi}}_{d}=({\boldsymbol{\xi}}_{d,1},{\boldsymbol{\xi}}_{d,2},{\boldsymbol{\xi}}_{d,3})\in[0:100]^{3}$
for $d=1,2$, which represent quantities for 3 items in each basket. There are
3 parameters ${\boldsymbol{\theta}}=(\rho,{\boldsymbol{\alpha}},u)$ in
dimension 5, whose prior distributions are respectively $\rho\sim Beta(1,1)$,
${\boldsymbol{\alpha}}=(\alpha_{1},\alpha_{2},\alpha_{3})\sim
Dirichlet(1,1,1)$ and $\log u\sim{\cal N}(1,3)$.

The model likelihood is given by the following model that uses a subjective
utility function $U$ and two hyperparameters $\epsilon=2^{-22}$ and
$\tau=0.005$,

$\displaystyle y$ $\displaystyle=f(\eta,\epsilon)$ where
$\displaystyle\eta\sim{\cal N}(\mu_{\eta},\sigma_{\eta}^{2})$ with
$\displaystyle\mu_{\eta}=(U({\boldsymbol{\xi}}_{1})-U({\boldsymbol{\xi}}_{2}))\;u$
$\displaystyle\sigma_{\eta}$
$\displaystyle=(1+\|{\boldsymbol{\xi}}_{1}-{\boldsymbol{\xi}}_{2}\|)\;\tau\;u$
$\displaystyle\mbox{For $d=1,2$,}\quad U({\boldsymbol{\xi}}_{d})$
$\displaystyle=\left(\alpha_{1}{\boldsymbol{\xi}}_{d,1}^{\rho}+\alpha_{2}{\boldsymbol{\xi}}_{d,2}^{\rho}+\alpha_{3}{\boldsymbol{\xi}}_{d,3}^{\rho}\right)^{1/\rho}$

where $f(\eta,\epsilon)$ takes it values in $[\epsilon,1-\epsilon]$ and is a
censored sigmoid defined by

$\displaystyle f(\eta,\epsilon)$ $\displaystyle=$ $\displaystyle
1-\epsilon\quad\mbox{if $\eta\geq logit(1-\epsilon)$}$ $\displaystyle=$
$\displaystyle\epsilon\quad\mbox{if $\eta\leq logit(\epsilon)$ }$
$\displaystyle=$ $\displaystyle(1+\exp(-\eta))^{-1}\quad\mbox{otherwise}$

with $logit(y)=\log(y/(1-y))$. In other words, $y$ is a censored logit-normal
distribution with parameters $\mu_{\eta}$ and $\sigma_{\eta}$. Its density is
a mixture

$p(y|{\boldsymbol{\theta}},{\boldsymbol{\xi}})=p_{0}\delta_{\epsilon}(y)+p_{1}\delta_{1-\epsilon}(y)+(1-p_{0}-p_{1})\;q(y|{\boldsymbol{\theta}},{\boldsymbol{\xi}})$

where
$q(y|{\boldsymbol{\theta}},{\boldsymbol{\xi}})=\frac{1}{\sigma_{\eta}\sqrt{2\pi}y(1-y)}\exp(\frac{(logit(y)-\mu_{\eta})^{2}}{2\sigma_{\eta}})$
is the density of a logit-normal distribution and $p_{0}$ and $p_{1}$ are
defined by the following logit-normal CDF values

$\displaystyle p_{0}$ $\displaystyle=$ $\displaystyle
q(y\leq\epsilon)=p(\eta\leq
logit(\epsilon))=F\left(\frac{logit(\epsilon)-\mu_{\eta}}{\sigma_{\eta}}\right)$
(21) $\displaystyle p_{1}$ $\displaystyle=$ $\displaystyle 1-q(y\leq
1-\epsilon)=1-p(\eta\leq
logit(1-\epsilon))=1-F\left(\frac{logit(1-\epsilon)-\mu_{\eta}}{\sigma_{\eta}}\right)$
(22)

with the last equalities involving the normal CDF values of variable $\eta$
and the standard normal CDF $F$ values. In practice, computing $\log p_{0}$
and $\log p_{1}$ may sometimes be numerically problematic when $p_{0}$ or
$p_{1}$ become too small. Computing $p_{0}$ or $p_{1}$ involves computing
lower and upper Gaussian tails. In this case, following [10], we use the
following first order asymptotic approximation of the standard normal CDF,
when $x$ is large,

$1-F(x)\approx\frac{1}{x\;\sqrt{2\pi}}\;\exp(-x^{2}/2).$

and when $x$ is small (negative)

$F(x)=1-F(-x)\approx\frac{1}{-x\sqrt{2\pi}}\exp(-x^{2}/2).$

Thus, denoting $f$ the pdf of the standard normal distribution, $\log
p_{0}\approx\log f(x)-log(-x)$ with
$x=\frac{logit(\epsilon)-\mu_{\eta}}{\sigma_{\eta}}$ and $\log
p_{1}\approx\log f(x)-log(x)$ with
$x=\frac{logit(1-\epsilon)-\mu_{\eta}}{\sigma_{\eta}}$ or to summarize both
approximations when $|x|$ is large, $\log f(x)-\log(|x|)$.

Implementation details, if not otherwise specified, are the same as for the
source location example. We plan $K=10$ successive design optimisations and
repeat the whole experiment 100 times for varying values of the true
parameters, for all methods, PASOA, SMC, RL-BOED, VPCE and the random design
baseline. This is overall a more challenging example as the objective function
has many suboptimal local maxima, and the stochastic gradient procedure may be
more sensitive to initialization. The number of gradient steps was set to
5000. For the SMC procedure, the ESS was set to 0.9, the Markov kernel is that
of a random walk Metropolis-Hasting scheme with prior transformations mapping
the parameters to $\mathbb{R}^{4}$. The transformations used are respectively,
$u^{\prime}=\log u$, $\rho^{\prime}=logit(\rho)=\log\frac{\rho}{1-\rho}$, and
$\alpha_{1}^{\prime}=\log\frac{\alpha_{1}}{\alpha_{3}},\alpha_{2}^{\prime}=\log\frac{\alpha_{2}}{\alpha_{3}}$,
with the inverse transformations being $u=\exp u^{\prime}$,
$\rho=\frac{\exp\rho^{\prime}}{1+\exp\rho^{\prime}}$,
$\alpha_{1}=\frac{\exp\alpha_{1}^{\prime}}{1+\exp\alpha_{1}^{\prime}+\exp\alpha_{2}^{\prime}}$,
$\alpha_{2}=\frac{\exp\alpha_{2}^{\prime}}{1+\exp\alpha_{1}^{\prime}+\exp\alpha_{2}^{\prime}}$
and $\alpha_{3}=\frac{1}{1+\exp\alpha_{1}^{\prime}+\exp\alpha_{2}^{\prime}}$.
We use then $L=100$ contrastive variables with each $N=500$ simulations.

Figure 4 in the paper and Figures 4 and 5 below show, with respect to $k$, the
median and standard error of the SPCE, SNMC and Wasserstein distances between
weighted particles and the true parameters. We observe for all methods more
variability in the repetitions for this example. In terms of total EIG, the
difference with RL-BOED is not as large as in the source location example, but
the difference remains large for Wasserstein distances. PASOA still shows
better performance in terms of information gain as measured by SPCE and SNMC.
In Figure 4 in the paper, we observe that in experiments 0-2, our approach
temporarily loses its advantage over RL-BOED due to insufficiently refined
particle approximations of the posteriors. However, this edge is regained in
subsequent experiments as more information from the posteriors becomes
available. Our better design sequences are also visible in the Wasserstein
distance plot presented in main-Figure 4. Figure 4 left shows on the same plot
the SPCE and SNMC curves. Without tempering, SMC gains an advantage only in
the latter steps 7-9 in terms of information gained, while, in the Wasserstein
distance plot presented in main-Figure 4, SMC shows better performance from
the start. A possible explanation is that, as shown in Figure 5 below, RL-BOED
performs better on parameter $\rho$ at the expense of sacrificing precision on
the others. Overall the Wasserstein distance for all parameters remains in
favor of our methods but it may be that a better precision on $\rho$ leads to
a slightly higher information gain (Figure 4 left).

In our current tempering implementation, the Markov kernel is fixed to a
standard Metropolis-Hastings scheme for all steps. It is out of the scope of
this paper but possible directions for improvement include using more
sophisticated kernels, such as Langevin or Hamiltonian Monte Carlo moves, as
suggested in the Tuning parameters section p.1591 of [8] and in references
therein. More generally, a number of recommendations, as reviewed in [8], have
been reported as efficient and could be investigated.

Figure 4: CES example. Median and standard error over 100 tollouts, with
respect to the number of experiments $k$. The [SPCE, SNMC] intervals
containing the totel EIG are plot, respectively with plain (SPCE lower bound)
and dashed (SNMC upper bound) lines. Left: SMC (blue) vs RL-BOED (green) and
VPCE (yellow). Right: SMC (blue) vs PASOA (red).

|  |
---|---|---

Figure 5: CES example. Median and standard error of Wasserstein distances for
each parameter $(\rho,{\boldsymbol{\alpha}},u)$ separately.

###  4.7 Non differentiable examples

When the model log-likelihood is not differentiable, either because the
gradient is not available or difficult to compute, or because the design space
is discrete, the stochastic gradient part of our method cannot be directly
applied. However, we can still use the other parts by replacing the
optimization step by either an exhaustive argmax, in the case of a finite
design space, or by Bayesian optimization [18, 11] which does not requires
gradients. This can be seen as an advantage of myopic solutions, which allow
such replacements to be easily performed. For each sequential optimization,
the search space remains of reasonable size or dimension and does not
increases exponentially. This is not the case for other policy-based approach,
e.g. [9], which would involve a challenging high-dimensional Bayesian
optimization of the policy parameters, or an exhaustive search in an
exponentially increasing space with $K$.

To illustrate this situation, another benchmark example used in [3, 16] is the
Prey population example. The design is a discrete variable. Instead of using
stochastic gradient descent to optimize the $I^{k}_{PCE}$ bound at each step
sequentially, we can compute it for every possible design and take the argmax.
The same can be done to adapt VPCE to this discrete design space, while RL-
BOED has the advantage to be applicable for both continuous and discrete
spaces. Without the gradient part, our approach is similar to that of [16] but
with an additional tempering, which was already reported to compare favorably
to RL-BOED in Figure 4 of [3].

## References

  * [1] L. Aitchison.  Tensor Monte Carlo: particle methods for the GPU era.  ArXiv, 2019.
  * [2] Igor Babuschkin, Kate Baumli, Alison Bell, Surya Bhupatiraju, Jake Bruce, Peter Buchlovsky, David Budden, Trevor Cai, Aidan Clark, Ivo Danihelka, Antoine Dedieu, Claudio Fantacci, Jonathan Godwin, Chris Jones, Ross Hemsley, Tom Hennigan, Matteo Hessel, Shaobo Hou, Steven Kapturowski, Thomas Keck, Iurii Kemaev, Michael King, Markus Kunesch, Lena Martens, Hamza Merzic, Vladimir Mikulik, Tamara Norman, George Papamakarios, John Quan, Roman Ring, Francisco Ruiz, Alvaro Sanchez, Rosalia Schneider, Eren Sezener, Stephen Spencer, Srivatsan Srinivasan, Wojciech Stokowiec, Luyu Wang, Guangyao Zhou, and Fabio Viola.  The DeepMind JAX Ecosystem, 2020.
  * [3] T. Blau, E. V. Bonilla, I. Chades, and A. Dezfouli.  Optimizing sequential experimental design with deep reinforcement learning.  In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 2107–2128. PMLR, 17–23 Jul 2022.
  * [4] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang.  JAX: composable transformations of Python+NumPy programs, 2020\.
  * [5] N. Chopin and O. Papaspiliopoulos.  An introduction to sequential Monte Carlo, volume 4.  Springer, 2020.
  * [6] D. Crisan and A. Doucet.  A survey of convergence results on particle filtering methods for practitioners.  IEEE Transactions on Signal Processing, 50(3):736–746, 2002.
  * [7] M. Cuturi, L. Meng-Papaxanthos, Y. Tian, C. Bunne, G. Davis, and O. Teboul.  Optimal Transport Tools (OTT): A JAX Toolbox for all things Wasserstein.  arXiv preprint arXiv:2201.12324, 2022.
  * [8] C. Dai, J. Heng, P. E. Jacob, and N. Whiteley.  An invitation to sequential Monte Carlo samplers.  Journal of the American Statistical Association, 117(539), 2022\.
  * [9] A. Foster, D. R. Ivanova, I. Malik, and T. Rainforth.  Deep Adaptive Design: Amortizing Sequential Bayesian Experimental Design.  International Conference on Machine Learning, 2021.
  * [10] A. Foster, M. Jankowiak, M. O’Meara, Y. W. Teh, and T. Rainforth.  A Unified Stochastic Gradient Approach to Designing Bayesian-Optimal Experiments.  In Proc. 23rd Int. Conf. Artificial Intelligence and Statistics, 2020.
  * [11] J. M. Henrández-Lobato, M. W. Hoffman, and Z. Ghahramani.  Predictive Entropy Search for Efficient Global Optimization of Black-Box Functions.  In Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 1, NIPS’14, page 918–926, 2014.
  * [12] J. Kuntz, F. Crucinio, and A.M. Johansen.  The divide-and-conquer sequential Monte Carlo algorithm: theoretical properties and limit theorems.  ArXiv, 2023.
  * [13] J. Kuntz, F.R. Crucinio, and A.M. Johansen.  Product-form estimators: exploiting independence to scale up Monte Carlo.  Statistics and Computing, (12), 2022.
  * [14] J. Lao and R. Louf.  Blackjax: A sampling library for JAX, 2020.
  * [15] F. Lindsten, A.M. Johansen, C.A. Naesseth, B. Kirkpatrick, T.B. Schön, J.A.D. Aston, and A. Bouchard-Côté.  Divide-and-Conquer With Sequential Monte Carlo.  Journal of Computational and Graphical Statistics, 26(2):445 – 458, 2017.
  * [16] H. Moffat, M. Hainy, N.E. Papanikolaou, and C. Drovandi.  Sequential experimental design for predator-prey functional response experiments.  Journal of Royal Society Interface, 166(17), 2020.
  * [17] P. Rebeschini and R. van Handel.  Can local particle filters beat the curse of dimensionality?  The Annals of Applied Probability, 25(5), 2015.
  * [18] J. Snoek, H. Larochelle, and R. P. Adams.  Practical Bayesian optimization of machine learning algorithms.  In Advances in neural information processing systems, pages 2951–2959, 2012.
  * [19] A. W. van der Vaart.  Asymptotic Statistics.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 1998.
  * [20] A. W. van der Vaart and J. A. Wellner.  Weak Convergence and Empirical Process: With Applications to Statistics.  springer, 1996.