$\displaystyle\geq\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{q^{\text{pos}}(\Theta|Y,\xi)}{p(\Theta)}\right]=\mathbb{E}_{Y,\Theta|\xi}[\log
q^{\text{pos}}(\Theta|Y,\xi)]+H(\Theta).$ (63)

Hence, for any approximation $q^{\text{pos}}(\theta|y,\xi)$ of the posterior
density, (63) is a lower bound on the mutual information (EIG). This bound is
sometime called the Barber–Agakov bound, after Barber and Agakov (2003). Note
that all expectations above, and specifically in (62) and (63), are with
respect to the true distribution $p_{Y,\Theta|\xi}$.

Evaluating the upper bound (62) requires the ability to evaluate the
likelihood function, and thus it does not apply to the implicit model setting
where the likelihood is intractable. Evaluating the lower bound, on the other
hand, only requires access to the prior density (or the differential entropy
of the prior) and the ability to find a tractable approximation
$q^{\text{pos}}(\theta|y,\xi)$.

Foster et al. (2019) were the first to suggest using these mutual information
bounds in OED, and in practice selected the approximations $q$ from simple
parametric families of densities (e.g., Gaussian, uniform) that were tailored
to the design problem at hand. Once such a family $\mathcal{Q}$ is specified,
the best member of the family—i.e., the density yielding the closest
approximation of the EIG—can be found by tightening the bound. Specifically,
for the marginal approximation, we seek (for any given design $\xi$)

$\displaystyle
q^{\text{mar},\ast}\in\operatornamewithlimits{argmax}_{q\in\mathcal{Q}}\mathbb{E}_{Y|\xi}\left[\log
q(Y|\xi)\right],$ (64)

which minimizes the upper bound (62), while for the posterior approximation we
seek

$\displaystyle
q^{\text{pos},\ast}\in\operatornamewithlimits{argmax}_{q\in\mathcal{Q}}\mathbb{E}_{Y,\Theta|\xi}\left[\log
q(\Theta|Y,\xi)\right],$ (65)

which maximizes the lower bound (63). In practice, the expectations in (64) or
(65) are approximated using samples from the model; for example, (64) becomes

$\displaystyle\hat{q}^{\text{mar}}\in\operatornamewithlimits{argmax}_{q\in\mathcal{Q}}\sum_{i=1}^{M}\log
q(y^{(i)}|\xi),\quad y^{(i)}\sim p(y|\xi),$ (66)

and analogously for (65). Inspecting (66), it is apparent that identifying a
member of the variational family in this way is none other than _maximum
likelihood estimation_ of either the marginal density $p(y|\xi)$ or
conditional density $p(\theta|y,\xi)$.333‘Maximum likelihood’ here refers to
the density estimation problem immediately at hand, i.e., estimating the best
$q\in\mathcal{Q}$ given samples, and should not be confused with the idea of
estimating $\theta$ by maximizing $\theta\mapsto p(y|\theta,\xi)$.

With this link in mind, we can immediately generalize the machinery used to
construct good density approximations. One powerful approach for estimating
both joint and conditional densities rests on transportation of measure
(Villani 2009, Marzouk et al. 2016, Spantini et al. 2018). Given some generic
target distribution $\pi$ on $\mathbb{R}^{d}$, the idea behind transport
methods is to find an invertible transformation
$S:\mathbb{R}^{d}\to\mathbb{R}^{d}$ that pushes forward $\pi$ to a simple
reference distribution $\rho$ on $\mathbb{R}^{d}$ whose density can be easily
evaluated. Abusing notation (by not distinguishing measures from densities),
we can write this as $S_{\sharp}\pi=\rho$, which means that
$\rho(A)=\pi\left(S^{-1}(A)\right)$ for any $\rho$-measurable set $A$.
Crucially, for any diffeomorphism $\tilde{S}$ on $\mathbb{R}^{d}$, the
distribution $\tilde{S}^{\sharp}\rho\coloneqq\tilde{S}^{-1}_{\sharp}\rho$,
called the pullback of $\rho$ under $\tilde{S}$, has a closed-form expression
for its density:

$\displaystyle\tilde{S}^{\sharp}\rho=(\rho\circ\tilde{S})\,\det\nabla\tilde{S},$
(67)

which is guaranteed to be positive and to integrate to one. Density estimation
can thus be recast as the problem of finding a map $\tilde{S}$ in some
suitable class such that $\tilde{S}^{\sharp}\rho$ is ‘close’ to $\pi$ (Wang
and Marzouk 2022).

In practice, these models can be quite expressive. Any measure absolutely
continuous with respect to the Lebesgue measure on $\mathbb{R}^{d}$ can be
represented as the pullback, under some map, of a Gaussian $\rho$ on
$\mathbb{R}^{d}$; in fact, for any pair of equivalent continuous measures
$(\rho,\pi)$, there exist infinitely many transport maps $S$ that achieve
$S^{\sharp}\rho=\pi$. Normalizing flows (Kobyzev et al. 2020, Papamakarios et
al. 2021) are a special case of this formulation, i.e., a particular class of
transport map parameterizations that guarantee invertibility,
differentiability, and easy evaluation of the Jacobian determinant $\det\nabla
S$. But many other representations are useful. Monotone triangular maps
(Bogachev et al. 2005), for instance, can represent arbitrary absolutely
continuous distributions and be parameterized in a way that endows the maximum
likelihood estimation problems (64)–(65) with optimization guarantees; see
Baptista et al. (2023) for details. Continuous optimal transport maps can also
be estimated by first solving a discrete optimal transport problem (i.e.,
between empirical measures) and smoothing the result (Manole et al. 2021,
Pooladian and Niles-Weed 2021) or by parameterizing a differentiable convex
potential (Huang et al. 2020).

In the context of OED, it is useful to employ _block-triangular_ maps
(Baptista et al. 2024b) (a class which includes strictly triangular maps but
infinitely many other choices), as they naturally capture conditional
densities. Consider the following block arrangement for a map
$S:\mathbb{R}^{n+p}\to\mathbb{R}^{n+p}$:

$\displaystyle S(y,\theta)=\begin{bmatrix}[l]S^{Y}(y)\\\
S^{\Theta}(y,\theta)\end{bmatrix},$ (68)

where $S^{Y}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and
$S^{\Theta}:\mathbb{R}^{n+p}\to\mathbb{R}^{p}$. Let the reference distribution
$\rho$ on $\mathbb{R}^{n+p}$ factor as $\rho=\rho_{n}\otimes\rho_{p}$, where
the subscripts denote the dimension of the factors. (The standard Gaussian
distribution on $\mathbb{R}^{n+p}$ naturally factorizes in this way, but any
distribution with appropriate block-independence would suffice.) Then, as
shown in Baptista et al. (2024b) and Marzouk et al. (2016), if
$S_{\sharp}p_{Y,\Theta}=\rho$, then $S^{Y}_{\sharp}p_{Y}=\rho_{n}$ and
$\bigl{(}\theta\mapsto
S^{\Theta}(y,\theta)\bigr{)}_{\sharp}p_{\Theta|Y=y}=\rho_{p}$ for any
$y\in\operatorname{supp}p_{Y}$. Thus the two component functions of the map
capture the evidence ($Y$-marginal) and the posterior, respectively. Note also
that (68) can easily be extended to depend on the design, that is,

$\displaystyle S(y,\theta;\xi)=\begin{bmatrix}[l]S^{Y}(\xi,y)\\\
S^{\Theta}(\xi,y,\theta)\end{bmatrix}.$ (69)

Suppose now that we obtain a maximum likelihood estimate $\widehat{S}$ of a
map of the form (69), given samples from $p(y,\theta|\xi)$ and some tractable
class of candidate block-triangular maps $\mathcal{S}$, for instance as
described in Baptista et al. (2023):

$\displaystyle\widehat{S}(\,\cdot\,\,;\xi)\in\arg\max_{S\in\mathcal{S}}\sum_{i=1}^{M}\log
S^{\sharp}\rho(y^{(i)},\theta^{(i)}).$ (70)

The map $\widehat{S}$ yields a plug-in estimate of the desired densities, via
its two component functions:

$\displaystyle\hat{q}^{\text{mar}}(y|\xi)$
$\displaystyle=\rho_{n}\bigl{(}\widehat{S}^{Y}(\xi,y)\bigr{)}\det\nabla_{y}\widehat{S}^{Y}(\xi,y),$
$\displaystyle\hat{q}^{\text{pos}}(\theta|y,\xi)$
$\displaystyle=\rho_{p}\bigl{(}\widehat{S}^{\Theta}(\xi,y,\theta)\bigr{)}\det\nabla_{\theta}\widehat{S}^{\Theta}(\xi,y,\theta).$

(For a statistical convergence analysis of transport-based density estimators,
in a general nonparametric setting, see Wang and Marzouk 2022.) In fact, since
the reference $\rho$ is a product distribution, (70) splits into two separate
optimization problems such that the component functions $S^{Y}$ and
$S^{\Theta}$ can be estimated separately, and hence only the component needed
for the desired variational bound needs to be learned. The resulting density
estimates $\hat{q}^{\text{mar}}$ or $\hat{q}^{\text{pos}}$ can then be
substituted into (62) or (63), respectively. The transport approach
effectively defines the approximating class of densities $\mathcal{Q}$ in (64)
or (65) as the set of all densities that can be expressed as $S^{\sharp}\rho$
for $S\in\mathcal{S}$.

Several recent instantiations of this transport approach in OED have
parameterized the maps as normalizing flows; see Kennamer et al. (2023),
Orozco et al. (2024), and Dong et al. (2024). In particular, the lower
component function $S^{\Theta}$ can be represented as a _conditional_
normalizing flow, which is essentially a structured invertible function of
$\theta$ that is parameterized by $y$. Another canonical choice of
$S^{\Theta}$ is a conditional Brenier map (Carlier et al. 2016), which can be
understood as a family of $L^{2}$-optimal transport maps from $p_{\Theta|Y=y}$
to $\rho_{p}$, parameterized by $y$. The essential requirements are to respect
the overall block structure of (68) in an invertible, differentiable map that
pushes forward $p_{Y,\Theta}$ to a block-independent reference distribution.

Figure 3.2, adapted from Li et al. (2024a), shows an application of transport-
based density estimators to the nonlinear Mössbauer spectroscopy example
described in Feng and Marzouk (2019, Section 4.2), for a fixed design $\xi$.
The orange violin plots and circles show repeated independent estimates of the
upper bound (62), while the blue violin plots and circles illustrate repeated
independent estimates of the lower bound (63). These estimates are produced by
learning the appropriate transport map component, $S^{Y}$ or $S^{\Theta}$, via
the adaptive semi-parametric procedure described in Baptista et al. (2023),
which naturally enlarges the family of maps being considered as the sample
size available for estimation increases. Here the maps $S$ are strictly
triangular, and hence approximations of the Knothe–Rosenblatt rearrangement;
see Rosenblatt (1952), Knothe (1957) and Santambrogio (2015, Section 2.3). The
horizontal axis shows the total number of independent samples drawn from the
model, comprising two batches: one to estimate the map and a separate batch to
estimate the outer expectations $\mathbb{E}_{Y,\Theta|\xi}$ in (62) and (63).
(Asymptotically optimal allocations of this sample budget are discussed in Li
et al. (2024a).) The upper and lower bounds ‘sandwich’ the true EIG as the
sample size increases.

It is useful to understand the source and nature of the randomness in these
results. Each estimate of a transport map yields a density estimate $\hat{q}$,
and plugging any such estimate into (62) and (63) yields a guaranteed upper or
lower bound, as appropriate. Approximating the expectations in (62) and (63)
with samples, however, yields an unbiased estimate of the bound. There is no
guarantee that each realization of the estimator will be above or below the
true EIG, though the figure suggests that no realizations cross this threshold
at the sample sizes considered here. Overall, however, fluctuation in the
estimates reflects randomness in both $\hat{q}$ and in the Monte Carlo
approximation of the outer expectation. The red dashed line in Figure 3.2 is
an NMC estimate of the EIG obtained with $1.58\times 10^{6}$ samples, which is
two to four orders of magnitude greater than the number of samples used in the
variational estimators. The NMC estimate appears to be positively biased, even
at this large sample size.

Figure 3.2: Variational upper (orange) and lower (blue) bounds on the EIG in a
nonlinear design problem, compared to a biased estimate obtained via NMC
(dashed red line). See the discussion in Section 3.2. Figure adapted from Li
et al. (2024a).

As noted earlier, estimating the posterior density to construct a lower bound
for EIG is well suited to implicit models, since (63) does not require the
ability to evaluate $p(y|\theta,\xi)$. But the flexibility of using such
density estimates also applies to other information-based design criteria
discussed in Section 2.2.2. For instance, the expected marginal information
gain objective in (28) can be bounded from below by estimating
$p(\theta_{1}|y,\xi)$, i.e., by replacing this marginal posterior density with
some $q(\theta_{1}|y,\xi)$. Baptista et al. (2024a) uses this approach to
assess the information value of different observables in a phase-field model
(a nonlinear PDE), where high-dimensional random initial conditions take the
role of the ‘nuisance parameters’ $\Theta_{2}$. If the prior marginal density
$p(\theta_{1})$ is not available in closed form, it can also be estimated from
samples, with the caveat that if $p(\theta_{1}|y,\xi)$ and $p(\theta_{1})$ are
both replaced with approximations, the resulting EIG approximation will be
neither an upper nor a lower bound for the true marginal EIG in general
(Foster et al. 2019).

In the case of goal-oriented OED (31)–(32), the need for learning both the
numerator and denominator in the associated density ratios is more acute: for
a generic $\Psi$, none of the densities in (31) or (32) may be readily
available. But as long as we can simulate from $p(z,y|\xi)$—which can be done
by drawing $\theta^{(i)}\sim p(\theta)$, drawing $y^{(i)}\sim
p(y|\theta^{(i)},\xi)$, and evaluating $z^{(i)}=\Psi(\theta^{(i)})$—the
density estimation techniques discussed above will apply. If we choose to use
transport in this setting, then two block-triangular maps will be relevant:

$\displaystyle S_{1}(y,z;\xi)=\begin{bmatrix}[l]S_{1}^{Y}(\xi,y)\\\\[6.0pt]
S_{1}^{Z}(\xi,y,z)\end{bmatrix}\quad\text{and}\quad
S_{2}(z,y;\xi)=\begin{bmatrix}[l]S_{2}^{Z}(\xi,z)\\\\[6.0pt]
S_{2}^{Y}(\xi,z,y)\end{bmatrix}.$ (71)

Approximating the densities in (31) would use the lower component of $S_{1}$
and the upper component of $S_{2}$; conversely, approximating the densities in
(32) would use the upper component of $S_{1}$ and the lower component of
$S_{2}$.

In closing, we note that transport approximations of the densities relevant to
an EIG calculation can be built by means other than maximum likelihood
estimation. For instance, Koval et al. (2024) first prescribe some density
over the continuous real-valued design variables $\xi$, and then build
triangular transport maps directly from unnormalized evaluations of the joint
density $p(y,\theta|\xi)p(\xi)$, using a functional tensor-train format (Cui
et al. 2023). This construction method does not correspond to tightening a
variational bound, in contrast to maximizing the likelihood, but the error in
EIG due to the map approximation can nonetheless be analyzed.

Finally we note that another use for the variational approximations discussed
in this subsection—specifically, the approximation
$q^{\text{pos}}(\theta|y,\xi)$—is as a biasing distribution in NMC estimation
of EIG. This idea, proposed in Foster et al. (2019), fits directly into (57)
by setting $q^{i,\xi}(\theta)=q^{\text{pos}}(\theta|y^{(i)},\xi)$ and recovers
the consistency guarantees of NMC.

###  3.3 More general variational bounds for mutual information

So far, we have discussed bounds on mutual information that follow from
parameterizing a density—or a transport map from which a density estimate can
be derived. These bounds are variational in the sense that they can be
tightened by solving an optimization problem over some class of densities or
transport maps. For the bounds discussed so far, if the class of densities or
transports is sufficiently rich, in principle there is a member of the class
that corresponds to the true density and the bound will attain the true mutual
information.

But there are many other variational bounds for mutual information, based on
different approximation formats and more general classes of functions (i.e.,
not just densities or invertible transport maps). Such bounds are in fact an
important topic in information theory, with many applications in machine
learning. And several have found recent use in OED. Poole et al. (2019)
provide a unifying framework for understanding many variational bounds, which
inspires some of the discussion below.

