Black-box Combinatorial Optimization using Models with Integer-valued Minima

When a black-box optimization objective can only be evaluated with costly or noisy measurements, most standard optimization algorithms are unsuited to find the optimal solution. Specialized algorithms that deal with exactly this situation make use of surrogate models. These models are usually continuous and smooth, which is beneficial for continuous optimization problems, but not necessarily for combinatorial problems. However, by choosing the basis functions of the surrogate model in a certain way, we show that it can be guaranteed that the optimal solution of the surrogate model is integer. This approach outperforms random search, simulated annealing and one Bayesian optimization algorithm on the problem of finding robust routes for a noise-perturbed traveling salesman benchmark problem, with similar performance as another Bayesian optimization algorithm, and outperforms all compared algorithms on a convex binary optimization problem with a large number of variables.


Introduction
Traditional optimization techniques make use of a known mathematical formulation of the objective function, for example by calculating the derivative or a lower bound. However, many objective functions in real-life situations have no known mathematical formulation. For example, smart grids or railways are complex networks where every decision influences the whole network. In such applications, we can observe the effect of decisions either in real life, or by running a simulation. Waiting for such a result can take some time, or may have some other cost associated with it. Furthermore, the outcome of two observations with the same decision variables may be different. Such problems have been approached using methods such as black-box or Bayesian optimization [1], simulation-based optimization [2], and derivative-free optimization [3].
Here, a model fits the relation between decision variables and objective function, and then standard optimization techniques are used on the model instead of the original objective. These so-called surrogate modeling techniques have been applied successfully to continuous optimization problems in signal processing [4], optics [4], machine learning [5], robotics [6], and more. However, it is still an on-going research question on how these techniques can be applied effectively to combinatorial optimization problems. A common approach is to simply round to the nearest integer, a method that is known to be sub-optimal in traditional optimization, and also in black-box optimization [7]. Another option is to use discrete surrogate models from machine learning such as regression trees [8] or linear model trees [9]. Although powerful, this makes both model fitting and optimization computationally expensive.
This work describes an approach where the surrogate model is still continuous, but where finding the optimum of the surrogate model gives an integer solution. The main contributions are as follows: • This surrogate modeling algorithm, called IDONE, with two variants (one with a basic and one with a more complex surrogate model).
• A proof that finding the optimum of the surrogate model gives an integer solution.
• Experimental results that show when IDONE outperforms random search, simulated annealing and Bayesian optimization.
Section 2 gives a general description of the problem and an overview of related work. Section 3 describes the IDONE algorithm and the proof. In Section 4, IDONE is compared to random search, simulated annealing and Bayesian optimization on two different problems: finding robust routes in a noise-perturbed traveling salesman benchmark problem, and a convex binary optimization problem. Finally, Section 5 contains conclusions and future work.

