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

1r""" 

2 

3Nonlinear solvers 

4----------------- 

5 

6.. currentmodule:: scipy.optimize 

7 

8This is a collection of general-purpose nonlinear multidimensional 

9solvers. These solvers find *x* for which *F(x) = 0*. Both *x* 

10and *F* can be multidimensional. 

11 

12Routines 

13~~~~~~~~ 

14 

15Large-scale nonlinear solvers: 

16 

17.. autosummary:: 

18 

19 newton_krylov 

20 anderson 

21 

22General nonlinear solvers: 

23 

24.. autosummary:: 

25 

26 broyden1 

27 broyden2 

28 

29Simple iterations: 

30 

31.. autosummary:: 

32 

33 excitingmixing 

34 linearmixing 

35 diagbroyden 

36 

37 

38Examples 

39~~~~~~~~ 

40 

41**Small problem** 

42 

43>>> def F(x): 

44... return np.cos(x) + x[::-1] - [1, 2, 3, 4] 

45>>> import scipy.optimize 

46>>> x = scipy.optimize.broyden1(F, [1,1,1,1], f_tol=1e-14) 

47>>> x 

48array([ 4.04674914, 3.91158389, 2.71791677, 1.61756251]) 

49>>> np.cos(x) + x[::-1] 

50array([ 1., 2., 3., 4.]) 

51 

52 

53**Large problem** 

54 

55Suppose that we needed to solve the following integrodifferential 

56equation on the square :math:`[0,1]\times[0,1]`: 

57 

58.. math:: 

59 

60 \nabla^2 P = 10 \left(\int_0^1\int_0^1\cosh(P)\,dx\,dy\right)^2 

61 

62with :math:`P(x,1) = 1` and :math:`P=0` elsewhere on the boundary of 

63the square. 

64 

65The solution can be found using the `newton_krylov` solver: 

66 

67.. plot:: 

68 

69 import numpy as np 

70 from scipy.optimize import newton_krylov 

71 from numpy import cosh, zeros_like, mgrid, zeros 

72 

73 # parameters 

74 nx, ny = 75, 75 

75 hx, hy = 1./(nx-1), 1./(ny-1) 

76 

77 P_left, P_right = 0, 0 

78 P_top, P_bottom = 1, 0 

79 

80 def residual(P): 

81 d2x = zeros_like(P) 

82 d2y = zeros_like(P) 

83 

84 d2x[1:-1] = (P[2:] - 2*P[1:-1] + P[:-2]) / hx/hx 

85 d2x[0] = (P[1] - 2*P[0] + P_left)/hx/hx 

86 d2x[-1] = (P_right - 2*P[-1] + P[-2])/hx/hx 

87 

88 d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy 

89 d2y[:,0] = (P[:,1] - 2*P[:,0] + P_bottom)/hy/hy 

90 d2y[:,-1] = (P_top - 2*P[:,-1] + P[:,-2])/hy/hy 

91 

92 return d2x + d2y - 10*cosh(P).mean()**2 

93 

94 # solve 

95 guess = zeros((nx, ny), float) 

96 sol = newton_krylov(residual, guess, method='lgmres', verbose=1) 

97 print('Residual: %g' % abs(residual(sol)).max()) 

98 

99 # visualize 

100 import matplotlib.pyplot as plt 

101 x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)] 

102 plt.pcolormesh(x, y, sol, shading='gouraud') 

103 plt.colorbar() 

104 plt.show() 

105 

106""" 

107# Copyright (C) 2009, Pauli Virtanen <pav@iki.fi> 

108# Distributed under the same license as SciPy. 

109 

110import sys 

111import numpy as np 

112from scipy.linalg import norm, solve, inv, qr, svd, LinAlgError 

113from numpy import asarray, dot, vdot 

114import scipy.sparse.linalg 

115import scipy.sparse 

116from scipy.linalg import get_blas_funcs 

117import inspect 

118from scipy._lib._util import getfullargspec_no_self as _getfullargspec 

119from .linesearch import scalar_search_wolfe1, scalar_search_armijo 

120 

121 

122__all__ = [ 

123 'broyden1', 'broyden2', 'anderson', 'linearmixing', 

124 'diagbroyden', 'excitingmixing', 'newton_krylov'] 

125 

126#------------------------------------------------------------------------------ 

127# Utility functions 

128#------------------------------------------------------------------------------ 

129 

130 

131class NoConvergence(Exception): 

132 pass 

133 

134 

135def maxnorm(x): 

136 return np.absolute(x).max() 

137 

138 

139def _as_inexact(x): 

140 """Return `x` as an array, of either floats or complex floats""" 

141 x = asarray(x) 

142 if not np.issubdtype(x.dtype, np.inexact): 

143 return asarray(x, dtype=np.float_) 

144 return x 

145 

146 

147def _array_like(x, x0): 

148 """Return ndarray `x` as same array subclass and shape as `x0`""" 

149 x = np.reshape(x, np.shape(x0)) 

150 wrap = getattr(x0, '__array_wrap__', x.__array_wrap__) 

151 return wrap(x) 

152 

153 

154def _safe_norm(v): 

155 if not np.isfinite(v).all(): 

156 return np.array(np.inf) 

157 return norm(v) 

158 

159#------------------------------------------------------------------------------ 

160# Generic nonlinear solver machinery 

161#------------------------------------------------------------------------------ 

162 

163 

164_doc_parts = dict( 

165 params_basic=""" 

166 F : function(x) -> f 

167 Function whose root to find; should take and return an array-like 

168 object. 

169 xin : array_like 

170 Initial guess for the solution 

171 """.strip(), 

172 params_extra=""" 

173 iter : int, optional 

174 Number of iterations to make. If omitted (default), make as many 

175 as required to meet tolerances. 

176 verbose : bool, optional 

177 Print status to stdout on every iteration. 

178 maxiter : int, optional 

179 Maximum number of iterations to make. If more are needed to 

180 meet convergence, `NoConvergence` is raised. 

181 f_tol : float, optional 

182 Absolute tolerance (in max-norm) for the residual. 

183 If omitted, default is 6e-6. 

184 f_rtol : float, optional 

185 Relative tolerance for the residual. If omitted, not used. 

186 x_tol : float, optional 

187 Absolute minimum step size, as determined from the Jacobian 

188 approximation. If the step size is smaller than this, optimization 

189 is terminated as successful. If omitted, not used. 

190 x_rtol : float, optional 

191 Relative minimum step size. If omitted, not used. 

192 tol_norm : function(vector) -> scalar, optional 

193 Norm to use in convergence check. Default is the maximum norm. 

194 line_search : {None, 'armijo' (default), 'wolfe'}, optional 

195 Which type of a line search to use to determine the step size in the 

196 direction given by the Jacobian approximation. Defaults to 'armijo'. 

197 callback : function, optional 

198 Optional callback function. It is called on every iteration as 

199 ``callback(x, f)`` where `x` is the current solution and `f` 

200 the corresponding residual. 

201 

202 Returns 

203 ------- 

204 sol : ndarray 

205 An array (of similar array type as `x0`) containing the final solution. 

206 

207 Raises 

208 ------ 

209 NoConvergence 

210 When a solution was not found. 

211 

212 """.strip() 

213) 

214 

215 

216def _set_doc(obj): 

