A layered multiple importance sampling scheme for focused optimal Bayesian experimental design
Chi Feng Youssef M. Marzouk


Information Geometric Complexity of Entropic Motion on Curved
Statistical Manifolds under Different Metrizations of Probability Spaces
    Steven Gassner and Carlo Cafaro
    December 30, 2023
======================================================================================================================================




  We develop a new computational approach for “focused” optimal
  Bayesian experimental design with nonlinear models, with the goal of
  maximizing expected information gain in targeted subsets of model
  parameters. Our approach considers uncertainty in the full set of
  model parameters, but employs a design objective that can exploit
  learning trade-offs among different parameter subsets.
  We introduce a new layered multiple importance sampling scheme
  that provides consistent estimates of expected information gain in
  this focused setting. This sampling scheme yields significant
  reductions in estimator bias and variance for a given computational
  effort, making optimal design more tractable for a wide range of
  computationally intensive problems. 





optimal experimental design, Bayesian inference, Monte Carlo methods,
multiple importance sampling, expected information gain, mutual
information








§ INTRODUCTION

Data acquired from experiments are essential for developing and refining models. Acquiring data can be a costly process, so it is important to carefully choose experimental conditions that maximize the value of the acquired data, as defined by some optimal experimental design criterion <cit.>. In general, the design criterion should depend on what one intends to do with the data acquired from an experiment.

In this paper, we consider the optimal design of experiments where the observed data are used to infer uncertain model parameters. This setting includes models that may contain a large number of parameters, but where the user may care only about inference of a targeted subset. For example, in the context of chemical kinetics, experimental outcomes might be influenced by multiple rate constants and thermodynamic parameters, but only the rates of certain subsets of reactions may be of interest. Another example is the Kennedy & O'Hagan model calibration framework <cit.>, where users may favor experimental designs that focus on reducing uncertainties in `physical' parameters over reducing uncertainty in `tuning' parameters. To this end, we employ a design criterion that favors experiments maximizing the expected information gain in user-selected parameters of interest; following <cit.>, we refer to this approach as focused optimal experimental design. Our key contribution is a new multiple importance sampling framework that makes focused experimental design tractable.

We formulate our experimental design problems in a Bayesian setting, where uncertain parameters are endowed with probability distributions and parameter inference can naturally be performed when the data are noisy and limited in number <cit.>. We adopt a decision-theoretic approach, following <cit.>, wherein an optimal design maximizes the expected utility of an experiment. As suggested by <cit.>, our utility function is chosen to be the information gain due to an experiment, which is equivalent to the Kullback–Leibler divergence from the prior to the posterior—or from selected prior marginals to the corresponding posterior marginals. This nonlinear Bayesian optimal design framework is a generalization of Bayesian 
D-optimal design, making no assumptions on the linearity of the forward model or the Gaussianity of the prior and noise distributions <cit.>.

For models where the observables depend nonlinearly on the parameters of interest, exact evaluation of this design criterion is in general intractable, as it requires integration over the possible outcomes of the experiment and the space of model parameters. A nested Monte Carlo estimator <cit.> of the expected information gain (jointly over all model parameters) is proposed in <cit.>, but has significant limitations due to computational expense. Alternatives to this nested Monte Carlo estimator include kernel estimators of mutual information 
as proposed in <cit.> and applied in <cit.>, but these approaches involve density estimation and thus require careful parameter choices to control bias and variance, along with large sample populations.
Laplace approximations to the posterior distribution <cit.> or computationally efficient surrogate models <cit.> can result in more tractable design criteria, but are unsuitable when these approximations break down, e.g., as the posterior distribution departs from normality or when a sufficiently accurate surrogate cannot be constructed. In general, these approximations induce asymptotic bias and do not yield consistent estimators of the expected information gain. Of course, it is possible to use Laplace approximations to construct consistent importance sampling estimators, as recently proposed in <cit.> and discussed below; these techniques are complementary to the Monte Carlo methods developed here. Further reviews of modern computational methods and approximations for decision-theoretic Bayesian design can be found in <cit.>. 

Other recent experimental design formulations avoid computing the expected information gain criterion altogether, and instead begin with the Bayesian A-optimality criterion <cit.> for linear problems, extending it to nonlinear settings by (i) using Laplace approximations of the posterior for tractability; and (ii) averaging over many possible realizations of the data <cit.>. These approaches are similar to nonlinear design criteria that average functionals of the expected or observed Fisher information matrix, as described in <cit.>. 
Though these alternatives provide computational scalability and can be rationalized in weakly nonlinear problems, they remain only approximations of a decision-theoretic expected utility. Also, while D-optimal design and its generalization to expected information gain are invariant under reparameterization of the underlying model, A-optimal design and its extensions discussed above are not <cit.>.


This paper proposes a consistent, yet tractable Monte Carlo estimator of the expected information gain, not limited by linear 
 approximations. Estimating the expected information gain requires estimating posterior normalizing constants for many realizations of the data. Yet estimating the posterior normalizing constant—even for a single realization of the data—is a challenging task. In addition, our criterion for focused experimental design introduces marginalization over nuisance parameters <cit.>, further reducing the efficiency and accuracy of existing Monte Carlo approaches. While many schemes <cit.> have been developed to estimate posterior normalizing constants, they only target a single posterior distribution. Estimating posterior normalizing constants over many realizations of the data, however, allows for the sharing of information between estimators. We use this idea to develop a layered multiple importance sampling (LMIS) estimator of the expected information gain, which iteratively improves the biasing distributions used to estimate posterior normalizing constants. This is accomplished by optimally combining samples and cached model evaluations from earlier iterations using multiple importance sampling <cit.>, yielding consistent and more accurate estimators.

Our approach has roots in the simpler importance sampling scheme of <cit.>. A useful generalization of the latter, recently proposed by <cit.>, involves computing a Laplace approximation of the posterior for each realization of the data, and using this approximation to construct a biasing distribution to estimate each posterior normalizing constant. The resulting expected information gain estimators are also consistent, but require solving an optimization problem and computing a Hessian for each realization. In contrast, the LMIS scheme developed here emphasizes recycling information from one realization of the data to the next. Within this framework, many different biasing schemes can be used. The particular biasing schemes employed here avoid repeated optimization and do not require gradients or Hessians of the forward model or log-likelihood. LMIS is thus entirely non-intrusive, in the sense of requiring only pointwise evaluations of the likelihood function.

We demonstrate using both linear and nonlinear numerical examples examples that our LMIS estimator offers multiple order-of-magnitude reductions in mean squared error over comparable Monte Carlo methods using the same number of model evaluations, while remaining asymptotically exact. This new estimator enables efficient estimation of expected information gain for both joint and the more challenging case of focused experimental design.


This paper is organized as follows. In <ref>, we introduce the Bayesian optimal experimental design problem and formulate the focused information gain criterion. We illustrate this formulation using a linear-Gaussian toy problem—which has closed-form expressions of the expected information gain, and a clear trade-off between inferring all the model parameters and inferring a subset thereof. In <ref> we derive a nested Monte Carlo estimator of our focused design criterion, allowing for a generic choice of importance sampling scheme in its inner loops; we then analyze the asymptotic bias and variance of this estimator. <Ref> describes our layered multiple importance
 sampling scheme. <Ref> evaluates the numerical performance of the LMIS estimator on linear problems of varying dimension, and in <ref> we apply the LMIS algorithm to the optimal design of a Mössbauer spectroscopy experiment, which is nonlinear in its parameters. 

















 


§ OPTIMAL EXPERIMENTAL DESIGN


The Bayesian paradigm treats unknown parameters as random variables. Let (Ω,ℱ,ℙ) be a probability space, where Ω is a sample space, ℱ is a σ-field, and ℙ is a probability measure on (Ω,ℱ). Let the real-valued random vector : Ω→^n_θ + n_η represent the uncertain parameters of a statistical model, and let it be partitioned coordinate-wise, (, ), into : Ω→^n_θ and : Ω→^n_η, which denote the parameters of interest and nuisance parameters, respectively. Both  and  are parameters to be conditioned on experimental data.
The random variable  is associated with a prior measure μ on ℝ^n_θ+n_η such that μ(A)=ℙ(^-1(A)) for A∈ℝ^n_θ + n_η. We can then define p() = dμ/d = p(, ) as the joint density of  and  with respect to the Lebesgue measure d. For the present purposes, we will assume that this density always exists. Similarly, the data  can be treated a real-valued random variable : Ω→^n_y. Finally, let ∈𝒟⊆^n_d denote the design variables, or experimental conditions. Hence, n_θ is the number of uncertain parameters of interest, n_η is the number of nuisance parameters, and n_d is the number of design variables. In the Bayesian setting, upon observing a realization of the data  obtained by performing an experiment under conditions , we can update our knowledge about the model parameters using Bayes' rule:

    p(, |,)_posterior=p( | , , )_likelihood p(, )_prior/p(| )_evidence.

Here p(, ) is the prior density, p(| , , ) is the likelihood function associated with our statistical model, p(,  | , ) is the posterior density, and p(| ) is the evidence or posterior normalizing constant. It is reasonable to assume that prior knowledge on  and  (i.e., the prior measure μ) does not depend on the experimental design, and thus we employ the simplification p(,| )=p(,).

The likelihood function often combines a deterministic forward model with a statistical model for measurement and/or model errors. For example, the observations might be modeled as additive Gaussian perturbations of some deterministic model predictions (, , ):

    = (, , ) + ,   ∼𝒩(0, Γ_ϵ).

Here,  is the forward model; we will assume that it is (in general) nonlinear and that its evaluations are the dominant computational cost of evaluating or simulating from the likelihood. Thus, while evaluations of (, , ) for new values of the parameters or design variables are expensive, it is cheap to evaluate p(|, , ) for different values of , given fixed values of , , and . This assumption is representative of many science and engineering problems, where  might represent the solution of a set of differential equations and where the noise model is relatively simpler.




 §.§ Expected utility framework

Following the decision-theoretic approach of Lindley <cit.>, an objective for optimal experimental design should have the following general form:

    U()    = ∫_𝒴∫_∫_u(,,,) p(,,|)   ,

where u(,,,) is a utility function representing the value of observing a particular outcome  of the experiment under design , when the parameters are  and . We do not know the outcome of the experiment a priori, however; nor do we know the true values of the parameters  and . So we instead consider the expected utility U(), where the expectation is taken over the joint distribution of , , and , and where ⊆ℝ^n_θ, ⊆ℝ^n_η and 𝒴⊆ℝ^n_y are the ranges of , , and , respectively.

To focus attention on inference of the parameters , we choose u(,,,) to be the Kullback–Leibler (KL) divergence (i.e., relative entropy) from the marginal prior on  to the marginal posterior of , where the marginalization is over the nuisance parameters . In other words, the divergence is from the marginal prior on the parameters of interest,

    p_θ() = ∫_p(,) d,

to the marginal posterior on the parameters of interest,

    p_θ(|,) =
      ∫_p(,| ,)d,

and hence our utility function is

    u(,,,)    (p_θ(| ,)   p_θ()) 
       = ∫_ p_θ(| ,) log[p_θ(| ,)/p_θ()]d  = u(, ).

Above,  is simply a dummy variable of integration. The final equality reflects the fact that the KL divergence does not actually depend on the values of  and . Thus we can write the expected utility (<ref>) of a design as an expectation over the data only:

    U()    =  ∫_∫_∫_u(,) p(,,|)ddd
       =  ∫_ u(,) p(| ) d.

Now, substituting the utility (<ref>) into the objective function (<ref>), we obtain

    U()   = ∫_∫_ p_θ(| ,) log[p_θ(| ,)/p_θ()]d p(| )d
       = _| [(p_θ(| ,)  p_θ())],

which means that the expected utility U() is the expected information gain in the parameters of interest , irrespective of the information gain in the nuisance parameters . Intuitively, the expected information gain describes the average reduction in entropy from the prior marginal of  to the posterior marginal of , where the average is taken over all possible experimental outcomes , weighed according to their prior predictive distribution.  U() is equivalent to the mutual information
between the parameters of interest  and the data , given a particular design  <cit.>.






 §.§ Linear–Gaussian example


To illustrate our formulation for focused optimal experimental design, we first apply it to a toy problem that can be solved analytically. Consider the following experiment with a linear model,

    [ y_1; y_2 ]=[   d   0;   0 1-d ][ θ; η ]+[ ϵ_1; ϵ_2 ],

additive independent noise ϵ_1,ϵ_2∼(0,σ_ϵ^2), and a design space d∈[0,1]. The prior distributions on θ and η are chosen to be independent standard normals, θ∼(0,1) and η∼(0,1). The posterior θ,η |,d ∼(_post, _post) is a multivariate normal with mean and covariance available in closed form.
The marginal posterior p( θ  | ,d) is also Gaussian and follows immediately.

The structure of this model encodes a simple trade-off between learning about θ and learning about η. As d increases, we obtain a larger signal-to-noise ratio in the measurement of θ while the signal-to-noise ratio associated with measuring η decreases. This trade-off can be visualized by observing the anisotropy of the posterior distribution when varying d for a fixed value of , as illustrated in <ref>. The marginals at the bottom of the figure suggest that the design minimizing uncertainty in θ is simply d=1.



Since we are working with a linear-Gaussian model, we can obtain analytical expressions for the 
expected information gain in (θ, η) jointly,

    U_θ,η(d)   =_|d[(p (θ,η | ,d ) ‖ p(θ,η))]=1/2log[((1-d)^2+^2)(d^2+^2)/^4],and in the parameter of interest θ alone,
    U_θ(d)   =_|d[(p(θ|,d) ‖ p(θ))]=1/2log[1+d^2/^2].

The expected information gains for both cases are plotted on the right of <ref>.

To find the optimal experimental design, the expected utility must be maximized over the design space 𝒟 = [0,1]. When the objective is to maximize information gain in the parameters of interest, we have U = U_θ and the optimal design is d_θ^∗ = 1; this is consistent with the θ marginals shown in <ref>. On the other hand, when the objective is to maximize information gain in all uncertain parameters, we have U = U_θ, η and the optimal design becomes d_θ,η^∗ = 1/2. Notice that maximizing the expected information gain in <ref> is equivalent to minimizing the determinant of the posterior covariance matrix, i.e., Bayesian D-optimal design <cit.>.
We will revisit more complex versions of this example in <ref>, where we benchmark the convergence of various Monte Carlo estimators of (<ref>).




 



§ MONTE CARLO ESTIMATION OF THE EXPECTED INFORMATION GAIN


In general, the expected information gain (EIG)—whether in the parameters  alone or jointly in (,)—does not have a closed form and must be approximated numerically, e.g., by a Monte Carlo estimate. In this section, we derive a nested Monte Carlo estimator of the focused EIG (<ref>), analyze its properties, and then propose a layered multiple importance sampling (LMIS) estimator that has better properties. To make the subsequent developments clear, we note that the nested estimators employed here comprise an outer loop that involves taking an expectation over the prior distribution and an inner loop that involves computing certain posterior normalizing constants. <Ref> presents a general format for nested estimators of (<ref>) that employ importance sampling in their inner loops, but deliberately leaves the choice of biasing distributions open. <Ref> analyzes the bias and variance of these estimators; this analysis also suggests what good biasing distributions ought to be. <Ref> discusses the convergence rate of the nested estimators and the optimal asymptotic scaling of inner loop and outer loop sample sizes.
<Ref> then presents a way of adaptively identifying biasing distributions via multiple importance sampling, which leads to our new LMIS estimator. LMIS can offer significant increases in performance over simpler and non-adaptive choices of biasing, while preserving asympotically exact estimation of the focused EIG.




 §.§ A nested importance sampling estimator


The EIG as written in (<ref>) contains an inner loop of integration over the posterior distribution, and also requires evaluating the normalized posterior density of ; both are computationally undesirable. By applying Bayes' rule to the terms involving the posterior density, we can rewrite (<ref>) so that integration is instead performed over the prior distribution of , , and :

    U()     = ∫_ ( ∫_log[p(,)/p()] p(,)  )  p() 
        =
        ∫_∫_∫_{log[p(,)]-log[p()]} p(,, ) p(, )   .

This form of the focused EIG can be approximated by the following Monte Carlo estimate:

    U() = 1/N∑_i=1^N log[p(^(i)^(i),)_conditional
    likelihood]-log[p(^(i))_marginal
    likelihood],

where the samples {^(i), ^(i)}_i=1^N are drawn i.i.d. from the prior as follows: for each i, draw (^(i), ^(i)) ∼ p(,) from the joint prior on the parameters; then draw ^(i)∼ p(^(i),^(i),). The sample ^(i) is necessary to realize a value of ^(i) but does not appear explicitly in (<ref>).

The summand in (<ref>) contains the posterior normalizing constant p(^(i)) for a given realization of the data ^(i) and design parameters , known as the marginal likelihood, along with a partial posterior normalizing constant p(^(i)^(i),), which we call the conditional likelihood. The marginal likelihood typically does not have an analytical form, but can be approximated by an importance sampling estimate for each ^(i):

    p(^(i))     =  ∫_×p(^(i),,) p(,) (̣,̣)̣
       ≈1/M_1∑_j=1^M_1 p(^(i)^(i,j),^(i,j),) w_marg^(i,j)p(^(i)),

where the samples {^(i,j), ^(i,j)}_j=1^M_2 are drawn i.i.d. from a suitable biasing distribution q_marg^(i)(, ) 
with corresponding importance weights:

    w_marg^(i,j) = p(^(i,j),^(i,j))/q_marg^(i)(^(i,j),^(i,j)).

Similarly, the conditional likelihood evaluated at ^(i), p(^(i)^(i),), can be estimated as follows:

    p(^(i)^(i),)     = ∫_ p(^(i)^(i),,) p(^(i)) 
       ≈1/M_2∑_k=1^M_2 p(^(i)^(i),^(i,k),) w_cond^(i,k)p(^(i)^(i),),

where samples {^(i,k)}_k=1^M_2 are drawn i.i.d. from a biasing distribution q_cond^(i)(), with importance weights

    w_cond^(i,k) = p(^(i,k)^(i))/q_cond^(i)(^(i,k)).

Substituting the above importance sampling estimates of the marginal and conditional likelihoods into (<ref>), we arrive at a nested Monte Carlo estimator for the EIG in :

    U()=1/N∑_i=1^N {log[ 1/M_2∑_k=1^M_2 p(^(i)^(i),^(i,k),) w_cond^(i,k)] 
        
            -log[ 1/M_1∑_j=1^M_1 p(^(i)^(i,j),^(i,j),) w_marg^(i,j)] }.


So far, we have not explicitly specified the biasing distributions { q_marg^(i)}_i=1^N and { q_cond^(i)}_i=1^N. In the next sections, we will identify optimal biasing distributions and describe methods to approximate them.




 §.§ Bias and variance of the nested estimator



Choosing an appropriate biasing distribution requires understanding the bias and variance of the focused EIG estimator (<ref>). Using the second-order delta method (see <ref>), we find that the bias of U() is, to leading order:

    [U()-U()]    ≈1/M_1_[σ^2_marg(,)/p^2(|)] - 1/M_2_,[σ^2_cond(,,)/p^2(|,)]
        = C_1()/M
            _1-C_2()/M_2,

and that its variance is, to leading order:

    [U()]    ≈1/N_,[lnp(|,)/p(|)]+1/NM_1_[σ^4_marg(,)/p^4(|)] + 1/NM_2_,[σ^4_cond(,,)/p^4(|,)]
        = D_3()/N+D_1()/NM_1+D_2()/NM_2.

The higher-order terms in the bias (<ref>) are 𝒪(1/M_1^2) and 𝒪(1/M_2^2), while the higher-order terms in the variance (<ref>) are 𝒪(1/(NM_1^2)) and 𝒪(1/(NM_2^2)); see <ref>  for details.  
Also, σ_marg^2(,) and σ_cond^2(,,) are variances of the summands in the marginal and conditional likelihood estimators (<ref>) and (<ref>), respectively:

    
    σ_marg^2(,)                                                                                                                                                                                               _,[p(|,,)p(,)/q_marg^()(,)],     (,)∼ q_marg^(),
    σ_cond^2(,,)                                                                                                                                                                                       _[p(|,,)p(|)/q_cond^(, )()],                                                           ∼ q_cond^(, ).



From (<ref>), we see that the bias is controlled by the variance of the inner marginal and conditional likelihood estimators, which can be reduced by choosing better biasing distributions q_marg^() and q_cond^(,) or by increasing the number of inner Monte Carlo samples M_1 and M_2.[
Here and in (<ref>)–(<ref>), we write the inner biasing distributions as q_marg^() and q_cond^(, ) to indicate that they can depend arbitrarily on the outer values of  and possibly . To recover the simpler notation of <ref>, where i indexes the outer loop, simply put q_marg^(i) =  q_cond^(^(i) ) and q_cond^(i) =  q_cond^ (^(i), ^(i) ).]
The bias can take both positive and negative values since it is the difference of two positive quantities. This property contrasts with the nested Monte Carlo estimator of unfocused EIG in <cit.>, where the bias only takes positive values. The latter corresponds to the case of no nuisance parameters (n_η = 0); then the conditional likelihood p(|,) (<ref>) can be evaluated exactly, so σ_cond^2(, , ) is effectively zero and (<ref>) is always positive.

The variance of U() decreases with the number of outer samples as 1/N. The first term, D_3(), depends only on the variability of the conditional and marginal likelihoods given (, ) drawn from the prior. The second and third terms, D_1() and D_2(), depend explicitly on the variance of the inner estimators. Thus, reducing the variance of the inner estimators will also reduce the variance of U().

In order to minimize the bias and variance of U(), we therefore need to find biasing distributions q_marg^() and q_cond^(,) that minimize the variance of the importance sampling estimates of the marginal and conditional likelihood. From <ref> and <ref>, we see that the using the joint and conditional posteriors as biasing distributions (with normalized densities) would result in zero-variance estimators. Current approaches for evaluating the marginal likelihood in this setting, such as <cit.>, instead use the prior distribution p(,) as a biasing distribution, i.e.,

    q_marg^()(,)=p(,)  and  
        q_cond^(,)()=p(|)  ,

so that the importance weights w_marg and w_cond in (<ref>) and (<ref>) are identically one. However, when the joint or conditional posterior distributions are more concentrated than their corresponding prior distributions, e.g., when the data ^(i) are informative, using the prior distribution as a biasing distribution can result in very high-variance  estimates of the marginal and conditional likelihoods, which ultimately lead to large bias (and variance) when estimating the expected information gain. (We will demonstrate this issue in <ref>.) To avoid this problem, a more sophisticated approach is needed to find biasing distributions that more closely resemble the joint and conditional posteriors.



 §.§ Sample sizes and scaling


The nested estimator (<ref>) is consistent (see <cit.>) and in particular, its bias and variance converge to zero as M_1, M_2,  N →∞. Since all three sample sizes must increase, it is natural to ask whether there is an optimal allocation of inner (M_1, M_2) and outer (N) samples for a given computational budget W = N(M_1 + M_2).

For simplicity, let M_1 = M_2 = M so that W = 2MN and define C (C_1- C_2)^2>0, where we drop dependence on  to reduce clutter.  Let α^2=2M/N be the ratio between the number of inner and outer Monte Carlo samples. The value of α and the corresponding ratio α^2 that minimize the leading-order MSE are

    α_∗=2  ( C/D_3 )^1/31/W^1/6  ⇒  α_∗^2= 4  ( C/D_3 )^2/31/W^1/3.

This result suggests that the ratio between the number of inner and outer samples should decrease slowly as the computational budget increases, in particular scaling as W^-1/3. The ratio C/D_3 reflects the relative magnitudes of leading-order terms in the bias and variance; if the bias dominates, one should choose a relatively larger number of inner loop samples. A detailed derivation of the optimal scaling and a numerical study of the asymptotic behavior are presented in <ref>.

It is important to emphasize, however, that all these results are asymptotic in M, N →∞. For finite M and N, higher-order terms in the bias and variance may come into play. Also, the constants in (<ref>) are in practice difficult to estimate, and an incorrect value—one that does not capture the relative magnitudes of the bias and the variance—may cause the α∼ W^-1/6 scaling to perform worse than a fixed M/N ratio for finite sample sizes. 



 §.§ Layered multiple importance sampling



Now we introduce our new algorithm. 



  §.§.§ Overall algorithm

To implement the focused EIG estimator (<ref>), one should identify good biasing distributions q_marg^(i) and  q_cond^(i) at each outer-loop iteration i=1, …, N. An ideal approach might involve directly sampling from the posterior corresponding to each ^(i), but doing so is computationally taxing—and, moreover, should only be done when evaluating the corresponding normalized posterior density (to avoid instabilities of the harmonic mean estimator <cit.>), which requires knowledge of the quantity we actually wish to estimate.
Instead, we will restrict our choice of biasing distributions to a parametric family with known density and use approximations of the joint and conditional posteriors, within this family, as biasing distributions. We describe our approach below using multivariate t-distributions; since these distributions are heavy-tailed, they are broadly useful in importance sampling <cit.>. As we shall emphasize later, however, our approach can easily be used with other parametric families.


To identify a multivariate t-distribution that is good for biasing relative to a given posterior, it is useful to estimate the posterior mean and covariance. Ideally, we would like to estimate these moments without drawing additional samples that require costly forward model evaluations. Our key idea is as follows: if we proceed sequentially  through the outer loop of (<ref>), indexed by i, then we can use the samples drawn from each inner-loop biasing distribution in the previous steps {1, …, i-1} to help characterize the posterior corresponding to the current ^(i).
To reuse these existing samples, we consider a multiple importance sampling approach <cit.>, which combines the previously-drawn samples into a single mixture distribution 
that is used to estimate the mean _post^(i) and covariance _post^(i) of the posterior distribution p(,|^(i),). To this end, we first introduce self-normalized importance sampling estimators of the posterior mean and covariance:

    _post^(i)    =∑_ℓ=1^L w_post^(i,ℓ)[ ^(i,ℓ); ^(i,ℓ) ] 
    _post^(i)    =∑_ℓ=1^L w_post^(i,ℓ)([ ^(i,ℓ); ^(i,ℓ) ]-_post^(i))([ ^(i,ℓ); ^(i,ℓ) ]-_post^(i))^

where the samples {(^(i,ℓ),^(i,ℓ))}_ℓ=1^L are drawn from a biasing distribution q_post^(i)(,) and the unnormalized importance weights are given by

    w_post^(i,ℓ)∝p(^(i,ℓ),^(i,ℓ)|^(i),)^joint posterior/q_post^(i)(^(i,ℓ),^(i,ℓ))∝p(^(i)|^(i,ℓ),^(i,ℓ),)^likelihood p(^(i,ℓ),^(i,ℓ))^prior/q_post^(i)(^(i,ℓ),^(i,ℓ)),

which can be rescaled to obtain the self-normalized importance weights

    w_post^(i,ℓ) = w_post^(i,ℓ)/∑_ℓ=1^Lw_post^(i,ℓ).


Now, let q_post^(i) be a mixture distribution whose components are distributions from which we have already drawn samples—e.g., the prior distribution (since we can easily draw all N outer loop samples in advance) and the inner-loop biasing distributions from earlier iterations:

    q_post^(i)(,)=N/Lp(,)+M_1/Lq_marg^(1)(,) + M_1/Lq_marg^(2)(,)+⋯+M_1/Lq_marg^(i-1)(,).

Instead of drawing new samples from this mixture distribution, we simply realize that the existing samples are a draw from q_post^(i):

    {(^(i,ℓ),^(i,ℓ))}_ℓ=1^L_from q_post^(i) = {(^(i),^(i))}_i=1^N_from prior   ⊔{(^(1,j),^(1,j))}_j=1^M_1_from q_marg^(1)⊔{(^(2,j),^(2,j))}_j=1^M_1_from q_marg^(2)
       ⊔…⊔{(^(i-1,j),^(i-1,j))}_j=1^M_1_from q_marg^(i-1)

Moreover, because the forward model  has already been evaluated each sample of (, ) in this empirical mixture, evaluating the likelihood at the current data point ^(i) will be very cheap. This process results in the incremental enrichment of the samples used to estimate the posterior mean _post^(i) and covariance _post^(i). Pseudocode for this process is given in Algorithm <ref>.

The inner estimators for the marginal and conditional likelihood still retain their original forms, (<ref>) and (<ref>), but now the biasing distributions q_marg^(i) and q_cond^(i) are set according to the estimated posterior moments. Thus the biasing distribution for estimating the marginal likelihood is

    q_marg^(i) (,)     =t_ν(_post^(i), _post^(i)).

Since it is helpful for biasing distributions to have heavier tails, we typically choose small ν, e.g., ν = 2.5. One possible choice of q_cond^(i) is another multivariate t-distribution derived from conditionals of the posterior Gaussian approximation:

    q_cond^(i)()     =t_ν(_+_^T_^-1(^(i)-_), _-_^T_^-1_).

where we used the block-matrix notation

    _post^(i)    = [ _; _ ],  _post^(i) = [   _   _; _^T   _ ].

Alternatively, one could use the conditional of q_marg^(i) directly; note that this is just another multivariate t distribution, very similar to (<ref>), but with a scaling of the scale matrix _|_-_^T_^-1_ and a modification of the degrees of freedom ν <cit.>. Another possible choice is simply the –marginal of  q_marg^(i). This last choice is more diffuse, since it does not depend on ^(i), but may be useful in situations where the preceding approximations to the true posterior conditional are poor.


Note, finally, that the importance sampling estimates for the marginal and conditional likelihood involve drawing new samples {(^(i,j),^(i,j))}_j=1^M_1 and {(^(i,k))}_k=1^M_2 from the biasing distributions just constructed. Since we know the normalized densities of these biasing distributions, the resulting importance sampling estimates are unbiased at finite M_1, M_2. 

We call this complete strategy
layered multiple importance sampling, with two layers:

    
  * A biased layer (with bias due to self-normalization) to estimate posterior moments, via multiple importance sampling with sample reuse;
    
  * An unbiased layer to estimate the marginal and conditional likelihoods, using biasing distributions selected from a parametric family according to the estimated posterior moments.




  §.§.§ Pruning the mixture biasing distribution

Although layered multiple importance sampling does not require additional model evaluations beyond a standard nested scheme, if we include all the biasing distributions from earlier iterations in (<ref>), the number of likelihood evaluations—using cached model evaluations—grows as O(N^2M_1^2), which can result in a nontrivial computational overhead. 
Instead, we can include only the most “useful” samples in the empirical mixture distribution. Since we do not know the posterior beforehand, we can use the (^(i),^(i)) which gave rise to the data ^(i) as our best guess for where the posterior p(,|^(i),) is most concentrated. This selection criterion allows us to define an index set

    𝒥^(i)={m :  q_marg^(m)(^(i),^(i)) > p(^(i),^(i)),  m<i}.

where we only include distributions that have a higher probability density on the sample (^(i),^(i)) than the prior. Then, the mixture density (<ref>) can be modified as

    q_post^(i)(,) = N/L p(,) + M_1/L∑_m∈𝒥^(i) q_marg^(m)(,),

and the samples drawn from this mixture distribution 
q_post^(i) can be written as a disjoint union of samples drawn previously from the prior and biasing distributions:

    {(^(i,ℓ),^(i,ℓ))}_ℓ=1^L_from q_post^(i) = {(^(i),^(i))}_i=1^N_from prior⊔_m∈𝒥^(i){(^(i,m,j),^(i,m,j))}_j=1^M_1_from q_marg^(m),

where the total number of samples is

    L=N+M_1|𝒥^(i)|.

The specification of an index set 𝒥^(i) allows some generality: instead of the strategy in (<ref>), one could include all the mixture components by setting 𝒥^(i)={1,2,…,i-1}, or one could avoid incremental enrichment by setting 𝒥^(i)={∅}, so that q_post^(i) becomes the prior distribution p(,) and we reuse only the N samples drawn from the prior.

Our overall strategy—layered multiple importance sampling with mixture pruning—is illustrated via the flowchart in <ref>.






  §.§.§ Ordering samples

Since the layered multiple importance sampling approach is sequential by construction, a natural question is whether the ordering of the prior samples {(^(i),^(i),^(i))}_i=1^N will affect the behavior of the overall estimator. Since we draw these samples before evaluating the inner estimators, we have the freedom to specify a particular ordering. Intuitively, a good ordering is one that maximizes incremental enrichment—such that samples drawn from earlier biasing distributions will be reused in many subsequent mixture distributions. This goal suggests an ordering that prioritizes the posterior distributions p(,|^(i),) that overlap the most with all other posterior distributions. However, since the N different posterior distributions are not characterized before the ordering is chosen, we suggest using the prior density of p(^(i),^(i)) as a surrogate. With sufficient data and under appropriate consistency conditions, the posterior distributions p(,|^(i),) will be concentrated near the ^(i),^(i) that gave rise to the ^(i); more generally, the centers of the posteriors will (collectively) be distributed similarly to the prior. Thus, ordering the samples by decreasing prior density should achieve a similar effect as ordering the posteriors by decreasing overlap.

Our complete algorithm, LMIS with mixture pruning and an initial sorting of samples, is described in Algorithm <ref>.




[htbp]
    

    funcFunction
    EIGExpectedInformationGain
    MISMomentsEstimatePosteriorMoments
    funcFunction
    textsl

    
    i  1 … N
        Sample (^(i),^(i)) from  the prior p(,) 
        Sample ^(i) from the likelihood p(|^(i),^(i),) 
    
    Rearrange {(^(i),^(i))}_i=1^N and ^(i) in decreasing order of p(^(i),^(i))  
    i  1 … N
    𝒥^(i){m:q_marg^(m)(^(i),^(i)) > p(^(i),^(i)), m<i} 
    _post^(i), _post^(i)^(i), 𝒥^(i) 
    Estimate the marginal likelihood p^(i)|
    q_marg^(i) t_ν(_post^(i),_post^(i)) 
    j 1… M_1
    Sample (^(i,j), ^(i,j)) from q_marg^(i)(^(i,j), ^(i,j)) 
    
    p(^(i)|)1/M_1∑_j=1^M_1p(^(i)|^(i,j),^(i,j),) p(^(i,j),^(i,j))/q_marg^(i)(^(i,j),^(i,j)) 
    Estimate the conditional likelihood p^(i)|^(i),
    Condition _post^(i) and _post^(i) on ^(i) to obtain _cond^(i) and _cond^(i) 
    q_cond^(i) t_ν(_cond^(i),_cond^(i)) 
    k  1 … M_2 
        Sample ^(i,k) from q_cond^(i)() 
    
    p(^(i)|^(i),)1/M_2∑_k=1^M_2p(^(i)|^(i),^(i,k),) p(^(i,k))/q_cond^(i)(^(i,k)) 
    
    U()1/N∑_i=1^Nlog[p(^(i)|^(i),)]-log[p(^(i)|)]
    
    Layered multiple importance sampling estimate of the expected information gain
    



[htbp]
    1.25
    textsl
    funcFunction
    EIGExpectedInformationGain
    MISMomentsEstimatePosteriorMoments
    funcFunction
    textsl

    ^(i), 𝒥^(i)


    Create mixture biasing distribution q_post^(i) according to the index set J^(i)
    Initialize {(^(i,ℓ),^(i,ℓ))}_ℓ=1^N{(^(i),^(i))}_i=1^N drawn earlier from the prior  
    L N 
    m∈𝒥^(i)
    {(^(i,ℓ),^(i,ℓ))}_ℓ=1^L+M_1{(^(i,ℓ),^(i,ℓ))}_ℓ=1^L⊔{(^(m,j), ^(m,j))}_j=1^M_1 
     
    L L+M_1
    

    q_post^(i)(,)N/L p(,) + M_1/L∑_m∈𝒥^(i)q_marg^(m)(,) 
    Compute importance weights and normalize them
    (^(i,ℓ),^(i,ℓ))∈{(^(i,ℓ),^(i,ℓ))}_ℓ=1^L
    w_post^(i,ℓ)p(^(i)|^(i,ℓ),^(i,ℓ),) p(^(i,ℓ),^(i,ℓ))q_post^(i)(^(i,ℓ),^(i,ℓ))
    
    w_post^(i,ℓ)w_post^(i,ℓ) / ∑_ℓ=1^L w_post^(i,ℓ) 
    Estimate posterior mean and covariance
    _post^(i)∑_ℓ=1^Lw_post^(i,ℓ)[[ ^(i,ℓ); ^(i,ℓ) ]] 
    _post^(i)∑_ℓ=1^Lw_post^(i,ℓ)([[ ^(i,ℓ); ^(i,ℓ) ]]-_post^(i))([[ ^(i,ℓ); ^(i,ℓ) ]]-_post^(i))^ 
    
    Multiple importance sampling estimate of posterior moments
    




§ NUMERICAL RESULTS


In this section, we examine two experimental design problems. The first is a linear Gaussian problem with a known, exact solution; this example is used to assess the performance and convergence properties of the layered multiple importance sampling algorithm, for varying dimensions. The second example is a more challenging nonlinear experimental design problem that occurs in Mössbauer spectroscopy, where we again illustrate the computational efficieny of the LMIS approach as well as the impact of the focused design objective.



 §.§ Higher-dimensional linear Gaussian example


We revisit the linear Gaussian problem that was first introduced in <ref>, but increase the dimension of the inference problem to n=4 and n=8 dimensions. The dimension of the parameter of interest remains the same, i.e., n_θ=1, so the dimension of the nuisance parameters is n_η = n-n_θ. 
In general, the observation model for our n-dimensional inference problem is:

    [ y_1; y_2;   ⋮; y_n ]=
        [     kd                    1;        k(1-d)              ;                    ⋱       ;      1               k(1-d) ]_(d,k)[     θ;   η_1;     ⋮; η_n-1 ] +
        [ ϵ_1; ϵ_2;   ⋮; ϵ_n ],    k>0.

The design matrix  is parameterized by the parameter d∈[0,1] and a gain parameter k>0.
The diagonal terms describe a trade-off between learning about the parameter of interest and learning about the nuisance parameters where d=1 maximizes the information gain in θ at the expense of information gain in the nuisance parameters (η_1,η_2,…,η_n-1). Furthermore, the nonzero off-diagonal terms in the design matrix introduce coupling between the parameter of interest and the nuisance parameters so that the resulting optimal design is not trivially d=1.

Each of the input parameters have standard normal priors, i.e., (θ,η_1,…,η_n-1)^T∼𝒩(0, ^pr), with ^pr=𝐈_n × n.
The observation noise is additive Gaussian with zero mean and variance σ_ϵ^2. i.e., ∼𝒩(0,^obs) where ^obs=σ_ϵ^2𝐈_n × n.

The EIG for all of the parameters is

    U_joint(d)=-1/2log[(^post)]

where ^post is the posterior covariance matrix, given by

    ^post=^pr-^pr^T(^pr^T+^obs)^-1^pr,

and the EIG for the marginal in the parameter of interest is

    U_marg(d)=-1/2log[γ^post_1,1]

where γ^post_1,1 is the  top left entry of the posterior covariance matrix ^post.



  §.§.§ Expected information gain profiles

For the following numerical examples in n=4 and n=8 dimensions, we set k=5 and σ_ϵ=0.4.
The EIG for the focused and unfocused settings are plotted in <ref>.
We see that the optimal designs for the two objectives are very different: the optimal design for learning about all the parameters is d_joint^∗=0 and the optimal design for learning only about the parameter of interest is d_marg^∗≈ 0.93.
The determinant of the posterior covariance matrix is plotted in <ref> as a function of the design parameter d, which illustrates the relationship between the posterior covariance and the EIG.



From the analysis of the EIG estimator in <ref> and <ref>, we expect that the bias and variance will be largest when the biasing distributions and the posterior distributions have the smallest overlap.
Thus, if we use the the prior distribution as the biasing distribution, the bias and variance of the EIG estimator should be large when there is significant posterior concentration.
Accordingly, in <ref> we observe that the volume (computed as the determinant of the posterior covariance matrix) of the posterior distribution varies over several orders of magnitude across the design space, and in <ref> and <ref> we observe large positive bias for designs that correspond to significant posterior concentration.

The correlation of estimator bias with posterior concentration is particlularly undesirable behavior—particularly in the context of focused experimental design, since bias that is correlated with the determinant of the joint posterior can mask the true shape of the objective function (which depends only on the marginal posterior).
We observe this behavior in <ref> where the location of local and global maxima correspond to sub-optimal designs that maximize posterior concentration instead of those that maximize information gain in the parameter of interest.
Furthermore, when using the prior as a biasing distribution, the large estimator variance may pose a challenge for stochastic optimization algorithms. On the other hand, LMIS yields estimates of the EIG with much smaller bias and variance while using the same number of model evaluations.
More importantly, the LMIS estimates of the EIG in both the focused and unfocused cases capture the correct locations of the global maxima.





  §.§.§ Sample allocations, bias, and variance

The allocation of samples between the outer and inner estimators can play an important role in the bias and variance of the overall estimator.
In <ref>, we investigate the effect of the ratio between the number of outer and inner samples on the mean squared error of the EIG estimator, both using the prior as a biasing distribution and using layered multiple importance sampling.
When using the prior as a biasing distribution, the bias dominates the mean squared error, so allocating more inner samples—at the expense of affording fewer outer samples—will reduce the overall mean squared error.
On the other hand, with LMIS, the biasing distributions closely approximate the posterior distributions, so fewer inner samples are required, and more samples can be allocated for the outer estimator, which not only decreases the overall variance, but also provides additional opportunities for incremental enrichment, further improving the approximation accuracy of the inner estimators.




Using favorable but fixed ratios between the number of outer and inner samples,[When using prior-based sampling, we set M=100N and when using LMIS, we set N=100M.] we can examine the mean squared error as a function of number of model evaluations.
In <ref>, we observe a dramatic reduction in the mean squared error when using LMIS. 
With 10^6 model evaluations, the resulting LMIS estimate has mean squared error comparable to using the exact posteriors as biasing distributions.[When using the exact joint and conditional posteriors as biasing distributions, the importance estimators of the marginal and conditional likelihoods are zero-variance estimators. In theory, a single inner sample is then sufficient to obtain an exact value of the marginal and conditional likelihood. However, for a fair comparison, we use the same number of inner and outer samples as LMIS, even if it is not the most efficient allocation when using the exact posterior.]

To better understand the observed convergence rates in <ref>, we plot the mean squared error versus the number of inner and outer samples in <ref>.
When using the prior as the biasing distribution, we observe the convergence rates predicted by <ref> and <ref>. The estimator bias decreases as 1/M and is unaffected by the number of outer samples, N. The variance also decreases with increasing number of inner and outer samples, which suggests that the 1/NM terms in <ref> are not negligible.

On the other hand, LMIS allows the estimator bias to decrease with the number of outer samples. 
With larger N, more samples can initially be used to estimate the posterior moments used to find optimal biasing distributions. 
Furthermore, the LMIS estimator applies incremental enrichment a total of N times, so as N becomes larger, more samples are made available for use when estimating the posterior moments. 
When using LMIS, the variances of the inner estimators are significantly reduced, so the 1/NM terms in the variance for the LMIS estimator are negligible in comparison to prior-based sampling. 
Thus, instead of the usual bias–variance trade-off between maximizing the number of inner and outer samples, the earlier analysis and the numerical results suggest that maximizing the number of outer samples is the best way to reduce both the bias and variance!
In fact, if the biasing distributions converge to the posterior distributions as N→∞ (e.g., if they are of the same family of distributions), it is conceivable that one could set M=1 in this limit, since importance samping using the normalized posterior results in a zero-variance estimator. However, for finite sample sizes and non-Gaussian posterior distributions, we still require M→∞ to ensure asymptotic consistency, since the inner estimators correct for the mismatch between the biasing distributions and the posterior distributions. 




  §.§.§ Importance sampling diagnostics and incremental enrichment
 
The exact value of the EIG is generally not available, so in practice, we must instead rely on importance sampling diagnostics to measure the performance of the estimator.
A useful “rule of thumb” for measuring importance sampling performance is to use the effective sample size (ESS) to measure how different the biasing distribution is from the target distribution <cit.>.
However, in our case, the target distribution is the prior distribution, so maximizing the ESS diagnostic would suggest choosing the prior as the biasing distribution, which we showed to be sub-optimal.
Instead, we consider the effective sample size customized for the target function—the marginal or conditional likelihood—as suggested in <cit.>. For a diagnostic that is specific to the importance sampling estimator in <ref>, we define

    w_ML^(i,j)   p(^(i)|^(i,j),^(i,j),) p(^(i,j),^(i,j))/q_ML^(i)(^(i,j),^(i,j))/∑_j=1^M_1p(|^(i,j),^(i,j),)p(^(i,j),^(i,j))/q_ML^(i)(^(i,j),^(i,j)),
    w_CL^(i,j)   p(^(i)|^(i),^(i,j),) p(^(i,j)^(i))/q_CL^(i)(^(i,j))/∑_j=1^M_1p(|^(i),^(i,j),)p(^(i,j)^(i))/q_CL^(i)(^(i,j)).

Then, the effective sample sizes customized to the marginal likelihood and the conditional likelihood are

    cESS(w_ML)=1/∑_j=1^M_1(w_ML^(i,j))^2,  cESS(w_CL)=1/∑_j=1^M_1(w_CL^(i,j))^2.

In <ref>, we can see that LMIS greatly increases the effective sample size customized to the marginal and conditional likelihood for both small N=320 and large N=1000 values of N and M.





 In <ref>, we explore the mechanism by which incremental enrichment can improve the quality of the biasing distributions for the inner estimators. Since the samples of (^(i),^(i)) and the corresponding ^(i) are sorted in decreasing order of prior density p(^(i),^(i)), we should expect that the effective sample size of the inner estimators decrease with increasing outer loop index i. This is particularly evident when using the prior as a biasing distribution, since as i increases, the ith posterior p(,^(i),) will be further out in the tail of the prior distribution p(,), and is reflected in the monotonic decrease of cESS with i seen in <ref>. On the other hand, when using LMIS, the cESS initially increases for small i, as a direct result of incremental enrichment. Then, the cESS decays more slowly for the bulk of the inner estimates than in the prior case. Finally, for large i, the quality of the biasing distributions deteriorates rapidly, since there are fewer available samples for reuse when estimating the moments of posterior distributions concentrated in the tail of the prior. However, as N becomes larger, the impact of the deterioration in the tails diminishes (cf. the decrease in probability mass in the lower-right corner of each subplot in <ref>). 
The deterioration in the tails is also reflected in <ref> where the bias of the LMIS does not decrease as 1/M for fixed N. This can be attributed to the observation that the biasing distributions q_CL^(i) for estimating the conditional likelihoods p(^(i),) are obtained by calculating conditionals of an estimated posterior covariance. The conditional mean and variance are very sensitive to the correlation in the estimated posterior covariance, which may not be a robust estimator in the small sample size limit. The result is that the variance of the inner estimator of the conditional likelihood may decrease more slowly than 1/M. 






 §.§ Nonlinear example: Mössbauer spectroscopy


The Mössbauer effect refers to recoil-free nuclear resonance fluorescence, which involves the resonant and recoil-free emission and absorption of gamma radiation by atomic nuclei bound in a solid <cit.>. Mössbauer spectroscopy is responsible for the discovery of the Mössbauer isomeric shift <cit.>, which carries 
valuable information regarding the electronic and chemical structure of the lattice in which emitting and absorbing nuclei are embedded. These shifts are incredibly small, usually on the order of 10^-8 eV, which is roughly 12 orders of magnitude smaller than the energy of the emitted gamma rays themselves. Spectroscopic experiments involve moving the gamma ray source relative to the absorber. Change in the energy of the emitted gamma ray from Doppler shifting can be precisely controlled by the relative velocity of the source and absorber. In practice, shifts of 10^-8 eV correspond to velocities on the order of mm/s.

The general experimental setup for Mössbauer spectroscopy involves a radioisotope source, a drive system to move the source (either with constant acceleration or constant velocity), an absorber, and a gamma ray detector.
The Mössbauer drive is often operated at constant velocity for precision measurements, so that integration time may be spent only at the velocities of interest, which, in general, are not uniformly or symmetrically distributed about zero <cit.>. If one uses long counting times and small velocity steps to obtain high precision, the total time of the run can take days, so careful selection of velocities at which to acquire data is critical for a successful experiment.



  §.§.§ Problem formulation

An experimental design consists of a set of n_d velocities for high precision measurement of the isomeric shift—the center of the absorption peak—which we denote by θ. 
We use the standard parameterization of the absorption peak as having a Lorentzian profile <cit.>, such that the the number of detector counts y_i at a velocity d_i can be described by the following relationship, also illustrated in <ref>:

    y_i(d_i; · ) = offset - heightwidth^2/width^2 + (center-d_i)^2 +ϵ_i,   i=1,2,…,n_d

The model parameters are assigned the following prior distributions,

    center∼𝒩(0,1),      exp(width)∼𝒩(0,0.3^2),
    exp(height)∼𝒩(0,0.3^2),      exp(offset)∼𝒩(1.0,0.2^2),

where log-normal priors are assigned to parameters which can only have positive values, and the mean and variance of the prior distributions are chosen to “standardize” the signal, which would be standard experimental practice before collecting data. The additive observation error is assumed to be Gaussian  where ϵ_i∼𝒩(0,0.1^2), which represents a signal-to-noise ratio of 10:1. 


In <ref>, we plot two “maps” of EIG 
estimated using LMIS, using sample sizes M_1 = M_2 = 100 and N= 1000. 
One map uses θ = `center' and lets η comprise the remaining parameters; the other map puts θ = `offset' and lets η comprise the rest. The optimal designs differ for these two scenarios: d_center^∗ = (-1.3,0,1.3) and d_offset^∗ = (-2,0,2). For visualization purposes, we restrict the design problem to three velocities, with one of the three velocities fixed to zero. To better understand the impact of these design choices, <Ref> plots the posterior distribution obtained using a data realization from each optimal design. Comparing the two focused designs, we observe a distinct tradeoff between learning about the `center' parameter and the `offset' parameter, especially in the lower-left pairwise marginal. Furthermore, the two-dimensional marginals clearly illustrate the non-Gaussianity of the posterior distribution.




