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.