Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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 from sksparse.cholmod import cholesky as cholmod 

32 from sksparse.cholmod import analyze as cholmod_analyze 

33except ImportError: 

34 has_cholmod = False 

35try: 

36 import scikits.umfpack # test whether to use factorized 

37except ImportError: 

38 has_umfpack = False 

39 

40 

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

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

43 """ 

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

45 solver. 

46 

47 Parameters 

48 ---------- 

49 M : 2-D array 

50 As defined in [4] Equation 8.31 

51 sparse : bool (default = False) 

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

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

54 lstsq : bool (default = False) 

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

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

57 needed as the solution is approached. 

58 sym_pos : bool (default = True) 

59 True if the system matrix is symmetric positive definite 

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

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

62 numerical difficulties. 

63 cholesky : bool (default = True) 

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

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

66 small or prone to numerical difficulties. 

67 permc_spec : str (default = 'MMD_AT_PLUS_A') 

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

69 

70 - ``NATURAL``: natural ordering. 

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

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

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

74 

75 See SuperLU documentation. 

76 

77 Returns 

78 ------- 

79 solve : function 

80 Handle to the appropriate solver function 

81 

82 """ 

83 try: 

84 if sparse: 

85 if lstsq: 

86 def solve(r, sym_pos=False): 

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

88 elif cholesky: 

89 try: 

90 # Will raise an exception in the first call, 

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

92 _get_solver.cholmod_factor.cholesky_inplace(M) 

93 except Exception: 

94 _get_solver.cholmod_factor = cholmod_analyze(M) 

95 _get_solver.cholmod_factor.cholesky_inplace(M) 

96 solve = _get_solver.cholmod_factor 

97 else: 

98 if has_umfpack and sym_pos: 

99 solve = sps.linalg.factorized(M) 

100 else: # factorized doesn't pass permc_spec 

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

102 

103 else: 

104 if lstsq: # sometimes necessary as solution is approached 

105 def solve(r): 

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

107 elif cholesky: 

108 L = sp.linalg.cho_factor(M) 

109 

110 def solve(r): 

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

112 else: 

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

114 # with multiple right hand sides is much faster 

115 def solve(r, sym_pos=sym_pos): 

116 return sp.linalg.solve(M, r, sym_pos=sym_pos) 

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

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

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

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

121 except KeyboardInterrupt: 

122 raise 

123 except Exception: 

124 return None 

125 return solve 

126 

127 

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

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

130 permc_spec='MMD_AT_PLUS_A'): 

131 """ 

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

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

134 algorithmic parameters ``gamma and ``eta; 

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

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

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

138 

139 Parameters 

140 ---------- 

141 As defined in [4], except: 

142 sparse : bool 

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

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

145 lstsq : bool 

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

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

148 needed as the solution is approached. 

149 sym_pos : bool 

150 True if the system matrix is symmetric positive definite 

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

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

153 numerical difficulties. 

154 cholesky : bool 

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

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

157 small or prone to numerical difficulties. 

158 pc : bool 

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

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

161 the solution of an additional linear system, the factorization 

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

163 number of algorithm iterations is typically reduced. 

164 ip : bool 

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

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

167 permc_spec : str (default = 'MMD_AT_PLUS_A') 

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

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

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

171 sparsity preservation. Acceptable values are: 

172 

173 - ``NATURAL``: natural ordering. 

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

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

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

177 

178 This option can impact the convergence of the 

179 interior point algorithm; test different values to determine which 

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

181 ``scipy.sparse.linalg.splu``. 

182 

183 Returns 

184 ------- 

185 Search directions as defined in [4] 

186 

187 References 

188 ---------- 

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

190 optimizer for linear programming: an implementation of the 

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

192 2000. 197-232. 

193 

194 """ 

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

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

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

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

199 n_x = len(x) 

200 

201 # [4] Equation 8.8 

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

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

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

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

206 

207 # Assemble M from [4] Equation 8.31 

208 Dinv = x / z 

209 

210 if sparse: 

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

212 else: 

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

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

215 

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

217 # In development this option could be turned off 

218 # but it always seems to improve performance substantially 

219 n_corrections = 1 if pc else 0 

220 

221 i = 0 

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