217 if obj.__doc__: 

218 obj.__doc__ = obj.__doc__ % _doc_parts 

219 

220 

221def nonlin_solve(F, x0, jacobian='krylov', iter=None, verbose=False, 

222 maxiter=None, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None, 

223 tol_norm=None, line_search='armijo', callback=None, 

224 full_output=False, raise_exception=True): 

225 """ 

226 Find a root of a function, in a way suitable for large-scale problems. 

227 

228 Parameters 

229 ---------- 

230 %(params_basic)s 

231 jacobian : Jacobian 

232 A Jacobian approximation: `Jacobian` object or something that 

233 `asjacobian` can transform to one. Alternatively, a string specifying 

234 which of the builtin Jacobian approximations to use: 

235 

236 krylov, broyden1, broyden2, anderson 

237 diagbroyden, linearmixing, excitingmixing 

238 

239 %(params_extra)s 

240 full_output : bool 

241 If true, returns a dictionary `info` containing convergence 

242 information. 

243 raise_exception : bool 

244 If True, a `NoConvergence` exception is raise if no solution is found. 

245 

246 See Also 

247 -------- 

248 asjacobian, Jacobian 

249 

250 Notes 

251 ----- 

252 This algorithm implements the inexact Newton method, with 

253 backtracking or full line searches. Several Jacobian 

254 approximations are available, including Krylov and Quasi-Newton 

255 methods. 

256 

257 References 

258 ---------- 

259 .. [KIM] C. T. Kelley, \"Iterative Methods for Linear and Nonlinear 

260 Equations\". Society for Industrial and Applied Mathematics. (1995) 

261 https://archive.siam.org/books/kelley/fr16/ 

262 

263 """ 

264 # Can't use default parameters because it's being explicitly passed as None 

265 # from the calling function, so we need to set it here. 

266 tol_norm = maxnorm if tol_norm is None else tol_norm 

267 condition = TerminationCondition(f_tol=f_tol, f_rtol=f_rtol, 

268 x_tol=x_tol, x_rtol=x_rtol, 

269 iter=iter, norm=tol_norm) 

270 

271 x0 = _as_inexact(x0) 

272 func = lambda z: _as_inexact(F(_array_like(z, x0))).flatten() 

273 x = x0.flatten() 

274 

275 dx = np.inf 

276 Fx = func(x) 

277 Fx_norm = norm(Fx) 

278 

279 jacobian = asjacobian(jacobian) 

280 jacobian.setup(x.copy(), Fx, func) 

281 

282 if maxiter is None: 

283 if iter is not None: 

284 maxiter = iter + 1 

285 else: 

286 maxiter = 100*(x.size+1) 

287 

288 if line_search is True: 

289 line_search = 'armijo' 

290 elif line_search is False: 

291 line_search = None 

292 

293 if line_search not in (None, 'armijo', 'wolfe'): 

294 raise ValueError("Invalid line search") 

295 

296 # Solver tolerance selection 

297 gamma = 0.9 

298 eta_max = 0.9999 

299 eta_treshold = 0.1 

300 eta = 1e-3 

301 

302 for n in range(maxiter): 

303 status = condition.check(Fx, x, dx) 

304 if status: 

305 break 

306 

307 # The tolerance, as computed for scipy.sparse.linalg.* routines 

308 tol = min(eta, eta*Fx_norm) 

309 dx = -jacobian.solve(Fx, tol=tol) 

310 

311 if norm(dx) == 0: 

312 raise ValueError("Jacobian inversion yielded zero vector. " 

313 "This indicates a bug in the Jacobian " 

314 "approximation.") 

315 

316 # Line search, or Newton step 

317 if line_search: 

318 s, x, Fx, Fx_norm_new = _nonlin_line_search(func, x, Fx, dx, 

319 line_search) 

320 else: 

321 s = 1.0 

322 x = x + dx 

323 Fx = func(x) 

324 Fx_norm_new = norm(Fx) 

325 

326 jacobian.update(x.copy(), Fx) 

327 

328 if callback: 

329 callback(x, Fx) 

330 

331 # Adjust forcing parameters for inexact methods 

332 eta_A = gamma * Fx_norm_new**2 / Fx_norm**2 

333 if gamma * eta**2 < eta_treshold: 

334 eta = min(eta_max, eta_A) 

335 else: 

336 eta = min(eta_max, max(eta_A, gamma*eta**2)) 

337 

338 Fx_norm = Fx_norm_new 

339 

340 # Print status 

341 if verbose: 

342 sys.stdout.write("%d: |F(x)| = %g; step %g\n" % ( 

343 n, tol_norm(Fx), s)) 

344 sys.stdout.flush() 

345 else: 

346 if raise_exception: 

347 raise NoConvergence(_array_like(x, x0)) 

348 else: 

349 status = 2 

350 

351 if full_output: 

352 info = {'nit': condition.iteration, 

353 'fun': Fx, 

354 'status': status, 

355 'success': status == 1, 

356 'message': {1: 'A solution was found at the specified ' 

357 'tolerance.', 

358 2: 'The maximum number of iterations allowed ' 

359 'has been reached.' 

360 }[status] 

361 } 

362 return _array_like(x, x0), info 

363 else: 

364 return _array_like(x, x0) 

365 

366 

367_set_doc(nonlin_solve) 

368 

369 

370def _nonlin_line_search(func, x, Fx, dx, search_type='armijo', rdiff=1e-8, 

371 smin=1e-2): 

372 tmp_s = [0] 

373 tmp_Fx = [Fx] 

374 tmp_phi = [norm(Fx)**2] 

375 s_norm = norm(x) / norm(dx) 

376 

377 def phi(s, store=True): 

378 if s == tmp_s[0]: 

379 return tmp_phi[0] 

380 xt = x + s*dx 

381 v = func(xt) 

382 p = _safe_norm(v)**2 

383 if store: 

384 tmp_s[0] = s 

385 tmp_phi[0] = p 

386 tmp_Fx[0] = v 

387 return p 

388 

389 def derphi(s): 

390 ds = (abs(s) + s_norm + 1) * rdiff 

391 return (phi(s+ds, store=False) - phi(s)) / ds 

392 

393 if search_type == 'wolfe': 

394 s, phi1, phi0 = scalar_search_wolfe1(phi, derphi, tmp_phi[0], 

395 xtol=1e-2, amin=smin) 

396 elif search_type == 'armijo': 

397 s, phi1 = scalar_search_armijo(phi, tmp_phi[0], -tmp_phi[0], 

398 amin=smin) 

399 

400 if s is None: 

401 # XXX: No suitable step length found. Take the full Newton step, 

402 # and hope for the best. 

403 s = 1.0 

404 

405 x = x + s*dx 

406 if s == tmp_s[0]: 

407 Fx = tmp_Fx[0] 

408 else: 

409 Fx = func(x) 

410 Fx_norm = norm(Fx) 

411 

412 return s, x, Fx, Fx_norm 

413 

414 

415class TerminationCondition(object): 

416 """ 

417 Termination condition for an iteration. It is terminated if 

418 

419 - |F| < f_rtol*|F_0|, AND 

420 - |F| < f_tol 

421 

422 AND 

423 

424 - |dx| < x_rtol*|x|, AND 

425 - |dx| < x_tol 

426 

427 """ 