To provide some context, we first recall a classical nonparametric estimator
of mutual information (which is not variational). The $k$-nearest-neighbor
estimator of Kraskov et al. (2004) (KSG) takes joint sample pairs
$(y^{(i)},\theta^{(i)})\sim p_{Y,\Theta}$ as input, and uses the statistics of
nearest-neighbor distances in both the joint space and in the marginal
directions to construct an estimator. Its construction is related to
$k$-nearest-neighbor estimators of the Shannon entropy (Kozachenko and
Leonenko 1987). The KSG estimator has been used for OED (Terejanu et al.
2012), but is known to scale poorly with dimension. Rigorous statistical
analyses of this estimator, e.g., proofs of consistency for fixed $k$ and
bounds on the rate of convergence of the MSE with sample size $N$, have
appeared relatively recently (Gao et al. 2018). For instance, Gao et al.
(2018, Corollary 2) establish an upper bound on the MSE that is, up to
polylogarithmic factors, $\mathcal{O}(N^{-2/(n+p)})$, where $n$ is the
dimension of $Y$ and $p$ is the dimension of $X$. They suggest that this
convergence rate cannot be refined even if the densities of interest are
assumed to be Hölder-smooth.

Nonparametric estimators of the entropy (rather than the mutual information)
have been used to construct lower bounds for mutual information in Ao and Li
(2020). Here the idea is to write
$\mathcal{I}(Y;\Theta|\xi)=H(Y|\xi)-H(Y|\Theta,\xi)$ and then to seek an upper
bound for the conditional entropy $H(Y|\Theta,\xi)$. The authors develop an
upper bound that requires only (unconditional) entropy estimation, performed
via a $k$-nearest-neighbors approach (Kozachenko and Leonenko 1987).

####  3.3.1 Variational lower bounds and approximate density ratios

Now we turn our discussion back to variational bounds, parameterized by
learnable functions. Dropping our earlier requirement that the function be a
properly normalized density or a smooth invertible transport map, Poole et al.
(2019) develop the following bound, which holds for any function
$f:\mathbb{R}^{p}\times\mathbb{R}^{n}\to\mathbb{R}$ and any positive function
$a:\mathbb{R}^{n}\to\mathbb{R}_{>0}$:

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{\Theta,Y}\left[f(\Theta,Y)\right]-\mathbb{E}_{Y}\left[\frac{\mathbb{E}_{\Theta}\left[e^{f(\Theta,Y)}\right]}{a(Y)}+\log
a(Y)-1\right]$ $\displaystyle\eqqcolon\mathcal{L}^{\text{TUBA}}(f,a)\,.$ (72)

Here ‘TUBA’ denotes ‘tractable unnormalized Barber–Agakov (BA),’ for reasons
we will explain shortly. The function $f$ is known as the ‘critic.’ This bound
is tightened by maximizing simultaneously over $f$ and $a$. Since the mutual
information will generally depend on the design $\xi$, the critic and the
function $a$ can both depend on the design as well, but here we have
temporarily dropped dependence on $\xi$ to lighten notation.444For the
remainder of this section, we will keep notation lighter by not explicitly
writing dependence of the mutual information, the critic, the likelihood, and
other quantities on $\xi$, with the understanding that everything is
appropriately conditioned on $\xi$—since we are considering the mutual
information at a particular design. We can link this bound to the BA bound in
(63), which uses variational approximation over a class of normalized
densities, by writing

$\displaystyle q(\theta|y)=\frac{p(\theta)}{Z(y)}e^{f(\theta,y)}$ (73)

and

$\displaystyle Z(y)=\mathbb{E}_{\Theta}[\exp f(\Theta,y)].$ (74)

Hence $f$ parameterizes an approximation ${q(\theta|y)}/{p(\theta)}$ of the
true posterior-to-prior density ratio ${p(\theta|y)}/{p(\theta)}$, and $Z$ is
the $y$-dependent normalizing constant. The bound
$\mathcal{L}^{\text{TUBA}}(f,a)$ is then tight when

$f(\theta,y)=\log p(y|\theta)+c(y)\quad\text{and}\quad a(y)=Z(y),$

where $c$ is any function of $y$. One specific case is when $c=0$, and hence
$f(\theta,y)=\log p(y|\theta)$ and $a(y)=p(y)$.

Poole et al. (2019) point out that the lower bound for mutual information
introduced earlier by Nguyen et al. (2010), often called the NWJ bound, is a
special case of (72) that follows from setting $a=e$:

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{\Theta,Y}\left[f(\Theta,Y)\right]-\frac{1}{e}\,\mathbb{E}_{Y}\left[Z(Y)\right]\eqqcolon\mathcal{L}^{\text{NWJ}}(f).$
(75)

In this case, the unique optimal critic is

$f^{\ast}(\theta,y)=1+\log\frac{p(\theta|y)}{p(\theta)}.$

Unbiased estimates of these lower bounds, for any $f$ and $a$, can be
constructed by simple Monte Carlo estimation of the expectations therein.
Similarly, one can produce unbiased estimates of the gradient of each bound
with respect to $f$ and $a$ (or some parameterization thereof). The modern
approach is to parameterize $f$ and $a$ with deep neural networks, and to
maximize over the network parameters, with the hope that the resulting class
of functions is rich enough to yield a good approximation of the optimal
critic. The statistical properties of these nonlinear M-estimators are not yet
fully understood, though there are some recent results. For instance,
Sreekumar and Goldfeld (2022) develop non-asymptotic error bounds for neural
estimators of statistical divergences, using shallow neural network
parameterizations of $f$ (i.e., with a single hidden layer) and under strong
smoothness assumptions, capturing the impact of sample size and neural network
size. Statistical convergence results have also been developed for the NWJ
estimator in the more classical setting where the critic is constrained to lie
in a reproducing kernel Hilbert space and parameterized with kernels (Nguyen
et al. 2010).

A more general observation is that while these variational bounds can in
principle become tight, their estimators exhibit high variance—especially when
the true value of the mutual information is large; see McAllester and Stratos
(2020) for a discussion.

For completeness, we should note that the classical Donsker–Varadhan bound
(Donsker and Varadhan 1983), a variational lower bound on KL divergence, can
be adapted to mutual information to yield

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{\Theta,Y}\left[f(\Theta,Y)\right]-\log\mathbb{E}_{Y}\left[Z(Y)\right]\eqqcolon\mathcal{L}^{\text{DV}}(f).$
(76)

When the critic $f$ is parameterized by a neural network, this is the ‘mutual
information neural estimation’ (MINE) approach of Belghazi et al. (2018). The
challenge in using this bound is that the second term on the right yields a
nested expectation, following the definition of $Z$ in (74). Hence, in
practice, when estimating the bound with Monte Carlo as in Belghazi et al.
(2018), we revert to the NMC setting of Section 3.1, where estimates are
biased at finite inner loop sample sizes. Then, as emphasized in Poole et al.
(2019), we will obtain (to leading order) a positively biased estimator of a
lower bound. The mean of this estimator is neither an upper nor a lower bound
on the desired $\mathcal{I}(Y;\Theta)$.

####  3.3.2 Multi-sample lower bounds

A different class of lower bounds are called ‘multi-sample’ in that they
involve not only an expectation over the joint distribution of $\Theta,Y\sim
p_{\Theta,Y}$ but simultaneously another expectation over
$\Theta_{2:M}\sim\prod_{j=2}^{M}r_{\Theta}$, where $r_{\Theta}$ could be
different from the marginal $p_{\Theta}$.

The simplest among these is the so-called ‘prior contrastive estimator’ (PCE)
proposed in Foster et al. (2020):

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{Y_{1}|\Theta_{1}}\mathbb{E}_{\Theta_{1}}\mathbb{E}_{\Theta_{2:M}}\left[\log\frac{p(Y_{1}|\Theta_{1})}{\frac{1}{M}\sum_{j=1}^{M}p(Y_{1}|\Theta_{j})}\right]\eqqcolon\mathcal{L}^{\text{PCE}}(M),$
(77)

where the expectation is over $\Theta_{1},Y_{1}\sim p_{\Theta,Y}$ and
$\Theta_{2:M}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}p_{\Theta}$. Keep in
mind that to estimate the right-hand side of (77), one must draw many
independent and identically distributed (i.i.d.) realizations of
$(Y_{1},\Theta_{1},\Theta_{2:M})$. Suppose that we draw $N$ such (outer-loop)
realizations, with sample index $i$. The result is very much like the NMC
estimator (53), with one crucial modification: each outer loop sample of the
parameters, $\Theta_{1}^{(i)}$, is used once more, within the inner-loop
estimate of the evidence for the corresponding $Y_{1}^{(i)}$. Doing so
guarantees that $\mathcal{L}^{\text{PCE}}(M)$ is a lower bound of the mutual
information for any $M>1$; moreover, like NMC, the bound becomes tight as
$M\to\infty$ (Foster et al. 2020, Theorem 1). Intuitively, recycling the
outer-loop parameter sample makes the denominator of (77) over-estimate the
evidence and hence under-estimate the mutual information. The samples
$\Theta_{2:M}^{(i)}$ are called contrastive samples, in that, unlike the
original sample $\Theta_{1}^{(i)}$, they are not responsible for
$Y_{1}^{(i)}$.

Just as the variance and bias of the NMC estimators in Section 3.1 were
improved by importance sampling, a $y$-dependent, properly normalized, biasing
distribution $q(\theta|y)$ can be used within a contrastive estimator as well.
(The biasing distribution will generally depend on $\xi$ too, but recall that
we are not explicitly writing dependence on $\xi$ in this section, to keep
notation simple.) This leads to the idea of an ‘adaptive contrastive
estimator’ (ACE) proposed by Foster et al. (2020),

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{Y_{1},\Theta_{1}}\mathbb{E}_{\Theta_{2:M}}\left[\log\frac{p(Y_{1}|\Theta_{1})}{\frac{1}{M}\sum_{j=1}^{M}\frac{p(Y_{1}|\Theta_{j})p(\Theta_{j})}{q(\Theta|Y_{1})}}\right]\eqqcolon\mathcal{L}^{\text{ACE}}(M,q),$
(78)

where now the expectation is over $\Theta_{1},Y_{1}\sim p_{\Theta,Y}$ and
$\Theta_{2:M}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}q_{\Theta|Y_{1}}$. The
optimal biasing distribution in this case is $q(\theta|y)=p(\theta|y)$, which
leads to zero-variance estimates of the evidence $p(Y_{1})$. Seeking this $q$
is exactly the idea described at the end of Section 3.2—combining a
variational approximation of the posterior density with nested Monte
Carlo—with only the addition of a contrastive recycling of the outer-loop
$\Theta_{1}$. With the optimal choice of $q$, this bound is tight for any
finite $M$; otherwise, for sub-optimal $q$, it becomes tight as $M\to\infty$.

Both (77) and (78) require the ability to evaluate the likelihood, and thus in
contrast with the variational bounds (63), (72), and (75), they are not suited
to implicit models. A further modification can make (77) likelihood-free. It
introduces a critic function $f$ as follows:

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{Y_{1},\Theta_{1}}\mathbb{E}_{\Theta_{2:M}}\left[\log\frac{\exp
f(\Theta_{1},Y_{1})}{\frac{1}{M}\sum_{j=1}^{M}\exp
f(\Theta_{j},Y_{1})}\right]\eqqcolon\mathcal{L}^{\text{NCE}}(M,f),$ (79)

where the expectation is again over $\Theta_{1},Y_{1}\sim p_{\Theta,Y}$ and
$\Theta_{2:M}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}p_{\Theta}$. This is
essentially the ‘information noise-contrastive estimation’ (InfoNCE) approach
proposed in van den Oord et al. (2018). $\mathcal{L}^{\text{NCE}}(M,f)$ is a
lower bound of the mutual information for all $f$ and $M$ (Poole et al. 2019).
An amalgam of (78) and (79), incorporating both a biasing distribution $q$ and
a critic $f$, is called ‘likelihood-free ACE’ in Foster et al. (2020). We note
that both $\mathcal{L}^{\text{PCE}}(M)$ and $\mathcal{L}^{\text{NCE}}(M,f)$
are bounded above by $\log(M)$—no matter how the critic is chosen in the
latter. Hence if $\mathcal{I}(Y;\Theta)>\log(M)$, these contrastive bounds
cannot become tight; caution is therefore needed when the mutual information
is high or the contrastive sample size $M$ is small. Note also that the
$\log(M)$ limit does not apply to (78) and to other variants that use
importance sampling, because in principle there exists a biasing distribution
that yields a zero-variance estimate of the evidence, and in such a case even
$M=1$ is sufficient.

An alternative way of estimating the PCE and InfoNCE bounds given only a
single set of i.i.d. sample pairs $(Y_{i},\Theta_{i})_{i=1}^{M}\sim
p_{Y,\Theta}$ is to rotate the role of the ‘data-generating sample’ and the
‘contrastive samples’ through the set. We can take the $i$th pair as the
‘data-generating sample’ and use the remaining $M-1$ samples as ‘contrastive
samples’ to obtain

$\log\frac{\exp f(\Theta_{i},Y_{i})}{\frac{1}{M}\sum_{j=1}^{M}\exp
f(\Theta_{j},Y_{i})}.$

Repeating this for $i=1,\ldots,M$ so that each pair has the opportunity to be
the ‘data-generating sample,’ and then taking the average, we obtain the
estimator,

$\displaystyle\frac{1}{M}\sum_{i=1}^{M}\log\frac{\exp
f(\Theta_{i},Y_{i})}{\frac{1}{M}\sum_{j=1}^{M}\exp f(\Theta_{j},Y_{i})}.$ (80)

As shown in Poole et al. (2019), the expectation of this estimator is a lower
bound of the true mutual information, for any critic $f$:

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{(Y_{i},\Theta_{i})_{i=1}^{M}}\left[\frac{1}{M}\sum_{i=1}^{M}\log\frac{\exp
f(\Theta_{i},Y_{i})}{\frac{1}{M}\sum_{j=1}^{M}\exp
f(\Theta_{j},Y_{i})}\right].$ (81)

The expectation is also bounded above by $\log M$, in the same way as (77) and
(79) above. Inequality (81) is in fact the form of InfoNCE bound originally
proposed in van den Oord et al. (2018). It is in essence the same as (79) but
using correlated, rather than independent, terms to estimate the expectation
of

$\log\Bigl{(}\exp f(\Theta_{1},Y_{1})\Big{/}\frac{1}{M}\sum_{j=1}^{M}\exp
f(\Theta_{j},Y_{1})\Bigr{)}.$

If the conditional density $p(y|\theta)$ is tractable, then we can replace $f$
with an optimal critic, $f(\theta,y)=\log p(y|\theta)$, such that (81) becomes

$\displaystyle\mathcal{I}(Y;\Theta)$
$\displaystyle\geq\mathbb{E}_{(Y_{i},\Theta_{i})_{i=1}^{M}}\left[\frac{1}{M}\sum_{i=1}^{M}\log\frac{p(Y_{i}|\Theta_{i})}{\frac{1}{M}\sum_{j=1}^{M}p(Y_{i}|\Theta_{j})}\right],$
(82)

which can be understood as a ‘sample re-use’ analogue of the PCE bound (77)
above. Interestingly, the corresponding estimator (i.e., (80) with the log-
likelihood substituting for $f$) was used in Huan and Marzouk (2013), before
it was formally understood to be a lower bound in expectation. There, the
rationale was to lower the cost of nested Monte Carlo by reducing the number
of expensive model evaluations from $\mathcal{O}(NM)$ to $\mathcal{O}(M)$, in
a setting where evaluations of the likelihood $p(y|\theta)$ are costly for
each new value of $\theta$, but cheap for new values of $y$ given a fixed
$\theta$.

####  3.3.3 Deploying variational bounds

The relative merits of the variational bounds discussed above are very much a
current subject of research, involving both empirical comparisons and some
theory (Czyż et al. 2023, Song and Ermon 2020, Letizia et al. 2023). In
general, the TUBA, NWJ, and DV bounds can become tight with increasingly good
choices of critic, but their Monte Carlo estimates generally have much higher
variance than estimates of multi-sample bounds such as PCE and InfoNCE. The
latter bounds, on the other hand, are hampered by the fact they cannot become
tight for arbitrarily large mutual information when using finite contrastive
sample sizes. We also emphasize that new variants of these variational bounds
continue to be proposed; we have attempted to cover the key constructions
relevant to OED, but do not claim to be comprehensive. We refer to Song and
Ermon (2020) and Czyż et al. (2023) for empirical evaluations of the mutual
information estimators (and estimators of mutual information bounds) discussed
here, along with many others. Important settings for these evaluations include
low versus high values of mutual information, heavy versus light tails,
various forms of non-Gaussianity, high dimensionality, sparse interactions,
and so forth. As emphasized by Czyż et al. (2023), evaluations limited to the
Gaussian setting, while convenient in that exact analytical results are
available for comparison, are not fully representative.

Deploying the bounds discussed in Sections 3.2 and 3.3 requires optimization.
Specifically, it is useful to have an _unbiased estimate_ of the gradient of
each bound with respect to parameters representing the critic $f$ (and
possibly $a$) or the density approximation $q$, along with an unbiased
estimate of the gradient of the bound with respect to the design $\xi$, for
any given critic or density approximation. Explicit forms of these gradient
estimates, for several bounds above, are detailed in Kleinegesse and Gutmann
(2021). Given these gradients, one can _simultaneously_ maximize the bound
with respect to both the variational parameters and the design, e.g., via
stochastic gradient ascent. We will return to this idea in Section 4.3.

###  3.4 Low-dimensional structure

So far, we have discussed a variety of methods that, essentially, approximate
the densities or density ratios appearing in expressions for the mutual
information (expected information gain). A different and complementary class
of approximation methods involve _dimension reduction_. Here we will discuss
methods for Bayesian OED that exploit possible low dimensionality in the
_update_ from prior to posterior; such methods rest on the idea that the
posterior is well-approximated by a distribution that departs from the prior
only on a subspace of dimension $r\ll p$. At the same time, we will consider
dimension reduction for the observation vector $y\in\mathbb{R}^{n}$; the
central idea is that conditioning on $y$ could be replaced by conditioning on
the _projection_ of $y$ onto a subspace of dimension $s\ll n$.

Both of these notions of low dimensionality appear quite often in Bayesian
inverse problems, where the data are related to the parameters $\theta$ by the
action of a forward operator that is somehow smoothing, e.g., in that it
suppresses information about finer scales (Cui et al. 2014). Low
dimensionality more generally results from limited informativeness of the data
in certain directions of $\bm{\Theta}$ and/or redundancy among the components
of $y$. There is additional potential for dimension reduction in _goal-
oriented_ problems, as in Bayesian $\text{D}_{\text{A}}$-optimality and
(31)–(32), where the experimenter wishes to learn about a specific quantity of
interest $Z=\Psi(\Theta)$. Identifying and exploiting this low dimensionality
can lead to much more computationally efficient approximations of the mutual
information, as we shall discuss below.

####  3.4.1 Linear models

We first consider dimension reduction, of both parameters and data, in the
context of finite-dimensional Bayesian linear-Gaussian models. We begin with
design for parameter inference and then consider goal-oriented design.

##### Optimal low-rank approximations for parameter inference.

Recall the setup of the Bayesian linear-Gaussian model from the start of
Section 2.1.2. The mutual information between parameters and observations in
this setting can already be written in closed form as a ratio of log-
determinants, as in (25) and (26). It is nonetheless useful to express this
quantity in a way that reveals the structure and _intrinsic dimensionality_ of
the linear-Gaussian model.

The central objects in such a construction are generalized eigenvalue problems
involving combinations of the prior covariance matrix $\Gamma_{\Theta}$; the
noise covariance $\Gamma_{Y|\Theta}$; the posterior covariance
$\Gamma_{\Theta|Y}$ (8); the negative Hessian of the log-likelihood,
$H\coloneqq G^{\top}\Gamma_{Y|\Theta}^{-1}G$; and the marginal covariance of
the data, $\Gamma_{Y}=G\Gamma_{\Theta}G^{\top}+\Gamma_{Y|\Theta}$. Several
different generalized eigenvalue problems involving these matrices can be
written, all closely related via simple transformations. Perhaps the most
direct is the symmetric definite generalized eigenproblem

$\displaystyle Hw_{j}=\lambda_{j}\Gamma_{\Theta}^{-1}w_{j}\,,$ (83)

where, since $H\succeq 0$ and $\Gamma^{-1}_{\Theta}\succ 0$, we have
$\lambda_{j}\geq 0$. This eigenproblem can be understood as balancing
informativeness of the likelihood with informativeness of the prior, via the
Rayleigh ratio $w^{\top}Hw/w^{\top}\Gamma_{\Theta}^{-1}w$: an eigendirection
$w_{j}\in\mathbb{R}^{p}$ associated with $\lambda_{j}>1$ is more constrained
by the likelihood than by the prior, and thus more important to the update
from prior to posterior. See Spantini et al. (2015) for thorough
interpretations and demonstrations of this eigenstructure. Closely related to
(83) are two other $p\times p$ parameter-space matrix pencils:
$(\Gamma_{\Theta}-\Gamma_{\Theta|Y},\Gamma_{\Theta|Y})$ and
$(\Gamma_{\Theta},\Gamma_{\Theta|Y})$. The first has eigenvalue-eigenvector
pairs $(\lambda_{j},v_{j})$, where $v_{j}=\Gamma_{\Theta}^{-1}w_{j}$, and the
second has eigenvalue-eigenvector pairs $(1+\lambda_{j},v_{j})$. See Jagalur-
Mohan and Marzouk (2021, Proposition 10) for a proof of this equivalence. The
leading eigendirections $v_{j}$ (ordered from largest to smallest
$\lambda_{j}$) are the directions along which the ratio of posterior to prior
variance is minimized; see Spantini et al. (2015, Section 3.1).

In the linear-Gaussian setting, there is a duality between spectral structure
in the parameter space and in the data space. Consider the $n\times n$ data-
space matrix pencils $(\Gamma_{Y}-\Gamma_{Y|\Theta},\Gamma_{Y|\Theta})$ and
$(\Gamma_{Y},\Gamma_{Y|\Theta})$. Jagalur-Mohan and Marzouk (2021, Proposition
10) show that $(\Gamma_{Y}-\Gamma_{Y|\Theta},\Gamma_{Y|\Theta})$ and
$(\Gamma_{\Theta}-\Gamma_{\Theta|Y},\Gamma_{\Theta|Y})$ have identical non-
trivial generalized eigenvalues, $\lambda_{j}>0$, and that
$(\Gamma_{Y},\Gamma_{Y|\Theta})$ and $(\Gamma_{\Theta},\Gamma_{\Theta|Y})$
have identical eigenvalues that are strictly greater than one, i.e.,
$1+\lambda_{j}>1$. The eigenvectors $u_{j}$ of the two data-space pencils are
identical, with $w_{j}=\frac{1}{\alpha_{j}}\Gamma_{\Theta}G^{\top}u_{j}$ for
all $j$ with $\lambda_{j}>0$, where $\alpha_{j}$ is a scaling parameter that
can be set to $\alpha_{j}=\sqrt{\lambda_{j}}$ to obtain $\langle
U_{i},U_{j}\rangle_{\Gamma_{Y|\Theta}}=\delta_{ij}$, given $\langle
W_{i},W_{j}\rangle_{\Gamma_{\Theta}^{-1}}=\delta_{ij}$.

A convenient way of simultaneously capturing spectral structure in the
parameter space and the data space is through the ‘whitened’ forward operator
$\Gamma_{Y|\Theta}^{-1/2}G\Gamma_{\Theta}^{1/2}$. The singular value
decomposition of this matrix is

$\Gamma_{Y|\Theta}^{-1/2}G\Gamma_{\Theta}^{1/2}=\sum_{j\geq
1}\sqrt{\lambda_{j}}a_{j}b_{j}^{\top},$

where the left singular vectors are related to the preceding data-space
eigendirections via $a_{j}=\Gamma_{Y|\Theta}^{1/2}u_{j}$ and the right
singular vectors similarly are $b_{j}=\Gamma_{\Theta}^{-1/2}w_{j}$. See also
Spantini et al. (2015, Remark 4). The whitened forward operator thus
summarizes much of the essential structure of the linear-Gaussian Bayesian
inverse problem and appears in many contexts. In the case of uncorrelated
observational errors $\Gamma_{Y|\Theta}=\sigma^{2}I$, for example, the adjoint
of the whitened forward operator, $\sigma^{-1}\Gamma_{\Theta}^{1/2}G^{\top}$,
is the central object in column subset selection approaches to finding
D-optimal designs (Eswar et al. 2024). We will comment further on these
methods in Section 4.

Spantini et al. (2015, Theorem 2.3) consider optimal approximations of the
posterior covariance $\Gamma_{\Theta|Y}$ that take the form of low-rank
updates of the prior covariance. Optimality is cast in terms of minimizing a
class of loss functions $\ell(M)$ over the manifold of all symmetric positive
definite (SPD) matrices $M$; this class includes the geodesic distance from
$M$ to $\Gamma_{\Theta|Y}$ on SPD manifold, and the KL divergence or Hellinger
distance from $\mathcal{N}(\mu,M)$ to $\mathcal{N}(\mu,\Gamma_{\Theta|Y})$,
where $\mu\in\mathbb{R}^{p}$ is an arbitrary (common) mean vector. The
approximations sought are of the form
$\widetilde{\Gamma}^{r}_{\Theta|Y}\in\\{\Gamma_{\Theta}-KK^{T}\succ 0,\
\text{rank}(K)\leq r\\}$, and the optimal approximation for any $r\leq n$,
simultaneously minimizing all loss functions in the class, is given by

$KK^{\top}=\sum_{i=1}^{r}\frac{\lambda_{i}}{1+\lambda_{i}}w_{i}w_{i}^{\top};$

hence it follows from the leading eigenpairs of the parameter-space pencil
$(H,\Gamma_{\Theta}^{-1})$. The decay of the eigenvalue sequence
$(\lambda_{j})_{j\geq 1}$, common to all of the matrix pencils above, captures
the intrinsic dimensionality of the Bayesian model.

These optimal approximation results are related to the approximation of a more
central object of interest in OED: the mutual information. Using the fact that
the generalized eigenvectors $v_{j}$ simultaneously diagonalize
$\Gamma_{\Theta}$ and $\Gamma_{\Theta|Y}$, or that the generalized
eigenvectors $u_{j}$ simultaneously diagonalize $\Gamma_{Y}$ and
$\Gamma_{Y|\Theta}$, the mutual information between parameters $\Theta$ and
data $Y$ (25)–(26) can be written as

$\displaystyle\mathcal{I}(Y;\Theta)=\frac{1}{2}\sum_{j}\log(1+\lambda_{j})\,.$
(84)

This expression for mutual information has appeared in many places in recent
literature, e.g., Alexanderian et al. (2016a) and Giraldi et al. (2018). We
note also that $\mathcal{I}(Y;\Theta)$ can be written using the squared
canonical correlation scores (Bach and Jordan 2002) between $Y$ and $\Theta$;
these scores follow from generalized eigenvalue problems that are slightly
different from those discussed above, but closely related. Truncating the sum
in (84) after $r<\min(n,p)$ terms yields a low-rank approximation of the
mutual information. From a computational perspective, since the MI is
dominated by the largest generalized eigenvalues (though the $\log$ function
slows this decay), problems with quickly decaying spectra and hence low
intrinsic dimension are easier to approximate: one must compute only the
leading eigenvalues of any of the matrix pencils described above.

Truncating (84) in this way in fact corresponds to the mutual information
obtained from an optimal $r$-dimensional projection of the data $Y$, of the
parameters $\Theta$, or of both simultaneously. Specifically, Giraldi et al.
(2018) show that the leading eigenvectors $U_{1:r}=[u_{1}\cdots u_{r}]$ of the
data-space matrix pencils define low-dimensional projections of the data that
are optimal at any given dimension $r\leq n$, in the sense of maximizing
$\mathcal{I}(A_{r}^{\top}Y;\Theta)$ over matrices $A\in\mathbb{R}^{n\times
r}$. The resulting mutual information is

$\mathcal{I}(U_{1:r}^{\top}Y;\Theta)=\frac{1}{2}\sum_{j=1}^{r}\log(1+\lambda_{j}).$

Similarly, replacing $\Theta$ with the $r$-dimensional projection
$V_{1:r}^{\top}\Theta$ yields

$\mathcal{I}(Y;V_{1:r}^{\top}\Theta)=\frac{1}{2}\sum_{j=1}^{r}\log(1+\lambda_{j}).$

Going further, one can show that the same, optimal value of mutual information
is also achieved by

$\mathcal{I}(U_{1:r}^{\top}Y;V_{1:r}^{\top}\Theta)=\frac{1}{2}\sum_{j=1}^{r}\log(1+\lambda_{j}).$

Note that this truncated mutual information is equivalent to the mutual
information that would be obtained if the forward operator $G$ in (4) were
replaced by certain ‘projected’ versions. Define

$\mathbb{P}_{\text{obs}}\coloneqq\Gamma_{Y|\Theta}U_{1:r}U_{1:r}^{\top}\quad\text{and}\quad\mathbb{P}_{\text{param}}\coloneqq\Gamma_{\Theta}V_{1:r}V_{1:r}^{\top}=W_{1:r}W_{1:r}^{\top}\Gamma_{\Theta}^{-1}.$

Then one can show that

$\displaystyle\mathbb{P}_{\text{obs}}G=G\mathbb{P}_{\text{param}}=\mathbb{P}_{\text{obs}}G\mathbb{P}_{\text{param}},$
(85)

and that these projected forward operators achieve the mutual information

$\frac{1}{2}\sum_{j=1}^{r}\log(1+\lambda_{j}).$

We emphasize, however, that the equivalence of these three projected models,
and more broadly the strict duality between parameter and observation
reduction discussed above, is specific to the linear-Gaussian setting.

The expression (84) naturally appears in infinite-dimensional formulations of
linear Bayesian inverse problems, as discussed in Section 2.3.1. Indeed,
finding the leading eigenpairs of (83) is equivalent to constructing a low-
rank approximation of the prior-preconditioned Hessian operator
$\widetilde{H}$, as discussed in Alexanderian et al. (2016a). To see this
intuitively in finite dimensions, note that
$\widetilde{H}=\Gamma_{\Theta}^{1/2}H\Gamma_{\Theta}^{1/2}$, and that its
(simple, not generalized) eigenvalue-eigenvector pairs are
$(\lambda_{j},\Gamma_{\Theta}^{-1/2}w_{j})$, where $(\lambda_{j},w_{j})$ are
the eigenpairs of $(H,\Gamma_{\Theta}^{-1})$. The algorithms proposed in
Alexanderian and Saibaba (2018) (see also Ghattas and Willcox 2021) compute
these low-rank approximations efficiently in a discretization-invariant
manner, when the forward operator $G$ is described by partial differential
equations.

##### Optimal low-rank approximations for goal-oriented design.

Now we turn to the _goal-oriented_ linear setting. Our interest here lies not
in learning the parameters $\Theta$ per se, but rather in informing a quantity
of interest $Z=A^{\top}\Theta$, where $A\in\mathbb{R}^{p\times s}$ has full
column rank $s<p$. Maximizing the mutual information $\mathcal{I}(Y;Z)$ is
equivalent to maximizing the Bayesian $\text{D}_{\text{A}}$-optimality
criterion mentioned in Section 2.1.2, and is a linear special case of
(31)–(32). Now, however, there is further potential for dimension reduction,
stemming from the interaction of $A$ with other elements of the problem.
Intuitively, the update from the prior predictive distribution to posterior
predictive distribution of $Z$ (and hence the mutual information) should be
well captured by directions in the parameter space $\mathbb{R}^{p}$ that are
simultaneously informed by the data (relative to the prior) _and_ influential
on $Z$. We can expect the latter consideration to help ‘screen out’ parameter
directions that are dampened by $A^{\top}$. (Consider, as a limiting example,
elements of the kernel of $A^{\top}$.)

This idea is formalized by Spantini et al. (2017). Beginning with the linear-
Gaussian model $Y=G\Theta+\mathcal{E}$ (4), a Gaussian prior distribution
$\Theta\sim\mathcal{N}(0,\Gamma_{\Theta})$ (zero mean for simplicity), and the
specification of $Z=A^{\top}\Theta$, Spantini et al. (2017, Lemma 2.2)
introduce an _equivalent_ linear-Gaussian model that directly relates $Y$ and
$Z$,

$\displaystyle Y=GA_{\dagger}Z+\Delta,$ (86)