Next we evaluate the performance of the LMIS estimator, focusing on the estimation of EIG in the marginal distribution of θ = `center.' Since we do not have a closed-form expression for this EIG, we resort to calculating a “reference” value by using a large computational budget: 10^8 model evaluations with the LMIS estimator. In <ref>, we plot the mean squared error with respect to this reference value for both the LMIS estimator and the prior estimator. We choose to evaluate the focused EIG for the design  _center^∗=(-1.3, 0, 1.3) since it presents the most challenging case—with the greatest amount of posterior concentration. Compared to the prior estimator, we observe a reduction of up to two orders of magnitude in the mean squared error when using LMIS. This result is consistent with the earlier linear-Gaussian example and suggests that the algorithm is able to identify high-quality biasing distributions even in the nonlinear case.


We can also use the cESS diagnostic <ref> to evaluate the performance of the LMIS scheme. The distribution of cESS is plotted in <ref>, and displays a marked increase over that obtained when using the prior as a biasing distribution, comparable to the improvement for the linear Gaussian example of the same dimension. 
Again, this result suggests that LMIS is able to find high-quality biasing distributions for estimating both the marginal and conditional likelihoods, even in the case of a nonlinear observational model and non-Gaussian posteriors. 









§ CONCLUSION


This paper introduces a new computational approach for focused optimal Bayesian experimental design, where the design objective is to maximize expected information gain in the marginal distribution of parameters of interest. In particular, we develop a consistent layered multiple importance sampling (LMIS) estimator of this objective, which uses incremental enrichment to construct a sequence of increasingly tailored biasing distributions.
We demonstrate the effectiveness of our approach on design problems of moderate dimension. Compared to previous algorithms that draw samples from the prior, the LMIS estimator reduces both bias and variance by several orders of magnitude, for the same number of model evaluations. The reduction in variance may improve the efficiency of optimization over the design space, but more importantly, the reduction in bias diminishes the confounding effect caused by link between estimator bias and posterior concentration. The latter is particularly important in the context of focused design, as posterior concentration in general has a different dependence on the design parameters than the information gain in selected marginals of interest.

The enrichment process of LMIS suggests a different balance of inner-loop and outer-loop samples than previous nested Monte Carlo methods. We also emphasize that LMIS does not require gradients of the likelihood or deterministic models embedded therein, and thus is well suited to problems where gradient information may not be available. While here we implemented LMIS with a particular family of inner-loop biasing distributions (multivariate t), the overall algorithm can be used with any parametric family or mixture, as suited to the form of the posterior.


Possible extensions of this work include using surrogate models to accelerate estimation; while previous work did so in a biased way <cit.>, the LMIS framework could incorporate surrogate models only in its outer layer, to help construct biasing distributions more cheaply. When gradients of the posterior density are readily available, one also could use Laplace approximations of the posterior distribution to construct suitable biasing distributions for estimation of the marginal <cit.> and conditional likelihoods. Moreover, surrogate models and approximations of the marginal likelihood could be combined with the full-model LMIS estimator using a control variate or multilevel Monte Carlo approach, to further reduce the number of full model evaluations while preserving consistency. Another clear extension is to pair LMIS estimators with optimization methods designed for stochastic objectives (e.g., <cit.>) to facilitate search over the design space. 

An interesting application of focused experimental design is for models that contain discrepancy terms <cit.> intended to capture model misspecification. In this setting, one could treat parameters of the discrepancy model as nuisance parameters (η), where the goal is to learn about other parameters of interest (θ) in the presence of model discrepancy. Alternatively, one could focus on learning about the discrepancy term to better capture the nature of the mismatch between the model and data.





§ ASYMPTOTIC BIAS AND VARIANCE OF THE EIG ESTIMATOR

In this section, we determine the asymptotic behavior of the bias and
variance of the expected information gain (EIG) estimator U() (<ref>) with respect to the number of samples in both the inner and outer sums.
Let A:ℝ^n_η→ℝ^+ and B:ℝ^n_η×ℝ^n_θ→ℝ^+ be scalar functions defined as

    A(;,,)    p(,,),
    
      B(,;,)    p(,,),

where underlined variables are fixed parameters. 
The expectations of these functions are equal to the conditional and
marginal likelihoods, respectively:

    [    μ_A(,;)            _|[A(;,,)]          =      p(,),;     μ_B(;)            _,[B(,;,)]          =       p(). ]

Let ={^(k)}_k=1^M_2 and ={(^(j),^(j))}_j=1^M_1 be samples drawn i.i.d. from the biasing distributions q_cond^(,)() and q_marg^()(,), respectively.
Then, the importance sampling estimates of the conditional and marginal likelihoods from <ref> can be written as

    A(,,;)    = 1/M_2∑_j=1^M_2 A(^(j);,,)p(^(k)^(k))/q_cond^(,)(^(k)),
    B(,;)    = 1/M_1∑_j=1^M_1 B(^(j),^(j);,)p(^(j),^(j))/q_marg^()(^(j),^(j)),

Now, define

    u(,,,;) lnA(,,;)/B(,;)

so that U() from (<ref>) can be written as

    U() = 1/N∑_i=1^N u(^(i),^(i),S_A^(i),S_B^(i);)

where ^(i)∼ p(^(i)^(i),^(i),), (^(i),^(i))∼ p(,),
and S_A^(i) and S_B^(i) are sets of i.i.d. samples drawn from q_cond^(^(i),^(i)) and q_marg^(^(i)), respectively.

From the law of iterated expectations, the expected value of U() can be written as

    [U()]   =_,[_,[u(,,,;)]],
      and from the law of total variance, the variance of U() can be written as[U()]   =1/N_,,,[u(,,,)]
       =1/N_,[_,[u(,,,;)]] + 1/N_,[_,[u(,,,;)]].

Using the delta method, we can obtain the second-order approximation of the expected value and variance of u(,,,;):

    _,[u(,,,;)]
       =lnμ_A(,)/μ_B()+1/2[1/M_1σ_B^2(;)/μ_B^2(;)-1/M_2σ_A^2(,;)/μ_A^2(,;)]
           +𝒪(1/M_1^2) + 𝒪(1/M_2^2),
    _,[u(,,,;)]
       =1/M_2σ_A^4(,;)/μ_A^4(;)+1/M_1σ_B^4(,;)/μ_B^4(;)
           +𝒪(1/M_1^2) + 𝒪(1/M_2^2),

where we have used the fact that A(,;) and B(;) are unbiased estimators, so their expected values are

    _,[A(,,;)]    = μ_A(,;), 
    _,[B(,;)]   =μ_B(;),

and their variances are

    _,[A(,,;)]    =
    1/M_2_∼ q_cond^(, )[A(;,,)p()/q_cond^(, )()]
    =σ_A^2(,;)/M_2,
    _,[B(,;)]    =
    1/M_1_(,)∼ q_marg^()[B(,;,)p(,)/q_marg^()(,)]
    =σ_B^2(;)/M_1.

Then, the expected value of U() from (<ref>) becomes

    [U()]    =
      _,[
        lnμ_A(,;)/μ_B(;)]_U()
       + 1/M_1_[
        σ_B^2(;)/2μ_B^2(;)]_C_1()
       - 1/M_2_,[
        σ_A^2(,;)/2μ_A^2(,;)]_C_2()
           +𝒪(1/M_1^2) + 𝒪(1/M_2^2)
       = U() + C_1()/M_1-C_2()/M_2+𝒪(1/M_1^2) + 𝒪(1/M_2^2),

And thus the bias of U() is

    [U()-U()] = C_1()/M_1 - C_2()/M_2+𝒪(1/M_1^2) + 𝒪(1/M_2^2).

Similarly, the variance of U() from (<ref>) becomes

    [U()]    ≈1/N(1/M_2_,[σ_A^4(,;)/μ_A^4(,;)]^D_1()+1/M_1_[σ_B^4(;)/μ_B^4(;)]^D_2()+𝒪(1/M_1^2) + 𝒪(1/M_2^2))+
           1/N(_,[lnμ_A(,;)/μ_B(;)]_D_3() +1/(M_1)^2_[σ_B^2(;)/2μ_B^2(;)]_D_4()+1/(M_2)^2_,[σ_A^2(,;)/2μ_A^2(,;)]_D_5()
                    +𝒪(1/M_1^4) + 𝒪(1/M_2^4))
       =D_3()/N + D_1()/NM_2  + D_2()/NM_1  + D_4()/N(M_2)^2  + D_5()/N(M_1)^2 + 𝒪(1NM_1^3)+𝒪(1/NM_2^3).







§ SAMPLE SIZES AND SCALING

The nested estimator (<ref>) is consistent (see <cit.>) and in particular, its bias and variance converge to zero as M_1, M_2,  N →∞. Since all three sample sizes must increase, it is natural to ask whether there is an optimal allocation of inner (M_1, M_2) and outer (N) samples for a given computational budget W = N(M_1 + M_2).

For simplicity, let M_1 = M_2 = M and define  C (C_1- C_2)^2>0 and D D_1+D_2, where we drop dependence on  to reduce clutter. Then W = 2MN and the bias and variance from (<ref>) and (<ref>) can be written as

    bias   ∼C^1/2/M,    variance∼D_3/N+D/NM,

yielding an expression for the mean squared error (MSE) of the estimator, to leading order, in terms of M and N:

    MSE    = 𝔼 [ |U - U |^2  ] = bias^2 + variance∼C/M^2+D_3/N+D/MN.

Let α^2=2M/N be the ratio between the number of inner and outer Monte Carlo samples. Then, by substituting M=Nα^2/2 and N=√(W)/α, 
we can rewrite the MSE in terms of the total computational work W and α:

    MSE(W,α)    ∼4C/α^2 W+α D_3/W^1/2+2D/W.

The value of α and the corresponding ratio α^2 that minimize the leading-order MSE are

    α_∗=2  ( C/D_3 )^1/31/W^1/6  ⇒  α_∗^2= 4  ( C/D_3 )^2/31/W^1/3.

This result suggests that the ratio between the number of inner and outer samples should decrease slowly as the computational budget increases, in particular scaling as W^-1/3. The ratio C/D_3 reflects the relative magnitudes of leading-order terms in the bias and variance; if the bias dominates, one should choose a relatively larger number of inner loop samples. The minimum MSE, realized by this value of α_∗, is

    MSE(W,α_∗) ∼2D/W+3C^1/3D_3^2/3/W^2/3.

Of course, if one does not know the constant in (<ref>) and instead just uses α = α_0/W^1/6, the work-scaling of the MSE remains the same: 

    MSE(W, α_0/W^1/6 ) ∼2D/W +  ( 4 C/α_0^2 + α_0 D_3  ) 1/W^2/3

Thus we can compare two scenarios. For a fixed ratio between M and N, i.e., if  α^2 were not allowed to depend on W, the leading-order term of (<ref>) would yield MSE = 𝒪(W^-1/2). If, instead, we adjust the ratio according to the optimal scaling α^2 ∝ W^-1/3, we obtain a faster decay rate of MSE = 𝒪(W^-2/3).



We use the 4-dimensional linear Gaussian problem to verify the results of the optimal scaling analysis in <ref>.
In <ref>, we plot the estimator MSE against the number of model evaluations using two scaling strategies: a fixed ratio (dashed line) and a scaled ratio (solid line).
In practice, the constant term (C/D_3)^2/3 in <ref> for the optimal ratio is not known.
In <ref>, we anchor the ratio at W=10^4 to be M=10N (left) and M=1000N (right) to investigate the influence of the constant term on the convergence rate for finite sample sizes.
We observe that when the MSE is dominated by the estimator bias (e.g., M=10N), the trade-off explored by the optimal scaling actually results in sub-optimal convergence in MSE.
This is because the benefits of a faster convergence rate in the estimator variance are overshadowed by the slower convergence rate in the bias.
On the other hand, when the bias and variance are of similar magnitude (e.g., M=1000N at W=10^4), we can observe benefits of a faster convergence rate.
Between W=10^6 and W=10^7, the MSE converges at the predicted rates: W^-1/2 when using a fixed ratio and W^-2/3 when using the optimal scaling.
In practice, without knowing the optimal value of the constant from <ref> or the magnitude of the estimator bias, a strategy that uses a fixed ratio between the number of samples in the inner and outer estimators may be appropriate for realistic computational budgets.  



siamplain

