Discussion: “Bayesian Optimal Design of Experiments for Inferring the Statistical Expectation of Expensive Black-Box Functions” (Pandita, P., Bilionis, I., and Panchal, J., 2019, ASME J. Mech. Des., 141(10), p. 101404)

In Ref. [1], the authors developed a sequential Bayesian optimal design framework to estimate the statistical expectation of a black-box function $f(x):Rd→R$⁠. Let x ∼ p(x) with p(x) the probability distribution of the input x, the statistical expectation is then defined as follows:The function f(x) is not known a priori but can be evaluated at arbitrary x with Gaussian noise of variance $σ2$⁠:

Based on the Gaussian process surrogate learned from the available samples D_n = {X_n, Y_n}, i.e., $f(x)|Dn∼GP(mn(x),kn(x,x′))$⁠, the next-best sample is chosen by maximizing the information-based acquisition $G(x~)$⁠:where $G(x~)$ computes the information gain of adding a sample $y~$ at $x~$⁠, i.e., the expected KL divergence between the current estimation p(q|D_n) and the hypothetical next-step estimation $p(q|Dn,x~,y~)$⁠:In Eq. (4), $y~$ is chosen based on the surrogate f(x)|D_n following a distribution of $N(mn(x~),kn(x~,x~)+σ2)$⁠); q is considered as a random variable with uncertainties coming from the (current and next-step) surrogates. It is noted that $G(x~)$ also depends on the hyperparameter $θ$ in the learned Gaussian process f(x)|D_n. We neglect this dependence for simplicity, which does not affect the main derivation.

As a major contribution of the discussed paper, the authors simplified the information-based acquisition as Eq. (30) in Ref. [1]:where $σ12$ and $σ22$ are, respectively, the variances of current estimation and hypothetical next-step estimation of q; $v(x~)=∫kn(x~,x)p(x)dx$ and $σn2(x~)=kn(x~,x~)$⁠. Furthermore, for numerical computation of Eq. (5), the authors developed analytical formulas for each involved quantity (important for high-dimensional computation) under uniform distribution of x.

To simplify Eq. (4), we first notice that q|D_n follows a Gaussian distribution with mean μ₁ and variance $σ12$⁠:After adding one hypothetical sample ${x~,y~}$⁠, the function follows an updated surrogate $f(x)|Dn,x~,y~∼GP(mn+1(x),kn+1(x,x′))$ withThe quantity $q|Dn,x~,y~$ can then be represented by another Gaussian with mean μ₂ and variance $σ22$⁠:

We note that Eqs. (7), (8), (12), and (13) are, respectively, intermediate steps of Eqs. (19), (21), (26), and (28) in the discussed paper. Substitute Eq. (6) and Eq. (11) into Eq. (4), one can obtain:where Eq. (14) is exactly Eq. (5) (or Eq. (30) in discussed paper). The fact that the last three terms of Eq. (14) sum up to zero is a direct result of Eq. (13).

The advantage of having a simplified form (15) is that the optimization (3) yields a much more intuitive physical interpretation. Since $σ1$ does not depend on $x~$⁠, Eq. (3) can be reformulated aswhich selects the next-best sample minimizing the expected variance of q. Similar optimization criterion is also used in Refs. [2,3] for the purpose of computing the extreme-event probability.

Another alternative interpretation can be obtained by writing (3) aswhere Eq. (17) is a result of Eq. (13), and $ρy|Dn(x~,x)=kn(x~,x)/((kn(x~,x~)+σ2)(kn(x,x)+σ2))1/2$ is the correlation of y for two inputs $x~$ and x. Equation (18) can be interpreted as to select the next sample which has overall most (weighted) correlation with all x.

We finally remark that the above derivation is for given hyperparamter values $θ$ in f(x)|D_n. This is consistent with the Bayesian approach where the optimal values of $θ$ are chosen from maximizing the likelihood function. However, the discussed paper used a different approach by sampling a distribution of $θ$ and computed $G(x~)$ as an average of the sampling. In the latter case, the above analysis should be likewise considered in a slightly different way, i.e., Eq. (16) should be considered as maximization of the multiplication of $σ22$ from all samples of ${θ(i)}i=1s$⁠:

In the computation of G(x) in the form of Eq. (17), the most heavy computation involved is the integral $∫kn(x~,x)p(x)dx$ (which is prohibitive in high-dimensional problem if direct integration is performed). Following the discussed paper, the integral can be reformulated aswherewith s and Λ involving hyperparameters of the kernel function (with either optimized values from training or selected values as in Ref. [1]).

The main computation is then Eq.(21), for which the authors of the discussed paper addressed the situation of uniform p(x). To generalize the formulation to arbitrary p(x), we can approximate p(x) with the Gaussian mixture model (as an universal approximator of distributions [4]) :Equation (21) can then be formulated as:which yields an analytical computation. In practice, the number of mixtures n_GMM is determined by the complexity of the input distributions, but any distribution of p(x) can be approximated in such a way.

Finally, in computing $G(x~)$ in the form of Eq. (19), computation of $σ12$ in Eq. (8) is necessary. This can also be generalized for arbitrary p(x) using the Gaussian mixture model as follows: