Optimal Experimental Design Using A Consistent Bayesian Approach

TL;DR

This paper introduces a novel optimal experimental design method based on a consistent Bayesian framework, which ensures the posterior distribution aligns with observed data and model predictions, improving data acquisition strategies.

Contribution

It develops a new OED approach utilizing the consistent Bayesian method, ensuring posterior-data consistency and computational efficiency in PDE-based models.

Findings

Outperforms classical Bayesian OED methods in numerical tests.

Ensures posterior consistency with observed data.

Effective for PDE-based model experiments.

Abstract

We consider the utilization of a computational model to guide the optimal acquisition of experimental data to inform the stochastic description of model input parameters. Our formulation is based on the recently developed consistent Bayesian approach for solving stochastic inverse problems which seeks a posterior probability density that is consistent with the model and the data in the sense that the push-forward of the posterior (through the computational model) matches the observed density on the observations almost everywhere. Given a set a potential observations, our optimal experimental design (OED) seeks the observation, or set of observations, that maximizes the expected information gain from the prior probability density on the model parameters. We discuss the characterization of the space of observed densities and a computationally efficient approach for rescaling observed densities to satisfy the fundamental assumptions of the consistent Bayesian approach. Numerical results are presented to compare our approach with existing OED methodologies using the classical/statistical Bayesian approach and to demonstrate our OED on a set of representative PDE-based models.

Figures

Tables

Table 1. Table 1: The top 5 experimental designs chosen using the full set of 5,000 samples. For each of these designs, we compute E ​ ( I Q z ) 𝐸 subscript 𝐼 superscript 𝑄 𝑧 E(I_{Q^{z}}) for 50, 200, 1,000 and 5,000 samples. Notice the change in E ​ ( I Q z ) 𝐸 subscript 𝐼 superscript 𝑄 𝑧 E(I_{Q^{z}}) for a given design decreases as we increase to 5,000 samples.

Design Location	50	200	1,000	5,000
$(0.558,0.571)$	2.758	2.767	2.815	2.826
$(0.561,0.546)$	2.752	2.762	2.809	2.820
$(0.582,0.574)$	2.729	2.736	2.782	2.793
$(0.549,0.570)$	2.728	2.735	2.781	2.792
$(0.593,0.596)$	2.726	2.733	2.779	2.790

Table 2. Table 2: The top 5 experimental designs chosen using the full set of 5,000 samples. For each of these designs, we compute E ​ ( I Q z ) 𝐸 subscript 𝐼 superscript 𝑄 𝑧 E(I_{Q^{z}}) for 50, 200, 1,000 and 5,000 samples. Notice the change in E ​ ( I Q z ) 𝐸 subscript 𝐼 superscript 𝑄 𝑧 E(I_{Q^{z}}) for a given design decreases as we increase to 5,000 samples.

Design Location	50	200	1,000	5,000
$(0,0)$	0.687	0.828	0.738	0.741
$(1,1)$	0.653	0.713	0.747	0.740
$(0,0.1)$	0.687	0.817	0.733	0.736
$(0.1,0)$	0.687	0.810	0.728	0.735
$(1,0.9)$	0.648	0.713	0.742	0.735

Table 3. Table 3: The top 5 experimental designs chosen using the full set of 1,000 samples. For each of these designs, we compute E ​ ( I Q z ) 𝐸 subscript 𝐼 superscript 𝑄 𝑧 E(I_{Q^{z}}) for 10, 50, 100 and 1,000 samples.

Design Location	10	50	100	1,000
$(-1.001,0.961)$	2.303	3.662	4.085	4.569
$(-1.050,1.035)$	2.303	3.391	3.811	4.261
$(1.041,0.959)$	2.303	3.437	3.839	4.256
$(-1.050,1.130)$	2.302	3.418	3.661	4.096
$(1.005,0.870)$	2.277	3.246	3.544	3.813

Table 4. Table 4: The top 5 experimental designs chosen using the full set of 10,000 samples. For each of these designs, we compute E ​ ( I Q z ) 𝐸 subscript 𝐼 superscript 𝑄 𝑧 E(I_{Q^{z}}) for 50, 100, 1,000 and 10,000 samples. Notice the change in E ​ ( I Q z ) 𝐸 subscript 𝐼 superscript 𝑄 𝑧 E(I_{Q^{z}}) for a given design decreases as we increase to 10,000 samples.

Design Location	50	100	1,000	10,000
$(0.48,1)$	1.966	2.001	1.992	2.008
$(0.5,0.98)$	1.952	1.963	1.993	2.007
$(0.46,0.98)$	1.930	1.985	1.986	2.006
$(0.52,0)$	1.777	1.915	1.999	2.006
$(0.56,0)$	1.751	1.863	1.996	2.006

Equations

$$Q^{-1}(A)=\left\{\lambda \in \Lambda |Q(\lambda )\in A\right\}\in {\mathcal{B}}_{\Lambda },\text{and}Q(Q^{-1}(A))=A.$$

$$P_{\mathcal{D}}^{Q(P_{\Lambda })}(A)=P_{\Lambda }(Q^{-1}(A)).$$

$$P_{\Lambda }(Q^{-1}(A))=P_{\mathcal{D}}^{Q(P_{\Lambda })}(A)=P_{\mathcal{D}}^{\text{obs}}(A),$$

$$P_{\mathcal{D}}^{Q(\text{prior})}(A)=P_{\Lambda }^{\text{prior}}(Q^{-1}(A)).$$

$$P_{\Lambda }^{\text{post}}(B):=\{\begin{array}{ll}{P}_{\Lambda }^{\text{prior}}(B)\frac{P_{\mathcal{D}}^{\text{obs}}(Q(B))}{P_{\mathcal{D}}^{Q(\text{prior})}(Q(B))},& \text{if }{P}_{\Lambda }^{\text{prior}}(B)>0,\\ 0,& \text{otherwise},\end{array}$$

$${\pi }_{\Lambda }^{\text{post}}(\lambda )={\pi }_{\Lambda }^{\text{prior}}(\lambda )\frac{{\pi }_{\mathcal{D}}^{\text{obs}}(Q(\lambda ))}{{\pi }_{\mathcal{D}}^{Q(\text{prior})}(Q(\lambda ))},\lambda \in \Lambda .$$

$$I_{Q}(\pi_{{\Lambda}}^{\text{prior}}:\pi_{{\Lambda}}^{\text{post}}):=\int_{{\Lambda}}\pi_{{\Lambda}}^{\text{post}}\log\Bigg{(}{\frac{\pi_{{\Lambda}}^{\text{post}}}{\pi_{{\Lambda}}^{\text{prior}}}}\Bigg{)}d\mu_{{\Lambda}}.$$

$$I_Q({\pi }_{\Lambda }^{\text{prior}}:{\pi }_{\Lambda }^{\text{post}})=I_Q({\pi }_{\Lambda }^{\text{post}}),$$

$$I_Q({\pi }_{\Lambda }^{\text{post}})=I_Q({\pi }_{\mathcal{D}}^{\text{obs}}).$$

$$I_Q({\pi }_{\Lambda }^{\text{post}})$$

$${\lambda }_1{x}_1^2+x_2^2=1$$

$$x_1^2-{\lambda }_2{x}_2^2=1$$

$$Q_2(\lambda )=x_1.$$

$$\mathcal{O}_{\mathbf{\mathcal{D}}}=\Bigg{\{}\hat{N}(q,\sigma^{2}):q\in\mathbf{\mathcal{D}}\Bigg{\}},$$

$$E(I_Q):={\int }_{\mathcal{D}}{I}_Q(q){\pi }_{\mathcal{D}}^{Q(\text{prior})}(q)d\mu ={\int }_{\mathcal{D}}{I}_Q(q)d{P}_{\mathcal{D}}^{Q(\text{prior})}.$$

$$E(I_{Q})=\int_{\mathbf{\mathcal{D}}}\int_{{\Lambda}}\pi_{{\Lambda}}^{\text{post}}(\lambda;q)\log\Bigg{(}{\frac{\pi_{{\Lambda}}^{\text{post}}(\lambda;q)}{\pi_{{\Lambda}}^{\text{prior}}(\lambda)}}\Bigg{)}d\mu_{{\Lambda}}dP_{\mathbf{\mathcal{D}}}^{Q(\text{prior})},$$

$$I_{Q}(q^{(j)})\approx\frac{\mu_{{\Lambda}}({\Lambda})}{N}\sum_{i=1}^{N}\frac{\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}(Q(\lambda^{(i)}))}{\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}(Q(\lambda^{(i)}))}\log\Bigg{(}\frac{\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}(Q(\lambda^{(i)}))}{\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}(Q(\lambda^{(i)}))}\Bigg{)}$$

$$Q^{\text{opt}}:=\arg \underset{Q^z\in \mathcal{Q}}{\max }E(I_{Q^z}).$$

$$\mathcal{O}_{\mathbf{\mathcal{D}}}=\Bigg{\{}\frac{\hat{N}(q,\sigma^{2})}{C_{q}}:q\in\mathbf{\mathcal{D}}\Bigg{\}},$$

$$C_q={\int }_{\mathcal{D}}\hat{N}(q,{\sigma }^2)d\mu .$$

$$P_{\Lambda }^{\text{post}}(\Lambda )=P_{\Lambda }^{\text{prior}}(\Lambda )\frac{P_{\mathcal{D}}^{\text{obs}}(Q(\Lambda ))}{P_{\mathcal{D}}^{Q(\text{prior})}(Q(\Lambda ))},$$

$$P_{\Lambda }^{\text{post}}(\Lambda )=P_{\mathcal{D}}^{\text{obs}}(Q(\Lambda )).$$

$${\pi }_{\Lambda }^{\text{post}}=\frac{{\tilde{\pi }}_{\Lambda }^{\text{post}}}{{\int }_{\Lambda }{\tilde{\pi }}_{\Lambda }^{\text{post}}d{\mu }_{\Lambda }}.$$

$$P_{\Lambda }^{\text{post}}(\Lambda )={\int }_{\Lambda }{\pi }_{\Lambda }^{\text{post}}d{\mu }_{\Lambda }={\int }_{\Lambda }\frac{{\pi }_{\mathcal{D}}^{\text{obs}}(Q(\lambda ))}{{\pi }_{\mathcal{D}}^{Q(\text{prior})}(Q(\lambda ))}d{P}_{\Lambda }^{\text{prior}}.$$

