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.

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: Refactor linear programming solvers #9145

Merged
merged 20 commits into from
Sep 11, 2018

Conversation

Kai-Striega
Copy link
Member

Early implementation of changes discussed in #8942 and #9076. It's far from complete but I think it's worth discussing the changes so far.
Major changes:

  • Move functions related to pre and post processing into _linrpog.py
  • Remove redundant code in _linrpog_simplex(…).
  • Update tests to account for pre-solve.
  • Move simplex solver into own module

The effects on the current issues:

For 5400 and 8174 the simplex method pivots with a tolerance value very close to zero. The solution then becomes (wrongly) infeasible with a status of 0 (success). A warning is raised during the simplex method if this occurred. An error should be raised instead. Post-process catches this but returns a status of 4 and an error message instead of an error, should this be changed to a ValueError?

There are more less important changes, most of which are added as comments:

  • References / docstrings need to be updated
  • How do we handle input for presolve
  • How to test for callback functions with the interior-point
  • etc..

@mdhaber again this is far from complete, but how does this look so far?

Move ``_clean_inputs(...)``, ``_presolve(...)``, ``_get_Abc(...)`` and
``_post_process(...)`` from module ``_linprog_ip`` into ``_linprog``.

Remove duplicate behavior from ``_linprog_simplex(...)`` and modify
return values to match those required in pre and postsolve functions.
@tylerjereddy tylerjereddy added scipy.optimize maintenance Items related to regular maintenance tasks labels Aug 15, 2018
@mdhaber
Copy link
Contributor

mdhaber commented Aug 15, 2018

Excited about this! I'll try to look today.

Copy link
Contributor

@mdhaber mdhaber left a comment

Choose a reason for hiding this comment

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

Cool! Tough to follow because of how GitHub detects additions/deletions but I think it all makes sense and doesn't look too dangerous. I wrote a bunch of comments about what I think was changed and why. These are just to confirm that I checked the changed and tried to make sense of them, but if any are wrong, would you let me know?

I didn't check the test changes very carefully. Anything substantial, besides relaxing some of the tests that previously were not applied to simplex?

Only thing I saw that I was a little concerned about was poping options. I don't understand why get isn't sufficient.

tol = options.get('tol', default_tol)

# This is an undocumented option for unit testing sparse presolve
_sparse_presolve = options.pop('_sparse_presolve', False)
Copy link
Contributor

Choose a reason for hiding this comment

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

Why pop it (instead of get it)?

Copy link
Member Author

Choose a reason for hiding this comment

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

_sparse_presolve applies only to this part of the program (changing A to sparse for testing). Previously this was done in _linprog_ip(...) with _sparse_presolve passed as a (undocumented) option. The conversion now occurs in _linprog(...). _linprog_ip(...) no longer requires _sparse_presolve to be passed as an option, and has been removed. If options still contains _sparse_presolve it will be passed to _check_unknown_options(...), not be recognized as an option (as it no longer is) and raise an error.


# Solve trivial problem, eliminate variables, tighten bounds, etc...
c0 = 0 # we might get a constant term in the objective
if options.pop('presolve', True):
Copy link
Contributor

Choose a reason for hiding this comment

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

Again why are we modifying options? This is a reference to the dictionary passed in by the user, not our own local copy, right? So this would modify the user's object, which we wouldn't want to do.

Copy link
Member Author

Choose a reason for hiding this comment

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

Same reasoning as for _sparse_presolve.

I considered making a local copy of options, however the user's dictionary would not be used anywhere else. Unless you think keeping it as a reference is worth doing I don't see a benefit to keeping a duplicate copy.


# Solve trivial problem, eliminate variables, tighten bounds, etc...
c0 = 0 # we might get a constant term in the objective
if options.pop('presolve', True):
Copy link
Contributor

Choose a reason for hiding this comment

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

Again why are we modifying options here and anywhere else rather than just reading it? This is a reference to the dictionary passed in by the user, not our own local copy, right? So this would modify the user's object, which we wouldn't want to do.

@@ -1754,11 +759,9 @@ def eta(g=gamma):

