Mixed-Variable Bayesian Optimization

TL;DR

MiVaBo is a novel Bayesian optimization algorithm designed for efficient optimization of mixed-variable functions with discrete constraints, demonstrating superior sample efficiency in hyperparameter tuning tasks.

Contribution

Introduces MiVaBo, the first BO method capable of handling complex discrete constraints in mixed-variable domains with convergence guarantees.

Findings

MiVaBo outperforms existing methods in hyperparameter tuning tasks.

Provides the first convergence analysis for mixed-variable BO.

Efficiently handles complex discrete constraints in optimization.

Abstract

The optimization of expensive to evaluate, black-box, mixed-variable functions, i.e. functions that have continuous and discrete inputs, is a difficult and yet pervasive problem in science and engineering. In Bayesian optimization (BO), special cases of this problem that consider fully continuous or fully discrete domains have been widely studied. However, few methods exist for mixed-variable domains and none of them can handle discrete constraints that arise in many real-world applications. In this paper, we introduce MiVaBo, a novel BO algorithm for the efficient optimization of mixed-variable functions combining a linear surrogate model based on expressive feature representations with Thompson sampling. We propose an effective method to optimize its acquisition function, a challenging problem for mixed-variable domains, making MiVaBo the first BO method that can handle complex constraints over the discrete variables. Moreover, we provide the first convergence analysis of a mixed-variable BO algorithm. Finally, we show that MiVaBo is significantly more sample efficient than state-of-the-art mixed-variable BO algorithms on several hyperparameter tuning tasks, including the tuning of deep generative models.

Tables

Table 1. Table 1 : Hyperparameters of the XGBoost algorithm. 10 parameters, 7 of which are continuous, and 3 of which are discrete.

Name	Type	Domain
booster	discr.	[’gbtree’, ’gblinear’]
nrounds	discr.	$[3,5000]$
alpha	contin.	$[0.000985,1009.209690]$
lambda	contin.	$[0.000978,999.020893]$
colsample_bylevel	contin.	$[0.046776,0.998424]$
colsample_bytree	contin.	$[0.062528,0.999640]$
eta	contin.	$[0.000979,0.995686]$
max_depth	discr.	$[1,15]$
min_child_weight	contin.	$[1.012169,127.041806]$
subsample	contin.	$[0.100215,0.999830]$

Table 2. Table 2 : Hyperparameters of the VAE. The architecture of the VAE (if all layers are enabled) is C1-C2-F1-F2-z-F3-F4-D1-D2 , with C denoting a convolutional (conv.) layer, F a fully-connected (fc.) layer, D a deconvolutional (deconv.) layer and z the latent space. Layers F2 and F3 have fixed sizes of 2 ​ d z 2 subscript 𝑑 𝑧 2d_{z} and d z subscript 𝑑 𝑧 d_{z} units respectively, where d z subscript 𝑑 𝑧 d_{z} denotes the dimensionality of the latent space z . The domain of the number of units of the fc. layers F1 and F4 is discretized with a step size of 64, i.e. [ 0 , 64 , 128 , … , 832 , 896 , 960 ] 0 64 128 … 832 896 960 [0,64,128,\ldots,832,896,960] , denoted by [ 0 ​ … ​ 960 ] delimited-[] 0 … 960 [0\ldots 960] in the table for brevity. For d z subscript 𝑑 𝑧 d_{z} , the domain [ 16 ​ … ​ 64 ] delimited-[] 16 … 64 [16\ldots 64] refers to all integers within that interval.

#	Name	Type	Domain	Bits
1	Number of conv. layers in encoder	discrete	[0,1,2]	2
	Parameters of C1
2	Number of channels of C1	discrete	[4,8,16,24]	2
3	Stride of C1	discrete	[1,2]	1
4	Filter size of C1	discrete	[3,5]	1
5	Padding of C1	discrete	[0,1,2,3]	2
	Parameters of C2
