Mixed-Variable Bayesian Optimization

1 Introduction

Bayesian optimization (BO) [Moˇckus, 1975] is a wellestablished paradigm to optimize costly-to-evaluate, complex, black-box objectives that has been successfully applied to many scientific domains. Most of the existing BO literature focuses on objectives that have purely continuous domains, such as those arising in tuning of continuous hyperparameters of machine learning algorithms, recommender systems, and preference learning [Shahriari et al., 2016]. More recently, problems with purely discrete domains, such as food safety control and model-sparsification in multi-component systems [Baptista and Poloczek, 2018] have been considered.

However, many real-world optimization problems in science and engineering are of mixed-variable nature, involving both continuous and discrete input variables, and exhibit complex constraints. For example, tuning the hyperparameters of a convolutional neural network involves both continuous variables, e.g., learning rate and momentum, and discrete ones, e.g., kernel size, stride and padding. Also, these hyperparameters impose validity constraints, as some combinations of kernel size, stride and padding define invalid networks. Further examples of mixed-variable, potentially constrained, optimization problems include sensor placement [Krause et al., 2008], drug discovery [Negoescu et al., 2011], optimizer configuration [Hutter et al., 2011] and many others. Nonetheless, only few BO methods can address the unconstrained version of such problem and no existing method can handle the constrained one. This work introduces the first algorithm that can efficiently optimize mixed-variable functions subject to known constraints with provable convergence guarantees.

∗Equal contribution. †Work done while at ETH Zurich.

Related Work. Extending continuous BO methods [Shahriari et al., 2016] to mixed inputs requires ad-hoc relaxation methods to map the problem to a fully continuous one and rounding methods to map the solution back. This ignores the original domain structure, makes the solution quality dependent on the relaxation and rounding methods, and makes it hard to handle discrete constraints. Extending discrete BO methods [Baptista and Poloczek, 2018; Oh et al., 2019] to mixed inputs requires a discretization of the continuous domain part, the granularity of which is crucial: If it is too small, the domain becomes prohibitively large; if it is too large, the domain may only contain poorly performing values of the continuous inputs. Few BO methods address the mixed-variable setting. SMAC [Hutter et al., 2011] uses a random forest surrogate model. However, its frequentist uncertainty estimates may be too inaccurate to steer the sampling. TPE [Bergstra et al., 2011] uses kernel density estimation to find inputs that will likely improve upon and unlikely perform worse than the incumbent solution. While SMAC and TPE can handle hierarchical constraints, they cannot handle more general constraints over the discrete variables, e.g., cardinality constraints. They also lack convergence guarantees. Hyperband (HB) [Li et al., 2018] uses cheap but less accurate approximations of the objective to dynamically allocate resources for function evaluations. BOHB [Falkner et al., 2018] is the model-based counterpart of HB, based on TPE. They thus extend existing mixed-variable methods to the multi-fidelity setting rather than proposing new ones, which is complementary to our approach, rather than in competition with it. [Garrido-Merch´an and Hern´andez-Lobato, 2018] propose a Gaussian process kernel to model discrete inputs without rounding bias. Their method lacks guarantees and cannot handle discrete constraints. We instead use discrete optimizers for the acquisition function, which avoid bias by only making integer evaluations. Finally, while [Hern´andez-Lobato et al., 2015; Gardner et al., 2014; Sui et al., 2015] extend continuous BO methods to handle unknown constraints, no method can handle known discrete constraints in a mixed-variable domain.

Contributions. We introduce MIVABO, the first BO algorithm for efficiently optimizing mixed-variable functions subject to known linear and quadratic integer constraints, encompassing many of the constraints present in real-world domains (e.g. cardinality, budget and hierarchical constraints). It relies on a linear surrogate model that decouples the continuous, discrete and mixed components of the function using an expressive feature expansion (Sec. 3.1). We exploit the ability of this model to efficiently draw samples from the posterior over the objective (Sec. 3.2) by combining it with Thompson sampling, and show how to optimize the resulting constrained acquisition function (Sec. 3.3). While in continuous BO, optimizing the acquisition function is difficult but has well-established solutions, this is not true for mixedvariable spaces and doing this efficiently and accurately is a key challenge that hugely impacts the algorithm's performance. We also provide the first convergence analysis of a mixed-variable BO algorithm (Sec. 3.5). Finally, we demonstrate the effectiveness of MIVABO on a set of complex hyperparameter tuning tasks, where it outperforms state-of-theart methods and is competitive with human experts (Sec. 4).

