∎

11institutetext: 30 Blue Ridge Lane, West Hartford, CT, USA
Tel.: +18605018265
11email<EMAIL_ADDRESS>
# Multimodal Information Gain in Bayesian Design of Experiments

Quan Long

(Received: date / Accepted: date)

###### Abstract

One of the well-known challenges in optimal experimental design is how to
efficiently estimate the nested integrations of the expected information gain.
The Gaussian approximation and associated importance sampling have been shown
to be effective at reducing the numerical costs. However, they may fail due to
the non-negligible biases and the numerical instabilities. A new approach is
developed to compute the expected information gain, when the posterior
distribution is multimodal - a situation previously ignored by the methods
aiming at accelerating the nested numerical integrations. Specifically, the
posterior distribution is approximated using a mixture distribution
constructed by multiple runs of global search for the modes and weighted local
Laplace approximations. Under any given probability of capturing all the
modes, we provide an estimation of the number of runs of searches, which is
dimension independent. It is shown that the novel global-local multimodal
approach can be significantly more accurate and more efficient than the other
existing approaches, especially when the number of modes is large. The methods
can be applied to the designs of experiments with both calibrated and
uncalibrated observation noises.

###### Keywords:

Weighted Laplace approximation Gaussian mixture Expected information gain
Optimal design of experimentsMachine learning

##  1 Introduction

Statistical optimal experimental design is a systematic approach to increase
data efficiency (i.e., actively acquiring as less data as possible to train a
model for given tolerance). In linear optimal experimental design, the
optimalities are the norms of the predictive covariance matrix, for example,
A-optimality, D-optimality and E-optimality atkinson2007optimum
atkinson1992optimum . A-optimality denotes the trace of the predictive
covariance matrix, D-optimality denotes the determinant of the predictive
covariance matrix and E-optimality refers to the matrix norm of the predictive
covariance matrix. The optimalities can also be correspondingly interpreted
from geometrical perspectives. For example, Titterington  studied the
D-optimality from such a point of view.

Nonlinear optimal experimental design is more complicated in that the
corresponding optimalities rely on the values of the unknown parameters. The
Bayesian approach provides a powerful framework for nonlinear designs in that
the prior knowledge of the unknown parameters is naturally incorporated in
terms of the relationship between the posterior distribution, the likelihood
function and the prior distribution. In a Bayesian framework, the optimal
design of an experiment can be obtained by maximizing the expected information
gain Chaloner1995 Verdinelli2000 . Specifically, the expected information gain
is the expected logarithmic ratio of the posterior probability density
function (pdf) and the prior pdf of the unknown random parameters kullback1951
. The expected information gain has been widely used as a measurement of
information in engineering and sciences, for example, seismology LONG2015123 ,
combustionBISETTI20160112 , tomographyLONG201324 ,
biologydehideniya2019synthetic , nonlinear dynamical systems Giovanni
HUAN2013288  and active learning Visser2008 .

Obtaining an accurate estimation of the expected information gain is
challenging. Using a direct nested sample average HUAN2013288 Ryan2003  leads
to a computational complexity of ${\cal{O}}(TOL^{-3})$, where $TOL$ is the
tolerance on the estimation error BECK2018523 . Many efforts have been made to
accelerate the computations involved in estimating the expected information
gain BECK2018523 BISETTI20160112 LONG2015123 LONG201324 LONG2015849 , where
the Laplace approximation Tierney1986 Tierney1989  has been used to
analytically complete the inner integration and to provide a proposal
distribution in an importance sampling strategy. The Laplace approximation
Tierney1986 Tierney1989  expands a posterior pdf at the mode and results in a
Gaussian integral which approximates the original integration. The discrepancy
between the Laplace approximation and the true integral diminishes as the
number of data increases LONG201324 Tierney1989 . It has also been extended to
the cases, where the posterior pdf concentrates on a submanifold in
LONG2015849 . The Laplace-based importance sampling has been used in the
sequential design of experiments in senarathne2019laplace . It is also note-
worthy that the multilevel Monte Carlo method and the Laplace-based importance
sampling can be combined to balance the number of samples in the inner and
outer integrals, leading to efficient computations of the expected information
gain Beck2019 Takashi2020 . The Laplace approximation has been used together
with the gradient descent method to tackle the continuous optimization
problems in optimal Bayesian experimental design in carlon2020nesterov
CHEN2019163 . Some authors used surrogate to approximate the expected
information gain against the design space Overstall2017 TARAKANOV2020 .
Additionally, it is note-worthy that consistent formulas have been derived for
Bayesian optimal experimental design based on infinite dimensional models, for
example, the partial differential equations (pdes) Alexanderian2014
Alexanderian2018 .

In the previous literature, the applicability of the Laplace approximation
relies on the assumption that a single mode dominates the overall shape of the
posterior pdf, with the exception of LONG2015849  which developed the Laplace
approximations on the data-informed submanifold. However, this assumption
clearly fails to be true, when the posterior pdf is multimodal (see Ihler2005
Lan2014 Shaw2007  for examples of multimodal posterior pdfs). Hence, the
imperceptive applications of the Laplace approximation and the importance
sampling scheme based on the Laplace approximation would induce irreducible
bias and numerical instability.

In this paper, a multimodal Laplace approximation is proposed to remove the
forementioned restrictions of the Laplace approximation. Under some mild
regularity constraints on the posterior pdf, the approach offers the following
advantages: First, it removes the bias of the conventional Laplace
approximation; second, it is numerically more stable, when used as a proposal
distribution in an importance sampling scheme; third, it avoids numerical
under-flow BECK2018523 , when the direct nested sampling is used in a big data
scenario.

The outline of the rest of the contents is as follows. Section 2 briefly
reviews Bayesian optimal experimental design, the expected information gain
and its numerical approximations. Section 3 derives and analyzes the
multimodal Laplace approximation. Given the probability of capturing all the
modes, a lower bound is provided for the number of independent mode searches.
Section 4 introduces an importance sampling scheme based on the multimodal
Laplace approximations presented in Section 3. The accuracy and efficiency of
the proposed methods are demonstrated using two numerical experiments in
Section 5.

##  2 Bayesian optimal experimental design and expected information gain

We firstly provide the context of our experimentation. In the current
scenario, we can take $m$ observations of the responses of an experiment under
the fixed experimental setup. Specifically, we assume the experimental
observations can be decomposed into a deterministic model and a measurement
noise:

$\bm{y}_{i}=\bm{g}(\bm{\theta}_{t},\bm{\xi})+\bm{\epsilon}_{i}\quad\text{with}\quad\bm{\epsilon}_{i}\sim{\cal{N}}(\bm{0},\bm{\Sigma}_{e})\quad\text{and}\quad
i=1,...,m\,,$

where $\bm{y}_{i}\in\mathbb{R}^{d_{m}}$ is the $i^{th}$ observation of the
measurement data, $\bm{\theta}_{t}\in{\mathbb{R}}^{d}$ is the vector of the
“true” parameters, $\bm{g}(\bm{\theta}_{t},\bm{\xi})$ is a deterministic
function of the system parameters and the experimental setup,
$\bm{\epsilon}_{i}$ is the observational noise, $\bm{\Sigma}_{e}$ is the
covariance matrix of the random measurement noise, $\bm{\xi}\in\mathbb{R}^{s}$
is the vector of the design parameters, for example, the locations of the
electrodes in a tomography experiment LONG201324  and the initial temperature
in a combustion experiment BISETTI20160112 . For the sake of conciseness, we
denote $\bar{\bm{y}}=[\bm{y}_{1}\,,\bm{y}_{2}\,,...,\bm{y}_{m}]^{\top}$ in the
rest of the paper. Note that we restrict our methodology to batch experimental
design and multiple different experimental setup can be considered by the
multi-dimensional $\bm{\xi}$ vector.

Consequently, the Bayes Theorem can be expressed as follows:

$\displaystyle
p(\bm{\theta}|\bar{\bm{y}},\bm{\xi})=\frac{p(\bar{\bm{y}}|\bm{\theta},\bm{\xi})p(\bm{\theta})}{p(\bar{\bm{y}}|\bm{\xi})}\,.$
(1)

where $p(\cdot)$ denotes a pdf. $p(\bar{\bm{y}}|\bm{\theta},\bm{\xi})$,
$p(\bm{\theta})$, $p(\bar{\bm{y}}|\bm{\xi})$ and
$p(\bm{\theta}|\bar{\bm{y}},\bm{\xi})$ denote the likelihood function, the
prior pdf, the marginal likelihood function (also called the evidence), and
the posterior pdf, respectively.

