Bayesian Optimization in High Dimensions via Random Embeddings

1 Introduction

Let f : X → R be a function on a compact subset X ⊆ RD . We address the following global optimization problem

We are particularly interested in objective functions f that may satisfy one or more of the following criteria: they do not have a closed-form expression, are expensive to evaluate, do not have easily available derivatives, or are non-convex. We treat f as a blackbox function that only allows us to query its function value at arbitrary x ∈ X. To address objectives of this challenging nature, we adopt the Bayesian optimization framework. There is a rich literature on Bayesian optimization, and we refer readers unfamiliar with the topic to more tutorial treatments [Brochu et al., 2009; Jones et al., 1998; Jones, 2001; Lizotte et al., 2011; Moˇckus, 1994; Osborne et al., 2009] and recent theoretical results [Srinivas et al., 2010; de Freitas et al., 2012].

In a nutshell, in order to optimize a blackbox function f, Bayesian optimization uses a prior distribution that captures our beliefs about the behavior of f, and updates this prior with sequentially acquired data. Specifically, it iterates the following phases: (1) use the prior to decide at which input x ∈ X to query f next; (2) evaluate f(x); and (3) update the prior based on the new data 〈x, f(x)〉. Step 1 uses a socalled acquisition function that quantifies the expected value of learning the value of f(x) for each x ∈ X. Bayesian optimization methods differ in their choice of prior and their choice of this acquisition function.

In recent years, the artificial intelligence community has increasingly used Bayesian optimization; see for example [Martinez-Cantin et al., 2009; Brochu et al., 2009; Srinivas et al., 2010; Hoffman et al., 2011; Lizotte et al., 2011; Azimi et al., 2012]. Despite many success stories, the approach is restricted to problems of moderate dimension, typically up to about 10; see for example the excellent and very recent overview in [Snoek et al., 2012]. Of course, for a great many problems this is all that is needed. However, to advance the state of the art, we need to scale Bayesian optimization to high-dimensional parameter spaces. This is a difficult problem: To ensure that a global optimum is found, we require good coverage of X, but as the dimensionality increases, the number of evaluations needed to cover X increases exponentially.

For linear bandits, Carpentier et al [2012] recently proposed a compressed sensing strategy to attack problems with a high degree of sparsity. Also recently, Chen et al [2012] made significant progress by introducing a two stage strategy for optimization and variable selection of high-dimensional GPs. In the first stage, sequential likelihood ratio tests with a couple of tuning parameters are used to select the relevant dimensions. This, however, requires the relevant dimensions to be axis-aligned with an ARD kernel. Chen et al provide empirical results only for synthetic examples (of up to 400 dimensions), but they provide key theoretical guarantees. Hutter et al [2011] used Bayesian optimization with random forests based on frequentist uncertainty estimates. Their method does not have theoretical guarantees for continuous optimization, but it achieved state-of-the-art performance for tuning up to 76 parameters of algorithms for solving combinatorial problems.

Many researchers have noted that for certain classes of problems most dimensions do not change the objective function significantly; examples include hyper-parameter optimization for neural networks and deep belief networks [Bergstra and Bengio, 2012] and automatic configuration of state-of-the-art algorithms for solving NP-hard problems [Hutter et al., 2013]. That is to say these problems have "low effective dimensionality". To take advantage of this property, [Bergstra and Bengio, 2012] proposed to simply use random search for optimization - the rationale being that points sampled uniformly at random in each dimension can densely cover each low-dimensional subspace. As such, random search can exploit low effective dimensionality without knowing which dimensions are important. In this paper, we exploit the same property, while still capitalizing on the strengths of Bayesian optimization. By combining randomization with Bayesian optimization, we are able to derive a new approach that outperforms each of the individual components.

Figure 1 illustrates our approach in a nutshell. Assume we know that a given D = 2 dimensional black-box function f(x1, x2) only has d = 1 important dimensions, but we do not know which of the two dimensions is the important one. We can then perform optimization in the embedded 1-dimensional subspace defined by x1 = x2 since this is guaranteed to include the optimum. This idea enables us to perform Bayesian optimization in a low-dimensional space to optimize a high-dimensional function with low intrinsic dimensionality. Importantly, it is not restricted to cases with axis-aligned intrinsic dimensions.