where $A_{\dagger}\coloneqq\Gamma_{\Theta}A\Gamma_{Z}^{-1}$,
$\Gamma_{Z}\coloneqq A^{\top}\Gamma_{\Theta}A$,
$Z\sim\mathcal{N}(0,\Gamma_{Z})$ is the prior induced on $Z$, and

$\displaystyle\Delta\sim\mathcal{N}(0,\Gamma_{\Delta})\quad\text{with}\quad\Gamma_{\Delta}\coloneqq\Gamma_{Y|\Theta}+G(\Gamma_{\Theta}-\Gamma_{\Theta}A\Gamma_{Z}^{-1}A^{\top}\Gamma_{\Theta})G^{\top}.$

Crucially, these choices render the ‘effective noise’ $\Delta$ _independent_
of $Z$. Now the low-rank approximation methods discussed earlier for the
$Y$–$\Theta$ system can be applied directly to the $Y$–$Z$ system above.

For example, as a goal-oriented counterpart to (83), one can compute the
leading eigenpairs of the matrix pencil
$(A^{\top}_{\dagger}G^{\top}\Gamma_{\Delta}^{-1}GA_{\dagger},\Gamma_{Z}^{-1})$.
Letting $(\tilde{\lambda}_{j})_{j\geq 1}$ denote the generalized eigenvalues
of this pencil, arranged in descending order, we can then write the mutual
information as
$\mathcal{I}(Y;Z)=\frac{1}{2}\sum_{j}\log(1+\tilde{\lambda}_{j})$ and truncate
this sum as desired to obtain a low-rank approximation; see Wu et al. (2023a)
for more details and a first application of this approach in OED. As
demonstrated in Spantini et al. (2017), $(\tilde{\lambda}_{j})$ can decay much
more quickly than the spectrum $(\lambda_{j})$ of (83), yielding a more
efficient approximation than would be obtained by first approximating the
posterior covariance matrix of $\Theta$ and then applying $A$. Moreover,
analogously to the discussion above, one can obtain optimal _goal-oriented_
projections of the data using the leading eigenvectors of
$(\Gamma_{Y},\Gamma_{\Delta})$. Each of these eigenproblems optimally balances
all the ingredients of the goal-oriented inference problem: the structure of
the forward operator, the scale and structure of the observational noise
covariance and prior covariance, and the ultimate prediction goals.

In closing, we should point out that (as in Section 3.3) we have not
explicitly noted dependence of the quantities above on the design $\xi$.
Recall, however, that our generic linear-Gaussian model (4) can be written
more explicitly as
$Y|\Theta,\xi\sim\mathcal{N}(G(\xi)\Theta,\Gamma_{Y|\Theta}(\xi))$. Hence all
of the eigenvalues, eigenspaces, and projectors described above may depend on
$\xi$ via $G$ and $\Gamma_{Y|\Theta}$. The details of this dependence are
problem-dependent, but it is of interest to understand how to exploit any
smoothness or other regularity that may be present—for instance, by
differentiating the eigenvalues $\lambda_{j}$ with respect to $\xi$. To our
knowledge, algorithmic approaches in this vein are in their infancy.

####  3.4.2 Nonlinear models

Extending the low-dimensional structure discussed in Section 3.4.1 to
nonlinear models is considerably more complex, and a subject of ongoing
research. The tools we will highlight here originate in dimension reduction
for Bayesian inverse problems and Bayesian statistical models more generally:
likelihood-informed subspaces (Cui et al. 2014, Cui and Tong 2022), certified
dimension reduction (Zahm et al. 2022, Li et al. 2024b), and related efforts.
Our perspective is similar to that taken in the linear case: find the subspace
of the parameter space where the posterior differs ‘most’ from the prior,
which is equivalent to finding the subspace that is best informed by the data;
and find the subspace of the data that is most informative of the parameters.

We focus here on the approximation of mutual information, highlighted in
Baptista et al. (2022), which makes the intuitive notions just stated more
precise. In the nonlinear case, it is generally difficult to find _optimal_
low-dimensional projections. Instead, one can develop gradient-based upper
bounds on the error induced by projecting $Y$ and $\Theta$ onto lower-
dimension subspaces, i.e., find upper bounds on the difference

$\mathcal{I}(Y;\Theta)-\mathcal{I}(U_{1:s}^{\top}Y;V_{1:r}^{\top}\Theta)$

for some matrices $U_{1:s}\in\mathbb{R}^{n\times s}$ and
$V_{1:r}\in\mathbb{R}^{p\times r}$, and then minimize these bounds over $U$
and $V$. (Note that the matrices $U$ and $V$ discussed in this section are
generally different from the $U$ and $V$ found in Section 3.4.1.)

The key assumption underlying this analysis is that the joint distribution
$p_{Y,\Theta}$ must satisfy a _subspace logarithmic Sobolev inequality_ ; see
Baptista et al. (2022, Definition 2) or Zahm et al. (2022, Theorem 2.10).
Letting $Z=(Y,\Theta)$ be a random variable taking values in
$\mathbb{R}^{n+p}$, and letting $W\in\mathbb{R}^{(n+p)\times(n+p)}$ be a
unitary matrix, the essence of this assumption is that any conditional
distribution $p(Z_{t}|Z_{\perp}=z_{\perp})$ defined by the decomposition
$W=[W_{t}\ W_{\perp}]$, with $W_{t}\in\mathbb{R}^{(n+p)\times t}$,
$Z_{t}=W_{t}^{\top}Z$, and $Z_{\perp}=W_{\perp}^{\top}Z$, satisfies a
logarithmic Sobolev inequality with constant
$C\left(p(Z_{t}|Z_{\perp}=z_{\perp})\right)$ bounded above by some
$\overline{C}<\infty$, for all $W$, $t$, and $z_{\perp}$. As shown in Zahm et
al. (2022), a sufficient condition for a distribution to satisfy the subspace
logarithmic Sobolev inequality is that it has convex support and that its log-
density is a bounded perturbation of a strongly convex function; see Zahm et
al. (2022, Examples 2.5–2.9) and Baptista et al. (2022, Examples 1–2). This
allows, for example, Gaussian mixtures (and hence certain multi-modal
distributions with strictly positive densities) but not distributions with
exponential tails.

Now define the diagnostic matrices

$\displaystyle H_{\Theta}$ $\displaystyle=\iint(\nabla_{\theta}\nabla_{y}\log
p_{Y|\Theta}(y|\theta))^{\top}(\nabla_{\theta}\nabla_{y}\log
p_{Y|\Theta}(y|\theta))p_{Y,\Theta}(y,\theta)\,\mathrm{d}\theta\,\mathrm{d}y,$
(87) $\displaystyle H_{Y}$ $\displaystyle=\iint(\nabla_{\theta}\nabla_{y}\log
p_{Y|\Theta}(y|\theta))(\nabla_{\theta}\nabla_{y}\log
p_{Y|\Theta}(y|\theta))^{\top}p_{Y,\Theta}(y,\theta)\,\mathrm{d}\theta\,\mathrm{d}y,$
(88)

where $\nabla_{\theta}\nabla_{y}\log p_{Y|\Theta}\in\mathbb{R}^{n\times p}$.
We then have the following theorem, adapted from Baptista et al. (2022).

###### Theorem 3.1

Let $p_{Y,\Theta}$ satisfy a subspace logarithmic Sobolev inequality with
constant $\overline{C}<\infty$. Then, for any unitary matrices $V=[V_{1:r}\
V_{r+1:p}]\in\mathbb{R}^{p\times p}$ and $U=[U_{1:s}\
U_{s+1:n}]\in\mathbb{R}^{n\times n}$, we have

$\displaystyle\mathcal{I}(Y;\Theta)-\mathcal{I}(U_{1:s}^{\top}Y;V_{1:r}^{\top}\Theta)$
$\displaystyle\leq\overline{C}^{2}\Bigl{(}\text{tr}\bigl{(}V_{r+1:p}^{\top}H_{\Theta}V_{r+1:p}\bigr{)}+\text{tr}\bigl{(}U_{s+1:n}^{\top}H_{Y}U_{s+1:n}\bigr{)}\Bigr{)}\,.$
(89)

Crucially, the right-hand side of (89) can be minimized for any dimensions
$r\leq p$ and $s\leq n$ by choosing $V_{1:r}$ to be the leading eigenvectors
of $H_{\Theta}$ and $U_{1:s}$ to be the leading eigenvectors of $H_{Y}$. This
choice yields the bound

$\displaystyle\mathcal{I}(Y;\Theta)-\mathcal{I}(U_{1:s}^{\top}Y;V_{1:r}^{\top}\Theta)$
$\displaystyle\leq\overline{C}^{2}\left(\sum_{i=r+1}^{p}\lambda_{i}(H_{\theta})+\sum_{i=s+1}^{n}\lambda_{i}(H_{Y})\right).$
(90)

Fast decay of the eigenvalues $\lambda_{i}$ of the two diagnostic matrices
thus provides more opportunity for dimension reduction, with controlled error.

Computationally, this dimension reduction approach can be viewed as a first
step towards approximating a mutual information design objective: first, find
low-dimensional projections of the data and parameters,
$Y_{s}=U_{1:s}^{\top}Y$ and $\Theta_{r}=V_{1:r}^{\top}\Theta$, for some $r$
and $s$, by finding the dominant eigenspaces of $H_{Y}$ and $H_{\Theta}$.
(Note that, unlike in the linear-Gaussian case, the spectra of these two
matrices are generally different.) Then, use any of the techniques discussed
in Sections 3.1–3.3 to estimate or bound $\mathcal{I}(Y_{s};\Theta_{r})$. By
the data-processing inequality (or the positivity of the right-hand side of
(89)) we always have $\mathcal{I}(Y_{s};\Theta_{r})\leq\mathcal{I}(Y;\Theta)$.
For example, Li et al. (2024a) use the transport-based density estimators of
Section 3.2 to bound the projected mutual information
$\mathcal{I}(Y_{s};\Theta_{r})$. The density estimation problems to be solved
are thus of reduced dimension, and involve estimating transport maps that act
on lower-dimensional spaces. Alternatively, Wu et al. (2023b) apply $V_{1:r}$
for parameter dimension reduction as proposed above, but use a simpler
dimension reduction scheme for the data (principal component analysis based on
the marginal covariance $\mathbb{C}\text{ov}(Y)$); both reductions are a
prelude to constructing a reduced-dimensional surrogate $\widetilde{G}$ for
the forward operator in a Bayesian inverse problem (see Section 3.1), and then
using this surrogate within a standard nested Monte Carlo estimator of the
mutual information.

In Theorem 3.1, for simplicity, we restricted the matrices $U$ and $V$ to be
unitary. This constraint can effectively be relaxed by preconditioning the
problem. Note that with any change of variables $\overline{\Theta}=A\Theta$
and $\overline{Y}=BY$, the left-hand side of (89) does not change, but the
right-hand side does—via the diagnostic matrices and the log-Sobolev constant.
It is unclear how to choose an optimal change of variables, as the log-Sobolev
constant is in general very hard to compute, but one heuristic suggested in
Baptista et al. (2022, Section 4), for the case of a Gaussian prior and
Gaussian additive noise, is to ‘whiten’ the problem so that
$\Gamma_{\overline{\Theta}}$ and $\Gamma_{\overline{Y}|\overline{\Theta}}$
become identity matrices. In the case of a linear-Gaussian model, this change
of variables leads the proposed method to recover the optimal subspaces of
Section 3.4.1. We also note that the change of variables need not be linear;
any bijective transformation will leave the mutual information terms on the
left-hand side of (89) unchanged. Thus, in principle one could compose a
nonlinear change of variables with linear dimension reduction to effectively
achieve nonlinear dimension reduction!

Another remark is that the bounds in Theorem 3.1 are generally not tight. One
can intuitively see this in the linear-Gaussian setting where, after the
linear preconditioning described above, we have
$\lambda_{i}=\lambda_{i}(H_{\Theta})=\lambda_{i}(H_{Y})$. Then we can compare
(84), where the exact truncation error involves logarithms of trailing
eigenvalues $\log(1+\lambda_{i})$, with (90), which sums the trailing
eigenvalues directly. These quantities are not the same, and the difference is
analyzed in Baptista et al. (2022, Section 4.3); see also Zahm et al. (2022,
Section 2.3). There is reason to believe that recent dimensional improvements
to log-Sobolev inequalities (Eskenazis and Shenfeld 2024) could help tighten
the bounds used here.

##  4 Design optimization methods

Now we will discuss methods for optimizing (generally, maximizing) OED
criteria over possible designs $\xi\in\Xi$. Different representations of
candidate designs and choices of design criterion can yield very different
classes of optimization problems. This section will therefore review a broad
array of optimization approaches, ranging from algorithms for combinatorial
optimization (including continuous relaxations thereof) to gradient-based
algorithms for intrinsically continuous problems. Section 4.1 focuses on exact
designs for linear models, as linearity lends the problem special structure;
here, we discuss various discrete algorithms, as well as convex continuous
relaxations that are followed by rounding. Section 4.2 discusses exact designs
for more general nonlinear problems, with design criteria treated as set
functions. Section 4.3 turns to designs that are parameterized by real
variables (e.g., sensor positions or times, dosages), and discusses how a
variety of derivative-free or gradient-based optimization approaches interact
with estimators or bounds for decision-theoretic design criteria.

###  4.1 Optimization methods for linear design

We state a prototypical linear design problem, by recalling the linear
regression model with features $f:\mathcal{X}\to\mathbb{R}^{p}$ and
uncorrelated observational errors. If the design $\xi$ is supported on $n$
equally weighted points $x_{i}\in\mathcal{X}\in\mathbb{R}^{d}$,
$i=1,\ldots,n$, then the Fisher information matrix of the model is given by
(11); as pointed out earlier, this expression follows from (5) by setting

$G=\left[f(x_{1})^{\top};f(x_{2})^{\top};\ldots;f(x_{n})^{\top}\right]\in\mathbb{R}^{n\times
p}$

and $\Gamma_{Y|\theta}=\sigma^{2}I_{n}$. The classical design objectives
discussed in Section 2.1 can be viewed as functions of the Fisher information
matrix.

We will use the index set $\mathcal{V}\coloneqq\\{1,\ldots,m\\}$ to represent
the support of the design space. It is instructive to think of $\mathcal{V}$
as a specific instance of $\mathcal{X}$ in the setting of a discrete design
domain. In the simplest setting, the task of finding an optimal exact design
amounts to selecting a set $\mathcal{S}\subseteq\mathcal{V}$ of specified
cardinality $n$ that maximizes the chosen objective, i.e., any of the design
criteria for linear models given in Section 2.1. More complex cases can
involve constraints on the overall budget and variable costs associated with
each individual design instance.

Strictly speaking, the set $\mathcal{S}$ is better characterized as a
multiset, where elements can be duplicated. In a conventional set, elements
are typically assumed to be different from one another. The notion of a
multiset is more general, and accommodates the ‘with repetition’ scenario
where a given design point can be chosen multiple times. (Repeated experiments
at the same design point are not necessarily redundant, since they may yield
different outcomes due to statistically independent noise.) The contrasting
scenario is ‘without repetition,’ where each point is selected at most once;
this is applicable, for instance, in the setting of deterministic computer
experiments.

####  4.1.1 Local search and exchange algorithms

Algorithms guided by local search and exchange heuristics are non-sequential
in nature. (We discuss sequential algorithms later in Section 4.1.3.) In the
non-sequential setting, we start with an initial design of the required size,
and at each iteration attempt to improve the quality of the design by
deleting, adding, or exchanging points as guided by certain rules. The non-
sequential nature of these algorithms implies that the resulting solution sets
are non-nested, meaning that optimal design set of cardinality $n$ may look
quite different from the optimal set of cardinality $n-1$ or $n+1$.

In the context of exchange algorithms, the approach of Fedorov (1972) has been
quite influential. Fedorov’s exchange method starts with an arbitrary initial
set $\mathcal{S}_{0}$ of $n$ candidate points, and at each iteration $t$
attempts to improve the objective by exchanging one of the points,

$\mathcal{S}_{t}\leftarrow\\{j\\}\cup\mathcal{S}_{t-1}\setminus\\{i\\},$

where $i\in\mathcal{S}_{t-1}$ and $j\in\mathcal{V}\setminus\mathcal{S}_{t-1}$.
The search at each iteration is exhaustive: the improvement resulting from
every possible exchange is calculated, and the best exchange is selected. The
process is continued as long as an exchange improves the design objective.