The current derivation of Laplace approximation is based on the additive model
of data and the noise is restricted to be Gaussian. That said, the validity of
Laplace approximation is not compromised, while other scenarios are
considered, for example, data $\bm{y}$ are discrete random variables, noises
in the responses are stemmed from intrinsic randomness Calvetti .

Note that we can similarly derive a Bayesian framework for the parameter
identification of time dependent problems, which has been demonstrated in
LONG2015123 . Additionally, the assumption of additive noise is a flexible
one, for example, a multiplicative model (i.e.,
$\bm{y}=\bm{g}\cdot\bm{\epsilon}$) can be converted to its additive form via a
logarithmic transformation of the involved variables.

We focus on the information efficiency of the experiments denoted by the setup
$\bm{\xi}$, which controls the generation of data $\bar{\bm{y}}$ in terms of
model $\bm{g}$. A widely used metric is the so-called expected information
gain kullback1951 Chaloner1995 . It measures the amount of new “information”
encoded in $p(\bm{\theta}|\bar{\bm{y}},\bm{\xi})$ on top of the legacy
“information” represented by $p(\bm{\theta})$. It is essentially an averaged
Kullback-Leibler (K-L) divergence kullback1951  (also called discrimination
information Ghosh1987 ) against all the possible values of the data. It can be
used as a utility function to quantify the information contents in the data
with respect to (w.r.t.) the unknown parameters and other focused unknown
quantities of interest (QoI) LONG201324 . The expected information gain
associated to (1) can be written as the following:

$\displaystyle
I=\int_{\cal{\bm{Y}}}\int_{{\bm{\Theta}}}log\left[\frac{p(\bm{\theta}|\bar{\bm{y}})}{p(\bm{\theta})}\right]p(\bm{\theta}|\bar{\bm{y}})p(\bar{\bm{y}})d\bm{\theta}d\bar{\bm{y}}\,,$
(2)

where we ignore the notation of $\bm{\xi}$ in the conditional distributions
for the sake of conciseness. To the best of our knowledge, Ryan2003  is the
first paper which noted the computational challenge of estimating the utility
function using the so called double loop Monte Carlo method (DLMC). The DLMC
sampler, which has been benchmarked against by many more efficient methods in
recent years, is a direct nested random sampling method:

$\displaystyle
I_{DLMC}=\frac{1}{M_{1}}\sum^{M_{1}}_{i=1}log\left[\frac{p(\bar{\bm{y}}_{i}|\bm{\theta}_{i})}{\frac{1}{N_{1}}\sum^{N_{1}}_{j=1}p(\bar{\bm{y}}_{i}|\bm{\theta}_{j})}\right]\,,$
(3)

where $M_{1}$ is the number of samples in the outer loop, and $N_{1}$ is the
number of samples in the inner loop. $\bm{\theta}_{i}\sim p(\bm{\theta})$,
$\bar{\bm{y}}_{i}\sim p(\bar{\bm{y}}|\bm{\theta}_{i})$, and
$\bm{\theta}_{j}\sim p(\bm{\theta})$. The variance and bias of this estimator
can be bounded by
$\mathop{\mathbb{V}}(I_{DLMC})={\cal{O}}\left(\frac{1}{M_{1}}\right)$ and
$\mathop{\mathbb{B}}(I_{DLMC})={\cal{O}}\left(\frac{1}{N_{1}}\right)$,
respectively (see Ryan2015 LONG201324 LONG2015849  for details of how the
errors are estimated ). The computational cost of this estimator can be
characterized by the number of likelihood computations (i.e.,
$W_{1}=N_{1}\times M_{1}$).

If a numerical tolerance is imposed on the mean square error:

$\mathop{\mathbb{V}}(I_{DLMC})+\mathop{\mathbb{B}}(I_{DLMC})^{2}=TOL^{2}\,,$

we can obtain the magnitudes of $N_{1}$ and $M_{1}$ in terms of the tolerance:

$M_{1}={\cal{O}}\left(TOL^{-2}\right)\quad\text{and}\quad
N_{1}={\cal{O}}\left(TOL^{-1}\right)\,.$

Consequently, the cost of the total computational cost is proportional to the
cubic of the reciprocal of the tolerance:

$W_{1}={\cal{O}}\left(TOL^{-3}\right)\,.$

Note that the bias caused by the inexact numerical approximations of the
forward model $\bm{g}(\bm{\theta})$ is not considered in this study.
Interested readers can refer to BECK2018523  from this perspective.

One of the more efficient approaches than the DLMC method is to invoke the
Laplace approximation of the inner integral. The Laplace approximation laplace
(also referred to as Laplace method) is conventionally used to convert the
integral of an exponential function to a Gaussian integral, which can be
written as follows:

$\displaystyle\int
e^{-mf(x)}dx=\sqrt{\frac{2\pi}{m|f^{\prime\prime}(\hat{x})|}}e^{-mf(\hat{x})}+{\cal{O}}\left(\frac{1}{m}\right),\,$
(4)

where $\hat{x}$ is the global optimum of $f(x)$. The asymptotic error,
${\cal{O}}\left(\frac{1}{m}\right)$, may have different rate w.r.t. $m$,
dependent on how the expansion of $f(x)$ is truncated at $\hat{x}$ Tierney1986
Tierney1989 schillings2020 . Expanding the log-posterior function at the
maximum a posterior (MAP) estimate and following a multivariate version of (4)
(see LONG201324  for details of the derivation procedure and error analysis),
the K-L divergence (the inner integral of (2) ) can be analytically obtained
using the Laplace approximation. Consequently, the expected information gain
can be approximated as follows:

$\displaystyle I=$
$\displaystyle\int_{\Theta}\int_{\cal{\bm{Y}}}\left[-\frac{1}{2}log\left(|\bm{\Sigma}|\right)-\frac{\bm{tr}(\bm{\Sigma}\bm{H}_{h}(\hat{\bm{\theta}}))}{2}\right]p(\bar{\bm{y}}|\bm{\theta}_{t})p(\bm{\theta}_{t})d\bar{\bm{y}}d\bm{\theta}_{t}$
$\displaystyle-\frac{d}{2}-\frac{d}{2}log(2\pi)+{\cal{O}}\left(\frac{1}{m}\right)\,,$
(5)

where the posterior covariance matrix can be approximated by the inverse of
the Hessian matrix of the log-posterior function at the MAP estimate (i.e.,
$\bm{\Sigma}\approx\bm{H}^{-1}(\hat{\bm{\theta}})\approx\left[\nabla_{\bm{\theta}}\bm{g}^{\top}(\hat{\bm{\theta}})\bm{\Sigma}^{-1}_{\epsilon}\nabla_{\bm{\theta}}\bm{g}(\hat{\bm{\theta}})\right]^{-1}$),
$\hat{\bm{\theta}}$ is the MAP estimate of $\bm{\theta}$ conditioned on the
synthetic data $\bar{\bm{y}}$, $\bm{tr}(\cdot)$ denotes matrix trace and
$\bm{H}_{h}(\hat{\bm{\theta}})=\nabla_{\bm{\theta}}\nabla_{\bm{\theta}}logp(\hat{\bm{\theta}})$.
The corresponding numerical discritization of (5) can be written as follows:

$\displaystyle I_{LA}=$
$\displaystyle\frac{1}{M_{2}}\sum^{M_{2}}_{i=1}\left[-\frac{1}{2}log\left(|\bm{\Sigma}_{i}|\right)-\frac{\bm{tr}(\bm{\Sigma}_{i}\bm{H}_{h}(\hat{\bm{\theta}}_{i}))}{2}\right]-\frac{d}{2}-\frac{d}{2}log(2\pi)\,,$
(6)

where the $\hat{\bm{\theta}}_{i}$ is the MAP estimate of $\bm{\theta}$
conditioned on the synthetic data $\bar{\bm{y}}_{i}\sim
p(\bar{\bm{y}}|{\bm{\theta}_{t}}^{i})$ and ${\bm{\theta}_{t}}^{i}\sim
p(\bm{\theta}_{t})$. The computational cost of estimating the inner integral
of (3) by running thousands times the likelihood functions can now be reduced
to one run of optimization to search for $\hat{\bm{\theta}}_{i}$.