428 def __init__(self, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None, 

429 iter=None, norm=maxnorm): 

430 

431 if f_tol is None: 

432 f_tol = np.finfo(np.float_).eps ** (1./3) 

433 if f_rtol is None: 

434 f_rtol = np.inf 

435 if x_tol is None: 

436 x_tol = np.inf 

437 if x_rtol is None: 

438 x_rtol = np.inf 

439 

440 self.x_tol = x_tol 

441 self.x_rtol = x_rtol 

442 self.f_tol = f_tol 

443 self.f_rtol = f_rtol 

444 

445 self.norm = norm 

446 

447 self.iter = iter 

448 

449 self.f0_norm = None 

450 self.iteration = 0 

451 

452 def check(self, f, x, dx): 

453 self.iteration += 1 

454 f_norm = self.norm(f) 

455 x_norm = self.norm(x) 

456 dx_norm = self.norm(dx) 

457 

458 if self.f0_norm is None: 

459 self.f0_norm = f_norm 

460 

461 if f_norm == 0: 

462 return 1 

463 

464 if self.iter is not None: 

465 # backwards compatibility with SciPy 0.6.0 

466 return 2 * (self.iteration > self.iter) 

467 

468 # NB: condition must succeed for rtol=inf even if norm == 0 

469 return int((f_norm <= self.f_tol 

470 and f_norm/self.f_rtol <= self.f0_norm) 

471 and (dx_norm <= self.x_tol 

472 and dx_norm/self.x_rtol <= x_norm)) 

473 

474 

475#------------------------------------------------------------------------------ 

476# Generic Jacobian approximation 

477#------------------------------------------------------------------------------ 

478 

479class Jacobian(object): 

480 """ 

481 Common interface for Jacobians or Jacobian approximations. 

482 

483 The optional methods come useful when implementing trust region 

484 etc., algorithms that often require evaluating transposes of the 

485 Jacobian. 

486 

487 Methods 

488 ------- 

489 solve 

490 Returns J^-1 * v 

491 update 

492 Updates Jacobian to point `x` (where the function has residual `Fx`) 

493 

494 matvec : optional 

495 Returns J * v 

496 rmatvec : optional 

497 Returns A^H * v 

498 rsolve : optional 

499 Returns A^-H * v 

500 matmat : optional 

501 Returns A * V, where V is a dense matrix with dimensions (N,K). 

502 todense : optional 

503 Form the dense Jacobian matrix. Necessary for dense trust region 

504 algorithms, and useful for testing. 

505 

506 Attributes 

507 ---------- 

508 shape 

509 Matrix dimensions (M, N) 

510 dtype 

511 Data type of the matrix. 

512 func : callable, optional 

513 Function the Jacobian corresponds to 

514 

515 """ 

516 

517 def __init__(self, **kw): 

518 names = ["solve", "update", "matvec", "rmatvec", "rsolve", 

519 "matmat", "todense", "shape", "dtype"] 

520 for name, value in kw.items(): 

521 if name not in names: 

522 raise ValueError("Unknown keyword argument %s" % name) 

523 if value is not None: 

524 setattr(self, name, kw[name]) 

525 

526 if hasattr(self, 'todense'): 

527 self.__array__ = lambda: self.todense() 

528 

529 def aspreconditioner(self): 

530 return InverseJacobian(self) 

531 

532 def solve(self, v, tol=0): 

533 raise NotImplementedError 

534 

535 def update(self, x, F): 

536 pass 

537 

538 def setup(self, x, F, func): 

539 self.func = func 

540 self.shape = (F.size, x.size) 

541 self.dtype = F.dtype 

542 if self.__class__.setup is Jacobian.setup: 

543 # Call on the first point unless overridden 

544 self.update(x, F) 

545 

546 

547class InverseJacobian(object): 

548 def __init__(self, jacobian): 

549 self.jacobian = jacobian 

550 self.matvec = jacobian.solve 

551 self.update = jacobian.update 

552 if hasattr(jacobian, 'setup'): 

553 self.setup = jacobian.setup 

554 if hasattr(jacobian, 'rsolve'): 

555 self.rmatvec = jacobian.rsolve 

556 

557 @property 

558 def shape(self): 

559 return self.jacobian.shape 

560 

561 @property 

562 def dtype(self): 

563 return self.jacobian.dtype 

564 

565 

566def asjacobian(J): 

567 """ 

568 Convert given object to one suitable for use as a Jacobian. 

569 """ 

570 spsolve = scipy.sparse.linalg.spsolve 

571 if isinstance(J, Jacobian): 

572 return J 

573 elif inspect.isclass(J) and issubclass(J, Jacobian): 

574 return J() 

575 elif isinstance(J, np.ndarray): 

576 if J.ndim > 2: 

577 raise ValueError('array must have rank <= 2') 

578 J = np.atleast_2d(np.asarray(J)) 

579 if J.shape[0] != J.shape[1]: 

580 raise ValueError('array must be square') 

581 

582 return Jacobian(matvec=lambda v: dot(J, v), 

583 rmatvec=lambda v: dot(J.conj().T, v), 

584 solve=lambda v: solve(J, v), 

585 rsolve=lambda v: solve(J.conj().T, v), 

586 dtype=J.dtype, shape=J.shape) 

587 elif scipy.sparse.isspmatrix(J): 

588 if J.shape[0] != J.shape[1]: 

589 raise ValueError('matrix must be square') 

590 return Jacobian(matvec=lambda v: J*v, 

591 rmatvec=lambda v: J.conj().T * v, 

592 solve=lambda v: spsolve(J, v), 

593 rsolve=lambda v: spsolve(J.conj().T, v), 

594 dtype=J.dtype, shape=J.shape) 

595 elif hasattr(J, 'shape') and hasattr(J, 'dtype') and hasattr(J, 'solve'): 

596 return Jacobian(matvec=getattr(J, 'matvec'), 

597 rmatvec=getattr(J, 'rmatvec'), 

598 solve=J.solve, 

599 rsolve=getattr(J, 'rsolve'), 

600 update=getattr(J, 'update'), 

601 setup=getattr(J, 'setup'), 

602 dtype=J.dtype, 

603 shape=J.shape) 

604 elif callable(J): 

605 # Assume it's a function J(x) that returns the Jacobian 

606 class Jac(Jacobian): 

607 def update(self, x, F): 

608 self.x = x 

609 

610 def solve(self, v, tol=0): 

611 m = J(self.x) 

612 if isinstance(m, np.ndarray): 

613 return solve(m, v) 

614 elif scipy.sparse.isspmatrix(m): 

615 return spsolve(m, v) 

616 else: 

617 raise ValueError("Unknown matrix type") 

618 

619 def matvec(self, v): 

620 m = J(self.x) 

621 if isinstance(m, np.ndarray): 

622 return dot(m, v) 

623 elif scipy.sparse.isspmatrix(m): 

624 return m*v 

625 else: 

626 raise ValueError("Unknown matrix type") 

627 

628 def rsolve(self, v, tol=0): 

629 m = J(self.x) 

630 if isinstance(m, np.ndarray): 

631 return solve(m.conj().T, v) 

632 elif scipy.sparse.isspmatrix(m): 

633 return spsolve(m.conj().T, v) 

634 else: 

635 raise ValueError("Unknown matrix type") 

636 

637 def rmatvec(self, v): 

638 m = J(self.x) 

639 if isinstance(m, np.ndarray): 

640 return dot(m.conj().T, v) 

641 elif scipy.sparse.isspmatrix(m): 

642 return m.conj().T * v 

643 else: 

644 raise ValueError("Unknown matrix type") 

645 return Jac() 

646 elif isinstance(J, str): 

647 return dict(broyden1=BroydenFirst, 

648 broyden2=BroydenSecond, 

649 anderson=Anderson, 

650 diagbroyden=DiagBroyden, 

651 linearmixing=LinearMixing, 

652 excitingmixing=ExcitingMixing, 

653 krylov=KrylovJacobian)[J]() 

654 else: 

655 raise TypeError('Cannot convert object to a Jacobian') 

656 

657 

658#------------------------------------------------------------------------------ 

659# Broyden 

660#------------------------------------------------------------------------------ 

661 

662class GenericBroyden(Jacobian): 

663 def setup(self, x0, f0, func): 

664 Jacobian.setup(self, x0, f0, func) 

665 self.last_f = f0 

666 self.last_x = x0 

667 

668 if hasattr(self, 'alpha') and self.alpha is None: 

669 # Autoscale the initial Jacobian parameter 

670 # unless we have already guessed the solution. 

671 normf0 = norm(f0) 

672 if normf0: 

673 self.alpha = 0.5*max(norm(x0), 1) / normf0 

674 else: 

675 self.alpha = 1.0 

676 

677 def _update(self, x, f, dx, df, dx_norm, df_norm): 

678 raise NotImplementedError 

679 

680 def update(self, x, f): 

681 df = f - self.last_f 

682 dx = x - self.last_x 

683 self._update(x, f, dx, df, norm(dx), norm(df)) 

684 self.last_f = f 

685 self.last_x = x 

686 

687 

688class LowRankMatrix(object): 

689 r""" 

690 A matrix represented as 

691 

692 .. math:: \alpha I + \sum_{n=0}^{n=M} c_n d_n^\dagger 

693 

694 However, if the rank of the matrix reaches the dimension of the vectors, 

695 full matrix representation will be used thereon. 

696 

697 """ 

698 

699 def __init__(self, alpha, n, dtype): 

700 self.alpha = alpha 

701 self.cs = [] 

702 self.ds = [] 

703 self.n = n 

704 self.dtype = dtype 

705 self.collapsed = None 

706 

707 @staticmethod 

708 def _matvec(v, alpha, cs, ds): 

709 axpy, scal, dotc = get_blas_funcs(['axpy', 'scal', 'dotc'], 

710 cs[:1] + [v]) 

711 w = alpha * v 

712 for c, d in zip(cs, ds): 

713 a = dotc(d, v) 

714 w = axpy(c, w, w.size, a) 

715 return w 

716 

717 @staticmethod 

718 def _solve(v, alpha, cs, ds): 

719 """Evaluate w = M^-1 v""" 

720 if len(cs) == 0: 

721 return v/alpha 

722 

723 # (B + C D^H)^-1 = B^-1 - B^-1 C (I + D^H B^-1 C)^-1 D^H B^-1 

724 

725 axpy, dotc = get_blas_funcs(['axpy', 'dotc'], cs[:1] + [v]) 

726 

727 c0 = cs[0] 

728 A = alpha * np.identity(len(cs), dtype=c0.dtype) 

729 for i, d in enumerate(ds): 

730 for j, c in enumerate(cs): 

731 A[i,j] += dotc(d, c) 

732 

733 q = np.zeros(len(cs), dtype=c0.dtype) 

734 for j, d in enumerate(ds): 

735 q[j] = dotc(d, v) 

736 q /= alpha 

737 q = solve(A, q) 

738 

739 w = v/alpha 

740 for c, qc in zip(cs, q): 

741 w = axpy(c, w, w.size, -qc) 

742 

743 return w 

744 

745 def matvec(self, v): 

746 """Evaluate w = M v""" 

747 if self.collapsed is not None: 

748 return np.dot(self.collapsed, v) 

749 return LowRankMatrix._matvec(v, self.alpha, self.cs, self.ds) 

750 

751 def rmatvec(self, v): 

752 """Evaluate w = M^H v""" 

753 if self.collapsed is not None: 

754 return np.dot(self.collapsed.T.conj(), v) 

755 return LowRankMatrix._matvec(v, np.conj(self.alpha), self.ds, self.cs) 

756 

757 def solve(self, v, tol=0): 

758 """Evaluate w = M^-1 v""" 

759 if self.collapsed is not None: 

760 return solve(self.collapsed, v) 

761 return LowRankMatrix._solve(v, self.alpha, self.cs, self.ds) 

762 

763 def rsolve(self, v, tol=0): 

764 """Evaluate w = M^-H v""" 

765 if self.collapsed is not None: 

766 return solve(self.collapsed.T.conj(), v) 

767 return LowRankMatrix._solve(v, np.conj(self.alpha), self.ds, self.cs) 

768 

769 def append(self, c, d): 

770 if self.collapsed is not None: 

771 self.collapsed += c[:,None] * d[None,:].conj() 

772 return 

773 

774 self.cs.append(c) 

775 self.ds.append(d) 

776 

777 if len(self.cs) > c.size: 

778 self.collapse() 

779 

780 def __array__(self): 

781 if self.collapsed is not None: 

782 return self.collapsed 

783 

784 Gm = self.alpha*np.identity(self.n, dtype=self.dtype) 

785 for c, d in zip(self.cs, self.ds): 

786 Gm += c[:,None]*d[None,:].conj() 

787 return Gm 

788 

789 def collapse(self): 

790 """Collapse the low-rank matrix to a full-rank one.""" 

791 self.collapsed = np.array(self) 

792 self.cs = None 

793 self.ds = None 

794 self.alpha = None 

795 

796 def restart_reduce(self, rank): 

797 """ 

798 Reduce the rank of the matrix by dropping all vectors. 

799 """ 

800 if self.collapsed is not None: 

801 return 

802 assert rank > 0 

803 if len(self.cs) > rank: 

804 del self.cs[:] 

805 del self.ds[:] 

806 

807 def simple_reduce(self, rank): 

808 """ 

809 Reduce the rank of the matrix by dropping oldest vectors. 

810 """ 

811 if self.collapsed is not None: 

812 return 

813 assert rank > 0 

814 while len(self.cs) > rank: 

815 del self.cs[0] 

816 del self.ds[0] 

817 

818 def svd_reduce(self, max_rank, to_retain=None): 

819 """ 

820 Reduce the rank of the matrix by retaining some SVD components. 

821 

822 This corresponds to the \"Broyden Rank Reduction Inverse\" 

823 algorithm described in [1]_. 

824 

825 Note that the SVD decomposition can be done by solving only a 

826 problem whose size is the effective rank of this matrix, which 

827 is viable even for large problems. 

828 

829 Parameters 

830 ---------- 

831 max_rank : int 

832 Maximum rank of this matrix after reduction. 

833 to_retain : int, optional 

834 Number of SVD components to retain when reduction is done 

835 (ie. rank > max_rank). Default is ``max_rank - 2``. 

836 

837 References 

838 ---------- 

839 .. [1] B.A. van der Rotten, PhD thesis, 

840 \"A limited memory Broyden method to solve high-dimensional 

841 systems of nonlinear equations\". Mathematisch Instituut, 

842 Universiteit Leiden, The Netherlands (2003). 

843 

844 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf 

845 

846 """ 

