RelativeGap() calculation inconsistency

[jd-lara]

The documentation for the attribute RelativeGap() shows that it is calculated as \\frac{|b-f|}{|f|}, where b is the best bound and f is the best feasible objective value. However, it seems that the implementation is \\frac{|b-f|}{|b|}.

See the result from Gurobi's log after defining the MIPGap = 0.8 for testing purposes.

Optimize a model with 1202 rows, 1810 columns and 5539 nonzeros
Variable types: 1624 continuous, 186 integer (186 binary)
Coefficient statistics:
  Matrix range     [2e-02, 2e+02]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 8e+01]
  RHS range        [1e+00, 5e+01]
Presolve removed 232 rows and 242 columns
Presolve time: 0.01s
Presolved: 970 rows, 1568 columns, 4559 nonzeros
Variable types: 1382 continuous, 186 integer (186 binary)

Root relaxation: objective 8.497417e-02, 818 iterations, 0.02 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    0.08497    0    8          -    0.08497      -     -    0s
     0     0    0.08497    0    8          -    0.08497      -     -    0s
H    0     0                      47.0000000    0.08497   100%     -    0s
     0     2    0.08497    0    8   47.00000    0.08497   100%     -    0s
* 1240  1094              34      13.0000000    0.47609  96.3%  14.1    1s
H11008  7007                      10.0000000    1.20863  87.9%  12.7    4s
 11907  7773 infeasible  134        10.00000    1.21516  87.8%  12.4    5s
 41214 28804 infeasible  129        10.00000    1.31800  86.8%  10.9   10s
 69165 48345    5.46379  137   37   10.00000    1.36625  86.3%  10.7   15s
 93785 65011    2.40770   69   20   10.00000    1.40558  85.9%  10.6   20s
 124505 85422    2.64401  117   28   10.00000    1.44142  85.6%  10.5   25s
 154518 105008    2.40117   80   21   10.00000    1.47464  85.3%  10.6   30s
 182972 122830    2.55937  102   26   10.00000    1.50595  84.9%  10.6   35s
 207747 137727    8.50196   96   24   10.00000    1.53065  84.7%  10.8   40s
 222404 146142    2.74201  126   33   10.00000    1.54609  84.5%  10.8   45s
 249731 162469    2.51607   92   26   10.00000    1.57528  84.2%  10.9   50s
 275416 177177    2.58414  117   26   10.00000    1.60650  83.9%  11.0   55s
 296684 188804 infeasible  132        10.00000    1.63516  83.6%  11.2   60s
 315551 198991    3.62581  118   26   10.00000    1.65984  83.4%  11.2   65s
 333312 208174    2.70836  134   32   10.00000    1.68626  83.1%  11.3   70s
 350528 217001    3.10813  137   39   10.00000    1.71566  82.8%  11.4   76s
 363226 222997    3.06461  129   38   10.00000    1.74139  82.6%  11.5   80s
 380600 231411    2.86165  136   32   10.00000    1.77978  82.2%  11.6   85s
 398193 239288    7.22327   51   18   10.00000    1.82048  81.8%  11.6   90s
 414707 246352 infeasible  142        10.00000    1.87331  81.3%  11.7   95s
 432218 253516    4.06421  139   33   10.00000    1.93321  80.7%  11.8  100s

Explored 448041 nodes (5288084 simplex iterations) in 104.38 seconds
Thread count was 8 (of 8 available processors)

Solution count 3: 10 13 47 

Optimal solution found (tolerance 8.00e-01)
Best objective 1.000000000000e+01, best bound 3.000000000000e+00, gap 70.0000%
104.429030 seconds (169.79 k allocations: 10.398 MiB)

When I ran the line MathOptInterface.get(model, MathOptInterface.RelativeGap()) the result I got was 2.333 corresponding to (10-3)/3 instead of 0.7.

[odow]

This isn't strictly a MOI issue. Here is the offending code in Gurobi.jl:
https://github.com/JuliaOpt/Gurobi.jl/blob/3e77b5af7935ab053786dca977fc9d52dfa5d280/src/MOIWrapper.jl#L472-L476

We probably need a test in MOI for this, but it's really hard to do...

For consistency, we should probably adopt Gurobi's definition and return Inf when no bound or incumbent solution has been found, or when the current incumbent is 0.

[mlubin]

I seem to remember also that SCIP and Gurobi use different definitions of the MIP gap. We should double check before going all in on Gurobi's definition.

[jd-lara]

@odow not sure if this goes here as well, but when I tested it for GLPK to see what happens when the solver can't find the Incumbent, got an error for the method. It seems GLPK.jl doesn't have the function implemented.

MethodError: no method matching get_relative_mip_gap(::GLPK.Optimizer)
Closest candidates are:
  get_relative_mip_gap(!Matched::Gurobi.Optimizer) at /Users/jdlara/.julia/packages/Gurobi/kaIvm/src/MOIWrapper.jl:473
  get_relative_mip_gap(!Matched::LinQuadOptInterface.MockLinQuadOptimizer) at /Users/jdlara/.julia/packages/LinQuadOptInterface/0IqWS/src/mockoptimizer.jl:823

[odow]

Here is what different solvers do for reference:

• CPLEX: |bound - incumbent| / (|incumbent| + 1e-10)
• GLPK: |bound - incumbent| / (|incumbent| + EPSILON)
• Gurobi: |bound - incumbent| / |incumbent|
• SCIP: |bound - incumbent|/MIN(|bound|,|incumbent|)

[jd-lara]

Would it make more sense to get the Gap directly from the solver's API?

[mlubin]

Would it make more sense to get the Gap directly from the solver's API?

It wouldn't hurt to make this accessible by some attribute if it's clearly documented as returning the gap using the solver's definition. But doesn't this only hide the issue? What if you want to compare the performance of different solvers?

[jd-lara]

@mlubin I think using the solver's definition makes everything more transparent and that's a value added. Even if it complicates the comparison of performance, the user interested in that can get the bound and solution using MathOptInterface.get to implement the comparison directly from the solver's result.

I would also say, from a user perspective, that general approach in MOI seems to be that one has to know which solver related parameters to use when implementing specific configurations or wanting to get some behavior. Hence, getting the bound from the solver directly would be consistent with that approach, even if it creates headaches from time to time :P

[blegat]

Why don't we define a RelativeGap(εmin, εmax) defined as |bound - incumbent| / (|incumbent| + ε) where ε should be between εmin and εmax ? We could define a result fallback and if εmin and εmax do not match what the solver returns then it the wrapper use the fallback with εmax.

[joaquimg]

Here is what different solvers do for reference:

• CPLEX: |bound - incumbent| / (|incumbent| + 1e-10)
• GLPK: |bound - incumbent| / (|incumbent| + EPSILON)
• Gurobi: |bound - incumbent| / |incumbent|
• SCIP: |bound - incumbent|/MIN(|bound|,|incumbent|)

for completeness:

• Xpress |bound - incumbent|/MAX(|bound|,|incumbent|)