Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_linprog_ip.py: 11%

1"""Interior-point method for linear programming 

2 

3The *interior-point* method uses the primal-dual path following algorithm 

4outlined in [1]_. This algorithm supports sparse constraint matrices and 

5is typically faster than the simplex methods, especially for large, sparse 

6problems. Note, however, that the solution returned may be slightly less 

7accurate than those of the simplex methods and will not, in general, 

8correspond with a vertex of the polytope defined by the constraints. 

9 

10 .. versionadded:: 1.0.0 

11 

12References 

13---------- 

14.. [1] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

15 optimizer for linear programming: an implementation of the 

16 homogeneous algorithm." High performance optimization. Springer US, 

17 2000. 197-232. 

18""" 

19# Author: Matt Haberland 

20 

21import numpy as np 

22import scipy as sp 

23import scipy.sparse as sps 

24from warnings import warn 

25from scipy.linalg import LinAlgError 

26from ._optimize import OptimizeWarning, OptimizeResult, _check_unknown_options 

27from ._linprog_util import _postsolve 

28has_umfpack = True 

29has_cholmod = True 

30try: 

31 import sksparse 

32 from sksparse.cholmod import cholesky as cholmod 

33 from sksparse.cholmod import analyze as cholmod_analyze 

34except ImportError: 

35 has_cholmod = False 

36try: 

37 import scikits.umfpack # test whether to use factorized 

38except ImportError: 

39 has_umfpack = False 

40 

41 

42def _get_solver(M, sparse=False, lstsq=False, sym_pos=True, 

43 cholesky=True, permc_spec='MMD_AT_PLUS_A'): 

44 """ 

45 Given solver options, return a handle to the appropriate linear system 

46 solver. 

47 

48 Parameters 

49 ---------- 

50 M : 2-D array 

51 As defined in [4] Equation 8.31 

52 sparse : bool (default = False) 

53 True if the system to be solved is sparse. This is typically set 

54 True when the original ``A_ub`` and ``A_eq`` arrays are sparse. 

55 lstsq : bool (default = False) 

56 True if the system is ill-conditioned and/or (nearly) singular and 

57 thus a more robust least-squares solver is desired. This is sometimes 

58 needed as the solution is approached. 

59 sym_pos : bool (default = True) 

60 True if the system matrix is symmetric positive definite 

61 Sometimes this needs to be set false as the solution is approached, 

62 even when the system should be symmetric positive definite, due to 

63 numerical difficulties. 

64 cholesky : bool (default = True) 

65 True if the system is to be solved by Cholesky, rather than LU, 

66 decomposition. This is typically faster unless the problem is very 

67 small or prone to numerical difficulties. 

68 permc_spec : str (default = 'MMD_AT_PLUS_A') 

69 Sparsity preservation strategy used by SuperLU. Acceptable values are: 

70 

71 - ``NATURAL``: natural ordering. 

72 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A. 

73 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A. 

74 - ``COLAMD``: approximate minimum degree column ordering. 

75 

76 See SuperLU documentation. 

77 

78 Returns 

79 ------- 

80 solve : function 

81 Handle to the appropriate solver function 

82 

83 """ 

84 try: 

85 if sparse: 

86 if lstsq: 

87 def solve(r, sym_pos=False): 

88 return sps.linalg.lsqr(M, r)[0] 

89 elif cholesky: 

90 try: 

91 # Will raise an exception in the first call, 

92 # or when the matrix changes due to a new problem 

93 _get_solver.cholmod_factor.cholesky_inplace(M) 

94 except Exception: 

95 _get_solver.cholmod_factor = cholmod_analyze(M) 

96 _get_solver.cholmod_factor.cholesky_inplace(M) 

97 solve = _get_solver.cholmod_factor 

98 else: 

99 if has_umfpack and sym_pos: 

100 solve = sps.linalg.factorized(M) 

101 else: # factorized doesn't pass permc_spec 

102 solve = sps.linalg.splu(M, permc_spec=permc_spec).solve 

103 

104 else: 

105 if lstsq: # sometimes necessary as solution is approached 

106 def solve(r): 

107 return sp.linalg.lstsq(M, r)[0] 

108 elif cholesky: 

109 L = sp.linalg.cho_factor(M) 

110 

111 def solve(r): 

112 return sp.linalg.cho_solve(L, r) 

113 else: 

114 # this seems to cache the matrix factorization, so solving 

115 # with multiple right hand sides is much faster 

116 def solve(r, sym_pos=sym_pos): 

117 if sym_pos: 

118 return sp.linalg.solve(M, r, assume_a="pos") 

119 else: 

120 return sp.linalg.solve(M, r) 

121 # There are many things that can go wrong here, and it's hard to say 

122 # what all of them are. It doesn't really matter: if the matrix can't be 

123 # factorized, return None. get_solver will be called again with different 

124 # inputs, and a new routine will try to factorize the matrix. 

125 except KeyboardInterrupt: 

126 raise 

127 except Exception: 

128 return None 

129 return solve 

130 

131 

