CMlargesymbols”00 CMlargesymbols”01 Optimal experimental design: Formulations
and computations–References

# Optimal experimental design:
Formulations and computations

Xun Huan
University of Michigan  1231 Beal Ave  Ann Arbor  MI 48109  USA
Email<EMAIL_ADDRESS>
Jayanth Jagalur
Center for Applied Scientific Computing  Lawrence Livermore National
Laboratory
7000 East Ave  Livermore  CA 94550  USA
Email<EMAIL_ADDRESS>
Youssef Marzouk
Massachusetts Institute of Technology
77 Massachusetts Ave  Cambridge  MA 02139  USA
Email<EMAIL_ADDRESS>
###### Abstract

Questions of ‘how best to acquire data’ are essential to modeling and
prediction in the natural and social sciences, engineering applications, and
beyond. Optimal experimental design (OED) formalizes these questions and
creates computational methods to answer them. This article presents a
systematic survey of modern OED, from its foundations in classical design
theory to current research involving OED for complex models. We begin by
reviewing criteria used to formulate an OED problem and thus to encode the
goal of performing an experiment. We emphasize the flexibility of the Bayesian
and decision-theoretic approach, which encompasses information-based criteria
that are well-suited to nonlinear and non-Gaussian statistical models. We then
discuss methods for estimating or bounding the values of these design
criteria; this endeavor can be quite challenging due to strong nonlinearities,
high parameter dimension, large per-sample costs, or settings where the model
is implicit. A complementary set of computational issues involves optimization
methods used to find a design; we discuss such methods in the discrete
(combinatorial) setting of observation selection and in settings where an
exact design can be continuously parameterized. Finally we present emerging
methods for sequential OED that build non-myopic design policies, rather than
explicit designs; these methods naturally adapt to the outcomes of past
experiments in proposing new experiments, while seeking coordination among all
experiments to be performed. Throughout, we highlight important open questions
and challenges.

2020 Mathematics Subject Classification: Primary 62-02, 62-08, 62B15, 62K05

Secondary 62L05, 94A17, 65M32

###### doi:

XXXXXXXX

###### Contents

  1. 1 Introduction
    1. 1.1 Scope and organization of this article
  2. 2 Optimal design criteria
    1. 2.1 Design criteria for the linear-Gaussian setting
      1. 2.1.1 Classical alphabetic optimality
      2. 2.1.2 Bayesian alphabetic optimality
      3. 2.1.3 Designs as probability measures
      4. 2.1.4 Challenges of nonlinear design
    2. 2.2 Decision- and information-theoretic formulations
      1. 2.2.1 Expected information gain in parameters
      2. 2.2.2 Other utility functions
    3. 2.3 Design criteria for infinite-dimensional problems
      1. 2.3.1 Inverse problems in the Bayesian setting
      2. 2.3.2 Gaussian process regression
    4. 2.4 Related problems and their distinctions
  3. 3 Numerical approximation of design criteria
    1. 3.1 Nested Monte Carlo estimators
    2. 3.2 Mutual information bounds via density approximations
    3. 3.3 More general variational bounds for mutual information
      1. 3.3.1 Variational lower bounds and approximate density ratios
      2. 3.3.2 Multi-sample lower bounds
      3. 3.3.3 Deploying variational bounds
    4. 3.4 Low-dimensional structure
      1. 3.4.1 Linear models
      2. 3.4.2 Nonlinear models
  4. 4 Design optimization methods
    1. 4.1 Optimization methods for linear design
      1. 4.1.1 Local search and exchange algorithms
      2. 4.1.2 Continuous convex relaxation approaches
      3. 4.1.3 Sequential algorithms for finding optimal designs
    2. 4.2 Combinatorial approaches for discrete design variables
      1. 4.2.1 Background
      2. 4.2.2 Set function attributes of design criteria
      3. 4.2.3 Greedy algorithms for cardinality-constrained designs
      4. 4.2.4 Beyond greedy: algorithms for more general design problems
    3. 4.3 Optimizing real-valued design variables
      1. 4.3.1 Derivative-free methods
      2. 4.3.2 Gradient-based methods
      3. 4.3.3 Unbiased EIG gradient estimators
      4. 4.3.4 Simultaneous bound tightening and design optimization
  5. 5 Sequential optimal experimental design
    1. 5.1 Background
    2. 5.2 Formulation as a Markov decision process
    3. 5.3 Problem statement
    4. 5.4 Solution approaches
      1. 5.4.1 Approximate dynamic programming (ADP-sOED)
      2. 5.4.2 Actor-critic policy gradient (PG-sOED)
      3. 5.4.3 Methods that leverage information bounds
    5. 5.5 Open questions in sequential optimal design
  6. 6 Outlook
    1. 6.1 Model misspecification
    2. 6.2 Risk-aware design criteria
    3. 6.3 OED in practice

##  1 Introduction

Acquiring data to inform models and guide decisions is an essential part of
scientific inquiry, engineering design, and even policy making. Seldom can we
construct a useful model in isolation from data. Rather, data must be used to
infer parameters of models, to assess whether models can provide useful
predictions, and to spur a wide variety of model improvements. In this
setting, it is natural to consider how to acquire data _efficiently_.
Experimental data and field measurements can be costly or time-consuming to
acquire; other information sources, similarly, may be expensive to query. Yet
we face a multitude of choices in designing such queries. Where to place a
sensor? What experimental conditions to impose? What quantity to observe? Do
measurements need to be very precise, or would a noisier measurement suffice?
When should measurements be made? And how much data should be collected? More
broadly, what _combination_ of observations would be most informative or
useful—and how should we precisely define notions of ‘informative’ or ‘useful’
in the first place? Crucially, these notions must allow experimental choices
to be made _before_ data are acquired.

_Optimal experimental design_ (OED) aims to answer these questions, through
mathematical formulations that formalize and tailor them to the ultimate goals
of acquiring data. A model developer may have many possible goals, and hence
there are many possible criteria for what comprises a good experimental
design. Another essential aspect of OED involves numerical algorithms, e.g.,
for evaluating suitable design criteria, for optimizing over possible
experimental configurations, and possibly doing so in a ‘closed-loop’
sequential fashion. Collectively, these endeavors lie at the intersection of
many fields: statistical inference and decision theory, information theory,
Monte Carlo methods, continuous and combinatorial optimization, dynamic
programming and stochastic control, and even reinforcement learning.

Every OED problem has two essential ingredients: an _experiment_ , which is
the source of data, and a mathematical _model_. The role of the model is to
simulate what might happen in candidate experiments, and to assess how the
results of such experiments might improve the model and its predictions.
Formally, the model is a statistical model; in some applications, evaluating
this statistical model also involves significant numerical simulation. An
underlying presumption of OED is that it is _worthwhile_ to perform many
calculations involving the model, in order to plan experiments that are more
efficient and effective. These calculations might be far less expensive (by
some metric) than performing experiments, or there might be other impediments
to experimentation that make finding an optimal design in advance, or building
an online optimal design policy for online settings, worth the effort.

The notion of an ‘experiment’ should be construed quite broadly, and certainly
not limited to laboratory experiments in the traditional sense. A high-
fidelity simulation of a complex model, used to produce data to calibrate the
parameters of a simpler model, constitutes an experiment. Arranging a network
of sensors, or planning the path of an airborne vehicle carrying instruments,
also constitutes designing an experiment. Application domains in which OED is
used are similarly vast. OED has long been an essential part of statistical
modeling, from the design of clinical trials (Berry et al. 2010) to the design
of computer experiments (Sacks et al. 1989); the latter is closely related to
the classical problem of design for regression (Elfving 1952, Kiefer and
Wolfowitz 1959, Kiefer 1961a). But OED can also be applied to problems
involving parameter inference in ordinary or partial differential equations
(Huan and Marzouk 2013), to a wide range of inverse problems (Haber et al.
2008, Alexanderian et al. 2016b, Ruthotto et al. 2018, Alexanderian 2021,
Helin et al. 2022), and to myriad other ‘complex’ models of data-generating
processes—in astronomy (Loredo 2011), systems biology (Liepe et al. 2013),
aeroelasticity (Riley et al. 2019), reliability testing (Weaver and Meeker
2021), nuclear physics (Melendez et al. 2021), and beyond.

The evolution of experimental design has a rich and fascinating history. Early
twentieth-century approaches to experimental design were motivated by
agricultural experiments and similar applications. Statistical methods in this
setting often involved hypothesis testing, and a good design was one that
maximized the sensitivity of the test. Notable works by Fisher and his
collaborators during this period introduced concepts such as balance,
orthogonality, blocking, and aliasing (Fisher 1936, Craig and Fisher 1936,
Yates 1933, 1937, 1940, Bose 1939, Bose and Nair 1939). Wald recognized that
these ideas were relevant to many other fields, and in a seminal paper (Wald
1943) introduced formal notions of the efficiency of a design. With this
framework, he was able to explain the success of designs based on Latin
squares, commonly used in agricultural experimentation. Following Wald’s work,
many core ideas of modern OED emerged towards the middle of the twentieth
century (Elfving 1952, Lindley 1956, Kiefer 1958, Stone 1959, Kiefer and
Wolfowitz 1959, Kiefer 1959, 1961a), some of them influenced by the emerging
discipline of decision theory. Kiefer distinguished various design criteria by
giving them meaningful names, and thus originated current ‘alphabetic
optimality’ terminology (Kiefer 1958). Much later, Kiefer also showed inter-
relationships among various design criteria by introducing more general
notions of ‘universal’ optimality (Kiefer 1974). Lindley’s work in the same
period (Lindley 1956) aligns more closely with Bayesian statistical thinking,
and deserves special mention in light of the current popularity of information
theoretic design criteria. These approaches remained largely intractable half
a century ago, but with advances in computing power and algorithms, these more
general and arguably more rigorous design criteria have become feasible to
implement. For a more detailed history of OED, we refer readers to Wynn (1984)
and to the texts by Fedorov (1972), Shah and Sinha (1989), and Pukelsheim
(2006).

In recent years, interest in OED has expanded from the statistics literature
into the uncertainty quantification and applied mathematics communities,
where, as noted earlier, an animating goal has been to advance OED
methodologies for ‘complex’ models. Here many forms of complexity are
relevant: high-dimensional parameters, computationally intensive likelihood
functions that involve the evaluation of ordinary or partial differential
equations, strong nonlinearity in the dependence of observables on parameters,
and the associated non-Gaussianity of posterior distributions in the Bayesian
setting. Additional complexities can arise due to data-generating processes
that evolve in time, which present the opportunity to design and implement
experiments adaptively, i.e., where the results of previous experiments
influence the next; this is the setting of sequential optimal experimental
design. Another thread relevant to modern OED has come from the computer
science literature, where much attention has been paid to the optimization of
_set functions_ ; methods here underpin a combinatorial approach to OED, where
the problem is cast as choosing a subset of a given ‘ground set’ of feasible
candidate experiments. And in recent years, tools from machine learning and in
particular deep learning have become quite useful for OED, for instance by
offering new ways of evaluating complex design criteria in high dimensional,
non-Gaussian settings. Of particular interest are expressive machine learning
models and learning algorithms for the underlying density (or density ratio)
estimation tasks. We will discuss all of these threads, and many more, in the
ensuing sections.

We mention here several other excellent surveys that may be of interest to the
reader. Steinberg and Hunter (1984) and Pukelsheim (2006) (updated from the
original 1993 version) provide comprehensive reviews of non-Bayesian OED for
linear models. Ford et al. (1989) discuss design for nonlinear models, while
DasGupta (1995) discusses Bayesian designs, mostly for linear models. Atkinson
et al. (2007) provide extensive coverage of linear optimal design theory, with
some extension to nonlinear and Bayesian methods. The paper of Chaloner and
Verdinelli (1995) is an influential review from a statistical perspective,
highlighting features of the Bayesian and decision-theoretic approach to OED.
We view the Bayesian approach as very natural in the setting of design, and
will largely adopt such a perspective here. Clyde (2001) presents an overview
of Bayesian OED formulations with various choices of utility.

More recent reviews have placed a greater emphasis on computation. Ryan et al.
(2016) discuss Bayesian formulations of OED and then survey computational
algorithms for realizing Bayesian designs, emphasizing Monte Carlo methods for
estimating design criteria and for searching through design space.
Alexanderian (2021) reviews OED for Bayesian inverse problems, emphasizing
formulations and algorithms for the infinite-dimensional (function space)
setting. Rainforth et al. (2023) provide a survey highlighting recent machine
learning and reinforcement learning tools in OED, including sequential design.
Strutz and Curtis (2024) present a review of variational OED methodologies and
their application to geophysical and in particular seismological problems.
There may be other recent reviews of which we are unaware, and we apologize
for such omissions.

###  1.1 Scope and organization of this article

This article aims to provide a broad, comprehensive survey of methodologies
for optimal experimental design. Our perspective covers both
formulations—i.e., the many ways of posing an OED problem—and
computations—i.e., numerical algorithms for finding optimal designs, or
tractable approximations of optimal designs, for a range of problem settings.
The second topic in particular is multi-faceted—we will draw upon Monte Carlo
methods, algorithms for inference and density estimation, and a variety of
optimization approaches—but naturally enjoys a close interplay with the first.
Our goal is to guide readers who are new to OED from the basic ideas to the
research frontier, and to illuminate open issues and challenges at that
frontier. Indeed, we believe that the present moment is ripe for a survey of
the field: the problems that remain are quite challenging, and a synthesis of
ideas and approaches from many different intellectual communities (some of
which have been rather disconnected) is needed to make progress.

The article is organized as follows. We begin in Section 2 by describing
possible design criteria for OED, each encoding a different notion of what
constitutes a ‘good’ experiment. The applicability of these criteria ranges
from the rather specific (e.g., linear-Gaussian problems) to the very general.
Many are rooted in information-theoretic and/or decision-theoretic
formulations. We also mention alternative design heuristics that have been
used in practice, and clarify distinctions between the OED problem and several
related but different problems, such as Bayesian optimization.

In Section 3 we turn to the first step of computation: numerical algorithms
for estimating or bounding the values of these design criteria. We survey
Monte Carlo schemes, as well as variational approximations and density
estimation methods, for this purpose. Some of these algorithms are applicable
to so-called ‘implicit’ Bayesian models, where evaluations of the likelihood
function or prior density may be unavailable. We also discuss challenges
associated with high-dimensional parameters and data, and dimension reduction
schemes that mitigate these challenges. In Section 4 we discuss strategies and
algorithms for efficiently optimizing design criteria in various settings. On
one hand, we address problems where an exact design of interest is represented
by continuous (real-valued) variables. On the other hand, we discuss settings
where the design optimization problem is discrete and combinatorial, e.g., a
form of subset selection, and continuous relaxations thereof.

Sections 2–4 focus on the all-at-once (‘batch’) design of experiments—that is,
‘fixed’ or static designs that cannot be adjusted as the outcomes of
experiments are realized. In contrast, Section 5 turns to the problem of
sequential OED, where design decisions are naturally adapted according to the
outcomes of previous experiments, while taking into consideration the
information to be gained from future experiments. We present sequential OED in
a fully Bayesian setting, leveraging the formalism of Markov decision
processes. We then highlight recent computational methods for solving this
challenging problem, which make use of dynamic programming, various
reinforcement learning techniques, and information bounds.

We close in Section 6 with a discussion of open questions and opportunities
for future work.

##  2 Optimal design criteria

We begin by addressing the central question of any OED formulation: In what
sense should one deem a candidate design to be ‘good?’ More specifically, by
what considerations should one candidate design be deemed better than another?
Answering these questions is essential to the notion of _optimal_ design. The
answers are formalized by choosing a quantitative design criterion—a function
that can then be maximized in order to identify the corresponding optimal
design. In this section, we will discuss a wide variety of design criteria and
the goals they encode.

