Year: 2021

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.

heading: 1 INTRODUCTION
text: The problem of optimizing an expensive black-box function f : Ω 7→ R over a 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-constrained settings, it is desirable to develop algorithms that exploit known combinatorial structure in Ω to search the space more efficiently.
Model-based Black-box Optimization (MBO), a popular paradigm that includes Bayesian Optimization as a special case, iteratively refines a function approximator f̂(x) ≈ f(x) and selects new points to query by optimizing an acquisition function a(x) derived from a point estimate or posterior distribution over f̂ (Section 2.1). This inner-loop optimization problem is assumed to be easier than the original, since, for example, a(x) is less expensive to query than f(x) or provides “white-box” properties such as gradients.
There is a vast literature addressing the challenges of applying MBO in practice. We focus on two of these: first, optimizing a(x) may itself be a computationally-difficult optimization problem; second, in many applications, practitioners are confronted by additional constraints on x. For example, in neural architecture search, x might represent a computation graph that must be both connected and acyclic. Due to the difficulty in optimizing the acquisition function over a combinatorial domain, most approaches resort to heuristic inner-loop solvers, which often need to be specialized to the problem at hand to ensure feasibility, e.g. evolutionary solvers with custom mutation operators.
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 a 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 f̂ (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
performance of MILP for optimizing such NNs with respect to their inputs (Anderson et al., 2020). For us, MILP is practical to use in the inner loop because the dimensionality of typical black-box optimization 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 use 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 a simple declarative language for problemspecific constraints.
• We show that NN+MILP matches or surpasses the performance of strong MBO baselines based on problem-specific evolutionary algorithms on a wide range of synthetic and real-world discrete black-box problems.
• We use the NAS-Bench-101 neural architecture search benchmark as a case study, presenting a novel MILP formulation of its graph-structured domain.

heading: 2 BACKGROUND AND RELATED WORK

heading: 2.1 MODEL-BASED BLACK-BOX OPTIMIZATION
text: Model-based Black-box Optimization (MBO) is a broad family of methods that includes Bayesian optimization as a special case (Mockus et al., 1978; Jones et al., 1998; Hutter et al., 2011; Snoek et al., 2012; Shahriari et al., 2015). As depicted in Algorithm 1, the method proposes xt at iteration t using three steps. First, the user performs inference over a surrogate model f̂ to approximate f using the data previously collected from the black-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 f̂(x) is posed that quantifies the quality of new points to query. Finally, xt is 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(x) may be orders of magnitude less expensive to evaluate or have a tractable functional form. Practitioners can encode prior knowledge about the structure of f via a choice of inductive bias for f̂ , e.g. a 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 to 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 closedform posterior inference, naturally adjust their expressivity as the dataset grows, and users can inject domain knowledge via a choice of kernel (Rasmussen & Williams, 2006; Oh et al., 2019). On the other hand, neural networks provide a practical alternative (Snoek et al., 2015; Hernández-Lobato et al., 2017), since they often scale more gracefully, either computationally or statistically, to large datasets or high-dimensional domains.

heading: 2.2 SOLVING THE DISCRETE MBO ACQUISITION PROBLEM
text: 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 limited guarantees (Hutter et al., 2011; Müller, 2016; Oh et al., 2019; Kandasamy et al., 2020). A key advantage of 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 (approximately) solved using specialized combinatorial solvers (Baptista & Poloczek, 2018; Deshwal et al., 2020), mixed-integer nonlinear programming (MINLP) solvers (Costa & Nannicini, 2018; Kim & Boukouvala, 2020), or continuous optimization solvers (Bliek et al., 2021).
Algorithm 1 MBO Input: hypothesis class F , budget N , initial dataset Dn = {xi, f(xi)}ni=1, optimization domain Ω for t = n+ 1 to t = N do P (f̂t)← fit(F ,Dt−1) a(x)← get acquis func(P (f̂t)) xt ← inner solver(a(x),Ω) Dt ← Dt−1 ∪ {xt, f(xt)}
end for return arg max(xt,yt)∈DN yt
Algorithm 2 NN+MILP Input: hypothesis classF , budgetN , initial dataset Dn = {xi, f(xi)}ni=1, MILP domain formulation MΩ for t = n+ 1 to t = N do f̂t ← fit(F ,Dt−1) (3.2) Mt ← build milp(f̂t, MΩ, Dt−1) (3.3) xt ← optimize(Mt) (generic MILP solver) Dt ← Dt−1 ∪ {xt, f(xt)}
end for return arg max(xt,yt)∈DN yt
Therefore, practitioners must decide between either introducing difficult-to-analyze approximations due to inexact heuristic solvers or using tractable surrogate models that may be mis-specified for the application domain. This serves as a key motivation for our work: we seek to enable practitioners to employ broad families of surrogate models and exactly solve the acquisition problem with reasonable computational overhead in practice.

heading: 2.3 CONSTRAINED MBO
text: In many applications, x is subject to non-trivial structural constraints. Prior work has largely focused on the case where determining whether x 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; Hernández-Lobato 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(x).
Problems with inexpensive white-box h(x) can be tackled using these approaches for black-box constraints, but doing so may lead to slower optimization and may query f(x) at invalid x, which can be unsafe when performing physical experiments (Berkenkamp et 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 customized for each family of constraints, a task usually left to the user. Prior work employing MINLP solvers addresses white-box constraints either by adding a penalty for constraint violation (Costa & Nannicini, 2018) or in small-scale settings (Kim & Boukouvala, 2020).

heading: 2.4 MIXED INTEGER LINEAR PROGRAMMING
text: Mixed Integer Linear Programming (MILP) seeks to maximize a linear function over a set of decision variables, some of which may be integral, subject to linear inequality constraints. Decades of development have allowed MILP to have a significant impact in a wide range of applications due to its better-than-expected computational performance (Jünger et al., 2010). Indeed, while MILP problems are computationally hard (NP-complete), they are routinely solved (to global or near-global optimality) in production environments thanks to state-of-the-art solvers that nearly double their machine-independent performance every year (Achterberg & Wunderling, 2013; Bixby, 2012).
A notable aspect of MILP is that it 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 to 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-linear functions (Huchette & Vielma, 2019; Vielma et al., 2010). Specifically, our work leverages MILP formulations for trained neural networks with piecewise-linear activation functions such as ReLUs (Anderson 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, 2017; Tjeng et al., 2019), reinforcement learning (Ryu et al., 2020; Delarue et al., 2020), and analysis and exact com-
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.

heading: 3 MILP FOR MBO
text: We propose the NN+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 a wide variety of models and leverage MILP’s versatile declarative language to incorporate constraints. This section describes various design choices to make the approach practical.

heading: 3.1 PROBLEM SETTING
text: Our goal is to find: x∗ = arg max
x∈Ω f(x), (1)
where f : Ω 7→ R is an expensive, noiseless black-box reward function and Ω ⊆ Ω1 × . . .×Ωn is a domain on n decision variables. We assume Ω can be described by an inexpensive function hΩ(x) indicating whether x is in Ω. Algorithms are allowed a fixed budget of N sequential queries to f . Xt := {xi}ti=1 refers to the set of sampled points by iteration t, and Dt := {xi, yi = f(xi)}ti=1 includes corresponding rewards. An algorithm’s performance is measured as the best reward inDN . Since f is noiseless, it is advantageous for algorithms to avoid repeated evaluations of the same x.
We choose to focus on finite discrete sets Ω as we believe this is the area where MILP can provide the greatest benefit. As noted in Section 2.4, there are many well-studied formulation techniques for Ω with combinatorial structure, such as directed graphs. More generally, such sets have a polynomially-sized MILP formulation whenever hΩ(x) can be evaluated in polynomial time (e.g., Yannakakis (1991)). 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.

heading: 3.2 SURROGATE MODEL AND ACQUISITION FUNCTION
text: For surrogate model f̂ , we allow any feedforward neural network with piecewise-linear activation functions, as they can be represented by MILP (Section 3.3). Though we focus on fully-connected networks using ReLUs, we note that a wide range of architectures (e.g., those including convolutional and max-pooling layers) are piecewise-linear and may place suitable inductive bias on f̂ .
Our experiments employ a simple acquisition function loosely motivated by Thompson sampling (Thompson, 1933; Hernández-Lobato et al., 2017; Kandasamy et al., 2018). At each optimization step, we train a new regressor f̂(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 et al., 2017). We discuss extensions to alternative acquisition functions in Section 6. We use a flattened one-hot encoding of x for the input layer, and train each network f̂t ∈ F on Dt−1 using `2 loss. Finally, we re-scale the observed rewards inDt−1 before training models, to aid both in training and optimization. Poorly-scaled data may result in slower performance or small inaccuracies in MILP solvers (Miltenberger et al., 2018).

heading: 3.3 MILP FORMULATION OF THE ACQUISITION PROBLEM
text: The inner-loop solver then seeks to find
xt = arg max x∈Ω\Xt−1 f̂t(x), (2)
where Ω is the feasible set for (1) and Xt−1 is the set of points where the noiseless f(x) has been queried already. The MILP formulation of (2) is denoted byMt and has the following three components:
Domain We use a one-hot encoding of decision variables x (unless they are already binary), defining the binary decision vector z with zij ≡ I{xi = j} for i ∈ [n], j ∈ Ωi, and subject to linear
constraints ∑
j∈Ωi zij = 1 ∀i. Integer domains with small range may be one-hot encoded; see Appendix D for a comparison between integer and one-hot encodings. Additional constraints due to Ω are added as necessary, with form dependent on the application at hand. We assume that these are MILP-representable, which as noted in Section 2.4 could include a wide range of combinatorial, logical, and polyhedral constraints. We useMΩ to denote the domain formulation itself. No-good Constraints A no-good constraint is one that eliminates one or more undesirable solutions from the domain. Here, we leverage the binary nature of z to exactly eliminate the set Xt−1 from Mt. For illustrative purposes, consider a single point x̄ ∈ Ω we wish to exclude from the acquisition problem’s domain, and let z̄ denote its one-hot encoding (or x̄ itself if the problem is binary). Then the constraint: ∑
i,j : z̄ij=0
zij + ∑
i,j : z̄ij=1
(1− zij) ≥ 1 (3)
enforces that any feasible z has a Hamming distance of at least 1 from z̄. As z are binary, this effectively eliminates just the single point z̄ from the feasible region. We therefore formulate Ω \ Xt−1 by including one such constraint for each x̄ ∈ Xt−1. Note that the right-hand side can be tightened to 2 for one-hot encodings, and these compact no-good constraints do not extend naturally to continuous x (Section 6).
Neural Network We formulate the neural network by introducing auxiliary decision variables encoding the activation of each neuron for a given z. We present here the formulation for a single ReLU, commonly used throughout the literature (Section 2.4), while noting that the full formulation is obtained by combining all ReLU formulations and matching their input and output variables according to the structure of the network. The overall MILP objective is the activation corresponding to the regressor’s output neuron.
A ReLU neuron with vector input x and scalar output y has the piecewise-linear form y = max(0, w>x + b), where w and b are its weights and bias respectively. At optimization time, w and b are fixed, while x 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 a binary decision variable α that indicates whether the ReLU is active or not. We then write the following set of constraints to enforce that y = 0 when α = 0 and y = w>x+ b when α = 1:
0 ≤ y ≤Mα (4) w>x+ b ≤ y ≤ w>x+ b+M(1− α) (5)
where M is a sufficiently large fixed value, such as an upper bound on the range of y. As w and b are fixed, values for M can be computed in advance of the optimization, e.g. by propagating bounds from Ω. Our experiments use a more advanced method to compute M , detailed in Appendix A.

heading: 3.4 OPTIMALITY GUARANTEES FOR MILP
text: The full acquisition problem formulation, denoted byMt, is passed to a generic MILP solver with fixed time limit. If the solver does not time out, it is guaranteed to have produced a global optimum of (2), 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 an optimal solution before finding the upper bound that guarantees its optimality, so timing out do not imply sub-optimality. Finally, innerloop 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 a useful empirical tool for understanding the impact of exact inner-loop optimization (Section 4).

heading: 4 EXPERIMENTS
text: This section presents experimental results on a 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 NN+MILP where the inner-loop solver is replaced by an inexact evolutionary alternative. Depending on the problem, we also include independent baselines tailored to the application domain.
In all experiments, we fix the surrogate model hypothesis class F to networks with a single, fullyconnected hidden layer of 16 neurons. Models are 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 et al., 2020) using default settings and a time limit of 500 seconds. We use standard CPU machines with ∼1G RAM and ≤ 10 cores.

heading: 4.1 BENCHMARKING TASKS
text: Unless otherwise stated, tasks’ domains consist of discrete decision vectors of length n, with a common alphabet A for all elements. We consider three families of black-box objectives:
• RandomMLP The output of a multi-layer perceptron operating on a one-hot encoding of the input. Notably, architectures have significantly more layers/parameters than the 16-neuron networks used as surrogates by NN+MILP.
• TfBind Binding strength of a length-8 DNA sequence to a 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 are treated as unordered and categorical.
We use parentheses after the family name to denote dimensionality of a problem, e.g. RandomMLP(10,5) refers to a RandomMLP objective over a discrete domain with n = 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 a 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 a given problem is assigned a score of one (resp. zero), and intermediate values express relative distance from these extremes. No hyper-parameter tuning was performed across problems for any algorithm.

heading: 4.2 UNCONSTRAINED OPTIMIZATION
text: Before considering problems with combinatorial white-box constraints, we first tackle simple problems with no additional constraints on the discrete domain, i.e., Ω = An. This allows us to compare against general-purpose algorithms for unconstrained discrete black-box optimization. We vary the problem sizes over 30 functions, consisting of eight RandomMLP(25,5), ten BBOB(10,10) and twelve TfBind(8,4) targets (Appendix B).
NN-MILP provides an analytical tool for understanding the relative impacts of the choice of surrogate 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 configuration details are provided in Appendix C.
• 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 et al. (2020), using an ensemble of linear and random forest regressors as the surrogate, where hyperparameters are dynamically selected at each iteration. The acquisition function is the ensemble mean and inner-loop optimization uses RegEvo.
• RBFOpt A competitive mixed-integer black-box optimization solver that uses the ‘Radial Basis Function method’ as a surrogate model (Costa & Nannicini, 2018).
Figure 1 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 NN+MILP 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. NN+RegEvo. The only difference between the two is the former’s stronger optimality guarantees when solving the acquisition problem. We observe that NN+MILP obtains a greater or equal score than its evolution-based counterpart in 22 of the 30 problems considered, and variance in its normalized scores is lower within a given objective family.
The comparison of NN+MILP and Ensemble+RegEvo 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 a single neural network surrogate and exact optimization is able to achieve comparable performance.

heading: 4.3 CONSTRAINED OPTIMIZATION
text: Next, problems are augmented with combinatorial constraints on the domain. We simulate finebalance constraints in observational studies (Zubizarreta et al., 2018; Bennett et al., 2020), where the same number of items must be selected from given sub-populations (e.g., groups sharing a common attribute). These simple, yet highly combinatorial, constraints allow for fair comparison with evolutionary 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 a binary alphabet A = {0, 1} to indicate whether each of n = 100 items is selected. These have been partitioned into given subsets S1, ..., S2k for some integer k and constraints enforce that the number of selected items is equal in pairs of subsets: ∑ i∈S2j−1 zi1 = ∑ i∈S2j zi1 for j ∈ [k]. We use equally-sized subsets, each of cardinality n2k , creating 30 problems over a grid of ten objectives and three values of k ∈ {5, 10, 25}. Further details can be found in Appendix B. RandomMLP(100,2) functions then simulate non-linear reward structures for a given selection. A different objective class based on the Ising model is also considered in Appendix E.
The following optimization approaches provide ablations to contrast declarative vs. procedural approaches to handling the constraints. Further configuration details are given in Appendix C.
• ConEvo RegEvo with a 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 as 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-random from the domain. By using rejection sampling, the solver can leverage the same declarative definition of the constraints as in NN+MILP.
• Ensemble+RegEvo Same as Section 4.2, using the constraint definition to assign negative rewards to infeasible proposed points, which in turn affect the surrogate model.
Figure 2a plots algorithms’ best observed reward as a function of iteration for a typical problem. Here, we observe that NN+MILP and NN+ConEvo perform similarly, both benefiting from the ability 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 NN+MILP and NN+ConEvo to be statistically significant in terms of optimization performance, they differ considerably in terms of ease of implementation. In particular, NN+MILP 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-representable constraints. Conversely, NN+ConEvo relied on the implementation of a custom mutator, tailored to the constraint structure at hand, and would likely require significant reworking if other constraints were added.

heading: 4.4 PRACTICALITY OF MILP
text: Figure 2b plots the distribution of MILP solve times for inner-loop optimization as function of iteration over all trials of all TfBind problems. Despite the computational complexity of the acquisition problem, MILP finds globally optimal solutions in seconds. Averaging across all unconstrained experiments (Section 4.2), the inner-loop optimization took 7.92 ± 4.23s for NN+MILP compared to 9.00±1.94s for NN+RegEvo (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 average solve times for lower-dimensional TfBind problems than for RandomMLP and BBOB. We also explore the impact of surrogate model network size on solve time in Appendix H. Solve times tend to increase over time, possibly due to the increasing number of no-good constraints and the nature of surrogate models that have been fit on more data.

heading: 5 NAS-BENCH-101 CASE STUDY
text: Finally, we use the NAS-Bench-101 (Ying et al., 2019) neural architecture search (NAS) benchmark to illustrate the power of MILP’s declarative constraint language in formulating complex combinatorial 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 a directed path, while the remaining nodes are each assigned to be 1x1 convolution, 3x3 convolution, or 3x3 max-pooling. Edges specify the flow of activations between nodes. The objective f(x) is out-of-sample image classification accuracy. More details can be found in Appendix G.

heading: 5.1 MILP FORMULATIONS FOR GRAPH SEARCH
text: We next introduce a novel MILP formulation that precisely characterizes the set of valid NASBench-101 cells. We use 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, we introduce a new “null” operation, allowing the MILP to represent DAGs with fewer than V nodes. Constraints enforce that all non-null nodes appear on a path from the input to output node, and that there exists at least one such path. A full formulation in terms of linear constraints appears in Appendix G.
A second formulation, nearly identical but enforcing so-called symmetry-breaking 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. An alternative approach uses an isomorphism-invariant surrogate (Wen et al., 2020), but MBO based on such a model is still vulnerable to proposing isomorphic x.

heading: 5.2 RESULTS
text: We refer to an optimizer using the first formulation as NN+MILP and the second as NN+MILP (w/ symmetry). Otherwise, we use the same configuration as in Section 4. Linear+MILP replaces the neural network surrogate with a linear model trained on Dt−1 with randomization provided by bootstrapping that training data. Regularized evolution (RE) and random search (RS) baselines are 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, NN+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 MILP-compatible models at each iteration.

heading: 6 CONCLUSION AND FUTURE WORK
text: In this work we propose the NN+MILP framework for discrete MBO, using neural networks with ReLU activations for surrogate modeling and MILP to solve the acquisition problem. A 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 optimization. Our experiments show that NN+MILP performs well on a range of discrete black-box problems with practical computational overhead using standard packages and hardware.
For future applications to continuous or 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 MILP-representable. For example, Expected Improvement defined over the posterior predictive distribution of an ensemble of neural networks is piecewise-linear.

heading: ETHICS STATEMENT
text: In this paper, we introduce a 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, it is crucial to formulate the objective and constraints of the problem to ensure they reflect ethical design principles. Users should understand whether optimizing for a given objective might implicitly result in biased/discriminatory solutions, or have unintended consequences for some unmodeled objective. For example, in the NAS case study of Section 5 we follow common practice by maximizing out-of-sample accuracy, which might result in larger networks that are environmentally costly to 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.

heading: REPRODUCIBILITY STATEMENT
text: We include detailed implementation sections in the Appendix to accompany all of our experiments. In particular, 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 hyper-parameters we used. Appendix C lists implementation details and hyper-parameter settings for NN+MILP, as well as all baseline algorithms from Sections 4.2 and 4.3 (RegEvo, NN+RegEvo, Ensemble+RegEvo, RBFOpt, ConEvo, NN+ConEvo, NN+RejSample). Finally, Appendix 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.

heading: A STRENGTHENING THE MILP FORMULATION FOR NEURAL NETWORKS
text: Here we discuss more advanced techniques for formulating the neural network surrogate model in the MILP problem. Recall the ReLU formulation constraints (4) and (5) from Section 3.3, except that we consider M separately for each constraint:
0 ≤ y ≤M0α (4’) w>x+ b ≤ y ≤ w>x+ b+M1(1− α), (5’)
Here, we require a nonnegative value M0 such that the right-hand side of (4’) is greater or equal than a valid upper bound on y when α = 1. Similarly, M1 must be a nonnegative value such that the right-hand side of (5’) is greater or equal than zero when α = 0. Therefore, we may choose M0 to be any upper bound of maxx∈Ω′ w>x+ b and M1 to be any upper bound of maxx∈Ω′ −(w>x+ b), where Ω′ 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 M0 or M1, then we may (in fact, must) replace the formulation by y = 0 or y = w>x + b respectively, since in these cases the ReLU is always inactive or active for any x ∈ Ω′. This replacement must be done because the formulation assumes nonnegative M0 and M1 for feasibility.
The simplest way to compute M0 and M1 is to start from bounds in Ω and propagate them via interval arithmetic. For example, if x ∈ [L,U ], thenM0 can be set to ∑ i:wi>0 wiUi + ∑ i:wi<0 wiLi +b
and M1 to −( ∑
i:wi>0 wiLi + ∑ i:wi<0
wiUi + b). However, despite being fast, the drawback of this simple approach is that it does not take into account constraints on Ω or one-hot and no-good constraints.
In our experiments, we compute M0 and M1 by solving the linear programming (LP) relaxations of maxx∈Ω′ w>x + b and maxx∈Ω′ −(w>x + 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 Ω 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 we do not add them. Future work could explore warm-starting the MILP solver using results from earlier MBO iterations.

heading: B BENCHMARKING TASKS
text: 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 fixedlength discrete vectors of length n, with each element drawn from an alphabet A of fixed size.

heading: B.1 TFBIND
text: The objective function is given by the binding affinity of a length-8 DNA sequence to a particular transcription factor, characterized experimentally in the dataset described by Barrera et al. (2016). The problem size is thus fixed by the application at hand, with n = 8 and |A| = 4 (each input element corresponding to a given DNA nucleotide). We min/max-normalize the binding affinity 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 R1, FOXC1 REF R1, GFI1B REF R1, HOXD13 Q325R R1, HOXD13 REF R1, NR1H4 C144R R1, NR1H4 REF R1, PAX4 REF R1, PAX4 REF R2, POU6F2 REF R1, SIX6 REF R1. Here, the 3 fields separated by underscores represent 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.

heading: B.2 RANDOMMLP
text: The objective function is given by the output of a multi-layer perceptron (MLP) with randomlysampled weights. Different functions are generated by varying the architecture type (described below) and random seed. All architectures employ a one-hot encoding of the inputs as the first layer. Weights are sampled using the default behavior of tf.keras.layers.Dense (glorot uniform).
We consider two architecture types, both utilizing more layers/parameters than the 16-neuron networks used by NN+MILP (Section 4). The RandomFCC architecture uses two fully-connected layers with 128 hidden units each, while the RandomCNN architecture uses two convolutional layers each with 64 hidden units each, a kernel width of 13 and stride size of 1. We use a linear activation function for the output and ReLU activations for all intermediate layers.
Unconstrained RandomMLP problems (Section 4.2) all have size n = 25 and |A| = 5 . Eight 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 n = 100 and |A| = 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).

heading: B.3 BBOB
text: The objective is given by a function from the continuous Black-Box Optimization Benchmarking library (Hansen et al., 2009). All BBOB functions are defined for a variable number of dimensions n and the search domain is given as [−5, 5]n, with the global optimum centered at zero. We normalize each function’s output range by evaluating it 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 a point in the grid. Concretely, we use a fixed alphabetA = {1, . . . ,m} for all coordinates, denoting the index of one of m allowed values for that coordinate. Allowed values for each coordinate are m equally-spaced points in the range [−5, 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 on 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 Sector (ATTRAC-
TIVE SECTOR) and Step Ellipsoidal (STEP ELLIPSOID). 3. Functions with high conditioning and unimodal: Discus (DISCUS) and Bent Cigar
(BENT CIGAR). 4. Multi-modal functions with adequate global structure: Weierstrass (WEIERSTRASS) and
Schaffers F7 (SCHAFFERS F7). 5. Multi-model functions with weak global structure: Schwefel (SCHWEFEL) and Gal-
lagher’s Gaussian 21-hi Peaks (GALLAGHER 21ME).
We set the dimension for all of these to n = 10 and discretize as 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 a given coordinate.

heading: C BASELINE OPTIMIZATION ALGORITHMS
text: In this section we describe implementation and configuration details for all baseline optimization algorithms described in Section 4.

heading: C.1 NN+MILP
text: 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., 2016), using the ADAM optimizer for 25K epochs with a batch size of 64 and no explicit regularization. We use a constant learning rate of α = 0.01 and default decay parameters (β1, β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 a time limit of 500 seconds. In order to increase the diversity of trained models, we train each model from scratch at each iteration of optimization instead of fine-tuning a model from an earlier iteration.

heading: C.2 REGEVO
text: 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 recombines two parent sequences. The algorithm proposes xt+1 by selecting two parent sequences from the existing population, recombining them and mutating them. Parents are chosen by tournament selection, taking the two best samples from a randomly-selected subset of size T of previously sampled points. The pool from which parents can be selected is limited to the D most recentlyproposed points (referred to as the “alive population”), to avoid high-reward points from early rounds dominating the process. The selected parent sequences are recombined by copying them left-to-right, starting a pointer at one parent at switching reading to the other parent with a fixed cross-over probability pc after each copy. The resulting sequence is finally mutated by changing each position to a different token from A with a fixed probability pm. In the unconstrained experiments (Section 4.2), we use RegEvo as the outer-loop optimization algorithm and set the tournament size to T = 10, the alive population size to D = 100, and the crossover/mutation probabilities to (pc, pm) = (0.1, 0.1).

heading: C.3 NN+REGEVO
text: This algorithm is an ablation of NN+MILP, with the only difference being the use of RegEvo in lieu of MILP to solve the acquisition problem at every iteration. A surrogate neural network f̂t ∈ F is trained as in NN+MILP, and the acquisition function is a(x) = f̂t(x). The problem of selecting xt+1 is posed as a batched optimization problem and solved by RegEvo.
More concretely, at iteration t, the acquisition function is evaluated for all points in the existing population Dt to generate the initial inner-loop population D̂t := {xi, a(xi)}. 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 blackbox function. That is, RegEvo generates b points by recombination/mutation of parents from D̂t, which are evaluated on the acquisition function and added to the inner-loop population. The process repeats until a 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 xt+1.
In the unconstrained experiments (Section 4.2), we use NN+RegEvo and set surrogate model hyperparameters exactly as in NN+MILP (Section C.1). 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 RegEvo optimizer’s hyper-parameters, defined in Section C.2, are set to T = 20, D = 1, 000 and (pc, pm) = (0.2, 0.01).

heading: C.4 ENSEMBLE+REGEVO
text: We recreate the MBO baseline of Angermueller et al. (2020). Here, surrogate modeling proceeds by optimizing the hyper-parameters of a diverse set of regressor models through randomized search. Regressors are trained using the scikit-learn libary (Pedregosa et al., 2011), drawing from the following model classes (randomized search parameters are listed in parentheses):
• LassoRegressor (alpha) • RidgeRegressor (alpha) • RandomForestRegressor (max depth, max features, n estimators) • LGBMRegressor (learning rate, n estimators)
Each model is evaluated by an explained variance score using five-fold cross validation on the training set. All models with a score ≥ 0.4 are used as an ensemble for the surrogate model, with their average prediction serving as 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 hyper-parameters, defined in Section C.2, are set to T = 20, D = 1, 000 and (pc, pm) = (0.2, 0.01).
We use Ensemble+RegEvo 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 are assigned a highly negative reward (worse than any observed). In this way, the surrogate model might be expected to implicitly model infeasibility with low predictions which should be avoided by the inner-loop optimizer.

heading: C.5 RBFOPT
text: RBFOpt (Costa & Nannicini, 2018) is a black-box optimization solver for mixed-integer unconstrained problems (i.e. with only bound constraints) that performs competitively with respect to other solvers of its type. It uses a Radial Basis Function as a surrogate model and includes a number of practical enhancements. It relies on a mixed-integer nonlinear programming (MINLP) solver, BONMIN (Bonami et al., 2008), to optimize the inner loop problems. The MINLP solver could in theory incorporate constraints in a similar fashion as in our work, although this is not offered by the open-source implementation (aside from manually penalizing the objective function) and we expect it to not scale as well as a MILP solver in practice since MINLP is a significantly more difficult problem class than MILP.
We use RBFOpt for our unconstrained experiments (Section 4.2), using the open-source implementation available at https://github.com/coin-or/rbfopt. We leave all settings at their defaults, including building the initial set of points. We note in particular that the API for this implementation uses an integer encoding for categorical variables (constraints are not supported, which precludes a 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.

heading: C.6 CONEVO
text: 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 a single parent sequence (using the same tournament procedure as RegEvo) and mutates it in a 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 applicationspecific mutator below.
Recall that the domain encodes the selection or not of each of n items using a binary alphabet A = {0, 1}. The items’ indices are partitioned into disjoint, equally-sized subsets S1, . . . , S2k for some k and the constraints enforce that the number of selected items should be the same in pairs of subsets; that is:
∑ i∈S2j−1 I{xi = 1} = ∑ i∈S2j I{xi = 1} ∀j ∈ [k]
where we have used indicator notation and the original decision variables x rather the one-hot encoding.
The mutator begins with a single parent sequence x, assumed feasible, and is given access to the item subsets S1, . . . , S2k. Each pair of subsets (S2j−1, S2j) is mutated concurrently to create the child sequence y, ensuring that mutations to one subset are counter-balanced 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 ∈ 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-1 conversions in positions I+. If c is positive (i.e. there were more 0-to-1 conversions than 1-to-0 conversions) then exactly c indices are chosen randomly from {i ∈ I− : xi = 0}, and also flipped in the child. If c is negative, then −c indices are selected randomly from {i ∈ I− : xi = 1} and flipped in the child. If c is zero, the positions in I− are left unchanged in the child. As a result, the subsets S2j−1 and S2j 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 algorithm setting the tournament size for selecting parent sequences to T = 20 and the mutation probability to pm = 0.05.

heading: C.7 NN+CONEVO
text: This algorithm is an ablation of NN+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 f̂t ∈ F is trained as in NN+MILP, and the acquisition function is a(x) = f̂t(x). Selecting xt+1 is posed as a batched optimization problem and solved by ConEvo. The batched optimization procedure is exactly as described for NN+RegEvo (Section C.3).
In the constrained experiments (Section 4.3), we use NN+ConEvo and set surrogate model hyperparameters exactly as in NN+MILP (Section C.1). 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 T = 20 and pm = 0.05.

heading: C.8 NN+REJSAMPLE
text: This algorithm is an ablation of NN+MILP used in Section 4.3, with the only difference being the use of random search in lieu of MILP to solve the constrained acquisition problem at every iteration. A surrogate neural network f̂t ∈ F is trained as in NN+MILP, and the acquisition function is a(x) = f̂t(x).
To select xt+1, we first generate B points uniformly-at-random from the constrained domain via rejection sampling. That is, candidates are drawn uniformly at random from An and rejected if they do no 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 as xt+1.
In the constrained experiments (Section 4.3) we use NN+RejSample and set surrogate model hyperparameters 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 a rejection rate of ≈ 91% for the subsetequality constraints in our experiments.

heading: D BINARY VS INTEGER VARIABLES
text: 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 our method frequently proposes solutions where several variable values are at either their lower bound or upper bound. As a result, our method would underexplore solutions away from the boundary. A possible explanation for this is that feedforward ReLU networks tend
to extrapolate linearly, and thus their optima may often lie on the boundary (Xu et al., 2021). In contrast, every feasible point of a binary problem lies on a corner of the 0-1 hypercube. A similar observation has been made in the context of IDONE (Bliek et al., 2021), which also uses a ReLU-based surrogate model: encoding the Rosenbrock problem using binary variables improves the performance of the IDONE algorithm, although the opposite happens for a Bayesian optimization algorithm (Karlsson et al., 2020).
We provide computational evidence of this behavior in Figure 3 for TfBind and BBOB instances. The binary variables are encoded as one-hot variables, whereas the integer variables follow an arbitrary ordering for TfBind and the problem ordering for BBOB.

heading: E CONSTRAINED OPTIMIZATION WITH AN ISING PROBLEM
text: To provide further evidence that the results from Section 4.3 can generalize beyond the RandomMLP function, Figure 4 shows the result for the same set of constrained experiments with the only difference that we switch the RandomMLP objective by another random function based on the Ising model. Here, the Ising model represents a fully-connected graph with nodes that can take binary values. Each edge is defined by a 2× 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 a 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+RegEvo with penalty as described in Section 4.3 since it does not perform well with these constraints. We observe that the results are qualitatively similar to those of Section 4.3: NN+MILP continues to perform significantly better than ConEvo, while NN+MILP and NN+ConEvo 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 a single value for the number of paired subsets in the constraints, namely k = 20 when n = 200 or k = 40 when n = 400 (in other experiments, different values of k resulted in qualitatively similar results), and we do not run NN+RejSample. We aggregate over 10 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 it is beneficial to solve the acquisition problem to optimality at larger scales.
Finally, in Figure 6 we examine the impact the surrogate model’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 of the 16 used before. Despite optimizing over a larger domain, neither NN+MILP nor NN+ConEvo show substantial improvements in performance when using a larger
surrogate network. This suggests that, in this case at least, the relatively small number of training points is a more significant bottleneck for approximation than the capacity of the surrogate.

heading: F CONSTRAINED BINARY QUADRATIC PROBLEMS FROM MINLPLIB
text: We study the performance of our method for linearly-constrained binary quadratic problems from the MINLPLib benchmark library (Vigerske, 2021). In practice, one would use a specialized mixedinteger quadratic programming solver to tackle these problems, but they serve well as a benchmark for our black-box optimization method since they are still harder to solve than MILP and include practically-motivated constraints, along with offering good feasible solutions to compare with. Many of these instances are of a slightly larger scale than typical black-box optimization problems (i.e. a few hundreds of variables), which helps us evaluate the method at larger scale and observe its scalability limitations.
We select the instances of type “BQP” from MINLPLib with at least one constraint. We discard the 11 instances prefixed by celar6-sub0, color lab, and max csp, which are large instances for which our method was unable to find a solution with primal gap at most 10% (see below). This
leaves us with the 50 instances listed in Table 1. For consistency with the remainder of the paper, we 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.1). One difference 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 an 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 1 the number of runs that found a 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): if v∗ is the best known value and v is the objective value given by our method, then the primal gap is |v−v
∗| max(|v|,|v∗|) , or 0 if |v| = |v ∗| = 0, or 1 if v 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 a number of instances that are considered to be large for blackbox 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 7, in which we observe how the method scales with the graph partitioning problems denoted by graphpart clique. For the smaller instance (with 60 variables), our method finds an optimal solution in relatively few steps, but it has 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. we have not attempted to tune the parameters or change the architecture of the surrogate model for larger instances).

heading: G NAS-BENCH-101 CASE STUDY

heading: G.1 BACKGROUND
text: In Section 5 we consider the NAS-Bench-101 (Ying et 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 acyclic
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 a directed path, while the remaining nodes are each assigned to be 1x1 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-ofsample 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, f(x) 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 to halt unpromising evaluations.
Despite these differences, we apply NN+MILP exactly as described in Section 3/Appendix C.1 (including no-good constraints which prevent repeated queries). We reimplement the Regularized Evolution (RE) and random search (RS) baselines from the original NAS-Bench-101 paper (Ying et al., 2019), using the same hyper-parameters settings specified therein.
The NAS-Bench-101 dataset contains pre-computed validation and test accuracies for three independently 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 budget of evaluations, we use a fixed budget of T = 5 × 106 seconds, and allow algorithms to query the objective until cumulative training time exceeds the budget. For evaluation purposes (e.g. Figure 2c) we plot the out-of-sample accuracy of the incumbent architecture as a function of cumulative architecture training time.

heading: G.2 DOMAIN FORMULATION
text: To formulate the NAS-Bench-101 domain, we first define a representation of cell architectures as fixed-length binary vectors. We split the representation into two components; one set 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 acyclic, we limit the edge variables to the strict upper triangle of the adjacency matrix, which implicitly enforces a 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 = {conv1x1,conv3x3,maxpool3x3}. To ensure a fixed-length set of decision variables while allowing for graphs with a 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
representations, constrained to yield feasible cell architectures. Denoting by V andM the maximum number of allowable nodes and edges respectively, the decision variables (all binary) are:
• mi,j for 1 ≤ i < j ≤ V , 1 if there is an edge from node i to node j, 0 otherwise. • wi,k for 1 ≤ i ≤ V, 1 ≤ k ≤ |S|, 1 if node i is assigned the k’th operation in S , 0
otherwise.
• zi for 1 < i < V , 1 if node i is assigned the null operation,0 otherwise.
The feasible set of cell architectures can then be given in terms of linear constraints as follows:
w1,k = wV,k = z1 = zV = 0 for 1 ≤ k ≤ |S| (1)
zi + S∑ k=1 wi,k = 1 for 1 < i < V (2)
V∑ i=1 V∑ j=i+1 mi,j ≤M (3)
mi,j ≤ 1− zj for 1 ≤ i < j ≤ V (4) mi,j ≤ 1− zi for 1 ≤ i < j ≤ V (5)
j−1∑ i=1 mi,j ≥ 1− zj for 1 ≤ j ≤ V (6)
V∑ j=i+1 mi,j ≥ 1− zi for 1 ≤ i ≤ V (7)
zi ≤ zi+1 for 1 < i < V − 1 (8)
Constraints 1 ensure that the input and output nodes are not assigned any operation from S or null, while 2 enforces the one-hot encoding of operations for intermediate nodes (including the possibility of a null operation). Constraint 3 imposes a limit on the number of edges in the graph, per the NASBench-101 specifications. Constraints 4 & 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-triangular adjacency matrix, these also ensure that there is always a path from the input to 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 a higher-indexed non-null node – and all non-null nodes (including the output) have at least oneincoming 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-breaking constraints. These assert that a 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, it introduces a “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 (Figure 2c). We hypothesize that this is due to a reduction in the exploration behaviour of NN+MILP, as the surrogate’s predictive distribution 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 augmentDt with isomorphic representations before training, e.g. by random reordering of nodes in the representations of sampled points.

heading: H ADDITIONAL PLOTS
text: We present in this section additional figures related to the experiments in the main text.

heading: H.1 UNCONSTRAINED OPTIMIZATION
text: While the best observed reward inDN (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, we 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, we min/max normalize algorithms’ AUC scores within each problem exactly as we did for the best observed reward (Section 4.1). That is, the best (resp. worst) on-average algorithm in terms of AUC is assigned a score of one (resp. zero) and intermediate values express relative distance from these extremes.
Figure 8 plots the distribution of algorithms’ 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 1). Figures 11 and 12 plot the individual reward curves as function of outer-loop iteration for each unconstrained problem.

heading: H.2 CONSTRAINED OPTIMIZATION
text: In Figure 9 we plot the distribution of algorithms’ normalized final reward scores across all constrained problems from Section 4.3 (paralleling Figure 1, top). The individual reward curves for all such problems can be found in Figures 13 and 14, which did not differ significantly.

heading: H.3 PRACTICALITY OF MILP
text: Here we present experiments to explore the impact of surrogate network size on MILP solve times. We varied NN+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 a simple Linear model (no hidden layer). Note that the first architecture, FC(16), was what was presented in the main results. Otherwise, we used the same training and optimization hyper-parameters for NN+MILP as described in Section C.1. We ran this experiment on all 12 unconstrained TfBind problems from Section 4.2, with 20 trials using different random initial data sets.
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 out 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.

Review summaries:
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.   The experiments cover both the cases of unconstrained and constrained optimization of the acquisition function.  The 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.  In 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. Finally, 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.

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.

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.          

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.