847 if self.collapsed is not None: 

848 return 

849 

850 p = max_rank 

851 if to_retain is not None: 

852 q = to_retain 

853 else: 

854 q = p - 2 

855 

856 if self.cs: 

857 p = min(p, len(self.cs[0])) 

858 q = max(0, min(q, p-1)) 

859 

860 m = len(self.cs) 

861 if m < p: 

862 # nothing to do 

863 return 

864 

865 C = np.array(self.cs).T 

866 D = np.array(self.ds).T 

867 

868 D, R = qr(D, mode='economic') 

869 C = dot(C, R.T.conj()) 

870 

871 U, S, WH = svd(C, full_matrices=False, compute_uv=True) 

872 

873 C = dot(C, inv(WH)) 

874 D = dot(D, WH.T.conj()) 

875 

876 for k in range(q): 

877 self.cs[k] = C[:,k].copy() 

878 self.ds[k] = D[:,k].copy() 

879 

880 del self.cs[q:] 

881 del self.ds[q:] 

882 

883 

884_doc_parts['broyden_params'] = """ 

885 alpha : float, optional 

886 Initial guess for the Jacobian is ``(-1/alpha)``. 

887 reduction_method : str or tuple, optional 

888 Method used in ensuring that the rank of the Broyden matrix 

889 stays low. Can either be a string giving the name of the method, 

890 or a tuple of the form ``(method, param1, param2, ...)`` 

891 that gives the name of the method and values for additional parameters. 

892 

893 Methods available: 

894 

895 - ``restart``: drop all matrix columns. Has no extra parameters. 

896 - ``simple``: drop oldest matrix column. Has no extra parameters. 

897 - ``svd``: keep only the most significant SVD components. 

898 Takes an extra parameter, ``to_retain``, which determines the 

899 number of SVD components to retain when rank reduction is done. 

900 Default is ``max_rank - 2``. 

901 

902 max_rank : int, optional 

903 Maximum rank for the Broyden matrix. 

904 Default is infinity (i.e., no rank reduction). 

905 """.strip() 

906 

907 

908class BroydenFirst(GenericBroyden): 

909 r""" 

910 Find a root of a function, using Broyden's first Jacobian approximation. 

911 

912 This method is also known as \"Broyden's good method\". 

913 

914 Parameters 

915 ---------- 

916 %(params_basic)s 

917 %(broyden_params)s 

918 %(params_extra)s 

919 

920 See Also 

921 -------- 

922 root : Interface to root finding algorithms for multivariate 

923 functions. See ``method=='broyden1'`` in particular. 

924 

925 Notes 

926 ----- 

927 This algorithm implements the inverse Jacobian Quasi-Newton update 

928 

929 .. math:: H_+ = H + (dx - H df) dx^\dagger H / ( dx^\dagger H df) 

930 

931 which corresponds to Broyden's first Jacobian update 

932 

933 .. math:: J_+ = J + (df - J dx) dx^\dagger / dx^\dagger dx 

934 

935 

936 References 

937 ---------- 

938 .. [1] B.A. van der Rotten, PhD thesis, 

939 \"A limited memory Broyden method to solve high-dimensional 

940 systems of nonlinear equations\". Mathematisch Instituut, 

941 Universiteit Leiden, The Netherlands (2003). 

942 

943 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf 

944 

945 """ 

946 

947 def __init__(self, alpha=None, reduction_method='restart', max_rank=None): 

948 GenericBroyden.__init__(self) 

949 self.alpha = alpha 

950 self.Gm = None 

951 

952 if max_rank is None: 

953 max_rank = np.inf 

954 self.max_rank = max_rank 

955 

956 if isinstance(reduction_method, str): 

957 reduce_params = () 

958 else: 

959 reduce_params = reduction_method[1:] 

960 reduction_method = reduction_method[0] 

961 reduce_params = (max_rank - 1,) + reduce_params 

962 

963 if reduction_method == 'svd': 

964 self._reduce = lambda: self.Gm.svd_reduce(*reduce_params) 

965 elif reduction_method == 'simple': 

966 self._reduce = lambda: self.Gm.simple_reduce(*reduce_params) 

967 elif reduction_method == 'restart': 

968 self._reduce = lambda: self.Gm.restart_reduce(*reduce_params) 

969 else: 

970 raise ValueError("Unknown rank reduction method '%s'" % 

971 reduction_method) 

972 

973 def setup(self, x, F, func): 

974 GenericBroyden.setup(self, x, F, func) 

975 self.Gm = LowRankMatrix(-self.alpha, self.shape[0], self.dtype) 

976 

977 def todense(self): 

978 return inv(self.Gm) 

979 

980 def solve(self, f, tol=0): 

981 r = self.Gm.matvec(f) 

982 if not np.isfinite(r).all(): 

983 # singular; reset the Jacobian approximation 

984 self.setup(self.last_x, self.last_f, self.func) 

985 return self.Gm.matvec(f) 

986 

987 def matvec(self, f): 

988 return self.Gm.solve(f) 

989 

990 def rsolve(self, f, tol=0): 

991 return self.Gm.rmatvec(f) 

992 

993 def rmatvec(self, f): 

994 return self.Gm.rsolve(f) 

995 

996 def _update(self, x, f, dx, df, dx_norm, df_norm): 

997 self._reduce() # reduce first to preserve secant condition 

998 

999 v = self.Gm.rmatvec(dx) 

1000 c = dx - self.Gm.matvec(df) 

1001 d = v / vdot(df, v) 

1002 

1003 self.Gm.append(c, d) 

1004 

1005 

1006class BroydenSecond(BroydenFirst): 

1007 """ 

1008 Find a root of a function, using Broyden\'s second Jacobian approximation. 

1009 

1010 This method is also known as \"Broyden's bad method\". 

1011 

1012 Parameters 

1013 ---------- 

1014 %(params_basic)s 

1015 %(broyden_params)s 

1016 %(params_extra)s 

1017 

1018 See Also 

1019 -------- 

1020 root : Interface to root finding algorithms for multivariate 

1021 functions. See ``method=='broyden2'`` in particular. 

1022 

1023 Notes 

1024 ----- 

1025 This algorithm implements the inverse Jacobian Quasi-Newton update 

1026 

1027 .. math:: H_+ = H + (dx - H df) df^\\dagger / ( df^\\dagger df) 

1028 

1029 corresponding to Broyden's second method. 

1030 

1031 References 

1032 ---------- 

1033 .. [1] B.A. van der Rotten, PhD thesis, 

1034 \"A limited memory Broyden method to solve high-dimensional 

1035 systems of nonlinear equations\". Mathematisch Instituut, 

1036 Universiteit Leiden, The Netherlands (2003). 

1037 

1038 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf 

1039 

1040 """ 

1041 

1042 def _update(self, x, f, dx, df, dx_norm, df_norm): 

1043 self._reduce() # reduce first to preserve secant condition 

1044 

1045 v = df 

1046 c = dx - self.Gm.matvec(df) 