To formulate these criteria, we must first specify a _statistical model_ for
the observations $Y$ obtained under a candidate design. We need such a model
in order to predict the outcomes of candidate experiments, and to relate these
outcomes to the ensuing estimation or prediction tasks that motivated the
acquisition of data in the first place. We will initially consider parametric
statistical models. Any such model is a family of probability distributions
for $Y$, indexed by (unknown) model parameters
$\theta\in\bm{\Theta}\subseteq\mathbb{R}^{p}$ and by the choice of design
$\xi\in\Xi$:

$\mathcal{M}=\\{\textit{{F}}_{Y;\theta,\xi}:\theta\in\bm{\Theta},\
\xi\in\Xi\\}.$ (1)

This model encodes, for any given value of $\theta$ and $\xi$, a complete
probabilistic description of the resulting observations, via the cumulative
distribution function $\textit{{F}}_{Y;\theta,\xi}$. Further, if every element
of $\mathcal{M}$ is absolutely continuous with respect to Lebesgue measure, we
can also write the statistical model as a family of conditional probability
density functions, i.e.,
$\mathcal{M}=\\{p_{Y|\theta,\xi}(y|\theta,\xi):\theta\in\bm{\Theta},\
\xi\in\Xi\\}$. We will make this simplifying assumption from here on, with the
understanding that the conditional density of $Y$ could be replaced by a
conditional probability measure $\mu(dy|\theta,\xi)$ as needed.

The parameters $\theta$ in the model are unknown and hence the object of
estimation or inference. On the other hand, $\xi$ can be controlled by the
experimenter. For instance, if candidate $Y$ corresponded to observations of a
spatiotemporal process, $\xi$ might represent the spatial and temporal
coordinates of a chosen set of observations. In a regression model, $\xi$
encodes the values of the covariates (i.e., independent variables) at which
observations are obtained. Different ways of representing $\xi$ will give rise
to different optimization problems, which we will discuss in Section 4.

The statistical model immediately yields a likelihood function $\theta\mapsto
p_{Y|\theta,\xi}(y|\theta,\xi)$, and the corresponding symmetric, positive
semi-definite $p\times p$ Fisher information matrix

$F(\theta,\xi)\coloneqq\mathbb{E}_{Y|\theta,\xi}[\nabla_{\theta}\log
p_{Y|\theta,\xi}(Y|\theta,\xi)\otimes\nabla_{\theta}\log
p_{Y|\theta,\xi}(Y|\theta,\xi)],$ (2)

which is a central object in estimation theory (Lehmann and Casella 1998).
Below we will discuss various design criteria based on the Fisher information
matrix.

Alternatively, one can take a Bayesian approach and endow the unknown
parameters, now denoted as $\Theta$, with a prior distribution. Let this
distribution have (Lebesgue) density $p_{\Theta}$ on $\mathbb{R}^{p}$. We will
always assume that the prior density is functionally independent of the design
$\xi$. The posterior density of $\Theta$ is then given by Bayes’ rule:

$p_{\Theta|Y,\xi}(\theta|y,\xi)=\frac{p_{Y|\Theta,\xi}(y|\theta,\xi)\,p_{\Theta}(\theta)}{p_{Y|\xi}(y|\xi)}.$
(3)

In the Bayesian paradigm, the prior and posterior distributions represent,
respectively, our states of knowledge before and after having observed $Y=y$.
The marginal density of the observations,

$\displaystyle p_{Y|\xi}(y|\xi)=\int
p_{Y|\Theta,\xi}(y|\theta,\xi)\,p_{\Theta}(\theta)\,\mathrm{d}\theta,$

appearing in the denominator of (3), is called the evidence or marginal
likelihood. We will also call this distribution the prior predictive, as it
reflects our probabilistic prediction of future values of $Y$ given only the
prior on $\Theta$, the statistical model, and a chosen design. Having observed
a particular value of the data, say $y^{\ast}$, at some chosen design $\xi$,
the posterior predictive density of the data $Y$ for a new design $\xi^{+}$ is

$\displaystyle p_{Y|\xi^{+},y^{\ast},\xi}(y|\xi^{+},y^{\ast},\xi)=\int
p_{Y|\Theta,\xi}(y|\theta,\xi^{+})\,p_{\Theta|Y,\xi}(\theta|y^{\ast},\xi)\,\mathrm{d}\theta.$

Many design criteria discussed below will explicitly take advantage of this
Bayesian formulation of the inference problem. Indeed, we will see that it is
useful to have the ability to incorporate prior information on $\Theta$—in
general, but especially in nonlinear design settings—and that the Bayesian
approach to OED is natural for decision-theoretic reasons as well.

###  2.1 Design criteria for the linear-Gaussian setting

A rich variety of design criteria, both Bayesian and non-Bayesian, have been
developed for linear-Gaussian models. This class of statistical models has
numerous practical applications: certainly linear regression, but also _linear
inverse problems_ , where observations depend _indirectly_ on the parameters
through the action of some linear operator. Canonical linear inverse problems
include deconvolution, computerized tomography, and source inversion, among
many others (Kaipio and Somersalo 2006). Design criteria in this setting are
often quite explicit and tractable, and also serve as a building block for
certain nonlinear design approaches.

We specify a general linear-Gaussian model as

$Y=G\theta+\mathcal{E},$ (4)

where $G\in\mathbb{R}^{n\times p}$ represents the linear ‘forward’ operator,
mapping parameters to data in $\mathbb{R}^{n}$, and $\mathcal{E}$ is a
Gaussian random variable with full-rank covariance matrix
$\Gamma_{Y|\Theta}\in\mathbb{R}^{n\times n}$ that does not depend on $\theta$.
We let $\mathcal{E}$ have mean zero; choosing otherwise would not affect the
developments below. In general, both $G$ and $\Gamma_{Y|\Theta}$ can depend on
the design $\xi$; that is, we have $G(\xi)$ and $\Gamma_{Y|\Theta}(\xi)$. The
linear-Gaussian model can be summarized as
$Y|\theta,\xi\,\sim\,\mathcal{N}\left(G(\xi)\theta,\Gamma_{Y|\Theta}(\xi)\right)$.

####  2.1.1 Classical alphabetic optimality

The Fisher information matrix associated with (4) is

$F(\theta,\xi)=F(\xi)=G(\xi)^{\top}\Gamma_{Y|\Theta}(\xi)^{-1}G(\xi).$ (5)

It is thus independent of the value of the parameters $\theta$; this property
of linear models greatly simplifies the construction of design criteria. Note
also that the inverse of the Fisher information matrix, $F^{-1}$, when it
exists, is the covariance of the least-squares estimate and hence the maximum
likelihood estimate, $\hat{\theta}(y)$, of $\theta$. Many classical design
criteria are therefore chosen to be scalar-valued functionals of the matrix
$F$. Indeed, we must somehow ‘scalarize’ $F$ to produce a useful optimization
objective, and various scalarizations encode different goals. We recall
several of these so-called ‘alphabetic optimality’ criteria as follows.

A-optimal design seeks

$\xi^{\ast}\in\arg\max_{\xi\in\Xi}\text{tr}\left(F(\xi)\right).$

When $F$ is invertible, which is invariably the situation of interest in
classical design and what we shall assume in the remainder of this subsection,
the optimization problem above is equivalent to
$\min_{\xi\in\Xi}\text{tr}\left(F^{-1}(\xi)\right)$, which can be interpreted
as minimizing the average variance of the $p$ components of $\hat{\theta}$.
D-optimal design, similarly, seeks

$\xi^{\ast}\in\arg\min_{\xi\in\Xi}\log\det\left(F^{-1}(\xi)\right),$

which minimizes the volume (in $\mathbb{R}^{p}$) of the smallest
$100(1-\alpha)\%$ confidence ellipsoid for $\theta$, for any confidence level
$1-\alpha$. A useful distinguishing feature of D-optimal designs is that they
are invariant under linear reparameterization (and hence rescaling) of
$\theta$; that is, if ${\theta}^{\prime}=M\theta$ for some invertible matrix
$M$, then a design that is D-optimal for $\theta$ is also D-optimal for
$\theta^{\prime}$. This is not true, in general, for other optimality
criteria.

While the A- and D-optimality criteria explicitly involve all the eigenvalues
of $F$, E-optimal designs minimize the maximum eigenvalue of $F^{-1}(\xi)$,
$\lambda_{\text{max}}(F^{-1}(\xi))$ (equivalently, maximize
$\lambda_{\text{min}}(F(\xi))=1/\lambda_{\text{max}}(F^{-1}(\xi))$). Such
designs thus minimize the variance of $q^{\top}\hat{\theta}$ among all
$q\in\mathbb{R}^{p}$ satisfying a norm constraint, that is, they minimize
$\max_{\|q\|=w}\mathbb{V}\text{ar}(q^{\top}\hat{\theta})=\max_{\|q\|=w}q^{\top}F(\xi)^{-1}q$
for any $w>0$.

A-, D-, and E-optimality are perhaps the essential building-block design
criteria for linear models, but there are numerous others. Some focus on the
estimation of a linear combination of the parameters $\theta$, or a subset of
the elements of $\theta$. Suppose, for example, that we are primarily
interested in $A^{\top}\theta$, where $A\in\mathbb{R}^{p\times s}$, $s<p$, and
$A$ has rank $s$. Then $\text{D}_{\text{A}}$-optimality generalizes
D-optimality by seeking

$\xi^{\ast}\in\arg\min_{\xi\in\Xi}\log\det\left(A^{\top}F^{-1}(\xi)A\right).$
(6)

This objective is justified by noting that $A^{\top}F^{-1}(\xi)A$ is the
covariance matrix of $A^{\top}\hat{\theta}$. If we put $A=[I_{s}\ \
0]^{\top}$, then the design criterion (6) focuses on the first $s$ elements of
$\theta$; this is called $\text{D}_{\text{S}}$-optimality.

An analogous generalization of A-optimality, for some matrix
$L\in\mathbb{R}^{p\times p}$, is called L-optimality (Atkinson et al. 2007):

$\xi^{\ast}\in\arg\min_{\xi\in\Xi}\text{tr}(F^{-1}(\xi)L)\,.$

If $L$ is symmetric and has rank $s\leq p$, then it can be expressed as
$L=AA^{\top}$ for $A\in\mathbb{R}^{p\times s}$. Then, using the cyclic
property of the trace, the L-optimality objective can be rewritten as
$\text{tr}(A^{\top}F^{-1}(\xi)A)$. If $A^{\top}$ is a row vector in this
setting (i.e., $s=1$), then the criterion seeks to minimize the variance of a
single linear combination of the parameters, and it is called c-optimality.

Other criteria instead seek to control the variance of predictions of the
linear model. Consider, specifically, a regression model on some compact
domain $\mathcal{X}$, where each row of $G$ is given by the evaluation of a
feature vector $f:\mathcal{X}\to\mathbb{R}^{p}$; that is, the $i$th row of $G$
is $G(i,:)=f^{\top}(x_{i})$ for some $x_{i}\in\mathcal{X}$ that is in the
support of the design $\xi$. The G-optimality criterion considers the variance
of the predicted response at any point $x\in\mathcal{X}$,
$f^{\top}(x)F^{-1}(\xi)f(x)$, and seeks a design $\xi$ that will minimize the
maximum value of this variance, that is,

$\xi^{\ast}\in\operatornamewithlimits{argmin}_{\xi\in\Xi}\max_{x\in\mathcal{X}}f^{\top}(x)F^{-1}(\xi)f(x).$
(7)

Variants of this criterion that target the average variance of the predicted
response over a region, rather than its maximum, are called I-optimality or
V-optimality. For a much more comprehensive discussion of classical alphabetic
optimality criteria and their properties, we refer to Hedayat (1981) and Shah
and Sinha (1989) as well as Atkinson et al. (2007, Chapter 10).

We should note here that the assumption of normality on $Y$ is generally not
needed for these criteria to apply, and for the statistical interpretations
given above to hold. Rather, we need only that $\mathbb{E}[Y]=G\theta$ and
$\mathbb{C}\text{ov}(Y)=\Gamma_{Y|\Theta}$. The best linear unbiased estimator
still follows from the least squares solution in this setting, and its
performance is bounded by the Fisher information matrix (cf. Cramér–Rao
bounds). In many situations (including the usual setting for G-optimality
described above), it is further assumed that $\Gamma_{Y|\Theta}=\sigma^{2}I$,
that is, the observational errors are uncorrelated and have constant variance.

####  2.1.2 Bayesian alphabetic optimality

In the Bayesian setting, the parameters $\theta$ of the linear model (4) are
endowed with a prior distribution. Since these parameters are now modeled as
random variables, we write them as uppercase $\Theta$ and require that
$\mathcal{E}$ and $\Theta$ be independent. Choosing a conjugate Gaussian
prior, $\Theta\sim\mathcal{N}(\mu_{\Theta},\Gamma_{\Theta})$, we obtain a
posterior distribution that is again Gaussian; it can be written in closed
form as
$\Theta|y,\xi\,\sim\mathcal{N}(\mu_{\Theta|Y}(y,\xi),\Gamma_{\Theta|Y}(\xi))$,
where

$\displaystyle\Gamma_{\Theta|Y}(\xi)$ $\displaystyle\coloneqq$
$\displaystyle\left(G(\xi)^{\top}\Gamma_{Y|\Theta}(\xi)^{-1}G(\xi)+\Gamma_{\Theta}^{-1}\right)^{-1}$
(8) $\displaystyle\mu_{\Theta|Y}(y,\xi)$ $\displaystyle\coloneqq$
$\displaystyle\Gamma_{\Theta|Y}(\xi)\left(G(\xi)^{\top}\Gamma_{Y|\Theta}(\xi)^{-1}y+\Gamma_{\Theta}^{-1}\mu_{\Theta}\right).$
(9)

The posterior mean therefore depends on the realized value of the data $y$,
but the posterior covariance matrix does not.

The Bayesian analogue to classical alphabetic optimality uses design criteria
that are functions of the _posterior covariance matrix_. Prior knowledge—and
more generally the ‘balance’ of information between the prior and the
likelihood, where the latter may be affected by the number of
observations—will therefore affect the optimal design, since the design
criteria depend on the dispersion (shape and scale) of the posterior.

For instance, Bayesian A-optimality seeks designs that minimize the trace of
the posterior covariance matrix,

$\xi^{\ast}\in\arg\min_{\xi\in\Xi}\text{tr}\left(\Gamma_{\Theta|Y}(\xi)\right)=\text{tr}\left((F(\xi)+\Gamma_{\Theta}^{-1})^{-1}\right).$

Note that, in contrast with classical A-optimality, this objective no longer
requires $F$ to be invertible, as long as the prior covariance
$\Gamma_{\Theta}$ is chosen to have full rank. The same is true of all other
Bayesian alphabetic optimality criteria, making these criteria well-suited to
designs with fewer than $p$ support points, and to ill-posed inverse problems
generally (Haber et al. 2008). Bayesian D-optimality seeks to minimize the
log-determinant of the posterior covariance matrix,

$\xi^{\ast}\in\arg\min_{\xi\in\Xi}\log\det\left(\Gamma_{\Theta|Y}(\xi)\right).$
(10)

Similarly, Bayesian $\text{D}_{\text{A}}$-optimality seeks to minimize the
log-determinant of the covariance of the posterior predictive distribution of
$A^{\top}\theta$, for some matrix $A\in\mathbb{R}^{p\times s}$ with rank
$s<p$; hence the goal is to minimize
$\log\det\left(A^{\top}\Gamma_{\Theta|Y}(\xi)A\right)$. See Attia et al.
(2018) for an application of this criterion, and of a Bayesian analogue of
classical L-optimality. Bayesian E-optimal design minimizes the maximum
eigenvalue of $\Gamma_{\Theta|Y}(\xi)$, and so on. For an extensive discussion
of Bayesian alphabetic optimality criteria and their interpretations, we refer
the reader to Chaloner and Verdinelli (1995) and DasGupta (1995). We will also
revisit several of these criteria from the more general viewpoint of decision
theory in Section 2.2; the decision-theoretic formulation lets us derive many
Bayesian alphabetic optimality criteria from specific utility functions and,
crucially, enables generalization to nonlinear models.

