Non-convex Quadratic Integer Programming

Here, I want to revisit a particular model from [1]:

Model 3: Quadratic Preemptive Model

\[\begin{align}\max\>&\color{darkred}z_{model3}=\sum_p \color{darkred}z_p \\ & \color{darkred}z_p = \sum_{p',g} \color{darkblue}{\mathit{pref}}_{p,p'}\cdot \color{darkred}{\mathit{assign}}_{p,g}\cdot\color{darkred}{\mathit{assign}}_{p',g} \\ & \color{darkblue}z_{model2}^* \le \color{darkred}z_p & \forall p\\ & \sum_g \color{darkred}{\mathit{assign}}_{p,g} = 1 & \forall p \\ & \sum_p \color{darkred}{\mathit{assign}}_{p,g} = \color{darkblue}{\mathit{groupSize}} & \forall g \\& \color{darkred}{\mathit{assign}}_{p,g} \in \{0,1\}\end{align}\]

Cplex was especially behaving strangely. This deserves a bit more investigation. \(\color{darkblue}z_{model2}^*\) is a constant (it is the objective function value of an earlier model; in our case \(\color{darkblue}z_{model2}^*=0\)). The model as specified here, where \(\color{darkblue}z_{model2}^* \le \color{darkred}z_p\) is implemented as a lower bound on \(\color{darkred}z_p\), is refused outright by Cplex. The message is: 

CPLEX Error 5002: 'zpdef(person1)' is not convex.

The name zpdef refers to the quadratic constraint. Indeed, it is not convex. It can be linearized, however. Cplex only wants to try linearizing quadratic terms when they are in the objective. Why? I am not sure. We can trick Cplex into accepting the model by replacing the lower bound with an explicit constraint. This allows Cplex to make the quadratic constraint convex. Which means: 

• Making it an inequality. Any quadratic equality constraint is non-convex.
• Convexify by making the quadratic coefficient matrix positive semi-definite. This can be done by perturbing the diagonal [2]. In this case: add \(v\) times the diagonal, where \(-v\) is the most negative eigenvalue. To compensate, we can add a linear term using \(-v\cdot x^2\), which is, in this special case, the same as \(-v \cdot x\). For binary variables, we have \(x_i=x^2_i\).

I am guessing a bit here about what Cplex is doing. But I think I am right. 

This convexification makes the model a MISOCP (Mixed-Integer Second-Order Cone Program), which unfortunately shows dismal performance.

----    182 PARAMETER report  performance statistics

                 MIQCP       MIQCP       MIQCP       MIQCP      MISOCP         MIP
                 Baron    Antigone        SCIP       Cplex       Cplex       Cplex

Variables      253.000     253.000     253.000     253.000     253.000    1279.000
Equations       79.000      79.000      79.000      79.000     115.000    3157.000
Nonzeros      1735.000    1735.000    1735.000    1735.000    1771.000    8713.000
Status         optimal      intsol     optimal  nosolution      intsol     optimal
Objective       79.000      79.000      79.000          NA      74.000      79.000
Gap%                         0.010                              99.986
Time            10.270    1901.000      42.000                7215.407       3.141
Nodes            1.000       3.000    1269.000              597571.000    6207.000
Iterations                   1.000  123268.000             6.176701E+7  200181.000

The first three columns show the performance of global MINLP solvers. The MIQCP/Cplex column shows that Cplex refused to solve this model. The MISOCP column is for the model where we replaced lower bounds with explicit constraints so that Cplex can convexify. We stop after a two-hour time limit. It still has a gap of 99% (this is not a typo). Finally, in the last column, we linearized by hand. Cplex is both the slowest and the fastest in this table, depending on the formulation.

It is interesting to see that global, non-convex MINLP solvers are doing a much better job on the original non-convex problem than Cplex on the convexified problem.

I believe Cplex never linearizes constraints, only the objective. Obviously, if Cplex had linearized instead of convexified, it would have done infinitely better on this model. 

Conclusions

• For some models, quadratic formulations are very natural. 
• It is interesting to see that global, non-convex MINLP solvers are doing a much better job on the original non-convex problem than Cplex on the convexified problem.
• I believe Cplex never linearizes constraints, only the objective. Obviously, if Cplex had linearized instead of convexified, it would have done infinitely better on this model. 
• Solvers should be clear about what they are doing: convexify or linearize. The difference in performance can be dramatic.
• The behavior and design decisions behind some solvers remain a mystery to me.
• This is a very simple model. The results are a bit confusing and it is difficult  to understand what solvers are doing with this model. For larger, more complex models, this may be even more difficult to understand. Quadratic optimization is not as straightforward as one would expect.

 

References

1. Equity in optimization models, https://yetanothermathprogrammingconsultant.blogspot.com/2024/10/equity-in-optimization-models.html
2. Quadratic programming with binary variables, https://yetanothermathprogrammingconsultant.blogspot.com/2018/11/quadratic-programming-with-binary.html

Appendix: GAMS Model

$onText     Investigate Cplex performance on model 3 of   https://yetanothermathprogrammingconsultant.blogspot.com/2024/10/equity-in-optimization-models.html    $offText   option mip=cplex, miqcp=cplex;   *----------------------------------------------------------------------- * problem size *-----------------------------------------------------------------------   Set    p 'persons' /person1*person36/    g 'groups'  /group1*group6/ ;   abort$(card(p)/card(g)<>round(card(p)/card(g)))     "number of persons should be a multiple of number of groups";   alias (p,p1,p2);   sets    ne(p1,p2) 'off diagonal'    lt(p1,p2) 'strict triangular' ; ne(p1,p2) = ord(p1)<>ord(p2); lt(p1,p2) = ord(p1)<ord(p2);   scalar groupsize 'number of persons in each group'; groupsize = round(card(p)/card(g));   *----------------------------------------------------------------------- * random sparse preferences *-----------------------------------------------------------------------   option seed = 45634; parameter pref(p1,p2) 'preferences of person p1'; pref(ne)$(uniform(0,1)<0.15) = uniformint(-3,3);   option pref:0:0:4; display pref;   *----------------------------------------------------------------------- * reporting macros *-----------------------------------------------------------------------   parameter report(*,*,*) 'performance statistics'; option report:3:1:2;   acronym optimal,nosolution,intsol;   $macro statistics(m,name,desc) \ assign.l(p,g) = round(assign.l(p,g));  \ option assign:0; display assign.l;     \ report('Variables',name,desc)  = m.numvar;  \ report('Equations',name,desc)  = m.numequ;  \ report('Nonzeros',name,desc)   = m.numnz;   \ report('Status',name,desc)     = m.modelstat;   \ report('Status',name,desc)$(m.modelstat=1) = optimal; \ report('Status',name,desc)$(m.modelstat=8) = intsol; \ report('Status',name,desc)$(m.modelstat=14) = nosolution; \ report('Objective',name,desc)  = m.objval;  \ report('Gap%',name,desc)$(m.modelstat=8)  = 100*abs(m.objval-m.objest)/abs(m.objval); \ report('Time',name,desc)       = m.resusd;  \ report('Nodes',name,desc)      = m.nodusd;  \ report('Iterations',name,desc) = m.iterusd; \ display report;   *----------------------------------------------------------------------- * Non-convex Model v1 *-----------------------------------------------------------------------   scalar z2 'previously found objective value for obj2' /0/;   variable    zp(p)  'sum preferences for each person'    z1     'obj value for obj1' ;   binary variable assign(p,g) 'assign person to group';   equations     obj1        'obj: maximize overall sum preferences'     zpdef(p)    'quadratic constraint: sum preferences each person'     eassign1(p) 'each person is assigned to exactly one group'     eassign2(g) 'size of group' ;   obj1.. z1 =e= sum(p,zp(p));   zpdef(p1)..    zp(p1) =e= sum((p2,g),pref(p1,p2)*assign(p1,g)*assign(p2,g));   * assignment constraints eassign1(p).. sum(g,assign(p,g)) =e= 1; eassign2(g).. sum(p,assign(p,g)) =e= groupsize;   * we want to enforce *  z2 = smin(p,zp(p)) * or *  z2 <= smin(p,zp(p)) * they are the same in practice zp.lo(p) = z2;   model m1/obj1,zpdef,eassign1,eassign2/;   *----------------------------------------------------------------------- * Solve with global MINLP solvers and Cplex *-----------------------------------------------------------------------   option miqcp = baron; solve m1 maximizing z1 using miqcp; statistics(m1,'MIQCP','Baron')   option miqcp = antigone; solve m1 maximizing z1 using miqcp; statistics(m1,'MIQCP','Antigone')   option miqcp = scip; solve m1 maximizing z1 using miqcp; statistics(m1,'MIQCP','SCIP')   option miqcp = Cplex; solve m1 maximizing z1 using miqcp; statistics(m1,'MIQCP','Cplex')     *----------------------------------------------------------------------- * Model v2 * replace lowerbound by explicit constraint  *-----------------------------------------------------------------------   * remove lowerbound zp.lo(p) = -INF;   * introduce constraints instead equation lobnd(p); lobnd(p).. z2 =l= zp(p);   model m2/m1,lobnd/;     *----------------------------------------------------------------------- * Solve with Cplex * assumption: convexified and solved as MISOCP *-----------------------------------------------------------------------   option miqcp = Cplex, reslim = 7200; solve m2 maximizing z1 using miqcp; statistics(m2,'MISOCP','Cplex')   *----------------------------------------------------------------------- * Model v3 * handcrafted linearization  *-----------------------------------------------------------------------   * reinstate lowerbound zp.lo(p) = z2;   positive variable assign2(p1,p2,g) 'assign combo (p1,p2) to group'; assign2.up(p1,p2,g) = 1;   equation    zpdef2(p) 'linearized version of zpdef'    lin1(p1,p2,g) 'linearization of assign(p1,g)*assign(p2,g)'    lin2(p1,p2,g) 'linearization of assign(p1,g)*assign(p2,g)'    lin3(p1,p2,g) 'linearization of assign(p1,g)*assign(p2,g)' ;   zpdef2(p1)..  zp(p1) =e= sum((p2,g),pref(p1,p2)*assign2(p1,p2,g));   set nz(p1,p2) 'pref<>0'; nz(p1,p2) = pref(p1,p2);   lin1(nz(p1,p2),g).. assign2(p1,p2,g) =l= assign(p1,g); lin2(nz(p1,p2),g).. assign2(p1,p2,g) =l= assign(p2,g); lin3(nz(p1,p2),g).. assign2(p1,p2,g) =g= assign(p1,g)+assign(p2,g)-1;   model m3 /obj1,zpdef2,lin1,lin2,lin3,eassign1,eassign2/;     *----------------------------------------------------------------------- * Solve with Cplex *-----------------------------------------------------------------------   option mip = Cplex; solve m3 maximizing z1 using mip; statistics(m3,'MIP','Cplex')