132def _get_delta(A, b, c, x, y, z, tau, kappa, gamma, eta, sparse=False, 

133 lstsq=False, sym_pos=True, cholesky=True, pc=True, ip=False, 

134 permc_spec='MMD_AT_PLUS_A'): 

135 """ 

136 Given standard form problem defined by ``A``, ``b``, and ``c``; 

137 current variable estimates ``x``, ``y``, ``z``, ``tau``, and ``kappa``; 

138 algorithmic parameters ``gamma and ``eta; 

139 and options ``sparse``, ``lstsq``, ``sym_pos``, ``cholesky``, ``pc`` 

140 (predictor-corrector), and ``ip`` (initial point improvement), 

141 get the search direction for increments to the variable estimates. 

142 

143 Parameters 

144 ---------- 

145 As defined in [4], except: 

146 sparse : bool 

147 True if the system to be solved is sparse. This is typically set 

148 True when the original ``A_ub`` and ``A_eq`` arrays are sparse. 

149 lstsq : bool 

150 True if the system is ill-conditioned and/or (nearly) singular and 

151 thus a more robust least-squares solver is desired. This is sometimes 

152 needed as the solution is approached. 

153 sym_pos : bool 

154 True if the system matrix is symmetric positive definite 

155 Sometimes this needs to be set false as the solution is approached, 

156 even when the system should be symmetric positive definite, due to 

157 numerical difficulties. 

158 cholesky : bool 

159 True if the system is to be solved by Cholesky, rather than LU, 

160 decomposition. This is typically faster unless the problem is very 

161 small or prone to numerical difficulties. 

162 pc : bool 

163 True if the predictor-corrector method of Mehrota is to be used. This 

164 is almost always (if not always) beneficial. Even though it requires 

165 the solution of an additional linear system, the factorization 

166 is typically (implicitly) reused so solution is efficient, and the 

167 number of algorithm iterations is typically reduced. 

168 ip : bool 

169 True if the improved initial point suggestion due to [4] section 4.3 

170 is desired. It's unclear whether this is beneficial. 

171 permc_spec : str (default = 'MMD_AT_PLUS_A') 

172 (Has effect only with ``sparse = True``, ``lstsq = False``, ``sym_pos = 

173 True``.) A matrix is factorized in each iteration of the algorithm. 

174 This option specifies how to permute the columns of the matrix for 

175 sparsity preservation. Acceptable values are: 

176 

177 - ``NATURAL``: natural ordering. 

178 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A. 

179 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A. 

180 - ``COLAMD``: approximate minimum degree column ordering. 

181 

182 This option can impact the convergence of the 

183 interior point algorithm; test different values to determine which 

184 performs best for your problem. For more information, refer to 

185 ``scipy.sparse.linalg.splu``. 

186 

187 Returns 

188 ------- 

189 Search directions as defined in [4] 

190 

191 References 

192 ---------- 

193 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

194 optimizer for linear programming: an implementation of the 

195 homogeneous algorithm." High performance optimization. Springer US, 

196 2000. 197-232. 

197 

198 """ 

199 if A.shape[0] == 0: 

200 # If there are no constraints, some solvers fail (understandably) 

201 # rather than returning empty solution. This gets the job done. 

202 sparse, lstsq, sym_pos, cholesky = False, False, True, False 

203 n_x = len(x) 

204 

205 # [4] Equation 8.8 

206 r_P = b * tau - A.dot(x) 

207 r_D = c * tau - A.T.dot(y) - z 

208 r_G = c.dot(x) - b.transpose().dot(y) + kappa 

209 mu = (x.dot(z) + tau * kappa) / (n_x + 1) 

210 

211 # Assemble M from [4] Equation 8.31 

212 Dinv = x / z 

213 

214 if sparse: 

215 M = A.dot(sps.diags(Dinv, 0, format="csc").dot(A.T)) 

216 else: 

217 M = A.dot(Dinv.reshape(-1, 1) * A.T) 

218 solve = _get_solver(M, sparse, lstsq, sym_pos, cholesky, permc_spec) 

219 

220 # pc: "predictor-corrector" [4] Section 4.1 

221 # In development this option could be turned off 

222 # but it always seems to improve performance substantially 

223 n_corrections = 1 if pc else 0 

224 

225 i = 0 

226 alpha, d_x, d_z, d_tau, d_kappa = 0, 0, 0, 0, 0 

227 while i <= n_corrections: 

228 # Reference [4] Eq. 8.6 

229 rhatp = eta(gamma) * r_P 

230 rhatd = eta(gamma) * r_D 

231 rhatg = eta(gamma) * r_G 

232 

233 # Reference [4] Eq. 8.7 

234 rhatxs = gamma * mu - x * z 

235 rhattk = gamma * mu - tau * kappa 

236 

237 if i == 1: 

238 if ip: # if the correction is to get "initial point" 

239 # Reference [4] Eq. 8.23 

240 rhatxs = ((1 - alpha) * gamma * mu - 

241 x * z - alpha**2 * d_x * d_z) 

242 rhattk = ((1 - alpha) * gamma * mu - 

243 tau * kappa - 

244 alpha**2 * d_tau * d_kappa) 