In the scenario where the posterior Hessian matrix $\bm{H}(\hat{\bm{\theta}})$
is low-rank LONG2015849  (i.e., $r<d$ with $r$ be the rank of
$\bm{H}(\hat{\bm{\theta}})$), new information is only gained in the subspace
of the parameters, which are sensitive to the perturbations of the observed
data. The non-intuitive relationship between the expected information gain and
the nonzero eigenvalues of the posterior Hessian matrix has been revealed in
LONG2015849 . The Laplace approximation of the expected information gain can
be extended to the following form:

$\displaystyle
I=\int_{\Theta}\int_{{\cal{Y}}}\left[-\frac{1}{2}log|\bm{\Sigma}_{p}|-\frac{r}{2}-\frac{r}{2}log(2\pi)-log\int_{\bm{T}}p_{\bm{s},\bm{t}}(\bm{0},\bm{t})d\bm{t}\right]p(\bar{\bm{y}}|\bm{\theta}_{t})p(\bm{\theta}_{t})d\bm{\theta}_{t}+{\cal{O}}\left(\frac{1}{m}\right)\,,$
(7)

where
$\bm{\Sigma}_{p}(\hat{\bm{\theta}})=\left[\bm{U}^{\top}(\hat{\bm{\theta}})\bm{H}(\hat{\bm{\theta}})\bm{U}(\hat{\bm{\theta}})\right]^{-1}$,
the columns of $\bm{U}$ are the basis spanning the orthogonal space of the
Jacobian kernal of $\bm{H}$,
$\int_{\bm{T}}p_{\bm{s},\bm{t}}(\bm{0},\bm{t})d\bm{t}$ is an integral on the
non-informative manifold $\bm{T}$, $\bm{t}$ is the new variable which
parameterizes $\bm{T}$ and $\bm{s}$ is the variable which parameterizes the
normal direction of the manifold, $p_{\bm{s},\bm{t}}(\bm{s},\bm{t})$ is the
prior pdf after a change of parameters (i.e.,
$p_{\bm{s},\bm{t}}(\bm{s},\bm{t})d\bm{s}d\bm{t}=p(\bm{\theta})d\bm{\theta}$).

The equation (7) can be discritized as

$\displaystyle
I_{LAS}=\frac{1}{M_{3}}\sum^{M_{3}}_{i=1}\left[-\frac{1}{2}log|\bm{\Sigma}_{p}^{i}|\right]-\frac{r}{2}-\frac{r}{2}log(2\pi)-log\int_{\bm{T}}p_{\bm{s},\bm{t}}(\bm{0},\bm{t})d\bm{t}\,,$
(8)

where
$\bm{\Sigma}_{p}^{i}(\hat{\bm{\theta}})=\left[\bm{U}^{\top}(\hat{\bm{\theta}}_{i})\bm{H}(\hat{\bm{\theta}}_{i})\bm{U}(\hat{\bm{\theta}}_{i})\right]^{-1}$.
Note that (7) and (8) become (5) and (6), respectively, when $r=d$.

Although Laplace approximation enables efficient estimation of the inner loop
of (3), it entails an asymptotic bias when the posterior distribution is non-
Gaussian. To achieve both acceleration and consistency, importance sampling
based on Laplace approximation has been developed in BECK2018523  and Ryan2015
for optimal experimental design. The specific scheme can be written as follows

$\displaystyle
I_{LAIS}=\frac{1}{M_{4}}\sum^{M_{4}}_{i=1}log\left[\frac{p(\bar{\bm{y}}_{i}|\bm{\theta}_{i})}{\frac{1}{N_{4}}\sum^{N_{4}}_{j=1}p(\bar{\bm{y}}_{i}|\bm{\theta}_{j})\alpha_{j}^{i}}\right]\,,$
(9)

where the likelihood ratio,
$\alpha_{j}^{i}=\frac{p(\bm{\theta}_{j})}{p_{g}^{i}(\bm{\theta}_{j})}$,
$p_{g}^{i}$ is a Gaussian pdf with mean $\hat{\bm{\theta}}_{i}$ and covariance
matrix $\bm{\Sigma}_{i}$,
$\bm{\theta}_{j}\sim{\cal{N}}(\hat{\bm{\theta}}_{i},\bm{\Sigma}_{i})$. On top
of being efficient, direct Laplace approximation and corresponding importance
sampling can both avoid the numerical underflow as demonstrated in BECK2018523
, which occurs with very high probability in DLMC, when a large amount of data
lead to a very concentrated posterior pdf. Note that a Student’s
t-distribution should be used as the proposal distribution in the importance
sampling scheme, if the posterior pdf has a heavy tail.

Method | Complexity
---|---
DLMC | $W_{1}=M_{1}\times(C_{l}+N_{1}\times C_{l})$
LA | $W_{2}=M_{2}\times C_{o}$
Manifold LA | $W_{3}=M_{3}\times(C_{o}+C_{e})$
LAIS | $W_{4}=M_{4}\times(C_{o}+C_{l}+N_{4}\times C_{l})$
Table 1: Computational costs of the numerical schemes in the literature.

The computational costs of the numerical schemes, which we have reviewed so
far, are listed in Table 1, where $C_{o}$ denotes the cost of a single run of
optimization algorithm (e.g. the gradient descent method) used to find the MAP
estimate and the associated Hessian matrix, $C_{l}$ is the cost of computing a
likelihood function including computing the forward model
$\bm{g}(\bm{\theta})$ and the associated pdf, $C_{e}$ is the cost of the
numerical eigenvalue decomposition of the Hessian matrix.

Number of samples | Estimation
---|---
$M_{1}$ | $\mathop{\mathbb{V}}\left(log\left[\frac{p(\bar{\bm{y}}|\bm{\theta})}{p(\bar{\bm{y}})}\right]\right)\times TOL^{-2}$
$M_{2}$ | $\mathop{\mathbb{V}}\left(\frac{1}{2}log|\bm{H}(\bm{\theta})|\right)\times TOL^{-2}$
$M_{3}$ | $\mathop{\mathbb{V}}\left(\frac{1}{2}log|\bm{H}_{p}(\bm{\theta})|\right)\times TOL^{-2}$
$M_{4}$ | $\mathop{\mathbb{V}}\left(log\left[\frac{p(\bar{\bm{y}}|\bm{\theta})}{p(\bar{\bm{y}})}\right]\right)\times TOL^{-2}$
$N_{1}$ | $\frac{1}{2}{\mathop{\mathbb{E}}}\left(\mathop{\mathbb{V}}\left[\frac{p(\bar{\bm{y}}|\bm{\theta})}{p(\bar{\bm{y}})}|\bar{\bm{y}}\right]\right)\times TOL^{-1}$
$N_{4}$ | $\frac{1}{2}{\mathop{\mathbb{E}}}\left(\mathop{\mathbb{V}}\left[\frac{p(\bar{\bm{y}}|\bm{\theta})\alpha(\bm{\theta})}{p(\bar{\bm{y}})}|\bar{\bm{y}}\right]\right)\times TOL^{-1}$
Table 2: Estimation of the numbers of samples.

We tabulate the estimates of the number of samples in all of the introduced
estimators in Table 2. Note that $N_{4}\ll N_{1}$ for a concentrated single
modal posterior pdf, because of the much reduced variance in the importance
sampling scheme of LAIS.

##  3 Multimodal Laplace approximation (MLA)

The numerical integrations based on the Laplace approximations are only valid
in the cases, where the posterior pdf has a single dominant mode. An example
of Gaussian mixture is used here to illustrate the potential erroneous
information gain caused by a single modal Laplace approximation. In Figure 1,
the solid curve denotes a Gaussian mixture pdf of three modes and the dashed
line denotes a Gaussian pdf matching the first mode of the mixture pdf. The
multimodal pdf involves three weighted modes:
$0.3\times{\cal{N}}(2,0.5)+0.3\times{\cal{N}}(5,0.2)+0.4\times{\cal{N}}(7,0.5)$.
Note that the error in considering a Gaussian pdf rather than a truncated
Gaussian pdf is negligible.

Figure 1: A Gaussian mixture pdf and its single modal approximation using the
Laplace approximation.

Its single modal approximation can be the Gaussian distribution:
${\cal{N}}(2,0.5)$. Using a uniform distribution, ${\cal{U}}(0,10)$, as the
prior pdf, and decomposing the K-L divergence into three integrals
corresponding to the disparate modes, the approximated information gain can be
analytically estimated by the following expression

