Coverage for /usr/lib/python3/dist-packages/scipy/optimize/_nonlin.py: 22%
631 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1# Copyright (C) 2009, Pauli Virtanen <pav@iki.fi>
2# Distributed under the same license as SciPy.
4import sys
5import numpy as np
6from scipy.linalg import norm, solve, inv, qr, svd, LinAlgError
7from numpy import asarray, dot, vdot
8import scipy.sparse.linalg
9import scipy.sparse
10from scipy.linalg import get_blas_funcs
11import inspect
12from scipy._lib._util import getfullargspec_no_self as _getfullargspec
13from ._linesearch import scalar_search_wolfe1, scalar_search_armijo
16__all__ = [
17 'broyden1', 'broyden2', 'anderson', 'linearmixing',
18 'diagbroyden', 'excitingmixing', 'newton_krylov',
19 'BroydenFirst', 'KrylovJacobian', 'InverseJacobian']
21#------------------------------------------------------------------------------
22# Utility functions
23#------------------------------------------------------------------------------
26class NoConvergence(Exception):
27 pass
30def maxnorm(x):
31 return np.absolute(x).max()
34def _as_inexact(x):
35 """Return `x` as an array, of either floats or complex floats"""
36 x = asarray(x)
37 if not np.issubdtype(x.dtype, np.inexact):
38 return asarray(x, dtype=np.float_)
39 return x
42def _array_like(x, x0):
43 """Return ndarray `x` as same array subclass and shape as `x0`"""
44 x = np.reshape(x, np.shape(x0))
45 wrap = getattr(x0, '__array_wrap__', x.__array_wrap__)
46 return wrap(x)
49def _safe_norm(v):
50 if not np.isfinite(v).all():
51 return np.array(np.inf)
52 return norm(v)
54#------------------------------------------------------------------------------
55# Generic nonlinear solver machinery
56#------------------------------------------------------------------------------
59_doc_parts = dict(
60 params_basic="""
61 F : function(x) -> f
62 Function whose root to find; should take and return an array-like
63 object.
64 xin : array_like
65 Initial guess for the solution
66 """.strip(),
67 params_extra="""
68 iter : int, optional
69 Number of iterations to make. If omitted (default), make as many
70 as required to meet tolerances.
71 verbose : bool, optional
72 Print status to stdout on every iteration.
73 maxiter : int, optional
74 Maximum number of iterations to make. If more are needed to
75 meet convergence, `NoConvergence` is raised.
76 f_tol : float, optional
77 Absolute tolerance (in max-norm) for the residual.
78 If omitted, default is 6e-6.
79 f_rtol : float, optional
80 Relative tolerance for the residual. If omitted, not used.
81 x_tol : float, optional
82 Absolute minimum step size, as determined from the Jacobian
83 approximation. If the step size is smaller than this, optimization
84 is terminated as successful. If omitted, not used.
85 x_rtol : float, optional
86 Relative minimum step size. If omitted, not used.
87 tol_norm : function(vector) -> scalar, optional
88 Norm to use in convergence check. Default is the maximum norm.
89 line_search : {None, 'armijo' (default), 'wolfe'}, optional
90 Which type of a line search to use to determine the step size in the
91 direction given by the Jacobian approximation. Defaults to 'armijo'.
92 callback : function, optional
93 Optional callback function. It is called on every iteration as
94 ``callback(x, f)`` where `x` is the current solution and `f`
95 the corresponding residual.
97 Returns
98 -------
99 sol : ndarray
100 An array (of similar array type as `x0`) containing the final solution.
102 Raises
103 ------
104 NoConvergence
105 When a solution was not found.
107 """.strip()
108)
111def _set_doc(obj):
112 if obj.__doc__:
113 obj.__doc__ = obj.__doc__ % _doc_parts
116def nonlin_solve(F, x0, jacobian='krylov', iter=None, verbose=False,
117 maxiter=None, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None,
118 tol_norm=None, line_search='armijo', callback=None,
119 full_output=False, raise_exception=True):
120 """
121 Find a root of a function, in a way suitable for large-scale problems.
123 Parameters
124 ----------
125 %(params_basic)s
126 jacobian : Jacobian
127 A Jacobian approximation: `Jacobian` object or something that
128 `asjacobian` can transform to one. Alternatively, a string specifying
129 which of the builtin Jacobian approximations to use:
131 krylov, broyden1, broyden2, anderson
132 diagbroyden, linearmixing, excitingmixing
134 %(params_extra)s
135 full_output : bool
136 If true, returns a dictionary `info` containing convergence
137 information.
138 raise_exception : bool
139 If True, a `NoConvergence` exception is raise if no solution is found.
141 See Also
142 --------
143 asjacobian, Jacobian
145 Notes
146 -----
147 This algorithm implements the inexact Newton method, with
148 backtracking or full line searches. Several Jacobian
149 approximations are available, including Krylov and Quasi-Newton
150 methods.
152 References
153 ----------
154 .. [KIM] C. T. Kelley, \"Iterative Methods for Linear and Nonlinear
155 Equations\". Society for Industrial and Applied Mathematics. (1995)
156 https://archive.siam.org/books/kelley/fr16/
158 """
159 # Can't use default parameters because it's being explicitly passed as None
160 # from the calling function, so we need to set it here.
161 tol_norm = maxnorm if tol_norm is None else tol_norm
162 condition = TerminationCondition(f_tol=f_tol, f_rtol=f_rtol,
163 x_tol=x_tol, x_rtol=x_rtol,
164 iter=iter, norm=tol_norm)
166 x0 = _as_inexact(x0)
167 def func(z):
168 return _as_inexact(F(_array_like(z, x0))).flatten()
169 x = x0.flatten()
171 dx = np.full_like(x, np.inf)
172 Fx = func(x)
173 Fx_norm = norm(Fx)
175 jacobian = asjacobian(jacobian)
176 jacobian.setup(x.copy(), Fx, func)
178 if maxiter is None:
179 if iter is not None:
180 maxiter = iter + 1
181 else:
182 maxiter = 100*(x.size+1)
184 if line_search is True:
185 line_search = 'armijo'
186 elif line_search is False:
187 line_search = None
189 if line_search not in (None, 'armijo', 'wolfe'):
190 raise ValueError("Invalid line search")
192 # Solver tolerance selection
193 gamma = 0.9
194 eta_max = 0.9999
195 eta_treshold = 0.1
196 eta = 1e-3
198 for n in range(maxiter):
199 status = condition.check(Fx, x, dx)
200 if status:
201 break
203 # The tolerance, as computed for scipy.sparse.linalg.* routines
204 tol = min(eta, eta*Fx_norm)
205 dx = -jacobian.solve(Fx, tol=tol)
207 if norm(dx) == 0:
208 raise ValueError("Jacobian inversion yielded zero vector. "
209 "This indicates a bug in the Jacobian "
210 "approximation.")
212 # Line search, or Newton step
213 if line_search:
214 s, x, Fx, Fx_norm_new = _nonlin_line_search(func, x, Fx, dx,
215 line_search)
216 else:
217 s = 1.0
218 x = x + dx
219 Fx = func(x)
220 Fx_norm_new = norm(Fx)
222 jacobian.update(x.copy(), Fx)
224 if callback:
225 callback(x, Fx)
227 # Adjust forcing parameters for inexact methods
228 eta_A = gamma * Fx_norm_new**2 / Fx_norm**2
229 if gamma * eta**2 < eta_treshold:
230 eta = min(eta_max, eta_A)
231 else:
232 eta = min(eta_max, max(eta_A, gamma*eta**2))
234 Fx_norm = Fx_norm_new
236 # Print status
237 if verbose:
238 sys.stdout.write("%d: |F(x)| = %g; step %g\n" % (
239 n, tol_norm(Fx), s))
240 sys.stdout.flush()
241 else:
242 if raise_exception:
243 raise NoConvergence(_array_like(x, x0))
244 else:
245 status = 2
247 if full_output:
248 info = {'nit': condition.iteration,
249 'fun': Fx,
250 'status': status,
251 'success': status == 1,
252 'message': {1: 'A solution was found at the specified '
253 'tolerance.',
254 2: 'The maximum number of iterations allowed '
255 'has been reached.'
256 }[status]
257 }
258 return _array_like(x, x0), info
259 else:
260 return _array_like(x, x0)
263_set_doc(nonlin_solve)
266def _nonlin_line_search(func, x, Fx, dx, search_type='armijo', rdiff=1e-8,
267 smin=1e-2):
268 tmp_s = [0]
269 tmp_Fx = [Fx]
270 tmp_phi = [norm(Fx)**2]
271 s_norm = norm(x) / norm(dx)
273 def phi(s, store=True):
274 if s == tmp_s[0]:
275 return tmp_phi[0]
276 xt = x + s*dx
277 v = func(xt)
278 p = _safe_norm(v)**2
279 if store:
280 tmp_s[0] = s
281 tmp_phi[0] = p
282 tmp_Fx[0] = v
283 return p
285 def derphi(s):
286 ds = (abs(s) + s_norm + 1) * rdiff
287 return (phi(s+ds, store=False) - phi(s)) / ds
289 if search_type == 'wolfe':
290 s, phi1, phi0 = scalar_search_wolfe1(phi, derphi, tmp_phi[0],
291 xtol=1e-2, amin=smin)
292 elif search_type == 'armijo':
293 s, phi1 = scalar_search_armijo(phi, tmp_phi[0], -tmp_phi[0],
294 amin=smin)
296 if s is None:
297 # XXX: No suitable step length found. Take the full Newton step,
298 # and hope for the best.
299 s = 1.0
301 x = x + s*dx
302 if s == tmp_s[0]:
303 Fx = tmp_Fx[0]
304 else:
305 Fx = func(x)
306 Fx_norm = norm(Fx)
308 return s, x, Fx, Fx_norm
311class TerminationCondition:
312 """
313 Termination condition for an iteration. It is terminated if
315 - |F| < f_rtol*|F_0|, AND
316 - |F| < f_tol
318 AND
320 - |dx| < x_rtol*|x|, AND
321 - |dx| < x_tol
323 """
324 def __init__(self, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None,
325 iter=None, norm=maxnorm):
327 if f_tol is None:
328 f_tol = np.finfo(np.float_).eps ** (1./3)
329 if f_rtol is None:
330 f_rtol = np.inf
331 if x_tol is None:
332 x_tol = np.inf
333 if x_rtol is None:
334 x_rtol = np.inf
336 self.x_tol = x_tol
337 self.x_rtol = x_rtol
338 self.f_tol = f_tol
339 self.f_rtol = f_rtol
341 self.norm = norm
343 self.iter = iter
345 self.f0_norm = None
346 self.iteration = 0
348 def check(self, f, x, dx):
349 self.iteration += 1
350 f_norm = self.norm(f)
351 x_norm = self.norm(x)
352 dx_norm = self.norm(dx)
354 if self.f0_norm is None:
355 self.f0_norm = f_norm
357 if f_norm == 0:
358 return 1
360 if self.iter is not None:
361 # backwards compatibility with SciPy 0.6.0
362 return 2 * (self.iteration > self.iter)
364 # NB: condition must succeed for rtol=inf even if norm == 0
365 return int((f_norm <= self.f_tol
366 and f_norm/self.f_rtol <= self.f0_norm)
367 and (dx_norm <= self.x_tol
368 and dx_norm/self.x_rtol <= x_norm))
371#------------------------------------------------------------------------------
372# Generic Jacobian approximation
373#------------------------------------------------------------------------------
375class Jacobian:
376 """
377 Common interface for Jacobians or Jacobian approximations.
379 The optional methods come useful when implementing trust region
380 etc., algorithms that often require evaluating transposes of the
381 Jacobian.
383 Methods
384 -------
385 solve
386 Returns J^-1 * v
387 update
388 Updates Jacobian to point `x` (where the function has residual `Fx`)
390 matvec : optional
391 Returns J * v
392 rmatvec : optional
393 Returns A^H * v
394 rsolve : optional
395 Returns A^-H * v
396 matmat : optional
397 Returns A * V, where V is a dense matrix with dimensions (N,K).
398 todense : optional
399 Form the dense Jacobian matrix. Necessary for dense trust region
400 algorithms, and useful for testing.
402 Attributes
403 ----------
404 shape
405 Matrix dimensions (M, N)
406 dtype
407 Data type of the matrix.
408 func : callable, optional
409 Function the Jacobian corresponds to
411 """
413 def __init__(self, **kw):
414 names = ["solve", "update", "matvec", "rmatvec", "rsolve",
415 "matmat", "todense", "shape", "dtype"]
416 for name, value in kw.items():
417 if name not in names:
418 raise ValueError("Unknown keyword argument %s" % name)
419 if value is not None:
420 setattr(self, name, kw[name])
422 if hasattr(self, 'todense'):
423 self.__array__ = lambda: self.todense()
425 def aspreconditioner(self):
426 return InverseJacobian(self)
428 def solve(self, v, tol=0):
429 raise NotImplementedError
431 def update(self, x, F):
432 pass
434 def setup(self, x, F, func):
435 self.func = func
436 self.shape = (F.size, x.size)
437 self.dtype = F.dtype
438 if self.__class__.setup is Jacobian.setup:
439 # Call on the first point unless overridden
440 self.update(x, F)
443class InverseJacobian:
444 def __init__(self, jacobian):
445 self.jacobian = jacobian
446 self.matvec = jacobian.solve
447 self.update = jacobian.update
448 if hasattr(jacobian, 'setup'):
449 self.setup = jacobian.setup
450 if hasattr(jacobian, 'rsolve'):
451 self.rmatvec = jacobian.rsolve
453 @property
454 def shape(self):
455 return self.jacobian.shape
457 @property
458 def dtype(self):
459 return self.jacobian.dtype
462def asjacobian(J):
463 """
464 Convert given object to one suitable for use as a Jacobian.
465 """
466 spsolve = scipy.sparse.linalg.spsolve
467 if isinstance(J, Jacobian):
468 return J
469 elif inspect.isclass(J) and issubclass(J, Jacobian):
470 return J()
471 elif isinstance(J, np.ndarray):
472 if J.ndim > 2:
473 raise ValueError('array must have rank <= 2')
474 J = np.atleast_2d(np.asarray(J))
475 if J.shape[0] != J.shape[1]:
476 raise ValueError('array must be square')
478 return Jacobian(matvec=lambda v: dot(J, v),
479 rmatvec=lambda v: dot(J.conj().T, v),
480 solve=lambda v: solve(J, v),
481 rsolve=lambda v: solve(J.conj().T, v),
482 dtype=J.dtype, shape=J.shape)
483 elif scipy.sparse.isspmatrix(J):
484 if J.shape[0] != J.shape[1]:
485 raise ValueError('matrix must be square')
486 return Jacobian(matvec=lambda v: J*v,
487 rmatvec=lambda v: J.conj().T * v,
488 solve=lambda v: spsolve(J, v),
489 rsolve=lambda v: spsolve(J.conj().T, v),
490 dtype=J.dtype, shape=J.shape)
491 elif hasattr(J, 'shape') and hasattr(J, 'dtype') and hasattr(J, 'solve'):
492 return Jacobian(matvec=getattr(J, 'matvec'),
493 rmatvec=getattr(J, 'rmatvec'),
494 solve=J.solve,
495 rsolve=getattr(J, 'rsolve'),
496 update=getattr(J, 'update'),
497 setup=getattr(J, 'setup'),
498 dtype=J.dtype,
499 shape=J.shape)
500 elif callable(J):
501 # Assume it's a function J(x) that returns the Jacobian
502 class Jac(Jacobian):
503 def update(self, x, F):
504 self.x = x
506 def solve(self, v, tol=0):
507 m = J(self.x)
508 if isinstance(m, np.ndarray):
509 return solve(m, v)
510 elif scipy.sparse.isspmatrix(m):
511 return spsolve(m, v)
512 else:
513 raise ValueError("Unknown matrix type")
515 def matvec(self, v):
516 m = J(self.x)
517 if isinstance(m, np.ndarray):
518 return dot(m, v)
519 elif scipy.sparse.isspmatrix(m):
520 return m*v
521 else:
522 raise ValueError("Unknown matrix type")
524 def rsolve(self, v, tol=0):
525 m = J(self.x)
526 if isinstance(m, np.ndarray):
527 return solve(m.conj().T, v)
528 elif scipy.sparse.isspmatrix(m):
529 return spsolve(m.conj().T, v)
530 else:
531 raise ValueError("Unknown matrix type")
533 def rmatvec(self, v):
534 m = J(self.x)
535 if isinstance(m, np.ndarray):
536 return dot(m.conj().T, v)
537 elif scipy.sparse.isspmatrix(m):
538 return m.conj().T * v
539 else:
540 raise ValueError("Unknown matrix type")
541 return Jac()
542 elif isinstance(J, str):
543 return dict(broyden1=BroydenFirst,
544 broyden2=BroydenSecond,
545 anderson=Anderson,
546 diagbroyden=DiagBroyden,
547 linearmixing=LinearMixing,
548 excitingmixing=ExcitingMixing,
549 krylov=KrylovJacobian)[J]()
550 else:
551 raise TypeError('Cannot convert object to a Jacobian')
554#------------------------------------------------------------------------------
555# Broyden
556#------------------------------------------------------------------------------
558class GenericBroyden(Jacobian):
559 def setup(self, x0, f0, func):
560 Jacobian.setup(self, x0, f0, func)
561 self.last_f = f0
562 self.last_x = x0
564 if hasattr(self, 'alpha') and self.alpha is None:
565 # Autoscale the initial Jacobian parameter
566 # unless we have already guessed the solution.
567 normf0 = norm(f0)
568 if normf0:
569 self.alpha = 0.5*max(norm(x0), 1) / normf0
570 else:
571 self.alpha = 1.0
573 def _update(self, x, f, dx, df, dx_norm, df_norm):
574 raise NotImplementedError
576 def update(self, x, f):
577 df = f - self.last_f
578 dx = x - self.last_x
579 self._update(x, f, dx, df, norm(dx), norm(df))
580 self.last_f = f
581 self.last_x = x
584class LowRankMatrix:
585 r"""
586 A matrix represented as
588 .. math:: \alpha I + \sum_{n=0}^{n=M} c_n d_n^\dagger
590 However, if the rank of the matrix reaches the dimension of the vectors,
591 full matrix representation will be used thereon.
593 """
595 def __init__(self, alpha, n, dtype):
596 self.alpha = alpha
597 self.cs = []
598 self.ds = []
599 self.n = n
600 self.dtype = dtype
601 self.collapsed = None
603 @staticmethod
604 def _matvec(v, alpha, cs, ds):
605 axpy, scal, dotc = get_blas_funcs(['axpy', 'scal', 'dotc'],
606 cs[:1] + [v])
607 w = alpha * v
608 for c, d in zip(cs, ds):
609 a = dotc(d, v)
610 w = axpy(c, w, w.size, a)
611 return w
613 @staticmethod
614 def _solve(v, alpha, cs, ds):
615 """Evaluate w = M^-1 v"""
616 if len(cs) == 0:
617 return v/alpha
619 # (B + C D^H)^-1 = B^-1 - B^-1 C (I + D^H B^-1 C)^-1 D^H B^-1
621 axpy, dotc = get_blas_funcs(['axpy', 'dotc'], cs[:1] + [v])
623 c0 = cs[0]
624 A = alpha * np.identity(len(cs), dtype=c0.dtype)
625 for i, d in enumerate(ds):
626 for j, c in enumerate(cs):
627 A[i,j] += dotc(d, c)
629 q = np.zeros(len(cs), dtype=c0.dtype)
630 for j, d in enumerate(ds):
631 q[j] = dotc(d, v)
632 q /= alpha
633 q = solve(A, q)
635 w = v/alpha
636 for c, qc in zip(cs, q):
637 w = axpy(c, w, w.size, -qc)
639 return w
641 def matvec(self, v):
642 """Evaluate w = M v"""
643 if self.collapsed is not None:
644 return np.dot(self.collapsed, v)
645 return LowRankMatrix._matvec(v, self.alpha, self.cs, self.ds)
647 def rmatvec(self, v):
648 """Evaluate w = M^H v"""
649 if self.collapsed is not None:
650 return np.dot(self.collapsed.T.conj(), v)
651 return LowRankMatrix._matvec(v, np.conj(self.alpha), self.ds, self.cs)
653 def solve(self, v, tol=0):
654 """Evaluate w = M^-1 v"""
655 if self.collapsed is not None:
656 return solve(self.collapsed, v)
657 return LowRankMatrix._solve(v, self.alpha, self.cs, self.ds)
659 def rsolve(self, v, tol=0):
660 """Evaluate w = M^-H v"""
661 if self.collapsed is not None:
662 return solve(self.collapsed.T.conj(), v)
663 return LowRankMatrix._solve(v, np.conj(self.alpha), self.ds, self.cs)
665 def append(self, c, d):
666 if self.collapsed is not None:
667 self.collapsed += c[:,None] * d[None,:].conj()
668 return
670 self.cs.append(c)
671 self.ds.append(d)
673 if len(self.cs) > c.size:
674 self.collapse()
676 def __array__(self):
677 if self.collapsed is not None:
678 return self.collapsed
680 Gm = self.alpha*np.identity(self.n, dtype=self.dtype)
681 for c, d in zip(self.cs, self.ds):
682 Gm += c[:,None]*d[None,:].conj()
683 return Gm
685 def collapse(self):
686 """Collapse the low-rank matrix to a full-rank one."""
687 self.collapsed = np.array(self)
688 self.cs = None
689 self.ds = None
690 self.alpha = None
692 def restart_reduce(self, rank):
693 """
694 Reduce the rank of the matrix by dropping all vectors.
695 """
696 if self.collapsed is not None:
697 return
698 assert rank > 0
699 if len(self.cs) > rank:
700 del self.cs[:]
701 del self.ds[:]
703 def simple_reduce(self, rank):
704 """
705 Reduce the rank of the matrix by dropping oldest vectors.
706 """
707 if self.collapsed is not None:
708 return
709 assert rank > 0
710 while len(self.cs) > rank:
711 del self.cs[0]
712 del self.ds[0]
714 def svd_reduce(self, max_rank, to_retain=None):
715 """
716 Reduce the rank of the matrix by retaining some SVD components.
718 This corresponds to the \"Broyden Rank Reduction Inverse\"
719 algorithm described in [1]_.
721 Note that the SVD decomposition can be done by solving only a
722 problem whose size is the effective rank of this matrix, which
723 is viable even for large problems.
725 Parameters
726 ----------
727 max_rank : int
728 Maximum rank of this matrix after reduction.
729 to_retain : int, optional
730 Number of SVD components to retain when reduction is done
731 (ie. rank > max_rank). Default is ``max_rank - 2``.
733 References
734 ----------
735 .. [1] B.A. van der Rotten, PhD thesis,
736 \"A limited memory Broyden method to solve high-dimensional
737 systems of nonlinear equations\". Mathematisch Instituut,
738 Universiteit Leiden, The Netherlands (2003).
740 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
742 """
743 if self.collapsed is not None:
744 return
746 p = max_rank
747 if to_retain is not None:
748 q = to_retain
749 else:
750 q = p - 2
752 if self.cs:
753 p = min(p, len(self.cs[0]))
754 q = max(0, min(q, p-1))
756 m = len(self.cs)
757 if m < p:
758 # nothing to do
759 return
761 C = np.array(self.cs).T
762 D = np.array(self.ds).T
764 D, R = qr(D, mode='economic')
765 C = dot(C, R.T.conj())
767 U, S, WH = svd(C, full_matrices=False)
769 C = dot(C, inv(WH))
770 D = dot(D, WH.T.conj())
772 for k in range(q):
773 self.cs[k] = C[:,k].copy()
774 self.ds[k] = D[:,k].copy()
776 del self.cs[q:]
777 del self.ds[q:]
780_doc_parts['broyden_params'] = """
781 alpha : float, optional
782 Initial guess for the Jacobian is ``(-1/alpha)``.
783 reduction_method : str or tuple, optional
784 Method used in ensuring that the rank of the Broyden matrix
785 stays low. Can either be a string giving the name of the method,
786 or a tuple of the form ``(method, param1, param2, ...)``
787 that gives the name of the method and values for additional parameters.
789 Methods available:
791 - ``restart``: drop all matrix columns. Has no extra parameters.
792 - ``simple``: drop oldest matrix column. Has no extra parameters.
793 - ``svd``: keep only the most significant SVD components.
794 Takes an extra parameter, ``to_retain``, which determines the
795 number of SVD components to retain when rank reduction is done.
796 Default is ``max_rank - 2``.
798 max_rank : int, optional
799 Maximum rank for the Broyden matrix.
800 Default is infinity (i.e., no rank reduction).
801 """.strip()
804class BroydenFirst(GenericBroyden):
805 r"""
806 Find a root of a function, using Broyden's first Jacobian approximation.
808 This method is also known as \"Broyden's good method\".
810 Parameters
811 ----------
812 %(params_basic)s
813 %(broyden_params)s
814 %(params_extra)s
816 See Also
817 --------
818 root : Interface to root finding algorithms for multivariate
819 functions. See ``method='broyden1'`` in particular.
821 Notes
822 -----
823 This algorithm implements the inverse Jacobian Quasi-Newton update
825 .. math:: H_+ = H + (dx - H df) dx^\dagger H / ( dx^\dagger H df)
827 which corresponds to Broyden's first Jacobian update
829 .. math:: J_+ = J + (df - J dx) dx^\dagger / dx^\dagger dx
832 References
833 ----------
834 .. [1] B.A. van der Rotten, PhD thesis,
835 \"A limited memory Broyden method to solve high-dimensional
836 systems of nonlinear equations\". Mathematisch Instituut,
837 Universiteit Leiden, The Netherlands (2003).
839 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
841 Examples
842 --------
843 The following functions define a system of nonlinear equations
845 >>> def fun(x):
846 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
847 ... 0.5 * (x[1] - x[0])**3 + x[1]]
849 A solution can be obtained as follows.
851 >>> from scipy import optimize
852 >>> sol = optimize.broyden1(fun, [0, 0])
853 >>> sol
854 array([0.84116396, 0.15883641])
856 """
858 def __init__(self, alpha=None, reduction_method='restart', max_rank=None):
859 GenericBroyden.__init__(self)
860 self.alpha = alpha
861 self.Gm = None
863 if max_rank is None:
864 max_rank = np.inf
865 self.max_rank = max_rank
867 if isinstance(reduction_method, str):
868 reduce_params = ()
869 else:
870 reduce_params = reduction_method[1:]
871 reduction_method = reduction_method[0]
872 reduce_params = (max_rank - 1,) + reduce_params
874 if reduction_method == 'svd':
875 self._reduce = lambda: self.Gm.svd_reduce(*reduce_params)
876 elif reduction_method == 'simple':
877 self._reduce = lambda: self.Gm.simple_reduce(*reduce_params)
878 elif reduction_method == 'restart':
879 self._reduce = lambda: self.Gm.restart_reduce(*reduce_params)
880 else:
881 raise ValueError("Unknown rank reduction method '%s'" %
882 reduction_method)
884 def setup(self, x, F, func):
885 GenericBroyden.setup(self, x, F, func)
886 self.Gm = LowRankMatrix(-self.alpha, self.shape[0], self.dtype)
888 def todense(self):
889 return inv(self.Gm)
891 def solve(self, f, tol=0):
892 r = self.Gm.matvec(f)
893 if not np.isfinite(r).all():
894 # singular; reset the Jacobian approximation
895 self.setup(self.last_x, self.last_f, self.func)
896 return self.Gm.matvec(f)
897 return r
899 def matvec(self, f):
900 return self.Gm.solve(f)
902 def rsolve(self, f, tol=0):
903 return self.Gm.rmatvec(f)
905 def rmatvec(self, f):
906 return self.Gm.rsolve(f)
908 def _update(self, x, f, dx, df, dx_norm, df_norm):
909 self._reduce() # reduce first to preserve secant condition
911 v = self.Gm.rmatvec(dx)
912 c = dx - self.Gm.matvec(df)
913 d = v / vdot(df, v)
915 self.Gm.append(c, d)
918class BroydenSecond(BroydenFirst):
919 """
920 Find a root of a function, using Broyden\'s second Jacobian approximation.
922 This method is also known as \"Broyden's bad method\".
924 Parameters
925 ----------
926 %(params_basic)s
927 %(broyden_params)s
928 %(params_extra)s
930 See Also
931 --------
932 root : Interface to root finding algorithms for multivariate
933 functions. See ``method='broyden2'`` in particular.
935 Notes
936 -----
937 This algorithm implements the inverse Jacobian Quasi-Newton update
939 .. math:: H_+ = H + (dx - H df) df^\\dagger / ( df^\\dagger df)
941 corresponding to Broyden's second method.
943 References
944 ----------
945 .. [1] B.A. van der Rotten, PhD thesis,
946 \"A limited memory Broyden method to solve high-dimensional
947 systems of nonlinear equations\". Mathematisch Instituut,
948 Universiteit Leiden, The Netherlands (2003).
950 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
952 Examples
953 --------
954 The following functions define a system of nonlinear equations
956 >>> def fun(x):
957 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
958 ... 0.5 * (x[1] - x[0])**3 + x[1]]
960 A solution can be obtained as follows.
962 >>> from scipy import optimize
963 >>> sol = optimize.broyden2(fun, [0, 0])
964 >>> sol
965 array([0.84116365, 0.15883529])
967 """
969 def _update(self, x, f, dx, df, dx_norm, df_norm):
970 self._reduce() # reduce first to preserve secant condition
972 v = df
973 c = dx - self.Gm.matvec(df)
974 d = v / df_norm**2
975 self.Gm.append(c, d)
978#------------------------------------------------------------------------------
979# Broyden-like (restricted memory)
980#------------------------------------------------------------------------------
982class Anderson(GenericBroyden):
983 """
984 Find a root of a function, using (extended) Anderson mixing.
986 The Jacobian is formed by for a 'best' solution in the space
987 spanned by last `M` vectors. As a result, only a MxM matrix
988 inversions and MxN multiplications are required. [Ey]_
990 Parameters
991 ----------
992 %(params_basic)s
993 alpha : float, optional
994 Initial guess for the Jacobian is (-1/alpha).
995 M : float, optional
996 Number of previous vectors to retain. Defaults to 5.
997 w0 : float, optional
998 Regularization parameter for numerical stability.
999 Compared to unity, good values of the order of 0.01.
1000 %(params_extra)s
1002 See Also
1003 --------
1004 root : Interface to root finding algorithms for multivariate
1005 functions. See ``method='anderson'`` in particular.
1007 References
1008 ----------
1009 .. [Ey] V. Eyert, J. Comp. Phys., 124, 271 (1996).
1011 Examples
1012 --------
1013 The following functions define a system of nonlinear equations
1015 >>> def fun(x):
1016 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
1017 ... 0.5 * (x[1] - x[0])**3 + x[1]]
1019 A solution can be obtained as follows.
1021 >>> from scipy import optimize
1022 >>> sol = optimize.anderson(fun, [0, 0])
1023 >>> sol
1024 array([0.84116588, 0.15883789])
1026 """
1028 # Note:
1029 #
1030 # Anderson method maintains a rank M approximation of the inverse Jacobian,
1031 #
1032 # J^-1 v ~ -v*alpha + (dX + alpha dF) A^-1 dF^H v
1033 # A = W + dF^H dF
1034 # W = w0^2 diag(dF^H dF)
1035 #
1036 # so that for w0 = 0 the secant condition applies for last M iterates, i.e.,
1037 #
1038 # J^-1 df_j = dx_j
1039 #
1040 # for all j = 0 ... M-1.
1041 #
1042 # Moreover, (from Sherman-Morrison-Woodbury formula)
1043 #
1044 # J v ~ [ b I - b^2 C (I + b dF^H A^-1 C)^-1 dF^H ] v
1045 # C = (dX + alpha dF) A^-1
1046 # b = -1/alpha
1047 #
1048 # and after simplification
1049 #
1050 # J v ~ -v/alpha + (dX/alpha + dF) (dF^H dX - alpha W)^-1 dF^H v
1051 #
1053 def __init__(self, alpha=None, w0=0.01, M=5):
1054 GenericBroyden.__init__(self)
1055 self.alpha = alpha
1056 self.M = M
1057 self.dx = []
1058 self.df = []
1059 self.gamma = None
1060 self.w0 = w0
1062 def solve(self, f, tol=0):
1063 dx = -self.alpha*f
1065 n = len(self.dx)
1066 if n == 0:
1067 return dx
1069 df_f = np.empty(n, dtype=f.dtype)
1070 for k in range(n):
1071 df_f[k] = vdot(self.df[k], f)
1073 try:
1074 gamma = solve(self.a, df_f)
1075 except LinAlgError:
1076 # singular; reset the Jacobian approximation
1077 del self.dx[:]
1078 del self.df[:]
1079 return dx
1081 for m in range(n):
1082 dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m])
1083 return dx
1085 def matvec(self, f):
1086 dx = -f/self.alpha
1088 n = len(self.dx)
1089 if n == 0:
1090 return dx
1092 df_f = np.empty(n, dtype=f.dtype)
1093 for k in range(n):
1094 df_f[k] = vdot(self.df[k], f)
1096 b = np.empty((n, n), dtype=f.dtype)
1097 for i in range(n):
1098 for j in range(n):
1099 b[i,j] = vdot(self.df[i], self.dx[j])
1100 if i == j and self.w0 != 0:
1101 b[i,j] -= vdot(self.df[i], self.df[i])*self.w0**2*self.alpha
1102 gamma = solve(b, df_f)
1104 for m in range(n):
1105 dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha)
1106 return dx
1108 def _update(self, x, f, dx, df, dx_norm, df_norm):
1109 if self.M == 0:
1110 return
1112 self.dx.append(dx)
1113 self.df.append(df)
1115 while len(self.dx) > self.M:
1116 self.dx.pop(0)
1117 self.df.pop(0)
1119 n = len(self.dx)
1120 a = np.zeros((n, n), dtype=f.dtype)
1122 for i in range(n):
1123 for j in range(i, n):
1124 if i == j:
1125 wd = self.w0**2
1126 else:
1127 wd = 0
1128 a[i,j] = (1+wd)*vdot(self.df[i], self.df[j])
1130 a += np.triu(a, 1).T.conj()
1131 self.a = a
1133#------------------------------------------------------------------------------
1134# Simple iterations
1135#------------------------------------------------------------------------------
1138class DiagBroyden(GenericBroyden):
1139 """
1140 Find a root of a function, using diagonal Broyden Jacobian approximation.
1142 The Jacobian approximation is derived from previous iterations, by
1143 retaining only the diagonal of Broyden matrices.
1145 .. warning::
1147 This algorithm may be useful for specific problems, but whether
1148 it will work may depend strongly on the problem.
1150 Parameters
1151 ----------
1152 %(params_basic)s
1153 alpha : float, optional
1154 Initial guess for the Jacobian is (-1/alpha).
1155 %(params_extra)s
1157 See Also
1158 --------
1159 root : Interface to root finding algorithms for multivariate
1160 functions. See ``method='diagbroyden'`` in particular.
1162 Examples
1163 --------
1164 The following functions define a system of nonlinear equations
1166 >>> def fun(x):
1167 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
1168 ... 0.5 * (x[1] - x[0])**3 + x[1]]
1170 A solution can be obtained as follows.
1172 >>> from scipy import optimize
1173 >>> sol = optimize.diagbroyden(fun, [0, 0])
1174 >>> sol
1175 array([0.84116403, 0.15883384])
1177 """
1179 def __init__(self, alpha=None):
1180 GenericBroyden.__init__(self)
1181 self.alpha = alpha
1183 def setup(self, x, F, func):
1184 GenericBroyden.setup(self, x, F, func)
1185 self.d = np.full((self.shape[0],), 1 / self.alpha, dtype=self.dtype)
1187 def solve(self, f, tol=0):
1188 return -f / self.d
1190 def matvec(self, f):
1191 return -f * self.d
1193 def rsolve(self, f, tol=0):
1194 return -f / self.d.conj()
1196 def rmatvec(self, f):
1197 return -f * self.d.conj()
1199 def todense(self):
1200 return np.diag(-self.d)
1202 def _update(self, x, f, dx, df, dx_norm, df_norm):
1203 self.d -= (df + self.d*dx)*dx/dx_norm**2
1206class LinearMixing(GenericBroyden):
1207 """
1208 Find a root of a function, using a scalar Jacobian approximation.
1210 .. warning::
1212 This algorithm may be useful for specific problems, but whether
1213 it will work may depend strongly on the problem.
1215 Parameters
1216 ----------
1217 %(params_basic)s
1218 alpha : float, optional
1219 The Jacobian approximation is (-1/alpha).
1220 %(params_extra)s
1222 See Also
1223 --------
1224 root : Interface to root finding algorithms for multivariate
1225 functions. See ``method='linearmixing'`` in particular.
1227 """
1229 def __init__(self, alpha=None):
1230 GenericBroyden.__init__(self)
1231 self.alpha = alpha
1233 def solve(self, f, tol=0):
1234 return -f*self.alpha
1236 def matvec(self, f):
1237 return -f/self.alpha
1239 def rsolve(self, f, tol=0):
1240 return -f*np.conj(self.alpha)
1242 def rmatvec(self, f):
1243 return -f/np.conj(self.alpha)
1245 def todense(self):
1246 return np.diag(np.full(self.shape[0], -1/self.alpha))
1248 def _update(self, x, f, dx, df, dx_norm, df_norm):
1249 pass
1252class ExcitingMixing(GenericBroyden):
1253 """
1254 Find a root of a function, using a tuned diagonal Jacobian approximation.
1256 The Jacobian matrix is diagonal and is tuned on each iteration.
1258 .. warning::
1260 This algorithm may be useful for specific problems, but whether
1261 it will work may depend strongly on the problem.
1263 See Also
1264 --------
1265 root : Interface to root finding algorithms for multivariate
1266 functions. See ``method='excitingmixing'`` in particular.
1268 Parameters
1269 ----------
1270 %(params_basic)s
1271 alpha : float, optional
1272 Initial Jacobian approximation is (-1/alpha).
1273 alphamax : float, optional
1274 The entries of the diagonal Jacobian are kept in the range
1275 ``[alpha, alphamax]``.
1276 %(params_extra)s
1277 """
1279 def __init__(self, alpha=None, alphamax=1.0):
1280 GenericBroyden.__init__(self)
1281 self.alpha = alpha
1282 self.alphamax = alphamax
1283 self.beta = None
1285 def setup(self, x, F, func):
1286 GenericBroyden.setup(self, x, F, func)
1287 self.beta = np.full((self.shape[0],), self.alpha, dtype=self.dtype)
1289 def solve(self, f, tol=0):
1290 return -f*self.beta
1292 def matvec(self, f):
1293 return -f/self.beta
1295 def rsolve(self, f, tol=0):
1296 return -f*self.beta.conj()
1298 def rmatvec(self, f):
1299 return -f/self.beta.conj()
1301 def todense(self):
1302 return np.diag(-1/self.beta)
1304 def _update(self, x, f, dx, df, dx_norm, df_norm):
1305 incr = f*self.last_f > 0
1306 self.beta[incr] += self.alpha
1307 self.beta[~incr] = self.alpha
1308 np.clip(self.beta, 0, self.alphamax, out=self.beta)
1311#------------------------------------------------------------------------------
1312# Iterative/Krylov approximated Jacobians
1313#------------------------------------------------------------------------------
1315class KrylovJacobian(Jacobian):
1316 r"""
1317 Find a root of a function, using Krylov approximation for inverse Jacobian.
1319 This method is suitable for solving large-scale problems.
1321 Parameters
1322 ----------
1323 %(params_basic)s
1324 rdiff : float, optional
1325 Relative step size to use in numerical differentiation.
1326 method : str or callable, optional
1327 Krylov method to use to approximate the Jacobian. Can be a string,
1328 or a function implementing the same interface as the iterative
1329 solvers in `scipy.sparse.linalg`. If a string, needs to be one of:
1330 ``'lgmres'``, ``'gmres'``, ``'bicgstab'``, ``'cgs'``, ``'minres'``,
1331 ``'tfqmr'``.
1333 The default is `scipy.sparse.linalg.lgmres`.
1334 inner_maxiter : int, optional
1335 Parameter to pass to the "inner" Krylov solver: maximum number of
1336 iterations. Iteration will stop after maxiter steps even if the
1337 specified tolerance has not been achieved.
1338 inner_M : LinearOperator or InverseJacobian
1339 Preconditioner for the inner Krylov iteration.
1340 Note that you can use also inverse Jacobians as (adaptive)
1341 preconditioners. For example,
1343 >>> from scipy.optimize import BroydenFirst, KrylovJacobian
1344 >>> from scipy.optimize import InverseJacobian
1345 >>> jac = BroydenFirst()
1346 >>> kjac = KrylovJacobian(inner_M=InverseJacobian(jac))
1348 If the preconditioner has a method named 'update', it will be called
1349 as ``update(x, f)`` after each nonlinear step, with ``x`` giving
1350 the current point, and ``f`` the current function value.
1351 outer_k : int, optional
1352 Size of the subspace kept across LGMRES nonlinear iterations.
1353 See `scipy.sparse.linalg.lgmres` for details.
1354 inner_kwargs : kwargs
1355 Keyword parameters for the "inner" Krylov solver
1356 (defined with `method`). Parameter names must start with
1357 the `inner_` prefix which will be stripped before passing on
1358 the inner method. See, e.g., `scipy.sparse.linalg.gmres` for details.
1359 %(params_extra)s
1361 See Also
1362 --------
1363 root : Interface to root finding algorithms for multivariate
1364 functions. See ``method='krylov'`` in particular.
1365 scipy.sparse.linalg.gmres
1366 scipy.sparse.linalg.lgmres
1368 Notes
1369 -----
1370 This function implements a Newton-Krylov solver. The basic idea is
1371 to compute the inverse of the Jacobian with an iterative Krylov
1372 method. These methods require only evaluating the Jacobian-vector
1373 products, which are conveniently approximated by a finite difference:
1375 .. math:: J v \approx (f(x + \omega*v/|v|) - f(x)) / \omega
1377 Due to the use of iterative matrix inverses, these methods can
1378 deal with large nonlinear problems.
1380 SciPy's `scipy.sparse.linalg` module offers a selection of Krylov
1381 solvers to choose from. The default here is `lgmres`, which is a
1382 variant of restarted GMRES iteration that reuses some of the
1383 information obtained in the previous Newton steps to invert
1384 Jacobians in subsequent steps.
1386 For a review on Newton-Krylov methods, see for example [1]_,
1387 and for the LGMRES sparse inverse method, see [2]_.
1389 References
1390 ----------
1391 .. [1] C. T. Kelley, Solving Nonlinear Equations with Newton's Method,
1392 SIAM, pp.57-83, 2003.
1393 :doi:`10.1137/1.9780898718898.ch3`
1394 .. [2] D.A. Knoll and D.E. Keyes, J. Comp. Phys. 193, 357 (2004).
1395 :doi:`10.1016/j.jcp.2003.08.010`
1396 .. [3] A.H. Baker and E.R. Jessup and T. Manteuffel,
1397 SIAM J. Matrix Anal. Appl. 26, 962 (2005).
1398 :doi:`10.1137/S0895479803422014`
1400 Examples
1401 --------
1402 The following functions define a system of nonlinear equations
1404 >>> def fun(x):
1405 ... return [x[0] + 0.5 * x[1] - 1.0,
1406 ... 0.5 * (x[1] - x[0]) ** 2]
1408 A solution can be obtained as follows.
1410 >>> from scipy import optimize
1411 >>> sol = optimize.newton_krylov(fun, [0, 0])
1412 >>> sol
1413 array([0.66731771, 0.66536458])
1415 """
1417 def __init__(self, rdiff=None, method='lgmres', inner_maxiter=20,
1418 inner_M=None, outer_k=10, **kw):
1419 self.preconditioner = inner_M
1420 self.rdiff = rdiff
1421 # Note that this retrieves one of the named functions, or otherwise
1422 # uses `method` as is (i.e., for a user-provided callable).
1423 self.method = dict(
1424 bicgstab=scipy.sparse.linalg.bicgstab,
1425 gmres=scipy.sparse.linalg.gmres,
1426 lgmres=scipy.sparse.linalg.lgmres,
1427 cgs=scipy.sparse.linalg.cgs,
1428 minres=scipy.sparse.linalg.minres,
1429 tfqmr=scipy.sparse.linalg.tfqmr,
1430 ).get(method, method)
1432 self.method_kw = dict(maxiter=inner_maxiter, M=self.preconditioner)
1434 if self.method is scipy.sparse.linalg.gmres:
1435 # Replace GMRES's outer iteration with Newton steps
1436 self.method_kw['restart'] = inner_maxiter
1437 self.method_kw['maxiter'] = 1
1438 self.method_kw.setdefault('atol', 0)
1439 elif self.method in (scipy.sparse.linalg.gcrotmk,
1440 scipy.sparse.linalg.bicgstab,
1441 scipy.sparse.linalg.cgs):
1442 self.method_kw.setdefault('atol', 0)
1443 elif self.method is scipy.sparse.linalg.lgmres:
1444 self.method_kw['outer_k'] = outer_k
1445 # Replace LGMRES's outer iteration with Newton steps
1446 self.method_kw['maxiter'] = 1
1447 # Carry LGMRES's `outer_v` vectors across nonlinear iterations
1448 self.method_kw.setdefault('outer_v', [])
1449 self.method_kw.setdefault('prepend_outer_v', True)
1450 # But don't carry the corresponding Jacobian*v products, in case
1451 # the Jacobian changes a lot in the nonlinear step
1452 #
1453 # XXX: some trust-region inspired ideas might be more efficient...
1454 # See e.g., Brown & Saad. But needs to be implemented separately
1455 # since it's not an inexact Newton method.
1456 self.method_kw.setdefault('store_outer_Av', False)
1457 self.method_kw.setdefault('atol', 0)
1459 for key, value in kw.items():
1460 if not key.startswith('inner_'):
1461 raise ValueError("Unknown parameter %s" % key)
1462 self.method_kw[key[6:]] = value
1464 def _update_diff_step(self):
1465 mx = abs(self.x0).max()
1466 mf = abs(self.f0).max()
1467 self.omega = self.rdiff * max(1, mx) / max(1, mf)
1469 def matvec(self, v):
1470 nv = norm(v)
1471 if nv == 0:
1472 return 0*v
1473 sc = self.omega / nv
1474 r = (self.func(self.x0 + sc*v) - self.f0) / sc
1475 if not np.all(np.isfinite(r)) and np.all(np.isfinite(v)):
1476 raise ValueError('Function returned non-finite results')
1477 return r
1479 def solve(self, rhs, tol=0):
1480 if 'tol' in self.method_kw:
1481 sol, info = self.method(self.op, rhs, **self.method_kw)
1482 else:
1483 sol, info = self.method(self.op, rhs, tol=tol, **self.method_kw)
1484 return sol
1486 def update(self, x, f):
1487 self.x0 = x
1488 self.f0 = f
1489 self._update_diff_step()
1491 # Update also the preconditioner, if possible
1492 if self.preconditioner is not None:
1493 if hasattr(self.preconditioner, 'update'):
1494 self.preconditioner.update(x, f)
1496 def setup(self, x, f, func):
1497 Jacobian.setup(self, x, f, func)
1498 self.x0 = x
1499 self.f0 = f
1500 self.op = scipy.sparse.linalg.aslinearoperator(self)
1502 if self.rdiff is None:
1503 self.rdiff = np.finfo(x.dtype).eps ** (1./2)
1505 self._update_diff_step()
1507 # Setup also the preconditioner, if possible
1508 if self.preconditioner is not None:
1509 if hasattr(self.preconditioner, 'setup'):
1510 self.preconditioner.setup(x, f, func)
1513#------------------------------------------------------------------------------
1514# Wrapper functions
1515#------------------------------------------------------------------------------
1517def _nonlin_wrapper(name, jac):
1518 """
1519 Construct a solver wrapper with given name and Jacobian approx.
1521 It inspects the keyword arguments of ``jac.__init__``, and allows to
1522 use the same arguments in the wrapper function, in addition to the
1523 keyword arguments of `nonlin_solve`
1525 """
1526 signature = _getfullargspec(jac.__init__)
1527 args, varargs, varkw, defaults, kwonlyargs, kwdefaults, _ = signature
1528 kwargs = list(zip(args[-len(defaults):], defaults))
1529 kw_str = ", ".join([f"{k}={v!r}" for k, v in kwargs])
1530 if kw_str:
1531 kw_str = ", " + kw_str
1532 kwkw_str = ", ".join([f"{k}={k}" for k, v in kwargs])
1533 if kwkw_str:
1534 kwkw_str = kwkw_str + ", "
1535 if kwonlyargs:
1536 raise ValueError('Unexpected signature %s' % signature)
1538 # Construct the wrapper function so that its keyword arguments
1539 # are visible in pydoc.help etc.
1540 wrapper = """
1541def %(name)s(F, xin, iter=None %(kw)s, verbose=False, maxiter=None,
1542 f_tol=None, f_rtol=None, x_tol=None, x_rtol=None,
1543 tol_norm=None, line_search='armijo', callback=None, **kw):
1544 jac = %(jac)s(%(kwkw)s **kw)
1545 return nonlin_solve(F, xin, jac, iter, verbose, maxiter,
1546 f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search,
1547 callback)
1548"""
1550 wrapper = wrapper % dict(name=name, kw=kw_str, jac=jac.__name__,
1551 kwkw=kwkw_str)
1552 ns = {}
1553 ns.update(globals())
1554 exec(wrapper, ns)
1555 func = ns[name]
1556 func.__doc__ = jac.__doc__
1557 _set_doc(func)
1558 return func
1561broyden1 = _nonlin_wrapper('broyden1', BroydenFirst)
1562broyden2 = _nonlin_wrapper('broyden2', BroydenSecond)
1563anderson = _nonlin_wrapper('anderson', Anderson)
1564linearmixing = _nonlin_wrapper('linearmixing', LinearMixing)
1565diagbroyden = _nonlin_wrapper('diagbroyden', DiagBroyden)
1566excitingmixing = _nonlin_wrapper('excitingmixing', ExcitingMixing)
1567newton_krylov = _nonlin_wrapper('newton_krylov', KrylovJacobian)