6	Number of channels of C2	discrete	[8,16,32,48]	2
7	Stride of C2	discrete	[1,2]	1
8	Filter size of C2	discrete	[3,5]	1
9	Padding of C2	discrete	[0,1,2,3]	2
10	Number of fc. layers in encoder	discrete	[0,1,2]	2
11	Number of units of F1	discrete	[0…960]	4
12	Dimensionality $d_z$ of z	discrete	[16…64]	6
13	Number of fc. layers in decoder	discrete	[0,1,2]	2
14	Number of units of F4	discrete	[0…960]	4
15	Number of deconv. layers in decoder	discrete	[0,1,2]	2
	Parameters of D1
16	Number of channels of D1	discrete	[8,16,32,48]	2
17	Stride of D1	discrete	[1,2]	1
18	Filter size of D1	discrete	[3,5]	1
19	Padding of D1	discrete	[0,1,2,3]	2
20	Output padding of D1	discrete	[0,1,2,3]	2
	Parameters of D2
21	Number of channels of D2	discrete	[4,8,16,24]	2
22	Stride of D2	discrete	[1,2]	1
23	Filter size of D2	discrete	[3,5]	1
24	Padding of D2	discrete	[0,1,2,3]	2
25	Output padding of D2	discrete	[0,1,2,3]	2
26	Learning rate	continuous	$[{10}^{-4},{10}^{-2}]$	-
27	Learning rate decay factor	continuous	$[0.5,1.0]$	-
28	Weight decay regularization	continuous	$[{10}^{-6},{10}^{-2}]$	-
	Total			50

Table 3. Table 3 : Mean plus/minus one standard deviation of the negative test log-likelihood over 8 random initializations, achieved by the best VAE configuration found by SMAC, TPE and GPyOpt after 16 BO iterations, for constraint violation penalties of 500, 250 and 125 nats. Performance values of MiVaBo and random search (which are not affected by the penalty) are included for reference.

Algorithm	Penalty (nats)
	500	250	125
SMAC	$113.0\pm 1.8$	$112.1\pm 1.8$	$111.1\pm 1.6$
TPE	$108.8\pm 1.2$	$108.1\pm 1.3$	$108.1\pm 1.3$
GPyOpt	$108.5\pm 1.1$	$108.5\pm 0.6$	$106.5\pm 1.4$
RS			$106.3\pm 0.9$
MiVaBo			$94.4\pm 0.8$

Table 4. Table 4 : Mean plus/minus one standard deviation of the number of constraint violations by SMAC, TPE, GPyOpt and random search within 16 BO iterations over 8 random initializations, for constraint violation penalties of 500, 250 and 125 nats.

Algorithm	Penalty (nats)
	500	250	125
SMAC	$37\pm 21.7$	$36\pm 21.9$	$28\pm 11.6$
TPE	$67\pm 21.3$	$68\pm 22.2$	$68\pm 22.2$
GPyOpt	$36\pm 19.3$	$32\pm 18.0$	$27\pm 10.4$
Random search	$71\pm 25.5$	$71\pm 25.5$	$71\pm 25.5$

Equations

$${\min }_{\mathbf{x}\in \mathcal{X}}f(\mathbf{x})\text{s.t. }{g}^c(\mathbf{x})\ge 0,g^d(\mathbf{x})\ge 0,$$

$$f(\mathbf{x})={\sum }_{j\in \{d,c,m\}}{{\mathbf{w}}^j}^{\top }{\mathbf{\varphi }}^j({\mathbf{x}}^j)$$

$$W_{\text{out}}^e=(W_{\text{in}}^e-F^e+P^e)/S^e+1\in \mathbb{N}$$

$$W_{\text{out1}}^e=(28-F_1^e+P_1^e)/S_1^e+1$$

$$W_{\text{out2}}^e=(W_{\text{out1}}^e-F_2^e+P_2^e)/S_2^e+1$$

$$W_{\text{out}}^d=(W_{\text{in}}^d-1)\times S^d+F^d-2P^d+O^d$$

$$W_{\text{out}}^d=(W_{\text{in}}^d-1)\times S_1^d+F_1^d-2P_1^d+O_1^d,$$

$$28=(W_{\text{out}}^d-1)\times S_2^d+F_2^d-2P_2^d+O_2^d,$$