1047 d = v / df_norm**2 

1048 self.Gm.append(c, d) 

1049 

1050 

1051#------------------------------------------------------------------------------ 

1052# Broyden-like (restricted memory) 

1053#------------------------------------------------------------------------------ 

1054 

1055class Anderson(GenericBroyden): 

1056 """ 

1057 Find a root of a function, using (extended) Anderson mixing. 

1058 

1059 The Jacobian is formed by for a 'best' solution in the space 

1060 spanned by last `M` vectors. As a result, only a MxM matrix 

1061 inversions and MxN multiplications are required. [Ey]_ 

1062 

1063 Parameters 

1064 ---------- 

1065 %(params_basic)s 

1066 alpha : float, optional 

1067 Initial guess for the Jacobian is (-1/alpha). 

1068 M : float, optional 

1069 Number of previous vectors to retain. Defaults to 5. 

1070 w0 : float, optional 

1071 Regularization parameter for numerical stability. 

1072 Compared to unity, good values of the order of 0.01. 

1073 %(params_extra)s 

1074 

1075 See Also 

1076 -------- 

1077 root : Interface to root finding algorithms for multivariate 

1078 functions. See ``method=='anderson'`` in particular. 

1079 

1080 References 

1081 ---------- 

1082 .. [Ey] V. Eyert, J. Comp. Phys., 124, 271 (1996). 

1083 

1084 """ 

1085 

1086 # Note: 

1087 # 

1088 # Anderson method maintains a rank M approximation of the inverse Jacobian, 

1089 # 

1090 # J^-1 v ~ -v*alpha + (dX + alpha dF) A^-1 dF^H v 

1091 # A = W + dF^H dF 

1092 # W = w0^2 diag(dF^H dF) 

1093 # 

1094 # so that for w0 = 0 the secant condition applies for last M iterates, i.e., 

1095 # 

1096 # J^-1 df_j = dx_j 

1097 # 

1098 # for all j = 0 ... M-1. 

1099 # 

1100 # Moreover, (from Sherman-Morrison-Woodbury formula) 

1101 # 

1102 # J v ~ [ b I - b^2 C (I + b dF^H A^-1 C)^-1 dF^H ] v 

1103 # C = (dX + alpha dF) A^-1 

1104 # b = -1/alpha 

1105 # 

1106 # and after simplification 

1107 # 

1108 # J v ~ -v/alpha + (dX/alpha + dF) (dF^H dX - alpha W)^-1 dF^H v 

1109 # 

1110 

1111 def __init__(self, alpha=None, w0=0.01, M=5): 

1112 GenericBroyden.__init__(self) 

1113 self.alpha = alpha 

1114 self.M = M 

1115 self.dx = [] 

1116 self.df = [] 

1117 self.gamma = None 

1118 self.w0 = w0 

1119 

1120 def solve(self, f, tol=0): 

1121 dx = -self.alpha*f 

1122 

1123 n = len(self.dx) 

1124 if n == 0: 

1125 return dx 

1126 

1127 df_f = np.empty(n, dtype=f.dtype) 

1128 for k in range(n): 

1129 df_f[k] = vdot(self.df[k], f) 

1130 

1131 try: 

1132 gamma = solve(self.a, df_f) 

1133 except LinAlgError: 

1134 # singular; reset the Jacobian approximation 

1135 del self.dx[:] 

1136 del self.df[:] 

1137 return dx 

1138 

1139 for m in range(n): 

1140 dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m]) 

1141 return dx 

1142 

1143 def matvec(self, f): 

1144 dx = -f/self.alpha 

1145 

1146 n = len(self.dx) 

1147 if n == 0: 

1148 return dx 

1149 

1150 df_f = np.empty(n, dtype=f.dtype) 

1151 for k in range(n): 

1152 df_f[k] = vdot(self.df[k], f) 

1153 

1154 b = np.empty((n, n), dtype=f.dtype) 

1155 for i in range(n): 

1156 for j in range(n): 

1157 b[i,j] = vdot(self.df[i], self.dx[j]) 

1158 if i == j and self.w0 != 0: 

1159 b[i,j] -= vdot(self.df[i], self.df[i])*self.w0**2*self.alpha 

1160 gamma = solve(b, df_f) 

1161 

1162 for m in range(n): 

1163 dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha) 

1164 return dx 

1165 

1166 def _update(self, x, f, dx, df, dx_norm, df_norm): 

1167 if self.M == 0: 

1168 return 

1169 

1170 self.dx.append(dx) 

1171 self.df.append(df) 

1172 

1173 while len(self.dx) > self.M: 

1174 self.dx.pop(0) 

1175 self.df.pop(0) 

1176 

1177 n = len(self.dx) 

1178 a = np.zeros((n, n), dtype=f.dtype) 

1179 

1180 for i in range(n): 

1181 for j in range(i, n): 

1182 if i == j: 

1183 wd = self.w0**2 

1184 else: 

1185 wd = 0 

1186 a[i,j] = (1+wd)*vdot(self.df[i], self.df[j]) 

1187 

1188 a += np.triu(a, 1).T.conj() 

1189 self.a = a 

1190 

1191#------------------------------------------------------------------------------ 

1192# Simple iterations 

1193#------------------------------------------------------------------------------ 

1194 

1195 

1196class DiagBroyden(GenericBroyden): 

1197 """ 

1198 Find a root of a function, using diagonal Broyden Jacobian approximation. 

1199 

1200 The Jacobian approximation is derived from previous iterations, by 

1201 retaining only the diagonal of Broyden matrices. 

1202 

1203 .. warning:: 

1204 

1205 This algorithm may be useful for specific problems, but whether 

1206 it will work may depend strongly on the problem. 

1207 

1208 Parameters 

1209 ---------- 

1210 %(params_basic)s 

1211 alpha : float, optional 

1212 Initial guess for the Jacobian is (-1/alpha). 

1213 %(params_extra)s 

1214 

1215 See Also 

1216 -------- 

1217 root : Interface to root finding algorithms for multivariate 

1218 functions. See ``method=='diagbroyden'`` in particular. 

1219 """ 

1220 

1221 def __init__(self, alpha=None): 

1222 GenericBroyden.__init__(self) 

1223 self.alpha = alpha 

1224 

1225 def setup(self, x, F, func): 

1226 GenericBroyden.setup(self, x, F, func) 

1227 self.d = np.full((self.shape[0],), 1 / self.alpha, dtype=self.dtype) 

1228 

1229 def solve(self, f, tol=0): 

1230 return -f / self.d 

1231 

1232 def matvec(self, f): 

1233 return -f * self.d 

1234 

1235 def rsolve(self, f, tol=0): 

1236 return -f / self.d.conj() 

1237 

1238 def rmatvec(self, f): 

1239 return -f * self.d.conj() 

1240 

1241 def todense(self): 

1242 return np.diag(-self.d) 

1243 

1244 def _update(self, x, f, dx, df, dx_norm, df_norm): 

1245 self.d -= (df + self.d*dx)*dx/dx_norm**2 

1246 

1247 

1248class LinearMixing(GenericBroyden): 

