Presolver incorrectly reports infeasibility in awkwardly scaled problem

[fuglede]

Here's one that came up during today's Advent of Code puzzle.

Take the following ILP:

Minimize
  3x + y
Subject to
  47x + 19y = 10000000002226
  23x + 57y = 10000000013254
General
  x y
End

If solved with --presolve off, HiGHS finds the correct optimal solution. With --presolve on it incorrectly reports infeasibility, although it does emit a warning suggesting to scale the bounds so they aren't quite as big. Doing so, however, is not sufficient: The presolver concludes that the one below is infeasible, even though it can be solved without the presolver.

Minimize
  3x + y
Subject to
  0.0000000047x + 0.0000000019y = 1000.0000002226
  0.0000000023x + 0.0000000057y = 1000.0000013254
General
  x y
End

[jajhall]

Although pathological, this is an instance that's easy to explore, so @jajhall and @fwesselm may learn something from it.

[fwesselm]

The problem is found to be infeasible in HPresolve::sparsify. I will have a look.

[fwesselm]

Sparsification adds the first row to the second row with a multiplier of -23/47, thus creating the row singleton:

(2242/47) * y = 10000000013254 - (23/47) * 10000000002226 = 5106382990888.0850

Singleton row reduction then deduces that

y = 5106382990888.0850 / (2242/47)

but due to floating-point roundoff we actually get

y = 107047279469.99998

(and not the needed integer 107047279470). Thus, presolve (incorrectly) deduces that the problem is infeasible.

I am wondering if rows with huge right-hand/left-hand sides should be skipped in HPresolve::sparsify.

[jajhall]

I see the issue, thanks - I've just worked through this with VScode.

• The equal bounds should be an integer, so that they are shifted by boundTol either side of that value, and ceil/floor yields a fixed column.
• Rather, the equal bounds are near-integer but fractional, and their magnitude is such that the action of boundTol doesn't change them. Hence ceil/floor yields inconsistent bounds on the column.

Working backwards, there will be round-off in computing y (in general), but the round-off needs to be a couple of orders of magnitude smaller than boundTol so that, when it's used to shift y, ceil/floor yields a fixed column.

So, for variable coefficient $v$, RHS value $r$, and primal feasibility tolerance $\epsilon$, we need

• For $|v|\ge1$, $(r/v)10^{-16} < 10^{-2} (\epsilon/v)$ so $R<10^{14}\epsilon$, which is $R<10^8$ for default $\epsilon=10^{-6}$
• For $|v|<1$, $(r/v)10^{-16} < 10^{-2} (\epsilon)$ so $R<10^{14}\epsilon v$, which is $R<10^8 v$ for default $\epsilon=10^{-6}$

so the latter case will not apply the sparsification according to the RHS value and variable coefficient $v$. But there is a limit on $v$ being small for sparsification already.

[jajhall]

Instance added as (failing) unit test in https://github.com/ERGO-Code/HiGHS/tree/fix-2084