Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

MAINT: Add warnings if pivot value is close to tolerance in linprog(method='simplex') #9081

Merged
merged 6 commits into from
Jul 29, 2018
Merged

MAINT: Add warnings if pivot value is close to tolerance in linprog(method='simplex') #9081

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
dacb823
Add contribution to THANKS.txt
Kai-Striega Jul 26, 2018
ed0cb85
Refactor pivot technique into new function.
Kai-Striega Jul 27, 2018
65559ac
Add warning if pivot value is close to the tolerance.
Kai-Striega Jul 27, 2018
84a73f0
Update tests to reflect warnings.
Kai-Striega Jul 27, 2018
f095ebb
Fix comment on expected behaviour of test.
Kai-Striega Jul 27, 2018
b081e9e
Slight changes to codestyle and warning message.
Kai-Striega Jul 28, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions THANKS.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,7 @@ Konrad Griessinger for the small sample Kendall test
Tony Xiang for improvements in scipy.sparse
Roy Zywina for contributions to scipy.fftpack.
Christian H. Meyer for bug fixes in subspace_angles.
Kai Striega for improvments to the numerical stability of the simplex method.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

*improvements
And I don't think your contributions are limited to numerical stability. I would say "improvements to the scipy.optimize.linprog simplex method."


Institutions
------------
Expand Down
65 changes: 45 additions & 20 deletions scipy/optimize/_linprog.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
"""
A top-level linear programming interface. Currently this interface only
solves linear programming problems via the Simplex Method.
A top-level linear programming interface. Currently this interface solves
linear programming problems via the Simplex and Interior-Point methods.

.. versionadded:: 0.15.0

Expand All @@ -18,7 +18,8 @@
from __future__ import division, print_function, absolute_import

import numpy as np
from .optimize import OptimizeResult, _check_unknown_options
from warnings import warn
from .optimize import OptimizeResult, OptimizeWarning, _check_unknown_options
from ._linprog_ip import _linprog_ip

__all__ = ['linprog', 'linprog_verbose_callback', 'linprog_terse_callback']
Expand Down Expand Up @@ -128,7 +129,6 @@ def linprog_terse_callback(xk, **kwargs):
(and this is the final call to callback), otherwise False.
"""
nit = kwargs["nit"]

if nit == 0:
print("Iter: X:")
print("{0: <5d} ".format(nit), end="")
Expand Down Expand Up @@ -218,6 +218,45 @@ def _pivot_row(T, basis, pivcol, phase, tol=1.0E-12, bland=False):
return True, min_rows[0]


def _apply_pivot(T, basis, pivrow, pivcol, tol=1e-12):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't thought about what this does, but I have compared it to the code it replaces. I'm not seeing any difference between calling this function and writing out all those (repeated lines) except that pivval is not available after the function is done. Apparently that is not needed, so I don't think this changes anything functionally - it's better style in a few ways.

Copy link
Member Author

@Kai-Striega Kai-Striega Jul 28, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The change was basically to make sure I didn't forget to change the warning message in both places. Other than that the code is copy and pasted from the two sections.

"""
Pivot the simplex tableau inplace on the element given by (pivrow, pivol).
The entering variable corresponds to the column given by pivcol forcing
the variable basis[pivrow] to leave the basis.

Parameters
----------
T : 2-D array
A 2-D array representing the simplex T to the corresponding
maximization problem.
basis : 1-D array
An array of the indices of the basic variables, such that basis[i]
contains the column corresponding to the basic variable for row i.
Basis is modified in place by _apply_pivot.
pivrow : int
Row index of the pivot.
pivcol : int
Column index of the pivot.
"""
basis[pivrow] = pivcol
pivval = T[pivrow, pivcol]
T[pivrow] = T[pivrow] / pivval
for irow in range(T.shape[0]):
if irow != pivrow:
T[irow] = T[irow] - T[pivrow] * T[irow, pivcol]

# The selected pivot should never lead to a pivot value less than the tol.
if np.isclose(pivval, tol, atol=0, rtol=1e4):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just confirming that I checked this: this ensures that:
abs(pivval - tol) <= rtol * abs(tol)
For example, if pivval = 1e-9, tol is the default 1e-12, and rtol is the default 1e4:
abs(1e-9 - 1e-12) <= 1e4 * 1e-12
~ 1e-9 <= 1e-8
Which is True. On the other hand, if pivval were 1e-7, this would be False. This is what we want.