2 Problem Statement

We consider the problem of optimizing an unknown, costlyto-evaluate function defined over a mixed-variable domain, accessible through noisy evaluations and subject to known linear and quadratic constraints. Formally, we aim to solve

minx∈X f(x) s.t. gc (x) ≥ 0, gd (x) ≥ 0, (1) where X ⊆ X c × X d with continuous subspace X c and discrete subspace X d . Both constraints gc (x) ≥ 0 over X c and gd (x) ≥ 0 over X d are known, and specifically gd (x) are linear or quadratic. We assume, that the domain of the continuous inputs is box-constrained and can thus, w.l.o.g., be scaled to the unit hypercube, X c = [0, 1]Dc . We further assume, w.l.o.g., that the discrete inputs are binary, i.e., vectors xd ∈ X d = {0, 1}Dd are vertices of the unit hypercube. This representation can effectively capture the domain of any discrete function. For example, a vector xd = [xdi ]Dd i=1 ∈ X d can encode a subset A of a ground set of Dd elements, such that xdi = 1 ⇔ ai ∈ A and xdi = 0 ⇔ ai /∈ A, yielding a set function. Alternatively, xd ∈ X d can be a binary encoding of integer variables, yielding a function defined over integers.

Background. BO algorithms are iterative black-box optimization methods which, at every step t, select an input xt ∈ X and observe a noise-perturbed output yt ≜ f(xt) + ϵ with ϵ iid∼ N(0, β-1 ), β > 0. As evaluating f is costly, the goal is to query inputs based on past observations to find a global minimizer x∗ ∈ argminx∈X f(x) as efficiently and accurately as possible. To this end, BO algorithms leverage two components: (i) a probabilistic function model (or surrogate), that encodes the belief about f based on the observations available, and (ii) an acquisition function α : X → R that expresses the informativeness of input x about the location of x∗, given the surrogate of f. Based on the model of f, we query the best input measured by the acquisition function, then update the model with the observation and repeat this procedure. The goal of the acquisition function is to simultaneously learn about inputs that are likely to be optimal and about poorly explored regions of the input space, i.e., to trade-off exploitation against exploration. Thus, BO reduces the original hard black-box optimization problem to a series of cheaper problems xt ∈ argmaxx∈X αt(x). However, in our case, these optimization problems involve mixed variables and exhibit linear and quadratic constraints and are thus still challenging. We now present MIVABO, an algorithm to efficiently solve the optimization problem in Eq. (1).

3 MIVABO Algorithm

We first introduce the linear model used to represent the objective (Sec. 3.1) and describe how to do inference with it (Sec. 3.2). We then show how to use Thompson sampling to query informative inputs (Sec. 3.3) and, finally, provide a bound on the regret incurred by MIVABO. (Sec. 3.5).

3.1 Model

We propose a surrogate model that accounts for both discrete and continuous variables in a principled way, while balancing two conflicting goals: Model expressiveness versus feasibility of Bayesian inference and of the constrained optimization of the mixed-variable acquisition function. Linear models defined over non-linear feature mappings, f (x) = w⊤φ(x), are a class of flexible parametric models that strike a good tradeoff between model capacity, interpretability and ease of use through the definition of features φ : X → RM . While the complexity of the model is controlled by the number of features, M, its capacity depends on their definition. Therefore, to make the design of a set of expressive features more intuitive, we treat separately the contribution to the objective f from the discrete part of the domain, from the continuous part of the domain, and from the interaction of the two,

where, for j ∈ {d, c,m}, φj (xj ) = [φji (xj )]Mj i=1 ∈ RMj and wj ∈ RMj are the feature and weight vector for the discrete, continuous and mixed function component, respectively.

