Goal Seeking Quadratic Unconstrained Binary Optimization

Quadratic Unconstrained Binary Optimization (QUBO) is a popular modeling framework wherein many combinatorial optimization problems can be recast into the form $min\;x^{\prime}Qx;x\in\{0,1\}$ where $Q$ is a symmetric matrix of size $n$ with integer or real components and $x$ is a binary vector (see [1] for a survey of the broad applicability of this modeling structure). Optimization software such as CPLEX and Gurobi find an exact solution that optimizes the objective function value. Constraint Programming (CP) finds feasible solutions satisfying a set of constraints. Thus, CP supports a decision maker who is interested in setting goals and generating a set of solutions to consider for implementation, for example, a target time to complete a task, multi-criteria profit, risk, or time-frame targets. In the mixed integer linear programming paradigm, the objective target is modeled as lower and upper bound constraints on the objective function. However for quadratic binary objectives, these constraints are also binary quadratic and are very difficult for conventional linear solvers such as CPLEX and Gurobi. In this short paper, we present two variants of a goal-seeking or target QUBO. In the first case, the target $t$ is exact and we are interested in finding binary vectors $x$ such that $x^{\prime}Qx=t$. In the second case, the target is uncertain but defined by an interval $[lb,ub]$.

Target QUBO belongs to the class of multi-criteria goal programming where the objective function $x^{\prime}Qx$ is provided a goal or target to be achieved and deviation from the target may be seen as the quadratic Taguchi quality loss function ([2]), although other measures such as absolute value are valid. Our implementation utilizes the squared deviation from the goal, which is minimized in the tabu-search heuristic. However, the importance of the magnitude of deviation from the target is typically determined by the decision-maker.

Our proposed model finds utility in four scenarios. First, the settings where we are interested in finding the input parameters that achieve a desired target level benefit from this analysis which is closely related to the “goal-seeking problem” defined in [3] and exemplified in the parameter selection method of Taguchi ([4]). For instance, near-optimal solutions in RNA folding prediction are often better predictors of RNA structure than the optimal, and the near-optimal solutions occur in the neighborhood around a target energy level ([5]). Second, in sensitivity analysis of a discrete event simulation (see [3] for more details). The third application arises in the context of the Constraint Satisfaction Problem modeled as a QUBO wherein a target number of constraints are specified for satisfaction.

Fourth, target QUBO can be utilized as a tool to solve multi-objective optimization problems. For instance, when the model appears as a subproblem for lexicographic optimization and the objectives are sorted from most important to least important ([6]). Next, we find the binary variables satisfying the inequality for the $i^{\prime}th$ objective $x^{\prime}Q_{i}x\leq f_{i}^{*}+\delta_{i}$ which guarantees that the search focuses on near-optimal function evaluation $f_{i}^{*}$ controlled by a tolerance level $\delta_{i}$ which translates to $x^{\prime}Q_{i}x\in[f_{i}^{*},f_{i}^{*}+\delta_{i}]$ in our modeling framework. Thus, the results are useful for lexicographic optimization wherein we are looking for solutions that are close to the optimal value. Moreover, various objectives can be combined into a weighted multi-objective function. Thus, the various solution techniques proposed for multi-objective QUBO (like [7, 8, 9, 10]) benefit from our analysis.

The paper is organized as follows. The goal-seeking QUBO model is presented in Section 2. The solution technique and the results are presented in Section 3 and 4, followed by conclusions and future research directions.

First, we present the exact goal-seeking QUBO for a deterministic target $t$. An exact model with quadratic penalties minimizing the squared deviation from the target is $min\;(x^{\prime}Qx-t)^{2}=min\;(x^{\prime}Qx)^{2}-2t(x^{\prime}Qx)+t^{2}$ wherein the first term has fourth-degree pseudo boolean functions of the form $x_{i}x_{j}x_{k}x_{l}$ significantly increasing the computational complexity. We could iteratively use Rosenberg quadratization substituting the quadratic subterm $x_{i}x_{j}$ with an additional binary variable $z_{ij}$ thereby introducing a penalty term $M(x_{i}x_{j}-2x_{i}z_{ij}-2x_{j}z_{ij}+3z_{ij})$ for a minimization problem. In this way, we could transform the fourth degree polynomial to a quadratic polynomial and utilize existing QUBO solvers. However, this reformulation leads to a large number of auxiliary variables and penalty terms ([11]).