$\displaystyle
D_{KL}\approx\sum^{3}_{k=1}D_{KL}(w_{k}{\cal{N}}(\mu_{k},\sigma^{2}_{k})\,||\,{\cal{U}}(0,10))=\sum^{3}_{k=1}\left[w_{k}log(w_{k})-0.5w_{k}log(\sigma^{2}_{k})-w_{k}h\right]\,,$
(10)

where $w_{k}$, $\mu_{k}$ and $\sigma^{2}_{k}$ are the weight, the mean and the
variance of the $k^{th}$ Gaussian mode, respectively, $h=-log10$ is the
logarithm of the prior pdf. Substituting the values of $w_{k}$, $\mu_{k}$ and
$\sigma^{2}_{k}$ into (10), we obtain ${D_{KL}}\approx 0.97$. To the contrary,
a conventional Laplace approximation taking the first mode as the single
dominant mode leads to $D_{KL}^{laplace}\approx
0.5log(\sigma_{1}^{2})+h(\mu_{1})\approx 1.23$.

In a special case, where $w_{k}=1/K$ and $\sigma^{2}_{k}=\sigma^{2}$,
$k=1,...,K$, the discrepancy between the conventional Laplace approximation
and the true value of the information gain can be found as the following:

$D_{KL}-D_{KL}^{laplace}=logK\,,$

which increases as the number of modes increases.

Next, we present the multidimensional multimodal approximation of the expected
information gain in Theorem 1.

###### Theorem 3.1

Assuming that the posterior distribution of $\bm{\theta}$ conditioning on the
synthetic data $\bar{\bm{y}}$ entails $K$ modes, data being modeled by
$\bm{y}_{i}=\bm{g}(\bm{\theta}_{t})+\bm{\epsilon}_{i}$, the third derivatives
of the model $\bm{g}(\bm{\theta}_{t})$ being bounded from above, we can
approximate the expected information gain as follows

$\displaystyle I=$
$\displaystyle\int_{\bm{\Theta}}\int_{{\cal{\bm{Y}}}}\sum^{K}_{k=1}[w_{k}log(w_{k})-0.5w_{k}log(|\bm{\Sigma}_{k}|)$
$\displaystyle-
w_{k}h(\hat{\bm{\theta}}_{k})]p(\bar{\bm{y}}|\bm{\theta}_{t})p(\bm{\theta}_{t})d\bar{\bm{y}}d\bm{\theta}_{t}-0.5log(2\pi)^{d}$
$\displaystyle-0.5d+{\cal{O}}\left(\frac{1}{m}\right)\,,$ (11)

where $\hat{\bm{\theta}}_{k}$ is the $k^{th}$ mode of the posterior pdf,
$w_{k}$ is the weight of the $k^{th}$ mode, $\bm{\Sigma}_{k}$ is the inverse
of the Hessian matrix of the negative logarithm of the posterior pdf at the
$k^{th}$ mode and $h(\hat{\bm{\theta}}_{k})=log(p(\hat{\bm{\theta}}_{k}))$.

###### Proof

We first express the integral of the K-L divergence as a summation of several
integrals on the sub-domains related to the multiple disparate modes:

$\displaystyle D(\bar{\bm{y}})=$
$\displaystyle\sum^{K}_{k=1}\int_{\Omega_{k}}log\left[\frac{p(\bm{\theta}|\bar{\bm{y}})}{p(\bm{\theta})}\right]p(\bm{\theta}|\bar{\bm{y}})d\bm{\theta}+{\cal{O}}\left(e^{-m}\right)\,,$
(12)

where $\Omega_{k}$ is a ball of radius $R_{k}$ covering a neighborhood of the
$k^{th}$ mode. According to the Laplace principle Dembo , most of the
probability mass is concentrated at the modes as $m$ increases, namely, a big-
data scenario schillings2020 . The integral outside
$\bm{\Omega}={\displaystyle\cup}_{k=1}^{K}\Omega_{k}$ adds up to $Ce^{-m}$
with the constant
$C=e^{ess\,inf_{\bm{\theta}\in\bm{\Theta}\setminus\Omega}\phi(\bm{\theta})}$,
$\phi(\bm{\theta})$ is the integrand function in (12) scaled by the number of
data:

$\phi(\bm{\theta})=\frac{1}{m}log\left[\frac{p(\bm{\theta}|\bar{\bm{y}})}{p(\bm{\theta})}\right]p(\bm{\theta}|\bar{\bm{y}})\,.$

Next, locally approximating $log[p(\bm{\theta}|\bar{\bm{y}})]$ in each
$\Omega_{k}$ by a truncated second order Taylor expansion leads to:

$\displaystyle
p(\bm{\theta}|\bar{\bm{y}})=\tilde{p}_{k}(\bm{\theta}|\bar{\bm{y}})+e^{\left[{\cal{O}}\left(|\bm{\theta}-\hat{\bm{\theta}}_{k}|^{3}\right)\right]}\,,$

where

$\displaystyle\tilde{p}_{k}(\bm{\theta}|\bar{\bm{y}})$
$\displaystyle=p(\hat{\bm{\theta}}_{k}|\bar{\bm{y}})exp\left[-\frac{1}{2}(\bm{\theta}-\hat{\bm{\theta}}_{k})^{\top}\bm{\Sigma}_{k}^{-1}(\bm{\theta}-\hat{\bm{\theta}}_{k})\right]$
$\displaystyle=w_{k}p^{g}_{k}(\bm{\theta})\,,$

with $p^{g}_{k}(\bm{\theta})$ be the pdf of the Gaussian distribution-
${\cal{N}}(\hat{\bm{\theta}}_{k},\bm{\Sigma}_{k})$, with
$\bm{\Sigma}_{k}^{-1}=-\nabla_{\bm{\theta}}\nabla_{\bm{\theta}}log\left[p(\hat{\bm{\theta}}_{k}|\bar{\bm{y}})\right]$
is the Hessian matrix of the negative logarithm of the posterior pdf at
$\hat{\bm{\theta}}_{k}$ and
$w_{k}=p(\hat{\bm{\theta}}_{k}|\bar{\bm{y}})\sqrt{2\pi}^{d}|\bm{\Sigma}_{k}|^{1/2}$.
Note that the first derivative of $log[p(\bm{\theta}|\bar{\bm{y}})]$ at
$\hat{\bm{\theta}}_{k}$ is zero (i.e.,
$\nabla_{\bm{\theta}}log\left[p(\hat{\bm{\theta}}_{k}|\bar{\bm{y}})\right]=0$).
Substituting $p(\bm{\theta}|\bar{\bm{y}})$ by
$\tilde{p}_{k}(\bm{\theta}|\bar{\bm{y}})$ in $\Omega_{k}$, we can approximate
the K-L divergence of (12) as follows:

$\displaystyle D(\bar{\bm{y}})=$
$\displaystyle\sum^{K}_{k=1}\int_{\Omega_{k}}log\left[\frac{\tilde{p}_{k}(\bm{\theta}|\bar{\bm{y}})}{p(\bm{\theta})}\right]\tilde{p}_{k}(\bm{\theta}|\bar{\bm{y}})d\bm{\theta}+{\cal{O}}\left(\frac{1}{m}\right)$
$\displaystyle=$
$\displaystyle\sum^{K}_{k=1}\int_{\Omega_{k}}log\left[\frac{w_{k}p^{g}_{k}(\bm{\theta})}{p(\bm{\theta})}\right]w_{k}p^{g}_{k}(\bm{\theta})d\bm{\theta}+{\cal{O}}\left(\frac{1}{m}\right)$
$\displaystyle=$
$\displaystyle\sum^{K}_{k=1}w_{k}log(w_{k})+\sum^{K}_{k=1}\int_{\Omega_{k}}w_{k}log\left[\frac{p^{g}_{k}(\bm{\theta})}{p(\bm{\theta})}\right]p^{g}_{k}(\bm{\theta})d\bm{\theta}+{\cal{O}}\left(\frac{1}{m}\right)\,.$
(13)

The second term in the last line of (13) can be obtained analytically by
expanding the prior pdf locally around $\hat{\bm{\theta}}_{k}$, and the
Laplace approximation of a K-L divergence, in the case of a multidimensional
multimodal posterior pdf, can be expressed as

$\displaystyle
D(\bar{\bm{y}})=\tilde{D}(\bar{\bm{y}})+{\cal{O}}\left(\frac{1}{m}\right)$

with

