On the solution of nonconvex cardinality Boolean quadratic programming problems: a computational study

Abstract

This paper addresses the solution of a cardinality Boolean quadratic programming problem using three different approaches. The first transforms the original problem into six mixed-integer linear programming (MILP) formulations. The second approach takes one of the MILP formulations and relies on the specific features of an MILP solver, namely using starting incumbents, polishing, and callbacks. The last involves the direct solution of the original problem by solvers that can accomodate the nonlinear combinatorial problem. Particular emphasis is placed on the definition of the MILP reformulations and their comparison with the other approaches. The results indicate that the data of the problem has a strong influence on the performance of the different approaches, and that there are clear-cut approaches that are better for some instances of the data. A detailed analysis of the results is made to identify the most effective approaches for specific instances of the data.

Appendices

Appendix 1: A note on formulation (F3)

The validity of formulation (F3) can be established as follows. The constraint $z_{ij} \ge x_{i}+ x_{j} - 1$ is not considered in (F3), which forces $z_{ij}=1$ when $x_{i}=1$ and $x_{j}=1$. In formulation (F3), this relation is enforced by constraint (10). To demonstrate this, consider that $M=3$, $x_{j}=1$, $x_{k}=1$, $x_{l}=1$, $x_{i}=0, \forall \, i \ne j,i\ne k,i\ne l$, and $j<k<l$ then the constraint in (10) results in

$$\begin{aligned}&\sum \limits _{s< j}{z_{sj}} + \sum \limits _{s>j}{z_{js}}= 2, \end{aligned}$$

(MIQP$_\textrm {F}$)

$$\begin{aligned}&\sum \limits _{s< k}{z_{sk}} + \sum \limits _{s>k}{z_{ks}}= 2, \end{aligned}$$