245 else: # if the correction is for "predictor-corrector" 

246 # Reference [4] Eq. 8.13 

247 rhatxs -= d_x * d_z 

248 rhattk -= d_tau * d_kappa 

249 

250 # sometimes numerical difficulties arise as the solution is approached 

251 # this loop tries to solve the equations using a sequence of functions 

252 # for solve. For dense systems, the order is: 

253 # 1. scipy.linalg.cho_factor/scipy.linalg.cho_solve, 

254 # 2. scipy.linalg.solve w/ sym_pos = True, 

255 # 3. scipy.linalg.solve w/ sym_pos = False, and if all else fails 

256 # 4. scipy.linalg.lstsq 

257 # For sparse systems, the order is: 

258 # 1. sksparse.cholmod.cholesky (if available) 

259 # 2. scipy.sparse.linalg.factorized (if umfpack available) 

260 # 3. scipy.sparse.linalg.splu 

261 # 4. scipy.sparse.linalg.lsqr 

262 solved = False 

263 while not solved: 

264 try: 

265 # [4] Equation 8.28 

266 p, q = _sym_solve(Dinv, A, c, b, solve) 

267 # [4] Equation 8.29 

268 u, v = _sym_solve(Dinv, A, rhatd - 

269 (1 / x) * rhatxs, rhatp, solve) 

270 if np.any(np.isnan(p)) or np.any(np.isnan(q)): 

271 raise LinAlgError 

272 solved = True 

273 except (LinAlgError, ValueError, TypeError) as e: 

274 # Usually this doesn't happen. If it does, it happens when 

275 # there are redundant constraints or when approaching the 

276 # solution. If so, change solver. 

277 if cholesky: 

278 cholesky = False 

279 warn( 

280 "Solving system with option 'cholesky':True " 

281 "failed. It is normal for this to happen " 

282 "occasionally, especially as the solution is " 

283 "approached. However, if you see this frequently, " 

284 "consider setting option 'cholesky' to False.", 

285 OptimizeWarning, stacklevel=5) 

286 elif sym_pos: 

287 sym_pos = False 

288 warn( 

289 "Solving system with option 'sym_pos':True " 

290 "failed. It is normal for this to happen " 

291 "occasionally, especially as the solution is " 

292 "approached. However, if you see this frequently, " 

293 "consider setting option 'sym_pos' to False.", 

294 OptimizeWarning, stacklevel=5) 

295 elif not lstsq: 

296 lstsq = True 

297 warn( 

298 "Solving system with option 'sym_pos':False " 

299 "failed. This may happen occasionally, " 

300 "especially as the solution is " 

301 "approached. However, if you see this frequently, " 

302 "your problem may be numerically challenging. " 

303 "If you cannot improve the formulation, consider " 

304 "setting 'lstsq' to True. Consider also setting " 

305 "`presolve` to True, if it is not already.", 

306 OptimizeWarning, stacklevel=5) 

307 else: 

308 raise e 

309 solve = _get_solver(M, sparse, lstsq, sym_pos, 

310 cholesky, permc_spec) 

311 # [4] Results after 8.29 

312 d_tau = ((rhatg + 1 / tau * rhattk - (-c.dot(u) + b.dot(v))) / 

313 (1 / tau * kappa + (-c.dot(p) + b.dot(q)))) 

314 d_x = u + p * d_tau 

315 d_y = v + q * d_tau 

316 

317 # [4] Relations between after 8.25 and 8.26 

318 d_z = (1 / x) * (rhatxs - z * d_x) 

319 d_kappa = 1 / tau * (rhattk - kappa * d_tau) 

320 

321 # [4] 8.12 and "Let alpha be the maximal possible step..." before 8.23 

322 alpha = _get_step(x, d_x, z, d_z, tau, d_tau, kappa, d_kappa, 1) 

323 if ip: # initial point - see [4] 4.4 

324 gamma = 10 

325 else: # predictor-corrector, [4] definition after 8.12 

326 beta1 = 0.1 # [4] pg. 220 (Table 8.1) 

327 gamma = (1 - alpha)**2 * min(beta1, (1 - alpha)) 

328 i += 1 

329 

330 return d_x, d_y, d_z, d_tau, d_kappa 

331 

332 

333def _sym_solve(Dinv, A, r1, r2, solve): 

334 """ 

335 An implementation of [4] equation 8.31 and 8.32 

336 

337 References 

338 ---------- 

339 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

340 optimizer for linear programming: an implementation of the 

341 homogeneous algorithm." High performance optimization. Springer US, 

342 2000. 197-232. 

343 

344 """ 

345 # [4] 8.31 

346 r = r2 + A.dot(Dinv * r1) 

347 v = solve(r) 

348 # [4] 8.32 

349 u = Dinv * (A.T.dot(v) - r1) 

350 return u, v 

351 

352 

353def _get_step(x, d_x, z, d_z, tau, d_tau, kappa, d_kappa, alpha0): 

354 """ 

355 An implementation of [4] equation 8.21 

356 

357 References 

358 ---------- 

359 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

360 optimizer for linear programming: an implementation of the 

361 homogeneous algorithm." High performance optimization. Springer US, 

362 2000. 197-232. 

363 

364 """ 