Fedorov’s original work focused exclusively on D-optimal design, although its
applicability to other design criteria has also been demonstrated. While the
exchange algorithm is relatively simple, it is accompanied by a sizable
computational overhead due to the exhaustive search at each step. Cook and
Nachtsheim (1980) proposed a modified Fedorov exchange procedure, where they
consider each design point in turn, perhaps in random order, carrying out any
beneficial exchange as soon as it is discovered. Johnson and Nachtsheim (1983)
suggested further improvements by focusing the search on fewer points from the
current design set, specifically those that have the lowest predictive
variance.

The underlying ideas of the exchange algorithm can be extended to the setting
where design points are selected from a continuous compact space, without an
underlying candidate set (Meyer and Nachtsheim 1995). A preliminary step here
is to specify the expected utility of changing one factor or coordinate at a
time, while holding the others constant. This step can be made more efficient
by using Gaussian process emulators, as demonstrated by Overstall and Woods
(2017). These ideas in principle are also applicable for nonlinear design
problems.

Although the original Fedorov exchange algorithm was proposed several decades
ago, rigorous analyses establishing approximation guarantees have only
recently emerged. Here, an approximation guarantee is a provable bound on the
ratio between the value of the objective function evaluated at the solution
returned by the algorithm, and the optimal solution. It is thus a bound on the
worst-case performance of the algorithm (in a relative sense). This is a
typical way to assess the performance of an algorithm that produces
approximate solutions to optimization problems, in particular NP-hard problems
(Hochba 1997, Vazirani 2001).

Madan et al. (2019) establish approximation guarantees for local search
algorithms for D-optimal design and A-optimal design; their results for the
case of A-optimality are restricted to the less general ‘with repetition’
setting. They show that the algorithms are asymptotically optimal when $n/p$
is large, and that in the case of A-optimal design there could be arbitrarily
bad local optima. They also propose approximate local search algorithms, where
exchanges are made only when the objective improves by a factor of $1+\delta$,
leading to improved running times with a slight degradation in the
approximation guarantees. Lau and Zhou (2022) improve upon these results using
novel analysis tools, and provide approximation guarantees for D/A/E-optimal
designs. In particular they show that Fedorov’s exchange method for A-optimal
design works well as long as there exists an near-optimal solution with a
well-conditioned design matrix.

####  4.1.2 Continuous convex relaxation approaches

The combinatorial linear design problem is strictly speaking an integer
program, whose exact solution is often intractable. One elegant solution
approach is to perform a continuous relaxation of the design variables into a
convex program that can be solved using established techniques (Boyd and
Vandenberghe 2004, Ben-Tal and Nemirovski 2001). A solution to the original
combinatorial problem is then obtained by appropriately rounding the convex
solution to an integer solution. Various polynomial-time (and possibly
randomized) rounding algorithms, which convert the convex solution to an
integer solution for the combinatorial problem, have been proposed.

Consider the linear design problem, where $\phi$ denotes an operator that acts
on the Fisher information matrix to form the design criterion (e.g.,
$\phi=-\text{tr}$ for A-optimality, $\phi=-\log\det$ for D-optimality). The
relaxed continuous optimization problem takes the form

$\displaystyle\operatornamewithlimits{argmin}_{\xi=(\xi_{1},\ldots,\xi_{n})}\phi\left(\sum_{i=1}^{m}\xi_{i}f(x_{i})f(x_{i})^{\top}\right)\quad\text{subject
to}\quad\lVert\xi\rVert_{1}\leq n.$ (91)

Written in this way, the traditional optimality criteria discussed in Section
2.1.1 are all convex functions of the Fisher information matrix. The convex
constraint $\lVert\xi\rVert_{1}\leq n$ ensures that only $n$ points are
selected. If we further specify $0\leq\xi_{i}\leq 1$, then we constrain each
point to be selected at most once, corresponding to the without-replacement
scenario. More complex constraints, e.g., encoding varied costs associated
with selecting each point and a ceiling on the cumulative cost (i.e., a
knapsack constraint), can easily be added to (91). The fractional solution
weights from this convex program have natural analogues to the notion of
designs as probability measures, specifically to continuous designs, as
discussed in Section 2.1.3. However the notion of a continuous design arose
organically from the work of early experimental design researchers,
independently of these modern continuous relaxation formulations.

There are a number of convex programming relaxation approaches for the
D/A/E-optimal design problems, differing from each other in how the rounding
algorithm provides integer solutions. Even though these design criteria share
attributes such as convexity, they behave differently in terms of
approximability. Allen-Zhu et al. (2017) connect experimental design to matrix
sparsification (Spielman and Srivastava 2008, Batson et al. 2009) and use
regret minimization methods (Allen-Zhu et al. 2015) to obtain approximate
designs for popular optimality criteria. Singh and Xie (2020) devise an
approximation algorithm for D-optimal design where the rounding strategy uses
approximate positively correlated distributions. Nikolov et al. (2022) develop
an approximation algorithm for D/A-optimal design using proportional volume
sampling; interestingly, they also show that the same approach will not work
for E-optimal design. Lau and Zhou (2022) modify the iterative randomized
rounding algorithm based on the regret minimization framework of Lau and Zhou
(2020) for D/A/E-optimal design problems with knapsack constraints.
Approximation algorithms using more involved relaxation techniques have
enabled D-optimal design under partition (Nikolov and Singh 2016) and matroid
(Madan et al. 2020) constraints.

####  4.1.3 Sequential algorithms for finding optimal designs

Sequential algorithms arrive at a solution set either by the gradual addition
of candidate points to a smaller design set, or by the sequential deletion of
candidate points from a larger design set. These are called ‘forward’ and
‘backward’ procedures, respectively (Atkinson et al. 2007). Typically the
greedy heuristic guides the selection of candidate points at each step. In
some cases, designs obtained using sequential procedures can be further
improved using the non-sequential techniques outlined in Section 4.1.1. For
the case of D/A-optimal designs, Madan et al. (2019) analyze the sequential
forward procedure, and provide approximation guarantees that retain a
specificity to the initialized set. The dependence on the initial set is
rather undesirable, which can be ameliorated by leveraging distributed
computing resources and running the algorithms with different initializations.
Sequential design algorithms and greedy heuristics have been more generally
studied using the language of set functions for combinatorial problems; we
will discuss these methods more deeply in Section 4.2.

###  4.2 Combinatorial approaches for discrete design variables

In this section we formulate experimental design problems as optimization of
set functions. The latter topic has been studied extensively (Wolsey and
Nemhauser 1999, Lovász 2007, Schrijver 2003, Papadimitriou and Steiglitz
1998), and has a rich mathematical structure. We discuss some fundamental
properties of set functions and their relevance to popular design criteria,
including criteria for nonlinear design. We also discuss algorithms for their
optimization.

####  4.2.1 Background

Previously we defined the candidate index set as
$\mathcal{V}\coloneqq\\{1,\ldots,m\\},\ m\in\mathbb{Z}_{>0}$. Let its power
set (i.e., the set of all subsets) be denoted by $2^{\mathcal{V}}$. An
important property of a certain class of set functions is submodularity: any
real-valued set function $\psi:2^{\mathcal{V}}\rightarrow\mathbb{R}$ such that
$\psi(\emptyset)=0$ is submodular (Fujishige 2005, Bach 2013) if and only if,
for all subsets $\mathcal{A},\mathcal{B}\subseteq\mathcal{V}$, we have

$\psi(\mathcal{A})+\psi(\mathcal{B})\geq\psi\left(\mathcal{A}\cup\mathcal{B}\right)+\psi\left(\mathcal{A}\cap\mathcal{B}\right).$
(92)

A function is _supermodular_ if its negation is submodular, and it is
_modular_ if it is both supermodular and submodular. Many objective functions
that arise in experimental design are naturally submodular. Consider, for
instance, sensor placement where each sensor has a certain coverage area, and
our interest is in optimizing the collective locations of a fixed number of
sensors. There is a natural _diminishing returns_ property that accompanies
such objectives, which is characteristic of submodular functions. An
alternative but equivalent (Fujishige 2005, Bach 2013) definition of
submodularity using first-order differences highlights this _diminishing
returns_ property and is often easier to demonstrate in practice: the set
function $\psi$ is submodular if and only if, for all
$\mathcal{A},\mathcal{B}\subseteq\mathcal{V}$ and $\nu\in\mathcal{V}$ such
that $\mathcal{A}\subseteq\mathcal{B}$ and $\nu\notin\mathcal{B}$, we have

$\psi(\mathcal{A}\cup\\{\nu\\})-\psi(\mathcal{A})\geq\psi(\mathcal{B}\cup\\{\nu\\})-\psi(\mathcal{B}).$

The term
$\rho_{\nu}(\mathcal{A})\coloneqq\psi(\mathcal{A}\cup\\{\nu\\})-\psi(\mathcal{A})$
is the incremental gain associated with the element $\nu$ when added to the
set $\mathcal{A}$. It refers to the amount of change in the objective when the
individual item $\nu$ is added to an existing pool of items in set
$\mathcal{A}$. The definition is applicable by extension when we add a set
$\mathcal{B}$ to set $\mathcal{A}$.

In certain cases when the function is not strictly submodular, it is helpful
to understand how much it deviates from submodularity, or in other
circumstances to quantify how close a set function is to being modular. Such
characterizations also prove useful in analyzing the performance of many
optimization algorithms. The notion of curvature introduced by Conforti and
Cornuéjols (1984) in the context of non-negative set-functions is a bound on
the intrinsic value of an item against its value in conjunction with all items
in the candidate set. Formally, curvature is defined as the scalar $c$:

$c\coloneqq\max_{\nu\in\mathcal{V}}\frac{\rho_{\nu}(\emptyset)-\rho_{\nu}(\mathcal{V}\setminus\\{\nu\\})}{\rho_{\nu}(\emptyset)}=1-\min_{\nu\in\mathcal{V}}\frac{\rho_{\nu}(\mathcal{V}\setminus\\{\nu\\})}{\rho_{\nu}(\emptyset)}.$

The classical notion of curvature (Conforti and Cornuéjols 1984) measures how
close a submodular set function is to being modular, while the notion of
generalized curvature (Bian et al. 2017) measures how close a set function—not
necessarily submodular—is to being supermodular. The submodularity ratio
introduced by Das and Kempe (2011) for a general non-negative set function is
a lower bound on the ratio of the sum of incremental gains associated with
elements in a set, to the incremental gain associated with the set itself.
Intuitively it captures ‘how close’ to submodularity is the function in
question. Formally it is defined as the scalar $\gamma$:

$\gamma\coloneqq\min_{\mathcal{B}\subseteq\mathcal{V},\mathcal{A}\cap\mathcal{B}=\emptyset}\frac{\sum_{\nu\in\mathcal{A}}\rho_{\nu}(\mathcal{B})}{\rho_{\mathcal{A}}(\mathcal{B})}.$

In Bian et al. (2017) the submodularity ratio and the generalized curvature
together quantify how close a set function is to being modular.

The counterpart of the submodularity ratio termed the supermodularity ratio
(Tzoumas et al. 2021, Bogunovic et al. 2018, Karaca and Kamgarpour 2018,
Jagalur-Mohan and Marzouk 2021) has also proved useful in analyzing many
algorithms. It is defined as the scalar $\eta$:

$\eta\coloneqq\min_{\mathcal{B}\subseteq\mathcal{V},\mathcal{A}\cap\mathcal{B}=\emptyset}\frac{\rho_{\mathcal{A}}(\mathcal{B})}{\sum_{\nu\in\mathcal{A}}\rho_{\nu}(\mathcal{B})}.$

In Jagalur-Mohan and Marzouk (2021), the product of supermodularity and
submodularity ratios measures deviation from modularity. We will revisit these
notions and their implications for performance guarantees in Section 4.2.3. We
now discuss set function attributes of widely used experimental design
criteria.

####  4.2.2 Set function attributes of design criteria

In Section 2.1.1 we considered several classical alphabetic optimality design
criteria. The D-optimality criterion defined using the log determinant is
submodular. Using strictly linear algebraic techniques, it can be shown that
the log determinant of a principal submatrix is a submodular function with
respect to the indices defining the submatrices (Gantmacher and Kreĭn 1960,
Kotelyanskiĭ 1950, Fan 1967, 1968, Kelmans and Kimelfeld 1983, Johnson and
Barrett 1985). Since the log determinant evaluations may be non-positive,
however, most existing algorithms for submodular maximization are not directly
applicable. The A-optimality criterion in general is not submodular, as shown
in Krause et al. (2008). We can, however, qualify its non-submodular nature by
bounding the submodularity ratio $\gamma$ and curvature $c$ for the Bayesian
A-optimal design objective (Bian et al. 2017).

In Section 2.2, we introduced a variety of information theoretic design
criteria. Specifically, we defined the expected information gain (EIG) in
parameters and showed that it is equivalent to the mutual information between
the parameters and the observations given the design,
$\mathcal{I}(Y;\Theta|\xi)$. Suppose that selecting a design corresponds to
choosing individual components of the $\mathbb{R}^{m}$-valued random variable
$Y$. Then we can recast the optimization problem as finding a selection
operator $\mathcal{P}\in\mathbb{R}^{m\times n}$, $n<m$,
$\mathcal{P}\coloneqq[e_{i_{1}},\ldots,e_{i_{n}}]$, where $e_{i_{j}}$ are
distinct canonical unit vectors in $\mathbb{R}^{m}$:

$\mathcal{P}_{\text{opt}}=\operatornamewithlimits{argmax}_{\mathcal{P}\in\mathbb{R}^{m\times
n}}\mathcal{I}\left(\mathcal{P}^{\top}Y;\Theta\right).$ (93)

Given a desired number of observations $n<m$, we seek a selection operator
$\mathcal{P}_{\text{opt}}$ such that the mutual information between the
inference parameter $\Theta$ and the selected observations
$Y_{\mathcal{P}_{\text{opt}}}\coloneqq\mathcal{P}_{\text{opt}}^{\top}Y$ is
maximized. Jagalur-Mohan and Marzouk (2021) show that the mutual information
between parameters $\Theta$ and subsets of data $Y_{\mathcal{P}}$ is a
submodular function if the observations are conditionally independent.
Interestingly, this property holds even when the underlying joint distribution
is non-Gaussian. In the setting of Bayesian inverse problems with additive
noise, correlated noise renders the mutual information design objective non-
submodular. When the inverse problem is linear, it is possible to quantify the
non-submodularity of the information-theoretic objective by bounding the
submodularity ratio $\gamma$ and supermodularity ratio $\eta$. It was shown in
Jagalur-Mohan and Marzouk (2021) that those parameters can be both lower-
bounded by ${\log\zeta_{\text{min}}}/{\log\zeta_{\text{max}}}$, where $\zeta$
is any generalized eigenvalue of the definite pair
$\left(\Gamma_{Y},\Gamma_{Y|\Theta}\right)$.

####  4.2.3 Greedy algorithms for cardinality-constrained designs

The case of cardinality-constrained optimization commonly arises in OED, with
the prototypical example being sensor placement. When the design objective is
monotone and non-negative, the greedy heuristic of successively picking the
candidate corresponding to the highest incremental gain performs well despite
its simplicity. Nemhauser et al. (1978) were the first to analyze the greedy
heuristic for the class of submodular functions, and showed that the algorithm
has a constant factor $1-1/e$ approximation guarantee. Nemhauser and Wolsey
(1978) further showed that the approximation guarantee cannot be improved in
general by any other polynomial-time algorithm.

In the above paragraph, the term ‘constant factor’ refers to the approximation
guarantee not depending on the particular instance of the function and only
requiring the function to be submodular. By incorporating parameters that
capture more specific attributes of the function, however, more expressive
results can be obtained and offer useful insights into performance of the
greedy heuristic in different scenarios. For instance, Conforti and Cornuéjols
(1984) proved the more refined guarantee $\frac{1}{c}(1-e^{-c})$ for
submodular functions, where $c\in[0,1]$ is the curvature. When the function is
known to have a small curvature, the improved performance of the greedy
heuristic is easily explained. Now suppose the function is not submodular but
has a submodularity ratio $\gamma\in[0,1]$; Das and Kempe (2011) proved that
the greedy heuristic has a $1-e^{-\gamma}$ approximation guarantee. Bian et
al. (2017) improved upon that result by incorporating the notion of a
generalized curvature $\alpha\in[0,1]$ to obtain the factor
$\frac{1}{\alpha}(1-e^{-\alpha\gamma})$. As we alluded to in Section 4.2.2,
these parameters can be concretely bounded for many experimental design
objectives.