1249 """ 

1250 Find a root of a function, using a scalar Jacobian approximation. 

1251 

1252 .. warning:: 

1253 

1254 This algorithm may be useful for specific problems, but whether 

1255 it will work may depend strongly on the problem. 

1256 

1257 Parameters 

1258 ---------- 

1259 %(params_basic)s 

1260 alpha : float, optional 

1261 The Jacobian approximation is (-1/alpha). 

1262 %(params_extra)s 

1263 

1264 See Also 

1265 -------- 

1266 root : Interface to root finding algorithms for multivariate 

1267 functions. See ``method=='linearmixing'`` in particular. 

1268 

1269 """ 

1270 

1271 def __init__(self, alpha=None): 

1272 GenericBroyden.__init__(self) 

1273 self.alpha = alpha 

1274 

1275 def solve(self, f, tol=0): 

1276 return -f*self.alpha 

1277 

1278 def matvec(self, f): 

1279 return -f/self.alpha 

1280 

1281 def rsolve(self, f, tol=0): 

1282 return -f*np.conj(self.alpha) 

1283 

1284 def rmatvec(self, f): 

1285 return -f/np.conj(self.alpha) 

1286 

1287 def todense(self): 

1288 return np.diag(np.full(self.shape[0], -1/self.alpha)) 

1289 

1290 def _update(self, x, f, dx, df, dx_norm, df_norm): 

1291 pass 

1292 

1293 

1294class ExcitingMixing(GenericBroyden): 

1295 """ 

1296 Find a root of a function, using a tuned diagonal Jacobian approximation. 

1297 

1298 The Jacobian matrix is diagonal and is tuned on each iteration. 

1299 

1300 .. warning:: 

1301 

1302 This algorithm may be useful for specific problems, but whether 

1303 it will work may depend strongly on the problem. 

1304 

1305 See Also 

1306 -------- 

1307 root : Interface to root finding algorithms for multivariate 

1308 functions. See ``method=='excitingmixing'`` in particular. 

1309 

1310 Parameters 

1311 ---------- 

1312 %(params_basic)s 

1313 alpha : float, optional 

1314 Initial Jacobian approximation is (-1/alpha). 

1315 alphamax : float, optional 

1316 The entries of the diagonal Jacobian are kept in the range 

1317 ``[alpha, alphamax]``. 

1318 %(params_extra)s 

1319 """ 

1320 

1321 def __init__(self, alpha=None, alphamax=1.0): 

1322 GenericBroyden.__init__(self) 

1323 self.alpha = alpha 

1324 self.alphamax = alphamax 

1325 self.beta = None 

1326 

1327 def setup(self, x, F, func): 

1328 GenericBroyden.setup(self, x, F, func) 

1329 self.beta = np.full((self.shape[0],), self.alpha, dtype=self.dtype) 

1330 

1331 def solve(self, f, tol=0): 

1332 return -f*self.beta 

1333 

1334 def matvec(self, f): 

1335 return -f/self.beta 

1336 

1337 def rsolve(self, f, tol=0): 

1338 return -f*self.beta.conj() 

1339 

1340 def rmatvec(self, f): 

1341 return -f/self.beta.conj() 

1342 

1343 def todense(self): 

1344 return np.diag(-1/self.beta) 

1345 

1346 def _update(self, x, f, dx, df, dx_norm, df_norm): 

1347 incr = f*self.last_f > 0 

1348 self.beta[incr] += self.alpha 

1349 self.beta[~incr] = self.alpha 

1350 np.clip(self.beta, 0, self.alphamax, out=self.beta) 

1351 

1352 

1353#------------------------------------------------------------------------------ 

1354# Iterative/Krylov approximated Jacobians 

1355#------------------------------------------------------------------------------ 

1356 

1357class KrylovJacobian(Jacobian): 

1358 r""" 

1359 Find a root of a function, using Krylov approximation for inverse Jacobian. 

1360 

1361 This method is suitable for solving large-scale problems. 

1362 

1363 Parameters 

1364 ---------- 

1365 %(params_basic)s 

1366 rdiff : float, optional 

1367 Relative step size to use in numerical differentiation. 

1368 method : {'lgmres', 'gmres', 'bicgstab', 'cgs', 'minres'} or function 

1369 Krylov method to use to approximate the Jacobian. 

1370 Can be a string, or a function implementing the same interface as 

1371 the iterative solvers in `scipy.sparse.linalg`. 

1372 

1373 The default is `scipy.sparse.linalg.lgmres`. 

1374 inner_maxiter : int, optional 

1375 Parameter to pass to the "inner" Krylov solver: maximum number of 

1376 iterations. Iteration will stop after maxiter steps even if the 

1377 specified tolerance has not been achieved. 

1378 inner_M : LinearOperator or InverseJacobian 

1379 Preconditioner for the inner Krylov iteration. 

1380 Note that you can use also inverse Jacobians as (adaptive) 

1381 preconditioners. For example, 

1382 

1383 >>> from scipy.optimize.nonlin import BroydenFirst, KrylovJacobian 

1384 >>> from scipy.optimize.nonlin import InverseJacobian 

1385 >>> jac = BroydenFirst() 

1386 >>> kjac = KrylovJacobian(inner_M=InverseJacobian(jac)) 

1387 

1388 If the preconditioner has a method named 'update', it will be called 

1389 as ``update(x, f)`` after each nonlinear step, with ``x`` giving 

1390 the current point, and ``f`` the current function value. 

1391 outer_k : int, optional 

1392 Size of the subspace kept across LGMRES nonlinear iterations. 

1393 See `scipy.sparse.linalg.lgmres` for details. 

1394 inner_kwargs : kwargs 

1395 Keyword parameters for the "inner" Krylov solver 

1396 (defined with `method`). Parameter names must start with 

1397 the `inner_` prefix which will be stripped before passing on 

1398 the inner method. See, e.g., `scipy.sparse.linalg.gmres` for details. 

1399 %(params_extra)s 

1400 

1401 See Also 

1402 -------- 

1403 root : Interface to root finding algorithms for multivariate 

1404 functions. See ``method=='krylov'`` in particular. 

1405 scipy.sparse.linalg.gmres 

1406 scipy.sparse.linalg.lgmres 

1407 

1408 Notes 

1409 ----- 

1410 This function implements a Newton-Krylov solver. The basic idea is 

1411 to compute the inverse of the Jacobian with an iterative Krylov 

1412 method. These methods require only evaluating the Jacobian-vector 

1413 products, which are conveniently approximated by a finite difference: 

1414 

1415 .. math:: J v \approx (f(x + \omega*v/|v|) - f(x)) / \omega 

1416 

1417 Due to the use of iterative matrix inverses, these methods can 

1418 deal with large nonlinear problems. 

1419 

1420 SciPy's `scipy.sparse.linalg` module offers a selection of Krylov 

1421 solvers to choose from. The default here is `lgmres`, which is a 

1422 variant of restarted GMRES iteration that reuses some of the 

1423 information obtained in the previous Newton steps to invert 

1424 Jacobians in subsequent steps. 

1425 

1426 For a review on Newton-Krylov methods, see for example [1]_, 

1427 and for the LGMRES sparse inverse method, see [2]_. 

1428 

1429 References 

1430 ---------- 

1431 .. [1] D.A. Knoll and D.E. Keyes, J. Comp. Phys. 193, 357 (2004). 

1432 :doi:`10.1016/j.jcp.2003.08.010` 

1433 .. [2] A.H. Baker and E.R. Jessup and T. Manteuffel, 

1434 SIAM J. Matrix Anal. Appl. 26, 962 (2005). 

1435 :doi:`10.1137/S0895479803422014` 

1436 

1437 """ 

