# Computational Comparison of Exact Solution Methods for 0-1 Quadratic Programs: Recommendations for Practitioners

Authors: Richard J. Forrester; Noah Hunt-Isaak
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

This paper is concerned with binary quadratic programs (BQPs), which are among the most well-studied classes of nonlinear integer optimization problems because of their wide variety of applications. While a number of different solution approaches have been proposed for tackling BQPs, practitioners need techniques that are both efficient and easy to implement. We revisit two of the most widely used linearization strategies for BQPs and examine the effectiveness of enhancements to these formulations that have been suggested in the literature. We perform a detailed large-scale computational study over five different classes of BQPs to compare these two linearizations with a more recent linear reformulation and direct submission of the nonlinear integer program to an optimization solver. The goal is to provide practitioners with guidance on how to best approach solving BQPs in an effective and easily implemented manner.


---

## Body

## 1. Introduction

Binary quadratic programs (BQPs) are one of the most well-studied classes of nonlinear integer optimization problems. These problems appear in a wide variety of applications (see [1, 2] for examples) and are known to be NP-hard. There are a number of different solution techniques that have been proposed for BQPs, including heuristics and exact solution methods. Given the difficulty of implementing and maintaining custom algorithms, the most commonly used exact solution method for BQPs involves linearizing the nonlinear problem and subsequently submitting the equivalent linear form to a standard mixed-integer linear programming (MILP) solver. Interestingly, the majority of commercial MILP solvers have been updated to handle the direct submission of a BQP, which provides an attractive solution alternative.In this paper, we revisit two of the most widely used linearization strategies for BQPs: the standard linearization [3] and Glover’s method [4]. There have been a number of proposed enhancements to these methods, but to the best of our knowledge, there has not been a comprehensive study to determine the optimal manner in which to apply these linearization techniques. In order to investigate these enhancements, we utilize five different classes of BQPs from the literature (the unconstrained BQP, the multidimensional quadratic knapsack problem, the k-item quadratic knapsack problem, the heaviest k-subgraph problem, and the quadratic semi-assignment problem) and three different MILP solvers (CPLEX, GUROBI, and XPRESS). We also make comparisons to the more recent linearization strategy of Sherali and Smith [5] and the direct submission of the BQP to the solver. The goal of this paper is to provide practitioners with suggestions for solving BQPs in an efficient and easily implemented manner. Note that our work is similar to that of [6, 7] who also perform computational studies of different linearization strategies for BQPs. However, our focus is not only on comparing different linearizations but also on how enhancements to the standard linearization and Glover’s method affect algorithmic performance.The general formulation of a BQP is given by(1)BQP:maximize∑i=1ncixi+∑i=1n∑j=1j≠inCijxixj:x∈X,xbinary,where cj and Cij are the linear and quadratic cost coefficients, respectively, and X is a polyhedral set representing the feasible region. For notational ease, we henceforth let the indices i and j run from 1 to n unless otherwise stated. Note that for each binary variable, xi2=xi, and thus, Cii can be incorporated into the linear term ci without loss of generality. Furthermore, note that the overall cost associated with the product xixj is Cij+Cji. Therefore, even if Cij and Cji are modified so that their sum is unchanged, the problem remains the same. The two most common quadratic coefficient representations in the literature are to either assume Cij=Cji for all i≠j (so that the quadratic coefficient matrix C is symmetric) or assume Cij=0 for all i≥j (so that the quadratic coefficient matrix C is upper triangular); both of which can be assumed without loss of generality. For our purposes, we will explicitly keep both terms Cijxixj and Cjixjxi. Note that there is no assumption that C is a negative-semidefinite matrix, and therefore, the continuous relaxation of problem BQP is not necessarily a convex optimization problem. However, with regards to the linearizations, the concavity of the objective function is irrelevant as it is rewritten into a linear form.


## 2. Linearizations

One of the most widely used techniques for optimizing a BQP is to utilize a linearization step that reformulates the nonlinear program into an equivalent linear form through the introduction of auxiliary variables and constraints. The linearized model can then be submitted to a standard MILP solver. While the size and continuous relaxation strength of a reformulation certainly play a role in the performance of a linearization when submitted to an MILP solver, it is not always possible to infer which formulation will have a better performance. Indeed, the use of preprocessing techniques, cutting planes, heuristics, and other enhancements makes it challenging to predict how two linearizations will compare computationally.In this section, we review two of the most well-known linearization strategies and consider modifications that can affect their computational performance. In addition, we describe the more recent linearization strategy of Sherali and Smith [5].


### 2.1. Standard Linearization

A standard method to linearize a BQP, due to Glover and Woolsey [3], is to replace each product xixj in the objective function with a continuous variable wij. In order to simplify the presentation, we first rewrite the objective function of BQP as
(2)∑i=1ncixi+∑i=1n−1∑j=i+1nCij′xixj,where Cij′=Cij+Cji for all i,j with i<j, which we can do without loss of generality. We now define the standard linearization of BQP below.
(3)STD:maximize∑j=1ncjxj+∑i=1n−1∑j=i+1nCij′wij.subject to
(4)wij≤xi∀i,j,i<j,(5)wij≤xj∀i,j,i<j,(6)wij≥xi+xj−1∀i,j,i<j,(7)wij≥0∀i,j,i<j,(8)x∈X,xbinary.Note that (4)–(7) enforce that wij=xixj for all binary x. While this linearization method is straightforward, it has the disadvantage that it requires the addition of nn−1/2 auxiliary variables and 4nn−1/2 auxiliary constraints.As noted in [8], we can reduce the number of auxiliary constraints by using the sign of Cij′. Let C−=i,j:Cij′<0: and C+=i,j:Cij′>0. Then, we can omit the constraints (4) and (5) bounding wij from above for i,j∈C−, and we can omit the constraints (6) and (7) bounding wij from below for i,j∈C+, as these will be implied at optimality. We refer to this reduced form as STD′ and note that the continuous relaxation of STD′ may potentially be weaker than that of STD′.This leads us to the first question that we would like to address:Question 1: When applying the standard linearization to a BQP, should you reduce the size of the formulation based on the sign of the quadratic objective coefficients? That is, which reformulation, STD or STD′, provides a better performance when submitted to a standard MILP solver?We will address this question in our computational study.


### 2.2. Glover’s Linearization

A more compact linearization method is due to Glover [4]. Given problem BQP, this method replaces each xj∑i=1,i≠jnCijxi) found in the objective function with a continuous variable zj and uses four linear restrictions to enforce that zj=xj∑i=1,i≠jnCijxi. Specifically, the Glover linearization is as follows.
(9)G1:maximize∑j=1ncjxj+∑j=1nzjsubject to
(10)zj≤Uj1xj∀j,(11)zj≥Lj1xj∀i,(12)zj≤∑i=1i≠jnCijxi−Lj01−xj∀j,(13)zj≥∑i=1i≠jnCijxi−Uj01−xj∀j,(14)x∈X,xbinaryFor eachj, Uj1 and Uj0 are upper bounds on ∑i=1,i≠jnCijxi, while Lj1xj and Lj0xj are lower bounds on ∑i=1,i≠jnCijxi. Problems BQP and G1 are equivalent in that, given any binary x, constraints (10)–(13) ensure that zj=xj∑i=1,i≠jnCijxi) for each j. Observe that this formulation only requires the addition of n (unrestricted) auxiliary continuous variables and 4n auxiliary constraints and is therefore considerably more compact than the standard linearization.As noted earlier, the two most common representations of the quadratic objective coefficientsCij in the literature are to either assume Cij=Cji for all i<j (so that the quadratic coefficient matrix C is symmetric) or assume Cij=0 for all i≥j (so that the quadratic coefficient matrix C is upper triangular); both of which can be assumed without loss of generality. While the continuous relaxation strength of STD and STD′ are not affected by the choice of objective function representation, the relaxation value of G1 is dependent on the manner in which the objective function of a BQP is expressed as shown in [9, 10]. This follows from the fact that the quadratic objective coefficients Cij do not appear in the auxiliary constraints of STD or STD′, whereas they do appear in those of G1. Moreover, if the quadratic coefficient matrix is upper triangular, then the variable z1 can be removed from G1 along with the associated constraints in (10)–(13). To see this, note that when C is upper triangular, (10)–(13) ensure that z1=x1∑i=1,i≠nnCi1xi=0 because Cij=0 for all i≥j. Thus, when the quadratic coefficient matrix is upper triangular, Glover’s formulation only requires the addition of n−1 auxiliary variables and 4n−1 auxiliary constraints. Furthermore, when C is upper triangular, the auxiliary constraints of G1 are less dense than when C is symmetric, which may provide a computational advantage when submitted to an MILP solver. This leads us to our next question, which we will address in our computational study.Question 2: When formulating Glover’s formulation, should you represent the quadratic objective coefficient matrix C in upper triangular or symmetric form?The bounds within G1 can be computed in a number of different ways. As originally suggested in [4], for each j, they can easily be computed as
(15)Lj1=Lj0=∑i=1,i≠jCij<0nCijandUj1=Uj0=∑i=1,i≠jCij>0nCij.As shown in [9], stronger bounds that take into consideration the feasible region can be computed as
(16)Ljp=min∑i=1i≠jnCijxi:x∈X,xj=pandUjp=max∑i=1i≠jnCijxi:x∈X,xj=pfor p∈0,1. These bounds could potentially be made even tighter by enforcing x binary so that
(17)Ljp=min∑i=1i≠jnCijxi:x∈X,xj=p,xbinaryandUjp=max∑i=1i≠jnCijxi:x∈X,xj=p,xbinary.While the continuous relaxation of G1 could be potentially tightened by strengthening the values of theUj1, Uj0, Lj1, and Lj0 bounds within the formulation, we need to take into consideration the amount of computational effort required to find the bounds. This leads us to our third question.Question 3: With regards to the overall computational effort to formulate and solve problem G1 to optimality using an MILP solver, should we compute the bounds Uj1, Uj0, Lj1, and Lj0 using (15), (16), or (17)?As noted in [11], we can reduce the size of G1 by removing the constraints (11) and (13) that bound zj from below because of the maximization objective and the fact that the zj terms do not appear elsewhere in the problem. Thus, problem G1 can be written more concisely as the modified Glover G2:
(18)G2:maximize∑j=1ncjxj+∑j=1nzjsubject to(19)5,7x∈X,xbinary.Interestingly, problem G2 only requires2n auxiliary constraints as opposed to the 4n required for G1, but G2 will have the same continuous relaxation strength as G1.Our next question to be addressed in the computational study is as follows.Question 4: How do formulations G1 and G2 compare when submitted to an MILP solver?It turns out that we can further reduce the number of structural constraints in G2 through the substitution of variablessj=Uj1xj−zj or sj=∑i≠jCijxi−Lj01−xj−zj for all j, which express the variables zj in terms of the slack variables to the inequalities (10) and (12), respectively (see [11]). Upon performing the first substitution, we obtain G2a:
(20)G2a:maximize∑j=1ncjxj+∑j=1nUj1xj−sjsubject to
(21)sj≥0∀j,sj≥Uj1−Lj0xj−∑i=1i≠jnCijxi+Lj0∀j,x∈X,xbinary.The second substitution yields G2b:(22)G2b:maximize∑j=1ncjxj+∑j=1n∑i=1i≠jnCijxi−Lj01−xj−sjsubject to
(23)sj≥∑i=1i≠jnCijxi−Uj1−Lj0xj−Lj0∀j,sj≥0∀j,x∈X,xbinary.Note that both G2a and G2b only requiren (nonnegative) auxiliary variables and n auxiliary constraints, yet they both have the same continuous relaxation strength as G1 and G2. This leads us to our next question:Question 5: Is it advantageous to perform the substitution of variables to reduce problem G2 to G2a or G2b when submitting the model to an MILP solver? That is, how do formulations G2, G2a, and G2b compare when submitted to an MILP solver?


