HiGHS returns wrong solution?

[jajhall]

issue-1259.mps.txt returns a wrong solution with HiGHS:

Solving report
Status Optimal
Primal bound 241.252481008
Dual bound 241.2348756
Gap 0.0073% (tolerance: 0.01%)
Solution status feasible
241.252481008 (objective)
0 (bound viol.)
0 (int. viol.)
0 (row viol.)
Timing 0.21 (total)
0.01 (presolve)
0.00 (postsolve)
Nodes 1
LP iterations 404 (total)
15 (strong br.)
19 (separation)
336 (heuristics)

So it says that the optimal solution is 241.252481008

However when I solve this with other solvers (Lindo, scip, lpsolve) they all report a cheaper solution:

Lindo: 241.234839463
scip: 241.234839463423
lpsolve: 241.23483946

Regards,

Peter

[fuglede]

Perhaps interestingly, lowering mip_rel_gap is enough to get a dual solution with an objective value of 241.252481008, but mip_abs_gap seems to do nothing at all! Edit: I guess the convergence criterion is "mip_rel_gap or mip_abs_gap", where I would have guessed that both of them would have to be satisfied.

$ highs issue-1259.mps | grep -A 4 report
Solving report
  Status            Optimal
  Primal bound      241.252481008
  Dual bound        241.2348756
  Gap               0.0073% (tolerance: 0.01%)

$ echo "mip_rel_gap=0.0000001" > options && highs --options_file options issue-1259.mps | grep -A 4 report
Solving report
  Status            Optimal
  Primal bound      241.252481008
  Dual bound        241.252481008
  Gap               0% (tolerance: 1e-05%)

$ echo "mip_abs_gap=0.0000000001" > options && highs --options_file options issue-1259.mps | grep -A 4 report
Solving report
  Status            Optimal
  Primal bound      241.252481008
  Dual bound        241.2348756
  Gap               0.0073% (tolerance: 0.01%)

But SCIP does report OP's result for both bounds, so something is amiss:

$ scip -f issue-1259.mps 2>/dev/null | grep -A 2 " Primal Bound"
  Primal Bound     : +2.41234839463423e+02   (in run 1, after 59 nodes, 0.04 seconds, depth 8, found by <relaxation>)
  Dual Bound       : +2.41234839463423e+02
  Gap              :       0.00 %

[jajhall]

I, too, reduced mip_rel_gap to zero and got this outcome. Indeed, the termination criterion is that at least one of the two tolerances is satisfied: just to give flexibility to users. We can't expect both to be satisfied, as problems with very large optimal objective values would "never" satisfy mip_abs_gap when mip_rel_gap is small (and positive).

There's clearly a bug in the MIP solver, as HiGHS gets 241.252481008, even when --presolve=off

[erdembanak]

I have debugged the issue. The value of C19_F163_I79 is returned as 9.99e-07 from the analytic center computation, which is the upper bound (1.00E-6); therefore, it is fixed to that value while it should have been 0. It is possible to check this claim by fixing the variable in the MPS file.

The problem occurs in the post-process (https://github.com/ERGO-Code/HiGHS/blob/master/src/lp_data/HighsSolution.cpp#L623). In the analytic center computation, the column's value is 4.02e-07, and the value of the dual is -9.19e-07. Then, the below conditions are met:

      if (dual_infeasibility <= dual_feasibility_tolerance) continue;
      if (residual < dual_infeasibility && !force_dual_feasibility) {

(dual_feasibility_tolerance = 9.99e-08 and dual_infeasibility = 9.19e-07 since it is abs(dual)) and the variable is fixed at 9.99e-07.

I have experimented with changing the optimality_tol_ of the IPM from 1e-08 to 1e-10. The iteration count of the IPM increased from 14 to 17, and the new column value became 4.24e-07 while the dual is -1.18e-10. Since the dual_infeasibility is now less than the tolerance, the problem is solved correctly.

There are different ways of tackling the issue (decreasing the optimality_tol_, or scaling the problem since the residual can at most be 9.99e-07 for this variable). I can experiment with different options and report the results. Based on your suggestions, I would like to open a PR to solve this.

Thanks,
Erdem

Edit: Scaling the residual might be an easier option. I have played around with something like this, and it solved the problem:

      double residual = std::fabs(std::max(lower - value, value - upper));
      double var_range = upper - lower;
      double residual_rate = residual/var_range;
      residual = std::max(residual, residual_rate);

(There is already a control that checks upper > lower before). This looks safe, but I didn't check the performance implications.

[jajhall]

Thanks, but this won't work when (for example) one of lower and upper is infinite.

I'll have a look, but changes at this level of the HiGHS solvers have to be made with extreme care

[erdembanak]

I proposed it mainly as a starting point, but it works with infinity since we have std::max(residual, residual_rate). If the range is large, we fall back to using residual. You are right about saying it should be done with extreme care; this might end up with slowing down HiGHS since less variables might be fixed after the analytic center computation. Also there might be some other extreme cases I haven't thought about.

If there is anything I can help with the tests, I would like to help. I can try out different options based on discussions with you. I want to propose using HiGHS as an open-source solver for my company; therefore, I would like to help out with the bugs.

Thanks

[peno64]

I have another example where the difference is even bigger. Even with presolve off and a relative mip tolerance of 1e-8, HiGHS reports 241.5540085 as optimal value while other solvers like scip and Gurobi report 241.43565000
I guess it's the same issue.
HighsNotOptimal.mps.txt

[erdembanak]

Yes, it's the same issue. The proposed fix solves this problem too.