def _linprog_ip(
c,
A_ub=None,
Copy link
Contributor

Choose a reason for hiding this comment

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

These deletions make sense because the problem is already in standard form.

A_eq=None,
b_eq=None,
bounds=None,
c0=0,
Copy link
Contributor

Choose a reason for hiding this comment

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

Makes sense because presolve might have added constant to objective - which is just needed in here for display, right?

Copy link
Member Author

Choose a reason for hiding this comment

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

Correct, c0 will only be used for display.

@@ -0,0 +1,496 @@
"""
Copy link
Contributor

Choose a reason for hiding this comment

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

_linprog_simplex, _solve_simplex, _apply_pivot, _pivot_row, and _pivot_col simply moved from _linprog.py

sup.filter(OptimizeWarning, "Solving system with option...")
res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
method=self.method, options=self.options)
method=self.method, options=self.options)
Copy link
Contributor

Choose a reason for hiding this comment

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

I think the indentation of this line is off? Should match with ( on previous line.

Copy link
Member Author

Choose a reason for hiding this comment

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

I'll fix it in the next commit

@@ -874,7 +885,7 @@ def cb(xk, **kwargs):
res = linprog(c, A_ub=A_ub, b_ub=b_ub, callback=cb, method=self.method)

assert_(callback_complete[0])
assert_allclose(last_xk[0], res.x)
assert_allclose(last_xk[0][:res.x.size], res.x)
Copy link
Contributor

Choose a reason for hiding this comment

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

Was there a bug in this test before?

Copy link
Member Author

Choose a reason for hiding this comment

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

The simplex method sees the artificial variables / variable substitutions from pre processing as "true" variables expecting them to appear in the solution. The variables are instead removed during _postprocess and therefore do not appear.

res = linprog(
c, A_ub, b_ub, A_eq, b_eq, method=self.method,
bounds=bounds, options={'tol': low_tol}
bounds=bounds, options=self.options
Copy link
Contributor

Choose a reason for hiding this comment

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

Not digging into this - can you explain why this changed?

if self.method == "simplex":
# Including the callback here ensures the solution can be
# calculated correctly.

Copy link
Contributor

Choose a reason for hiding this comment

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

Cool that several of these tests pass with simple now.

@mdhaber
Copy link
Contributor

mdhaber commented Aug 16, 2018

CI failure is due to a doctest in _linprog.py line 1322. Just that the number of iterations is different. Any idea why that would have changed?

AppVeyor finding some assertion errors in simplex. Not sure what's going on there.

@mdhaber
Copy link
Contributor

mdhaber commented Aug 16, 2018

Maybe also check whether test_issue_6139 and most of BaseTestLinprogIP can be moved into the common tests?

@Kai-Striega
Copy link
Member Author

Kai-Striega commented Aug 16, 2018

I think the number of iterations are different because of using presolve, at least that seems logical.

The relevant output from AppVoyeur:

    assert_(isinstance(basis, np.ndarray))
> assert_(basis.dtype == np.int_)
E  AssertionError
    basis      = array([5, 6, 7], dtype=int64)

on my system:

>>> import numpy as np

>>> a = np.array([5, 6, 7], dtype=np.int64)
>>> a.dtype == np.int_
True

So it shouldn't really be happening. Maybe it has to do with only being on python 2.7?

Copy link
Member

@ilayn ilayn left a comment

Choose a reason for hiding this comment

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

I've tried to make it at least halfway. I would really recommend checking out previous code to stay coherent with the remaining parts of the library.

And mostly for the logical comparisons we can be a bit more concise.

from warnings import warn
from .optimize import OptimizeResult, OptimizeWarning, _check_unknown_options
from scipy.optimize._remove_redundancy import _remove_redundancy
Copy link
Member

Choose a reason for hiding this comment

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

You might want to group these imports with ( on the first one and ) on the last over multiple lines.

regardless of magnitude).
c : array_like
Coefficients of the linear objective function to be minimized.
A_ub : array_like, optional
Copy link
Member

Choose a reason for hiding this comment

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

This can simply use the [...] such that A_ub @ x gives the upper bound constraints or something similar. Also - is not required for dimension notation to be consistent with the rest.

Copy link
Member Author

Choose a reason for hiding this comment

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

This was simply copied from _linporg_ip. The description has NOT been consistent between additions. Standardising the descriptions will need to be done at some point, preferably before anything new is added.

1-D array of values representing the RHS of each equality constraint
(row) in ``A_eq``.
bounds : sequence, optional
``(min, max)`` pairs for each element in ``x``, defining
Copy link
Member

Choose a reason for hiding this comment

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

Does it have to be tuple? list of tuples work? Also the explanation might use some explicit example :

[...] [(a1, b1), (a2, b2) ...] defining upper bounds on the parameter. By default all as are set to 0 and all bs to None hence modeling the interval [0, inf]

Copy link
Member Author

Choose a reason for hiding this comment

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

I agree these SHOULD be a different form. This is how they were implemented as standard @mdhaber and myself have already noted that it should be changed, but in a separate PR.


c : 1-D array
Coefficients of the linear objective function to be minimized.
A_ub : 2-D array
Copy link
Member

Choose a reason for hiding this comment

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

Same as above with dimension dashes and A_ub@x. Why do you explain the outputs if they are going to be checked anyways.

bland : bool
If True, use Bland's rule for selection of the row (if more than one
row can be used, choose the one with the lowest variable index).
try:
Copy link
Member

Choose a reason for hiding this comment

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

This looks pretty nonstandard array checking. You can already use the misc.utils functions for checking valid input arraylikes. You don't want to chase exception types but trust on np.arrayability. Those functions will complain anyways.

Copy link
Contributor

@mdhaber mdhaber Aug 16, 2018

Choose a reason for hiding this comment

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

Didn't know about those when I wrote it. It's quite thorough and works very well, so I hope we don't have to overhaul this right now. See history of such comments in #7123.

ub_newub = ubs[ub_some]
n_bounds = np.count_nonzero(ub_some)
A_ub = vstack((A_ub, zeros((n_bounds, A_ub.shape[1]))))
b_ub = np.concatenate((b_ub, np.zeros(n_bounds)))
Copy link
Member

Choose a reason for hiding this comment

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

hstack maybe ? Since you have vstack elsewhere?


# add slack variables
A2 = vstack([eye(A_ub.shape[0]), zeros((A_eq.shape[0], A_ub.shape[0]))])
A = hstack([A1, A2])
Copy link
Member

Choose a reason for hiding this comment

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

ah you have hstack afterall

where: -inf <= x[0] <= inf
Parameters
----------
x : 1-D array
Copy link
Member

Choose a reason for hiding this comment

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

Same comments as above.

1 : Iteration limit reached
2 : Problem appears to be infeasible
3 : Problem appears to be unbounded
4 : Serious numerical difficulties encountered
Copy link
Member

Choose a reason for hiding this comment

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

What is not "serious" in terms of difficulties?

Copy link
Contributor

Choose a reason for hiding this comment

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

Difficulties that have been resolved using a more robust, albeit less efficient, solver.

optimal objective value for original problem
slack: 1-D array
The (non-negative) slack in the upper bound constraints, that is,
``b_ub - A_ub * x``
Copy link
Member

Choose a reason for hiding this comment

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

A_ub @ x ?

Copy link
Contributor

Choose a reason for hiding this comment

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

You know I just learned about that a few weeks ago?

@mdhaber
Copy link
Contributor

mdhaber commented Aug 16, 2018

@ilayn Thanks for all these checks but can we wait for a separate PR to fix existing issues in the code? Most of the stuff that shows up in green is not new; it is just moved from elsewhere.

@Kai-Striega
Copy link
Member Author

@ilayn Thanks for having such a thorough look through it! I'm going to respond to most of the reviews within a single comment. I have not written most of the code. The PR moves existing code from _linporg_ip.py into _linprog.py and creates a new module for the simplex method in an attempt to remove code duplication between simplex and interior-point solvers. The existing presolve code has already been discussed elsewhere (see @mdhaber's comment on #9076) and will be addressed once the two have been merged.

@ilayn
Copy link
Member

ilayn commented Aug 16, 2018

Ah I was going through only the green lines and didn't pay too much attention to the existing code. In that case indeed let's keep it as only the refactoring.

But we at least need the stacklevel arguments while we are at it. They turn out to be quite important for CI blocking issues.

@Kai-Striega
Copy link
Member Author

Based on the discussion with @mdhaber and @ilayn, here is what I think we should be looking at for merging the solvers:

Modifying options:

The user may provide an options storing any additional solver options. My change modifies options by using pop to access presolve and _sparse_presolve.

These parameters previously only applied to the interior-point method. Now presolving occurs in the top-level linear programming module. Causing me to remove them from _linprog_ip. presolve and _sparse_presolve are now unknown options and will cause _check_unknown_options(...) to fail.

I agree changing the user provided options is not the best way to go. Further the docs define options to be:

options : dict, optional

    A dictionary of solver options. All methods accept the following generic options:

        maxiter : int

            Maximum number of iterations to perform.
        disp : bool

            Set to True to print convergence messages.

    For method-specific options, see show_options('linprog').

Which is now no longer correct. Maybe chaning presolve and _sparse_presolve to keyword arguments will provide a more accurate approach?

Tests:

  • The interior-point solver does not implement callback functions. This is currently tested for. However linprog should only raise this error if the interior-point method is actually used i.e. ``method='interior-point' and the problem is not solved in presolve.

  • Check if more tests from BaseTestLinprogIP and issue 6139 and can be moved into common tests.

Documentation

  • Edit parameter descriptions to be consistent between different functions.

  • Move docstrings describing presolve from _linprog_ip into linprog

  • Change reference numbers to be consistent between ALL linear programming modules (maybe?)

@mdhaber
Copy link
Contributor

mdhaber commented Aug 18, 2018

@Kai-Striega more to come, but why don't we just change it so that the options are not unknown? We should also add presolve to the documentation. They should be options and not keyword arguments.

@mdhaber
Copy link
Contributor

mdhaber commented Aug 18, 2018

Also the thoughts about docstrings sound great.

@Kai-Striega
Copy link
Member Author

Is there a way to add presolve as an option without adding to _linprog_simplex and _linprog_ip? I don't like passing variables to a function if they are not going to be used anywhere.

@mdhaber
Copy link
Contributor

mdhaber commented Aug 18, 2018

@Kai-Striega It doesn't bother me. If it really bothers you, you could copy the dictionary.

options = {"a":1, "b":2, "presolve":10} # example
used_options = {"rr", "presolve", "_sparse_presolve"}
remaining_options= {key:val for key, val in options.items() if key not in used_options} #unfortunately options.items() creates a list in Python 2, but I don't think it's worth it to have separate commands for Python 2/3
print(remaining_options)

Or perhaps someone else has a more clever solution. But I think it's more trouble than it's worth.

Changes to the ``simplex`` method now provide a more robust simplex
solver. Tests previously only passing for ``interior-point`` now also
pass for ``simplex``. The tests from ``BaseTestLinprogIP`` have been
moved into ``LinprogCommonTests`` to reflect this.
``test_issue_6139`` incorrectly failed for ``TestLinprogIPDense`` due
to the tolerance being set aribtrarily low. The commit increases the
``tol`` value inside ``_postprocess(...)`` by a factor of 10.
By default linprog attempts to presolve a problem. If the problem is
solved linprog avoids using a solver at all. However the test
expects a NotImplementedError to be raised regardless of if a callback
is actually required. The test now expects a NotImplementedError only
if a callback would actually have been used.
@mdhaber
Copy link
Contributor

mdhaber commented Aug 25, 2018

Hey Kai, this is pretty close? Let me know when you'd like me to take another look.

@Kai-Striega
Copy link
Member Author

Kai-Striega commented Aug 25, 2018

@mdhaber I'm getting close. I've been busy with uni giving me less time to work on this.

The doc string for linprog(...) now contains the section describing presolve with the references standardised between solver (maybe double check these?). LinprogCommonTests now tests for all tests from BaseTestLinprogIP and all pass.

I still haven't changed how linprog handles the option values and to get the doc tests passing. I've also left most the improvements by @ilayn (for now). They are not directly related to the merging, although I want to keep them in the discussion as I think they are worthwhile but should be a new PR.

Ideally I'll have time to get the above working today or tomorrow. If you would like to have a look at it those changes will be independent of the changes so far.

@rgommers
Copy link
Member

rgommers commented Sep 3, 2018

For now it would just have the sort of information as the one returned by linprog but it is flexible enough to add more information later. How about that?

that sounds good to me.

@Kai-Striega
Copy link
Member Author

Kai-Striega commented Sep 3, 2018

@mdhaber @rgommers In summary, the signature for callback should be changed to require an OptimizeResult object only. I could change this fairly quickly, however, I would prefer any further changes to be a separate PR.

EDIT: Adding only an OptimizeResult will break backwards compatibility. The current callback(...) function requires the current tableau and is used in linprog_verbose_callback(...). An OptimizeResult object does not contain the tableau.

- Add ``# may vary`` to output slack arguments
- Add unit test for docstring example
- Extend docstring for status 4 to be more descriptive
@mdhaber
Copy link
Contributor

mdhaber commented Sep 4, 2018

@Kai-Striega Please review the comments. As you can see, we were reluctant to ask for this change, as we know this PR was intended only to refactor code. However, the PR as it stood broke backwards compatibility silently by changing the meaning of the information passed to the callback function without changing the signature. Unfortunatley, this was considered unacceptable. We iterated on how the situation should be handled, and decided that the best course of action - easiest for now and sufficiently flexible for the future - would be to pass an OptimizeResult into the callback function to ensure that callback functions relying on the old callback signature would raise an error. This OptimizeResult should contain, at a minimum, an array x representing the current solution to the original problem. Thoughts on obtaining that (eliminating any phase one artificial variables and calling _postprocess) are above.

You are welcome to modify linprog_verbose_callback as necessary to avoid reliance on the tableau. We might decide in the future to provide the tableau to the current problem - whenever we are ready to also provide information about how the variables of the current problem correspond with those in the original problem. We can always add these things to an OptimizeResult later.

Thanks!

Recent changes to ``linprog`` silently break the (undefined) behaviour
of ``callback``. The commit changes the callback signature to
require an ``OptimizeResult`` rather than the existing ``xk``,
``kwargs`` by refactoring ``_postprocess`` into multiple, smaller,
functions in a new module ``_linprog_util``.
Copy link
Contributor

@mdhaber mdhaber left a comment

Choose a reason for hiding this comment

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

@Kai-Striega Looks great! I don't see very much I'd change. I really like the callback test's way of verifying that the callback function is receiving the information it is supposed to. It is good that problem chosen for the callback test has inequality constraints because then it is known to have variables added in conversion to standard form, so checking, for instance, that the values of x inside the callback are equal to those reported at completion is meaningful.

I might have checked all the fields of the OptimizeResult and chosen a problem that would also have been modified in presolve, and I think that the name of _T_o/T_o is descriptive enough.

The only thing I have to ask for is that the documentation of the callback function argument for linprog and linprog(method='simplex') be modified. They still specify the old signature and reference the tableau.

@mdhaber
Copy link
Contributor

mdhaber commented Sep 9, 2018

@rgommers See anything else that needs to be changed, or can this be merged once the documentation is updated? Regarding the documentation, can you suggest how to call attention to the compatibility break (in 1.2, presumably) and the rationale behind it?

@pv
Copy link
Member

pv commented Sep 9, 2018

Notices of backward-compatibility issues should be added to the release notes.
The draft release notes can be edited here: https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.2.0 --- the notice should be added there before release.

@mdhaber
Copy link
Contributor

mdhaber commented Sep 10, 2018

Draft release notes edited to reflect this change. @Kai-Striega You're welcome to edit if they could be more clear; just wanted to help.

@Kai-Striega
Copy link
Member Author

Kai-Striega commented Sep 10, 2018

@mdhaber I've edited the docstrings as requested, plus added more tests for fun returned by the callback. Your release notes cover everything I can think of and I appreciate the help =).

I've done some further refactoring by moving the private functions from _linprog into the linprog_util module. _linprog's module doc string states A top-level linear programming interface. I feel that if we already implement a utility module including the private functions seems appropriate.

Setting the solver_options, _clean_inputs(...) and sparsity checking are now wrapped in _parse_linprog. As mentioned by @ilayn this should help make the internals more concrete.

@mdhaber
Copy link
Contributor

mdhaber commented Sep 10, 2018

@Kai-Striega Took a quick look; it appears that you've updated the docstrings and moved things around in a way that makes sense. Looking forward to the checks passing! If so, are you done?

@Kai-Striega
Copy link
Member Author

Yes. This should be it. I think the PR is getting too unwieldly anyway so any extra changes should really be in a new PR

@mdhaber
Copy link
Contributor

mdhaber commented Sep 10, 2018

@rgommers Think this is ready? Do the draft release notes adequately describe the change?

@rgommers
Copy link
Member

Release note update looks good to me. Go ahead and merge if code is good to go.

@mdhaber mdhaber merged commit f3e0dcf into scipy:master Sep 11, 2018
@rgommers rgommers added this to the 1.2.0 milestone Sep 11, 2018
@rgommers
Copy link
Member

Awesome!

Could we close the issues that are now solved? Not quite sure if there's anything left to do in, e.g., gh-8942 or gh-9076.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
maintenance Items related to regular maintenance tasks scipy.optimize
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants