Coverage for /usr/lib/python3/dist-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 # noqa: F401 

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 

728 while go: 

729 

730 iteration += 1 

731 

732 if ip: # initial point 

733 # [4] Section 4.4 

734 gamma = 1 

735 

736 def eta(g): 

737 return 1 

738 else: 

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

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

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

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

743 # [4] Section 4.1 

744 

745 def eta(g=gamma): 

746 return 1 - g 

747 

748 try: 

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

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

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

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

753 

754 if ip: # initial point 

755 # [4] 4.4 

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

757 # take it negative 

758 alpha = 1.0 

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

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

761 d_z, d_tau, d_kappa, alpha) 

762 x[x < 1] = 1 

763 z[z < 1] = 1 

764 tau = max(1, tau) 

765 kappa = max(1, kappa) 

766 ip = False # done with initial point 

767 else: 

768 # [4] Section 4.3 

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

770 d_tau, kappa, d_kappa, alpha0) 

771 # [4] Equation 8.9 

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

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

774 

775 except (LinAlgError, FloatingPointError, 

776 ValueError, ZeroDivisionError): 

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

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

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

780 status = 4 

781 message = _get_message(status) 

782 break 

783 

784 # [4] 4.5 

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

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

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

788 

789 if disp: 

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

791 if callback is not None: 

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

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

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

795 'complete': False, 'status': 0, 

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

797 callback(res) 

798 

799 # [4] 4.5 

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

801 max(1, kappa)) 

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

803 if inf1 or inf2: 

804 # [4] Lemma 8.4 / Theorem 8.3 

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

806 status = 2 

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

808 status = 3 

809 message = _get_message(status) 

810 break 

811 elif iteration >= maxiter: 

812 status = 1 

813 message = _get_message(status) 

814 break 

815 

816 x_hat = x / tau 

817 # [4] Statement after Theorem 8.2 

818 return x_hat, status, message, iteration 

819 

820 

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

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

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

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

825 r""" 

826 Minimize a linear objective function subject to linear 

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

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

829 of the following form: 

830 

831 Minimize:: 

832 

833 c @ x 

834 

835 Subject to:: 

836 

837 A @ x == b 

838 x >= 0 

839 

840 User-facing documentation is in _linprog_doc.py. 

841 

842 Parameters 

843 ---------- 

844 c : 1-D array 

845 Coefficients of the linear objective function to be minimized. 

846 c0 : float 

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

848 variables. (Purely for display.) 

849 A : 2-D array 

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

851 constraints at ``x``. 

852 b : 1-D array 

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

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

855 callback : callable, optional 

856 Callback function to be executed once per iteration. 

857 postsolve_args : tuple 

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

859 problem into the solution to the original problem. 

860 

861 Options 

862 ------- 

863 maxiter : int (default = 1000) 

864 The maximum number of iterations of the algorithm. 

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

866 Termination tolerance to be used for all termination criteria; 

867 see [4]_ Section 4.5. 

868 disp : bool (default = False) 

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

870 to the console each iteration. 

871 alpha0 : float (default = 0.99995) 

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

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

874 beta : float (default = 0.1) 

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

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

877 sparse : bool (default = False) 

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

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

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

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

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

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

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

885 lstsq : bool (default = False) 

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

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

888 numerical difficulties are encountered. Leave this at the default 

889 unless you receive a warning message suggesting otherwise. 

890 sym_pos : bool (default = True) 

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

892 symmetric positive definite normal equation matrix 

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

894 a warning message suggesting otherwise. 

895 cholesky : bool (default = True) 

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

897 Cholesky decomposition followed by explicit forward/backward 

898 substitution. This is typically faster for problems 

899 that are numerically well-behaved. 

900 pc : bool (default = True) 

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

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

903 ip : bool (default = False) 

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

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

906 depends on the problem. 

907 permc_spec : str (default = 'MMD_AT_PLUS_A') 

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

909 True``, and no SuiteSparse.) 

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

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

912 sparsity preservation. Acceptable values are: 

913 

914 - ``NATURAL``: natural ordering. 

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

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

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

918 

919 This option can impact the convergence of the 

920 interior point algorithm; test different values to determine which 

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

922 ``scipy.sparse.linalg.splu``. 

923 unknown_options : dict 

924 Optional arguments not used by this particular solver. If 

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

926 unused options. 

927 

928 Returns 

929 ------- 

930 x : 1-D array 

931 Solution vector. 

932 status : int 

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

934 

935 0 : Optimization terminated successfully 

936 1 : Iteration limit reached 

937 2 : Problem appears to be infeasible 

938 3 : Problem appears to be unbounded 

939 4 : Serious numerical difficulties encountered 

940 

941 message : str 

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

943 iteration : int 

944 The number of iterations taken to solve the problem. 

945 

946 Notes 

947 ----- 

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

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

950 

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

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

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

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

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

956 formulation, which provides certificates of infeasibility or unboundedness 

957 where applicable. 

958 

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

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

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

962 point can be calculated according to the additional recommendations of 

963 [4]_ Section 4.4. 

964 

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

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

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

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

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

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

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

972 solving the normal equations rather than 8.25 directly is that the 

973 matrices involved are symmetric positive definite, so Cholesky 

974 decomposition can be used rather than the more expensive LU factorization. 

975 

976 With default options, the solver used to perform the factorization depends