Constrained Bayesian Optimization with Noisy Experiments

Randomized experiments are the gold standard for evaluating the effects of changes to real-world systems, including Internet services. Data in these tests may be difficult to collect and outcomes may have high variance, resulting in potentially large measurement error. Bayesian optimization is a promising technique for optimizing multiple continuous parameters for field experiments, but existing approaches degrade in performance when the noise level is high. We derive an exact expression for expected improvement under greedy batch optimization with noisy observations and noisy constraints, and develop a quasi-Monte Carlo approximation that allows it to be efficiently optimized. Experiments with synthetic functions show that optimization performance on noisy, constrained problems outperforms existing methods. We further demonstrate the effectiveness of the method with two real experiments conducted at Facebook: optimizing a production ranking system, and optimizing web server compiler flags.


Introduction
Many real-world systems, including software, distributed systems, and machine learning infrastructure, have multiple continuous parameters that affect outcomes of interest and are subject to various trade-offs. These parameters often have complex interactions that make it impossible to know a priori what the best configuration is. Measuring the outcomes of any given set of parameters often requires time-consuming experiments, commonly referred to as A/B tests in the Internet industry. As a result, many systems used in practice involve various constants that have been chosen with a limited amount of manual tuning.
Bayesian optimization is a powerful tool for solving black-box global optimization problems with expensive function evaluations . Most commonly, this process begins by evaluating a small number of randomly selected function values, and fitting a Gaussian process (GP) regression model to the results. The GP posterior provides an estimate of the function value at each point, as well as the uncertainty in that estimate. We then choose a new point at which to evaluate the function by balancing exploration (high uncertainty) and exploitation (best estimated function value). This is done by optimizing an acquisition function, which encodes the value of potential points in the optimization and defines the balance between exploration and exploitation. A common choice for the acquisition function is expected improvement (EI), which measures the expected value of the improvement at each point over the best observed point. Optimization then continues sequentially, at each iteration updating the model to include all past observations. Bayesian optimization has recently become an important tool for optimizing machine learning hyperparameters (Snoek et al. 2012), where in each iteration a machine learning model is fit to data and prediction quality is observed. Our work is motivated by a need to develop robust algorithms for optimizing via field experiments. There are three aspects of A/B tests that differ from the usual hyperparameter optimization paradigm. First, there are typically high noise levels when measuring performance of systems. Extensions of Bayesian optimization to handle noisy observations use heuristics to simplify the acquisition function that can have pathological behaviors for high noise levels. In particular, we show that they lead to over-exploitation. Second, there are almost always trade-offs involved in optimizing real systems: improving the quality of images may result in increased data usage; increasing cache sizes may improve the speed of a mobile application, but decrease reliability on some devices; optimizing click-through rates may result in decreases in quality. Practitioners have stressed the importance of considering multiple outcomes (Deng and Shi 2016), and noisy constraints must be incorporated into the optimization. Finally, it is often straightforward to run multiple A/B tests in parallel, with limited wall time in which to complete the optimization. Methods for batch optimization have been developed in the noiseless case; here we unify the approach for handling noise and batches.
This paper provides three main contributions.
• We develop an exact expression for the expected improvement under noisy observations and noisy constraints. We show with synthetic problems that in highnoise settings the heuristics that have previously been developed lead to overexploitation and perform poorly compared to our approach. • We derive a quasi-Monte Carlo approximation that allows us to efficiently optimize expected improvement without relying on simplifying heuristics. Our approach also provides a substantial speed-up over current approaches for noiseless batch optimization that use Monte Carlo integration. • We demonstrate the applicability of Bayesian optimization to A/B testing with real examples: optimizing a production ranking system, and optimizing web server performance.

Prior work on expected improvement
The EI acquisition function was introduced by Jones et al. (1998) for efficient optimization of expensive black-box functions. They considered an unconstrained problem , we first use GP regression to estimate f . Let g(x|D f ) be the GP posterior at x and f * = min i f (x i ) the current best observation. The EI of a candidate x is the expectation of its improvement over f * : The GP posterior g(x|D f ) is normally distributed with known mean µ f (x) and variance σ 2 f (x), so this expectation has an elegant closed form in terms of the Gaussian density and distribution functions: This function is easy to implement, easy to optimize, has strong theoretical guarantees (Bull 2011), and performs well in practice (Snoek et al. 2012).

