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.

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 -item quadratic knapsack problem, the heaviest -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 where and are the linear and quadratic cost coefficients, respectively, and is a polyhedral set representing the feasible region. For notational ease, we henceforth let the indices and run from 1 to unless otherwise stated. Note that for each binary variable, , and thus, can be incorporated into the linear term without loss of generality. Furthermore, note that the overall cost associated with the product is . Therefore, even if and 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 for all (so that the quadratic coefficient matrix is symmetric) or assume for all (so that the quadratic coefficient matrix is upper triangular); both of which can be assumed without loss of generality. For our purposes, we will explicitly keep both terms and . Note that there is no assumption that 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 in the objective function with a continuous variable . In order to simplify the presentation, we first rewrite the objective function of BQP as where for all with , which we can do without loss of generality. We now define the standard linearization of BQP below. subject to

Note that (4)–(7) enforce that for all binary . While this linearization method is straightforward, it has the disadvantage that it requires the addition of auxiliary variables and auxiliary constraints.

As noted in [8], we can reduce the number of auxiliary constraints by using the sign of . Let : and . Then, we can omit the constraints (4) and (5) bounding from above for , and we can omit the constraints (6) and (7) bounding from below for , as these will be implied at optimality. We refer to this reduced form as and note that the continuous relaxation of may potentially be weaker than that of .

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 ) found in the objective function with a continuous variable and uses four linear restrictions to enforce that . Specifically, the Glover linearization is as follows. subject to

For each , and are upper bounds on , while and are lower bounds on . Problems BQP and G1 are equivalent in that, given any binary , constraints (10)–(13) ensure that ) for each . Observe that this formulation only requires the addition of (unrestricted) auxiliary continuous variables and auxiliary constraints and is therefore considerably more compact than the standard linearization.