In the Bayesian setting, it is also natural to consider linear models with
unknown variance parameters, e.g., $\Gamma_{Y|\Theta}=\sigma^{2}I$ with
unknown $\sigma^{2}$, and to endow these variance parameters with suitable
priors. In general, this extension modifies the optimality criteria discussed
above—with some exceptions, such as Bayesian A-optimality using a conjugate
inverse Gamma prior for $\sigma^{2}$; for more information, see Chaloner and
Verdinelli (1995).

####  2.1.3 Designs as probability measures

So far, we have deliberately remained somewhat non-specific in describing how
the design $\xi$ enters the statistical models (1) or (4), other than to think
of $\xi$ as representing all the chosen ‘coordinates’ or locations of the
observations on some continuous domain $\mathcal{X}$, or the indices of
observations selected from some countable set of candidates. One rather
elegant way of formalizing these examples is to cast the design as a
_probability measure_ on some domain $\mathcal{X}$. This viewpoint,
originating in Kiefer and Wolfowitz (1959), is widely adopted in the classical
literature on optimal design.

To explain this perspective, let us first consider the discrete case, where
$\mathcal{X}$ is a countable and perhaps finite set of observation indices,
$\mathcal{X}=\\{1,2,3,\ldots\\}$. Suppose that we wish to choose $n$
observations in total. If $n_{i}$ observations are taken at each point
$i=1,2,3,\ldots$ and $\sum_{i}n_{i}=n$, then we can write $\xi_{i}=n_{i}/n$
and consider $(\xi_{i})_{i}$ to be a probability measure (Chaloner and
Verdinelli 1995). This class of ‘quantized’ measures, with integer $n_{i}$ and
hence weights that are multiples of $1/n$, is called an exact design. If each
point can only be observed once—i.e., if the selection is without
replacement—then we further require $n_{i}\in\\{0,1\\}$. It is often useful to
relax the integer constraint, such that $\xi$ is any probability measure over
the set of candidate indices; in this case, $\xi$ is called a continuous or
approximate design. Then $\Xi$ is the set of all probability mass functions
over $\mathcal{X}$.

Now consider the setting of continuous observation indices, on a compact set
$\mathcal{X}\subseteq\mathbb{R}^{d}$, for $d\geq 1$. This setting allows
candidate observations to be indexed by some continuous coordinates, e.g.,
angles for a tomography problem, real-valued spatial coordinates for a sensor
placement problem, or any other covariates in a generic regression problem.
The notion of a continuous _design_ extends naturally: $\xi$ is simply a
probability measure on $\mathcal{X}$, and $\Xi$ is the set of all such
probability measures. An exact design here would be a finite mixture of Dirac
measures with quantized weights—that is, a measure supported on a finite
collection of points in $\mathcal{X}$ with mixture weights that are multiples
of $1/n$, i.e., $\xi=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$, for
$x_{i}\in\mathcal{X}$.

A nice consequence of this general viewpoint is that many quantities relevant
to the preceding design criteria can be written as integrals with respect to
$\xi$. Consider a 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}$, $i=1,\ldots,n$, then the Fisher information matrix of
the model is

$\displaystyle
F(\xi)=\frac{1}{\sigma^{2}}\sum_{i=1}^{n}f(x_{i})f(x_{i})^{\top}.$ (11)

This expression follows from (5) by setting

$\displaystyle
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}$. If the design has continuous support,
then we simply have instead

$F(\xi)=\frac{n}{\sigma^{2}}\int_{\mathcal{X}}f(x)f(x)^{\top}\xi(\mathrm{d}x),$

where the $n$ above ensures that the scaling of the integral is consistent
with (11).

It is natural to wonder how to reconcile this continuous-support viewpoint
with the existence of classical optimal designs supported on a finite set of
points in $\mathcal{X}$ (Atwood 1969). In other words, when optimizing some
design criterion over the set of all probability measures on the infinite set
$\mathcal{X}$, why should a minimizer be supported only on a finite number of
locations? In fact, as explained in Atwood (1969) and Kiefer (1961b), under
conditions satisfied by any of the design criteria in Section 2.1.1, there
exists an optimal design supported on at most $p(p+1)/2$ points. Intuition for
the result follows from Carathéodory’s theorem, compactness of $\mathcal{X}$,
and the fact that $F$ is a $p\times p$ symmetric matrix: the optimal
information matrix $F$ can always be expressed as a convex combination of at
most $p(p+1)/2$ rank-one matrices, each produced by a single-point design. The
need for $F$ to be full rank, on the other hand, imposes a lower bound of $p$
on the number of points in an optimal design. In the case of Bayesian
alphabetic optimality criteria for linear models, similar upper bounds for the
number of support points in an optimal design have also been derived (Chaloner
1984). In the nonlinear setting, however, the Fisher information matrix
depends on the parameter $\theta$ (see Section 2.1.4). A common approach,
discussed below, is then to average a local design criterion over a
distribution on $\theta$. Because now (infinitely) many information matrices
$F(\theta,\xi)$ are involved, upper bounds on the number of support points do
not in general hold (Atkinson et al. 2007, Chaloner and Larntz 1989).

An important theme in classical design theory has been the construction of so-
called ‘equivalence theorems’ for continuous designs. The first such result
was due to Kiefer and Wolfowitz (1960), who showed that any continuous
D-optimal design is also G-optimal, and vice-versa. This result was
substantially generalized by Kiefer (1974) and Whittle (1973). The resulting
‘general equivalence theorem’ relies on the fact that the design criteria to
be minimized are convex functionals $\phi$ of the probability measure $\xi$.
Under this condition, with some further assumptions on the regularity of
$\phi$, the equivalence theorem provides multiple equivalent conditions for
the optimality of a design $\xi^{\ast}$, some of which correspond to verifying
that the directional derivatives of $\phi$ (with respect to feasible designs
$\xi$) are zero at $\xi^{\ast}$. The latter are useful to check optimality of
a given design measure (in the continuous case), and have been employed in
algorithms (Wynn 1972). Variants of the general equivalence theorem have been
established for Bayesian alphabetic optimality in linear models (Chaloner
1984, Pilz 1991) and for certain local optimality criteria averaged over the
prior in nonlinear models (Chaloner and Larntz 1989). There is a vast array of
results along this theme, which we will not attempt to survey here. Instead we
refer the reader to the comprehensive framework in Pukelsheim (2006), which
tackles design optimality for linear models using tools of convex analysis,
and the historical perspective in Wynn (1984).

Since our interest is largely in nonlinear problems (as well as infinite-
dimensional linear problems), we will generally resort to numerical methods
for optimizing $\xi$, which must in any case employ _tractable_ finite-
dimensional parameterizations of candidate designs. Moreover, _only an exact
design can be realized in an experiment_ , and hence some notion of ‘rounding’
is usually needed if a problem is initially solved from a continuous design
perspective. We will discuss these considerations further in Section 4.

####  2.1.4 Challenges of nonlinear design

In problems where dependence on the model parameters $\theta$ is nonlinear,
the Fisher information matrix $F$ will generally vary with $\theta$. Since
$\theta$ is unknown, it is then unclear where to evaluate $F(\theta,\xi)$ and
any of the associated design criteria.

A relatively crude approach is simply to choose some reference or ‘best-guess’
parameter value $\theta_{0}$ and proceed; this is known as ‘local’ design, but
it clearly ignores parameter uncertainty and the impact of nonlinearity, and
may have sharp dependence on the choice of $\theta_{0}$. One could iterate
this approach, by estimating $\theta$ after collecting data from a first local
design, and then using the estimated parameter value to produce a new local
design, and so on (Korkel et al. 1999). A more principled alternative is a
minimax formulation (Fedorov and Hackl 1997, King and Wong 2000, Berger and
Wong 2009). If $\phi(\theta,\xi)$ is any ‘local’ design criterion (containing
the parameter-dependent Fisher information matrix $F(\theta,\xi)$, for
instance, $\phi(\theta,\xi)=\log\det F^{-1}(\theta,\xi)$), then we seek

$\xi^{\ast}\in\operatornamewithlimits{argmin}_{\xi\in\Xi}\max_{\theta\in\bm{\Theta}}\phi(\theta,\xi).$

An interpretation of this objective is that it seeks the design with best
performance for the worst-case parameter value, i.e., the value of $\theta$
that is most challenging to estimate (in the sense of $\phi$). This idea has
been explored by Sun and Yeh (2007) and Siade et al. (2017), among others, but
generally leads to a rather difficult optimization problem. Typically $\phi$
is chosen to be a classical alphabetic optimality criterion (as in the example
of D-optimality above), and thus the only ‘prior’ information on $\theta$ in
these formulations is via the set $\bm{\Theta}$.

Another natural alternative is to introduce a prior distribution $p_{\Theta}$
for $\theta$ and to average any local design criterion $\phi(\theta,\xi)$ over
this prior. This approach is widely adopted, due in no small part to its
computational tractability; see Pronzato and Walter (1985) and Fedorov and
Hackl (1997). If $\phi$ is solely based on the Fisher information, however,
then such a formulation is only ‘pseudo-Bayesian’ (Atkinson et al. 2007, Ryan
et al. 2016). Indeed, in some such works, the prior is used as a means of
handling the parameter dependence arising from nonlinearity but then discarded
for subsequent analysis. Ryan et al. (2016) argue that any ‘fully Bayesian’
design criterion must be a functional of the posterior distribution.
Interestingly, however, Walker (2016) shows that prior expectation of the
trace of the Fisher information matrix,

$\displaystyle\int\text{tr}(F(\theta,\xi))p_{\Theta}(\theta)\,\mathrm{d}\theta,$

has an information-theoretic interpretation. Overstall (2022) and Prangle et
al. (2023) both explore using this quantity as a design criterion in nonlinear
problems, and show that it has a decision-theoretic justification as well.

As a step in the ‘fully Bayesian’ direction, given a nonlinear model and a
Gaussian prior $p_{\Theta}$ with full rank covariance matrix
$\Gamma_{\Theta}$, one could seek, as proposed in Chaloner and Verdinelli
(1995),

$\displaystyle\xi^{\ast}\in\arg\max_{\xi\in\Xi}\int\log\det\left(F(\theta,\xi)+\Gamma_{\Theta}^{-1}\right)p_{\Theta}(\theta)\,\mathrm{d}\theta.$
(12)

This can be loosely interpreted as minimizing the average, over the prior
predictive distribution of $Y$, of the log-determinant of the covariance of a
Gaussian approximation of each resulting posterior. Yet this interpretation is
rather imprecise, as for a general nonlinear model,
$F(\theta^{\prime},\xi)+\Gamma_{\Theta}^{-1}$ can be very far from the
precision matrix of the posterior distribution that results from a realization
of the data $y\sim p_{Y|\theta,\xi}(\cdot|\theta^{\prime},\xi)$. Criteria such
as this are best viewed as approximations of an expected utility function
arising from a more principled decision-theoretic formulation, which we
discuss next.

###  2.2 Decision- and information-theoretic formulations

The decision-theoretic approach to OED was formalized by Lindley (1956) (see
also Stone 1959, Raiffa and Schlaifer 1961) and has two primary ingredients: a
utility function $u$, chosen to reflect the purpose of the experiment; and the
Bayesian principle of averaging over what is uncertain (Berger 1985). In this
framework, any design criterion takes the form of an expected utility, to be
maximized with respect to $\xi$:111Beginning in this section, we will drop
subscripts from probability density functions when the arguments are explicit,
letting these arguments make the choice of density clear.

$\displaystyle U(\xi)=\mathbb{E}_{Y,\Theta|\xi}[u(\xi,Y,\Theta)]=\iint
p(y,\theta|\xi)\,u(\xi,y,\theta)\,\mathrm{d}\theta\,\mathrm{d}y.$ (13)

Here, the utility function $u(\xi,y,\theta)$ quantifies the value of an
experiment performed with a design $\xi$ that yields observations $y$, if the
true parameter value is $\theta$. Since the outcome of the experiment is not
known when selecting $\xi$, and since the parameter value is also uncertain,
we average over the joint prior distribution of $Y$ and $\Theta$. This process
yields the expected utility $U(\xi)$. Many choices of utility function $u$
have been proposed and explored in the literature. We review some of the
possibilities below.

####  2.2.1 Expected information gain in parameters

The influential paper of Lindley (1956) proposed using the expected gain in
Shannon information, from prior to posterior, as an optimal experimental
design criterion. It is evocative to think of this quantity in at least two
ways, namely

$\displaystyle U_{\text{KL}}(\xi)$
$\displaystyle=\mathbb{E}_{Y|\xi}\left[D_{\text{KL}}(p_{\Theta|Y,\xi}||p_{\Theta})\right]$
(14) $\displaystyle=H(\Theta)-H(\Theta|Y,\xi)$ (15)

where $D_{\text{KL}}$ denotes the Kullback–Leibler (KL) divergence, or
relative entropy, from prior to posterior,

$D_{\text{KL}}(p_{\Theta|y,\xi}||p_{\Theta})=\int
p(\theta|y,\xi)\log\frac{p(\theta|y,\xi)}{p(\theta)}\,\mathrm{d}\theta,$ (16)

and the two terms in (15) are the entropy and conditional entropy,
respectively:

$\displaystyle H(\Theta)$ $\displaystyle=-\int p(\theta)\log
p(\theta)\,\mathrm{d}\theta,$ (17) $\displaystyle H(\Theta|Y,\xi)$
$\displaystyle=-\iint p(y,\theta|\xi)\log
p(\theta|y,\xi)\,\mathrm{d}\theta\,\mathrm{d}y.$ (18)

Equality of the expressions (14) and (15) is easily verified. Note that
$U_{\text{KL}}$ is always non-negative: conditioning reduces entropy on
average (not necessarily for each realization $Y=y$, but when averaging over
values of $y$), with zero entropy reduction $H(\Theta)=H(\Theta|Y,\xi)$ if and
only if $\Theta$ and $Y|\xi$ are independent. Similarly the KL divergence,
whose expectation yields (14), is always non-negative and reaches zero if and
only if the two distributions being compared are identical (Cover and Thomas
2006). Some intuition for maximizing this criterion is that the design $\xi$
yielding data $Y$ that most reduce Shannon entropy is, in this information-
theoretic sense, the most informative. Another intuition is that an optimal
experiment should maximize the ‘change’ (here, quantified by the KL
divergence) from the prior to the posterior, averaged over the prior
predictive distribution of the data $Y$.

Lindley’s original rationale for the design criterion $U_{\text{KL}}$ was not
rooted in decision theory, but the criterion was later given a decision-
theoretic justification by Bernardo (1979); see also the discussion in Prangle
et al. (2023). Bernardo frames the task of inference as a decision problem,
where making a decision amounts to returning a probability density function
for the parameters of interest $\Theta$. He argues that the utility function
for this decision should be a _proper, local scoring rule_ (Gneiting and
Raftery 2007), and that these desiderata in turn dictate that $u$ must
specifically be a _logarithmic_ scoring rule. In the language of (13) above,
this means that one should choose

$\displaystyle u^{\text{score}}(\xi,y,\theta)=\log p(\theta|y,\xi)-\log
p(\theta).$ (19)

Substituting this utility into (13) immediately yields (14) and (15).

Note also that choosing the utility to be the KL divergence from prior to
posterior, which depends explicitly on $y$ and $\xi$ but not on $\theta$,

$\displaystyle
u^{\text{div}}(\xi,y,\theta)=D_{\text{KL}}(p_{\Theta|y,\xi}||p_{\Theta})=u^{\text{div}}(\xi,y)$
(20)

and substituting this utility into (13), yields the same expected utility
$U_{\text{KL}}$ (14). We also call $U_{\text{KL}}$ the expected information
gain (EIG) in $\Theta$, from prior to posterior. It is useful for subsequent
computations (see Section 3) to write the EIG more explicitly as follows:

$\displaystyle U_{\text{KL}}(\xi)$ $\displaystyle=\iint
p(y,\theta|\xi)\log\frac{p(\theta|y,\xi)}{p(\theta)}\,\mathrm{d}\theta\,\mathrm{d}y$
(21) $\displaystyle=\iint
p(y,\theta|\xi)\log\frac{p(y|\theta,\xi)}{p(y|\xi)}\,\mathrm{d}\theta\,\mathrm{d}y$
(22) $\displaystyle=\iint
p(y,\theta|\xi)\log\frac{p(y,\theta|\xi)}{p(y|\xi)p(\theta)}\,\mathrm{d}\theta\,\mathrm{d}y$
(23) $\displaystyle\eqqcolon\mathcal{I}(Y;\Theta|\xi).$

In all of these expressions, the joint prior predictive density of $Y$ and
$\Theta$ can be factored as $p(y,\theta|\xi)=p(y|\theta,\xi)p(\theta)$, i.e.,
as a product of likelihood and prior. Moving from (21) to (22) is an
application of Bayes’ rule (3). Both (21) and (22) make clear that some kind
of posterior calculation is necessary: the former involves the _normalized_
posterior density $p(\theta|y,\xi)$, while the latter involves the posterior
normalizing constant $p(y|\xi)$. Moreover, these expressions must be evaluated
for a range of $y$ values—i.e., for ‘all possible’ posterior distributions—via
the outer expectation. The last expression above, (23), shows that EIG is
equivalent to the _mutual information_ (MI) between the parameters and
observations given the design, $\mathcal{I}(Y;\Theta|\xi)$. Henceforth, we
will use the terms EIG and MI interchangeably.

Expanding (22) into two terms also shows that the EIG is a difference of
entropies of the data, paralleling (15),

$\displaystyle U_{\text{KL}}(\xi)$ $\displaystyle=H(Y|\xi)-H(Y|\Theta,\xi).$
(24)

For some statistical models, $H(Y|\Theta,\xi)$ is a constant function of
$\xi$. One example is the case of a nonlinear model with additive noise,
$Y=G(\Theta,\xi)+\mathcal{E}$, where $\mathcal{E}$ is independent of $\Theta$
and its distribution does not depend on $\xi$; then the entropy
$H(Y|\Theta=\theta,\xi)=H(\mathcal{E})$ and hence does not depend on $\theta$
or $\xi$. Maximizing EIG then specializes to maximizing the entropy of the
prior marginal distribution of $Y$. This design strategy is called ‘maximum
entropy sampling’; see Shewry and Wynn (1987) and Sebastiani and Wynn (2000).

The EIG objective has additional desirable properties. For one, it is
invariant under bijective transformations of $\theta$; this property follows
from the invariance of the KL divergence to such transformations, and thus
includes rescalings of the parameters as well as more complex
reparameterizations. We also emphasize that there are _no assumptions of
linearity or Gaussianity_ in the motivation for the EIG objective, and in any
of its expressions above.

In the linear-Gaussian case, however, maximizing EIG is equivalent to seeking
a Bayesian D-optimal design (10). To see this, note that when $\Theta$ and
$Y|\xi$ are jointly Gaussian (as in Section 2.1.2), the EIG or MI can be
written in closed form as

$\displaystyle\mathcal{I}(Y;\Theta|\xi)$
$\displaystyle=\frac{1}{2}\left(\log\det\Gamma_{\Theta}-\log\det\Gamma_{\Theta|Y}(\xi)\right)$
(25)
$\displaystyle=\frac{1}{2}\left(\log\det\Gamma_{Y}(\xi)-\log\det\Gamma_{Y|\Theta}(\xi)\right),$
(26)

where $\Gamma_{Y}$ is the marginal (prior predictive) covariance of $Y$,

$\Gamma_{Y}(\xi)=G(\xi)\Gamma_{\Theta}G(\xi)^{\top}+\Gamma_{Y|\Theta}(\xi).$

From (25), it is apparent that maximizing EIG with respect to $\xi$ is
equivalent to minimizing the log-determinant of the posterior covariance. If,
additionally, the observational error covariance is independent of the design,
then $\Gamma_{Y|\Theta}(\xi)=\Gamma_{Y|\Theta}$ and an equivalent goal, via
(26), is to maximize the log-determinant of the marginal covariance
$\Gamma_{Y}$; this is a specific (linear-Gaussian) case of the maximum entropy
sampling described above.

It is interesting to note that the Bayesian D-optimality criterion for linear-
Gaussian models can be derived from other utility functions, besides (19) and
(20); see Chaloner and Verdinelli (1995, Section 2.2) for details.

####  2.2.2 Other utility functions

Having defined an information- and decision-theoretic design criterion for
inference of the model parameters $\Theta$, it is natural to extend this
construction to other goals.

Suppose that we are interested in only a subset of the model parameters.
Partitioning $\Theta$ as $\Theta=(\Theta_{1},\Theta_{2})$, information gain in
$\Theta_{1}$ is captured by the KL divergence from its prior marginal
distribution to its posterior marginal distribution:

$\displaystyle u(\xi,y)$
$\displaystyle=D_{\text{KL}}(p_{\Theta_{1}|y,\xi}||p_{\Theta_{1}}),$ (27)

where

$p(\theta_{1})=\int p(\theta_{1},\theta_{2})\,\mathrm{d}\theta_{2}\ \ \
\text{and}\ \ \ p(\theta_{1}|y,\xi)=\int
p(\theta_{1},\theta_{2}|y,\xi)\,\mathrm{d}\theta_{2}.$

The outer expectation over $Y$ in (13), yielding the _expected_ information
gain in $\Theta_{1}$, must still account for uncertainty in both blocks of
$\Theta$. The optimal design criterion in this case becomes

$\displaystyle U(\xi)$ $\displaystyle=\iiint
p(y|\theta_{1},\theta_{2},\xi)p(\theta_{1},\theta_{2})\,\log\frac{p(\theta_{1}|y,\xi)}{p(\theta_{1})}\,\mathrm{d}\theta_{1}\,\mathrm{d}\theta_{2}\,\mathrm{d}y$
(28) $\displaystyle=\iiint
p(y|\theta_{1},\theta_{2},\xi)p(\theta_{1},\theta_{2})\,\log\frac{p(y|\theta_{1},\xi)}{p(y|\xi)}\,\mathrm{d}\theta_{1}\,\mathrm{d}\theta_{2}\,\mathrm{d}y.$
(29)

The second expression is analogous to (22) in that it uses densities for the
data $Y$, but now even the numerator of the density ratio involves
marginalization, as

$\displaystyle p(y|\theta_{1},\xi)$ $\displaystyle=\int
p(y|\theta_{1},\theta_{2},\xi)p(\theta_{2}|\theta_{1})\,\mathrm{d}\theta_{2}.$
(30)

Compared to the EIG in $\Theta$, evaluating this objective therefore requires
an additional integration over $\Theta_{2}$. Note also that (28) and (29) are
equivalent to the MI, $\mathcal{I}(Y;\Theta_{1}|\xi)$. In this formulation,
$\Theta_{2}$ could represent a variety of possible ‘nuisance’ parameters in
the statistical model, i.e., any parameters that are uncertain but simply not
the modeler’s immediate object of interest (Feng and Marzouk 2019,
Alexanderian et al. 2021). Special examples include the parameters of a
discrepancy model (Kennedy and O’Hagan 2001) designed to capture model error,
or the background medium in an inverse scattering problem (Borges and Biros
2018).

A generalization of the preceding formulation is to consider the EIG in some
(generally nonlinear) function of the parameters $\Theta$, $Z=\Psi(\Theta)$,
where $\Psi:\bm{\Theta}\to\mathbb{R}^{q}$ for some $q\leq p$. This can be
thought of as a ‘goal-oriented’ objective (just like Bayesian
$\text{D}_{\text{A}}$-optimality in the linear-Gaussian case) where the
function $\Psi$ encodes the true quantity of interest. Now we seek to maximize
the expected KL divergence from the prior predictive distribution of $Z$ to
its posterior predictive distribution:

$\displaystyle U(\xi)$
$\displaystyle=\mathbb{E}_{Y|\xi}\left[D_{\text{KL}}(p_{Z|Y,\xi}||p_{Z})\right]$
$\displaystyle=\iint
p(y,z|\xi)\,\log\frac{p(z|y,\xi)}{p(z)}\,\mathrm{d}z\,\mathrm{d}y$ (31)
$\displaystyle=\iint
p(y,z|\xi)\,\log\frac{p(y|z,\xi)}{p(y|\xi)}\,\mathrm{d}z\,\mathrm{d}y$ (32)
$\displaystyle=\mathcal{I}(Y;Z|\xi).$

Notably, this objective is also the original object of interest in Bernardo
(1979). While it is straightforward to write down, it raises significant
computational challenges, beyond those associated with calculating EIG in the
parameters alone. For generic $\Psi$, $p(\theta)$, and $p(y|\theta,\xi)$, we
do not have simple expressions for the prior density $p(z)$ or the posterior
density $p(z|y,\xi)$ of $Z$ (even up to a normalizing constant) appearing in
(31). Nor do we have easy access to the marginal likelihood $p(y|z,\xi)$ to
instead evaluate (32). Numerical approximations are needed, involving density
estimation or approximate Bayesian computation; see Section 3. Of course, in
specific cases (such as linear $\Psi$ with Gaussian priors) some aspects of
the expressions above become more tractable. We note also that if $\Psi$ is
bijective, then EIG in $Z$ is identical to EIG in $\Theta$, since (as noted
earlier) the information gain objective is invariant under bijective
transformations; otherwise the EIG in $Z$ is smaller (Bernardo 1979, Theorem
1).

The optimal design obtained by maximizing the EIG in some $Z$ can differ
drastically from the design maximizing EIG in $\Theta$. Figure 2.1 illustrates
these contrasts for sensor placement in a time-dependent advection-diffusion
problem (Zhong et al. 2024). The parameter $\Theta$ is the unknown source
location, endowed with a uniform prior on $[0,1]^{2}$. The source emits a
scalar quantity that diffuses and is advected towards the top-right of the
$[0,1]^{2}$ domain, with the advection velocity increasing linearly in time,
from a value of zero at $t=0$. The design entails placing a single sensor,
with coordinates $\xi=(\xi_{1},\xi_{2})\in[0,1]^{2}$, that measures the
concentration of the scalar at time $t_{1}>0$. Figure
2.1LABEL:sub@f:KL_nonGOOED shows a map of EIG in $\Theta$ as function of
sensor location (estimated numerically; see Section 3). We see that the
optimal measurement location is not unique, but lies roughly 0.2 units of
distance from the center of the domain; the slight asymmetry is due to the
direction of advection. Figure 2.1LABEL:sub@f:KL_GOOED, in contrast, shows
maps of EIG for four different choices of $Z$; each $Z$ is the predicted
concentration of the passive scalar at a future time $t_{2}>t_{1}$ and at the
specific location marked by a red star in each panel. We see that the optimal
sensor locations, maximizing these goal-oriented EIG criteria, are markedly
different from those in Figure 2.1LABEL:sub@f:KL_nonGOOED.

(a) EIG in $\Theta$

(b) EIG in various $Z$

Figure 2.1: Optimal sensor placement in a time-dependent advection-diffusion
problem. Each figure shows a map of expected information gain (EIG) in a
chosen quantity, as a function of the sensor location $\xi\in[0,1]^{2}$.
Measurements are made at time $t_{1}>0$, and advection is towards the top
right. (a) EIG in the unknown source location $\Theta$. (b) EIG for different
quantities of interest $Z$, where each $Z$ is the predicted concentration at
some future time $t_{2}>t_{1}$, at the location marked by the red star. The
optimal designs, maximizing EIG in each case, differ significantly.

If the quantity of interest $Z$ is random given $\theta$, e.g., if it is
described by a conditional density $p(z|\theta)$, then the formulation above
still applies. Computations may actually be easier: the problem of estimating
$p(z)$ or $p(z|y,\xi)$ is _smoothed_ by the kernel $p(z|\theta)$, and the
restriction that $q\leq p$ is lifted as long as the conditional density
$p(z|\theta)$ on $\mathbb{R}^{q}$ exists. A special case is when
$Z=Y|\xi^{+}$, that is, we wish to maximize information gain in the model
prediction for some given design $\xi^{+}$. We call this future prediction
$Y^{+}$, to distinguish it from the potential outcomes $Y$ of the experiment
being currently designed. In this case, the utility $u$ can be set to

$\displaystyle u(\xi,y)$
$\displaystyle=D_{\text{KL}}(p_{Y^{+}|y,\xi,\xi^{+}}||p_{Y^{+}|\xi^{+}}),$
(33)

which yields, as an expected utility, the expected information gain in
$Y^{+}$:

$\displaystyle U(\xi)$ $\displaystyle=\iint
p(y,y^{+}|\xi,\xi^{+})\,\log\frac{p(y^{+}|y,\xi,\xi^{+})}{p(y^{+}|\xi^{+})}\,\mathrm{d}y^{+}\,\mathrm{d}y$
(34) $\displaystyle=\iiint
p(y|\theta,\xi)p(y^{+}|\theta,\xi^{+})p(\theta)\,\log\frac{p(y^{+}|y,\xi,\xi^{+})}{p(y^{+}|\xi^{+})}\,\mathrm{d}y^{+}\,\mathrm{d}y\,\mathrm{d}\theta$
(35) $\displaystyle=\mathcal{I}(Y;Y^{+}|\xi,\xi^{+}).$ (36)

The penultimate line above reflects the fact that $Y$ and $Y^{+}$ are
conditionally independent given $\theta$, and hence their joint prior
predictive can be expanded and factored by introducing the parameter $\theta$
explicitly.

A final information-theoretic utility that we will consider arises from
problems of _model discrimination_ in the Bayesian setting (Myung and Pitt
2009). Suppose we have a countable set of models $\mathcal{M}_{m}$,
$m=1,2,\ldots$, each with its own parameters,
$\theta_{m}\in\bm{\Theta}_{m}\subseteq\mathbb{R}^{p_{m}}$, and its own prior
on parameters, $p(\theta_{m})$. Suppose there is also a (discrete) prior
distribution over the model indicators $m$. As a utility, we choose the
relative entropy from this prior to the posterior distribution over model
indicators,

$\displaystyle u(\xi,y)=\sum_{m}P(m|y,\xi)\log\frac{P(m|y,\xi)}{P(m)}.$ (37)

Following (13), we must take an expectation over the prior predictive
distribution of $Y$ to obtain an expected utility. Now, however, because there
are multiple possible models, the prior predictive distribution is itself a
mixture of the prior predictives of each model:

$\displaystyle
p(y|\xi)=\sum_{m}P(m)p(y|m,\xi)=\sum_{m}P(m)\int_{\bm{\Theta}_{m}}p(y|\theta_{m},\xi)p(\theta_{m})\,\mathrm{d}\theta_{m},$
(38)

where conditioning on $\theta_{m}$ also implies conditioning on $m$ at the
same time. The EIG in the model indicator $m$ follows by combining (37) and
(38): $U(\xi)=\int u(\xi,y)p(y|\xi)\,\mathrm{d}y$. As noted in Ryan et al.
(2016), this design approach applies in the $\mathcal{M}$-closed framework for
model selection (see Bernardo and Smith 2000, Chapter 6); that is, the true
(data-generating) model is assumed to be within the set of models considered,
and one must assign a prior weight $P(m)$ to the event that each model is
true. Examples of Bayesian OED for model discrimination can be found in Myung
and Pitt (2009), Cavagnaro et al. (2010), McGree et al. (2012), Drovandi et
al. (2014), Aggarwal et al. (2016), and Hainy et al. (2022).

While most of the preceding discussion has focused on information-theoretic
utilities, the decision-theoretic framework described at the start of Section
2.2 is certainly not limited to utility functions of this kind. As described
in Ryan et al. (2016) and Chaloner and Verdinelli (1995), another natural
utility is a quadratic loss, motivated by the desire to extract a point
estimate of $\theta$ from the posterior. For instance, let
$\bar{\theta}(y,\xi)\coloneqq\mathbb{E}[\Theta|y,\xi]$ denote the posterior
mean. Then we can write

$\displaystyle
u(\xi,y,\theta)=-\left(\theta-\bar{\theta}(y,\xi)\right)^{\top}B\left(\theta-\bar{\theta}(y,\xi)\right)$
(39)

for some symmetric positive semi-definite matrix $B\in\mathbb{R}^{p\times p}$.
The expected utility is then

$\displaystyle
U(\xi)=-\iint\left(\theta-\bar{\theta}(y,\xi)\right)^{\top}B\left(\theta-\bar{\theta}(y,\xi)\right)p(y,\theta|\xi)\,\mathrm{d}\theta\,\mathrm{d}y,$
(40)

which is the negative Bayes risk of the posterior mean under a _weighted
squared error_ loss. Maximizing the expected utility over $\xi$ thus minimizes
this risk. As noted in Chaloner and Verdinelli (1995), in the case of a
linear-Gaussian model, this formulation reverts to Bayesian A-optimal design
with weight matrix $B$, i.e.,
$\min_{\xi}\text{tr}\left(B\,\Gamma_{\Theta|Y}(\xi)\right)$.

Another family of non-information-theoretic utilities involves scalar
functionals of the posterior covariance matrix; these are not strictly
motivated by point estimation, but rather can be seen as a more
computationally tractable alternative to information-theoretic utilities that
require calculation of posterior normalizing constants. Ryan et al. (2016)
specifically suggest using as a utility the determinant of the posterior
precision matrix:

$\displaystyle
u(\xi,y)=\frac{1}{\det\left(\mathbb{C}\text{ov}(\Theta|y,\xi)\right)}$ (41)

and then, as usual, averaging this quantity over the prior predictive of $Y$
to obtain an expected utility:

$\displaystyle U(\xi)$
$\displaystyle=\int\left(\det(\mathbb{C}\text{ov}(\Theta|y,\xi))\right)^{-1}p(y|\xi)\,\mathrm{d}y$
(42)
$\displaystyle=\iint\left(\det(\mathbb{C}\text{ov}(\Theta|y,\xi))\right)^{-1}p(y|\theta,\xi)p(\theta)\,\mathrm{d}\theta\,\mathrm{d}y.$
(43)

We emphasize that this criterion is intended for nonlinear/non-Gaussian
problems, and thus calculation of the posterior covariance for different
realizations of $Y$ is not a computationally trivial undertaking. It is
instructive to compare this criterion to the similar but cruder heuristic
(12), which is motivated by a series of Gaussian approximations as described
in Chaloner and Verdinelli (1995).

To close this section, we point the reader to a more general formalism for
what comprises a ‘valid’ notion of information gain from a statistical
experiment, due to Ginebra (2007). In this formalism, an information measure
must satisfy a _minimal_ set of requirements: (i) it is real-valued; (ii) it
returns zero for a ‘totally non-informative experiment,’ where $Y$ is
independent of $\Theta$; and (iii) it satisfies sufficiency ordering
(Blackwell 1951, 1953, Le Cam 1964). The last requirement can be understood as
follows. Let $Y|\theta,\xi_{1}$ and $Y|\theta,\xi_{2}$ be the outcomes of two
different experiments, for the same parameter value $\theta$, and let $\eta$
be an independent random variable with fixed and known distribution,
introducing auxiliary randomness. If there exists a function $W$ such that
$W(Y,\eta)|\theta,\xi_{1}$ has the same distribution as $Y|\theta,\xi_{2}$ for
all $\theta$, then the experiment with design $\xi_{1}$ is said to be
‘sufficient for’ or ‘always at least as informative as’ the experiment with
design $\xi_{2}$. That is to say, the data from $\xi_{1}$ can generate data
from $\xi_{2}$ with an additional randomization mechanism and without knowing
$\theta$. In such a situation, $\xi_{1}$ is preferred. This generalized notion
of an information measure broadly encompasses several commonly used objectives
in OED, including mutual information. We refer readers to Ginebra (2007) for
an extended discussion, including connections to likelihood ratio and
posterior-to-prior ratio statistics.

###  2.3 Design criteria for infinite-dimensional problems

Infinite-dimensional statistical models arise in the Bayesian approach to
inverse problems (Stuart 2010, Dashti and Stuart 2017, Knapik et al. 2011)
and, more broadly, in nonparametric estimation and nonparametric Bayesian
procedures (Giné and Nickl 2021). These problems can be understood as
estimation or inference of functions; in other words, the parameter $\theta$
of the statistical model now belongs to a function space, rather than to
$\mathbb{R}^{p}$ for finite $p$. Application domains of such models are vast,
and we will not attempt to review them here. Instead we will focus on two
classes of problems where the integration of OED with infinite-dimensional
models has proved to be particularly fruitful.

####  2.3.1 Inverse problems in the Bayesian setting

The infinite-dimensional setting is natural for inverse problems involving
partial differential equations, where the parameter to be learned is typically
an initial condition, a source term, a boundary condition, or a heterogeneous
coefficient—and thus a function of space and/or time. An important research
theme, at the intersection of applied mathematics and statistics, has been to
create statistical formulations of inverse problems that are well-defined in
the infinite-dimensional setting (Stuart 2010). This is necessary, for
instance, to create consistent Bayesian models for inverse problems—Bayesian
models that have a well-defined limit as the discretization of the underlying
functions is refined (Bui-Thanh et al. 2013). Another important practical
result of these efforts is algorithms with discretization-invariant (and hence
dimension-independent) performance, for example Markov chain Monte Carlo
methods whose sampling efficiency does not deteriorate with grid refinement
(Hairer et al. 2011, Cotter et al. 2013, Cui et al. 2016, Rudolf and Sprungk
2018, Villa et al. 2021).

OED for infinite-dimensional Bayesian inverse problems is well explored in the
setting of Gaussian priors, particularly for linear inverse problems
(Alexanderian et al. 2014), and summarized in the recent review of
Alexanderian (2021). To explain these developments, we first briefly sketch
the setting and refer the reader to Stuart (2010) and Dashti and Stuart (2017)
for full details and precise results. The parameter $\Theta$ is modeled as a
random variable taking values in an infinite-dimensional separable Hilbert
space $(\mathcal{H},\langle\cdot,\cdot\rangle_{\mathcal{H}})$. $\Theta$ is
further assumed to be Gaussian, which means that the real scalar-valued random
variable $\langle v,\Theta\rangle_{\mathcal{H}}$ is Gaussian for any
$v\in\mathcal{H}$. The mean of $\Theta$ can be defined as an element
$\bar{\theta}$ of $\mathcal{H}$ satisfying
$\langle\bar{\theta},v\rangle_{\mathcal{H}}=\mathbb{E}[\langle\Theta,v\rangle_{\mathcal{H}}]$
for all $v\in\mathcal{H}$. The covariance operator of $\Theta$ is the
positive, self-adjoint, and compact linear operator
$C:\mathcal{H}\to\mathcal{H}$ defined through $\langle
v,Cw\rangle_{\mathcal{H}}=\mathbb{E}[\langle\Theta-\bar{\theta},v\rangle_{\mathcal{H}}\langle\Theta-\bar{\theta},w\rangle_{\mathcal{H}}]$
for all $v,w\in\mathcal{H}$. Because $\Theta$ takes values in $\mathcal{H}$,
the trace of $C$ is finite; it is then said that $C$ is of trace class. We can
write this Gaussian (prior) measure of $\Theta$ as
$\mu=\mathcal{N}(\bar{\theta},C)$, and will assume that $C$ is strictly
positive.

A common assumption on the statistical model for the observations
$Y\in\mathbb{R}^{n}$ is that they result from the action of a (possibly
nonlinear) ‘forward’ operator $G_{\xi}:\mathcal{H}\to\mathbb{R}^{n}$ perturbed
with additive Gaussian noise: $Y=G_{\xi}(\Theta)+\mathcal{E}$, where
$\mathcal{E}\sim\mathcal{N}(0,\Gamma_{Y|\Theta})$ is independent of $\Theta$.
For any fixed $Y=y$, this in turn defines a likelihood function that is
proportional to

$\mathcal{L}^{y}_{\xi}:\theta\mapsto-\frac{1}{2}\exp\Bigl{(}\bigl{(}y-G_{\xi}(\theta)\bigr{)}^{\top}\Gamma_{Y|\Theta}^{-1}\bigl{(}y-G_{\xi}(\theta)\bigr{)}\Bigr{)}.$

Under appropriate conditions on $G_{\xi}$ and the prior $\mu$, detailed in
Stuart (2010), the posterior distribution of $\Theta$, i.e., the distribution
of $\Theta$ given $Y=y$, denoted by $\mu^{y}_{\xi}$, is well defined and
dominated by the prior measure, $\mu^{y}_{\xi}\ll\mu$. One can then write
$\mu^{y}_{\xi}$ in terms of its Radon–Nikodym derivative with respect to
$\mu$:

$\displaystyle\frac{\mathrm{d}\mu^{y}_{\xi}}{\mathrm{d}\mu}(\theta)\propto\mathcal{L}^{y}_{\xi}(\theta).$
(44)

The KL divergence from prior to posterior,
$D_{\text{KL}}(\mu^{y}_{\xi}||\mu)$, can also be defined under these
conditions.

With this background in hand, we can summarize several design criteria that
have been proposed for infinite-dimensional Bayesian inverse problems. When
the forward operator $G$ is linear, the posterior measure $\mu^{y}_{\xi}$ is
again Gaussian, with a covariance operator $C_{\text{pos}}(\xi)$ that is
independent of $y$. In this setting, Alexanderian et al. (2014) propose
minimizing the trace of the posterior covariance operator with respect to
$\xi$: $\min_{\xi\in\Xi}\text{tr}(C_{\text{pos}}(\xi))$ . This is the
infinite-dimensional version of Bayesian A-optimality; the objective is well
defined because the posterior covariance $C_{\text{pos}}(\xi)$ is also of
trace class, under the conditions noted above. A marginalized version of
infinite-dimensional Bayesian A-optimality, focused on the covariance of a
subset of variables of interest, was used for design in Alexanderian et al.
(2021). This can be compared to L-optimality in the finite-dimensional
setting.

The analog of Bayesian D-optimality is somewhat less straightforward, as the
eigenvalues of the posterior covariance operator $C_{\text{pos}}(\xi)$
accumulate at zero, and hence minimizing the log-determinant of this operator
is not meaningful (Alexanderian 2021). Instead, Alexanderian et al. (2016a)
use the correspondence between D-optimality and maximizing EIG in linear
problems (see the discussion at the end of Section 2.2.1) to derive an
alternative objective for Bayesian D-optimality in the linear setting.
Specifically, they start with the EIG—taking advantage of the fact that the KL
divergence from prior to posterior is well defined, under conditions
summarized above. Specializing this quantity and its expectation over $Y$ to
the linear-Gaussian setting, the objective thus obtained is
$\frac{1}{2}\log\det(\text{Id}+\widetilde{H})$, where $\widetilde{H}$ is the
prior-preconditioned Hessian operator of the negative log-likelihood
(Alexanderian et al. 2016a, Theorem 1). This expression coincides with the
log-determinant of the posterior covariance in the finite-dimensional case.
Efficient ways to estimate this objective, leveraging low rank structure, are
discussed in Alexanderian and Saibaba (2018). We note also that
$\widetilde{H}$ is a central quantity in dimension reduction for Bayesian
inverse problems, and will appear again in Section 3.

For nonlinear forward operators, a full treatment of EIG in the infinite-
dimensional setting has not (to our knowledge) been used as a design
criterion. This may be largely due to computational tractability, though some
theoretical gaps (e.g., checking that the mutual information between $Y$ and
the infinite-dimensional $\Theta$ is well defined; see Duncan 1970) may
remain. Instead, researchers have focused on simpler design criteria and their
further approximations. For instance, Alexanderian (2021), motivated by the
finite-dimensional approach in Haber et al. (2009), discusses design that
minimizes the Bayes risk of the posterior mode $\hat{\theta}(y)$ under the
squared error loss defined by the inner product on $\mathcal{H}$, i.e.,
$\|\cdot\|^{2}\equiv\langle\cdot,\cdot\rangle_{\mathcal{H}}$. This is
analogous to (39) but using the posterior mode (also called the maximum a
posteriori (MAP) estimate (Dashti et al. 2013)) $\hat{\theta}(y)$ rather than
the posterior mean, as the former is typically more computationally tractable.
Another closely related heuristic suggested in Alexanderian (2021) is to
minimize the trace of the posterior covariance operator, in expectation over
the data $Y$:

$\displaystyle\min_{\xi\in\Xi}\int_{\mathcal{H}}\int_{\mathbb{R}^{n}}\text{tr}\left(C_{\text{pos}}(y,\xi)\right)\mathcal{L}^{y}_{\xi}(\theta)\,\mathrm{d}y\,\mu(\mathrm{d}\theta).$
(45)

This is essentially the ‘Bayesian A-posterior precision’ expected utility
described in Ryan et al. (2016, Section 3.1.2). As the posterior covariance
operator is difficult to approximate in nonlinear inverse problems (not to
mention many posterior covariances, one for each realization of the data used
to evaluate the integral in (45)), Alexanderian et al. (2016b) instead propose
replacing $\text{tr}(C_{\text{pos}}(y,\xi))$ above with the trace of the
inverse Hessian of a Laplace approximation of the posterior at the MAP point
$\hat{\theta}(y)$. This approximation is reasonable when the posterior is
‘close’ to Gaussian (Schillings et al. 2020, Helin and Kretschmann 2022,
Spokoiny 2023).

Computing any of these design criteria in the setting of infinite-dimensional
Bayesian inverse problems is a computationally challenging undertaking, due to
the high discretization dimension used to represent the parameter $\theta$ in
practice, as well as the cost of forward operator evaluations. Considerable
computational ingenuity is required; an essential step towards mitigating the
impact of high discretization dimension is to take advantage of low-rank
structure in the prior-preconditioned Hessian, and to exploit randomized
numerical linear algebra methods for computing eigendecompositions, estimating
the trace, and so on.

####  2.3.2 Gaussian process regression

Gaussian process (GP) regression, also known as kriging, is a ubiquitous tool
in spatial statistics, time series modeling, machine learning, engineering
design, surrogate modeling, and countless other applications. Rasmussen and
Williams (2006) and Gramacy (2020) provide excellent expositions of both
applications and some theoretical foundations. Experimental design for GP
regression has thus received considerable attention.

From the perspective of the previous section, GP regression can be viewed as a
linear Bayesian ‘inverse problem’ on function space, with a trivial forward
operator: a selection operator. The underlying true function,
$\theta^{\ast}:\mathcal{X}\to\mathbb{R}$, for
$\mathcal{X}\subseteq\mathbb{R}^{d}$, is observed directly through evaluation
at a finite collection of points $x_{i}\in\mathcal{X}$, $i=1,\ldots,n$,
perhaps with additive Gaussian noise, e.g.,
$Y_{i}=\theta^{*}(x_{i})+\mathcal{E}_{i}$ with
$\mathcal{E}_{i}\sim\mathcal{N}(0,\sigma^{2})$. Here we let the prior model
for the true function be a Gaussian process $\Theta$ on $\mathcal{X}$, with
mean function $m(x)\coloneqq\mathbb{E}[\Theta(x)]$ and positive semi-definite
covariance function
$c(x,x^{\prime})\coloneqq\mathbb{E}[(\Theta(x)-m(x))(\Theta(x^{\prime})-m(x^{\prime}))]$,
which defines the prior covariance operator $C$ via

$(Cv)(x):x\mapsto\int_{\mathcal{X}}c(x,x^{\prime})v(x^{\prime})\,\mathrm{d}x^{\prime},$

for functions $v\in L^{2}(\mathcal{X})$. Given a collection of observations
$y_{s}\coloneqq(y_{i})_{i=1}^{n}\in\mathbb{R}^{n}$, taken at corresponding
covariate values $x_{s}\coloneqq(x_{i})_{i=1}^{n}$, performing GP regression
entails conditioning $\Theta$ on these data. The posterior distribution,
describing this conditioned process, remains Gaussian,

$\displaystyle\Theta|(x_{s},y_{s})\sim\mathcal{N}(m_{\text{pos}},C_{\text{pos}}),$
(46)