Noisy observations
Suppose that we do not observe f (x i ), rather we observe y i = f (x i ) + i , where i is the observation noise, for the purposes of GP regression assumed to be i ∼ N (0, τ 2 i ).

Given noisy observations with uncertainty estimates
, GP regression proceeds similarly to the noiseless case and we obtain the GP posterior g(x|D f ).
Computing EI with observation noise is challenging because we no longer know the function value of the current best point, f * . Gramacy and Lee (2011) recognize this problem and propose replacing f * in (1) with the GP mean estimate of the best function value: g * = min x µ f (x). With this substitution, EI can be computed and optimized in a similar way as in the noiseless case. Picheny et al. (2013) extend this idea by substituting a GP quantile in the place of the mean, and then optimize expected improvement of that quantile.
Measuring EI relative to the GP mean can be a reasonable heuristic, but when noise levels are high it can have odd behavior. Unlike in the noiseless case, the expected improvement at the current best is generally positive: in (1), at the current best z = 0 and so α EI = σ f (x)/ √ 2π, which is positive. Fig. 1 illustrates how this can lead to clumping of candidates, using a GP fit to noisy observations of sin(x). The point that maximizes EI over the best GP mean actually repeats an earlier observation. Effectively, the Bayesian optimization has switched to a pure exploit strategy. We show in the experiments of Section 5 that this over-exploiting can hurt performance. Schonlau et al. (1998) extend EI to solve noiseless constrained optimization problems of the form min

Constraints
where the constraint functions c j (x) are also black-box functions that are observed together with f . As with f , we give each c j a GP prior and denote its posterior mean and variance as µ cj (x) and σ 2 cj (x). Let f * c be the value of the best feasible observation. The choice that maximizes EI relative to the GP mean is to repeat the best observation, thus being purely exploitative. Schonlau et al. (1998) define the improvement of a candidate x over f * c to be 0 if x is infeasible, and otherwise to be the usual improvement. Assuming independence between f and each c j given x, the expected improvement with constraints is then As with unconstrained EI, this quantity is easy to optimize and works well in practice (Gardner et al. 2014). When the observations of the constraint functions are noisy, a similar challenge arises as with noisy observations of f : We may not know which observations are feasible, and so cannot compute the best feasible value f * c . Gelbart et al. (2014) propose using the best GP mean value that satisfies each constraint c j (x) with probability at least 1 − δ j , for a user-specified threshold δ j (0.05 in their experiments). If there is no x that satisfies the constraints with the required probability, then they select the candidate that maximizes the probability of feasibility, regardless of the objective value. In a high-noise setting, this heuristic for setting f * c can lead to over exploitation in a similar way as Fig. 1. Suppose δ = 0.05 and the best observation has probability 0.94 of being feasible. This point will not be used for f * c and so will have positive expected improvement and can be duplicated as in Fig. 1.

Batch optimization
EI can be used for batch or asynchronous optimization by iteratively maximizing EI integrated over pending outcomes (Ginsbourger et al. 2011). Let x b 1 , . . . , x b m be m candidates whose observations are pending, and ] the corresponding unobserved outcomes at those points. Candidate m + 1 is chosen as the point that maximizes Because of the GP prior on f , the conditional posterior f b |D f has a multivariate normal distribution with known mean and covariance matrix. The integral in (3) does not have an analytic expression, but we can sample from p(f b |D f ) and so can use a Monte Carlo approximation of the integral. Snoek et al. (2012) describe this approach to batch optimization, and show that despite the Monte Carlo integration it is efficient enough to be practically useful for optimizing machine learning hyperparameters. Batch optimization of EI has not previously been studied in noisy settings.

Exact EI with noise
We now derive the exact form of EI under noisy observations and constraints without heuristics, and will see that it extends immediately to handle asynchronous optimization. The result will be an integral similar to that of (3), but in Section 4 we develop a substantially more efficient estimate than has previously been used for batch optimization.

