BAYESIAN EXPERIMENTAL DESIGN VIA CONTRASTIVE DIFFUSIONS

## Abstract

Abstract Bayesian Optimal Experimental Design (BOED) is a powerful tool to reduce the cost of running a sequence of experiments.
When based on the Expected Information Gain (EIG), design optimization corresponds to the maximization of some intractable expected contrast between prior and posterior distributions.
Scaling this maximization to high dimensional and complex settings has been an issue due to BOED inherent computational complexity.
In this work, we introduce an expected posterior distribution with cost-effective sampling properties and provide a tractable access to the EIG contrast maximization via a new EIG gradient expression. Diffusion-based samplers are used to compute the dynamics of the expected posterior and ideas from bi-level optimization are leveraged to derive an efficient joint sampling-optimization loop, without resorting to lower bound approximations of the EIG. The resulting efficiency gain allows to extend BOED to the well-tested generative capabilities of diffusion models.
By incorporating generative models into the BOED framework, we expand its scope and its use in scenarios that were previously impractical. Numerical experiments and comparison with state-of-the-art methods show the potential of the approach.

## 1 Introduction

Designing optimal experiments can be critical in numerous applied contexts where experiments are constrained in terms of resources or more generally costly and limited. In this work, design is assumed to be characterized by some continuous parameters ${\boldsymbol{\xi}}\leavevmode\nobreak\ \in\leavevmode\nobreak\ {\cal E}\subset
\mathbb{R}^{d}$, which
refers to the experimental part, such as the choice of a measurement location, that can be controlled to optimize the experimental outcome.
We consider a Bayesian setting in which the parameters of interest is ${\boldsymbol{\theta}}\in{\boldsymbol{\Theta}}\subset\mathbb{R}^{m}$ and design is optimized to maximize the information gain on ${\boldsymbol{\theta}}$. Bayesian optimal experimental design (BOED) is not a new topic in statistics, see e.g. Chaloner and Verdinelli, (); Sebastiani and Wynn, (); Amzal et al., () but has recently gained new interest with the use of machine learning techniques, see Rainforth et al., (); Huan et al., () for recent reviews.
The most common approach consists of maximizing the so-called expected information
gain (EIG), which is a mutual information criterion that accounts for information via the Shannon’s entropy.
Let $p({\boldsymbol{\theta}})$ denote a prior probability distribution and $p({\mathb{y}}|{\boldsymbol{\theta}},{\boldsymbol{\xi}})$ a likelihood defining the observation ${\mathb{y}}\in{\cal Y}$ generating process. The prior is assumed
to be independent on ${\boldsymbol{\xi}}$ and $p({\mathb{y}}|{\boldsymbol{\theta}},{\boldsymbol{\xi}})$ available in closed-form.
To our knowledge, all previous BOED approaches also assume that the prior is available in closed-form, a setting that we will refer to as density-based BOED. In this work, we will also consider data-based BOED, a setting where only samples from the prior are available.
The EIG, denoted below by $I$, admits several equivalent expressions, see e.g. Foster et al., (). It can be written as the expected loss in entropy when accounting for an observation ${\mathb{y}}$ at ${\boldsymbol{\xi}}$ (eq. (1)) or as a mutual information (MI) or expected Kullback-Leibler (KL) divergence (eq. (2)). Denoting $p_{\boldsymbol{\xi}}({\boldsymbol{\theta}},{\mathb{y}})=p({\boldsymbol{\theta}
},{\mathb{y}}|{\boldsymbol{\xi}})$ the joint distribution of $({\boldsymbol{\theta}},{\mathb{Y}})$ and using $p({\boldsymbol{\theta}},{\mathb{y}}|{\boldsymbol{\xi}})=p({\boldsymbol{\theta}
}|{\mathb{y}},{\boldsymbol{\xi}})p({\mathb{y}}|{\boldsymbol{\xi}})=p({\mathb{y
}}|{\boldsymbol{\theta}},{\boldsymbol{\xi}})p({\boldsymbol{\theta}})$, it comes,

$$ $\displaystyle I({\boldsymbol{\xi}})$ $\displaystyle=\mathbb{E}_{p({\mathb{y}}|{\boldsymbol{\xi}})}[\text{H}(p({ \boldsymbol{\theta}}))-\text{H}(p({\boldsymbol{\theta}}|{\mathb{Y}},{ \boldsymbol{\xi}})]$ (1) $\displaystyle=\mathbb{E}_{p({\mathb{y}}|{\boldsymbol{\xi}})}\left[\text{KL}(p( {\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}}),p({\boldsymbol{\theta}}) )\right]\;=\;\text{MI}(p_{\boldsymbol{\xi}})\;,$ (2) $$

where random variables are indicated with uppercase letters, $\mathbb{E}_{p(\cdot)}[\cdot]$ or $\mathbb{E}_{p}[\cdot]$ denotes the expectation with respect to $p$ and $\text{H}(p({\boldsymbol{\theta}}))=-\mathbb{E}_{p({\boldsymbol{\theta}})}[\log
p
({\boldsymbol{\theta}})]$
is the entropy of $p$. The joint distribution $p_{\boldsymbol{\xi}}$ completely determines all other distributions, marginal (prior) and conditional (posterior) distributions, so that the mutual information, which is the KL between the joint and the product of its marginal distributions, can be written as a function of $p_{\boldsymbol{\xi}}\in{\cal P}({\boldsymbol{\Theta}}\times{\cal Y})$ only. In the following ${\cal P}({\boldsymbol{\Theta}}\times{\cal Y})$, resp. ${\cal P}({\boldsymbol{\Theta}})$, resp. ${\cal P}({\cal Y})$, denotes the set of probability measures on ${\boldsymbol{\Theta}}\times{\cal Y}$, resp. ${\boldsymbol{\Theta}}$, resp. ${\cal Y}$.

In BOED, we look for ${\boldsymbol{\xi}}^{*}$ satisfying

$$ ${\boldsymbol{\xi}}^{*}\in\arg\max_{{\boldsymbol{\xi}}\in\mathbb{R}^{d}}I({ \boldsymbol{\xi}})=\arg\max_{{\boldsymbol{\xi}}\in\mathbb{R}^{d}}\text{MI}(p_{ \boldsymbol{\xi}})\;.$ (3) $$

The above optimization is usually referred to as static design optimization.
The main challenge in EIG-based BOED is that both the EIG and its gradient with respect to ${\boldsymbol{\xi}}$ are doubly intractable. Their respective expressions involve an expectation of an intractable integrand over a posterior distribution which is itself not straightforward to sample from. The posterior distribution is generally only accessible through an iterative algorithm providing approximate samples.
In practice, the inference problem is further complicated as design optimization is considered in a sequential context, in which a series of experiments is planned sequentially and each successive design has to be accounted for.
In order to remove the integrand intractability issue, solutions have been proposed which optimize an EIG lower bound (Foster et al.,, ). This lower bound can be expressed as an expectation of a tractable integrand and becomes tight with increased simulation budgets. The remaining posterior sampling issue has then been solved in different ways. A set of approaches consists of approximating the problematic posterior distribution, either with variational techniques (Foster et al.,, ) or with efficient sequential Monte Carlo (SMC) sampling (Iollo et al.,, ; Drovandi et al.,, ). Other approaches avoid posterior estimation, using reinforcement learning (RL) and off-line policy learning to bypass the need for sampling (Foster et al.,, ; Ivanova et al.,, ; Blau et al.,, ).
However, some studies have shown that estimating the posterior was beneficial, e.g. Iollo et al., () and Ivanova et al., (), which improves on Foster et al., () by introducing posterior estimation steps in order to refine the learned policy. In addition, one should keep in mind that posterior inference is central in BOED as the ultimate goal is not design per se but to gain information on the parameter of interest. This is challenging especially in a sequential context. Previous attempts that provide both candidate design and estimates of the posterior distribution, such as Foster et al., (); Iollo et al., (), are thus essentially in 2 alternating stages, approximate design optimization being dependent on approximate posterior sampling and vice-versa.

In this work, we propose a novel 1-stage approach which leverages a sampling-as-optimization setting (Korba and Salim,, ; Marion et al.,, ) where sampling is seen as an optimization task over the space of probability distributions.
We introduce a new EIG gradient expression (Section 3), which highlights the EIG gradient as a function of both the design and some sampling outcome. This fits into a bi-level optimization framework adapted to BOED in Section 4.
So doing, at each step, both an estimation of the optimal design and samples from the current posterior distribution can be provided in a single loop described in Section 5.
It results an efficient procedure that can handle both traditional density-based samplers and data-based samplers such as provided by the highly successful diffusion-based generative models. The resulting efficiency gain enables BOED applications at significantly larger scales than previously feasible. Notably, this opens the first access to data-based BOED, beyond density-based approaches. Figure 1 is an illustration on a $28\times 28$ image ${\boldsymbol{\theta}}$ reconstruction problem from $7\times 7$ sub-images centered at locations ${\boldsymbol{\xi}}$ to be selected, details in Section 6.
For simpler notation, we first present our approach in the static design case. Adaptation to the sequential case is specified in Section 6 and all numerical examples are in the sequential setting.

Figure: Figure 1: $28\times 28$ Image ${\boldsymbol{\theta}}$ (1st column) reconstruction from seven $7\times 7$ sub-images ${\mathb{y}}={\mathb{A}}_{\boldsymbol{\xi}}{\boldsymbol{\theta}}+{\boldsymbol{ \eta}}$ centered at seven central pixels ${\boldsymbol{\xi}}$ (designs) selected sequentially. Optimized vs. random designs: measured outcome ${\mathb{y}}$ (2nd vs. 3rd column) and parameter ${\boldsymbol{\theta}}$ estimates (reconstruction) with highest weights (upper vs. lower sub-row).
Refer to caption: extracted/5929226/img/mnist_3/comparison.png

## 2 Related work

We focus on gradient-based BOED for continuous problems. Applying a first-order method to solve (3) requires computing gradients of the EIG $I$, which are no more tractable than $I$ itself.
Gradient-based BOED is generally based on stochastic gradient-type algorithms (see Section 4.3.2. in Huan et al., ()). This requires in principle unbiased gradient estimators, although stochastic approximation solutions using biased oracles have also been investigated, see e.g. Demidovich et al., (); Liu and Tajbakhsh, (). To meet this requirement, most stochastic gradient-based approaches start from an EIG lower bound that yields tractable unbiased gradient estimators. More specifically, EIG lower bounds have usually the advantage to remove the nested expectation issue, see e.g. Foster et al., (). In contrast, very few approaches focus on direct EIG gradient estimators. To our knowledge, this is only the case in Goda et al., () and Ao and Li, ().
Goda et al., () propose an unbiased estimator of the EIG gradient using a randomized version of a multilevel
nested Monte Carlo (MLMC) estimator from Rhee and Glynn, ().
A different estimator is proposed by Ao and Li, (), who use MCMC samplers leading to biased estimators, for which the authors show empirically that the bias could be made negligible.
In this work, we first show, in Section 3, that their two apparently different solutions actually only differ in the way the intractable posterior distribution is approximated. We then propose a third way to compute EIG gradients that is more computationally efficient and scales better to larger data volumes and sequential design contexts.
This new expression makes use of a distribution that we introduce and name the expected posterior distribution.
This latter distribution has interesting sampling features that allows us to leverage score-based sampling techniques and connect to the so-called implicit diffusion framework of Marion et al., (). Our single loop procedure in Section 5 is inspired by Marion et al., () and
other recent developments in bi-level optimization (Yang et al.,, ; Dagréou et al.,, ; Hong et al.,, ). However, these latter settings do not cover doubly intractable objectives such as the EIG, which requires both appropriate gradient estimators and sampling operators, see our Sections 3 and 4. In BOED, efficient single loop procedures have been proposed by Foster et al., () but they rely heavily on variational approximations, which may limit accuracy in scenarios
with complex posterior distributions. Additionally, they focus on lower-bound maximization which might sacrifice precision.

## 3 Expected-posterior estimation of the EIG gradient

Efficient EIG gradient estimators are central for accurate scalable BOED.
Gradients derived from the reparameterization trick are often preferred, over the ones obtained with score-based techniques,
as they have been reported to exhibit
lower variance (Xu et al.,, ).

##### EIG gradient via a reparametrization trick.

Assuming $p({\mathb{y}}|{\boldsymbol{\theta}},{\boldsymbol{\xi}})$ is such that ${\mathb{Y}}$ can be rewritten as ${\mathb{Y}}=T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})$ with $T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}$ invertible so that ${\mathb{U}}=T^{-1}_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{Y}})$ and
${\mathb{U}}$ is a random variable independent on ${\boldsymbol{\theta}}$ and ${\boldsymbol{\xi}}$ with a tractable distribution $p_{U}({\mathb{U}})$.
The existence of $T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}$ is straightforward if the direct model corresponds to an additive Gaussian noise as the transformation is then linear in ${\mathb{U}}$. Results exist to guarantee the existence of such a transformation in more general situations (Papamakarios et al.,, ). Using this change of variable,
two expressions of the EIG gradient, (5) and (6) below, can be derived. Detailed steps are given in Appendix A. With $p_{\boldsymbol{\xi}}$ denoting the joint distribution $p({\boldsymbol{\theta}},{\mathb{y}}|{\boldsymbol{\xi}})$,
$g$ a quantity related to the score $g({\boldsymbol{\xi}},{\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{
\prime})=\nabla_{\boldsymbol{\xi}}\log p(T_{{\boldsymbol{\xi}},{\boldsymbol{
\theta}}}({\mathb{u}})|{\boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})_{|{
\mathb{u}}=T^{-1}_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{y}})}\;$ and denoting
$h({\boldsymbol{\xi}},{\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{
\prime})=\nabla_{\boldsymbol{\xi}}p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}
}}({\mathb{u}})|{\boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})_{|{\mathb{u
}}=T^{-1}_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{y}})}$,
a first expression is

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})=$ $\displaystyle\mathbb{E}_{p_{\boldsymbol{\xi}}}\left[g({\boldsymbol{\xi}},{ \mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}})-\frac{\mathbb{E}_{p({ \boldsymbol{\theta}}^{\prime})}{\left[h({\boldsymbol{\xi}},{\mathb{Y}},{ \boldsymbol{\theta}},{\boldsymbol{\theta}}^{\prime})\right]}}{\mathbb{E}_{p({ \boldsymbol{\theta}}^{\prime})}\left[p({\mathb{Y}}|{\boldsymbol{\theta}}^{ \prime},{\boldsymbol{\xi}})\right]}\right]\;.$ (4) $$

