1 Introduction
We consider optimization of composite objective functions, i.e., of the form , where is a blackbox expensivetoevaluate vectorvalued function, and is a realvalued function that can be cheaply evaluated. We assume evaluations are noisefree. These problems arise, for example, in calibration of simulators to realworld data (vrugt2001calibration; cullick2006improved; schultz2018bayesian); in materials and drug design (kapetanovic2008computer; frazier2016bayesian) when seeking to design a compound with a particular set of physical or chemical properties; when finding maximum a posteriori estimators with expensivetoevaluate likelihoods (bliznyuk2008bayesian); and in constrained optimization (gardner14; hernandez_constrained) when seeking to maximize one expensivetoevaluate quantity subject to constraints on others (See Section 2 for a more detailed description of these problems.).
One may ignore the composite structure of the objective and solve such problems using Bayesian optimization (BO) (brochu2010tutorial), which has been shown to perform well compared with other generalpurpose optimization methods for blackbox derivativefree expensivetoevaluate objectives (snoek2012practical). In the standard BO approach, one would build a Gaussian process (GP) prior over based on past observations of , and then choose points at which to evaluate by maximizing an acquisition function computed from the posterior. This approach would not use observations of or knowledge of .
In this paper, we describe a novel BO approach that leverages the structure of composite objectives to optimize them more efficiently. This approach builds a multioutput GP on , and uses the expected improvement (jones1998efficient) under the implied statistical model on as its acquisition function. This implied statistical model is typically nonGaussian when is nonlinear. We refer to the resulting acquisition function as expected improvement for composite functions (EICF) to distinguish it from the classical expected improvement (EI) acquisition function evaluated on a GP posterior on .
Intuitively, the above approach can substantially outperform standard BO when observations of provide information relevant to optimization that is not available from observations of alone. As one example, suppose and are both onedimensional and . If is continuous, , and , then our approach knows that there is a global minimum in the interval , while the standard approach does not. This informational benefit is compounded further when is vectorvalued.
While EICF
is simply the expected improvement under a different statistical model, unlike the classical EI acquisition function, it lacks a closedform analytic expression and must be evaluated through simulation. We provide a simulationbased method for computing unbiased estimators of the gradient of the
EICF acquisition function, which we use within multistart stochastic gradient ascent to allow efficient maximization. We also show that optimizing using EICF is asymptotically consistent under suitable regularity conditions, in the sense that the best point found converges to the global maximum of as the number of samples grows to infinity.In numerical experiments comparing with standard BO benchmarks, EICF provides immediate regret that is several orders of magnitude smaller, and reaches their final solution quality using less than 1/4 the function evaluations.
2 Related Work
2.1 Related Methodological Literature
We work within the Bayesian optimization framework, whose origins date back to the seminal work of movckus1975bayesian
, and which has recently become popular due to its success in hyperparameter optimization of machine learning algorithms
(snoek2012practical; swersky2013multi).Optimizing composite functions has been studied in first and secondorder optimization (shapiro2003class; drusvyatskiy2016efficiency). This literature differs from our paper in that it assumes derivatives are available, and also often assumes convexity and that evaluations are inexpensive. In this setting, leveraging the structure of the objective has been found to improve performance, just as we find here in the setting of derivativefree optimization. However, to the best of our knowledge, ours is the first paper to study composite objective functions within the BO framework and also the first within the more general context of optimization of blackbox derivativefree expensivetoevaluate functions.
Our work is related to marque2017efficient
, which proposes a framework for efficient sequential experimental design for GPbased modeling of nested computer codes. In contrast with our work, that work’s goal is not to optimize a composite function, but instead to learn it as accurately as possible within a limited evaluation budget. A predictive variance minimization sampling policy is proposed and methods for efficient computation are provided. Moreover, it is assumed that both the inner (
) and outer () functions are realvalued and expensivetoevaluate blackbox functions, while our method uses the easeofevaluation of the outer function for substantial benefit.Our work is also similar in spirit to overstall2013strategy
, which proposes to model an expensivetoevaluate vector of parameters of a posterior probability density function using a multioutput GP instead of the function directly using a singleoutput GP. The surrogate model is then used to perform Bayesian inference.
Constrained optimization is a special case of optimization of a composite objective. To see this, take to be the objective to be maximized and take , for , to be the constraints, all of which are constrained to be nonnegative without loss of generality. Then, we recover the original constrained optimization problem by setting
Moreover, when specialized to this particular setting, our EICF acquisition function is equivalent to the expected improvement for constrained optimization as proposed by schonlau1998global and gardner14.
Within the constrained BO literature, our work also shares several methodological similarities with picheny2016bayesian, which considers an augmented Lagrangian and models its components as GPs instead of it directly as a GP. As in our work, the expected improvement under this statistical model is used as acquisition function. Moreover, it is shown that this approach outperforms the standard BO approach.
Our method for optimizing the EICF acquisition function uses an unbiased estimator of the gradient of EICF within a multistart stochastic gradient ascent framework. This technique is structurally similar to methods developed for optimizing acquisition functions in other BO settings without composite objectives, including the parallel expected improvement (wang2016parallel) and the parallel knowledgegradient (wu2016parallel).
2.2 Related Application Literature
Optimization of composite blackbox derivativefree expensivetoevaluate functions arises in a number of application settings in the literature, though this literature does not leverage the composite structure of the objective to optimize it more efficiently.
In materials design, it arises when the objective is the combination of multiple material properties via a performance index (ashby1993materials), e.g., the specific stiffness, which is the ratio of Young’s modulus and the density, or normalization (jahan2015state). Here, is the set of material properties that results from a particular chemical composition or set of processing conditions, , and is given by the performance index or normalization method used. Evaluating the material properties, , that result from a materials design typically requires running expensive physical or computational experiments that do not provide derivative information, for which BO is appropriate (kapetanovic2008computer; ueno2016combo; ju2017designing; griffiths2017constrained).
Optimization of composite functions also arises in calibration of expensive blackbox simulators (vrugt2001calibration; cullick2006improved; schultz2018bayesian), where the goal is to find input parameters, , to the simulator such that its vectorvalued output, , most closely matches a vector data observed in the real world, . Here, the objective to be minimized is , where is often the norm, norm, or some monotonic transformation of the likelihood of observation errors.
If one has a prior probability density
on , and the loglikelihood of some realworld observation error, , is proportional to(as it would be, for example, with independent normally distributed errors taking
to be the norm), then, finding the maximum a posteriori estimator of (bliznyuk2008bayesian) is an optimization problem with a composite objective: the logposterior is equal to the sum of a constant and (In this example, is actually a function of both and . Our framework extends easily to this setting as long as remains a cheaptoevaluate function.).3 Problem Description and Standard Approach
As described above, we consider optimization of objectives of the form , where is a blackbox expensivetoevaluate continuous function whose evaluations do not provide derivatives, is a function that can be cheaply evaluated, and . As is common in BO, we assume that is not too large () and that projections onto can be efficiently computed. We also place the technical condition that , where is an variate standard normal random vector. The problem to be solved is
(1) 
As discussed before, one can solve problem (1) by applying standard BO to the objective function, . This approach models
as drawn from a GP prior probability distribution. Then, iteratively, indexed by
, this approach would choose the point at which to evaluate next by optimizing an acquisition function, such as the EI acquisition function (movckus1975bayesian; jones1998efficient). This acquisition function would be calculated from the posterior distribution given , which is itself a GP, and would quantify the value of an evaluation at a particular point. Although would be observed as part of this standard approach, these evaluations would be ignored when calculating the posterior distribution and acquisition function.4 Our Approach
We now describe our approach, which like the standard BO approach is comprised of a statistical model and an acquisition function. Unlike standard BO, however, our approach leverages the additional information in evaluations of , along with knowledge of . We argue below and demonstrate in our numerical experiments that this additional information can substantially reduce the number of evaluations required to find good approximate global optima.
Briefly, our statistical model is a multioutput Gaussian process on (alvarez2012kernels) (Section 4.1), and our acquisition function, EICF, is the expected improvement under this statistical model (Section 4.2). This acquisition function, unfortunately, cannot be computed in closed form for most functions . In Section 4.3, under mild regularity conditions, we provide a technique for efficiently maximizing EICF. We also provide a theoretical analysis showing that EICF is asymptotically consistent (Section 4.4). Finally, we conclude this section by discussing the computational complexity of our approach (Section 4.5).
4.1 Statistical Model
We model as drawn from a multioutput GP distribution (alvarez2012kernels), , where is the mean function, is the covariance function, and is the cone of positive definite matrices. Analogously to the singleoutput case, after observing evaluations of , , the posterior distribution on is again a multioutput GP, , where and can be computed in closed form in terms of and (liu2018remarks).
4.2 Expected Improvement for Composite Functions
We define the expected improvement for composite functions analogously to the classical expected improvement, but where our posterior on is given by the composition of and the normally distributed posterior distribution on :
(2) 
where is the maximum value across the points that have been evaluated so far, , indicates the conditional expectation given the available observations at time , , and is the positive part function.
When is scalarvalued and is the identity function, is equivalent to the classical expected improvement computed directly from a GP prior on , and has an analytic expression that makes it easy to compute and optimize. For general nonlinear functions , however, cannot be computed in closed form. Despite this, as we shall see next, under mild regularity conditions, can be efficiently maximized.
Figure 6 illustrates the EICF and classical EI acquisition functions in a setting where is scalarvalued, , we have evaluated and at four points, and we wish to minimize . The righthand column shows the posterior distribution on
and EI acquisition function using the standard approach: posterior credible intervals have 0 width at points where we have evaluated (since evaluations are free from noise), and become wider as we move away from them. The classical expected improvement is largest near the right limit of the domain, where the posterior mean computed using observations of
alone is relatively small and has large variance.The lefthand column shows the posterior distribution on , computed using a GP (singleoutput in this case, since is scalarvalued), the resulting posterior distribution on , and the resulting EICF acquisition function. The posterior distribution on (which is not normally distributed) has support only on nonnegative values, and places higher probability on small values of in the regions , which creates a larger value for EICF in these regions.
Examining past observations of , the points with high EICF () seem substantially more valuable to evaluate than the point with largest EI (). Indeed, for the region , we know that is below near the left limit, and is above near the right limit. The continuity of implies that is 0 at some point in this region, which in turn implies that has a global optimum in this region. Similarly, is also quite likely to have a global optimum in . EICF takes advantage of this knowledge in its sampling decisions while classical EI does not.
4.3 Computation and Maximization of EiCf
We now describe computation and efficient maximization of EICF. For any fixed , the time posterior distribution on is multivariate normal. (By the “time posterior distribution”, we mean the conditional distribution given .) We let denote its mean vector and denote its covariance matrix. We also let denote the lower Cholesky factor of . Then, we can express , where is a variate standard normal random vector under the time posterior distribution, and thus
(3) 
Thus, we can compute via Monte Carlo, as summarized in Algorithm 1. We note that (3) and the condition imply that is finite for all .
In principle, the above is enough to maximize using a derivativefree global optimization algorithm (for nonexpensive noisy functions). However, such methods typically require a large number of samples, and optimization can be typically performed with much greater efficiency if derivative information is available (jamieson2012query; swisher2000survey). The following proposition describes a simulationbased procedure for generating such derivative information. A formal statement and proof can be found in Appendix A.
Proposition 1.
Under mild regularity conditions, is differentiable almost everywhere, and its gradient, when it exists, is given by
(4) 
where
(5) 
Thus, provides an unbiased estimator of . To construct such an estimator, we would draw an independent standard normal random vector and then compute using (5), optionally averaging across multiple samples as in Algorithm 1. To optimize , we then use this gradient estimation procedure within stochastic gradient ascent, using multiple restarts. The final iterate from each restart is an approximate stationary point of the . We then use Algorithm 1 to select the stationary point with the best value of .
4.4 Theoretical Analysis
Here we present two results giving insight into the properties of the expected improvement for composite functions. The first result simply states that, when is linear, EICF has a closed form analogous to the one of the classical EI.
Proposition 2.
Suppose that is given by for some fixed . Then,
where , , and and
are the standard normal probability density function and cumulative distribution function, respectively.
This result can be easily verified by noting that, since the time posterior distribution of is variate normal with mean vector and covariance matrix , the time posterior distribution of is normal with mean and variance . Proposition 2 does not, however, mean that our approach is equivalent to the classical one when g is linear. This is because, in general, the posterior distribution given observations of is different from the one given observations of . We refer the reader to Appendix B for a discussion.
Our second result states that, under suitable conditions, our acquisition function is asymptotically consistent, i.e., the solution found by our method converges to the global optimum when the number of evaluations goes to infinity. An analogous result for the classical expected improvement was proved by vazquez2010convergence.
Theorem 1.
Let be the sequence of evaluated points and suppose there exists such that for all ,
Then, under suitable regularity conditions and as ,
4.5 Computational Complexity of Posterior Inference
The computation required to maximize the classical EI acquisition function is dominated by the computation of the posterior mean and variance and thus in principle scales as (with a precomputation of complexity ) with respect to the number of evaluations (shahriari2016taking). However, recent advances on approximate fast GP training and prediction may considerably reduce the computational burden (pleiss2018constant).
In our approach, the computational cost is again dominated by the computation of the posterior mean and covariance matrix, and , respectively. When the outputs of are modeled independently, the components of and can be computed separately ( is diagonal in this case) and thus computation of the posterior mean and covariance scales as . This allows our approach to be used even if has a relatively large number of outputs. However, in general, if correlation between components of is modeled, these computations scale as . Therefore, in principle there is a tradeoff between modeling correlation between components of , which presumably allows for a faster learning of , and retaining tractability in the computation of the acquisition function.
5 Numerical Experiments
We compare the performance of three acquisition functions: expected improvement (EI), probability of improvement (PI) (kushner1964new), and the acquisition function that chooses points uniformly at random (Random), both under our proposed statistical model and the standard one, i.e., modeling using a multioutput GP and modeling directly using a singleoutput GP, respectively. We refer the reader to Appendix C for a formal definition of the probability of improvement under our statistical model, and a discussion of how we maximize this acquisition function in our numerical experiments. To distinguish each acquisition function under our proposed statistical model from its standard version, we append ”CF” to its abbreviation; so if the classical expected improvement acquisition function is denoted EI, then the expected improvement under our statistical model is denoted EICF, as previously defined. We also include as a benchmark the predictive entropy search (PES) acquisition function (hernandez2014predictive) under the standard statistical model, i.e., modeling directly using a singleoutput GP. For all problems and methods, an initial stage of evaluations is performed using points chosen uniformly at random over .
For EICF, PICF, and RandomCF, the outputs of are modeled using independent GP prior distributions. All GP distributions involved, including those used by the standard BO methods (EI, PI, Random, and PES), have a constant mean function and ARD squared exponential covariance function; the associated hyperparameters are estimated under a Bayesian approach. As proposed in snoek2012practical, for all methods we use an averaged version of the acquisition function, obtained by first drawing 10 samples of the GP hyperparameters, computing the acquisition function conditioned on each of these hyperparameters, and then averaging the results.
We calculate each method’s simple regret at the point it believes to be the best based on evaluations observed thus far. We take this point to be the point with the largest (or smallest, if minimizing) posterior mean. For EICF, PICF, and RandomCF, we use the posterior mean on implied by the GP posterior on , and for EI, PI, Random, and PES we use the GP posterior on
. Error bars in the plots below show the mean of the base10 logarithm of the simple regret plus and minus 1.96 times the standard deviation divided by the square root of the number of replications. Each experiment was replicated 100 times.
Our code is available at bocf.
5.1 GPGenerated Test Problems
The first two problems used functions generated at random from GPs. Each component of was generated by sampling on a uniform grid from independent GP distributions with different fixed hyperparameters and then taking the resulting posterior mean as a proxy; the hyperparameters were not known to any of the algorithms. The details of each problem, including the function used, are summarized in Table 1.
Results are shown on a logarithmic scale in Figures 7 and 8, where the horizontal axis indicates the number of samples following the initial stage. EICF outperforms the other methods significantly. Regret is smaller than the best of the standard BO benchmarks throughout and by several orders of magnitude after 50 evaluations (5 orders of magnitude smaller in test 1, and 2 in test 2). It also requires substantially fewer evaluations beyond the first stage to reach the regret achieved by the best of the standard BO benchmarks in 100 evaluations: approximately 30 in test 1, and 10 in test 2. RandomCF performs surprisingly well in type2 GPgenerated problems, suggesting that a substantial part of the benefit provided by our approach is the value of the additional information available from observing . In type1 problems it does not perform as well, suggesting that highquality decisions about where to sample are also important.
Problem  

1  5  
2  4 
5.2 Standard Global Optimization Test Problems
We assess our approach’s performance on two standard benchmark functions from the global optimization literature: the Langermann (langermann) and Rosenbrock (rosenbrock) functions. We refer the reader to Appendix D for a description of how these functions are adapted to our setting.
Results of applying our method and benchmarks to these problems are shown on a logarithmic scale in Figures 9 and 10. As before, EICF outperforms competing methods with respect to the final achieved regret. PICF and RandomCF also perform well compared to methods other than EICF. Moreover, in the Langermann test problem, PICF outperforms EICF during the first 20 evaluations.
5.3 Environmental Model Function
The environmental model function was originally proposed by bliznyuk2008bayesian and is now a wellknown test problem in the literature of Bayesian calibration of expensive computer models. It models a chemical accident that has caused a pollutant to spill at two locations into a long and narrow holding channel, and is based on a firstorder approach to modeling the concentration of substances in such channels under the assumption that the channel can be approximated by an infinitely long onedimensional system with diffusion as the only method of transport. This leads to the concentration representation
where is the mass of pollutant spilled at each location, is the diffusion rate in the channel, is the location of the second spill, and is the time of the second spill.
We observe in a grid of values; specifically, we observe , where , , and are the underlying true values of these parameters. Since we assume noiseless observations, the calibration problem reduces to finding so that the observations at the grid minimize the sum of squared errors, i.e., our goal is to minimize
In our experiment, we take , , and . The search domain is , , and .
Results from this experiment are shown in Figure 11. As above, EICF performs best, with PICF and RandomCF also significantly outperforming benchmarks that do not leverage the composite structure.
6 Conclusion and Future Work
We have proposed a novel Bayesian optimization approach for objective functions of the form , where is a blackbox expensivetoevaluate vectorvalued function, and is a realvalued function that can be cheaply evaluated. Our numerical experiments show that this approach may substantially outperform standard Bayesian optimization, while retaining computational tractability.
There are several relevant directions for future work. Perhaps the most evident is to understand whether other wellknown acquisition functions can be generalized to our setting in a computationally tractable way. We believe this to be true for predictive entropy search (hernandez2014predictive) and knowledge gradient (scott2011correlated). Importantly, these acquisition functions would allow noisy and decoupled evaluations of the components of , thus increasing the applicability of our approach. However, in the standard Bayesian optimization setting, they are already computationally intensive and thus a careful analysis is required to make them computationally tractable in our setting.
Acknowledgements
This work was partially supported by NSF CAREER CMMI1254298, NSF CMMI1536895 and AFOSR FA95501510038. The authors also thank Eytan Bakshy for helpful comments.
References
Appendix A Unbiased Estimator of the Gradient of EiCf
In this section we prove that, under mild regularity conditions, is differentiable and an unbiased estimator of its gradient can be efficiently computed. More concretely, we prove the following.
Proposition A.1.
Suppose that is differentiable and let be an open subset of so that and are differentiable on and there exists a measurable function satisfying