223 while i <= n_corrections: 

224 # Reference [4] Eq. 8.6 

225 rhatp = eta(gamma) * r_P 

226 rhatd = eta(gamma) * r_D 

227 rhatg = eta(gamma) * r_G 

228 

229 # Reference [4] Eq. 8.7 

230 rhatxs = gamma * mu - x * z 

231 rhattk = gamma * mu - tau * kappa 

232 

233 if i == 1: 

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

235 # Reference [4] Eq. 8.23 

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

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

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

239 tau * kappa - 

240 alpha**2 * d_tau * d_kappa) 

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

242 # Reference [4] Eq. 8.13 

243 rhatxs -= d_x * d_z 

244 rhattk -= d_tau * d_kappa 

245 

246 # sometimes numerical difficulties arise as the solution is approached 

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

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

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

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

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

252 # 4. scipy.linalg.lstsq 

253 # For sparse systems, the order is: 

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

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

256 # 3. scipy.sparse.linalg.splu 

257 # 4. scipy.sparse.linalg.lsqr 

258 solved = False 

259 while(not solved): 

260 try: 

261 # [4] Equation 8.28 

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

263 # [4] Equation 8.29 

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

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

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

267 raise LinAlgError 

268 solved = True 

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

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

271 # there are redundant constraints or when approaching the 

272 # solution. If so, change solver. 

273 if cholesky: 

274 cholesky = False 

275 warn( 

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

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

278 "occasionally, especially as the solution is " 

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

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

281 OptimizeWarning, stacklevel=5) 

282 elif sym_pos: 

283 sym_pos = False 

284 warn( 

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

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

287 "occasionally, especially as the solution is " 

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

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

290 OptimizeWarning, stacklevel=5) 

291 elif not lstsq: 

292 lstsq = True 

293 warn( 

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

295 "failed. This may happen occasionally, " 

296 "especially as the solution is " 

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

298 "your problem may be numerically challenging. " 

299 "If you cannot improve the formulation, consider " 

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

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

302 OptimizeWarning, stacklevel=5) 

303 else: 

304 raise e 

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

306 cholesky, permc_spec) 

307 # [4] Results after 8.29 

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

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

310 d_x = u + p * d_tau 

311 d_y = v + q * d_tau 

312 

313 # [4] Relations between after 8.25 and 8.26 

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

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

316 

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

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

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

320 gamma = 10 

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

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

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

324 i += 1 

325 

326 return d_x, d_y, d_z, d_tau, d_kappa 

327 

328 

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

330 """ 

331 An implementation of [4] equation 8.31 and 8.32 

332 

333 References 

334 ---------- 

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

336 optimizer for linear programming: an implementation of the 

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

338 2000. 197-232. 

339 

340 """ 

341 # [4] 8.31 

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

343 v = solve(r) 

344 # [4] 8.32 

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

346 return u, v 

347 

348 

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

350 """ 

351 An implementation of [4] equation 8.21 

352 

353 References 

354 ---------- 

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

356 optimizer for linear programming: an implementation of the 

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

358 2000. 197-232. 

359 

360 """ 

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

362 # same step is taken in primal and dual spaces 

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

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

365 i_x = d_x < 0 

366 i_z = d_z < 0 

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

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

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

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

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

372 return alpha 

373 

374 

375def _get_message(status): 

376 """ 

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

378 

379 Parameters 

380 ---------- 

381 status : int 

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

383 

384 0 : Optimization terminated successfully 

385 1 : Iteration limit reached 

386 2 : Problem appears to be infeasible 

387 3 : Problem appears to be unbounded 

388 4 : Serious numerical difficulties encountered 

389 

390 Returns 

391 ------- 

392 message : str 

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

394 

395 """ 

396 messages = ( 

397 ["Optimization terminated successfully.", 

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

399 "The algorithm terminated successfully and determined that the " 

400 "problem is infeasible.", 

401 "The algorithm terminated successfully and determined that the " 

402 "problem is unbounded.", 

403 "Numerical difficulties were encountered before the problem " 

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

405 "independence of linear equality constraints, and reasonable " 

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

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

408 ]) 

409 return messages[status] 

410 

411 

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

413 """ 

414 An implementation of [4] Equation 8.9 

415 

416 References 

417 ---------- 

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

419 optimizer for linear programming: an implementation of the 

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

421 2000. 197-232. 

422 

423 """ 

424 x = x + alpha * d_x 

425 tau = tau + alpha * d_tau 

426 z = z + alpha * d_z 

427 kappa = kappa + alpha * d_kappa 

428 y = y + alpha * d_y 

429 return x, y, z, tau, kappa 

430 

431 

432def _get_blind_start(shape): 

433 """ 

434 Return the starting point from [4] 4.4 

435 

436 References 

437 ---------- 

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

439 optimizer for linear programming: an implementation of the 

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

441 2000. 197-232. 

442 

443 """ 

444 m, n = shape 

445 x0 = np.ones(n) 

446 y0 = np.zeros(m) 

447 z0 = np.ones(n) 

448 tau0 = 1 

449 kappa0 = 1 

450 return x0, y0, z0, tau0, kappa0 

451 

452 

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

454 """ 

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

456 the status of optimization. 

457 

458 References 

459 ---------- 

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

461 optimizer for linear programming: an implementation of the 

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

463 2000. 197-232. 

464 

465 """ 

466 

467 # residuals for termination are relative to initial values 

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

469 

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

471 def r_p(x, tau): 

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

473 

474 def r_d(y, z, tau): 

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

476 

477 def r_g(x, y, kappa): 

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

479 

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

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

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

483 

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

485 

486 def norm(a): 

487 return np.linalg.norm(a) 

488 

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

490 r_p0 = r_p(x0, tau0) 

491 r_d0 = r_d(y0, z0, tau0) 

492 r_g0 = r_g(x0, y0, kappa0) 

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

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

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

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

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

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

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

500 

501 

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

503 """ 

504 Print indicators of optimization status to the console. 

505 

506 Parameters 

507 ---------- 

508 rho_p : float 

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

510 rho_d : float 

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

512 rho_g : float 

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

514 alpha : float 

515 The step size, see [4] 4.3 

516 rho_mu : float 

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

518 obj : float 

519 The objective function value of the current iterate 

520 header : bool 

521 True if a header is to be printed 

522 

523 References 

524 ---------- 

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

526 optimizer for linear programming: an implementation of the 

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

528 2000. 197-232. 

529 

530 """ 

531 if header: 

532 print("Primal Feasibility ", 

533 "Dual Feasibility ", 

534 "Duality Gap ", 

535 "Step ", 

536 "Path Parameter ", 

537 "Objective ") 

538 

539 # no clue why this works 

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

541 print(fmt.format( 

542 float(rho_p), 

543 float(rho_d), 

544 float(rho_g), 

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

546 float(rho_mu), 

547 float(obj))) 

548 

549 

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

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

552 r""" 

553 Solve a linear programming problem in standard form: 

554 

555 Minimize:: 

556 

557 c @ x 

558 

559 Subject to:: 

560 

561 A @ x == b 

562 x >= 0 

563 

564 using the interior point method of [4]. 

565 

566 Parameters 

567 ---------- 

568 A : 2-D array 

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

570 constraints at ``x``. 

571 b : 1-D array 

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

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

574 c : 1-D array 

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

576 standard form problem). 

577 c0 : float 

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

579 variables. (Purely for display.) 

580 alpha0 : float 

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

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

583 beta : float 

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

585 maxiter : int 

586 The maximum number of iterations of the algorithm. 

587 disp : bool 

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

589 to the console each iteration. 

590 tol : float 

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

592 sparse : bool 

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

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

595 (dense) arrays rather than sparse matrices. 

596 lstsq : bool 

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

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

599 numerical difficulties are frequently encountered, and a better option 

600 would be to improve the formulation of the problem. 

601 sym_pos : bool 

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

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

604 cholesky : bool 

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

606 Cholesky decomposition followed by explicit forward/backward 

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

608 that are numerically well-behaved. 

609 pc : bool 

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

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

612 ip : bool 

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

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

615 permc_spec : str (default = 'MMD_AT_PLUS_A') 

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

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

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

619 sparsity preservation. Acceptable values are: 

620 

621 - ``NATURAL``: natural ordering. 

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

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

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

625 

626 This option can impact the convergence of the 

627 interior point algorithm; test different values to determine which 

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

629 ``scipy.sparse.linalg.splu``. 

630 callback : callable, optional 

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

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

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

634 

635 x : 1-D array 

636 Current solution vector 

637 fun : float 

638 Current value of the objective function 

639 success : bool 

640 True only when an algorithm has completed successfully, 

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

642 only while the algorithm is still iterating. 

643 slack : 1-D array 

644 The values of the slack variables. Each slack variable 

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

646 the corresponding constraint is active. 

647 con : 1-D array 

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

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

650 phase : int 

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

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

653 status : int 

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

655 status is detected, the algorithm terminates. 

656 nit : int 

657 The number of iterations performed. 

658 message : str 

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

660 postsolve_args : tuple 

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

662 problem into the solution to the original problem. 

663 

664 Returns 

665 ------- 

666 x_hat : float 

667 Solution vector (for standard form problem). 

668 status : int 

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

670 

671 0 : Optimization terminated successfully 

672 1 : Iteration limit reached 

673 2 : Problem appears to be infeasible 

674 3 : Problem appears to be unbounded 

675 4 : Serious numerical difficulties encountered 

676 

677 message : str 

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

679 iteration : int 

680 The number of iterations taken to solve the problem 

681 

682 References 

683 ---------- 

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

685 optimizer for linear programming: an implementation of the 

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

687 2000. 197-232. 

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

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

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

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

692 

693 """ 

694 

695 iteration = 0 

696 

697 # default initial point 

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

699 

700 # first iteration is special improvement of initial point 

701 ip = ip if pc else False 

702 

703 # [4] 4.5 

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

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

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

707 

708 if disp: 

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

710 if callback is not None: 

711 x_o, fun, slack, con, _ = _postsolve(x/tau, postsolve_args, 

712 tol=tol) 

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

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

715 'complete': False, 'status': 0, 

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

717 callback(res) 

718 

719 status = 0 

720 message = "Optimization terminated successfully." 

721 

722 if sparse: 

723 A = sps.csc_matrix(A) 

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

725 # Redefine it to avoid calculating again 

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

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 tol=tol) 

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

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

796 'complete': False, 'status': 0, 

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

798 callback(res) 

799 

800 # [4] 4.5 

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

802 max(1, kappa)) 

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

804 if inf1 or inf2: 

805 # [4] Lemma 8.4 / Theorem 8.3 

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

807 status = 2 

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

809 status = 3 

810 message = _get_message(status) 

811 break 

812 elif iteration >= maxiter: 

813 status = 1 

814 message = _get_message(status) 

815 break 

816 

817 x_hat = x / tau 

818 # [4] Statement after Theorem 8.2 

819 return x_hat, status, message, iteration 

820 

821 

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

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

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

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