The greedy algorithm has a $\mathcal{O}(mn)$ complexity, where $m$ is the size
of the candidate set and $n$ is the desired cardinality. Closely related
variants of the greedy heuristic can be better choices depending on the
context and needs. Robertazzi and Schwartz (1989) explored an accelerated
version wherein the computed incremental gains are stored and exploited in the
successive step, possibly reducing the overall number of function evaluations.
This is much like the modified Fedorov algorithm we discussed previously in
Section 4.1.1. To reduce the complexity further and make it independent of the
cardinality constraint, Mirzasoleiman et al. (2015) analyze a randomized
version of the greedy heuristic, termed stochastic greedy. This algorithm
achieves, in expectation, a ($1-1/e-\epsilon$) approximation guarantee
relative to the optimum solution. The number of function evaluations does not
depend on the cardinality constraint, but linearly on the size of the
candidate set, thus reducing the complexity substantially.

To benefit from modern HPC platforms, Mirzasoleiman et al. (2013) proposed a
two-stage parallelized version which reduces the number of function
evaluations per parallel process. The approximation guarantee for the
algorithm, however, in general depends on the size of the candidate set and
cardinality constraint, which can only be overcome in special cases. Jagalur-
Mohan and Marzouk (2021) analyzed a batch version of stochastic and
distributed greedy algorithms with applicability to non-submodular objectives.
The heuristic involved choosing multiple candidates in each step but relying
solely on the incremental gains associated with individual candidates. This
reduces the computational overhead but with reduced approximation guarantees.

####  4.2.4 Beyond greedy: algorithms for more general design problems

As mentioned earlier in passing, many experimental design problems involve
constraints more complex than simple cardinality bounds. Consider, for
instance, sensor placement with non-uniform costs $c(\nu)$ associated with
placing the sensors. Here the goal could be to maximize the design objective
$\psi(\cdot)$ while ensuring that the cumulative cost is within a budget
$\sum_{\nu\in\mathcal{S}}c(\nu)\leq b$:

$\max_{\mathcal{S}\subset\mathcal{V}}\psi(\mathcal{S})\quad\text{subject
to}\quad\sum_{\nu\in\mathcal{S}}c(\nu)\leq b.$

This is the well known knapsack problem, where the goal is to maximize the set
function $\psi$ subject to a non-negative modular constraint. A natural
modification of the the greedy algorithm in this setting is to pick the
element maximizing the benefit-to-cost ratio, at each step. Surprisingly, for
the case when $\psi$ is submodular, it can be shown that either the output of
the standard greedy algorithm or the index set returned by the cost-benefit
greedy algorithm will be within ${(1-1/e)}/{2}$ of the optimal solution
(Leskovec et al. 2007). A stronger result ($1-1/e$) is possible via a more
computationally involved algorithm that exhaustively enumerates all subsets of
size 3, and augments them using the cost-benefit greedy algorithm (Sviridenko
2004).

Although we exclusively focused on non-negative monotone functions in Section
4.2.3, several design criteria can be non-monotone. In Section 2.3.2 we
discussed an information-theoretic design objective arising in sensor
placement that is symmetric submodular but strictly speaking non-monotone:
$\mathcal{I}\left(\Theta(x_{s});\Theta(x_{s^{c}})\right)$. For the
maximization of such submodular objectives with cardinality constraints,
Buchbinder et al. (2014) have analyzed discrete random greedy and continuous
double greedy algorithms. We refer the interested readers to this work for
more details on such methods.

An important tool in the field of submodular optimization is the use of
extensions, particularly the multilinear extension (Vondrák 2008) in the
context of maximization. The multilinear extension of the set function $\psi$
is the function $\Psi:[0,1]^{m}\rightarrow\mathbb{R}$ defined as

$\Psi(\mathbf{\xi})=\sum_{\mathcal{S}\subseteq\mathcal{V}}\psi(\mathcal{S})\prod_{i\in\mathcal{S}}\xi_{i}\prod_{i\in\mathcal{V}\setminus\mathcal{S}}(1-\xi_{i}).$

An intuitive interpretation of this extension is to think of the original set
function as defined over the corners of the hypercube $\\{0,1\\}^{m}$, while
the extension is valid over the entire unit cube $[0,1]^{m}$. $\Psi(\xi)$ is
the expected value of $\psi$ over sets, where for any set $\mathcal{S}$, each
element $i$ is included independently with probability $\xi_{i}$, and not
included with probability $1-\xi_{i}$. If we write
$\mathcal{S}\sim\mathbf{\xi}$ to indicate that $\mathcal{S}$ is the random
subset sampled according to $\mathbf{\xi}$, then the multilinear extension is
simply
$\Psi(\mathbf{\xi})=\mathbb{E}_{\mathcal{S}\sim\mathbf{\xi}}[\psi(\mathcal{S})]$.
The multilinear extension (Vondrák 2008) is quite different from the Lovász
extension (Lovász 1983); the latter maps any discrete submodular function to
its continuous convex counterpart while the former maps it to its continuous
concave counterpart. Note that these continuous extensions are quite different
from the convex relaxation formulations we described in Section 4.1.2! Many of
the more general approaches to maximizing submodular functions under a wide
class of constraints rely on the multilinear extension. The approach here
involves first approximately solving the problem

$\max\Psi(\mathbf{\xi})\quad\text{subject
to}\quad\mathbf{\xi}\in\mathcal{X}\subseteq[0,1]^{m}$

and then rounding the continuous solution to obtain a near-optimal set
(Vondrák 2010, Chekuri et al. 2014, Vondrák et al. 2011, Calinescu et al.
2011, Sviridenko et al. 2017). We refer the interested reader to the cited
works for details of the algorithms. For a broader survey on submodular
maximization, see Krause and Golovin (2014).

###  4.3 Optimizing real-valued design variables

In this section we consider OED problems where the design is naturally
parameterized by real-valued coordinates. These coordinates could be, for
example, the spatial locations of a finite collection of sensors, times at
which measurements should be taken, the values of certain experimental
controls (pressure, temperature), and so on. In Section 2.1.3 we formalized
this continuous-parameter case by letting these coordinates be elements of a
set $\mathcal{X}\subseteq\mathbb{R}^{d}$ and considering the design $\xi$ to
be a probability measure on $\mathcal{X}$. An exact $n$-point design can then
be understood as a mixture of Dirac measures,
$\xi=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$, where each
$x_{i}\in\mathcal{X}$. Here we will focus on continuous optimization methods
for such exact designs, and assume the number of support points $n$ to be
prescribed. The design measure $\xi$ is thus entirely parameterized by
$(x_{1},\ldots,x_{n})\in\bigtimes_{i=1}^{n}\mathcal{X}\subseteq\mathbb{R}^{D}$,
where $D=nd$. To keep notation intuitive in this setting, we will elide the
notions of the design measure and its parameterization and simply write
$\xi=(x_{1},\ldots,x_{n})$. In other words, throughout this section we will
consider $\xi$ to be a vector in
$\mathcal{X}^{n}=\bigtimes_{i=1}^{n}\mathcal{X}=\Xi$. This way, the design
remains exact, and gradients with respect to $\xi$ are straightforward to
define.

Of course, one could discretize $\mathcal{X}^{n}$ and employ approaches from
Section 4.2 to find a design within this discretized space. This approach may
be convenient when dealing with complicated constraints on the set of feasible
designs or when the number of desired experiments $n$ is not fixed (as in the
case of knapsack constraints). For any naïve discretization of
$\mathcal{X}^{n}$, however, the number of candidate designs would grow
exponentially with $nd$. Instead, we will focus here on optimization
algorithms that make use of a continuous (Euclidean) parameterization of
$\xi$.

When the design criterion can be evaluated exactly (i.e., in closed form) for
any given $\xi$, the problem of optimizing over a real-valued $\xi$ in some
sense becomes standard—though of course it inherits the natural difficulties
associated with how the parameterized design enters the objective, or with the
geometry of the feasible domain $\mathcal{X}$ (e.g., nonlinear constraints,
non-convexity). These issues are generally quite problem dependent, however,
and are not specific to OED.

The decision-theoretic OED objectives formulated in Section 2.2, on the other
hand, do raise certain cross-cutting challenges for optimization. The expected
utility $U$, which is the objective function to be maximized, typically must
be estimated using Monte Carlo techniques. The optimal design problem can thus
be written as

$\displaystyle\xi^{\ast}\in\arg\max_{\xi\in\Xi}\mathbb{E}_{W|\xi}\big{[}\widehat{U}(\xi,W)\big{]},$
(94)

where $\widehat{U}$ is an estimator of $U$, and the random variable $W$ is the
source of stochasticity in this estimator. For example, the Monte Carlo
estimator of a general expected utility $U$ given in (52),

$\displaystyle\widehat{U}(\xi,W)=\frac{1}{N}\sum_{i=1}^{N}u(\xi,Y^{(i)},\Theta^{(i)}),$
(95)

with $\Theta^{(i)}\sim p(\theta)$ and $Y^{(i)}\sim p(y|\theta^{(i)},\xi)$, has
$W=(Y^{(i)},\Theta^{(i)})_{i=1}^{N}$. If, at any given design $\xi$, all we
can compute are estimates $\widehat{U}(\xi,W)$ for different realizations of
$W\sim p(W|\xi)$, the optimization problem (94) is essentially a _stochastic
approximation_ (SA) problem.

Implicit in writing (94) is the assumption that we are content to maximize
$\mathbb{E}[\widehat{U}]$: that is, either $\widehat{U}(\xi,W)$ is an unbiased
estimator of the desired utility $U(\xi)$ at any $\xi$ or the bias
$\mathbb{E}[\widehat{U}]-U$ is small and acceptable to the experimenter. In
Sections 4.3.3 and 4.3.4 we will discuss additional techniques that can be
deployed when this is not the case, as in the case of mutual information
maximization with biased nested Monte Carlo estimators or variational bounds.
Starting below, however, we will use the notation
$\bar{U}(\xi)\coloneqq\mathbb{E}_{W|\xi}[\widehat{U}(\xi,W)]$ for the
objective in (94) to make the possible presence of bias clear.

Many nonlinear optimization methods can be applied to (94), and we do not
attempt to survey this vast literature here. Instead, we will briefly recall a
number of derivative-free and gradient-based approaches that have been used in
the context of OED with continuous design variables. We call a method
‘derivative-free’ if it only requires evaluations of $\widehat{U}$ (even if
these evaluations are used to estimate derivatives) and we call a method
‘gradient-based’ if it requires additional derivative information as an input,
e.g., unbiased estimates of the gradient of $\bar{U}$ with respect to $\xi$.
In both cases, we assume that the objective $\bar{U}$ is differentiable, but
we make no assumptions on convexity or other structure. We also note that
estimating $\nabla_{\xi}\bar{U}$ may, in many cases, require evaluating
derivatives of the log-likelihood or an underlying simulation model with
respect to $\xi$, and that this task may not be straightforward for PDE-based
forward operators or intractable likelihoods.

####  4.3.1 Derivative-free methods

Larson et al. (2019) provide a comprehensive survey of derivative-free
optimization methods, focusing on local optimization, with methods for solving
stochastic problems of the form (94) specifically discussed in Larson et al.
(2019, Section 6). A key desideratum for algorithms for this setting is that
they perform well with noisy objective evaluations. Below we briefly point out
several classes of derivative-free optimization methods that do exactly this,
and that are therefore useful for OED. We make no attempt to be comprehensive,
and instead refer the reader to the preceding survey for more information.

##### Bayesian optimization.

We introduced Bayesian optimization (BO) in Section 2.4 to elucidate its
differences from OED for Gaussian processes. But BO can be quite useful as a
tool within OED—specifically, as a means of solving (94) when $\xi$ is of
moderate dimension. BO (Močkus 1975, Jones et al. 1998, Wang et al. 2023) is
essentially a derivative-free algorithm for global optimization, well suited
to ‘black-box’ objective functions that are expensive to evaluate and
potentially noisy. BO uses Gaussian process regression to construct and refine
a model for the objective (in the case of OED, the function $\bar{U}$); this
regression naturally handles noisy pointwise evaluations and smooths the
underlying function estimate. Both the level of noise and the smoothness of
the reproducing kernel Hilbert space (RKHS) containing this estimate (the
posterior mean of the Gaussian process) can be adjusted by learning
hyperparameters of the covariance function of the Gaussian process. Examples
of the application of BO to OED include, among others, works by Weaver et al.
(2016), Overstall and Woods (2017), Xu and Liao (2020), and Zhong et al.
(2024). For more information on BO, we refer to the surveys by Shahriari et
al. (2016) and Frazier (2018).

##### Stochastic approximation methods based on finite differences.

A prototypical deterministic optimization method is gradient ascent:

$\displaystyle\xi_{k+1}=\xi_{k}+\alpha_{k}g(\xi_{k}),$ (96)

where $g(\xi_{k})\coloneqq\nabla_{\xi}\bar{U}(\xi_{k})$ is the gradient of the
objective $\bar{U}(\xi)=\mathbb{E}_{W|\xi}[\widehat{U}(\xi,W)]$ evaluated at
$\xi_{k}$, and $\\{\alpha_{k}\\}$ is a sequence of scalar step sizes with
$\alpha_{k}>0$. The Kiefer–Wolfowitz algorithm (Kiefer and Wolfowitz 1952)
estimates $g(\xi_{k})$ by applying centered differences to unbiased estimates
of $\bar{U}$. Specifically, for the $i$th component of $g$, where
$i=1,\ldots,D$, we have

$\displaystyle
g_{i}(\xi_{k})\approx\frac{\widehat{U}(\xi_{k}+e_{i}\Delta_{k},w_{i}^{+})-\widehat{U}(\xi_{k}-e_{i}\Delta_{k},w_{i}^{-})}{2\Delta_{k}},$
(97)

where $w_{i}^{+}$ and $w_{i}^{-}$ are independent draws from $p_{W|\xi_{k}}$,
$e_{i}$ is the unit vector in coordinate $i$, and the scalar $\Delta_{k}>0$ is
a difference parameter. Hence we need to evaluate the estimator $\widehat{U}$
$2D$ times to compute each gradient estimate. Under appropriate conditions on
the step size sequence $\\{\alpha_{k}\\}$ and difference sequence
$\\{\Delta_{k}\\}$, the objective function $\bar{U}$ and the noise $W$,
Kiefer–Wolfowitz iterations converge to a first-order critical point of
$\bar{U}$ almost surely; see, e.g., Blum (1954) and Bhatnagar et al. (2013).
Simultaneous perturbation stochastic approximation (SPSA) (Spall 1998b, a) is
similar in form to Kiefer–Wolfowitz, except that it always uses a single
centered difference at each iteration, the direction of which is chosen from a
particular probability distribution. SPSA thus uses only _two_ independent
realizations of $\widehat{U}$ at each step $k$, independent of the dimension
of $\xi$. An intuitive justification for SPSA is that error in the gradient
produced by restricting the direction of the finite differences ‘averages out’
over a large number of iterations Spall (1998a).

##### Model-based methods.

Model-based methods for derivative-free optimization use pointwise evaluations
of the objective function (and constraints, if applicable) to form local
approximations (called models) that can be analyzed to produce the next
optimization step. There is an extensive literature on these methods, with
many effective algorithms; for thorough reviews, see Larson et al. (2019) and
Conn et al. (2009). Local models can take the form of low-order polynomials or
radial basis function approximations, constructed via interpolation or
regression on carefully chosen point sets. Updating of these models is often
set within a trust region framework that carefully manages the point sets, the
trust region radius, and the acceptance of each optimization step. Rigorous
convergence guarantees have been developed for most prevalent algorithms
(Larson et al. 2019).

Model-based derivative-free optimization has been applied in settings where
only noisy estimates $\widehat{U}$ of the objective are available. The models
fit by these methods (e.g., via local regression) are stochastic, and hence
certain modifications to the algorithms and certainly to their convergence
analyses are required, as explained in Larson et al. (2019, Section 6.3). On
the algorithmic side, one modification that has been proposed is to introduce
adaptive schemes for Monte Carlo sampling that balance noise in $\widehat{U}$
with other errors (Shashaani et al. 2018).

##### Direct search methods.