Problem description and related work
Consider the problem of minimizing an objective f : R d → R with integer and bound constraints: These bounds are also assumed to be integer. It is assumed that f does not have a known mathematical formulation, and can only be accessed via noisy measurements y = f (x) + , y ∈ R, with ∈ R zero-mean noise with finite variance. Furthermore, taking a measurement y is computationally expensive or is expensive due to a long measuring time, human decision making, deteriorating biological samples, or other reasons. Examples are hyper-parameter optimization in deep learning [10], contamination control of a food supply chain [11], and structure design in material science [12].
Although many standard optimization methods are unfit for this problem, there exists a vast number of methods that were designed with most of the above assumptions in mind. For example, local search heuristics [13] such as hill-climbing, simulated annealing, or taboo search are general enough to be applied to this problem, and have the advantage of being easy to implement. These heuristics are considered as the baseline in this work, together with random search (simply assign the variables completely random values and keep the best results). Population-based heuristics such as genetic algorithms [14], particle swarm optimization [15], and ant colony optimization [16] operate in the same way as local search algorithms, but keep track of multiple candidate solutions that together form a population. These algorithms have been applied successfully in many applications, but are unfit for the problem described in this paper since evaluating a whole population of candidate solutions is not practical if each measurement is expensive. The same holds for algorithms that filter out the noise by averaging, such as COMPASS [17], since they evaluate one candidate solution multiple times. Typical integer programming algorithms such as branch and bound [18] work well on standard combinatorial problems, but when the objective function is unknown and observations of single values are noisy or expensive, it is very difficult to obtain reasonable bounds.
Surrogate modeling techniques operate in a different way from the above methods: past measurements are used to fit a model, which is then used to select a next candidate solution. Bayesian optimization algorithms [1,19], for example, have been successfully applied in many different fields. These methods use probability theory to determine the most promising candidate point according to the surrogate model. However, when the variables are discrete, the typical approach is to relax the integer constraints, which often leads to sub-optimal solutions [7]. The authors in [7] tackled this problem by modifying the covariance function used in the surrogate model. Another approach, based on semi-definite programming, is given in [11]. And the Hy-perOpt algorithm [20] takes yet a different approach by using a Tree-structured Parzen Estimator as the surrogate model, which is discrete in case the variables are discrete. HyperOpt is considered the main contender in this paper.
A downside of many Bayesian optimization algorithms is that the computation time per iteration scales quadratically with the number of measurements taken up until that point. This causes these methods to become slower over time, and after a certain number of iterations they may even violate the assumption that the bottleneck in the problem is the cost of evaluating the objective. This downside is recognized and overcome in two similar algorithms: COMBO [12] and DONE [4]. Both algorithms use the power of random features [21] to get a fixed computation time every iteration, but COMBO is designed for combinatorial optimization while DONE is designed for continuous optimization. A disadvantage of COMBO is that it evaluates the surrogate model at every possible candidate point, a number that grows exponentially with the input dimension d. Though evaluating the surrogate model takes very little time (compared to evaluating the original objective f ), this still makes the algorithm unfit for problems where the input dimension d is large. A variant of DONE named CDONE has been applied to a mixed-integer problem, where the integer constraints were relaxed [10], but as mentioned earlier, this can lead to sub-optimal solutions.
However, the downside of having to relax the integer constraints can be circum-vented. By choosing the basis functions in a certain way, we show how a model can be constructed for which it is known beforehand that the minima of the model lie exactly in integer points. This makes it possible to apply the algorithm to combinatorial problems, as explained in the next section.

IDONE algorithm
The DONE algorithm [4] and its variants are able to solve problem (1) without the integer constraint by making use of a surrogate model. Every time a new measurement y = f (x) + comes in, the surrogate model is updated, the minimum of the surrogate model is found, and an exploration step is performed. To make the algorithm suitable for combinatorial optimization we propose a variant of DONE called IDONE (integer-DONE), where the surrogate model is guaranteed to have integer-valued minima.

Piece-wise linear surrogate model
The proposed surrogate model g : R d → R is a linear combination of rectified linear units ReLU(z) = max(0, z), a basis function that is commonly used in the deep learning community [22]: Unlike what is common practice in the deep learning community, the parameters w k and b k remain fixed in this surrogate model. This makes the model linear in its parameters (c k ), allowing it to be trained via linear regression instead of iterative methods. This is explained in Section 3.2.
Because of the choice of basis functions, the surrogate model is actually piece-wise linear, which causes its local minima to lie in one of its corner points: The reverse of this theorem is not necessarily true: ifx satisfies z k (x) = 0 for d linearly independent z k , then it depends on the parameters c k of the model whetherx is actually a local minimum or not.
The number of local minima and their locations depend on the parameters w k and b k . In this work, we provide two options for choosing these parameters in such a way that the local minima are always found in integer solutions. In the first case, the functions z k are simply chosen to have zeros on hyper-planes that together form an integer lattice: Definition 1 (Basic model ). Let g be as in (2). The parameters w k , b k of the basic model are chosen according to Algorithm 1. That is, every function z k is zero on a Algorithm 1 Basic model parameters

An example of a basis function in this model
The 1 comes from the model bias, a basis function that is equal to 1 everywhere. This allows the model to be shifted up or down.
Since all the basis functions depend only on one variable, this basic model might not be fit for problems where the decision variables have complex interactions. Therefore, in the advanced model, we use the same basis functions, but we also add basis functions that depend on two variables: Definition 2 (Advanced model). Let g be as in (2). The parameters w k and b k of the advanced model are chosen according to Algorithm 2. That is, every function z k different from the ones in the basic model is zero on a

This model has
more basis functions than the basic model. The added functions z k are zero on diagonal lines through two variables, see Figure 1.
We now show one of our main contributions.
(III) If x * is a non-strict local minimum of the basic model, it holds that the model retains the same value when going from x * to the nearest pointx that satisfiesx ∈ Z d and l i ≤x i ≤ u k , ∀i = 1, . . . , d.