$\displaystyle\tilde{D}(\bar{\bm{y}})=$
$\displaystyle\sum^{K}_{k=1}\left[w_{k}log(w_{k})-0.5w_{k}log\left[|\bm{\Sigma}_{k}(\hat{\bm{\theta}}_{k})|\right]-w_{k}h(\hat{\bm{\theta}}_{k})\right]$
$\displaystyle-0.5log(2\pi)^{d}-0.5d\,.$ (14)

Consequently, the corresponding expected information gain can be expressed as

$\displaystyle I=$
$\displaystyle\int_{{\cal{\bm{Y}}}}D(\bar{\bm{y}})p(\bar{\bm{y}})\bar{\bm{y}}=\int_{\bm{\Theta}}\int_{{\cal{\bm{Y}}}}D(\bar{\bm{y}})p(\bar{\bm{y}}|\bm{\theta}_{t})p(\bm{\theta}_{t})d\bar{\bm{y}}d\bm{\theta}_{t}$
$\displaystyle=$
$\displaystyle\int_{\bm{\Theta}}\int_{{\cal{\bm{Y}}}}\tilde{D}(\bar{\bm{y}})p(\bar{\bm{y}}|\bm{\theta}_{t})p(\bm{\theta}_{t})d\bar{\bm{y}}d\bm{\theta}_{t}+{\cal{O}}\left(\frac{1}{m}\right)$
(15)

$\square$

The posterior modes are commonly unknown and need to be obtained numerically.
We find the modes using runs of optimization starting from a few randomized
initial points. Algorithm 1 shows the procedure of computing the expected
information gain via multimodal Laplace approximation.

Algorithm 1 Multimodal Laplace approximation

1: inputs: $p(\bm{\theta})$, $g(\bm{\theta})$, $\bm{\xi}$, $n$

2: draw a sample from the prior distribution $\bm{\theta}^{i}_{t}\sim
p(\bm{\theta})$

3: draw a sample of data from the likelihood function $\bar{\bm{y}}_{i}\sim
p(\bar{\bm{y}}|\bm{\theta}^{i}_{t})$

4: solve the minimization problem:
$\arg\min_{\bm{\theta}}[-log(p(\bm{\theta}|\bar{\bm{y}}_{i}))]$ $n$ times with
distinct initial points, the distinct local optimal solutions are
$\\{\hat{\bm{\theta}}^{i}_{k},k=1,...,K\\}$ (note that $K$ can also vary
w.r.t. $i$)

5: compute $\tilde{D}_{i}$ of (14) using the Hessian
(${\bm{\Sigma}^{i}_{k}}^{-1}$) and the normalized likelihood value
($w_{k}^{i}$)

6: Repeat steps 2-5 $M_{5}$ times and take the sample average of
$\tilde{D}_{i}$: $I_{MLA}=\frac{1}{M_{5}}\sum_{i=1}^{M_{5}}\tilde{D}_{i}$

The cost of computing the expected information gain using Algorithm 1 can be
estimated as the following:

$W_{5}=M_{5}\times(n\times C_{o}+C_{l})\,,$

where $M_{5}$ is the number of samples of $\bm{\theta}_{t}$ drawn from the
prior distribution, $n$ is the number of optimization runs.

In addition, $n$ can be estimated for any given probability that all modes are
captured using Algorithm 1. Such probability can be expressed as the
following:

$\displaystyle 1-P(\cup e_{i})$
$\displaystyle=1-\left[\sum_{i=1}^{K}P(e_{i})-\sum_{i>j}P(e_{i}\cap
e_{j})+\sum_{i>j>k}P(e_{i}\cap e_{j}\cap e_{k})-...\right]$ $\displaystyle\geq
1-\sum_{i}^{K}P(e_{i})\,,$

where $Pr(\cdot)$ is a probability function, $Pr(\cup e_{i})$ is the
probability that at least one mode is missed by the $n$ runs of optimization,
$e_{i}$ is the set of events where the $i^{th}$ mode is missed in all the $n$
runs of optimization. The probability of $e_{i}$ can be written as

$Pr(e_{i})=(1-p_{i})^{n}$

, where $p_{i}$ is the probability that an i.i.d. initial value of $\theta$
converges to the $i^{th}$ mode. It is straightforward to see that the
following relationship exists:

$\displaystyle 1-\sum_{i}^{K}Pr(e_{i})\leq 1-K(1-p)^{n}\,,$ (16)

where $p=min(p_{i})$. Let $\beta=\sum_{i}^{K}Pr(e_{i})$, we can obtain the
following inequality from (16):

$\beta\geq K(1-p)^{n}\,.$

Subsequently, we have the following lower bound of $n$:

$\displaystyle n\geq\frac{log\beta-logK}{log(1-p)}\,,$ (17)

which is independent to the dimension of $\bm{\theta}$. Note that $p$ is
associated to the smallest size of the basins of attractions of all the modes.

##  4 Multimodal nested importance sampling (MNIS)

The method in the previous section approximates a multimodal posterior pdf
using a mixture Gaussian pdf and replaces the inner integration in (2) using
multiple weighted Gaussian integrals. It provides a remedy to the problems
caused a blindly applied conventional Laplace approximation. However, its
error term, ${\cal{O}}\left(\frac{1}{m}\right)$, may not be negligible, when
$m$ is not so large.

Importance sampling is an ideal method to marry the powers of the Laplace
approximation and random sampling. The Laplace approximation based importance
sampling has been demonstrated with salient efficiency and accuracy, when it
is used to calculate a posterior expectation, for example, the expected
information gains used in Bayesian optimal experimental design BECK2018523
Ryan2015 schillings2020 . A multimodal importance sampling method based on the
multimodal Laplace approximations is presented in this section to relieve the
numerical instability (i.e., value of the likelihood ratio be close to
infinity, when samples are taken from the tail of the proposal pdf) induced by
inappropriately using a Gaussian proposal. Theorem 2 summarizes the results of
the change of measure.

###### Theorem 4.1

Assuming that the posterior distribution of $\bm{\theta}$, conditioning on the
synthetic data $\bar{\bm{y}}$, entails $K$ modes, data being modeled by
$\bm{y}_{i}=\bm{g}(\bm{\theta}_{t})+\bm{\epsilon}_{i}$, the third derivatives
of the model $\bm{g}(\bm{\theta}_{t})$ being bounded from above, for any
confidence level $\alpha>0$, there exist $M_{0}$ and $N_{0}$, for all
$M_{6}>M_{0}$ and $N_{6}>N_{0}$, we have

$\displaystyle
Pr\left(\left|I-\frac{1}{M_{6}}\sum^{M_{6}}_{i=1}log\left[\frac{p(\bar{\bm{y}}_{i}|\bm{\theta}_{i})}{\frac{1}{N_{6}}\sum^{N_{6}}_{j=1}p(\bar{\bm{y}}_{i}|\bm{\theta}_{j})\beta^{i}_{j}}\right]\right|>\frac{C_{1}}{\sqrt{M_{6}}}+\frac{C_{2}}{N_{6}}\right)<\alpha\,,$
(18)

with

$\displaystyle
C_{1}=\sqrt{\mathop{\mathbb{V}}\left(log\left[\frac{p(\bar{\bm{y}}|\bm{\theta})}{p(\bar{\bm{y}})}\right]\right)}\,,\,$

and

$\displaystyle
C_{2}=\frac{1}{2}{\mathop{\mathbb{E}}}\left(\mathop{\mathbb{V}}\left[\frac{p(\bar{\bm{y}}|\bm{\theta})\beta(\bm{\theta})}{p(\bar{\bm{y}})}|\bar{\bm{y}}\right]\right)\,,$

where $Pr(\cdot)$ denotes a probability function,
$\beta^{i}_{j}=\frac{p(\bm{\theta}_{j})}{p^{i}_{gm}(\bm{\theta}_{j})}$,
$p^{i}_{gm}$ is the pdf of the Gaussian mixture
$\sum^{K}_{k=1}m_{k}{\cal{N}}\left(\hat{\bm{\theta}}^{i}_{k},\bm{\Sigma}^{i}_{k}\right)$,
$M_{6}$ is the number of samples of the unknown parameters drawn from the
prior distribution, $N_{6}$ is the number of samples of the unknown parameters
drawn from the Gaussian mixture.

###### Proof

The theorem is a direct entailment of Proposition 1 in BECK2018523 . $\square$

Algorithm 3 lists the steps of computing the expected information gain via a
change of Gaussian mixture measure, with the errors stated in Theorem 4.1.