Direct search methods (Torczon 1997, Larson et al. 2019) are also good choices
for (94). In general, these methods do not explicitly build local
approximations of the objective or attempt to approximate its gradient, but
rather compute optimization steps based on the relative ordering of the
evaluated function values. A classical pattern search approach is the
Nelder–Mead nonlinear simplex (Nelder and Mead 1965), which has been explored
for OED by Huan and Marzouk (2013). Modifications to Nelder–Mead for strongly
noisy objectives include re-sampling schemes to account for the impact of
noise on rankings (Chang 2012). Other direct search methods, such as
generalized pattern search methods (Audet and Dennis 2002, Audet 2004), have
variants for stochastic objectives as well (Sriver et al. 2009).

####  4.3.2 Gradient-based methods

##### Robbins–Monro.

The Robbins–Monro (RM) algorithm (Robbins and Monro 1951)—the progenitor of
modern stochastic gradient descent algorithms—follows (96) but requires access
to an unbiased estimator $\hat{g}$ of the gradient $\nabla_{\xi}\bar{U}$. That
is, we need a random $\hat{g}$ satisfying

$\mathbb{E}[\hat{g}(\xi)]=g(\xi)\coloneqq\nabla_{\xi}\bar{U}(\xi).$

One way of constructing such an estimator is to reparameterize $\widehat{U}$
such that the distribution over which we take the expectation is functionally
independent of the design variable $\xi$. Such a reparameterization is always
feasible (Mohamed et al. 2020). For example, in (95) where
$W=(Y^{(i)},\Theta^{(i)})_{i=1}^{N}$, we initially have
$\mathbb{E}_{W|\xi}=\mathbb{E}_{Y^{(1:N)}|\theta^{(1:N)},\xi}\,\mathbb{E}_{\Theta^{(1:N)}}$,
but it is easy to replace the design-dependent observations $Y$ in this
expectation with random variables that do not depend on $\xi$. If the
observations are described as $Y=G_{\xi}(\theta)+\mathcal{E}$, for instance,
with an observation noise $\mathcal{E}$ whose distribution is independent of
the design, this goal is immediately achieved: we can replace the original
expectation with
$\mathbb{E}_{\mathcal{E}^{(1:N)}}\,\mathbb{E}_{\Theta^{(1:N)}}$ and
reparameterize the expected utility estimator accordingly. More fundamentally,
however, the source of randomness in any statistical model is always a
transformation of some _independent_ random input. Call this random input
$\check{W}\equiv\left(\check{W}^{(i)}\right)_{i=1}^{N}\sim p_{\check{W}}$. (In
the example we just raised,
$\check{W}^{(i)}=(\mathcal{E}^{(i)},\Theta^{(i)})$.) After reparameterization,
the estimator (95) of the expected utility is rewritten as

$\displaystyle\check{U}(\xi,\check{W})=\frac{1}{N}\sum_{i=1}^{N}\check{u}\left(\xi,\check{W}^{(i)}\right),$
(98)

and satisfies
$\mathbb{E}_{\check{W}}[\check{U}(\xi,\check{W})]=\mathbb{E}_{W|\xi}[\widehat{U}(\xi,W)]=\bar{U}(\xi)$.
Since $\check{W}$ is independent of $\xi$, we can then write

$\nabla_{\xi}\bar{U}(\xi)=\nabla_{\xi}\mathbb{E}_{\check{W}}[\check{U}(\xi,\check{W})]=\mathbb{E}_{\check{W}}[\nabla_{\xi}\check{U}(\xi,\check{W})]=\mathbb{E}_{\check{W}}[\hat{g}(\xi,\check{W})],$

where

$\displaystyle\hat{g}(\xi,\check{W})\coloneqq\nabla_{\xi}\check{U}(\xi,\check{W})=\frac{1}{N}\sum_{i=1}^{N}\nabla_{\xi}\check{u}\left(\xi,\check{W}^{(i)}\right)$
(99)

is hence an unbiased estimator of the gradient.

Under appropriate conditions on the objective, the noise, and the step size
sequence $\\{\alpha_{k}\\}$, it can be shown that RM converges almost surely
to the global optimum of $\bar{U}$ in convex problems (Kushner and Yin 2003)
or to the first-order critical set in certain non-convex problems
(Mertikopoulos et al. 2020). Many improvements on the original RM algorithm
have been devised and widely deployed, including Polyak–Ruppert averaging
(Ruppert 1988, Polyak and Juditsky 1992), which can improve robustness and
convergence rates by averaging iterates along the optimization path, and
Nesterov acceleration (Nesterov 1983). RM-type algorithms have seen extensive
use in OED, for example by Huan and Marzouk (2014), Foster et al. (2019),
Carlon et al. (2020), Foster et al. (2020, 2021), Kleinegesse and Gutmann
(2020), Ivanova et al. (2021), Shen and Huan (2021), Zhang et al. (2021), and
Shen and Huan (2023).

##### Sample-average approximation.

Sample-average approximation (SAA) (also referred to as the retrospective
method or the sample-path method) (Shapiro 1991, Healy and Schruben 1991,
Gürkan et al. 1994, Kleywegt et al. 2002) takes a different strategy: it seeks
to reduce the stochastic objective to a deterministic one by fixing the
randomness throughout the entire optimization process. This again requires
reparameterizing the distribution over which we take the expectation (i.e.,
$\mathbb{E}_{W|\xi}$) to be independent of the design variable $\xi$, as in
the Monte Carlo estimator $\check{U}$ of (98). The SAA problem then becomes

$\displaystyle\xi^{\ast}_{\text{SAA}}\in\arg\max_{\xi\in\Xi}\check{U}(\xi,\check{W}=w_{s}),$
(100)

where $w_{s}$ is a realization of $\check{W}$ that is _fixed_ across $\Xi$.
SAA can also be viewed as an application of common random numbers, where
holding the sample of $\check{W}$ fixed essentially correlates the estimates
of $\bar{U}$ across the design space and yields a smoother objective surface.
(For further discussion of the relationships between common random numbers,
stochastic finite difference approximations to the gradient as in (97), and
stochastic gradient estimates as in (99), see Asmussen and Glynn (2007,
Chapter VII.2).) The optimization problem (100) is _deterministic_ , and hence
any standard deterministic optimization algorithm can be adopted. For
instance, Huan and Marzouk (2014) applied SAA to a nonlinear OED problem with
an EIG objective, and used the Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-
Newton scheme (Nocedal and Wright 2006, Chapter 6) to find optimal designs.
Stochastic bounds on the optimality gap
$\bar{U}(\xi^{\ast})-\bar{U}(\xi^{\ast}_{\text{SAA}})$ can also be obtained
Norkin et al. (1998), Mak et al. (1999).

####  4.3.3 Unbiased EIG gradient estimators

As discussed near the start of Section 4.3, the algorithms described thus far
focus on the objective $\bar{U}$; in other words, they assume that we have an
unbiased estimator of the desired objective at any $\xi$ (or, practically
speaking, that the bias is ‘small enough’ to be ignored). Section 4.3.2 then
describes how to obtain an unbiased estimator of the gradient of $\bar{U}$.
But we have not yet addressed what to do when $\bar{U}$ departs from the true
design objective $U$.

For the mutual information (EIG) objectives discussed throughout Section 3, we
do not have unbiased estimators; this is due to the presence of nested
expectations and hence nested Monte Carlo estimators, or the use of density or
density-ratio estimators within a nonlinear function (i.e., the logarithm).
Yet directly applying any of the optimization methods that we have just
discussed to a biased estimator of EIG can lead to arbitrarily large
departures from the true maximizer (recall the example of Figure 3.1), as the
objectives at hand may have multiple local maxima over $\Xi$ and the bias
itself may vary significantly as a function of $\xi$. Similarly, stochastic
gradient-based methods will converge to critical points specified by the
biased gradient, i.e., points where $\nabla_{\xi}\bar{U}(\xi)=0$, which again
differ from points where $\nabla_{\xi}U(\xi)=0$ if we cannot enforce
$\nabla_{\xi}\bar{U}(\xi)=\nabla_{\xi}U(\xi)$ for all $\xi$.

To remedy this situation for gradient-based optimization methods, Goda et al.
(2022) propose an unbiased estimator of the gradient of the EIG (14) with
respect to $\xi$, i.e., of $\nabla_{\xi}U_{\text{KL}}(\xi)$. Their approach
builds on the multilevel nested Monte Carlo estimator of Goda et al. (2020) by
combining it with the random truncation approach of Rhee and Glynn (2015),
which generally aims to ‘de-bias’ a consistent sequence of estimators. Goda et
al. (2022) use the Rhee and Glynn scheme to randomly truncate the infinite
telescoping sum in the inner loop of nested Monte Carlo, in a way that
guarantees unbiasedness of the overall gradient estimator, allowing it to
satisfy the requirements of the RM algorithm.

A different estimator of $\nabla_{\xi}U_{\text{KL}}(\xi)$ is proposed by Ao
and Li (2024). Here, the authors first use the reparameterization trick to
express the observations $Y$ as a function of the parameters $\Theta$ and
additional design-independent random variables $\mathcal{E}$, i.e., as
$Y=h(\Theta,\mathcal{E},\xi)$, and then apply the following identity to
rewrite the gradient of the log-evidence term in (60),

$\displaystyle\nabla_{\xi}\log
p_{Y|\xi}(y^{(i)}|\xi)=-\mathbb{E}_{\Theta|y^{(i)},\xi}\bigl{[}\nabla_{\xi}\log
p_{Y|\Theta,\xi}\bigl{(}h(\theta^{(i)},\mathcal{E}^{(i)},\xi)\,|\,\Theta,\xi\bigr{)}\bigr{]},$

where $y^{(i)}=h(\theta^{(i)},\mathcal{E}^{(i)},\xi)$ is any realization of
$Y$. MCMC sampling from the posterior $p(\theta|y^{(i)},\xi)$ yields an
estimate of $\nabla_{\xi}\log p(y^{(i)}|\xi)$, which is then embedded in an
outer loop of standard Monte Carlo sampling over the joint distribution of
$\Theta$ and $\mathcal{E}$. (One must also estimate the gradient of the
expected log-likelihood term in (60), but this is straightforward to do in an
unbiased way using the reparameterization trick.) Insofar as the inner
posterior samples are exact, the resulting nested estimator of
$\nabla_{\xi}U_{\text{KL}}(\xi)$ is unbiased. Of course, finite-length MCMC
chains produce biased estimates of the associated expectations, but the
authors demonstrate empirically that this bias can be made negligible.
Moreover, other posterior sampling schemes could be used instead.

####  4.3.4 Simultaneous bound tightening and design optimization

Focusing again on design objectives that involve EIG in parameters,
predictions, or some other quantity of interest, an important theme of
Sections 3.2 and 3.3 was the construction of variational lower bounds for the
EIG, i.e., functions $\mathcal{L}(f,\xi)\leq\text{EIG}(\xi)$. Most of these
lower bounds (with some exceptions, e.g., $\mathcal{L}^{\text{PCE}}$ (77)) are
parameterized by learnable functions $f$—which could take the form of
probability densities, transport maps, or more generic ‘critic’ functions.
Moreover, most of these learnable bounds can become tight for an appropriate
choice of $f$ (again with some exceptions, e.g., certain contrastive multi-
sample bounds in Section 3.3.2). In these cases, it is then natural to
maximimize _simultaneously_ over $f$ and $\xi$, as proposed in Foster et al.
(2020) and Kleinegesse and Gutmann (2020, 2021):

$\displaystyle\xi^{\ast}\in\operatornamewithlimits{argmax}_{\xi\in\Xi}\max_{f}\mathcal{L}(f,\xi)\,.$
(101)

This approach simultaneously seeks to tighten the lower bound and to find a
good design. The maximizer should in principle be a maximizer of the true EIG.

Recall that all of the variational bounds $\mathcal{L}$ were expressed as
expectations. To apply a stochastic gradient technique to (101), unbiased
estimates of gradients $\nabla_{\xi}\mathcal{L}(f,\xi)$ and
$\nabla_{f}\mathcal{L}(f,\xi)$ are required. As noted in Section 3.3.3, these
estimates are generally available through simple Monte Carlo estimation,
without the need for nested Monte Carlo or other complexities. Simultaneous
maximization (101) via stochastic gradient techniques has been used in the
setting of batch OED by Foster et al. (2020), Kleinegesse and Gutmann (2020,
2021), Zhang et al. (2021) and Orozco et al. (2024), and in the setting of
sequential experimental design by Ivanova et al. (2021) and Shen et al.
(2023), which we will discuss in the next section.

A related approach to leveraging bounds on EIG is proposed by Zheng et al.
(2020). Here, the authors take advantage of the ability to obtain _both_ lower
and upper bounds (in expectation) for the EIG—e.g., the prior-contrastive
estimator (77) as the lower bound and the standard NMC estimator (53) as the
upper bound. Both of these bounds can be refined and made tight by increasing
the inner-loop sample size $M$. Zheng et al. (2020) thus develop an adaptive
refinement strategy, in which the bounds are tightened over the course of
optimization, that achieves regret-style guarantees.

##  5 Sequential optimal experimental design

Sequential experimental design is concerned with the planning of multiple
experiments that are conducted in a sequence, where the _results_ of previous
experiments in the sequence can inform the design of subsequent experiments.
This is the crucial difference between sequential design and the ‘batch’ (also
called ‘static’) design approaches that were the focus of previous sections.
Here we will focus on a Bayesian approach to sequential experimental design.

One straightforward sequential design procedure is to apply the batch optimal
experimental design (OED) framework and methods from Sections 2–4 to one
experiment, or subset of experiments, at a time: optimally choose the next
design, perform that experiment, update the prior to the posterior based on
the outcome of the experiment, and repeat the process for the next. Such a
procedure is called _greedy_ or _myopic_ , because it does not take into
account future experiments when finding the (immediate) next experiment; it
involves no ‘lookahead.’ Yet greedy design is conceptually simple to
implement, especially if computational tools for batch design already exist,
and it is flexible for situations where the total number of desired
experiments is unknown. As a result, a large body of sequential experimental
design research has been based on some form of greedy design (Box 1992, Dror
and Steinberg 2008, Cavagnaro et al. 2010, Solonen et al. 2012, Drovandi et
al. 2013, 2014, Kim et al. 2014, Hainy et al. 2016, Kleinegesse et al. 2021).

Of course, to improve coordination among the experiments, one could design
_all_ the experiments simultaneously and thus revert to an overall batch
design. But doing so would forgo the opportunity to adapt to new observations;
in other words, it would allow for coordination among the experiments, but _no
feedback_. A hallmark of sequential experimental design is precisely the idea
of feedback.

We note that the greedy design approach described just now, and more generally
the sequential design approaches to be detailed in this section, differ from
the sequential algorithms for solving batch design problems discussed in
Sections 4.1 and 4.2. In particular, the latter are not for sequential
experimental design, as they do not involve observing the results of
experiments between design decisions, and thus do not incorporate any
feedback. Rather, they still seek a static design for a single batch of
experiments, but break the search into stages to control the dimension of the
design space and avoid combinatorial scaling of computational complexity.

In the remainder of this section, we will present a sequential optimal
experimental design formulation that includes both lookahead and feedback;
see, e.g., Müller et al. (2007), Von Toussaint (2011, VII.G), and Huan (2015,
Chapter 3). We refer to this formulation as _sequential OED_ (sOED). sOED
generalizes both greedy and batch design (Shen and Huan 2023, Section 2.3).
The key ideas of sOED are that (i) each design should be selected while taking
into consideration _all remaining experiments_ to be performed, and (ii)
designs should be given in the form of _functions_ , called policies, that
adaptively specify the next experiment in the sequence given the current state
of information. We will formalize these ideas using the framework of Markov
decision processes, to be presented shortly.

As we shall see, the numerical solution of the sOED problem is rather
challenging, and there have been relatively few attempts to solve it in great
generality. For example, Carlin et al. (1998), Gautier and Pronzato (2000),
Pronzato and Thierry (2002), Brockwell and Kadane (2003), Christen and
Nakamura (2003), Murphy (2003), Wathen and Christen (2006), Müller et al.
(2022), and Tec et al. (2023) have all made advances largely limited to
discrete settings, or did not employ a Bayesian framework with information-
theoretic design criteria. In this section, we will survey recent progress
towards realizing fully Bayesian sOED, including methods that make use of
dynamic programming Huan (2015), Huan and Marzouk (2016), reinforcement
learning Blau et al. (2022), Shen and Huan (2023), and information bounds
Foster et al. (2021), Ivanova et al. (2021), Shen et al. (2023).