Algorithm 2 Advanced model parameters
To arrive at a contradiction, suppose that ∃s ∈ {1, . . . , d} such that x * s ∈ Z. Then by (3), x * i ∈ Z for all i ∈ {1, . . . , d}. However, this is only possible if none of the z k have the form z k (x) = ±(x i − j), and all d of the z k have the form z k (x) = ±(x i − x i−1 − j), for d different i. But by construction, there are only d − 1 of these last ones available, see Algorithm 2 (the for-loop starts at 2). Therefore, it is not true that ∃s ∈ {1, . . . , d} such that x * s ∈ Z. Hence, x * i ∈ Z ∀i = 1, . . . , d, which is what the theorem claims.
(III) Let x * be a non-strict local minimum of the basic model and suppose x * s ∈ Z ∪ [l s , u s ] for some s ∈ {1, . . . , d}. Let L be the line segment from x * s to the nearest pointx s in the set Z ∪ [l s , u s ], without including that point. Since the only z k functions that depend on x s have the form z k (x) = ±(x s − j), j = l s , . . . , u s , it follows that z k (x) = 0 does not happen on L for any z k that depends on x s . Therefore, model g is linear on this line segment, and since x * is a non-strict local minimum and g is continuous, g retains the same value when replacing x * s byx s . This can be repeated for all s for which x * s ∈ Z ∪ [l s , u s ], which proves the claim. (IV) Let x * be a non-strict local minimum of the advanced model and suppose x * ∈ Z. We first show that rounding x * to the nearest integer does not change the sign of any z k . Note that all the z k of the advanced model have the form , for some i = 1, . . . , d and some integer j. Letx i denote rounding x i to the nearest integer. Then we have (because j is integer): Since the sign of none of the z k change when rounding, and model g is only nonlinear when going from z k (x) < 0 to z k (x) > 0 for some k = 1, . . . , D, it follows that g is linear on the line segment from x * to the nearest integer. Together with the fact that x * is a non-strict local minimum, it follows that g retains the same value on this line segment. Finally, the claim is valid because g is continuous.

Fitting the model
Because the surrogate model g is linear in its parameters c k , fitting the model can be done with linear regression. Given input-output pairs (x i , y i , i = 1, . . . , N ), this is done by solving the regularized linear least squares problem with regularization parameter λ and initial weights c 0 . The regularization part is added to overcome ill-conditioning, noise, and model over-fitting. Furthermore, by choosing c 0 = [0, 1, . . . , 1] T , it is ensured that the surrogate model is convex before the first iteration [10]. In this work, λ = 0.001 has been chosen.
To prevent having to solve this problem at every iteration (with runtime O(N 3 )), (4) is solved with the recursive least squares algorithm [23]. This algorithm has runtime O(D 2 ) per iteration, with D the number of basis functions used by the model. This implies that the computation time per iteration does not depend on the number of measurements, which is a big advantage compared to Bayesian optimization algorithms (which usually have complexity O(N 2 ) per iteration). The memory is also O(D 2 ), because a D × D covariance matrix needs to be stored. Since D scales linearly with the input dimension d and with the lower and upper bounds, the computational complexity of fitting the surrogate model is O(p 2 d 2 ), with p = max i (u i − l i ).

Model visualization
To visualize the surrogate model used by the IDONE algorithm, the fitting procedure is applied to a simple traveling salesman problem with four cities. The distance matrix for the cities is shown in Table 1. The decision variables are chosen as follows: the route starts at city 1, then variable x 1 ∈ {1, 2, 3} determines which of the three remaining cities is visited, then variable x 2 ∈ {1, 2} determines which of the two remaining cities is visited; then the one remaining city is visited, then city 1 is visited again. This problem has two optimal solutions: x = [1, 2] T (route 1-2-4-3-1) and x = [2, 2] T (route 1-3-4-2-1), both with a total distance of 80. All other solutions have a total distance of 95. Figure 2 shows what the surrogate model looks like after taking measurements in all possible data points for this problem, which is possible due to the low number of possibilities. It can be observed that this model is piece-wise linear and that any local minimum retains the same value when rounding to the nearest integer. Furthermore, the diagonal lines (see also Figure 1) make the advanced model more accurate.