Algorithm 2 Importance sampling based on multimodal Laplace approximation

1: inputs: $p(\bm{\theta})$, $g(\bm{\theta})$, $\bm{\xi}$, $n$

2: draw a sample from the prior distribution $\bm{\theta}^{i}_{t}\sim
p(\bm{\theta})$

3: draw a sample of data from the likelihood function $\bar{\bm{y}}_{i}\sim
p(\bar{\bm{y}}|\bm{\theta}^{i}_{t})$

4: solve problem
$\arg\min_{\bm{\theta}}[-log(p(\bm{\theta}|\bar{\bm{y}}_{i}))]$ $n$ times with
distinct initial points, the distinct local optimal solutions are
$\\{\hat{\bm{\theta}}^{i}_{k},k=1,...,K\\}$ (note that $K$ can also vary
w.r.t. $i$)

5: draw $N_{6}$ i.i.d. samples from Gaussian mixture:
$\bm{\theta}_{j}\sim\sum^{K}_{k=1}w_{k}{\cal{N}}(\hat{\bm{\theta}}^{i}_{k},\bm{\Sigma}^{i}_{k})$

6: compute the inner sample average in (18):
$L_{i}=\frac{1}{N_{6}}\sum^{N_{6}}_{j=1}p(\bar{\bm{y}}_{i}|\bm{\theta}_{j})\beta^{i}_{j}$

7: Repeat steps 2-6 $M_{6}$ times and compute the outer sample average in
(18):
$I_{MNIS}=\frac{1}{M_{6}}\sum^{M_{6}}_{i=1}log\left[\frac{p(\bar{\bm{y}}_{i}|\bm{\theta}_{i})}{L_{i}}\right]$

The computational complexity of Algorithm 2 can be expressed as follows:

$W_{6}=M_{6}\times[n\times C_{o}+C_{l}+N_{6}\times C_{l}]\,.$

Thanks to the change of measure, $C_{2}$ in Theorem 2 is very small and only a
small $N_{6}=C_{2}/(\gamma TOL)$ is needed to control the bias of this
multimodal importance sampler. $\gamma TOL$ is the tolerance on the bias.
Consequently, the following inequality can be true in practice

$N_{6}\times C_{l}\ll n\times C_{o}+C_{l}\,.$

Therefore, the computational cost of the MNIS method is similar to the cost of
the direct multimodal Laplace approximation. We will demonstrate the efficacy
of the proposed methods using numerical examples in the next section.

##  5 Numerical validation

###  5.1 A simple model of multiple factors

First, we prove the concept using a multidimensional model consisting of
quadratic monomials:

$\displaystyle\bm{y}=\bm{g}+\bm{\epsilon}\quad\text{with}\quad\bm{g}=[\xi\theta_{1}^{2}\quad(1-0.5\times\xi)\theta_{2}^{2}\quad\theta_{3}^{2}]^{\top}\,,$
(19)

where $\bm{y}\in\bm{R}^{3}$ is the vector of the synthetic data, $\bm{g}$ is a
three-dimensional deterministic model, whose components consist of distinct
monomials in each dimension,
$\bm{\epsilon}\sim{\cal{N}}\left(\bm{0},\bm{\Sigma}_{e}\right)$ is a three-
dimensional Gaussian random vector with mean zero and covariance matrix
$\bm{\Sigma}_{e}=\bm{I}\sigma^{2}_{e}$. The unknown parameters are uniformly
distributed between $-10$ and $10$ (i.e.,
$\theta_{i}\sim{\cal{U}}(-10,10),\,i=1,2,3$).

Figure 2 demonstrates the superior performance of the MLA and MNIS methods
against the DLMC method (the baseline). The bias of the MLA method is
negligible even with only $10$ sampling points (i.e., $M_{5}=10$ in Algorithm
1). To the contrary, the error of DLMC method is dominated by a bias, which is
controlled by the number of samples in the inner loop (i.e., $N_{1}$ in
equation (3)). It is shown that at least $10^{4}$ inner samples are needed to
reduce the relative error below $10\%$. Two curves of the MLA method produced
by $10$ and $20$ distinct optimization runs are plotted on the left of Figure
2.

An estimation of the required optimization runs using (17) is
$\frac{log(0.1)-log(8)}{log(7/8)}\approx 30$, where we let $\beta=0.1$, $K=8$
and $p=1/K$. It is observed that the MLA method converges to the DLMC method,
when $n=20$. , When $n=10$, the converged result of the MLA method slightly
deviates from the baseline. Note that we always used $M_{1}=10^{3}$, which
effectively removed the statistical variance of the DLMC method. It took
around $390$ s on a MacBook Air to obtain the MLA result for $n=20$ and
$M_{5}=10^{3}$. In comparison, it took around $2000$ s to obtain comparable
result of DLMC for $M_{1}=10^{3}$ and $N_{1}=10^{5}$.

The MNIS results on the right of Figure 2 were produced using $20$ independent
optimization runs. The horizontal axis represents $N_{1}$ and $N_{6}$ of the
DLMC and MNIS methods, respectively. With few number of samples, the MNIS
method is able to eliminate the bias of the approximation, hence,
significantly outperforms the DLMC method. The CPU time of the MNIS method
with $M_{6}=10^{3}$ and $N_{6}=10^{3}$, is very close to that of the MLA
method with $M_{5}=10^{3}$, which produces similar results.

Figure 2: The left figure shows the convergences of the estimated expected
information gain using the MLA method and the DLMC method. $\xi=1$. The right
figure shows the convergence of the MNIS method and the DLMC method. $\xi=1$.

The initial starting point in each run of optimization is obtained according
to the Latin hypercube sampling stein1987large . Figure 3 demonstrates the
effect of the number of optimization runs on the expected information gain and
the averaged number of modes. The expected information gains were computed
using the MLA method with $M_{5}=10^{4}$. It is observed that the relative
error is around $2\%$, when $10$ optimization runs are carried out. The
relative error is smaller than $0.5\%$, when we used more than $20$
optimization runs.

Figure 3: The left figure shows the estimated expected information gains
w.r.t. various numbers of optimization runs. $\xi=1$. The right figure shows
the average number of obtained modes w.r.t. various numbers of optimization
runs. $\xi=1$.

We plot the expected information gains of various experimental setup (i.e.,
$\xi\in[0.05,0.55]$) in Figure 4. The variance of measurement noise is
$\sigma_{e}^{2}=4$ on the left of Figure 4. It is clear that the MLA method is
biased in this case due to the error term in (11), while the MNIS method can
exactly recover the baseline.

Additionally, we consider a scenario, where the observation noise is not
calibrated. Specifically, $\sigma_{e}$ is a random variable. By a little abuse
of notation, the posterior pdf,
$p(\bm{\theta}|\bar{\bm{y}})=\int_{\sigma_{e}}p(\bm{\theta}|\bar{\bm{y}},\sigma_{e})p(\sigma_{e})d\sigma_{e}$,
is now a marginal pdf against the unknown standard deviation of the Gaussian
measurement noise. It can be estimated using
$p(\bm{\theta}|\bar{\bm{y}})=\frac{1}{L}\sum_{l=1}^{L}p(\bm{\theta}|\bar{\bm{y}},\sigma_{e}^{l})\,\text{with}\,\sigma_{e}^{l}\sim
p(\sigma_{e})$. With few modifications to Algorithm 2, the expected
information gain can then be estimated using Algorithm 3.

Algorithm 3 Importance sampling based on multimodal Laplace approximation for
uncalibrated observation noises

1: inputs: $p(\bm{\theta})$, $g(\bm{\theta})$, $\bm{\xi}$, $n$

2: draw two samples from the prior distributions $\bm{\theta}^{i}_{t}\sim
p(\bm{\theta})$ and $\sigma_{e}^{i}\sim p(\sigma_{e})$

3: draw a sample of data from the likelihood function $\bar{\bm{y}}_{i}\sim
p(\bar{\bm{y}}|\bm{\theta}^{i}_{t},\sigma^{i}_{e})$

4: solve problem
$\arg\min_{\bm{\theta}}[-log(\frac{1}{L}\sum_{l=1}^{L}p(\bm{\theta}|\bar{\bm{y}}_{i},\sigma_{e}^{l}))]$
$n$ times with distinct initial points, the distinct local optimal solutions
are $\\{\hat{\bm{\theta}}^{i}_{k},k=1,...,K\\}$ (note that $K$ can also vary
w.r.t. $i$)