### 2.3. Sherali-Smith Linear Formulation

A more recent linearization was introduced by Sherali and Smith [5]. This method converts BQP to the linear form below.
(24)SS:maximize∑i=1nsi+∑i=1nci+Lixisubject to(25)yi=∑j=1nCijxj−si−Li∀i,(26)yi≤Ui−Li1−xi∀i,(27)si≤Ui−Lixi∀i,(28)yi≥0∀i,(29)si≥0∀i,(30)x∈X,xbinary.For eachi, Ui and Li are upper and lower bounds, respectively, on ∑j=1,j≠inCijxj, and can be computed in a similar fashion as in (15), (16), or (17). However, based upon preliminary computational results, we will construct SS using the bounds computed as in (16). Note that the authors of [5] actually introduced three formulations, called BP, BP¯, and BP-strong, which all consider instances of BQP that include quadratic constraints. However, each of these three formulations is equivalent to SS for BQP. After performing the substitutions suggested by constraints (25), SS increases the size of the problem by adding an additional n nonnegative auxiliary variables and 3n auxiliary constraints (26)–(29).


## 2.1. Standard Linearization

A standard method to linearize a BQP, due to Glover and Woolsey [3], is to replace each product xixj in the objective function with a continuous variable wij. In order to simplify the presentation, we first rewrite the objective function of BQP as
(2)∑i=1ncixi+∑i=1n−1∑j=i+1nCij′xixj,where Cij′=Cij+Cji for all i,j with i<j, which we can do without loss of generality. We now define the standard linearization of BQP below.
(3)STD:maximize∑j=1ncjxj+∑i=1n−1∑j=i+1nCij′wij.subject to
(4)wij≤xi∀i,j,i<j,(5)wij≤xj∀i,j,i<j,(6)wij≥xi+xj−1∀i,j,i<j,(7)wij≥0∀i,j,i<j,(8)x∈X,xbinary.Note that (4)–(7) enforce that wij=xixj for all binary x. While this linearization method is straightforward, it has the disadvantage that it requires the addition of nn−1/2 auxiliary variables and 4nn−1/2 auxiliary constraints.As noted in [8], we can reduce the number of auxiliary constraints by using the sign of Cij′. Let C−=i,j:Cij′<0: and C+=i,j:Cij′>0. Then, we can omit the constraints (4) and (5) bounding wij from above for i,j∈C−, and we can omit the constraints (6) and (7) bounding wij from below for i,j∈C+, as these will be implied at optimality. We refer to this reduced form as STD′ and note that the continuous relaxation of STD′ may potentially be weaker than that of STD′.This leads us to the first question that we would like to address:Question 1: When applying the standard linearization to a BQP, should you reduce the size of the formulation based on the sign of the quadratic objective coefficients? That is, which reformulation, STD or STD′, provides a better performance when submitted to a standard MILP solver?We will address this question in our computational study.


## 2.2. Glover’s Linearization

A more compact linearization method is due to Glover [4]. Given problem BQP, this method replaces each xj∑i=1,i≠jnCijxi) found in the objective function with a continuous variable zj and uses four linear restrictions to enforce that zj=xj∑i=1,i≠jnCijxi. Specifically, the Glover linearization is as follows.
(9)G1:maximize∑j=1ncjxj+∑j=1nzjsubject to
(10)zj≤Uj1xj∀j,(11)zj≥Lj1xj∀i,(12)zj≤∑i=1i≠jnCijxi−Lj01−xj∀j,(13)zj≥∑i=1i≠jnCijxi−Uj01−xj∀j,(14)x∈X,xbinaryFor eachj, Uj1 and Uj0 are upper bounds on ∑i=1,i≠jnCijxi, while Lj1xj and Lj0xj are lower bounds on ∑i=1,i≠jnCijxi. Problems BQP and G1 are equivalent in that, given any binary x, constraints (10)–(13) ensure that zj=xj∑i=1,i≠jnCijxi) for each j. Observe that this formulation only requires the addition of n (unrestricted) auxiliary continuous variables and 4n auxiliary constraints and is therefore considerably more compact than the standard linearization.As noted earlier, the two most common representations of the quadratic objective coefficientsCij in the literature are to either assume Cij=Cji for all i<j (so that the quadratic coefficient matrix C is symmetric) or assume Cij=0 for all i≥j (so that the quadratic coefficient matrix C is upper triangular); both of which can be assumed without loss of generality. While the continuous relaxation strength of STD and STD′ are not affected by the choice of objective function representation, the relaxation value of G1 is dependent on the manner in which the objective function of a BQP is expressed as shown in [9, 10]. This follows from the fact that the quadratic objective coefficients Cij do not appear in the auxiliary constraints of STD or STD′, whereas they do appear in those of G1. Moreover, if the quadratic coefficient matrix is upper triangular, then the variable z1 can be removed from G1 along with the associated constraints in (10)–(13). To see this, note that when C is upper triangular, (10)–(13) ensure that z1=x1∑i=1,i≠nnCi1xi=0 because Cij=0 for all i≥j. Thus, when the quadratic coefficient matrix is upper triangular, Glover’s formulation only requires the addition of n−1 auxiliary variables and 4n−1 auxiliary constraints. Furthermore, when C is upper triangular, the auxiliary constraints of G1 are less dense than when C is symmetric, which may provide a computational advantage when submitted to an MILP solver. This leads us to our next question, which we will address in our computational study.Question 2: When formulating Glover’s formulation, should you represent the quadratic objective coefficient matrix C in upper triangular or symmetric form?The bounds within G1 can be computed in a number of different ways. As originally suggested in [4], for each j, they can easily be computed as
(15)Lj1=Lj0=∑i=1,i≠jCij<0nCijandUj1=Uj0=∑i=1,i≠jCij>0nCij.As shown in [9], stronger bounds that take into consideration the feasible region can be computed as
(16)Ljp=min∑i=1i≠jnCijxi:x∈X,xj=pandUjp=max∑i=1i≠jnCijxi:x∈X,xj=pfor p∈0,1. These bounds could potentially be made even tighter by enforcing x binary so that
(17)Ljp=min∑i=1i≠jnCijxi:x∈X,xj=p,xbinaryandUjp=max∑i=1i≠jnCijxi:x∈X,xj=p,xbinary.While the continuous relaxation of G1 could be potentially tightened by strengthening the values of theUj1, Uj0, Lj1, and Lj0 bounds within the formulation, we need to take into consideration the amount of computational effort required to find the bounds. This leads us to our third question.Question 3: With regards to the overall computational effort to formulate and solve problem G1 to optimality using an MILP solver, should we compute the bounds Uj1, Uj0, Lj1, and Lj0 using (15), (16), or (17)?As noted in [11], we can reduce the size of G1 by removing the constraints (11) and (13) that bound zj from below because of the maximization objective and the fact that the zj terms do not appear elsewhere in the problem. Thus, problem G1 can be written more concisely as the modified Glover G2:
(18)G2:maximize∑j=1ncjxj+∑j=1nzjsubject to(19)5,7x∈X,xbinary.Interestingly, problem G2 only requires2n auxiliary constraints as opposed to the 4n required for G1, but G2 will have the same continuous relaxation strength as G1.Our next question to be addressed in the computational study is as follows.Question 4: How do formulations G1 and G2 compare when submitted to an MILP solver?It turns out that we can further reduce the number of structural constraints in G2 through the substitution of variablessj=Uj1xj−zj or sj=∑i≠jCijxi−Lj01−xj−zj for all j, which express the variables zj in terms of the slack variables to the inequalities (10) and (12), respectively (see [11]). Upon performing the first substitution, we obtain G2a:
(20)G2a:maximize∑j=1ncjxj+∑j=1nUj1xj−sjsubject to
(21)sj≥0∀j,sj≥Uj1−Lj0xj−∑i=1i≠jnCijxi+Lj0∀j,x∈X,xbinary.The second substitution yields G2b:(22)G2b:maximize∑j=1ncjxj+∑j=1n∑i=1i≠jnCijxi−Lj01−xj−sjsubject to
(23)sj≥∑i=1i≠jnCijxi−Uj1−Lj0xj−Lj0∀j,sj≥0∀j,x∈X,xbinary.Note that both G2a and G2b only requiren (nonnegative) auxiliary variables and n auxiliary constraints, yet they both have the same continuous relaxation strength as G1 and G2. This leads us to our next question:Question 5: Is it advantageous to perform the substitution of variables to reduce problem G2 to G2a or G2b when submitting the model to an MILP solver? That is, how do formulations G2, G2a, and G2b compare when submitted to an MILP solver?


