**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 |

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 |

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 |

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 |

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 |

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 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 |

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

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

## No comments:

## Post a Comment