Considering importance sampling formulations for the second term
of
(4), with an importance distribution $q\in{\cal P}({\boldsymbol{\Theta}})$, potentially depending on ${\mathb{y}}$, ${\boldsymbol{\theta}}$ and ${\boldsymbol{\xi}}$, further leads to

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})$ $\displaystyle=\mathbb{E}_{p_{\boldsymbol{\xi}}}\left[g({\boldsymbol{\xi}},{ \mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}})-\frac{\mathbb{E}_{q({ \boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{ \xi}})}{\left[\frac{p({\boldsymbol{\theta}}^{\prime})}{q({\boldsymbol{\theta}} ^{\prime}|{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\xi}})}\;h({ \boldsymbol{\xi}},{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{ \prime})\right]}}{\mathbb{E}_{q({\boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{ \boldsymbol{\theta}},{\boldsymbol{\xi}})}\left[\frac{p({\boldsymbol{\theta}}^{ \prime})}{q({\boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{\theta}},{ \boldsymbol{\xi}})}\;p({\mathb{Y}}|{\boldsymbol{\theta}}^{\prime},{\boldsymbol {\xi}})\right]}\right]\;.$ (5) $$

In Goda et al., (), this latter expression is used in a randomized MLMC procedure with $q$ set to a Laplace approximation of the posterior distribution, without justification for this specific choice of $q$.
It results an estimator which is not unbiased
but can be de-biased following Rhee and Glynn, ().
Alternatively, a second expression of the EIG gradient is the starting point of Ao and Li, (),

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})=$ $\displaystyle\mathbb{E}_{p_{\boldsymbol{\xi}}}\left[g({\boldsymbol{\xi}},{ \mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}})-\mathbb{E}_{p({ \boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{\xi}})}\left[g({ \boldsymbol{\xi}},{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{ \prime})\right]\right]\;.$ (6) $$

It follows a nested Monte Carlo estimator (30) given in Appendix A, using samples $\{({\mathb{y}}_{i},{\boldsymbol{\theta}}_{i})\}_{i=1:N}$ from the joint $p_{\boldsymbol{\xi}}$ and for each ${\mathb{y}}_{i}$, samples $\{{\boldsymbol{\theta}}^{\prime}_{i,j}\}_{j=1:M}$ from an MCMC procedure approximating the intractable posterior $p({\boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$.
Interestingly, expression (6) can also be recovered by setting the importance proposal
$q({\boldsymbol{\theta}}^{\prime}|{\mathb{y}},{\boldsymbol{\theta}},{
\boldsymbol{\xi}})$ to $p({\boldsymbol{\theta}}^{\prime}|{\mathb{y}},{\boldsymbol{\xi}})$ in (5), which provides a clear justification of why the choice of $q$ made in Goda et al., () is relevant.
Approaches by Goda et al., () and Ao and Li, () thus mainly differ in their choice of approximations for the posterior distribution. Using a Laplace approximation as in Goda et al., () is relevant only if the posterior is unimodal, which may not be the case in practice. The MCMC version of Ao and Li, () is then potentially more general but also more costly as it requires running $N$ times a MCMC sampler, targeting each time a different posterior $p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$.
In the next paragraph, we introduce the expected posterior distribution and derive another, more computationally efficient, gradient expression.

##### Importance sampling EIG gradient estimator with an expected posterior proposal.

In their work, Ao and Li, () consider only static design, which hides the fact that for more realistic sequential design contexts, their solution is not tractable due to its computational complexity.
Their solution faces the standard issue of nested estimation
(Rainforth et al.,, ).
To avoid this issue we propose to use an importance sampling expression for the second term in (6), which has the advantage to move the dependence on ${\mathb{y}}$ (and ${\boldsymbol{\theta}}$) from the sampling part to the integrand part. We consider a proposal distribution $q\in{\cal P}({\boldsymbol{\theta}})$ that does not depend on ${\mathb{Y}}$ nor ${\boldsymbol{\theta}}$. It comes,

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})$ $\displaystyle=$ $\displaystyle\mathbb{E}_{p_{\boldsymbol{\xi}}}\left[g({\boldsymbol{\xi}},{ \mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}})-\mathbb{E}_{q({ \boldsymbol{\theta}}^{\prime}|{\boldsymbol{\xi}})}\left[\frac{p({\boldsymbol{ \theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{\xi}})}{q({\boldsymbol{\theta}}^{ \prime}|{\boldsymbol{\xi}})}\;g({\boldsymbol{\xi}},{\mathb{Y}},{\boldsymbol{ \theta}},{\boldsymbol{\theta}}^{\prime})\right]\right],$ (7) $$

and an approximate gradient can be obtained as

$$ $\frac{1}{N}\sum_{i=1}^{N}\left[g({\boldsymbol{\xi}},{\mathb{y}}_{i},{ \boldsymbol{\theta}}_{i},{\boldsymbol{\theta}}_{i})-\mathbb{E}_{q({\boldsymbol {\theta}}^{\prime}|{\boldsymbol{\xi}})}\left[\frac{p({\boldsymbol{\theta}}^{ \prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}})}{q({\boldsymbol{\theta}}^{\prime}| {\boldsymbol{\xi}})}\;g({\boldsymbol{\xi}},{\mathb{y}}_{i},{\boldsymbol{\theta }}_{i},{\boldsymbol{\theta}}^{\prime})\right]\right]\;.$ (8) $$

The second term in (8) still requires $N$ importance sampling approximations whose quality depends on the choice of the proposal distribution $q$. The ideal proposal $q$ is easy to simulate, with computable weights at least up to a constant, and so that $q$ and the multiple target distributions $p(\cdot|{\mathb{y}}_{i},{\boldsymbol{\xi}})$ are not too far apart.
Given $N$ samples $\{({\boldsymbol{\theta}}_{i},{\mathb{y}}_{i})\}_{i=1:N}$ from $p_{\boldsymbol{\xi}}$, we propose thus to take $q=q_{{\boldsymbol{\xi}},N}$ where $q_{{\boldsymbol{\xi}},N}$ is the following logarithmic pooling or geometric mixture, with $\sum_{i=1}^{N}\nu_{i}=1$,

$$ $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}})\propto\prod_{i=1}^{N}p({ \boldsymbol{\theta}}|{\mathb{y}}_{i},{\boldsymbol{\xi}})^{\nu_{i}}\propto p({ \boldsymbol{\theta}})\prod_{i=1}^{N}p({\mathb{y}}_{i}|{\boldsymbol{\theta}},{ \boldsymbol{\xi}})^{\nu_{i}}\;.$ (9) $$

We refer to $q_{{\boldsymbol{\xi}},N}$ as the expected posterior distribution, defined in a more general way in (13). It allows to assess the effect of a candidate design ${\boldsymbol{\xi}}$ on samples from the prior. It differs from a standard posterior as no real data ${\mathb{y}}$ obtained by running the experiment ${\boldsymbol{\xi}}$ is available during the optimization. We only have access to samples $\{({\boldsymbol{\theta}}_{i},{\mathb{y}}_{i})\}_{i=1:N}$ from the joint $p_{\boldsymbol{\xi}}$.
The expected posterior can be seen as a distribution that takes into account all possible outcomes of a candidate experiment ${\boldsymbol{\xi}}$ given the samples $\{{\boldsymbol{\theta}}_{i}\}_{i=1:N}$ from the prior.
This choice of $q_{{\boldsymbol{\xi}},N}$ is justified in Appendix B, using Lemma 2, proved therein.
Lemma 2 shows that, for $\sum_{i=1}^{N}\!\nu_{i}=1$, $q_{{\boldsymbol{\xi}},N}$ is the distribution $q$ that minimizes the weighted sum of the KLs against each posterior $p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$, i.e. $\sum_{i=1}^{N}\!\nu_{i}\text{KL}(q,p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{
\boldsymbol{\xi}}))$, leading to an efficient importance sampling proposal.
It follows our new gradient estimator,

$$ $\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})\approx\frac{1}{N}\sum_{i=1}^{N} \left[g({\boldsymbol{\xi}},{\mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{ \boldsymbol{\theta}}_{i})-\frac{1}{M}\sum_{j=1}^{M}w_{i,j}\;g({\boldsymbol{\xi }},{\mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{\boldsymbol{\theta}}^{\prime}_{j })\right]\;,$ (10) $$

where $\{({\boldsymbol{\theta}}_{i},{\mathb{y}}_{i})\}_{i=1:N}$ follow $p_{\boldsymbol{\xi}}$, $\{{\boldsymbol{\theta}}^{\prime}_{i}\}_{j=1:M}$ follow $q_{{\boldsymbol{\xi}},N}$ and $w_{i,j}=\frac{p({\boldsymbol{\theta}}^{\prime}_{j}|{\mathb{y}}_{i},{
\boldsymbol{\xi}})}{q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}}^{\prime}_{j
})}$ denotes the importance sampling weight. When this fraction can only be evaluated up to a constant, we consider self normalized importance sampling (SNIS) using $\tilde{p}$, $\tilde{q}_{{\boldsymbol{\xi}},N}$ the unnormalized versions of $p$ and $q_{{\boldsymbol{\xi}},N}$,

$$ $\tilde{w}_{i,j}=\frac{\tilde{p}({\boldsymbol{\theta}}^{\prime}_{j}|{\mathb{y}} _{i},{\boldsymbol{\xi}})}{\tilde{q}_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta }}^{\prime}_{j})}=\frac{p({\mathb{y}}_{i}|{\boldsymbol{\theta}}^{\prime}_{j},{ \boldsymbol{\xi}})}{\prod_{\ell=1}^{N}p({\mathb{y}}_{\ell}|{\boldsymbol{\theta }}^{\prime}_{j},{\boldsymbol{\xi}})^{\nu_{\ell}}}\qquad\text{and}\qquad w_{i,j }=\frac{\tilde{w}_{i,j}}{\sum_{j=1}^{M}\tilde{w}_{i,j}}\;.$ (11) $$

Although with a reduced computational cost, computing gradients with (10) still requires an iterative sampling algorithm ideally run for a large number of iterations to reach satisfying approximations of the joint $p_{\boldsymbol{\xi}}$ and the expected posterior $q_{{\boldsymbol{\xi}},N}$.
In static design, sampling from the joint is not generally difficult as the prior and the likelihood are assumed available but this becomes problematic in sequential design, as further detailed in Section 6.1.
Sequential design is the setting to be kept in mind in this paper and
in practice, the exact distributions are rarely reached. To assess the impact on gradient approximations,
it is convenient to introduce, as in Marion et al., (),
gradient operators. In the next section, we show how to adapt the formalism of Marion et al., () to our BOED task.

## 4 EIG optimization through sampling

The intractability of the EIG and its gradient links optimization and sampling steps in a nested and costly way.
Marion et al., (), in a different context, propose to adopt a bi-level optimization point of view to integrate sampling as an additional optimization step and derive simpler procedures.

##### Estimation of gradients through sampling.

Denote $\Gamma$ a function from
${\cal P}({\boldsymbol{\Theta}}\times{\cal Y})\times{\cal P}({\boldsymbol{
\Theta}})\times\mathbb{R}^{d}$ to $\mathbb{R}^{d}$, defined as,

$$ $\displaystyle\Gamma(p,q,{\boldsymbol{\xi}})$ $\displaystyle=\mathbb{E}_{p}\left[g({\boldsymbol{\xi}},{\mathb{Y}},{ \boldsymbol{\theta}},{\boldsymbol{\theta}})-\mathbb{E}_{q}\left[\frac{p({ \boldsymbol{\theta}}^{\prime}|{\mathb{Y}})}{q({\boldsymbol{\theta}}^{\prime})} g({\boldsymbol{\xi}},{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{ \prime})\right]\right]\;.$ (12) $$

Expression (7) shows that
$\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})=\Gamma(p_{\boldsymbol{\xi}},q,{
\boldsymbol{\xi}})$,
where $q$ is a distribution $q({\boldsymbol{\theta}}^{\prime}|{\boldsymbol{\xi}})$ on ${\boldsymbol{\theta}}^{\prime}$ possibly depending on ${\boldsymbol{\xi}}$.
The gradient estimator (10) corresponds then to
$\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})\approx\Gamma(\hat{p}_{
\boldsymbol{\xi}},\hat{q}_{{\boldsymbol{\xi}},N},{\boldsymbol{\xi}})$,
where $\hat{p}_{\boldsymbol{\xi}}=\sum_{i=1}^{N}\delta_{({\boldsymbol{\theta}}_{i},{
\mathb{y}}_{i})}$ and
$\hat{q}_{{\boldsymbol{\xi}},N}=\sum_{j=1}^{M}\delta_{{\boldsymbol{\theta}}^{
\prime}_{j}}$.
In general, sampling from $p_{\boldsymbol{\xi}}$, or its sequential counterpart, and $q_{{\boldsymbol{\xi}},N}$ is challenging and only possible through an iterative procedure.
However, an interesting feature of our expected posterior is that it does not add additional sampling difficulties.

##### Expected posterior distribution.

More generally (details in Appendix B), we define,

$$ $\displaystyle q_{{\boldsymbol{\xi}},\rho}({\boldsymbol{\theta}})$ $\displaystyle\propto\exp\left(\mathbb{E}_{\rho}[\log p({\boldsymbol{\theta}}|{ \mathb{Y}},{\boldsymbol{\xi}})]\right)$ (13) $$

where $\rho$ is a measure on ${\cal Y}$. When $\rho({\mathb{y}})\!=\!\sum_{i=1}^{N}\!\nu_{i}\delta_{{\mathb{y}}_{i}}\!({
\mathb{y}})$ with $\sum_{i=1}^{N}\!\nu_{i}\!=\!1$, we recover $q_{{\boldsymbol{\xi}},\rho}({\boldsymbol{\theta}})\!=\!q_{{\boldsymbol{\xi}},N
}({\boldsymbol{\theta}})$ in (9).
The special structure of the expected posterior allows to sample from it using the same algorithmic structure to sample from a single posterior $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$. Indeed, the score of
$q_{{\boldsymbol{\xi}},\rho}({\boldsymbol{\theta}})$ is linked to the posterior score,

$$ $\displaystyle\nabla_{\boldsymbol{\theta}}\log q_{{\boldsymbol{\xi}},\rho}({ \boldsymbol{\theta}})$ $\displaystyle=\mathbb{E}_{\rho}[\nabla_{\boldsymbol{\theta}}\log p({ \boldsymbol{\theta}}|{\mathb{Y}},{\boldsymbol{\xi}})]\;,$ (14) $$

which for $q_{{\boldsymbol{\xi}},N}$ simplifies into $\sum_{i=1}^{N}\nu_{i}\nabla_{\boldsymbol{\theta}}\log p({\boldsymbol{\theta}}|
{\mathb{y}}_{i},{\boldsymbol{\xi}})$.
In practice, we consider the operation of sampling as the output of a
stochastic process iterating a so-called sampling operator.

##### Iterative sampling operators.

Iterative sampling operators, as introduced in Marion et al., (), are mappings from a space, to be specified below, to a space of probabilities.
In our BOED setting, we consider two such operators.
The first one is defined, for each ${\boldsymbol{\xi}}$, through a sequence over $s$ of parameterized functions from ${\cal P}({\boldsymbol{\Theta}}\times{\cal Y})$ to ${\cal P}({\boldsymbol{\Theta}}\times{\cal Y})$ and denoted by
$\Sigma_{s}^{{\mathb{Y}},{\boldsymbol{\theta}}}(p,{\boldsymbol{\xi}})$. Sampling is defined as the outcome in the limit $s\rightarrow\infty$ or for some finite $s=S$ of the following process starting from $p^{(0)}\in{\cal P}({\boldsymbol{\Theta}}\times{\cal Y})$ and iterating

$$ $\displaystyle p^{(s+1)}=\Sigma_{s}^{{\mathb{Y}},{\boldsymbol{\theta}}}(p^{(s)} ,{\boldsymbol{\xi}})\;.$ (15) $$

This formalism includes normalizing flows (Papamakarios et al.,, ) for which $p^{(s)}$ is explicit. More generally, we also consider situations in which $p^{(s)}$ is represented by a random variable ${\mathb{X}}_{s}\sim p^{(s)}$. For example, in density-based BOED, when considering a Langevin diffusion converging to $p_{\boldsymbol{\xi}}$, for a given ${\boldsymbol{\xi}}$,
$p^{(s+1)}=\Sigma_{s}^{{\mathb{Y}},{\boldsymbol{\theta}}}(p^{(s)},{\boldsymbol{
\xi}})$
is a shorthand for its Euler discretization,

$$ $\displaystyle{\mathb{y}}^{(s+1)}$ $\displaystyle={\mathb{y}}^{(s)}-\gamma_{s}\nabla_{\mathb{y}}V({\mathb{y}}^{(s) },{\boldsymbol{\theta}}^{(s)},{\boldsymbol{\xi}})+\sqrt{2\gamma_{s}}{\mathb{B} }_{{\mathb{y}},s}$ $\displaystyle{\boldsymbol{\theta}}^{(s+1)}$ $\displaystyle={\boldsymbol{\theta}}^{(s)}-\gamma_{s}\nabla_{\boldsymbol{\theta }}V({\mathb{y}}^{(s)},{\boldsymbol{\theta}}^{(s)},{\boldsymbol{\xi}})+\sqrt{2 \gamma_{s}}{\mathb{B}}_{{\boldsymbol{\theta}},s}$ $$

where ${\mathb{B}}_{{\mathb{y}},s}$ and ${\mathb{B}}_{{\boldsymbol{\theta}},s}$ are realizations of independent standard Gaussian variables, $\gamma_{s}$ is a step-size and $V({\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\xi}})=-\log p({\boldsymbol{
\theta}})-\log p({\mathb{y}}|{\boldsymbol{\theta}},{\boldsymbol{\xi}})$ is the $p_{\boldsymbol{\xi}}$ potential,
$p_{\boldsymbol{\xi}}({\mathb{y}},{\boldsymbol{\theta}})\propto\exp\left(-V({
\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\xi}})\right)$.
The dynamics induced lead to samples from $p_{\boldsymbol{\xi}}$ for $s\rightarrow\infty$.
In the following, we will thus use the notation $\Sigma_{s}^{{\mathb{Y}},{\boldsymbol{\theta}}}$ to mean that we have access to samples from $p^{(s+1)}$, which is equivalent to apply $\Sigma_{s}^{{\mathb{Y}},{\boldsymbol{\theta}}}$ to an empirical version of $p^{(s)}$ built from samples $\{({\mathb{y}}^{(s)}_{i},{\boldsymbol{\theta}}^{(s)}_{i})\}_{i=1:N}$.
Similarly for a given posterior $p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$, we can produce samples $\{{\boldsymbol{\theta}}^{\prime(s+1)}_{j}\}_{j=1:M}$, with a fixed ${\mathb{y}}_{i}$, using the following updating,

$$ $\displaystyle{\boldsymbol{\theta}}^{\prime(s+1)}$ $\displaystyle={\boldsymbol{\theta}}^{\prime(s)}-\gamma^{\prime}_{s}\nabla_{ \boldsymbol{\theta}}V({\mathb{y}}_{i},{\boldsymbol{\theta}}^{\prime(s)},{ \boldsymbol{\xi}})+\sqrt{2\gamma^{\prime}_{s}}{\mathb{B}}_{{\boldsymbol{\theta }}^{\prime},s}\;,$ (16) $$

while the expected posterior $q_{{\boldsymbol{\xi}},N}$ can be sampled with

$$ $\displaystyle{\boldsymbol{\theta}}^{\prime(s+1)}$ $\displaystyle={\boldsymbol{\theta}}^{\prime(s)}-\gamma^{\prime}_{s}\sum_{i=1}^ {N}\nu_{i}\nabla_{\boldsymbol{\theta}}V({\mathb{y}}_{i},{\boldsymbol{\theta}}^ {\prime(s)},{\boldsymbol{\xi}})+\sqrt{2\gamma^{\prime}_{s}}{\mathb{B}}_{{ \boldsymbol{\theta}}^{\prime},s}\;.$ (17) $$

For a sampling operator of the expected posterior general form (13), we need to extend the definition in Marion et al., () by adding a dependence on some distribution $\rho\in{\cal P}({\cal Y})$ for the conditioning part.
The second sampling operator is defined, for some given ${\boldsymbol{\xi}}$ and $\rho$, through a sequence over $s$ of parameterized functions from ${\cal P}({\boldsymbol{\Theta}})$ to ${\cal P}({\boldsymbol{\Theta}})$ and denoted by
$\Sigma_{s}^{{\boldsymbol{\theta}}^{\prime}}(q,{\boldsymbol{\xi}},\rho)$. The sampling operator is defined as the outcome
of the following process starting from $q^{(0)}\in{\cal P}({\boldsymbol{\Theta}})$ and iterating

$$ $\displaystyle q^{(s+1)}$ $\displaystyle=\Sigma_{s}^{{\boldsymbol{\theta}}^{\prime}}(q^{(s)},{\boldsymbol {\xi}},\rho)\;.$ (18) $$

For instance, when $\rho=\sum_{i=1}^{N}\nu_{i}\delta_{{\mathb{y}}_{i}}$, $q^{(s+1)}=\Sigma_{s}^{{\boldsymbol{\theta}}^{\prime}}(q^{(s)},{\boldsymbol{\xi
}},\rho)$ can then be a shorthand for (17).

Figure: Figure 2: Source localisation example. Prior (left) and expected posterior (right) samples at experiment $k$. Final ${\boldsymbol{\xi}}_{k}^{*}$ (orange cross) at the end of the optimization sequence ${\boldsymbol{\xi}}_{0},\cdot,{\boldsymbol{\xi}}_{T}$ (blue crosses). This optimization "contrasts" the two distributions by making the expected posterior "as different as possible" from the prior.
Refer to caption: extracted/5929226/img/expected/expected_3.png

## 5 Single loop contrastive EIG optimization

The perspective of optimization through sampling
leads naturally to a nested loop procedure.
An inner loop is performed, iterating sampling operators, to reach good approximations of $p_{\boldsymbol{\xi}}$ and $q_{{\boldsymbol{\xi}},N}$ via samples used to evaluate gradient (10).
It results a nested loop algorithm in which the design is updated only after completion of the inner loop, which in our BOED context involves two samplers, as specified in Section 4 and summarized in Algorithm 1.
More efficient procedures consist of iterating alternately between the sampling operator,
e.g. one Langevin step, and the design update, via one iteration of stochastic gradient descent (SGD) or another optimizer.
In particular, Marion et al., () point out the link to bi-level optimization when sampling can be cast as an optimization over the space of distributions.
For instance, the Langevin diffusion dynamics follow a gradient flow of the KL with respect to the Wasserstein-2 distance (Jordan et al.,, ).
This point of view suggests to transform the nested loop Algorithm 1 into a single loop summarized in Algorithm 2. At each optimization step over ${\boldsymbol{\xi}}$, the sampling operators are applied only once using the current ${\boldsymbol{\xi}}$, which is updated in turn, etc.
Algorithms 1 and 2 are in a general form. In this section, we specify the steps that are further illustrated in Section 6.
Regarding the choice of sampling operators,
this approach can handle both traditional density-based sampling where an expression of the target distribution is available and data-based sampling where only training samples are available and a generative model is used to generate additional similar samples.
Examples of both cases are specified below.

Figure: Algorithm 1 ​Nested-loop optimization

Figure: Algorithm 2 Single loop optimization

##### Sampling operators examples for BOED.

Among density-based samplers, we can mention score-based MCMC samplers, including Langevin dynamics via the Unadjusted Langevin Algorithm (ULA) and Metropolis Adjusted Langevin Algorithm (MALA) (Roberts and Tweedie,, ), Hamiltonian Monte Carlo (HMC) samplers (Hoffman and Gelman,, ).
In Section 6.2, an illustration is given with Langevin and sequential Monte Carlo (SMC) to handle a sequential density-based BOED task.

Examples of data-based samplers are generative models such as denoising diffusions and others (Sohl-Dickstein et al.,, ; Ho et al.,, ; Song et al.,, ) that use score-matching techniques (Hyvärinen,, ; Vincent,, ).
Diffusion models, over a time index $t$, learn a score function $\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}}^{(t)})$ using a neural network $s_{\phi}({\boldsymbol{\theta}}^{(t)},t)$ with parameters $\phi$. Given this learned approximation, it is possible to generate new samples of ${\boldsymbol{\theta}}$ by running a backward diffusion process using $s_{\phi}({\boldsymbol{\theta}}^{(t)},t)$ instead of the true score.
In BOED, we are interested in sampling from a conditional distribution $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$ with the following objectives. First we need to sample from the expected posterior $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}})$ which requires conditioning on the observation ${\mathb{y}}$. Second, in sequential design problems (see section 6.1 and Appendix D), we need to condition on the history of observations ${\mathb{D}}_{k-1}$ and produce samples from $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$. Both these issues can be tackled in the framework of denoising diffusions with linear inverse problems (Daras et al.,, ).
When the likelihood corresponds to a measurement ${\mathb{Y}}$ with ${\mathb{Y}}={\mathb{A}}_{\boldsymbol{\xi}}{\boldsymbol{\theta}}+{\boldsymbol{
\eta}}$ and ${\boldsymbol{\eta}}\!\sim\!\mathcal{N}(\mathbf{0},{\boldsymbol{\Sigma}})$, inspiring recent attempts, such as (Corenflos et al.,, ; Cardoso et al.,, ), have addressed the problem of sampling efficiently from $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$ using only the learned score $s_{\phi}({\boldsymbol{\theta}}^{(t)},t)$, without the need for any kind of retraining.
Denoising diffusions include several variants but some formulations consist of running a stochastic differential equation (SDE), see Appendix C. For
conditional sampling, this would mean running an SDE with a conditional score
$\nabla_{\boldsymbol{\theta}}\log p({\boldsymbol{\theta}}^{(t)}|{\mathb{y}},{
\boldsymbol{\xi}})$, which is intractable. See (42) and more details in Appendix C.2.
As a solution, Dou and Song, () propose a method named FPS that approximates $p({\boldsymbol{\theta}}^{(t)}|{\mathb{y}},{\boldsymbol{\xi}})$ by $p({\boldsymbol{\theta}}^{(t)}|{\mathb{y}}^{(t)},{\boldsymbol{\xi}})$ with ${\mathb{y}}^{(t)}$ the noised observation at time $t$. As the score $\nabla_{\boldsymbol{\theta}}\!\log p({\boldsymbol{\theta}}^{(t)}|{\mathb{y}}^{
(t)},{\boldsymbol{\xi}})$ can be written as $\nabla_{\boldsymbol{\theta}}\!\log p({\boldsymbol{\theta}}^{(t)})\!+\!\nabla_{
\boldsymbol{\theta}}\!\log p({\mathb{y}}^{(t)}|{\boldsymbol{\theta}}^{(t)},{
\boldsymbol{\xi}})$, we can leverage the learned score $s_{\phi}({\boldsymbol{\theta}}^{(t)},t)$ and the closed form of $\nabla_{\boldsymbol{\theta}}\log p({\mathb{y}}^{(t)}|{\boldsymbol{\theta}}^{(t
)},{\boldsymbol{\xi}})$ to sample approximately from $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$ using a backward SDE with the approximate score, see (44) in Appendix C.2.
This allows to sample efficiently from $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$
and, using $\nabla_{\boldsymbol{\theta}}\log q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}
}^{\prime})=\sum_{i=1}^{N}\nu_{i}\nabla_{\boldsymbol{\theta}}\log p({
\boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$,
from the expected posterior with the extension of (44) below, where
${\mathb{y}}_{i}^{(t)}$ is the noised ${\mathb{y}}_{i}$ at time $t$,

$$ $d{\boldsymbol{\theta}}^{{}^{\prime}(t)}=\left[-\frac{\beta(t)}{2}{\boldsymbol{ \theta}}^{{}^{\prime}(t)}-\beta(t)\sum_{i=1}^{N}\nu_{i}\nabla_{\boldsymbol{ \theta}}\log p({\boldsymbol{\theta}}^{{}^{\prime}(t)}|{\mathb{y}}_{i}^{(t)},{ \boldsymbol{\xi}})\right]dt+\sqrt{\beta(t)}d{\mathb{B}}_{t}\;.$ (19) $$

An additional resampling SMC-like step can also be added as explained in Appendix E. The approach allows to handle new sequential data-based BOED tasks as illustrated in Section 6.3.

##### Contrastive Optimization.

Optimizing ${\boldsymbol{\xi}}$ using the gradient expression (10) encourages to select a ${\boldsymbol{\xi}}$ that gives either high probability $p({\mathb{y}}_{i}|{\boldsymbol{\theta}}_{i},{\boldsymbol{\xi}})$ to samples $({\boldsymbol{\theta}}_{i},{\mathb{y}}_{i})$ from $p_{\boldsymbol{\xi}}$ or low probability $p({\mathb{y}}_{j}|{\boldsymbol{\theta}}^{\prime}_{j},{\boldsymbol{\xi}})$ to samples ${\boldsymbol{\theta}}^{\prime}_{j}$ from $q_{{\boldsymbol{\xi}},N}$.
This contrastive behaviour is also visible in (2) where the EIG is defined as the mean over the experiment outcomes of the KL between posterior and prior distributions.
The expected posterior $q_{{\boldsymbol{\xi}},N}$ is then used as a proxy to the intractable posterior, to perform this contrastive optimization. Figure 2 provides a visualization of this contrastive behavior in the source localization example of Section 6.2. It corresponds to set the next design ${\boldsymbol{\xi}}$ to a value that eliminates the most parameter ${\boldsymbol{\theta}}$ values (right plot) among the possible ones a priori (left plot).
This is analogous to Noise Constrastive Estimation (Gutmann and Hyvärinen,, ) methods where model parameters are computed so that the data samples are as different as possible from the noise samples.
Additional illustrations are given in Appendix Figure 5.

## 6 Numerical experiments

Two sequential density-based (Section 6.2) and data-based (Section 6.3) BOED examples are considered to illustrate that our method extends to the sequential case in both settings.

### 6.1 Sequential Bayesian experimental design

In the sequential setting, a sequence of $K$ experiments is planned while gradually accounting for the successively collected data. At step $k$, we wish to pick the best design ${\boldsymbol{\xi}}_{k}$ given previous outcomes ${\mathb{D}}_{k-1}=\{({\mathb{y}}_{1},{\boldsymbol{\xi}}_{1}),\dots,({\mathb{y}
}_{k-1},{\boldsymbol{\xi}}_{k-1})\}$. The expected information gain in this scenario is given by:

$$ $I_{k}({\boldsymbol{\xi}},{\mathb{D}}_{k-1})=\mathbb{E}_{p({\mathb{y}}|{ \boldsymbol{\xi}},{\mathb{D}}_{k-1})}\left[\mathrm{KL}(p({\boldsymbol{\theta}} |{\mathb{y}},{\boldsymbol{\xi}},{\mathb{D}}_{k-1}),p({\boldsymbol{\theta}}|{ \mathb{D}}_{k-1}))\right]\;,$ $$

where $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ and $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}},{\mathb{D}}_{k-1})$ act respectively as prior and posterior analogues to the static case (2). See Appendix D for more detailed explanations.
The main difference is that we no longer have direct access to samples from the step $k$ prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$. However, as $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})\propto p({\boldsymbol{\theta}})
\prod_{n=1}^{k-1}p({\mathb{y}}_{n}|{\boldsymbol{\theta}},{\boldsymbol{\xi}}_{n})$ and $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}},{\mathb{D}}_{k-1})
\propto p({\boldsymbol{\theta}})p({\mathb{y}}|{\boldsymbol{\theta}},{
\boldsymbol{\xi}})\prod_{n=1}^{k-1}p({\mathb{y}}_{n}|{\boldsymbol{\theta}},{
\boldsymbol{\xi}}_{n})$ we can still compute the score of these distributions and run sampling operators similar to (15) and (18). To emphasize their dependence on ${\mathb{D}}_{k-1}$, they are denoted by $\Sigma_{s}^{{\mathb{Y}},{\boldsymbol{\theta}}|{\mathb{D}}_{k-1}}(p^{(s)},{
\boldsymbol{\xi}})$ and
$\Sigma_{s}^{{\boldsymbol{\theta}}^{\prime}|{\mathb{D}}_{k-1}}(q^{(s)},{
\boldsymbol{\xi}},\rho)$. Examples of these operators are provided in (20) and (21) in Section 6.2.