$$\{\begin{array}{ll}-D{\nabla }^{2}u+\nabla \cdot (vu)=S,& x\in \Omega ,\\ \nabla u\cdot \mathbf{n}=0,& x\in {\Gamma }_N\subset \partial \Omega ,\\ u=0,& x\in {\Gamma }_D\subset \partial \Omega ,\end{array}$$

$$S(x)=A\exp\Big{(}-\frac{||x_{src}-x||^{2}}{2h^{2}}\Big{)}$$

$$\{\begin{array}{ll}\frac{\partial u}{\partial t}={\nabla }^{2}u+S,& x\in \Omega ,t>0,\\ \nabla u\cdot \mathbf{n}=0,& x\in \partial \Omega ,t>0,\\ u=0,& x\in \Omega ,t=0.\end{array}$$

$$S(x)=\begin{cases}\frac{s}{2\pi h^{2}}\exp\Big{(}-\frac{||x_{src}-x||^{2}}{2h^{2}}\Big{)},&\text{if }0\leq t<\tau,\\ 0,&\text{if }t\geq\tau,\end{cases}$$

$${\sigma }_i=0.1+0.1|q_i|\text{for}i=1\dots M,$$

$$\mathcal{O}_{\mathbf{\mathcal{D}}^{z}}=\Bigg{\{}\frac{\hat{N}(q,(\sigma(q))^{2})}{C_{q}}:q\in\mathbf{\mathcal{D}}^{z}\Bigg{\}}.$$

$$\{\begin{array}{ll}-\nabla \cdot \sigma (\mathbf{u})=0,& x\in \Omega =[-5,5]\times [0,2],\\ \mathbf{u}=\mathbf{g},& x\in {\Gamma }_D=\left\{(x,y)\in \Omega |x=0\right\},\\ \sigma (\mathbf{u})\mathbf{n}=\mathbf{t},& x\in {\Gamma }_N=\partial \Omega \backslash {\Gamma }_D,\end{array}$$

Full text

1 Introduction

Experimental data is often used to infer valuable information about parameters for models of physical systems. However, the collection of experimental data can be costly and time consuming. For example, exploratory drilling can reveal valuable information about subsurface hydrocarbon reservoirs, but each well can cost upwards of tens of millions of US dollars. In such situations we can only afford to gather some limited number of experimental data, however not all experiments provide the same amount of information about the processes they are helping inform. Consequently, it is important to design experiments in an optimal way, i.e., to choose some limited number of experimental data to maximize the value of each experiment.

The first experimental design methods employed mainly heuristics, based on concepts such as space-filling and blocking, to select field experiments [14, 15, 32, 33, 38]. While these methods can perform well in some situations, these methods can be improved upon by incorporating any knowledge of the underlying physical processes being inferred or measured. Using physical models to guide experiment selection has been shown to drastically improve the cost effectiveness of experimental designs for a variety of models based on ordinary differential equations [26, 12, 6], partial differential equations [21] and differential algebraic equations [4]. When model observables are linear with respect to the model parameters the alphabetic optimality criteria are often used [19, 1, 17]. For example $A$-optimality to minimize the average variance of parameter estimates, $D$-optimality to maximize the differential Shannon entropy, or $G$-optimality to minimize the maximum variance of model predictions. These criteria have been developed in both Bayesian and non-Bayesian settings [2, 1, 19, 3, 11, 28].

In this paper we focus attention on Bayesian methods for OED that can be applied to both linear and nonlinear models [30, 31, 35]. Specifically we pursue OEDs which are optimal for inferring model parameters on finite-dimensional spaces from experimental data observed at a set of sensor locations. In the context of OED for inference, analogues of the alphabetic criterion, for linear models have also been applied to nonlinear models [2, 29, 18, 21]. In certain situations, for example infinite-dimensional problems (random variables are random fields) or problems with computational expensive models, OED based upon linearizations of the model response and Laplace (Gaussian) approximations of the posterior distribution have been necessary [2, 29]. In other settings non-Gaussian approximations of the posterior have also been pursued [30, 31, 35, 28].

This manuscript presents a new approach for OED based upon consistent Bayesian inference, introduced in [8]. We adopt an approach for OED similar to the approach in [22] and seek an OED that maximizes the expected information gain from the prior to the posterior over the set of possible observational densities. Although our OED framework is Bayesian in nature, this approach is fundamentally different from the statistical Bayesian methods mentioned above. The aforementioned Bayesian OED methods use what we will refer to as the classical/statistical Bayesian approach for stochastic inference (see e.g., [36]) to characterize posterior densities that reflects an assumed error model. In contrast, consistent Bayesian inference assumes a probability density on the observations is given and produces a posterior density that is consistent with the model and the data in the sense that the push-forward of the posterior (through the computational model) matches the observed density almost everywhere. We direct the interested reader to [8] for a discussion on the differences between the consistent and statistical Bayesian approaches. Consistent Bayesian inference has some connections with measure-theoretic inference [7], which was used for OED in [9], but the two approaches make different assumptions and therefore typically give different solutions to the stochastic inverse problem.

The consistent Bayesian approach is appealing for OED since it can be used in an offline-online mode. Consistent Bayesian inference requires an estimate of the push-forward of the prior, which although expensive can be computed offline or obtained from archival simulation data. Once the push-forward of the prior is constructed, the posterior density can be approximated cheaply. Moreover, this push-forward of the prior does not depend on the density on the observations which enables a computationally efficient approach for solving multiple stochastic inverse problems for different densities on the observations. This can significantly reduce the cost of computing the expected information gain if the set of candidate observation is known a priori.

The main objectives in this paper are to derive an OED formulation using the consistent Bayesian framework and to present a computational strategy to estimate the expected information gained for an experimental design. The pursuit of a computationally efficient approach for coupling our OED method with continuous optimization techniques is an intriguing topic that we leave for future work. Here, we consider batch design over a discrete set of possible experiments. Batch design, also known as open-loop design, involves selecting a set of experiments concurrently such that the outcome of any experiment does not effect the selection of the other experiments. Such an approach is often necessary when one cannot wait for the results of one experiment before starting another, but is limited in terms of the number of observations we can consider.

The remainder of this paper is outlined as follows. In Section 2 we summarize the consistent Bayesian method for solving stochastic inverse problems. In Section 3 we discuss the information content of an experiment, and present our OED formulation based upon expected information gain. During the process of defining the expected information gain of a given experimental design, care must be taken to ensure that the model can predict all of our potential observed data. In Section 4 we discuss situations for which this assumption is violated and means for avoiding these situations. Numerical examples are presented in Section 5 and concluding remarks are provided in Section 6.

2 A Consistent Bayes formulation for stochastic inverse

problems

We are interested in experimental designs which are optimal for inferring model parameters from experimental data. Inferring model parameters for a single design and realization of experimental data is a fundamental component of producing such optimal designs. In this section we summarize the consistent Bayes method for parametric inference, originally presented in [8]. Although Bayesian in nature, the consistent Bayesian approach differs significantly from its classical Bayesian counterpart [36, 34, 25] which was used for OED in [22, 24, 10, 1, 2, 29, 28, 23]. We refer the interested reader to [8] for a full discussion of these differences.

2.1 Notation, Assumptions, and a Stochastic Inverse Problem

Let $M(Y,\lambda)$ denote a deterministic model with solution $Y(\lambda)$ that is an implicit function of model parameters $\lambda\in{\Lambda}\subset\mathbb{R}^{n}$. The set ${\Lambda}$ represents the largest physically meaningful domain of parameter values, and, for simplicity, we assume that ${\Lambda}$ is compact. In practice, modelers are often only concerned with computing a relatively small set of quantities of interest (QoI), $\{Q_{i}(Y)\}_{i=1}^{m}$, where each $Q_{i}$ is a real-valued functional dependent on the model solution $Y$. Since $Y$ is a function of parameters $\lambda$, so are the QoI and we write $Q_{i}(\lambda)$ to make this dependence explicit. Given a set of QoI, we define the QoI map $Q(\lambda):=(Q_{1}(\lambda),\cdots,Q_{m}(\lambda))^{\top}:{\Lambda}\to\mathbf{\mathcal{D}}\subset\mathbb{R}^{m}$ where $\mathbf{\mathcal{D}}:=Q({\Lambda})$ denotes the range of the QoI map.

Assume $({\Lambda},\mathcal{B}_{{\Lambda}},\mu_{{\Lambda}})$ and $(\mathbf{\mathcal{D}},\mathcal{B}_{\mathbf{\mathcal{D}}},\mu_{\mathbf{\mathcal{D}}})$ are measure spaces. We assume $\mathcal{B}_{{\Lambda}}$ and $\mathcal{B}_{\mathbf{\mathcal{D}}}$ are the Borel $\sigma$-algebras inherited from the metric topologies on $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$, respectively. The measures $\mu_{{\Lambda}}$ and $\mu_{\mathbf{\mathcal{D}}}$ are volume measures.

We assume that the QoI map $Q$ is at least piecewise smooth implying that $Q$ is a measurable map between the measurable spaces $({\Lambda},\mathcal{B}_{{\Lambda}})$ and $(\mathbf{\mathcal{D}},\mathcal{B}_{\mathbf{\mathcal{D}}})$. For any $A\in\mathcal{B}_{\mathbf{\mathcal{D}}}$, we then have

[TABLE]

Furthermore, $B\subseteq Q^{-1}(Q(B))$ for any $B\in\mathcal{B}_{{\Lambda}}$, although in most cases $B\neq Q^{-1}(Q(B))$ even when $n=m$.

Finally, we assume that an observed probability measure, $P_{\mathbf{\mathcal{D}}}^{\text{obs}}$, is given on $(\mathbf{\mathcal{D}},\mathcal{B}_{\mathbf{\mathcal{D}}})$ and is absolutely continuous with respect to $\mu_{\mathbf{\mathcal{D}}}$, which implies it can be described in terms of an observed probability density, $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$. The stochastic inverse problem is then defined as determining a probability measure, $P_{\Lambda}$, described as a probability density, $\pi_{\Lambda}$, such that, the push-forward measure agrees with $P_{\mathbf{\mathcal{D}}}^{\text{obs}}$. We use $P^{Q(P_{\Lambda})}_{\mathbf{\mathcal{D}}}$ to denote the push-forward of $P_{\Lambda}$ through $Q(\lambda)$, i.e.,

[TABLE]

for all $A\in\mathcal{B}_{\mathbf{\mathcal{D}}}$. Using this notation, a solution to the stochastic inverse problem is defined formally as follows:

Definition 1** (Consistency).**

Given a probability measure $P_{\mathbf{\mathcal{D}}}^{\text{obs}}$ on $(\mathbf{\mathcal{D}},\mathcal{B}_{\mathbf{\mathcal{D}}})$ that is absolutely continuous with respect $\mu_{\mathbf{\mathcal{D}}}$ and admits a density $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$, the stochastic inverse problem seeks a probability measure $P_{\Lambda}$ on $({\Lambda},\mathcal{B}_{{\Lambda}})$ that is absolutely continuous with respect to $\mu_{{\Lambda}}$ and admits a probability density $\pi_{\Lambda}$, such that the subsequent push-forward measure induced by the map, $Q(\lambda)$, satisfies

[TABLE]

for any $A\in\mathcal{B}_{\mathbf{\mathcal{D}}}$. We refer to any probability measure $P_{\Lambda}$ that satisfies (1) as a consistent solution to the stochastic inverse problem.

Clearly, a consistent solution may not be unique, i.e., there may be multiple probability measures that are consistent in the sense of Definition 1. This is analogous to a deterministic inverse problem where multiple sets of parameters may produce the observed data. A unique solution may be obtained by imposing additional constraints or structure on the stochastic inverse problem. In this paper, such structure is obtained by incorporating prior information to construct a unique Bayesian solution to the stochastic inverse problem.

2.2 A Bayesian solution to the stochastic inverse problem

Following the Bayesian philosophy [37], we introduce a prior probability measure $P_{{\Lambda}}^{\text{prior}}$ on $({\Lambda},\mathcal{B}_{{\Lambda}})$ that is absolutely continuous with respect to $\mu_{{\Lambda}}$ and admits a probability density $\pi_{{\Lambda}}^{\text{prior}}$. The prior probability measure encapsulates the existing knowledge about the uncertain parameters.

Assuming that $Q$ is at least measurable, then the prior probability measure on ${\Lambda}$, $P_{{\Lambda}}^{\text{prior}}$, and the map, $Q$, induce a push-forward measure $P_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$ on $\mathbf{\mathcal{D}}$, which is defined for all $A\in\mathcal{B}_{\mathbf{\mathcal{D}}}$,

[TABLE]

We utilize the following expression for the posterior,

[TABLE]

which we describe in terms of a probability density given by

[TABLE]

We note that if $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}=\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$, i.e., if the prior solves the stochastic inverse problem, then the posterior density will be equal to the prior density.

It was recently shown in [8] that the posterior given by (3) defines a consistent probability measure using a contour $\sigma$-algebra. When interpreted as a particular iterated integral of (4), the posterior defines a probability measure on $({\Lambda},\mathcal{B}_{{\Lambda}})$ in the sense of Definition 1, i.e., the push-forward of the posterior matches the observed probability density. Approximating the posterior density using the consistent Bayesian approach only requires an approximation of the push-forward of the prior probability on the model parameters, which is fundamentally a forward propagation of uncertainty. While numerous approaches have been developed in recent years to improve the efficiency and accuracy of the forward propagation of uncertainty using computational models, in this paper we only consider the most basic of methods, namely Monte Carlo sampling, to sample from the prior. We evaluate the computational model for each of the samples from the prior and use a standard kernel density estimator [40] to approximate the push-forward of the prior.

Given the approximation of the push-forward of the prior, we can evaluate the posterior at any point $\lambda\in{\Lambda}$ if we compute $Q(\lambda)$. This provides several possibilities for interogating the posterior. In Section 3.2, we compute $Q(\lambda)$ on a uniform grid of points to visualize the posterior after we compute the push-forward of the prior. This does require additional model evaluations, but visualizing the posterior is rarely required and only useful for illustrative purposes in 1 or 2 dimensions. More often, we are interested in obtaining samples from the posterior. This is also demonstrated in Section 3.2 where the samples from the prior are either accepted or rejected using a standard rejection sampling procedure. For a given $\lambda$, we compute the ratio $\pi_{{\Lambda}}^{\text{post}}(\lambda)/(M\pi_{{\Lambda}}^{\text{prior}}(\lambda))$, where $M$ is an estimate of the maximum of the ratio over ${\Lambda}$, and compare this value with a sample, $\eta$, drawn from a uniform distribution on (0,1). If the ratio is larger than $\eta$, then we accept the sample. We apply the accept-reject algorithm to the samples from the prior and therefore the samples from the posterior are a subset of the samples used to compute the push-forward of the prior. Since we have already computed $Q(\lambda)$ for each of these samples, the computational cost to select a subset of the samples for the posterior is minimal. However, in the context of OED we are primarily interested in computing the information gained from the prior to the posterior which only involves integrating with respect to the prior (see Section 3.1) and does not require additional model evaluations or rejection sampling.

In practice, we prefer to use data that is sensitive to the parameters since otherwise it is difficult to infer useful information about the uncertain parameters. Specifically, if $m\leq n$ and the Jacobian of $Q$ is defined a.e. in ${\Lambda}$ and is full rank a.e., then the push-forward volume measure $\mu_{\mathbf{\mathcal{D}}}$ is absolutely continuous with respect to the Lebesgue measure [8].

For the rest of this work we maintain the following assumptions needed to produce a unique consistent solution to the stochastic inverse problem:

(A1) We have a mathematical model and a description of our prior knowledge about the model input parameters,

(A2) The data exhibits sensitivity to the parameters a.e. in ${\Lambda}$, hence, we use the Lebesgue measure $\mu$ as the volume measure on the data space,

(A3) The observed density is absolutely continuous with respect to the push-forward of the prior.

The assumption concerning the absolute continuity of the observed density with respect to the prior is essential to define a solution to the stochastic inverse problem [8]. While this assumption may appear rather abstract, it simply assures that the prior and the model can predict, with non-zero probability, any event that we have observed. Since the observed density and the model are assumed to be fixed, this is only an assumption on the prior.

In the remainder of this work, we focus on quantifying the value of these posterior densities. We use the Kullback-Leibler divergence [27, 39], to measure the information gained about the parameters from the prior to the posterior. We compute the expected information gain of a given set of QoI (a given experimental design), and then determine the OED to deploy in the field.

3 The information content of an

experiment

We are interested in finding the OED for inferring model input parameters. Conceptually, a design is informative if the posterior distribution of the model parameters is significantly different from the prior. To quantify the information gain of a design we use the Kullback-Leibler (KL) divergence [39] as a measure of the difference between a prior and posterior distribution. While the KL divergence is by no means the only way to compare two probability densities, it does provide a reasonable measure of the information gained in the sense of Shannon information [13] and is commonly used in Bayesian OED [24]. In this section we discuss how to compute the KL divergence and define our OED formulation based upon expected information gain over a specific space of possible observed densities.

3.1 Information gain: Kullback-Leibler divergence

Suppose we are given a description of the uncertainty on the observed data in terms of a probability density $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$. This produces a unique solution to the stochastic inverse problem $P_{{\Lambda}}^{\text{post}}$ that is absolutely continuous with respect to the Lebesgue measure $\mu_{{\Lambda}}$ [8] and admits a probability density, $\pi_{{\Lambda}}^{\text{post}}$. The KL divergence of the posterior from the prior (information gain), denoted $I_{Q}$, is given by

[TABLE]

Note that because $\pi_{{\Lambda}}^{\text{prior}}$ is fixed, $I_{Q}$ is simply a function of the posterior

[TABLE]

and from Eq. (4) the posterior is a function of the observed density. Therefore, we write $I_{Q}$ as a function of the observed density,

[TABLE]

The observation that $I_{Q}$ is a function of only $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ allows us to define the expected information gain in Section 3.3 based on a specific space of observed densities.

Given a high dimensional parameter space, it may be computationally infeasible to accurately approximate the integral in Eq. (5). For example, a multi-variate normal density with unit variance in 100-dimensions has a maximum value of $(1/\sqrt{2\pi})^{100}\approx 1\times 10^{-40}$. However, we may write this integral in terms of densities on the data space evaluated at $Q(\lambda)$ as follows

[TABLE]

where the second equality comes from a simple substitution using Eq. 4. Given a set of samples from the prior, we only need to compute the push-forward of the prior in the data space to approximate $I_{Q}$. This observation provides an efficient method for approximating $I_{Q}$ given a high dimensional parameter space and a low dimensional data space. In fact, we found it convenient to use (8) whenever the prior is not uniform. In the consistent Bayesian formulation, we evaluate the model at the samples generated from the prior to estimate the push-forward of the prior. It is a computational advantage to also use these samples to integrate with respect to the prior rather than integrating with respect to the volume measure which would require additional model evaluations.

3.2 A motivating nonlinear system

Consider the following 2-component nonlinear system of equations with two parameters introduced in [7]:

[TABLE]

The first QoI is the second component, i.e., $Q_{1}(\lambda)=x_{2}(\lambda)$. The parameter ranges are given by $\lambda_{1}\in[0.79,0.99]$ and $\lambda_{2}\in[1-4.5\sqrt{0.1},1+4.5\sqrt{0.1}]$ which are chosen as in [7] to induce an interesting variation in the QoI. We assume the observed density on $Q_{1}$ is a truncated normal distribution with mean 0.3 and standard deviation of 0.01, see Figure 1 (right).

We generate 40,000 samples from the uniform prior and use a kernel density estimator (KDE) to construct an approximation to the resulting push-forward density, see Figure 1 (right). Then we use Eq. (4) to construct an approximation to the posterior density using the same 40,000 samples, see Figure 1 (left), and a simple accept/reject algorithm to generate a set of samples from the posterior, see Figure 1 (middle). We propagate this set of samples from the posterior through the model and approximate the resulting push-forward of the posterior density using a KDE. In Figure 1 (right) we see the push-forward of the posterior agrees quite well with the observed density. Notice the support of the posterior lies in a relatively small region of the parameter space. The information gain from this posterior is $I_{Q_{1}}(\pi_{\mathbf{\mathcal{D}}}^{\text{obs}})\approx 2.015$.

Next, we consider a different QoI to use in the inverse problem, and compare the support of its posterior to the one we just observed. Specifically consider,

[TABLE]

We assume the observed density on $Q_{2}$ is a truncated normal distribution with mean 1.015 and standard deviation of 0.01. We approximate the push-forward density and the posterior using the same 40,000 samples and again generate a set of samples from the posterior and propagate these samples through the model to approximate the push-forward of the posterior, see Figure 2.

Although both $Q_{1}$ and $Q_{2}$ have the same standard deviation in their observed densities, clearly the two QoI produce very different posterior densities. The posterior corresponding to data from $Q_{2}$ has a much larger region of support within the parameter space compared to that of the posterior corresponding to $Q_{1}$. This is quantified with the information gain from this posterior $I_{Q_{2}}(\pi_{\mathbf{\mathcal{D}}}^{\text{obs}})\approx 0.466$. Given these two maps, $Q_{1}$ and $Q_{2}$, and the specified observed data on each of these data spaces, the data $Q_{1}$ is more informative of the parameters than the data $Q_{2}$.

Next, we consider using the data from both $Q_{1}$ and $Q_{2}$, $Q:{\Lambda}\to(Q_{1},Q_{2})$, with the same means and standard deviations as specified above. Again, we approximate the push-forward density and the posterior using the same 40,000 samples, see Figure 3. With the information from both $Q_{1}$ and $Q_{2}$ we see a substantial decrease in the support of the posterior density. Intuitively, the support of the posterior using both $Q_{1}$ and $Q_{2}$ is the support of the posterior using $Q_{1}$ intersected with the support of the posterior using $Q_{2}$. This is quantified in the information gain of this posterior $I_{Q}(\pi_{\mathbf{\mathcal{D}}}^{\text{obs}})\approx 2.98$.

In the scenario in which we can afford to gather data on both $Q_{1}$ and $Q_{2}$, we benefit greatly in terms of reducing the uncertainties on the model input parameters. However, suppose we could only afford to gather one of these QoI in the field. Based on the information gain from each posterior, $Q_{1}$ is more informative about the parameters than $Q_{2}$. However, consider a scenario in which the observed data has different means in both $Q_{1}$ and $Q_{2}$. Due to the nonlinearities of the maps, it is not necessarily true that $Q_{1}$ is still more informative than $Q_{2}$. If we do not know the mean of the data for either $Q_{1}$ or $Q_{2}$, then we want to determine which of these QoI we expect to produce the most informative posterior.

3.3 Expected information gain

Optimal experimental design must select a design before experimental data becomes available. In the absence of data we use the simulation model to quantify the expected information gain of a given experimental design. Let $\mathcal{O}$ denote the space of densities over $\mathbf{\mathcal{D}}$. We want to define the expected information gain as some kind of average over this density space in a meaningful way. However, this is far too general of a space to use to define the expected information gain. This space includes densities that are unlikely to be observed in reality. Therefore, we restrict $\mathcal{O}$ to be a space more representative of densities that may be observed in reality.

With no experimental data available to specify an observed density on a single QoI, we assume the density is a truncated Gaussian with a standard deviation determined by some estimate of the measurement instrument error. With Gaussians of (possibly) varying standard deviations specified for each QoI, this defines the shape of the observed densities we consider. We let $\mathcal{O}_{\mathbf{\mathcal{D}}}$ denote the space of all densities of this shape centered in $\mathbf{\mathcal{D}}=Q({\Lambda})$,

[TABLE]

where $\hat{N}(q,\sigma^{2})$ is a truncated Gaussian function with mean $q$ and standard deviation $\sigma$. More details of this definition of $\mathcal{O}_{\mathbf{\mathcal{D}}}$ are addressed in Section 4. We can easily generalize our description of $\mathcal{O}$. For example, we could also consider the standard deviation of the observed data to be uncertain, in which case we would also average over some interval of possible values for $\sigma$. However, in this work we only vary the center of the Gaussian densities.

Remark 1**.**

We can restrict $\mathcal{O}$ in other ways as well. For example, if we expect the uncertainty in each QoI to be described by a uniform density, then we define the restriction on $\mathcal{O}$ accordingly. This choice of characterization of the observed density space is largely dependent on the application. The only limitation is that we require the measure specified on the observed density space to be defined in terms of the push-forward measure, $P_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$, as described below. In Section 5.2 we describe one approach for defining a restricted observed density space where the observed density of each QoI has a Gaussian profile and the standard deviations are functions of the magnitudes of each QoI.

The restriction of possible $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ to this specific space of densities allows us to represent each density uniquely with a single point $q\in\mathbf{\mathcal{D}}$. Based on our prior knowledge of the parameters and the sensitivities of the map $Q$, the model informs us that some data are more likely to be observed than other data, this is seen in the plot of $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$ in Figure 3 (upper left). This implies we do not want to average over $\mathbf{\mathcal{D}}$ with respect to $\mu$ or $\mu_{\mathbf{\mathcal{D}}}$, but rather with respect to the push-forward of the prior on $D$, $P_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$. This respects the prior knowledge of the parameters and the sensitivity information provided by the model. We define the expected information gain, denoted $E(I_{Q})$, as just described,

[TABLE]

From Eq. (5), $I_{Q}$ itself is defined in terms of an integral. The expanded form for $E(I_{Q})$ is then an iterated integral,

[TABLE]

where we make explicit that $\pi_{{\Lambda}}^{\text{post}}$ is a function of the observed density and, by our restriction of the space of observed densities in Eq. (9), therefore a function of $q\in\mathbf{\mathcal{D}}$. We utilize Monte Carlo sampling to approximate the integral in Eq. (10) as described in Algorithm 1.

Algorithm 1 appears to be a computationally expensive procedure since it requires solving $M$ stochastic inverse problems and, as noted in [8], approximating $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$ can be expensive. In [8] and in this paper we use kernel density estimation techniques to approximate $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$ which does not scale well as the dimension of $\mathbf{\mathcal{D}}$ increases [40]. On the other hand, for a given experimental design, we only need to compute this approximation once, as each $I_{Q}$ in Step 4 of Algorithm 1 is computed using the same prior and map $Q$ and, therefore, the same $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$. In other words, the fact that the consistent Bayes method only requires approximating the push-forward of the prior implies that this information can be used to approximate posteriors for different observed densities without requiring additional model evaluations. This significantly improves the computational efficiency of the consistent Bayesian approach in the context of OED. We leverage this computational advantage throughout this paper by considering a discrete set of designs which allows us to compute the push-forward for all of the candidate designs simultaneously. Utilizing a continuous design space might require computing the push-forward for each iteration of the optimization algorithm, since the designs (locations of the observations) are not known a priori. The additional model simulations required to compute the push-forward of the prior at new design points might be intractable if the number of iterations is large, however the need for new simulations may be avoided, if the new observations can be extracted from archived state-space data. For example, if one stores the finite element solutions of a PDE at all samples of the prior at the first iteration of the design optimization, one can evaluate obsevations at new design locations, which are functionals of this PDE solution, via interpolation using the finite element basis.

3.4 Defining the OED

We are now in a position for define our OED formulation. Recall that our experimental design is defined as the set of QoI computed from the model and we seek the optimal set of QoI to deploy in the field. Given a physics based model, prior information on the model parameters, a space of potential experimental designs, and a generic description of the uncertainties for each QoI, we define our OED as follows.

Definition 2** (OED).**

Let $\mathcal{Q}$ represent the design space, i.e., the space of all possible experimental designs, and $Q^{z}\in\mathcal{Q}$ be a specific design. Then the OED is the $Q^{z}\in\mathcal{Q}$ that maximizes the expected information gain,

[TABLE]

As previously mentioned, the focus in this paper is on the utilization of the consistent Bayesian methodology within the OED framework, so we do not explore different approaches for solving the optimization problem given by Definition 2 and simply find the optimal design over a discrete set of candidate designs.

Remark 2**.**

Consistent Bayesian inference is potentially well suited to finding OED in continuous design spaces. Typically OED based upon statistical Bayesian methods uses Markov Chain Monte Carlo (MCMC) methods to characterize the posterior distribution. MCMC methods do not provide a functional form for the posterior but rather only provide samples from the posterior. Consequently, gradient-free or stochastic gradient-based optimization methods must be used to find the optimal design. In contrast consistent Bayesian inference provides a functional form for the posterior which allows the use of more efficient gradient based optimizers. Exploring the use of more efficient continuous optimization procedures will be the subject of future work.

4 Infeasible data

The OED procedure proposed in this manuscript is based upon consistent Bayesian inference which requires that the observed measure is absolutely continuous with respect to the push-forward measure induced by the prior and the model (assumption A3). In other words, any event that we observe with non-zero probability will be predicted using the model and prior with non-zero probability. During the process of computing $E(I_{Q})$, it is possible that we violate this assumption. Specifically, depending on the mean and variance of the observational density we may encounter $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}\in\mathcal{O}_{\mathbf{\mathcal{D}}}$ such that $\int_{\mathbf{\mathcal{D}}}\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}d\mu\ <1$, i.e., support of $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ extends beyond the range of the map $Q$, see Figure 5 (upper right). In this section we discuss the causes of infeasible data and options for avoiding infeasible data when estimating an optimal experimental design.