While the number of auxiliary variables required for the transformation of a fourth degree goal-seeking QUBO is O($n^{2}$), calculating the impact on the achievement function $(x^{\prime}Qx-t)^{2}$ as part of a one-flip search routine is fast given that the incremental function evaluation of $x^{\prime}Qx$ could be done in $O(1)$ ([12]). Thus, selecting the best variable out of one-flip neighborhood is $O(n)$. Only binary vectors that meet the target are collected as part of the search routine and presented to the decision-maker.

A common feature of integer optimization problems is the presence of alternate optima ([13]) and while it may not be practical to enumerate all alternate optima, our solution approach generates a list of binary solution vectors meeting the target value or range. This feature is important to multi-criteria decision-making, wherein the decision-maker selects the preferred solution from the list of pareto optimal solution vectors.

The landscape theory provides us techniques to quantify the impact on the number of local optima due to width of the target interval. According to [14], the number of local optima $M$ in a combinatorial landscape is given by the following expression:

where $|X|$ denotes the total size of the solution space and $X(x_{o},l)$ represents the set of solutions accessible from the vector $x_{o}$ in $l$ or less local moves. Note that the width of the interval around $x_{o}$ increases as $l$ increases. An estimate of $M$ for QUBO of size $n$ provided by [15] is given by:

Using this estimate, we could derive the following related to Lemma 2 in [16]:

Second, we describe the interval goal-seeking QUBO that seeks to find solution vectors that satisfy $x^{\prime}Qx\in[lb,ub]$. This is attained by optimizing the achievement function $min\;(x^{\prime}Qx-lb)(x^{\prime}Qx-ub)=min\;(x^{\prime}Qx)^{2}-(lb+ub)(x^{\prime}Qx)+lb*ub$. Again we face computational difficulties with solving this exact fourth-degree pseudo boolean function involving a large number of auxiliary variables and penalty terms. Alternatively, lower and upper bound constraints on the objective function can be used, but most solvers do not work well with quadratic binary constraints. Thus we integrate the target constraint into the solution method as described in Section 3.

Because of the quadratic nature of the achievement function, the optimal value is achieved at $x^{\prime}Qx=(lb+ub)/2$. However, in general the achievement function is non-positive as long as $x^{\prime}Qx\in[lb,ub]$. Moreover, the binary vectors outside this interval lead to a higher penalty term, which increases with the squared distance from the interval bounds. This incentivizes the solver to limit the search to the solution landscape around the interval. Any binary vectors with non-positive function evaluation are easily collected during the search routine. Such vectors satisfying the interval constraint are presented to the decision-maker in descending order according to their objective function.

The following simple lemma captures the relationship between satisficing vectors and achievement function for either of the variants of goal-seeking QUBO.

Our solution method is summarized in Algorithm 1. We utilize a tabu-search (TS) heuristic for finding the binary vectors that meet the user-defined goals. A neighborhood move is defined by a one-flip i.e. setting $x_{i}$ to $1-x_{i}$ for a specific variable $i$. The best impact on the objective function $x^{\prime}Qx$ due to one-flip is calculated in an incremental manner building on a previous solution requiring $O(n)$ ([12]) to select the best one-flip from the $n$ variables in the neighborhood. We start from the all-zero solution and iteratively choose the non-tabu move that leads to the best improvement in the achievement function ($AF$). In case an improvement in $AF$ is not possible i.e., we reach a local optima, we chose the non-improving and non-tabu move having the least impact on $AF$. As discussed earlier, $AF$ is defined as $(x^{\prime}Qx-t)^{2}$ for deterministic target and $(x^{\prime}Qx-lb)(x^{\prime}Qx-ub)$ for interval target.

To help prevent cycling to previously visited solutions, we define a tabu list consisting of the tabu tenure of the variables that were flipped. Thus, a non-tabu variable $x_{i}$ is added to the tabu list with maximum tabu tenure if it was the best choice for one-flip. The variable is eventually removed from the tabu list after a fixed number of iterations set by the tabu tenure parameter. Note that this parameter has an impact on the computational performance of the heuristic. In general, as the tabu tenure increases, the chances for cycling reduces while a too large tabu tenure parameter inhibits the ability to reach a target.

Compared to a traditional heuristic wherein we are interested in finding the best solution in a limited runtime, we focus on finding multiple solutions that meet the user-defined objective goal. In this way, we can provide various choices to the decision-maker. The termination criteria in our TS heuristic is the solution runtime and our algorithm might find duplicate solutions due to cycling, so the end of the routine, we run a post-processing method that identifies the unique solutions among the solution set.