2 Bayesian Optimization

Bayesian optimization has two ingredients that need to be specified: The prior and the acquisition function. In this work, we adopt GP priors. We review GPs very briefly and refer the interested reader to [Rasmussen and Williams, 2006]. A GP is a distribution over functions specified by its mean function m(·) and covariance k(·, ·). More specifically, given a set of points x1:t, with xi ⊆ RD , we have

where K(x1:t, x1:t)i,j = k(xi, xj) serves as the covariance matrix. A common choice of k is the squared exponential function, but many other choices are possible depending on our degree of belief about the smoothness of the objective function.

An advantage of using GPs lies in their analytical tractability. In particular, given observations x1:n with corresponding values f1:t, where fi = f(xi), and a new point x∗ , the joint distribution is given by:

For simplicity, we assume that m(x1:t) = 0. Using the Sherman-Morrison-Woodbury formula, one can easily arrive at the posterior predictive distribution:

with data Dt = {x1:t, f1:t}, mean µ(x∗ |Dt) = k(x∗ , x1:t)K(x1:t, x1:t)-1 f1:t and variance σ(x∗ |Dt) = k(x∗ , x∗ ) - k(x∗ , x1:t)K(x1:t, x1:t)-1k(x1:t, x∗ ). That is, we can compute the posterior predictive mean µ(·) and variance σ(·) exactly for any point x∗ .

At each iteration of Bayesian optimization, one has to re-compute the predictive mean and variance. These two quantities are used to construct the second ingredient of Bayesian optimization: The acquisition function. In this work, we report results for the expected improvement acquisition function u(x|Dt) = E(max{0, ft+1(x) - f (x+)}|Dt) [Moˇckus, 1982; Bull, 2011]. In this definition, x+ = argmaxx∈{x1:t} f(x) is the element with the best objective value in the first t steps of the optimization process. The next query is: xt+1 = argmaxx∈X u(x|Dt). Note that this utility favors the selection of points with high variance (points in regions not well explored) and points with high mean value (points worth exploiting). We also experimented with the UCB acquisition function [Srinivas et al., 2010; de Freitas et al., 2012] and found it to yield similar results. The optimization of the closed-form acquisition function can be carried out by off-the-shelf numerical optimization procedures. The Bayesian optimization procedure is shown in Algorithm 1.

Algorithm 1 Bayesian Optimization

• 1: for t = 1, 2, . . . do
• 2: Find xt+1 ∈ RD by optimizing the acquisition function u: xt+1 = argmaxx∈X u(x|Dt).
• 3: Augment the data Dt+1 = {Dt, (xt+1, f(xt+1))}
• 4: end for

3 REMBO

Before introducing our new algorithm and its theoretical properties, we need to define what we mean by effective dimensionality formally.

Definition 1. A function f : RD → R is said to have effective dimensionality de, with de < D, if there exists a linear subspace T of dimension de such that for all x⊤ ∈ T ⊂ RD and x⊥ ∈ T ⊥ ⊂ RD , we have f(x) = f(x⊤ + x⊥) = f(x⊤), where T ⊥ denotes the orthogonal complement of T . We call T the effective subspace of f and T ⊥ the constant subspace.

This definition simply states that the function does not change along the coordinates x⊥, and this is why we refer to T ⊥ as the constant subspace. Given this definition, the following theorem shows that problems of low effective dimensionality can be solved via random embedding.

Theorem 2. Assume we are given a function f : RD → R with effective dimensionality de and a random matrix A ∈ RD×d with independent entries sampled according to N(0, 1) and d ≥ de. Then, with probability 1, for any x ∈ RD , there exists a y ∈ Rd such that f(x) = f(Ay).

Proof. Since f has effective dimensionality de, there exists an effective subspace T ⊂ RD , such that rank(T ) = de. Furthermore, any x ∈ RD decomposes as x = x⊤ + x⊥, where x⊤ ∈ T and x⊥ ∈ T ⊥. Hence, f(x) = f(x⊤+x⊥) = f(x⊤). Therefore, without loss of generality, it will suffice to show that for all x⊤ ∈ T , there exists a y ∈ Rd such that f (x⊤) = f (Ay).