As noted earlier, the two most common representations of the quadratic objective coefficients in the literature are to either assume for all (so that the quadratic coefficient matrix is symmetric) or assume for all (so that the quadratic coefficient matrix 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 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 can be removed from G1 along with the associated constraints in (10)–(13). To see this, note that when is upper triangular, (10)–(13) ensure that because for all . Thus, when the quadratic coefficient matrix is upper triangular, Glover’s formulation only requires the addition of auxiliary variables and auxiliary constraints. Furthermore, when is upper triangular, the auxiliary constraints of G1 are less dense than when 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 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 , they can easily be computed as

As shown in [9], stronger bounds that take into consideration the feasible region can be computed as for . These bounds could potentially be made even tighter by enforcing x binary so that

While the continuous relaxation of G1 could be potentially tightened by strengthening the values of the , , , and 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 , , , and 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 from below because of the maximization objective and the fact that the terms do not appear elsewhere in the problem. Thus, problem G1 can be written more concisely as the modified Glover G2: subject to

Interestingly, problem G2 only requires auxiliary constraints as opposed to the 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 variables or for all , which express the variables in terms of the slack variables to the inequalities (10) and (12), respectively (see [11]). Upon performing the first substitution, we obtain G2a: subject to

The second substitution yields G2b: subject to

Note that both G2a and G2b only require n (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.

subject to

For each , and are upper and lower bounds, respectively, on , 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, , 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 nonnegative auxiliary variables and 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 that is defined as

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 from a collection of noisy measurements. Traditionally, this problem is modeled as subject to where and are given. By ignoring the constant term , we can rewrite this formulation in our notation and in maximization form by setting , , and subsequently defining for , for , and for all .

We randomly generated instances of BLSP as in [1]. Specifically, we set and generated a random matrix with elements from the standard normal distribution and a binary vector with elements from the uniform distribution . The vector is then constructed as where is a random noise term with elements from . We generated 10 instances for each value of , 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. subject to

Within QMKP, the , , , and coefficients are typically nonnegative scalars. We assume for every that for each since otherwise variables can be fixed to 0 and that for the th 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 wherein for all , so that the problem reduces to linear form. The QKP retains the quadratic objective terms but has only a single knapsack constraint defining 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 with and 10 knapsack constraints. The coefficients are integers taken from a uniform distribution over the interval [1,50], and is an integer from a uniform distribution between 50 and .

For the first test bed, the objective coefficients were all nonnegative. Specifically, the (nonzero) objective coefficients for all and for with are integers from a uniform distribution over the interval [1,100]. To assess the effects of the density of nonzero and on CPU time, each of these coefficients is nonzero with some predetermined probability . We considered instances with probabilities (densities) For each value of , we generated ten instances for each value of n, which varied from 80 to 110 in increments of 10 when and varied from 30 to 70 in increments of 10 when or 10.

For the second test bed, the objective coefficients were mixed in sign. Specifically, the (nonzero) objective coefficients for all and for with are integers from a uniform distribution over the interval . As with the first test bed, the density of the objective coefficients was varied from – in increments of . For each value of , we generated 10 instances for each value of , which varied from 20 to 50 in increments of 10 for , and 10. We utilized smaller sizes as objective coefficients that are mixed in sign are significantly more difficult to solve.

3.3. The Heaviest -Subgraph Problem

The heaviest -subgraph problem (HSP) is concerned with determining a block of 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 subject to

The HSP is also known under the name of the cardinality constrained quadratic binary program, the densest -subgraph problem, the -dispersion-sum problem, and the -cluster problem. For each number of nodes, , we randomly generated 10 instances with three different graph densities, . Specifically, for each density and any pair of indexes such that , we assigned with probability and otherwise. We set for all of our tests.

3.4. -Item Quadratic Knapsack Problem

The -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: subject to

Here, we assume that for each i since otherwise variables can be fixed to 0 and that for the constraint to be restrictive. Let us denote by the largest number of items which can be filled in the knapsack, that is, the largest number of the smallest whose sum does not exceed . We can assume that .

Note that QKP 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 from the discrete uniform distribution over the range . For the first test bed, the objective coefficients were all nonnegative, where the (nonzero) objective coefficients for all and for with are integers from a uniform distribution over the interval [1,100]. The density of the objective coefficients was varied from – in increments of 25%. For each value of , we generated ten instances with variables. For the second test bed, the objective coefficients were mixed in sign, where the objective coefficients for all and for with are integers from a uniform distribution over the interval . As with the first test bed, the density of the objective coefficients was varied from in increments of 25%. For each value of , we generated 10 instances for each value of , where 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 is to be allocated to a set of processors . The execution costs for assigning task to processor is denoted by and the communication costs between on and on is denoted by . The constraints limit every task to exactly one processor. The specific formulation is as follows. subject to

For our tests, we considered two values of (4 and 5) and three values of (10, 12, 15). For each of the six pairs, we generated 10 instances where the coefficients and 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 using performance 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 set of formulations and a set of problem instances. For each problem and formulation , we define

We compare the performance of formulation on problem with the best performance by any formulation on this problem using the performance ratio

Then, the cumulative performance profile is defined as so that is the probability for formulation that a performance ratio is within a factor of the best possible ratio. Thus, for example, is the probability of formulation to be the best formulation to solve any given problem, while is the probability of formulation 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 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?

Figure 1 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 QKP, 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, as we are comparing the two formulations STD and STD, while sizes . 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 is the probability formulation 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%.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(a)
(b)
(c)
(d)
(e)
(f)
(g)

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.

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 Figure 2. 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 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 in the objective function of BQP with defined as follows. To obtain an upper triangular objective coefficient matrix, we set for all with and for . To obtain a symmetric objective coefficient matrix, we set ) for all with and for all .

Figure 3 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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(a)
(b)
(c)
(d)
(e)
(f)
(g)

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.

The results in Figure 3 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 QKP 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 () 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 Figure 4. 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.

Figure 4 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, Figures 3 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. Figure 5 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).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(a)
(b)
(c)
(d)
(e)
(f)
(g)

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.

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 Figure 6, 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 variables 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.

(a)

(b)

(c)

(d)

(e)

(f)

(a)
(b)
(c)
(d)
(e)
(f)

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.

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 Table 1.

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 matrix C 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 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 identity when is binary. The quadratic portion of the objective function can therefore be replaced by . A number of strategies for selecting a so that 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 Sections 4.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, was varied from 110 to 140 in increments of 10 when and varied from 50 to 80 in increments of 10 when or 10. For instances of QMKP with objective coefficients of mixed sign, was varied from 40 to 70 in increments of 10 for . We varied from 20 to 50 in increments of 10 for Problem HSP. For instances of kQKP with positive objective coefficients, was varied from 60 to 90 in increments of 10, while for those instances with objective coefficients of mixed sign, we varied from 80 to 110 in increments of 10. With regards to Problem QSAP, we considered ten pairs where and . For BLSP, we maintained the same sizes as before since this set was already challenging.

The performance profiles of our results are presented in Figure 9, where we once again aggregate across all objective coefficient densities and sizes. Note that in order to ensure that the number of problem instances 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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(a)
(b)
(c)
(d)
(e)
(f)
(g)

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.

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 Figures 9 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 Figure 9 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 QKP 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 QKP 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 QKP with objective coefficients of mixed sign that were 100% dense.

In terms of the average CPU times in Figure 10, 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.