Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/optimize/nonlin.py : 21%

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"""
3Nonlinear solvers
4-----------------
6.. currentmodule:: scipy.optimize
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.
12Routines
13~~~~~~~~
15Large-scale nonlinear solvers:
17.. autosummary::
19 newton_krylov
20 anderson
22General nonlinear solvers:
24.. autosummary::
26 broyden1
27 broyden2
29Simple iterations:
31.. autosummary::
33 excitingmixing
34 linearmixing
35 diagbroyden
38Examples
39~~~~~~~~
41**Small problem**
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.])
53**Large problem**
55Suppose that we needed to solve the following integrodifferential
56equation on the square :math:`[0,1]\times[0,1]`:
58.. math::
60 \nabla^2 P = 10 \left(\int_0^1\int_0^1\cosh(P)\,dx\,dy\right)^2
62with :math:`P(x,1) = 1` and :math:`P=0` elsewhere on the boundary of
63the square.
65The solution can be found using the `newton_krylov` solver:
67.. plot::
69 import numpy as np
70 from scipy.optimize import newton_krylov
71 from numpy import cosh, zeros_like, mgrid, zeros
73 # parameters
74 nx, ny = 75, 75
75 hx, hy = 1./(nx-1), 1./(ny-1)
77 P_left, P_right = 0, 0
78 P_top, P_bottom = 1, 0
80 def residual(P):
81 d2x = zeros_like(P)
82 d2y = zeros_like(P)
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
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
92 return d2x + d2y - 10*cosh(P).mean()**2
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())
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()
106"""
107# Copyright (C) 2009, Pauli Virtanen <pav@iki.fi>
108# Distributed under the same license as SciPy.
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
122__all__ = [
123 'broyden1', 'broyden2', 'anderson', 'linearmixing',
124 'diagbroyden', 'excitingmixing', 'newton_krylov']
126#------------------------------------------------------------------------------
127# Utility functions
128#------------------------------------------------------------------------------
131class NoConvergence(Exception):
132 pass
135def maxnorm(x):
136 return np.absolute(x).max()
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
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)
154def _safe_norm(v):
155 if not np.isfinite(v).all():
156 return np.array(np.inf)
157 return norm(v)
159#------------------------------------------------------------------------------
160# Generic nonlinear solver machinery
161#------------------------------------------------------------------------------
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.
202 Returns
203 -------
204 sol : ndarray
205 An array (of similar array type as `x0`) containing the final solution.
207 Raises
208 ------
209 NoConvergence
210 When a solution was not found.
212 """.strip()
213)
216def _set_doc(obj):
217 if obj.__doc__:
218 obj.__doc__ = obj.__doc__ % _doc_parts
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.
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:
236 krylov, broyden1, broyden2, anderson
237 diagbroyden, linearmixing, excitingmixing
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.
246 See Also
247 --------
248 asjacobian, Jacobian
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.
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/
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)
271 x0 = _as_inexact(x0)
272 func = lambda z: _as_inexact(F(_array_like(z, x0))).flatten()
273 x = x0.flatten()
275 dx = np.inf
276 Fx = func(x)
277 Fx_norm = norm(Fx)
279 jacobian = asjacobian(jacobian)
280 jacobian.setup(x.copy(), Fx, func)
282 if maxiter is None:
283 if iter is not None:
284 maxiter = iter + 1
285 else:
286 maxiter = 100*(x.size+1)
288 if line_search is True:
289 line_search = 'armijo'
290 elif line_search is False:
291 line_search = None
293 if line_search not in (None, 'armijo', 'wolfe'):
294 raise ValueError("Invalid line search")
296 # Solver tolerance selection
297 gamma = 0.9
298 eta_max = 0.9999
299 eta_treshold = 0.1
300 eta = 1e-3
302 for n in range(maxiter):
303 status = condition.check(Fx, x, dx)
304 if status:
305 break
307 # The tolerance, as computed for scipy.sparse.linalg.* routines
308 tol = min(eta, eta*Fx_norm)
309 dx = -jacobian.solve(Fx, tol=tol)
311 if norm(dx) == 0:
312 raise ValueError("Jacobian inversion yielded zero vector. "
313 "This indicates a bug in the Jacobian "
314 "approximation.")
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)
326 jacobian.update(x.copy(), Fx)
328 if callback:
329 callback(x, Fx)
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))
338 Fx_norm = Fx_norm_new
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
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)
367_set_doc(nonlin_solve)
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)
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
389 def derphi(s):
390 ds = (abs(s) + s_norm + 1) * rdiff
391 return (phi(s+ds, store=False) - phi(s)) / ds
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)
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
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)
412 return s, x, Fx, Fx_norm
415class TerminationCondition(object):
416 """
417 Termination condition for an iteration. It is terminated if
419 - |F| < f_rtol*|F_0|, AND
420 - |F| < f_tol
422 AND
424 - |dx| < x_rtol*|x|, AND
425 - |dx| < x_tol
427 """
428 def __init__(self, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None,
429 iter=None, norm=maxnorm):
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
440 self.x_tol = x_tol
441 self.x_rtol = x_rtol
442 self.f_tol = f_tol
443 self.f_rtol = f_rtol
445 self.norm = norm
447 self.iter = iter
449 self.f0_norm = None
450 self.iteration = 0
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)
458 if self.f0_norm is None:
459 self.f0_norm = f_norm
461 if f_norm == 0:
462 return 1
464 if self.iter is not None:
465 # backwards compatibility with SciPy 0.6.0
466 return 2 * (self.iteration > self.iter)
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))
475#------------------------------------------------------------------------------
476# Generic Jacobian approximation
477#------------------------------------------------------------------------------
479class Jacobian(object):
480 """
481 Common interface for Jacobians or Jacobian approximations.
483 The optional methods come useful when implementing trust region
484 etc., algorithms that often require evaluating transposes of the
485 Jacobian.
487 Methods
488 -------
489 solve
490 Returns J^-1 * v
491 update
492 Updates Jacobian to point `x` (where the function has residual `Fx`)
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.
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
515 """
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])
526 if hasattr(self, 'todense'):
527 self.__array__ = lambda: self.todense()
529 def aspreconditioner(self):
530 return InverseJacobian(self)
532 def solve(self, v, tol=0):
533 raise NotImplementedError
535 def update(self, x, F):
536 pass
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)
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
557 @property
558 def shape(self):
559 return self.jacobian.shape
561 @property
562 def dtype(self):
563 return self.jacobian.dtype
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')
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
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")
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")
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")
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')
658#------------------------------------------------------------------------------
659# Broyden
660#------------------------------------------------------------------------------
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
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
677 def _update(self, x, f, dx, df, dx_norm, df_norm):
678 raise NotImplementedError
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
688class LowRankMatrix(object):
689 r"""
690 A matrix represented as
692 .. math:: \alpha I + \sum_{n=0}^{n=M} c_n d_n^\dagger
694 However, if the rank of the matrix reaches the dimension of the vectors,
695 full matrix representation will be used thereon.
697 """
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
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
717 @staticmethod
718 def _solve(v, alpha, cs, ds):
719 """Evaluate w = M^-1 v"""
720 if len(cs) == 0:
721 return v/alpha
723 # (B + C D^H)^-1 = B^-1 - B^-1 C (I + D^H B^-1 C)^-1 D^H B^-1
725 axpy, dotc = get_blas_funcs(['axpy', 'dotc'], cs[:1] + [v])
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)
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)
739 w = v/alpha
740 for c, qc in zip(cs, q):
741 w = axpy(c, w, w.size, -qc)
743 return w
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)
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)
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)
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)
769 def append(self, c, d):
770 if self.collapsed is not None:
771 self.collapsed += c[:,None] * d[None,:].conj()
772 return
774 self.cs.append(c)
775 self.ds.append(d)
777 if len(self.cs) > c.size:
778 self.collapse()
780 def __array__(self):
781 if self.collapsed is not None:
782 return self.collapsed
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
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
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[:]
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]
818 def svd_reduce(self, max_rank, to_retain=None):
819 """
820 Reduce the rank of the matrix by retaining some SVD components.
822 This corresponds to the \"Broyden Rank Reduction Inverse\"
823 algorithm described in [1]_.
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.
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``.
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).
844 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
846 """
847 if self.collapsed is not None:
848 return
850 p = max_rank
851 if to_retain is not None:
852 q = to_retain
853 else:
854 q = p - 2
856 if self.cs:
857 p = min(p, len(self.cs[0]))
858 q = max(0, min(q, p-1))
860 m = len(self.cs)
861 if m < p:
862 # nothing to do
863 return
865 C = np.array(self.cs).T
866 D = np.array(self.ds).T
868 D, R = qr(D, mode='economic')
869 C = dot(C, R.T.conj())
871 U, S, WH = svd(C, full_matrices=False, compute_uv=True)
873 C = dot(C, inv(WH))
874 D = dot(D, WH.T.conj())
876 for k in range(q):
877 self.cs[k] = C[:,k].copy()
878 self.ds[k] = D[:,k].copy()
880 del self.cs[q:]
881 del self.ds[q:]
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.
893 Methods available:
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``.
902 max_rank : int, optional
903 Maximum rank for the Broyden matrix.
904 Default is infinity (i.e., no rank reduction).
905 """.strip()
908class BroydenFirst(GenericBroyden):
909 r"""
910 Find a root of a function, using Broyden's first Jacobian approximation.
912 This method is also known as \"Broyden's good method\".
914 Parameters
915 ----------
916 %(params_basic)s
917 %(broyden_params)s
918 %(params_extra)s
920 See Also
921 --------
922 root : Interface to root finding algorithms for multivariate
923 functions. See ``method=='broyden1'`` in particular.
925 Notes
926 -----
927 This algorithm implements the inverse Jacobian Quasi-Newton update
929 .. math:: H_+ = H + (dx - H df) dx^\dagger H / ( dx^\dagger H df)
931 which corresponds to Broyden's first Jacobian update
933 .. math:: J_+ = J + (df - J dx) dx^\dagger / dx^\dagger dx
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).
943 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
945 """
947 def __init__(self, alpha=None, reduction_method='restart', max_rank=None):
948 GenericBroyden.__init__(self)
949 self.alpha = alpha
950 self.Gm = None
952 if max_rank is None:
953 max_rank = np.inf
954 self.max_rank = max_rank
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
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)
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)
977 def todense(self):
978 return inv(self.Gm)
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)
987 def matvec(self, f):
988 return self.Gm.solve(f)
990 def rsolve(self, f, tol=0):
991 return self.Gm.rmatvec(f)
993 def rmatvec(self, f):
994 return self.Gm.rsolve(f)
996 def _update(self, x, f, dx, df, dx_norm, df_norm):
997 self._reduce() # reduce first to preserve secant condition
999 v = self.Gm.rmatvec(dx)
1000 c = dx - self.Gm.matvec(df)
1001 d = v / vdot(df, v)
1003 self.Gm.append(c, d)
1006class BroydenSecond(BroydenFirst):
1007 """
1008 Find a root of a function, using Broyden\'s second Jacobian approximation.
1010 This method is also known as \"Broyden's bad method\".
1012 Parameters
1013 ----------
1014 %(params_basic)s
1015 %(broyden_params)s
1016 %(params_extra)s
1018 See Also
1019 --------
1020 root : Interface to root finding algorithms for multivariate
1021 functions. See ``method=='broyden2'`` in particular.
1023 Notes
1024 -----
1025 This algorithm implements the inverse Jacobian Quasi-Newton update
1027 .. math:: H_+ = H + (dx - H df) df^\\dagger / ( df^\\dagger df)
1029 corresponding to Broyden's second method.
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).
1038 https://web.archive.org/web/20161022015821/http://www.math.leidenuniv.nl/scripties/Rotten.pdf
1040 """
1042 def _update(self, x, f, dx, df, dx_norm, df_norm):
1043 self._reduce() # reduce first to preserve secant condition
1045 v = df
1046 c = dx - self.Gm.matvec(df)
1047 d = v / df_norm**2
1048 self.Gm.append(c, d)
1051#------------------------------------------------------------------------------
1052# Broyden-like (restricted memory)
1053#------------------------------------------------------------------------------
1055class Anderson(GenericBroyden):
1056 """
1057 Find a root of a function, using (extended) Anderson mixing.
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]_
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
1075 See Also
1076 --------
1077 root : Interface to root finding algorithms for multivariate
1078 functions. See ``method=='anderson'`` in particular.
1080 References
1081 ----------
1082 .. [Ey] V. Eyert, J. Comp. Phys., 124, 271 (1996).
1084 """
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 #
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
1120 def solve(self, f, tol=0):
1121 dx = -self.alpha*f
1123 n = len(self.dx)
1124 if n == 0:
1125 return dx
1127 df_f = np.empty(n, dtype=f.dtype)
1128 for k in range(n):
1129 df_f[k] = vdot(self.df[k], f)
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
1139 for m in range(n):
1140 dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m])
1141 return dx
1143 def matvec(self, f):
1144 dx = -f/self.alpha
1146 n = len(self.dx)
1147 if n == 0:
1148 return dx
1150 df_f = np.empty(n, dtype=f.dtype)
1151 for k in range(n):
1152 df_f[k] = vdot(self.df[k], f)
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)
1162 for m in range(n):
1163 dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha)
1164 return dx
1166 def _update(self, x, f, dx, df, dx_norm, df_norm):
1167 if self.M == 0:
1168 return
1170 self.dx.append(dx)
1171 self.df.append(df)
1173 while len(self.dx) > self.M:
1174 self.dx.pop(0)
1175 self.df.pop(0)
1177 n = len(self.dx)
1178 a = np.zeros((n, n), dtype=f.dtype)
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])
1188 a += np.triu(a, 1).T.conj()
1189 self.a = a
1191#------------------------------------------------------------------------------
1192# Simple iterations
1193#------------------------------------------------------------------------------
1196class DiagBroyden(GenericBroyden):
1197 """
1198 Find a root of a function, using diagonal Broyden Jacobian approximation.
1200 The Jacobian approximation is derived from previous iterations, by
1201 retaining only the diagonal of Broyden matrices.
1203 .. warning::
1205 This algorithm may be useful for specific problems, but whether
1206 it will work may depend strongly on the problem.
1208 Parameters
1209 ----------
1210 %(params_basic)s
1211 alpha : float, optional
1212 Initial guess for the Jacobian is (-1/alpha).
1213 %(params_extra)s
1215 See Also
1216 --------
1217 root : Interface to root finding algorithms for multivariate
1218 functions. See ``method=='diagbroyden'`` in particular.
1219 """
1221 def __init__(self, alpha=None):
1222 GenericBroyden.__init__(self)
1223 self.alpha = alpha
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)
1229 def solve(self, f, tol=0):
1230 return -f / self.d
1232 def matvec(self, f):
1233 return -f * self.d
1235 def rsolve(self, f, tol=0):
1236 return -f / self.d.conj()
1238 def rmatvec(self, f):
1239 return -f * self.d.conj()
1241 def todense(self):
1242 return np.diag(-self.d)
1244 def _update(self, x, f, dx, df, dx_norm, df_norm):
1245 self.d -= (df + self.d*dx)*dx/dx_norm**2
1248class LinearMixing(GenericBroyden):
1249 """
1250 Find a root of a function, using a scalar Jacobian approximation.
1252 .. warning::
1254 This algorithm may be useful for specific problems, but whether
1255 it will work may depend strongly on the problem.
1257 Parameters
1258 ----------
1259 %(params_basic)s
1260 alpha : float, optional
1261 The Jacobian approximation is (-1/alpha).
1262 %(params_extra)s
1264 See Also
1265 --------
1266 root : Interface to root finding algorithms for multivariate
1267 functions. See ``method=='linearmixing'`` in particular.
1269 """
1271 def __init__(self, alpha=None):
1272 GenericBroyden.__init__(self)
1273 self.alpha = alpha
1275 def solve(self, f, tol=0):
1276 return -f*self.alpha
1278 def matvec(self, f):
1279 return -f/self.alpha
1281 def rsolve(self, f, tol=0):
1282 return -f*np.conj(self.alpha)
1284 def rmatvec(self, f):
1285 return -f/np.conj(self.alpha)
1287 def todense(self):
1288 return np.diag(np.full(self.shape[0], -1/self.alpha))
1290 def _update(self, x, f, dx, df, dx_norm, df_norm):
1291 pass
1294class ExcitingMixing(GenericBroyden):
1295 """
1296 Find a root of a function, using a tuned diagonal Jacobian approximation.
1298 The Jacobian matrix is diagonal and is tuned on each iteration.
1300 .. warning::
1302 This algorithm may be useful for specific problems, but whether
1303 it will work may depend strongly on the problem.
1305 See Also
1306 --------
1307 root : Interface to root finding algorithms for multivariate
1308 functions. See ``method=='excitingmixing'`` in particular.
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 """
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
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)
1331 def solve(self, f, tol=0):
1332 return -f*self.beta
1334 def matvec(self, f):
1335 return -f/self.beta
1337 def rsolve(self, f, tol=0):
1338 return -f*self.beta.conj()
1340 def rmatvec(self, f):
1341 return -f/self.beta.conj()
1343 def todense(self):
1344 return np.diag(-1/self.beta)
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)
1353#------------------------------------------------------------------------------
1354# Iterative/Krylov approximated Jacobians
1355#------------------------------------------------------------------------------
1357class KrylovJacobian(Jacobian):
1358 r"""
1359 Find a root of a function, using Krylov approximation for inverse Jacobian.
1361 This method is suitable for solving large-scale problems.
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`.
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,
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))
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
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
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:
1415 .. math:: J v \approx (f(x + \omega*v/|v|) - f(x)) / \omega
1417 Due to the use of iterative matrix inverses, these methods can
1418 deal with large nonlinear problems.
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.
1426 For a review on Newton-Krylov methods, see for example [1]_,
1427 and for the LGMRES sparse inverse method, see [2]_.
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`
1437 """
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)
1451 self.method_kw = dict(maxiter=inner_maxiter, M=self.preconditioner)
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)
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
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)
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
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
1503 def update(self, x, f):
1504 self.x0 = x
1505 self.f0 = f
1506 self._update_diff_step()
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)
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)
1519 if self.rdiff is None:
1520 self.rdiff = np.finfo(x.dtype).eps ** (1./2)
1522 self._update_diff_step()
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)
1530#------------------------------------------------------------------------------
1531# Wrapper functions
1532#------------------------------------------------------------------------------
1534def _nonlin_wrapper(name, jac):
1535 """
1536 Construct a solver wrapper with given name and Jacobian approx.
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`
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)
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"""
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
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)