365 # [4] 4.3 Equation 8.21, ignoring 8.20 requirement 

366 # same step is taken in primal and dual spaces 

367 # alpha0 is basically beta3 from [4] Table 8.1, but instead of beta3 

368 # the value 1 is used in Mehrota corrector and initial point correction 

369 i_x = d_x < 0 

370 i_z = d_z < 0 

371 alpha_x = alpha0 * np.min(x[i_x] / -d_x[i_x]) if np.any(i_x) else 1 

372 alpha_tau = alpha0 * tau / -d_tau if d_tau < 0 else 1 

373 alpha_z = alpha0 * np.min(z[i_z] / -d_z[i_z]) if np.any(i_z) else 1 

374 alpha_kappa = alpha0 * kappa / -d_kappa if d_kappa < 0 else 1 

375 alpha = np.min([1, alpha_x, alpha_tau, alpha_z, alpha_kappa]) 

376 return alpha 

377 

378 

379def _get_message(status): 

380 """ 

381 Given problem status code, return a more detailed message. 

382 

383 Parameters 

384 ---------- 

385 status : int 

386 An integer representing the exit status of the optimization:: 

387 

388 0 : Optimization terminated successfully 

389 1 : Iteration limit reached 

390 2 : Problem appears to be infeasible 

391 3 : Problem appears to be unbounded 

392 4 : Serious numerical difficulties encountered 

393 

394 Returns 

395 ------- 

396 message : str 

397 A string descriptor of the exit status of the optimization. 

398 

399 """ 

400 messages = ( 

401 ["Optimization terminated successfully.", 

402 "The iteration limit was reached before the algorithm converged.", 

403 "The algorithm terminated successfully and determined that the " 

404 "problem is infeasible.", 

405 "The algorithm terminated successfully and determined that the " 

406 "problem is unbounded.", 

407 "Numerical difficulties were encountered before the problem " 

408 "converged. Please check your problem formulation for errors, " 

409 "independence of linear equality constraints, and reasonable " 

410 "scaling and matrix condition numbers. If you continue to " 

411 "encounter this error, please submit a bug report." 

412 ]) 

413 return messages[status] 

414 

415 

416def _do_step(x, y, z, tau, kappa, d_x, d_y, d_z, d_tau, d_kappa, alpha): 

417 """ 

418 An implementation of [4] Equation 8.9 

419 

420 References 

421 ---------- 

422 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

423 optimizer for linear programming: an implementation of the 

424 homogeneous algorithm." High performance optimization. Springer US, 

425 2000. 197-232. 

426 

427 """ 

428 x = x + alpha * d_x 

429 tau = tau + alpha * d_tau 

430 z = z + alpha * d_z 

431 kappa = kappa + alpha * d_kappa 

432 y = y + alpha * d_y 

433 return x, y, z, tau, kappa 

434 

435 

436def _get_blind_start(shape): 

437 """ 

438 Return the starting point from [4] 4.4 

439 

440 References 

441 ---------- 

442 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

443 optimizer for linear programming: an implementation of the 

444 homogeneous algorithm." High performance optimization. Springer US, 

445 2000. 197-232. 

446 

447 """ 

448 m, n = shape 

449 x0 = np.ones(n) 

450 y0 = np.zeros(m) 

451 z0 = np.ones(n) 

452 tau0 = 1 

453 kappa0 = 1 

454 return x0, y0, z0, tau0, kappa0 

455 

456 

457def _indicators(A, b, c, c0, x, y, z, tau, kappa): 

458 """ 

459 Implementation of several equations from [4] used as indicators of 

460 the status of optimization. 

461 

462 References 

463 ---------- 

464 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

465 optimizer for linear programming: an implementation of the 

466 homogeneous algorithm." High performance optimization. Springer US, 

467 2000. 197-232. 

468 

469 """ 

470 

471 # residuals for termination are relative to initial values 

472 x0, y0, z0, tau0, kappa0 = _get_blind_start(A.shape) 

473 

474 # See [4], Section 4 - The Homogeneous Algorithm, Equation 8.8 

475 def r_p(x, tau): 

476 return b * tau - A.dot(x) 

477 

478 def r_d(y, z, tau): 

479 return c * tau - A.T.dot(y) - z 

480 

481 def r_g(x, y, kappa): 

482 return kappa + c.dot(x) - b.dot(y) 

483 

484 # np.dot unpacks if they are arrays of size one 

485 def mu(x, tau, z, kappa): 

486 return (x.dot(z) + np.dot(tau, kappa)) / (len(x) + 1) 

487 

488 obj = c.dot(x / tau) + c0 

489 

490 def norm(a): 

491 return np.linalg.norm(a) 

492 

493 # See [4], Section 4.5 - The Stopping Criteria 

494 r_p0 = r_p(x0, tau0) 

495 r_d0 = r_d(y0, z0, tau0) 

496 r_g0 = r_g(x0, y0, kappa0) 

497 mu_0 = mu(x0, tau0, z0, kappa0) 

498 rho_A = norm(c.T.dot(x) - b.T.dot(y)) / (tau + norm(b.T.dot(y))) 

499 rho_p = norm(r_p(x, tau)) / max(1, norm(r_p0)) 

500 rho_d = norm(r_d(y, z, tau)) / max(1, norm(r_d0)) 

501 rho_g = norm(r_g(x, y, kappa)) / max(1, norm(r_g0)) 

502 rho_mu = mu(x, tau, z, kappa) / mu_0 

503 return rho_p, rho_d, rho_A, rho_g, rho_mu, obj 

504 

505 

506def _display_iter(rho_p, rho_d, rho_g, alpha, rho_mu, obj, header=False): 

507 """ 

508 Print indicators of optimization status to the console. 

509 

510 Parameters 

511 ---------- 

512 rho_p : float 

513 The (normalized) primal feasibility, see [4] 4.5 

514 rho_d : float 

515 The (normalized) dual feasibility, see [4] 4.5 

516 rho_g : float 

517 The (normalized) duality gap, see [4] 4.5 

518 alpha : float 

519 The step size, see [4] 4.3 

520 rho_mu : float 

521 The (normalized) path parameter, see [4] 4.5 

522 obj : float 

523 The objective function value of the current iterate 

524 header : bool 

525 True if a header is to be printed 

526 

527 References 

528 ---------- 

529 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

530 optimizer for linear programming: an implementation of the 

531 homogeneous algorithm." High performance optimization. Springer US, 

532 2000. 197-232. 

533 

534 """ 

535 if header: 

536 print("Primal Feasibility ", 

537 "Dual Feasibility ", 

538 "Duality Gap ", 

539 "Step ", 

540 "Path Parameter ", 

541 "Objective ") 

542 

543 # no clue why this works 

544 fmt = '{0:<20.13}{1:<20.13}{2:<20.13}{3:<17.13}{4:<20.13}{5:<20.13}' 

545 print(fmt.format( 

546 float(rho_p), 

547 float(rho_d), 

548 float(rho_g), 

549 alpha if isinstance(alpha, str) else float(alpha), 

550 float(rho_mu), 

551 float(obj))) 

552 

553 

554def _ip_hsd(A, b, c, c0, alpha0, beta, maxiter, disp, tol, sparse, lstsq, 

555 sym_pos, cholesky, pc, ip, permc_spec, callback, postsolve_args): 