$$\text{MAP}(\mathbf{\theta })=\underset{\mathbf{x}}{\max }\left\{{\theta }^d({\mathbf{x}}^d)+{\theta }^c({\mathbf{x}}^c)+\sum _{f\in F}{\theta }_f^m({\mathbf{x}}_f^d,{\mathbf{x}}_f^c)\right\}$$

$$\underset{\mathbf{x},{\mathbf{x}}^F}{\max }\left\{{\theta }^d({\mathbf{x}}^d)+{\theta }^c({\mathbf{x}}^c)+\sum _{f\in F}{\theta }_f^m({\mathbf{x}}_f^{df},{\mathbf{x}}_f^{cf})\right\}$$

$$x_i^{df}=x_i^d,\forall f\in F,i\in f^d$$

$$x_j^{cf}=x_j^c,\forall f\in F,j\in f^c$$

$$L(\mathbf{\lambda },\mathbf{x},{\mathbf{x}}^F)={\theta }^d({\mathbf{x}}^d)+{\theta }^c({\mathbf{x}}^c)+\sum _{f\in F}{\theta }_f^m({\mathbf{x}}_f^{df},{\mathbf{x}}_f^{cf})$$

$$+\sum _{f\in F}\sum _{i\in f^d}{\lambda }_{fi}^d\left(x_i^d-x_i^{df}\right)+\sum _{f\in F}\sum _{j\in f^c}{\lambda }_{fj}^c\left(x_j^c-x_j^{cf}\right)$$

$$=\left({\theta }^d({\mathbf{x}}^d)+\sum _{f\in F}\sum _{i\in f^d}{\lambda }_{fi}^d{x}_i^d\right)$$

$$+\left({\theta }^c({\mathbf{x}}^c)+\sum _{f\in F}\sum _{j\in f^c}{\lambda }_{fj}^c{x}_j^c\right)$$

$$+\sum _{f\in F}\left({\theta }_f^m({\mathbf{x}}_f^{df},{\mathbf{x}}_f^{cf})-\sum _{i\in f^d}{\lambda }_{fi}^d{x}_i^{df}-\sum _{j\in f^c}{\lambda }_{fj}^c{x}_j^{cf}\right).$$

$$\underset{\mathbf{x},{\mathbf{x}}^F}{\max }L(\mathbf{\lambda },\mathbf{x},{\mathbf{x}}^F)$$

$$x_i^{df}=x_i^d,\forall f\in F,i\in f^d$$

$$x_j^{cf}=x_j^c,\forall f\in F,j\in f^c$$

$$L(\mathbf{\lambda })=\underset{\mathbf{x},{\mathbf{x}}^F}{\max }L(\mathbf{\lambda },\mathbf{x},{\mathbf{x}}^F)$$

$$=\underset{{\mathbf{x}}^d}{\max }\left({\theta }^d({\mathbf{x}}^d)+\sum _{f\in F}\sum _{i\in f^d}{\lambda }_{fi}^d{x}_i^d\right)$$

$$+\underset{{\mathbf{x}}^c}{\max }\left({\theta }^c({\mathbf{x}}^c)+\sum _{f\in F}\sum _{j\in f^c}{\lambda }_{fj}^c{x}_j^c\right)$$

$$+\sum _{f\in F}\underset{{\mathbf{x}}_f^{df},{\mathbf{x}}_f^{cf}}{\max }\left({\theta }_f^m({\mathbf{x}}_f^{df},{\mathbf{x}}_f^{cf})-\sum _{i\in f^d}{\lambda }_{fi}^d{x}_i^{df}-\sum _{j\in f^c}{\lambda }_{fj}^c{x}_j^{cf}\right)$$

$$L(\mathbf{\lambda })=\underset{{\mathbf{x}}^d}{\max }\left({\theta }^d({\mathbf{x}}^d)+\sum _{f\in F}\sum _{i\in f^d}{\lambda }_{fi}^d{x}_i^d\right)$$

$$+\underset{{\mathbf{x}}^c}{\max }\left({\theta }^c({\mathbf{x}}^c)+\sum _{f\in F}\sum _{j\in f^c}{\lambda }_{fj}^c{x}_j^c\right)$$