4.1 Infeasible data and consistent Bayesian inference

When inferring model parameters using consistent Bayesian inference the most common cause for infeasible data is that the model being used to estimate the OED is inadequate. That is, the deviation between the computational model and reality is large enough to prohibit the model from predicting all of the observational data. The deviation between the model prediction and the observational data is often referred to as model structure error and can often be a major source of uncertainty. This is an issue of most if not all inverse parameter estimation problems [25]. Recently there has been a number of attempts to quantify this error (see e.g., [34]) however such approaches are beyond the scope of this paper. In the following we will assume that the model structure error does not prevent the model from predicting all the observational data.

4.2 Infeasible data and OED

To estimate an approximate OED we must quantify the expected information gain of a given experimental design (see Section 3.3). The expectation is over all possible normal observation densities with mean $q\in\mathcal{D}$ and variance $\sigma$, defined by the space (9). When the support of $\mathcal{D}$ is bounded these densities may produce infeasible data. The effect of this violation increases as $q$ approaches the boundary of $\mathcal{D}$.

To remedy this violation of (A3) we must modify the set of observational densities. In this paper we choose to normalize $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ over $\mathbf{\mathcal{D}}$. We redefine the observed density space $\mathcal{O}_{\mathbf{\mathcal{D}}}$ so that (A3) holds for each density in the space,