Finding the minimum of the model
After fitting the model g at iteration N , the algorithm proceeds to find a local minimum using the new weights c N : The BFGS method [24] with a relaxation on the integer constraint was used to solve the above problem, with a provided analytical derivative of g. In this work, the derivative of the basis function ReLU(z) = max(0, z) has been chosen to be 0.5 at z = 0. The optimal solution was rounded to the nearest integer per Theorem 2.

Exploration
After fitting the model and finding its minimum, a new point x N +1 needs to be chosen to evaluate the function f . As in DONE [4], a random perturbation δ is added to the found minimum: x N +1 = x * + δ, but instead of a continuous random variable, δ ∈ {−1, 0, 1} d is a discrete random variable with the following probabilities: In this work, p = 1/d has been chosen (d is the number of variables).

IDONE algorithm
The IDONE algorithm iterates over three phases: updating the surrogate model with recursive least squares, finding the minimum of the model, and performing the exploration step. The pseudocode for the algorithm is shown in Algorithm 3. Depending on which subroutine is used in the first line, we refer to this algorithm as either IDONEbasic (using the basic model) or IDONE-advanced (using the advanced model).

Experimental results
To determine how the IDONE algorithm compares to other black-box optimization algorithms in terms of convergence speed and scalability, it has been applied to two problems: finding a robust route for a noise-perturbed asymmetric traveling salesman benchmark problem with 17 cities, and an artificial convex binary optimization problem. The first problem gives a first indication of the algorithm's performance on an objective function that follows from a simulation where there is a network structure. The second problem shows an easier and more tangible situation -due to the convexity and the fact that we know the global optimum -which makes it easier to interpret results. The algorithm is compared with several different black-box optimization algorithms: random search (RS), simulated annealing (SA), and two different Bayesian optimization algorithms: Matlab's bayesopt function [25] (BO), and the Python library HyperOpt [20] (HypOpt). The RS, IDONE-advanced and HypOpt algorithms are implemented in Python and run on a cluster (32 Intel Xeon E5-2650 2.0 GHz CPUs), while the SA, BO and IDONE-basic algorithms are implemented in Matlab and run on a laptop (Intel Core i7-6600U 2.6 GHz CPU with 8 GB RAM), which is about 34% faster 1 according to https://www.cpubenchmark.net/singleThread.html. Whenever we compare the runtimes of algorithms evaluated in these two different machine environments we will be more careful with our conclusions. For BO and HypOpt, we used the default settings. It should be noted that BO and HypOpt are both aimed at minimizing black-box functions using as few function evaluations as possible. For SA, the settings are explained below.
In the context of the IDONE algorithm, the SA algorithm essentially consists of just the exploration step of the IDONE algorithm (see Section 3.4), coupled with a probability of returning to the previous candidate solution. Suppose the current best solution is (x b , y b ), and that the exploration step as defined in Section 3.4 gives a new candidate solution (x c , y c ). If y c < y b , then x c is accepted as the new best solution. Else, there is a probability that x c is still accepted as the new best solution. This probability is equal to e (y b −yc)/T , with T a so-called temperature. In this work, the simulated annealing algorithm starts out with a starting temperature T = T 0 , and the temperature is multiplied with a factor T f every iteration. This strategy is called a cooling schedule. For the asymmetric traveling salesman problem, T 0 = 4.48 and T f = 0.996 have been chosen. For the convex binary optimization problem, T 0 = 1 and T f = 0.95 have been chosen.

Robust routes for an asymmetric traveling salesman problem (17 cities)
Consider the asymmetric TSP benchmark called BR17. This benchmark was taken from the TSPLIB website [26], a library of sample instances for the traveling salesman problem. While there exist specific solvers developed for this problem, these solvers are not adequate if the objective to be minimized is perturbed by noise. Here, noise ∈ [0, 1], with a uniform distribution, was added to the distances between any two cities (for distances other than 0 or infinity, which both occurred once per city, the mean distance between cities is 16.43 for this instance). Furthermore, every time a sequence of cities has been chosen, we evaluate this route 100 times, with different noise samples. The objective is the worst-case (largest) route length of the chosen sequence of cities. Minimizing this objective should then result in a route that is robust to the noise in the problem. For the variables the same encoding as in Section 3.2.1 has been used, giving 15 integer variables in total.
All algorithms were run 5 times on this problem, and the results are shown in Figure 3. The differences between the run times of the algorithms are significantly larger than between the two machine environments, so we ignore this subtlety here. The BO algorithm was not included as it took over 80 hours per run. It can be seen that IDONE-advanced achieves similar results as HyperOpt, although its computation time is over ten times larger. Both HyperOpt and IDONE-advanced outperform the simpler benchmark methods. It seems IDONE-basic is unable to deal with the complex interaction between the variables due to the basic structure of the model, as it performs worse than the baseline algorithms.