Let Φ ∈ RD×de be a matrix, whose columns form an orthonormal basis for T . Hence, for each x⊤ ∈ T , there exists a c ∈ Rde such that x⊤ = Φc. Let us for now assume that ΦT A has rank de. If ΦT A has rank de, there exists a y such that (ΦT A)y = c. The orthogonal projection of Ay onto T is given by ΦΦT Ay = Φc = x⊤. Thus Ay = x⊤ + x′ for some x′ ∈ T ⊥ since x⊤ is the projection Ay onto T . Consequently, f(Ay) = f(x⊤ + x′ ) = f(x⊤).

It remains to show that, with probability one, the matrix ΦT A has rank de. Let Ae ∈ RD×de be a submatrix of A consisting of any de columns of A, which are i.i.d. samples distributed according to N(0, I). Then, ΦT ai are i.i.d. samples from N(0,ΦT Φ) = N(0de , Ide×de ), and so we have ΦT Ae, when considered as an element of Rd2e, is a sample from N(0d2e , Id2e×d2e ). On the other hand, the set of singular matrices in Rd2e has Lebesgue measure zero, since it is the zero set of a polynomial (i.e. the determinant function) and polynomial functions are Lebesgue measurable. Moreover, the Normal distribution is absolutely continuous with respect to the Lebesgue measure, so our matrix ΦT Ae is almost surely non-singular, which means that it has rank de and so the same is true of ΦT A, whose columns contain the columns of ΦT Ae.

Theorem 2 says that given any x ∈ RD and a random matrix A ∈ RD×d , with probability 1, there is a point y ∈ Rd such that f(x) = f(Ay). This implies that for any optimizer x⋆ ∈ RD , there is a point y⋆ ∈ Rd with f(x⋆ ) = f(Ay⋆ ). Therefore, instead of optimizing in the high dimensional space, we can optimize the function g(y) = f(Ay) in the lower dimensional space. This observation gives rise to our new algorithm Bayesian Optimization with Random Embedding (REMBO), described in Algorithm 2. REMBO first draws a random embedding (given by A) and then performs Bayesian optimization in this embedded space.

An important detail is how REMBO chooses the bounded region Y, inside which it performs Bayesian optimization. This is important because its effectiveness depends on the size of Y. Locating the optimum within Y is easier if Y is small, but if we set Y too small it may not actually contain the global optimizer (see Figure 2). In the following theorem, we show that we can choose Y in a way that only depends on the effective dimensionality de such that the optimizer of the original problem is contained in the low-dimensional space Algorithm 2 REMBO: Bayesian Optimization with Random Embedding

