Using Nonlinear Mixed-integer Programming

In [1] we explored some linear mixed-integer programming formulations for the Least Median of Squares regression problem. Now let’s look at some nonlinear formulations.

Quadratic model

A mixed-integer quadratically constrained formulation can look like [1]:

\[\begin{align}\min\>&z\\&z\ge r^2_i-M\delta_i\\ &\sum_i \delta_i = n-h\\&r_i = y_i-\sum_j X_{i,j} \beta_j\\&\delta_i \in \{0,1\}\end{align}\]

This MIQCP model is convex and can be solved with top-of-the-line commercial solvers like Cplex and Gurobi. Some timings I see are:

----    145 PARAMETER report                    obj        time       nodes big-M MIP       0.262      30.594  160058.000 MIQCP           0.069     712.156 1783661.000

The objectives look different but they are actually correct: in the MIP formulation the objective reflects \(|r|_{(h)}\) while the quadratic formulation uses \(r^2_{(h)}\).

This performance differential is not unusual. If there are good MIP formulations available, in most cases quadratic formulations are not competitive.

MINLP model

The MINLP model from [1]:

\[\begin{align}\min\>&z\\&z\ge \delta_i \cdot r^2_i\\ &\sum_i \delta_i = h\\&r_i = y_i-\sum_j X_{i,j} \beta_j\\&\delta_i \in \{0,1\}\end{align}\]

is more problematic. This is no longer convex, so we need a global solver. When we try to solve this model as is we are in a heap of problems:

----    193 PARAMETER report                        obj        time       nodes   modelstat MINLP/Couenne       0.452       5.562     261.000       1.000 (optimal) MINLP/Baron         0.225    1000.000   71811.000       8.000 (integer solution, exceeded time limit)

Obviously, the Couenne solution is not correct, and Baron is not doing so great either. Global solvers have troubles when variables do not have proper bounds. It would be good it all kind of alarm bells start ringing, but they don’t.

Let’s add the following bounds:

\[\begin{align}&-100\le \beta_j \le 100\\&-100\le r_i \le 100\end{align}\]

Of course in practice you may not have good bounds available. These bounds yield much better results:

----    193 PARAMETER report                        obj        time       nodes   modelstat MINLP/Couenne       0.069     108.875   10827.000       1.000 MINLP/Baron         0.069     170.890      77.000       1.000

This is actually quite good compared to the MIQCP model.

Initial solutions

Recently I was asked about the effect of an initial point on the behavior of Couenne. Let’s try the experiment where we do two solves in sequence. The result is:

----    181 PARAMETER report                    obj        time       nodes   modelstat Couenne.1       0.069     108.469   10827.000       1.000 Couenne.2       0.069     113.359   11782.000       1.000

This is surprising: the second solve is actually more expensive! I have no good explanation for this. May be the initial point only affected the first sub-problem and in a bad way. Now try Baron:

----    181 PARAMETER report                  obj        time       nodes   modelstat Baron.1       0.069     185.600      77.000       1.000 Baron.2       0.069      94.980      29.000       1.000

This looks much better. In addition the Baron log gives some clear messages about the initial point being used as an incumbent:

Starting solution is feasible with a value of    0.686748259896E-001   Doing local search   Solving bounding LP   Starting multi-start local search   Done with local search ===========================================================================    Iteration    Open nodes         Time (s)    Lower bound      Upper bound            1             1             8.25      0.00000         0.686748E-01            3+            3            36.98      0.00000         0.686748E-01            7+            3            66.42      0.00000         0.686748E-01           29             0            94.98     0.686747E-01     0.686747E-01 Cleaning up

We can see the effect of the initial point: the upper bound is immediately equal to the objective value of the initial point.

Setting a bound

What about instead of providing an initial point, just set a bound on \(z\):

\[z\le 0.069\]

This looks like a good strategy to help the solver.

----    193 PARAMETER report                  obj        time       nodes   modelstat Couenne       0.069     130.125   10104.000       1.000 Baron            NA     115.390      15.000      19.000 (infeasible, no solution)

Bummer: it does not help at all.

Solving global MINLP models is not yet as “automatic” as solving an LP. You need to be aware: solvers for these type of models are a little bit more fragile. A little bit of a sobering conclusion.

References

1. Integer Programming and Least Median of Squares Regression, http://yetanothermathprogrammingconsultant.blogspot.com/2017/11/integer-programming-and-least-median-of.html