Evaluation metrics and comparison.
We refer to our method as CoDiff. In Section 6.2, comparison is provided with other recent approaches, namely a reinforcement learning-based approach RL-BOED from Blau et al., (), the variational prior contrastive estimation VPCE of Foster et al., () and a recent approach named PASOA (Iollo et al.,, ) based on tempered sequential Monte Carlo samplers.
We also compare with a non tempered version of this latter approach (SMC) and with a random baseline, where the observations $\{{\mathb{y}}_{1},\cdot,{\mathb{y}}_{K}\}$ are simulated with designs generated randomly.
More details about these methods are given in Appendix F.2.
To compare methods in terms of information gains, we use
the sequential prior contrastive estimation (SPCE) and sequential nested Monte Carlo (SNMC) bounds introduced in Foster et al., () and used in Blau et al., ().
These quantities allow to compare methods on the produced design sequences only, via their [SPCE, SNMC] intervals which contain the total EIG.
Their expressions are given in Appendix F.1.
We also provide the L2 Wasserstein distance between the produced samples and the true parameter ${\boldsymbol{\theta}}$.
For methods that do not provide posterior estimations or poor quality ones (RL-BOED, VPCE, Random), we compute Wasserstein distances on posterior samples obtained by using tempered SMC on their design and observation sequences. In contrast, SMC Wasserstein distances are computed on the SMC posterior samples.
In Section 6.3, our evaluation is mainly qualitative. The previous methods do not apply and we are not aware of existing attempts that could handle such a generative setting.

### 6.2 Sources location finding

We present a source localization example inspired by Foster et al., (); Blau et al., ().
The setup involves $C$ sources in $\mathbb{R}^{2}$, with unknown positions ${\boldsymbol{\theta}}=\{{\boldsymbol{\theta}}_{1},\dots,{\boldsymbol{\theta}}_
{C}\}$.
The challenge is to determine optimal measurement locations to accurately infer the sources positions. When a measurement is taken at location ${\boldsymbol{\xi}}\in\mathbb{R}^{2}$, the signal strength is defined as $\mu({\boldsymbol{\theta}},{\boldsymbol{\xi}})=b+\sum_{c=1}^{C}\frac{\alpha_{c}
}{m+\|{\boldsymbol{\theta}}_{c}-{\boldsymbol{\xi}}\|_{2}^{2}}$ where $\alpha_{c}$, $b$, and $m$ are predefined constants.
We assume a standard Gaussian prior for each source location, ${\boldsymbol{\theta}}_{c}\sim\mathcal{N}(0,{\mathb{I}}_{2})$, and model the likelihood as log-normal:
$(\log{\mathb{y}}\mid{\boldsymbol{\theta}},{\boldsymbol{\xi}})\sim\mathcal{N}(
\log\mu({\boldsymbol{\theta}},{\boldsymbol{\xi}}),\sigma)$,
with $\sigma$ representing the standard deviation. For this experiment, we set $C=2$, $\alpha_{1}=\alpha_{2}=1$, $m=10^{-4}$, $b=10^{-1}$, $\sigma=0.5$, and plan $K=30$ sequential design optimizations.
In the notation of the single loop Algorithm 2, we consider $\Sigma_{t}^{{\mathb{Y}},{\boldsymbol{\theta}}|{\mathb{D}}_{k-1}}(p^{(t)},{
\boldsymbol{\xi}}_{t})$ and
$\Sigma_{t}^{{\boldsymbol{\theta}}^{\prime}|{\mathb{D}}_{k-1}}(q^{(t)},{
\boldsymbol{\xi}}_{t},\hat{\rho}_{t+1})$ operators
that correspond respectively to the update of batch samples
of size $N=200$ and $M=200$ $\{({\mathb{y}}_{i}^{(t)},{\boldsymbol{\theta}}_{i}^{(t)})\}_{i=1:N}$ and $\{{\boldsymbol{\theta}}_{j}^{{}^{\prime}(t)}\}_{j=1:M}$ with $\hat{\rho}_{t+1}=\sum_{i=1}^{N}\nu_{i}\delta_{{\mathb{y}}_{i}^{(t+1)}}$ using Langevin diffusions. Making use of the availability of the likelihood in this example, sampling from it, is straightforward and sampling operator iterations simplify into, for $i=1:N$ and $j=1:M$

$$ $\displaystyle{\boldsymbol{\theta}}^{(t+1)}_{i}$ $\displaystyle={\boldsymbol{\theta}}^{(t)}_{i}-\gamma_{t}\nabla_{\boldsymbol{ \theta}}\log p({\boldsymbol{\theta}}^{(t)}_{i}|{\mathb{D}}_{k-1})+\sqrt{2 \gamma_{t}}{\mathb{B}}_{{\boldsymbol{\theta}},t}\quad\text{and}\quad{\mathb{y} }^{(t+1)}_{i}\sim p({\mathb{y}}|{{\boldsymbol{\theta}}}^{(t+1)}_{i},{ \boldsymbol{\xi}}_{t})$ (20) $\displaystyle{{\boldsymbol{\theta}}^{\prime}}^{(t+1)}_{j}$ $\displaystyle={{\boldsymbol{\theta}}^{\prime}}^{(t)}_{j}-\gamma^{\prime}_{t} \sum_{i=1}^{N}\nu_{i}\nabla_{\boldsymbol{\theta}}\log p({{\boldsymbol{\theta}} ^{\prime}}^{(t)}_{j}|{\mathb{y}}^{(t+1)}_{i},{\boldsymbol{\xi}}_{t},{\mathb{D} }_{k-1})+\sqrt{2\gamma^{\prime}_{t}}{\mathb{B}}_{{\boldsymbol{\theta}}^{\prime },t}\;.$ (21) $$

In practice, Langevin diffusion can get trapped in local minima and causes the sampling to be too slow to keep pace with the optimization process. To address this, we augment the Langevin diffusion with the Diffusive Gibbs (DiGS) MCMC kernel proposed by Chen et al., ().
DiGS is an auxiliary variable MCMC method where the auxiliary variable $\tilde{{\mathb{X}}}$ is a noisy version of the original variable ${\mathb{X}}$. DiGS enhances mixing and helps escape local modes by alternately sampling from the distributions $p(\tilde{x}|x)$, which introduces noise via Gaussian convolution, and $p(x|\tilde{x})$, which denoises the sample back to the original space using a score-based update (here a Langevin diffusion).
With 400 total samples, each measurement step takes 2.9 s. This number of samples is insightful as it is usually the amount of samples one can afford to compute in the diffusion models of Section 6.3.
The whole experiment is repeated 100 times with random source locations each time.
Figure 3 shows, with respect to $k$, the median for SPCE, the L2 Wasserstein distances between weighted samples and the true source locations and SNMC.
CoDiff clearly outperforms all other methods, with
significant improvement, both in terms of information gain and posterior estimation. It improves by $30\%$ the non-myopic RL-BOED results on SPCE and provides much higher SNMC. The L2 Wasserstein distance is two order of magnitude lower, suggesting the higher quality of our measurements.

Figure: Figure 3: Source location. Median over 100 rollouts for SPCE, L2 Wasserstein distance (logarithmic scale) and SNMC with respect to the number of experiments $k$. Number of samples $N+M=400$
Refer to caption: extracted/5929226/img/sources/spce.png

### 6.3 Image reconstruction with Diffusion models

We build an artificial experimental design task to illustrate the ability of our method to handle design parameters related to inverse problems where the parameter ${\boldsymbol{\theta}}$ of interest is high dimensional. More specifically, we consider the task that consists of recovering an hidden image
from only partial observations of its pixels.
The image to be recovered is denoted by ${\boldsymbol{\theta}}$. An experiment corresponds to the choice of a pixel ${\boldsymbol{\xi}}$ around which an observation mask is centered and the image becomes visible. The measured observation ${\mathb{y}}$ is then a masked version of ${\boldsymbol{\theta}}$.
The likelihood derives from the model ${\mathb{Y}}={\mathb{A}}_{\boldsymbol{\xi}}{\boldsymbol{\theta}}+{\boldsymbol{
\eta}}$ where ${\mathb{A}}_{\boldsymbol{\xi}}$ is a square mask centered at ${\boldsymbol{\xi}}$ and ${\boldsymbol{\eta}}$ some Gaussian variable. For the image prior, we consider a diffusion model trained for generation of the MNIST dataset (LeCun et al.,, ). The goal is thus to select sequentially the best central pixel locations for $7\times 7$ masks so as to reconstruct an entire $28\times 28$ MNIST image in the smallest number of experiments. The smaller the mask the more interesting it becomes to optimally select the mask centers.
Algorithm 2 is used with diffusion-based sampling operators specified in Appendix F.2.2. The gain in optimizing the mask placements is illustrated in Figure 1 and Appendix Figure 8. Progressive image reconstructions are shown in Figure 4 and Appendix Figures 6 and 7.
The digit to be recovered is shown in the 1st column. The successively selected masks are shown (red line squares) in the 2nd column with the resulting gradually discovered part of the image.
The reconstruction per se can be estimated from the posterior samples shown in the last 16 columns. At each experiment, the upper sub-row shows the 16 most-likely reconstructed images, while the lower sub-row shows the 16 less probable ones.
As the number of experiments increases the posterior samples gradually concentrate on the right digit.

Figure: Figure 4: Image reconstruction. First 6 experiments (rows): image ground truth, measurement at experiment $k$, samples from current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$, with best (resp. worst) weights in upper (resp. lower) sub-row. The samples incorporate past measurement information as the procedure advances. Each design steps takes $\sim 7.3$s
Refer to caption: extracted/5929226/img/mnist_3/samples_0.png

## 7 Conclusion

We presented a new approach to gradient-based BOED that allows very efficient implementations.
Introducing the concept of expected posterior distribution, we proposed a new EIG gradient estimator,
which allows the EIG maximization to be performed in a bi-level optimization framework and without resorting to lower bounds.
The performance was illustrated in a traditional density-based setting with superior accuracy and lower computational cost compared to state-of-the-art methods.
In addition, the possibility of our method to also handle data-based sampling represents, to our knowledge, the first extension of BOED to diffusion-based generative models.
By integrating the highly successful framework of diffusion models for our sampling operators, we were able to optimize a design parameter ${\boldsymbol{\xi}}$ concurrently with the diffusion process.
This was illustrated in a new application for BOED involving high dimensional image parameters.
The foundation of our approach lies on conditional diffusion models and their application to inverse problems.
Thanks to this advancement, there are as many new potential applications of BOED as there are trained diffusion models for specific inverse problem tasks. Currently, most models assume a linear forward model.
However, the non-linear setting is an active field of research, and advancements in this area could be directly applied to our framework.
The applicability of our method could also be extended by considering settings with no explicit expression of the likelihood and investigating simulation-based inference such as developed by Ivanova et al., (); Kleinegesse and Gutmann, (); Kleinegesse et al., ().
In addition, although in density-based BOED, we have shown that greedy approaches could outperform long-sighted reinforcement learning procedures, in a data-based setting, it would be interesting to investigate an extension to non myopic approaches such as Iqbal et al., ().

## Appendix A Two expressions for the EIG gradient

Both approaches presented below, that of Goda et al., () and Ao and Li, (), start from EIG gradient expressions derived using a reparameterization trick.
Using the change of variable ${\mathb{Y}}=T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})$, we can derive the following expression for the EIG gradient,

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})=$ $\displaystyle\mathbb{E}_{p_{U}({\mathb{u}})p({\boldsymbol{\theta}})}\left[ \nabla_{\boldsymbol{\xi}}\log p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({ \mathb{U}})|{\boldsymbol{\theta}};{\boldsymbol{\xi}})\right]-\mathbb{E}_{p_{U} ({\mathb{u}})p({\boldsymbol{\theta}})}\left[\nabla_{\boldsymbol{\xi}}\log p(T_ {{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\xi}}) \right]\;.$ (22) $$

The first term in (22) involves only the known likelihood and is generally not problematic. For the second term, we can use,

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}\log p(T_{{\boldsymbol{\xi}},{ \boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\xi}})$ $\displaystyle=\frac{\nabla_{\boldsymbol{\xi}}p(T_{{\boldsymbol{\xi}},{ \boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\xi}})}{p(T_{{\boldsymbol{\xi} },{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\xi}})}$ (23) $$