1438 

1439 def __init__(self, rdiff=None, method='lgmres', inner_maxiter=20, 

1440 inner_M=None, outer_k=10, **kw): 

1441 self.preconditioner = inner_M 

1442 self.rdiff = rdiff 

1443 self.method = dict( 

1444 bicgstab=scipy.sparse.linalg.bicgstab, 

1445 gmres=scipy.sparse.linalg.gmres, 

1446 lgmres=scipy.sparse.linalg.lgmres, 

1447 cgs=scipy.sparse.linalg.cgs, 

1448 minres=scipy.sparse.linalg.minres, 

1449 ).get(method, method) 

1450 

1451 self.method_kw = dict(maxiter=inner_maxiter, M=self.preconditioner) 

1452 

1453 if self.method is scipy.sparse.linalg.gmres: 

1454 # Replace GMRES's outer iteration with Newton steps 

1455 self.method_kw['restrt'] = inner_maxiter 

1456 self.method_kw['maxiter'] = 1 

1457 self.method_kw.setdefault('atol', 0) 

1458 elif self.method is scipy.sparse.linalg.gcrotmk: 

1459 self.method_kw.setdefault('atol', 0) 

1460 elif self.method is scipy.sparse.linalg.lgmres: 

1461 self.method_kw['outer_k'] = outer_k 

1462 # Replace LGMRES's outer iteration with Newton steps 

1463 self.method_kw['maxiter'] = 1 

1464 # Carry LGMRES's `outer_v` vectors across nonlinear iterations 

1465 self.method_kw.setdefault('outer_v', []) 

1466 self.method_kw.setdefault('prepend_outer_v', True) 

1467 # But don't carry the corresponding Jacobian*v products, in case 

1468 # the Jacobian changes a lot in the nonlinear step 

1469 # 

1470 # XXX: some trust-region inspired ideas might be more efficient... 

1471 # See e.g., Brown & Saad. But needs to be implemented separately 

1472 # since it's not an inexact Newton method. 

1473 self.method_kw.setdefault('store_outer_Av', False) 

1474 self.method_kw.setdefault('atol', 0) 

1475 

1476 for key, value in kw.items(): 

1477 if not key.startswith('inner_'): 

1478 raise ValueError("Unknown parameter %s" % key) 

1479 self.method_kw[key[6:]] = value 

1480 

1481 def _update_diff_step(self): 

1482 mx = abs(self.x0).max() 

1483 mf = abs(self.f0).max() 

1484 self.omega = self.rdiff * max(1, mx) / max(1, mf) 

1485 

1486 def matvec(self, v): 

1487 nv = norm(v) 

1488 if nv == 0: 

1489 return 0*v 

1490 sc = self.omega / nv 

1491 r = (self.func(self.x0 + sc*v) - self.f0) / sc 

1492 if not np.all(np.isfinite(r)) and np.all(np.isfinite(v)): 

1493 raise ValueError('Function returned non-finite results') 

1494 return r 

1495 

1496 def solve(self, rhs, tol=0): 

1497 if 'tol' in self.method_kw: 

1498 sol, info = self.method(self.op, rhs, **self.method_kw) 

1499 else: 

1500 sol, info = self.method(self.op, rhs, tol=tol, **self.method_kw) 

1501 return sol 

1502 

1503 def update(self, x, f): 

1504 self.x0 = x 

1505 self.f0 = f 

1506 self._update_diff_step() 

1507 

1508 # Update also the preconditioner, if possible 

1509 if self.preconditioner is not None: 

1510 if hasattr(self.preconditioner, 'update'): 

1511 self.preconditioner.update(x, f) 

1512 

1513 def setup(self, x, f, func): 

1514 Jacobian.setup(self, x, f, func) 

1515 self.x0 = x 

1516 self.f0 = f 

1517 self.op = scipy.sparse.linalg.aslinearoperator(self) 

1518 

1519 if self.rdiff is None: 

1520 self.rdiff = np.finfo(x.dtype).eps ** (1./2) 

1521 

1522 self._update_diff_step() 

1523 

1524 # Setup also the preconditioner, if possible 

1525 if self.preconditioner is not None: 

1526 if hasattr(self.preconditioner, 'setup'): 

1527 self.preconditioner.setup(x, f, func) 

1528 

1529 

1530#------------------------------------------------------------------------------ 

1531# Wrapper functions 

1532#------------------------------------------------------------------------------ 

1533 

1534def _nonlin_wrapper(name, jac): 

1535 """ 

1536 Construct a solver wrapper with given name and Jacobian approx. 

1537 

1538 It inspects the keyword arguments of ``jac.__init__``, and allows to 

1539 use the same arguments in the wrapper function, in addition to the 

1540 keyword arguments of `nonlin_solve` 

1541 

1542 """ 

1543 signature = _getfullargspec(jac.__init__) 

1544 args, varargs, varkw, defaults, kwonlyargs, kwdefaults, _ = signature 

1545 kwargs = list(zip(args[-len(defaults):], defaults)) 

1546 kw_str = ", ".join(["%s=%r" % (k, v) for k, v in kwargs]) 

1547 if kw_str: 

1548 kw_str = ", " + kw_str 

1549 kwkw_str = ", ".join(["%s=%s" % (k, k) for k, v in kwargs]) 

1550 if kwkw_str: 

1551 kwkw_str = kwkw_str + ", " 

1552 if kwonlyargs: 

1553 raise ValueError('Unexpected signature %s' % signature) 

1554 

1555 # Construct the wrapper function so that its keyword arguments 

1556 # are visible in pydoc.help etc. 

1557 wrapper = """ 

1558def %(name)s(F, xin, iter=None %(kw)s, verbose=False, maxiter=None, 

1559 f_tol=None, f_rtol=None, x_tol=None, x_rtol=None, 

1560 tol_norm=None, line_search='armijo', callback=None, **kw): 

1561 jac = %(jac)s(%(kwkw)s **kw) 

1562 return nonlin_solve(F, xin, jac, iter, verbose, maxiter, 

1563 f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, 

1564 callback) 

1565""" 

1566 

1567 wrapper = wrapper % dict(name=name, kw=kw_str, jac=jac.__name__, 

1568 kwkw=kwkw_str) 

1569 ns = {} 

1570 ns.update(globals()) 

1571 exec(wrapper, ns) 

1572 func = ns[name] 

1573 func.__doc__ = jac.__doc__ 

1574 _set_doc(func) 

1575 return func 

1576 

1577 

1578broyden1 = _nonlin_wrapper('broyden1', BroydenFirst) 

1579broyden2 = _nonlin_wrapper('broyden2', BroydenSecond) 

1580anderson = _nonlin_wrapper('anderson', Anderson) 

1581linearmixing = _nonlin_wrapper('linearmixing', LinearMixing) 

1582diagbroyden = _nonlin_wrapper('diagbroyden', DiagBroyden) 

1583excitingmixing = _nonlin_wrapper('excitingmixing', ExcitingMixing) 

1584newton_krylov = _nonlin_wrapper('newton_krylov', KrylovJacobian)