diff --git a/THANKS.txt b/THANKS.txt index defef3374499..3c1895223336 100644 --- a/THANKS.txt +++ b/THANKS.txt @@ -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 improvements to the scipy.optimize.linprog simplex method. Institutions ------------ diff --git a/scipy/optimize/_linprog.py b/scipy/optimize/_linprog.py index 8d639f47c5db..ce4ed9e8135f 100644 --- a/scipy/optimize/_linprog.py +++ b/scipy/optimize/_linprog.py @@ -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 @@ -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'] @@ -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="") @@ -218,6 +218,46 @@ 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): + """ + 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): + message = ( + "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. " + "Removing redundant constraints, changing the pivot strategy " + "via Bland's rule or increasing the tolerance 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): """ @@ -337,14 +377,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: @@ -384,14 +417,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 diff --git a/scipy/optimize/tests/test_linprog.py b/scipy/optimize/tests/test_linprog.py index ebc38eb08d19..0cffc1ffd765 100644 --- a/scipy/optimize/tests/test_linprog.py +++ b/scipy/optimize/tests/test_linprog.py @@ -746,6 +746,62 @@ def test_bug_8663(self): desired_x=[0, 6./7], desired_fun=5*6./7) + 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, 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): + 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" @@ -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" @@ -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],