Infeasibility in the noiseless setting
We build up from the noiseless case, where both objective and constraints are observed exactly. We begin by defining a utility function that gives the utility after n iterations of optimization. To correctly deal with noisy constraints later, we must explicitly consider the case where no observations are feasible. Let S = {i : c j (x i ) ≤ 0 ∀j} be the set of feasible observations. The utility function is Here M is the cost of not having a feasible solution. 1 As before, f * c is the objective value of the best feasible point after n iterations. The improvement in utility from iteration n to iteration n + 1 is We choose x n+1 to maximize the expected improvement under the posterior distributions of f (x) and c j (x). For convenience, let f n = [f (x 1 ), . . . , f (x n )] be the objective values at the observations, c n j = [c j (x 1 ), . . . , c j (x n )] the values for each constraint, and c n = [c n 1 , . . . , c n J ] all constraint observations. In the noiseless setting, f n and c n are known, the best feasible value f * c can be computed, and the EI with infeasibility is This extends the constrained EI of (2) to explicitly handle the case where there are no feasible observations. Without a feasible best, this acquisition function balances the expected objective value with the probability of feasibility, according to the cost M . As M gets large, it approaches the strategy of Gelbart et al. (2014) and maximizes the probability of feasibility. For finite M , however, given two points with the same probability of being feasible, this acquisition function will choose the one with the better objective value.

Noisy EI
We now extend the expectation in (5) to noisy observations and noisy constraints. This is done exactly by iterating the expectation over the posterior distributions of f n and c n given their noisy observations. Let D cj be the noisy observations of the constraint functions, potentially with heteroscedastic noise. Then, by their GP priors and assumed independence, These are the GP posteriors for the true (noiseless) values of the objective and constraints at the observed points. The means and covariance matrices of these posterior distributions have closed forms in terms of the GP kernel function and the observed data (Rasmussen and Williams 2006). Let D = {D f , D c1 , . . . , D c J } denote the full set of data. Noisy expected improvement (NEI) is then: This acquisition function does not have an analytic expression, but we will show in the next section that both it and its gradient can be efficiently estimated, and so it can be optimized. This approach extends directly to allow for batch or asynchronous optimization with noise and constraints following the approach of Section 2.3. The objective values at the observed points, f n , and at the earlier points in the batch, f b , are jointly normally distributed with known mean and covariance. The integral in (6) is over the true values of all previously sampled points. For batch optimization, we simply extend that integral to be over both the previously sampled points and over any pending observations. Replacing f n in (6) with [f n , f b ] and making the corresponding replacement for c n yields the formula for batch optimization.

Efficient quasi-Monte Carlo integration of noisy EI
For batch optimization in the noiseless unconstrained case, the integral in (3) is estimated with Monte Carlo (MC) sampling. To solve the higher dimensionality NEI integral in (6), we benefit from a more efficient integral approximation, for which we turn to quasi-Monte Carlo (QMC) methods.
QMC methods provide an efficient approximation of high-dimensional integrals on the unit cube as a sum of function evaluations: When t k are chosen from a uniform distribution on [0, 1] d , this is MC integration. The Central Limit Theorem provides a convergence rate of O(1/ √ N ) (Caflisch 1998). QMC methods provide faster convergence and lower error by using a better choice of t k . For the purposes of integration, random samples can be wasteful because they tend to clump; a point that is very close to another provides little additional information about a smooth f . QMC methods replace random samples for t k with a deterministic sequence that is constructed to be low-discrepancy, or space-filling. There are a variety of such sequences, and here we use Sobol sequences (Owen 1998). Theoretically, QMC methods achieve a convergence rate of O((log N ) d /N ), and typically achieve much faster convergence in practice (Dick et al. 2013). The main theoretical result for QMC integration is the Koksma-Hlawka theorem, which provides a deterministic bound on the integration error in terms of the smoothness of f and the discrepancy of t k (Caflisch 1998).
To use QMC integration to estimate the NEI in (6), we must transform that integral to the unit cube.
Proposition 1 (Dick et al. 2013). Let p(x|µ, Σ) be the multivariate normal density function and choose A such that Σ = AA . Then, The matrix A can be the Cholesky decomposition of Σ. We now apply this result to the NEI integral in (6).
α EIx (x|f n (u),c n (u))du. QMC methods thus provide an estimate for the NEI integral according to The transform AΦ −1 (u)+µ is the typical way that multivariate normal random samples are generated from uniform random samples u (Gentle 2009). Thus when each t k is chosen uniformly at random, this corresponds exactly to Monte Carlo integration using draws from the GP posterior. Using a quasirandom sequence {t 1 , . . . , t N } provides faster convergence, and so reduces the number of samples N required for optimization. As an illustration, Fig. 2 shows random draws from a multivariate normal alongside quasirandom "draws" from the same distribution, generated by applying the transform of Proposition 1 to a scrambled Sobol sequence. The quasirandom samples have better coverage of the distribution and will provide lower integration error.
The algorithm for computing NEI is summarized in Algorithm 1. In essence, we draw QMC samples from the posteriors for the true values of the noisy observations, and for each sampled "true" value, we compute noiseless EI using (5). The expensive steps in Algorithm 1 are kernel inference in line 1 and constructing the noiseless GP models in line 8. For the noiseless GP models we reuse the kernel hyperparameters from line 1, but must still invert each of their covariance matrices. Lines 1-8 (the QMC sampling and constructing the noiseless models for each sample) are independent of the candidate x. In practice, we do these steps once at the beginning of the optimization and cache the models. When we wish to evaluate the expected improvement at any point x during the optimization, we evaluate the GP posteriors at x for each of these cached models and compute EI (lines 10-13). This allows NEI to be quickly computed and optimized. For asynchronous or batch optimization, the posteriors in line 2 are those of both completed and pending observations, and all other steps remain the same.
The gradient of α EIx can be computed analytically, and so the gradient of (7) is available analytically and NEI can be optimized with standard nonlinear optimization methods. Besides the increased dimensionality of the integral, it is no harder to optimize (7) than it is to optimize (3), which has been shown to be efficient enough for practical use. Optimizing (3) for batch EI requires sampling from the GP posterior and fitting conditional models for each sample just as in Algorithm 1. We now show that