826 r""" 

827 Minimize a linear objective function subject to linear 

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

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

830 of the following form: 

831 

832 Minimize:: 

833 

834 c @ x 

835 

836 Subject to:: 

837 

838 A @ x == b 

839 x >= 0 

840 

841 Parameters 

842 ---------- 

843 c : 1-D array 

844 Coefficients of the linear objective function to be minimized. 

845 c0 : float 

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

847 variables. (Purely for display.) 

848 A : 2-D array 

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

850 constraints at ``x``. 

851 b : 1-D array 

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

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

854 callback : callable, optional 

855 Callback function to be executed once per iteration. 

856 postsolve_args : tuple 

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

858 problem into the solution to the original problem. 

859 

860 Options 

861 ------- 

862 maxiter : int (default = 1000) 

863 The maximum number of iterations of the algorithm. 

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

865 Termination tolerance to be used for all termination criteria; 

866 see [4]_ Section 4.5. 

867 disp : bool (default = False) 

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

869 to the console each iteration. 

870 alpha0 : float (default = 0.99995) 

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

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

873 beta : float (default = 0.1) 

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

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

876 sparse : bool (default = False) 

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

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

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

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

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

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

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

884 lstsq : bool (default = False) 

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

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

887 numerical difficulties are encountered. Leave this at the default 

888 unless you receive a warning message suggesting otherwise. 

889 sym_pos : bool (default = True) 

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

891 symmetric positive definite normal equation matrix 

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

893 a warning message suggesting otherwise. 

894 cholesky : bool (default = True) 

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

896 Cholesky decomposition followed by explicit forward/backward 

897 substitution. This is typically faster for problems 

898 that are numerically well-behaved. 

899 pc : bool (default = True) 

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

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

902 ip : bool (default = False) 

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

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

905 depends on the problem. 

906 permc_spec : str (default = 'MMD_AT_PLUS_A') 

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

908 True``, and no SuiteSparse.) 

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

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

911 sparsity preservation. Acceptable values are: 

912 

913 - ``NATURAL``: natural ordering. 

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

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

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

917 

918 This option can impact the convergence of the 

919 interior point algorithm; test different values to determine which 

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

921 ``scipy.sparse.linalg.splu``. 

922 unknown_options : dict 

923 Optional arguments not used by this particular solver. If 

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

925 unused options. 

926 

927 Returns 

928 ------- 

929 x : 1-D array 

930 Solution vector. 

931 status : int 

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

933 

934 0 : Optimization terminated successfully 

935 1 : Iteration limit reached 

936 2 : Problem appears to be infeasible 

937 3 : Problem appears to be unbounded 

938 4 : Serious numerical difficulties encountered 

939 

940 message : str 

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

942 iteration : int 

943 The number of iterations taken to solve the problem. 

944 

945 Notes 

946 ----- 

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

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

949 

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

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

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

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

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

955 formulation, which provides certificates of infeasibility or unboundedness 

956 where applicable. 

957 

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

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

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

961 point can be calculated according to the additional recommendations of 

962 [4]_ Section 4.4. 

963 

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

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

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

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

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

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

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

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

972 matrices involved are symmetric positive definite, so Cholesky 

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

974 

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

976 on third-party software availability and the conditioning of the problem. 

977 

978 For dense problems, solvers are tried in the following order: 

979 

980 1. ``scipy.linalg.cho_factor`` 

981 

982 2. ``scipy.linalg.solve`` with option ``sym_pos=True`` 

983 

984 3. ``scipy.linalg.solve`` with option ``sym_pos=False`` 

985 

986 4. ``scipy.linalg.lstsq`` 

987 

988 For sparse problems: 

989 

990 1. ``sksparse.cholmod.cholesky`` (if scikit-sparse and SuiteSparse are installed) 

991 

992 2. ``scipy.sparse.linalg.factorized`` (if scikit-umfpack and SuiteSparse are installed) 

993 

994 3. ``scipy.sparse.linalg.splu`` (which uses SuperLU distributed with SciPy) 

995 

996 4. ``scipy.sparse.linalg.lsqr`` 

997 

998 If the solver fails for any reason, successively more robust (but slower) 

999 solvers are attempted in the order indicated. Attempting, failing, and 

1000 re-starting factorization can be time consuming, so if the problem is 

1001 numerically challenging, options can be set to bypass solvers that are 

1002 failing. Setting ``cholesky=False`` skips to solver 2, 

1003 ``sym_pos=False`` skips to solver 3, and ``lstsq=True`` skips 

1004 to solver 4 for both sparse and dense problems. 

1005 

1006 Potential improvements for combatting issues associated with dense 

1007 columns in otherwise sparse problems are outlined in [4]_ Section 5.3 and 

1008 [10]_ Section 4.1-4.2; the latter also discusses the alleviation of 

1009 accuracy issues associated with the substitution approach to free 

1010 variables. 

1011 

1012 After calculating the search direction, the maximum possible step size 

1013 that does not activate the non-negativity constraints is calculated, and 

1014 the smaller of this step size and unity is applied (as in [4]_ Section 

1015 4.1.) [4]_ Section 4.3 suggests improvements for choosing the step size. 

1016 

1017 The new point is tested according to the termination conditions of [4]_ 

1018 Section 4.5. The same tolerance, which can be set using the ``tol`` option, 

1019 is used for all checks. (A potential improvement would be to expose 

1020 the different tolerances to be set independently.) If optimality, 

1021 unboundedness, or infeasibility is detected, the solve procedure 

1022 terminates; otherwise it repeats. 

1023 

1024 The expected problem formulation differs between the top level ``linprog`` 

1025 module and the method specific solvers. The method specific solvers expect a 

1026 problem in standard form: 

1027 

1028 Minimize:: 

1029 

1030 c @ x 

1031 

1032 Subject to:: 

1033 

1034 A @ x == b 

1035 x >= 0 

1036 

1037 Whereas the top level ``linprog`` module expects a problem of form: 

1038 

1039 Minimize:: 

1040 

1041 c @ x 

1042 

1043 Subject to:: 

1044 

1045 A_ub @ x <= b_ub 

1046 A_eq @ x == b_eq 

1047 lb <= x <= ub 

1048 

1049 where ``lb = 0`` and ``ub = None`` unless set in ``bounds``. 

1050 

1051 The original problem contains equality, upper-bound and variable constraints 

1052 whereas the method specific solver requires equality constraints and 

1053 variable non-negativity. 

1054 

1055 ``linprog`` module converts the original problem to standard form by 

1056 converting the simple bounds to upper bound constraints, introducing 

1057 non-negative slack variables for inequality constraints, and expressing 

1058 unbounded variables as the difference between two non-negative variables. 

1059 

1060 

1061 References 

1062 ---------- 

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

1064 optimizer for linear programming: an implementation of the 

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

1066 2000. 197-232. 

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

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

1069 March 2004. Available 2/25/2017 at 

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

1071 .. [8] Andersen, Erling D., and Knud D. Andersen. "Presolving in linear 

1072 programming." Mathematical Programming 71.2 (1995): 221-245. 

1073 .. [9] Bertsimas, Dimitris, and J. Tsitsiklis. "Introduction to linear 

1074 programming." Athena Scientific 1 (1997): 997. 

1075 .. [10] Andersen, Erling D., et al. Implementation of interior point methods 

1076 for large scale linear programming. HEC/Universite de Geneve, 1996. 

1077 

1078 """ 

1079 

1080 _check_unknown_options(unknown_options) 

1081 

1082 # These should be warnings, not errors 

1083 if (cholesky or cholesky is None) and sparse and not has_cholmod: 

1084 if cholesky: 

1085 warn("Sparse cholesky is only available with scikit-sparse. " 

1086 "Setting `cholesky = False`", 

1087 OptimizeWarning, stacklevel=3) 

1088 cholesky = False 

1089 

1090 if sparse and lstsq: 

1091 warn("Option combination 'sparse':True and 'lstsq':True " 

1092 "is not recommended.", 

1093 OptimizeWarning, stacklevel=3) 

1094 

1095 if lstsq and cholesky: 

1096 warn("Invalid option combination 'lstsq':True " 

1097 "and 'cholesky':True; option 'cholesky' has no effect when " 

1098 "'lstsq' is set True.", 

1099 OptimizeWarning, stacklevel=3) 

1100 

1101 valid_permc_spec = ('NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', 'COLAMD') 

1102 if permc_spec.upper() not in valid_permc_spec: 

1103 warn("Invalid permc_spec option: '" + str(permc_spec) + "'. " 

1104 "Acceptable values are 'NATURAL', 'MMD_ATA', 'MMD_AT_PLUS_A', " 

1105 "and 'COLAMD'. Reverting to default.", 

1106 OptimizeWarning, stacklevel=3) 

1107 permc_spec = 'MMD_AT_PLUS_A' 

1108 

1109 # This can be an error 

1110 if not sym_pos and cholesky: 

1111 raise ValueError( 

1112 "Invalid option combination 'sym_pos':False " 

1113 "and 'cholesky':True: Cholesky decomposition is only possible " 

1114 "for symmetric positive definite matrices.") 

1115 

1116 cholesky = cholesky or (cholesky is None and sym_pos and not lstsq) 

1117 

1118 x, status, message, iteration = _ip_hsd(A, b, c, c0, alpha0, beta, 

1119 maxiter, disp, tol, sparse, 

1120 lstsq, sym_pos, cholesky, 

1121 pc, ip, permc_spec, callback, 

1122 postsolve_args) 

1123 

1124 return x, status, message, iteration