with $p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\xi}}
)=\mathbb{E}_{p({\boldsymbol{\theta}}^{\prime})}\left[p(T_{{\boldsymbol{\xi}},
{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\theta}}^{\prime},{
\boldsymbol{\xi}})\right]$ and

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}p(T_{{\boldsymbol{\xi}},{\boldsymbol{ \theta}}}({\mathb{U}})|{\boldsymbol{\xi}})$ $\displaystyle=\nabla_{\boldsymbol{\xi}}\mathbb{E}_{p({\boldsymbol{\theta}}^{ \prime})}\left[p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{ \boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})\right]$ $\displaystyle=\mathbb{E}_{p({\boldsymbol{\theta}}^{\prime})}\left[\nabla_{ \boldsymbol{\xi}}\;p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}}) |{\boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})\right]$ (24) $\displaystyle=\mathbb{E}_{p({\boldsymbol{\theta}}^{\prime})}\left[p(T_{{ \boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\theta}}^{ \prime},{\boldsymbol{\xi}})\;\nabla_{\boldsymbol{\xi}}\log p(T_{{\boldsymbol{ \xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\theta}}^{\prime},{ \boldsymbol{\xi}})\right]\;.$ (25) $$

Subsequently, two expressions of the EIG gradient can be derived depending on which of (24) or (25) is used.
Using (24) and
$p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\xi}}
)=\mathbb{E}_{p({\boldsymbol{\theta}}^{\prime})}\left[p(T_{{\boldsymbol{\xi}},
{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\theta}}^{\prime},{
\boldsymbol{\xi}})\right],$
it comes

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})=$ $\displaystyle\mathbb{E}_{p_{U}({\mathb{u}})p({\boldsymbol{\theta}})}\left[ \nabla_{\boldsymbol{\xi}}\log p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({ \mathb{U}})|{\boldsymbol{\theta}};{\boldsymbol{\xi}})-\frac{\mathbb{E}_{p({ \boldsymbol{\theta}}^{\prime})}{\left[\nabla_{\boldsymbol{\xi}}\;p(T_{{ \boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\theta}}^{ \prime},{\boldsymbol{\xi}})\right]}}{\mathbb{E}_{p({\boldsymbol{\theta}}^{ \prime})}\left[p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{ \boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})\right]}\right]$ $\displaystyle=\mathbb{E}_{p_{\boldsymbol{\xi}}}\left[g({\boldsymbol{\xi}},{ \mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}})-\frac{\mathbb{E}_{p({ \boldsymbol{\theta}}^{\prime})}{\left[h({\boldsymbol{\xi}},{\mathb{Y}},{ \boldsymbol{\theta}},{\boldsymbol{\theta}}^{\prime})\right]}}{\mathbb{E}_{p({ \boldsymbol{\theta}}^{\prime})}\left[p({\mathb{Y}}|{\boldsymbol{\theta}}^{ \prime},{\boldsymbol{\xi}})\right]}\right]\;.$ (26) $$

with

$$ $g({\boldsymbol{\xi}},{\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{ \prime})\!=\!\nabla_{\boldsymbol{\xi}}\!\log p(T_{{\boldsymbol{\xi}},{ \boldsymbol{\theta}}}({\mathb{u}})|{\boldsymbol{\theta}}^{\prime},\!{ \boldsymbol{\xi}})_{|{\mathb{u}}=T^{-1}_{{\boldsymbol{\xi}},{\boldsymbol{ \theta}}}({\mathb{y}})}$ $$

and

$$ $h({\boldsymbol{\xi}},{\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{ \prime})\!=\!\nabla_{\boldsymbol{\xi}}p(T_{{\boldsymbol{\xi}},{\boldsymbol{ \theta}}}({\mathb{u}})|{\boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})_{|{ \mathb{u}}=T^{-1}_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{y}})}.$ $$

Considering, in the second term, an additional importance distribution $q({\boldsymbol{\theta}}^{\prime}|{\mathb{y}},{\boldsymbol{\theta}},{
\boldsymbol{\xi}})$ leads to
the expression used in Goda et al., (),

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})$ $\displaystyle=\mathbb{E}_{p_{\boldsymbol{\xi}}}\left[g({\boldsymbol{\xi}},{ \mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}})-\frac{\mathbb{E}_{q({ \boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{ \xi}})}{\left[\frac{p({\boldsymbol{\theta}}^{\prime})}{q({\boldsymbol{\theta}} ^{\prime}|{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\xi}})}h({\boldsymbol {\xi}},{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{\prime})\right ]}}{\mathbb{E}_{q({\boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{ \theta}},{\boldsymbol{\xi}})}\left[\frac{p({\boldsymbol{\theta}}^{\prime})}{q( {\boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{ \xi}})}p({\mathb{Y}}|{\boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})\right] }\right]\;.$ (27) $$

It can be used to
derive estimators of the form,

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})$ $\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\left[g({\boldsymbol{\xi}},{ \mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{\boldsymbol{\theta}}_{i})-\frac{ \frac{1}{M}\sum_{j=1}^{M}\frac{p({\boldsymbol{\theta}}^{\prime}_{i,j})}{q({ \boldsymbol{\theta}}^{\prime}_{i,j}|{\mathb{y}}_{i},{\boldsymbol{\theta}}_{i}, {\boldsymbol{\xi}})}h({\boldsymbol{\xi}},{\mathb{y}}_{i},{\boldsymbol{\theta}} _{i},{\boldsymbol{\theta}}^{\prime}_{i,j})}{\frac{1}{M}\sum_{j=1}^{M}\frac{p({ \boldsymbol{\theta}}^{\prime}_{i,j})}{q({\boldsymbol{\theta}}^{\prime}_{i,j}|{ \mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{\boldsymbol{\xi}})}p({\mathb{y}}_{i} |{\boldsymbol{\theta}}^{\prime}_{i,j},{\boldsymbol{\xi}})}\right]\;,$ (28) $$

where $\{({\mathb{y}}_{i},{\boldsymbol{\theta}}_{i})\}_{i=1:N}$ are simulated from the joint distribution $p_{\boldsymbol{\xi}}$ and for each $i=1:N$,
$\{{\boldsymbol{\theta}}^{\prime}_{i,j}\}_{j=1:M}$ is a sample from $q(\cdot|{\mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{\boldsymbol{\xi}})$.
Goda et al., () use (28) with $N=1$.
Even with perfect sampling, this estimator is not unbiased due to the ratio in the second term but can be de-biased following Rhee and Glynn, ().
The randomized MLMC procedure of Rhee and Glynn, () is a post-hoc general procedure that can be more generally applied to de-bias a sequence of possibly biased estimators, provided the estimators are consistent.

Alternatively, using (25) instead, another expression of the EIG gradient can be derived. Replacing (25) in (23), it comes,

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}\log p(T_{{\boldsymbol{\xi}},{ \boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\xi}})$ $\displaystyle=\mathbb{E}_{p({\boldsymbol{\theta}}^{\prime})}\left[\frac{p(T_{{ \boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\theta}}^{ \prime},{\boldsymbol{\xi}})\;\nabla_{\boldsymbol{\xi}}\log p(T_{{\boldsymbol{ \xi}},{\boldsymbol{\theta}}}({\mathb{U}})|{\boldsymbol{\theta}}^{\prime},{ \boldsymbol{\xi}})}{p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U}} )|{\boldsymbol{\xi}})}\right]$ $\displaystyle=\mathbb{E}_{p({\boldsymbol{\theta}}^{\prime}|T_{{\boldsymbol{\xi }},{\boldsymbol{\theta}}}({\mathb{U}}),{\boldsymbol{\xi}})}\left[\nabla_{ \boldsymbol{\xi}}\log p(T_{{\boldsymbol{\xi}},{\boldsymbol{\theta}}}({\mathb{U }})|{\boldsymbol{\theta}}^{\prime},{\boldsymbol{\xi}})\right]\;,$ $$

which, with the definition of $g$ above, leads to

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})=$ $\displaystyle\mathbb{E}_{p_{\boldsymbol{\xi}}}\left[g({\boldsymbol{\xi}},{ \mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}})-\mathbb{E}_{p({ \boldsymbol{\theta}}^{\prime}|{\mathb{Y}},{\boldsymbol{\xi}})}\left[g({ \boldsymbol{\xi}},{\mathb{Y}},{\boldsymbol{\theta}},{\boldsymbol{\theta}}^{ \prime})\right]\right]\;.$ (29) $$

This alternative expression (29) is the starting point of Ao and Li, (), who subsequently use the following estimator,

$$ $\displaystyle\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})\approx$ $\displaystyle\frac{1}{N}\sum_{i=1}^{N}\left[g({\boldsymbol{\xi}},{\mathb{y}}_{ i},{\boldsymbol{\theta}}_{i},{\boldsymbol{\theta}}_{i})-\mathbb{E}_{q({ \boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}})}\left[g({ \boldsymbol{\xi}},{\mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{\boldsymbol{ \theta}}^{\prime})\right]\right]\;,$ $$

where $\{({\mathb{y}}_{i},{\boldsymbol{\theta}}_{i})\}_{i=1:N}$ is as before a sample from the joint distribution $p_{\boldsymbol{\xi}}$ and where for each ${\mathb{y}}_{i}$, $q({\boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$ is a tractable approximation of the intractable posterior $p({\boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$. More specifically, Ao and Li, () propose to approximate each posterior distribution by $q({\boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}})=\frac{1}{
M}\sum_{j=1}^{M}\delta_{{\boldsymbol{\theta}}^{\prime}_{i,j}}$, using a
sample $\{{\boldsymbol{\theta}}^{\prime}_{i,j}\}_{j=1:M}$ from an MCMC procedure. It follows the nested Monte Carlo estimator below,

$$ $\nabla_{\boldsymbol{\xi}}I({\boldsymbol{\xi}})\approx\frac{1}{N}\sum_{i=1}^{N} \left[g({\boldsymbol{\xi}},{\mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{ \boldsymbol{\theta}}_{i})-\frac{1}{M}\sum_{j=1}^{M}g({\boldsymbol{\xi}},{ \mathb{y}}_{i},{\boldsymbol{\theta}}_{i},{\boldsymbol{\theta}}_{i,j})\right]\;.$ (30) $$

## Appendix B Logarithmic pooling as a good importance sampling proposal

When considering importance sampling with a proposal distribution $q$ and a target distribution $p$,
Chatterjee and Diaconis, () proved that under certain conditions, the number of simulation draws required for both importance sampling and self normalized importance sampling (SNIS) estimators to have small L1 error with high probability was roughly $\exp(\text{KL}(p,q))$, see Theorem 1.2 in Chatterjee and Diaconis, () for SNIS.
Similarly, selecting a proposal distribution which minimizes the importance sampling estimator variance is equivalent to finding a distribution with small $\chi^{2}$-distance to $p$, see e.g. Appendix E of Minka, ().
More generally, finding a good proposal $q$ is linked to the problem of minimizing $\alpha$-divergences or $f$-divergence between $p$ and $q$, which are jointly convex in $p$ and $q$, see Minka, ().
In this work, we consider $\text{KL}(q,p)$ as a measure of proximity between $p$ and $q$. This choice is ultimately arbitrary but has the advantage of leading to an interpretable proposal with interesting sampling properties.
To justify the expected posterior $q_{{\boldsymbol{\xi}},N}$ in (9) and its use in (10), we then use Lemma 2 below to show that
for $\sum_{i=1}^{N}\nu_{i}=1$, the
distribution $q^{*}$ that minimizes the weighted sum of the KL against each posterior $p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$, i.e. $\sum_{i=1}^{N}\nu_{i}\text{KL}(q,p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{
\boldsymbol{\xi}}))$ is

$$ $\displaystyle q^{*}({\boldsymbol{\theta}})$ $\displaystyle\propto$ $\displaystyle p({\boldsymbol{\theta}})\prod_{i=1}^{N}p({\mathb{y}}_{i}|{ \boldsymbol{\theta}},{\boldsymbol{\xi}})^{\nu_{i}}$ (31) $\displaystyle\propto$ $\displaystyle\prod_{i=1}^{N}p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{ \boldsymbol{\xi}})^{\nu_{i}}\;,$ (32) $$

which is the logarithmic pooling (or geometric mixture) of the respective posterior distributions $p({\boldsymbol{\theta}}|{\mathb{y}}_{i},{\boldsymbol{\xi}})$.
Lemma 2 results from an application of a lemma mentioned by Alquier, () (Lemma 2.2 therein), and recalled below in Lemma 1.
This Lemma 1 has been known since Kullback (Kullback,, ) in the case of a finite parameter space ${\boldsymbol{\Theta}}$,
but the general case is due to Donsker and Varadhan (Donsker and Varadhan,, ).
Recall that ${\cal P}({\boldsymbol{\Theta}})$ denotes the set of probability measures on ${\boldsymbol{\Theta}}$ and $p$ a given probability measure in ${\cal P}({\boldsymbol{\Theta}})$.

###### Lemma 1 (Donsker and Varadhan’s variational formula)

For any measurable, bounded
function $f:{\boldsymbol{\Theta}}\rightarrow\mathbb{R}$, the supremum with respect to $q\in{\cal P}({\boldsymbol{\Theta}})$ of

$$ $\mathbb{E}_{q}\left[f({\boldsymbol{\theta}})\right]-\text{KL}(q,p)$ $$

is the following Gibbs measure $p_{f}$ defined by its density with respect to $p$,

$$ $dp_{f}=\frac{\exp(f({\boldsymbol{\theta}}))}{\mathbb{E}_{p}\left[\exp(f({ \boldsymbol{\theta}}))\right]}dp\;.$ $$

The following Lemma 2 is an application of Lemma 1.

###### Lemma 2

For a given probability measure $p\in{\cal P}({\boldsymbol{\Theta}})$ and a measure $\rho$ on ${\cal Y}$ (not necessarily a probability measure), define for any probability measure $q\in{\cal P}({\boldsymbol{\Theta}})$

$$ $\displaystyle\ell(q)$ $\displaystyle=$ $\displaystyle\mathbb{E}_{q}\left[\mathbb{E}_{{\mathb{Y}}\sim\rho}\left[\log p( {\mathb{Y}}|{\boldsymbol{\theta}})\right]\right]-\text{KL}(q,p)\;.$ (33) $$

It results from the Donsker and Varadhan’s variational formula Lemma 1 that the supremum of $\ell(q)$ with respect to $q$ is reached for the Gibbs
measure $q^{*}$ defined by its density with respect to $p$,

$$ $q^{*}({\boldsymbol{\theta}})\propto p({\boldsymbol{\theta}})\;\exp(\mathbb{E}_ {Y\sim\rho}[\log p({\mathb{Y}}|{\boldsymbol{\theta}})])\;.$ $$

In addition maximizing $\ell$ is equivalent to minimising

$$ $\mathbb{E}_{{\mathb{Y}}\sim\rho}\left[\text{KL}(q({\boldsymbol{\theta}}),p({ \boldsymbol{\theta}}|{\mathb{Y}}))\right]$ $$

which means that $q^{*}$ is the measure that minimizes the KL to each $p({\boldsymbol{\theta}}|{\mathb{y}})$ on average with respect to ${\mathb{y}}$.

###### Lemma 1 (Donsker and Varadhan’s variational formula)

###### Lemma 2

##### Proof of Lemma 2

The expression of $q^{*}$ results from a direct application of Lemma 1 to $f({\boldsymbol{\theta}})=\mathbb{E}_{{\mathb{Y}}\sim\rho}\left[\log p({\mathb{
Y}}|{\boldsymbol{\theta}})\right]$ assuming it is measurable and bounded as a function of ${\boldsymbol{\theta}}$ (to be checked in practice).
The second part results from rewriting $\ell$ as

$$ $\displaystyle\ell(q)$ $\displaystyle=$ $\displaystyle\mathbb{E}_{{\mathb{Y}}\sim\rho}\left[\mathbb{E}_{{\boldsymbol{ \theta}}\sim q}\left[\frac{\log(p({\mathb{Y}}|{\boldsymbol{\theta}})p({ \boldsymbol{\theta}}))}{\log q({\boldsymbol{\theta}})}\right]\right]$ $\displaystyle=$ $\displaystyle-\mathbb{E}_{{\mathb{Y}}\sim\rho}\left[\text{KL}(q,p({\boldsymbol {\theta}}|{\mathb{Y}})\right]+\mathbb{E}_{{\mathb{Y}}\sim\rho}\left[\log p({ \mathb{Y}})\right]\;.$ $$

Example: as already mentioned our expected posterior $q_{{\boldsymbol{\xi}},N}$ corresponds to the application of this result to $\rho=\sum_{i=1}^{N}\nu_{i}\delta_{{\mathb{y}}_{i}}$ with $\sum_{i=1}^{N}\nu_{i}=1$.

Remark 1: If $\rho=\delta_{\mathb{y}}$, or $\rho=\sum_{i=1}^{N}\delta_{{\mathb{y}}_{i}}$, we recover the standard variational formulation of the posterior distribution (see e.g. Table 1 in (Knoblauch et al.,, )). The posterior distribution $p({\boldsymbol{\theta}}|{\mathb{y}}_{1},\ldots,{\mathb{y}}_{N})$ differs from the logarithmic pooling (for which the weights $\nu_{i}$ sum to 1) in the relative weight given to the prior. The result is valid for very general $\ell$ not necessarily expressed as an expectation.

Remark 2: Regarding logarithmic pooling, the result is similar to a result in Carvalho et al., () (Remark 3.1 therein) by showing that, in the case of the sum of the KL, $\sum_{i=1}^{N}\text{KL}(q,p({\boldsymbol{\theta}}|{\mathb{y}}_{i}))$, the optimal pooling weights are equal, $\nu_{i}=\frac{1}{N}$.

Remark 3: The expected posterior distribution can also be recovered as a constrained mean field solution. Indeed, it is easy to show that $q^{*}$ is also the measure that minimizes the KL between the joint distribution and a product form approximation where one of the factor is fixed to $\rho({\mathb{y}})$,

$$ $q^{*}=\arg\min_{q\in{\cal P}({\boldsymbol{\Theta}})}\text{KL}(q({\boldsymbol{ \theta}})\rho({\mathb{y}}),p({\boldsymbol{\theta}},{\mathb{Y}}))\;.$ $$

## Appendix C Diffusion-based generative models

### C.1 Denoising diffusion models

Given a distribution $p_{0}\in{\cal P}({\boldsymbol{\Theta}})$ only available through a set of samples of ${\boldsymbol{\theta}}$, diffusion models are based on the addition of noise to the available samples in such a manner that allows to learn the reverse process that "denoises" the samples. This learned process can then be exploited to generate new samples by denoising random noise samples until we get back to the original data distribution.
As an appropriate noising process, in our experiments we ran the Variance Preserving SDE from Dhariwal and Nichol, ():

$$ $d\tilde{{\boldsymbol{\theta}}}^{(t)}=-\frac{\beta(t)}{2}\tilde{{\boldsymbol{ \theta}}}^{(t)}dt+\sqrt{\beta(t)}d\tilde{{\mathb{B}}}_{t}$ (34) $$

where $\beta(t)>0$ is a linear noise schedule that controls the amount of noise added at time $t$.
Solving SDE (34) leads to

$$ $(\tilde{{\boldsymbol{\theta}}}^{(t)}|\tilde{{\boldsymbol{\theta}}}^{(0)})\sim \mathcal{N}\!\left(\tilde{{\boldsymbol{\theta}}}^{(0)}\exp(-\frac{1}{2}\int_{0 }^{t}\beta(s)ds),\;(1-\exp(-\int_{0}^{t}\beta(s)ds)){\mathb{I}}\right)\;,$ (35) $$

which can be written as

$$ $\tilde{{\boldsymbol{\theta}}}^{(t)}=\sqrt{\bar{\alpha}_{t}}\tilde{{\boldsymbol {\theta}}}^{(0)}+\sqrt{1-\bar{\alpha}_{t}}{\boldsymbol{\epsilon}}\quad\text{ with}\quad\bar{\alpha}_{t}=\exp(-\int_{0}^{t}\beta(s)ds)\quad\text{and}$ (36) $$

where ${\boldsymbol{\epsilon}}\sim\mathcal{N}(\mathbf{0},{\mathb{I}})$ is a standard Gaussian random variable.
Samples from $p_{0}$ are transformed to samples from a standard Gaussian distribution after some large time $T$.

The reverse denoising process can then be written as the reverse of the diffusion process (34), which as stated by
Anderson, () is:

$$ $d{\boldsymbol{\theta}}^{(t)}=\left[-\frac{\beta(t)}{2}{\boldsymbol{\theta}}^{( t)}-\beta(t)\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}}^{(t)} )\right]dt+\sqrt{\beta(t)}d{{\mathb{B}}}_{t}\;,$ (37) $$

