Statistical problems often involve learning about an intractable, real-valued function f (x) that can be evaluated at given values of p continuous input variables x = (x 1 , . . . , x p ) ∈ [0, 1] p . For example, physical experiments are often infeasible in many engineering problems so computer experiments are performed instead through evaluations of a computationally expensive simulator. Santner, Williams, and Notz (2003) overviewed the design and analysis of computer experiments and apply the methodology to understand the evolution of wildfires (Berk et al. 2002), to design a prosthesis device (Chang et al. 2001), and to optimize a helicopter blade design across 31 input variables (Booker et al. 1999). Jala et al. (2016) recently used computer experiments to assess the impact of electromagnetic exposure on fetuses. We focus our attention on optimization of a continuous, infinitely differentiable function, f , that is, (1) expensive to evaluate; (2) depends on a moderate to large number of variables, and ( 3) is measured with error.
Due to the assumed cost of evaluation, we desire an optimization strategy that requires few evaluations. For expensive f , a sequential design approach is commonly employed to find χ = arg max x f (x). The approach begins with the evaluation of f at an initial design of input settings. The resulting observations are modeled with a surrogate function, f , often taken to be a Gaussian process (GP) model, and χ is estimated from f . To improve this estimation, a new design point x * is chosen based on an acquisition function that assigns a numeric value to each potential design point, which is related to the point's expected ability to improve estimation of χ if it were added to the initial design. The function, f , is evaluated at x * , and f and χ are updated.
The sequential design process requires estimation of two optima at each step, that for f and the acquisition function. Although these functions are more tractable than f , they are still difficult to optimize in high dimensions (Kandasamy, Schneider, and Poczos 2015). Indeed, acquisition functions are often multimodal and contain regions where both the functions and their gradients are essentially zero, which is problematic for gradientbased optimization methods (Lizotte, Greiner, and Schuurmans 2012). Dimension reduction techniques are commonly employed to improve performance of maximizer estimation. Regis (2016) reduces the optimization space to a trust region centered at the current estimator χ . Djolonga, Krause, and Cevher (2013) assumed that f (x) = g(Ax) for some smooth function g(•) : R q → R and row-orthogonal matrix A ∈ R q×p with q < p. Their SI-BO algorithm uses low-rank approximation techniques to identify the subspace that supports f with a Bayesian bandit framework for optimization with respect to g. Wang et al. (2016) proposed the REMBO algorithm which uses a similar dimension reduction technique but identifies maximizers within randomly generated embeddings z = Ax where A is randomly generated.
A special case of the SI-BO and REMBO algorithms could require Ax to simply produce a selection of q < p input variables, that is, to have the algorithms remove variables from consideration. For example, the "importance" of each variable may be quantified through a sensitivity analysis that assesses the variability of f as x changes over each dimension (Shan and Wang 2010). If that variability is reasonably large (small), then the variable is called globally active (inactive). To this end, Linkletter et al. (2006) specified mixture priors on the GP parameters and use Markov chain Monte Carlo (MCMC) to determine the posterior probabilities of each variable being globally active. The globally inactive variables are removed from the design and analysis, and the resulting lower-dimensional space is easier to optimize across.
Even after employing global variable selection, further dimension reduction is possible if we were to focus our attention on a localized region of the input space, similar to the idea of a trust region. For example, local sensitivity examines the partial derivatives of f evaluated at a particular input x * (Oakley and O'Hagan 2004), say at χ . Bai et al. (2014) proposed two approaches for local variable selection. The first assumes a local linear model around some input x, assesses variable importance using local sensitivity, and implements a penalized LASSO framework to perform local variable selection. The second approach uses a forward/backward stepwise approach to choose the set of locally active variables around x using local linear estimators. Zhao et al. (2018) offered a generalization of the earlier algorithms and demonstrate set convergence (of the locally active variables) as well as parameter convergence. These articles, however, do not consider using the localized information to improve estimation of χ .
While the ultimate goal is to identify χ = arg max x f (x), one may not know how many additional iterations are affordable nor how many are needed to meet this goal. Therefore, we desire a sequential design approach that consistently increases f ( χ ) following each sequential run. In this article, we develop a Bayesian sequential design framework called sequential optimization in locally important dimensions (SOLID) that accomplishes this by performing global variable selection and localized variable selection around χ to optimize f . In Section 2, we review Bayesian estimation of a GP, Bayesian global variable selection for response surfaces, and two common acquisition functions, expected improvement (EI) and augmented EI. We introduce in Section 3 a new measure of local variable importance, based on local changes in f near χ after perturbing the posterior GP parameters. In Section 4, we detail the SOLID algorithm and illustrate it on a toy example. In Section 5, we compare SOLID with standard sequential optimization methods on multiple test functions and in Section 6, demonstrate SOLID's effectiveness on a robotics dataset from Vijayakumar and Schaal (2000). We find that SOLID provides larger values of f ( χ) in the first few evaluations of f , whereas standard sequential methods require more evaluations of f to obtain comparable values of f ( χ ). In Section 7, we discuss the advantages and disadvantages of using SOLID to sequentially optimize an expensive black-box function and propose some areas of further development.
Let X 0 denote the n 0 × p initial design matrix whose rows are the p input settings of the n 0 initial runs. The success of a sequential design strongly depends on the initial design (Crombecq, Laermans, and Dhaene 2011) and the statistical model used to make predictions. Space-filling designs, in which the inputs are "spread out" across the entire design space, are a popular choice for initial designs (Kleijnen et al. 2005) because they maximize the possibility of identifying potential regions that contain the optimum when we have no prior information about the function. In this article, we use maximin LHS designs (Joseph and Hung 2008) because their projection properties provide useful information for performing variable selection.
Evaluating f at each row of X 0 produces a response vector, y, from the model Y(x) = f (x) + where ∼ N(0, τ 2 ). A surrogate model, f , is constructed from y and is used to make predictions for an arbitrary input x. The surrogate model considered in this article assumes that f is a realization of a GP with mean function E[f (x)] = μ(x) and covariance function cov[f (x), f (x )] = σ 2 K(x, x ) for any two inputs x and x . Following Welch et al. (1992), we set μ(x) ≡ μ for all x. There are numerous choices for correlation functions, including Matérn, nonstationary correlation functions, and BSS-ANOVA (Reich, Storlie, and Bondell 2009). Although a nonstationary correlation function could be more appropriate, they can require a large number of design points for proper estimation, which we cannot afford for our problem of interest. Instead, we choose the squared exponential correlation function (Sacks, Schiller, and Welch 1989)
where γ 1 , . . . , γ p ≥ 0 are the correlation range parameters.
If γ k = 0, then varying x k across [0, 1] has no effect on the response.
The covariance function for f induces a covariance function for Y(x), which includes a nugget term τ 2 to account for random variation. Even for deterministic functions where Y(x) = f (x), including a nugget effect can protect against violations of model assumptions (Gramacy and Lee 2012). Letting Y i (x) denote the ith observation at x, we then have
Let V X be the n × n covariance matrix of y from design matrix X and let v(x) be the n × 1 vector of covariances between y and new observation Y(x). The prediction for f (x), conditional on y, is a Gaussian random variable with mean and variance
where denotes the vector of GP parameters (Gelman et al. 2004). Of course, needs to be estimated from the available data, which we do following Linkletter et al. (2006) which incorporates global variable selection, described next.
Each input variable x k influences f through its corresponding range parameter γ k in K(x, x ), where γ k = 0 implies that the input variable is globally inactive. Linkletter et al. (2006) placed positive mass on γ k = 0 through a mixture prior such that
) where u k is independent of b k , and θ ∼ Beta(a θ , b θ ) is the probability of each variable being globally active. More details on parameter priors are available in Appendix A.1.
The decision to declare an input variable globally active is based on the posterior probability Pr(b k = 1 | y) = Pr(γ k > 0 | y) ≡ bk . Variable k is declared globally inactive if bk < g where g ∈ (0, 1) is some threshold. Following Linkletter et al. (2006), a data-driven estimate of g may be found by augmenting the design with one or more random inputs and setting g to be the estimated probability of those variables being active. In this article, once variable k is deemed globally inactive, the kth column of X is permanently removed from future consideration and the remaining GP parameters are re-estimated. Henceforth, p will always reference the current number of variables that are deemed globally active at the current sequential step.
Each posterior draw t , t = 1, . . . , M, results in a new prediction surface ft = f (•| t ) and χ t , that is estimated from ft . We will also make use of the marginal prediction surface f = M -1 t ft and define the estimated global maximizer to be χ = arg max x f .
( 5 )
Note this estimator may differ from the alternative estimator M -1 t χt , the average of the maximizer posterior draws.
To determine a new design point to help identify χ , Jones, Schonlau, and Welch (1998) introduced the efficient global optimization (EGO) algorithm, which balances exploring the design space and honing in on areas likely containing χ. As introduced by Močkus (1975), the improvement at any x is
, where x opt is the row of X where y(x opt ) is the largest observed response in y. Since f is unknown, the EGO algorithm instead uses f to compute the expected improvement, EI( Jones, Schonlau, and Welch (1998) showed that EI can be written as
where
, and (•) and φ(•) are the CDF and PDF of a standard normal distribution, respectively. The next input is x * ≡ arg max x EI(x).
The EGO algorithm was built for deterministic computer simulations, where τ 2 = 0. The augmented EI criterion (Huang et al. 2006), or AEI, is more appropriate for nondeterministic functions 2006) stated that the x opt design point is chosen to reflect the user's degree of risk aversion, where ν = 1 represents a "willingness to trade 1 unit of predicted objective value for 1 unit of the standard deviation of prediction uncertainty. " See Brochu et al. (2010) for a discussion of other acquisition functions for identifying χ .
There are numerous algorithms to optimize AEI. Picheny and Ginsbourger (2014) optimized AEI through genetic optimization with derivatives, developed by Mebane and Sekhon (2011). Kleijnen (2015) constructed a space-filling design of candidate points C ⊂ [0, 1] p and sets the next design point to be x * = arg max x∈C AEI(x). These approaches are not immediately appealing for the problem at hand because (1) they may fail to find the true maximizer of AEI due to its multimodal nature in moderate to high dimensions and ( 2) they may encourage initial exploration of the design space, leading to poor initial improvement over the current χ. Sections 3 and 4 describe how we address these issues using a local variable selection algorithm and adaptive candidate sets to improve optimization of AEI and f .
Optimizing AEI around the current χ is appealing for multiple reasons. For one, it limits the possibility of the next design point to explore unobserved regions of the design space having high uncertainty under the surrogate model and instead encourages identification of a local optimum in an area that has been estimated to contain χ . Hence, it is more likely to lead to an updated χ with a larger f ( χ ) than if we chose a design point by globally optimizing AEI. Another reason is that, even when a variable is determined to be globally active, and hence has γk > 0, it may be that the variable is not important in a localized region of interest. Employing local variable selection can help to optimize AEI and update χ by significantly reducing the dimensionality of the optimization problem. This localized strategy would be especially beneficial for expensive functions with a potentially limited number of additional evaluations.
We define here a new measure of local importance defined on some region of the input space and, in the next section, develop a flexible algorithm that uses local variable selection to identify the maximizer for AEI. An appealing aspect of the proposed measure is its avoidance of expensive gradient evaluations. To motivate the measure, consider the two-dimensional toy example in Figure 1.
Both x 1 and x 2 are needed to describe the function globally, but there are areas that would require only one of the variables for optimization. Focusing on the localized, rectangular region labeled as (b), we make a baseline predicted surface using our current posterior estimates of γ 1 and γ 2 .
The global parameter γ 2 is likely greater than 0, but x 2 clearly does not substantially affect f in this region. Consider the alternative predicted surface restricted to this region, where γ 2 is temporarily set to 0 and γ 1 is the same as in the baseline surface. If this alternative predicted surface is similar to the baseline predicted surface, we would conclude that x 2 is locally inactive. Figure 2 shows the baseline and alternative predicted surfaces for this rectangular region and demonstrates how one would reach the conclusion that x 2 is locally inactive. Our approach assesses local variable importance within a neighborhood of each maximizer posterior draw χ t by comparing the baseline predicted surface, ft to each of the p alternative predicted surfaces, denoted by f k t , generated by temporarily fixing γ k = 0. To this end, we first generate q prediction points Q t from a truncated multivariate normal distribution
where δ controls how far the prediction points are spread from χ t , and the truncation keeps Q t within the [0, 1] p design space. We calculate the baseline and alternative predictions at the points in Q t , denoted ft (Q t ) and f k t (Q t ), respectively. We compare the baseline and alternative predictions using the squared correlation
If R 2 kt is close to one, then setting γ k = 0 did not greatly affect the predictions, offering evidence that x k is locally inactive.
It is possible that the χt 's will be dispersed across the input space and different variables are likely to be locally active with respect to different χ t (Bai et al. 2014). Anticipating this possibility, we average over the R 2 kt values and define the local importance L k of input k as
Then L k is an averaged measure of local importance across the posterior distribution of χ. We declare a variable to be locally active if L k ≥ ρ for 0 < ρ < 1 and let A denote the set of locally active variables. Algorithm 1 summarizes this procedure, including an additional step to perform the above calculations on only m < M of the posterior draws for computational reasons. The choice of the m draws should be done carefully to be representative of the entire posterior distribution. Declaring a variable to be locally active does not necessarily mean the function exhibits nonstationary behavior. For a function generated from a stationary GP, the parameters γ k describe the correlation with respect to the entire input space. As we focus our attention to a smaller region of interest, variables having γ k > 0 will start to appear unimportant. The larger γ k is, the smaller the region needs to be for this to happen. We allow the uncertainty of χ to dictate the size of the region. Even if f is generated from a nonstationary GP, our use of a surrogate Algorithm 1 Identifying locally active variables ). 9: return A = {k : L k ≥ ρ | χ } the set of locally active variables function assuming a stationary GP is necessary given our cost assumptions of evaluating f . To further support this, our toy example and numerical studies involve functions that exhibit nonstationary behavior.
In practice, finding the optimal AEI often involves reducing the optimization space, such as with trust regions. If χ is near χ , further exploration would be unnecessary, and restricting the search for the AEI maximizer to a small neighborhood of χ would be advantageous. Local variable selection can further reduce the dimension of this localized search space. We detail here the SOLID algorithm that uses global and local variable importance measures to improve estimation of χ .
In SOLID, rather than restrict the AEI search space for the kth locally active variable to be within some neighborhood centered at χ , we restrict the space to be
using the kth coordinate of the m χt 's and δ implemented in Algorithm 1. For the jth locally inactive variable, we set R δ j := χj , the jth coordinate of χ which is calculated from f . Let R δ denote the corresponding restricted search space.
Unlike the trust region in Regis (2016), the range of the R δ search is guided directly by the estimates of χk,t and explores only the locally active variables. Moreover, (11) incorporates the uncertainty of χ into the search space R δ . The inclusion of the δ parameter allows us to further expand the region if the distribution of the χk,1 , . . . , χk,m may be too narrow (perhaps due to selecting a smaller m for better computational performance). One could also use a different parameter than the δ used in Algorithm 1.
To search for the AEI optimum in R δ , we construct an |A|dimensional maximin LHS design L δ ⊂ R δ with c settings to evaluate AEI. For example, if only the first a < p variables are locally active, then the set of restricted candidate points C δ ⊆ R δ would be
It is possible that R δ is still too restrictive so we also consider a slightly larger space R A with R A k := [0, 1] for each locally active variable k, and R A j := χ j for each locally inactive variable j. This allows us to consider exploration of unobserved locations, but only within the locally active dimension. We again use a maximin LHS design L of dimension |A|, within R A for the c candidate points. These unrestricted candidate points C A are constructed in the same manner as in ( 12), where the column of each locally inactive variable j ∈ A c is χj 1 c . Note that R δ ⊆ R A but that R δ is more densely concentrated around the χt 's.
While it is possible to combine both C δ and C A into one large candidate set, we have found that the candidate points with the greatest AEI often all reside in one of the two sets. Whichever set has the largest AEI becomes the final set of candidate points C. Conceptually, this helps us see if SOLID is honing-in on a restricted space R δ or exploring the larger space R A . Using the |A|-dimensional gradient of AEI (see Picheny and Ginsbourger 2014), we conduct line searches from the five most promising candidates in C, restricting the search to lie within a ball of radius δ (as specified in ( 8)). After the line searches are complete, the one with the largest AEI is chosen as the next design point.
Thus far, we have described using the local variable selection results only for optimizing AEI, but they could also apply for the estimate of χ . One may be skeptical of doing this since the proposed local variable selection procedure uses the posterior draws χt which are calculated without local variable selection. By estimating χ across all globally active variables, we are likely to observe larger variation of χ t because we are optimizing in higher dimensions. Modifying the χ optimization space to R A should provide a more stable estimator.
We have described, up to this point, a single iteration of a sequential optimization algorithm detailing the global and local variable selection procedure, the localized AEI optimization, and the localized estimation of χ. The next iteration starts with evaluation of the recommended design point from the localized AEI optimization. We then re-estimate the GP parameters using MCMC with priors from the initial iteration. Recall, if at least one variable is deemed globally inactive at this next step then we would again re-estimate the GP parameters and perform global variable selection using only the globally active variables. Next we perform our Bayesian local variable selection across all globally active variables and maximize AEI in the new localized region. Note that previous results of the local variable selection algorithm are ignored. This way, any misclassification of locally active/inactive variables will not have long-lasting consequences. This flexibility also allows assessment of local importance to recalibrate when χ changes. The SOLID procedure is summarized in Algorithm 2.
We demonstrate SOLID in Figure 3 on the toy function (Figure 1) that includes a third, unimportant variable x 3 . We set the global and local thresholds to be g = 0.50 and ρ = 0.30, and set δ = 0.15. We considered c = 300 candidate points when optimizing AEI. Observations were generated with noise, τ 2 = 0.08. For simplicity, the marginal surfaces were built using m = 25 random draws using MCMC chains of length M = 500. We start with an initial maximin LHS design with n 0 = 10, p = 3, shown as in Figure 3 in the upper left panel.
In the first iteration, y(x opt ) = 7.71 and all three variables were deemed globally active. Local importance was assessed around the m = 25 posterior draws of χ t , shown as open circles in Figure 3. All three variables were also deemed locally active, with L 1 = 0.52, L 2 = 0.81, and L 3 = 0.76. It follows then that R A = [0, 1] 3 , the entire input space, and R δ (visualized in just the important dimensions as the shaded rectangle in the upper left panel of Figure 3) was bounded by 0.66 ≤ x 1 ≤ 1.00, 0.00 ≤ x 2 ≤ 0.59 and 0.44 ≤ x 3 ≤ 0.94. The localized optimum estimation step was not technically localized since R A = [0, 1] 3 and we determined χ0 = (0.88, 0.42, 0.57).
Next we optimized AEI to determine the next design point. The c = 300 candidate points were generated in R A and R δ , producing C A and C δ , respectively. The set C δ contained the point with the largest AEI and a line search optimization algorithm contained in this restricted space determined the next design point to be x * = (1.00, 0.31, 0.67). This point was added to X and we then evaluated Y(x * ).
In the second run, all variables were again found to be globally and locally active. The updated optimum was χ1 = (0.90, 0.40, 0.58). Here the unrestricted candidate set C A = [0, 1] 3 was preferred for AEI optimization and the next selected point was x * = (0.99, 0.47, 0.55), close to χ 1 .
In the third run, we found b3 = 0.49 < g = 0.5 so variable x 3 was permanently removed and the remaining GP parameters were re-estimated. Both x 1 and x 2 were still deemed globally active, as they should be. Their respective local importance measures were L 1 = 0.85 and L 2 = 0.29. With L 2 < ρ = 0.3, x 2 was declared locally inactive, and so R A and R δ were entirely contained within the x 1 dimension, both fixing x 2 at Algorithm 2 Summary of SOLID 1: Set n 0 , N (maximum number of evaluations), g, δ, ρ, M, m, c 2: Create an initial maximin LHS(n 0 , p) design, X 3: Generate y from Y(X) 4: for step i ∈ {0, . . . , N} do 5:
Obtain M posterior draws of t and χ t (Section 2.2); calculate f and χ .
Global variable selection: Remove variables with bk < g from X; if variables removed, repeat step (5) with new X 7:
Local variable selection: Implement Algorithm 1 with δ, ρ, and m < M χ t 's; store A 8:
Define restricted R δ and unrestricted R A search spaces 9:
Localized optimum estimation: Update estimate χ in R A using f ; store as χ i . Localized AEI estimation: Perform line search optimization to identify x * = arg max x∈C AEI(x) 13: Augment x * to X; generate Y(x * ) and add to y 14: return { χ0 , . . . , χ N } 0.41. Maximizing AEI, the next design point x * = (0.88, 0.41) was added to X. Note that the setting for the third input variable was no longer considered. If a value for that variable were required for the function to be evaluated, one could choose the corresponding coordinate from the previous optimum estimate, which in this case was 0.58 for x 3 .
For the remaining runs, x 1 and x 2 were always found to be globally active but x 2 was found locally inactive in runs 4 and 9. As design points were added and parameters were updated, f ( χ ) rose steadily: 7. 42, 9.41, 9.89, 9.92, 9.93, and 10.00, with 10.00 being the largest possible value.
We conducted a simulation study to evaluate the effects of global and local variable selection on sequential optimization. We compared four approaches: (1) GVS conducts global variable selection only; (2) SOLID, as described in Section 4; (3) Oracle uses only the known globally active variables (without performing any variable selection); and (4) None uses all variables. Within each simulation run, all four approaches used the same initial design, a maximin Latin hypercube design, and the same vector of initial responses. The responses were measured with error, τ 2 = 0.05.
We compared results for three different test functions in p = 15 dimensions, named Beach, Drum, and Simba. Although these are not conventionally high-dimensional functions, they are still sufficiently large enough to be important for practical concerns. All three functions have 6 truly globally active variables. The names for each function come from their visualization in the x 1 and x 2 subspace with all other variables set to χ j , their values in the global maximizer, visualized in the top row of Figure 4. The Beach function resembles a sandy beach along a pink sea; Drum resembles an oval shaped drum; and Simba is reminiscent of the scene from Disney's "The Lion King" (Hahn, Allers, and Minkoff 1994), where Rafiki holds Simba high up on Pride Rock, against the rolling hills and surrounding plains. The Beach function is constructed to have a local mode in a region far away from χ. Four variables are locally active around this local mode, but only x 1 , x 2 and x 3 are locally active around χ. The Drum function is primarily influenced by x 3 , but around χ, x 1 through x 5 are all locally active. The Simba function is especially challenging to optimize, since it has a large number of local modes involving all six globally active variables. Around χ, however, only x 1 , x 2 and x 3 are locally active. Results for three additional test functions commonly found in the optimization literature are shown in Section 3 of the supplementary materials.
It was important to choose n 0 to be large enough to construct a reasonable f without being so large that χ would easily be known. To that end, the initial designs for Beach and Drum had n 0 = 70 observations, and Simba had n 0 = 80. For all methods we used noninformative priors σ μ = 100, a η = b η = 0.1, a θ = b θ = 1. We set a u = 1 and b u = 10 so that E(u k ) = 10 and var(u k ) = 100 (see ( 4)). We ran MCMC chains of length M = 1000, of which m = 100 posterior draws were used for the marginal surfaces. We set δ = 0.30 and chose conservative global and local variable selection thresholds, g = 0.05 and ρ = 0.02. Of primary interest was determining the response of the true function f evaluated at χ across N = 25 additional runs. Because each initial designs gave a different χ 0 , our performance metric was relative improvement
To summarize across all N = 25 runs, we defined overall improvement 1 25
. Averaging across 100 simulated initial designs, we present the mean relative improvement for each added sequential design point, for each approach and test function in Figure 5. As expected, Oracle performed the best on Beach and Drum, largely because it optimized across only the 6 globally active variables. Across all test functions, None performed the worst, since it always optimized over a 15-dimensional space. SOLID had higher mean relative improvement than GVS for each of the first 10 runs on the Drum, and each of the first 20 runs on the Beach. For the Simba function, SOLID had higher mean improvement than GVS and Oracle for each of the first 6 runs. In Table 1, SOLID had significantly higher overall improvement than GVS (p-value < 0.001) in terms on Beach and Simba, and even outperformed Oracle on Simba (p = 0.003).
SOLID was able to achieve its enhanced performance, not only by permanently removing variables through global variable selection, but by honing in on more promising lowerdimensional subspaces. Figure 6 shows that our proposed measure of local importance successfully captured the locally active variables. Table 2 shows that across all three test functions, SOLID was optimizing over fewer variables than GVS, as defined by the number of variables explored by the AEI function.
We also compared the methods in terms of computational costs. Oracle, which knows the set of globally active variables and does not perform any variable selection, took 1.8 hr to add 25 design points, averaged across all test functions. For every hour that Oracle took to obtain 25 new design points, None took 1.4 hr, GVS took 4.0 hr, and SOLID took 5.3 hr.
Although computationally more expensive, if each evaluation of f is expensive, SOLID would still be preferable to GVS and None, since it required fewer evaluations of f to obtain equivalent or better estimates of χ . To see this, we compared the mean improvement value that Oracle achieves after 7 runs with the number of runs the other approaches needed to achieve at least that value (see Figure 5). On the Beach function, SOLID needed 10 runs, whereas GVS required 13, and none required 15. Similar patterns held for the Drum and Simba functions. As evidenced in Figure 5, the mean improvement value for GVS eventually met (on the Beach and Simba functions) or exceeded (on the Drum function) the value obtained by SOLID. The difference in how these two methods selected inputs at which to evaluate f next could explain this result. By selecting inputs whose values vary in all p globally active dimensions, GVS may be better able to identify truly globally inactive variables and correctly remove them from the design matrix. By removing more globally inactive variables than SOLID, GVS could eventually experience comparable or better mean improvement values.
In the first few sequential runs, however, the results favored SOLID.
The Sarcos robot dataset (Vijayakumar and Schaal 2000) consists of n = 44,484 observations and p = 21 input variables, available at www.gaussianprocess.org/gpml/data. The input variables are the positions, velocities, and accelerations of seven different points on a robot arm as it draws a figure eight (Vijayakumar, D'Souza, and Schaal 2005). We transform the inputs such that x ∈ [0, 1] 21 . The response variable Y(x) is the first of seven joint torque measurements (Parker 2015). Sequential optimization requires being able to evaluate the response surface at arbitrary input values, but this is not possible with the discrete Sarcos data. Therefore, for illustration purposes, we generated data assuming the true response surface is a kernel smoothed function. For any input x, we have
where S = {1, . . . , n}, x i are the observed inputs in the Sarcos dataset, and the kernel smoother is K
Based on 5-fold cross-validation minimizing the out-of-sample prediction MSE, the best bandwidth was h = 0.08272.
With the fields package in R, we randomly selected inputs that led to space-filling designs. We included ten times as many initial design points as dimensions (Loeppky, Sacks, and Welch 2009). We considered only 15 sequential evaluations in this analysis due to the computational demands. Using only the GVS and SOLID approaches, we evaluated the improvement (13) at each run i ∈ {1, . . . , 15}. We set g = 0.15 and ρ = 0.01 to provide a moderate amount of variable selection, and we set δ = 0.20 to emphasize local searches. τ 2 = 0.05. We used the same priors and number of MCMC chains and posterior samples as in Section 5.
We present results for 100 simulated datasets in Figure 7. SOLID achieved greater improvement than GVS over 100 simulated datasets at nearly every run of the sequential design. Comparing overall improvement, SOLID was significantly better than GVS (p-value < 0.001). SOLID consistently used fewer variables for optimization than GVS. We found that neither method removed any variables based on global variable selection. However, by the final run, SOLID used 15.88 variables during its optimization and design selection, compared to 21 for GVS. Figure 7 shows the proportion of datasets with globally and locally active variables at the final run. Of the 21 input variables, SOLID identified several as locally active. At the last run, variables Position 1, Acceleration 1, and Acceleration 4 were identified as locally active in 93%, 100%, and 96% of the simulated datasets.
When optimizing a function f (•) where each evaluation is expensive, one goal is to obtain the largest f ( χ ) in as few evaluations of f as possible. To that end, we proposed SOLID, a new method that measures local variable importance around χ and uses this information to optimize f in a sequential design. Whereas global variable selection permanently removes globally inactive variables, our local variable selection approach is flexible, adapting to the uncertainty of χ . We tailored local variable selection to optimize the search for both the maximizer of the AEI acquisition function and the global maximizer. Rather than exploring the entire p-dimensional space, SOLID examines only the locally active variables. In a simulation study, we found that our definition of local importance successfully captured the subset of locally active variables across multiple test functions. By reducing the optimization dimension global and local variable selection, our SOLID algorithm achieved higher f ( χ ) values compared to the standard methods, for a fixed number of sequential evaluations.
One reviewer pointed out that the presence of locally active variables could lead to nonstationary behavior in the response surface. Depending on the nature of nonstationarity, this could result in the estimated GP spatial range parameters, which assume a stationary covariance, as poor indicators of a variable's global importance. It is our intention that any variable that influences f anywhere in the input space be classified as globally active. Using a GP model with a nonstationary covariance function is the next step, though computational costs would increase and our definition of globally active would need to be modified. On the other hand, a nonstationary covariance function may not be necessary. The three test functions in Section 5 exhibit nonstationary behavior, yet the SOLID behaves well and rarely drops a globally active variable. It is also possible that a stationary covariance function would still able to predict f in a subregion near χ , which is one reason why we chose our local selection criterion to be based on prediction comparisons instead of directly inspecting the estimates of the spatial range parameters.
There are several ways that SOLID could be further improved. In our implementation, we fixed the global and local selection threshold parameters (i.e., g and ρ) but an adaptive approach could lead to improved performance. This is also true for the δ parameter used to declare local importance. Ideally this parameter would shrink as the design space is filled. Finally, there needs to be further exploration of the design impact on SOLID's performance, both in terms of the design size and how well the design supports estimation of the γ k . Following Linkletter et al. (2006), we used a maximin Latin hypercube design as our initial design. We have started to explore the distancedistributed designs discussed in Zhang, Cole, and Gramacy (2019) as they can better estimate the γ k .
One limitation of SOLID is its required computations to estimate local importance and its utilization of MCMC to estimate f . In instances where the underlying function is inexpensive to evaluate, it would be faster to use conventional sequential design approaches. Additionally, for experiments involving, say, more than 50 variables, the MCMC algorithm presented here for global and local variable selection could become excessively slow and other methods would be preferable, such as the random embedding approach (Wang et al. 2016) or by specifying an additive model (Kandasamy, Schneider, and Poczos 2015). That said, SOLID's local variable selection could be used for many functions f , without needing to know which or how many variables are locally active. An area of future work would be to incorporate aspects of these other methods within the SOLID framework, perhaps to perform fast initial screening. radius for local importance, we considered 100 simulations using a two-factor crossed design. Because so few variables were declared to be globally inactive, we disabled the global variable selection feature in the sensitivity analysis. We set the initial design to have size n = 210 and considered the improvement after 7 runs, limiting the number of runs due to computational costs. Results in Table A.1 show that a larger radius (δ = 0.60) and conservative threshold (ρ = 0.01 or ρ = 0.05) provided for the best performance. As long as a conservative local variable selection threshold was chosen, the results were not too sensitive to ρ.