5: draw $N_{6}$ i.i.d. samples from Gaussian mixture:
$\bm{\theta}_{j}\sim\sum^{K}_{k=1}w_{k}{\cal{N}}(\hat{\bm{\theta}}^{i}_{k},\bm{\Sigma}^{i}_{k})$

6: compute the inner sample average in (18):
$L_{i}=\frac{1}{N_{6}}\sum^{N_{6}}_{j=1}p(\bar{\bm{y}}_{i}|\bm{\theta}_{j})\beta^{i}_{j}$

7: Repeat steps 2-6 $M_{6}$ times and compute the outer sample average in
(18):
$I_{MNIS}=\frac{1}{M_{6}}\sum^{M_{6}}_{i=1}log\left[\frac{p(\bar{\bm{y}}_{i}|\bm{\theta}_{i})}{L_{i}}\right]$

We assume $\sigma_{e}\sim U(2,4)$ and compute the expected information gain of
the uncalibrated experiments using Algorithm 3. The results are plotted on the
right of Figure 4. The uncalibrated experiments provide less expected
information gains than the calibrated experiments do.

Figure 4: The left figure demonstrates the expected information gains of
various designs. $10^{4}$ and $10^{5}$ samples are used in the outer and inner
loops of the DLMC estimator. $10^{4}$ samples are used in both the outer loop
and the inner loop of the MNIS method. $10^{4}$ samples are used in the MLA
method. The right figure demonstrates the expected information gains, when the
observation noises are uncalibrated.

###  5.2 Identification of sensor networks

We infer the locations of the sensors in a two-dimensional space from some
measurements of the distances between them. The posterior distributions of the
locations are known to have multiple modes Tak2018 Ihler2005 Lan2014 . We
demonstrate the proposed methodologies using two scenarios, where there are
four and six sensors, respectively.

In the first case, the coordinates ($x_{1}$, $z_{1}$, $x_{2}$ and $z_{2}$) of
the first and the second sensors are unknown random variables distributed
uniformly between zero and one, (i.e., $x_{i}\sim{\cal{U}}(0,1)$ and
$z_{i}\sim{\cal{U}}(0,1)$, $i=1,2$). The coordinates of the third and forth
sensors are fixed at $x_{3}=0.5$, $z_{3}=0.3$, $x_{4}=0.3$ and $z_{4}=0.5$,
respectively. Figure 5 shows a set of realizations of the four sensors on the
left.

Figure 5: The left figure shows a network of four sensors. The blue circles
denote the fixed sensors. The black circles denote the sensors of unknown
locations. The arrows connecting sensors indicate measurements of distances.
The right figure shows a network of six sensors.

The data model is

$\bm{y}=\bm{M}\bm{D}+\bm{\epsilon}\,,$

where $\bm{D}$ is a six-dimensional vector:

$\bm{D}=[D_{12},\,D_{13},\,D_{14},\,D_{23},\,D_{24},\,D_{34}]^{\top}\,,$

with $D_{ij}=\sqrt{(x_{i}-x_{j})^{2}+(z_{i}-z_{j})^{2}}$. $\bm{M}$ is a matrix
of $d_{m}$ rows and six columns. It is filled by $0$s and $1$s. We measure the
distances between the first and third sensors, the first and fourth sensors,
the second and third sensors, and the second and fourth sensors. Accordingly,
the matrix $\bm{M}$ reads

$\bm{M}=\begin{bmatrix}0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&1&0&0\\\
0&0&0&0&1&0\end{bmatrix}\,.$

The measurement noises are independent multivariate Gaussian:
$\bm{\epsilon}\sim{\cal{N}}(\bm{0},\bm{\Sigma}_{e}),\bm{\Sigma}_{e}=\bm{I}\sigma^{2}_{e},\sigma^{2}_{e}=0.0005^{2}$.
We then estimate the expected information gain using the fore-mentioned
methods, namely, the MNIS, DLMC and MLA methods. The convergences w.r.t. the
number of samples are plotted in Figure 6. We used $10^{3}$ samples in the
outer loop of the DLMC method and the MNIS method. The DLMC method takes ten
million samples in the inner loop to reach a negligible bias. To the contrary,
the MNIS and MLA methods achieved very accurate results, when few samples
(around $10^{2}$) are used. We used $20$ optimization runs to search for the
modes in the MNIS and MLA methods. It took around $74.28$ s on a MacBook Air
to obtain the MLA result for $n=20$ and $M_{5}=10^{2}$. In comparison, it took
around $2.25\times 10^{6}$ s to obtain comparable result of DLMC for
$M_{1}=10^{3}$ and $N_{1}=10^{7}$. Note that all the methods can be equally
accelerated by parallel computing. Note that the horizontal axis of Figure 6
means differently for different methods (i.e., $N_{1}$ of the DLMC method,
$N_{6}$ of the MNIS methods, and $M_{5}$ of the MLA method). It is a fair
comparison in that we used identical number of samples in the outer loop of
both the DLMC and MNIS methods (i.e., $M_{1}=M_{6}$).

Figure 6: The convergences of the expected information gain estimated using
the MNIS , DLMC and MLA methods, respectively. $\sigma_{e}^{2}=0.0005^{2}$ in
the four-sensor scenario on the left. $\sigma^{2}_{e}=0.16$ in the six-sensor
scenario on the right.

In the second scenario shown on the right of Figure 5, we increased the number
of sensors to six, out of which the first four sensors are independently
distributed according to a multivariate Gaussian distribution (i.e.,
$x_{i}\sim{\cal{N}}(0,1)$, and $z_{i}\sim{\cal{N}}(0,1)$, $i=1,...,4$). The
locations of the fifth and sixth sensors are fixed at $x_{5}=0.5$,
$z_{5}=0.3$, $x_{6}=0.3$ and $z_{6}=0.7$. $\bm{D}$ is a vector of $14$
components, which reads:

$\bm{D}=[D_{12},\,D_{13},\,D_{14},\,D_{15},\,D_{16},\,D_{23},\,D_{24},\,D_{25},\,D_{26},\,D_{34},\,D_{35},\,D_{36},\,D_{45},\,D_{46}]^{\top}\,,$

where $D_{ij}=\sqrt{(x_{i}-x_{j})^{2}+(z_{i}-z_{j})^{2}}$. The matrix $\bm{M}$
is the following:

$\bm{M}=\setcounter{MaxMatrixCols}{12}\begin{bmatrix}1&0&0&0&0&0&0&0&0&0&0&0\\\
0&0&1&0&0&0&0&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0&0&0&0&0\\\
0&0&0&0&1&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&1&0&0&0&0&0\\\
0&0&0&0&0&0&0&1&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0&1&0\\\
0&0&0&0&0&0&0&0&0&0&0&1\end{bmatrix}\,.$

The variance of measurement noise is $\sigma_{e}^{2}=0.16^{2}$. We compare the
convergences of the MNIS, DLMC and MLA methods on the right of Figure 6.
Again, we fixed the number of samples in the outer loop to be $1000$ for the
DLMC and MNIS methods. It is observed that the DLMC method has a relative
error of $10\%$ even with $N_{1}=10^{7}$. Then, we set
$\sigma_{e}^{2}=0.005^{2}$. The MNIS and the MLA methods can achieve accurate
estimations with high confidence using just $100$ samples, while the DLMC
method is erroneous as shown on the left of Figure 7.

Figure 7: The left figure shows the expected information gains computed using
MNIS, DLMC and MLA methods, respectively. $\sigma^{2}_{e}=0.005$. The right
figure shows the maximum expected information gain of various number of
measurements.

Figure 7 on the right shows the maximum expected information gain w.r.t.
various numbers of measurements (i.e., $d_{m}=1,...,14$). It is noted that the
rate of information gain saliently slows down after eight of the fourteen
unknown distances have been measured. Some of the best configurations for
various numbers of measurements are plotted in Figure 8.

Figure 8: The best designs w.r.t. various numbers of measurements.

##  6 Conclusion

Two novel multimodal approaches have been devleoped to estimate the expected
information gain in Bayesian optimal design of experiments. The methods are
developed to consider both calibrated and uncalibrated observation noises. A
multivariate Gaussian mixture distribution is constructed by few direct
optimization runs with randomized initial guesses and weighted local Laplace
approximations. In the first approach, the Gaussian mixture distribution is
used to approximate the posterior pdf. Hence, the K-L divergence can be
computed analytically. And, we can efficiently remove the nested integral in
the DLMC method. In the second approach, the Gaussian mixture is used as the
proposal pdf in an importance sampling scheme, hence, we can remove the
possible Laplace bias in the first approach. It is shown both theoretically
and numerically that the methods have similar computational costs. They can be
several magnitudes more efficient than the DLMC methods. Additionally, the
proposed methodologies are used to efficiently control the error of the
expected information gain, therefore, the error in the numerical solution of
the optimal design is also restrained. We consider applying the novel
multimodal experimental design approach to models based on complex numerical
PDEs as one future research direction.

###### Acknowledgements.

The author would like to thank the associate editor and anonymous referees for
their careful reading and constructive comments.

## References

  * (1) Alexanderian, A., Petra, N., Stadler, G., Ghattas, O.: A-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems with regularized l0-sparsification.  SIAM Journal on Scientific Computing 36(5), A2122–A2148 (2014)
  * (2) Alexanderian, A., Saibaba, A.K.: Efficient D-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems.  SIAM Journal on Scientific Computing 40(5), A2956–A2985 (2018)
  * (3) Atkinson, A., Donev, A., Tobias, R., et al.: Optimum experimental designs, with SAS, vol. 34.  Oxford University Press (2007)
  * (4) Atkinson, A.C., Donev, A.N.: Optimum experimental designs.  Oxford University Press (1992)
  * (5) Beck, J., Dia, B.M., Espath, L.F., Long, Q., Tempone, R.: Fast Bayesian experimental design: Laplace-based importance sampling for the expected information gain.  Computer Methods in Applied Mechanics and Engineering 334, 523 – 553 (2018)
  * (6) Beck, J., Dia, B.M., Espath, L.F., Tempone, R.: Multilevel double loop monte carlo and stochastic collocation methods with importance sampling for Bayesian optimal experimental design.  International Journal for Numerical Methods in Engineering (2020)
  * (7) Bisetti, F., Kim, D., Knio, O., Long, Q., Tempone, R.: Optimal Bayesian experimental design for priors of compact support with application to shock–tube experiments for combustion kinetics.  International Journal for Numerical Methods in Engineering 108, 136 – 155 (2016)
  * (8) Busetto, A.G., Ong, C.S., Buhmann, J.M.: Optimized expected information gain for nonlinear dynamical systems.  In: Proceedings of the 26th Annual International Conference on Machine Learning, ICML ’09, p. 97–104. Association for Computing Machinery, New York, NY, USA (2009)
  * (9) Calvetti, D., Somersalo, E.: Introduction to Bayesian Scientific Computing.  Springer (2007)
  * (10) Carlon, A.G., Dia, B.M., Espath, L., Lopez, R.H., Tempone, R.: Nesterov-aided stochastic gradient methods using Laplace approximation for Bayesian design optimization.  Computer Methods in Applied Mechanics and Engineering 363, 112909 (2020)
  * (11) Chaloner, K., Verdinelli, I.: Bayesian experimental design: A review.  Statistical Science 10(3), 273–304 (1995)
  * (12) Chen, P., Villa, U., Ghattas, O.: Taylor approximation and variance reduction for pde-constrained optimal control under uncertainty.  Journal of Computational Physics 385, 163 – 186 (2019)
  * (13) Dehideniya, M., Overstall, A.M., Drovandi, C.C., McGree, J.M.: A synthetic likelihood-based Laplace approximation for efficient design of biological processes.  arXiv:1903.04168 (2019)
  * (14) Dembo, A., Zeitouni, O.: Large Deviations Techniques and Applications.  Springer (2010)
  * (15) Ghosh, S., Burnham, K.P., Laubscher, N.F., Dallal, G.E., Wilkinson, L., Morrison, D.F., Loyer, M.W., Eisenberg, B., Kullback, S., Jolliffe, I.T., Simonoff, J.S.: Letters to the editor.  The American Statistician 41(4), 338–341 (1987)
  * (16) Goda, T., Hironaka, T., Iwamoto, T.: Multilevel monte carlo estimation of expected information gains.  Stochastic Analysis and Applications 38(4), 581–600 (2020)
  * (17) Huan, X., Marzouk, Y.M.: Simulation-based optimal Bayesian experimental design for nonlinear systems.  Journal of Computational Physics 232(1), 288 – 317 (2013)
  * (18) Ihler, A.T., Fisher, J.W., Moses, R.L., Willsky, A.S.: Nonparametric belief propagation for self-localization of sensor networks.  IEEE Journal on Selected Areas in Communications 23(4), 809–819 (2005)
  * (19) Kullback, S., Leibler, R.A.: On information and sufficiency.  Annals of Mathematical Statistics 22(1), 79–86 (1951)
  * (20) Lan, S., Streets, J., Shahbaba, B.: Wormhole Hamiltonian Monte Carlo.  In: Proceedings of the Twenty-Eighth AAAI Conference on Artificial Intelligence, AAAI’14, pp. 1953–1959. AAAI Press (2014)
  * (21) Laplace, P.S.: Memoir on the probability of the causes of events.  Statistical Science 1(3), 364–378 (1986)
  * (22) Long, Q., Motamed, M., Tempone, R.: Fast Bayesian optimal experimental design for seismic source inversion.  Computer Methods in Applied Mechanics and Engineering 291, 123 – 145 (2015)
  * (23) Long, Q., Scavino, M., Tempone, R., Wang, S.: Fast estimation of expected information gains for Bayesian experimental designs based on Laplace approximations.  Computer Methods in Applied Mechanics and Engineering 259, 24 – 39 (2013)
  * (24) Long, Q., Scavino, M., Tempone, R., Wang, S.: A Laplace method for under-determined Bayesian optimal experimental designs.  Computer Methods in Applied Mechanics and Engineering 285, 849 – 876 (2015)
  * (25) Overstall, A.M., Woods, D.C.: Bayesian design of experiments using approximate coordinate exchange.  Technometrics 59(4), 458–470 (2017)
  * (26) Ryan, E., Drovandi, C., Pettitt, A.: Fully Bayesian experimental design for pharmacokinetic studies.  Entropy 17, 1063–1089 (2015)
  * (27) Ryan, K.J.: Estimating expected information gains for experimental designs with application to the random fatigue-limit model.  Journal of Computational and Graphical Statistics 12(3), 585–603 (2003)
  * (28) Schillings, C., Sprungk, B., Wacker, P.: On the convergence of the Laplace approximation and noise-level-robustness of Laplace-based monte carlo methods for Bayesian inverse problems.  Numerische Mathematik pp. 1–57 (2020)
  * (29) Senarathne, S., Drovandi, C., McGree, J.: A Laplace-based algorithm for Bayesian adaptive design.  Statistics and Computing 30, 1183–1208 (2020)
  * (30) Shaw, J.R., Bridges, M., Hobson, M.P.: Efficient Bayesian inference for multimodal problems in cosmology.  Monthly Notices of the Royal Astronomical Society 378(4), 1365–1370 (2007)
  * (31) Stein, M.: Large sample properties of simulations using Latin hypercube sampling.  Technometrics 29(2), 143–151 (1987)
  * (32) Tak, H., Meng, X.L., van Dyk, D.A.: A repelling-attracting Metropolis algorithm for multimodality.  Journal of Computational and Graphical Statistics 27(3), 479–490 (2018)
  * (33) Tarakanov, A., Elsheikh, A.H.: Optimal bayesian experimental design for subsurface flow problems.  Computer Methods in Applied Mechanics and Engineering 370, 113208 (2020)
  * (34) Tierney, L., Kadane, J.B.: Accurate approximations for posterior moments and marginal densities.  Journal of the American Statistical Association 81(393), 82–86 (1986)
  * (35) Tierney, L., Kass, R.E., Kadane, J.B.: Fully exponential Laplace approximations to expectations and variances of nonpositive functions.  Journal of the American Statistical Association 84(407), 710–716 (1989)
  * (36) Titterington, D.M.: Optimal design: Some geometrical aspects of D-optimality.  Biometrika 62(2), 313–320 (1975)
  * (37) Verdinelli, I.: A note on Bayesian design for the normal linear model with unknown error variance.  Biometrika 87(1), 222–227 (2000)
  * (38) Visser, A., Slamet, B.A.: Including communication success in the estimation of information gain for multi-robot exploration.  In: 2008 6th International Symposium on Modeling and Optimization in Mobile, Ad Hoc, and Wireless Networks and Workshops, pp. 680–687 (2008)