556 r""" 

557 Solve a linear programming problem in standard form: 

558 

559 Minimize:: 

560 

561 c @ x 

562 

563 Subject to:: 

564 

565 A @ x == b 

566 x >= 0 

567 

568 using the interior point method of [4]. 

569 

570 Parameters 

571 ---------- 

572 A : 2-D array 

573 2-D array such that ``A @ x``, gives the values of the equality 

574 constraints at ``x``. 

575 b : 1-D array 

576 1-D array of values representing the RHS of each equality constraint 

577 (row) in ``A`` (for standard form problem). 

578 c : 1-D array 

579 Coefficients of the linear objective function to be minimized (for 

580 standard form problem). 

581 c0 : float 

582 Constant term in objective function due to fixed (and eliminated) 

583 variables. (Purely for display.) 

584 alpha0 : float 

585 The maximal step size for Mehrota's predictor-corrector search 

586 direction; see :math:`\beta_3`of [4] Table 8.1 

587 beta : float 

588 The desired reduction of the path parameter :math:`\mu` (see [6]_) 

589 maxiter : int 

590 The maximum number of iterations of the algorithm. 

591 disp : bool 

592 Set to ``True`` if indicators of optimization status are to be printed 

593 to the console each iteration. 

594 tol : float 

595 Termination tolerance; see [4]_ Section 4.5. 

596 sparse : bool 

597 Set to ``True`` if the problem is to be treated as sparse. However, 

598 the inputs ``A_eq`` and ``A_ub`` should nonetheless be provided as 

599 (dense) arrays rather than sparse matrices. 

600 lstsq : bool 

601 Set to ``True`` if the problem is expected to be very poorly 

602 conditioned. This should always be left as ``False`` unless severe 

603 numerical difficulties are frequently encountered, and a better option 

604 would be to improve the formulation of the problem. 

605 sym_pos : bool 

606 Leave ``True`` if the problem is expected to yield a well conditioned 

607 symmetric positive definite normal equation matrix (almost always). 

608 cholesky : bool 

609 Set to ``True`` if the normal equations are to be solved by explicit 

610 Cholesky decomposition followed by explicit forward/backward 

611 substitution. This is typically faster for moderate, dense problems 

612 that are numerically well-behaved. 

613 pc : bool 

614 Leave ``True`` if the predictor-corrector method of Mehrota is to be 

615 used. This is almost always (if not always) beneficial. 

616 ip : bool 

617 Set to ``True`` if the improved initial point suggestion due to [4]_ 

618 Section 4.3 is desired. It's unclear whether this is beneficial. 

619 permc_spec : str (default = 'MMD_AT_PLUS_A') 

620 (Has effect only with ``sparse = True``, ``lstsq = False``, ``sym_pos = 

621 True``.) A matrix is factorized in each iteration of the algorithm. 

622 This option specifies how to permute the columns of the matrix for 

623 sparsity preservation. Acceptable values are: 

624 

625 - ``NATURAL``: natural ordering. 

626 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A. 

627 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A. 

628 - ``COLAMD``: approximate minimum degree column ordering. 

629 

630 This option can impact the convergence of the 

631 interior point algorithm; test different values to determine which 

632 performs best for your problem. For more information, refer to 

633 ``scipy.sparse.linalg.splu``. 

634 callback : callable, optional 

635 If a callback function is provided, it will be called within each 

636 iteration of the algorithm. The callback function must accept a single 

637 `scipy.optimize.OptimizeResult` consisting of the following fields: 

638 

639 x : 1-D array 

640 Current solution vector 

641 fun : float 

642 Current value of the objective function 

643 success : bool 

644 True only when an algorithm has completed successfully, 

645 so this is always False as the callback function is called 

646 only while the algorithm is still iterating. 

647 slack : 1-D array 

648 The values of the slack variables. Each slack variable 

649 corresponds to an inequality constraint. If the slack is zero, 

650 the corresponding constraint is active. 

651 con : 1-D array 

652 The (nominally zero) residuals of the equality constraints, 

653 that is, ``b - A_eq @ x`` 

654 phase : int 

655 The phase of the algorithm being executed. This is always 

656 1 for the interior-point method because it has only one phase. 

657 status : int 

658 For revised simplex, this is always 0 because if a different 

659 status is detected, the algorithm terminates. 

660 nit : int 

661 The number of iterations performed. 

662 message : str 

663 A string descriptor of the exit status of the optimization. 

664 postsolve_args : tuple 

665 Data needed by _postsolve to convert the solution to the standard-form 

666 problem into the solution to the original problem. 

667 

668 Returns 

669 ------- 

670 x_hat : float 

671 Solution vector (for standard form problem). 

672 status : int 

673 An integer representing the exit status of the optimization:: 

674 

675 0 : Optimization terminated successfully 

676 1 : Iteration limit reached 

677 2 : Problem appears to be infeasible 

678 3 : Problem appears to be unbounded 

679 4 : Serious numerical difficulties encountered 

680 

681 message : str 

682 A string descriptor of the exit status of the optimization. 

683 iteration : int 

684 The number of iterations taken to solve the problem 

685 

686 References 

687 ---------- 

688 .. [4] Andersen, Erling D., and Knud D. Andersen. "The MOSEK interior point 

689 optimizer for linear programming: an implementation of the 

690 homogeneous algorithm." High performance optimization. Springer US, 

691 2000. 197-232. 

692 .. [6] Freund, Robert M. "Primal-Dual Interior-Point Methods for Linear 

693 Programming based on Newton's Method." Unpublished Course Notes, 

694 March 2004. Available 2/25/2017 at: 

695 https://ocw.mit.edu/courses/sloan-school-of-management/15-084j-nonlinear-programming-spring-2004/lecture-notes/lec14_int_pt_mthd.pdf 

696 

697 """ 

698 

699 iteration = 0 

700 

701 # default initial point 

702 x, y, z, tau, kappa = _get_blind_start(A.shape) 

703 

704 # first iteration is special improvement of initial point 

705 ip = ip if pc else False 

706 

707 # [4] 4.5 

708 rho_p, rho_d, rho_A, rho_g, rho_mu, obj = _indicators( 

709 A, b, c, c0, x, y, z, tau, kappa) 

710 go = rho_p > tol or rho_d > tol or rho_A > tol # we might get lucky : ) 

711 

712 if disp: 

713 _display_iter(rho_p, rho_d, rho_g, "-", rho_mu, obj, header=True) 

714 if callback is not None: 

715 x_o, fun, slack, con = _postsolve(x/tau, postsolve_args) 

716 res = OptimizeResult({'x': x_o, 'fun': fun, 'slack': slack, 

717 'con': con, 'nit': iteration, 'phase': 1, 

718 'complete': False, 'status': 0, 

719 'message': "", 'success': False}) 

720 callback(res) 

721 

722 status = 0 

723 message = "Optimization terminated successfully." 

724 

725 if sparse: 

726 A = sps.csc_matrix(A) 

727 A.T = A.transpose() # A.T is defined for sparse matrices but is slow 

728 # Redefine it to avoid calculating again 

729 # This is fine as long as A doesn't change 

730 

731 while go: 

732 

733 iteration += 1 

734 

735 if ip: # initial point 

736 # [4] Section 4.4 

737 gamma = 1 

738 

739 def eta(g): 

740 return 1 

741 else: 

742 # gamma = 0 in predictor step according to [4] 4.1 

743 # if predictor/corrector is off, use mean of complementarity [6] 

744 # 5.1 / [4] Below Figure 10-4 

745 gamma = 0 if pc else beta * np.mean(z * x) 

746 # [4] Section 4.1 

747 

748 def eta(g=gamma): 

749 return 1 - g 

750 

751 try: 