with the posterior mean $m_{\text{pos}}$ and the posterior covariance function
$c_{\text{pos}}$ (yielding $C_{\text{pos}}$) expressible in closed form:

$m_{\text{pos}}(x)=m(x)+\alpha^{\top}c(x_{s},x)\ \ \text{and}\ \
c_{\text{pos}}(x,x^{\prime};x_{s})=c(x,x^{\prime})-c(x_{s},x)^{\top}Rc(x_{s},x^{\prime}),$

where the matrix $R^{-1}\in\mathbb{R}^{n\times n}$ has entries
$[R^{-1}]_{ij}=c(x_{i},x_{j})+\delta_{ij}\sigma^{2}$, the coefficient vector
$\alpha\in\mathbb{R}^{n}$ has entries $\alpha_{i}=R{[i,:]}(y_{s}-m(x_{s}))$,
and
$c(x_{s},x)\coloneqq\left(c(x_{1},x),\ldots,c(x_{n},x)\right):\mathcal{X}\to\mathbb{R}^{n}$
(Rasmussen and Williams 2006, Gramacy 2020).

Since the purpose of GP regression is generally to make predictions about
$\theta^{\ast}$ at unseen values of the covariates $x$, most design criteria
involve the (posterior) predictive variance, i.e., the variance of
$\Theta|(x_{s},y_{s})$. In this setting, however, the ‘parameter’ is the
process $\Theta$, and hence the boundary between parameters and predictions is
rather blurred.

Another perspective on GP regression follows intuitively from finite-
dimensional distributions of the process $\Theta$, which are always
multivariate Gaussian (both before and after conditioning on the data). For a
finite number of sites
$x_{\mathcal{S}}\coloneqq(x_{i})_{i=1}^{m}\in\mathcal{X}$,
$\Theta(x_{\mathcal{S}})\coloneqq(\Theta(x_{i}))_{i\in\mathcal{S}}$ is simply
a multivariate normal random vector. We can observe some components of this
vector, and we wish to use these observations to predict other components.

Maximum entropy sampling (Shewry and Wynn 1987) originates with this
discretized perspective. Let $x_{s}\subset x_{\mathcal{S}}$ denote the $n<m$
distinct locations selected for a candidate design and let
$x_{s^{c}}=x_{\mathcal{S}}\setminus x_{s}$ denote its complement. Then the
chain rule for entropy yields

$\displaystyle
H\left(\Theta(x_{\mathcal{S}})\right)=H\left(\Theta(x_{s})\right)+H\left(\Theta(x_{s^{c}})|\Theta(x_{s})\right).$
(47)

Since the left-hand side of (47) is fixed, minimizing entropy in the
predictions at unobserved sites $x_{s^{c}}$ given the observations (the second
term on the right) can be accomplished by maximizing the first term on the
right. Thus finding an optimal design is cast as maximizing the entropy of the
model predictions (the _joint_ entropy of these predictions) at the observed
locations, $H\left(\Theta(x_{s})\right)$. (Recall that we also discussed
maximum entropy sampling for general parametric statistical models in Section
2.2.1.) More explicitly, the optimization problem is typically posed with some
cardinality constraint on the number of observations, e.g., $|x_{s}|\leq n$:

$\displaystyle\operatornamewithlimits{argmax}_{x_{s}\subset{x_{\mathcal{S}}},\
|x_{s}|\leq n}H\left(\Theta(x_{s})\right).$ (48)

The objective above can be understood as a set function, i.e., a function of
all subsets of $x_{\mathcal{S}}$. In the present case, since $\Theta(x_{s})$
is a Gaussian vector, closed-form expressions for the entropy are immediately
available. The problem is equivalent to finding the principal submatrix of
$\mathbb{C}\text{ov}(\Theta(x_{\mathcal{S}}))$ that has largest determinant.
This problem is NP-hard, but many practical algorithms have been developed to
tackle it (Ko et al. 1995).

A crude approximation to maximum entropy sampling is to choose the elements of
$x_{s}$ one at a time, in a greedy fashion: beginning with $x_{s}=\emptyset$
and $x_{s^{c}}=x_{\mathcal{S}}$, at each iteration select from $x_{s^{c}}$ the
point with maximum predictive variance,222Recall that the entropy of a
univariate Gaussian random variable is an increasing function of its variance.
by ranking the diagonal elements of
$\mathbb{C}\text{ov}(\Theta(x_{s^{c}})|\Theta(x_{s}))$. Then add this point to
$x_{s}$ and repeat. Seo et al. (2000) and subsequent papers call this approach
‘active learning MacKay’, after MacKay (1992). Its performance can be far from
optimal, however, as the entropy objective is not submodular (see Section 4).

An alternative design approach, advocated by Krause et al. (2008) (see also
Caselton and Zidek (1984)) is to maximize MI, rather than entropy.
Specifically, the problem is posed as

$\displaystyle\operatornamewithlimits{argmax}_{x_{s}\subset x_{\mathcal{S}},\
|x_{s}|\leq n}\mathcal{I}\left(\Theta(x_{s});\Theta(x_{s^{c}})\right).$ (49)

Greedy approaches are typically applied to this problem, for reasons of
computational tractability. Section 4 will discuss these algorithmic
considerations in much more detail. Here, however, we will note that greedy
approaches tend to work far better for MI maximization than for entropy
maximization, as MI is submodular (Krause et al. 2008, Nemhauser et al. 1978,
Fisher et al. 1978). Krause et al. (2008, Section 4.1) provide some useful
intuition contrasting greedy selection via MI and greedy selection via
predictive entropy. A key consideration in the setup above is to ensure also
that the objective is monotone increasing for $|x_{s}|\leq n$, which is
required for optimization guarantees to hold. Krause et al. (2008) shows that
$\mathcal{I}\left(\Theta(x_{s});\Theta(x_{s^{c}})\right)$ is approximately
monotone in this regime as long as the discretization of the underlying domain
$\mathcal{X}$, via $x_{\mathcal{S}}$, is sufficiently fine. Beck and Guillas
(2016) present improvements to greedy MI maximization tailored to computer
model emulation.

A rather different class of design approaches is more rooted in the continuous
view of GPs, seeking the observation locations that minimize the resulting
posterior predictive variance, integrated over the domain of the process,
$\mathcal{X}$. Letting $x_{s}\subset\mathcal{X}$ denote a finite collection of
observation locations (not necessarily chosen from a finite candidate set),
the objective to be minimized, over feasible $x_{s}$, can be written as

$\displaystyle\int_{\mathcal{X}}c_{\text{pos}}(x,x;x_{s})\,\mathrm{d}x.$ (50)

This objective has appeared in many papers (Sacks et al. 1989, Seo et al.
2000, Santner et al. 2018, Gorodetsky and Marzouk 2016) and has been variously
called the integrated mean-squared error (IMSE) criterion, the integrated
mean-squared prediction error (IMSPE) criterion, or the integrated variance
(IVAR) criterion. With discrete candidate sets and a greedy one-point-at-a-
time approach to constructing $x_{s}$ (see below), it is also called ‘active
learning Cohn’ (ALC), after Cohn et al. (1996). In the language of classical
design, minimizing (50) can be understood as a kind of (Bayesian)
V-optimality, in that one minimizes the predictive variance integrated over a
region. We should also note that (50) is precisely
$\text{tr}(C_{\text{pos}})$, i.e., the trace of the posterior covariance
operator, which is the infinite-dimensional notion of A-optimality discussed
in Section 2.3.1. In practice, the integral (50) is approximated by a large
set of points chosen uniformly over $\mathcal{X}$, or perhaps non-uniformly to
reflect some desired weight. Both discrete selection methods (see (Seo et al.
2000)) and methods that optimize over the continuous coordinates of $n$ points
$x_{s}=(x_{1},\ldots,x_{n})$ have been explored in the literature. For the
latter, see Sacks et al. (1989) and Gorodetsky and Marzouk (2016). Here, as
with maximum entropy sampling, one can also optimize for all $n$ elements of
$x_{s}$ simultaneously (a ‘full batch’ design procedure), or proceed in a
greedier but sub-optimal fashion: select a subset of design points to minimize
(50), ‘freeze’ these points and update $c_{\text{pos}}$ accordingly, and
repeat for the next subset. Batch approaches are more computationally
demanding, but generally yield better performance; see demonstrations in
Gorodetsky and Marzouk (2016). We will discuss related optimization issues
further in Section 4.

Finally, we mention several approaches that rely on spectral decompositions of
the GP, specifically the Karhunen–Loève representation of
$\Theta\sim\mathcal{N}(m,C)$ (Karhunen 1947, Loève 1948). The idea is to find
the leading eigenvalues and eigenfunctions $(\lambda_{i},\phi_{i})$ of the
prior covariance operator $C$, and to write the GP as

$\displaystyle\Theta(x)=m(x)+\sum_{i=1}^{\infty}\sqrt{\lambda_{i}}\phi_{i}(x)\zeta_{i},$
(51)

where the scalar-valued random variables $\zeta_{i}$ are standard Gaussian and
mutually independent. If the eigensystem is truncated to $r<\infty$ eigenpairs
$(\lambda_{i},\phi_{i})_{i=1}^{r}$, then GP regression is reduced to
parametric regression with coefficients $(\zeta_{1},\ldots,\zeta_{r})$. Then
any standard Bayesian alphabetic optimality criterion can be applied; see
Fedorov and Flanagan (1997), Fedorov and Müller (2007), and Harari and
Steinberg (2014). See also Spöck (2012) for a related approach based on the
polar spectral representation of $\Theta$, assuming that the process is
stationary and isotropic. The truncated eigendecomposition can also be used to
approximate the integrated posterior variance objective (50); see Fedorov
(1996).

One outstanding issue in model-based design for GP regression is that
hyperparameters of the prior covariance function (controlling, e.g., the scale
of the prior variance, correlation lengths, and smoothness) must often be
learned from data as well. The methods discussed above all take the prior
covariance as fixed, and thus ignore the ‘outer’ hyperparameter learning
process—as well as the impact that uncertainty in the hyperparameters has on
the predictive distribution. In some applications, moreover, the main interest
is not in prediction of an unknown function, but rather in learning the
parameters $\gamma$ of the covariance function itself (Pardo-Igúzquiza 1998);
it is natural to expect that optimal designs for this purpose should differ
from optimal design for prediction. A common way of learning $\gamma$ is an
‘empirical Bayes’ approach, which provides a point estimate of the covariance
parameters by maximizing the log-marginal likelihood,
$\hat{\gamma}\in\operatornamewithlimits{argmax}_{\gamma}\log
p(y_{s}|x_{s},\gamma)$, where

$\displaystyle\log
p(y_{s}|x_{s},\gamma)=-\frac{1}{2}y_{s}^{\top}R_{\gamma}y_{s}-\frac{1}{2}\log\det
R_{\gamma}^{-1}-\frac{n}{2}\log 2\pi,$

and the notation $R_{\gamma}$ emphasizes that $\gamma$ controls the matrix
$R\in\mathbb{R}^{n\times n}$ defined earlier via
$[R^{-1}]_{ij}=c(x_{i},x_{j})+\delta_{ij}\sigma^{2}$. Here, $\gamma$ may
include parameters of the covariance function $c$ as well as the noise
variance $\sigma^{2}$. Solving this optimization problem is easier,
computationally, than treating $\gamma$ in a fully Bayesian way. That said,
there are many papers and software packages (Gramacy 2022) that do the latter,
endowing $\gamma$ with a prior distribution and inferring it jointly with
$\Theta$; usually, this task requires Markov chain Monte Carlo (MCMC) or
sequential Monte Carlo methods, as the joint distribution of $\Theta$ and
$\gamma$ is non-Gaussian. These fully Bayesian approaches are therefore
computationally demanding—even more so in the setting of design. In principle,
one could use the joint posterior distribution of $(\Theta,\gamma)$ to find
designs that maximize EIG in $\Theta$, $\gamma$, or both, following the
criteria developed in Sections 2.2.1–2.2.2. We are not aware of methods that
completely realize this approach. Instead, a variety of more practical
computational schemes have been devised for experimental design in GP
regression when the prior covariance parameters are uncertain.

Zhu and Stein (2005) focus on design for parameters of the covariance function
only, using the Fisher information matrix derived from the marginal
likelihood. Since dependence on these parameters $\gamma$ is generally
nonlinear, the authors propose using either ‘local’ D-optimal design, a
minimax approach, or a prior-averaged D-optimality criterion (cf. Section
2.1.4). It is interesting to note that the resulting designs involve points
that are non-uniformly spaced over $\mathcal{X}$; intuitively, such point sets
are useful for learning correlation lengths in the prior covariance function
(Gramacy 2020). Zhu and Stein (2006) then suggest a two-stage design process,
where some fraction of the design points are chosen according to a criterion
focused on estimation of the covariance parameters $\gamma$, while the rest
are chosen to improve prediction of $\Theta$; uncertainty in the covariance
parameters, given some asymptotic approximations, is accounted for in the
design criterion for the latter. Spöck and Pilz (2010) cast design for
prediction, using a spectral decomposition of the GP, within a minimax
formulation over a compact set of parameterized covariance functions. On the
other hand, the local, sequential design scheme of Gramacy and Apley (2015)
interleaves design for prediction with local Fisher-infomation matrix-based
design for the length-scale parameter in $\gamma$. And simpler sequential
schemes, where batches of standard (e.g., IMSE) design for $\Theta$ are
interleaved with maximum likelihood estimation of covariance parameters, are
pursued in Harari and Steinberg (2014) and Gorodetsky and Marzouk (2016).
Further perspective on such sequential schemes is given in Gramacy (2020,
Section 6.2). Indeed, it is very natural to interleave updates of the
covariance function parameters with actual observations, in a sequential
design fashion. In the purely Gaussian case (i.e., with a fixed covariance
function), the realized values of $y_{s}$ would have _no_ impact on subsequent
designs, since the posterior covariance and entropy depend only on $x_{s}$,
not on $y_{s}$. When parameters of the covariance function must also be
inferred, however, we are in the nonlinear design setting and hence there is
value to feedback: $y_{s}$ informs the covariance parameters, and in turn
these parameters reshape the predictive uncertainty of the Gaussian process
for the next stage.

We will discuss closed-loop sequential design much more systematically in
Section 5. Here we will mention just a few more instantiations in the setting
of GP regression. Riis et al. (2022) perform myopic sequential design using
criteria based on the marginal posterior of $\Theta$, where marginalization
over the kernel hyperparameters is performed with MCMC samples. Hoang et al.
(2014) develop a non-myopic sequential design policy that can be understood as
approximating a sequential variant of maximum entropy sampling; the authors
argue that this policy naturally balances effort between informing covariance
hyperparameters and directly reducing uncertainty in the prediction of
$\Theta$ itself.

###  2.4 Related problems and their distinctions

To help orient the reader, here we discuss several classes of problems in the
broader literature that have some conceptual overlap with optimal experimental
design, but also some essential differences.

##### Space-filling and other non-model based designs.

Space-filling designs, as the name suggests, spread design points throughout
the domain $\mathcal{X}$ so that one can reasonably assess variations of a
generic response or the parameters of an associated statistical model. In
their simplest form, these methods do not attempt to exploit the structure of
a statistical model for the response, or make any assumptions on such a model;
they are hence _not_ model-based designs, in contrast to the focus of this
article. Spread in the design space is achieved by formulating an optimization
problem, cases of which are primarily distinguished by whether the focus is
exclusively on the distance among the design points in $\xi$ (maximin distance
design), or the distance to all points in the ground set $\mathcal{X}$
(minimax distance design) (Johnson et al. 1990). These notions have clear
analogies to the well-known problem of sphere packing (Zong 1999). Pure space-
filling designs tend to have poor projection properties, failing to retain
their optimality properties when viewed in subspaces. This is undesirable in
circumstances when the response is insensitive to one or more of the design
variables. In such cases, Latin hypercube designs (McKay et al. 1979), which
ensure that any of their projections along a single coordinate axis yields a
maximin distance design, are a suitable alternative. Space-filling properties
in larger subspaces can be induced through orthogonal array extensions of
Latin hypercube design (Owen 1992, Tang 1993). Further, ‘maximum projection’
(MaxPro) designs have been developed to achieve space-filling properties on
_all_ possible subsets of factors (Joseph et al. 2015, 2020).

Other approaches using entropy maximization have also been suggested for
space-filling design (Jourdan and Franco 2010). Note that this is not akin to
the maximum entropy sampling methods we previously discussed in Section 2.3.2.
Here, ‘entropy’ is that of the empirical distribution of design points on
$\mathcal{X}$, and is used as a design criterion to be maximized. The core
idea is to relate the space-filling quality of the design to the uniform
distribution on $\mathcal{X}$, motivated by the fact that the uniform
distribution has maximum entropy among all distributions with prescribed
finite support. A similar idea is pursued through designs that seek to
minimize the discrepancy (Niederreiter 1992) of a set of design points. This
notion relates to the much broader topic of low-discrepancy sequences and
quasi-Monte Carlo methods for integration (Caflisch 1998, Dick et al. 2013).
We refer the reader also to Pronzato and Müller (2012) and Santner et al.
(2018) for a more comprehensive discussion of space-filling designs targeting
computer experiments.

Other standard non-model-based design strategies include factorial designs,
with blocking and fractional variants, as well as composite designs; see Wu
and Hamada (2011) and Atkinson et al. (2007, Chapter 7) for more.

##### Active learning.

Active learning is a term originating in the computer science and machine
learning communities, referring to a diverse array of algorithms for choosing
which data to ‘label,’ usually in a supervised learning setting (Dasgupta
2011). In statistical terms, we can understand this setting as regression with
real or discrete-valued outcomes $y_{i}$ (where the latter case is
classification). In so-called ‘pool-based active learning,’ there is a large
pool of candidate covariates or feature values $\\{x_{i}\\}_{i}$, referred to
as the unlabeled data (Schein and Ungar 2007). To label a chosen data point is
to obtain its associated outcome, thus creating the pair
$(x_{i^{\ast}},y_{i^{\ast}})$ for some chosen index $i^{\ast}$. We can thus
think of this problem as OED with a countable and even finite design space
$\mathcal{X}$, corresponding to which indices should be chosen from the
unlabeled pool. Other learning scenarios might select covariates $x_{i}$ from
an infinite set (e.g., with $\mathcal{X}$ now a region of $\mathbb{R}^{d}$);
alternatively, one might be presented with a stream of successive $x_{i}$ and
be required to choose, on-the-fly, whether to label the current value (Settles
2009). In any of these cases, active learning usually refers to a _sequential_
version of the OED problem where data to be labeled are selected one at a time
or in a batch, then labeled, then used to update a model, and then the process
repeats. Most often, these iterations take a _greedy_ approach (see Section
5).

An important point of differentiation among active learning methods is _by
what mechanism_ this selection occurs. One class of selection methods, known
as uncertainty sampling, selects the unlabeled data point(s) for which the
model’s current predictions are most uncertain. The notions of uncertainty
used here vary widely, and can include many heuristics that do not have a
statistical justification. This is an important distinction from OED. If the
notion of uncertainty is tied to the posterior predictive distribution of a
Bayesian model, however, then we can recover maximum entropy sampling,
discussed in Section 2.2.1 and 2.3.2 above. Indeed, predictive entropy is
commonly used to rank candidate points in active learning (Lewis 1995).
Another widely used class of selection methods is ‘query-by-committee’ (Freund
et al. 1997), where disagreement among an ensemble of models is used to rank
candidate points, such that points with greater disagreement are chosen. Yet
another class of criteria ranks candidate points by how much their label would
reduce uncertainty in the predictions of the model being trained. An example
of the latter, where uncertainty is captured by integrated variance, is the
ALC approach discussed in Section 2.3.2. Other selection methods focus on more
tailored prediction goals; for instance, Blanchard and Sapsis (2021) describe
selection criteria favoring regions of the input space that yield unusual
output values, to help build regression models capable of predicting extreme
events.

We should also emphasize that many active learning methods are not based on
explicit design criteria, but rather on other heuristics (Settles 2009).
Moreover, even criterion-based active learning methods usually focus on
reducing uncertainty in predictions, rather than on improving estimation or
inference of the parameters of a statistical model. Again, this is an
important distinction from many OED problems.

##### Bayesian optimization.

Bayesian optimization (BO) (Močkus 1975, Jones et al. 1998, Wang et al. 2023)
is widely used in applications from engineering design to machine learning, to
name just a few. It is essentially a derivative-free optimization method, used
to maximize ‘black-box’ (and often computationally expensive) objective
functions—that is, functions which can be evaluated pointwise, but whose
derivatives cannot be directly evaluated. Gaussian process regression is a key
ingredient of modern BO. GP regression is used to build an approximation of
the objective function (via the mean of the GP) and an estimate of uncertainty
in the predicted value of this objective (via the variance of the GP), and
both are refined over the course of the optimizaton iterations. A design-type
question then arises in choosing where (i.e., at what points in the input
domain $\mathcal{X}$) to evaluate the objective. In BO, this question is
typically resolved by defining a real-valued, easy-to-evaluate ‘acquisition
function’ over $\mathcal{X}$ and finding its maxima. A point at which the
acquisition function is maximized is then taken to be the next evaluation
point for the objective. The realized value of the objective at each iteration
affects the acquisition function at the next stage, and hence the design
process is sequential and adaptive. The acquisition process is generally
formulated in a myopic way, though there are a few exceptions (Lam and Willcox
2017, Wu and Frazier 2019).

Many acquisition functions have been proposed in the literature, beginning
with Jones et al. (1998) and in the decades since; these functions generally
balance a notion of ‘exploration’ (learning the objective in unseen places,
where the GP model has large predictive variance) with ‘exploitation’
(evaluating the objective at points where it is expected to be larger than the
current best value). A prominent choice is the ‘expected improvement’ function
and its many variants (Močkus 1975, Zhan and Xing 2020). Yet many other
choices are possible, with batch/parallel (Chevalier and Ginsbourger 2013,
Wang et al. 2020) and even multi-fidelity (Song et al. 2019) schemes. BO is a
large and active field which we cannot hope to survey here; instead we point
the reader to a few recent reviews and tutorials (Shahriari et al. 2016,
Frazier 2018, Wang et al. 2023). To set it in context relative to OED,
however, we emphasize that BO and OED are essentially different: the goal of
BO is to _maximize_ a function, not to predict the value of the function over
all of $\mathcal{X}$ or even some a priori chosen subset of $\mathcal{X}$. The
resulting sets of evaluation points thus differ, both in their configuration
and in their purpose, from the GP regression designs discussed in Section
2.3.2.

##### Data summarization.

Reducing or somehow ‘summarizing’ large data sets is a problem of frequent
practical interest, motivated by considerations of both storage and
computation. The cost of most Bayesian inference algorithms, for example,
scales at least linearly with the size of the data. Random subsampling of data
is of course an option, but a more effective approach is to identify a small
weighted subset of the data that is somehow representative of the full
dataset. This is the notion of a coreset, which originated in computational
geometry and computer science (Agarwal et al. 2005, Feldman and Langberg
2011). It has since been formulated in a statistical setting, and notably the
Bayesian setting (Huggins et al. 2016, Campbell and Broderick 2018, 2019,
Campbell and Beronov 2019), where the idea is to find weighted data subset of
given size that least changes the likelihoood or the posterior distribution
from its original full-data version. Ostensibly this problem seems similar to
OED, but a crucial difference is that coresets are generally identified after
the data $Y$ are realized, and depend on the realized values of data. In
optimal design, on the other hand, a design must be chosen before $Y|\xi$ are
observed.

##  3 Numerical approximation of design criteria

Now we turn to one of the central computational questions of optimal design:
how to approximate, numerically, the value of a chosen design criterion at any
candidate design? This issue is not particularly vexing for standard
alphabetic optimality criteria and linear models, where closed-form
expressions are generally available. Evaluating these expressions can become
costly when parameters $\Theta$ are high- or infinite-dimensional, however,
and we discuss dimension reduction methods relevant to this setting in Section
3.4. The information-theoretic design criteria introduced in Section 2.2, on
the other hand—which have become a mainstay of modern OED due to their
flexibility and their applicability to complex nonlinear models—can be very
challenging to evaluate, even when the parameter dimension is low. This
section will discuss a variety of computational approaches for approximating
such objectives, resting on nested Monte Carlo estimation (Section 3.1),
approximation of the relevant densities within tractable families (Section
3.2), and more general methods for constructing variational bounds (Section
3.3). Some of these methods differ with regard to which aspects of the
underlying Bayesian model are assumed to be computationally accessible—i.e.,
is the problem in the ‘standard’ setting where the likelihood function can be
evaluated, or is it in the ‘implicit model’ setting where likelihood
evaluations and possibly prior density evaluations are unavailable? Throughout
this section, we will comment on the applicability of the methods being
discussed to either setting.

###  3.1 Nested Monte Carlo estimators

For a generic expected utility design criterion, as formulated in (13), a
standard Monte Carlo estimator of the expectation employs pairs of samples
$(y^{(i)},\theta^{(i)})$ drawn from the joint prior of parameters and
observations given the design $\xi$, $p(y,\theta|\xi)$:

$\displaystyle
U(\xi)=\mathbb{E}_{Y,\Theta|\xi}[u(\xi,Y,\Theta)]\approx\frac{1}{N}\sum_{i=1}^{N}u(\xi,y^{(i)},\theta^{(i)}).$
(52)

These samples are typically obtained by first drawing $\theta^{(i)}\sim
p(\theta)$ from the prior and then drawing $y^{(i)}\sim p(y|\theta^{(i)},\xi)$
from the conditional density of the observations. The utility function $u$
must then be evaluated at the samples generated:
$u(\xi,y^{(i)},\theta^{(i)})$. This may not be easy to do. For example,
evaluating either of the utility functions $u^{\text{score}}$ (19) or
$u^{\text{div}}$ (20), which render $U$ equal to the expected information gain
(EIG) in $\Theta$, requires evaluating a normalized posterior density, where
the normalizing constant may change for each realization of $y$ and each value
of $\xi$. One way forward is to estimate these normalizing constants by
another Monte Carlo simulation, which gives rise to nested Monte Carlo (NMC)
estimators.

As a canonical/representative case, we describe NMC approaches to estimating
the EIG in parameters $\Theta$, from prior to posterior—i.e., the
$U_{\text{KL}}$ defined in (14) and the expressions thereafter. A widely used
estimator, proposed in Ryan (2003), employs the form of $U_{\text{KL}}$ given
in (22): we simply replace the posterior normalizing constant $p(y|\xi)$ with
a Monte Carlo estimate, as follows,

$\displaystyle U_{\text{KL}}(\xi)$ $\displaystyle=\iint
p(y,\theta|\xi)\log\frac{p(y|\theta,\xi)}{p(y|\xi)}\,\mathrm{d}\theta\,\mathrm{d}y$
$\displaystyle=\iint p(y,\theta|\xi)\log\frac{p(y|\theta,\xi)}{\int
p(y|\tilde{\theta},\xi)p(\tilde{\theta})\,\text{d}\tilde{\theta}}\,\mathrm{d}\theta\,\mathrm{d}y$
$\displaystyle\approx\frac{1}{N}\sum_{i=1}^{N}\left(\log
p(y^{(i)}|\theta^{(i)},\xi)-\log\left[\frac{1}{M}\sum_{j=1}^{M}p(y^{(i)}|\tilde{\theta}^{(i,j)},\xi)\right]\right)\eqqcolon\widehat{U}_{\text{KL}}^{N,M}(\xi).$
(53)

Here the ‘outer loop’ sample pairs $\\{(y^{(i)},\theta^{(i)})\\}_{i=1}^{N}$
are drawn from $p(y,\theta|\xi)$ as before, but we also require an
_independent_ collection of samples from the prior,

$\\{\tilde{\theta}^{(i,j)}\\}_{i=1,j=1}^{i=N,j=M}\sim p(\theta),$

amounting to $M$ independent samples for each outer-loop iteration. Overall,
evaluating $\widehat{U}_{\text{KL}}^{N,M}$ requires (i) the ability to sample
from the prior, (ii) the ability to sample from the statistical model for
$Y|\xi$, and (iii) the ability to evaluate the likelihood function.

Properties of the estimator $\widehat{U}_{\text{KL}}^{N,M}$ have been analyzed
in Ryan (2003), Beck et al. (2018), and Rainforth et al. (2018). It is biased
at finite $M$, but asymptotically unbiased and consistent. More precisely, the
leading order terms of its bias and variance are given by

$\displaystyle\mathbb{E}\left[\widehat{U}_{\text{KL}}^{N,M}(\xi)\right]-U_{\text{KL}}(\xi)$
$\displaystyle=\frac{C_{1}(\xi)}{M}+\mathcal{O}\left(\frac{1}{M^{2}}\right),$
(54)
$\displaystyle\mathbb{V}\text{ar}\left[\widehat{U}_{\text{KL}}^{N,M}(\xi)\right]$
$\displaystyle=\frac{C_{2}(\xi)}{N}+\frac{C_{3}(\xi)}{NM}+\mathcal{O}\left(\frac{1}{NM^{2}}\right),$
(55)

where $C_{1},C_{2},C_{3}$ are design-dependent constants. The constant $C_{1}$
is always positive (see Beck et al. 2018, Proposition 1), and thus the NMC
estimator $\widehat{U}_{\text{KL}}^{N,M}$ is, to leading order, positively
biased. Intuitively, bias arises because the Monte Carlo estimator of
$p(y|\xi)$ (which is unbiased) is transformed by a nonlinear function.

It is useful then to ask how to optimally allocate the sample sizes $M$ and
$N$ to minimize the mean-square error

$\text{MSE}=\mathbb{E}\big{[}(\widehat{U}_{\text{KL}}^{N,M}-U_{\text{KL}})^{2}\big{]}$

for any given budget. The total number of samples drawn is $W=(M+1)N$;
similarly, if cost lies in evaluating the statistical model $p(y|\theta,\xi)$
for any new value of $\theta$, then computing $\widehat{U}_{\text{KL}}^{N,M}$
incurs $W=(M+1)N$ evaluations. Letting $\alpha^{2}=M/N$ denote the ratio of
inner-to-outer loop sample sizes, one can show that the optimal value of this
ratio for any given $W$ scales as $\alpha_{\ast}^{2}=\mathcal{O}(W^{-1/3})$
(Beck et al. 2018, Feng and Marzouk 2019, Rainforth et al. 2018). In other
words, the ratio should decrease slowly as the computational budget increases.
This scaling translates to setting $M=\mathcal{O}(\sqrt{N})$, and an optimal
convergence rate of $\text{MSE}=\mathcal{O}(W^{-2/3})$. Significantly, this is
_slower_ than the standard Monte Carlo rate! Put another way, the
computational effort required to achieve an MSE of $\epsilon^{2}$ is
$\mathcal{O}(\epsilon^{-3})$, rather than $\mathcal{O}(\epsilon^{-2})$.

Many improvements upon the ‘vanilla’ NMC estimator (53) have been proposed. A
straightforward idea is to use importance sampling to estimate the evidence
$p(y|\xi)$ in the inner loop, i.e., to sample from some possibly $y$-dependent
biasing distribution rather than the prior $p_{\Theta}$. To make this
explicit, we write only the outer loop of (53) and replace the log-evidence
term with a plugin estimate,

$\displaystyle U_{\text{KL}}(\xi)$
$\displaystyle\approx\widehat{U}_{\text{KL}}^{N,\text{is}}(\xi)\coloneqq\frac{1}{N}\sum_{i=1}^{N}\left(\log
p(y^{(i)}|\theta^{(i)},\xi)-\log\widehat{p}(y^{(i)}|\xi)\right),$ (56)

where