The solution quality of target QUBO is assessed by the number of unique solutions obtained by the solver in a limited runtime. Because the tabu tenure parameter has an impact on the solution quality, we conducted preliminary experiments to determine a good setting. For this purpose, we utilized the $2500$ node ORLIB instances ([18]). The deterministic target was set at $80\%$ of the best known solutions of the problems (obtained from [19]). In general, smaller tabu tenures lead to quick solutions, but tend to cycle, and the solver is often stuck in local optima. On the other hand, large tabu tenure restricts the neighborhood visited by the solver. Based on preliminary experiments, we observed that tabu tenure between $10$ and $20$ provides good results on the ORLIB problems and we set the tabu tenure to $10$ for the problem set in the following section.

For computational experiments, we use the QUBO instances with $2500,3000$ and $4000$ variables (presented in [18] and [20]). The solution technique was implemented using C. We also utilize the IBM Constraint Programming (CP) Optimizer for benchmarking with the exact model. The exact non-linear model for the two variants namely $x^{\prime}Qx=t$ and $x^{\prime}Qx>=lb\;\&\;x^{\prime}Qx<=ub$ was input to the CP solver using Python DOCPLEX API. The target level is set as a higher percentage of best known solutions ($bks$) namely $80\%,85\%,90\%$ and $95\%$. The experiments were performed on a 3.40 GHz Intel Core i7 processor with 16 GB RAM running 64 bit Windows 7 OS.

In terms of benchmarking, the IBM CP solver was unable to find any feasible solution for either of the problem variants in $100$ seconds and for bigger instances, the size of the exact model was too large in terms of memory requirements. We also tried the default portfolio of solvers provided by MiniZinc ([21]) for the non-linear models, and those also yielded no solutions. Thus, these problems are computationally challenging.

First, we present the results using our heuristic for the exact target QUBO problem $x^{\prime}Qx=t$ in Table 1. It is evident from the results that as we get closer to $bks$, the goal-setting problem becomes more challenging and we conclude that the solution landscape for these problems are distributed such that the number of feasible solutions decrease as the target approaches $bks$.

Second, we investigate the performance of our heuristic on interval targets i.e. $x^{\prime}Qx\in[lb,ub]$ in Table 2. For this case, we choose the intervals as $(80\%,85\%)$,$(85\%,90\%)$,$(90\%,95\%)$ and $(95\%,100\%)$ of $bks$. Note that $(lb,ub)$ is equivalent to $[lb+1,ub-1]$ for functions with integral coefficients. We again observe from the results that the intervals close to $bks$ do not contain as many feasible solutions. Moreover, the number of unique feasible solutions found while satisfying an interval goal as well as the deterministic goal are dependent on the solution landscape which is problem specific.

The diversity attributes for the set of solution vectors ($S$) output by the heuristic for a specific target were also analyzed. For this purpose, the hamming distance $h(x,y)$ between two solution vectors $x$ and $y$ given by the number of difference bits was computed for each pair of solution vectors. Thereafter, we calculated the mean hamming distance given by $\sum_{x,y\in S}^{x\neq y}h(x,y)/{|S|\choose 2}$. The denominator represents the number of unique solution pairs. The mean hamming distance for all output sets was between $1$ and $2$. The median hamming distance is $1$ or $2$ while the maximum hamming distance is $2$ or $3$ bits. Thus the solutions found were not diverse. In other words, alternate optima for a specific target near $bks$ are located in close proximity in the solution landscape. Based on this finding, it is more likely to find the next satisficing binary vector near the current one.

We present a novel model for two variants of goal-seeking QUBO involving either deterministic or uncertain targets. A number of alternate satisficing solutions were obtained using a greedy one-flip tabu search routine. Testing on benchmark instances clearly established the efficacy of our approach compared to state-of-the-art CP solvers which were unable to find solutions to any of the problems in the test set.

Regarding more detailed future research, our technique is easily extendable to multiple targets. For instance, achieving target levels of $t_{1}$ or $t_{2}$ would require parallelized instances of the routine searching for the two target levels. These threads could operate independently. Moreover, two different objective goals $t_{1}$ and $t_{2}$ with priorities $w_{1}$ and $w_{2}$ could be handled by changing the achievement function to $w_{1}(x^{\prime}Q_{1}x-t_{1})^{2}+w_{2}(x^{\prime}Q_{2}x-t_{2})^{2}$. In this way, we construct the efficient solution frontier by enumerating different choices of $w_{1}$ and $w_{2}$. Herein, the weights are determined by relative importance according to the decision-maker and is closely related to weighted or non-preemptive goal programming ([22]).