## Tuesday, November 21, 2017

### Fun with indicator constraints

In [1] a simple model is presented for the Least Median of Squares regression problem:

 \begin{align}\min\>&z\\&\delta_i=1 \Rightarrow \> –z \le r_i \le z\\&\sum_i \delta_i = h\\&r_i = y_i - \sum_j X_{i,j} \beta_j\\&\delta_i \in \{0,1\} \end{align}

Here $$X$$ and $$y$$ are data. Modern solvers like Cplex and Gurobi support indicator constraints. This allows us to directly pass on the above model to a solver without resorting to other formulations (such as binary variables with big-M constraints or SOS1 structures).

This particular model is interesting: it has an unbounded LP relaxation. E.g. the start of the Cplex log shows:

 Nodes                                         Cuts/    Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap       0     0     unbounded                                          0               0     2     unbounded                                          0         Elapsed time = 0.03 sec. (1.27 ticks, tree = 0.01 MB, solutions = 0) *    44    44      integral     0        0.4886       -0.0000       63  100.00% Found incumbent of value 0.488611 after 0.05 sec. (3.02 ticks) *    45    43      integral     0        0.4779       -0.0000       66  100.00%. . .

Gurobi has some real problems on this model:

 Optimize a model with 48 rows, 97 columns and 188 nonzeros Model has 94 general constraints Variable types: 50 continuous, 47 integer (47 binary) Coefficient statistics:   Matrix range     [1e+00, 5e+00]   Objective range  [1e+00, 1e+00]   Bounds range     [1e+00, 1e+00]   RHS range        [4e+00, 2e+01] Presolve added 47 rows and 47 columns Presolve time: 0.00s Presolved: 95 rows, 144 columns, 423 nonzeros Presolved model has 94 SOS constraint(s) Variable types: 97 continuous, 47 integer (47 binary) Root relaxation: unbounded, 77 iterations, 0.00 seconds     Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time      0     0  postponed    0               -          -      -     -    0s     0     0  postponed    0               -          -      -     -    0s     0     2  postponed    0               -          -      -     -    0s *  136   136              68       0.9658407          -      -   6.1    0s *  138   136              69       0.9135526          -      -   6.4    0s H  363   353                       0.8212338          -      -   5.1    0s H  997   786                       0.6097619          -      -   5.5    0s H 1414   964                       0.5006481          -      -   6.9    0s H 1564   892                       0.3806481          -      -   8.7    0s H 1713   931                       0.3611765          -      -   8.8    0s...139309129 6259738    0.02538   65   23    0.26206          -      -  28.3 28775s 139337799 6262622     cutoff   67         0.26206          -      -  28.3 28780s 139367156 6266141     cutoff   68         0.26206          -      -  28.3 28785s 139395590 6268935     cutoff   67         0.26206          -      -  28.3 28790s 139424342 6272187  postponed   64         0.26206          -      -  28.3 28795s

This model does not seem to finish (the model has only 47 discrete variables so this log indicates something is seriously wrong: this model should solve in a manner of seconds). Gurobi was never able to establish a lower bound (column BestBd). Notice that Gurobi translates the implications to SOS1 variables (see [1] for the SOS1 version of this model).

There is a simple work around: add the bound $$z \ge 0$$. Now things solve very fast. We have a normal bound and no more any of these “postponed” nodes.

 Root relaxation: objective 0.000000e+00, 52 iterations, 0.00 seconds     Nodes    |    Current Node    |     Objective Bounds      |     Work Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time      0     0    0.00000    0   24          -    0.00000      -     -    0s H    0     0                       0.7556195    0.00000   100%     -    0s     0     2    0.00000    0   24    0.75562    0.00000   100%     -    0s H   54    30                       0.5603846    0.00000   100%  19.3    0s. . .

The bound $$z\ge 0$$ can be deduced from the constraints $$–z \le r_i \le z$$ (we know $$h$$ of them have to hold). Presolvers are not able to find this out automatically.

See [2] for some other interesting observations on indicator constraints.

Update: Gurobi has identified the problem and it has been fixed (available in next version).

##### References
1. Integer Programming and Least Median of Squares Regression, http://yetanothermathprogrammingconsultant.blogspot.com/2017/11/integer-programming-and-least-median-of.html
2. Belotti, P., Bonami, P., Fischetti, M. et al., On handling indicator constraints in mixed integer programming, Comput Optim Appl (2016) 65: 545.