1""" 2Unified interfaces to root finding algorithms. 3 4Functions 5--------- 6- root : find a root of a vector function. 7""" 8__all__ = ['root'] 9 10import numpy as np 11 12ROOT_METHODS = ['hybr', 'lm', 'broyden1', 'broyden2', 'anderson', 13 'linearmixing', 'diagbroyden', 'excitingmixing', 'krylov', 14 'df-sane'] 15 16from warnings import warn 17 18from .optimize import MemoizeJac, OptimizeResult, _check_unknown_options 19from .minpack import _root_hybr, leastsq 20from ._spectral import _root_df_sane 21from . import nonlin 22 23 24def root(fun, x0, args=(), method='hybr', jac=None, tol=None, callback=None, 25 options=None): 26 """ 27 Find a root of a vector function. 28 29 Parameters 30 ---------- 31 fun : callable 32 A vector function to find a root of. 33 x0 : ndarray 34 Initial guess. 35 args : tuple, optional 36 Extra arguments passed to the objective function and its Jacobian. 37 method : str, optional 38 Type of solver. Should be one of 39 40 - 'hybr' :ref:`(see here) <optimize.root-hybr>` 41 - 'lm' :ref:`(see here) <optimize.root-lm>` 42 - 'broyden1' :ref:`(see here) <optimize.root-broyden1>` 43 - 'broyden2' :ref:`(see here) <optimize.root-broyden2>` 44 - 'anderson' :ref:`(see here) <optimize.root-anderson>` 45 - 'linearmixing' :ref:`(see here) <optimize.root-linearmixing>` 46 - 'diagbroyden' :ref:`(see here) <optimize.root-diagbroyden>` 47 - 'excitingmixing' :ref:`(see here) <optimize.root-excitingmixing>` 48 - 'krylov' :ref:`(see here) <optimize.root-krylov>` 49 - 'df-sane' :ref:`(see here) <optimize.root-dfsane>` 50 51 jac : bool or callable, optional 52 If `jac` is a Boolean and is True, `fun` is assumed to return the 53 value of Jacobian along with the objective function. If False, the 54 Jacobian will be estimated numerically. 55 `jac` can also be a callable returning the Jacobian of `fun`. In 56 this case, it must accept the same arguments as `fun`. 57 tol : float, optional 58 Tolerance for termination. For detailed control, use solver-specific 59 options. 60 callback : function, optional 61 Optional callback function. It is called on every iteration as 62 ``callback(x, f)`` where `x` is the current solution and `f` 63 the corresponding residual. For all methods but 'hybr' and 'lm'. 64 options : dict, optional 65 A dictionary of solver options. E.g., `xtol` or `maxiter`, see 66 :obj:`show_options()` for details. 67 68 Returns 69 ------- 70 sol : OptimizeResult 71 The solution represented as a ``OptimizeResult`` object. 72 Important attributes are: ``x`` the solution array, ``success`` a 73 Boolean flag indicating if the algorithm exited successfully and 74 ``message`` which describes the cause of the termination. See 75 `OptimizeResult` for a description of other attributes. 76 77 See also 78 -------- 79 show_options : Additional options accepted by the solvers 80 81 Notes 82 ----- 83 This section describes the available solvers that can be selected by the 84 'method' parameter. The default method is *hybr*. 85 86 Method *hybr* uses a modification of the Powell hybrid method as 87 implemented in MINPACK [1]_. 88 89 Method *lm* solves the system of nonlinear equations in a least squares 90 sense using a modification of the Levenberg-Marquardt algorithm as 91 implemented in MINPACK [1]_. 92 93 Method *df-sane* is a derivative-free spectral method. [3]_ 94 95 Methods *broyden1*, *broyden2*, *anderson*, *linearmixing*, 96 *diagbroyden*, *excitingmixing*, *krylov* are inexact Newton methods, 97 with backtracking or full line searches [2]_. Each method corresponds 98 to a particular Jacobian approximations. See `nonlin` for details. 99 100 - Method *broyden1* uses Broyden's first Jacobian approximation, it is 101 known as Broyden's good method. 102 - Method *broyden2* uses Broyden's second Jacobian approximation, it 103 is known as Broyden's bad method. 104 - Method *anderson* uses (extended) Anderson mixing. 105 - Method *Krylov* uses Krylov approximation for inverse Jacobian. It 106 is suitable for large-scale problem. 107 - Method *diagbroyden* uses diagonal Broyden Jacobian approximation. 108 - Method *linearmixing* uses a scalar Jacobian approximation. 109 - Method *excitingmixing* uses a tuned diagonal Jacobian 110 approximation. 111 112 .. warning:: 113 114 The algorithms implemented for methods *diagbroyden*, 115 *linearmixing* and *excitingmixing* may be useful for specific 116 problems, but whether they will work may depend strongly on the 117 problem. 118 119 .. versionadded:: 0.11.0 120 121 References 122 ---------- 123 .. [1] More, Jorge J., Burton S. Garbow, and Kenneth E. Hillstrom. 124 1980. User Guide for MINPACK-1. 125 .. [2] C. T. Kelley. 1995. Iterative Methods for Linear and Nonlinear 126 Equations. Society for Industrial and Applied Mathematics. 127 <https://archive.siam.org/books/kelley/fr16/> 128 .. [3] W. La Cruz, J.M. Martinez, M. Raydan. Math. Comp. 75, 1429 (2006). 129 130 Examples 131 -------- 132 The following functions define a system of nonlinear equations and its 133 jacobian. 134 135 >>> def fun(x): 136 ... return [x[0] + 0.5 * (x[0] - x[1])**3 - 1.0, 137 ... 0.5 * (x[1] - x[0])**3 + x[1]] 138 139 >>> def jac(x): 140 ... return np.array([[1 + 1.5 * (x[0] - x[1])**2, 141 ... -1.5 * (x[0] - x[1])**2], 142 ... [-1.5 * (x[1] - x[0])**2, 143 ... 1 + 1.5 * (x[1] - x[0])**2]]) 144 145 A solution can be obtained as follows. 146 147 >>> from scipy import optimize 148 >>> sol = optimize.root(fun, [0, 0], jac=jac, method='hybr') 149 >>> sol.x 150 array([ 0.8411639, 0.1588361]) 151 152 """ 153 if not isinstance(args, tuple): 154 args = (args,) 155 156 meth = method.lower() 157 if options is None: 158 options = {} 159 160 if callback is not None and meth in ('hybr', 'lm'): 161 warn('Method %s does not accept callback.' % method, 162 RuntimeWarning) 163 164 # fun also returns the Jacobian 165 if not callable(jac) and meth in ('hybr', 'lm'): 166 if bool(jac): 167 fun = MemoizeJac(fun) 168 jac = fun.derivative 169 else: 170 jac = None 171 172 # set default tolerances 173 if tol is not None: 174 options = dict(options) 175 if meth in ('hybr', 'lm'): 176 options.setdefault('xtol', tol) 177 elif meth in ('df-sane',): 178 options.setdefault('ftol', tol) 179 elif meth in ('broyden1', 'broyden2', 'anderson', 'linearmixing', 180 'diagbroyden', 'excitingmixing', 'krylov'): 181 options.setdefault('xtol', tol) 182 options.setdefault('xatol', np.inf) 183 options.setdefault('ftol', np.inf) 184 options.setdefault('fatol', np.inf) 185 186 if meth == 'hybr': 187 sol = _root_hybr(fun, x0, args=args, jac=jac, **options) 188 elif meth == 'lm': 189 sol = _root_leastsq(fun, x0, args=args, jac=jac, **options) 190 elif meth == 'df-sane': 191 _warn_jac_unused(jac, method) 192 sol = _root_df_sane(fun, x0, args=args, callback=callback, 193 **options) 194 elif meth in ('broyden1', 'broyden2', 'anderson', 'linearmixing', 195 'diagbroyden', 'excitingmixing', 'krylov'): 196 _warn_jac_unused(jac, method) 197 sol = _root_nonlin_solve(fun, x0, args=args, jac=jac, 198 _method=meth, _callback=callback, 199 **options) 200 else: 201 raise ValueError('Unknown solver %s' % method) 202 203 return sol 204 205 206def _warn_jac_unused(jac, method): 207 if jac is not None: 208 warn('Method %s does not use the jacobian (jac).' % (method,), 209 RuntimeWarning) 210 211 212def _root_leastsq(fun, x0, args=(), jac=None, 213 col_deriv=0, xtol=1.49012e-08, ftol=1.49012e-08, 214 gtol=0.0, maxiter=0, eps=0.0, factor=100, diag=None, 215 **unknown_options): 216 """ 217 Solve for least squares with Levenberg-Marquardt 218 219 Options 220 ------- 221 col_deriv : bool 222 non-zero to specify that the Jacobian function computes derivatives 223 down the columns (faster, because there is no transpose operation). 224 ftol : float 225 Relative error desired in the sum of squares. 226 xtol : float 227 Relative error desired in the approximate solution. 228 gtol : float 229 Orthogonality desired between the function vector and the columns 230 of the Jacobian. 231 maxiter : int 232 The maximum number of calls to the function. If zero, then 233 100*(N+1) is the maximum where N is the number of elements in x0. 234 epsfcn : float 235 A suitable step length for the forward-difference approximation of 236 the Jacobian (for Dfun=None). If epsfcn is less than the machine 237 precision, it is assumed that the relative errors in the functions 238 are of the order of the machine precision. 239 factor : float 240 A parameter determining the initial step bound 241 (``factor * || diag * x||``). Should be in interval ``(0.1, 100)``. 242 diag : sequence 243 N positive entries that serve as a scale factors for the variables. 244 """ 245 246 _check_unknown_options(unknown_options) 247 x, cov_x, info, msg, ier = leastsq(fun, x0, args=args, Dfun=jac, 248 full_output=True, 249 col_deriv=col_deriv, xtol=xtol, 250 ftol=ftol, gtol=gtol, 251 maxfev=maxiter, epsfcn=eps, 252 factor=factor, diag=diag) 253 sol = OptimizeResult(x=x, message=msg, status=ier, 254 success=ier in (1, 2, 3, 4), cov_x=cov_x, 255 fun=info.pop('fvec')) 256 sol.update(info) 257 return sol 258 259 260def _root_nonlin_solve(fun, x0, args=(), jac=None, 261 _callback=None, _method=None, 262 nit=None, disp=False, maxiter=None, 263 ftol=None, fatol=None, xtol=None, xatol=None, 264 tol_norm=None, line_search='armijo', jac_options=None, 265 **unknown_options): 266 _check_unknown_options(unknown_options) 267 268 f_tol = fatol 269 f_rtol = ftol 270 x_tol = xatol 271 x_rtol = xtol 272 verbose = disp 273 if jac_options is None: 274 jac_options = dict() 275 276 jacobian = {'broyden1': nonlin.BroydenFirst, 277 'broyden2': nonlin.BroydenSecond, 278 'anderson': nonlin.Anderson, 279 'linearmixing': nonlin.LinearMixing, 280 'diagbroyden': nonlin.DiagBroyden, 281 'excitingmixing': nonlin.ExcitingMixing, 282 'krylov': nonlin.KrylovJacobian 283 }[_method] 284 285 if args: 286 if jac: 287 def f(x): 288 return fun(x, *args)[0] 289 else: 290 def f(x): 291 return fun(x, *args) 292 else: 293 f = fun 294 295 x, info = nonlin.nonlin_solve(f, x0, jacobian=jacobian(**jac_options), 296 iter=nit, verbose=verbose, 297 maxiter=maxiter, f_tol=f_tol, 298 f_rtol=f_rtol, x_tol=x_tol, 299 x_rtol=x_rtol, tol_norm=tol_norm, 300 line_search=line_search, 301 callback=_callback, full_output=True, 302 raise_exception=False) 303 sol = OptimizeResult(x=x) 304 sol.update(info) 305 return sol 306 307def _root_broyden1_doc(): 308 """ 309 Options 310 ------- 311 nit : int, optional 312 Number of iterations to make. If omitted (default), make as many 313 as required to meet tolerances. 314 disp : bool, optional 315 Print status to stdout on every iteration. 316 maxiter : int, optional 317 Maximum number of iterations to make. If more are needed to 318 meet convergence, `NoConvergence` is raised. 319 ftol : float, optional 320 Relative tolerance for the residual. If omitted, not used. 321 fatol : float, optional 322 Absolute tolerance (in max-norm) for the residual. 323 If omitted, default is 6e-6. 324 xtol : float, optional 325 Relative minimum step size. If omitted, not used. 326 xatol : float, optional 327 Absolute minimum step size, as determined from the Jacobian 328 approximation. If the step size is smaller than this, optimization 329 is terminated as successful. If omitted, not used. 330 tol_norm : function(vector) -> scalar, optional 331 Norm to use in convergence check. Default is the maximum norm. 332 line_search : {None, 'armijo' (default), 'wolfe'}, optional 333 Which type of a line search to use to determine the step size in 334 the direction given by the Jacobian approximation. Defaults to 335 'armijo'. 336 jac_options : dict, optional 337 Options for the respective Jacobian approximation. 338 alpha : float, optional 339 Initial guess for the Jacobian is (-1/alpha). 340 reduction_method : str or tuple, optional 341 Method used in ensuring that the rank of the Broyden 342 matrix stays low. Can either be a string giving the 343 name of the method, or a tuple of the form ``(method, 344 param1, param2, ...)`` that gives the name of the 345 method and values for additional parameters. 346 347 Methods available: 348 349 - ``restart`` 350 Drop all matrix columns. Has no 351 extra parameters. 352 - ``simple`` 353 Drop oldest matrix column. Has no 354 extra parameters. 355 - ``svd`` 356 Keep only the most significant SVD 357 components. 358 359 Extra parameters: 360 361 - ``to_retain`` 362 Number of SVD components to 363 retain when rank reduction is done. 364 Default is ``max_rank - 2``. 365 max_rank : int, optional 366 Maximum rank for the Broyden matrix. 367 Default is infinity (i.e., no rank reduction). 368 """ 369 pass 370 371def _root_broyden2_doc(): 372 """ 373 Options 374 ------- 375 nit : int, optional 376 Number of iterations to make. If omitted (default), make as many 377 as required to meet tolerances. 378 disp : bool, optional 379 Print status to stdout on every iteration. 380 maxiter : int, optional 381 Maximum number of iterations to make. If more are needed to 382 meet convergence, `NoConvergence` is raised. 383 ftol : float, optional 384 Relative tolerance for the residual. If omitted, not used. 385 fatol : float, optional 386 Absolute tolerance (in max-norm) for the residual. 387 If omitted, default is 6e-6. 388 xtol : float, optional 389 Relative minimum step size. If omitted, not used. 390 xatol : float, optional 391 Absolute minimum step size, as determined from the Jacobian 392 approximation. If the step size is smaller than this, optimization 393 is terminated as successful. If omitted, not used. 394 tol_norm : function(vector) -> scalar, optional 395 Norm to use in convergence check. Default is the maximum norm. 396 line_search : {None, 'armijo' (default), 'wolfe'}, optional 397 Which type of a line search to use to determine the step size in 398 the direction given by the Jacobian approximation. Defaults to 399 'armijo'. 400 jac_options : dict, optional 401 Options for the respective Jacobian approximation. 402 403 alpha : float, optional 404 Initial guess for the Jacobian is (-1/alpha). 405 reduction_method : str or tuple, optional 406 Method used in ensuring that the rank of the Broyden 407 matrix stays low. Can either be a string giving the 408 name of the method, or a tuple of the form ``(method, 409 param1, param2, ...)`` that gives the name of the 410 method and values for additional parameters. 411 412 Methods available: 413 414 - ``restart`` 415 Drop all matrix columns. Has no 416 extra parameters. 417 - ``simple`` 418 Drop oldest matrix column. Has no 419 extra parameters. 420 - ``svd`` 421 Keep only the most significant SVD 422 components. 423 424 Extra parameters: 425 426 - ``to_retain`` 427 Number of SVD components to 428 retain when rank reduction is done. 429 Default is ``max_rank - 2``. 430 max_rank : int, optional 431 Maximum rank for the Broyden matrix. 432 Default is infinity (i.e., no rank reduction). 433 """ 434 pass 435 436def _root_anderson_doc(): 437 """ 438 Options 439 ------- 440 nit : int, optional 441 Number of iterations to make. If omitted (default), make as many 442 as required to meet tolerances. 443 disp : bool, optional 444 Print status to stdout on every iteration. 445 maxiter : int, optional 446 Maximum number of iterations to make. If more are needed to 447 meet convergence, `NoConvergence` is raised. 448 ftol : float, optional 449 Relative tolerance for the residual. If omitted, not used. 450 fatol : float, optional 451 Absolute tolerance (in max-norm) for the residual. 452 If omitted, default is 6e-6. 453 xtol : float, optional 454 Relative minimum step size. If omitted, not used. 455 xatol : float, optional 456 Absolute minimum step size, as determined from the Jacobian 457 approximation. If the step size is smaller than this, optimization 458 is terminated as successful. If omitted, not used. 459 tol_norm : function(vector) -> scalar, optional 460 Norm to use in convergence check. Default is the maximum norm. 461 line_search : {None, 'armijo' (default), 'wolfe'}, optional 462 Which type of a line search to use to determine the step size in 463 the direction given by the Jacobian approximation. Defaults to 464 'armijo'. 465 jac_options : dict, optional 466 Options for the respective Jacobian approximation. 467 468 alpha : float, optional 469 Initial guess for the Jacobian is (-1/alpha). 470 M : float, optional 471 Number of previous vectors to retain. Defaults to 5. 472 w0 : float, optional 473 Regularization parameter for numerical stability. 474 Compared to unity, good values of the order of 0.01. 475 """ 476 pass 477 478def _root_linearmixing_doc(): 479 """ 480 Options 481 ------- 482 nit : int, optional 483 Number of iterations to make. If omitted (default), make as many 484 as required to meet tolerances. 485 disp : bool, optional 486 Print status to stdout on every iteration. 487 maxiter : int, optional 488 Maximum number of iterations to make. If more are needed to 489 meet convergence, ``NoConvergence`` is raised. 490 ftol : float, optional 491 Relative tolerance for the residual. If omitted, not used. 492 fatol : float, optional 493 Absolute tolerance (in max-norm) for the residual. 494 If omitted, default is 6e-6. 495 xtol : float, optional 496 Relative minimum step size. If omitted, not used. 497 xatol : float, optional 498 Absolute minimum step size, as determined from the Jacobian 499 approximation. If the step size is smaller than this, optimization 500 is terminated as successful. If omitted, not used. 501 tol_norm : function(vector) -> scalar, optional 502 Norm to use in convergence check. Default is the maximum norm. 503 line_search : {None, 'armijo' (default), 'wolfe'}, optional 504 Which type of a line search to use to determine the step size in 505 the direction given by the Jacobian approximation. Defaults to 506 'armijo'. 507 jac_options : dict, optional 508 Options for the respective Jacobian approximation. 509 510 alpha : float, optional 511 initial guess for the jacobian is (-1/alpha). 512 """ 513 pass 514 515def _root_diagbroyden_doc(): 516 """ 517 Options 518 ------- 519 nit : int, optional 520 Number of iterations to make. If omitted (default), make as many 521 as required to meet tolerances. 522 disp : bool, optional 523 Print status to stdout on every iteration. 524 maxiter : int, optional 525 Maximum number of iterations to make. If more are needed to 526 meet convergence, `NoConvergence` is raised. 527 ftol : float, optional 528 Relative tolerance for the residual. If omitted, not used. 529 fatol : float, optional 530 Absolute tolerance (in max-norm) for the residual. 531 If omitted, default is 6e-6. 532 xtol : float, optional 533 Relative minimum step size. If omitted, not used. 534 xatol : float, optional 535 Absolute minimum step size, as determined from the Jacobian 536 approximation. If the step size is smaller than this, optimization 537 is terminated as successful. If omitted, not used. 538 tol_norm : function(vector) -> scalar, optional 539 Norm to use in convergence check. Default is the maximum norm. 540 line_search : {None, 'armijo' (default), 'wolfe'}, optional 541 Which type of a line search to use to determine the step size in 542 the direction given by the Jacobian approximation. Defaults to 543 'armijo'. 544 jac_options : dict, optional 545 Options for the respective Jacobian approximation. 546 547 alpha : float, optional 548 initial guess for the jacobian is (-1/alpha). 549 """ 550 pass 551 552def _root_excitingmixing_doc(): 553 """ 554 Options 555 ------- 556 nit : int, optional 557 Number of iterations to make. If omitted (default), make as many 558 as required to meet tolerances. 559 disp : bool, optional 560 Print status to stdout on every iteration. 561 maxiter : int, optional 562 Maximum number of iterations to make. If more are needed to 563 meet convergence, `NoConvergence` is raised. 564 ftol : float, optional 565 Relative tolerance for the residual. If omitted, not used. 566 fatol : float, optional 567 Absolute tolerance (in max-norm) for the residual. 568 If omitted, default is 6e-6. 569 xtol : float, optional 570 Relative minimum step size. If omitted, not used. 571 xatol : float, optional 572 Absolute minimum step size, as determined from the Jacobian 573 approximation. If the step size is smaller than this, optimization 574 is terminated as successful. If omitted, not used. 575 tol_norm : function(vector) -> scalar, optional 576 Norm to use in convergence check. Default is the maximum norm. 577 line_search : {None, 'armijo' (default), 'wolfe'}, optional 578 Which type of a line search to use to determine the step size in 579 the direction given by the Jacobian approximation. Defaults to 580 'armijo'. 581 jac_options : dict, optional 582 Options for the respective Jacobian approximation. 583 584 alpha : float, optional 585 Initial Jacobian approximation is (-1/alpha). 586 alphamax : float, optional 587 The entries of the diagonal Jacobian are kept in the range 588 ``[alpha, alphamax]``. 589 """ 590 pass 591 592def _root_krylov_doc(): 593 """ 594 Options 595 ------- 596 nit : int, optional 597 Number of iterations to make. If omitted (default), make as many 598 as required to meet tolerances. 599 disp : bool, optional 600 Print status to stdout on every iteration. 601 maxiter : int, optional 602 Maximum number of iterations to make. If more are needed to 603 meet convergence, `NoConvergence` is raised. 604 ftol : float, optional 605 Relative tolerance for the residual. If omitted, not used. 606 fatol : float, optional 607 Absolute tolerance (in max-norm) for the residual. 608 If omitted, default is 6e-6. 609 xtol : float, optional 610 Relative minimum step size. If omitted, not used. 611 xatol : float, optional 612 Absolute minimum step size, as determined from the Jacobian 613 approximation. If the step size is smaller than this, optimization 614 is terminated as successful. If omitted, not used. 615 tol_norm : function(vector) -> scalar, optional 616 Norm to use in convergence check. Default is the maximum norm. 617 line_search : {None, 'armijo' (default), 'wolfe'}, optional 618 Which type of a line search to use to determine the step size in 619 the direction given by the Jacobian approximation. Defaults to 620 'armijo'. 621 jac_options : dict, optional 622 Options for the respective Jacobian approximation. 623 624 rdiff : float, optional 625 Relative step size to use in numerical differentiation. 626 method : {'lgmres', 'gmres', 'bicgstab', 'cgs', 'minres'} or function 627 Krylov method to use to approximate the Jacobian. 628 Can be a string, or a function implementing the same 629 interface as the iterative solvers in 630 `scipy.sparse.linalg`. 631 632 The default is `scipy.sparse.linalg.lgmres`. 633 inner_M : LinearOperator or InverseJacobian 634 Preconditioner for the inner Krylov iteration. 635 Note that you can use also inverse Jacobians as (adaptive) 636 preconditioners. For example, 637 638 >>> jac = BroydenFirst() 639 >>> kjac = KrylovJacobian(inner_M=jac.inverse). 640 641 If the preconditioner has a method named 'update', it will 642 be called as ``update(x, f)`` after each nonlinear step, 643 with ``x`` giving the current point, and ``f`` the current 644 function value. 645 inner_tol, inner_maxiter, ... 646 Parameters to pass on to the "inner" Krylov solver. 647 See `scipy.sparse.linalg.gmres` for details. 648 outer_k : int, optional 649 Size of the subspace kept across LGMRES nonlinear 650 iterations. 651 652 See `scipy.sparse.linalg.lgmres` for details. 653 """ 654 pass 655