[TABLE]

where $\hat{N}(q,\sigma^{2})$ is a truncated Gaussian function with mean $q$ and standard deviation $\sigma$, and $C_{q}$ is the integral of $\hat{N}(q,\sigma^{2})$ over $\mathbf{\mathcal{D}}$ with respect to the Lebesgue measure on $\mathbf{\mathcal{D}}$,

[TABLE]

A similar approach for normalizing Gaussian densities over compact domains was taken in [5].

4.3 A nonlinear model with infeasible data

In this section, we use the nonlinear model introduced in Section 3.2 to demonstrate that infeasible data can arise from relatively benign assumptions. Suppose the observed density on $Q_{1}$ is a truncated normal distribution with mean 0.3 and standard deviation of 0.04. In this one dimensional data space, this observed density is absolutely continuous with respect to the push-forward of the prior on $Q_{1}$, see Figure 4 (left). Next, suppose the observed density on $Q_{2}$ is a truncated normal distribution with mean 0.982 and standard deviation of 0.01. Again, in this new one dimensional data space, this observed density is absolutely continuous with respect to the push-forward of the prior on $Q_{2}$, see Figure 4 (right). Both of these observe densities are dominated by their corresponding push-forward densities, i.e., the model can reach all of the observed data in each case.

However, consider the data space defined by both $Q_{1}$ and $Q_{2}$ and the corresponding push-forward and observed densities on this space, see Figure 5. The non-rectangular shape of the combined data space is induced by the nonlinearity in the model and the correlations between $Q_{1}$ and $Q_{2}$. As we see in Figure 5, the observed density using the product of the 1-dimensional Gaussian densities is not absolutely continuous with respect to the push-forward density on $(Q_{1},Q_{2})$, i.e., the support of $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ extends beyond the support of $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$. Referring to Eq. (13), we normalize this observed density over $\mathbf{\mathcal{D}}$, see Figure 5 (right). Now that the new observed density obeys the assumptions needed, we could solve the stochastic inverse problem as described in Section 2.

