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.