for all ,

, where is a variate standard normal random vector.
Further, suppose that for almost every the set is countable. Then, is differentiable on and its gradient is given by
where the expectation is with respect to and
Proof.
Since is differentiable and and are differentiable on , for any fixed the function is differentiable on as well. This in turn implies that the function is continuous on and differentiable at every such that , with gradient equal to . From our assumption that for almost every the set is countable, it follows that for almost every the function is continuous on and differentiable on all , except maybe on a countable subset. Using this, along with conditions 1 and 2, and Theorem 1 in l1990unified, the desired result follows. ∎
We end this section by making a few remarks.

If and are differentiable on , then one can show that and are differentiable on .

If one imposes the stronger condition , then
has finite second moment, and thus this unbiased estimator of
can be used within stochastic gradient ascent to find a stationary point of (bottou1998online). 
In Proposition A.1, the condition that for almost every the set is countable, can be weakened to the following more technical condition: for almost every , every and every , there exists such that the set is countable, where denotes the th canonical vector in .
Appendix B EICF and EI Do Not Coincide When Is Linear
Recall the following result that was stated in the main paper.
Proposition B.1.
Suppose that is given by for some fixed . Then,
The resemblance of the above expression to the classical EI acquisition function may make one think that, in the above case, EICF coincides, in some sense, with the classical EI under an appropriate choice of the prior distribution.
Indeed, suppose that we set a singleoutput GP prior with mean and covariance function of (and fix its hyperparameters), then
However, if we condition on rather than in the lefthand side, then the equality is no longer true, even if the values on which we condition satisfy :
Thus, even if we initiate optimization using EICF and a parallel optimization using EI with a singleoutput Gaussian process as described above, their acquisition functions will cease to agree once we condition on the results of an evaluation.
Appendix C Probability of Improvement for Composite Functions
In this section, we formally define the probability of improvement for composite functions (PICF) acquisition function and specify its implementation details used within our experimental setup.
Analogously to EICF, PICF is simply defined as the probability of improvement evaluated with respect to the implied posterior distribution on when we model as a multioutput GP:
where denotes the conditional probability given the available observations at time , , and is a parameter to be specified. As we did with EICF, we can express as
where is a variate standard normal random vector under the time posterior distribution.
We can further rewrite using an indicator function as
which implies that PICF can be computed with arbitrary precision following a Monte Carlo approach as well:
where are draws of an variate standard normal random vector. However, an unbiased estimator of the gradient of PICF cannot be computed following an analogous approach to the one used with EICF. In fact, at those points for which the function is differentiable. Thus, even if exists, in general
unless .
In our experiments, we adopt a sample average approximation (kim2015guide) scheme for approximately maximizing PICF. At each iteration we fix and choose the next point to evaluate as
where and . We solve the above optimization problem using the derivativefree optimization algorithm, CMAES (hansen2016cma).
Appendix D Description of Langermann and Rosenbrock Test Problems
The following pair of test problems are standard benchmarks in the global optimization literature. In this section, we describe in detail how they are adapted our setting, i.e., how we express them as composite functions.
d.1 Langermann Function
The Langermann function (langermann_supp) is defined by where
In our experiment we set , , ,
and .
d.2 Rosenbrock Function
The Rosenbrock function (rosenbrock_supp) is
We adapt this problem to our framework by taking and defining and by
Appendix E Asymptotic Consistency of the Expected Improvement for Composite Functions
e.1 Basic Definitions and Assumptions
In this section we prove that, under suitable conditions, the expected improvement sampling policy is asymptotically consistent in our setting. In the standard Bayesian optimization setting, this was first proved under quite general conditions by vazquez2010convergence_supp. Later, bull2011convergence provided convergence rates for several expected improvementtype policies both with fixed hyperparameters and hyperparameters estimated from the data in suitable way. Here, we restrict to prove asymptotic consistency, under fixed hyperparameters, following a similar approach to vazquez2010convergence_supp. In particular, we provide a generalization of the NoEmptyBall (NEB) condition, under which the expected improvement sampling policy is guaranteed to be asymptotically consistent in our setting. In the reminder of this work denotes the sequence of points at which is evaluated, which is not necessarily given by the expected improvement acquisition function, unless explicitly stated.
Definition E.1 (GeneralizedNoEmptyBall property).
We shall say that a kernel, , satisfies the GeneralizedNoEmptyBall (GNEB) property if, for all sequences in and all , the following assertions are equivalent:

is a limit point of .

There exists a subsequence of converging to a singular matrix.
We highlight that, if is diagonal, i.e. if the output components are independent of each other, the GNEB property holds provided that at least one of its components satisfies the standard NEB property. In particular, the following result is a corollary of Proposition 10 in vazquez2010convergence_supp.
Corollary E.2.
Suppose is diagonal and at least one of its components has a spectral density whose inverse has at most polynomial growth. Then, satisfies the GNEB property.
Thus, the GNEB property holds, in particular, if is diagonal and at least one of its components is a Matérn kernel (stein2012interpolation).
Now we introduce some additional notation. We denote by to the the reproducing kernel Hilbert space associated with (alvarez2012kernels_supp). As is standard in Bayesian optimization, we make a slight abuse of notation and denote by both a fixed element of and a random function distributed according to a Gaussian process with mean and kernel (below we assume ); we shall explicitly state whenever is held fixed. As before, we denote by . Finally, we make the following standing assumptions.

is a compact subset of , for some .

The prior mean function is identically .

is continuous, positive definite, and satisfies the GNEB property.

is continuous.

For any bounded sequences and , , where the expectation is over and is a dimensional standard normal random vector.
The assumption that is continuous guarantees that is continuous, provided that is continuous as well. Moreover, in this case, since is compact, attains its maximum value in ; we shall denote this maximum value by , i.e., .
e.2 Preliminary Results
Before proving asymptotic consistency, we prove several auxiliary results. We begin by proving that is continuous.
Proposition E.3.
For any , the function defined by
where the expectation is over and is a dimensional standard normal random vector, is continuous.
Proof.
Let be a convergent sequence with limit . Since is continuous, and are both continuous functions of , and thus and as . Moreover, since is continuous too, it follows by the continuous mapping theorem (billingsley2013convergence) that
almost surely as .
Now observe that
Moreover, the sequences and are both convergent (with finite limits) and thus are bounded. Hence, the above inequality, along with assumption 5 and the dominated convergence theorem (williams1991), imply that
as , i.e., . Hence, is continuous. ∎
Lemma E.4.
Let and be two sequences in . Assume that is convergent, and denote by its limit. Then, each of the following conditions implies the next one:

is a limit point of .

as .

For any fixed , as .
Proof.
First we prove that 1 implies 2. If is an element of , say , then, for , we have
where we use Lemma F.3 and the fact that is continuous. Now assume is not an element of . Let be a subsequence of converging to and let . Then, by Lemmas F.1 and F.2 we obtain
Finally, since is not an element of , , and it follows from the continuity of that
and thus .
Now we prove that 2 implies 3. Using the CauchySchwarz inequality in , we obtain
Thus,
since is continuous. ∎
Lemma E.5.
Let . Then, for all , .
Proof.
Fix and let be the sequence of points generated by the expected improvement policy, i.e., . Let be a limit point of and let be any subsequence converging to . Consider the sequence given by for all , . Clearly, , and thus Lemma E.4 implies that and . In particular, and , i.e., and . Moreover, is a bounded increasing sequence, and thus has a finite limit, , which satisfies as is a limit point of and is continuous.
The sequences and are convergent and thus bounded. Hence, from assumption 5 and the dominated convergence theorem we obtain that
but
and thus the desired conclusion follows. ∎
e.3 Proof of the Main Result
We are now in position to prove that the expected improvement acquisition function is asymptotically consistent in the composite functions setting.
Theorem E.6 (Asymptotic consistency of EICF).
Assume that the covariance function, , satisfies the GNEB property. Then, for any fixed and , any (measurable) sequence