(QP$_{\textrm {Fu)

$$\begin{aligned}&\sum \limits _{s< l}{z_{sl}} + \sum \limits _{s>l}{z_{ls}}= 2, \end{aligned}$$

(MIQP$_{\text {Fu)

$$\begin{aligned}&\sum \limits _{s<i}{z_{is}} + \sum \limits _{s>i}{z_{is}}= 0, \qquad \forall \, i\ne j,i\ne k,i\ne l, \end{aligned}$$

(15)

where in each (15), (16), and (17) there are only two positive variables z, because the fourth set of equations forces all the other variables to zero, leading to $z_{jk}=1$, $z_{jl}=1$, $z_{kl}=1$, which is a solution of the original problem. Note that in (18), we have $x_{i}=0$ and then the $rhs=0$.

Appendix 2: Formulation (F6)

Formulation (F6) compared with the other formulations has some additional elements that require some clarification. Therefore, in this appendix, we derive formulation (F6), and we explain the meaning of the variable $t_{i}$, and the parameters $L_{i}$ and $U_{i}$. The derivation of the formulation is based on the linearization technique proposed in [43], and on one of the formulations derived in [9].

We start by defining the CBQP problem of interest as

$$\begin{aligned} \begin{array}{rl} \text{ min }&{} \sum \limits _{i}{q_{ii}x_{i}} + \sum \limits _{i}{\sum \limits _{j>i}{(q_{ij} + q_{ji})x_{i}x_{j}}} \\ \text{ subject } \text{ to } &{} \sum \limits _{i}{x_{i}} = M \\ &{} x_{i} \in \{0,1\}, \qquad \forall \, i. \end{array} \end{aligned}$$

(16)

Following Glover [43] we define a new variable $z_{i}$ that is only indexed by i as

$$\begin{aligned} z_{i} = x_{i}\left[ \sum \limits _{j<i}{(q_{ij} + q_{ji})x_{j}} + \sum \limits _{j>i}{(q_{ij} + q_{ji})x_{j}}\right] , \qquad \forall \, i, \end{aligned}$$

(17)

which is used to replace the quadratic term in the objective function of (19). Next the linearization technique is applied to the product between $x_{i}$ and the term within the square brackets in (20). For this linearization we need to define the lower and upper bounds for the term inside of the square brackets in (20):

$$\begin{aligned} L_{i} \le \sum \limits _{j<i}{(q_{ij} + q_{ji})x_{j}} + \sum \limits _{j>i}{(q_{ij} + q_{ji})x_{j}} \le U_{i}, \qquad \forall \, i. \end{aligned}$$

(18)

In the context of the problem (19), because of the cardinality constraint, $L_{i}$ is equal to the sum over j of the $M-1$ smallest coefficients $q_{ij} + q_{ji}$, and $U_{i}$ is equal to the sum over j of the $M-1$ greatest coefficients $q_{ij} + q_{ji}$, see [9] for additional details. Using the linearization technique from Glover [43] applied to (20), problem (19) can be re-written as

$$\begin{aligned} \begin{array}{rl} \text{ min }&{} \sum \limits _{i}{q_{ii}x_{i}} + \frac{1}{2} \sum \limits _{i}{z_{i}} \\ \text{ subject } \text{ to } &{} \sum \limits _{i}{x_{i}} = M \\ &{} z_{i} \ge \sum \limits _{j<i}{(q_{ij} + q_{ji})x_{j}} + \sum \limits _{j>i}{(q_{ij} + q_{ji})x_{j}} - U_{i}(1 - x_{i}), \qquad \forall \, i \\ &{} z_{i} \ge L_{i}x_{i}, \qquad \forall \, i \\ &{} x_{i} \in \{0,1\}, \qquad \forall \, i. \end{array} \end{aligned}$$

(19)

Note that the fraction 1 / 2 in the objective function of (22) is needed to prevent the double accounting of $(q_{ij} + q_{ji})$ in $z_{i}$ and $z_{j}$ when $x_{i}=1$ and $x_{j}=1$. Note also that because we are minimizing, we do not need to include in (22) the other two constraints derived in the linearization technique proposed in [43]:

$$\begin{aligned} z_{i}\le & {} \sum \limits _{j<i}{(q_{ij} + q_{ji})x_{j}} + \sum \limits _{j>i}{(q_{ij} + q_{ji})x_{j}} - L_{i}(1 - x_{i}), \qquad \forall \, i \end{aligned}$$

(20)

$$\begin{aligned} z_{i}\le & {} U_{i}x_{i}, \qquad \forall \, i. \end{aligned}$$

(21)

Defining $t_{i} = z_{i} - U_{i} x_{i}$ and replacing $z_{i}$ by $t_{i} + U_{i} x_{i}$ in problem (22), this problem can be re-written as

$$\begin{aligned} \begin{array}{rl} \min &{} \sum \limits _{i}{q_{ii}x_{i}} + \frac{1}{2} \sum \limits _{i}^{}{t_{i}} + \frac{1}{2} \sum \limits _{i}^{}{U_{i} x_{i}}\\ \text{ subject } \text{ to } &{} \sum \limits _{i}{x_{i}} = M \\ &{} t_{i} \ge \sum \limits _{j<i}{\left( q_{ji}+q_{ij}\right) x_{j}}+ \sum \limits _{j>i}{\left( q_{ij}+q_{ji}\right) x_{j}} - U_{i}, \qquad \forall \, i \\ &{} t_{i} \ge \left( L_{i} - U_{i}\right) x_{i}, \qquad \forall \, i \\ &{} x_{i} \in \{0,1\}, \qquad \forall \, i, \end{array} \end{aligned}$$

(22)

which leads to formulation (F6).

Appendix 3: Performance profiles

The performance profiles presented are based on the work of Dolan and Moré [52], who have proposed an alternative graphical comparison methodology to support the benchmarking of optimization software. This type of profiles are used in several benchmark tests, and they are an alternative to other comparison techniques based on averages and standard deviations, or other statistical indicators.

The performance profiles represent for a given set of results, P, the cumulative distribution function, $p_{s}(\tau )$, of the performance ratios between the computational time of the instance p solved with the solver s and the minimum of the computational times taken by the other solvers for the same instance. Thus the performance ratio is given by

$$\begin{aligned} \rho _{p,s} = \frac{t_{p,s}}{\min \{t_{p,s}: 1 \le s \le n_{s}\}}\;, \quad \forall \, p,s, \end{aligned}$$

(23)

and the cumulative performance profile is defined as

$$\begin{aligned} p_{s}(\tau ) = \frac{1}{n_{p}}\;\text {size} \left\{ p\in P: \rho _{p,s} \le \tau \right\} , \quad \forall \, s \end{aligned}$$

(24)

where $t_{p,s}$ is the computational time that solver s needed to solve problem p, $n_{s}$ is the number of solvers, $n_{p}$ is the number of instances, and $\tau$ denotes the ratio of the time factor of interest. The performance ratios with $\tau =1$ correspond to the solvers with the minimum computational time, and therefore, $p_s(1)$ is the probability of the solver s to have the best performance on any given problem [52].

Fig. 13

Example of performance profile built over the CPU time

As an example consider the performance profiles based on the CPU time in Figure 1, where the profiles for the MILP formulations are presented. These performance profiles are built over the formulations, and not over the solvers. This means that each formulation corresponds to an index s (denoted by solver in the nomenclature above), and thus $n_{s}=6$, and $n_{p}= 80\, \mathrm { instances} \times 3\, \mathrm {MILP\, solvers} = 240$. Each instance of these 240 is then solved for each formulation. For each instance the minimum computational time of the six formulations is calculated, and then the computational time of each formulation is divided by the minimum time, this is the performance ratio. Based on the calculation of the performance ratios for all instances, the cumulative profile is then calculated for a set of fixed ratios of time factors. $p_s(1)$ is the probability of the formulation s to be the best formulation to solve any given problem independently of the solver, and of the parameters of the 80 instances defined for matrix Q. For $\tau =4$, $p_s(4)$ is the probability of any formulation to be the best within a time factor of 4 of the best formulation. The performance profiles provide also information about the formulations that fail to solve to optimality a set of problems. This is given by $1 - p_{s}(\tau )$, and thus, the profiles that do not meet $p_{s}(\tau )= 1$, mean that there is a probability that they will not solve a subset of problems to optimality within the maximum CPU time set. Figure 13 illustrates a performance profile based on the computational times.

The performance profiles based on the absolute gap [53] are easier to understand and construct than the profiles based on the computational times describe above. A cumulative performance profile using absolute gaps for each solver s is defined as

$$\begin{aligned} p_{s}(Gap (\%)) = \frac{1}{n_{p}}\;\text {size} \left\{ p\in P: \gamma _{p,s} \le Gap (\%) \right\} , \quad \forall \, s, \end{aligned}$$

(F6)

where $Gap (\%)$ is the optimality gap of interest, and $\gamma _{p,s}$ is the optimality gap that solver s obtained for problem p. As an example consider the performance profiles based on the optimality gap in Fig. 1, where the profiles for the MILP formulations are presented. These performance profiles are built over the formulations, and not over the solvers. This means that each formulation corresponds to an index s (denoted by solver in the nomenclature above), and thus $n_{s}=6$, and $n_{p}= 80\, \mathrm { instances} \times 3\, \mathrm {MILP\, solvers} = 240$. For each formulation the cumulative profile of the optimality gaps of the 240 optimization runs are determined. $p_s(0)$ is the probability of the formulation s to solve any given problem to optimality independently of the solver, and of the parameters of the 80 instances defined for matrix Q. Figure 14 depicts an absolute profile for optimality gaps.

Fig. 14

Example of a performance profile based on the absolute gap