4.4 Computational considerations

The main computational challenge in the consistent Bayesian approach is the approximation of the push-forward of the prior. Following [8], we use Monte Carlo sampling for the forward propagation of uncertainty. While the rate of convergence is independent of the number of parameters (dimension of ${\Lambda}$), the accuracy in the statistics for the QoI may be relatively poor unless a large number of samples can be taken. Alternative approaches based on surrogate models can significantly improve the accuracy, but are generally limited to small number of parameters. We also employ kernel density estimation techniques to construct a non-parametric approximation of the push-forward density, but it is well-known that these techniques do not scale well with the number of observations (dimension of $\mathbf{\mathcal{D}}$) [40].

Next, we address the computational issue of normalizing $\hat{N}(q,\sigma^{2})$, i.e., $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$, over $\mathbf{\mathcal{D}}$. From the plot of $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$ in Figure 5 (left) it is clear the data space may be a complex region. Normalizing $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$, as in Figure 5 (right), over $\mathbf{\mathcal{D}}$ would be computationally expensive. Fortunately, the consistent Bayesian approach provides a means to avoid this expense. Note that from Eq. (3) we have,

[TABLE]

where $P_{{\Lambda}}^{\text{prior}}({\Lambda})=P_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}(Q({\Lambda}))=1$ which implies,

[TABLE]

Therefore, normalizing $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ over $\mathbf{\mathcal{D}}$ is equivalent to solving the inverse problem and then normalizing $\tilde{\pi}_{{\Lambda}}^{\text{post}}$ (where we use the tilde over $\pi$ to indicate this function does not integrate to 1 because we have violated (A3)) over ${\Lambda}$. Although ${\Lambda}$ may not always be a generalized rectangle, (A1) implies we have a clear definition of ${\Lambda}$ and therefore can efficiently integrate $\tilde{\pi}_{{\Lambda}}^{\text{post}}$ over ${\Lambda}$ and then normalize $\tilde{\pi}_{{\Lambda}}^{\text{post}}$ by

[TABLE]

In fact, this normalization factor can be estimated without additional model evaluations and without using the values of the prior or the posterior, which may not be usable in high-dimensional spaces. We observe that

[TABLE]

Thus, we can use the values of $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ and $\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$ computed for the samples generated from the prior, which were used to estimate the push-forward of prior, to integrate $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}/\pi_{\mathbf{\mathcal{D}}}^{Q(\text{prior})}$ with respect to the prior.

5 Numerical examples

In this section we consider several models of physical systems. First, we consider a stationary convection-diffusion model with a single uncertain parameter controlling the magnitude of the source term. Next, we consider a transient transport model with a two dimensional parameter space determining the location of the source of a contaminant. Then, we consider a inclusion problem in computational mechanics where two uncertain parameters control the shape of the inclusion. Finally, we consider a high-dimensional example of single-phase incompressible flow in porous media where the uncertain permeability field is given by a Karhunen-Loeve expansion [42].

In each example, we have a parameter space ${\Lambda}$, a set of possible QoI, and a specified number of QoI we can afford to gather during the experiment. This in turn defines a design space $\mathcal{Q}$ and we let $Q^{z}\in\mathcal{Q}$ represent a single experimental design and $\mathbf{\mathcal{D}}^{z}=Q^{z}({\Lambda})$ the corresponding data space. For each experimental design, we let $\sigma^{z}$ represent the standard deviations defined by the uncertainties in each QoI that compose $Q^{z}$ and $\mathcal{O}_{\mathbf{\mathcal{D}}^{z}}$ represent the observed density space.

All of these examples have continuous design spaces, so we approximate the OED by selecting the OED from a large set of candidate designs. This approach was chosen because it is much more efficient to perform the forward propagation of uncertainty using random sampling only once and to compute all of the candidate measurements for each of these random samples. Alternatively, one could pursue a continuous optimization formulation which would require a full forward propagation of uncertainty for each new design. As mentioned in Section 3.4, one could limit the number of designs using a gradient-based or Newton-based optimization approach, but this is beyond the scope of this paper.

5.1 Stationary convection-diffusion: uncertain source amplitude

In this section we consider a convection-diffusion problem with a single uncertain parameter controlling the magnitude of a source term. This example serves to demonstrate that the OED formulation gives intuitive results for simple problems.

5.1.1 Problem setup

Consider a stationary convection diffusion model on a square domain:

[TABLE]

with

[TABLE]

where $\Omega=[0,1]^{2}$, $u$ is the concentration field, the diffusion coefficient $D=0.01$, the convection vector $v=[1,1]$, and $S$ is a Gaussian source with the following parameters: $x_{src}$ is the location, $A$ is the amplitude, $h$ is the width. We impose homogeneous Neumann boundary conditions on $\Gamma_{N}$ (right and top boundaries) and homogeneous Dirichlet conditions on $\Gamma_{D}$ (left and bottom boundaries). For this problem, we choose $x_{src}=[0.5,0.5]$, and $h=0.05$. We let $A$ be uncertain within $[50,150]$, thus the parameter space for this problem is ${\Lambda}=[50,150]$. Hence, our goal is to gather some limited amount of data that provides the best information about the amplitude of the source, i.e., reduces our uncertainty in $A$. To approximate solutions to the PDE in Eq. 18 given a source amplitude $A$, we use a finite element discretization with continuous piecewise bilinear basis functions defined on a uniform ($25\times 25$) spatial grid.

5.1.2 Results

We assume that we have limited resources for gathering experimental data, specifically, we can only afford to place one sensor in the domain to gather a single concentration measurement. Our goal is to place this single sensor in $\Omega$ to maximize the expected information gained about the amplitude of the source. We discretize $\Omega$ using 2,000 uniform random points which produces a design space with 2,000 possible experimental designs. For this problem, we let the uncertainty in each QoI be described by a truncated Gaussian profile with a fixed standard deviation of 0.1. This produces observed density spaces, $\mathcal{O}_{\mathbf{\mathcal{D}}^{z}}$, as described in Eq. 13.

We generate 5,000 uniform samples from the prior and simulate measurements of each QoI for each of these 5,000 samples. We consider approximate solutions to the OED problem using subsets of the 5,000 samples of size 50, 200, 1,000 and 5,000. For each experimental design, we calculate $E(I_{Q^{z}})$ using Algorithm 1 and plot $E(I_{Q^{z}})$ as a function of the discretized design space in Figure 6. Notice the expected information gain is greatest near the center of the domain (near the location of the source) and in the direction of the convection vector away from the source. This result matches intuition, as we expect data gathered in regions of the domain that exhibit sensitivity to the parameters to produce high expected information gains.

We note that, for this example, a sufficiently accurate approximation to the design space and the OED is obtained using only 50 samples corresponding to 50 model evaluations. In Table 1 we show the top 5 experimental designs (computed using the full set of 5,000 samples) and corresponding $E(I_{Q^{z}})$ for each set of samples.

5.2 Time dependent diffusion: uncertain source location

In this section, we compare results from a statistical Bayesian formulation of OED to the formulation described in this paper. Specifically, we consider the model in [22] where the author uses a classical Bayesian framework for OED to determine the optimal placement of a single sensor that maximizes the expected information about the location of a contaminant source.

5.2.1 Problem setup

Consider a contaminant transport model on a square domain:

[TABLE]

with

[TABLE]

where $\Omega=[0,1]^{2}$, $u$ is the space-time concentration field, we impose homogeneous Neumann boundary conditions along with a zero initial condition, and $S$ is a Gaussian source with the following parameters: $x_{src}$ is the location, $s$ is the intensity, $h$ is the width, and $\tau$ is the shutoff time.

Our goal is to gather some limited amount of data that provides the best information about the location of the source, i.e., reduces our uncertainty in $x_{src}$. For this problem, we choose $s=2.0,\,h=0.05,$ and $\tau=0.3$ and let $x_{src}$ be uncertain within $[0,1]^{2}$ such that ${\Lambda}=[0,1]^{2}$. To approximate solutions to the PDE in Eq. 19 given a location of $S$, i.e., a given $x_{src}$, we use a finite element discretization with continuous piecewise bilinear basis functions defined on a uniform ($25\times 25$) spatial grid and backward Euler time integration with a step size $\Delta t=0.004$ (100 time steps).

5.2.2 Results

We assume that we have limited resources for gathering experimental data, specifically, we can only afford to place one sensor in the domain and can only gather a single concentration measurement at time $t=0.24$. Our goal is to place this single sensor in $\Omega$ to maximize the expected information gained about the location of the contaminant source. For simplicity, we discretize $\Omega$ using an $11\times 11$ regular grid of points which produces a design space with 121 possible experimental designs. We let the uncertainty in each QoI be described by a Gaussian profile with a standard deviation that is a function of the magnitude of the QoI, i.e.,

[TABLE]

where $M$ is the dimension of the data space. This produces observed density spaces, $\mathcal{O}_{\mathbf{\mathcal{D}}^{z}}$, that consist of truncated Gaussian functions with varying standard deviations,

[TABLE]

We generate 5,000 uniform samples from the prior and simulate measurements of each QoI for each of these 5,000 samples. We consider approximate solutions to the OED problem using subsets of the 5,000 samples of size 50, 200, 1,000 and 5,000. For each experimental design, we use this data to calculate $E(I_{Q^{z}})$ using Algorithm 1 and plot $E(I_{Q^{z}})$ as a function of the discretized design space in Figure 7. Notice the expected information gain is greatest near the corners of the domain and smallest near the center, this is consistent with [22]. In Table 2 we show the top 5 experimental designs, approximated using the full set of 5,000 samples, and corresponding $E(I_{Q^{z}})$ for each set of samples.

In Figure 8 we consider three different posteriors computed using data from the OED approximated using 5,000 samples, i.e., data gathered by a sensor placed in the bottom left corner of the domain, where each posterior corresponds to a different possible location of the source. We see varying levels of information gain in these three scenarios, reiterating the point that we choose the OED based on the average of these information gains, $E(I_{Q})$.

Remark 3**.**

*Although many of the results in this section seem to match our intuition about which measurement locations should produce high expected information gains, this may not always be the case. In particular, we have found that our results can depend on our choice of the variance in the the observed densities $\sigma$. If $\sigma$ is chosen to be large relative to the range of a data space, then the posteriors produced as we average over $\mathcal{O}_{\mathbf{\mathcal{D}}}$ are all nearly the same and potentially produce unusually high information gains when the observed densities have substantial support over regions of the data space with very small probability (very small values of the push-forward of the prior). Another way to think of this is the push-forward densities have high entropy and because $\sigma$ is large $\pi_{\mathbf{\mathcal{D}}}^{\text{obs}}$ is very close to uniform and this produces posterior densities with high information gains. If $\sigma$ is chosen to be small relative to the range of the data space, i.e., if we expect the experiments to be informative, we do not encounter this issue because we are integrating over $\mathbf{\mathcal{D}}$ with respect to the push-forward measure so most of our potential observed data lies in high probability regions of the data space. *

5.3 A Parameterized Inclusion

In this section, we consider a simple problem in computational mechanics where the precise boundary of an inclusion is uncertain. We parameterize the inclusion and seek to determine the location to place a sensor that will maximize the information gained regarding the shape of the inclusion. We use a linear elastic formulation to model the response of the media to surface forces and measure horizontal stress at each sensor location. We assume that the material properties (Poisson ratio and Young’s modulus) are different inside the inclusion and that these properties are known a prior.

5.3.1 Problem setup

Consider a linear elastic plane strain model,

[TABLE]

where $\mathbf{\sigma(u)}$ is given by the linear elastic constitutive relation,

[TABLE]

We express this relation in terms of the Lamé parameters, $\lambda$ and $\mu$, which are related to the Poisson ratio, $\nu$, and Young’s modulus, $E$, via the following expressions,

[TABLE]

Now assume that there is an inclusion within the media defined by an ellipse

[TABLE]

where $x_{0}=y_{0}=0$ and $\alpha$ is uniformly distributed on $[0.5,1]$ and $\beta$ is uniformly distributed on $[0.25,0.5]$. The material properties are assumed to be known and are given by,

[TABLE]

These material properties were not chosen to emulate any particular materials, just to demonstrate the proposed OED formulation.

