Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_linprog_rs.py: 9%
194 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Revised simplex method for linear programming
3The *revised simplex* method uses the method described in [1]_, except
4that a factorization [2]_ of the basis matrix, rather than its inverse,
5is efficiently maintained and used to solve the linear systems at each
6iteration of the algorithm.
8.. versionadded:: 1.3.0
10References
11----------
12.. [1] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear
13 programming." Athena Scientific 1 (1997): 997.
14.. [2] Bartels, Richard H. "A stabilization of the simplex method."
15 Journal in Numerische Mathematik 16.5 (1971): 414-434.
17"""
18# Author: Matt Haberland
20import numpy as np
21from numpy.linalg import LinAlgError
23from scipy.linalg import solve
24from ._optimize import _check_unknown_options
25from ._bglu_dense import LU
26from ._bglu_dense import BGLU as BGLU
27from ._linprog_util import _postsolve
28from ._optimize import OptimizeResult
31def _phase_one(A, b, x0, callback, postsolve_args, maxiter, tol, disp,
32 maxupdate, mast, pivot):
33 """
34 The purpose of phase one is to find an initial basic feasible solution
35 (BFS) to the original problem.
37 Generates an auxiliary problem with a trivial BFS and an objective that
38 minimizes infeasibility of the original problem. Solves the auxiliary
39 problem using the main simplex routine (phase two). This either yields
40 a BFS to the original problem or determines that the original problem is
41 infeasible. If feasible, phase one detects redundant rows in the original
42 constraint matrix and removes them, then chooses additional indices as
43 necessary to complete a basis/BFS for the original problem.
44 """
46 m, n = A.shape
47 status = 0
49 # generate auxiliary problem to get initial BFS
50 A, b, c, basis, x, status = _generate_auxiliary_problem(A, b, x0, tol)
52 if status == 6:
53 residual = c.dot(x)
54 iter_k = 0
55 return x, basis, A, b, residual, status, iter_k
57 # solve auxiliary problem
58 phase_one_n = n
59 iter_k = 0
60 x, basis, status, iter_k = _phase_two(c, A, x, basis, callback,
61 postsolve_args,
62 maxiter, tol, disp,
63 maxupdate, mast, pivot,
64 iter_k, phase_one_n)
66 # check for infeasibility
67 residual = c.dot(x)
68 if status == 0 and residual > tol:
69 status = 2
71 # drive artificial variables out of basis
72 # TODO: test redundant row removal better
73 # TODO: make solve more efficient with BGLU? This could take a while.
74 keep_rows = np.ones(m, dtype=bool)
75 for basis_column in basis[basis >= n]:
76 B = A[:, basis]
77 try:
78 basis_finder = np.abs(solve(B, A)) # inefficient
79 pertinent_row = np.argmax(basis_finder[:, basis_column])
80 eligible_columns = np.ones(n, dtype=bool)
81 eligible_columns[basis[basis < n]] = 0
82 eligible_column_indices = np.where(eligible_columns)[0]
83 index = np.argmax(basis_finder[:, :n]
84 [pertinent_row, eligible_columns])
85 new_basis_column = eligible_column_indices[index]
86 if basis_finder[pertinent_row, new_basis_column] < tol:
87 keep_rows[pertinent_row] = False
88 else:
89 basis[basis == basis_column] = new_basis_column
90 except LinAlgError:
91 status = 4
93 # form solution to original problem
94 A = A[keep_rows, :n]
95 basis = basis[keep_rows]
96 x = x[:n]
97 m = A.shape[0]
98 return x, basis, A, b, residual, status, iter_k
101def _get_more_basis_columns(A, basis):
102 """
103 Called when the auxiliary problem terminates with artificial columns in
104 the basis, which must be removed and replaced with non-artificial
105 columns. Finds additional columns that do not make the matrix singular.
106 """
107 m, n = A.shape
109 # options for inclusion are those that aren't already in the basis
110 a = np.arange(m+n)
111 bl = np.zeros(len(a), dtype=bool)
112 bl[basis] = 1
113 options = a[~bl]
114 options = options[options < n] # and they have to be non-artificial
116 # form basis matrix
117 B = np.zeros((m, m))
118 B[:, 0:len(basis)] = A[:, basis]
120 if (basis.size > 0 and
121 np.linalg.matrix_rank(B[:, :len(basis)]) < len(basis)):
122 raise Exception("Basis has dependent columns")
124 rank = 0 # just enter the loop
125 for i in range(n): # somewhat arbitrary, but we need another way out
126 # permute the options, and take as many as needed
127 new_basis = np.random.permutation(options)[:m-len(basis)]
128 B[:, len(basis):] = A[:, new_basis] # update the basis matrix
129 rank = np.linalg.matrix_rank(B) # check the rank
130 if rank == m:
131 break
133 return np.concatenate((basis, new_basis))
136def _generate_auxiliary_problem(A, b, x0, tol):
137 """
138 Modifies original problem to create an auxiliary problem with a trivial
139 initial basic feasible solution and an objective that minimizes
140 infeasibility in the original problem.
142 Conceptually, this is done by stacking an identity matrix on the right of
143 the original constraint matrix, adding artificial variables to correspond
144 with each of these new columns, and generating a cost vector that is all
145 zeros except for ones corresponding with each of the new variables.
147 A initial basic feasible solution is trivial: all variables are zero
148 except for the artificial variables, which are set equal to the
149 corresponding element of the right hand side `b`.
151 Runnning the simplex method on this auxiliary problem drives all of the
152 artificial variables - and thus the cost - to zero if the original problem
153 is feasible. The original problem is declared infeasible otherwise.
155 Much of the complexity below is to improve efficiency by using singleton
156 columns in the original problem where possible, thus generating artificial
157 variables only as necessary, and using an initial 'guess' basic feasible
158 solution.
159 """
160 status = 0
161 m, n = A.shape
163 if x0 is not None:
164 x = x0
165 else:
166 x = np.zeros(n)
168 r = b - A@x # residual; this must be all zeros for feasibility
170 A[r < 0] = -A[r < 0] # express problem with RHS positive for trivial BFS
171 b[r < 0] = -b[r < 0] # to the auxiliary problem
172 r[r < 0] *= -1
174 # Rows which we will need to find a trivial way to zero.
175 # This should just be the rows where there is a nonzero residual.
176 # But then we would not necessarily have a column singleton in every row.
177 # This makes it difficult to find an initial basis.
178 if x0 is None:
179 nonzero_constraints = np.arange(m)
180 else:
181 nonzero_constraints = np.where(r > tol)[0]
183 # these are (at least some of) the initial basis columns
184 basis = np.where(np.abs(x) > tol)[0]
186 if len(nonzero_constraints) == 0 and len(basis) <= m: # already a BFS
187 c = np.zeros(n)
188 basis = _get_more_basis_columns(A, basis)
189 return A, b, c, basis, x, status
190 elif (len(nonzero_constraints) > m - len(basis) or
191 np.any(x < 0)): # can't get trivial BFS
192 c = np.zeros(n)
193 status = 6
194 return A, b, c, basis, x, status
196 # chooses existing columns appropriate for inclusion in initial basis
197 cols, rows = _select_singleton_columns(A, r)
199 # find the rows we need to zero that we _can_ zero with column singletons
200 i_tofix = np.isin(rows, nonzero_constraints)
201 # these columns can't already be in the basis, though
202 # we are going to add them to the basis and change the corresponding x val
203 i_notinbasis = np.logical_not(np.isin(cols, basis))
204 i_fix_without_aux = np.logical_and(i_tofix, i_notinbasis)
205 rows = rows[i_fix_without_aux]
206 cols = cols[i_fix_without_aux]
208 # indices of the rows we can only zero with auxiliary variable
209 # these rows will get a one in each auxiliary column
210 arows = nonzero_constraints[np.logical_not(
211 np.isin(nonzero_constraints, rows))]
212 n_aux = len(arows)
213 acols = n + np.arange(n_aux) # indices of auxiliary columns
215 basis_ng = np.concatenate((cols, acols)) # basis columns not from guess
216 basis_ng_rows = np.concatenate((rows, arows)) # rows we need to zero
218 # add auxiliary singleton columns
219 A = np.hstack((A, np.zeros((m, n_aux))))
220 A[arows, acols] = 1
222 # generate initial BFS
223 x = np.concatenate((x, np.zeros(n_aux)))
224 x[basis_ng] = r[basis_ng_rows]/A[basis_ng_rows, basis_ng]
226 # generate costs to minimize infeasibility
227 c = np.zeros(n_aux + n)
228 c[acols] = 1
230 # basis columns correspond with nonzeros in guess, those with column
231 # singletons we used to zero remaining constraints, and any additional
232 # columns to get a full set (m columns)
233 basis = np.concatenate((basis, basis_ng))
234 basis = _get_more_basis_columns(A, basis) # add columns as needed
236 return A, b, c, basis, x, status
239def _select_singleton_columns(A, b):
240 """
241 Finds singleton columns for which the singleton entry is of the same sign
242 as the right-hand side; these columns are eligible for inclusion in an
243 initial basis. Determines the rows in which the singleton entries are
244 located. For each of these rows, returns the indices of the one singleton
245 column and its corresponding row.
246 """
247 # find indices of all singleton columns and corresponding row indicies
248 column_indices = np.nonzero(np.sum(np.abs(A) != 0, axis=0) == 1)[0]
249 columns = A[:, column_indices] # array of singleton columns
250 row_indices = np.zeros(len(column_indices), dtype=int)
251 nonzero_rows, nonzero_columns = np.nonzero(columns)
252 row_indices[nonzero_columns] = nonzero_rows # corresponding row indicies
254 # keep only singletons with entries that have same sign as RHS
255 # this is necessary because all elements of BFS must be non-negative
256 same_sign = A[row_indices, column_indices]*b[row_indices] >= 0
257 column_indices = column_indices[same_sign][::-1]
258 row_indices = row_indices[same_sign][::-1]
259 # Reversing the order so that steps below select rightmost columns
260 # for initial basis, which will tend to be slack variables. (If the
261 # guess corresponds with a basic feasible solution but a constraint
262 # is not satisfied with the corresponding slack variable zero, the slack
263 # variable must be basic.)
265 # for each row, keep rightmost singleton column with an entry in that row
266 unique_row_indices, first_columns = np.unique(row_indices,
267 return_index=True)
268 return column_indices[first_columns], unique_row_indices
271def _find_nonzero_rows(A, tol):
272 """
273 Returns logical array indicating the locations of rows with at least
274 one nonzero element.
275 """
276 return np.any(np.abs(A) > tol, axis=1)
279def _select_enter_pivot(c_hat, bl, a, rule="bland", tol=1e-12):
280 """
281 Selects a pivot to enter the basis. Currently Bland's rule - the smallest
282 index that has a negative reduced cost - is the default.
283 """
284 if rule.lower() == "mrc": # index with minimum reduced cost
285 return a[~bl][np.argmin(c_hat)]
286 else: # smallest index w/ negative reduced cost
287 return a[~bl][c_hat < -tol][0]
290def _display_iter(phase, iteration, slack, con, fun):
291 """
292 Print indicators of optimization status to the console.
293 """
294 header = True if not iteration % 20 else False
296 if header:
297 print("Phase",
298 "Iteration",
299 "Minimum Slack ",
300 "Constraint Residual",
301 "Objective ")
303 # :<X.Y left aligns Y digits in X digit spaces
304 fmt = '{0:<6}{1:<10}{2:<20.13}{3:<20.13}{4:<20.13}'
305 try:
306 slack = np.min(slack)
307 except ValueError:
308 slack = "NA"
309 print(fmt.format(phase, iteration, slack, np.linalg.norm(con), fun))
312def _display_and_callback(phase_one_n, x, postsolve_args, status,
313 iteration, disp, callback):
314 if phase_one_n is not None:
315 phase = 1
316 x_postsolve = x[:phase_one_n]
317 else:
318 phase = 2
319 x_postsolve = x
320 x_o, fun, slack, con = _postsolve(x_postsolve,
321 postsolve_args)
323 if callback is not None:
324 res = OptimizeResult({'x': x_o, 'fun': fun, 'slack': slack,
325 'con': con, 'nit': iteration,
326 'phase': phase, 'complete': False,
327 'status': status, 'message': "",
328 'success': False})
329 callback(res)
330 if disp:
331 _display_iter(phase, iteration, slack, con, fun)
334def _phase_two(c, A, x, b, callback, postsolve_args, maxiter, tol, disp,
335 maxupdate, mast, pivot, iteration=0, phase_one_n=None):
336 """
337 The heart of the simplex method. Beginning with a basic feasible solution,
338 moves to adjacent basic feasible solutions successively lower reduced cost.
339 Terminates when there are no basic feasible solutions with lower reduced
340 cost or if the problem is determined to be unbounded.
342 This implementation follows the revised simplex method based on LU
343 decomposition. Rather than maintaining a tableau or an inverse of the
344 basis matrix, we keep a factorization of the basis matrix that allows
345 efficient solution of linear systems while avoiding stability issues
346 associated with inverted matrices.
347 """
348 m, n = A.shape
349 status = 0
350 a = np.arange(n) # indices of columns of A
351 ab = np.arange(m) # indices of columns of B
352 if maxupdate:
353 # basis matrix factorization object; similar to B = A[:, b]
354 B = BGLU(A, b, maxupdate, mast)
355 else:
356 B = LU(A, b)
358 for iteration in range(iteration, maxiter):
360 if disp or callback is not None:
361 _display_and_callback(phase_one_n, x, postsolve_args, status,
362 iteration, disp, callback)
364 bl = np.zeros(len(a), dtype=bool)
365 bl[b] = 1
367 xb = x[b] # basic variables
368 cb = c[b] # basic costs
370 try:
371 v = B.solve(cb, transposed=True) # similar to v = solve(B.T, cb)
372 except LinAlgError:
373 status = 4
374 break
376 # TODO: cythonize?
377 c_hat = c - v.dot(A) # reduced cost
378 c_hat = c_hat[~bl]
379 # Above is much faster than:
380 # N = A[:, ~bl] # slow!
381 # c_hat = c[~bl] - v.T.dot(N)
382 # Can we perform the multiplication only on the nonbasic columns?
384 if np.all(c_hat >= -tol): # all reduced costs positive -> terminate
385 break
387 j = _select_enter_pivot(c_hat, bl, a, rule=pivot, tol=tol)
388 u = B.solve(A[:, j]) # similar to u = solve(B, A[:, j])
390 i = u > tol # if none of the u are positive, unbounded
391 if not np.any(i):
392 status = 3
393 break
395 th = xb[i]/u[i]
396 l = np.argmin(th) # implicitly selects smallest subscript
397 th_star = th[l] # step size
399 x[b] = x[b] - th_star*u # take step
400 x[j] = th_star
401 B.update(ab[i][l], j) # modify basis
402 b = B.b # similar to b[ab[i][l]] =
404 else:
405 # If the end of the for loop is reached (without a break statement),
406 # then another step has been taken, so the iteration counter should
407 # increment, info should be displayed, and callback should be called.
408 iteration += 1
409 status = 1
410 if disp or callback is not None:
411 _display_and_callback(phase_one_n, x, postsolve_args, status,
412 iteration, disp, callback)
414 return x, b, status, iteration
417def _linprog_rs(c, c0, A, b, x0, callback, postsolve_args,
418 maxiter=5000, tol=1e-12, disp=False,
419 maxupdate=10, mast=False, pivot="mrc",
420 **unknown_options):
421 """
422 Solve the following linear programming problem via a two-phase
423 revised simplex algorithm.::
425 minimize: c @ x
427 subject to: A @ x == b
428 0 <= x < oo
430 User-facing documentation is in _linprog_doc.py.
432 Parameters
433 ----------
434 c : 1-D array
435 Coefficients of the linear objective function to be minimized.
436 c0 : float
437 Constant term in objective function due to fixed (and eliminated)
438 variables. (Currently unused.)
439 A : 2-D array
440 2-D array which, when matrix-multiplied by ``x``, gives the values of
441 the equality constraints at ``x``.
442 b : 1-D array
443 1-D array of values representing the RHS of each equality constraint
444 (row) in ``A_eq``.
445 x0 : 1-D array, optional
446 Starting values of the independent variables, which will be refined by
447 the optimization algorithm. For the revised simplex method, these must
448 correspond with a basic feasible solution.
449 callback : callable, optional
450 If a callback function is provided, it will be called within each
451 iteration of the algorithm. The callback function must accept a single
452 `scipy.optimize.OptimizeResult` consisting of the following fields:
454 x : 1-D array
455 Current solution vector.
456 fun : float
457 Current value of the objective function ``c @ x``.
458 success : bool
459 True only when an algorithm has completed successfully,
460 so this is always False as the callback function is called
461 only while the algorithm is still iterating.
462 slack : 1-D array
463 The values of the slack variables. Each slack variable
464 corresponds to an inequality constraint. If the slack is zero,
465 the corresponding constraint is active.
466 con : 1-D array
467 The (nominally zero) residuals of the equality constraints,
468 that is, ``b - A_eq @ x``.
469 phase : int
470 The phase of the algorithm being executed.
471 status : int
472 For revised simplex, this is always 0 because if a different
473 status is detected, the algorithm terminates.
474 nit : int
475 The number of iterations performed.
476 message : str
477 A string descriptor of the exit status of the optimization.
478 postsolve_args : tuple
479 Data needed by _postsolve to convert the solution to the standard-form
480 problem into the solution to the original problem.
482 Options
483 -------
484 maxiter : int
485 The maximum number of iterations to perform in either phase.
486 tol : float
487 The tolerance which determines when a solution is "close enough" to
488 zero in Phase 1 to be considered a basic feasible solution or close
489 enough to positive to serve as an optimal solution.
490 disp : bool
491 Set to ``True`` if indicators of optimization status are to be printed
492 to the console each iteration.
493 maxupdate : int
494 The maximum number of updates performed on the LU factorization.
495 After this many updates is reached, the basis matrix is factorized
496 from scratch.
497 mast : bool
498 Minimize Amortized Solve Time. If enabled, the average time to solve
499 a linear system using the basis factorization is measured. Typically,
500 the average solve time will decrease with each successive solve after
501 initial factorization, as factorization takes much more time than the
502 solve operation (and updates). Eventually, however, the updated
503 factorization becomes sufficiently complex that the average solve time
504 begins to increase. When this is detected, the basis is refactorized
505 from scratch. Enable this option to maximize speed at the risk of
506 nondeterministic behavior. Ignored if ``maxupdate`` is 0.
507 pivot : "mrc" or "bland"
508 Pivot rule: Minimum Reduced Cost (default) or Bland's rule. Choose
509 Bland's rule if iteration limit is reached and cycling is suspected.
510 unknown_options : dict
511 Optional arguments not used by this particular solver. If
512 `unknown_options` is non-empty a warning is issued listing all
513 unused options.
515 Returns
516 -------
517 x : 1-D array
518 Solution vector.
519 status : int
520 An integer representing the exit status of the optimization::
522 0 : Optimization terminated successfully
523 1 : Iteration limit reached
524 2 : Problem appears to be infeasible
525 3 : Problem appears to be unbounded
526 4 : Numerical difficulties encountered
527 5 : No constraints; turn presolve on
528 6 : Guess x0 cannot be converted to a basic feasible solution
530 message : str
531 A string descriptor of the exit status of the optimization.
532 iteration : int
533 The number of iterations taken to solve the problem.
534 """
536 _check_unknown_options(unknown_options)
538 messages = ["Optimization terminated successfully.",
539 "Iteration limit reached.",
540 "The problem appears infeasible, as the phase one auxiliary "
541 "problem terminated successfully with a residual of {0:.1e}, "
542 "greater than the tolerance {1} required for the solution to "
543 "be considered feasible. Consider increasing the tolerance to "
544 "be greater than {0:.1e}. If this tolerance is unnaceptably "
545 "large, the problem is likely infeasible.",
546 "The problem is unbounded, as the simplex algorithm found "
547 "a basic feasible solution from which there is a direction "
548 "with negative reduced cost in which all decision variables "
549 "increase.",
550 "Numerical difficulties encountered; consider trying "
551 "method='interior-point'.",
552 "Problems with no constraints are trivially solved; please "
553 "turn presolve on.",
554 "The guess x0 cannot be converted to a basic feasible "
555 "solution. "
556 ]
558 if A.size == 0: # address test_unbounded_below_no_presolve_corrected
559 return np.zeros(c.shape), 5, messages[5], 0
561 x, basis, A, b, residual, status, iteration = (
562 _phase_one(A, b, x0, callback, postsolve_args,
563 maxiter, tol, disp, maxupdate, mast, pivot))
565 if status == 0:
566 x, basis, status, iteration = _phase_two(c, A, x, basis, callback,
567 postsolve_args,
568 maxiter, tol, disp,
569 maxupdate, mast, pivot,
570 iteration)
572 return x, status, messages[status].format(residual, tol), iteration