In many real-world domains, a large set of features can be discarded a priori to simplify the design space. It is common practice in high-dimensional BO to assume that only low-order interactions between the variables contribute significantly to the objective, which was shown for many practical problems [Rolland et al., 2018; Mutn´y and Krause, 2018], including deep neural network hyperparameter tuning [Hazan et al., 2017]. Similarly, we focus on features defined over small subsets of the inputs. Formally, we consider

φ(x) = [φk(xk)]Mk=1 , where xk is a subvector of x containing exclusively continuous or discrete variables or a mix of both. Thus, the objective f(x) can be decomposed into a sum of low-dimensional functions fk(xk) ≜ wkφk(xk) defined over subspaces Xk ⊆ X with dim(Xk) ≪ dim(X). This defines a generalized additive model [Rolland et al., 2018; Hastie, 2017], where the same variable can be included in multiple subvectors/features. The complexity of this model is controlled by the effective dimensionality (ED) of the subspaces, which is crucial under limited computational resources. In particular, let ¯Dd ≜ maxk∈[M] dim(X d k ) denote the ED of the discrete component in Eq. (2), i.e. the dimensionality of the largest subspace that exclusively contains discrete variables. Analogously, ¯Dc and ¯Dm denote the EDs of the continuous and mixed component, respectively. Intuitively, the ED corresponds to the maximum order of the variable interactions present in f. Then, the number of features M ∈ O(D ¯Dd d +D ¯Dc c +(Dd +Dc) ¯Dm ) scales exponentially in the EDs only (as modeling up to L-th order interactions of N inputs requires ∑Ll=0 (Nl ) ∈ O(N L) terms), which are usually small, even if the true dimensionality is large.

Discrete Features φd . We aim to define features φd that can effectively represent the discrete component of Eq. (2) as a linear function, which should generally be able to capture arbitrary interactions between the discrete variables. To this end, we consider all subsets S of the discrete variables in X d (or, equivalently, all elements S of the powerset 2Xd of Xd) and define a monomial ∏j∈S xdj for each subset S (where for S = ∅, ∏j∈∅ xdj = 1). We then form a weighted sum of all monomials to yield the multi-linear polynomial wd⊤φd (xd ) = ∑S∈2Xd wS ∏j∈S xdj . This functional representation corresponds to the Fourier expansion of a pseudoBoolean function (PBF) [Boros and Hammer, 2002]. In practice, an exponential number of features can be prohibitively expensive and may lead to high-variance estimators as in BO one typically does not have access to enough data to robustly fit a large model. Alternatively, [Baptista and Poloczek, 2018; Hazan et al., 2017] empirically found that a second-order polynomial in the Fourier basis provides a practical balance between expressiveness and efficiency, even when the true function is of higher order. In our model, we also consider quadratic PBFs, wd⊤φd (xd ) = w∅ + ∑ni=1 w{i}xdi + ∑1≤i<j≤n w{i,j}xdi xdj , which induces the discrete feature representation φd (xd ) ≜ [1, {xdi }Dd i=1 , {xdi xdj }1≤i<j≤Dd ]⊤ and reduces the number of model weights to Md ∈ O(D2 d ).

Continuous Features φc . In BO over continuous spaces, most approaches are based on Gaussian process (GP) models [Williams and Rasmussen, 2006] due to their flexibility and ability to capture large classes of continuous functions. To fit our linear model formulation, we leverage GPs' expressiveness by modeling the continuous part of our model in Eq. (2) using feature expansions φc (xc ) that result in a finite linear approximation of a GP. One simple, yet theoretically sound, choice is the class of Random Fourier Features (RFFs) [Rahimi and Recht, 2008], which use Monte Carlo integration for a randomized approximation of a GP. Alternatively, one can use Quadrature Fourier Features [Mutn´y and Krause, 2018], which instead use numerical integration for a deterministic approximation, which is particularly effective for problems with low effective dimensionality. Both feature classes were successfully used in BO [Jenatton et al., 2017; Mutn´y and Krause, 2018]. In our experiments, we use RFFs approximating a GP with a squared exponential kernel, which we found to best trade off complexity vs. accuracy in practice.

Mixed Features φm. The mixed term should capture as rich and realistic interactions between the discrete and continuous variables as possible, while keeping model inference and acquisition function optimization efficient. To this end, we stack products of all pairwise combinations of features of the two variable types, i.e. φm(xd , xc ) ≜ [φdi (xd ) · φcj (xc )]⊤1≤i≤Md,1≤j≤Mc . This formulation provides a good trade-off between modeling accuracy and computational complexity. In particular, it allows us to reduce φm to the discrete feature representation φd when conditioned on a fixed assignment of continuous variables φc (and vice versa). This property is crucial for optimizing the acquisition function, as it allows us to optimize the mixed term of our model by leveraging the tools for optimizing the discrete and continuous parts individually. The proposed representation contains MdMc features, resulting in a total of M = Md + Mc + MdMc. To reduce model complexity, prior knowledge about the problem can be incorporated into the construction of the mixed features. In particular, one may consider the following approaches. Firstly, one can exploit a known interaction structure between variables, e.g., in form of a dependency graph, and ignore the features that are known to be irrelevant. Secondly, one can start by including all of the proposed pairwise feature combinations and progressively discard not-promising ones. Finally, for high-dimensional problems, one can do the opposite and progressively add pairwise feature combinations, starting from the empty set.

3.2 Model Inference

Let X1:t ∈ Rt×D be the matrix whose ith row contains the input xi ∈ X queried at iteration i, dimX = D, and let y1:t = [y1, . . . , yt]⊤ ∈ Rt be the array of the corresponding noisy function observations. Also, let Φ1:t ∈ Rt×M be the matrix whose ith row contains the featurized input φ(xi) ∈ RM . The formulation of f in Eq. (2) and the noisy observation model induce the Gaussian likelihood p(y1:t|X1:t, w) = N(Φ1:tw, β-1 I). To reflect our a priori belief about the weight vector w and thus f, we specify a prior distribution over w. A natural choice for this is a zero-mean isotropic Gaussian prior p(w|α) = N (0, α-1 I), with precision α > 0, which encourages w to be uniformly small, so that the final predictor is a sum of all features, each giving a small, nonzero contribution. Given the likelihood and prior, we infer the posterior p(w|X1:t, y1:t, α, β) ∝ p(y1:t|X1:t, w, β)p(w|α), which due to conjugacy is Gaussian, p(w|X1:t, y1:t) = N(m, S-1 ), with mean m = βS-1Φ⊤1:ty1:t ∈ RM and precision S = αI + βΦ⊤1:tΦ1:t ∈ RM×M [Williams and Rasmussen, 2006]. This simple analytical treatment of the posterior distribution over w is a main benefit of this model, which can be viewed as a GP with a linear kernel in feature space.

3.3 Acquisition Function

We propose to use Thompson sampling (TS) [Thompson, 1933], which samples weights ˜w ∼ p(w|X1:t, y1:t, α, β) from the posterior and chooses the next input by solving ̂x ∈ arg minx∈X ˜w⊤φ(x). TS intuitively focuses on inputs that are plausibly optimal and has previously been successfully applied in both discrete and continuous domains [Baptista and Poloczek, 2018; Mutn´y and Krause, 2018].

TS requires solving ̂x ∈ arg minx∈X ˜w⊤ t φ(x), which is a challenging mixed-variable optimization problem. However, as ˜w⊤ t φ(x) decomposes as in Eq. (2), we can naturally use an alternating optimization scheme which iterates between optimizing the discrete variables xd conditioned on a particular setting of the continuous variables xc and vice versa, until convergence to some local optimum. While this scheme provides no theoretical guarantees, it is simple and thus widely and effectively applied in many contexts where the objective is hard to optimize. In particular, we iteratively solve ̂xd ∈ arg minxd∈X d ( ˜wd⊤φd (xd ) + ˜wm⊤φm(xd , xc = ̂xc )), ̂xc ∈ arg minxc∈X c ( ˜wc⊤φc (xc ) + ˜wm⊤φm(xd = ̂xd , xc )). Importantly, using the mixed features proposed in Sec. 3.1, these problems can be optimized by purely discrete and continuous optimizers, respectively. This also holds in the presence of mixed constraints gm(x) ≥ 0 if those decompose accordingly into discrete and continuous constraints.