$\displaystyle\widehat{p}(y^{(i)}|\xi)=\frac{1}{M}\sum_{j=1}^{M}p(y^{(i)}|\tilde{\theta}^{(i,j)},\xi)w^{(i,j)},\quad\tilde{\theta}^{(i,j)}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}q^{i,\xi},\quad
w^{(i,j)}=\frac{p(\tilde{\theta}^{(i,j)})}{q^{i,\xi}(\tilde{\theta}^{(i,j)})}.$
(57)

Here the superscripts on the density $q^{i,\xi}$ of the biasing distribution
emphasize that it can depend on the outer-loop index $i$, and thus on the
point $y^{(i)}$ where the evidence is being evaluated. The biasing
distribution generally depends on the design $\xi$ as well. The question then
becomes to how to choose this biasing distribution, for each summand of the
outer loop.

One natural choice, proposed in Beck et al. (2018), is to use a Laplace
approximation of the posterior distribution $p(\theta|y^{(i)},\xi)$ associated
with each outer-loop sample. The Laplace approximation seeks the point of
highest posterior density,
$\theta_{\text{map}}^{i}\equiv\theta_{\text{map}}(y^{(i)},\xi)\in\operatornamewithlimits{argmax}_{\theta}\log
p(\theta|y^{(i)},\xi)$, and builds a Gaussian approximation centered at this
point, with covariance equal to the negative Hessian of $\log
p(\theta|y^{(i)},\xi)$ evaluated at $\theta_{\text{map}}^{i}$. In general, we
define

$\displaystyle\Sigma_{\text{map}}(y,\xi)=\bigl{(}\left.\nabla_{\theta}^{2}\log
p\left(\theta|y,\xi\right)\right|_{\theta=\theta_{\text{map}}(y,\xi)}\bigr{)}^{-1}.$

The biasing distribution suggested by Beck et al. (2018) is then

$\displaystyle
q^{i,\xi}=\mathcal{N}\bigl{(}\theta_{\text{map}}^{i},\Sigma_{\text{map}}(y^{(i)},\xi)\bigr{)}.$
(58)

This choice aims to focus samples of the biasing distribution onto regions of
higher posterior density, thus reducing the variance of the evidence estimate
$\widehat{p}(y^{(i)}|\xi)$ and hence both the bias and the variance of
estimates of $U_{\text{KL}}$, as demonstrated in Beck et al. (2018). As is
typically the case with importance sampling, choosing a better biasing
distribution reduces the magnitudes of the constants $C_{1},C_{2},C_{3}$ in
(54)–(55), but does not change the rates of convergence of the estimator with
$M,N$.

Computing a Laplace approximation efficiently, which must be done here $N$
times (once for each outer-loop sample), generally requires gradients of the
log-posterior density, $\nabla_{\theta}\log p(\theta|y^{(i)},\xi)$, and a good
approximation of the Hessian of this log-density Bui-Thanh et al. (2013),
Schillings and Schwab (2016). We also note that, in general,
$\theta_{\text{map}}^{i}\neq\theta^{(i)}$; that is, the posterior mode does
not coincide with the data-generating value of the parameter. Englezou et al.
(2022) propose a simplification of Laplace-based importance sampling that
instead centers the biasing distribution at the data-generating value of the
parameter, drawn from the outer loop. Specifically, both the mean of the
Gaussian and the position of the Hessian evaluation in (58) are set to
$\theta^{(i)}$, rather than $\theta_{\text{map}}^{i}$. Doing so avoids the
cost of numerical optimization to find $\theta_{\text{map}}^{i}$, and may
yield only a modest increase in the MSE of the estimator at finite sample
sizes (Englezou et al. 2022). An (approximate) Hessian of the log-posterior,
however, is still required. We should also emphasize that choosing a Gaussian
approximation as a biasing distribution can become unstable when the posterior
is strongly non-Gaussian, and particularly if the posterior has heavy tails.
Here other choices of biasing—for instance, using the mean and covariance
matrix of the Laplace approximation to parameterize a heavier-tailed
multivariate-$t$ distribution—could be more robust Owen (2013).

Another approach to importance sampling, developed in Feng and Marzouk (2019),
is derivative-free. It is a multiple importance sampling scheme that proceeds
iteratively over the outer-loop index $i$: at any given $i$, past inner-loop
samples $\\{(\tilde{\theta}^{(k,j)})_{j=1}^{M}\\}_{k<i}$ and their associated
likelihood evaluations are used to create a _mixture_ biasing distribution
$q^{i,\xi}$ tailored to estimating the current evidence $p(y^{(i)}|\xi)$. To
make this process efficient and maximize reuse of information, the outer-loop
iterations are ordered from largest to smallest prior density of
$\theta^{(i)}$, $p(\theta^{(i)})$. Again, this approach can substantially
decrease the bias and variance of EIG estimates relative to a vanilla NMC
scheme.

A different way of accelerating nested Monte Carlo involves multilevel
formulations (Giles 2015). Goda et al. (2020) introduce a multilevel Monte
Carlo estimator of $U_{\text{KL}}$, where the level controls the number of
inner-loop samples. They also develop an antithetic coupling for the inner-
loop estimates that reduces variance, and show that the overall construction
reduces the computational complexity required to achieve an MSE of
$\epsilon^{2}$ to $\mathcal{O}(\epsilon^{-2})$, improving on the
$\mathcal{O}(\epsilon^{-3})$ complexity of the standard NMC estimator. In
other words, we recover the optimal Monte Carlo rate. Goda et al. (2020) also
show how importance sampling—for instance, the Laplace approximation-based
importance scheme discussed above—can be incorporated within their multilevel
Monte Carlo estimator to reduce constants of the error terms. Beck et al.
(2020) also introduce a multilevel scheme for EIG, varying not only the number
of inner-loop samples but also the discretization/approximation of some
underlying partial differential equation (PDE) model used to define the
likelihood.

All of the estimators discussed so far are consistent, that is, they converge
in probability to the true EIG as the relevant sample sizes (e.g., $M$ and
$N$, or their analogues in a multilevel scheme) are sent to infinity. Other
schemes proposed in the literature, with the goal of improving computational
efficiency, are not consistent. For instance, Long et al. (2013) propose
replacing the utility $u^{\text{div}}(\xi,y)$ (20), whose expectation yields
the EIG, with an approximation of the Kullback–Leibler (KL) divergence from
the prior to a Laplace approximation
$\mathcal{N}\bigl{(}\theta_{\text{map}}(y,\xi),\Sigma_{\text{map}}(y,\xi)\bigr{)}$
of the posterior $p(\theta|y,\xi)$. Their construction introduces a series of
additional approximations—replacing the true MAP estimate
$\theta_{\text{map}}(y,\xi)$ with the data-generating value of the parameter
$\theta$, using the Gauss–Newton approximation of the Hessian at this point
rather than at the MAP (the inverse of which yields a ‘covariance’ matrix
$\widetilde{\Sigma}(\theta,\xi)$ that does not depend on $y$), and taking
further asymptotic approximations of the relevant integrals—to obtain an
approximation of EIG that involves only integration over the prior,

$\displaystyle\text{EIG}(\xi)\approx\int\left(-\frac{1}{2}\log\det\widetilde{\Sigma}(\theta,\xi)-\frac{p}{2}\left(\log
2\pi+1\right)-\log p(\theta)\right)p(\theta)\,\text{d}\theta.$ (59)

As a design criterion, we note that this EIG approximation is similar in
structure to (12). The error of this approximation can be related to the
Gaussianity of the posterior _and_ the number of independent repeated trials
$N_{\text{tr}}$ of the experiments specified by $\xi$, and is bounded by
$\mathcal{O}(1/N_{\text{tr}})$ in probability (Long et al. 2013). Intuitively,
repeated trials cause the posterior to concentrate around its mode, and this
concentration controls the error of integral approximations leading to (59).
In practice, for small $N_{\text{tr}}$ and a non-Gaussian problem, this error
can be large. An alternative and perhaps more straightforward way of using the
Laplace approximation is proposed by Overstall et al. (2018); here the
approach is simply to replace the estimate of the log-evidence term $\log
p(y^{(i)}|\xi)$ in (56) with the log-evidence of the standard Laplace
approximation of the posterior $p(\theta|y^{(i)})$. The explicit likelihood
term and the outer Monte Carlo sum in (56) are unchanged. The bias of this EIG
approximation increases from zero as the posterior departs from Gaussianity,
but compared to the scheme in Long et al. (2013), it does not rely directly on
posterior concentration.

A different family of approximations follows by replacing computationally
expensive aspects of statistical model $p(y|\theta,\xi)$ with a
computationally cheaper ‘surrogate.’ For instance, suppose that the data arise
from a nonlinear forward operator $G_{\xi}:\mathbb{R}^{p}\to\mathbb{R}^{n}$
perturbed with additive noise $\mathcal{E}$, i.e.,
$Y=G_{\xi}(\theta)+\mathcal{E}$. This setup corresponds to a discretization of
the Bayesian inverse problems discussed in Section 2.3.1. Here, evaluating the
function $G_{\xi}$ often involves solving a set of PDEs or integral equations,
and it is natural to replace this solution with a cheaper approximation
$\widetilde{G}_{\xi}\approx G_{\xi}$, which in turn induces an approximation
$\widetilde{p}(y|\theta,\xi)$ of the likelihood and an approximation
$\widetilde{p}(\theta|y,\xi)$ of the posterior distribution via (3). Inserting
any of these approximations into (15) and equivalent expressions thus yields
an approximate EIG. By the same token, using evaluations of
$\widetilde{p}(y|\theta,\xi)$ and samples $\tilde{y}^{(i)}$ drawn from this
approximate model in any of the NMC estimators discussed above will yield
consistent estimates of this approximate EIG.

Huan and Marzouk (2013) proposed such a procedure, using polynomial
approximations $\widetilde{G}_{\xi}$ of $G_{\xi}$ built via sparse quadrature.
But a wide variety of other approximation schemes and formats are possible:
direct function approximation, whether via polynomials, Gaussian processes, or
neural networks (Herrmann et al. 2020); but also reduced-order models (Benner
et al. 2015) or even coarser numerical discretizations of the model giving
rise to the parameter-to-observable map $G_{\xi}$. The impact of such
approximations on the posterior distribution is by now reasonably well
understood, especially when one can build a family of approximations
$\widetilde{G}_{\xi}^{\ell}$, indexed by $\ell$, that converge (in some
appropriate sense) to $G_{\xi}$ as $\ell\to\infty$ (Marzouk and Xiu 2009,
Stuart 2010, Stuart and Teckentrup 2018, Sprungk 2020). The impact of such
approximations on the EIG has only recently been analyzed, however. Duong et
al. (2023) show that the difference between the EIG and its approximation is
controlled by the prior expectation of likelihood perturbations under the KL
divergence. As a consequence, in the setting of Gaussian likelihoods that we
sketched here, Duong et al. (2023, Theorem 4.4) show that closeness of
$\widetilde{G}_{\xi}^{\ell}(\theta)$ to $G_{\xi}(\theta)$ in a prior-weighted
$L^{2}$ sense, uniformly over designs $\xi\in\Xi$, guarantees uniform control
over the error in the approximate EIG $U^{\ell}_{\text{KL}}(\xi)$ and
convergence of the maximizers of $U^{\ell}_{\text{KL}}$ as $\ell\to\infty$.

We close this section with a caution. Even the consistent NMC estimators of
$U_{\text{KL}}(\xi)$ discussed so far are biased at finite inner loop sample
sizes. This bias can be significant but, more importantly, will vary with
$\xi$ in general. Figure 3.1, adapted from Feng and Marzouk (2019),
illustrates this phenomenon for a four-dimensional linear-Gaussian problem
(thus allowing comparison with exact solutions). EIG in a subset of the model
parameters (28)–(29) is estimated using an adaptation of the vanilla NMC
estimator $\widehat{U}_{\text{KL}}^{N,M}$ (53) (in red), and with the multiple
importance sampling scheme developed in Feng and Marzouk (2019) (in green).
Sample sizes are _fixed_ for all $\xi$, and the true EIG is shown in black. We
see that the bias of the vanilla estimator actually obscures the location of
the true maximum. The bias is largest where posterior concentration is
maximized, as this is where the prior-weighted estimates of the evidence have
greatest variance. The multiple importance scheme fares better, but there is
generally no guarantee regarding stability of the maxima for any finite sample
size. Care is thus needed to adjust the approximation of EIG in conjunction
with the optimization procedure. We will revisit this issue in Section 4.

Figure 3.1: Estimated EIG as a function of a scalar design parameter $\xi$ for
a linear-Gaussian model, using vanilla NMC (red) or an improved multiple
importance sampling scheme (green), compared to the true EIG (black). Shaded
areas represent the interval containing 95% of 2000 independent estimates of
EIG at each $\xi$; red dashed and solid green lines are the means of these
estimates. Figure adapted from Feng and Marzouk (2019).

###  3.2 Mutual information bounds via density approximations

As anticipated in the discussion following (21)–(22), one of the core
computational tasks of the NMC estimators discussed in Section 3.1 is to
estimate the posterior normalizing constant $p(y|\xi)$ across many different
values of $y$. The vanilla NMC approach does this entirely independently for
each value of $y$, as does the Laplace-based importance sampling method. (The
adaptive importance sampling scheme of Feng and Marzouk (2019), on the other
hand, could be said to ‘borrow’ information from other values of $y$ to create
a local biasing distribution.) A rather different way of approaching this
problem is to approximate the marginal density $p(y|\xi)$, or similarly, the
normalized posterior density $p(\theta|y,\xi)$, directly, e.g., in some
parametric family of densities.

Let us recall (21)–(23) in a more concise form:

$\displaystyle\mathcal{I}(Y;\Theta|\xi)$
$\displaystyle=\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{p(Y|\Theta,\xi)}{p(Y|\xi)}\right]$
(60)
$\displaystyle=\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{p(\Theta|Y,\xi)}{p(\Theta)}\right].$
(61)

In either (60) or (61), one could in principle seek to approximate the density
in the numerator, the density in the denominator, or both. Suppose that we
replace $p(y|\xi)$ in (60) with some approximating probability density
function $q^{\text{mar}}(y|\xi)$ (where the superscript stands for
‘marginal’). Then, as noted in Barber and Agakov (2003), Poole et al. (2019),
and Foster et al. (2019),

$\displaystyle\mathcal{I}(Y;\Theta|\xi)$
$\displaystyle=\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{p(Y|\Theta,\xi)\,q^{\text{mar}}(Y|\xi)}{q^{\text{mar}}(Y|\xi)\,p(Y|\xi)}\right]$
$\displaystyle=\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{p(Y|\Theta,\xi)}{q^{\text{mar}}(Y|\xi)}\right]-D_{\text{KL}}(p_{Y|\xi}||q^{\text{mar}}_{Y|\xi})$
$\displaystyle\leq\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{p(Y|\Theta,\xi)}{q^{\text{mar}}(Y|\xi)}\right]=\mathbb{E}_{\Theta}\left[D_{\text{KL}}(p_{Y|\Theta,\xi}||q^{\text{mar}}_{Y|\xi})\right],$
(62)

where the inequality follows from the non-negativity of the KL divergence.
Hence, for any approximation $q^{\text{mar}}(y|\xi)$ of the marginal density
of $Y|\xi$, (62) is an upper bound on the mutual information (EIG).

Similarly, if we replace $p(\theta|y,\xi)$ in (61) with some approximating
probability density function $q^{\text{pos}}(\theta|y,\xi)$ (where the
superscript denotes ‘posterior’), we obtain

$\displaystyle\mathcal{I}(Y;\Theta|\xi)$
$\displaystyle=\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{q^{\text{pos}}(\Theta|Y,\xi)\,p(\Theta|Y,\xi)}{p(\Theta)\,q^{\text{pos}}(\Theta|Y,\xi)}\right]$
$\displaystyle=\mathbb{E}_{Y,\Theta|\xi}\left[\log\frac{q^{\text{pos}}(\Theta|Y,\xi)}{p(\Theta)}\right]+\mathbb{E}_{Y|\xi}\left[D_{\text{KL}}(p_{\Theta|Y,\xi}||q^{\text{pos}}_{\Theta|Y,\xi})\right]$