752 # Solve [4] 8.6 and 8.7/8.13/8.23 

753 d_x, d_y, d_z, d_tau, d_kappa = _get_delta( 

754 A, b, c, x, y, z, tau, kappa, gamma, eta, 

755 sparse, lstsq, sym_pos, cholesky, pc, ip, permc_spec) 

756 

757 if ip: # initial point 

758 # [4] 4.4 

759 # Formula after 8.23 takes a full step regardless if this will 

760 # take it negative 

761 alpha = 1.0 

762 x, y, z, tau, kappa = _do_step( 

763 x, y, z, tau, kappa, d_x, d_y, 

764 d_z, d_tau, d_kappa, alpha) 

765 x[x < 1] = 1 

766 z[z < 1] = 1 

767 tau = max(1, tau) 

768 kappa = max(1, kappa) 

769 ip = False # done with initial point 

770 else: 

771 # [4] Section 4.3 

772 alpha = _get_step(x, d_x, z, d_z, tau, 

773 d_tau, kappa, d_kappa, alpha0) 

774 # [4] Equation 8.9 

775 x, y, z, tau, kappa = _do_step( 

776 x, y, z, tau, kappa, d_x, d_y, d_z, d_tau, d_kappa, alpha) 

777 

778 except (LinAlgError, FloatingPointError, 

779 ValueError, ZeroDivisionError): 

780 # this can happen when sparse solver is used and presolve 

781 # is turned off. Also observed ValueError in AppVeyor Python 3.6 

782 # Win32 build (PR #8676). I've never seen it otherwise. 

783 status = 4 

784 message = _get_message(status) 

785 break 

786 

787 # [4] 4.5 

788 rho_p, rho_d, rho_A, rho_g, rho_mu, obj = _indicators( 

789 A, b, c, c0, x, y, z, tau, kappa) 

790 go = rho_p > tol or rho_d > tol or rho_A > tol 

791 

792 if disp: 

793 _display_iter(rho_p, rho_d, rho_g, alpha, rho_mu, obj) 

794 if callback is not None: 

795 x_o, fun, slack, con = _postsolve(x/tau, postsolve_args) 

796 res = OptimizeResult({'x': x_o, 'fun': fun, 'slack': slack, 

797 'con': con, 'nit': iteration, 'phase': 1, 

798 'complete': False, 'status': 0, 

799 'message': "", 'success': False}) 

800 callback(res) 

801 

802 # [4] 4.5 

803 inf1 = (rho_p < tol and rho_d < tol and rho_g < tol and tau < tol * 

804 max(1, kappa)) 

805 inf2 = rho_mu < tol and tau < tol * min(1, kappa) 

806 if inf1 or inf2: 

807 # [4] Lemma 8.4 / Theorem 8.3 

808 if b.transpose().dot(y) > tol: 

809 status = 2 

810 else: # elif c.T.dot(x) < tol: ? Probably not necessary. 

811 status = 3 

812 message = _get_message(status) 

813 break 

814 elif iteration >= maxiter: 

815 status = 1 

816 message = _get_message(status) 

817 break 

818 

819 x_hat = x / tau 

820 # [4] Statement after Theorem 8.2 

821 return x_hat, status, message, iteration 

822 

823 

824def _linprog_ip(c, c0, A, b, callback, postsolve_args, maxiter=1000, tol=1e-8, 

825 disp=False, alpha0=.99995, beta=0.1, sparse=False, lstsq=False, 

826 sym_pos=True, cholesky=None, pc=True, ip=False, 

827 permc_spec='MMD_AT_PLUS_A', **unknown_options): 