This scheme leverages independent subroutines for discrete and continuous optimization: For the discrete part, we exploit the fact that optimizing a second-order pseudoBoolean function is equivalent to a binary integer quadratic program (IQP) [Boros and Hammer, 2002], allowing us to exploit commonly-used efficient and robust solvers such as Gurobi or CPLEX. While solving general binary IQPs is NP-hard [Boros and Hammer, 2002], these optimizers are in practice very efficient for the dimensionalities we consider (i.e., Dd < 100). This approach allows us to use any functionality offered by these tools, such as the ability to optimize objectives subject to linear constraints Axd ≤ b, A ∈ RK×Dd , b ∈ RK or quadratic constraints xd⊤Qxd + q⊤xd ≤ b, Q ∈ RDd×Dd , q ∈ RDd , b ∈ R. For the continuous part, one can use optimizers commonly used in continuous BO, such as L-BFGS or DIRECT. In our experiments, we use Gurobi as the discrete and L-BFGS as the continuous solver within the alternating optimization scheme, which we always run until convergence.

3.4 Model Discussion

BO algorithms are comprised of three major design choices: the surrogate model to estimate the objective, the acquisition function to measure informativeness of the inputs and the acquisition function optimizer to select queries. Due to the widespread availability of general-purpose optimizers for continuous functions, continuous BO is mostly concerned with the first two design dimensions. However, this is different for mixed-variable constrained problems. We show in Sec. 4 that using a heuristic optimizer for the acquisition function optimization leads to poor queries and, therefore, poor performance of the BO algorithm. Therefore, the tractability of the acquisition function optimization influences and couples the other design dimensions. In particular, the following considerations make the choice of a linear model and TS the ideal combination of surrogate and acquisition function for our problem. Firstly, the linear model is preferable to a GP with a mixed-variable kernel as the latter would complicate the acquisition function optimization for two reasons: (i) the posterior samples would be arbitrary nonlinear functions of the discrete variables and (ii) it would be non-trivial to evaluate them at arbitrary points in the domain. In contrast, our explicit feature expansion solves both problems, while second order interactions provide a valid discrete function representation [Baptista and Poloczek, 2018; Hazan et al., 2017] and lead to tractable quadratic MIPs with capacity for complex discrete constraints. Moreover, Random Fourier Features approximate common GP kernels arbitrarily well, and inference in MIVABO scales linearly with the number of data points, making it applicable in cases where GP inference, which scales cubically with the number of data points, would be prohibitive. Secondly, TS induces a simple relation between the surrogate and the resulting optimization problem for the acquisition function, allowing to trade off model expressiveness and optimization tractability, which is a key challenge in mixed-variable domains. Finally, the combination of TS and the linear surrogate facilitates the convergence analysis described in Sec. 3.5, making MIVABO the first mixed-variable BO method with theoretical guarantees.

3.5 Convergence Analysis

Using a linear model and Thompson sampling, we can leverage convergence analysis from linearly parameterized multiarmed bandits, a well-studied class of methods for solving structured decision making problems [Abeille et al., 2017]. These also assume the objective to be linear in features φ(x) ∈ RM with a fixed but unknown weight vector w ∈ RM , i.e. E[f(x)|φ(x)] = w⊤φ(x), and aim to minimize the total regret up to time T: R(T) = ∑Tt=1 (f(x∗) - f(xt)). We obtain the following regret bound for MIVABO:

Proposition 1. Assume that the following assumptions hold in every iteration t = 1, . . . , T of the MIVABO algorithm:

1. ˜wt∼N(m, 24M lnT ln 1δ S-1 ), i.e. with scaled variance.
2. xt = argminx ˜w⊤φ(x) is selected exactly. 1

Prop. 1 follows from Theorem 1 in [Abeille et al., 2017] and works for infinite arms x ∈ X, |X| = ∞. In our setting, both the discrete and continuous Fourier features (and, thus, the mixed features) satisfy the standard boundedness assumption, such that the proof indeed holds. Prop. 1 implies no-regret, limT→∞R(T)/T = 0, i.e., convergence to the global minimum, since the minimum found after T iterations is no further away from f (x∗) than the mean regret R(T )/T . To our knowledge, MIVABO is the first mixed-variable BO algorithm for which such a guarantee is known to hold.

