linprog returns 'infeasible' for simple unbounded case #11617
Description
The linprog with interior-point method can return 'infeasible' for cases where the problem is unbounded.
Reproducing code example:
If you enter the simple LP
min x1
s.t. x1 <= -1
then pre-solve correctly detects that it's unbounded:
c = [1.]
A_ub = [[1.]]
b_ub = [-1.]
scipy.optimize.linprog(c, A_ub, b_ub, None, None, bounds=(None, None), method='interior-point') con: array([], dtype=float64)
fun: -inf
message: 'The problem is (trivially) unbounded because there are no non-trivial constraints and a) at least one decision variable is unbounded above and its corresponding cost is negative, or b) at least one decision variable is unbounded below and its corresponding cost is positive. '
nit: 0
slack: array([inf])
status: 3
success: False
x: array([-inf])
But if you convert that LP to standard form
min x1 - x3
s.t. x1 + x2 - x3 = -1
x1, x2, x3 >= 0
you find that linprog (interior-point only) considers it infeasible:
c = [1., 0., -1.]
A_eq = [[1., 1., -1.]]
b_eq = [-1.]
scipy.optimize.linprog(c, None, None, A_eq, b_eq, method='interior-point') con: array([-33.89879])
fun: -1735502760.3390775
message: 'The algorithm terminated successfully and determined that the problem is infeasible.'
nit: 3
slack: array([], dtype=float64)
status: 2
success: False
x: array([1.40062e+09, 1.73550e+09, 3.13612e+09])
The bug is very sensitive to coefficients: changing to c = [.9, 0., -.9] correctly returns 'unbounded.'
Possible cause:
The cause seems to be the logic below:
scipy/scipy/optimize/_linprog_ip.py
Lines 806 to 810 in 0538284
As far as I can tell, the bug occurs because, even though -c.T @ x (0.8669123426) is a large positive value and thus kappa > 0 is due to unboundedness rather than infeasibility, we still have b.T @ y (0.0000000151) slightly greater than tol (0.0000000100), possibly due to numerical inaccuracies. The above logic falsely reports 'infeasible' in such cases rather than unbounded.
Scipy/Numpy/Python version information:
scipy=1.4.1 (or 1.5.0.dev0+9beb4c5)
numpy=1.18.1 (or 1.19.0.dev0+f1a247f)
sys.version_info(major=3, minor=8, micro=1, releaselevel='final', serial=0)