message = (
"The pivot operation produces a pivot value of:{0: .1e}. "
"Being only slightly greater than the set tolerance{1: .1e}. "
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"The pivot operation produces a pivot value of:{0: .1e}, "
"which is only slightly greater than the specified tolerance{1: .1e}. "

"This may lead to issues regarding the numerical stability of "
"the simplex method. Increasing the tolerance, changing the pivot "
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will increasing the tolerance actually allow some problems to be solved that won't otherwise be solved correctly?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, in theory. Having this as the first option probably wasn't the best decision and increasing the tolerance should be avoided anyway. I've change the order of the suggestions to:

  • Removing redundant constraints.
  • Using Bland's rule.
  • Increasing the tolerance

"strategy via Bland's rule or removing redundant constraints may "
"help reduce the issue.".format(pivval, tol))
warn(message, OptimizeWarning)


def _solve_simplex(T, n, basis, maxiter=1000, phase=2, callback=None,
tol=1.0E-12, nit0=0, bland=False):
"""
Expand Down Expand Up @@ -337,14 +376,7 @@ def _solve_simplex(T, n, basis, maxiter=1000, phase=2, callback=None,
if abs(T[pivrow, col]) > tol]
if len(non_zero_row) > 0:
pivcol = non_zero_row[0]
# variable represented by pivcol enters
# variable in basis[pivrow] leaves
basis[pivrow] = pivcol
pivval = T[pivrow][pivcol]
T[pivrow, :] = T[pivrow, :] / pivval
for irow in range(T.shape[0]):
if irow != pivrow:
T[irow, :] = T[irow, :] - T[pivrow, :]*T[irow, pivcol]
_apply_pivot(T, basis, pivrow, pivcol)
nit += 1

if len(basis[:m]) == 0:
Expand Down Expand Up @@ -384,14 +416,7 @@ def _solve_simplex(T, n, basis, maxiter=1000, phase=2, callback=None,
status = 1
complete = True
else:
# variable represented by pivcol enters
# variable in basis[pivrow] leaves
basis[pivrow] = pivcol
pivval = T[pivrow][pivcol]
T[pivrow, :] = T[pivrow, :] / pivval
for irow in range(T.shape[0]):
if irow != pivrow:
T[irow, :] = T[irow, :] - T[pivrow, :]*T[irow, pivcol]
_apply_pivot(T, basis, pivrow, pivcol)
nit += 1

return nit, status
Expand Down
152 changes: 122 additions & 30 deletions scipy/optimize/tests/test_linprog.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -746,6 +746,62 @@ def test_bug_8663(self):
desired_x=[0, 6./7],
desired_fun=5*6./7)

def test_bug_5400(self):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've looked through this. Looks faithful to the original post.

# https://github.com/scipy/scipy/issues/5400
bounds = [
(0, None),
(0, 100), (0, 100), (0, 100), (0, 100), (0, 100), (0, 100),
(0, 900), (0, 900), (0, 900), (0, 900), (0, 900), (0, 900),
(0, None), (0, None), (0, None), (0, None), (0, None), (0, None)]

f = 1 / 9
g = -1e4
h = -3.1
A_ub = np.array([
[1, -2.99, 0, 0, -3, 0, 0, 0, -1, -1, 0, -1, -1, 1, 1, 0, 0, 0, 0],
[1, 0, -2.9, h, 0, -3, 0, -1, 0, 0, -1, 0, -1, 0, 0, 1, 1, 0, 0],
[1, 0, 0, h, 0, 0, -3, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, 1, 1],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 1.99, -1, -1, 0, 0, 0, -1, f, f, 0, 0, 0, g, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 2, -1, -1, 0, 0, 0, -1, f, f, 0, g, 0, 0, 0, 0],
[0, -1, 1.9, 2.1, 0, 0, 0, f, -1, -1, 0, 0, 0, 0, 0, g, 0, 0, 0],
[0, 0, 0, 0, -1, 2, -1, 0, 0, 0, f, -1, f, 0, 0, 0, g, 0, 0],
[0, -1, -1, 2.1, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, 0, 0, g, 0],
[0, 0, 0, 0, -1, -1, 2, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, g]])

b_ub = np.array([
0.0, 0, 0, 100, 100, 100, 100, 100, 100, 900, 900, 900, 900, 900,
900, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

c = np.array([-1.0, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 0, 0, 0, 0, 0, 0])

if self.method == 'simplex':
with pytest.warns(OptimizeWarning):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a warning assertion according to https://docs.pytest.org/en/latest/warnings.html
Looks fine.

res = linprog(c, A_ub, b_ub, bounds=bounds,
method=self.method, options=self.options)
elif self.method == 'interior-point':
res = linprog(c, A_ub, b_ub, bounds=bounds,
method=self.method, options=self.options)
_assert_success(res, desired_fun=-106.63507541835018)


class TestLinprogSimplex(LinprogCommonTests):
method = "simplex"
Expand Down Expand Up @@ -822,6 +878,72 @@ def test_issue_6139(self):
res2, desired_fun=14.95, desired_x=np.array([5, 4.95, 5])
)

def test_issue_7237_passes_with_bland(self):
# https://github.com/scipy/scipy/issues/7237
# The simplex method sometimes "explodes" if the pivot value is very
# close to zero. Bland's rule provides an alternative pivot selection
# and produces a valid result.

c = np.array([-1., 0., 0., 0., 0., 0., 0., 0., 0.])
A_ub = np.array([
[ 1., -724., 911., -551., -555., -896., 478., -80., -293.],
[ 1., 566., 42., 937., 233., 883., 392., -909., 57.],
[ 1., -208., -894., 539., 321., 532., -924., 942., 55.],
[ 1., 857., -859., 83., 462., -265., -971., 826., 482.],
[ 1., 314., -424., 245., -424., 194., -443., -104., -429.],
[ 1., 540., 679., 361., 149., -827., 876., 633., 302.],
[ 0., -1., -0., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -1., -0., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -1., -0., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -1., -0., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -1., -0., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -1., -0., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -1., -0.],
[ 0., -0., -0., -0., -0., -0., -0., -0., -1.],
[ 0., 1., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 1., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 1., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 1., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 1., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 1., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0., 1.]
])
b_ub = np.array([
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 1., 1., 1., 1., 1., 1., 1., 1.])
A_eq = np.array([[ 0., 1., 1., 1., 1., 1., 1., 1., 1.]])
b_eq = np.array([[ 1.]])
bounds = [(None, None)] * 9

# Should warn if Bland's rule is not used.
with pytest.warns(OptimizeWarning):
res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
bounds=bounds, method=self.method, options=self.options)

o = self.options.copy()
o['bland'] = True
res_bland = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
bounds=bounds, method=self.method, options=o)
_assert_success(res_bland, desired_fun=108.568535)

def test_issue_8174_warns_if_pivval_near_tol(self):
# https://github.com/scipy/scipy/issues/8174
# The simplex method sometimes "explodes" if the pivot value is very
# close to zero.
A_ub = np.array([
[ 22714., 1008., 13380., -2713.5, -1116. ],
[ -4986., -1092., -31220., 17386.5, 684. ],
[ -4986., 0., 0., -2713.5, 0. ],
[ 22714., 0., 0., 17386.5, 0. ]])
b_ub = np.zeros(A_ub.shape[0])
c = -np.ones(A_ub.shape[1])
bounds = [(0,1)] * A_ub.shape[1]

with pytest.warns(OptimizeWarning):
linprog(c=c, A_ub=A_ub, b_ub=b_ub, bounds=bounds,
options=self.options, method=self.method)


class BaseTestLinprogIP(LinprogCommonTests):
method = "interior-point"
Expand All @@ -844,36 +966,6 @@ def test_bounds_equal_but_infeasible2(self):
method=self.method, options=self.options)
_assert_infeasible(res)

def test_bug_5400(self):
# https://github.com/scipy/scipy/issues/5400
bounds = [
(0, None),
(0, 100), (0, 100), (0, 100), (0, 100), (0, 100), (0, 100),
(0, 900), (0, 900), (0, 900), (0, 900), (0, 900), (0, 900),
(0, None), (0, None), (0, None), (0, None), (0, None), (0, None)]

f = 1 / 9
g = -1e4
h = -3.1
A_ub = np.array([
[1, -2.99, 0, 0, -3, 0, 0, 0, -1, -1, 0, -1, -1, 1, 1, 0, 0, 0, 0],
[1, 0, -2.9, h, 0, -3, 0, -1, 0, 0, -1, 0, -1, 0, 0, 1, 1, 0, 0],
[1, 0, 0, h, 0, 0, -3, -1, -1, 0, -1, -1, 0, 0, 0, 0, 0, 1, 1],
[0, 1.99, -1, -1, 0, 0, 0, -1, f, f, 0, 0, 0, g, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 2, -1, -1, 0, 0, 0, -1, f, f, 0, g, 0, 0, 0, 0],
[0, -1, 1.9, 2.1, 0, 0, 0, f, -1, -1, 0, 0, 0, 0, 0, g, 0, 0, 0],
[0, 0, 0, 0, -1, 2, -1, 0, 0, 0, f, -1, f, 0, 0, 0, g, 0, 0],
[0, -1, -1, 2.1, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, 0, 0, g, 0],
[0, 0, 0, 0, -1, -1, 2, 0, 0, 0, f, f, -1, 0, 0, 0, 0, 0, g]])

b_ub = np.array([0.0, 0, 0, 0, 0, 0, 0, 0, 0])
c = np.array([-1.0, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 0, 0, 0, 0, 0, 0])

res = linprog(c, A_ub, b_ub, bounds=bounds,
method=self.method, options=self.options)
_assert_success(res, desired_fun=-106.63507541835018)

def test_empty_constraint_1(self):
# detected in presolve?
res = linprog([-1, 1, -1, 1],
Expand Down