1 To this end, one can use more expensive but theoretically backed optimization methods instead of the alternating one, such as the powerful and popular dual decomposition [Sontag et al., 2011].

4 Experiments

We present experimental results on tuning the hyperparameters of two machine learning algorithms, namely gradient boosting and a deep generative model, on multiple datasets.

Experimental Setup. For MIVABO 2 , we set the prior variance α, observation noise variance β, and kernel bandwidth σ to 1.0, and scale the variance as stated in Prop. 1. We compare against SMAC, TPE, random search, and the popular GPyOpt BO package. GPyOpt uses a GP model with the upper confidence bound acquisition function [Srinivas et al., 2010], and accounts for mixed variables by relaxing discrete variables to be continuous and later rounding them to the nearest discrete neighbor. To separate the influence of model choice and acquisition function optimization, we also consider the MIVABO model optimized by simulated annealing (SA) (MIVABO-SA) and the GP approach optimized by SA (GP-SA). We compare against the SA-based variants only in constrained settings, using more principled methods in unconstrained ones. To handle constraints, SA assigns high energy values to invalid inputs, making the probability of moving there negligible. We use SMAC, TPE and GPyOpt and SA with their respective default settings.

4.1 Gradient Boosting Tuning

The OpenML database [Vanschoren et al., 2014] contains evaluations for various machine learning methods trained on several datasets with many hyperparameter settings. We consider extreme gradient boosting (XGBoost) [Chen and Guestrin, 2016], one of the most popular OpenML benchmarks, and tune its ten hyperparameters - three are discrete and seven continuous - to minimize the classification error on a held-out test set (without any constraints). We use two datasets, each containing more than 45000 hyperparameter settings. To evaluate hyperparameter settings for which no data is available, we use a surrogate modeling approach based on nearest neighbor [Eggensperger et al., 2015], meaning that the objective returns the error of the closest (w.r.t. Euclidean distance) setting available in the dataset. Fig. 1 shows that MIVABO achieves performance which is either significantly stronger than (left dataset) or competitive with (right dataset) the state-of-the-art mixed-variable BO algorithms on this challenging task. GPyOpt performs poorly, likely because it cannot account for discrete variables in a principled way. As compared to TPE and SMAC, MIVABO seems to benefit from more sophisticated uncertainty estimation.

4.2 Deep Generative Model (DGM) Tuning

DGMs recently received considerable attention in the machine learning community. Despite their popularity and importance, effectively tuning their hyperparameters is a major challenge. We consider tuning the hyperparameters of a variational autoencoder (VAE) [Kingma and Welling, 2014] composed of a convolutional encoder and a deconvolutional decoder [Salimans et al., 2015]. The VAEs are evaluated on stochastically binarized MNIST, as in [Burda et al., 2016], and FashionMNIST. They are trained on 60000 images for 32 epochs, using Adam with a mini-batch size of 128. We report the negative log-likelihood (NLL; in nats) achieved by the VAEs on a held-out test set of 10000 images, as estimated via importance sampling using 32 samples per test point. To our knowledge, no other BO paper considered DGM tuning.

2 We provide a Python implementation of MIVABO at https: //github.com/edaxberger/mixed_variable_bo.

VAE tuning is difficult due to the high-dimensional and structured nature of its hyperparameter space, and, in particular, due to constraints arising from dependencies between some of its parameters. We tune 25 discrete parameters defining the model architecture, e.g. the number of convolutional layers, their stride, padding and filter size, the number and width of fully-connected layers, and the latent space dimensionality. We further tune three continuous parameters for the optimizer and regularization. Crucially, mutual dependencies between the discrete parameters result in complex constraints, as certain combinations of stride, padding and filter size lead to invalid architectures. Particularly, for the encoder, the shapes of all layers must be integral, and for the decoder, the output shape must match the input data shape, i.e., one channel of size 28 × 28 for {Fashion}MNIST. The latter constraint is especially challenging, as only a small number of decoder configurations yield the required output shape. Thus, even for rather simple datasets such as {Fashion}MNIST, tuning such a VAE is significantly more challenging than, say, tuning a convolutional neural network for classification.

While MIVABO can conveniently capture these restrictions via linear and quadratic constraints, the competing methods cannot. To enable a comparison that is as fair as possible, we thus use the following sensible heuristic to incorporate the knowledge about the constraints into the baselines: If a method tries to evaluate an invalid parameter configuration, we return a penalty error value, which will discourage a model-based method to sample this (or a similar) setting again. However, for fairness, we only report valid observations and ignore all configurations that violated a constraint. We set the penalty value to 500 nats, which is the error incurred by a uniformly random generator. We investigated the impact of the penalty value (e.g., we also tried 250 and 125 nats) and found that it does not qualitatively affect the results.

Fig. 3 shows that MIVABO significantly outperforms the competing methods on this task, both on MNIST (left) and FashionMNIST (right). This is because MIVABO can naturally encode the constraints and thus directly optimize over the feasible region in parameter space, while TPE, SMAC and GPyOpt need to learn the constraints from data. They fail to do so and get stuck in bad local optima early on. The modelbased approaches likely struggle due to sharp discontinuities in hyperparameter space induced by the constraint violation penalties (i.e., as invalid configurations may lie close to wellperforming configurations). In contrast, random search is agnostic to these discontinuities, and thus notably outperforms the model-based methods. Lastly, GP-SA and MIVABO-SA struggle as well, suggesting that while SA can avoid invalid inputs, the effective optimization of complex constrained objectives crucially requires more principled approaches for acquisition function optimization, such as the one we propose. This shows that all model choices for MIVABO (as discussed in Sec. 3.4) are necessary to achieve such strong results.

Although log-likelihood scores allow for a quantitative comparison, they are hard to interpret for humans. Thus, for a qualitative comparison, Fig. 2 visualizes the reconstruction quality achieved on MNIST by the best VAE configuration found by all methods after 32 BO iterations. The VAEs were trained for 32 epochs each, as in Fig. 3. The likelihoods seem to correlate with the quality of appearance, and the model found by MIVABO arguably produces the visually most appealing reconstructions among all models. Note that while MIVABO requires more time than the baselines (see Fig. 4), this is still negligible compared to the cost of a function evaluation, which involves training a deep generative model. Finally, the best VAE found by MIVABO achieves 84.25 nats on MNIST when trained for 3280 epochs and using 5000 importance samples for log-likelihood estimation, i.e. the setting used in [Burda et al., 2016] (see Fig. 4). This is comparable to the performance of 82-83 nats achieved by human expert tuned models, e.g. as reported in [Salimans et al., 2015] (which use even more convolutional layers and a more sophisticated inference method), highlighting MIVABO's effectiveness in tuning complex deep neural network architectures.

5 Conclusion

We propose MIVABO, the first method for efficiently optimizing expensive mixed-variable black-box functions subject to linear and quadratic discrete constraints. MIVABO combines a linear model of expressive features with Thompson sampling, making it simple yet effective. Moreover, it is highly flexible due to the modularity of its components, i.e., the mixed-variable features, and the optimization oracles for the acquisition procedure. This allows practitioners to tailor MIVABO to specific objectives, e.g. by incorporating prior knowledge in the feature design or by leveraging optimizers handling specific types of constraints. We show that MIVABO enjoys theoretical convergence guarantees that competing methods lack. Finally, we empirically demonstrate that MIVABO significantly improves optimization performance as compared to state-of-the-art methods for mixed-variable optimization on complex hyperparameter tuning tasks.

Acknowledgements

This research has been partially supported by SNSF NFP75 grant 407540 167189. Matteo Turchetta was supported through the ETH-MPI Center for Learning Systems. Erik Daxberger was supported through the EPSRC and Qualcomm. The authors thank Josip Djolonga, Mojm´ır Mutn´y, Johannes Kirschner, Alonso Marco Valle, David R. Burt, Ross Clarke, Wolfgang Roth as well as the anonymous reviewers of an earlier version of this paper for their helpful feedback.