## 2.3. Sherali-Smith Linear Formulation

A more recent linearization was introduced by Sherali and Smith [5]. This method converts BQP to the linear form below.
(24)SS:maximize∑i=1nsi+∑i=1nci+Lixisubject to(25)yi=∑j=1nCijxj−si−Li∀i,(26)yi≤Ui−Li1−xi∀i,(27)si≤Ui−Lixi∀i,(28)yi≥0∀i,(29)si≥0∀i,(30)x∈X,xbinary.For eachi, Ui and Li are upper and lower bounds, respectively, on ∑j=1,j≠inCijxj, and can be computed in a similar fashion as in (15), (16), or (17). However, based upon preliminary computational results, we will construct SS using the bounds computed as in (16). Note that the authors of [5] actually introduced three formulations, called BP, BP¯, and BP-strong, which all consider instances of BQP that include quadratic constraints. However, each of these three formulations is equivalent to SS for BQP. After performing the substitutions suggested by constraints (25), SS increases the size of the problem by adding an additional n nonnegative auxiliary variables and 3n auxiliary constraints (26)–(29).


## 3. Binary Quadratic Program Classifications

In this section, we introduce the five different classes of BQPs that we consider in our computational study and discuss the manner in which our problem instances are generated.


### 3.1. Unconstrained 0-1 Quadratic Program: Boolean Least Squares Problem

The first family of 0-1 quadratic problems that we investigate is the unconstrained BQP (UBQP), so thatX is defined as
(31)X≡x∈ℝn:0≤xi≤1fori=1,⋯,n.While simplistic, the UBQP is notable for its ability to represent a wide range of applications (see [2] for a recent survey). For our experiments, we decided to focus on the Boolean least squares problem (BLSP), which is a basic problem in digital communication where the objective is to identify a binary signal x from a collection of noisy measurements. Traditionally, this problem is modeled as
(32)minimizeDx−d2=xTDTDx−2dTDx+dTdsubject to
(33)X≡x∈ℝn:0≤xi≤1fori=1,⋯,n,where D∈ℝm×n and d∈ℝm are given. By ignoring the constant term dTd, we can rewrite this formulation in our notation and in maximization form by setting Q=DTD, q=−2dTD, and subsequently defining Cij=−Qij for i≠j, Cij=0 for i=j, and ci=2qi+Qii for all i.We randomly generated instances of BLSP as in [1]. Specifically, we set m=n and generated a random D matrix with elements from the standard normal distribution N0,1 and a binary vector y∈ℝn with elements from the uniform distribution U0,1. The vector d is then constructed as d=Dy+ε where ε is a random noise term with elements from N0,1. We generated 10 instances for each value of n, which we varied from 30 to 70 in increments of 10 to provide a range of progressively more challenging problems.


### 3.2. Quadratic Multidimensional Knapsack Problem

The next class of problems that we consider is the quadratic knapsack problem with multiple constraints, also known as the quadratic multidimensional knapsack problem (QMKP) (see [10, 12]). These problems have the following form.
(34)QMKP:maximize∑j=1ncjxj+∑j=1n∑i=1i≠jnCijxixjsubject to
(35)∑i=1nakixi≤bkfork=1,⋯m,xbinary.Within QMKP, thecj, aki, Cij, and bk coefficients are typically nonnegative scalars. We assume for every k that aki≤bk for each i since otherwise variables can be fixed to 0 and that bk<∑i=1naki for the kth constraint to be restrictive.Problem QMKP is a generalization of both the multidimensional knapsack problem and the quadratic knapsack problem (QKP). The multidimensional knapsack problem is that case of QMKP whereinCij=0 for all i,j, so that the problem reduces to linear form. The QKP retains the quadratic objective terms but has only a single knapsack constraint defining Xm=1. This latter problem is one of the most extensively studied nonlinear discrete optimization programs in the literature (see [13] for an excellent survey).We randomly generated two different test beds; both of which contained instances withk=1,5 and 10 knapsack constraints. The coefficients aki are integers taken from a uniform distribution over the interval [1,50], and bk is an integer from a uniform distribution between 50 and ∑i=1naki.For the first test bed, the objective coefficients were all nonnegative. Specifically, the (nonzero) objective coefficientscj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval [1,100]. To assess the effects of the density of nonzero cj and Cij on CPU time, each of these coefficients is nonzero with some predetermined probability Δ. We considered instances with probabilities (densities) Δ=25%,50%,75%,and100%.For each value of Δ, we generated ten instances for each value of n, which varied from 80 to 110 in increments of 10 when k=1 and varied from 30 to 70 in increments of 10 when k=5 or 10.For the second test bed, the objective coefficients were mixed in sign. Specifically, the (nonzero) objective coefficientscj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval −100,100. As with the first test bed, the density of the objective coefficients was varied from Δ=25–100% in increments of 25%. For each value of Δ, we generated 10 instances for each value of n, which varied from 20 to 50 in increments of 10 for k=1,5, and 10. We utilized smaller sizes as objective coefficients that are mixed in sign are significantly more difficult to solve.


### 3.3. The Heaviestk-Subgraph Problem

The heaviestk-subgraph problem (HSP) is concerned with determining a block of k nodes of a weighted graph such that the total edge weight within the subgraph induced by the block is maximized [14]. The HSP can be formulated as
(36)HSP:maximize∑j=1n∑i=1i≠jnCijxixjsubject to
(37)∑i=1nxi=k,xbinary.The HSP is also known under the name of the cardinality constrained quadratic binary program, the densestk-subgraph problem, the p-dispersion-sum problem, and the k-cluster problem. For each number of nodes, n=10,20,30,and40, we randomly generated 10 instances with three different graph densities, Δ=25%,50%,and75%. Specifically, for each density Δ and any pair of indexes i,j such that i<j, we assigned Cij=1 with probability Δ and Cij=0 otherwise. We set k=0.5n for all of our tests.


### 3.4.k-Item Quadratic Knapsack Problem

Thek-item quadratic knapsack problem (kQKP), introduced in [15, 16], consists of maximizing a quadratic function subject to a linear knapsack constraint with an additional equality cardinality constraint:
(38)kQKP:maximize∑j=1ncjxj+∑j=1n∑i=1i≠jnCijxixjsubject to
(39)∑i=1naixi≤b,(40)∑i=1nxi=k,(41)xbinary.Here, we assume thatai≤b for each i since otherwise variables can be fixed to 0 and that b<∑i=1nai for the constraint to be restrictive. Let us denote by kmax the largest number of items which can be filled in the knapsack, that is, the largest number of the smallest ai whose sum does not exceed b. We can assume that k∈2,⋯,kmax.Note thatkQKP includes two classical subproblems, the HSP by dropping constraint (39) and the QKP by dropping constraint (40). We randomly generated two different tests beds; both of which contained instances with randomly chosen values of k from the discrete uniform distribution over the range 2,n/4. For the first test bed, the objective coefficients were all nonnegative, where the (nonzero) objective coefficients cj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval [1,100]. The density of the objective coefficients was varied from Δ=25–100% in increments of 25%. For each value of Δ, we generated ten instances with n=50,60,and70 variables. For the second test bed, the objective coefficients were mixed in sign, where the objective coefficients cj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval −100,100. As with the first test bed, the density of the objective coefficients was varied from Δ=25–100% in increments of 25%. For each value of Δ, we generated 10 instances for each value of n, where n varied from 50 to 90 in increments of 10.


### 3.5. Quadratic Semiassignment Problem: Task Allocation Problem

Here, we consider a specific instance of the quadratic semiassignment problem (QSAP) called the task allocation problem [17]. A set of tasks T1,⋯,Tm is to be allocated to a set of processors P1,⋯,Pn. The execution costs for assigning task Ti to processor Pk is denoted by eik and the communication costs between Ti on Pk and Tj on Pl is denoted by Cikjl. The constraints limit every task to exactly one processor. The specific formulation is as follows.
(42)QSAP:minimize∑i=1m∑k=1neikxik+∑i=1m−1∑j=i+1m∑k=1n∑l=1nCikjlxikxjlsubject to
(43)∑k=1nxik=1fori=1,⋯m,xbinary.For our tests, we considered two values ofn (4 and 5) and three values of m (10, 12, 15). For each of the six m,n pairs, we generated 10 instances where the coefficients eik and Cikjl are integers from a uniform distribution over the interval [−50,50] as in [1]. Note that in our tests, we changed the sense of the objective function of QSAP to maximization.


## 3.1. Unconstrained 0-1 Quadratic Program: Boolean Least Squares Problem

The first family of 0-1 quadratic problems that we investigate is the unconstrained BQP (UBQP), so thatX is defined as
(31)X≡x∈ℝn:0≤xi≤1fori=1,⋯,n.While simplistic, the UBQP is notable for its ability to represent a wide range of applications (see [2] for a recent survey). For our experiments, we decided to focus on the Boolean least squares problem (BLSP), which is a basic problem in digital communication where the objective is to identify a binary signal x from a collection of noisy measurements. Traditionally, this problem is modeled as
(32)minimizeDx−d2=xTDTDx−2dTDx+dTdsubject to
(33)X≡x∈ℝn:0≤xi≤1fori=1,⋯,n,where D∈ℝm×n and d∈ℝm are given. By ignoring the constant term dTd, we can rewrite this formulation in our notation and in maximization form by setting Q=DTD, q=−2dTD, and subsequently defining Cij=−Qij for i≠j, Cij=0 for i=j, and ci=2qi+Qii for all i.We randomly generated instances of BLSP as in [1]. Specifically, we set m=n and generated a random D matrix with elements from the standard normal distribution N0,1 and a binary vector y∈ℝn with elements from the uniform distribution U0,1. The vector d is then constructed as d=Dy+ε where ε is a random noise term with elements from N0,1. We generated 10 instances for each value of n, which we varied from 30 to 70 in increments of 10 to provide a range of progressively more challenging problems.


## 3.2. Quadratic Multidimensional Knapsack Problem

The next class of problems that we consider is the quadratic knapsack problem with multiple constraints, also known as the quadratic multidimensional knapsack problem (QMKP) (see [10, 12]). These problems have the following form.
(34)QMKP:maximize∑j=1ncjxj+∑j=1n∑i=1i≠jnCijxixjsubject to
(35)∑i=1nakixi≤bkfork=1,⋯m,xbinary.Within QMKP, thecj, aki, Cij, and bk coefficients are typically nonnegative scalars. We assume for every k that aki≤bk for each i since otherwise variables can be fixed to 0 and that bk<∑i=1naki for the kth constraint to be restrictive.Problem QMKP is a generalization of both the multidimensional knapsack problem and the quadratic knapsack problem (QKP). The multidimensional knapsack problem is that case of QMKP whereinCij=0 for all i,j, so that the problem reduces to linear form. The QKP retains the quadratic objective terms but has only a single knapsack constraint defining Xm=1. This latter problem is one of the most extensively studied nonlinear discrete optimization programs in the literature (see [13] for an excellent survey).We randomly generated two different test beds; both of which contained instances withk=1,5 and 10 knapsack constraints. The coefficients aki are integers taken from a uniform distribution over the interval [1,50], and bk is an integer from a uniform distribution between 50 and ∑i=1naki.For the first test bed, the objective coefficients were all nonnegative. Specifically, the (nonzero) objective coefficientscj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval [1,100]. To assess the effects of the density of nonzero cj and Cij on CPU time, each of these coefficients is nonzero with some predetermined probability Δ. We considered instances with probabilities (densities) Δ=25%,50%,75%,and100%.For each value of Δ, we generated ten instances for each value of n, which varied from 80 to 110 in increments of 10 when k=1 and varied from 30 to 70 in increments of 10 when k=5 or 10.For the second test bed, the objective coefficients were mixed in sign. Specifically, the (nonzero) objective coefficientscj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval −100,100. As with the first test bed, the density of the objective coefficients was varied from Δ=25–100% in increments of 25%. For each value of Δ, we generated 10 instances for each value of n, which varied from 20 to 50 in increments of 10 for k=1,5, and 10. We utilized smaller sizes as objective coefficients that are mixed in sign are significantly more difficult to solve.


## 3.3. The Heaviestk-Subgraph Problem

The heaviestk-subgraph problem (HSP) is concerned with determining a block of k nodes of a weighted graph such that the total edge weight within the subgraph induced by the block is maximized [14]. The HSP can be formulated as
(36)HSP:maximize∑j=1n∑i=1i≠jnCijxixjsubject to
(37)∑i=1nxi=k,xbinary.The HSP is also known under the name of the cardinality constrained quadratic binary program, the densestk-subgraph problem, the p-dispersion-sum problem, and the k-cluster problem. For each number of nodes, n=10,20,30,and40, we randomly generated 10 instances with three different graph densities, Δ=25%,50%,and75%. Specifically, for each density Δ and any pair of indexes i,j such that i<j, we assigned Cij=1 with probability Δ and Cij=0 otherwise. We set k=0.5n for all of our tests.


## 3.4.k-Item Quadratic Knapsack Problem

Thek-item quadratic knapsack problem (kQKP), introduced in [15, 16], consists of maximizing a quadratic function subject to a linear knapsack constraint with an additional equality cardinality constraint:
(38)kQKP:maximize∑j=1ncjxj+∑j=1n∑i=1i≠jnCijxixjsubject to
(39)∑i=1naixi≤b,(40)∑i=1nxi=k,(41)xbinary.Here, we assume thatai≤b for each i since otherwise variables can be fixed to 0 and that b<∑i=1nai for the constraint to be restrictive. Let us denote by kmax the largest number of items which can be filled in the knapsack, that is, the largest number of the smallest ai whose sum does not exceed b. We can assume that k∈2,⋯,kmax.Note thatkQKP includes two classical subproblems, the HSP by dropping constraint (39) and the QKP by dropping constraint (40). We randomly generated two different tests beds; both of which contained instances with randomly chosen values of k from the discrete uniform distribution over the range 2,n/4. For the first test bed, the objective coefficients were all nonnegative, where the (nonzero) objective coefficients cj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval [1,100]. The density of the objective coefficients was varied from Δ=25–100% in increments of 25%. For each value of Δ, we generated ten instances with n=50,60,and70 variables. For the second test bed, the objective coefficients were mixed in sign, where the objective coefficients cj for all j and Cij=Cji for i,j with i<j are integers from a uniform distribution over the interval −100,100. As with the first test bed, the density of the objective coefficients was varied from Δ=25–100% in increments of 25%. For each value of Δ, we generated 10 instances for each value of n, where n varied from 50 to 90 in increments of 10.


## 3.5. Quadratic Semiassignment Problem: Task Allocation Problem

Here, we consider a specific instance of the quadratic semiassignment problem (QSAP) called the task allocation problem [17]. A set of tasks T1,⋯,Tm is to be allocated to a set of processors P1,⋯,Pn. The execution costs for assigning task Ti to processor Pk is denoted by eik and the communication costs between Ti on Pk and Tj on Pl is denoted by Cikjl. The constraints limit every task to exactly one processor. The specific formulation is as follows.
(42)QSAP:minimize∑i=1m∑k=1neikxik+∑i=1m−1∑j=i+1m∑k=1n∑l=1nCikjlxikxjlsubject to
(43)∑k=1nxik=1fori=1,⋯m,xbinary.For our tests, we considered two values ofn (4 and 5) and three values of m (10, 12, 15). For each of the six m,n pairs, we generated 10 instances where the coefficients eik and Cikjl are integers from a uniform distribution over the interval [−50,50] as in [1]. Note that in our tests, we changed the sense of the objective function of QSAP to maximization.


## 4. Computational Study

In this section, we attempt to answer the questions posed in the previous sections. In addition, we compare the standard linearization and Glover’s method with both the linearization strategy of Sherali and Smith and direct submission of the BQP to the solver. All tests were implemented in Python and executed on a Dell Precision Workstation equipped with dual 2.4 GHz Xeon processors and 64 GB of RAM running 64bit Windows 10 Professional. In order to better understand the performance of a linear formulation independent of the solver, we utilized three different MILP optimization packages, CPLEX 12.8.0, GUROBI 8.0.0, and XPRESS 8.5.0, which were the most current versions of the software available when we performed the computational tests. We utilized the default settings on all the MILP solvers with a time limit of 3,600 seconds (1 hr) per instance. Note that for the sake of brevity, we only present a subset of the data generated during our study. However, more complete data is available upon request to the corresponding author.We present our results usingperformance profiles, which is a graphical comparison methodology introduced by Dolan and Moré [18] to benchmark optimization software. Rather than analyzing the data using average solve times that tend to be skewed by a small number of difficult problem instances, performance profiles provide the cumulative distribution function of the ratio of the solve time of a particular formulation versus the best time of all formulations as the performance metric. This effectively eliminates the influence of any outliers and provides a graphical presentation of the expected performance difference among the different formulations.In order to describe our specific implementation of the performance profiles, suppose that we have a setF of nf formulations and a set P of np problem instances. For each problem p and formulation f, we define tp,f=computingtimerequiredtosolveproblempusingformulationf.We compare the performance of formulationf on problem p with the best performance by any formulation on this problem using the performance ratio
(44)rp,f=tp,fmintp,f:f∈F.Then, thecumulative performance profile is defined as
(45)ρfτ=1npsizep∈P:rp,f≤τ,so that ρfτ is the probability for formulation f∈F that a performance ratio rp,f is within a factor τ of the best possible ratio. Thus, for example, ρf1 is the probability of formulation f to be the best formulation to solve any given problem, while ρf3 is the probability of formulation f to be the best within a time factor of 3 of the best formulation. Note that the performance profiles also provide information about the formulations that fail to reach optimality within the time limit. In particular, those profiles that do not meet ρfτ=1 indicate that there is a probability that they will not solve a subset of problems to optimality within the time limit.


### 4.1. Results for the Standard Linearization

We begin with our first question raised: When implementing the standard linearization, should you utilize the full formulation or use the formulation that is based on the sign of the quadratic objective coefficients? That is, how does STD compare with STD′ when submitted to an MILP solver?Figure1 shows the performance obtained for STD and STD′ for each class of BQP, where we aggregated all of the sizes and densities for each problem class, along with the three different MILP solvers. Note that since STD and STD′ vary depending on the sign of the objective coefficients, we provide two performance profiles for both QMKP and kQKP, one for instances with positive quadratic coefficients and one for instances with quadratic coefficients of mixed sign. For the reader unfamiliar with performance profiles, let us examine the results for the BLSP presented in Figure 1(a) in detail. Here, nf=2 as we are comparing the two formulations STD and STD′, while np=5 sizes n=30,40,50,60,70×10instances×3MILPsolvers=150. Each of these 150 instances is then solved using STD and STD′. For each instance, the minimum computational time of the two formulations is determined, and then the computational time of each formulation is divided by the minimum time, providing the performance ratio. Using the calculation of the performance ratios for all instances, the cumulative profile is then calculated for a set of fixed ratios of time factors. Recall that ρf1 is the probability formulation f will have the fastest solve time for any given problem. Thus, the probability that STD′ will have the fastest solve time for BLSP is approximately 86%, whereas the probability that STD will be the fastest formulation is approximately 14%.Figure 1
Performance profiles based on relative CPU time (s) for Question 1. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The performance profiles clearly indicate that STD′ is more likely to have the fastest solve time for all of the BQPs under consideration. While not presented in Figure 1, when aggregated across all problem classes, the probability that STD′ is the fastest formulation was 92%. We also note that a more detailed analysis showed that STD′ was the superior formulation regardless of objective coefficient density, problem size, or MILP solver—we omit the details for brevity.While performance profiles provide information on how likely a particular formulation will have a faster solve time, they do not indicate the actual differences in solve time. We report the average CPU times to solve the formulations to proven optimality for the different BQP classes in Figure2. These average CPU times were computed after removing all instances where at least one of the formulations did not solve within the time limit. That is, the averages only include instances in which both formulations were solved to optimality within the 1-hour time limit. Note that the performance profiles indicate that the number of instances that were not solved within the time limit was fairly small. To see this, recall that those profiles that do not meet ρfτ=1 indicate that there is a probability that they will not solve a subset of problems to optimality within the time limit. For example, note that approximately 5% of the instances of HSP formulated with STD did not solve within the 1-hour limit. The data in Figure 2 indicates that STD′ is significantly faster than STD on average for all BQP classes considered. Indeed, for 5 of the 7 classifications, the average CPU time of STD′ is less than half of that of STD.Figure 2
Average CPU times (s) for Question 1.Thus, our results clearly indicate that when utilizing the standard linearization, it is best to use the sign-based version STD′ regardless of BQP class, objective function density, problem size, or MILP solver. It is interesting to note that STD has a continuous relaxation that is at least as tight as STD′ but comes at the expense of a larger problem size. Our tests clearly indicate that the stronger continuous relaxation of STD does not provide an overall computational advantage.


### 4.2. Results for Glover’s Formulation

In this section, we investigate the four questions posed regarding Glover’s formulation. While ideally we would investigate these questions simultaneously, doing so would require a prohibitive number of test cases. Therefore, we investigate the questions in an incremental manner to make the task more tractable.We begin by focusing on answering the second and fourth questions raised. That is, when formulating Glover’s model, should we record the quadratic objective coefficient matrix in a symmetric or upper triangular form, and should we reduce the size of G1 by removing constraints (11) and (13) to obtain the modified Glover formulation G2? In order to compare the use of a symmetric quadratic coefficient matrix C with that of one in upper triangular form when formulating G1 and G2, we replaced Cij in the objective function of BQP with Cij′ defined as follows. To obtain an upper triangular objective coefficient matrix, we set Cij′=Cij+Cji for all i,j with i<j and Cij′=0 for i≥j. To obtain a symmetric objective coefficient matrix, we set Cij′=Cji′=1/2Cij+Cji) for all i,j with i<j and Cii′=0 for all i.Figure3 shows the performance of the four different formulations under consideration: G1 and G2 formulated with both symmetric and upper triangular objective coefficient matrices. We once again aggregate across all objective coefficient densities, sizes, and MILP solvers. Note that when formulating G1 and G2, we computed the bounds within the formulations as in (16) since preliminary computational experience showed this method provides the best performance. However, we will formally compare the bounds of (15), (16), and (17) in another set of tests. We also note that when computing solution times, we included the time needed to find the bounds within G1 and G2.Figure 3
Performance profiles based on relative CPU time (s) for Questions 2 and 4. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The results in Figure3 provide a fairly definitive answer to question 4—Problem G2 is superior to G1 when compared using the same objective coefficient representation. Indeed, G2 with objective coefficients in upper triangular form had the highest probability of being the fastest formulation for all classes of BQPs, while G2 with objective coefficients in symmetric form had the second highest probability of being the fastest formulation for all classes except BLSP and QSAP. With regards to Question 2, note that both G1 and G2 formulated with objective coefficients in upper triangular form were more likely to be faster than their respective versions formulated with symmetric objective coefficients, with the exception of instances of kQKP with positive objective coefficients. Interestingly, we note that a more detailed analysis (not presented here) showed that G2 with symmetric objective coefficients was also superior for higher density (Δ=75%and100%) instances of QMKP with positive coefficients. Note that the performance profiles also indicate that some of the formulations struggled with being able to solve all of the instances within the time limit. For example, G1 and G2 formulated with symmetric objective coefficients were unable to solve approximately 50% of the instances of QSAP, while none of the formulations were able to solve all of the instances of HSP.We report the average CPU times to solve the formulations to proven optimality for the different BQP classes in Figure4. As before, these average CPU times were computed after removing all instances where at least one of the formulations did not solve within the time limit.Figure 4
Average CPU times (s) for Questions 2 and 4.Figure4 provides some additional information not captured in the performance profiles. First, note that while the performance profiles in Figure 3 indicate that G2 was superior to G1, Figure 4 indicates that the actual difference in average solution times between G1 and G2 is fairly negligible when they are formulated with the same objective representation. However, there are stark differences between the symmetric and upper triangular coefficient representations for both G1 and G2. In particular, instances of G1 and G2 formulated with symmetric objective coefficients typically take considerably longer to solve on average than their respective models formulated with upper triangular coefficients. There are two exceptions—G2 with symmetric objective coefficients is superior for HSP and instances of kQKP with positive coefficients.While the results are not completely definitive, Figures3 and 4 indicate that in general, constraints (11) and (13) should be dropped from G1 to obtain G2 and the quadratic objective coefficients C should be represented in upper triangular form. We reiterate that while the average CPU time differences between G1 and G2 are fairly small, those differences are rather large for the two different objective function representations.We now attempt to answer Question 5, which is concerned with whether it is advantageous to perform the substitution of variables that reduces G2 to either G2a or G2b. Figure5 shows the performance profiles of G2, G2a, and G3b where we aggregate across all objective coefficient densities, sizes, and MILP solvers. As with our previous test cases, we computed the bounds within the formulations as in (16).Figure 5
Performance profiles based on relative CPU time (s) for Question 5. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Note that, typically, G2 was the slowest formulation while G2a and G2b had fairly comparable likelihoods of being the fastest. Although, a closer examination reveals that G2a was slightly more likely to be the fastest formulation for all problem classes except instances of QMKP with positive objective coefficients. We report the average CPU times in Figure6, which were computed as before.Figure 6
Average CPU times (s) for Question 5.While the average CPU times for G2, G2a, and G2b are all fairly similar for each problem classification, G2a tended to be the fastest formulation. Therefore, based on our tests, we recommend that G2 should be reduced in size via the substitution of variablessj=Uj1xj−zj to obtain G2a. While the results were certainly not definitive, they do suggest that G2a has the highest probability of being the fastest formulation.We now address Question 3 regarding how the bounds within Glover’s formulation should be computed. Within our tests, we refer to the bounds computed as in (15) as weak, those computed as in (16) as tight, and those as in (17) as tightest. Based upon our prior results, our tests were all performed using G2a formulated with the three different methods to compute the bounds. We again reiterate that our solution times include not only the time to solve G2a to integer optimality but also the time required to compute the bounds within the formulation. Our results are presented in Figure 7 where we once again aggregate across all objective coefficient densities, sizes, and MILP solvers. Note that we do not include the BLSP instances because for these problems, the bounds (15), (16), and (17) are all the same.Figure 7
Performance profiles based on relative CPU time (s) for Question 3. (a) QMKP with positive coefficients, (b) QMKP with mixed-sign coefficients, (c)kQKP with positive coefficients, (d) kQKP with mixed coefficients, (e) HSP, (f) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)Interestingly, note that G2a formulated with the weak bounds of (15) tended to be the fastest formulation. We found this result surprising as preliminary experience showed that the tight bounds of (16) were superior. A more detailed analysis revealed that for smaller-sized instances, the extra time needed to find the tighter bounds of (16) was not warranted. However, once the instances start to become more challenging, the additional strength afforded by the tight bounds (16) was advantageous even though the time to find the bounds increases. In order to better understand how these bounds affect CPU time, let us examine the average CPU times presented in Figure 8.Figure 8
Average CPU times (s) for Question 3.Note that in terms of average CPU times, G2a formulated with the tight bounds (16) was the best formulation for all problem classifications except instances of kQKP with positive objective coefficients. Moreover, not only did the weak bounds (15) perform poorly in terms of average CPU times, they typically performed significantly worse than the tighter bounds. This stems from the fact that the continuous relaxation strength of G2a formulated with (15) for anything other than smaller-sized instances is too weak for the linearization to be effective. Thus, our general recommendation is to formulate Glover’s model with the tight bounds of (16), especially for larger-sized instances.We summarize the answers and recommendations from Questions 1 through 5 in Table1.Table 1
Summary of answers to Questions 1 through 5.
Question
Answer/recommendation
1: When applying the standard linearization to a BQP, should you reduce the size of the formulation based on the sign of the quadratic objective coefficients?
Our recommendation is to use the sign-based formulation STD′ when using the standard linearization.
2: When applying Glover’s formulation, should you represent the quadratic objective coefficient matrixC in upper triangular or symmetric form?
Our recommendation is to represent the quadratic objective coefficient matrixC in upper triangular form when using Glover’s method.
3: With regards to the overall computational effort to formulate and solve problem G1 to optimality using a MILP solver, should we compute the boundsUj1, Uj0, L1j, and L0j using (15), (16), or (17)?
In general, we recommend using the bounds of (16) when applying Glover’s method. However, for smaller-sized instances it may be beneficial to use the weaker bounds of (15).
4: How do formulations G1 and G2 compare when submitted to a MILP solver?
In general, formulation G2 is superior to G1.
5: Is it advantageous to perform the substitution of variables to reduce Problem G2 to G2a or G2b when submitting the model to a MILP solver?
While formulations G2, G2a, and G2b are fairly comparable in terms of performance, our general recommendation is to use G2a when applying Glover’s method.In conclusion, our recommendation based on a detailed computational study across five different classes of BQPs and three different MILP solvers is to implement Glover’s formulation using the concise version G2a with an upper triangular objective coefficient matrixC where the bounds within the formulation are computed as in (16).


### 4.3. Comparison of Exact Solution Methods

In this section, we compare our recommended versions of the standard linearization and Glover’s method, STD′, and G2a formulated with an upper triangular quadratic coefficient matrix C and the bounds of (16), with the more recent linearization of Sherali and Smith and direct submission to an optimization solver. When formulating the Sherali-Smith model, we utilized the bounds of (16) and an upper triangular quadratic coefficient matrix.As noted earlier, many MILP solvers have recently been updated to directly handle nonconvex BQPs. Before examining the results of our study, we review the two approaches typically utilized by commercial MILP solvers to solve nonconvex BQPs directly. The MILP solvers of CPLEX, GUROBI, and XPRESS employ one of two different strategies for solving nonconvex BQPs. The first strategy consists of transforming the nonconvex problem into an equivalent convex BQP, which can then be solved using standard convex quadratic programming techniques. This is accomplished by using the identityx=xTIx when x is binary. The quadratic portion of the objective function xTCx can therefore be replaced by xTC+ρIx−ρx. A number of strategies for selecting a ρ so that C+ρI will be negative semidefinite have been proposed, such as those by Billionet et al. [19]. The second strategy is to employ a linearization technique to transform the nonconvex BQP into an equivalent linear form, such as the standard linearization.For this second part of our computational study, we randomly generated a completely new set of test problems from those used in Sections4.1 and 4.2. In addition, we also adjusted the sizes of the problem instances to include more challenging problems (although all other aspects of the instances were generated as presented in Section 3). With regards to instances of Problem QMKP with positive objective coefficients, n was varied from 110 to 140 in increments of 10 when k=1 and varied from 50 to 80 in increments of 10 when k=5 or 10. For instances of QMKP with objective coefficients of mixed sign, n was varied from 40 to 70 in increments of 10 for k=1,5,or10. We varied n from 20 to 50 in increments of 10 for Problem HSP. For instances of kQKP with positive objective coefficients, n was varied from 60 to 90 in increments of 10, while for those instances with objective coefficients of mixed sign, we varied n from 80 to 110 in increments of 10. With regards to Problem QSAP, we considered ten m,n pairs where n∈4,5 and m∈10,12,15,18,20. For BLSP, we maintained the same sizes as before since this set was already challenging.The performance profiles of our results are presented in Figure9, where we once again aggregate across all objective coefficient densities and sizes. Note that in order to ensure that the number of problem instances np is the same for each of the formulations as required to construct performance profiles, we only solved the linearizations STD′, G2a, and SS with CPLEX, as opposed to all three solvers as in Sections 4.1 and 4.2. We report the average CPU times in Figure 10, which were computed as before. That is, these average CPU times were computed after removing all instances where at least one of the formulations did not solve within the time limit. Thus, the averages only include instances in which all formulations were solved to optimality within the 1-hour time limit.Figure 9
Performance profiles based on relative CPU time (s) for all approaches. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 10
Average CPU times (s) for all approaches.While we had hoped that a single strategy would emerge as the best approach for solving all BQPs regardless of problem class, the results presented in Figures9 and 10 clearly indicate that this is not the case. Indeed, different strategies appear to be superior for each of the different problem classes. We begin by noting that in general, direct submission of the BQP does not appear to be competitive with submission of any of the linear reformulations. In fact, the average CPU times are typically much higher when using the quadratic solvers of CPLEX, GUROBI, and XPRESS as opposed to solving the linearizations. Moreover, the quadratic solvers were more likely to not be able to solve all instances. For example, the XPRESS solver was unable to solve approximately 55% of the instances of QMKP with objective coefficients of mixed sign. Of the three solvers, CPLEX performed the best for direct submission, which is not surprising given that the nonconvex BQP solver of CPLEX is more mature than those of GUROBI or XPRESS (CPLEX has had the ability to solve nonconvex BQP well before GUROBI or XPRESS). Overall, our tests suggest that while it is convenient to be able to directly submit a BQP to the solver, solution times may be significantly reduced by reformulating the problem into linear form by the user before applying the solver.We now examine the performance profiles of Figure9 to compare the linearizations STD′, G2a, and SS. The standard linearization STD′ had the highest probability of being the fastest formulation for instances of QMKP with objective coefficients of mixed sign and instances of HSP. Glover’s linearization G2a had the highest probability of being the best formulation to solve instances of BLSP and instances of kQKP with both positive objective coefficients and objective coefficients of mixed sign. The linearization of Sherali and Smith SS had the largest probability of being the fastest formulation for instances of QMKP with positive objective coefficients, instances of kQKP with objective coefficients of mixed sign, and instances of QSAP. When examining the data a little more closely, G2a was the top performing formulation for instances of HSP that were 75% dense, instances of QMKP with objective coefficients of mixed sign that had densities of 75% and 100%, and instances of kQKP with objective coefficients of mixed sign that were 100% dense.In terms of the average CPU times in Figure10, there were no significant differences between the three linearizations. All of the average solution times were fairly close with the exception of STD′ applied to instances of QSAP. Indeed, STD′ performed poorly for this class of BQP.We end our analysis by aggregating the test results for all of the classes of BQPs into a single data set. The performance profiles of the three linearizations are presented in Figure11, while the average CPU times are presented in Figure 12.Figure 11
Performance profiles based on relative CPU time (s) over all problem classes.Figure 12
Average CPU times (s) over all problem classes.We can see that G2a had the fastest solution time approximately 26% of the time, STD′ was the fastest reformulation approximately 38% of the time, while SS was the superior formulation approximately 34% of the time. Interestingly, in terms of average overall CPU time, we see a different pattern. Note that STD′ required an average overall CPU time of 70 seconds, while G2a only required an average of 28 seconds. While we omit the details for brevity, a closer examination of the data suggests that STD′ is better suited for smaller-sized instances of BQPs, likely because of the large number of auxiliary variables and constraints. For larger-sized instances of BQPs, the compactness of G2a and SS seems to provide an advantage when submitted to an MILP solver.


## 4.1. Results for the Standard Linearization

We begin with our first question raised: When implementing the standard linearization, should you utilize the full formulation or use the formulation that is based on the sign of the quadratic objective coefficients? That is, how does STD compare with STD′ when submitted to an MILP solver?Figure1 shows the performance obtained for STD and STD′ for each class of BQP, where we aggregated all of the sizes and densities for each problem class, along with the three different MILP solvers. Note that since STD and STD′ vary depending on the sign of the objective coefficients, we provide two performance profiles for both QMKP and kQKP, one for instances with positive quadratic coefficients and one for instances with quadratic coefficients of mixed sign. For the reader unfamiliar with performance profiles, let us examine the results for the BLSP presented in Figure 1(a) in detail. Here, nf=2 as we are comparing the two formulations STD and STD′, while np=5 sizes n=30,40,50,60,70×10instances×3MILPsolvers=150. Each of these 150 instances is then solved using STD and STD′. For each instance, the minimum computational time of the two formulations is determined, and then the computational time of each formulation is divided by the minimum time, providing the performance ratio. Using the calculation of the performance ratios for all instances, the cumulative profile is then calculated for a set of fixed ratios of time factors. Recall that ρf1 is the probability formulation f will have the fastest solve time for any given problem. Thus, the probability that STD′ will have the fastest solve time for BLSP is approximately 86%, whereas the probability that STD will be the fastest formulation is approximately 14%.Figure 1
Performance profiles based on relative CPU time (s) for Question 1. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The performance profiles clearly indicate that STD′ is more likely to have the fastest solve time for all of the BQPs under consideration. While not presented in Figure 1, when aggregated across all problem classes, the probability that STD′ is the fastest formulation was 92%. We also note that a more detailed analysis showed that STD′ was the superior formulation regardless of objective coefficient density, problem size, or MILP solver—we omit the details for brevity.While performance profiles provide information on how likely a particular formulation will have a faster solve time, they do not indicate the actual differences in solve time. We report the average CPU times to solve the formulations to proven optimality for the different BQP classes in Figure2. These average CPU times were computed after removing all instances where at least one of the formulations did not solve within the time limit. That is, the averages only include instances in which both formulations were solved to optimality within the 1-hour time limit. Note that the performance profiles indicate that the number of instances that were not solved within the time limit was fairly small. To see this, recall that those profiles that do not meet ρfτ=1 indicate that there is a probability that they will not solve a subset of problems to optimality within the time limit. For example, note that approximately 5% of the instances of HSP formulated with STD did not solve within the 1-hour limit. The data in Figure 2 indicates that STD′ is significantly faster than STD on average for all BQP classes considered. Indeed, for 5 of the 7 classifications, the average CPU time of STD′ is less than half of that of STD.Figure 2
Average CPU times (s) for Question 1.Thus, our results clearly indicate that when utilizing the standard linearization, it is best to use the sign-based version STD′ regardless of BQP class, objective function density, problem size, or MILP solver. It is interesting to note that STD has a continuous relaxation that is at least as tight as STD′ but comes at the expense of a larger problem size. Our tests clearly indicate that the stronger continuous relaxation of STD does not provide an overall computational advantage.


## 4.2. Results for Glover’s Formulation

In this section, we investigate the four questions posed regarding Glover’s formulation. While ideally we would investigate these questions simultaneously, doing so would require a prohibitive number of test cases. Therefore, we investigate the questions in an incremental manner to make the task more tractable.We begin by focusing on answering the second and fourth questions raised. That is, when formulating Glover’s model, should we record the quadratic objective coefficient matrix in a symmetric or upper triangular form, and should we reduce the size of G1 by removing constraints (11) and (13) to obtain the modified Glover formulation G2? In order to compare the use of a symmetric quadratic coefficient matrix C with that of one in upper triangular form when formulating G1 and G2, we replaced Cij in the objective function of BQP with Cij′ defined as follows. To obtain an upper triangular objective coefficient matrix, we set Cij′=Cij+Cji for all i,j with i<j and Cij′=0 for i≥j. To obtain a symmetric objective coefficient matrix, we set Cij′=Cji′=1/2Cij+Cji) for all i,j with i<j and Cii′=0 for all i.Figure3 shows the performance of the four different formulations under consideration: G1 and G2 formulated with both symmetric and upper triangular objective coefficient matrices. We once again aggregate across all objective coefficient densities, sizes, and MILP solvers. Note that when formulating G1 and G2, we computed the bounds within the formulations as in (16) since preliminary computational experience showed this method provides the best performance. However, we will formally compare the bounds of (15), (16), and (17) in another set of tests. We also note that when computing solution times, we included the time needed to find the bounds within G1 and G2.Figure 3
Performance profiles based on relative CPU time (s) for Questions 2 and 4. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)The results in Figure3 provide a fairly definitive answer to question 4—Problem G2 is superior to G1 when compared using the same objective coefficient representation. Indeed, G2 with objective coefficients in upper triangular form had the highest probability of being the fastest formulation for all classes of BQPs, while G2 with objective coefficients in symmetric form had the second highest probability of being the fastest formulation for all classes except BLSP and QSAP. With regards to Question 2, note that both G1 and G2 formulated with objective coefficients in upper triangular form were more likely to be faster than their respective versions formulated with symmetric objective coefficients, with the exception of instances of kQKP with positive objective coefficients. Interestingly, we note that a more detailed analysis (not presented here) showed that G2 with symmetric objective coefficients was also superior for higher density (Δ=75%and100%) instances of QMKP with positive coefficients. Note that the performance profiles also indicate that some of the formulations struggled with being able to solve all of the instances within the time limit. For example, G1 and G2 formulated with symmetric objective coefficients were unable to solve approximately 50% of the instances of QSAP, while none of the formulations were able to solve all of the instances of HSP.We report the average CPU times to solve the formulations to proven optimality for the different BQP classes in Figure4. As before, these average CPU times were computed after removing all instances where at least one of the formulations did not solve within the time limit.Figure 4
Average CPU times (s) for Questions 2 and 4.Figure4 provides some additional information not captured in the performance profiles. First, note that while the performance profiles in Figure 3 indicate that G2 was superior to G1, Figure 4 indicates that the actual difference in average solution times between G1 and G2 is fairly negligible when they are formulated with the same objective representation. However, there are stark differences between the symmetric and upper triangular coefficient representations for both G1 and G2. In particular, instances of G1 and G2 formulated with symmetric objective coefficients typically take considerably longer to solve on average than their respective models formulated with upper triangular coefficients. There are two exceptions—G2 with symmetric objective coefficients is superior for HSP and instances of kQKP with positive coefficients.While the results are not completely definitive, Figures3 and 4 indicate that in general, constraints (11) and (13) should be dropped from G1 to obtain G2 and the quadratic objective coefficients C should be represented in upper triangular form. We reiterate that while the average CPU time differences between G1 and G2 are fairly small, those differences are rather large for the two different objective function representations.We now attempt to answer Question 5, which is concerned with whether it is advantageous to perform the substitution of variables that reduces G2 to either G2a or G2b. Figure5 shows the performance profiles of G2, G2a, and G3b where we aggregate across all objective coefficient densities, sizes, and MILP solvers. As with our previous test cases, we computed the bounds within the formulations as in (16).Figure 5
Performance profiles based on relative CPU time (s) for Question 5. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Note that, typically, G2 was the slowest formulation while G2a and G2b had fairly comparable likelihoods of being the fastest. Although, a closer examination reveals that G2a was slightly more likely to be the fastest formulation for all problem classes except instances of QMKP with positive objective coefficients. We report the average CPU times in Figure6, which were computed as before.Figure 6
Average CPU times (s) for Question 5.While the average CPU times for G2, G2a, and G2b are all fairly similar for each problem classification, G2a tended to be the fastest formulation. Therefore, based on our tests, we recommend that G2 should be reduced in size via the substitution of variablessj=Uj1xj−zj to obtain G2a. While the results were certainly not definitive, they do suggest that G2a has the highest probability of being the fastest formulation.We now address Question 3 regarding how the bounds within Glover’s formulation should be computed. Within our tests, we refer to the bounds computed as in (15) as weak, those computed as in (16) as tight, and those as in (17) as tightest. Based upon our prior results, our tests were all performed using G2a formulated with the three different methods to compute the bounds. We again reiterate that our solution times include not only the time to solve G2a to integer optimality but also the time required to compute the bounds within the formulation. Our results are presented in Figure 7 where we once again aggregate across all objective coefficient densities, sizes, and MILP solvers. Note that we do not include the BLSP instances because for these problems, the bounds (15), (16), and (17) are all the same.Figure 7
Performance profiles based on relative CPU time (s) for Question 3. (a) QMKP with positive coefficients, (b) QMKP with mixed-sign coefficients, (c)kQKP with positive coefficients, (d) kQKP with mixed coefficients, (e) HSP, (f) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)Interestingly, note that G2a formulated with the weak bounds of (15) tended to be the fastest formulation. We found this result surprising as preliminary experience showed that the tight bounds of (16) were superior. A more detailed analysis revealed that for smaller-sized instances, the extra time needed to find the tighter bounds of (16) was not warranted. However, once the instances start to become more challenging, the additional strength afforded by the tight bounds (16) was advantageous even though the time to find the bounds increases. In order to better understand how these bounds affect CPU time, let us examine the average CPU times presented in Figure 8.Figure 8
Average CPU times (s) for Question 3.Note that in terms of average CPU times, G2a formulated with the tight bounds (16) was the best formulation for all problem classifications except instances of kQKP with positive objective coefficients. Moreover, not only did the weak bounds (15) perform poorly in terms of average CPU times, they typically performed significantly worse than the tighter bounds. This stems from the fact that the continuous relaxation strength of G2a formulated with (15) for anything other than smaller-sized instances is too weak for the linearization to be effective. Thus, our general recommendation is to formulate Glover’s model with the tight bounds of (16), especially for larger-sized instances.We summarize the answers and recommendations from Questions 1 through 5 in Table1.Table 1
Summary of answers to Questions 1 through 5.
Question
Answer/recommendation
1: When applying the standard linearization to a BQP, should you reduce the size of the formulation based on the sign of the quadratic objective coefficients?
Our recommendation is to use the sign-based formulation STD′ when using the standard linearization.
2: When applying Glover’s formulation, should you represent the quadratic objective coefficient matrixC in upper triangular or symmetric form?
Our recommendation is to represent the quadratic objective coefficient matrixC in upper triangular form when using Glover’s method.
3: With regards to the overall computational effort to formulate and solve problem G1 to optimality using a MILP solver, should we compute the boundsUj1, Uj0, L1j, and L0j using (15), (16), or (17)?
In general, we recommend using the bounds of (16) when applying Glover’s method. However, for smaller-sized instances it may be beneficial to use the weaker bounds of (15).
4: How do formulations G1 and G2 compare when submitted to a MILP solver?
In general, formulation G2 is superior to G1.
5: Is it advantageous to perform the substitution of variables to reduce Problem G2 to G2a or G2b when submitting the model to a MILP solver?
While formulations G2, G2a, and G2b are fairly comparable in terms of performance, our general recommendation is to use G2a when applying Glover’s method.In conclusion, our recommendation based on a detailed computational study across five different classes of BQPs and three different MILP solvers is to implement Glover’s formulation using the concise version G2a with an upper triangular objective coefficient matrixC where the bounds within the formulation are computed as in (16).