where $p_{t}$ is the distribution of $\tilde{{\boldsymbol{\theta}}}^{(t)}$ from
(34).
Solving this reverse SDE, the distribution of ${\boldsymbol{\theta}}^{(T)}$ is closed to $p_{0}$ for large $T$, which allows approximate sampling from $p_{0}$.

The score function $\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}}^{(t)})$ of the noisy data distribution at time $t$ is intractable and is then estimated by learning a neural network $s_{\phi}({\boldsymbol{\theta}},t)$ with parameters $\phi$. Score matching (Hyvärinen,, ) is a method to train $s_{\phi}$ by minimizing the following loss:

$$ $\mathbb{E}_{p_{t}({\boldsymbol{\theta}})}\left[||s_{\phi}({\boldsymbol{\theta} },t)-\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}})||^{2}\right ]\;.$ (38) $$

As $p_{t}({\boldsymbol{\theta}})$ is still unknown and only samples from $p_{t}({\boldsymbol{\theta}}|\tilde{{\boldsymbol{\theta}}}^{(0)})$ in (35) are available, Song et al., () rewrite this loss function as:

$$ $\mathbb{E}_{t\sim U[0,T]}\mathbb{E}_{p_{0}({\boldsymbol{\theta}}^{(0)})} \mathbb{E}_{p_{t}({\boldsymbol{\theta}}|{\boldsymbol{\theta}}^{(0)})}\left[ \lambda(t)||s_{\phi}({\boldsymbol{\theta}},t)-\nabla_{\boldsymbol{\theta}}\log p _{t}({\boldsymbol{\theta}}|{\boldsymbol{\theta}}^{(0)})||^{2}\right]$ (39) $$

where $\lambda(t)>0$ is a weighting function that allows to focus more on certain timesteps than others. It is common to take $\lambda(t)$ inversely proportional to the variance of $(\ref{eq:vp_noisy_sample})$ at time $t$.

Once the neural network $s_{\phi}$ has been trained by minimizing (39), it can be used to generate new samples approximately distributed as the target distribution $p_{0}$ by running a numerical scheme on the reverse SDE (37). By running for example the Euler-Maruyama scheme on (37), we get the following update step for the reverse process:

$$ ${\boldsymbol{\theta}}^{(t-\Delta t)}={\boldsymbol{\theta}}^{(t)}+\frac{\beta(t )}{2}{\boldsymbol{\theta}}^{(t)}\Delta t+\beta(t)s_{\phi}({\boldsymbol{\theta} }^{(t)},t)\Delta t+\sqrt{\beta(t)\Delta t}{\boldsymbol{\epsilon}}\;.$ (40) $$

We can then generate samples approximately from $p_{0}$ by running the reverse process (40) with a small enough $\Delta t$.

### C.2 Conditional Diffusion Models

Conditional diffusion models arise when, for some measurement ${\mathb{y}}$, we want to produce samples from some conditional distribution $p_{0}({\boldsymbol{\theta}}|{\mathb{y}})$.
Sampling from conditional distributions is a problem that arises in inverse problems. When using diffusion models, numerous solutions have been investigated as mentioned in a very recent review (Daras et al.,, ). We specify in this section the approach adopted for our applications.
With the application to experimental design in mind, we assume here that

$$ ${\mathb{Y}}={\mathb{A}}_{\boldsymbol{\xi}}{\boldsymbol{\theta}}+{\boldsymbol{ \eta}}$ (41) $$

where ${\boldsymbol{\eta}}\sim\mathcal{N}(\mathbf{0},\sigma^{2}{\mathb{I}})$ is the measurement noise, ${\mathb{A}}_{\boldsymbol{\xi}}$ is the operator that represents the experiment at ${\boldsymbol{\xi}}$.

Sampling from the conditional distribution $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$ can be done by running the reverse diffusion process on the conditional SDE:

$$ $d{\boldsymbol{\theta}}^{(t)}=\left[-\frac{\beta(t)}{2}{\boldsymbol{\theta}}^{( t)}-\beta(t)\nabla_{\boldsymbol{\theta}}\log p({\boldsymbol{\theta}}^{(t)}|{ \mathb{y}},{\boldsymbol{\xi}})\right]dt+\sqrt{\beta(t)}d{\mathb{B}}_{t},$ (42) $$

with the usual score $\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}}^{(t)})$ replaced by the conditionnal score $\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}}^{(t)}|{\mathb{y}}
,{\boldsymbol{\xi}})$. The main objective of conditional SDE is to generate samples from the conditional distribution $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$ without retraining a new neural network $s_{\phi}$ for the new conditional score.
By noting that the conditional score can be written as:

$$ $\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}}^{(t)}|{\mathb{y}} ,{\boldsymbol{\xi}})=\nabla_{\boldsymbol{\theta}}\log p_{t}({\mathb{y}}|{ \boldsymbol{\theta}}^{(t)},{\boldsymbol{\xi}})+\nabla_{\boldsymbol{\theta}} \log p_{t}({\boldsymbol{\theta}}^{(t)})$ (43) $$

we can leverage a pre-computed neural network $s_{\phi}$ that was trained to estimate the score $\nabla_{\boldsymbol{\theta}}\log p_{t}({\boldsymbol{\theta}}^{(t)})$ in the unconditional case. If we know how to evaluate the first term $\nabla_{\boldsymbol{\theta}}\log p_{t}({\mathb{y}}|{\boldsymbol{\theta}}^{(t)}
,{\boldsymbol{\xi}})$, we can then run the reverse process (42) to generate samples from the conditional distribution $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$.
Unfortunately this term does not have a closed form expression. As a solution, Dou and Song, () propose to approximate the intractable $\nabla_{\boldsymbol{\theta}}\log p_{t}({\mathb{y}}|{\boldsymbol{\theta}}^{(t)}
,{\boldsymbol{\xi}})$ by the tractable $\nabla_{\boldsymbol{\theta}}\log p_{t}({\mathb{y}}^{(t)}|{\boldsymbol{\theta}}
^{(t)},{\boldsymbol{\xi}})$ where ${\mathb{y}}^{(t)}$ is a noisy version of ${\mathb{y}}$ at time $t$. Then, the following backward SDE can be run to generate samples from the conditional distribution $p({\boldsymbol{\theta}}|{\mathb{y}},{\boldsymbol{\xi}})$:

$$ $d{\boldsymbol{\theta}}^{(t)}=\left[-\frac{\beta(t)}{2}{\boldsymbol{\theta}}^{( t)}-\beta(t)\nabla_{\theta}\log p_{t}({\mathb{y}}^{(t)}|{\boldsymbol{\theta}}^ {(t)},{\boldsymbol{\xi}})-\beta(t)\nabla_{\theta}\log p_{t}({\boldsymbol{ \theta}}^{(t)})\right]dt+\sqrt{\beta(t)}d{\mathb{B}}_{t}$ (44) $$

The sequence of noisy ${\mathb{y}}^{(t)}$ can be generated with a noising process like (36),

$$ ${\mathb{y}}^{(t)}=\sqrt{\bar{\alpha}_{t}}{\mathb{y}}+\sqrt{1-\bar{\alpha}_{t}} {\mathb{A}}_{\boldsymbol{\xi}}{\boldsymbol{\epsilon}}\quad\text{with}\quad\bar {\alpha}_{t}=\exp(-\int_{0}^{t}\beta(s)ds)\;,$ (45) $$

which using the forward model (41) can be written as:

$$ ${\mathb{y}}^{(t)}={\mathb{A}}_{\boldsymbol{\xi}}{\boldsymbol{\theta}}^{(t)}+ \sqrt{\bar{\alpha}_{t}}{\boldsymbol{\eta}}\;.$ (46) $$

We can then evaluate $\nabla_{\boldsymbol{\theta}}\log p_{t}({\mathb{y}}^{(t)}|{\boldsymbol{\theta}}
^{(t)},{\boldsymbol{\xi}})$ as :

$$ $\nabla_{\boldsymbol{\theta}}\log p_{t}({\mathb{y}}^{(t)}|{\boldsymbol{\theta}} ^{(t)},{\boldsymbol{\xi}})=\frac{1}{\sigma^{2}\bar{\alpha}_{t}}{\mathb{A}}_{ \boldsymbol{\xi}}^{T}({\mathb{y}}^{(t)}-{\mathb{A}}_{\boldsymbol{\xi}}{ \boldsymbol{\theta}}^{(t)})\;.$ (47) $$

## Appendix D Sequential Bayesian experimental design

In this framework, experimental conditions are determined sequentially, making use of measurements that are gradually made. This sequential view is referred to as sequential or iterated design.
In a sequential setting, we assume that we plan a sequence of $K$ experiments. For each experiment, we wish to pick the best ${\boldsymbol{\xi}}_{k}$ using the data that has already been observed
${\mathb{D}}_{k-1}=\{({\mathb{y}}_{1},{\boldsymbol{\xi}}_{1}),\ldots,({\mathb{y
}}_{k-1},{\boldsymbol{\xi}}_{k-1})\}$. Given this design, we conduct an experiment using ${\boldsymbol{\xi}}_{k}$ and obtain
outcome ${\mathb{y}}_{k}$. Both ${\boldsymbol{\xi}}_{k}$ and ${\mathb{y}}_{k}$ are then added to ${\mathb{D}}_{k-1}$ for a new set ${\mathb{D}}_{k}={\mathb{D}}_{k-1}\cup({\mathb{y}}_{k},{\boldsymbol{\xi}}_{k})$. After each step, our belief about ${\boldsymbol{\theta}}$ is updated and summarised by the current posterior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k})$, which acts as the next prior at step $k+1$. When the observations are assumed conditionally independent, it comes,

$$ $\displaystyle p({\boldsymbol{\theta}}|{\mathb{D}}_{k})$ $\displaystyle\propto p({\boldsymbol{\theta}})\;\prod_{n=1}^{k}p({\mathb{y}}_{n }|{\boldsymbol{\theta}},{\boldsymbol{\xi}}_{n})$ (48) $$

and

$$ $\displaystyle p({\mathb{y}},{\boldsymbol{\theta}}|{\boldsymbol{\xi}},{\mathb{D }}_{k-1})$ $\displaystyle\propto p({\boldsymbol{\theta}})\;p({\mathb{y}}|{\boldsymbol{ \theta}},{\boldsymbol{\xi}})\;\prod_{n=1}^{k-1}p({\mathb{y}}_{n}|{\boldsymbol{ \theta}},{\boldsymbol{\xi}}_{n})\;.$ (49) $$

A greedy design can be seen as choosing each design ${\boldsymbol{\xi}}_{k}$ as if it was the last one. This means that ${\boldsymbol{\xi}}_{k}$ is chosen as ${\boldsymbol{\xi}}_{k}^{*}$ the value that maximizes

$$ $\displaystyle{\boldsymbol{\xi}}_{k}^{*}$ $\displaystyle=\arg\max_{\boldsymbol{\xi}}I_{k}({\boldsymbol{\xi}},{\mathb{D}}_ {k-1})$ $$

where

$$ $\displaystyle I_{k}({\boldsymbol{\xi}},{\mathb{D}}_{k-1})$ $\displaystyle=\mathbb{E}_{p_{\boldsymbol{\xi}}^{k}}\left[\log\frac{p_{ \boldsymbol{\xi}}^{k}({\boldsymbol{\theta}},{\mathb{Y}})}{p({\mathb{Y}}|{ \boldsymbol{\xi}},{\mathb{D}}_{k-1})\;p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1 })}\right]\;=\;\text{MI}(p_{\boldsymbol{\xi}}^{k})$ (50) $$

with $p_{\boldsymbol{\xi}}^{k}$ denoting the joint distribution $p({\mathb{y}},{\boldsymbol{\theta}}|{\boldsymbol{\xi}},{\mathb{D}}_{k-1})=p({
\mathb{y}}|{\boldsymbol{\theta}},{\boldsymbol{\xi}})\;p({\boldsymbol{\theta}}|
{\mathb{D}}_{k-1})$. Distribution
$p_{\boldsymbol{\xi}}^{k}$ involves the current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$, which is not available in closed-form and is not straightforward to sample from.
Distribution $p_{\boldsymbol{\xi}}^{k}$ can be written as a Gibbs distribution by defining the potential $V_{k}$ as

$$ $\displaystyle p_{\boldsymbol{\xi}}^{k}({\mathb{y}},{\boldsymbol{\theta}})$ $\displaystyle\propto\exp\left(-V_{k}({\mathb{y}},{\boldsymbol{\theta}},{ \boldsymbol{\xi}})\right)$ $\displaystyle\mbox{with }V_{k}({\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{ \xi}})$ $\displaystyle=-\log p({\boldsymbol{\theta}})-\log p({\mathb{y}}|{\boldsymbol{ \theta}},{\boldsymbol{\xi}})-\sum_{n=1}^{k-1}\log p({\mathb{y}}_{n}|{ \boldsymbol{\theta}},{\boldsymbol{\xi}}_{n})$ $\displaystyle=V({\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\xi}})+\tilde{V }_{k}({\boldsymbol{\theta}})\;,$ $$

where $V({\mathb{y}},{\boldsymbol{\theta}},{\boldsymbol{\xi}})$ has been already defined in Section
4.
Note that the marginal in ${\boldsymbol{\theta}}$ of $p_{\boldsymbol{\xi}}^{k}$ is the posterior at step $k-1$ or equivalently the current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ and the marginal in ${\mathb{y}}$ is

$$ $p({\mathb{y}}|{\boldsymbol{\xi}},{\mathb{D}}_{k-1})=\mathbb{E}_{p({\boldsymbol {\theta}}|{\mathb{D}}_{k-1})}\left[p({\mathb{y}}|{\boldsymbol{\theta}},{\mathb {D}}_{k-1})\right].$ $$

Once a new ${\boldsymbol{\xi}}_{k}$ is computed and a new observation ${\mathb{y}}_{k}$ is performed,
the posterior at step $k$ is $p({\boldsymbol{\theta}}|{\mathb{y}}_{k},{\boldsymbol{\xi}}_{k},{\mathb{D}}_{k-
1})$ which is the conditional distribution of $p({\mathb{y}}_{k},{\boldsymbol{\theta}}|{\boldsymbol{\xi}}_{k},{\mathb{D}}_{k-
1})$.

## Appendix E Sequential Monte Carlo (SMC)-style resampling

SMC is an essential addition when dealing with sequential BOED.
In density-based BOED, it has been already exploited in the sequential context showing a real improvement in the quality of the generated samples (Iollo et al.,, ).
SMC is also useful in simpler static cases as it can improve the quality of the generated ${\boldsymbol{\theta}}$ and contrastive ${\boldsymbol{\theta}}^{\prime}$ samples, that in turn improves the accuracy of the gradient estimator (30).
A particularly central step in SMC is the resampling step, first recalled below in the density-based case.
Using our framework, is it also possible to derive a SMC-style resampling scheme in the data-based case. This is becoming a popular strategy in the context of generative models (Dou and Song,, ; Cardoso et al.,, ).

##### Density-based BOED.

In static density-based BOED, the prior $p({\boldsymbol{\theta}})$ and the likelihood $p({\mathb{y}}|{\boldsymbol{\theta}},{\boldsymbol{\xi}})$ are available in closed-form. In the sequential experiment context, we want to generate $N$ samples ${\boldsymbol{\theta}}_{1},\ldots,{\boldsymbol{\theta}}_{N}$ from the current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ and $M$ samples ${\boldsymbol{\theta}}^{\prime}_{1},\ldots,{\boldsymbol{\theta}}^{\prime}_{M}$ from the expected posterior $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}}^{\prime})$. As both $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ and $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}}^{\prime})\propto\prod_{i=1}^{N}
p({\boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}},{\mathb{D}
}_{k-1})^{\nu_{i}}$ can be evaluated up to a normalizing constant, it is straightforward to extend the sampling operators of Section 4 and add a resampling step to the samples ${\boldsymbol{\theta}}_{1},\ldots,{\boldsymbol{\theta}}_{N}$ and ${\boldsymbol{\theta}}^{\prime}_{1},\ldots,{\boldsymbol{\theta}}^{\prime}_{M}$ with weights $w_{i}$ and $w^{\prime}_{j}$:

$$ $\displaystyle w_{i}$ $\displaystyle=\frac{\tilde{w_{i}}}{\sum_{i=1}^{N}\tilde{w_{i}}}\quad\text{with }\quad\tilde{w_{i}}=\tilde{p}({\boldsymbol{\theta}}_{i})$ $\displaystyle w^{\prime}_{j}$ $\displaystyle=\frac{\tilde{w^{\prime}_{j}}}{\sum_{j=1}^{M}\tilde{w^{\prime}_{j }}}\quad\text{with}\quad\tilde{w^{\prime}_{j}}=\tilde{q}({\boldsymbol{\theta}} ^{\prime}_{j})$ $$

where $\tilde{p}$ and $\tilde{q}$ are the unnormalized versions of $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ and $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}}^{\prime})$ respectively.

##### Data-based BOED.

In the setting of data-based BOED, we assume access to a conditional diffusion model that allows to generate samples from $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ and $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}}^{\prime})$. The resampling scheme proposed in (Dou and Song,, ) can be used as is, to improve the quality of the samples from $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ as this is an usual conditional distribution.
The resampling scheme is based on the FPS update:

$$ $p({\boldsymbol{\theta}}^{(t-1)}|{\boldsymbol{\theta}}^{(t)},{\mathb{y}}^{(t-1) },{\boldsymbol{\xi}},{\mathb{D}}_{k-1})\propto p_{t}({\boldsymbol{\theta}}^{(t -1)}|{\boldsymbol{\theta}}^{(t)},{\mathb{D}}_{k-1})\;p({\mathb{y}}^{(t-1)}|{ \boldsymbol{\theta}}^{(t-1)},{\boldsymbol{\xi}})$ (51) $$

where $p_{t}({\boldsymbol{\theta}}^{(t-1)}|{\boldsymbol{\theta}}^{(t)},{\mathb{D}}_{k
-1})$ is given in closed form by the unconditional diffusion model and $p({\mathb{y}}^{(t-1)}|{\boldsymbol{\theta}}^{(t-1)},{\boldsymbol{\xi}})$ is given by (46).
As both these distributions are Gaussian, $p({\boldsymbol{\theta}}^{(t-1)}|{\boldsymbol{\theta}}^{(t)},{\mathb{y}}^{(t-1)
},{\boldsymbol{\xi}},{\mathb{D}}_{k-1})$ can be written in closed form and resampling weights can be written as:

$$ $w_{i}=\frac{\tilde{w_{i}}}{\sum_{i=1}^{N}\tilde{w_{i}}}\quad\text{with}\quad \tilde{w_{i}}=p({\mathb{y}}^{(t-1)}|{\boldsymbol{\theta}}_{i}^{(t)},{ \boldsymbol{\xi}})$ (52) $$

where $p({\mathb{y}}^{(t-1)}|{\boldsymbol{\theta}}_{i}^{(t)},{\boldsymbol{\xi}})$ is tractable (see Dou and Song, () for more details).

For the expected posterior $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}}^{\prime})\propto\prod_{i=1}^{N}
p({\boldsymbol{\theta}}^{\prime}|{\mathb{y}}_{i},{\boldsymbol{\xi}},{\mathb{D}
}_{k-1})^{\nu_{i}}$, update (51) takes the form:

$$ $\displaystyle\prod_{i=1}^{N}p({\boldsymbol{\theta}}^{(t-1)}|{\boldsymbol{ \theta}}^{(t)},{\mathb{y}}^{(t-1)}_{i},{\boldsymbol{\xi}},{\mathb{D}}_{k-1})^{ \nu_{i}}$ $\displaystyle\propto p_{t}({\boldsymbol{\theta}}^{(t-1)}|{\boldsymbol{\theta}} ^{(t)},{\mathb{D}}_{k-1})\;\prod_{i=1}^{N}p({\mathb{y}}^{(t-1)}_{i}|{ \boldsymbol{\theta}}^{(t-1)},{\boldsymbol{\xi}})^{\nu_{i}}$ $\displaystyle\propto\prod_{i=1}^{N}\left(p({\boldsymbol{\theta}}^{(t-1)}|{ \boldsymbol{\theta}}^{(t)},{\mathb{D}}_{k-1})\;p({\mathb{y}}^{(t-1)}_{i}|{ \boldsymbol{\theta}}^{(t-1)},{\boldsymbol{\xi}})\right)^{\nu_{i}}$ (53) $$

which leads to the following resampling weights:

$$ $w^{\prime}_{j}=\frac{\tilde{w^{\prime}_{j}}}{\sum_{j=1}^{M}\tilde{w^{\prime}_{ j}}}\quad\text{with}\quad\tilde{w^{\prime}_{j}}=\prod_{i=1}^{N}p({\mathb{y}}^{ (t-1)}_{i}|{\boldsymbol{\theta}}^{\prime(t)}_{j},{\boldsymbol{\xi}})^{\nu_{i}}\;.$ $$

## Appendix F Numerical experiments

### F.1 Sequential prior contrastive estimation (SPCE) and Sequential nested Monte Carlo (SNMC) criteria