$$+\sum _{f\in F}\underset{{\mathbf{x}}_f^d,{\mathbf{x}}_f^c}{\max }\left({\theta }_f^m({\mathbf{x}}_f^d,{\mathbf{x}}_f^c)-\sum _{i\in f^d}{\lambda }_{fi}^d{x}_i^d-\sum _{j\in f^c}{\lambda }_{fj}^c{x}_j^c\right)$$

$${\bar{\theta }}_d^{\mathbf{\lambda }}({\mathbf{x}}^d)$$

$$={\theta }^d({\mathbf{x}}^d)+\sum _{f\in F}{\mathbf{\lambda }}_f^d{\mathbf{x}}_{|f}^d$$

$${\bar{\theta }}_c^{\mathbf{\lambda }}({\mathbf{x}}^c)$$

$$={\theta }^c({\mathbf{x}}^c)+\sum _{f\in F}{\mathbf{\lambda }}_f^c{\mathbf{x}}_{|f}^c$$

$${\bar{\theta }}_f^{\mathbf{\lambda }}({\mathbf{x}}_f^d,{\mathbf{x}}_f^c)$$

$$={\theta }_f^m({\mathbf{x}}_f^d,{\mathbf{x}}_f^c)-{\mathbf{\lambda }}_f^d{\mathbf{x}}_f^d-{\mathbf{\lambda }}_f^c{\mathbf{x}}_f^c$$

$$L(\mathbf{\lambda })=\underset{{\mathbf{x}}^d}{\max }{\bar{\theta }}_d^{\mathbf{\lambda }}({\mathbf{x}}^d)+\underset{{\mathbf{x}}^c}{\max }{\bar{\theta }}_c^{\mathbf{\lambda }}({\mathbf{x}}^c)+\sum _{f\in F}\underset{{\mathbf{x}}_f^d,{\mathbf{x}}_f^c}{\max }{\bar{\theta }}_f^{\mathbf{\lambda }}({\mathbf{x}}_f^d,{\mathbf{x}}_f^c).$$

$$\text{MAP}(\mathbf{\theta })\le L(\mathbf{\lambda }).$$

$$\underset{\mathbf{\lambda }}{\min }L(\mathbf{\lambda })$$

$$s^d(\mathbf{\lambda })$$

$$s^c(\mathbf{\lambda })$$

$$s^f(\mathbf{\lambda })$$

$$L(\mathbf{\lambda })=s^d(\mathbf{\lambda })+s^c(\mathbf{\lambda })+\sum _{f\in F}{s}^f(\mathbf{\lambda }).$$

Full text

1 Introduction

Bayesian optimization (BO) movckus1975 is a well-established 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 shahriari16. More recently, problems with purely discrete domains, such as food safety control and model-sparsification in multi-component systems baptista2018 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 krause2008, drug discovery negoescu2011, optimizer configuration hutter2011 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.

Related Work.

Extending continuous BO methods shahriari16 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 baptista2018; oh2019 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 hutter2011 uses a random forest surrogate model. However, its frequentist uncertainty estimates may be too inaccurate to steer the sampling. TPE bergstra2011a 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) li2016 uses cheap but less accurate approximations of the objective to dynamically allocate resources for function evaluations. BOHB falkner2018 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. garrido2018 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 hernandez2015c; gardner2014; sui2015safe 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 mixed-variable 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-the-art methods and is competitive with human experts (Sec. 4).

2 Problem Statement

We consider the problem of optimizing an unknown, costly-to-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

[TABLE]