###  5.1 Background

We focus on the design of a finite number of experiments, indexed by
$k=0,1,\ldots,N-1$. We assume $N$ is known and fixed. The Bayesian update for
the $k$th experiment then becomes

$\displaystyle p(\theta|y_{k},\xi_{k},{\sf
i}_{k})=\frac{p(y_{k}|\theta,\xi_{k},{\sf i}_{k})\,p(\theta|\kern
0.50003pt{\sf i}_{k})}{p(y_{k}|\xi_{k},{\sf i}_{k})},$ (102)

where $\xi_{k}$ and $y_{k}$ are, respectively, the design and the value of the
observation realized in the $k$th experiment, and ${\sf
i}_{k}\coloneqq\left[\xi_{0},y_{0},\dots,\xi_{k-1},y_{k-1}\right]$ is a
realization of the random ‘background information’ sequence
$I_{k}\coloneqq\left[\xi_{0},Y_{0},\dots,\xi_{k-1},Y_{k-1}\right]$ composed of
the history of designs and observations from all experiments preceding the
current one, with $I_{0}=\emptyset$. The prior density $p(\theta|\kern
0.50003pt{\sf i}_{k})$ represents the state of knowledge about the uncertain
parameters $\Theta$ before the $k$th experiment, and the posterior density
$p(\theta|y_{k},\xi_{k},{\sf i}_{k})$ represents the updated state of
knowledge after having observed the outcome of the $k$th experiment. Equation
(102) uses the simplification $p(\theta|\xi_{k},{\sf i}_{k})=p(\theta|\kern
0.50003pt{\sf i}_{k})$, since the prior density does not depend on the pending
choice of design. The posterior after the $k$th experiment
$p(\theta|y_{k},\xi_{k},{\sf i}_{k})=p(\theta|\kern 0.50003pt{\sf i}_{k+1})$
then becomes the prior for the $(k+1)$th experiment, and (102) can be applied
recursively. In (102) we present the general setting where the density of the
current data $p(y_{k}|\theta,\xi_{k},{\sf i}_{k})$ may depend on the design
and observations from previous experiments, i.e., on ${\sf i}_{k}$. In many
situations, however, $Y_{k}$ is conditionally independent of the other
observations in the sequence given $\theta$ and $\xi_{k}$, in which case its
density simplifies to $p(y_{k}|\theta,\xi_{k},{\sf
i}_{k})=p(y_{k}|\theta,\xi_{k})$.

We note that the experiments in the sequence do not need be of the same type,
as long as they share the same parameters $\Theta$. The spaces of possible
designs $\Xi_{k}\ni\xi_{k}$ and observations $\mathcal{Y}_{k}\ni y_{k}$ can
differ from one experiment to the next, and even change dimension. Similarly,
the conditional density of the data in (102) is, in more explicit notation,
$p_{Y_{k}|\Theta,\xi_{k},I_{k}}(y_{k}|\theta,\xi_{k},{\sf i}_{k})$ and thus
has a $k$-dependence as well. For example, in learning the properties of a
fluid, a first experiment might entail selecting the shape of an obstacle to
place in the flow and observing the velocity in its wake, while the next
experiment might involve choosing where on the surface of this obstacle to
place a pressure sensor.

###  5.2 Formulation as a Markov decision process

Sequential experimental design can be modeled through a Markov decision
process (MDP) specified by a tuple
$\bigl{(}\mathcal{S},\\{\mathcal{A}_{k}\\}_{k},s_{0},\\{r_{k}(\cdot)\\}_{k},\\{T_{k}(\cdot)\\}_{k}\bigr{)}$.
$\mathcal{S}$ is the state space of the Markov process, where the state
$S_{k}$ is a random element with values in $\mathcal{S}$ and realizations of
this state are denoted by $s_{k}\in\mathcal{S}$; the action spaces
$\mathcal{A}_{k}$ comprise possible actions (which here are designs)
$\xi_{k}\in\mathcal{A}_{k}$;555$\mathcal{A}_{k}\equiv\Xi_{k}$, i.e., spaces of
candidate designs indexed by $k$, but we use $\mathcal{A}_{k}$ in this section
to be consistent with the usual notation in the MDP literature. $s_{0}$ is an
initial state; $r_{k}(s_{k},\xi_{k},y_{k})$ are scalar-valued reward functions
that evaluate the instantaneous reward when taking action $\xi_{k}$ and
observing $y_{k}$ at state $s_{k}$; and
$T_{k}(\mathcal{S}^{\prime}|s_{k},\xi_{k})$ are state transition kernels that
evaluate the probability of transitioning to any set of states
$\mathcal{S}^{\prime}\subseteq\mathcal{S}$ at stage $k+1$ having taken action
$\xi_{k}$ at state $s_{k}$. In the context of experimental design, the action
being taken is the selection of a design; thus we use the terms ‘action’ and
‘design’ interchangeably.

State. The state of the system before the $k$th experiment is described by
$S_{k}=\bigl{(}S_{k}^{b},S_{k}^{p}\bigr{)}$, a quantity that summarizes all
information deemed relevant to future design decisions. We split $S_{k}$ into
a ‘belief state’ $S_{k}^{b}$, representing the state of knowledge/uncertainty
in the parameters $\Theta$, and a ‘physical state’ $S_{k}^{p}$ comprising any
other deterministic variables that may be relevant to the design process.
Since $\Theta$ is not directly observed and can only be inferred from
observations $Y_{k}$, this setup can be viewed as a partially observed MDP
(POMDP) for $\Theta$ or a belief MDP with (belief) state $S_{k}$ (Kaelbling et
al. 1998).

A realization of the belief state, $s_{k}^{b}$, is simply the posterior
distribution of $\Theta$ given a set of past experimental designs and realized
observations, ${\sf i}_{k}$. For $\Theta$ taking values on $\mathbb{R}^{p}$,
$s_{k}^{b}$ is thus a probability distribution on $\mathbb{R}^{p}$.
Numerically, it can be represented by, for example, a probability density
function approximation, or a set of weighted particles, or by tracking ${\sf
i}_{k}$ directly. Tracking ${\sf i}_{k}$ is easiest to implement since it does
not require additional calculations translating $(\xi_{i})_{i<k}$ and
$(y_{i})_{i<k}$ to another representation, but the dimension of ${\sf i}_{k}$
grows with $k$—though it is bounded for finite $N$. In general, the set of
possible posterior distributions that can be realized is uncountably
infinite,666Unless both $\xi_{k}$ and $Y_{k}$ are discrete. and hence this
setting differs from a discrete or finite-state system.

Maintaining only a belief state, in the form of the posterior, does not
suffice to preserve the Markov property of the system if the likelihood
depends on the history of past experiments ${\sf i}_{k}$, as presented in
(102). This can be fixed by introducing a physical state. With regard to the
${\sf i}_{k}$-dependence in $p(y_{k}|\theta,\xi_{k},{\sf i}_{k})$, the
physical state essentially extracts and tracks relevant features from ${\sf
i}_{k}$ that allow the likelihood to be evaluated or the observations $Y_{k}$
simulated. An example of the physical state would be the known position of a
mobile sensor platform, which might evolve from one stage of the design
sequence to the next. Note that if $I_{k}$ is adopted as the belief state,
then information about the physical state is already contained in $I_{k}$,
even if only implicitly.

Action (design) and policy. Sequential experimental design is adaptive in
nature. Whereas a batch design problem seeks a single design $\xi$ (see
Section 2), sequential design now looks for a strategy, called a _policy_ ,
describing _how to choose_ the design depending on the current state. The
policy is a collection of functions
$\pi=\\{\mu_{k}:\mathcal{S}\to\mathcal{A}_{k},k=0,\ldots,N-1\\}$, where the
policy function $\mu_{k}$ returns the design for the $k$th experiment given
the current state, $\xi_{k}=\mu_{k}(s_{k})$. In general, $\mu_{k}$ differs
from experiment to experiment for finite $N$.777In contrast, in the infinite-
horizon setting, the policy must be stationary, i.e., independent of $k$. Thus
only a single $\mu$ needs to be found, making infinite-horizon problems
generally easier to tackle. Some intuition for this fact follows by
considering that even when starting from the same belief and physical state, a
‘good’ design can be quite different depending on how many experiments remain
in the overall sequence. We focus on deterministic policies in this
formulation, although stochastic policies that evaluate the probability of
choosing different candidate designs can also be adopted.

State transition dynamics. When an experiment is performed, the state changes
according to a transition kernel $T_{k}(\mathcal{S}^{\prime}|s_{k},\xi_{k})$
describing the probability of transitioning from the current state $s_{k}$,
having chosen design $\xi_{k}$ and observed the outcome of resulting
experiment, to any set of states at stage $k+1$,
$\mathcal{S}^{\prime}\subseteq\mathcal{S}$. This kernel is generally
intractable to evaluate, but it can instead be _simulated_ , by sampling from
the prior predictive distribution of $Y_{k}$ given the design $\xi_{k}$ and
then applying Bayes’ rule (102). Of course, both of these steps make use of
the statistical model for $Y_{k}$. We denote the latter transition dynamics by
$s_{k+1}=\mathcal{F}_{k}(s_{k},\xi_{k},y_{k})$; the function $\mathcal{F}_{k}$
encapsulates the transition from prior to posterior, given values of $\xi_{k}$
and the realized data $y_{k}$, following (102).

The physical state, if present, evolves according to a model for the relevant
physical process.

Reward (utility). Here $r_{k}(s_{k},\xi_{k},y_{k})$ denotes the real-valued
reward immediately obtained after performing the $k$th experiment, which may
depend on the state, design, and observation values. For example, the reward
can reflect the cost of carrying out the $k$th experiment and/or the
information gained in $\Theta$ as a result. On the other hand, $r_{N}(s_{N})$
denotes the terminal reward, which reflects any rewards that can only be
quantified after all experiments are completed. The choice of reward functions
$(r_{k})_{k=0}^{N}$ is quite flexible and can be based on the wide range of
optimal design utilities discussed in Section 2.2. In the next section we will
highlight the reward functions corresponding to $u^{\text{div}}$ (20), which
reflects information gain in the model parameters $\Theta$.

The overall progression of an MDP for the sequential design of $N$ experiments
is depicted in Figure 5.1, adapted from Shen and Huan (2023).

Figure 5.1: In the MDP progression of sequential experimental design, we start
with an initial state $s_{0}$ (i.e., initial prior and physical state),
evaluate the policy function at the state $\mu_{0}(s_{0})$ to obtain the
design $\xi_{0}$ for experiment 0, conduct the experiment to obtain its
outcome $y_{0}$ and immediate reward $r_{0}(s_{0},\xi_{0},y_{0})$, update the
state to the new state via the transition dynamics
$s_{1}=\mathcal{F}_{0}(s_{0},\xi_{0},y_{0})$ (i.e., updated posterior and
physical state), and repeat for the next experiments. Once the last experiment
$N-1$ is completed, the terminal state
$s_{N}=\mathcal{F}_{N-1}(s_{N-1},\xi_{N-1},y_{N-1})$ can be computed along
with the corresponding terminal reward $r_{N}(s_{N})$. Figure adapted from
Shen and Huan (2023).

###  5.3 Problem statement

The sOED problem entails finding a design policy $\pi^{\ast}$ that maximizes
the expected utility $U(\pi)$:

$\displaystyle\pi^{\ast}\in\operatornamewithlimits{argmax}_{\pi=\\{\mu_{0},\ldots,\mu_{N-1}\\}}$
$\displaystyle\ \
\left\\{U(\pi)\coloneqq\mathbb{E}_{Y_{0:N-1}|\pi,s_{0}}\left[\sum_{k=0}^{N-1}r_{k}(S_{k},\xi_{k},Y_{k})+r_{N}(S_{N})\right]\right\\}$
(103) subject to $\displaystyle\ \ \
\xi_{k}=\mu_{k}(S_{k})\in\mathcal{A}_{k},$ $\displaystyle\ \ \
S_{k+1}=\mathcal{F}_{k}(S_{k},\xi_{k},Y_{k}),\quad\text{for}\,\,k=0,\dots,N-1.$

Here $Y_{0:N-1}=Y_{0},Y_{1},\ldots,Y_{N-1}$. The initial state $s_{0}$ is
assumed known; if it is not known, then another expectation can be taken over
the random initial state $S_{0}$. Optionally, a discount factor
$\gamma\in(0,1]$ can multiply the reward at each successive stage (hence
yielding terms $\gamma^{k}r_{k}$) to artificially reduce the value of rewards
obtained further in the future. The expected utility here is also known as the
_expected return_ (or expected total reward) in MDP terminology; we will use
these terms interchangeably.

In the field of reinforcement learning (Kaelbling et al. 1996, Sutton and
Barto 2018), problem (103) corresponds to a model-based planning problem,
described by a finite-horizon belief-state MDP with continuous action and
observation spaces. Embedded in each transition $\mathcal{F}_{k}$ is a step of
Bayesian inference, which can be quite expensive to perform, especially in
nonlinear and non-Gaussian settings with computationally intensive
likelihoods.

###  5.4 Solution approaches

In general, problem (103) cannot be solved analytically. Different numerical
strategies must be adopted to find an approximate solution. Before we
highlight representative approaches from recent literature, we introduce two
types of value function that are central to the MDP formulation, as many
solution strategies are built on approximating these functions.

The _action-value functions_ (or _Q-functions_) corresponding to a policy
$\pi$ are

$\displaystyle Q_{k}^{\pi}(s_{k},\xi_{k})$
$\displaystyle=\mathbb{E}_{Y_{k:(N-1)}|\pi,s_{k},\xi_{k}}\left[r_{k}(s_{k},\xi_{k},Y_{k})+\sum_{t=k+1}^{N-1}r_{t}(S_{t},\mu_{t}(S_{t}),Y_{t})+r_{N}(S_{N})\right]$
(104)
$\displaystyle=\mathbb{E}_{Y_{k}|\pi,s_{k},\xi_{k}}\left[r_{k}(s_{k},\xi_{k},Y_{k})+Q^{\pi}_{k+1}(S_{k+1},\mu_{k+1}(S_{k+1}))\right],$
(105)

for $k=0,\ldots,N-1$ and subject to
$S_{t+1}=\mathcal{F}_{t}(S_{t},\xi_{t},Y_{t})$, with
$Q_{N}^{\pi}(s_{N})=r_{N}(s_{N})$. The value of the Q-function
$Q_{k}^{\pi}(s_{k},\xi_{k})$ is the expected remaining cumulative reward
(i.e., the expected sum of all remaining rewards) for performing the $k$th
experiment at design $\xi_{k}$ from state $s_{k}$ and thereafter following
policy $\pi$.

The _state-value functions_ (or _V-functions_) corresponding to a policy $\pi$
are

$\displaystyle V_{k}^{\pi}(s_{k})$
$\displaystyle=\mathbb{E}_{Y_{k:(N-1)}|\pi,s_{k}}\left[r_{k}(s_{k},\mu_{k}(s_{k}),Y_{k})+\sum_{t=k+1}^{N-1}r_{t}(S_{t},\mu_{t}(S_{t}),Y_{t})+r_{N}(S_{N})\right]$
(106)
$\displaystyle=\mathbb{E}_{Y_{k}|\pi,s_{k}}\left[r_{k}(s_{k},\mu_{k}(s_{k}),Y_{k})+V^{\pi}_{k+1}(S_{k+1})\right],$
(107)

for $k=0,\ldots,N-1$ and subject to
$S_{t+1}=\mathcal{F}_{t}(S_{t},\xi_{t},Y_{t})$, with
$V_{N}^{\pi}(s_{N})=r_{N}(s_{N})$. The value of the V-function
$V_{k}^{\pi}(s_{k})$ is the expected remaining cumulative reward starting from
a given state $s_{k}$ and following policy $\pi$ for all remaining
experiments. The expected utility in (103) can be succinctly written as
$U(\pi)=V_{0}^{\pi}(s_{0})$.

The V-function and Q-function are related to each other via

$\displaystyle V_{k}^{\pi}(s_{k})=Q_{k}^{\pi}(s_{k},\mu_{k}(s_{k})),$ (108)

and thus often only one of the two value functions needs to be solved for. We
note that both value functions can also be expressed in the recursive forms
(105) and (107). When the policy is the optimal policy $\pi^{\ast}$ from
(103), the recursive relations become the well-known _Bellman optimality
equations_ :

$\displaystyle V_{k}^{\ast}(s_{k})$