12
Use this GP posterior to compute EI as in the noiseless setting, α EIx in (5).

13
Increment α NEI = α NEI + 1 N α EIx . 14 return α NEI the QMC integration allows us to handle the increased dimensionality of the integral and makes NEI practically useful.

Synthetic problems
We use synthetic problems to provide a rigorous study of two aspects of our contribution. The first is to compare the performance of QMC integration to the MC approach used to estimate (3). We show that QMC integration allows the use of many fewer samples to achieve the same integration error and optimization performance, thus allowing us to efficiently optimize NEI. Second, we compare the optimization performance of NEI to that of several baseline approaches: a baseline using the heuristics described in Section 2, Spearmint EI (Snoek et al. 2012), and Spearmint constrained predictive entropy search (PESC) (Hernández-Lobato et al. 2015). Across many replicates of the simulation, the heuristic approach was over exploitative and did not perform as well. NEI also significantly outperformed both Spearmint EI and PESC.
We used four synthetic problems for our study. The equations and visualizations for each problem are given in the supplement. The first problem comes from Gramacy et al. (2016), and has two parameters and two constraints. The second is a constrained version of the Hartmann 6 problem, with six parameters and one constraint. The third problem is a constrained Branin problem used by Gelbart et al. (2014) and the fourth is a problem given by Gardner et al. (2014); these both have two parameters and one constraint. We simulated noisy objective and constraint observations by adding normally distributed noise to evaluations of the objective and constraints. Noise variances for each problem are given in the supplement.
In the experiments here and in Section 6, heteroscedastic GP regression was done using a Matérn 5/2 kernel, and posterior distributions for the kernel hyperparameters were inferred using the NUTS sampler (Hoffman and Gelman 2014). GP predictions were made using the posterior mean value for the hyperparameters. NEI was optimized using random restarts of the Scipy SLSQP optimizer. In a typical field experiment, including those of Section 6, we observe both the mean estimate and its standard error. All four methods (NEI, EI+heuristics, Spearmint EI, and PESC) were thus given the true noise variance of each observation.