Next let us impose homogeneous Dirichlet boundary conditions on the bottom boundary and stress free boundary conditions on the sides, and impose a uniform traction in the y-direction along the top boundary ($\mathbf{t}_{\text{top}}=(0,-1)^{T}$). And finally assume that we can probe the media and measure the horizontal stress at a given sensor location. We do not want to puncture the inclusion, so we only consider sensor locations outside the bounds on the inclusion. Equation 22 was solved using a finite element discretization with piecewise linear basis functions defined on a uniform $400\times 80$ mesh resulting in a system with 64,962 degrees of freedom. The computational model is implemented using the Trilinos toolkit [20] and each realization of the model requires approximately 1 second using 8 processors.

5.3.2 Results

As in previous examples, we assume that we have limited resources for gathering experimental data, specifically, we can only afford to place one sensor in the domain to gather a single stress measurement. Our goal is to place this single sensor to maximize the expected information gained about the shape of the inclusion. We select 2,000 random sensor locations (outside the inclusion bounds) which produces a design space with 2,000 possible experimental designs. For this problem, we let the probability density for the QoI be described by a truncated Gaussian profile with a fixed standard deviation of 0.001. We generate 1,000 uniform samples from the prior and compute the horizontal stress at each sensor location for each of these 1,000 samples.

First, we compare the posterior densities for two sensor locations, $(3.5294,1.3049)$ and $(1.3902,1.2100)$, under the assumption that we have already gathered data at these sensor locations. The purpose here is to demonstrate that we obtain different posterior densities and therefore gain different information from each sensor. The first sensor is further from the inclusion so we expect that the data from the second sensor will constrain the posterior more than the data from the first. In Figures 10 and 11, we plot the samples from the posterior and the corresponding kernel density estimate of the posterior for the first and second sensor locations respectively. It is clear that measuring the horizontal stress closer to the inclusion increases the information gained from the prior to the posterior.

We consider approximate solutions to the OED problem using subsets of the 1,000 samples of size 10, 50, 100 and 1,000. For each experimental design, we use this data to calculate $E(I_{Q^{z}})$ using Algorithm 1 and plot $E(I_{Q^{z}})$ as a function of the discretized design space in Figure 12. Notice the expected information gain is greatest near the bottom of the domain near the inclusion and is reasonably symmetric around the inclusion. Also note that the expected information gain is relatively large in the bottom corners of the domain. This is due to the choice of boundary conditions for the model which induces a large amount of stress in these corners. In Table 2 we show the top 5 experimental designs, approximated using the full set of 1,000 samples, and corresponding $E(I_{Q^{z}})$ for each set of samples.

5.4 A Higher-Dimensional Porous Media Example with Uncertain Permeability

In this section, we consider an example of single-phase incompressible flow in porous media with a Karhunen-Loéve expansion of the uncertain permeability field. The purpose of this example is to demonstrate the OED formulation on a problem with a high-dimensional parameter space and more than one sensor.

5.4.1 Problem setup

Consider a single-phase incompressible flow model:

[TABLE]

Here, $p$ is the pressure field and $K$ is the permeability field which we assume is a scalar field given by a Karhunen-Loéve expansion of the log transformation, $Y=\log{K}$, with

[TABLE]

where $\overline{Y}$ is the mean field. We assume the mean removed random media is given by a Gaussian process which implies that the $\xi_{i}$ are mutually uncorrelated random variables with zero mean and unit variance [16, 41]. The eigenvalues, $\eta_{i}$, and eigenfunctions, $f_{i}$, are computed numerically using the following covariance function,

[TABLE]

where $\sigma_{Y}$ and $\eta_{i}$ denote the variance and correlation length in the i${}^{\text{th}}$ spatial direction respectively. We assume a correlation length of $0.01$ in each spatial direction and truncate the expansion at 100 terms. This choice of truncation is purely for the sake of demonstration. In practice, the expansion is truncated once a sufficient fraction of the energy in the eigenvalues is retained [42, 16]. This truncation gives 100 uncorrelated random variables, $\xi_{1},\ldots,\xi_{1}00$, with zero mean and unit variance which implies ${\Lambda}=\mathbb{R}^{1}00$. To approximate solutions to the PDE in Eq. 23 we use a finite element discretization with continuous piecewise bilinear basis functions defined on a uniform ($50\times 50$) spatial grid.

5.4.2 Results

In this section, we present approximate solutions to several different design problems. We begin with the familiar problem of choosing a single sensor location within the physical domain. Then, we consider approximating the optimal location of a second sensor given the location of the first sensor. In this way, we solve the greedy OED problem and determine the greedy optimal locations of 1-8 sensors within the physical domain. We then consider solving the ehaustive OED problem where we limit the sensors to 25 locations and consider determining the optimal location of 5 available sensors.

First, assume that we have limited resources for gathering experimental data, specifically, we can only afford to place one sensor in the domain to gather a single pressure measurement. Our goal is to place this single sensor in $\Omega$ to maximize the expected information gained about the amplitude of the source. We discretize $\Omega$ using 1,301 pionts on a grid which produces a design space with 1,301 possible experimental designs. For this problem, we let the uncertainty in each QoI be described by a truncated Gaussian profile with a fixed standard deviation of 0.01. This produces observed density spaces, $\mathcal{O}_{\mathbf{\mathcal{D}}^{z}}$, as described in Eq. 13.

We generate 10,000 samples from the prior and simulate measurements of each QoI. We consider approximate solutions to the OED problem using subsets of the 10,000 samples of size 50, 100, 1,000 and 10,000. For each experimental design, we calculate $E(I_{Q^{z}})$ using Algorithm 1 and plot $E(I_{Q^{z}})$ as a function of the discretized design space in Figure 13. Notice the expected information gain is greatest near the top and bottom of the domain away from the left and right edges. This result matches intuition, as we expect data gathered near the left and right edges to be less informative given the Dirichlet boundary condition imposed on those boundaries. We note that, for this example, a sufficiently accurate approximation to the design space and the OED is obtained using only 1,000 samples corresponding to 1,000 model evaluations. In Table 1 we show the top 5 experimental designs (computed using the full set of 10,000 samples) and corresponding $E(I_{Q^{z}})$ for each set of samples.

Next, we consider the greedy OED problem of placing 8 sensors within the physical domain. We choose to use all of the available 10,000 samples to solve this problem. In Figure 14, we see the design space as a function of the previously determined locations of placed sensors. We observe a strong symmetry to this problem, as is expected due to the symmetry of the physical process defined on this domain with the given boundary conditions. In the bottom right of Figure 14, notice the very small range of the color bar indicating the possible values of the expected information gain. This suggests, for this example, we expect there is a limit on the number of useful sensor locations for informing likely parameter values.

Lastly, we consider the exhaustive OED problem of placing 5 sensors within the physical domain and, for computational feasibility, restrict the possible locations of these 5 sensors to 25 points in the physical domain, see Figure 15. We choose to use 1,000 samples to solve this problem. In Figure 15, we plot the design space for a single sensor location using these 1,000 samples and show the optimal location of 1, 2, 3, 4 and 5 sensors. The results are quite similar to the greedy results previously described.

6 Conclusion

In this manuscript, we developed an OED formulation based on the recently developed consistent Bayesian approach for solving stochastic inverse problems. We used the Kullback-Leibler divergence and the posterior obtained using consistent Bayesian inference to measure the information gain of a design and present a discrete optimization procedure for choosing the optimal experimental design that maximizes the expected information gain. The optimization procedure presented in this paper is limited in terms of the number of observations we can consider, but was chosen to focus attention on the definition and approximation of the expected information gained. More efficient strategies, utilizing gradient-based methods on continuous design spaces, will be pursued in a future work. We discussed a characterization of the space of observed densities needed to compute the expected information gain and a computationally efficient approach for rescaling observed densities to satisfy the requirements of the consistent Bayesian approach. Numerical examples were given to highlight the properties and utility of our approach.

7 Acknowledgments

J.D. Jakeman’s work was supported by DARPA EQUIPS. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525.