• 1: Generate a randommatrix A
• 2: Choose the set Y
• 3: for t = 1, 2, . . . do
• 4: Find yt+1 ∈ Rd by optimizing the acquisition function u: yt+1 = argmaxy∈Y u(y|Dt).
• 5: Augment the data Dt+1 = {Dt, (yt+1, f(Ayt+1)}
• 6: Update the kernel hyper-parameters.
• 7: end for

with constant probability. (If Ay is outside the box X, it is projected onto X .)

Theorem 3. Suppose we want to optimize a function f : RD → R with effective dimension de ≤ d subject to the box constraint X ⊂ RD , where X is centered around 0. Let us denote one of the optimizers by x⋆ . Suppose further that the effective subspace T of f is such that T is the span of de basis vectors. Let x⋆⊤ ∈ T ∩ X be an optimizer of f inside T . If A is a D×d random matrix with independent standard Gaussian entries, there exists an optimizer y⋆ ∈ Rd such that f(Ay⋆ ) = f(x⋆⊤) and ‖y⋆‖2 ≤ √de ϵ ‖x⋆⊤‖2 with probability at least 1- ϵ.

Proof. Since X is a box constraint, by projecting x⋆ to T we get x⋆⊤ ∈ T ∩ X . Also, since x⋆ = x⋆⊤ + x⊥ for some x⊥ ∈ T ⊥, we have f(x⋆ ) = f(x⋆⊤). Hence, x⋆⊤ is an optimizer. By using the same argument as appeared in Theorem 2, it is easy to see that with probability 1 ∀x ∈ T ∃y ∈ Rd such that Ay = x + x⊥ where x⊥ ∈ T ⊥. Let Φ be the matrix whose columns form a standard basis for T . Without loss of generality, we can assume that Φ = [Ide 0]T . Then, as shown in Proposition 2, there exists a y⋆ ∈ Rd such that ΦΦT Ay⋆ = x⋆⊤. Note that for each column of A, we have

Therefore ΦΦT Ay⋆ = x⋆⊤ is equivalent to By⋆ = ¯x⋆⊤ where B ∈ Rde×de is a random matrix with independent standard Gaussian entries and ¯x⋆⊤ is the vector that contains the first de entries of x⋆⊤ (the rest are 0's). By Theorem 3.4 of [Sankar et al., 2003], we have

Thus, with probability at least 1 - ϵ, ‖y⋆‖ ≤ ‖B-1‖2‖¯x⋆⊤‖2 = ‖B-1‖2‖x⋆⊤‖2 ≤ √de ϵ ‖x⋆⊤‖2.

Theorem 3 says that if the set X in the original space is a box constraint, then there exists an optimizer x⋆⊤ ∈ X that is de-sparse such that with probability at least 1 - ϵ, ‖y⋆‖2 ≤ √de ϵ ‖x⋆⊤‖2 where f(Ay⋆ ) = f(x⋆⊤). If the box constraint is X = [-1, 1]D (which is always achievable through rescaling), we have with probability at least 1 - ϵ that

Hence, to choose Y, we just have to make sure that the ball of radius de/ϵ satisfies (0, de ϵ ) ⊆ Y. Inmost practical scenarios, we found that the optimizer does not fall on the boundary which implies that ‖x⋆⊤‖2 < de. Thus setting Y too big may be unnecessarily wasteful. In all our experiments we set Y to be [-√d, √d]d .

Theorem 3 only guarantees that Y contains the optimum with probability at least 1 - ϵ; with probability δ ≤ ϵ the optimizer lies outside of Y. There are several ways to guard against this problem. One is to simply run REMBO multiple times with different independently drawn random embeddings. Since the probability of failure with each embedding is δ, the probability of the optimizer not being included in the considered space of k independently drawn embeddings is δk . Thus, the failure probability vanishes exponentially quickly in the number of REMBO runs, k. Note also that these independent runs can be trivially parallelized to harness the power of modern multi-core machines and large compute clusters.

3.1 Choice of Kernel

We begin our discussion on kernels with the definition of the squared exponential kernel between two points, y(1) , y(2) , on Y ⊆ Rd . Given a length scale r > 0, we define the corresponding squared exponential kernel as

It is possible to work with two variants of this kernel. First, we can use kd ℓ (y1 , y2 ) as above. We refer to this kernel as the low-dimensional kernel. We can also adopt an implicitly defined high-dimensional kernel on X :

where pX : RD → RD is the standard projection operator for our box-constraint: pX(y) = argminz∈X‖z - y‖2; see Figure 2.

Note that when using the high-dimensional kernel, we are fitting the GP in D dimensions. However, the search space is no longer the box X, but it is instead given by the much smaller subspace {pX(Ay) : y ∈ Y}. Importantly, in practice it is easier to maximize the acquisition function in this subspace. Our experiment on automatic algorithm configuration will show that indeed optimizing the acquisition function in low dimensions leads to significant gains of REMBO over standard Bayesian optimization.

The low-dimensional kernel has the benefit of only having to construct a GP in the space of intrinsic dimensionality d, whereas the high-dimensional kernel has to construct the GP in a space of extrinsic dimensionality D. However, the choice of kernel also depends on whether our variables are continuous, integer or categorical. The categorical case is important because we often encounter optimization problems that contain discrete choices. We define our kernel for categorical variables as:

where y(1) , y(2) ∈ Y ⊂ Rd and g defines the distance between 2 vectors. The function s maps continuous vectors to discrete vectors. In more detail, s(x) first projects x to [-1, 1]D to generate ¯x. For each dimension ¯xi of ¯x, s then maps ¯xi to the corresponding discrete parameters by scaling and rounding. In our experiments, following [Hutter, 2009], we defined g(x(1) , x(2) ) = |{i : x(1) i ̸= x(2) i }| so as not to impose an artificial ordering between the values of categorical parameters. In essence, we measure the distance between two points in the low-dimensional space as the distance between their mappings in the high-dimensional space.

Our demonstration of REMBO, in the domain of algorithm configuration, will use the high-dimensional kernel because the parameters in need of tuning are categorical. However, to provide the reader with a taste of the potential gains that the low-dimensional kernel could offer in continuous spaces, we will use a synthetic optimization problem. This synthetic example will also allow us to easily discuss different properties of the algorithm, including rotational invariance and dependency on d.

Finally, when using the high-dimensional kernel, the regret bounds of [Srinivas et al., 2010; Bull, 2011; de Freitas et al., 2012] apply. For the low-dimensional case, we have derived regret bounds that only depend on the intrinsic dimensionality. For lack of space, these are covered in a longer technical report version of this paper [Wang et al., 2013].

4 Experiments

For all our experiments, we used a single robust version of REMBO that automatically sets its GP's length scale parameter using a variant of maximum likelihood (see the technical report [Wang et al., 2013] for details). For each optimization of the acquisition function, this version runs both DIRECT [Jones et al., 1993] and CMA-ES [Hansen and Ostermeier, 2001] and uses the result of the better of the two. Source code for REMBO, as well as all data used in our experiments is publicly available at https://github.com/ziyuw/rembo.

4.1 Bayesian Optimization in a Billion Dimensions

The experiments in this section employ a standard de = 2dimensional benchmark function for Bayesian optimization, embedded in a D-dimensional space. That is, we add D - 2 additional dimensions which do not affect the function at all. More precisely, the function whose optimum we seek is f(x1:D) = g(xi, xj), where g is the Branin function (for its exact formula, see [Lizotte, 2008]) and where i and j are selected once using a random permutation. To measure the performance of each optimization method, we used the optimality gap: the difference of the best function value it found and the optimal function value.

We evaluate REMBO using a fixed budget of 500 function evaluations that is spread across multiple interleaved runs - for example, when using k = 4 interleaved REMBO runs, each of them was only allowed 125 function evaluations. We study the choices of k and d by considering several combinations of these values. The results in Table 1 demonstrate that interleaved runs helped improve REMBO's performance. We note that in 13/50 REMBO runs, the global optimum was indeed not contained in the box Y REMBO searched with d = 2; this is the reason for the poor mean performance of REMBO with d = 2 and k = 1. However, the remaining 37 runs performed very well, and REMBO thus performed well when using multiple interleaved runs: with a failure rate of 13/50=0.26 per independent run, the failure rate using k = 4 interleaved runs is only 0.264 ≈ 0.005. One could easily achieve an arbitrarily small failure rate by using many independent parallel runs. Using a larger d is also effective in increasing the probability of the optimizer falling into REMBO's box Y but at the same time slows down REMBO's convergence (such that interleaving several short runs loses its effectiveness).

Next, we compared REMBO to standard Bayesian optimization (BO) and to random search, for an extrinsic dimensionality of D = 25. Standard BO is well known to perform well in low dimensions, but to degrade above a tipping point of about 15-20 dimensions. Our results for D = 25 (see Figure 3, left) confirm that BO performed rather poorly just above this critical dimensionality (merely tying with random search). REMBO, on the other hand, still performed very well in 25 dimensions.

One important advantage of REMBO is that - in contrast to the approach of [Chen et al., 2012] - it does not require the effective dimension to be coordinate aligned. To demonstrate this fact empirically, we rotated the embedded Branin function by an orthogonal rotation matrix R ∈ RD×D . That is, we replaced f(x) by f(Rx). Figure 3 (middle) shows that REMBO's performance is not affected by this rotation. Finally, since REMBO is independent of the extrinsic dimensionality D as long as the intrinsic dimensionality de is small, it performed just as well in D = 1000000000 dimensions (see Figure 3, right). To the best of our knowledge, the only other existing method that can be run in such high dimensionality is random search.

For reference, we also evaluated the method of [Chen et al., 2012] for these functions, confirming that it does not handle rotation gracefully: while it performed best in the non-rotated case for D = 25, it performed worst in the rotated case. It could not be used efficiently for more than D = 1, 000. Based on a Mann-Whitney U test with Bonferroni multipletest correction, all performance differences were statistically significant, except Random vs. standard BO.

4.2 Automatic Configuration of aMixed Integer Linear Programming Solver

State-of-the-art algorithms for solving hard computational problems tend to parameterize several design choices in order to allow a customization of the algorithm to new problem domains. Automated methods for algorithm configuration have recently demonstrated that substantial performance gains of state-of-the-art algorithms can be achieved in a fully automated fashion [Moˇckus et al., 1999; Hutter et al., 2010; Vallati et al., 2011; Bergstra et al., 2011; Wang and de Freitas, 2011]. These successes have led to a paradigm shift in algorithm development towards the active design of highly parameterized frameworks that can be automatically customized to particular problem domains using optimization [Hoos, 2012; Bergstra et al., 2012].

It has recently been demonstrated that many algorithm configuration problems have low dimensionality [Hutter et al., 2013]. Here, we demonstrate that REMBO can exploit this low dimensionality even in the discrete spaces typically encountered in algorithm configuration. We use a configuration problem obtained from [Hutter et al., 2010], aiming to configure the 40 binary and 7 categorical parameters of lpsolve 1 , a popular mixed integer linear programming solver that has been downloaded over 40000 times in the last year. The objective is to minimize the optimality gap lpsolve can obtain in a time limit of five seconds for a mixed integer programming (MIP) encoding of a wildlife corridor problem from computational sustainability [Gomes et al., 2008]. Algorithm configuration aims to improve performance for a representative set of problem instances, and effective methods need to solve two orthogonal problems: searching the parameter space effectively and deciding how many resources to spend in each evaluation (to trade off computational overhead and over-fitting). Our contribution is for the first of these problems; to focus on how effectively the different methods search the parameter space, we only consider configuration on a single problem instance (i.e., a deterministic blackbox optimization problem; further work is required to tackle general algorithm configuration problems with randomized algorithms and distributions of problem instances).

1 http://lpsolve.sourceforge.net/

Due to the discrete nature of this optimization problem, we could only apply REMBO using the high-dimensional kernel for categorical variables kD λ (y(1) , y(2) ) described in Section 3.1. While we have not proven any theoretical guarantees for discrete optimization problems, REMBO appears to effectively exploit the low effective dimensionality of at least this particular optimization problem.

Figure 4 (top) compares BO, REMBO, and the baseline random search against ParamILS (which was used for all configuration experiments in [Hutter et al., 2010]) and SMAC [Hutter et al., 2011]. ParamILS and SMAC were specifically designed for the configuration of algorithms with many discrete parameters and define the current state of the art for this problem. Nevertheless, here SMAC and our vanilla REMBO method performed best. Based on a Mann-Whitney U test with Bonferroni multiple-test correction, they both yielded statistically significantly better results than both Random and standard BO; no other performance differences were significant. The figure only shows REMBO with d = 5 to avoid clutter, but we did not optimize this parameter; the only other value we tried (d = 3) resulted in indistinguishable .

As in the synthetic experiment, REMBO's performance could be further improved by using multiple interleaved runs. However, it is known that multiple independent runs also benefit the other procedures, especially ParamILS [Hutter et al., 2012]. Thus, to be fair, we re-evaluated all approaches using interleaved runs. Figure 4 (bottom) shows that ParamILS and REMBO benefitted most from interleaving k = 4 runs. However, the statistical test results did not change, still showing that SMAC and REMBO outperformed Random and BO, with no other significant performance differences.

5 Conclusion

This paper has shown that it is possible to use random embeddings in Bayesian optimization to optimize functions of high extrinsic dimensionality D, provided that they have low intrinsic dimensionality de. The new algorithm, REMBO, only requires a simple modification of the original Bayesian optimization algorithm; namely multiplication by a random matrix. We confirmed REMBO's independence of D empirically by optimizing low-dimensional functions embedded in high dimensions. Finally, we demonstrated that REMBO achieves excellent performance for optimizing the 47 discrete parameters of a popular mixed integer programming solver, thereby providing further evidence for the observation (already put forward by Bergstra, Hutter and colleagues) that, for many problems of great practical interest, the number of important dimensions indeed appears to be much lower than their extrinsic dimensionality. Of course, we do not yet know how many practical optimization problems fall within the class of problems where REMBO applies. With time and the release of the code, we will find more evidence. For the time being, the success achieved in the examples presented in this paper is very encouraging.