# Bayesian Optimization of Composite Functions

We consider optimization of composite objective functions, i.e., of the form f(x)=g(h(x)), where h is a black-box derivative-free expensive-to-evaluate function with vector-valued outputs, and g is a cheap-to-evaluate real-valued function. While these problems can be solved with standard Bayesian optimization, we propose a novel approach that exploits the composite structure of the objective function to substantially improve sampling efficiency. Our approach models h using a multi-output Gaussian process and chooses where to sample using the expected improvement evaluated on the implied non-Gaussian posterior on f, which we call expected improvement for composite functions (). Although cannot be computed in closed form, we provide a novel stochastic gradient estimator that allows its efficient maximization. We also show that our approach is asymptotically consistent, i.e., that it recovers a globally optimal solution as sampling effort grows to infinity, generalizing previous convergence results for classical expected improvement. Numerical experiments show that our approach dramatically outperforms standard Bayesian optimization benchmarks, reducing simple regret by several orders of magnitude.

Comments

There are no comments yet.

## Authors

• 2 publications
• 20 publications
• ### A hierarchical expected improvement method for Bayesian optimization

Expected improvement (EI) is one of the most popular Bayesian optimizati...
11/17/2019 ∙ by Zhehui Chen, et al. ∙ 0

read it

• ### A Tutorial on Bayesian Optimization

Bayesian optimization is an approach to optimizing objective functions t...
07/08/2018 ∙ by Peter I. Frazier, et al. ∙ 12

read it

• ### Stratified Bayesian Optimization

We consider derivative-free black-box global optimization of expensive n...
02/07/2016 ∙ by Saul Toscano-Palmerin, et al. ∙ 0

read it

• ### Efficient Nonmyopic Bayesian Optimization via One-Shot Multi-Step Trees

Bayesian optimization is a sequential decision making framework for opti...
06/29/2020 ∙ by Shali Jiang, et al. ∙ 0

read it

• ### Efficient batch-sequential Bayesian optimization with moments of truncated Gaussian vectors

We deal with the efficient parallelization of Bayesian global optimizati...
09/09/2016 ∙ by Sébastien Marmin, et al. ∙ 0

read it

• ### Model Guided Sampling Optimization for Low-dimensional Problems

Optimization of very expensive black-box functions requires utilization ...
08/31/2015 ∙ by Lukas Bajer, et al. ∙ 0

read it

• ### Integrals over Gaussians under Linear Domain Constraints

Integrals of linearly constrained multivariate Gaussian densities are a ...
10/21/2019 ∙ by Alexandra Gessner, et al. ∙ 0

read it

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

We consider optimization of composite objective functions, i.e., of the form , where is a black-box expensive-to-evaluate vector-valued function, and is a real-valued function that can be cheaply evaluated. We assume evaluations are noise-free. These problems arise, for example, in calibration of simulators to real-world 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 expensive-to-evaluate likelihoods (bliznyuk2008bayesian); and in constrained optimization (gardner14; hernandez_constrained) when seeking to maximize one expensive-to-evaluate 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 general-purpose optimization methods for black-box derivative-free expensive-to-evaluate 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 multi-output GP on , and uses the expected improvement (jones1998efficient) under the implied statistical model on as its acquisition function. This implied statistical model is typically non-Gaussian when is non-linear. We refer to the resulting acquisition function as expected improvement for composite functions (EI-CF) 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 one-dimensional 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 vector-valued.

While EI-CF

is simply the expected improvement under a different statistical model, unlike the classical EI acquisition function, it lacks a closed-form analytic expression and must be evaluated through simulation. We provide a simulation-based method for computing unbiased estimators of the gradient of the

EI-CF acquisition function, which we use within multi-start stochastic gradient ascent to allow efficient maximization. We also show that optimizing using EI-CF 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, EI-CF 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 second-order 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 derivative-free 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 black-box derivative-free expensive-to-evaluate functions.

Our work is related to marque2017efficient

, which proposes a framework for efficient sequential experimental design for GP-based 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 real-valued and expensive-to-evaluate black-box functions, while our method uses the ease-of-evaluation of the outer function for substantial benefit.

Our work is also similar in spirit to overstall2013strategy

, which proposes to model an expensive-to-evaluate vector of parameters of a posterior probability density function using a multi-output GP instead of the function directly using a single-output 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 non-negative without loss of generality. Then, we recover the original constrained optimization problem by setting

 g(y)={y1, if yi≥0 for all i>1,−∞, otherwise.

Moreover, when specialized to this particular setting, our EI-CF 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 EI-CF acquisition function uses an unbiased estimator of the gradient of EI-CF 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 knowledge-gradient (wu2016parallel).

### 2.2 Related Application Literature

Optimization of composite black-box derivative-free expensive-to-evaluate 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 black-box simulators (vrugt2001calibration; cullick2006improved; schultz2018bayesian), where the goal is to find input parameters, , to the simulator such that its vector-valued 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 log-likelihood of some real-world 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 log-posterior 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 cheap-to-evaluate function.).

## 3 Problem Description and Standard Approach

As described above, we consider optimization of objectives of the form , where is a black-box expensive-to-evaluate 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

 maxx∈Xg(h(x)). (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 multi-output Gaussian process on (alvarez2012kernels) (Section 4.1), and our acquisition function, EI-CF, 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 EI-CF. We also provide a theoretical analysis showing that EI-CF 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 multi-output GP distribution (alvarez2012kernels), , where is the mean function, is the covariance function, and is the cone of positive definite matrices. Analogously to the single-output case, after observing evaluations of , , the posterior distribution on is again a multi-output 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 :

 EI-CFn(x)=En[{g(h(x))−f∗n}+], (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 scalar-valued 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 EI-CF and classical EI acquisition functions in a setting where is scalar-valued, , we have evaluated and at four points, and we wish to minimize . The right-hand 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 left-hand column shows the posterior distribution on , computed using a GP (single-output in this case, since is scalar-valued), the resulting posterior distribution on , and the resulting EI-CF acquisition function. The posterior distribution on (which is not normally distributed) has support only on non-negative values, and places higher probability on small values of in the regions , which creates a larger value for EI-CF in these regions.

Examining past observations of , the points with high EI-CF () 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 . EI-CF takes advantage of this knowledge in its sampling decisions while classical EI does not.

### 4.3 Computation and Maximization of Ei-Cf

We now describe computation and efficient maximization of EI-CF. 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

 EI-CFn(x)=En[{g(μn(x)+Cn(x)Z)−f∗n}+]. (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 derivative-free global optimization algorithm (for non-expensive 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 simulation-based 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

 ∇EI-CFn(x)=En[γn(x,Z)], (4)

where

 γn(x,Z)={0, if g(μn(x)+Cn(x)Z)≤f∗n,∇g(μn(x)+Cn(x)Z), otherwise. (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, EI-CF has a closed form analogous to the one of the classical EI.

###### Proposition 2.

Suppose that is given by for some fixed . Then,

 EI-CFn(x)=Δn(x)Φ(Δn(x)σn(x))+σn(x)φ(Δn(x)σn(x)),

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 ,

 xn+1∈argmaxx∈XEI-CFn(x).

Then, under suitable regularity conditions and as ,

 f∗n→maxx∈Xf(x).

A formal statement and proof of Theorem 1 can be found in Appendix E.

### 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 pre-computation 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 trade-off 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 multi-output GP and modeling directly using a single-output 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 EI-CF, 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 single-output GP. For all problems and methods, an initial stage of evaluations is performed using points chosen uniformly at random over .

For EI-CF, PI-CF, and Random-CF, 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 EI-CF, PI-CF, and Random-CF, 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 base-10 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 GP-Generated 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. EI-CF 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. Random-CF performs surprisingly well in type-2 GP-generated problems, suggesting that a substantial part of the benefit provided by our approach is the value of the additional information available from observing . In type-1 problems it does not perform as well, suggesting that high-quality decisions about where to sample are also important.

### 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, EI-CF outperforms competing methods with respect to the final achieved regret. PI-CF and Random-CF also perform well compared to methods other than EI-CF. Moreover, in the Langermann test problem, PI-CF outperforms EI-CF during the first 20 evaluations.

### 5.3 Environmental Model Function

The environmental model function was originally proposed by bliznyuk2008bayesian and is now a well-known 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 first-order approach to modeling the concentration of substances in such channels under the assumption that the channel can be approximated by an infinitely long one-dimensional system with diffusion as the only method of transport. This leads to the concentration representation

 c(s,t;M,D,L,τ)= M√4πDtexp(−s24Dt)+ I{t>τ}M√4πD(t−τ)exp(−(s−L)24D(t−τ)),

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

 ∑(s,t)∈S×T(c(s,t;M0,D0,L0,τ0)−c(s,t;M,D,L,τ))2.

In our experiment, we take , , and . The search domain is , , and .

Results from this experiment are shown in Figure 11. As above, EI-CF performs best, with PI-CF and Random-CF 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 black-box expensive-to-evaluate vector-valued function, and is a real-valued 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 well-known 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 CMMI-1254298, NSF CMMI-1536895 and AFOSR FA9550-15-1-0038. The authors also thank Eytan Bakshy for helpful comments.

## Appendix A Unbiased Estimator of the Gradient of Ei-Cf

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

1. for all ,

2. , 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

 ∇EI-CFn(x)=En[γn(x,Z)],

where the expectation is with respect to and

 γn(x,z)={∇g(μn(x)+Cn(x)z), if g(μn(x)+Cn(x)z)>f∗n,0, otherwise.
###### 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 EI-CF and EI Do Not Coincide When g 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,

 EI-CFn(x)=Δn(x)Φ(Δn(x)σn(x))+σn(x)φ(Δn(x)σn(x))

The resemblance of the above expression to the classical EI acquisition function may make one think that, in the above case, EI-CF coincides, in some sense, with the classical EI under an appropriate choice of the prior distribution.

Indeed, suppose that we set a single-output GP prior with mean and covariance function of (and fix its hyperparameters), then

 E[{w⊤h(x)−f∗n}+∣xi,w⊤h(xi)=yi:i=1,…n]=E[{f(x)−f∗n}+∣xi,f(xi)=yi:i=1,…n].

However, if we condition on rather than in the left-hand side, then the equality is no longer true, even if the values on which we condition satisfy :

Thus, even if we initiate optimization using EI-CF and a parallel optimization using EI with a single-output 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 (PI-CF) acquisition function and specify its implementation details used within our experimental setup.

Analogously to EI-CF, PI-CF is simply defined as the probability of improvement evaluated with respect to the implied posterior distribution on when we model as a multi-output GP:

 PI-CF(x)=Pn(g(h(x))≥f∗n+δ),

where denotes the conditional probability given the available observations at time , , and is a parameter to be specified. As we did with EI-CF, we can express as

 PI-CF(x)=Pn(g(μn(x)+Cn(x)Z)≥f∗n+δ),

where is a -variate standard normal random vector under the time- posterior distribution.

We can further rewrite using an indicator function as

 PI-CF(x)=En[I{g(μn(x)+Cn(x)Z)≥f∗n+δ}],

which implies that PI-CF can be computed with arbitrary precision following a Monte Carlo approach as well:

 PI-CF(x)≈1LL∑ℓ=1I{g(μn(x)+Cn(x)Z(ℓ))≥f∗n+δ},

where are draws of an -variate standard normal random vector. However, an unbiased estimator of the gradient of PI-CF cannot be computed following an analogous approach to the one used with EI-CF. In fact, at those points for which the function is differentiable. Thus, even if exists, in general

 ∇PI-CF(x)≠En[∇I{g(μn(x)+Cn(x)Z)≥f∗n+δ}],

unless .

In our experiments, we adopt a sample average approximation (kim2015guide) scheme for approximately maximizing PI-CF. At each iteration we fix and choose the next point to evaluate as

 xn+1∈argmaxx∈X1LL∑ℓ=1I{g(μn(x)+Cn(x)Z(ℓ))≥f∗n+δ},

where and . We solve the above optimization problem using the derivative-free optimization algorithm, CMA-ES (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

 hj(x) =d∑i=1(xi−Aij)2, j=1,…m, g(y) =−m∑j=1cjexp(−yj/π)cos(πyj).

In our experiment we set , , ,

 A=(3521752149),

and .

### d.2 Rosenbrock Function

The Rosenbrock function (rosenbrock_supp) is

 f(x)=−d−1∑j=1100(xj+1−x2j)2+(xj−1)2

We adapt this problem to our framework by taking and defining and by

 hj(x) =xj+1−x2j, j=1,…,4, hj+4(x) =xj, j=1,…,4, g(y) =−4∑j=1100y2j+(yj+4−1)2.

## 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 improvement-type 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 No-Empty-Ball (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 (Generalized-No-Empty-Ball property).

We shall say that a kernel, , satisfies the Generalized-No-Empty-Ball (GNEB) property if, for all sequences in and all , the following assertions are equivalent:

1. is a limit point of .

2. 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.

1. is a compact subset of , for some .

2. The prior mean function is identically .

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

4. is continuous.

5. 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

 EI-CFn(x)=E[{g(μn(x)+Cn(x)Z)−f∗n}+],

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

 {g(μn(x′k)+Cn(x′k)Z)−f∗n}+→{g(μn(x′∞)+Cn(x′∞)Z)−f∗n}+

almost surely as .

Now observe that

 {g(μn(x′k)+Cn(x′k)Z)−f∗n}+≤supk|g(μn(x′k)+Cn(x′k)Z)|+|f⋆n|.

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

 E[{g(μn(x′k)+Cn(x′k)Z)−f∗n}+]→E[{g(μn(x′∞)+Cn(x′∞)Z)−f∗n}+],

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:

1. is a limit point of .

2. as .

3. For any fixed , as .

###### Proof.

First we prove that 1 implies 2. If is an element of , say , then, for , we have

 Kn(x′n)≲Kn0(x′n)→Kn0(x′∞)=Kn0(xn0)=0,

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

 Kn(x′n)=Cov(h(x′n)−μn(x′n))≲Cov(h(x′n)−h(xmn)).

Finally, since is not an element of , , and it follows from the continuity of that

 Cov(h(x′n)−h(xmn))=K(x′n,x′n)+K(xmn,xmn)−2K(x′n,xmn)→0,

and thus .

Now we prove that 2 implies 3. Using the Cauchy-Schwarz inequality in , we obtain

 ∥h(x′n)−μn(x′n)∥2≤∥Kn(x′n)∥122∥h∥H,

Thus,

 ∥h(x′∞)−μn(x′n)∥2 ≤∥h(x′∞)−h(x′n)∥2+|h(x′n)−μn(x′n)∥2 ≤∥h(x′∞)−h(x′n)∥2+∥Kn(x′n)∥122∥h∥H→0

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

 E[{g(μkn−1(xkn)+Ckn−1(xkn)Z)−f∗kn−1}+] →E[{g(h(~x))−f∗∞}+] =E[{f(~x)−f∗∞}+]=0,

but

 νkn−1=E[{g(μkn−1(xkn)+Ckn−1(xkn)Z)−f∗kn−1}+],

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 EI-CF).

Assume that the covariance function, , satisfies the GNEB property. Then, for any fixed and