Evaluating QMC performance
The first set of simulations analyze the performance of the QMC estimate in (7). We simulated computing NEI in a noisy, asynchronous setting by generating 10 points from a scrambled Sobol sequence, evaluating observations at the first 5 and then treating the next 5 as pending observations. We then evaluated the NEI integral of (6) using regular MC draws from the posterior, and using QMC draws as in Algorithm 1.
The dimensionality of the integral does not depend on the dimensionality of the optimization problem, rather it depends only on how many observations (completed and pending) there are. For the Gramacy problem, for instance, with 5 completed observations, 5 pending observations, and 2 constraints, the integral has a dimensionality of 30. The supplement shows the NEI surface for each problem's initialization. We computed a ground-truth for the NEI at its global optimum by estimating (6) with 10 4 regular MC samples drawn from the GP posterior.
We estimated NEI using a small number of samples (ranging from 4 to 50) with both standard MC and with the QMC method in (7) and measured the absolute error compared to the ground-truth. This was repeated 500 times for each number of samples. Fig. 3 shows that for the Gramacy problem, QMC reliably required half as many samples as MC to achieve the same integration error. Typically we are not interested in the actual value of NEI, rather we only want to find the optimizer. For 100 of the replicates, we optimized NEI using the small-sample approximation, and measured the Euclidean distance between the found optimizer and the ground-truth optimizer. Fig.  3 shows that the lower integration error leads to better optimization performance: 16 QMC samples achieved the same optimization performance as 50 MC samples. This same simulation was done for the other three problems, and similar results are shown in the supplement for them.

Optimization performance compared to heuristics and other methods
We compare optimization performance of NEI to using the heuristics of Section 2 to handle the noise in observations and constraints. Specifically, in the EI+heuristics method, we measure expected improvement relative to the best GP mean of points that satisfy the constraints in expectation. Batch optimization is done as described in Section 2.3, but using MC draws from a GP that includes the observation noise. The EI+heuristics method uses the same GP models and optimization routines as the NEI method, with the only difference being the use of heuristics in computing EI. In particular, the methods are identical in the absence of observation noise. In addition to the heuristics baseline, we also compare to two commonly used Bayesian optimization methods from the Spearmint package: Spearmint EI, and Spearmint PESC. Spearmint EI uses similar heuristics as EI+heuristics to handle noise, but also uses a different approach for GP estimation, different optimization routines, and other techniques like input warping . Spearmint PESC implements predictive entropy search -an entirely different class of acquisition function from EI which tries to minimize uncertainty about the location of the optimum. There are a number of other available packages for Bayesian optimization, however only Spearmint currently supports constraints and so our comparison is limited to these methods. Each optimization was begun from the same batch of 5 Sobol sequence points, after which Bayesian optimization was performed in 9 batches of 5 points each, for a total of 50 iterations. After each batch, noisy observations of the points in the batch were incorporated into the model. This simulation was repeated 100 times for each of the four problems, each with independent observation noise added to function and constraint evaluations. Fig. 4 shows the value of the best feasible point at each iteration of the optimization, for all four problems. NEI consistently performed the best of all of the methods. Compared to EI+heuristics, NEI was able to find better solutions with fewer iterations. Without noise, these two methods are identical; the improved performance comes entirely from correctly handling observation noise. PESC had equal performance as NEI on the Gardner problem, but performed worse even than EI+heuristics on the other problems. Computation time was similar for the four methods, all requiring around 10s per iteration; running times for each problem are given in the supplement.
The heuristic approaches have the potential to over exploit, and Gramacy problem that in these high-noise simulations there was substantial clumping in the EI+heuristics proposals, compared to those with NEI. Similar figures for the other three problems are given in the supplement. Being more exploitative in a noisy setting could potentially be advantageous by allowing the model to more accurately identify the best feasible solution. This was not the case in these simulations. Fig. 6 shows for the Gramacy problem that when the GP was used to select the best feasible solution after optimization, the best solution identified with NEI was better than that identified by EI+heuristics. Similar results are shown for the other three problems in the supplement.

Bayesian optimization with real field experiments
We present two case studies of how Bayesian optimization with NEI works in practice with real experiments at Facebook: an online field experiment and a server performance optimization problem. Both experiments involved tuning many continuous parameters simultaneously via noisy objectives and noisy constraints.