The SPCE introduced by Foster et al., () is a tractable quantity to assess the design sequence quality.
For a number $K$ of experiments, ${\mathb{D}}_{K}=\{({\mathb{y}}_{1},{\boldsymbol{\xi}}_{1}),\cdot,({\mathb{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({\mathb{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({\mathb{y}}_{k}|{ \boldsymbol{\theta}}_{0},{\boldsymbol{\xi}}_{k})}{\frac{1}{L+1}\sum\limits_{ \ell=0}^{L}\prod\limits_{k=1}^{K}p({\mathb{y}}_{k}|{\boldsymbol{\theta}}_{\ell },{\boldsymbol{\xi}}_{k})}\right]\;.$ (54) $$

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 $\{{\mathb{y}}_{k}\}_{k=1:K}$ are simulated using $p({\mathb{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({\mathb{y}}_{k}|{\boldsymbol{\theta}}^{*},{\boldsymbol{\xi}}_{k})}{ \frac{1}{L+1}\left(\prod\limits_{k=1}^{K}p({\mathb{y}}_{k}|{\boldsymbol{\theta }}^{*},{\boldsymbol{\xi}}_{k})+\sum\limits_{\ell=1}^{L}\prod\limits_{k=1}^{K}p ({\mathb{y}}_{k}|{\boldsymbol{\theta}}^{i}_{\ell},{\boldsymbol{\xi}}_{k}) \right)}\right\}\;.$ $$

Similarly, an upper bound on the total EIG has also been introduced by Foster et al., ()
and named the Sequential nested Monte Carlo (SNMC) criterion,

$$ $\displaystyle SNMC({\boldsymbol{\xi}}_{1},\cdot,{\boldsymbol{\xi}}_{K})$ $\displaystyle=\mathbb{E}_{\prod\limits_{k=1}^{K}p({\mathb{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({\mathb{y}}_{k}|{ \boldsymbol{\theta}}_{0},{\boldsymbol{\xi}}_{k})}{\frac{1}{L}\sum\limits_{\ell =1}^{L}\prod\limits_{k=1}^{K}p({\mathb{y}}_{k}|{\boldsymbol{\theta}}_{\ell},{ \boldsymbol{\xi}}_{k})}\right]\;.$ $$

As shown in Foster et al., () (Appendix A), SPCE increases with $L$ to reach the total EIG
when $L\rightarrow\infty$ at a rate ${\cal O}(L^{-1})$ of convergence.
It is also shown in Foster et al., () that for a given $L$, SPCE is bounded by $\log(L+1)$ while the upper bound SNMC below is potentially unbounded. As in Blau et al., (), 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.

### F.2 Implementation details

#### F.2.1 Source example

For VPCE (Foster et al.,, ) and RL-BOED (Blau et al.,, ), we
used the code available at [github.com/csiro-mlai/RL-BOED](https://github.com/csiro-mlai/RL-BOED/tree/master), using the settings recommended therein to reproduce the results in the respective papers.
VPCE optimizes an EIG lower bound in a myopic manner estimating posterior distributions with variational approximations. RL-BOED is a non-myopic approach which does not provide posterior distributions.
From the obtained sequences of observations and design values, we computed SPCE and SNMC as explained above and retrieved the same results as in their respective papers.
For PASOA and SMC procedures, we used the code available at [github.com/iolloj/pasoa](https://github.com/iolloj/pasoa).
PASOA is a myopic approach, optimizing an EIG lower bound using sequential Monte Carlo (SMC) samplers and tempering to also provide posterior estimations. The method refered to as SMC is a variant without tempering.

For CoDiff, the $\nu_{i}$’s in the expected posterior distribution were set to $\nu_{i}=\frac{1}{N}$. The current prior and posterior distributions at experimental step $k$ were initialized using respectively the prior and posterior samples at step $k-1$. Design optimization was performed using the Adam optimizer with an exponential learning rate decay schedule with initial learning rate $10^{-2}$ and decay rate $0.98$. The Langevin step-size in the DiGS method Chen et al., () was set to $10^{-2}$. The joint optimization-sampling loop was run for 5000 steps.

Figure: Figure 5: Several source localisation examples. Prior (left) and expected posterior (right) samples at experiment $k$. Final ${\boldsymbol{\xi}}_{k}^{*}$ value (orange cross) at the end of the optimization sequence ${\boldsymbol{\xi}}_{0},\cdot,{\boldsymbol{\xi}}_{T}$ (blue crosses). This optimization "contrasts" the two distributions by making the expected posterior "as different as possible" from the prior.
Refer to caption: extracted/5929226/img/expected/expected_1.png

For this example, the additional Figure 5 shows samples from the current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$ at some intermediate step $k$, which gradually concentrates around the true sources. Samples from the expected posterior distribution are shown in comparison to illustrate its contrastive nature.

#### F.2.2 MNIST example

For the numerical example of Section 6.3, we used the MNIST dataset (LeCun et al.,, ), the time varying SDE (34) with a noise schedule $(b_{max}-b_{min})/(T-t_{0})+(b_{min}-b_{max})t_{0}/(T-t_{0})$ (with $b_{max}=5$, $b_{min}=0.2,t_{0}=0,T=2$).
The training of the usual score matching was done for 3000 epochs with a batch size of 256 and using Adam optimizer Kingma and Ba, (). We used gradient clipping and the training was done on a single A100 GPU.

Update equations for the sampling operators were derived from SDE (19) for the contrastive samples of the expected posterior $q_{{\boldsymbol{\xi}},N}({\boldsymbol{\theta}}^{\prime})$ and (44) for samples from the current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$, where ${\mathb{D}}_{k-1}$ can be added in the conditioning part without difficulty.
Those updates are equivalent to (53) and (51) respectively. The resampling weights were computed as in Section E.

Figures 6 and 7 show additional image reconstruction processes. The digit to be recovered is shown in the first column. The successively selected masks are shown (red line squares) in the second column with the resulting gradually discovered part of the image.
The reconstruction per se can be estimated from the posterior samples shown in the last 16 columns. At each experiment, the upper sub-row shows the 16 most-likely reconstructed images, while the lower sub-row shows the 16 less-probable ones.
As the number of experiments increases the posterior samples gradually concentrate on the right digit.

Figure 8 then shows that design optimization is effective by showing better outcomes when masks locations are optimized (second column) than when masks are selected at random centers (third column). The highest posterior weight samples in the last 14 columns also clearly show more resemblance with the true digit in the optimized case.

Figure: Figure 6: Image reconstruction. First 7 experiments (rows): image ground truth, measurement at experiment $k$, samples from current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$, with best (resp. worst) weights in upper (resp. lower) sub-row. The samples incorporate past measurement information as the procedure advances.
Refer to caption: extracted/5929226/img/mnist_44/samples_0.png

Figure: Figure 7: Image reconstruction. First 7 experiments (rows): image ground truth, measurement at experiment $k$, samples from current prior $p({\boldsymbol{\theta}}|{\mathb{D}}_{k-1})$, with best (resp. worst) weights in upper (resp. lower) sub-row. The samples incorporate past measurement information as the procedure advances.
Refer to caption: extracted/5929226/img/mnist_444/samples_0.png

Figure: Figure 8: Image ${\boldsymbol{\theta}}$ (1st column) reconstruction from 7 sub-images ${\mathb{y}}={\mathb{A}}_{\boldsymbol{\xi}}{\boldsymbol{\theta}}+{\boldsymbol{ \eta}}$ selected sequentially at 7 central pixel ${\boldsymbol{\xi}}$. Optimized vs. random designs: measured outcome ${\mathb{y}}$ (2nd vs. 3rd column) and parameter ${\boldsymbol{\theta}}$ estimates (reconstruction) with highest weights (upper vs. lower sub-row).
Refer to caption: extracted/5929226/img/mnist_44/comparison.png

### F.3 Hardware details

The source example 6.2 can be run locally. It was tested on an Apple M1 Pro 16Gb chip but faster running times can be achieved on GPU. The MNIST example 6.3 was run on a single A100 80Gb GPU.

### F.4 Software details

Our code is implemented in Jax Bradbury et al., () and uses Flax as a Neural Network library and Optax as optimization one Babuschkin et al., ().

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: See the "Weaknesses" section.
soundness: 3
strengths: The paper advances the field of BOED by making it more computationally efficient, whilst broadening the scope of problems that can be tackled by BOED by extending it to diffusion-based generative models. The idea of performing joint sampling and optimization using bi-level optimization seems novel to, although I am not very familiar with this literature. The experimental results are significant and demonstrate the efficacy of CoDiff compared to other baselines. The method is presented clearly along with discussions w.r.t. to the related works. Overall, this is a well-written and technically solid paper, and I am happy to recommend its acceptance.
summary: This paper presents a new expression for the gradient of the expected information gain (EIG)---a central quantity in Bayesian optimal experimental design (BOED). Their method, named CoDiff, combines the two steps of maximizing the EIG and then sampling from the posterior in BOED into a single joint step by leveraging the sampling-as-optimization setting of e.g. Marion et al. (2024). Doing so results in a more efficient BOED method, as demonstrated by numerical experiments.
weaknesses: Some minor points:
* It would be nice to see some uncertainty bars around the results in Figure 3.
* Parts of Section 4 and 5 were a bit difficult for me to follow, but that could just be because I am not familiar with the relevant literature. However, if possible, it might be a good exercise for the authors to think about making these parts a bit more accessible to the BOED community who are not familiar with the sampling-as-optimization literature.
* It would have been nice to empirically test the performance of CoDiff as e.g. the dimensions of $\theta$ or the hyperparameters (if any) vary.
* Please include some limitations of CoDiff in the conclusion section.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: - Can the authors provide an intuition as to why the proposed expected posterior is preferable to the true expected posterior = prior? Is this purely for the sequential setting?
- how are the $\nu_i$ values determined in the expected posterior?
- what stop other methods from being used on the MNIST example? Is it just the computational cost from multiple nested loops and the sequential use case?
soundness: 3
strengths: - rigorously cites and discusses related ideas
- merging multiple nested loops into a single loop is intuitive idea and has been applied in multiple separate fields (training a VAE is effectively the EM algorithm was simplified into a single step, in Bayesian optimzation few-shot Knowledge Gradient)
- the writing is clear and each idea is easy enough to follow.
summary: This paper is a little outside of my expertise and I apologize in advance for any misunderstanding I may have and look forward to any correction from the authors. I have set my confidence to low accordingly.

The authors consider the problem of experimental design, given a statistical model of the form $P(y | \xi, \theta)$, we can choose values of $\xi$ and observe a corresponding $y$, and the goal is to infer $\theta$. In the first example in the paper, $\xi$ are locations on a 2D plane, $y$ are the measurements at the $\xi$ locations, and $\theta$ represents the source locations of contamination in the plane.

Assuming we have chosen the $\xi$ values, we can use the $p(y|\xi)$ to simulate $y_i$ values, one hallucinated rollout into the future, then infer the hallucinated posterior $p(\theta'|y_i, \xi)$ and measure the posterior distribution's entropy, $H(p(\theta'|y_i,\xi))$, typically requiring Monte-Carlo over the $\theta'$. We may repeat this process many many times for $i=1,...,N$, and compute the average posterior entropy, wrapping the inner Monte-Carlo over $\theta'$ in another outer Monte-Carlo over $y$. This average posterior entropy is the Mutual Information (MI) between $y$ and $\theta$ up to an additive constant (that is the prior entropy $H(p(\theta))$). This MI quantifies the benefit of sampling $y$ at $\xi$, we can therefore change $\xi$ in a way to maximise MI and this is a standard method for experimental design known as Expected Information Gain, or EIG.

However as evaluating EIG requires performing nested Monte-Carlo, and optimizing $\xi$ requires the gradient of EIG, much work has focused on various methods to construct efficient EIG gradient estimators. In past works and this works, the second inner integration over $\theta'$ is performed with MCMC or by Monte Carlo with importance sampling using a proposal $q(\theta'|y_i)$ that must be updated for every outer MC iteration $i$.

This work proposes to use a single proposal distribution across all outer MC iterations $i$ called the "expected posterior" which is a geometric mixture of all the hallucinated posteriors
$$
q(\theta') \propto \prod_i p(\theta|y_i, \xi)^{\nu_i}
$$

The paper proposes to Langevin dynamics to generate samples for the outer Monte Carlo $\\{y_i, \theta_i\\}$ from the joint distribution. As these samples are updated in a iterative manner, and the goal is to iteratively optimize $\xi$, these nested iterative procedures can be merged and updated together, if $\xi$ only incrementally moves, it is wasteful ti discard all the old samples $\\{y_i, \theta_i\\}$ but keep incrementing them so they're up-to-date.

Finally the proposed method is applied to two sequential design applications, the first is the source location example in which the $p(\theta)$, $p(y|\theta, \xi)$ have known analytical forms and the proposed method outperforms baselines. The second is where $p(\theta)$ cannot be evaluated but can be sampled, specifically a diffusion model trained on MNIST images, $\theta$ is a full image,  $\xi$ are pixel locations, and $y$ are the 7 X 7 patch of pixel values around $\xi$. The goal is to infer the full image by cherry-picking image patches.
weaknesses: ### Technical Comments
- __expected posterior__ the paper states that they aim to find one proposal distribution that covers all possible future posteriors and proposes a distribution that is a geometric mixing of hallucinated posteriors $q(\theta) \propto \prod_i p(\theta|y_i, \xi)^{\nu_i}$. The prior is the true expected posterior $p(\theta) = \mathbb{E}_Y[p(\theta|y)]$, exactly the average density of all possible posterior densities.  Why not use the prior $p(\theta)$ as the one global proposal distribution for integration over all the hallucinated posteriors? In the sequential case (introduced at the end of the paper)

### Writing Comments
- __Writing Density__ while the writing is thankfully clear, the paper cites and quickly introduces many ideas from prior works in quick succession making the paper very hard work to read for me, a lot of context switching making the paper very tiring. I would suggest reducing the details of referenced works, and providing higher level descriptions where possible and elaborating more on the novel parts of the proposed method.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: My questions are contained in the Weaknesses section.
soundness: 3
strengths: - The authors make a great connection to optimization by sampling via diffusion models. 
- Thorough gradient analysis.
- The connection to inverse problems e.g. inpainting in MNIST images as a BOED problem is an interesting new task in the BOED setting and could see BOED objectives (really mutual information objectives) gain further adoption in generative modeling.
summary: This paper makes an interesting connection between Bayesian optimal experimental design (BOED) and optimization by sampling using a diffusion model. While an interesting connection that may prove fruitful with some refinement, the authors seem to misunderstand the basis of BOED and maximization of a lower bound of mutual information using contrastive learning, as well as miss prior work that might better-inform this method. This disconnect, as well as confusing notation regrading how the diffusion model is used in multi-round BOED and lackluster empirical evaluation, makes me think this work is too early and could benefit from refinement in writing and additional experimental evaluations of the method.
weaknesses: While the authors make a great connection to optimization by sampling, theoretical analysis of the convergence of the desired goal distribution is missing. This would e.g. especially be helpful in the multi-step BOED setting to show how to improve bounds of the mutual information. I list out my critiques and questions below:

- The authors miss that they are also optimizing a lower bound of the mutual information by the InfoNCE bound, thus their comments that their approximation is exact is incorrect.
- The definition of information in Equation 2 does not make sense. The KL divergence is written in a non-standard way - do you mean $E_{p(y|\xi)}KL(p(\theta|y, \xi) || p(\theta))$?
- The authors miss the connection to Fotster et al. 2019b of LF-ACE that uses a posterior estimate in the InfoNCE bound to refine samples. e.g. in equation 5. 
- The notation for the expected posterior in equation 9, it's not clear what is the proposal distribution $q(\theta)$ in the paper. Are you using likelihood evaluations? Also, how are you choosing the weighting of $\nu_i$?
- I would prefer to see the notation for iterative sampling operators similar to Marion et al. 2024, where the Y, and $\Theta$ are variables of the operator instead of in the superscript. I think that would make the optimization problem more clear. 
- After equation 15, are you using two generative models to define the observed distribution $y$ and the parameter distribution $\theta$?
- in terms of the location finding experiment, usually we would expect to see the designs and final posterior distribution. the way it's portrayed in the paper is confusing. I'm not sure which samples are shown, which design round the designs come from, and why the seeming optimal design isn't on the ground truth values. It looks like the chosen design is quite far from the true theta_0.
- While the EIG seems to improve over others in Figure 3, the upper bound is _much_ worse. Why is that? 
- The MNIST experiment is interesting but lacks quantitative evaluation. Maybe including a classifier and using the number of successfully classified images would help? Also, there seems to be opportunity to change the design size to see how that influences outcomes.

In summary, I think the notation is a little bloated. The authors tend to repeat themselves in the text and could be more concise. Concision would help make space for exposition of more salient points that would improve this paper. Additionally, it would help the authors to explicitly lay out what are the likelihood, posterior, and priors they use in experiments. Also, it would be nice to see calibration of the posterior that results of the final design round using e.g. the l-C2ST method or other simulation based calibration methods. Most applications of amortized posteriors using BOED will be interested in this, but it can be computationally expensive, which is why I think just on the final round is sufficient. This is an interesting method that I think needs cleaner exposition.

Meta-review: The paper address the problem of Bayesian optimal experimental design (BOED), where the expected information gain EIG is the central quantity of interest. The paper presents CoDiff, a method that leverages the sampling-as-optimization approach to do the maximization of the EIG and sampling from the posterior. The paper improves the computational efficiency of the BOED and broadens its scope by incorporating diffusion models.   All the reviewers recommend acceptance and I agree.

Decision: Accept (Spotlight)