## 4.3. Comparison of Exact Solution Methods

In this section, we compare our recommended versions of the standard linearization and Glover’s method, STD′, and G2a formulated with an upper triangular quadratic coefficient matrix C and the bounds of (16), with the more recent linearization of Sherali and Smith and direct submission to an optimization solver. When formulating the Sherali-Smith model, we utilized the bounds of (16) and an upper triangular quadratic coefficient matrix.As noted earlier, many MILP solvers have recently been updated to directly handle nonconvex BQPs. Before examining the results of our study, we review the two approaches typically utilized by commercial MILP solvers to solve nonconvex BQPs directly. The MILP solvers of CPLEX, GUROBI, and XPRESS employ one of two different strategies for solving nonconvex BQPs. The first strategy consists of transforming the nonconvex problem into an equivalent convex BQP, which can then be solved using standard convex quadratic programming techniques. This is accomplished by using the identityx=xTIx when x is binary. The quadratic portion of the objective function xTCx can therefore be replaced by xTC+ρIx−ρx. A number of strategies for selecting a ρ so that C+ρI will be negative semidefinite have been proposed, such as those by Billionet et al. [19]. The second strategy is to employ a linearization technique to transform the nonconvex BQP into an equivalent linear form, such as the standard linearization.For this second part of our computational study, we randomly generated a completely new set of test problems from those used in Sections4.1 and 4.2. In addition, we also adjusted the sizes of the problem instances to include more challenging problems (although all other aspects of the instances were generated as presented in Section 3). With regards to instances of Problem QMKP with positive objective coefficients, n was varied from 110 to 140 in increments of 10 when k=1 and varied from 50 to 80 in increments of 10 when k=5 or 10. For instances of QMKP with objective coefficients of mixed sign, n was varied from 40 to 70 in increments of 10 for k=1,5,or10. We varied n from 20 to 50 in increments of 10 for Problem HSP. For instances of kQKP with positive objective coefficients, n was varied from 60 to 90 in increments of 10, while for those instances with objective coefficients of mixed sign, we varied n from 80 to 110 in increments of 10. With regards to Problem QSAP, we considered ten m,n pairs where n∈4,5 and m∈10,12,15,18,20. For BLSP, we maintained the same sizes as before since this set was already challenging.The performance profiles of our results are presented in Figure9, where we once again aggregate across all objective coefficient densities and sizes. Note that in order to ensure that the number of problem instances np is the same for each of the formulations as required to construct performance profiles, we only solved the linearizations STD′, G2a, and SS with CPLEX, as opposed to all three solvers as in Sections 4.1 and 4.2. We report the average CPU times in Figure 10, which were computed as before. That is, these average CPU times were computed after removing all instances where at least one of the formulations did not solve within the time limit. Thus, the averages only include instances in which all formulations were solved to optimality within the 1-hour time limit.Figure 9
Performance profiles based on relative CPU time (s) for all approaches. (a) BLSP, (b) QMKP with positive coefficients, (c) QMKP with mixed-sign coefficients, (d)kQKP with positive coefficients, (e) kQKP with mixed coefficients, (f) HSP, (g) QSAP.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 10
Average CPU times (s) for all approaches.While we had hoped that a single strategy would emerge as the best approach for solving all BQPs regardless of problem class, the results presented in Figures9 and 10 clearly indicate that this is not the case. Indeed, different strategies appear to be superior for each of the different problem classes. We begin by noting that in general, direct submission of the BQP does not appear to be competitive with submission of any of the linear reformulations. In fact, the average CPU times are typically much higher when using the quadratic solvers of CPLEX, GUROBI, and XPRESS as opposed to solving the linearizations. Moreover, the quadratic solvers were more likely to not be able to solve all instances. For example, the XPRESS solver was unable to solve approximately 55% of the instances of QMKP with objective coefficients of mixed sign. Of the three solvers, CPLEX performed the best for direct submission, which is not surprising given that the nonconvex BQP solver of CPLEX is more mature than those of GUROBI or XPRESS (CPLEX has had the ability to solve nonconvex BQP well before GUROBI or XPRESS). Overall, our tests suggest that while it is convenient to be able to directly submit a BQP to the solver, solution times may be significantly reduced by reformulating the problem into linear form by the user before applying the solver.We now examine the performance profiles of Figure9 to compare the linearizations STD′, G2a, and SS. The standard linearization STD′ had the highest probability of being the fastest formulation for instances of QMKP with objective coefficients of mixed sign and instances of HSP. Glover’s linearization G2a had the highest probability of being the best formulation to solve instances of BLSP and instances of kQKP with both positive objective coefficients and objective coefficients of mixed sign. The linearization of Sherali and Smith SS had the largest probability of being the fastest formulation for instances of QMKP with positive objective coefficients, instances of kQKP with objective coefficients of mixed sign, and instances of QSAP. When examining the data a little more closely, G2a was the top performing formulation for instances of HSP that were 75% dense, instances of QMKP with objective coefficients of mixed sign that had densities of 75% and 100%, and instances of kQKP with objective coefficients of mixed sign that were 100% dense.In terms of the average CPU times in Figure10, there were no significant differences between the three linearizations. All of the average solution times were fairly close with the exception of STD′ applied to instances of QSAP. Indeed, STD′ performed poorly for this class of BQP.We end our analysis by aggregating the test results for all of the classes of BQPs into a single data set. The performance profiles of the three linearizations are presented in Figure11, while the average CPU times are presented in Figure 12.Figure 11
Performance profiles based on relative CPU time (s) over all problem classes.Figure 12
Average CPU times (s) over all problem classes.We can see that G2a had the fastest solution time approximately 26% of the time, STD′ was the fastest reformulation approximately 38% of the time, while SS was the superior formulation approximately 34% of the time. Interestingly, in terms of average overall CPU time, we see a different pattern. Note that STD′ required an average overall CPU time of 70 seconds, while G2a only required an average of 28 seconds. While we omit the details for brevity, a closer examination of the data suggests that STD′ is better suited for smaller-sized instances of BQPs, likely because of the large number of auxiliary variables and constraints. For larger-sized instances of BQPs, the compactness of G2a and SS seems to provide an advantage when submitted to an MILP solver.


## 5. Conclusions and Recommendations

In this paper, we revisited the standard linearization [3] and Glover’s method [4] for reformulating a BQP into linear form. Using a large-scale computational study over five different classes of BQPs, we evaluated a number of enhancements that have been proposed for these two linearizations. Based on our analysis, we recommend that when reformulating a BQP using the standard linearization, the sign-based version STD′ should be utilized. When applying Glover’s linearization, we recommend using the concise version G2a with an upper triangular objective coefficient matrix C where the bounds within the formulation are computed as in (16).In the second part of our study, we compared the STD′ and G2a formulations with the more recent linearization of Sherali and Smith and direct submission of the BQP to the solver. Our analysis showed that while it is convenient to submit a BQP directly to the solver, in general, it is advantageous for the user to perform a linear reformulation themselves beforehand. That being said, the more mature nonconvex BQP solver of CPLEX fared acceptably well on a number of different problem classes. Among the linear reformulations studied, no linearization was superior over all problem classes in terms of the probability of having the fastest solution times. In terms of overall average CPU time across all classes of BQPs, G2a and SS yielded the best times. Based upon our analysis, we recommend using STD′ for smaller-sized instances of BQPs, while for larger-sized instances, we recommend using either G2a or SS.



---

*Source: 5974820-2020-04-27.xml*