Constrained Discrete Black-Box Optimization using Mixed-Integer Programming
Authors: Theodore Papalexopoulos
Christian Tjandraatmadja
Ross Anderson
Juan Pablo Vielma
David Benjamin Belanger
Venue: ICLR 2022 Submitted (2022)
Keywords: discrete blackbox optimization
mixed integer programming

Abstract: Discrete black-box optimization problems are challenging for model-based optimization (MBO) algorithms, such as Bayesian optimization, due to the size of the search space and the need to satisfy combinatorial constraints.  In particular, these methods require repeatedly solving a complex discrete global optimization problem in the inner loop, where popular heuristic inner-loop solvers introduce approximations and are difficult to adapt to combinatorial constraints. In response, we propose NN+MILP, a general discrete MBO framework using piecewise-linear neural networks as surrogate models and mixed-integer linear programming (MILP) to optimize the acquisition function. MILP provides optimality guarantees and a versatile declarative language for domain-specific constraints. We test our approach on a range of unconstrained and constrained problems, including DNA binding and the NAS-Bench-101 neural architecture search benchmark. NN+MILP surpasses or matches the performance of algorithms tailored to the domain at hand, with global optimization of the acquisition problem running in a few minutes using only standard software packages and hardware.

Under review as a conference paper at ICLR 2022
CONSTRAINED DISCRETE BLACK-BOX OPTIMIZA- TION USING MIXED-INTEGER PROGRAMMING
Anonymous authors Paper under double-blind review
ABSTRACT
Discrete black-box optimization problems are challenging for model-basee opti- mization (MBO) algorithms, such as Bayesian optimization, due the size of the search space and the need to satisfy combinatorial constraints. In partic- ular, these methods require repeatedly solving a complex discrete global opti- mization problem in the inner loop, where popular heuristic inner-loop solvers introduce approximations and are difficult to adapt to combinatorial constraints. In response, we propose NN+MILP, general discrete MBO framework using piecewise- linear neural networks as surrogate models and mixed-integer linear programming (MILP) to optimize the acquisition function. MILP provides op- timality guarantees and a versatile declarative language for domain-specific con- straints. We test our approach lems, including DNA binding and the NAS-Bench- 101 neural architecture search benchmark. NN+MILP surpasses or matches the performance of algorithms tai- lored to the domain at hand, with global optimization of the acquisition problem running in a few minutes using only standard software packages and hardware.
range of unconstrained and constrained prob-
INTRODUCTION
The problem of optimizing expensive black-box function f R over discrete, constrained domain arises in numerous application domains, e.g. neural architecture search (Zoph & Le 2017).
program synthesis (Summers) 1977 Biermann 1978), small-molecule design (Elton et al. 2019),
and protein design (Yang et al. 2019). In such resource-constraine settings, it is desirable to de- velop algorithms that exploit known combinatorial structure in Ω search the space more efficiently.
Model-based Black-box Optimization (MBO), a popular paradigm that includes Bayesian Optimiza- tion as special case, iteratively refines function approximator fxx = f(x) and selects new points to query by optimizing an acquisition function axx) derived from a point estimate or posterior dis- tribution over f (Section 2.1) This inner-loop optimization problem assumed to be easier than the original, since, for example, a(x) is less expensive to query than f(x) provides "white-box" properties such as gradients.
There is vast literature addressing the challenges of applying MBO in practice. We focus two of these: first, optimizing a(x) may itself be a compurationally-diffcall optimization problem; second, in many applications, practitioners neural architecture search, x might represent a computation graph that must be both connected and acyclic. Due the difficulty in optimizing the acquisition function over combinatorial domain, most approaches resort to heuristic inner-loop solvers, which often need to be specialized to the problem at hand ensure feasibility, e.g. evolutionary solvers with custom mutation operators.
confronted by additional constraints on x. For example, in
To address the challenge of inner-loop optimization, we introduce a new MBO framework, NN+MILP, that exactly solves the acquisition problem using mixed-integer linear programming (MILP). Crucially, by framing the inner-loop optimization as an MILP, our approach can flexibly incorporate wide variety of logical, combinatorial, and polyhedral constraints, which need only be provided in a declarative sense. Using MILP in the inner loop does restrict the functional form of (or the acquisition function based on it), but it supports any piecewise linear function. In particular, we employ the class of neural network (NN) approximators with ReLU activation functions due to their scalability and accuracy in practice, and because we can draw on recent work improving the

---

Under review as a conference paper at ICLR 2022
performance of MILP for optimizing such NNs with respect to their inputs (Anderson al. 2020
blxck-box For us, MILP is practical to use in the inner loop because the dimensionality of typical black box op- timization problems is orders of magnitude smaller than those usually considered by MILP solvers. Our contributions are as follows:
We introduce NN+MILP, an MBO framework for discrete black-box problems with NN surro- gates and exact optimality guarantees for solving the acquisition problem.
We observe in our experiments that the runtime of MILP is practical for with black-box prob- lems of real-world scale, often solving the inner acquisition problem in seconds using standard packages and hardware. Moreover, MILP provides simple declarative language for problem- specific constraints.
We show that MN+MILP matches or surpasses the performance of strong MBO baselines based on croblem-specifee evolutionary algorithms on wide range of synthetic and real-world discrete black-box problems.
We use the novel MILP formulation of its graph domain.
101 neural architecture search benchmark case study, presenting
BACKGROUND AND RELATED WORK
2.1 MODEL-BASEE BLACK-BAR BBx OPTIMIZATION
Model-based Black-box Optimization (MBO) is broad family of methods that includes Bayesian optimization as a special case (Mockus et al. 1978 Dones et al.) 1998] |Hutter et al., 2011; Snoek et al. 2012, Shahriari et al. 2015). As depicted in Algorithm| 1 the method proposes at iteration t surrogate model to approximate f using
using three steps. First, the user performs inference the data previously collected from the black-bor box function. Here, fit() may return a point estimate for f. a posterior distribution over f. or a posterior predictive distribution. Next, an acquisition function a(x) based on fax: is posed that quantifies the quality of new points to query. Finally, æi selected as the best point found by solving the acquisition problem, where an inner-loop solver (approximately) optimizes a(x). The acquisition problem is typically designed such that it is more approachable than directly solving the original problem. For example, a(z) may be orders of mag- nitude less expensive to evaluate have tractable functional form. Practitioners can encode prior knowledge about the structure of f via a choice of inductive bias for f. e.g. suitable Gaussian Process kernel or neural-network architecture.
Bayesian optimization performs Bayesian inference over f and employs an acquisition function that accounts for uncertainty in f. Doing so provides principled mechanisms for balancing exploration and exploitation (Mockus et al.) 1978) Srinivas et al. 2010) and is particularly important in early rounds of optimization when models are fit on limited data. We refer our method as an instance of MBO, not Bayesian optimization, because it does not assume formal Bayesian inference for f. Gaussian processes (GPs) are often used for f in Bayesian optimization, since they provide closed- form posterior inference, naturally adjust their expressivity the dataset grows, and users inject domain knowledge via a choice of kernel (Rasmussen & Williams) 2006 Oh et al. 2019). On the other hand, neural networks provide practical alternative (Snoek et al. 2015 Hernández-LLoato et al. 2017), since they often scale more gracefully, either computationally or statistically, to large datasets or nigh-dimensional domains.
2.2 SOLVING THE DISCRETE MBO ACQUISITION PROBLEM
In general, the inner-loop problem is itself a non-trivial global optimization problem. Prior work on discrete MBO has mainly employed local search solvers, such as evolutionary search, with lim- ited guarantees (Hutter et al.) 2011 Müller] 2016 Oh et al., 2019, [andasamy et al.} (020). A key advantage such solvers is that they treat a(x) as a black box, which provides practitioners with freedom when designing application-specific surrogate models. On the other hand, particular choices of surrogate model and acquisition function lead to acquisition problems that can be (ap- proximately) solved using specialized combinatorial solvers (Baptista & Poloczek, 2018 Deshwall et al., (020)). mixed-intege nonlinear programming (MINLP) solvers (Costa & Nannicini 2018 Kim & Boukouvala] 2020), continuous optimization solvers (Bliek et al. 2021).
2

---

Under review as a conference paper at ICLR 2022
Algorithm 1 MBO
Algorithm 2 NN+MILP
Input: hypothesis class F, budget N, initial dataset D {ai f(x)) op-
Input: hypothesis class F, budget N, initial dataset D, {z,, f(z) MILP domain formulation Ma
timization domain !
for t = n 1 N do
for t = 1 do
P(f.) t Eit(F,D.-))
ft t- fii((FD(--))
3.2
a(x) get acquis func(P(fi)
M, build mi 1p(pf Ma. D,-1)
3.3
3t - inner sol v)
Jt + opt.ppimmii.. meeimiieeMMe (Me) (generic MILP solver)
D, < D. U ( f(x))
D. < D.-1U U fae. / (xx)
end for return arg (1) Ve
end for return arg (ax 1/2
Therefore, practitioners must decide between either introducing difficult-to-ane due to inexact heuristic solvers or using tractable surrogate models that may be mis-specifies for the application domain. This serves key motivation for our work: we seek enable practition- ers to employ broad families of surrogate models and exactly solve the acquisition problem with reasonable computational overhead in practice.
approximations
2.3 CONSTRAINED MBO
In many applications, is subject to non-trivial structural constraints. Prior work has largely focused on the case where determining whether e is feasible requires evaluating an expensive, perhaps noisy,
black-box function h(x) with cost comparable to f(x)) (Schonlau et al. 1998 Gelbart et al. 2014
Hermnnezz-Lobal et al. 2016 Ariafar et al. 2019 Letham et al. 2019 Here, standard acquisition
functions can be extended to account for an additional classifier h(x) trained to predict h(z).
Problems with inexpensive white-box h(z) can be tackled using these approaches for black-box con- straints, but doing so may lead to slower optimization and may query f(x) at invalid æ, which can be unsafe when performing physical experiments (Berkenkamp al. 2016). Instead, the inner-loop solver can be modified directly to guarantee feasibility, e.g., by using rejection sampling (Shi et al.) 2020 Kandasamy et al. 2020). If using local search algorithms, the solver would need to be cus- tomized for each family of constraints, a task usually left to the user. Prior work employing MINLP solvers addresses white-bon boxte-box&ae constraints either by adding penalty for constraint violation (Costa & Nannicini 2018) or in small-scale settings (Kim & Boukouvala 2020).
2.4 MIXED INTEGER LINEAR PROGRAMMING
Mixed Integer Linear Programming (MILP) seeks to maximize a linear function over a set of deci- sion variables, some of which may be integral, subject linear inequality constraints. Decades of development have allowed MILP to have a significant impact in a wide range of applications due to its bettee-than-expeeceed computational performance (Jünger et al. 2010). Indeed, while MILP prob- lems are computationally hard (NP-complete), they are routinely solved (to global or near-global optimality) in production environments thanks to of-te-of-the-ar solvers that nearly double their machine-independent performance every year (Achterberg & Wunderling 2013 Bixby 2012).
A notable aspect of MILP that provides a simple yet extremely versatile declarative language for white-box constraints. It is well known that linear inequalities over integer variables can be used easily build pure-integer formulations for logical constraints and combinatorial optimization problems (Williams) 2013) |Schrijver 2003: Wolsey & Nemhauser] 1999). In addition, using both integer and continuous variables leads to mixed-integer formulations that can combine polyhedral
and logical constraints (Jeroslow 1989 Pochet & Wolsey 2006 Vielma 2015).
Particularly interesting to our proposed approach are MILP formulations for piecewise-lined func- tions (Huchette & Vielmal 2019 Vielma et al. 2010). Specifically, our work leverages MILP formu- activation functions such as ReLUs (An- derson et al.) 2020). Optimizing over trained ReLU networks with MILP has been done in contexts such as neural network verification (Cheng et al. 2017] Lomuscio & Maganti [007]] Tieng et al.) 2019), reinforcement learning (Ryu al. 2020
lations for trained neural networks with piecewise-lings
Delarue et al. 2020), and analysis and exact com-
3

---