828 r""" 

829 Minimize a linear objective function subject to linear 

830 equality and non-negativity constraints using the interior point method 

831 of [4]_. Linear programming is intended to solve problems 

832 of the following form: 

833 

834 Minimize:: 

835 

836 c @ x 

837 

838 Subject to:: 

839 

840 A @ x == b 

841 x >= 0 

842 

843 User-facing documentation is in _linprog_doc.py. 

844 

845 Parameters 

846 ---------- 

847 c : 1-D array 

848 Coefficients of the linear objective function to be minimized. 

849 c0 : float 

850 Constant term in objective function due to fixed (and eliminated) 

851 variables. (Purely for display.) 

852 A : 2-D array 

853 2-D array such that ``A @ x``, gives the values of the equality 

854 constraints at ``x``. 

855 b : 1-D array 

856 1-D array of values representing the right hand side of each equality 

857 constraint (row) in ``A``. 

858 callback : callable, optional 

859 Callback function to be executed once per iteration. 

860 postsolve_args : tuple 

861 Data needed by _postsolve to convert the solution to the standard-form 

862 problem into the solution to the original problem. 

863 

864 Options 

865 ------- 

866 maxiter : int (default = 1000) 

867 The maximum number of iterations of the algorithm. 

868 tol : float (default = 1e-8) 

869 Termination tolerance to be used for all termination criteria; 

870 see [4]_ Section 4.5. 

871 disp : bool (default = False) 

872 Set to ``True`` if indicators of optimization status are to be printed 

873 to the console each iteration. 

874 alpha0 : float (default = 0.99995) 

875 The maximal step size for Mehrota's predictor-corrector search 

876 direction; see :math:`\beta_{3}` of [4]_ Table 8.1. 

877 beta : float (default = 0.1) 

878 The desired reduction of the path parameter :math:`\mu` (see [6]_) 

879 when Mehrota's predictor-corrector is not in use (uncommon). 

880 sparse : bool (default = False) 

881 Set to ``True`` if the problem is to be treated as sparse after 

882 presolve. If either ``A_eq`` or ``A_ub`` is a sparse matrix, 

883 this option will automatically be set ``True``, and the problem 

884 will be treated as sparse even during presolve. If your constraint 

885 matrices contain mostly zeros and the problem is not very small (less 

886 than about 100 constraints or variables), consider setting ``True`` 

887 or providing ``A_eq`` and ``A_ub`` as sparse matrices. 

888 lstsq : bool (default = False) 

889 Set to ``True`` if the problem is expected to be very poorly 

890 conditioned. This should always be left ``False`` unless severe 

891 numerical difficulties are encountered. Leave this at the default 

892 unless you receive a warning message suggesting otherwise. 

893 sym_pos : bool (default = True) 

894 Leave ``True`` if the problem is expected to yield a well conditioned 

895 symmetric positive definite normal equation matrix 

896 (almost always). Leave this at the default unless you receive 

897 a warning message suggesting otherwise. 

898 cholesky : bool (default = True) 

899 Set to ``True`` if the normal equations are to be solved by explicit 

900 Cholesky decomposition followed by explicit forward/backward 

901 substitution. This is typically faster for problems 

902 that are numerically well-behaved. 

903 pc : bool (default = True) 

904 Leave ``True`` if the predictor-corrector method of Mehrota is to be 

905 used. This is almost always (if not always) beneficial. 

906 ip : bool (default = False) 

907 Set to ``True`` if the improved initial point suggestion due to [4]_ 

908 Section 4.3 is desired. Whether this is beneficial or not 

909 depends on the problem. 

910 permc_spec : str (default = 'MMD_AT_PLUS_A') 

911 (Has effect only with ``sparse = True``, ``lstsq = False``, ``sym_pos = 

912 True``, and no SuiteSparse.) 

913 A matrix is factorized in each iteration of the algorithm. 

914 This option specifies how to permute the columns of the matrix for 

915 sparsity preservation. Acceptable values are: 

916 

917 - ``NATURAL``: natural ordering. 

918 - ``MMD_ATA``: minimum degree ordering on the structure of A^T A. 

919 - ``MMD_AT_PLUS_A``: minimum degree ordering on the structure of A^T+A. 

920 - ``COLAMD``: approximate minimum degree column ordering. 

921 

922 This option can impact the convergence of the 

923 interior point algorithm; test different values to determine which 

924 performs best for your problem. For more information, refer to 

925 ``scipy.sparse.linalg.splu``. 

926 unknown_options : dict 

927 Optional arguments not used by this particular solver. If 

928 `unknown_options` is non-empty a warning is issued listing all 

929 unused options. 

930 

931 Returns 

932 ------- 

933 x : 1-D array 

934 Solution vector. 

935 status : int 

936 An integer representing the exit status of the optimization:: 

937 

938 0 : Optimization terminated successfully 

939 1 : Iteration limit reached 

940 2 : Problem appears to be infeasible 

941 3 : Problem appears to be unbounded 

942 4 : Serious numerical difficulties encountered 

943 

944 message : str 

945 A string descriptor of the exit status of the optimization. 

946 iteration : int 

947 The number of iterations taken to solve the problem. 

948 

949 Notes 

950 ----- 

951 This method implements the algorithm outlined in [4]_ with ideas from [8]_ 

952 and a structure inspired by the simpler methods of [6]_. 

953 

954 The primal-dual path following method begins with initial 'guesses' of 

955 the primal and dual variables of the standard form problem and iteratively 

956 attempts to solve the (nonlinear) Karush-Kuhn-Tucker conditions for the 

957 problem with a gradually reduced logarithmic barrier term added to the 

958 objective. This particular implementation uses a homogeneous self-dual 

959 formulation, which provides certificates of infeasibility or unboundedness 

960 where applicable. 

961 

962 The default initial point for the primal and dual variables is that 

963 defined in [4]_ Section 4.4 Equation 8.22. Optionally (by setting initial 

964 point option ``ip=True``), an alternate (potentially improved) starting 

965 point can be calculated according to the additional recommendations of 

966 [4]_ Section 4.4. 

967 

968 A search direction is calculated using the predictor-corrector method 

969 (single correction) proposed by Mehrota and detailed in [4]_ Section 4.1. 

970 (A potential improvement would be to implement the method of multiple 

971 corrections described in [4]_ Section 4.2.) In practice, this is 

972 accomplished by solving the normal equations, [4]_ Section 5.1 Equations 

973 8.31 and 8.32, derived from the Newton equations [4]_ Section 5 Equations 

974 8.25 (compare to [4]_ Section 4 Equations 8.6-8.8). The advantage of