where $\mathcal{X}\subseteq\mathcal{X}^{c}\times\mathcal{X}^{d}$ with continuous subspace $\mathcal{X}^{c}$ and discrete subspace $\mathcal{X}^{d}$. Both constraints $g^{c}(\mathbf{x})\geq 0$ over $\mathcal{X}^{c}$ and $g^{d}(\mathbf{x})\geq 0$ over $\mathcal{X}^{d}$ are known, and specifically $g^{d}(\mathbf{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, $\mathcal{X}^{c}=[0,1]^{D_{c}}$. We further assume, w.l.o.g., that the discrete inputs are binary, i.e., vectors $\mathbf{x}^{d}\in\mathcal{X}^{d}=\{0,1\}^{D_{d}}$ are vertices of the unit hypercube. This representation can effectively capture the domain of any discrete function. For example, a vector $\mathbf{x}^{d}=[x^{d}_{i}]_{i=1}^{D_{d}}\in\mathcal{X}^{d}$ can encode a subset $A$ of a ground set of $D_{d}$ elements, such that $x^{d}_{i}=1\Leftrightarrow a_{i}\in A$ and $x^{d}_{i}=0\Leftrightarrow a_{i}\notin A$, yielding a set function. Alternatively, $\mathbf{x}^{d}\in\mathcal{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 $\mathbf{x}_{t}\in\mathcal{X}$ and observe a noise-perturbed output $y_{t}\triangleq f(\mathbf{x}_{t})+\epsilon$ with $\epsilon\overset{\text{iid}}{\sim}\mathcal{N}(0,\beta^{-1})$, $\beta>0$. As evaluating $f$ is costly, the goal is to query inputs based on past observations to find a global minimizer $\mathbf{x}_{*}\in\operatorname{arg\,min}_{\mathbf{x}\in\mathcal{X}}f(\mathbf{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 $\alpha:\mathcal{X}\rightarrow\mathbb{R}$ that expresses the informativeness of input $\mathbf{x}$ about the location of $\mathbf{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 $\mathbf{x}_{t}\in\operatorname{arg\,max}_{\mathbf{x}\in\mathcal{X}}\alpha_{t}(\mathbf{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(\mathbf{x})=\mathbf{w}^{\top}\boldsymbol{\phi}(\mathbf{x})$, are a class of flexible parametric models that strike a good trade-off between model capacity, interpretability and ease of use through the definition of features $\boldsymbol{\phi}:\mathcal{X}\rightarrow\mathbb{R}^{M}$. 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,

[TABLE]

where, for $j\in\{d,c,m\}$, $\boldsymbol{\phi}^{j}(\mathbf{x}^{j})=[\phi^{j}_{i}(\mathbf{x}^{j})]_{i=1}^{M_{j}}\in\mathbb{R}^{M_{j}}$ and $\mathbf{w}^{j}\in\mathbb{R}^{M_{j}}$ are the feature and weight vector for the $d$iscrete, $c$ontinuous and $m$ixed 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 rolland2018; mutny2018, including deep neural network hyperparameter tuning hazan2017. Similarly, we focus on features defined over small subsets of the inputs. Formally, we consider $\boldsymbol{\phi}(\mathbf{x})=[\phi_{k}(\mathbf{x}_{k})]_{k=1}^{M}$, where $\mathbf{x}_{k}$ is a subvector of $\mathbf{x}$ containing exclusively continuous or discrete variables or a mix of both. Thus, the objective $f(\mathbf{x})$ can be decomposed into a sum of low-dimensional functions $f_{k}(\mathbf{x}_{k})\triangleq w_{k}\phi_{k}(\mathbf{x}_{k})$ defined over subspaces $\mathcal{X}_{k}\subseteq\mathcal{X}$ with $\text{dim}(\mathcal{X}_{k})\ll\text{dim}(\mathcal{X})$. This defines a generalized additive model rolland2018; hastie2017, 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 $\bar{D}_{d}\triangleq\max_{k\in[M]}\text{dim}(\mathcal{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, $\bar{D}_{c}$ and $\bar{D}_{m}$ 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\in\mathcal{O}\big{(}D_{d}^{\bar{D}_{d}}+D_{c}^{\bar{D}_{c}}+(D_{d}+D_{c})^{\bar{D}_{m}}\big{)} scales exponentially in the EDs only (as modeling up to $L$-th order interactions of $N$ inputs requires $\sum_{l=0}^{L}\binom{N}{l}\in\mathcal{O}(N^{L})$ terms), which are usually small, even if the true dimensionality is large.

Discrete Features $\boldsymbol{\phi}^{d}$.

We aim to define features $\phi^{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 $\mathcal{X}^{d}$ (or, equivalently, all elements $S$ of the powerset $2^{\mathcal{X}_{d}}$ of $\mathcal{X}_{d}$) and define a monomial $\prod_{j\in S}x^{d}_{j}$ for each subset $S$ (where for $S=\emptyset$, $\prod_{j\in\emptyset}x^{d}_{j}=1$). We then form a weighted sum of all monomials to yield the multi-linear polynomial ${\mathbf{w}^{d}}^{\top}\boldsymbol{\phi}^{d}(\mathbf{x}^{d})=\sum_{S\in 2^{\mathcal{X}_{d}}}w_{S}\prod_{j\in S}x^{d}_{j}$. This functional representation corresponds to the Fourier expansion of a pseudo-Boolean function (PBF) boros2002. 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, baptista2018; hazan2017 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, ${\mathbf{w}^{d}}^{\top}\boldsymbol{\phi}^{d}(\mathbf{x}^{d})=w_{\emptyset}+\sum_{i=1}^{n}w_{\{i\}}x^{d}_{i}+\sum_{1\leq i<j\leq n}w_{\{i,j\}}x^{d}_{i}x^{d}_{j}$, which induces the discrete feature representation $\boldsymbol{\phi}^{d}(\mathbf{x}^{d})\triangleq[1,\{x^{d}_{i}\}_{i=1}^{D_{d}},\{x^{d}_{i}x^{d}_{j}\}_{1\leq i<j\leq D_{d}}]^{\top}$ and reduces the number of model weights to $M_{d}\in\mathcal{O}(D_{d}^{2})$.

Continuous Features $\boldsymbol{\phi}^{c}$.

In BO over continuous spaces, most approaches are based on Gaussian process (GP) models williams2006gaussian 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 $\boldsymbol{\phi}^{c}(\mathbf{x}^{c})$ that result in a finite linear approximation of a GP. One simple, yet theoretically sound, choice is the class of Random Fourier Features (RFFs) rahimi2008, which use Monte Carlo integration for a randomized approximation of a GP. Alternatively, one can use Quadrature Fourier Features mutny2018, 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 jenatton2017; mutny2018. 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 $\boldsymbol{\phi}^{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. $\boldsymbol{\phi}^{m}(\mathbf{x}^{d},\mathbf{x}^{c})\triangleq[\phi^{d}_{i}(\mathbf{x}^{d})\cdot\phi^{c}_{j}(\mathbf{x}^{c})]^{\top}_{1\leq i\leq M_{d},1\leq j\leq M_{c}}$. This formulation provides a good trade-off between modeling accuracy and computational complexity. In particular, it allows us to reduce $\boldsymbol{\phi}^{m}$ to the discrete feature representation $\boldsymbol{\phi}^{d}$ when conditioned on a fixed assignment of continuous variables $\boldsymbol{\phi}^{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 $M_{d}M_{c}$ features, resulting in a total of $M=M_{d}+M_{c}+M_{d}M_{c}$. 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 $\mathbf{X}_{1:t}\in\mathbb{R}^{t\times D}$ be the matrix whose $i^{\text{th}}$ row contains the input $\mathbf{x}_{i}\in\mathcal{X}$ queried at iteration $i$, $\dim\mathcal{X}=D$, and let $\mathbf{y}_{1:t}=[y_{1},\dots,y_{t}]^{\top}\in\mathbb{R}^{t}$ be the array of the corresponding noisy function observations. Also, let $\mathbf{\Phi}_{1:t}\in\mathbb{R}^{t\times M}$ be the matrix whose $i^{\text{th}}$ row contains the featurized input $\boldsymbol{\phi}(\mathbf{x}_{i})\in\mathbb{R}^{M}$. The formulation of $f$ in Eq. (2) and the noisy observation model induce the Gaussian likelihood $p(\mathbf{y}_{1:t}|\mathbf{X}_{1:t},\mathbf{w})=\mathcal{N}(\mathbf{\Phi}_{1:t}\mathbf{w},\beta^{-1}\mathbf{I})$. To reflect our a priori belief about the weight vector $\mathbf{w}$ and thus $f$, we specify a prior distribution over $\mathbf{w}$. A natural choice for this is a zero-mean isotropic Gaussian prior $p(\mathbf{w}|\alpha)=\mathcal{N}(\mathbf{0},\alpha^{-1}\mathbf{I})$, with precision $\alpha>0$, which encourages $\mathbf{w}$ to be uniformly small, so that the final predictor is a sum of all features, each giving a small, non-zero contribution. Given the likelihood and prior, we infer the posterior $p(\mathbf{w}|\mathbf{X}_{1:t},\mathbf{y}_{1:t},\alpha,\beta)\propto p(\mathbf{y}_{1:t}|\mathbf{X}_{1:t},\mathbf{w},\beta)p(\mathbf{w}|\alpha)$, which due to conjugacy is Gaussian, $p(\mathbf{w}|\mathbf{X}_{1:t},\mathbf{y}_{1:t})=\mathcal{N}(\mathbf{m},\mathbf{S}^{-1})$, with mean $\mathbf{m}=\beta\mathbf{S}^{-1}\mathbf{\Phi}_{1:t}^{\top}\mathbf{y}_{1:t}\in\mathbb{R}^{M}$ and precision $\mathbf{S}=\alpha\mathbf{I}+\beta\mathbf{\Phi}_{1:t}^{\top}\mathbf{\Phi}_{1:t}\in\mathbb{R}^{M\times M}$ williams2006gaussian. This simple analytical treatment of the posterior distribution over $\mathbf{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) thompson1933, which samples weights $\widetilde{\mathbf{w}}\sim p(\mathbf{w}|\mathbf{X}_{1:t},\mathbf{y}_{1:t},\alpha,\beta)$ from the posterior and chooses the next input by solving $\widehat{\mathbf{x}}\in\operatorname{arg\ min}_{\mathbf{x}\in\mathcal{X}}\widetilde{\mathbf{w}}^{\top}\boldsymbol{\phi}(\mathbf{x})$. TS intuitively focuses on inputs that are plausibly optimal and has previously been successfully applied in both discrete and continuous domains baptista2018; mutny2018.

TS requires solving $\widehat{\mathbf{x}}\in\operatorname{arg\ min}_{\mathbf{x}\in\mathcal{X}}\widetilde{\mathbf{w}}_{t}^{\top}\boldsymbol{\phi}(\mathbf{x})$, which is a challenging mixed-variable optimization problem. However, as $\widetilde{\mathbf{w}}_{t}^{\top}\boldsymbol{\phi}(\mathbf{x})$ decomposes as in Eq. (2), we can naturally use an alternating optimization scheme which iterates between optimizing the discrete variables $\mathbf{x}^{d}$ conditioned on a particular setting of the continuous variables $\mathbf{x}^{c}$ 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 \widehat{\mathbf{x}}^{d}\in\operatorname{arg\ min}_{\mathbf{x}^{d}\in\mathcal{X}^{d}}\big{(}{{}\widetilde{\mathbf{w}}^{d}}^{\top}\boldsymbol{\phi}^{d}(\mathbf{x}^{d})+{{}\widetilde{\mathbf{w}}^{m}}^{\top}\boldsymbol{\phi}^{m}(\mathbf{x}^{d},\mathbf{x}^{c}=\widehat{\mathbf{x}}^{c})\big{)}, \widehat{\mathbf{x}}^{c}\in\operatorname{arg\ min}_{\mathbf{x}^{c}\in\mathcal{X}^{c}}\big{(}{{}\widetilde{\mathbf{w}}^{c}}^{\top}\boldsymbol{\phi}^{c}(\mathbf{x}^{c})+{{}\widetilde{\mathbf{w}}^{m}}^{\top}\boldsymbol{\phi}^{m}(\mathbf{x}^{d}=\widehat{\mathbf{x}}^{d},\mathbf{x}^{c})\big{)}. 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 $g^{m}(\mathbf{x})\geq 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 pseudo-Boolean function is equivalent to a binary integer quadratic program (IQP) boros2002, allowing us to exploit commonly-used efficient and robust solvers such as Gurobi or CPLEX. While solving general binary IQPs is NP-hard boros2002, these optimizers are in practice very efficient for the dimensionalities we consider (i.e., $D_{d}<100$). This approach allows us to use any functionality offered by these tools, such as the ability to optimize objectives subject to linear constraints $\mathbf{A}\mathbf{x}^{d}\leq\mathbf{b}$, $\mathbf{A}\in\mathbb{R}^{K\times D_{d}},\mathbf{b}\in\mathbb{R}^{K}$ or quadratic constraints ${\mathbf{x}^{d}}^{\top}\mathbf{Q}\mathbf{x}^{d}+\mathbf{q}^{\top}\mathbf{x}^{d}\leq b$, $\mathbf{Q}\in\mathbb{R}^{D_{d}\times D_{d}},\mathbf{q}\in\mathbb{R}^{D_{d}},b\in\mathbb{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 baptista2018; hazan2017 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 multi-armed bandits, a well-studied class of methods for solving structured decision making problems abeille2017. These also assume the objective to be linear in features $\boldsymbol{\phi}(\mathbf{x})\in\mathbb{R}^{M}$ with a fixed but unknown weight vector $\mathbf{w}\in\mathbb{R}^{M}$, i.e. $\mathbb{E}[f(\mathbf{x})|\boldsymbol{\phi}(\mathbf{x})]=\mathbf{w}^{\top}\boldsymbol{\phi}(\mathbf{x})$, and aim to minimize the total regret up to time $T$: $\mathcal{R}(T)=\sum_{t=1}^{T}(f(\mathbf{x}_{*})-f(\mathbf{x}_{t}))$. We obtain the following regret bound for MiVaBo:

Proposition 1**.**

Assume that the following assumptions hold in every iteration $t=1,\ldots,T$ of the MiVaBo algorithm:

$\widetilde{\mathbf{w}}_{t}\sim\mathcal{N}(\mathbf{m},24M\ln T\ln\frac{1}{\delta}\mathbf{S}^{-1})$, i.e. with scaled variance. 2. 2.

$\mathbf{x}_{t}=\arg\min_{\mathbf{x}}\widetilde{\mathbf{w}}^{\top}\boldsymbol{\phi}(\mathbf{x})$* is selected exactly.222To 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 sontag2011.* 3. 3.

\|\widetilde{\mathbf{w}}_{t}\|_{2}\leq c,\|\boldsymbol{\phi}(\mathbf{x}_{t})\|_{2}\leq c,\|f(\mathbf{x}_{*})-f(\mathbf{x}_{t})\|_{2}\leq c\, $c\in\mathbb{R}^{+}$.

Then, $\mathcal{R}(T)\leq\tilde{\mathcal{O}}\left(M^{3/2}\sqrt{T}\ln\frac{1}{\delta}\right)$ with probability $1-\delta$.

Prop. 1 follows from Theorem 1 in abeille2017 and works for infinite arms $\mathbf{x}\in\mathcal{X},|\mathcal{X}|=\infty$. 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, $\lim_{T\rightarrow\infty}\mathcal{R}(T)/T=0$, i.e., convergence to the global minimum, since the minimum found after $T$ iterations is no further away from $f(\mathbf{x}_{*})$ than the mean regret $\mathcal{R}(T)/T$. To our knowledge, MiVaBo is the first mixed-variable BO algorithm for which such a guarantee is known to hold.

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 MiVaBo333We provide a Python implementation of MiVaBo at https://github.com/edaxberger/mixed_variable_bo., we set the prior variance $\alpha$, observation noise variance $\beta$, and kernel bandwidth $\sigma$ 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 srinivas09, 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 vanschoren2014 contains evaluations for various machine learning methods trained on several datasets with many hyperparameter settings. We consider extreme gradient boosting (XGBoost) chen2016, 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 eggensperger2015, meaning that the objective returns the error of the closest (w.r.t. Euclidean distance) setting available in the dataset. Section 4 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) kingma2013 composed of a convolutional encoder and a deconvolutional decoder salimans2014. The VAEs are evaluated on stochastically binarized MNIST, as in burda2015, 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.

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\times 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. 4 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 model-based 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 well-performing 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. 4. 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 burda2015 (see Fig. 4). This is comparable to the performance of 82-83 nats achieved by human expert tuned models, e.g. as reported in salimans2014 (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ý, 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.