-
-
Notifications
You must be signed in to change notification settings - Fork 5.1k
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.
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
MAINT: Add warnings if pivot value is close to tolerance in linprog(method='simplex') #9081
Conversation
The pivoting procedure may produce pivot values very close to zero. Dividing by these values can cause the simplex method "blow up". Adding the warning notifes the user if this is likely to occur and provides some hints to prevent this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@Kai-Striega I don't see anything wrong. Is there something specific you want me to double-check?
There are some PEP8 issues to fix, though.
THANKS.txt
Outdated
@@ -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. |
There was a problem hiding this comment.
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."
scipy/optimize/_linprog.py
Outdated
if np.isclose(pivval, tol, atol=0, rtol=1e4): | ||
message = ( | ||
"The pivot operation produces a pivot value of:{0: .1e}. " | ||
"Being only slightly greater than the set tolerance{1: .1e}. " |
There was a problem hiding this comment.
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}. "
scipy/optimize/_linprog.py
Outdated
"The pivot operation produces a pivot value of:{0: .1e}. " | ||
"Being only slightly greater than the set tolerance{1: .1e}. " | ||
"This may lead to issues regarding the numerical stability of " | ||
"the simplex method. Increasing the tolerance, changing the pivot " |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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
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): |
There was a problem hiding this comment.
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.
@@ -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): |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
@@ -746,6 +746,62 @@ def test_bug_8663(self): | |||
desired_x=[0, 6./7], | |||
desired_fun=5*6./7) | |||
|
|||
def test_bug_5400(self): |
There was a problem hiding this comment.
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.
1, 1, 1, 1, 0, 0, 0, 0, 0, 0]) | ||
|
||
if self.method == 'simplex': | ||
with pytest.warns(OptimizeWarning): |
There was a problem hiding this comment.
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.
Removes additional whitespace from tests in test_linprog. Changes the order of suggestions in the previously added warning message. The original message firstly suggested trying to increase the tolerance. This may work in theory, but shouldn't be the first choice.
Yeah if there was nothing specific you wanted me to double check. You mentioned something felt uneasy but it looked good and straightforward to me. |
@mdhaber Looking over it again and it looks good to me, I think I simply spent too long over analysing it. |
Implementation of the first steps discussed in #9076.
I've implemented the reporting as warnings instead of a status/message, in theory (or depending on how the tolerance / problem is formed) it may still be on the tableau and maintain feasibility. The definition of "close" to tolerance has been set to be within 4 orders of magnitude (rtol of 1e4) this is really a guess. It catches #5400, #7237, #8174 comfortably and could possibly be fine tuned.
I'm not sure if the implementation of my tests are as well structured as they could be. All fail before and pass after the changes but just feel a bit off. #7237 returns the correct solution when using Bland's rule and has been added as a test.