Convex binary optimization
To gain a better understanding of the different algorithms, the second experiment is done on a function with a known mathematical formulation. Consider the function with A a random positive semi-definite matrix, and x * ∈ {0, 1} d a randomly chosen vector, with d the number of binary variables. The optimal solution x * or the structure of the function is not given to the different algorithms, only the number of variables and the fact that they are binary. Starting from a matrix U where each element is randomly generated from a uniform [0, 1] distribution, matrix A is constructed as with I the identity matrix. The function f can only be accessed via noisy measurements y = f (x) + , with ∈ [0, 1] a uniform random variable. We ran 100 experiments with this function, with IDONE and the other black-box optimization algorithms. For each run, A and x * were randomly generated, as well as the initial guess x 0 . All algorithms were stopped after taking 1000 function evaluations, and the best found objective value was stored at each iteration. Figure 4 shows a convergence plot for the case d = 100. It can be seen that the two variants of IDONE have the fastest convergence. The large number of variables is too much for a pure random search, but also for HyperOpt. Simulated annealing still gives decent results on this problem. Figure 5 shows the final objective value and computation time after 1000 iterations for the same problem for different values of d. The number of variables d was varied between 5 and 150. The Bayesian optimization implementation of Matlab was only evaluated for d = 5 due to its large computation time. As can be seen, IDONE (both versions) is the only algorithm that consistently gives a solution at or close to the optimal solution (which has an objective value between 0 and 1) for the highest dimensions. Where all algorithms get at or close to the optimal solution for problems with 10 variables or less, the difference between the algorithms becomes more distinguishable when 20 or more variables are considered. The strengths of HyperOpt, such as dealing with different types of variables that can have complex interactions, are not relevant for this particular problem, and the Parzen estimator surrogate model does not seem to scale well to higher dimensions compared to the piece-wise linear model used by IDONE.

Conclusions and future work
The IDONE algorithm is a black-box optimization algorithm that is designed for combinatorial problems with binary or integer constraints, and had shown to be useful in particular when the objective can only be accessed via costly and noisy evaluations. By using a surrogate model that is designed in such a way that its optimum lies in an integer solution, the algorithm can be applied to combinatorial optimization problems without having to resort to rounding in the objective function. IDONE has a fixed computation time per iteration that scales quadratically with the number of variables but is not influenced by the number of times the function has been evaluated, which is an advantage compared to Bayesian optimization algorithms. One variant of the proposed algorithm, IDONE-advanced, has been shown to outperform random search and simulated annealing on the problem of finding robust routes in a noise-perturbed traveling salesman benchmark problem, and on a convex binary optimization problem with up to 150 variables. The other variant of the algorithm, IDONE-basic, is a lot faster, but only performed well in the second experiment. HyperOpt, a popular Bayesian optimization algorithm, performs similar as IDONE-advanced on the traveling salesman benchmark problem, but does not scale as well on the binary optimization problem. These results show that there is room for improvement in the use of surrogate models for black-box combinatorial optimization, and that using continuous models with integer-valued local minima is a new and promising way forward.
In future work, the special structure of the surrogate model will be further exploited to provide a faster implementation, and the algorithm will be tested on real-life applications of combinatorial optimization with expensive objective functions. The question also arises whether this algorithm would perform well in situations where the objective function is not expensive to evaluate, or does not contain noise. Population-based methods perform particularly well on cheap black-box objective functions, so it would be interesting to see if they could be combined with the surrogate model used in this paper. As for the noiseless case, it is known that for continuous variables it becomes easy in this case to estimate the gradient and use more traditional gradient-based methods, but in the case of discrete variables the traditional combinatorial optimization methods might still benefit from IDONE's piece-wise linear surrogate model. Where surrogatebased optimization techniques have had great success in continuous optimization problems from many different fields, we hope that this work opens up the route to success of these techniques for the plenty of open combinatorial problems in these fields.