Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_nonlin.py: 22%
630 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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 func = lambda z: _as_inexact(F(_array_like(z, x0))).flatten()
168 x = x0.flatten()
170 dx = np.full_like(x, np.inf)
171 Fx = func(x)
172 Fx_norm = norm(Fx)
174 jacobian = asjacobian(jacobian)
175 jacobian.setup(x.copy(), Fx, func)
177 if maxiter is None:
178 if iter is not None:
179 maxiter = iter + 1
180 else:
181 maxiter = 100*(x.size+1)
183 if line_search is True:
184 line_search = 'armijo'
185 elif line_search is False:
186 line_search = None
188 if line_search not in (None, 'armijo', 'wolfe'):
189 raise ValueError("Invalid line search")
191 # Solver tolerance selection
192 gamma = 0.9
193 eta_max = 0.9999
194 eta_treshold = 0.1
195 eta = 1e-3
197 for n in range(maxiter):
198 status = condition.check(Fx, x, dx)
199 if status:
200 break
202 # The tolerance, as computed for scipy.sparse.linalg.* routines
203 tol = min(eta, eta*Fx_norm)
204 dx = -jacobian.solve(Fx, tol=tol)
206 if norm(dx) == 0:
207 raise ValueError("Jacobian inversion yielded zero vector. "
208 "This indicates a bug in the Jacobian "
209 "approximation.")
211 # Line search, or Newton step
212 if line_search:
213 s, x, Fx, Fx_norm_new = _nonlin_line_search(func, x, Fx, dx,
214 line_search)
215 else:
216 s = 1.0
217 x = x + dx
218 Fx = func(x)
219 Fx_norm_new = norm(Fx)
221 jacobian.update(x.copy(), Fx)
223 if callback:
224 callback(x, Fx)
226 # Adjust forcing parameters for inexact methods
227 eta_A = gamma * Fx_norm_new**2 / Fx_norm**2
228 if gamma * eta**2 < eta_treshold:
229 eta = min(eta_max, eta_A)
230 else:
231 eta = min(eta_max, max(eta_A, gamma*eta**2))
233 Fx_norm = Fx_norm_new
235 # Print status
236 if verbose:
237 sys.stdout.write("%d: |F(x)| = %g; step %g\n" % (
238 n, tol_norm(Fx), s))
239 sys.stdout.flush()
240 else:
241 if raise_exception:
242 raise NoConvergence(_array_like(x, x0))
243 else:
244 status = 2
246 if full_output:
247 info = {'nit': condition.iteration,
248 'fun': Fx,
249 'status': status,
250 'success': status == 1,
251 'message': {1: 'A solution was found at the specified '
252 'tolerance.',
253 2: 'The maximum number of iterations allowed '
254 'has been reached.'
255 }[status]
256 }
257 return _array_like(x, x0), info
258 else:
259 return _array_like(x, x0)
262_set_doc(nonlin_solve)
265def _nonlin_line_search(func, x, Fx, dx, search_type='armijo', rdiff=1e-8,
266 smin=1e-2):
267 tmp_s = [0]
268 tmp_Fx = [Fx]
269 tmp_phi = [norm(Fx)**2]
270 s_norm = norm(x) / norm(dx)
272 def phi(s, store=True):
273 if s == tmp_s[0]:
274 return tmp_phi[0]
275 xt = x + s*dx
276 v = func(xt)
277 p = _safe_norm(v)**2
278 if store:
279 tmp_s[0] = s
280 tmp_phi[0] = p
281 tmp_Fx[0] = v
282 return p
284 def derphi(s):
285 ds = (abs(s) + s_norm + 1) * rdiff
286 return (phi(s+ds, store=False) - phi(s)) / ds
288 if search_type == 'wolfe':
289 s, phi1, phi0 = scalar_search_wolfe1(phi, derphi, tmp_phi[0],
290 xtol=1e-2, amin=smin)
291 elif search_type == 'armijo':
292 s, phi1 = scalar_search_armijo(phi, tmp_phi[0], -tmp_phi[0],
293 amin=smin)
295 if s is None:
296 # XXX: No suitable step length found. Take the full Newton step,
297 # and hope for the best.
298 s = 1.0
300 x = x + s*dx
301 if s == tmp_s[0]:
302 Fx = tmp_Fx[0]
303 else:
304 Fx = func(x)
305 Fx_norm = norm(Fx)
307 return s, x, Fx, Fx_norm
310class TerminationCondition:
311 """
312 Termination condition for an iteration. It is terminated if
314 - |F| < f_rtol*|F_0|, AND
315 - |F| < f_tol
317 AND
319 - |dx| < x_rtol*|x|, AND
320 - |dx| < x_tol
322 """
323 def __init__(self, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None,
324 iter=None, norm=maxnorm):
326 if f_tol is None:
327 f_tol = np.finfo(np.float_).eps ** (1./3)
328 if f_rtol is None:
329 f_rtol = np.inf
330 if x_tol is None:
331 x_tol = np.inf
332 if x_rtol is None:
333 x_rtol = np.inf
335 self.x_tol = x_tol
336 self.x_rtol = x_rtol
337 self.f_tol = f_tol
338 self.f_rtol = f_rtol
340 self.norm = norm
342 self.iter = iter
344 self.f0_norm = None
345 self.iteration = 0
347 def check(self, f, x, dx):
348 self.iteration += 1
349 f_norm = self.norm(f)
350 x_norm = self.norm(x)
351 dx_norm = self.norm(dx)
353 if self.f0_norm is None:
354 self.f0_norm = f_norm
356 if f_norm == 0:
357 return 1
359 if self.iter is not None:
360 # backwards compatibility with SciPy 0.6.0
361 return 2 * (self.iteration > self.iter)
363 # NB: condition must succeed for rtol=inf even if norm == 0
364 return int((f_norm <= self.f_tol
365 and f_norm/self.f_rtol <= self.f0_norm)
366 and (dx_norm <= self.x_tol
367 and dx_norm/self.x_rtol <= x_norm))
370#------------------------------------------------------------------------------
371# Generic Jacobian approximation
372#------------------------------------------------------------------------------
374class Jacobian:
375 """
376 Common interface for Jacobians or Jacobian approximations.
378 The optional methods come useful when implementing trust region
379 etc., algorithms that often require evaluating transposes of the
380 Jacobian.
382 Methods
383 -------
384 solve
385 Returns J^-1 * v
386 update
387 Updates Jacobian to point `x` (where the function has residual `Fx`)
389 matvec : optional
390 Returns J * v
391 rmatvec : optional
392 Returns A^H * v
393 rsolve : optional
394 Returns A^-H * v
395 matmat : optional
396 Returns A * V, where V is a dense matrix with dimensions (N,K).
397 todense : optional
398 Form the dense Jacobian matrix. Necessary for dense trust region
399 algorithms, and useful for testing.
401 Attributes
402 ----------
403 shape
404 Matrix dimensions (M, N)
405 dtype
406 Data type of the matrix.
407 func : callable, optional
408 Function the Jacobian corresponds to
410 """
412 def __init__(self, **kw):
413 names = ["solve", "update", "matvec", "rmatvec", "rsolve",
414 "matmat", "todense", "shape", "dtype"]
415 for name, value in kw.items():
416 if name not in names:
417 raise ValueError("Unknown keyword argument %s" % name)
418 if value is not None:
419 setattr(self, name, kw[name])
421 if hasattr(self, 'todense'):
422 self.__array__ = lambda: self.todense()
424 def aspreconditioner(self):
425 return InverseJacobian(self)
427 def solve(self, v, tol=0):
428 raise NotImplementedError
430 def update(self, x, F):
431 pass
433 def setup(self, x, F, func):
434 self.func = func
435 self.shape = (F.size, x.size)
436 self.dtype = F.dtype
437 if self.__class__.setup is Jacobian.setup:
438 # Call on the first point unless overridden
439 self.update(x, F)
442class InverseJacobian:
443 def __init__(self, jacobian):
444 self.jacobian = jacobian
445 self.matvec = jacobian.solve
446 self.update = jacobian.update
447 if hasattr(jacobian, 'setup'):
448 self.setup = jacobian.setup
449 if hasattr(jacobian, 'rsolve'):
450 self.rmatvec = jacobian.rsolve
452 @property
453 def shape(self):
454 return self.jacobian.shape
456 @property
457 def dtype(self):
458 return self.jacobian.dtype
461def asjacobian(J):
462 """
463 Convert given object to one suitable for use as a Jacobian.
464 """
465 spsolve = scipy.sparse.linalg.spsolve
466 if isinstance(J, Jacobian):
467 return J
468 elif inspect.isclass(J) and issubclass(J, Jacobian):
469 return J()
470 elif isinstance(J, np.ndarray):
471 if J.ndim > 2:
472 raise ValueError('array must have rank <= 2')
473 J = np.atleast_2d(np.asarray(J))
474 if J.shape[0] != J.shape[1]:
475 raise ValueError('array must be square')
477 return Jacobian(matvec=lambda v: dot(J, v),
478 rmatvec=lambda v: dot(J.conj().T, v),
479 solve=lambda v: solve(J, v),
480 rsolve=lambda v: solve(J.conj().T, v),
481 dtype=J.dtype, shape=J.shape)
482 elif scipy.sparse.isspmatrix(J):
483 if J.shape[0] != J.shape[1]:
484 raise ValueError('matrix must be square')
485 return Jacobian(matvec=lambda v: J*v,
486 rmatvec=lambda v: J.conj().T * v,
487 solve=lambda v: spsolve(J, v),
488 rsolve=lambda v: spsolve(J.conj().T, v),
489 dtype=J.dtype, shape=J.shape)
490 elif hasattr(J, 'shape') and hasattr(J, 'dtype') and hasattr(J, 'solve'):
491 return Jacobian(matvec=getattr(J, 'matvec'),
492 rmatvec=getattr(J, 'rmatvec'),
493 solve=J.solve,
494 rsolve=getattr(J, 'rsolve'),
495 update=getattr(J, 'update'),
496 setup=getattr(J, 'setup'),
497 dtype=J.dtype,
498 shape=J.shape)
499 elif callable(J):
500 # Assume it's a function J(x) that returns the Jacobian
501 class Jac(Jacobian):
502 def update(self, x, F):
503 self.x = x
505 def solve(self, v, tol=0):
506 m = J(self.x)
507 if isinstance(m, np.ndarray):
508 return solve(m, v)
509 elif scipy.sparse.isspmatrix(m):
510 return spsolve(m, v)
511 else:
512 raise ValueError("Unknown matrix type")
514 def matvec(self, v):
515 m = J(self.x)
516 if isinstance(m, np.ndarray):
517 return dot(m, v)
518 elif scipy.sparse.isspmatrix(m):
519 return m*v
520 else:
521 raise ValueError("Unknown matrix type")
523 def rsolve(self, v, tol=0):
524 m = J(self.x)
525 if isinstance(m, np.ndarray):
526 return solve(m.conj().T, v)
527 elif scipy.sparse.isspmatrix(m):
528 return spsolve(m.conj().T, v)
529 else:
530 raise ValueError("Unknown matrix type")
532 def rmatvec(self, v):
533 m = J(self.x)
534 if isinstance(m, np.ndarray):
535 return dot(m.conj().T, v)
536 elif scipy.sparse.isspmatrix(m):
537 return m.conj().T * v
538 else:
539 raise ValueError("Unknown matrix type")
540 return Jac()
541 elif isinstance(J, str):
542 return dict(broyden1=BroydenFirst,
543 broyden2=BroydenSecond,
544 anderson=Anderson,
545 diagbroyden=DiagBroyden,
546 linearmixing=LinearMixing,
547 excitingmixing=ExcitingMixing,
548 krylov=KrylovJacobian)[J]()
549 else:
550 raise TypeError('Cannot convert object to a Jacobian')
553#------------------------------------------------------------------------------
554# Broyden
555#------------------------------------------------------------------------------
557class GenericBroyden(Jacobian):
558 def setup(self, x0, f0, func):
559 Jacobian.setup(self, x0, f0, func)
560 self.last_f = f0
561 self.last_x = x0
563 if hasattr(self, 'alpha') and self.alpha is None:
564 # Autoscale the initial Jacobian parameter
565 # unless we have already guessed the solution.
566 normf0 = norm(f0)
567 if normf0:
568 self.alpha = 0.5*max(norm(x0), 1) / normf0
569 else:
570 self.alpha = 1.0
572 def _update(self, x, f, dx, df, dx_norm, df_norm):
573 raise NotImplementedError
575 def update(self, x, f):
576 df = f - self.last_f
577 dx = x - self.last_x
578 self._update(x, f, dx, df, norm(dx), norm(df))
579 self.last_f = f
580 self.last_x = x
583class LowRankMatrix:
584 r"""
585 A matrix represented as
587 .. math:: \alpha I + \sum_{n=0}^{n=M} c_n d_n^\dagger
589 However, if the rank of the matrix reaches the dimension of the vectors,
590 full matrix representation will be used thereon.
592 """
594 def __init__(self, alpha, n, dtype):
595 self.alpha = alpha
596 self.cs = []
597 self.ds = []
598 self.n = n
599 self.dtype = dtype
600 self.collapsed = None
602 @staticmethod
603 def _matvec(v, alpha, cs, ds):
604 axpy, scal, dotc = get_blas_funcs(['axpy', 'scal', 'dotc'],
605 cs[:1] + [v])
606 w = alpha * v
607 for c, d in zip(cs, ds):
608 a = dotc(d, v)
609 w = axpy(c, w, w.size, a)
610 return w
612 @staticmethod
613 def _solve(v, alpha, cs, ds):
614 """Evaluate w = M^-1 v"""
615 if len(cs) == 0:
616 return v/alpha
618 # (B + C D^H)^-1 = B^-1 - B^-1 C (I + D^H B^-1 C)^-1 D^H B^-1
620 axpy, dotc = get_blas_funcs(['axpy', 'dotc'], cs[:1] + [v])
622 c0 = cs[0]
623 A = alpha * np.identity(len(cs), dtype=c0.dtype)
624 for i, d in enumerate(ds):
625 for j, c in enumerate(cs):
626 A[i,j] += dotc(d, c)
628 q = np.zeros(len(cs), dtype=c0.dtype)
629 for j, d in enumerate(ds):
630 q[j] = dotc(d, v)
631 q /= alpha
632 q = solve(A, q)
634 w = v/alpha
635 for c, qc in zip(cs, q):
636 w = axpy(c, w, w.size, -qc)
638 return w
640 def matvec(self, v):
641 """Evaluate w = M v"""
642 if self.collapsed is not None:
643 return np.dot(self.collapsed, v)
644 return LowRankMatrix._matvec(v, self.alpha, self.cs, self.ds)
646 def rmatvec(self, v):
647 """Evaluate w = M^H v"""
648 if self.collapsed is not None:
649 return np.dot(self.collapsed.T.conj(), v)
650 return LowRankMatrix._matvec(v, np.conj(self.alpha), self.ds, self.cs)
652 def solve(self, v, tol=0):
653 """Evaluate w = M^-1 v"""
654 if self.collapsed is not None:
655 return solve(self.collapsed, v)
656 return LowRankMatrix._solve(v, self.alpha, self.cs, self.ds)
658 def rsolve(self, v, tol=0):
659 """Evaluate w = M^-H v"""
660 if self.collapsed is not None:
661 return solve(self.collapsed.T.conj(), v)
662 return LowRankMatrix._solve(v, np.conj(self.alpha), self.ds, self.cs)
664 def append(self, c, d):
665 if self.collapsed is not None:
666 self.collapsed += c[:,None] * d[None,:].conj()
667 return
669 self.cs.append(c)
670 self.ds.append(d)
672 if len(self.cs) > c.size:
673 self.collapse()
675 def __array__(self):
676 if self.collapsed is not None:
677 return self.collapsed
679 Gm = self.alpha*np.identity(self.n, dtype=self.dtype)
680 for c, d in zip(self.cs, self.ds):
681 Gm += c[:,None]*d[None,:].conj()
682 return Gm
684 def collapse(self):
685 """Collapse the low-rank matrix to a full-rank one."""
686 self.collapsed = np.array(self)
687 self.cs = None
688 self.ds = None
689 self.alpha = None
691 def restart_reduce(self, rank):
692 """
693 Reduce the rank of the matrix by dropping all vectors.
694 """
695 if self.collapsed is not None:
696 return
697 assert rank > 0
698 if len(self.cs) > rank:
699 del self.cs[:]
700 del self.ds[:]
702 def simple_reduce(self, rank):
703 """
704 Reduce the rank of the matrix by dropping oldest vectors.
705 """
706 if self.collapsed is not None:
707 return
708 assert rank > 0
709 while len(self.cs) > rank:
710 del self.cs[0]
711 del self.ds[0]
713 def svd_reduce(self, max_rank, to_retain=None):
714 """
715 Reduce the rank of the matrix by retaining some SVD components.
717 This corresponds to the \"Broyden Rank Reduction Inverse\"
718 algorithm described in [1]_.
720 Note that the SVD decomposition can be done by solving only a
721 problem whose size is the effective rank of this matrix, which
722 is viable even for large problems.
724 Parameters
725 ----------
726 max_rank : int
727 Maximum rank of this matrix after reduction.
728 to_retain : int, optional
729 Number of SVD components to retain when reduction is done
730 (ie. rank > max_rank). Default is ``max_rank - 2``.
732 References
733 ----------
734 .. [1] B.A. van der Rotten, PhD thesis,
735 \"A limited memory Broyden method to solve high-dimensional
736 systems of nonlinear equations\". Mathematisch Instituut,
737 Universiteit Leiden, The Netherlands (2003).
739 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
741 """
742 if self.collapsed is not None:
743 return
745 p = max_rank
746 if to_retain is not None:
747 q = to_retain
748 else:
749 q = p - 2
751 if self.cs:
752 p = min(p, len(self.cs[0]))
753 q = max(0, min(q, p-1))
755 m = len(self.cs)
756 if m < p:
757 # nothing to do
758 return
760 C = np.array(self.cs).T
761 D = np.array(self.ds).T
763 D, R = qr(D, mode='economic')
764 C = dot(C, R.T.conj())
766 U, S, WH = svd(C, full_matrices=False)
768 C = dot(C, inv(WH))
769 D = dot(D, WH.T.conj())
771 for k in range(q):
772 self.cs[k] = C[:,k].copy()
773 self.ds[k] = D[:,k].copy()
775 del self.cs[q:]
776 del self.ds[q:]
779_doc_parts['broyden_params'] = """
780 alpha : float, optional
781 Initial guess for the Jacobian is ``(-1/alpha)``.
782 reduction_method : str or tuple, optional
783 Method used in ensuring that the rank of the Broyden matrix
784 stays low. Can either be a string giving the name of the method,
785 or a tuple of the form ``(method, param1, param2, ...)``
786 that gives the name of the method and values for additional parameters.
788 Methods available:
790 - ``restart``: drop all matrix columns. Has no extra parameters.
791 - ``simple``: drop oldest matrix column. Has no extra parameters.
792 - ``svd``: keep only the most significant SVD components.
793 Takes an extra parameter, ``to_retain``, which determines the
794 number of SVD components to retain when rank reduction is done.
795 Default is ``max_rank - 2``.
797 max_rank : int, optional
798 Maximum rank for the Broyden matrix.
799 Default is infinity (i.e., no rank reduction).
800 """.strip()
803class BroydenFirst(GenericBroyden):
804 r"""
805 Find a root of a function, using Broyden's first Jacobian approximation.
807 This method is also known as \"Broyden's good method\".
809 Parameters
810 ----------
811 %(params_basic)s
812 %(broyden_params)s
813 %(params_extra)s
815 See Also
816 --------
817 root : Interface to root finding algorithms for multivariate
818 functions. See ``method='broyden1'`` in particular.
820 Notes
821 -----
822 This algorithm implements the inverse Jacobian Quasi-Newton update
824 .. math:: H_+ = H + (dx - H df) dx^\dagger H / ( dx^\dagger H df)
826 which corresponds to Broyden's first Jacobian update
828 .. math:: J_+ = J + (df - J dx) dx^\dagger / dx^\dagger dx
831 References
832 ----------
833 .. [1] B.A. van der Rotten, PhD thesis,
834 \"A limited memory Broyden method to solve high-dimensional
835 systems of nonlinear equations\". Mathematisch Instituut,
836 Universiteit Leiden, The Netherlands (2003).
838 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
840 Examples
841 --------
842 The following functions define a system of nonlinear equations
844 >>> def fun(x):
845 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
846 ... 0.5 * (x[1] - x[0])**3 + x[1]]
848 A solution can be obtained as follows.
850 >>> from scipy import optimize
851 >>> sol = optimize.broyden1(fun, [0, 0])
852 >>> sol
853 array([0.84116396, 0.15883641])
855 """
857 def __init__(self, alpha=None, reduction_method='restart', max_rank=None):
858 GenericBroyden.__init__(self)
859 self.alpha = alpha
860 self.Gm = None
862 if max_rank is None:
863 max_rank = np.inf
864 self.max_rank = max_rank
866 if isinstance(reduction_method, str):
867 reduce_params = ()
868 else:
869 reduce_params = reduction_method[1:]
870 reduction_method = reduction_method[0]
871 reduce_params = (max_rank - 1,) + reduce_params
873 if reduction_method == 'svd':
874 self._reduce = lambda: self.Gm.svd_reduce(*reduce_params)
875 elif reduction_method == 'simple':
876 self._reduce = lambda: self.Gm.simple_reduce(*reduce_params)
877 elif reduction_method == 'restart':
878 self._reduce = lambda: self.Gm.restart_reduce(*reduce_params)
879 else:
880 raise ValueError("Unknown rank reduction method '%s'" %
881 reduction_method)
883 def setup(self, x, F, func):
884 GenericBroyden.setup(self, x, F, func)
885 self.Gm = LowRankMatrix(-self.alpha, self.shape[0], self.dtype)
887 def todense(self):
888 return inv(self.Gm)
890 def solve(self, f, tol=0):
891 r = self.Gm.matvec(f)
892 if not np.isfinite(r).all():
893 # singular; reset the Jacobian approximation
894 self.setup(self.last_x, self.last_f, self.func)
895 return self.Gm.matvec(f)
896 return r
898 def matvec(self, f):
899 return self.Gm.solve(f)
901 def rsolve(self, f, tol=0):
902 return self.Gm.rmatvec(f)
904 def rmatvec(self, f):
905 return self.Gm.rsolve(f)
907 def _update(self, x, f, dx, df, dx_norm, df_norm):
908 self._reduce() # reduce first to preserve secant condition
910 v = self.Gm.rmatvec(dx)
911 c = dx - self.Gm.matvec(df)
912 d = v / vdot(df, v)
914 self.Gm.append(c, d)
917class BroydenSecond(BroydenFirst):
918 """
919 Find a root of a function, using Broyden\'s second Jacobian approximation.
921 This method is also known as \"Broyden's bad method\".
923 Parameters
924 ----------
925 %(params_basic)s
926 %(broyden_params)s
927 %(params_extra)s
929 See Also
930 --------
931 root : Interface to root finding algorithms for multivariate
932 functions. See ``method='broyden2'`` in particular.
934 Notes
935 -----
936 This algorithm implements the inverse Jacobian Quasi-Newton update
938 .. math:: H_+ = H + (dx - H df) df^\\dagger / ( df^\\dagger df)
940 corresponding to Broyden's second method.
942 References
943 ----------
944 .. [1] B.A. van der Rotten, PhD thesis,
945 \"A limited memory Broyden method to solve high-dimensional
946 systems of nonlinear equations\". Mathematisch Instituut,
947 Universiteit Leiden, The Netherlands (2003).
949 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
951 Examples
952 --------
953 The following functions define a system of nonlinear equations
955 >>> def fun(x):
956 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
957 ... 0.5 * (x[1] - x[0])**3 + x[1]]
959 A solution can be obtained as follows.
961 >>> from scipy import optimize
962 >>> sol = optimize.broyden2(fun, [0, 0])
963 >>> sol
964 array([0.84116365, 0.15883529])
966 """
968 def _update(self, x, f, dx, df, dx_norm, df_norm):
969 self._reduce() # reduce first to preserve secant condition
971 v = df
972 c = dx - self.Gm.matvec(df)
973 d = v / df_norm**2
974 self.Gm.append(c, d)
977#------------------------------------------------------------------------------
978# Broyden-like (restricted memory)
979#------------------------------------------------------------------------------
981class Anderson(GenericBroyden):
982 """
983 Find a root of a function, using (extended) Anderson mixing.
985 The Jacobian is formed by for a 'best' solution in the space
986 spanned by last `M` vectors. As a result, only a MxM matrix
987 inversions and MxN multiplications are required. [Ey]_
989 Parameters
990 ----------
991 %(params_basic)s
992 alpha : float, optional
993 Initial guess for the Jacobian is (-1/alpha).
994 M : float, optional
995 Number of previous vectors to retain. Defaults to 5.
996 w0 : float, optional
997 Regularization parameter for numerical stability.
998 Compared to unity, good values of the order of 0.01.
999 %(params_extra)s
1001 See Also
1002 --------
1003 root : Interface to root finding algorithms for multivariate
1004 functions. See ``method='anderson'`` in particular.
1006 References
1007 ----------
1008 .. [Ey] V. Eyert, J. Comp. Phys., 124, 271 (1996).
1010 Examples
1011 --------
1012 The following functions define a system of nonlinear equations
1014 >>> def fun(x):
1015 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
1016 ... 0.5 * (x[1] - x[0])**3 + x[1]]
1018 A solution can be obtained as follows.
1020 >>> from scipy import optimize
1021 >>> sol = optimize.anderson(fun, [0, 0])
1022 >>> sol
1023 array([0.84116588, 0.15883789])
1025 """
1027 # Note:
1028 #
1029 # Anderson method maintains a rank M approximation of the inverse Jacobian,
1030 #
1031 # J^-1 v ~ -v*alpha + (dX + alpha dF) A^-1 dF^H v
1032 # A = W + dF^H dF
1033 # W = w0^2 diag(dF^H dF)
1034 #
1035 # so that for w0 = 0 the secant condition applies for last M iterates, i.e.,
1036 #
1037 # J^-1 df_j = dx_j
1038 #
1039 # for all j = 0 ... M-1.
1040 #
1041 # Moreover, (from Sherman-Morrison-Woodbury formula)
1042 #
1043 # J v ~ [ b I - b^2 C (I + b dF^H A^-1 C)^-1 dF^H ] v
1044 # C = (dX + alpha dF) A^-1
1045 # b = -1/alpha
1046 #
1047 # and after simplification
1048 #
1049 # J v ~ -v/alpha + (dX/alpha + dF) (dF^H dX - alpha W)^-1 dF^H v
1050 #
1052 def __init__(self, alpha=None, w0=0.01, M=5):
1053 GenericBroyden.__init__(self)
1054 self.alpha = alpha
1055 self.M = M
1056 self.dx = []
1057 self.df = []
1058 self.gamma = None
1059 self.w0 = w0
1061 def solve(self, f, tol=0):
1062 dx = -self.alpha*f
1064 n = len(self.dx)
1065 if n == 0:
1066 return dx
1068 df_f = np.empty(n, dtype=f.dtype)
1069 for k in range(n):
1070 df_f[k] = vdot(self.df[k], f)
1072 try:
1073 gamma = solve(self.a, df_f)
1074 except LinAlgError:
1075 # singular; reset the Jacobian approximation
1076 del self.dx[:]
1077 del self.df[:]
1078 return dx
1080 for m in range(n):
1081 dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m])
1082 return dx
1084 def matvec(self, f):
1085 dx = -f/self.alpha
1087 n = len(self.dx)
1088 if n == 0:
1089 return dx
1091 df_f = np.empty(n, dtype=f.dtype)
1092 for k in range(n):
1093 df_f[k] = vdot(self.df[k], f)
1095 b = np.empty((n, n), dtype=f.dtype)
1096 for i in range(n):
1097 for j in range(n):
1098 b[i,j] = vdot(self.df[i], self.dx[j])
1099 if i == j and self.w0 != 0:
1100 b[i,j] -= vdot(self.df[i], self.df[i])*self.w0**2*self.alpha
1101 gamma = solve(b, df_f)
1103 for m in range(n):
1104 dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha)
1105 return dx
1107 def _update(self, x, f, dx, df, dx_norm, df_norm):
1108 if self.M == 0:
1109 return
1111 self.dx.append(dx)
1112 self.df.append(df)
1114 while len(self.dx) > self.M:
1115 self.dx.pop(0)
1116 self.df.pop(0)
1118 n = len(self.dx)
1119 a = np.zeros((n, n), dtype=f.dtype)
1121 for i in range(n):
1122 for j in range(i, n):
1123 if i == j:
1124 wd = self.w0**2
1125 else:
1126 wd = 0
1127 a[i,j] = (1+wd)*vdot(self.df[i], self.df[j])
1129 a += np.triu(a, 1).T.conj()
1130 self.a = a
1132#------------------------------------------------------------------------------
1133# Simple iterations
1134#------------------------------------------------------------------------------
1137class DiagBroyden(GenericBroyden):
1138 """
1139 Find a root of a function, using diagonal Broyden Jacobian approximation.
1141 The Jacobian approximation is derived from previous iterations, by
1142 retaining only the diagonal of Broyden matrices.
1144 .. warning::
1146 This algorithm may be useful for specific problems, but whether
1147 it will work may depend strongly on the problem.
1149 Parameters
1150 ----------
1151 %(params_basic)s
1152 alpha : float, optional
1153 Initial guess for the Jacobian is (-1/alpha).
1154 %(params_extra)s
1156 See Also
1157 --------
1158 root : Interface to root finding algorithms for multivariate
1159 functions. See ``method='diagbroyden'`` in particular.
1161 Examples
1162 --------
1163 The following functions define a system of nonlinear equations
1165 >>> def fun(x):
1166 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0,
1167 ... 0.5 * (x[1] - x[0])**3 + x[1]]
1169 A solution can be obtained as follows.
1171 >>> from scipy import optimize
1172 >>> sol = optimize.diagbroyden(fun, [0, 0])
1173 >>> sol
1174 array([0.84116403, 0.15883384])
1176 """
1178 def __init__(self, alpha=None):
1179 GenericBroyden.__init__(self)
1180 self.alpha = alpha
1182 def setup(self, x, F, func):
1183 GenericBroyden.setup(self, x, F, func)
1184 self.d = np.full((self.shape[0],), 1 / self.alpha, dtype=self.dtype)
1186 def solve(self, f, tol=0):
1187 return -f / self.d
1189 def matvec(self, f):
1190 return -f * self.d
1192 def rsolve(self, f, tol=0):
1193 return -f / self.d.conj()
1195 def rmatvec(self, f):
1196 return -f * self.d.conj()
1198 def todense(self):
1199 return np.diag(-self.d)
1201 def _update(self, x, f, dx, df, dx_norm, df_norm):
1202 self.d -= (df + self.d*dx)*dx/dx_norm**2
1205class LinearMixing(GenericBroyden):
1206 """
1207 Find a root of a function, using a scalar Jacobian approximation.
1209 .. warning::
1211 This algorithm may be useful for specific problems, but whether
1212 it will work may depend strongly on the problem.
1214 Parameters
1215 ----------
1216 %(params_basic)s
1217 alpha : float, optional
1218 The Jacobian approximation is (-1/alpha).
1219 %(params_extra)s
1221 See Also
1222 --------
1223 root : Interface to root finding algorithms for multivariate
1224 functions. See ``method='linearmixing'`` in particular.
1226 """
1228 def __init__(self, alpha=None):
1229 GenericBroyden.__init__(self)
1230 self.alpha = alpha
1232 def solve(self, f, tol=0):
1233 return -f*self.alpha
1235 def matvec(self, f):
1236 return -f/self.alpha
1238 def rsolve(self, f, tol=0):
1239 return -f*np.conj(self.alpha)
1241 def rmatvec(self, f):
1242 return -f/np.conj(self.alpha)
1244 def todense(self):
1245 return np.diag(np.full(self.shape[0], -1/self.alpha))
1247 def _update(self, x, f, dx, df, dx_norm, df_norm):
1248 pass
1251class ExcitingMixing(GenericBroyden):
1252 """
1253 Find a root of a function, using a tuned diagonal Jacobian approximation.
1255 The Jacobian matrix is diagonal and is tuned on each iteration.
1257 .. warning::
1259 This algorithm may be useful for specific problems, but whether
1260 it will work may depend strongly on the problem.
1262 See Also
1263 --------
1264 root : Interface to root finding algorithms for multivariate
1265 functions. See ``method='excitingmixing'`` in particular.
1267 Parameters
1268 ----------
1269 %(params_basic)s
1270 alpha : float, optional
1271 Initial Jacobian approximation is (-1/alpha).
1272 alphamax : float, optional
1273 The entries of the diagonal Jacobian are kept in the range
1274 ``[alpha, alphamax]``.
1275 %(params_extra)s
1276 """
1278 def __init__(self, alpha=None, alphamax=1.0):
1279 GenericBroyden.__init__(self)
1280 self.alpha = alpha
1281 self.alphamax = alphamax
1282 self.beta = None
1284 def setup(self, x, F, func):
1285 GenericBroyden.setup(self, x, F, func)
1286 self.beta = np.full((self.shape[0],), self.alpha, dtype=self.dtype)
1288 def solve(self, f, tol=0):
1289 return -f*self.beta
1291 def matvec(self, f):
1292 return -f/self.beta
1294 def rsolve(self, f, tol=0):
1295 return -f*self.beta.conj()
1297 def rmatvec(self, f):
1298 return -f/self.beta.conj()
1300 def todense(self):
1301 return np.diag(-1/self.beta)
1303 def _update(self, x, f, dx, df, dx_norm, df_norm):
1304 incr = f*self.last_f > 0
1305 self.beta[incr] += self.alpha
1306 self.beta[~incr] = self.alpha
1307 np.clip(self.beta, 0, self.alphamax, out=self.beta)
1310#------------------------------------------------------------------------------
1311# Iterative/Krylov approximated Jacobians
1312#------------------------------------------------------------------------------
1314class KrylovJacobian(Jacobian):
1315 r"""
1316 Find a root of a function, using Krylov approximation for inverse Jacobian.
1318 This method is suitable for solving large-scale problems.
1320 Parameters
1321 ----------
1322 %(params_basic)s
1323 rdiff : float, optional
1324 Relative step size to use in numerical differentiation.
1325 method : str or callable, optional
1326 Krylov method to use to approximate the Jacobian. Can be a string,
1327 or a function implementing the same interface as the iterative
1328 solvers in `scipy.sparse.linalg`. If a string, needs to be one of:
1329 ``'lgmres'``, ``'gmres'``, ``'bicgstab'``, ``'cgs'``, ``'minres'``,
1330 ``'tfqmr'``.
1332 The default is `scipy.sparse.linalg.lgmres`.
1333 inner_maxiter : int, optional
1334 Parameter to pass to the "inner" Krylov solver: maximum number of
1335 iterations. Iteration will stop after maxiter steps even if the
1336 specified tolerance has not been achieved.
1337 inner_M : LinearOperator or InverseJacobian
1338 Preconditioner for the inner Krylov iteration.
1339 Note that you can use also inverse Jacobians as (adaptive)
1340 preconditioners. For example,
1342 >>> from scipy.optimize import BroydenFirst, KrylovJacobian
1343 >>> from scipy.optimize import InverseJacobian
1344 >>> jac = BroydenFirst()
1345 >>> kjac = KrylovJacobian(inner_M=InverseJacobian(jac))
1347 If the preconditioner has a method named 'update', it will be called
1348 as ``update(x, f)`` after each nonlinear step, with ``x`` giving
1349 the current point, and ``f`` the current function value.
1350 outer_k : int, optional
1351 Size of the subspace kept across LGMRES nonlinear iterations.
1352 See `scipy.sparse.linalg.lgmres` for details.
1353 inner_kwargs : kwargs
1354 Keyword parameters for the "inner" Krylov solver
1355 (defined with `method`). Parameter names must start with
1356 the `inner_` prefix which will be stripped before passing on
1357 the inner method. See, e.g., `scipy.sparse.linalg.gmres` for details.
1358 %(params_extra)s
1360 See Also
1361 --------
1362 root : Interface to root finding algorithms for multivariate
1363 functions. See ``method='krylov'`` in particular.
1364 scipy.sparse.linalg.gmres
1365 scipy.sparse.linalg.lgmres
1367 Notes
1368 -----
1369 This function implements a Newton-Krylov solver. The basic idea is
1370 to compute the inverse of the Jacobian with an iterative Krylov
1371 method. These methods require only evaluating the Jacobian-vector
1372 products, which are conveniently approximated by a finite difference:
1374 .. math:: J v \approx (f(x + \omega*v/|v|) - f(x)) / \omega
1376 Due to the use of iterative matrix inverses, these methods can
1377 deal with large nonlinear problems.
1379 SciPy's `scipy.sparse.linalg` module offers a selection of Krylov
1380 solvers to choose from. The default here is `lgmres`, which is a
1381 variant of restarted GMRES iteration that reuses some of the
1382 information obtained in the previous Newton steps to invert
1383 Jacobians in subsequent steps.
1385 For a review on Newton-Krylov methods, see for example [1]_,
1386 and for the LGMRES sparse inverse method, see [2]_.
1388 References
1389 ----------
1390 .. [1] C. T. Kelley, Solving Nonlinear Equations with Newton's Method,
1391 SIAM, pp.57-83, 2003.
1392 :doi:`10.1137/1.9780898718898.ch3`
1393 .. [2] D.A. Knoll and D.E. Keyes, J. Comp. Phys. 193, 357 (2004).
1394 :doi:`10.1016/j.jcp.2003.08.010`
1395 .. [3] A.H. Baker and E.R. Jessup and T. Manteuffel,
1396 SIAM J. Matrix Anal. Appl. 26, 962 (2005).
1397 :doi:`10.1137/S0895479803422014`
1399 Examples
1400 --------
1401 The following functions define a system of nonlinear equations
1403 >>> def fun(x):
1404 ... return [x[0] + 0.5 * x[1] - 1.0,
1405 ... 0.5 * (x[1] - x[0]) ** 2]
1407 A solution can be obtained as follows.
1409 >>> from scipy import optimize
1410 >>> sol = optimize.newton_krylov(fun, [0, 0])
1411 >>> sol
1412 array([0.66731771, 0.66536458])
1414 """
1416 def __init__(self, rdiff=None, method='lgmres', inner_maxiter=20,
1417 inner_M=None, outer_k=10, **kw):
1418 self.preconditioner = inner_M
1419 self.rdiff = rdiff
1420 # Note that this retrieves one of the named functions, or otherwise
1421 # uses `method` as is (i.e., for a user-provided callable).
1422 self.method = dict(
1423 bicgstab=scipy.sparse.linalg.bicgstab,
1424 gmres=scipy.sparse.linalg.gmres,
1425 lgmres=scipy.sparse.linalg.lgmres,
1426 cgs=scipy.sparse.linalg.cgs,
1427 minres=scipy.sparse.linalg.minres,
1428 tfqmr=scipy.sparse.linalg.tfqmr,
1429 ).get(method, method)
1431 self.method_kw = dict(maxiter=inner_maxiter, M=self.preconditioner)
1433 if self.method is scipy.sparse.linalg.gmres:
1434 # Replace GMRES's outer iteration with Newton steps
1435 self.method_kw['restart'] = inner_maxiter
1436 self.method_kw['maxiter'] = 1
1437 self.method_kw.setdefault('atol', 0)
1438 elif self.method in (scipy.sparse.linalg.gcrotmk,
1439 scipy.sparse.linalg.bicgstab,
1440 scipy.sparse.linalg.cgs):
1441 self.method_kw.setdefault('atol', 0)
1442 elif self.method is scipy.sparse.linalg.lgmres:
1443 self.method_kw['outer_k'] = outer_k
1444 # Replace LGMRES's outer iteration with Newton steps
1445 self.method_kw['maxiter'] = 1
1446 # Carry LGMRES's `outer_v` vectors across nonlinear iterations
1447 self.method_kw.setdefault('outer_v', [])
1448 self.method_kw.setdefault('prepend_outer_v', True)
1449 # But don't carry the corresponding Jacobian*v products, in case
1450 # the Jacobian changes a lot in the nonlinear step
1451 #
1452 # XXX: some trust-region inspired ideas might be more efficient...
1453 # See e.g., Brown & Saad. But needs to be implemented separately
1454 # since it's not an inexact Newton method.
1455 self.method_kw.setdefault('store_outer_Av', False)
1456 self.method_kw.setdefault('atol', 0)
1458 for key, value in kw.items():
1459 if not key.startswith('inner_'):
1460 raise ValueError("Unknown parameter %s" % key)
1461 self.method_kw[key[6:]] = value
1463 def _update_diff_step(self):
1464 mx = abs(self.x0).max()
1465 mf = abs(self.f0).max()
1466 self.omega = self.rdiff * max(1, mx) / max(1, mf)
1468 def matvec(self, v):
1469 nv = norm(v)
1470 if nv == 0:
1471 return 0*v
1472 sc = self.omega / nv
1473 r = (self.func(self.x0 + sc*v) - self.f0) / sc
1474 if not np.all(np.isfinite(r)) and np.all(np.isfinite(v)):
1475 raise ValueError('Function returned non-finite results')
1476 return r
1478 def solve(self, rhs, tol=0):
1479 if 'tol' in self.method_kw:
1480 sol, info = self.method(self.op, rhs, **self.method_kw)
1481 else:
1482 sol, info = self.method(self.op, rhs, tol=tol, **self.method_kw)
1483 return sol
1485 def update(self, x, f):
1486 self.x0 = x
1487 self.f0 = f
1488 self._update_diff_step()
1490 # Update also the preconditioner, if possible
1491 if self.preconditioner is not None:
1492 if hasattr(self.preconditioner, 'update'):
1493 self.preconditioner.update(x, f)
1495 def setup(self, x, f, func):
1496 Jacobian.setup(self, x, f, func)
1497 self.x0 = x
1498 self.f0 = f
1499 self.op = scipy.sparse.linalg.aslinearoperator(self)
1501 if self.rdiff is None:
1502 self.rdiff = np.finfo(x.dtype).eps ** (1./2)
1504 self._update_diff_step()
1506 # Setup also the preconditioner, if possible
1507 if self.preconditioner is not None:
1508 if hasattr(self.preconditioner, 'setup'):
1509 self.preconditioner.setup(x, f, func)
1512#------------------------------------------------------------------------------
1513# Wrapper functions
1514#------------------------------------------------------------------------------
1516def _nonlin_wrapper(name, jac):
1517 """
1518 Construct a solver wrapper with given name and Jacobian approx.
1520 It inspects the keyword arguments of ``jac.__init__``, and allows to
1521 use the same arguments in the wrapper function, in addition to the
1522 keyword arguments of `nonlin_solve`
1524 """
1525 signature = _getfullargspec(jac.__init__)
1526 args, varargs, varkw, defaults, kwonlyargs, kwdefaults, _ = signature
1527 kwargs = list(zip(args[-len(defaults):], defaults))
1528 kw_str = ", ".join(["%s=%r" % (k, v) for k, v in kwargs])
1529 if kw_str:
1530 kw_str = ", " + kw_str
1531 kwkw_str = ", ".join(["%s=%s" % (k, k) for k, v in kwargs])
1532 if kwkw_str:
1533 kwkw_str = kwkw_str + ", "
1534 if kwonlyargs:
1535 raise ValueError('Unexpected signature %s' % signature)
1537 # Construct the wrapper function so that its keyword arguments
1538 # are visible in pydoc.help etc.
1539 wrapper = """
1540def %(name)s(F, xin, iter=None %(kw)s, verbose=False, maxiter=None,
1541 f_tol=None, f_rtol=None, x_tol=None, x_rtol=None,
1542 tol_norm=None, line_search='armijo', callback=None, **kw):
1543 jac = %(jac)s(%(kwkw)s **kw)
1544 return nonlin_solve(F, xin, jac, iter, verbose, maxiter,
1545 f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search,
1546 callback)
1547"""
1549 wrapper = wrapper % dict(name=name, kw=kw_str, jac=jac.__name__,
1550 kwkw=kwkw_str)
1551 ns = {}
1552 ns.update(globals())
1553 exec(wrapper, ns)
1554 func = ns[name]
1555 func.__doc__ = jac.__doc__
1556 _set_doc(func)
1557 return func
1560broyden1 = _nonlin_wrapper('broyden1', BroydenFirst)
1561broyden2 = _nonlin_wrapper('broyden2', BroydenSecond)
1562anderson = _nonlin_wrapper('anderson', Anderson)
1563linearmixing = _nonlin_wrapper('linearmixing', LinearMixing)
1564diagbroyden = _nonlin_wrapper('diagbroyden', DiagBroyden)
1565excitingmixing = _nonlin_wrapper('excitingmixing', ExcitingMixing)
1566newton_krylov = _nonlin_wrapper('newton_krylov', KrylovJacobian)