Under review as a conference paper at ICLR 2022
pression of neural networks (Serra et al. 2018 2021 In particular, MILP has also been used to optimize ReLU network surrogates of simulation based constraints (Grimstad & Andersson 2019). although their approach optimizes a single surrogate model once, unlike in ours.
MILP FOR MBO
We propose the MN+MILP framework (Algorithm 2). which uses neural network surrogate models and solves the acquisition problem using MILP at every step. This provides practitioners with the flexibility to use wide variety of models and leverage MILP's versatile declarative language in- corporate constraints. This section describes various design choices to make the approach practical.
3.1 PROBLEM SETTING
Our goal to find:
x² = arg max f()), sE
(1)
where f : 12 IR is expensive. noiseless black-box reward function and 52 S 9i x... Ωn is
domain on n decision variables. We assume fi can be described by an inexpensive function (2(r) indicating whether
in Ω. Algorithms are allowed a fixed budget of N sequential queries f- X := {x; refers to the set of sampled points by iteration 6, and D; := { M/ f(x;) 1 includes corresponding rewards. An algorithm's 's performance is measured as the best reward DN- Since is noiseless, it is advantageous for algorithms to avoid repeated evaluations of the same
We choose to focus on finite discrete sets Ω as we believe this is the area where MILP can pro- vide the greatest benefit. As noted in Section 2.4 there are many well-studied formulation tech- niques for with combinatorial structure, such as directed graphs. More generally, such sets have polynomially-sizee MILP formulation whenever ha(x)) can be evaluated in polynomial time (e.g., Yannakakis (199T). Continuous and mixed-integer domains could be incorporated in our approach with some modifications (Section (6, although they are outside the scope of this paper.
3.2 SURROGATE MODEL AND ACQUISITION FUNCTION
For surrogate model f. we allow any feedforward neural network with piecewise-linee activation functions, as they can be represented by MILP (Section (.3). Though we focus on fully-connected networks using ReLUs, we note that wide range of architectures (e.g., those including convolu- tional and max-pooling layers) are piecewise-liness linear and may place suitable inductive bias on f.
Our experiments employ a simple acquisition function loosely motivated by Thompson sam- pling (Tommsson] 1933 Hernández-Lobato et al., 2017 Kandasamy et al. 2018). At each opti-
mization step, we train new regressor ((x) from scratch and set a(x) f(x). This relies on stochastic gradient descent training and random parameter initialization to increase the variability in surrogate models across iterations (Lakshminarayanan alternative acquisition functions in Section|6) We use a flattened one-hot encoding of x for the input layer, and train each network fe E F on Dr-1 using & loss. Finally, we re-scale the observed rewards D,-1 before training models, to aid both training and optimization. Poorly-scaled data may result in slower performance small inaccuracies in MILP solvers (Miltenberger et al. 2018).
al. (017). We discuss extensions to
3.3 MILP FORMULATION THE ACQUISITION PROBLEM
The inner-loop solver then seeks to find
f((x))
arg max XE
where Ω is the feasible set for 11) and X,_ is the of points where the noiseless f(x) has been queried already. The MILP formulation (2) is denoted by M2 and has the following three com- ponents:
Domain We use a hot encoding of decision variables x (unless they are already binary), defin- ing the binary decision vector with zij 1fa j} for [n],] E 2r and subject to linear
4

---

Under review as a conference paper at ICLR 2022
constraints Σ JEn, Zij = vi. Integer domains with small range may be one-hot encoded; Appendix D for a comparison between integer and one-hot encodings. Additional constraints due to St are added as necessary, with form dependent on the application hand. We assume that these are MILP-rrereed relp-rrepreeet.ale, which noted in Section 2.4 could include a wide range of combinatorial, logical, and polyhedral constraints. We use Ma to denote the domain formulation itself.
No-good Constraints A no-good constraint is that eliminates from the domain. Here, we leverage the binary nature of & to exactly eliminate the set X; from Mr. For illustrative purposes, consider single point E 5i wish to exclude from the acquisition problem's 's domain, and let ž denote its one-hot encoding (or 7 itself the problem is binary). Then the constraint:
undesirable solutions
Σ 2 1. :,y
(1-2)2)
(3)
J 5,=0
enforces that any feasible has a Hamming distance of at least from Σ. As are binary, this effec- tively eliminates just the single point from the feasible region. We therefore formulate \ &2-1 by including one such constraint for each ō E 1 Note that the right-hand side can be tightened to for one-hot encodings, and these compact no-good constraints do not extend naturally to continuous (Section (
Neural Network We formulate the neural network by introducing auxiliary decision variables en- coding the activation of each neuron for given We present here the formulation for single ReLU, commonly used throughout the literature (Section[.4). while noting that the full formulation is obtained by combining all ReLU formulations and matching their input and output variables cording to the structure of the network. The overall MILP objective is the activation corresponding the regressor's output neuron.
A ReLU neuron with vector input I and scalar output y/ has the piecewise-lineal form max(0, Tx b). where w and b are its weights and bias respectively. At optimization time, w and b are fixed, while æ and y are represented by decision variables (also used as the inputs and outputs of other ReLUs according to the feedforward structure). To handle the non-linearity, we add binary decision variable that indicates whether the ReLU is active or not. We then write the following set of constraints to enforce that y = when a 0 and y = w Jx b when a 1:
S S Ma
(4)
W x +b b < w bb+ +
(5)
where M is sufficiently large fixed value, such as an upper bound on the range of y. As w and are fixed, values for M can computed in advance of the optimization, e.g. by propagating bounds from Ω. Our experiments use more advanced method to compute M, detailed in Appendix A
3.4 OPTIMALITY GUARANTEES FOR MILP
The full acquisition problem formulation, denoted by M,, is passed generic MILP solver with fixed time limit. If the solver does not time out, it is guaranteed to have produced global optimum of ) and in Section 4.4 we find that our solver typically terminates within a practical time budget. Even if the solver times out, it will return the best feasible solution it found, plus an upper bound on the global optimal value. This bound can be used to evaluate the level of potential sub-optimality of the feasible solution. Note that solvers often find optimal solution before finding the upper bound that guarantees its optimality, so timing out do not imply sub-optimality. Finally, inner- loop optimality guarantees do not translate into guarantees for the overall black-box optimization, particularly when f(x does not belong to F. However, they do provide useful empirical tool for understanding the impact of exact inner-loop optimization (Section
EXPERIMENTS
This section presents experimental results wide range of discrete black-box problems, with and without combinatorial constraints. We focus primarily on analyzing the effect of global optimization of the acquisition function, by including controlled ablations of MN+MILP where the inner-loop solver is replaced by an inexact evolutionary alternative. Depending on the problem, we also include independent baselines tailored the application domain.
5

---

Under review as a conference paper at ICLR 2022
B80B(10,10)
RandoniMLP225.5
T1Bind[8,4)
Normalized Score
Normalized Score
1.00
Normalized Score
1.00
0.75
0.50
0.50
0.25
0.25
0.25
0.00
Encombie Reggvo
N8+
B0 Ve
RegEvo
NN+
Reg-Re
Ensembile RegEvo
RBFopt
Regivoo
Ensembile RegEvo
NAV WWrr
NW RegEro
Rarapt
Regivvo
8808(10,10) Example
RandomMLP/25 Example
TrBind(8 Example
Best Observee Reward
Best Observee Reward
Best Observed Reward
1.00
0.75
2.0
-2
MN MILP
0.50
M
MLP
MV-RegEvo
RN+RegEvo
NN+ReggEva RegFve
RegEvv Ensemble + RegFvo
0.25
RegEva Ensembles RegEvo
1.0
-4
Ensemble RegEvo
Opt
0.00
0.5
200
400
B0D
1000
200
800
1000
200
600
800
1000
Figure (Top) Distribution of algorithms' normalized scores (Section 4.1 on unconstrained problems split by objective class. Higher is better. NN+MILP matches or outperforms NN+RegEvo on 22/30 problems. An alternate plot where scores correspond to area under the best-observed reward curve (AUC) can be found in Appendix H (Bottom) Best observed reward function of iteration for an example problem in each class, averaged over 20 trials (bands indicate = 1sd). Dashed grey lines in the first 50 steps indicate the initial randomly sampled dataset, common to all methods except RBFOpt, which performs its own initialization.
In all experiments, we fix the surrogate model hypothesis class F to networks with a single, fully- connected hidden layer of 16 neurons. Models trained with TensorFlow (Abadi et al. 2016). using the ADAM optimizer. No hyper-parameter tuning is performed across problems. The MILP acquisition problem is solved with the Mixed-Integer Programming solver SCIP 7.0.1 (Gamrath ctt al. 2020) using default settings and a time limit of 500 seconds. We use standard CPU machines with IG RAM and < 10 cores.
4.1 BENCHMARKING TASKS
Unless otherwise stated, tasks' domains consist of discrete decision vectors of length n, with common alphabet for all elements. We consider three families of black-box objectives:
RandomMLP The output of multi-layer perceptron operating one-hot encoding of the in- put. Notably, architectures have significantly more layers/parameters than the 16-neuron networks used as surrogates by NN+MILP.
TfBind Binding strength of length-8 DNA sequence given transcription factor Barrera et al. 2016).
BBOB Non-linear function from the continuous Black-Box Optimization Benchmarking library (Hansen et al. 2009), where each coordinate has been uniformly discretized along its range. Despite the underlying continuous structure, inputs treated as unordered and categorical.
We use parentheses after the family name to denote dimensionality of problem, e.g. Ran- domMLP(10.5) refers a RandomMLP objective over discrete domain with 10 and A 5. Appendix B lists all functions considered, and provides further details on the BBOB discretization.
Algorithms are evaluated in terms of the best reward observed after 1000 queries, averaged over 20 trials per problem. Algorithms' performance is significantly influenced by the set of x that are proposed in early iterations. Therefore, to reduce variance when comparing algorithms, we initialize each of the 20 trials with different fixed dataset of 50 random points. To facilitate comparison across problems with different reward scales, the algorithms'' average final rewards are min/max normalized within each problem. That is, the best (resp. worst) on-average algorithm for given problem is assigned score one (resp. zero), and intermediate values express relative distance from these extremes. No parameer tuning was performed across problems for any algorithm.

---

Under review as a conference paper at ICLR 2022
4.2 UNCONSTRAINED OPTIMIZATION
Before considering problems with combinatorial white-box constraints, we first tackle simple prob- lems with additional constraints on the discrete domain, i.e., Ω = AA. This allows us compare against general-purpose algorithms for unconstrained discrete black-box optimization. We vary the problem sizes over 30 functions, consisting of eight RandomMLP125.5), ten BBOB(10.10) and twelve TjBind(8, 8,4) targets (Appendix BB)
NN-MILP provides an analytical tool for understanding the relative impacts of the choice of surro- gate model and whether the acquisition problem is solved to optimality. Doing so requires ablations that vary along two axes: the family of surrogate models and the inner-loop solver. Further configu- ration details are provided in Appendix
RegEvo Local evolutionary search (Real et al. 2019) using pointwise mutations of single parent
sequences and crossover recombination of two parent sequences.
NN + RegEvo An ablation of NN+MILP, with the only difference being the use of RegEvo in lieu of MILP for solving the acquisition problem. Here, the inner-loop solver is allowed 10k queries of the acquisition function batched over 100 rounds, and proposes the point it has visited with the highest acquisition function value. The surrogate model is fit exactly as in NN+MILP.
Ensemble RegEvo A re-implementation of the 'MBO' baseline from Angermueller (2020), using an ensemble of linear and random forest regressors as the surrogate, where hyper- parameters are dynamically selected at each iteration. The acquisition function is the ensemble mean and inner-loop optimization uses RegEvo.
al.
RBFOpt A competitive mixed-integer black-box optimization solver that uses the 'Radial Basis Function method' as a surrogate model (Costa & Nannicini 2018).
Figure I plots the distribution of algorithms' scores for all unconstrained problems and an example reward curve from each class. We omit RBFOpt from the BBOB problems since it proposes the integer midpoint (rounded down) as part of its initialization, which is close to optimal by design (see Appendix B.3). We observe that relative performance of algorithms varies significantly by objective family, with MNLP performing well across the board. In particular, we wish to highlight the empirical benefits of global optimization of the acquisition function, as illustrated by the improved performance of NN+MILP vs. RegRRe The only difference between the two is the former's stronger optimality guarantees when solving the acquisition problem. We observe that MMLP obtains a greater or equal score than its evolution- counterpart in 22 of the 30 problems con- sidered, and variance in its normalized scores is lower within given objective family.
The comparison of MILP and Ensemble+ Reggemo solver is also instructive. Here, the primary difference is the hypothesis class F. The strong performance of Ensemble + RegEvo on TfBind, and to a lesser extent BBOB, suggests that ensembles of linear and tree-based regressors are better suited to approximate those black-box objectives. However, the combination of single neural network surrogate and exact optimization is able to achieve comparable performance.
4.3 CONSTRAINED OPTIMIZATION
Next, problems are augmented with combinatorial constraints on the domain. We simulate fine- balance constraints in observational studies (Zubizarreta et al.) 2018 Bennett et al.) (020)) where the same number of items must be selected from given sub-populatee pobouuouultions (c.g., groups sharing a com- mon attribute). These simple, yet highly combinatorial, constraints allow for fair comparison with cvolutionary algorithms that are designed to maintain feasibility with every mutation. Appendix F includes additional results on binary quadratic optimization problems with more complex constraints (e.g. graph partitioning or quadratic assignment) from the MINLPLib benchmark (Vigerske) 2021).
We use binary alphabet A {0, 1} to indicate whether each of 100 items is selected. These have been partitioned into given subsets Si,..., S2 for some integer k and constraints enforce that the number of selected items is equal in pairs of subsets: Eiesay-1 2n1 = Eies,y zn for
B j E [k). We use equally-sized subsets, each of cardinality 2 creating 30 problems over a grid of ten objectives and three values of & E {5, 10, 25). Further details can be found in Appendix B RandomMLP(10002 functions then simulate non-lineat reward structures for a given selection. A different objective class based the Ising model is also considered in Appendix E
7

---

Under review as a conference paper at ICLR 2022
Best Observed Reward
MILP solver runtime (s)
20
942%
94.0%
Mi-HLP P
Mt-Contve
CFAR. 114--10
93.8%
Con
Exasmble Reggee
93.4%
200
400
600
800
1000
200
800
Step
Step
Oumdalive Architecture Trairing Time
(a)
(b)
(C)
Figure 2: (a) Best observed reward as a function of iteration for typical constrained problem (Sec-
tion averaged over 20 trials (bands indicate 1sd). Initial randomly sampled set of 50 points is omitted. Distribution of normalized final scores and more examples can be found in Appendix H function of iteration. Line and bands
(b) Distribution of MILP acquisition problem solve times show the median and 5th/95th percentile range over all trials of all TfBind(8,4) problems. (c) Test accuracy of algorithms' incumbent architecture as a function of cumulative training time on NAS- Bench-101, averaged over 100 trials. Bands indicate 95% confidence interval for the mean.
The following optimization approaches provide ablations to contrast declarative vs. procedural ap- proaches to handling the constraints. Further configuration details are given in Appendix C
ConEvo RegEvo with custom mutator that procedurally maintains feasibility. Paired subsets are mutated jointly, ensuring that the number of changes in each pair is the same.
NN+ConEvo An ablation of NN+MILP where ConEvo replaces MILP the inner-loop solver. The inner-loop solver is allowed 10k queries of the acquisition function, batched over 100 rounds.
NN+RejSample An ablation of NN+MILP where the inner-loop solver samples 10k feasible points uniformly-at-randomm from the domain. By using rejection sampling, the solver can lever- age the same declarative definition of the constraints in MNLMILP.
Ensemble+RegEvo Same as Section 4.2 using the constraint definition to assign negative re- wards to infeasible proposed points, which in turn affect the surrogate model.
Figure 2a plots algorithms' best observed reward as function of iteration for a typical problem. Here, we observe that MN+MILP and NonConEvo perform similarly, both benefiting from the abil- ity to model the objective with a surrogate, unlike ConEvo. The poor performance of NN+RejSample and Ensemble+RegEvo highlight the importance of exploiting combinatorial structure. Moreover, we obtain qualitatively similar results for an objective based on the Ising model (Appendix E
Finally, while we did not find the difference of MN+MILP and NN+ConEV to be statistically sig- nificant in terms of optimization performance, they differ considerably in terms of case of imple- mentation. In particular, MN+M required very few lines of extra code to add subset-equality constraints to the existing MILP formulation, and could have just as easily been extended to other, possibly interacting, MILP reppesealee constraints. Conversely, NN+ ConEvo relied on the imple- mentation of a custom mutator, tailored to the constraint structure at hand, and would likely require significant reworking if other constraints were added.
4.4 PRACTICALITY OF MILP
Figure 2b plots the distribution of MILP solve times for inner-loop optimization as function of iter- ation over all trials of all TiBind problems. Despite the computational complexity the acquisition problem, MILP finds globally optimal solutions in seconds. Averaging across all unconstrained ex- periments (Section 4.2). . the inner-loop optimization took 7.92 + 4.23s for MN+MIPP compared to 1.00 for RegReo (avg. + sd). We observed no MILP solves that exceeded the 500s time budget, and thus they were all provably optimal. As is often the case with MILP, the relationship between problem size and runtime can be unpredictable. For example, we encountered higher aver- age solve times for pover-dimensional TfBind problems than for RandomMLP and BBOB. We also explore the impact of surrogate model network size solve time in Appendix H Solve times tend to increase over time, possibly due to the increasing number of no-good constrains and the nature of surrogate models that have been fit more data.
8

---

Under review as a conference paper at ICLR 2022
5 NAS-BENCH-ION CASE STUDY
Finally, we the NAS-Bench-101 (Ying al. 2019) neural architecture search (NAS) benchmark to illustrate the power of MILP's declarative constraint language in formulating complex combinato- rial domains. The optimization domain consists of directed acyclic graphs (DAGs) representing the cell in a neural architecture. Two nodes represent the input and output, and must be connected by directed path, while the remaining nodes are each assigned be 1x1 convolution, 3x3 convolution, or 3x3 max-poo max-pooling. Edges specify the flow of activations between nodes. The objective fzz is out-of-sample image classification accuracy. More details can be found in Appendix
MILP FORMULATIONS FOR GRAPH SEARCH
We next introduce a novel MILP formulation that precisely characterizes the set of valid NAS- Bench-101 cells. We two sets of decision variables; the first set are binary and encode the upper-triangular adjacency matrix of a DAG with exactly V nodes. The second set are a one-hot binary encoding of nodes' operations. Crucially, the MILP to represent DAGs with fewer than nodes. Constraints enforce that all non-null nodes appear on path from the input to output node, and that there exists least one such path. A full formulation in terms of linear constraints appears in Appendix G
introduce new "null" operation, allowing
A second formulation, nearly identical but enforcing so-called symmetrry-breakine constraints, seeks to address isomorphisms in our representation. The fixed ordering of nodes above results in distinct representations for isomorphic graphs, which may be redundantly proposed despite the no-good constraints. In response, we add constraints to the MILP enforcing that null nodes must occur after any non-null nodes in the linear ordering of x, which removes some, but not all, isomorphisms. alternative approach uses an isomorphism-invan -nnvariaan surrogate (Wen et al. 2020), but MBO based on such a model still vulnerable to proposing isomorphic x.
5.2 RESULTS
We refer to an optimizer using the first formulation as NN+MILP and the second as NN+MILP (w/ symmetry). Otherwise, we the same configuration in Section Linear+MILP replaces the neural network surrogate with a linear model trained on D.-1 with randomization provided by bootstrapping that training data. Regularized evolution (RE) and random search (RS) baselines from Ying et al. (2019).
Figure 2c plots the out-of-sample accuracy of the proposed architecture with the highest observed validation accuracy (the "incumbent" architecture) vs. the cumulative architecture training time. NN+MILP, despite its more general design, significantly outperforms RE. Surprisingly, we observe that both methods using no symmetry constraints, MN+MILP and Linear+MILP, perform better than NN+MILP (w/ symmetry), despite optimizing over a larger search space. Finally, note that Linear+MILP outperforms NN+MILP in early iterations, but is eventually overtaken. Future work could select among compatcleeee models at each iteration.
CONCLUSION AND FUTURE WORK
In this work we propose the NN+MILP framework for discrete MBO, using neural networks with ReLU activations for surrogate modeling and MILP solve the acquisition problem. major advantage of our method is its generality, using MILP's versatile declarative constraint language to address domains that might otherwise require specialized search algorithms for inner-loop op- timization. Our experiments show that MN+MILP performs well range of discrete black-box problems with practical computational overhead using standard packages and hardware.
For future applications to continuous mixed-integer domains, the question arises as to how to best avoid redundant proposals and encourage exploration given that no-good constraints cannot be applied as stated. More complex acquisition functions could also be considered, as long as they remain rel.-reeresentable. For example, Expected Improvement defined over the posterior predictive distribution of ensemble of neural networks is piecewise-linear.

---

Under review as a conference paper at ICLR 2022
ETHICS STATEMENT
In this paper, we introduce general algorithm for constrained discrete black-box optimization that could be applied in any number of application domains, e-g. engineering system design, neural architecture search, or drug design. As with any optimization problem, is crucial to formulate the objective and constraints of the problem ensure they reflect ethical design principles. Users should understand whether optimizing for given objective might implicitly result in biased/diseriminatory solutions, or have unintended consequences for some unmodeled objective. For example, the NAS case study of Section 5 we follow common practice by maximizing out-of-sample accuracy, which might result larger networks that are environmentally costly train. To address this, users of our algorithm might want to incorporate environmental or ethical considerations in the (black-box) objective or constraints.
The authors do not have any conflicts of interest to disclose.
REPRODUCIBILITY STATEMENT
We include detailed implementation sections in the Appendix to accompany all of our experiments. In partieular, Appendix [A] describes the MILP formulation and solution techniques used for solving the NN+MILP acquisition problem. Appendix B contains descriptions of the RandomMLP, TfBind and BBOB benchmarks we consider in Sections 4.2 and 4.3 including how we selected problems and what settings for NN+MILP, as well as all baseline algorithms from Sections 442 and 4.3 (RegEvo, NN+RegEvo, Ensemble+ Heammble++gggee RegEvo, RBFOpt, ConEvo, NN+ConEve, NN+RejSample). Finally, Ap- pendix[G) contains the full MILP formulation for the NAS-Bench-101 graph search domain, as well as relevant details on the problem setting and algorithm evaluation.
we used. Appendix C lists implementation details and
REFERENCES
Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al. Tensorflow: system for large- scale machine learning. In 12th USENIX symposium on operating systems design and implemen- tution (OSDI 16), pp. 265-283, 2016.
T. Achterberg and R. Wunderling. Mixed integer programming: Analyzing 12 years of progress. In M. Jünger and G. Reinelt (eds.). Facets of Combinatorial Optimization: Festschrift for Martin Grötschel, pp. 449 491. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.
Ross Anderson, Joey Huchette, Will Ma. Christian Tjandraatmadja, and Juan Pablo Vielma. Strong mixed-integed programming formulations for trained neural networks. Mathematical Program- ming, pp. 1-37, 2020.
Christof Angermueller, David Belanger, Andreea Gane, Zelda Mariet, David Dohan, Kevin Murphy, Lucy Colwell, and D Sculley. Population-bar based black-box optimization for biological sequence design. In International Conference on Machine Learning, pp. 324-334. PMLR, 2020.
Setareh Ariafar, Jaume Coll-Font, Dana H Brooks, and Jennifer G Dy. ADMMBO: Bayesian op- timization with unknown constraints using ADMM. Journal of Machine Learning Research, 20 (123):1-26, 2019.
Ricardo Baptista and Matthias Poloczek. Bayesian optimization of combinatorial structures. In International Conference on Machine Learning, pp. 462-471. PMLR, 2018.
Luis A Barrera, Anastasia Vedenko, Jesse V Kurland, Julia M Rogers, Stephen S Gisselbrecht, Elizabeth J Rossin, Jaie Woodard, Luca Mariani, Kian Hong Kock, Sachi Inukai, al. Survey of variation in human transcription factors reveals prevalent DNA binding changes. Science, 351 (6280):1450-144 2016.
Magdalena Bennett, Juan-Pable Vielma, and Jose R. Zubizarreta. Building representative matched samples with valued treatments in large observational studies. Journal of Computational and Graphical Statistics, 29:744-757, 2020.
10

---

Under review as a conference paper at ICLR 2022
Felix Berkenkamp, Angela P Schoellig, and Andreas Krause. Safe controller optimization for quadrotors with gaussian processes. 2016 IEEE International Conference on Robotics and Automation (ICRA), pp. 491-496. IEEE, 2016.
Timo Berthold. Measuring the impact of primal heuristics. Operations Research Letters, 41(6): 611-614, 2013.
Alan W Biermann. The inference of regular LISP programs from examples. IEEE Transactions on Systems, Man, and Cybernetics, 8(8):585-600, 1978.
Robert E. Bixby. A brief history of linear and mixed- mixed-integer programming computation. Documenta Mathematica, pp. 107-121, 2012.
Laurens Bliek, Sicco Verwer, and Mathijs de Weerdt. Black-box combinatorial optimization using models with integer-valued minima. Annals of Mathematics and Artificial Intelligence, 89(7): 639-653, 2021.
Pierre Bonami, Lorenz T Biegler, Andrew R Conn, Gérard Cornuéjols, Ignacio E Grossmann, Carl D Laird, Jon Lee, Andrea Lodi, François Margot, Nicolas Sawaya, et al. An algorithmic framework for convex mixed integer nonlinear programs. Discrete Optimization, 5(2):186-204, 2008
Chih-Hong Cheng, Georg Nührenberg, and Harald Ruess. Maximum resilience of artificial neural networks. International Symposium on Automated Technology for Verification and Analysis, pp. 251-268. Springer, 2017.
Alberto Costa and Giacomo Nannicini. RBFOpt: an open-source library for black-box optimiza- tion with costly function evaluations. Mathematical Programming Computation, 10(4):597-629, 2018.
Arthur Delarue, Ross Anderson, and Christian Tjandraatmadja, Reinforcement learnnn with combi- natorial actions: An application to vehicle routing. Proceedings of the 34th Conference on Neural Information Processing Systems (NeurIPS 2020), arXiv:2010. 12001, 2020.
Aryan Deshwal, Syrine Belakaria, and Janardhan Rao Doppa. Mercer features for efficient combi- natorial Bayesian optimization. arXiv preprint arXiv:2012.07762. 2020.
Daniel C Elton, Zois Boukouvalas, Mark D Fuge, and Peter W Chung. Deep learning for molecular design a review of the state of the art. Molecular Systems Design & Engineering. 4(4):828-843 2019.
Gerald Gamrath, Daniel Anderson, Ksenia Bestuzheva, Wei-Kun Chen, Leon Eifler, Maxime Gasse, Patrick Gemander, Ambros Gleixner, Leona Gottwald, Katrin Halbig, Gregor Hendel, Christo- pher Hojny, Thorsten Koch, Pierre Le Bodic. Stephen J. Maher, Frederic Matter, Matthias Mil- tenberger, Erik Mühmer, Benjamin Müller, Marc F. Pfetsch, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Christine Tawfik, Stefan Vigerske, Fabian Wegscheider, Dieter Weninger, and Jakob Witzig. The SCIP Optimization Suite 7.0. ZIB-Report. Zuse Institut Berlin, 2020.
Michael A Gelbart, Jasper Snoek, and Ryan P Adams. Bayesian optimization with unknown con- straints. In Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, pp- 250-259, 2014.
Bjarne Grimstad and Henrik Andersson. ReLU networks surrogate models in mixed-integee linear programs. Computers & Chemical Engineering, 131:106580, 2019.
Nikolaus Hansen, Steffen Finck, Raymond Ros, and Anne Auger. perlmmerrmeter black-box opti- mization benchmarking 2009: Noiseless functions definitions. Research Report RR-6829, INRIA, 2009.
José Miguel Hemández-Lobato, Michael A Gelbart, Ryan P Adams. Matthew W Hoffman, and Zoubin Ghahramani. A general framework for constrained Bayesian optimization using information based search. Journal of Machine Learning Research, 17:5549-6601, 2016.

---

Under review as a conference paper at ICLR 2022
José Miguel Hrinândez-Lobatee James Requeima, Edward O Knapp,, and Alán Aspuru- Aspuru-Guzik. Parallel and distributed Thompson sampling for large-scale accelerated exploration of chemical space. In International Conference Machine Learning. pp. 1470-1479. PMLR, 2017.
Joey Huchette and Juan Pablo Vielma. Nonconvex piecewise linear functions: Advanced formula- tions and simple modeling tools. To appear in Operations Research, 10Xiv:770.00000 2019.
Frank Hutter, Holger H Hoos, and Kevin Leyton-Brown. Sequential model-based optimization for general algorithm configuration. In International conference on learning and intelligent optimiza- tion, pp. 507-523. Springer. 2011.
Robert G Jeroslow. Logic- decision support: Mixed integer model formulation. Elsevier, 1989.
Donald R Jones, Matthias Schonlau, and William J Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 13(4):455-492, 1998.
Michael Jünger, Thomas M. Liebling, Denis Naddef, George L. Nemhauser, William R. Pulley- blank, Gerhard Reinelt, Giovanni Rinaldi, and Laurence A. Wolsey (eds.). 50 Years of Integer Programming 1958-2008 From the Early Years to the State-of-the-Ae Springer, 2010. ISBN 978-3-540-68274-5.
Kirthevasan Kandasamy, Akshay Krishnamurthy, Jeff Schneider, and Barnabás Póczos. Parallelised Bayesian optimisation via Thompson sampling. In International Conference on Artificial Intelli- gence and Statistics, pp. 133-142. PMLR, 2018.
Kirthevasan Kandasamy, Karun Raju Vysyaraju, Willie Neiswanger, Biswajit Paria, Christopher R Collins, Jeff Schneider, Barnabas Poczos, and Eric P Xing. Tuning hyperparameters without grad students: Scalable and robust Bayesian optimisation with Dragonfly. Journal of Machine Learning Research, 21(81) ((11):1-22 2020.
Rickard Karlsson, Laurens Blick, Sicco Verwer, and Mathijs de Weerdt. Continuous surrogate-based optimization algorithms are well-suited for expensive discrete problems. In Benelux Conference on Artificial Intelligence, pp. 48 63. Springer, 2020.
Sun Hye Kim and Fani Boukouvala. Surrogate- based optimization for mixed-intege nonlinear problems. Computers & Chemical Engineering, 140:106847, 2020.
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. In Proceedings of the 31st International Conference on Neural Information Processing Systems, pp. 6405-6416, 2017.
Benjamin Letham, Brian Karrer, Guilherme Ottoni, Eytan Bakshy, et al. Constrained Bayesian optimization with noisy experiments. Bayesian Analysis, 5492)495-559, 2019.
Alessio Lomuscio and Lalit Maganti. An approach to reachability analysis for herdffr ReLU neural networks. arXiv preprint arXiv: 2017.
Matthias Miltenberger, Ted Ralphs, and Daniel E. Steffy. Exploring the numerics of branch-and-cut for mixed integer linear optimization. In Operations Research Proceedings 2017, pp. 151-157. Springer, 2018.
Jonas Mockus, Vytautas Tiesis, and Antanas Zilinskas. The application of bayesian methods for seeking the extremum. Towards global optimization, 2(117-129):2, 1978.
Juliane Müller. Miso: mixed-integer surrogate optimization framework. Optimization and Engi-
neering, 17(1):177-033, 2016.
Changyong Oh, Jakub M Tomezak, Efstratios Gavves, and Max Welling. Combinatorial Bayesian optimization using the graph Cartesian product. In Neural Information Processing Systems, 2019.
Fabian Pedregosa, Gaël Varoquaux. Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Seikit-leam: Machine learning in python. the Journal of machine Learning research, 12:2825-2830. 2011.
12

---

Under review as a conference paper at ICLR 2022
Yves Pochet and Laurence A Wolsey. Production planning by mixed integer programming. Springer Science & Business Media, 2006.
Carl E. Rasmussen and Christopher K.I. Williams. Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, USA, January 2006.
Esteban Real, Alok Aggarwal, Yanping Huang, and Quoe V Le. Regularized evolution for image classifier architecture search. In Proceedings of the AAAI conference on artificial intelligence, volume 33, pp- 4780-4789, 2019.
Moonkyung Ryu, Yinlam Chow, Ross Michael Anderson, Christian Tjandraatmadjj, and Craig Boutilier. CAQL: Continuous Action Q-Learning. In Proceedings the Eighth International Conference on Learning Representations (ICLR-20), Addis Ababa, Ethiopia, 2020.
Matthias Schonlau, William Welch, Donald Jones, et al. Global versus local search in constrained optimization of computer models. In New developments and applications in experimental design, pp. 11-25. Institute of Mathematical Statistics, 1998.
Alexander Schrijver. Combinatorial optimization: polyhedra and efficiency. Springer Science & Business Media, 2003.
Thiago Serra, Christian Tjandraatmadja, and Srikumar Ramalingam. Bounding and counting linear regions of deep neural networks. In International Conference on Machine Learning, pp. 4558- 4566. PMLR, 2018.
Thiago Serra, Abhinav Kumar, and Srikumar Ramalingam. Scaling up exact neural network com- pression by ReLU stability, arXiv preprint arXiv:2102.07804 2021.
Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas. Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104(1):
148-175. 2015.
Zhan Shi, Chirag Sakhuja, Milad Hashemi, Kevin Swersky, and Calvin Lin. Learned hard- ware/software co-design of neural accelerators. arXiv preprint arXiv:2010.02075, 2020.
Jasper Snoek, Hugo Larochelle, and Ryan Prescott Adams. Practical bayesian optimization of ma- chine learning algorithms. Advances in Neural Information Processing Systems, 2012.
Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mr Prabhat, and Ryan Adams. Scalable Bayesian optimization using deep neural networks. In International Conference on Machine Learning, pp. 2171-2180. PMLR, 2015.
Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. Gaussian process opti- mization in the bandit setting: No regret and experimental design. In Proceedings of the Interna- fional Conference on Machine Learning, 2010, 2010.
Phillip D Summers. A methodology for LISP program construction from examples. Journal of the ACM (JACM), 24(1):161-17,, 1977.
William R Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3/4):255-944 1933.
Vincent Tjeng, Kai Xiao, and Russ Tedrake. Evaluating robustness of neural networks with mixed integer programming. In International Conference on Learning Representations, 2019.
Juan Pablo Vielma. Mixed integer linear programming formulation techniques. SIAM Review, 57: 57,7, 2015.
Juan Pablo Vielma, Shabbir Ahmed, and George L Nemhauser. arable piecewise linear optimization: unifying framework and extensions. Operations Research, 58:303-315, 2010. 2010.
models for nonsep-
13

---

Under review as a conference paper at ICLR 2022
Stefan Vigerske. MINLPLib: A library of mixed-integer and continuous nonlinear programming instances, 2021. URL fattpa://Me http///wwwnnpppb.... org Accessed October 1, 2021 (git hash: 827fla2d).
Wei Wen, Hanxiao Liu, Yiran Chen, Hai Li, Gabriel Bender, and Pieter-Jan Kindermans. Neural predictor for neural architecture search. In European Conference on Computer Vision, pp. 660- 676. Springer, 2020.
H. Paul Williams. Model building in mathematical programming. John Wiley & Sons, 2013.
Laurence Wolsey and George Nemhauser. Integer and combinatorial optimization, volume 55. John Wiley & Sons, 1999.
Keyulu Xu, Mozhi Zhang, Jingling Li, Simon $ Du, Ken-ichi Kawarabayashi, and Stefanie Jegelka. How neural networks extrapolate: From feedforward to graph neural networks. In International Conference on Learning Representations, 2021.
Kevin K Yang, Zachary Wu, and Frances Arool. Machine-learning-uudee directed evolution for protein engineering. Nature methods, 16(8):687-694, 2019.
Mihalis Yannakakis. Expressing combinatorial optimization problems by linear programs. Journal of Computer and System Sciences, 43(3):441-466 1991.
Chris Ying, Aaron Klein, Eric Christiansen, Esteban Real, Kevin Murphy, and Frank Hutter. NAS- Bench- 101: Towards reproducible neural architecture search. In International Conference on Machine Learning, pp. 7105-7114, PMLR, 2019.
Barret Zoph and Quoe V Le. Neural architecture search with reinforcement learning. In Interna- tional Conference on Learning Representations, 2017.
Jose R. Zubizarreta, Cinar Kilcioglu, and Juan Pablo Vielma. designmatch: Matched samples that balanced and representative by design. R package version 0.3, 2018.
14

---

Under review as a conference paper at ICLR 2022
A STRENGTHENING THE MILP FORMULATION FOR NEURAL NETWORKS
Here we discuss more advanced techniques for formulating the neural network surrogate model in the MILP problem. Recall the ReLU formulation constraints (4) and from Section 3.3 except that we consider M separately for each constraint:
0 v S Mya
w + b S S w +b+ -a),
Here, we require nonnegative value M1 such that the right-hand side of (4) greater or equal than valid upper bound on when a = 1. Similarly, M1 must be nonnegative value such that the right-hand side of (5)) is greater or equal than zero when a = O. Therefore, we may choose Mg to be any upper bound of max w x + and M1 to be any upper bound of max -(w" x + b), where Q2' is the domain of the inputs of this ReLU, which depends on Ω. The tighter these bounds are, the better the MILP performs.
Moreover, if we find negative Mo or M,. then we may (in fact, must) replace the formulation by y 0 or y w + respectively, since in these cases the ReLU is always inactive or active for any x E Q2' This replacement must be done because the formulation assumes nonnegative My and M, for feasibility.
The simplest way to compute Mo and M, is to start from bounds in Q and propagate them via inter- val arithmetic. For example, if E [T_, U). then Mo can be set Eis,, W;U;; Eiseeo W; T. bbb and Mi to - -(inm,>0 W/ W Eiw W;U, b). However, despite being fast, the drawback of this simple approach is that it does not take into account constraints on S or one-hot and no-good constraints.
In our experiments, we compute Mo and M₁ by solving the linear programming (LP) relaxations of max w To x + b and max - To b) respectively (i.e. without integrality constraints). We remark that for neurons in the same layer these LPs have the same constraints but different objectives, and thus we may take advantage of the warm starting functionality in LP solvers. While this requires solving two LPs per neuron, taking into account the constraints from Q into the bounds often enable the overall MILP to be solved much faster.
The formulation can also be strengthened with cutting plane techniques (Anderson et al., 2020), but they are not particularly beneficial for the small network sizes considered in this paper (at most a single layer with 16 ReLUs) and thus do not add them. Future work could explore warm-starting the MILP solver using results from earlier MBO iterations.
B BENCHMARKING TASKS
This section details the black-box objective functions considered in both unconstrained (Section 4.2) and constrained (Section 4.3) experiments. Recall that all objective functions are defined over fixed- length discrete vectors of length n, with each element drawn from an alphabet A of fixed size.
B. TFBIND
The objective function is given by the binding affinity of length-8 DNA sequence particular transcription factor, characterized experimentally in the dataset described by Barrera et al. (2016). The problem size thus fixed by the application at hand, with n = and [ element corresponding to given DNA nucleotide). We mannmaxx-normalin values for each factor to the zero-one interval. We create 12 unconstrained problems (Section 4.2) using the following datasets: CRX R90W R1, CRX REF.N RI, FOXCI _EERR_ERR GFIIB REF RI, HOXD13 200133_332RREEE R1. HOXD13 REF REF.RI, DRRH44CC144 PAX4_REF RR,, POU6F2 REFFLRRER_RL SIX6_REF_R Here, the 3 fields separated by underscores rep- resent the transcription factor id, any mutations that have been made to the transcription factor, and the id of the experimental replicate used when collecting data.
(each input the binding affinity
R1, NRIH4 REF R1, PAX4 RI,
15

---

Under review as a conference paper at ICLR 2022
B.2 RANDOMMLP
The objective function is given by the output of a multi-layer perceptron (MLP) with randomly- sampled weights. Different functions are generated by varying the architecture type (described below) and random seed. All architectures employ one-hot encoding of the inputs as the first layer. Weights are sampled using the default behavior of lerass Dense iglorot.uniforms
We consider two architecture types, both utilizing more layers/parameters than the 16-neuron net- works used by MN+M (Section 4 The RandomFCC architecture uses two fully-connected lay- ers with 128 hidden units each, while the RandomCNN architecture uses two convolutional layers each with 64 hidden units each, kernel width of 13 and stride size of 1. We a linear activation function for the output and ReLU activations for all intermediate layers.
Unconstrained RandomMLP problems (Section 4.2 all have size objective functions are created by varying the architecture type (FCC or CNN) and random seed (0, 13, 42, 77). Constrained RandomMLP problems (Section 4.3) all have size 100 and |AA| 2. Thirty problems are created by varying the architecture type (FCC or CNN), random seed (0, 13, 42, 77, 100) and the number of paired subsets (5, 10, 15) in the subset-equality constraints (defined in Section 4.3
25 and |A| 5 Eight
B.3 BBOB
The objective is given by function from the continuous Black-Box Optimization Benchmarking library (Hansen et al. 2009). All BBOB functions are defined for a variable number of dimensions and the search domain is given as [-5,5] 5|"" with the global optimum centered at zero. We normalize each function's 's output range by evaluating at 30 fixed points and dividing outputs by the median absolute deviation in those points' values.
We discretize functions for our setting (Section 3.1 by defining a grid over the continuous search domain, adjusted so that the optimal solution exactly corresponds to point in the grid. Concretely, we use a fixed alphabet A = {L..... m } for all coordinates, denoting the index of one of m allowed values for that coordinate. Allowed values for each coordinate are equally-spaced points in the range [-5,55 5), except for a point lying closest to zero which is overwritten to exactly equal that value. In this way, the optimum is guaranteed to lie the discretized grid. Note that, despite the underlying continuous structure, all algorithms treat each dimension as an unordered, categorical variable.
For unconstrained BBOB problems (Section 4.2), we select a diverse set of objectives by taking two functions from each of the five categories defined by the BBOB library:
1. Separable functions: Sphere (SPHERE) and Ellipsoidal (ELLIPSOID SEPARABLE.
2. Functions with low or moderate conditioning: Attractive TIVE_SECTOR) and Step Ellipsoidal ELLIISODD))
Sector (ATTRAC-
3. Functions with high conditioning and unimodal: Diseus (DISCUS) and Bent Cigar (BENT.CIGAR).
4. Multi-modal functions with adequate global structure: Weierstrass (WEJERSTRASS) and Schaffers F7 (SCHAFFERS FF)..
5. Multi-model functions with weak global structure: Schwefel (SCHWEFEL) and Gal- lagher's Gaussian 21 -hi Peaks (GALLAGHER.2IME)
We set the dimension for all of these to 10 and discretize described above, using an alphabet of size |A| = 10 for all coordinates. We purposefully use a relatively large alphabet to ensure that the discretization does not obscure any inherent variance across given coordinate.
C BASELINE OPTIMIZATION ALGORITTMS
In this section we describe implementation and configuration details for all baseline optimization algorithms described in Section 4
16

---

Under review as a conference paper at ICLR 2022
C.I NN+MILP
For our experiments, we implement our main algorithm (Section 3) as follows: we use a fixed surrogate model hypothesis class F of networks with a single, fully-connected hidden layer of 16 neurons. Models are trained with TensorFlow (Abadi et al. 20T6), using the ADAM optimizer for 25K epochs with a batch size of 64 and no explicit regularization. We use a constant learning rate of a 0.01 and default decay parameters (81, (2) (0.9, 0.999). No hyper-parameter tuning is performed across problems. Model training is randomized due to the random example ordering of SGD training and random parameter initialization. The MILP acquisition problem is solved with the Mixed-Integer Programming solver SCIP 7.0.1 (Gamrath et al.) 2020) using default settings and time limit of 500 seconds. In order to increase the diversity of trained models, train each model from scratch at each iteration of optimization instead of fine-tuning a model from an earlier iteration.
C.2 REGEVO
We re-implement the local evolutionary search algorithm of Real et al. (2019), and extend the set of mutation operators from just pointwise mutators to also include a crossover operation that re- combines parent sequences. The algorithm proposes xx+1 by selecting two parent sequences from the existing population, recombining them and mutating them. Parents are chosen by tourna- ment selection, taking the two best samples from randomly-selected subset of size T of previously sampled points. The pool from which parents can be selected is limited to the D most recently- proposed points (referred rounds dominating the process. The selected parent sequences are recombined by copying them left-to-right, starting a pointer at one parent switching reading to the other parent with a fixed cross-over probability Pe after each copy. The resulting sequence is finally mutated by changing each position to a different token from AA with a fixed probability Pm-
the "alive population"), to avoid high-reward points from early
use RegEvo as the outer-loop optimization al- gorithm and set the tournament size to T 10, the alive population size to D 100, and the crossover/mutation probabilities to (Pe, Pm (0.1.
In the unconstrained experiments (Section 4.2
C.3 NN+REGEVO
This algorithm an ablation of MN+MILP, with the only difference being the use of RegEvo lieu of MILP to solve the acquisition problem at every iteration. A surrogate neural network fc E F is trained as in NN+MILP, and the acquisition function is a(x) = f((x). The problem of selecting Je-1 posed as a batched optimization problem and solved by RegEvo.
More concretely, at iteration f, the acquisition function is evaluated for all points in the existing population D. to generate the initial inner-loop population D. := {æ,, a(x;)) This population is iteratively extended by generating candidate proposals with RegEvo in batches of size b, and with rewards now corresponding to the value of the acquisition function rather than the original black- box function. That is, RegEvo generates b points by recombination/mutation of parents from D., which are evaluated on the acquisition function and added to the inner-loop population. The process repeats until total of B candidates have been generated, at which point the one with the highest acquisition function value (excluding any points already proposed) is proposed as Tt+1-
In the unconstrained experiments (Section 4.2 we use NN+RegEvo and set surrogate model hyper- parameters exactly as in NN+MILP (Section C.I For the inner-loop optimizer, we set the total
number of acquisition function evaluations to B optimizer's 's hyper-parameters, defined in Section (0.2, 0.01).
10, 000 and batch size to b 100. The RegEvo are set to T = 20, = 1, 000 and (pe, Pm) =
C.4 ENSEMBLE+
We recreate the MBO baseline of Angermueller et al. (2020). Here, surrogate modeling proceeds by optimizing the dyper-parameterd of a diverse set of regressor models through randomized search. Regressors are trained using the scikit meiit-lears -eaann libary (Pedregosa et al. (01T). drawing from the following model classes (randomized search parameters listed in parentheses):
17

---

Under review as a conference paper at ICLR 2022
LassoRegressor (alpha)
RidgeRegressor (alpha)
RandomnForesRRerresoon (max.depthh max featurrs, n.estimators)
LGBMRegressor (learning rate, n e.sssimators)
Each model is evaluated by explained variance score using five-fold cross validation the train- ing set. All models with a score 2 0.4 are used as an ensemble for the surrogate model, with their average prediction serving
the acquisition function. The acquisition problem is solved by batched RegEvo with a total of B = 12, 500 acquisition function evaluations and a batch size of b 25. The optimizer's syper-parameters, defined in Section C.2 and (Pe-Pm) = (0.2, 0.01).
are set to T 20, D 1, 000
We use ResembbleRRggEie in both the unconstrained (Section 4.2 and constrained (Section 4.3 experiments. In the latter case, we use the algorithm as a baseline that makes use of the declarative definition of constraints; during training of the ensemble, infeasible points assigned highly negative reward (worse than any observed). In this way. the surrogate model might be expected implicitly model infeasibility with low predictions which should be avoided by the inner-loop optimizer.
C.5 RBFOPT
RBFOpr (Costa & Nannicini 2018) is a black-box optimization solver for mixed-integer uncon- strained problems (i.e. with only bound constraints) that performs competitively with respect to other solvers of its type. It uses Radial Basis Function as surrogate model and includes a number of practical enhancements.
relies on a mixed-intege -ixeee nonlinear programming (MINLP) solver, BONMIN (Bonami et al. 2208), to optimize the inner loop problems. The MINLP solver could in theory incorporate constraints in similar fashion as in our work, although this not offered by the open-source implementation (aside from manually penalizing the objective function) and expect it to not scale as well as a MILP solver in practice since MINLP is a significantly more difficult problem class than MIL.P
We use RBFOpt for our unconstrained experiments (Section 4.2) using the open-source implemen- tation available tht / /gitbbbbb thub. com/ http:/oooooooooooeee rbt.ttt We leave all settings their defaults, including building the initial set of points. We note in particular that the API for this imple- mentation uses integer encoding for categorical variables (constraints are not supported, which precludes one-hot representation). As we note in the main text, we omit the RBFOpt results for BBOB because RBFOpt proposes the midpoint of this integer representation (rounded down) as part of its initialization, which is close to the optimal solution.
C.6 CONEVO
In Section 4.3 we introduce ConEvo, a local evolutionary search algorithm that exploits the known combinatorial structure of the subset-equality constraints considered therein. The method selects just single parent sequence (using the same tournament procedure RegEvo) and mutates it in way that guarantees feasibility of the child sequence. We do not implement recombination of multiple parent sequences since they are not likely to maintain feasibility. We described the application- specific mutator below.
Recall that the domain encodes the selection or not of each of items using a binary alphabet A = {0, 1}. The items' indices are partitioned into disjoint, equally-sized subsets Si,... S2x for some & and the constraints enforce that the number of selected items should be the same in pairs of subsets; that is:
Σ [1== 1ES2,11
11
Σ 1x:
=1
Vj E [A]
SES2
where we have used indicator notation and the original decision variables x rather the one-hot en- coding.
18

---

Under review as a conference paper at ICLR 2022
The mutator begins with single parent sequence z, assumed feasible, and is given access to the item subsets Si...., S2x- Each pair of subsets (S2g-1, S2y) is mutated concurrently to create the child sequence y, ensuring that mutations to subset counter- maaate-d by mutations to the second. Concretely, one of the two subsets is chosen randomly to be the "independent" mutatee with equal probability. We denote the selected subset I+, and the other subset in the pair by I Each position i E I of the child sequence is flipped from its parent value with some fixed probability Pm- We compute c, the net number of 0-to-I conversions in positions I+ If is positive (i.e. there were more 0-to-1 conversions than 1-to-0 conversions) then exactly indices are chosen randomly from {i E I 24 0), and also flipped in the child. If is negative, then -c indices are selected randomly from {i E I I; = 1] and flipped in the child. If è is zero, the positions in I are left unchanged the child. As a result, the subsets S2j-1 and Sa retain exactly the same number of selected items in the mutated sequence.
In the constrained experiments (Section 4.3). we use ConEvo as the outer-loop optimization al- gorithm setting the tournament size for selecting parent sequences to T 20 and the mutation probability to Pm = 0.05.
C.7 NN+CONEVO
This algorithm is an ablation of MN+MILP used in Section 4.3 with the only difference being the use of ConEvo in lieu of MILP to solve the constrained acquisition problem at every iteration. A surrogate neural network fr. E F is trained as in NN+MILP, and the acquisition function is a(æ) = ft(r). Selecting I++1 is posed as batched optimization problem and solved by ConEvo. The batched optimization procedure is exactly as described for NN+RegEvo (Section
In the constrained experiments (Section 4.3 parameters exactly as in NN+MILP (Section
we use NN+ConEvo and set surrogate model hyper-
C.I
For the inner-loop optimizer, we set the total
number of acquisition function evaluations to B
10, 000 and batch size to b = 100. The ConEvo
optimizer's hyper-parameters, defined in Section
C.6 are set to 20 and pm 0.05.
C.8 NN+REJSAMPLE
This algorithm is an ablation of NN+MILP used in Section4.3. with the only difference being the use of random search lieu of MILP to solve the constrained acquisition problem every iteration. surrogate neural network f. E FF trained as in NN+MILP, and the acquisition function is a(x) = fx)).
To select #z+1, first generate B points mniformly-at-ranee rejection sampling. That is, candidates are drawn uniformly at random from A and rejected they do satisfy the subset equality constraints (Section 4.3) The process continues until B feasible points have been generated, and the one with the highest observed acquisition function value (excluding any points already proposed) is proposed Tz+1-
from the constrained domain via
In the constrained experiments (Section 4.3 we use RejSamme and set surrogate model hyper- parameters exactly as in NN+MILP (Section 4.. We set the number of points to generate at each optimization step to B = 10, 000. Empirically, we observe rejection rate of z= 91% for the subset- equality constraints our experiments.
D BINARY VS INTEGER VARIABLES
In this work, we focus on problems with binary domain formulations (e.g. one-hot encoding of categorical domains), and even problems with integer variables such as the discretized BBOB are binarized. Part of the reason is to allow no-good constraints as described in Section 3.3 but in addition we have experimentally observed that the method performs better when using a binary encoding instead of an integer one.
When running this algorithm for unconstrained (bounded) integer or continuous problems, we have informally observed that method frequently proposes solutions where several variable values at either their lower bound or upper bound. As result, our method would underexplore solutions away from the boundary. A possible explanation for this is that feedforward ReL.U networks tend
19

---

Under review as a conference paper at ICLR 2022
to extrapolate linearly, and thus their optima may often lie the boundary (Xu al. 2021). In
contrast, every feasible point of a binary problem lies on a corner of the 0-1 hypercube. sim- ilar observation has been made in the context of IDONE (Blick et al., 202T). which also uses
ReI.U-hasee surrogate model: encoding the Rosenbrock problem using binary variables improves the performance of the IDONE algorithm, although the opposite happens for Bayesian optimiza- tion algorithm (Karlsson et al. 2020).
We provide computational evidence of this behavior in Figure for TYBind and BBOB instances.
The binary variables are encoded as one-hot variables, whereas the integer variables follow an arbi- trary ordering for TfBind and the problem ordering for BBOB.
(8Bind(8,4)
BBOB((0,10)
1.0
Best Reward
0.8
Best Reward
-1
0.6
-2
0.4
MILM Binary MN+MILP Integer
MN+M Binary
N Integer
0.2
200
400
600
800
1000
200
400
600
800
1000
Step
Step
Figure 3: Best observed reward as a function of iteration for all TfBind (left) and BBOB (right) instances, comparing the use of binary and integer variables. For TfBind, the categorical variables are transformed to integer with an arbitrary ordering, and for BBOB, we use the given ordering of the problem. Note that the error region is large here since we aggregate all of the instances of each class.
E CONSTRAINED OPTIMIZATION WITH AN ISING PROBLEM
To provide further evidence that the results from Section 4.3 can generalize beyond the RandomMLPP function, Figure 4 shows the result for the same set of constrained experiments with the only dif- ference that we switch the RandomMLP objective by another random function based on the Ising model. Here, the Ising model represents fully-connected graph with nodes that can take binary values. Each edge is defined by a 2 x 2 table of scores for each possible configuration of the nodes that are connected by the edge. Each value in the tables is drawn from normal distribution and the overall function is the sum of the scores over all edges. We aggregate over 10 random instances and 20 replications with different random seeds. We omit Ensemble+Reggxv with penalty described in Section 4.3 since it does not perform well with these constraints. We observe that the results are qualitatively similar to those Section 4.3 NN+MILP continues to perform significantly better than ConEvo, while NN+MILP and NNnConEvo have similar performance.
We also experiment with the same class of constrained problems but defined over larger domains (n = 200 and 400 binary decision variables). Given the larger scale, we set single value for the number of paired subsets the constraints, namely k 20 when
200 or k 40 when = 400 (in other experiments, different values of l: resulted in qualitatively similar results), and we not run NN+RejSample. We aggregate over random instances and 10 replications for each problem size. In Figure 5 we observe that NN+MILP significantly outperforms NN+ConEvo for the n = 400 instances, suggesting that beneficial to solve the acquisition problem to optimality larger scales.
Finally, in Figure we examine the impact the surrogate model's 's capacity in this larger-scale setting. We include ablations of both NN+MILP and NN+ConEvo where the surrogate model has 32 neurons in the hidden layer, instead the 16 used before. Despite optimizing over a larger domain, neither NN+MILE nor NN+ConEVo show substantial improvements in performance when using larger
20

---

Under review as a conference paper at ICLR 2022
1000))
1.0
Normalizes allizd Score
0.8
Best Reward
0.6
0.4
400
NN+MIL MILP
0.2
No
200
RejRejSample ConEvo
0.0
NN+MILP
ConEvo
ConEvo
NN+ RejSammea
200
400
600
800
1000
Step
Figure 4: (Left) Best observed reward function of iteration for random instances based the Ising model with the subset equality constraints described Section 4.3) (Right) Distribution of the algorithms' normalized scores on the same constrained Ising model problems. Higher is better.
surrogate network. This suggests that, in this case at least, the relatively small number of training points is more significant bottleneck for approximation than the capacity of the surrogate.
Ising (200,,) (200,2)
lsingModel(400,2
MN+M
2000
5000
NN-Con
Best Reward
Best Reward
4000
ConEvo
1500
3000
1000
2000
NN-MILP
500
Evv ConEva
1000
200
400
600
800
1000
200
400
600
800
1000
Step
Step
Figure 5: Best observed reward as function of iteration for random instances based on the Ising model with 200 (left) and 400 (right) binary variables and the subset equality constraints described in Section 4.3 except that the number of pairs of subsets is k 20 (left) and k: = 40 (right).
F CONSTRAINED BINARY QUADRATIC PROBLEMS FROM MINLPLIB
We study the performance of our method for linearly cenerrlycon binary quadratic problems from the MINLPLib benchmark library (Vigerske 2021). In practice, one would use a specialized mixed- integer quadratic programming solver to tackle these problems, but they serve well as bench- mark for our black-box optimization method since they are still harder to solve than MILP and include practically- motiivated constraints, along with offering good feasible solutions to compare with. Many of these instances are of slightly larger scale than typical black-box optimization problems (i.e. few hundreds of variables), which helps evaluate the method at larger scale and observe its scalability limitations.
We select the instances of type "BQP" from MINLPLib with at least constraint. We discard the 11 instances prefixed by c subb, color lab, and max csp, which are large instances for which our method was unable to find solution with primal gap at most 10% (see below). This
21

---

Under review as a conference paper at ICLR 2022
IsingModel (200,2)
(2000)
Model (400,2)
MN-MILP FC116)
2000
5000
MN--HILP FC232)
Best Reward
1500
Best Reward
4000
ConEvo
ConEvo- FC(82)
3000
1000
-N-NILP-
2000
CNnvon FC116)
500
NN-ConEve FC132)
NN-NILP FCE22
1000
200
400
600
800
1000
200
400
600
800
1000
Step
Step
Figure 6: Best observed reward as a function of iteration for random instances based on the Ising model with 200 (left) and 400 (right) binary variables and the subset equality constraints described in Section 4.3 except that the number of pairs of subsets is restricted to k = 10. FC(16) and FC(32) represent the runs where the surrogate neural network has single layer of 16 and 32 ReL.Us respectively.
leaves us with the 50 instances listed in Table | For consistency with the remainder of the paper, turn the problems into maximization problems by negating the objective function.
We run NN+MILP with the same settings as previous experiments (see Appendix C.T One dif- ference here is that the feasible set may be too small to sample from using rejection sampling. Therefore, we build our initial set of 50 points by randomly choosing an objective direction and solving MILP under the constraints of the problem, which is practically feasible since the scale of these problems is small in the context of MILP. In our experiments, we compare with the best known primal feasible solution from the MINLPLib benchmark itself (as of October 1, 2021).
Out of the 20 runs for each instance, we show in Table | the number of runs that found solution in 1000 steps with value equal to the best known value from MINLPLib, or with a primal gap of at most 1 % or 10% with respect to that value. We use the definition of primal gap from Berthold (2013): v* is the best known value and the objective value given by our method, then the primal gap is |e-v"| max(jul,l=-)>> or if [o| = |00 = 0, or if and v" have different signs.
Interestingly, in 20 out of the 50 instances, we match the best known objective value in at least one of the runs. This includes number of instances that considered to be large for black- box optimization, such as the general quadratic assignment problem pb302095 which has 600 variables.
On the other hand, we also observe that the method has difficulties in finding a good solution for larger instances. This is more clearly illustrated by Figure in which we observe how the method scales with the graph partitioning problems denoted by graphpart cl que..e. For the smaller instance (with 60 variables), our method finds an optimal solution in relatively few steps, difficulties in reaching the best known solution for the larger instance (with 180 variables) with the same constraint structure, That said, this does not mean that this method cannot be extended to scale further (e.g. 5. have not attempted to tune the parameters or change the architecture of the surrogate model for larger instances).
has
G NAS-BENCH-101 CASE STUDY
G. BACKGROUND
In Section 5 we consider the NAS-Bench-101 (Ying al., 2019) neural architecture search (NAS) benchmark as a case study, to illustrate the power of MILP's declarative constraint language in modeling complex combinatorial domains. The optimization domain consists of directed acyelic
22

---

Under review as a conference paper at ICLR 2022
Table 1: Results for subset of binary quadratic problems from the MINLPLib benchmark, indi- cating the number of runs out of 20 for which the primal gap with respect the best known primal feasible solution is at most 0%, 1,, and 10%.
Number of runs with solution at
Instance name
# variables
0% gap
1%
gap
10% gap
cardqp_iniy
50
5
14
20
cardqp-iqp
50
14
20
crossdock.15x7
210
0
0
20
crossdock 15x8
240
0
20
graphpart 2g-0044-1601
48
17
17
20
graphpart 20-0055-02
75
19
graphpaat.22--0066-00666
108
17
graphpputt22g00077000777
147
0
0
7
graphpart 2g-0088-0088
192
7
paapppart.2g-0099-222
243
0
2
graphpart2g--010--622
300
0
graphhaart22pm-00044-004
48
20
20
20
graphpaat22pn-0055-0055
75
13
13
graphpaatt2pmn-006660006666
108
15
graphpatt22mm0007700777
147
7
graphpart 20m-0088-0888
192
0
5
graphpart 22000999999
243
3
graphpart.3gg0033440234
72
18
graphpart_33g002444-224
96
17
sraphpart_3g--033--333
81
19
graphpart.3g-00334-0333
108
O
16
graphpart.3g-0344-044
144
8
"ruphpatt.g--04444400
192
0
0
8
graphpart
3pm-0024-24
72
graphpart-3pm-02244-02444
96
0
17
graphpart33pm-0333--:33.
81
1
graphpaatt.3pm-03344--23
108
0
0
11
graphpart33pm-03444-5444
144
0
0
5
graphpart.3pm-04444-04
192
6
graphpart.clique-20
60
20
20
20
graphpart..lique-30
90
20
20
20
graphpart_clique-40
120
13
13
18
graphpart_clique-50
150
0
graphpart clique-60
1 180
0
0
graphpart 70
210
0
pb302035
600
0
0
0
pb302055
600
20
pb302075
600
5
5
20
pb302095
13
20
20
pb351535
525
15
pb351555
525
1
20
pb351575
525
10
20
pb351595
525
6
19
20
qap
225
7
(Spp.0_10.00111001
180
5
20
0spp.0.11.0.1110.1
220
20
20
qspp.0.11.0011110.1
264
13
13
20
qspp 0 13.00.1111
312
19
20
qspp 0 1.100..1
364
0
16
20
q5pp.0.15.0.1110..1
420
12
13
20
23

---

Under review as a conference paper at ICLR 2022
graphpart_cliquez
graphpart clique-40
graphpart dllque-60
Best Observed Reward
-150
Best Observed Reward
-1000
Best Observed Reward
-4000
-1500
-250
2000
NN MILP
-8000
Best Known Value
-2500
-300
NN-MILP Optimal
NN++IIPP Optimal
-10000
-350
0
200
400
600
800
1000
200
400
600
800
1000
200
400
6.00
800
1000
Step
Step
Step
Figure 7: Best observed reward as function of iteration for three graph partitioning instances from MINLPLib (negated for maximization), with 60, 120, and 180 binary variables respectively. Black lines show the best known feasible solution to the problem (as of October I, 2021). Colored lines show the average over 20 trials, while bands indicate +1 sd. Note that bands that exceed the black line are an artifact of the symmetric nature of standard deviation, and do not necessarily mean trial found an improved solution.
graphs (DAGs) with a maximum of V = 7 nodes and M = 9 edges, representing the cell in a neural architecture. The overall model is obtained by stacking multiple copies of the cell. Two nodes represent the input and output, and must be connected by directed path, while the remaining nodes are each assigned to be Ixl convolution, 3x3 convolution, or 3x3 max-pooling. Edges specify the flow of activations between nodes. The goal is to find the cell architecture that maximizes out-of- sample accuracy on a given image classification task.
NAS differs from the problem setting in Section 3.1 in three key ways. First, algorithms do not have access to the true objective (out-of-sample accuracy), but instead a correlated proxy (validation accuracy). Second, far) is noisy due to the stochasticity of classifier training, and thus algorithms may benefit from repeated queries of the same point. Finally, algorithms may benefit by leveraging the validation accuracy at early epochs as a proxy halt unpromising evaluations.
Despite these differences, we apply MN+MILP exactly as described in Section 3 Appendix C.I (including no-good constraints which prevent repeated queries). We reimplement the Regularized Evolution (RE) and random search (RS) baselines from the original NAS-Bench-10 paper (Ying al. 2 using the same hyper-parameterr settings specified therein.
The NAS-Bench-10 101 dataset contains pre-computed validation and test accuracies for three inde- pendently trained replications of each architecture, as well as the training time of each. To simulate NAS, algorithms' observed reward after proposing an architecture is the validation accuracy of a randomly sampled replication from said architecture. This defines the notion of an "incumbent" proposal, namely the proposed architecture with the highest (observed) validation accuracy, which may not in fact be the best (unobserved) test accuracy. Instead of allowing algorithms a fixed bud- get of evaluations, the objective until cumulative training time exceeds the budget. For evaluation purposes (c.g. Fig- ure 2c) we plot the out-of-sampp sample accuracy of the incumbent architecture as a function of cumulative architecture training time.
fixed budget of T
10 seconds, and allow algorithms to query
G.2 DOMAIN FORMULATION
To formulate the MAS-Bench-101 domain, we first define a representation of cell architectures as fixed-length binary vectors. We split the representation into two components; one of variables encodes the presence or absence of each graph edge, while the second is a one-hot encoding of nodes' assigned operations. As all valid cell graphs are directed and acyelic, we limit the edge variables to the strict upper triangle of the adjacency matrix, which implicitly enforces topological ordering of the nodes in any feasible solution and ensures acyclicity. The first- and last-indexed nodes are always assigned the input and output operations respectively, while intermediates nodes can be assigned any operation from the set S = fconvixi, conv3x3, maxpool 133)).
To ensure fixed-length set of decision variables while allowing for graphs with variable number of nodes, we introduce a new null operation. Nodes assigned the null operation are not considered part of the computational graph of the cell. The algorithm then searches over the space of binary
24

---

Under review as a conference paper at ICLR 2022
representations, constrained to yield feasible cell architectures. Denoting by V and M the maximum number of allowable nodes and edges respectively, the decision variables (all binary) are:
S V, if there is edge from node to node j. 0 otherwise. may for 1 S i
Wi.k for otherwise.
S V, 1 < k [S], if node is assigned the k'th operation in S, 0
for 1 < V, if node is assigned the null operation, otherwise.
The feasible set of cell architectures can then be given in terms of linear constraints as follows:
W1.k
wv.k = 21 = zv =0
for S S [s]
8
24 +
Wi,k =
for < <<<<>
(2)
k=1
V
Σ Σ
mij S M
(3)
1=1==
Mi. S 2j
for <i<jsV
(4)
Mij $ 1 -
for 1 <i<jsv
(5)
majj 1 11-2 2j
for f
Σ $=1+1
Mi,j 2 1 2;
for SisV
S 2111
for 1 <i<<vv-
(8)
Constraints D ensure that the input and output nodes are not assigned any operation from S null, while[Zenforces the one-hot encoding of operations for intermediate nodes (including the possibility of null operation). Constraint
3 imposes limit on the number of edges in the graph, per the NAS-
Bench-101 specifications. Constraints J & 5 assert that null nodes have no incoming or outgoing
edges respectively, effectively disconnecting them from the remaining graph. Conversely, 6 7 assert that non-null nodes have at least one ingoing and one outgoing edge. Crucially, due to the implicit topological sorting of nodes by the upper-triant trianulaa adjacency matrix, these also ensure that there is always path from the input the output node using only non-null nodes. Intuitively, all non-null nodes (including the input) have at least one outgoing edge which necessarily leads to higher-indexed non-null node and all non-null nodes (including the output) have least one- incoming edge - which necessarily comes from a lower-indexed non-null node. The flow exiting the input node, must eventually enter the output node.
Finally, we focus on Constraints 8 which we refer to as symmetry--reaking constraints. These assert that node can only be assigned The null operation if its topological successor has also been assigned it. While not necessary for feasibility, this constraint serves to eliminate symmetry by ensuring that all null nodes are topologically sorted after any non-null nodes. In essence, introduces "canonical" labeling of null vs. non-null nodes, whose isomorphic representations are excluded from the feasible region.
In our experiments, we found that including symmetry-breaking constraints actually resulted in worse overall performance for the outer optimization problem (FigureRc). We hypothesize that this is due to a reduction in the exploration behaviour of MN+MILP, as the surrogate's 's predictive distri- bution was more uncertain in the larger search space and the inner-loop optimizer thus more likely to propose points in unexplored areas. One possible future line of work could be to augment D, with isomorphic representations before training, e.g. by random reordering of nodes in the representations of sampled points.
25

---

Under review as a conference paper at ICLR 2022
H
ADDITIONAL PLOTS
We present in this section additional figures related to the experiments in the main text.
H.1
UNCONSTRAINED OPTIMIZATION
8808/10,10)
RandomMLP125.3
TTBind
Nomalizee score
1.00
Nomalized Score
1.00
Normaized
0.25
0.25
0.00
0.00
BB B
NIV+
Enzombles EntombleeRRgEse RegEvo
RegEva
NN+
MV-RegEvo
Encemble+ RegEvo
RegiEvo
EnsemblesAegtt
MILL
NW HN--RegEvv
RegiEvo
Figure 8: Distribution of algorithms' normalized AUC scores (Section H.I on unconstrained prob- lems split by objective function class. Higher better. Relative performance of algorithms in terms of AUC is similar as best-observed reward (Figure I).
While the best observed reward in D (i.e. after all evaluations) is the primary metric of comparison for algorithms per Section 3.1 it is also instructive to consider a measure of how fast algorithms converge to their best observed reward. To this end, define an AUC metric that computes the area under the best observed reward curve; higher values indicate that an algorithm found better points in earlier iterations. To faciliate comparison across problems, AUC scores within each problem exactly as we did for the best observed reward (Section 4.J). That is, the best (resp. worst) on-average algorithm in terms of AUC is assigned score of (resp. zero) and intermediate values express relative distance from these extremes.
min/max normalize algorithms'
Figure 8 plots the distribution of algorithmss' AUC scores over all unconstrained problems, split by objective function class. The relative performance of algorithms in terms of this new AUC metric does not differ significantly from what we found for final reward (Section 4.2 Figure П Figures 11 and 12 plot the individual reward curves function of outer-loop iteration for each unconstrained problem.
H.2 CONSTRAINED OPTIMIZATION
In Figure 9
we plot the distribution of algorithmss' normalized final reward scores across all con- strained problems from Section 4.3 (paralleling Figure top). The individual reward curves for all
such problems can be found Higures 13 and 14] which did not differ significantly.
Normalized Score
0.6
0.2
NN UN+MMLPP MILP
ConEvo
ConEvo
NN+ Rejjamlee
Figure 9: Distribution of algorithms' normalized scores (Section 4.1
constrained problems.
Higher is better. MN+MILP and NN+ConEvo perform similarly across all problems. Reward pro- gression individual problems can be found in Figure 13 and
26

---

Under review as a conference paper at ICLR 2022
H.3 PRACTICALITY OF MILP
Table 2:
Distribution of per-step MILP immer-opiimizttion solve times in seconds for TfBind8 bench- marks when using different surrogate network architectures. The Network column denotes the num- ber of ReLUs in each fully-connected hidden layer. Runs were given a time limit (TL) of 300s. (*)
means that the time limit was hit.
Network
min
med
95%
99%
max
%TL
Linear
0.004
0.4
1.4
2.9
16.9
0%%
FC(16)
0.02
2.2
8.0
15.5
60.8
0%
FC(32)
0.04
11.7
49.2
85.5
300
0.1%
FC(16,16)
0.40
12.2
55.6
109.1
300
2.1%
Here we present experiments to explore the impact of surrogate network size on MILP solve times. We varied MILP's surrogate network architecture using fully-connected networks with different numbers of hidden layers and neurons, FC(16), FC(32), FC(16,16), as well as simple Linear model (no hidden layer). Note that the first architecture, FC(16), was what was presented in the main used the same training and optimization syper-parametery for MN+MILP as described in Section C.T We ran this experiment on all 12 unconstrained 7fBind problems from Section 4.2 with 20 trials using different random initial data sets.
results. Otherwise,
Table 2 shows aggregate distribution statistics of inner-loop optimization runtime for the different architectures, across all steps of all trials of all problems. We note that solve times increase as the network size increases, but even for the largest network (two layers with 16 neurons each), the solver rarely times and almost always terminates within a practical time limit. Furthermore, for larger networks we can improve scaling using advanced formulation techniques (e.g. Appendix [A]. We also note that, in these experiments, there was no single architecture that consistently produced better optimization across different instances (though Linear was almost always outperformed by the rest).
We also include Figure 10 which plots the distribution of MILP inner-loop solver runtimes for all constrained and unconstrained problems as a function of outer-loop iteration (paralleling Figure 2b which included only TfBind). As noted in Section 4.4 the relationship between problem size and runtime can be unpredictable. For example, the lower-dimensional TfBind problems showed the highest mean and variance in MILP inner-loop solve times among all constrained and unconstrained objective classes. Encouragingly, across all problem classes we observe only a linear (roughly) increase in solve time per iteration, presumably due to the increasing number of no-good constraints.
(803(15,10)
3undumMuri25.5)
TIB 41
Condra nes P11
MILP SAie Time
10
MILF Silee
100
200
800
250
830
200
000
Figure 10: Distribution of MILP acquisition problem solve times as function of iteration split by objective class for unconstrained problems (Section 4.2) and for all constrained problems (Sec- tion 4.3). Line and bands show the median and 5th/95th percentile range over all trials of all prob- lems in class.
27

---

Under review as a conference paper at ICLR 2022
BBOB(10.
SECTOR
8608(10.10)
1e-5 8608(10.10) DISCUS
0.0
0.00
Best Reward
Best Reward
-0.2
-0.4
-0.10
-1.0
-0.15
200
400
600
800
1000
200
400
600
800
1000
200
400
6.0
800
1000
Stop
Stop
Stop
B808/10.100
8808(10,10) GALLAGHER 21ME
8B08/10,10) BBOB/10.10LSCHAFFERS.. SCHAFFERS_
Best Reward
Reward
Beit Reward
-0.25
-0,025
-0,050
-0.50
000075
-0.75
-0.100
250
500
750
1000
200
400
600
800
1000
200
400
6:0
800
1000
Stop
Step
Step
B808(10,10) 10LSHWWEFEL SCHWEFLL
B8OB(1 0, SI_SPE
(188810010) STEP ELLIPSOID
0.0
Reward
-0.2
Reward
-0.2
Best Reward
-0.1
-0.4
0..4
-0.6
-0.6
0.8
200
600
1000
1000
Step
Step
Step
0808(10,10)
SOOBELL...LLWEERRSTTAASSS WEIERSTRASS
RandonCCNW/2555..eeee seed-:0
RandomCNN/25.0 (ttp:/NNN///553eddeeada seede11
3.5
Best Reward
Reward
Best Reward
2.0
-4
0.0
2.0
Stap
RandomCNN:25.5 51_ 51_sedmm42
RandomCNN/25.35 517dd=77
RandomFCN125..0 (050ii2255))__dddd
3.0
Roward
Rewarrd
Best Reward
2.0
2.0
1.5
1.5
1.0
Stop
Stop
Step
MILP
RegEvo
Ensemble RegEvo
Figure 11: Best observed reward as function of iteration for the first half of all unconstrained problems (Section[4.2), averaged over 20 iterations (bands indicate = 1sd), Dashed grey lines in the first 50 steps indicate the initial randomly sampled dataset, common to all methods except RBFOpt, which performs its own initialization.
28

---

Under review as a conference paper at ICLR 2022
Randomm seedm13
RandomFFNN225.0
RandomFFN1255
3.0
Best Reward
2.5
Best Reward
2.0
2.0
1.5
1.3
1.0
1.0
200
400
600
800
1000
2.00
400
600
800
1000
200
400
600
800
1000
Stop
Stop
Stop
TTBind(8,4) CRX ROOW
TIBind(B,4)_CRL REF
T1Bindr8.4), (0F_0) _RFF R1
1.0
1.0
1.0
Best Reward
Reward
Best Reward
0.8
0.8
0.7
0.4
0.6
200
400
600
800
1000
0.2
200
400
600
800
1000
200
400
600
800
1000
Stop
Stop
Stop
TrBin 4J11 REF R1
TrBind(8 41 HDXD13 Q3255_ R1
TrBind[8.4] T1Bindj8,4_HOXD113 HDXD13 REF
Reward
Reward
0.8
Best Reward
0.8
0.6
0.6
400
2.00
1000
800
Step
Stop
Step
T/Bindt8. 4) NR,H4 C144R, R1
T/Bind(8.4) NR,HH REF R1
TIBinot8 TIBino(8.4)PAAKC 4))PAXX4_AEEE R1
1.0
Reward
0.9
Reward
Best Bewaard
0.7
Stop
T[Bindl8. PAX4 PAX4_REF A2
TBBind88.4) 40_GFP0 REF_AA
TIBind 4)_SIX6__EEE R1
1.0
1.0
Best Rowwers
0.8
Reward
Best Reward
0.8
0.8
0.6
0.4
0.2
Stop
Stop
Stop
MILP
RegEvo
Ensemble RegEvo
Figure 12: Best observed reward function of iteration for the second half of all unconstrained problems (Section 4.2), averaged over 20 iterations (bands indicate + 1sdd, Dashed grey lines in the first 50 steps indicate the initial randomly sampled dataset, common to all methods except RBFOpt, which performs its own initialization.
29

---

Under review as a conference paper at ICLR 2022
RandomCNN:1100.2
RandomCNNi100.22
Random
Best Reward
Best. Reward
Best
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
Stop
Stop
Stop
RandomCNN(100.21,
Random/CNN1100.7 k=10
RandomCCWWr10..2) seed=13 k525
Best Reward
Reward
Best Reward
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
Step
Stop
Step
RandomCNM(1002) 40 k-5
send-44
RandomCNN/1000 500d-42 k-10
RandomCXVN(10072) 42 k-25
Best Best Rewardd
Best Rerward
400
1000
400
1000
Step
Stop
Step
RandomCNN(110222 seed- seed=
k=5
RandomCCNW1100.2 seed=77 7-000000 k=10
Random-CNN(100.2) 7eed=77 77.7777-111 k-
Best Reward
Rewwrd
Best Reward
RandomCNN/1000 (ttp://i000002sseddmm00 100_k=5 k-5
(NW1000CNW::00.22 (0lts_000)___eddeeee 100 100_k=10
RandomCNN1100022 00
Best Roward
Best Roward
Besi Reward
1000
Stop
Stop
Step
NN-MILP
NN-ReeSSample
ConEvo
Ensembler RegEvo
Figure 13: Best observed reward as function of iteration for the first half of all constrained prob- lems (.cc).on443). averaged over 20 iterations (bands indicate #1sd), Initial randomly sampled set of 50 points is omitted.
30

---

Under review as a conference paper at ICLR 2022
FCN/1002N
RandomFCN1100 20
RandomnFNN:
3.0
Best Reward
Rewarrg
Begt Reward
1.5
1.0
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
Stop
Stop
Stop
RandomFCN/100 seedd]]
RardoomFEN(100.2__eeed
k-10
RandomFCN:100.22) seed-]] k-25
Bost Reward
Reward
Best Rewardd
200
400
600
800
1000
200
400
600
800
1000
200
400
600
800
1000
Step
Stop
Stop
RandomFCN/100 2] 50 k-5
BardommFCM110022)eeedddd k-10
RandemFCN/100,2) 50cd-42 k-25
Rewarrd
Reward
Best Reward
400
1000
400
1000
600
Step
Stop
Stop
RandomFCNN1000 21 km5
BardommFCN/100222sseeeee 77_kk-7711110 k=10
(ardemmFGW/10022,sseeee k=25
Rowarrd
Best Reward
RandomFCN(100021, 20e002222eeeee000 km5
Rendom k=10
FandomFCN110002__eeed 51050dd00000- 20
Best Roward
Best Reward
-1
-2
-2
Stap
Stop
Stop
NN-MILP
KN-RejSample
ConEvo
Ensembles RegEvo
Figure 14: Best observed reward function of iteration for the second half of all constrained problems (.ection[4.3). averaged over 20 iterations (bands indicate 11sd). Initial randomly sampled set of 50 points is omitted.
31


Reviews:
[{"review_id": "go2CFvHRCzC", "reviewer": "Reviewer_j85k", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the paper is well-written and suggests an interesting and relevant approach to address black-box optimization problems equipped with discrete domains. I found the basic framework to be well thought-out and, in my view, of potential value for a large array of settings within this field.\n\nMy major concern, however, is that many design choices are somewhat unclear, and the paper often feels that lacks depth in more specific areas. In particular:\n\n(1) Optimizing inputs with neural networks is challenging with MILPs, often requiring more sophisticated implementations as that of Anderson et al., 2020 (as nicely emphasized by the authors). I wonder whether the fact that the authors were limited to a very simple network, with only a single layer, did not hinder some key insights on the numerical evaluation? For instance, wouldn't the function approximation be quite poor for larger instances? Also, is this approach really scalable?\n\nPerhaps my suggestion here is for authors to consider encodings that are much more efficient/scalable for MILPs (and other model-based approaches). For example, many black-box functions have a natural discrete structure, and authors could consider a decision tree, which is much \"simpler\" to optimize.\n\nSimilar reasoning applies to the \"no-good\" constraints, which in this case are quite simple and only eliminate one solution at a time. This could really impact MILP performance (e.g., similar to combinatorial Benders cuts). Would there be cases where it is possible to eliminate several points from the acquisition domain? For example, suppose \"x\" are binary and represent subsets. If you eventually learn all subsets of size <= N, you could eventually replace those constraints by \\sum_i x_i >= N + 1.\n\n(2) The numerical experiments do not seem to reflect well the benefits of the approach. To the best of my understanding, the final conclusion is that problems are \"easier\" to model with NN+MILP, but performance improvements are marginal (if any). I believe this is indeed the case, but the paper lacks more concrete evidence of this statement. \n\nIn particular, when comparing NN+MILP and NN+ConEvo, why not consider a problem class that the MILP has clear benefits? For example, a set-packing/set-covering acquisition domain, perhaps with other side constraints, or a scheduling feasible set (e.g., \"f\" could encode a weighted completion time with unknown weights, and the feasible set are the valid sequences). The authors could compare with any global optimizer (as opposed to RejSample and ConEvo)  because it is well-known that MILP is one of the state-of-the-art techniques for these problems.", "summary_of_the_paper": "The paper develops a technique to optimize an unknown, black-box function \"f\" by leveraging a combination of neural networks with mixed-integer linear programming (MILP) methodology. More specifically, authors encode an approximation of \"f\" using a neural network with piecewise-linear activation functions, which is optimized using its associated reformulation as an MILP with no-good cuts. Numerical results evaluate the approach with respect to other baselines and neural network optimization mechanisms. ", "summary_of_the_review": "The ideas are novel and significant, especially given the modeling expressivity provided by MILPs. However, the paper lacks some justification concerning the scalability of the approach and the fact that neural networks had quite a limited size. Moreover, in my view, the numerical experiments do not explore well the benefits of their approach.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "e2YgquJ5UT", "reviewer": "Reviewer_qsFh", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Reasons to Accept\n\n+ This paper is well-written and well-organized.\n+ It solves a very interesting problem that contains constraints on discrete search space.\n\n### Reasons to Reject\n\n- I do not think that a piecewise neural network models an uncertainty of unknown objective appropriately.\n- Following the point described above, even though we train a neural network every iteration, it tends not to reflect a factor for exploration, which implies that regression results are almost same where the same observations are given.\n- GP is a popular choice for a surrogate function, since it has sufficient expressiveness, which is defined on a RKHS space. I am curious that the surrogate function used in this paper is sufficiently expressive of unknown function.\n\n### Questions to Authors\n\nPlease answer the comments described in Reasons to Accept and Reasons to Reject.\n\n1. Can I ask what the difference between no-good constraints and generic equality constraints is? If they are similar, it can be considered in a continuous search space?\n1. Is there any specific reason not to compare the proposed method to COMBO? In my experience, COMBO can be thought of as the state-of-the-art model now. To cope with a no-good constraint, you can just apply a rejection strategy in a step of acquisition function optimization (i.e., local search with rejection sampling).", "summary_of_the_paper": "This paper solves a constrained discrete black-box optimization problem that employs a surrogate model in modeling an unknown objective function. Unlike the formulation of standard Bayesian optimization, it constructs a surrogate model using a piecewise linear neural network. Under the assumption that a randomly-initialized neural network is able to produce an uncertainty for exploration (against exploitation), the proposed method optimizes a surrogate model directly with mixed-integer linear programming. It follows a spirit of Thompson sampling. Finally the authors conduct their method on several experimental circumstances and show the validity of their method.", "summary_of_the_review": "It solves a very interesting topic which is defined on a discrete search space with constraints. However, the choice of surrogate model is not convincing and the baseline is missing. Thus, I would like to recommend rejection.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "_1U-aQ3Zdfd", "reviewer": "Reviewer_5edm", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Strong points\n- The declarative nature of the proposed method eliminates the need for developing an algorithm for solving the inner loop optimization (as long as the structural constraints can be represented in MIP formalsim). It also facilitates to adpating an existing method to similar variations by simply modifying the constraints. \n\n### Aspects to be improved\n- The choice of a piecewise linear neural network for both the surrogate model and acquisition function is not properly motivated. In particular, it is not clear how the proposed approach maintains an exploration-exploitation balance. \n\n- The proposed approach does not demonstrate significant improvements in empirical evaluations. In the unconstrained setting, it is outperformed by the baseline, and in the constrained setting, solving the inner loop problem to optimality does not seem to provide an advantage. The argument of \"ease of implementation\" would have been convincing if implementing the alternative (evolutionary) approach was prohibitively difficult. But this does not seem to be the case. \n\n- The proposed approach relies on solving a MIP problem repeatedly. This does not allow the method to be applied to problems with more than a certain number of variables. The problems studied in the  experiments have few number of variables. In particular, *TfBind(8, 4)* can be probably solved by simple enumeration. It is not clear to what extent the proposed method is applicable to larger problems. \n\n### Question\n\n- When comparing different methods, each algorithm is evaluated in terms of the best reward observed after 1000 queries. At each iteration, the MIP solver is given a 500 seconds timeout, and the evolutionary aglorithms is given a budget of 10k queries. Does this give the inner loop algorithms equal opportunities for finding good solutions? Isn't it more fair to compare the methods subject to an equal overall time budget?\n", "summary_of_the_paper": "The paper proposes a method using piecewise-linear neural networks as the surrogate model (and the acquisition function) in a black-box optimization framework. At each step, the acquisition function is optimized by casting the learned neural network (and a set of constraints excluding the already-visited data points) into a mixed integer linear program, which can be then solved using off-the-shelf solver. The proposed method is empirically evaluated on a number of unconstrained and constrained tasks. \n\n\n\n\n\n\n\n\n", "summary_of_the_review": "The paper presents an interesting direction for using a declarative framework (i.e. mixed integer linear programming) to encode combinatorial structures in black-box optimization. The proposed method lacks a theoretical motivation, and the empirical evaluation does not demonstrate significant advantages over existing methods. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "Tng2d6P_20p", "reviewer": "Reviewer_vgXk", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "# Strengths\nThis paper tackles an interesting problem of optimizing expensive blackbox functions even in the case of complex combinatorially constrained domains, by using a MILP formulation of a piecewise linear neural network. Even though such formulations have been used before (the authors mention neural network certification as an example), it offers a very flexible way of easily integrating combinatorial constraints, without handcrafting a local search method to maintain feasibility of solutions to the acquisition problem. The experiments show, that\n1) solving the acquisition problem to global optimality can provide benefits over local search methods even in case of unconstrained domains\n2) performance of handcrafted methods to ensure proposal feasibility in the constrained case can either be matched (NN+ConEvo in RandomMLP with subset constraints) or surpassed (RE in NAS-101), with a simpler way to implement the constrained domain in a general framework.\nWith the flexibility of the framework this can be used for various other future problems that involve MBO in the presence of combinatorially constrained domains.\n\n# Weaknesses\n### Limitations of surrogate model\nThe surrogate model is limited to a) neural networks with piecewise linear activation functions, this is inherent to the MILP formulation and restricts the framework from using potentially better-suited surrogate hypothesis classes such as Random Forests.\nThe surrogate model is also limited to b) a relatively small number of neurons in the network due to runtime limitations of the MILP. The dimensionality and number of constraints in the MILP formulation scales linearly with the number of nonlinearities in the network and as mentioned in the paper, the runtime scaling of the MILP solver is then often unpredictable. \nThis means that the runtime comparison between methods is a crucial factor in the comparison, which leads to some questions.\n\n1) In section 4.4 the authors reported the inner-loop optimization runtimes for the unconstrained case. To me it is very surprising that the local evolutionary search RegEvo has a larger average runtime than the global MILP optimization. It would be great if the authors could provide intuition on why this is the case. Additionally, is this also the case when comparing MILP to ConEvo in the constrained case, where those two methods are the main competitors? As the constrained case is the one this framework is designed for, this would be the more important comparison than in the case of an unconstrained domain.\n2) In the experiments a single hidden layer with 16 neurons is used, and ablations with two slightly bigger architectures (32 and 16+16 neurons) are provided. While the slightly bigger architectures are still within the computational limits, the network probably cannot be scaled to much larger (deep) networks. In the ablations it is stated that for the TFbind8 experiment larger architectures did not improve performance, which suggests that the small surrogate model is already expressive enough in this case. It would be interesting to know whether this is also the case for a larger experiment such as the NAS case study. \n\nIn any case, the runtime restriction could be a limiting factor in future applications, where more complex surrogates could be required, giving faster but less optimal local search methods an advantage over the global optimization.\n\n### Integration of continuous variables\nIn the current stage the framework doesn’t support mixed-integer domains because the no-good constraints do not naturally extend to this case, as stated in the conclusion. Therefore, in its current stage, the method should be called NN+ILP instead of NN+MILP, as the name otherwise suggests that mixed-integer domains can be used.\n\n\n## Typos\nThese do not affect my rating.\n1) Before section 2.2: ‘are’ should be removed", "summary_of_the_paper": "The authors present NN+MILP, a framework for the optimization of an expensive to evaluate blackbox fuction with a discrete combinatorially constrained domain. The acquisition problem of finding the surrogate minimum is solved to global optimality by solving a MILP formulation of the acquisition problem. The MILP formulation limits the considered neural network surrogate class to networks with piecewise linear activation functions. However, it provides a simple declarative language for integrating the problem-specific constraints. \n\nThe experiments cover both the cases of unconstrained and constrained optimization of the acquisition function. \nThe unconstrained case compares NN+MILP to general purpose algorithms for unconstrained discrete blackbox optimization. It shows that the global optimization often achieves better results than a local search evolution based method for solving the acquisition problem, while also solving the inner loop problem faster. Additionally, it demonstrates that even with the restrictive function class of a single layer neural network, comparable performance to better-suited surrogate hypothesis classes can be achieved, thanks to the global optimization. \nIn the constrained case artificial subset-equality constraints are used to create a combinatorial domain for the blackbox function. NN+MILP performs similar to NN+ConEvo, a manually adjusted local search method that ensures feasibility of the proposals at every step. Other methods employing random search for inner loop optimization or local search with a different surrogate hypothesis class perform much worse.\nFinally, a case study for the NAS-Bench-101 benchmark is provided. A novel MILP formulation for cells of a valid architecture design is described, which consists of an MILP formulation for directed acyclic graphs with added null operations to allow DAGs with a reduced number of cells. Despite the generality of NN+MILP it outperforms a strong evolution based baseline.", "summary_of_the_review": "Strengths: The presented framework is flexible and easier to implement than other methods in the presence of combinatorially constrained discrete domains. The experiments show that the performance is also competitive with other approaches.\n\nWeaknesses: The MILP formulation limits the hypothesis class of models to neural networks with linear activation functions, runtime limitations restrict the number of neurons in the network. Mixed-integer domains are currently not supported by the framework, even though the name suggests otherwise.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper considers the problem of black-box optimization and proposes a discrete MBO framework using piecewise-linear neural networks as surrogate models and mixed-integer linear programming. The reviewers generally agree that the paper suggests an interesting approach but they also raised several concerns in their initial reviews. The response from the authors addressed a number of these concerns, for instance regarding scalability and expressivity of the model. However, some of these concerns remained after the discussion period, including doubts about the usefulness for typical applications in discrete black-box optimization and some concerns about the balance between exploration and exploration. 

Overall the paper falls below the acceptance bar for now but the direction taken by the authors has some potential. I encourage the authors to address the problems discussed in the reviews before resubmitting.
Decision: Reject