Optimizing production ranking systems
Advances in modeling, feature engineering, and hyperparameter optimization are critical to the performance of the models that make up a production machine learning system. However, the performance of a production machine learning system also depends on the inputs to the model, which often come from many interconnected retrieval and ranking systems, each of which is controlled by many tuning parameters (Bendersky et al. 2010, Covington et al. 2016. For example, an indexer may retrieve a subset of items which are then fed into a high-precision ranking algorithm. The indexer has parameters such as the number of items to retrieve at each stage and how different items are valued (Rodriguez et al. 2012). These parameters can have large effects on latency, relevance, and quality, and can often lead to increases in user-level outcomes that rival those of modeling improvements in the ranking algorithm. The proportion of the model estimated best-feasible points that were actually feasible were similar for both methods.
While Bayesian optimization has proven to be an effective tool for optimizing the performance of machine learning models operating in isolation (Snoek et al. 2012), the evaluation of an entire production system requires live A/B testing. Since outcomes directly affected by machine learning systems such as click-through rates and other types of interaction are heavily skewed (Kohavi et al. 2014), measurement error is on the same order as the effect size itself.
We used NEI to optimize a real ranking system. This system consists of an indexer that aggregates content from various sources and identifies items to be sent to a model for ranking. The indexer has six parameters that control the aggregation, and our goal was to optimize these parameters to improve the overall performance of the ranking system. We maximized an objective metric (e.g., interaction) subject to a lower bound on another metric (e.g., quality).
NEI is ideally suited for this type of field experiment: noise levels are significant relative to the effect size, multiple variants are tested simultaneously in a batch fashion, and there are constraints that must be satisfied (e.g., measures of quality). The experiment was conducted in two batches: a quasirandom initial batch of 31 configurations selected with a scrambled Sobol sequence, and a second batch which used NEI to propose 3 configurations. Fig. 7 shows the results of the experiment as change relative to baseline, with axes scaled by the largest effect. In this experiment, the objective and constraint were highly negatively correlated (ρ = 0.78). NEI proposed candidates near the constraint boundary, and with only three points was able to find a feasible configuration that improved over both the baseline and anything from the initial batch.

Objective metric
Quasi-random NEI Figure 7: Posterior GP predictions (means and 2 standard deviations) from an A/B test using NEI to generate a batch of 3 candidates. The goal was to maximize the objective, subject to a lower bound on the constraint. The shaded region is infeasible. NEI found a feasible point with significantly better objective value than both the baseline and the quasirandom initialization.

Optimizing web server performance
We applied Bayesian optimization with NEI to improve the performance of the web servers that power Facebook. Facebook is written in a mix of the PHP and Hack programming languages, and it uses the HipHop Virtual Machine (HHVM)  to execute the PHP/Hack code in order to serve web requests. HHVM is an open-source virtual machine containing a just-in-time (JIT) compiler to translate the PHP/Hack code into Intel x86 machine code at runtime so it can be executed. During the compilation process, HHVM's JIT compiler performs a large number of code optimizations aimed at improving the performance of the final machine code. For example, code layout optimization splits the hot and cold code paths in order to improve the effectiveness of the instruction cache by increasing the chances of the hot code remaining in the cache. How often a code block is executed to be considered hot is a tunable parameter inside the JIT compiler. As another example, function inlining eliminates the overhead of calling and returning from a function, with tunable parameters determining which kinds of functions should be inlined.
Tuning compiler parameters can be very challenging for a number of reasons. First, even seemingly unrelated compiler optimizations, such as function inlining and code layout, can interfere with one another by affecting performance of the processor's instruction cache. Second, there are often additional constraints that limit the viable optimization space. Function inlining, for example, can drastically increase code size and, as a result, memory usage. Third, accurate modeling of all the factors inside a processor is so difficult that the only reasonable way to compare the performance of two different configurations is by running A/B tests. Facebook uses a system called Perflab for running A/B tests of server configurations (Bakshy and Frachtenberg 2015). At a high-level, a Perflab experiment assigns two isolated sets of machines to utilize the two configurations. It then replays a representative sample of user traffic against these hosts at high load, while measuring performance metrics including CPU time, memory usage, and database fetches, among other things. Crucially, Perflab provides confidence intervals on these noisy measurements, characterizing the noise level and allowing for rigorous comparison of the configurations. The system is described in detail in Bakshy and Frachtenberg (2015). Each A/B traffic experiment takes several hours to complete, however we had access to several machines on which to run these experiments, and so could use asynchronous optimization to run typically 3 function evaluations in parallel.
We tuned 7 numeric run-time compiler flags in HHVM that control inlining and code layout optimizations. This was a real experiment that we conducted, and the results were incorporated into the mainstream open-source HHVM (Ottoni 2016). Parameter names and their ranges are given in the supplement. Some parameters were integers -these values were rounded after optimization for each proposal. The goal of the optimization was to reduce CPU time with a constraint of not increasing peak memory usage on the web server.
We initialized with 30 configurations that were generated via scrambled Sobol sequences and then ran 70 more traffic experiments whose configurations were selected using NEI. Fig. 8 shows the CPU time and probability of feasibility across iterations. In the quasirandom initialization, CPU time and memory usage were only weakly correlated (ρ = 0.21). CPU times shown were scaled by the maximum observed difference. The optimization succeeded in finding a better parameter configuration, with experiment 83 providing the most reduction in CPU time while also not increasing peak memory. Nearly all of the NEI candidates provided a reduction of CPU time relative to baseline, while also being more likely to be feasible: the median probability of feasibility in the initialization was 0.77, which increased to 0.89 for the NEI candidates.

Discussion
Properly handling noisy observations and noisy constraints is essential when tuning parameters of a system via sequential experiments with measurement error. If the measurement error is small relative to the effect size, Bayesian optimization using a heuristic EI can be successful. However, when the measurement noise is high, previously suggested heuristics are exploitative and their proposals clump. We showed this overexploitation in synthetic problems, and a corresponding degradation in optimization performance.
We derived an exact expression for expected improvement in the presence of noisy observations and constraints, which extends directly to handle batch optimization. Our NEI requires solving a higher dimensional integral than has previously been used for batch optimization, but we developed a QMC integration technique allowed the integral to be estimated efficiently enough for optimization. Even in the noiseless case, the QMC approach that we developed here could be used to speed up the batch optimization strategy of Snoek et al. (2012).
For simplicity, here we assumed independence of the constraints. Extensions for multi-task GPs follow easily from the method. We expect that by leveraging these types of models, additional efficiency can be gained.
There are a number of other acquisition functions for Bayesian optimization besides EI. Methods such as quantile EI (Picheny et al. 2013) and knowledge gradient (Scott Predictive entropy search has been discussed in the context of noise, constraints, and batch optimization (Hernández-Lobato et al. 2014, Shah and Ghahramani 2015. Unlike EI, predictive entropy search does not lead to nice closed forms or integrals that can be directly estimated, rather it requires involved approximations that are described as being difficult to implement (Hernández-Lobato et al. 2015). We are interested in production optimization systems that are used and maintained by teams, and so the straightforward implementation of NEI is especially valuable. We found that not only did NEI generally outperform PESC, but even EI+heuristics outperformed PESC in three of the four experiments. PESC has been compared to Spearmint EI on these same problems before, but in settings more similar to hyperparameter optimization than our noisy experiments setting. Hernández-Lobato et al.
(2014) evaluated PESC on unconstrained Branin and Hartmann6 problems, but with a very low noise level: 0.03, whereas in our experiments the noise standard deviation was 5 for Branin and 0.2 for Hartmann6. Hernández-Lobato et al. (2015) evaluated PESC on the Gramacy problem, but with no observation noise. These previous experiments were also fully sequential, whereas ours required producing batches of 5 proposals before updating the model. Shah and Ghahramani (2015) evaluated predictive entropy search on unconstrained Branin and Hartmann6 problems with no noise, but with batches of size 3. They found for both of these problems that Spearmint EI outperformed predictive entropy search. Metzen (2016) showed that entropy search can perform worse than EI because it does not take into account the correlations in the observed function values. This can cause it to be over-exploitative, and is an issue that would be exacerbated by high observation noise. NEI explicitly incorporates the correlations between observations by integrating over joint samples from the GP posterior (line 2 of Algorithm 1).
Spearmint EI performed worse than EI+heuristics, despite also being an implementation of EI with heuristics. The most significant difference between the two is the way in which the constraint heuristic was implemented. EI+heuristics measured EI relative to the best point that was feasible in expectation. Spearmint EI requires the incumbent best to be feasible with probability at least 0.99 for each constraint. In our experiments with relatively noisy constraints, there were many iterations in which there were no observations with a probability of feasibility above 0.99, in which case Spearmint EI ignores the objective and proposes points that maximize the probability of feasibility. The sensitivity of the results to the way in which the heuristics are implemented provides additional motivation for ending our reliance on them with NEI.
We demonstrated the efficacy of our method to improve the performance of machine learning infrastructure and a JIT compiler. Our method is widely applicable to many other empirical settings which naturally produce measurement error, both in online and offline contexts.