Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_optimize.py: 6%
1351 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1#__docformat__ = "restructuredtext en"
2# ******NOTICE***************
3# optimize.py module by Travis E. Oliphant
4#
5# You may copy and use this module as you see fit with no
6# guarantee implied provided you keep this notice in all copies.
7# *****END NOTICE************
9# A collection of optimization algorithms. Version 0.5
10# CHANGES
11# Added fminbound (July 2001)
12# Added brute (Aug. 2002)
13# Finished line search satisfying strong Wolfe conditions (Mar. 2004)
14# Updated strong Wolfe conditions line search to use
15# cubic-interpolation (Mar. 2004)
18# Minimization routines
20__all__ = ['fmin', 'fmin_powell', 'fmin_bfgs', 'fmin_ncg', 'fmin_cg',
21 'fminbound', 'brent', 'golden', 'bracket', 'rosen', 'rosen_der',
22 'rosen_hess', 'rosen_hess_prod', 'brute', 'approx_fprime',
23 'line_search', 'check_grad', 'OptimizeResult', 'show_options',
24 'OptimizeWarning']
26__docformat__ = "restructuredtext en"
28import warnings
29import sys
30from numpy import (atleast_1d, eye, argmin, zeros, shape, squeeze,
31 asarray, sqrt, Inf, asfarray)
32import numpy as np
33from scipy.sparse.linalg import LinearOperator
34from ._linesearch import (line_search_wolfe1, line_search_wolfe2,
35 line_search_wolfe2 as line_search,
36 LineSearchWarning)
37from ._numdiff import approx_derivative
38from ._hessian_update_strategy import HessianUpdateStrategy
39from scipy._lib._util import getfullargspec_no_self as _getfullargspec
40from scipy._lib._util import MapWrapper, check_random_state
41from scipy.optimize._differentiable_functions import ScalarFunction, FD_METHODS
44# standard status messages of optimizers
45_status_message = {'success': 'Optimization terminated successfully.',
46 'maxfev': 'Maximum number of function evaluations has '
47 'been exceeded.',
48 'maxiter': 'Maximum number of iterations has been '
49 'exceeded.',
50 'pr_loss': 'Desired error not necessarily achieved due '
51 'to precision loss.',
52 'nan': 'NaN result encountered.',
53 'out_of_bounds': 'The result is outside of the provided '
54 'bounds.'}
57class MemoizeJac:
58 """ Decorator that caches the return values of a function returning `(fun, grad)`
59 each time it is called. """
61 def __init__(self, fun):
62 self.fun = fun
63 self.jac = None
64 self._value = None
65 self.x = None
67 def _compute_if_needed(self, x, *args):
68 if not np.all(x == self.x) or self._value is None or self.jac is None:
69 self.x = np.asarray(x).copy()
70 fg = self.fun(x, *args)
71 self.jac = fg[1]
72 self._value = fg[0]
74 def __call__(self, x, *args):
75 """ returns the function value """
76 self._compute_if_needed(x, *args)
77 return self._value
79 def derivative(self, x, *args):
80 self._compute_if_needed(x, *args)
81 return self.jac
84def _indenter(s, n=0):
85 """
86 Ensures that lines after the first are indented by the specified amount
87 """
88 split = s.split("\n")
89 indent = " "*n
90 return ("\n" + indent).join(split)
93def _float_formatter_10(x):
94 """
95 Returns a string representation of a float with exactly ten characters
96 """
97 if np.isposinf(x):
98 return " inf"
99 elif np.isneginf(x):
100 return " -inf"
101 elif np.isnan(x):
102 return " nan"
103 return np.format_float_scientific(x, precision=3, pad_left=2, unique=False)
106def _dict_formatter(d, n=0, mplus=1, sorter=None):
107 """
108 Pretty printer for dictionaries
110 `n` keeps track of the starting indentation;
111 lines are indented by this much after a line break.
112 `mplus` is additional left padding applied to keys
113 """
114 if isinstance(d, dict):
115 m = max(map(len, list(d.keys()))) + mplus # width to print keys
116 s = '\n'.join([k.rjust(m) + ': ' + # right justified, width m
117 _indenter(_dict_formatter(v, m+n+2, 0, sorter), m+2)
118 for k, v in sorter(d)]) # +2 for ': '
119 else:
120 # By default, NumPy arrays print with linewidth=76. `n` is
121 # the indent at which a line begins printing, so it is subtracted
122 # from the default to avoid exceeding 76 characters total.
123 # `edgeitems` is the number of elements to include before and after
124 # ellipses when arrays are not shown in full.
125 # `threshold` is the maximum number of elements for which an
126 # array is shown in full.
127 # These values tend to work well for use with OptimizeResult.
128 with np.printoptions(linewidth=76-n, edgeitems=2, threshold=12,
129 formatter={'float_kind': _float_formatter_10}):
130 s = str(d)
131 return s
134class OptimizeResult(dict):
135 """ Represents the optimization result.
137 Attributes
138 ----------
139 x : ndarray
140 The solution of the optimization.
141 success : bool
142 Whether or not the optimizer exited successfully.
143 status : int
144 Termination status of the optimizer. Its value depends on the
145 underlying solver. Refer to `message` for details.
146 message : str
147 Description of the cause of the termination.
148 fun, jac, hess: ndarray
149 Values of objective function, its Jacobian and its Hessian (if
150 available). The Hessians may be approximations, see the documentation
151 of the function in question.
152 hess_inv : object
153 Inverse of the objective function's Hessian; may be an approximation.
154 Not available for all solvers. The type of this attribute may be
155 either np.ndarray or scipy.sparse.linalg.LinearOperator.
156 nfev, njev, nhev : int
157 Number of evaluations of the objective functions and of its
158 Jacobian and Hessian.
159 nit : int
160 Number of iterations performed by the optimizer.
161 maxcv : float
162 The maximum constraint violation.
164 Notes
165 -----
166 `OptimizeResult` may have additional attributes not listed here depending
167 on the specific solver being used. Since this class is essentially a
168 subclass of dict with attribute accessors, one can see which
169 attributes are available using the `OptimizeResult.keys` method.
170 """
172 def __getattr__(self, name):
173 try:
174 return self[name]
175 except KeyError as e:
176 raise AttributeError(name) from e
178 __setattr__ = dict.__setitem__
179 __delattr__ = dict.__delitem__
181 def __repr__(self):
182 order_keys = ['message', 'success', 'status', 'fun', 'funl', 'x', 'xl',
183 'col_ind', 'nit', 'lower', 'upper', 'eqlin', 'ineqlin']
184 # 'slack', 'con' are redundant with residuals
185 # 'crossover_nit' is probably not interesting to most users
186 omit_keys = {'slack', 'con', 'crossover_nit'}
188 def key(item):
189 try:
190 return order_keys.index(item[0].lower())
191 except ValueError: # item not in list
192 return np.inf
194 def omit_redundant(items):
195 for item in items:
196 if item[0] in omit_keys:
197 continue
198 yield item
200 def item_sorter(d):
201 return sorted(omit_redundant(d.items()), key=key)
203 if self.keys():
204 return _dict_formatter(self, sorter=item_sorter)
205 else:
206 return self.__class__.__name__ + "()"
208 def __dir__(self):
209 return list(self.keys())
212class OptimizeWarning(UserWarning):
213 pass
216def _check_unknown_options(unknown_options):
217 if unknown_options:
218 msg = ", ".join(map(str, unknown_options.keys()))
219 # Stack level 4: this is called from _minimize_*, which is
220 # called from another function in SciPy. Level 4 is the first
221 # level in user code.
222 warnings.warn("Unknown solver options: %s" % msg, OptimizeWarning, 4)
225def is_finite_scalar(x):
226 """Test whether `x` is either a finite scalar or a finite array scalar.
228 """
229 return np.size(x) == 1 and np.isfinite(x)
232_epsilon = sqrt(np.finfo(float).eps)
235def vecnorm(x, ord=2):
236 if ord == Inf:
237 return np.amax(np.abs(x))
238 elif ord == -Inf:
239 return np.amin(np.abs(x))
240 else:
241 return np.sum(np.abs(x)**ord, axis=0)**(1.0 / ord)
244def _prepare_scalar_function(fun, x0, jac=None, args=(), bounds=None,
245 epsilon=None, finite_diff_rel_step=None,
246 hess=None):
247 """
248 Creates a ScalarFunction object for use with scalar minimizers
249 (BFGS/LBFGSB/SLSQP/TNC/CG/etc).
251 Parameters
252 ----------
253 fun : callable
254 The objective function to be minimized.
256 ``fun(x, *args) -> float``
258 where ``x`` is an 1-D array with shape (n,) and ``args``
259 is a tuple of the fixed parameters needed to completely
260 specify the function.
261 x0 : ndarray, shape (n,)
262 Initial guess. Array of real elements of size (n,),
263 where 'n' is the number of independent variables.
264 jac : {callable, '2-point', '3-point', 'cs', None}, optional
265 Method for computing the gradient vector. If it is a callable, it
266 should be a function that returns the gradient vector:
268 ``jac(x, *args) -> array_like, shape (n,)``
270 If one of `{'2-point', '3-point', 'cs'}` is selected then the gradient
271 is calculated with a relative step for finite differences. If `None`,
272 then two-point finite differences with an absolute step is used.
273 args : tuple, optional
274 Extra arguments passed to the objective function and its
275 derivatives (`fun`, `jac` functions).
276 bounds : sequence, optional
277 Bounds on variables. 'new-style' bounds are required.
278 eps : float or ndarray
279 If `jac is None` the absolute step size used for numerical
280 approximation of the jacobian via forward differences.
281 finite_diff_rel_step : None or array_like, optional
282 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
283 use for numerical approximation of the jacobian. The absolute step
284 size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``,
285 possibly adjusted to fit into the bounds. For ``method='3-point'``
286 the sign of `h` is ignored. If None (default) then step is selected
287 automatically.
288 hess : {callable, '2-point', '3-point', 'cs', None}
289 Computes the Hessian matrix. If it is callable, it should return the
290 Hessian matrix:
292 ``hess(x, *args) -> {LinearOperator, spmatrix, array}, (n, n)``
294 Alternatively, the keywords {'2-point', '3-point', 'cs'} select a
295 finite difference scheme for numerical estimation.
296 Whenever the gradient is estimated via finite-differences, the Hessian
297 cannot be estimated with options {'2-point', '3-point', 'cs'} and needs
298 to be estimated using one of the quasi-Newton strategies.
300 Returns
301 -------
302 sf : ScalarFunction
303 """
304 if callable(jac):
305 grad = jac
306 elif jac in FD_METHODS:
307 # epsilon is set to None so that ScalarFunction is made to use
308 # rel_step
309 epsilon = None
310 grad = jac
311 else:
312 # default (jac is None) is to do 2-point finite differences with
313 # absolute step size. ScalarFunction has to be provided an
314 # epsilon value that is not None to use absolute steps. This is
315 # normally the case from most _minimize* methods.
316 grad = '2-point'
317 epsilon = epsilon
319 if hess is None:
320 # ScalarFunction requires something for hess, so we give a dummy
321 # implementation here if nothing is provided, return a value of None
322 # so that downstream minimisers halt. The results of `fun.hess`
323 # should not be used.
324 def hess(x, *args):
325 return None
327 if bounds is None:
328 bounds = (-np.inf, np.inf)
330 # ScalarFunction caches. Reuse of fun(x) during grad
331 # calculation reduces overall function evaluations.
332 sf = ScalarFunction(fun, x0, args, grad, hess,
333 finite_diff_rel_step, bounds, epsilon=epsilon)
335 return sf
338def _clip_x_for_func(func, bounds):
339 # ensures that x values sent to func are clipped to bounds
341 # this is used as a mitigation for gh11403, slsqp/tnc sometimes
342 # suggest a move that is outside the limits by 1 or 2 ULP. This
343 # unclean fix makes sure x is strictly within bounds.
344 def eval(x):
345 x = _check_clip_x(x, bounds)
346 return func(x)
348 return eval
351def _check_clip_x(x, bounds):
352 if (x < bounds[0]).any() or (x > bounds[1]).any():
353 warnings.warn("Values in x were outside bounds during a "
354 "minimize step, clipping to bounds", RuntimeWarning)
355 x = np.clip(x, bounds[0], bounds[1])
356 return x
358 return x
361def rosen(x):
362 """
363 The Rosenbrock function.
365 The function computed is::
367 sum(100.0*(x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0)
369 Parameters
370 ----------
371 x : array_like
372 1-D array of points at which the Rosenbrock function is to be computed.
374 Returns
375 -------
376 f : float
377 The value of the Rosenbrock function.
379 See Also
380 --------
381 rosen_der, rosen_hess, rosen_hess_prod
383 Examples
384 --------
385 >>> import numpy as np
386 >>> from scipy.optimize import rosen
387 >>> X = 0.1 * np.arange(10)
388 >>> rosen(X)
389 76.56
391 For higher-dimensional input ``rosen`` broadcasts.
392 In the following example, we use this to plot a 2D landscape.
393 Note that ``rosen_hess`` does not broadcast in this manner.
395 >>> import matplotlib.pyplot as plt
396 >>> from mpl_toolkits.mplot3d import Axes3D
397 >>> x = np.linspace(-1, 1, 50)
398 >>> X, Y = np.meshgrid(x, x)
399 >>> ax = plt.subplot(111, projection='3d')
400 >>> ax.plot_surface(X, Y, rosen([X, Y]))
401 >>> plt.show()
402 """
403 x = asarray(x)
404 r = np.sum(100.0 * (x[1:] - x[:-1]**2.0)**2.0 + (1 - x[:-1])**2.0,
405 axis=0)
406 return r
409def rosen_der(x):
410 """
411 The derivative (i.e. gradient) of the Rosenbrock function.
413 Parameters
414 ----------
415 x : array_like
416 1-D array of points at which the derivative is to be computed.
418 Returns
419 -------
420 rosen_der : (N,) ndarray
421 The gradient of the Rosenbrock function at `x`.
423 See Also
424 --------
425 rosen, rosen_hess, rosen_hess_prod
427 Examples
428 --------
429 >>> import numpy as np
430 >>> from scipy.optimize import rosen_der
431 >>> X = 0.1 * np.arange(9)
432 >>> rosen_der(X)
433 array([ -2. , 10.6, 15.6, 13.4, 6.4, -3. , -12.4, -19.4, 62. ])
435 """
436 x = asarray(x)
437 xm = x[1:-1]
438 xm_m1 = x[:-2]
439 xm_p1 = x[2:]
440 der = np.zeros_like(x)
441 der[1:-1] = (200 * (xm - xm_m1**2) -
442 400 * (xm_p1 - xm**2) * xm - 2 * (1 - xm))
443 der[0] = -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0])
444 der[-1] = 200 * (x[-1] - x[-2]**2)
445 return der
448def rosen_hess(x):
449 """
450 The Hessian matrix of the Rosenbrock function.
452 Parameters
453 ----------
454 x : array_like
455 1-D array of points at which the Hessian matrix is to be computed.
457 Returns
458 -------
459 rosen_hess : ndarray
460 The Hessian matrix of the Rosenbrock function at `x`.
462 See Also
463 --------
464 rosen, rosen_der, rosen_hess_prod
466 Examples
467 --------
468 >>> import numpy as np
469 >>> from scipy.optimize import rosen_hess
470 >>> X = 0.1 * np.arange(4)
471 >>> rosen_hess(X)
472 array([[-38., 0., 0., 0.],
473 [ 0., 134., -40., 0.],
474 [ 0., -40., 130., -80.],
475 [ 0., 0., -80., 200.]])
477 """
478 x = atleast_1d(x)
479 H = np.diag(-400 * x[:-1], 1) - np.diag(400 * x[:-1], -1)
480 diagonal = np.zeros(len(x), dtype=x.dtype)
481 diagonal[0] = 1200 * x[0]**2 - 400 * x[1] + 2
482 diagonal[-1] = 200
483 diagonal[1:-1] = 202 + 1200 * x[1:-1]**2 - 400 * x[2:]
484 H = H + np.diag(diagonal)
485 return H
488def rosen_hess_prod(x, p):
489 """
490 Product of the Hessian matrix of the Rosenbrock function with a vector.
492 Parameters
493 ----------
494 x : array_like
495 1-D array of points at which the Hessian matrix is to be computed.
496 p : array_like
497 1-D array, the vector to be multiplied by the Hessian matrix.
499 Returns
500 -------
501 rosen_hess_prod : ndarray
502 The Hessian matrix of the Rosenbrock function at `x` multiplied
503 by the vector `p`.
505 See Also
506 --------
507 rosen, rosen_der, rosen_hess
509 Examples
510 --------
511 >>> import numpy as np
512 >>> from scipy.optimize import rosen_hess_prod
513 >>> X = 0.1 * np.arange(9)
514 >>> p = 0.5 * np.arange(9)
515 >>> rosen_hess_prod(X, p)
516 array([ -0., 27., -10., -95., -192., -265., -278., -195., -180.])
518 """
519 x = atleast_1d(x)
520 Hp = np.zeros(len(x), dtype=x.dtype)
521 Hp[0] = (1200 * x[0]**2 - 400 * x[1] + 2) * p[0] - 400 * x[0] * p[1]
522 Hp[1:-1] = (-400 * x[:-2] * p[:-2] +
523 (202 + 1200 * x[1:-1]**2 - 400 * x[2:]) * p[1:-1] -
524 400 * x[1:-1] * p[2:])
525 Hp[-1] = -400 * x[-2] * p[-2] + 200*p[-1]
526 return Hp
529def _wrap_scalar_function(function, args):
530 # wraps a minimizer function to count number of evaluations
531 # and to easily provide an args kwd.
532 ncalls = [0]
533 if function is None:
534 return ncalls, None
536 def function_wrapper(x, *wrapper_args):
537 ncalls[0] += 1
538 # A copy of x is sent to the user function (gh13740)
539 fx = function(np.copy(x), *(wrapper_args + args))
540 # Ideally, we'd like to a have a true scalar returned from f(x). For
541 # backwards-compatibility, also allow np.array([1.3]), np.array([[1.3]]) etc.
542 if not np.isscalar(fx):
543 try:
544 fx = np.asarray(fx).item()
545 except (TypeError, ValueError) as e:
546 raise ValueError("The user-provided objective function "
547 "must return a scalar value.") from e
548 return fx
550 return ncalls, function_wrapper
553class _MaxFuncCallError(RuntimeError):
554 pass
557def _wrap_scalar_function_maxfun_validation(function, args, maxfun):
558 # wraps a minimizer function to count number of evaluations
559 # and to easily provide an args kwd.
560 ncalls = [0]
561 if function is None:
562 return ncalls, None
564 def function_wrapper(x, *wrapper_args):
565 if ncalls[0] >= maxfun:
566 raise _MaxFuncCallError("Too many function calls")
567 ncalls[0] += 1
568 # A copy of x is sent to the user function (gh13740)
569 fx = function(np.copy(x), *(wrapper_args + args))
570 # Ideally, we'd like to a have a true scalar returned from f(x). For
571 # backwards-compatibility, also allow np.array([1.3]),
572 # np.array([[1.3]]) etc.
573 if not np.isscalar(fx):
574 try:
575 fx = np.asarray(fx).item()
576 except (TypeError, ValueError) as e:
577 raise ValueError("The user-provided objective function "
578 "must return a scalar value.") from e
579 return fx
581 return ncalls, function_wrapper
584def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None,
585 full_output=0, disp=1, retall=0, callback=None, initial_simplex=None):
586 """
587 Minimize a function using the downhill simplex algorithm.
589 This algorithm only uses function values, not derivatives or second
590 derivatives.
592 Parameters
593 ----------
594 func : callable func(x,*args)
595 The objective function to be minimized.
596 x0 : ndarray
597 Initial guess.
598 args : tuple, optional
599 Extra arguments passed to func, i.e., ``f(x,*args)``.
600 xtol : float, optional
601 Absolute error in xopt between iterations that is acceptable for
602 convergence.
603 ftol : number, optional
604 Absolute error in func(xopt) between iterations that is acceptable for
605 convergence.
606 maxiter : int, optional
607 Maximum number of iterations to perform.
608 maxfun : number, optional
609 Maximum number of function evaluations to make.
610 full_output : bool, optional
611 Set to True if fopt and warnflag outputs are desired.
612 disp : bool, optional
613 Set to True to print convergence messages.
614 retall : bool, optional
615 Set to True to return list of solutions at each iteration.
616 callback : callable, optional
617 Called after each iteration, as callback(xk), where xk is the
618 current parameter vector.
619 initial_simplex : array_like of shape (N + 1, N), optional
620 Initial simplex. If given, overrides `x0`.
621 ``initial_simplex[j,:]`` should contain the coordinates of
622 the jth vertex of the ``N+1`` vertices in the simplex, where
623 ``N`` is the dimension.
625 Returns
626 -------
627 xopt : ndarray
628 Parameter that minimizes function.
629 fopt : float
630 Value of function at minimum: ``fopt = func(xopt)``.
631 iter : int
632 Number of iterations performed.
633 funcalls : int
634 Number of function calls made.
635 warnflag : int
636 1 : Maximum number of function evaluations made.
637 2 : Maximum number of iterations reached.
638 allvecs : list
639 Solution at each iteration.
641 See also
642 --------
643 minimize: Interface to minimization algorithms for multivariate
644 functions. See the 'Nelder-Mead' `method` in particular.
646 Notes
647 -----
648 Uses a Nelder-Mead simplex algorithm to find the minimum of function of
649 one or more variables.
651 This algorithm has a long history of successful use in applications.
652 But it will usually be slower than an algorithm that uses first or
653 second derivative information. In practice, it can have poor
654 performance in high-dimensional problems and is not robust to
655 minimizing complicated functions. Additionally, there currently is no
656 complete theory describing when the algorithm will successfully
657 converge to the minimum, or how fast it will if it does. Both the ftol and
658 xtol criteria must be met for convergence.
660 Examples
661 --------
662 >>> def f(x):
663 ... return x**2
665 >>> from scipy import optimize
667 >>> minimum = optimize.fmin(f, 1)
668 Optimization terminated successfully.
669 Current function value: 0.000000
670 Iterations: 17
671 Function evaluations: 34
672 >>> minimum[0]
673 -8.8817841970012523e-16
675 References
676 ----------
677 .. [1] Nelder, J.A. and Mead, R. (1965), "A simplex method for function
678 minimization", The Computer Journal, 7, pp. 308-313
680 .. [2] Wright, M.H. (1996), "Direct Search Methods: Once Scorned, Now
681 Respectable", in Numerical Analysis 1995, Proceedings of the
682 1995 Dundee Biennial Conference in Numerical Analysis, D.F.
683 Griffiths and G.A. Watson (Eds.), Addison Wesley Longman,
684 Harlow, UK, pp. 191-208.
686 """
687 opts = {'xatol': xtol,
688 'fatol': ftol,
689 'maxiter': maxiter,
690 'maxfev': maxfun,
691 'disp': disp,
692 'return_all': retall,
693 'initial_simplex': initial_simplex}
695 res = _minimize_neldermead(func, x0, args, callback=callback, **opts)
696 if full_output:
697 retlist = res['x'], res['fun'], res['nit'], res['nfev'], res['status']
698 if retall:
699 retlist += (res['allvecs'], )
700 return retlist
701 else:
702 if retall:
703 return res['x'], res['allvecs']
704 else:
705 return res['x']
708def _minimize_neldermead(func, x0, args=(), callback=None,
709 maxiter=None, maxfev=None, disp=False,
710 return_all=False, initial_simplex=None,
711 xatol=1e-4, fatol=1e-4, adaptive=False, bounds=None,
712 **unknown_options):
713 """
714 Minimization of scalar function of one or more variables using the
715 Nelder-Mead algorithm.
717 Options
718 -------
719 disp : bool
720 Set to True to print convergence messages.
721 maxiter, maxfev : int
722 Maximum allowed number of iterations and function evaluations.
723 Will default to ``N*200``, where ``N`` is the number of
724 variables, if neither `maxiter` or `maxfev` is set. If both
725 `maxiter` and `maxfev` are set, minimization will stop at the
726 first reached.
727 return_all : bool, optional
728 Set to True to return a list of the best solution at each of the
729 iterations.
730 initial_simplex : array_like of shape (N + 1, N)
731 Initial simplex. If given, overrides `x0`.
732 ``initial_simplex[j,:]`` should contain the coordinates of
733 the jth vertex of the ``N+1`` vertices in the simplex, where
734 ``N`` is the dimension.
735 xatol : float, optional
736 Absolute error in xopt between iterations that is acceptable for
737 convergence.
738 fatol : number, optional
739 Absolute error in func(xopt) between iterations that is acceptable for
740 convergence.
741 adaptive : bool, optional
742 Adapt algorithm parameters to dimensionality of problem. Useful for
743 high-dimensional minimization [1]_.
744 bounds : sequence or `Bounds`, optional
745 Bounds on variables. There are two ways to specify the bounds:
747 1. Instance of `Bounds` class.
748 2. Sequence of ``(min, max)`` pairs for each element in `x`. None
749 is used to specify no bound.
751 Note that this just clips all vertices in simplex based on
752 the bounds.
754 References
755 ----------
756 .. [1] Gao, F. and Han, L.
757 Implementing the Nelder-Mead simplex algorithm with adaptive
758 parameters. 2012. Computational Optimization and Applications.
759 51:1, pp. 259-277
761 """
762 _check_unknown_options(unknown_options)
763 maxfun = maxfev
764 retall = return_all
766 x0 = asfarray(x0).flatten()
768 if adaptive:
769 dim = float(len(x0))
770 rho = 1
771 chi = 1 + 2/dim
772 psi = 0.75 - 1/(2*dim)
773 sigma = 1 - 1/dim
774 else:
775 rho = 1
776 chi = 2
777 psi = 0.5
778 sigma = 0.5
780 nonzdelt = 0.05
781 zdelt = 0.00025
783 if bounds is not None:
784 lower_bound, upper_bound = bounds.lb, bounds.ub
785 # check bounds
786 if (lower_bound > upper_bound).any():
787 raise ValueError("Nelder Mead - one of the lower bounds is greater than an upper bound.")
788 if np.any(lower_bound > x0) or np.any(x0 > upper_bound):
789 warnings.warn("Initial guess is not within the specified bounds",
790 OptimizeWarning, 3)
792 if bounds is not None:
793 x0 = np.clip(x0, lower_bound, upper_bound)
795 if initial_simplex is None:
796 N = len(x0)
798 sim = np.empty((N + 1, N), dtype=x0.dtype)
799 sim[0] = x0
800 for k in range(N):
801 y = np.array(x0, copy=True)
802 if y[k] != 0:
803 y[k] = (1 + nonzdelt)*y[k]
804 else:
805 y[k] = zdelt
806 sim[k + 1] = y
807 else:
808 sim = np.asfarray(initial_simplex).copy()
809 if sim.ndim != 2 or sim.shape[0] != sim.shape[1] + 1:
810 raise ValueError("`initial_simplex` should be an array of shape (N+1,N)")
811 if len(x0) != sim.shape[1]:
812 raise ValueError("Size of `initial_simplex` is not consistent with `x0`")
813 N = sim.shape[1]
815 if retall:
816 allvecs = [sim[0]]
818 # If neither are set, then set both to default
819 if maxiter is None and maxfun is None:
820 maxiter = N * 200
821 maxfun = N * 200
822 elif maxiter is None:
823 # Convert remaining Nones, to np.inf, unless the other is np.inf, in
824 # which case use the default to avoid unbounded iteration
825 if maxfun == np.inf:
826 maxiter = N * 200
827 else:
828 maxiter = np.inf
829 elif maxfun is None:
830 if maxiter == np.inf:
831 maxfun = N * 200
832 else:
833 maxfun = np.inf
835 if bounds is not None:
836 sim = np.clip(sim, lower_bound, upper_bound)
838 one2np1 = list(range(1, N + 1))
839 fsim = np.full((N + 1,), np.inf, dtype=float)
841 fcalls, func = _wrap_scalar_function_maxfun_validation(func, args, maxfun)
843 try:
844 for k in range(N + 1):
845 fsim[k] = func(sim[k])
846 except _MaxFuncCallError:
847 pass
848 finally:
849 ind = np.argsort(fsim)
850 sim = np.take(sim, ind, 0)
851 fsim = np.take(fsim, ind, 0)
853 ind = np.argsort(fsim)
854 fsim = np.take(fsim, ind, 0)
855 # sort so sim[0,:] has the lowest function value
856 sim = np.take(sim, ind, 0)
858 iterations = 1
860 while (fcalls[0] < maxfun and iterations < maxiter):
861 try:
862 if (np.max(np.ravel(np.abs(sim[1:] - sim[0]))) <= xatol and
863 np.max(np.abs(fsim[0] - fsim[1:])) <= fatol):
864 break
866 xbar = np.add.reduce(sim[:-1], 0) / N
867 xr = (1 + rho) * xbar - rho * sim[-1]
868 if bounds is not None:
869 xr = np.clip(xr, lower_bound, upper_bound)
870 fxr = func(xr)
871 doshrink = 0
873 if fxr < fsim[0]:
874 xe = (1 + rho * chi) * xbar - rho * chi * sim[-1]
875 if bounds is not None:
876 xe = np.clip(xe, lower_bound, upper_bound)
877 fxe = func(xe)
879 if fxe < fxr:
880 sim[-1] = xe
881 fsim[-1] = fxe
882 else:
883 sim[-1] = xr
884 fsim[-1] = fxr
885 else: # fsim[0] <= fxr
886 if fxr < fsim[-2]:
887 sim[-1] = xr
888 fsim[-1] = fxr
889 else: # fxr >= fsim[-2]
890 # Perform contraction
891 if fxr < fsim[-1]:
892 xc = (1 + psi * rho) * xbar - psi * rho * sim[-1]
893 if bounds is not None:
894 xc = np.clip(xc, lower_bound, upper_bound)
895 fxc = func(xc)
897 if fxc <= fxr:
898 sim[-1] = xc
899 fsim[-1] = fxc
900 else:
901 doshrink = 1
902 else:
903 # Perform an inside contraction
904 xcc = (1 - psi) * xbar + psi * sim[-1]
905 if bounds is not None:
906 xcc = np.clip(xcc, lower_bound, upper_bound)
907 fxcc = func(xcc)
909 if fxcc < fsim[-1]:
910 sim[-1] = xcc
911 fsim[-1] = fxcc
912 else:
913 doshrink = 1
915 if doshrink:
916 for j in one2np1:
917 sim[j] = sim[0] + sigma * (sim[j] - sim[0])
918 if bounds is not None:
919 sim[j] = np.clip(
920 sim[j], lower_bound, upper_bound)
921 fsim[j] = func(sim[j])
922 iterations += 1
923 except _MaxFuncCallError:
924 pass
925 finally:
926 ind = np.argsort(fsim)
927 sim = np.take(sim, ind, 0)
928 fsim = np.take(fsim, ind, 0)
929 if callback is not None:
930 callback(sim[0])
931 if retall:
932 allvecs.append(sim[0])
934 x = sim[0]
935 fval = np.min(fsim)
936 warnflag = 0
938 if fcalls[0] >= maxfun:
939 warnflag = 1
940 msg = _status_message['maxfev']
941 if disp:
942 warnings.warn(msg, RuntimeWarning, 3)
943 elif iterations >= maxiter:
944 warnflag = 2
945 msg = _status_message['maxiter']
946 if disp:
947 warnings.warn(msg, RuntimeWarning, 3)
948 else:
949 msg = _status_message['success']
950 if disp:
951 print(msg)
952 print(" Current function value: %f" % fval)
953 print(" Iterations: %d" % iterations)
954 print(" Function evaluations: %d" % fcalls[0])
956 result = OptimizeResult(fun=fval, nit=iterations, nfev=fcalls[0],
957 status=warnflag, success=(warnflag == 0),
958 message=msg, x=x, final_simplex=(sim, fsim))
959 if retall:
960 result['allvecs'] = allvecs
961 return result
964def approx_fprime(xk, f, epsilon=_epsilon, *args):
965 """Finite difference approximation of the derivatives of a
966 scalar or vector-valued function.
968 If a function maps from :math:`R^n` to :math:`R^m`, its derivatives form
969 an m-by-n matrix
970 called the Jacobian, where an element :math:`(i, j)` is a partial
971 derivative of f[i] with respect to ``xk[j]``.
973 Parameters
974 ----------
975 xk : array_like
976 The coordinate vector at which to determine the gradient of `f`.
977 f : callable
978 Function of which to estimate the derivatives of. Has the signature
979 ``f(xk, *args)`` where `xk` is the argument in the form of a 1-D array
980 and `args` is a tuple of any additional fixed parameters needed to
981 completely specify the function. The argument `xk` passed to this
982 function is an ndarray of shape (n,) (never a scalar even if n=1).
983 It must return a 1-D array_like of shape (m,) or a scalar.
985 .. versionchanged:: 1.9.0
986 `f` is now able to return a 1-D array-like, with the :math:`(m, n)`
987 Jacobian being estimated.
989 epsilon : {float, array_like}, optional
990 Increment to `xk` to use for determining the function gradient.
991 If a scalar, uses the same finite difference delta for all partial
992 derivatives. If an array, should contain one value per element of
993 `xk`. Defaults to ``sqrt(np.finfo(float).eps)``, which is approximately
994 1.49e-08.
995 \\*args : args, optional
996 Any other arguments that are to be passed to `f`.
998 Returns
999 -------
1000 jac : ndarray
1001 The partial derivatives of `f` to `xk`.
1003 See Also
1004 --------
1005 check_grad : Check correctness of gradient function against approx_fprime.
1007 Notes
1008 -----
1009 The function gradient is determined by the forward finite difference
1010 formula::
1012 f(xk[i] + epsilon[i]) - f(xk[i])
1013 f'[i] = ---------------------------------
1014 epsilon[i]
1016 Examples
1017 --------
1018 >>> import numpy as np
1019 >>> from scipy import optimize
1020 >>> def func(x, c0, c1):
1021 ... "Coordinate vector `x` should be an array of size two."
1022 ... return c0 * x[0]**2 + c1*x[1]**2
1024 >>> x = np.ones(2)
1025 >>> c0, c1 = (1, 200)
1026 >>> eps = np.sqrt(np.finfo(float).eps)
1027 >>> optimize.approx_fprime(x, func, [eps, np.sqrt(200) * eps], c0, c1)
1028 array([ 2. , 400.00004198])
1030 """
1031 xk = np.asarray(xk, float)
1032 f0 = f(xk, *args)
1034 return approx_derivative(f, xk, method='2-point', abs_step=epsilon,
1035 args=args, f0=f0)
1038def check_grad(func, grad, x0, *args, epsilon=_epsilon,
1039 direction='all', seed=None):
1040 """Check the correctness of a gradient function by comparing it against a
1041 (forward) finite-difference approximation of the gradient.
1043 Parameters
1044 ----------
1045 func : callable ``func(x0, *args)``
1046 Function whose derivative is to be checked.
1047 grad : callable ``grad(x0, *args)``
1048 Jacobian of `func`.
1049 x0 : ndarray
1050 Points to check `grad` against forward difference approximation of grad
1051 using `func`.
1052 args : \\*args, optional
1053 Extra arguments passed to `func` and `grad`.
1054 epsilon : float, optional
1055 Step size used for the finite difference approximation. It defaults to
1056 ``sqrt(np.finfo(float).eps)``, which is approximately 1.49e-08.
1057 direction : str, optional
1058 If set to ``'random'``, then gradients along a random vector
1059 are used to check `grad` against forward difference approximation
1060 using `func`. By default it is ``'all'``, in which case, all
1061 the one hot direction vectors are considered to check `grad`.
1062 If `func` is a vector valued function then only ``'all'`` can be used.
1063 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
1064 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
1065 singleton is used.
1066 If `seed` is an int, a new ``RandomState`` instance is used,
1067 seeded with `seed`.
1068 If `seed` is already a ``Generator`` or ``RandomState`` instance then
1069 that instance is used.
1070 Specify `seed` for reproducing the return value from this function.
1071 The random numbers generated with this seed affect the random vector
1072 along which gradients are computed to check ``grad``. Note that `seed`
1073 is only used when `direction` argument is set to `'random'`.
1075 Returns
1076 -------
1077 err : float
1078 The square root of the sum of squares (i.e., the 2-norm) of the
1079 difference between ``grad(x0, *args)`` and the finite difference
1080 approximation of `grad` using func at the points `x0`.
1082 See Also
1083 --------
1084 approx_fprime
1086 Examples
1087 --------
1088 >>> import numpy as np
1089 >>> def func(x):
1090 ... return x[0]**2 - 0.5 * x[1]**3
1091 >>> def grad(x):
1092 ... return [2 * x[0], -1.5 * x[1]**2]
1093 >>> from scipy.optimize import check_grad
1094 >>> check_grad(func, grad, [1.5, -1.5])
1095 2.9802322387695312e-08 # may vary
1096 >>> rng = np.random.default_rng()
1097 >>> check_grad(func, grad, [1.5, -1.5],
1098 ... direction='random', seed=rng)
1099 2.9802322387695312e-08
1101 """
1102 step = epsilon
1103 x0 = np.asarray(x0)
1105 def g(w, func, x0, v, *args):
1106 return func(x0 + w*v, *args)
1108 if direction == 'random':
1109 _grad = np.asanyarray(grad(x0, *args))
1110 if _grad.ndim > 1:
1111 raise ValueError("'random' can only be used with scalar valued"
1112 " func")
1113 random_state = check_random_state(seed)
1114 v = random_state.normal(0, 1, size=(x0.shape))
1115 _args = (func, x0, v) + args
1116 _func = g
1117 vars = np.zeros((1,))
1118 analytical_grad = np.dot(_grad, v)
1119 elif direction == 'all':
1120 _args = args
1121 _func = func
1122 vars = x0
1123 analytical_grad = grad(x0, *args)
1124 else:
1125 raise ValueError("{} is not a valid string for "
1126 "``direction`` argument".format(direction))
1128 return np.sqrt(np.sum(np.abs(
1129 (analytical_grad - approx_fprime(vars, _func, step, *_args))**2
1130 )))
1133def approx_fhess_p(x0, p, fprime, epsilon, *args):
1134 # calculate fprime(x0) first, as this may be cached by ScalarFunction
1135 f1 = fprime(*((x0,) + args))
1136 f2 = fprime(*((x0 + epsilon*p,) + args))
1137 return (f2 - f1) / epsilon
1140class _LineSearchError(RuntimeError):
1141 pass
1144def _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval,
1145 **kwargs):
1146 """
1147 Same as line_search_wolfe1, but fall back to line_search_wolfe2 if
1148 suitable step length is not found, and raise an exception if a
1149 suitable step length is not found.
1151 Raises
1152 ------
1153 _LineSearchError
1154 If no suitable step size is found
1156 """
1158 extra_condition = kwargs.pop('extra_condition', None)
1160 ret = line_search_wolfe1(f, fprime, xk, pk, gfk,
1161 old_fval, old_old_fval,
1162 **kwargs)
1164 if ret[0] is not None and extra_condition is not None:
1165 xp1 = xk + ret[0] * pk
1166 if not extra_condition(ret[0], xp1, ret[3], ret[5]):
1167 # Reject step if extra_condition fails
1168 ret = (None,)
1170 if ret[0] is None:
1171 # line search failed: try different one.
1172 with warnings.catch_warnings():
1173 warnings.simplefilter('ignore', LineSearchWarning)
1174 kwargs2 = {}
1175 for key in ('c1', 'c2', 'amax'):
1176 if key in kwargs:
1177 kwargs2[key] = kwargs[key]
1178 ret = line_search_wolfe2(f, fprime, xk, pk, gfk,
1179 old_fval, old_old_fval,
1180 extra_condition=extra_condition,
1181 **kwargs2)
1183 if ret[0] is None:
1184 raise _LineSearchError()
1186 return ret
1189def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf,
1190 epsilon=_epsilon, maxiter=None, full_output=0, disp=1,
1191 retall=0, callback=None, xrtol=0):
1192 """
1193 Minimize a function using the BFGS algorithm.
1195 Parameters
1196 ----------
1197 f : callable ``f(x,*args)``
1198 Objective function to be minimized.
1199 x0 : ndarray
1200 Initial guess.
1201 fprime : callable ``f'(x,*args)``, optional
1202 Gradient of f.
1203 args : tuple, optional
1204 Extra arguments passed to f and fprime.
1205 gtol : float, optional
1206 Terminate successfully if gradient norm is less than `gtol`
1207 norm : float, optional
1208 Order of norm (Inf is max, -Inf is min)
1209 epsilon : int or ndarray, optional
1210 If `fprime` is approximated, use this value for the step size.
1211 callback : callable, optional
1212 An optional user-supplied function to call after each
1213 iteration. Called as ``callback(xk)``, where ``xk`` is the
1214 current parameter vector.
1215 maxiter : int, optional
1216 Maximum number of iterations to perform.
1217 full_output : bool, optional
1218 If True, return ``fopt``, ``func_calls``, ``grad_calls``, and
1219 ``warnflag`` in addition to ``xopt``.
1220 disp : bool, optional
1221 Print convergence message if True.
1222 retall : bool, optional
1223 Return a list of results at each iteration if True.
1224 xrtol : float, default: 0
1225 Relative tolerance for `x`. Terminate successfully if step
1226 size is less than ``xk * xrtol`` where ``xk`` is the current
1227 parameter vector.
1229 Returns
1230 -------
1231 xopt : ndarray
1232 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1233 fopt : float
1234 Minimum value.
1235 gopt : ndarray
1236 Value of gradient at minimum, f'(xopt), which should be near 0.
1237 Bopt : ndarray
1238 Value of 1/f''(xopt), i.e., the inverse Hessian matrix.
1239 func_calls : int
1240 Number of function_calls made.
1241 grad_calls : int
1242 Number of gradient calls made.
1243 warnflag : integer
1244 1 : Maximum number of iterations exceeded.
1245 2 : Gradient and/or function calls not changing.
1246 3 : NaN result encountered.
1247 allvecs : list
1248 The value of `xopt` at each iteration. Only returned if `retall` is
1249 True.
1251 Notes
1252 -----
1253 Optimize the function, `f`, whose gradient is given by `fprime`
1254 using the quasi-Newton method of Broyden, Fletcher, Goldfarb,
1255 and Shanno (BFGS).
1257 See Also
1258 --------
1259 minimize: Interface to minimization algorithms for multivariate
1260 functions. See ``method='BFGS'`` in particular.
1262 References
1263 ----------
1264 Wright, and Nocedal 'Numerical Optimization', 1999, p. 198.
1266 Examples
1267 --------
1268 >>> import numpy as np
1269 >>> from scipy.optimize import fmin_bfgs
1270 >>> def quadratic_cost(x, Q):
1271 ... return x @ Q @ x
1272 ...
1273 >>> x0 = np.array([-3, -4])
1274 >>> cost_weight = np.diag([1., 10.])
1275 >>> # Note that a trailing comma is necessary for a tuple with single element
1276 >>> fmin_bfgs(quadratic_cost, x0, args=(cost_weight,))
1277 Optimization terminated successfully.
1278 Current function value: 0.000000
1279 Iterations: 7 # may vary
1280 Function evaluations: 24 # may vary
1281 Gradient evaluations: 8 # may vary
1282 array([ 2.85169950e-06, -4.61820139e-07])
1284 >>> def quadratic_cost_grad(x, Q):
1285 ... return 2 * Q @ x
1286 ...
1287 >>> fmin_bfgs(quadratic_cost, x0, quadratic_cost_grad, args=(cost_weight,))
1288 Optimization terminated successfully.
1289 Current function value: 0.000000
1290 Iterations: 7
1291 Function evaluations: 8
1292 Gradient evaluations: 8
1293 array([ 2.85916637e-06, -4.54371951e-07])
1295 """
1296 opts = {'gtol': gtol,
1297 'norm': norm,
1298 'eps': epsilon,
1299 'disp': disp,
1300 'maxiter': maxiter,
1301 'return_all': retall}
1303 res = _minimize_bfgs(f, x0, args, fprime, callback=callback, **opts)
1305 if full_output:
1306 retlist = (res['x'], res['fun'], res['jac'], res['hess_inv'],
1307 res['nfev'], res['njev'], res['status'])
1308 if retall:
1309 retlist += (res['allvecs'], )
1310 return retlist
1311 else:
1312 if retall:
1313 return res['x'], res['allvecs']
1314 else:
1315 return res['x']
1318def _minimize_bfgs(fun, x0, args=(), jac=None, callback=None,
1319 gtol=1e-5, norm=Inf, eps=_epsilon, maxiter=None,
1320 disp=False, return_all=False, finite_diff_rel_step=None,
1321 xrtol=0, **unknown_options):
1322 """
1323 Minimization of scalar function of one or more variables using the
1324 BFGS algorithm.
1326 Options
1327 -------
1328 disp : bool
1329 Set to True to print convergence messages.
1330 maxiter : int
1331 Maximum number of iterations to perform.
1332 gtol : float
1333 Terminate successfully if gradient norm is less than `gtol`.
1334 norm : float
1335 Order of norm (Inf is max, -Inf is min).
1336 eps : float or ndarray
1337 If `jac is None` the absolute step size used for numerical
1338 approximation of the jacobian via forward differences.
1339 return_all : bool, optional
1340 Set to True to return a list of the best solution at each of the
1341 iterations.
1342 finite_diff_rel_step : None or array_like, optional
1343 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
1344 use for numerical approximation of the jacobian. The absolute step
1345 size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,
1346 possibly adjusted to fit into the bounds. For ``method='3-point'``
1347 the sign of `h` is ignored. If None (default) then step is selected
1348 automatically.
1349 xrtol : float, default: 0
1350 Relative tolerance for `x`. Terminate successfully if step size is
1351 less than ``xk * xrtol`` where ``xk`` is the current parameter vector.
1352 """
1353 _check_unknown_options(unknown_options)
1354 retall = return_all
1356 x0 = asarray(x0).flatten()
1357 if x0.ndim == 0:
1358 x0.shape = (1,)
1359 if maxiter is None:
1360 maxiter = len(x0) * 200
1362 sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
1363 finite_diff_rel_step=finite_diff_rel_step)
1365 f = sf.fun
1366 myfprime = sf.grad
1368 old_fval = f(x0)
1369 gfk = myfprime(x0)
1371 k = 0
1372 N = len(x0)
1373 I = np.eye(N, dtype=int)
1374 Hk = I
1376 # Sets the initial step guess to dx ~ 1
1377 old_old_fval = old_fval + np.linalg.norm(gfk) / 2
1379 xk = x0
1380 if retall:
1381 allvecs = [x0]
1382 warnflag = 0
1383 gnorm = vecnorm(gfk, ord=norm)
1384 while (gnorm > gtol) and (k < maxiter):
1385 pk = -np.dot(Hk, gfk)
1386 try:
1387 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
1388 _line_search_wolfe12(f, myfprime, xk, pk, gfk,
1389 old_fval, old_old_fval, amin=1e-100, amax=1e100)
1390 except _LineSearchError:
1391 # Line search failed to find a better solution.
1392 warnflag = 2
1393 break
1395 sk = alpha_k * pk
1396 xkp1 = xk + sk
1398 if retall:
1399 allvecs.append(xkp1)
1400 xk = xkp1
1401 if gfkp1 is None:
1402 gfkp1 = myfprime(xkp1)
1404 yk = gfkp1 - gfk
1405 gfk = gfkp1
1406 if callback is not None:
1407 callback(xk)
1408 k += 1
1409 gnorm = vecnorm(gfk, ord=norm)
1410 if (gnorm <= gtol):
1411 break
1413 # See Chapter 5 in P.E. Frandsen, K. Jonasson, H.B. Nielsen,
1414 # O. Tingleff: "Unconstrained Optimization", IMM, DTU. 1999.
1415 # These notes are available here:
1416 # http://www2.imm.dtu.dk/documents/ftp/publlec.html
1417 if (alpha_k*vecnorm(pk) <= xrtol*(xrtol + vecnorm(xk))):
1418 break
1420 if not np.isfinite(old_fval):
1421 # We correctly found +-Inf as optimal value, or something went
1422 # wrong.
1423 warnflag = 2
1424 break
1426 rhok_inv = np.dot(yk, sk)
1427 # this was handled in numeric, let it remaines for more safety
1428 # Cryptic comment above is preserved for posterity. Future reader:
1429 # consider change to condition below proposed in gh-1261/gh-17345.
1430 if rhok_inv == 0.:
1431 rhok = 1000.0
1432 if disp:
1433 print("Divide-by-zero encountered: rhok assumed large")
1434 else:
1435 rhok = 1. / rhok_inv
1437 A1 = I - sk[:, np.newaxis] * yk[np.newaxis, :] * rhok
1438 A2 = I - yk[:, np.newaxis] * sk[np.newaxis, :] * rhok
1439 Hk = np.dot(A1, np.dot(Hk, A2)) + (rhok * sk[:, np.newaxis] *
1440 sk[np.newaxis, :])
1442 fval = old_fval
1444 if warnflag == 2:
1445 msg = _status_message['pr_loss']
1446 elif k >= maxiter:
1447 warnflag = 1
1448 msg = _status_message['maxiter']
1449 elif np.isnan(gnorm) or np.isnan(fval) or np.isnan(xk).any():
1450 warnflag = 3
1451 msg = _status_message['nan']
1452 else:
1453 msg = _status_message['success']
1455 if disp:
1456 print("%s%s" % ("Warning: " if warnflag != 0 else "", msg))
1457 print(" Current function value: %f" % fval)
1458 print(" Iterations: %d" % k)
1459 print(" Function evaluations: %d" % sf.nfev)
1460 print(" Gradient evaluations: %d" % sf.ngev)
1462 result = OptimizeResult(fun=fval, jac=gfk, hess_inv=Hk, nfev=sf.nfev,
1463 njev=sf.ngev, status=warnflag,
1464 success=(warnflag == 0), message=msg, x=xk,
1465 nit=k)
1466 if retall:
1467 result['allvecs'] = allvecs
1468 return result
1471def fmin_cg(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf, epsilon=_epsilon,
1472 maxiter=None, full_output=0, disp=1, retall=0, callback=None):
1473 """
1474 Minimize a function using a nonlinear conjugate gradient algorithm.
1476 Parameters
1477 ----------
1478 f : callable, ``f(x, *args)``
1479 Objective function to be minimized. Here `x` must be a 1-D array of
1480 the variables that are to be changed in the search for a minimum, and
1481 `args` are the other (fixed) parameters of `f`.
1482 x0 : ndarray
1483 A user-supplied initial estimate of `xopt`, the optimal value of `x`.
1484 It must be a 1-D array of values.
1485 fprime : callable, ``fprime(x, *args)``, optional
1486 A function that returns the gradient of `f` at `x`. Here `x` and `args`
1487 are as described above for `f`. The returned value must be a 1-D array.
1488 Defaults to None, in which case the gradient is approximated
1489 numerically (see `epsilon`, below).
1490 args : tuple, optional
1491 Parameter values passed to `f` and `fprime`. Must be supplied whenever
1492 additional fixed parameters are needed to completely specify the
1493 functions `f` and `fprime`.
1494 gtol : float, optional
1495 Stop when the norm of the gradient is less than `gtol`.
1496 norm : float, optional
1497 Order to use for the norm of the gradient
1498 (``-np.Inf`` is min, ``np.Inf`` is max).
1499 epsilon : float or ndarray, optional
1500 Step size(s) to use when `fprime` is approximated numerically. Can be a
1501 scalar or a 1-D array. Defaults to ``sqrt(eps)``, with eps the
1502 floating point machine precision. Usually ``sqrt(eps)`` is about
1503 1.5e-8.
1504 maxiter : int, optional
1505 Maximum number of iterations to perform. Default is ``200 * len(x0)``.
1506 full_output : bool, optional
1507 If True, return `fopt`, `func_calls`, `grad_calls`, and `warnflag` in
1508 addition to `xopt`. See the Returns section below for additional
1509 information on optional return values.
1510 disp : bool, optional
1511 If True, return a convergence message, followed by `xopt`.
1512 retall : bool, optional
1513 If True, add to the returned values the results of each iteration.
1514 callback : callable, optional
1515 An optional user-supplied function, called after each iteration.
1516 Called as ``callback(xk)``, where ``xk`` is the current value of `x0`.
1518 Returns
1519 -------
1520 xopt : ndarray
1521 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1522 fopt : float, optional
1523 Minimum value found, f(xopt). Only returned if `full_output` is True.
1524 func_calls : int, optional
1525 The number of function_calls made. Only returned if `full_output`
1526 is True.
1527 grad_calls : int, optional
1528 The number of gradient calls made. Only returned if `full_output` is
1529 True.
1530 warnflag : int, optional
1531 Integer value with warning status, only returned if `full_output` is
1532 True.
1534 0 : Success.
1536 1 : The maximum number of iterations was exceeded.
1538 2 : Gradient and/or function calls were not changing. May indicate
1539 that precision was lost, i.e., the routine did not converge.
1541 3 : NaN result encountered.
1543 allvecs : list of ndarray, optional
1544 List of arrays, containing the results at each iteration.
1545 Only returned if `retall` is True.
1547 See Also
1548 --------
1549 minimize : common interface to all `scipy.optimize` algorithms for
1550 unconstrained and constrained minimization of multivariate
1551 functions. It provides an alternative way to call
1552 ``fmin_cg``, by specifying ``method='CG'``.
1554 Notes
1555 -----
1556 This conjugate gradient algorithm is based on that of Polak and Ribiere
1557 [1]_.
1559 Conjugate gradient methods tend to work better when:
1561 1. `f` has a unique global minimizing point, and no local minima or
1562 other stationary points,
1563 2. `f` is, at least locally, reasonably well approximated by a
1564 quadratic function of the variables,
1565 3. `f` is continuous and has a continuous gradient,
1566 4. `fprime` is not too large, e.g., has a norm less than 1000,
1567 5. The initial guess, `x0`, is reasonably close to `f` 's global
1568 minimizing point, `xopt`.
1570 References
1571 ----------
1572 .. [1] Wright & Nocedal, "Numerical Optimization", 1999, pp. 120-122.
1574 Examples
1575 --------
1576 Example 1: seek the minimum value of the expression
1577 ``a*u**2 + b*u*v + c*v**2 + d*u + e*v + f`` for given values
1578 of the parameters and an initial guess ``(u, v) = (0, 0)``.
1580 >>> import numpy as np
1581 >>> args = (2, 3, 7, 8, 9, 10) # parameter values
1582 >>> def f(x, *args):
1583 ... u, v = x
1584 ... a, b, c, d, e, f = args
1585 ... return a*u**2 + b*u*v + c*v**2 + d*u + e*v + f
1586 >>> def gradf(x, *args):
1587 ... u, v = x
1588 ... a, b, c, d, e, f = args
1589 ... gu = 2*a*u + b*v + d # u-component of the gradient
1590 ... gv = b*u + 2*c*v + e # v-component of the gradient
1591 ... return np.asarray((gu, gv))
1592 >>> x0 = np.asarray((0, 0)) # Initial guess.
1593 >>> from scipy import optimize
1594 >>> res1 = optimize.fmin_cg(f, x0, fprime=gradf, args=args)
1595 Optimization terminated successfully.
1596 Current function value: 1.617021
1597 Iterations: 4
1598 Function evaluations: 8
1599 Gradient evaluations: 8
1600 >>> res1
1601 array([-1.80851064, -0.25531915])
1603 Example 2: solve the same problem using the `minimize` function.
1604 (This `myopts` dictionary shows all of the available options,
1605 although in practice only non-default values would be needed.
1606 The returned value will be a dictionary.)
1608 >>> opts = {'maxiter' : None, # default value.
1609 ... 'disp' : True, # non-default value.
1610 ... 'gtol' : 1e-5, # default value.
1611 ... 'norm' : np.inf, # default value.
1612 ... 'eps' : 1.4901161193847656e-08} # default value.
1613 >>> res2 = optimize.minimize(f, x0, jac=gradf, args=args,
1614 ... method='CG', options=opts)
1615 Optimization terminated successfully.
1616 Current function value: 1.617021
1617 Iterations: 4
1618 Function evaluations: 8
1619 Gradient evaluations: 8
1620 >>> res2.x # minimum found
1621 array([-1.80851064, -0.25531915])
1623 """
1624 opts = {'gtol': gtol,
1625 'norm': norm,
1626 'eps': epsilon,
1627 'disp': disp,
1628 'maxiter': maxiter,
1629 'return_all': retall}
1631 res = _minimize_cg(f, x0, args, fprime, callback=callback, **opts)
1633 if full_output:
1634 retlist = res['x'], res['fun'], res['nfev'], res['njev'], res['status']
1635 if retall:
1636 retlist += (res['allvecs'], )
1637 return retlist
1638 else:
1639 if retall:
1640 return res['x'], res['allvecs']
1641 else:
1642 return res['x']
1645def _minimize_cg(fun, x0, args=(), jac=None, callback=None,
1646 gtol=1e-5, norm=Inf, eps=_epsilon, maxiter=None,
1647 disp=False, return_all=False, finite_diff_rel_step=None,
1648 **unknown_options):
1649 """
1650 Minimization of scalar function of one or more variables using the
1651 conjugate gradient algorithm.
1653 Options
1654 -------
1655 disp : bool
1656 Set to True to print convergence messages.
1657 maxiter : int
1658 Maximum number of iterations to perform.
1659 gtol : float
1660 Gradient norm must be less than `gtol` before successful
1661 termination.
1662 norm : float
1663 Order of norm (Inf is max, -Inf is min).
1664 eps : float or ndarray
1665 If `jac is None` the absolute step size used for numerical
1666 approximation of the jacobian via forward differences.
1667 return_all : bool, optional
1668 Set to True to return a list of the best solution at each of the
1669 iterations.
1670 finite_diff_rel_step : None or array_like, optional
1671 If `jac in ['2-point', '3-point', 'cs']` the relative step size to
1672 use for numerical approximation of the jacobian. The absolute step
1673 size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``,
1674 possibly adjusted to fit into the bounds. For ``method='3-point'``
1675 the sign of `h` is ignored. If None (default) then step is selected
1676 automatically.
1677 """
1678 _check_unknown_options(unknown_options)
1680 retall = return_all
1682 x0 = asarray(x0).flatten()
1683 if maxiter is None:
1684 maxiter = len(x0) * 200
1686 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
1687 finite_diff_rel_step=finite_diff_rel_step)
1689 f = sf.fun
1690 myfprime = sf.grad
1692 old_fval = f(x0)
1693 gfk = myfprime(x0)
1695 k = 0
1696 xk = x0
1697 # Sets the initial step guess to dx ~ 1
1698 old_old_fval = old_fval + np.linalg.norm(gfk) / 2
1700 if retall:
1701 allvecs = [xk]
1702 warnflag = 0
1703 pk = -gfk
1704 gnorm = vecnorm(gfk, ord=norm)
1706 sigma_3 = 0.01
1708 while (gnorm > gtol) and (k < maxiter):
1709 deltak = np.dot(gfk, gfk)
1711 cached_step = [None]
1713 def polak_ribiere_powell_step(alpha, gfkp1=None):
1714 xkp1 = xk + alpha * pk
1715 if gfkp1 is None:
1716 gfkp1 = myfprime(xkp1)
1717 yk = gfkp1 - gfk
1718 beta_k = max(0, np.dot(yk, gfkp1) / deltak)
1719 pkp1 = -gfkp1 + beta_k * pk
1720 gnorm = vecnorm(gfkp1, ord=norm)
1721 return (alpha, xkp1, pkp1, gfkp1, gnorm)
1723 def descent_condition(alpha, xkp1, fp1, gfkp1):
1724 # Polak-Ribiere+ needs an explicit check of a sufficient
1725 # descent condition, which is not guaranteed by strong Wolfe.
1726 #
1727 # See Gilbert & Nocedal, "Global convergence properties of
1728 # conjugate gradient methods for optimization",
1729 # SIAM J. Optimization 2, 21 (1992).
1730 cached_step[:] = polak_ribiere_powell_step(alpha, gfkp1)
1731 alpha, xk, pk, gfk, gnorm = cached_step
1733 # Accept step if it leads to convergence.
1734 if gnorm <= gtol:
1735 return True
1737 # Accept step if sufficient descent condition applies.
1738 return np.dot(pk, gfk) <= -sigma_3 * np.dot(gfk, gfk)
1740 try:
1741 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
1742 _line_search_wolfe12(f, myfprime, xk, pk, gfk, old_fval,
1743 old_old_fval, c2=0.4, amin=1e-100, amax=1e100,
1744 extra_condition=descent_condition)
1745 except _LineSearchError:
1746 # Line search failed to find a better solution.
1747 warnflag = 2
1748 break
1750 # Reuse already computed results if possible
1751 if alpha_k == cached_step[0]:
1752 alpha_k, xk, pk, gfk, gnorm = cached_step
1753 else:
1754 alpha_k, xk, pk, gfk, gnorm = polak_ribiere_powell_step(alpha_k, gfkp1)
1756 if retall:
1757 allvecs.append(xk)
1758 if callback is not None:
1759 callback(xk)
1760 k += 1
1762 fval = old_fval
1763 if warnflag == 2:
1764 msg = _status_message['pr_loss']
1765 elif k >= maxiter:
1766 warnflag = 1
1767 msg = _status_message['maxiter']
1768 elif np.isnan(gnorm) or np.isnan(fval) or np.isnan(xk).any():
1769 warnflag = 3
1770 msg = _status_message['nan']
1771 else:
1772 msg = _status_message['success']
1774 if disp:
1775 print("%s%s" % ("Warning: " if warnflag != 0 else "", msg))
1776 print(" Current function value: %f" % fval)
1777 print(" Iterations: %d" % k)
1778 print(" Function evaluations: %d" % sf.nfev)
1779 print(" Gradient evaluations: %d" % sf.ngev)
1781 result = OptimizeResult(fun=fval, jac=gfk, nfev=sf.nfev,
1782 njev=sf.ngev, status=warnflag,
1783 success=(warnflag == 0), message=msg, x=xk,
1784 nit=k)
1785 if retall:
1786 result['allvecs'] = allvecs
1787 return result
1790def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
1791 epsilon=_epsilon, maxiter=None, full_output=0, disp=1, retall=0,
1792 callback=None):
1793 """
1794 Unconstrained minimization of a function using the Newton-CG method.
1796 Parameters
1797 ----------
1798 f : callable ``f(x, *args)``
1799 Objective function to be minimized.
1800 x0 : ndarray
1801 Initial guess.
1802 fprime : callable ``f'(x, *args)``
1803 Gradient of f.
1804 fhess_p : callable ``fhess_p(x, p, *args)``, optional
1805 Function which computes the Hessian of f times an
1806 arbitrary vector, p.
1807 fhess : callable ``fhess(x, *args)``, optional
1808 Function to compute the Hessian matrix of f.
1809 args : tuple, optional
1810 Extra arguments passed to f, fprime, fhess_p, and fhess
1811 (the same set of extra arguments is supplied to all of
1812 these functions).
1813 epsilon : float or ndarray, optional
1814 If fhess is approximated, use this value for the step size.
1815 callback : callable, optional
1816 An optional user-supplied function which is called after
1817 each iteration. Called as callback(xk), where xk is the
1818 current parameter vector.
1819 avextol : float, optional
1820 Convergence is assumed when the average relative error in
1821 the minimizer falls below this amount.
1822 maxiter : int, optional
1823 Maximum number of iterations to perform.
1824 full_output : bool, optional
1825 If True, return the optional outputs.
1826 disp : bool, optional
1827 If True, print convergence message.
1828 retall : bool, optional
1829 If True, return a list of results at each iteration.
1831 Returns
1832 -------
1833 xopt : ndarray
1834 Parameters which minimize f, i.e., ``f(xopt) == fopt``.
1835 fopt : float
1836 Value of the function at xopt, i.e., ``fopt = f(xopt)``.
1837 fcalls : int
1838 Number of function calls made.
1839 gcalls : int
1840 Number of gradient calls made.
1841 hcalls : int
1842 Number of Hessian calls made.
1843 warnflag : int
1844 Warnings generated by the algorithm.
1845 1 : Maximum number of iterations exceeded.
1846 2 : Line search failure (precision loss).
1847 3 : NaN result encountered.
1848 allvecs : list
1849 The result at each iteration, if retall is True (see below).
1851 See also
1852 --------
1853 minimize: Interface to minimization algorithms for multivariate
1854 functions. See the 'Newton-CG' `method` in particular.
1856 Notes
1857 -----
1858 Only one of `fhess_p` or `fhess` need to be given. If `fhess`
1859 is provided, then `fhess_p` will be ignored. If neither `fhess`
1860 nor `fhess_p` is provided, then the hessian product will be
1861 approximated using finite differences on `fprime`. `fhess_p`
1862 must compute the hessian times an arbitrary vector. If it is not
1863 given, finite-differences on `fprime` are used to compute
1864 it.
1866 Newton-CG methods are also called truncated Newton methods. This
1867 function differs from scipy.optimize.fmin_tnc because
1869 1. scipy.optimize.fmin_ncg is written purely in Python using NumPy
1870 and scipy while scipy.optimize.fmin_tnc calls a C function.
1871 2. scipy.optimize.fmin_ncg is only for unconstrained minimization
1872 while scipy.optimize.fmin_tnc is for unconstrained minimization
1873 or box constrained minimization. (Box constraints give
1874 lower and upper bounds for each variable separately.)
1876 References
1877 ----------
1878 Wright & Nocedal, 'Numerical Optimization', 1999, p. 140.
1880 """
1881 opts = {'xtol': avextol,
1882 'eps': epsilon,
1883 'maxiter': maxiter,
1884 'disp': disp,
1885 'return_all': retall}
1887 res = _minimize_newtoncg(f, x0, args, fprime, fhess, fhess_p,
1888 callback=callback, **opts)
1890 if full_output:
1891 retlist = (res['x'], res['fun'], res['nfev'], res['njev'],
1892 res['nhev'], res['status'])
1893 if retall:
1894 retlist += (res['allvecs'], )
1895 return retlist
1896 else:
1897 if retall:
1898 return res['x'], res['allvecs']
1899 else:
1900 return res['x']
1903def _minimize_newtoncg(fun, x0, args=(), jac=None, hess=None, hessp=None,
1904 callback=None, xtol=1e-5, eps=_epsilon, maxiter=None,
1905 disp=False, return_all=False,
1906 **unknown_options):
1907 """
1908 Minimization of scalar function of one or more variables using the
1909 Newton-CG algorithm.
1911 Note that the `jac` parameter (Jacobian) is required.
1913 Options
1914 -------
1915 disp : bool
1916 Set to True to print convergence messages.
1917 xtol : float
1918 Average relative error in solution `xopt` acceptable for
1919 convergence.
1920 maxiter : int
1921 Maximum number of iterations to perform.
1922 eps : float or ndarray
1923 If `hessp` is approximated, use this value for the step size.
1924 return_all : bool, optional
1925 Set to True to return a list of the best solution at each of the
1926 iterations.
1927 """
1928 _check_unknown_options(unknown_options)
1929 if jac is None:
1930 raise ValueError('Jacobian is required for Newton-CG method')
1931 fhess_p = hessp
1932 fhess = hess
1933 avextol = xtol
1934 epsilon = eps
1935 retall = return_all
1937 x0 = asarray(x0).flatten()
1938 # TODO: add hessp (callable or FD) to ScalarFunction?
1939 sf = _prepare_scalar_function(
1940 fun, x0, jac, args=args, epsilon=eps, hess=hess
1941 )
1942 f = sf.fun
1943 fprime = sf.grad
1944 _h = sf.hess(x0)
1946 # Logic for hess/hessp
1947 # - If a callable(hess) is provided, then use that
1948 # - If hess is a FD_METHOD, or the output fom hess(x) is a LinearOperator
1949 # then create a hessp function using those.
1950 # - If hess is None but you have callable(hessp) then use the hessp.
1951 # - If hess and hessp are None then approximate hessp using the grad/jac.
1953 if (hess in FD_METHODS or isinstance(_h, LinearOperator)):
1954 fhess = None
1956 def _hessp(x, p, *args):
1957 return sf.hess(x).dot(p)
1959 fhess_p = _hessp
1961 def terminate(warnflag, msg):
1962 if disp:
1963 print(msg)
1964 print(" Current function value: %f" % old_fval)
1965 print(" Iterations: %d" % k)
1966 print(" Function evaluations: %d" % sf.nfev)
1967 print(" Gradient evaluations: %d" % sf.ngev)
1968 print(" Hessian evaluations: %d" % hcalls)
1969 fval = old_fval
1970 result = OptimizeResult(fun=fval, jac=gfk, nfev=sf.nfev,
1971 njev=sf.ngev, nhev=hcalls, status=warnflag,
1972 success=(warnflag == 0), message=msg, x=xk,
1973 nit=k)
1974 if retall:
1975 result['allvecs'] = allvecs
1976 return result
1978 hcalls = 0
1979 if maxiter is None:
1980 maxiter = len(x0)*200
1981 cg_maxiter = 20*len(x0)
1983 xtol = len(x0) * avextol
1984 update = [2 * xtol]
1985 xk = x0
1986 if retall:
1987 allvecs = [xk]
1988 k = 0
1989 gfk = None
1990 old_fval = f(x0)
1991 old_old_fval = None
1992 float64eps = np.finfo(np.float64).eps
1993 while np.add.reduce(np.abs(update)) > xtol:
1994 if k >= maxiter:
1995 msg = "Warning: " + _status_message['maxiter']
1996 return terminate(1, msg)
1997 # Compute a search direction pk by applying the CG method to
1998 # del2 f(xk) p = - grad f(xk) starting from 0.
1999 b = -fprime(xk)
2000 maggrad = np.add.reduce(np.abs(b))
2001 eta = np.min([0.5, np.sqrt(maggrad)])
2002 termcond = eta * maggrad
2003 xsupi = zeros(len(x0), dtype=x0.dtype)
2004 ri = -b
2005 psupi = -ri
2006 i = 0
2007 dri0 = np.dot(ri, ri)
2009 if fhess is not None: # you want to compute hessian once.
2010 A = sf.hess(xk)
2011 hcalls = hcalls + 1
2013 for k2 in range(cg_maxiter):
2014 if np.add.reduce(np.abs(ri)) <= termcond:
2015 break
2016 if fhess is None:
2017 if fhess_p is None:
2018 Ap = approx_fhess_p(xk, psupi, fprime, epsilon)
2019 else:
2020 Ap = fhess_p(xk, psupi, *args)
2021 hcalls = hcalls + 1
2022 else:
2023 if isinstance(A, HessianUpdateStrategy):
2024 # if hess was supplied as a HessianUpdateStrategy
2025 Ap = A.dot(psupi)
2026 else:
2027 Ap = np.dot(A, psupi)
2028 # check curvature
2029 Ap = asarray(Ap).squeeze() # get rid of matrices...
2030 curv = np.dot(psupi, Ap)
2031 if 0 <= curv <= 3 * float64eps:
2032 break
2033 elif curv < 0:
2034 if (i > 0):
2035 break
2036 else:
2037 # fall back to steepest descent direction
2038 xsupi = dri0 / (-curv) * b
2039 break
2040 alphai = dri0 / curv
2041 xsupi = xsupi + alphai * psupi
2042 ri = ri + alphai * Ap
2043 dri1 = np.dot(ri, ri)
2044 betai = dri1 / dri0
2045 psupi = -ri + betai * psupi
2046 i = i + 1
2047 dri0 = dri1 # update np.dot(ri,ri) for next time.
2048 else:
2049 # curvature keeps increasing, bail out
2050 msg = ("Warning: CG iterations didn't converge. The Hessian is not "
2051 "positive definite.")
2052 return terminate(3, msg)
2054 pk = xsupi # search direction is solution to system.
2055 gfk = -b # gradient at xk
2057 try:
2058 alphak, fc, gc, old_fval, old_old_fval, gfkp1 = \
2059 _line_search_wolfe12(f, fprime, xk, pk, gfk,
2060 old_fval, old_old_fval)
2061 except _LineSearchError:
2062 # Line search failed to find a better solution.
2063 msg = "Warning: " + _status_message['pr_loss']
2064 return terminate(2, msg)
2066 update = alphak * pk
2067 xk = xk + update # upcast if necessary
2068 if callback is not None:
2069 callback(xk)
2070 if retall:
2071 allvecs.append(xk)
2072 k += 1
2073 else:
2074 if np.isnan(old_fval) or np.isnan(update).any():
2075 return terminate(3, _status_message['nan'])
2077 msg = _status_message['success']
2078 return terminate(0, msg)
2081def fminbound(func, x1, x2, args=(), xtol=1e-5, maxfun=500,
2082 full_output=0, disp=1):
2083 """Bounded minimization for scalar functions.
2085 Parameters
2086 ----------
2087 func : callable f(x,*args)
2088 Objective function to be minimized (must accept and return scalars).
2089 x1, x2 : float or array scalar
2090 Finite optimization bounds.
2091 args : tuple, optional
2092 Extra arguments passed to function.
2093 xtol : float, optional
2094 The convergence tolerance.
2095 maxfun : int, optional
2096 Maximum number of function evaluations allowed.
2097 full_output : bool, optional
2098 If True, return optional outputs.
2099 disp : int, optional
2100 If non-zero, print messages.
2101 0 : no message printing.
2102 1 : non-convergence notification messages only.
2103 2 : print a message on convergence too.
2104 3 : print iteration results.
2107 Returns
2108 -------
2109 xopt : ndarray
2110 Parameters (over given interval) which minimize the
2111 objective function.
2112 fval : number
2113 The function value evaluated at the minimizer.
2114 ierr : int
2115 An error flag (0 if converged, 1 if maximum number of
2116 function calls reached).
2117 numfunc : int
2118 The number of function calls made.
2120 See also
2121 --------
2122 minimize_scalar: Interface to minimization algorithms for scalar
2123 univariate functions. See the 'Bounded' `method` in particular.
2125 Notes
2126 -----
2127 Finds a local minimizer of the scalar function `func` in the
2128 interval x1 < xopt < x2 using Brent's method. (See `brent`
2129 for auto-bracketing.)
2131 References
2132 ----------
2133 .. [1] Forsythe, G.E., M. A. Malcolm, and C. B. Moler. "Computer Methods
2134 for Mathematical Computations." Prentice-Hall Series in Automatic
2135 Computation 259 (1977).
2136 .. [2] Brent, Richard P. Algorithms for Minimization Without Derivatives.
2137 Courier Corporation, 2013.
2139 Examples
2140 --------
2141 `fminbound` finds the minimizer of the function in the given range.
2142 The following examples illustrate this.
2144 >>> from scipy import optimize
2145 >>> def f(x):
2146 ... return (x-1)**2
2147 >>> minimizer = optimize.fminbound(f, -4, 4)
2148 >>> minimizer
2149 1.0
2150 >>> minimum = f(minimizer)
2151 >>> minimum
2152 0.0
2153 >>> minimizer = optimize.fminbound(f, 3, 4)
2154 >>> minimizer
2155 3.000005960860986
2156 >>> minimum = f(minimizer)
2157 >>> minimum
2158 4.000023843479476
2159 """
2160 options = {'xatol': xtol,
2161 'maxiter': maxfun,
2162 'disp': disp}
2164 res = _minimize_scalar_bounded(func, (x1, x2), args, **options)
2165 if full_output:
2166 return res['x'], res['fun'], res['status'], res['nfev']
2167 else:
2168 return res['x']
2171def _minimize_scalar_bounded(func, bounds, args=(),
2172 xatol=1e-5, maxiter=500, disp=0,
2173 **unknown_options):
2174 """
2175 Options
2176 -------
2177 maxiter : int
2178 Maximum number of iterations to perform.
2179 disp: int, optional
2180 If non-zero, print messages.
2181 0 : no message printing.
2182 1 : non-convergence notification messages only.
2183 2 : print a message on convergence too.
2184 3 : print iteration results.
2185 xatol : float
2186 Absolute error in solution `xopt` acceptable for convergence.
2188 """
2189 _check_unknown_options(unknown_options)
2190 maxfun = maxiter
2191 # Test bounds are of correct form
2192 if len(bounds) != 2:
2193 raise ValueError('bounds must have two elements.')
2194 x1, x2 = bounds
2196 if not (is_finite_scalar(x1) and is_finite_scalar(x2)):
2197 raise ValueError("Optimization bounds must be finite scalars.")
2199 if x1 > x2:
2200 raise ValueError("The lower bound exceeds the upper bound.")
2202 flag = 0
2203 header = ' Func-count x f(x) Procedure'
2204 step = ' initial'
2206 sqrt_eps = sqrt(2.2e-16)
2207 golden_mean = 0.5 * (3.0 - sqrt(5.0))
2208 a, b = x1, x2
2209 fulc = a + golden_mean * (b - a)
2210 nfc, xf = fulc, fulc
2211 rat = e = 0.0
2212 x = xf
2213 fx = func(x, *args)
2214 num = 1
2215 fmin_data = (1, xf, fx)
2216 fu = np.inf
2218 ffulc = fnfc = fx
2219 xm = 0.5 * (a + b)
2220 tol1 = sqrt_eps * np.abs(xf) + xatol / 3.0
2221 tol2 = 2.0 * tol1
2223 if disp > 2:
2224 print(" ")
2225 print(header)
2226 print("%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)))
2228 while (np.abs(xf - xm) > (tol2 - 0.5 * (b - a))):
2229 golden = 1
2230 # Check for parabolic fit
2231 if np.abs(e) > tol1:
2232 golden = 0
2233 r = (xf - nfc) * (fx - ffulc)
2234 q = (xf - fulc) * (fx - fnfc)
2235 p = (xf - fulc) * q - (xf - nfc) * r
2236 q = 2.0 * (q - r)
2237 if q > 0.0:
2238 p = -p
2239 q = np.abs(q)
2240 r = e
2241 e = rat
2243 # Check for acceptability of parabola
2244 if ((np.abs(p) < np.abs(0.5*q*r)) and (p > q*(a - xf)) and
2245 (p < q * (b - xf))):
2246 rat = (p + 0.0) / q
2247 x = xf + rat
2248 step = ' parabolic'
2250 if ((x - a) < tol2) or ((b - x) < tol2):
2251 si = np.sign(xm - xf) + ((xm - xf) == 0)
2252 rat = tol1 * si
2253 else: # do a golden-section step
2254 golden = 1
2256 if golden: # do a golden-section step
2257 if xf >= xm:
2258 e = a - xf
2259 else:
2260 e = b - xf
2261 rat = golden_mean*e
2262 step = ' golden'
2264 si = np.sign(rat) + (rat == 0)
2265 x = xf + si * np.maximum(np.abs(rat), tol1)
2266 fu = func(x, *args)
2267 num += 1
2268 fmin_data = (num, x, fu)
2269 if disp > 2:
2270 print("%5.0f %12.6g %12.6g %s" % (fmin_data + (step,)))
2272 if fu <= fx:
2273 if x >= xf:
2274 a = xf
2275 else:
2276 b = xf
2277 fulc, ffulc = nfc, fnfc
2278 nfc, fnfc = xf, fx
2279 xf, fx = x, fu
2280 else:
2281 if x < xf:
2282 a = x
2283 else:
2284 b = x
2285 if (fu <= fnfc) or (nfc == xf):
2286 fulc, ffulc = nfc, fnfc
2287 nfc, fnfc = x, fu
2288 elif (fu <= ffulc) or (fulc == xf) or (fulc == nfc):
2289 fulc, ffulc = x, fu
2291 xm = 0.5 * (a + b)
2292 tol1 = sqrt_eps * np.abs(xf) + xatol / 3.0
2293 tol2 = 2.0 * tol1
2295 if num >= maxfun:
2296 flag = 1
2297 break
2299 if np.isnan(xf) or np.isnan(fx) or np.isnan(fu):
2300 flag = 2
2302 fval = fx
2303 if disp > 0:
2304 _endprint(x, flag, fval, maxfun, xatol, disp)
2306 result = OptimizeResult(fun=fval, status=flag, success=(flag == 0),
2307 message={0: 'Solution found.',
2308 1: 'Maximum number of function calls '
2309 'reached.',
2310 2: _status_message['nan']}.get(flag, ''),
2311 x=xf, nfev=num, nit=num)
2313 return result
2316class Brent:
2317 #need to rethink design of __init__
2318 def __init__(self, func, args=(), tol=1.48e-8, maxiter=500,
2319 full_output=0, disp=0):
2320 self.func = func
2321 self.args = args
2322 self.tol = tol
2323 self.maxiter = maxiter
2324 self._mintol = 1.0e-11
2325 self._cg = 0.3819660
2326 self.xmin = None
2327 self.fval = None
2328 self.iter = 0
2329 self.funcalls = 0
2330 self.disp = disp
2332 # need to rethink design of set_bracket (new options, etc.)
2333 def set_bracket(self, brack=None):
2334 self.brack = brack
2336 def get_bracket_info(self):
2337 #set up
2338 func = self.func
2339 args = self.args
2340 brack = self.brack
2341 ### BEGIN core bracket_info code ###
2342 ### carefully DOCUMENT any CHANGES in core ##
2343 if brack is None:
2344 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
2345 elif len(brack) == 2:
2346 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
2347 xb=brack[1], args=args)
2348 elif len(brack) == 3:
2349 xa, xb, xc = brack
2350 if (xa > xc): # swap so xa < xc can be assumed
2351 xc, xa = xa, xc
2352 if not ((xa < xb) and (xb < xc)):
2353 raise ValueError(
2354 "Bracketing values (xa, xb, xc) do not"
2355 " fulfill this requirement: (xa < xb) and (xb < xc)"
2356 )
2357 fa = func(*((xa,) + args))
2358 fb = func(*((xb,) + args))
2359 fc = func(*((xc,) + args))
2360 if not ((fb < fa) and (fb < fc)):
2361 raise ValueError(
2362 "Bracketing values (xa, xb, xc) do not fulfill"
2363 " this requirement: (f(xb) < f(xa)) and (f(xb) < f(xc))"
2364 )
2366 funcalls = 3
2367 else:
2368 raise ValueError("Bracketing interval must be "
2369 "length 2 or 3 sequence.")
2370 ### END core bracket_info code ###
2372 return xa, xb, xc, fa, fb, fc, funcalls
2374 def optimize(self):
2375 # set up for optimization
2376 func = self.func
2377 xa, xb, xc, fa, fb, fc, funcalls = self.get_bracket_info()
2378 _mintol = self._mintol
2379 _cg = self._cg
2380 #################################
2381 #BEGIN CORE ALGORITHM
2382 #################################
2383 x = w = v = xb
2384 fw = fv = fx = fb
2385 if (xa < xc):
2386 a = xa
2387 b = xc
2388 else:
2389 a = xc
2390 b = xa
2391 deltax = 0.0
2392 iter = 0
2394 if self.disp > 2:
2395 print(" ")
2396 print(f"{'Func-count':^12} {'x':^12} {'f(x)': ^12}")
2397 print(f"{funcalls:^12g} {x:^12.6g} {fx:^12.6g}")
2399 while (iter < self.maxiter):
2400 tol1 = self.tol * np.abs(x) + _mintol
2401 tol2 = 2.0 * tol1
2402 xmid = 0.5 * (a + b)
2403 # check for convergence
2404 if np.abs(x - xmid) < (tol2 - 0.5 * (b - a)):
2405 break
2406 # XXX In the first iteration, rat is only bound in the true case
2407 # of this conditional. This used to cause an UnboundLocalError
2408 # (gh-4140). It should be set before the if (but to what?).
2409 if (np.abs(deltax) <= tol1):
2410 if (x >= xmid):
2411 deltax = a - x # do a golden section step
2412 else:
2413 deltax = b - x
2414 rat = _cg * deltax
2415 else: # do a parabolic step
2416 tmp1 = (x - w) * (fx - fv)
2417 tmp2 = (x - v) * (fx - fw)
2418 p = (x - v) * tmp2 - (x - w) * tmp1
2419 tmp2 = 2.0 * (tmp2 - tmp1)
2420 if (tmp2 > 0.0):
2421 p = -p
2422 tmp2 = np.abs(tmp2)
2423 dx_temp = deltax
2424 deltax = rat
2425 # check parabolic fit
2426 if ((p > tmp2 * (a - x)) and (p < tmp2 * (b - x)) and
2427 (np.abs(p) < np.abs(0.5 * tmp2 * dx_temp))):
2428 rat = p * 1.0 / tmp2 # if parabolic step is useful.
2429 u = x + rat
2430 if ((u - a) < tol2 or (b - u) < tol2):
2431 if xmid - x >= 0:
2432 rat = tol1
2433 else:
2434 rat = -tol1
2435 else:
2436 if (x >= xmid):
2437 deltax = a - x # if it's not do a golden section step
2438 else:
2439 deltax = b - x
2440 rat = _cg * deltax
2442 if (np.abs(rat) < tol1): # update by at least tol1
2443 if rat >= 0:
2444 u = x + tol1
2445 else:
2446 u = x - tol1
2447 else:
2448 u = x + rat
2449 fu = func(*((u,) + self.args)) # calculate new output value
2450 funcalls += 1
2452 if (fu > fx): # if it's bigger than current
2453 if (u < x):
2454 a = u
2455 else:
2456 b = u
2457 if (fu <= fw) or (w == x):
2458 v = w
2459 w = u
2460 fv = fw
2461 fw = fu
2462 elif (fu <= fv) or (v == x) or (v == w):
2463 v = u
2464 fv = fu
2465 else:
2466 if (u >= x):
2467 a = x
2468 else:
2469 b = x
2470 v = w
2471 w = x
2472 x = u
2473 fv = fw
2474 fw = fx
2475 fx = fu
2477 if self.disp > 2:
2478 print(f"{funcalls:^12g} {x:^12.6g} {fx:^12.6g}")
2480 iter += 1
2481 #################################
2482 #END CORE ALGORITHM
2483 #################################
2485 self.xmin = x
2486 self.fval = fx
2487 self.iter = iter
2488 self.funcalls = funcalls
2490 def get_result(self, full_output=False):
2491 if full_output:
2492 return self.xmin, self.fval, self.iter, self.funcalls
2493 else:
2494 return self.xmin
2497def brent(func, args=(), brack=None, tol=1.48e-8, full_output=0, maxiter=500):
2498 """
2499 Given a function of one variable and a possible bracket, return
2500 the local minimum of the function isolated to a fractional precision
2501 of tol.
2503 Parameters
2504 ----------
2505 func : callable f(x,*args)
2506 Objective function.
2507 args : tuple, optional
2508 Additional arguments (if present).
2509 brack : tuple, optional
2510 Either a triple (xa,xb,xc) where xa<xb<xc and func(xb) <
2511 func(xa), func(xc) or a pair (xa,xb) which are used as a
2512 starting interval for a downhill bracket search (see
2513 `bracket`). Providing the pair (xa,xb) does not always mean
2514 the obtained solution will satisfy xa<=x<=xb.
2515 tol : float, optional
2516 Relative error in solution `xopt` acceptable for convergence.
2517 full_output : bool, optional
2518 If True, return all output args (xmin, fval, iter,
2519 funcalls).
2520 maxiter : int, optional
2521 Maximum number of iterations in solution.
2523 Returns
2524 -------
2525 xmin : ndarray
2526 Optimum point.
2527 fval : float
2528 Optimum value.
2529 iter : int
2530 Number of iterations.
2531 funcalls : int
2532 Number of objective function evaluations made.
2534 See also
2535 --------
2536 minimize_scalar: Interface to minimization algorithms for scalar
2537 univariate functions. See the 'Brent' `method` in particular.
2539 Notes
2540 -----
2541 Uses inverse parabolic interpolation when possible to speed up
2542 convergence of golden section method.
2544 Does not ensure that the minimum lies in the range specified by
2545 `brack`. See `fminbound`.
2547 Examples
2548 --------
2549 We illustrate the behaviour of the function when `brack` is of
2550 size 2 and 3 respectively. In the case where `brack` is of the
2551 form (xa,xb), we can see for the given values, the output need
2552 not necessarily lie in the range (xa,xb).
2554 >>> def f(x):
2555 ... return x**2
2557 >>> from scipy import optimize
2559 >>> minimum = optimize.brent(f,brack=(1,2))
2560 >>> minimum
2561 0.0
2562 >>> minimum = optimize.brent(f,brack=(-1,0.5,2))
2563 >>> minimum
2564 -2.7755575615628914e-17
2566 """
2567 options = {'xtol': tol,
2568 'maxiter': maxiter}
2569 res = _minimize_scalar_brent(func, brack, args, **options)
2570 if full_output:
2571 return res['x'], res['fun'], res['nit'], res['nfev']
2572 else:
2573 return res['x']
2576def _minimize_scalar_brent(func, brack=None, args=(), xtol=1.48e-8,
2577 maxiter=500, disp=0,
2578 **unknown_options):
2579 """
2580 Options
2581 -------
2582 maxiter : int
2583 Maximum number of iterations to perform.
2584 xtol : float
2585 Relative error in solution `xopt` acceptable for convergence.
2586 disp: int, optional
2587 If non-zero, print messages.
2588 0 : no message printing.
2589 1 : non-convergence notification messages only.
2590 2 : print a message on convergence too.
2591 3 : print iteration results.
2592 Notes
2593 -----
2594 Uses inverse parabolic interpolation when possible to speed up
2595 convergence of golden section method.
2597 """
2598 _check_unknown_options(unknown_options)
2599 tol = xtol
2600 if tol < 0:
2601 raise ValueError('tolerance should be >= 0, got %r' % tol)
2603 brent = Brent(func=func, args=args, tol=tol,
2604 full_output=True, maxiter=maxiter, disp=disp)
2605 brent.set_bracket(brack)
2606 brent.optimize()
2607 x, fval, nit, nfev = brent.get_result(full_output=True)
2609 success = nit < maxiter and not (np.isnan(x) or np.isnan(fval))
2611 if success:
2612 message = ("\nOptimization terminated successfully;\n"
2613 "The returned value satisfies the termination criteria\n"
2614 f"(using xtol = {xtol} )")
2615 else:
2616 if nit >= maxiter:
2617 message = "\nMaximum number of iterations exceeded"
2618 if np.isnan(x) or np.isnan(fval):
2619 message = f"{_status_message['nan']}"
2621 if disp:
2622 print(message)
2624 return OptimizeResult(fun=fval, x=x, nit=nit, nfev=nfev,
2625 success=success, message=message)
2628def golden(func, args=(), brack=None, tol=_epsilon,
2629 full_output=0, maxiter=5000):
2630 """
2631 Return the minimum of a function of one variable using golden section
2632 method.
2634 Given a function of one variable and a possible bracketing interval,
2635 return the minimum of the function isolated to a fractional precision of
2636 tol.
2638 Parameters
2639 ----------
2640 func : callable func(x,*args)
2641 Objective function to minimize.
2642 args : tuple, optional
2643 Additional arguments (if present), passed to func.
2644 brack : tuple, optional
2645 Triple (a,b,c), where (a<b<c) and func(b) <
2646 func(a),func(c). If bracket consists of two numbers (a,
2647 c), then they are assumed to be a starting interval for a
2648 downhill bracket search (see `bracket`); it doesn't always
2649 mean that obtained solution will satisfy a<=x<=c.
2650 tol : float, optional
2651 x tolerance stop criterion
2652 full_output : bool, optional
2653 If True, return optional outputs.
2654 maxiter : int
2655 Maximum number of iterations to perform.
2657 See also
2658 --------
2659 minimize_scalar: Interface to minimization algorithms for scalar
2660 univariate functions. See the 'Golden' `method` in particular.
2662 Notes
2663 -----
2664 Uses analog of bisection method to decrease the bracketed
2665 interval.
2667 Examples
2668 --------
2669 We illustrate the behaviour of the function when `brack` is of
2670 size 2 and 3, respectively. In the case where `brack` is of the
2671 form (xa,xb), we can see for the given values, the output need
2672 not necessarily lie in the range ``(xa, xb)``.
2674 >>> def f(x):
2675 ... return x**2
2677 >>> from scipy import optimize
2679 >>> minimum = optimize.golden(f, brack=(1, 2))
2680 >>> minimum
2681 1.5717277788484873e-162
2682 >>> minimum = optimize.golden(f, brack=(-1, 0.5, 2))
2683 >>> minimum
2684 -1.5717277788484873e-162
2686 """
2687 options = {'xtol': tol, 'maxiter': maxiter}
2688 res = _minimize_scalar_golden(func, brack, args, **options)
2689 if full_output:
2690 return res['x'], res['fun'], res['nfev']
2691 else:
2692 return res['x']
2695def _minimize_scalar_golden(func, brack=None, args=(),
2696 xtol=_epsilon, maxiter=5000, disp=0,
2697 **unknown_options):
2698 """
2699 Options
2700 -------
2701 xtol : float
2702 Relative error in solution `xopt` acceptable for convergence.
2703 maxiter : int
2704 Maximum number of iterations to perform.
2705 disp: int, optional
2706 If non-zero, print messages.
2707 0 : no message printing.
2708 1 : non-convergence notification messages only.
2709 2 : print a message on convergence too.
2710 3 : print iteration results.
2711 """
2712 _check_unknown_options(unknown_options)
2713 tol = xtol
2714 if brack is None:
2715 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, args=args)
2716 elif len(brack) == 2:
2717 xa, xb, xc, fa, fb, fc, funcalls = bracket(func, xa=brack[0],
2718 xb=brack[1], args=args)
2719 elif len(brack) == 3:
2720 xa, xb, xc = brack
2721 if (xa > xc): # swap so xa < xc can be assumed
2722 xc, xa = xa, xc
2723 if not ((xa < xb) and (xb < xc)):
2724 raise ValueError(
2725 "Bracketing values (xa, xb, xc) do not"
2726 " fulfill this requirement: (xa < xb) and (xb < xc)"
2727 )
2728 fa = func(*((xa,) + args))
2729 fb = func(*((xb,) + args))
2730 fc = func(*((xc,) + args))
2731 if not ((fb < fa) and (fb < fc)):
2732 raise ValueError(
2733 "Bracketing values (xa, xb, xc) do not fulfill"
2734 " this requirement: (f(xb) < f(xa)) and (f(xb) < f(xc))"
2735 )
2736 funcalls = 3
2737 else:
2738 raise ValueError("Bracketing interval must be length 2 or 3 sequence.")
2740 _gR = 0.61803399 # golden ratio conjugate: 2.0/(1.0+sqrt(5.0))
2741 _gC = 1.0 - _gR
2742 x3 = xc
2743 x0 = xa
2744 if (np.abs(xc - xb) > np.abs(xb - xa)):
2745 x1 = xb
2746 x2 = xb + _gC * (xc - xb)
2747 else:
2748 x2 = xb
2749 x1 = xb - _gC * (xb - xa)
2750 f1 = func(*((x1,) + args))
2751 f2 = func(*((x2,) + args))
2752 funcalls += 2
2753 nit = 0
2755 if disp > 2:
2756 print(" ")
2757 print(f"{'Func-count':^12} {'x':^12} {'f(x)': ^12}")
2759 for i in range(maxiter):
2760 if np.abs(x3 - x0) <= tol * (np.abs(x1) + np.abs(x2)):
2761 break
2762 if (f2 < f1):
2763 x0 = x1
2764 x1 = x2
2765 x2 = _gR * x1 + _gC * x3
2766 f1 = f2
2767 f2 = func(*((x2,) + args))
2768 else:
2769 x3 = x2
2770 x2 = x1
2771 x1 = _gR * x2 + _gC * x0
2772 f2 = f1
2773 f1 = func(*((x1,) + args))
2774 funcalls += 1
2775 if disp > 2:
2776 if (f1 < f2):
2777 xmin, fval = x1, f1
2778 else:
2779 xmin, fval = x2, f2
2780 print(f"{funcalls:^12g} {xmin:^12.6g} {fval:^12.6g}")
2782 nit += 1
2783 # end of iteration loop
2785 if (f1 < f2):
2786 xmin = x1
2787 fval = f1
2788 else:
2789 xmin = x2
2790 fval = f2
2792 success = nit < maxiter and not (np.isnan(fval) or np.isnan(xmin))
2794 if success:
2795 message = ("\nOptimization terminated successfully;\n"
2796 "The returned value satisfies the termination criteria\n"
2797 f"(using xtol = {xtol} )")
2798 else:
2799 if nit >= maxiter:
2800 message = "\nMaximum number of iterations exceeded"
2801 if np.isnan(xmin) or np.isnan(fval):
2802 message = f"{_status_message['nan']}"
2804 if disp:
2805 print(message)
2807 return OptimizeResult(fun=fval, nfev=funcalls, x=xmin, nit=nit,
2808 success=success, message=message)
2811def bracket(func, xa=0.0, xb=1.0, args=(), grow_limit=110.0, maxiter=1000):
2812 """
2813 Bracket the minimum of the function.
2815 Given a function and distinct initial points, search in the
2816 downhill direction (as defined by the initial points) and return
2817 new points xa, xb, xc that bracket the minimum of the function
2818 f(xa) > f(xb) < f(xc). It doesn't always mean that obtained
2819 solution will satisfy xa<=x<=xb.
2821 Parameters
2822 ----------
2823 func : callable f(x,*args)
2824 Objective function to minimize.
2825 xa, xb : float, optional
2826 Bracketing interval. Defaults `xa` to 0.0, and `xb` to 1.0.
2827 args : tuple, optional
2828 Additional arguments (if present), passed to `func`.
2829 grow_limit : float, optional
2830 Maximum grow limit. Defaults to 110.0
2831 maxiter : int, optional
2832 Maximum number of iterations to perform. Defaults to 1000.
2834 Returns
2835 -------
2836 xa, xb, xc : float
2837 Bracket.
2838 fa, fb, fc : float
2839 Objective function values in bracket.
2840 funcalls : int
2841 Number of function evaluations made.
2843 Examples
2844 --------
2845 This function can find a downward convex region of a function:
2847 >>> import numpy as np
2848 >>> import matplotlib.pyplot as plt
2849 >>> from scipy.optimize import bracket
2850 >>> def f(x):
2851 ... return 10*x**2 + 3*x + 5
2852 >>> x = np.linspace(-2, 2)
2853 >>> y = f(x)
2854 >>> init_xa, init_xb = 0, 1
2855 >>> xa, xb, xc, fa, fb, fc, funcalls = bracket(f, xa=init_xa, xb=init_xb)
2856 >>> plt.axvline(x=init_xa, color="k", linestyle="--")
2857 >>> plt.axvline(x=init_xb, color="k", linestyle="--")
2858 >>> plt.plot(x, y, "-k")
2859 >>> plt.plot(xa, fa, "bx")
2860 >>> plt.plot(xb, fb, "rx")
2861 >>> plt.plot(xc, fc, "bx")
2862 >>> plt.show()
2864 """
2865 _gold = 1.618034 # golden ratio: (1.0+sqrt(5.0))/2.0
2866 _verysmall_num = 1e-21
2867 fa = func(*(xa,) + args)
2868 fb = func(*(xb,) + args)
2869 if (fa < fb): # Switch so fa > fb
2870 xa, xb = xb, xa
2871 fa, fb = fb, fa
2872 xc = xb + _gold * (xb - xa)
2873 fc = func(*((xc,) + args))
2874 funcalls = 3
2875 iter = 0
2876 while (fc < fb):
2877 tmp1 = (xb - xa) * (fb - fc)
2878 tmp2 = (xb - xc) * (fb - fa)
2879 val = tmp2 - tmp1
2880 if np.abs(val) < _verysmall_num:
2881 denom = 2.0 * _verysmall_num
2882 else:
2883 denom = 2.0 * val
2884 w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom
2885 wlim = xb + grow_limit * (xc - xb)
2886 if iter > maxiter:
2887 raise RuntimeError("Too many iterations.")
2888 iter += 1
2889 if (w - xc) * (xb - w) > 0.0:
2890 fw = func(*((w,) + args))
2891 funcalls += 1
2892 if (fw < fc):
2893 xa = xb
2894 xb = w
2895 fa = fb
2896 fb = fw
2897 return xa, xb, xc, fa, fb, fc, funcalls
2898 elif (fw > fb):
2899 xc = w
2900 fc = fw
2901 return xa, xb, xc, fa, fb, fc, funcalls
2902 w = xc + _gold * (xc - xb)
2903 fw = func(*((w,) + args))
2904 funcalls += 1
2905 elif (w - wlim)*(wlim - xc) >= 0.0:
2906 w = wlim
2907 fw = func(*((w,) + args))
2908 funcalls += 1
2909 elif (w - wlim)*(xc - w) > 0.0:
2910 fw = func(*((w,) + args))
2911 funcalls += 1
2912 if (fw < fc):
2913 xb = xc
2914 xc = w
2915 w = xc + _gold * (xc - xb)
2916 fb = fc
2917 fc = fw
2918 fw = func(*((w,) + args))
2919 funcalls += 1
2920 else:
2921 w = xc + _gold * (xc - xb)
2922 fw = func(*((w,) + args))
2923 funcalls += 1
2924 xa = xb
2925 xb = xc
2926 xc = w
2927 fa = fb
2928 fb = fc
2929 fc = fw
2930 return xa, xb, xc, fa, fb, fc, funcalls
2933def _line_for_search(x0, alpha, lower_bound, upper_bound):
2934 """
2935 Given a parameter vector ``x0`` with length ``n`` and a direction
2936 vector ``alpha`` with length ``n``, and lower and upper bounds on
2937 each of the ``n`` parameters, what are the bounds on a scalar
2938 ``l`` such that ``lower_bound <= x0 + alpha * l <= upper_bound``.
2941 Parameters
2942 ----------
2943 x0 : np.array.
2944 The vector representing the current location.
2945 Note ``np.shape(x0) == (n,)``.
2946 alpha : np.array.
2947 The vector representing the direction.
2948 Note ``np.shape(alpha) == (n,)``.
2949 lower_bound : np.array.
2950 The lower bounds for each parameter in ``x0``. If the ``i``th
2951 parameter in ``x0`` is unbounded below, then ``lower_bound[i]``
2952 should be ``-np.inf``.
2953 Note ``np.shape(lower_bound) == (n,)``.
2954 upper_bound : np.array.
2955 The upper bounds for each parameter in ``x0``. If the ``i``th
2956 parameter in ``x0`` is unbounded above, then ``upper_bound[i]``
2957 should be ``np.inf``.
2958 Note ``np.shape(upper_bound) == (n,)``.
2960 Returns
2961 -------
2962 res : tuple ``(lmin, lmax)``
2963 The bounds for ``l`` such that
2964 ``lower_bound[i] <= x0[i] + alpha[i] * l <= upper_bound[i]``
2965 for all ``i``.
2967 """
2968 # get nonzero indices of alpha so we don't get any zero division errors.
2969 # alpha will not be all zero, since it is called from _linesearch_powell
2970 # where we have a check for this.
2971 nonzero, = alpha.nonzero()
2972 lower_bound, upper_bound = lower_bound[nonzero], upper_bound[nonzero]
2973 x0, alpha = x0[nonzero], alpha[nonzero]
2974 low = (lower_bound - x0) / alpha
2975 high = (upper_bound - x0) / alpha
2977 # positive and negative indices
2978 pos = alpha > 0
2980 lmin_pos = np.where(pos, low, 0)
2981 lmin_neg = np.where(pos, 0, high)
2982 lmax_pos = np.where(pos, high, 0)
2983 lmax_neg = np.where(pos, 0, low)
2985 lmin = np.max(lmin_pos + lmin_neg)
2986 lmax = np.min(lmax_pos + lmax_neg)
2988 # if x0 is outside the bounds, then it is possible that there is
2989 # no way to get back in the bounds for the parameters being updated
2990 # with the current direction alpha.
2991 # when this happens, lmax < lmin.
2992 # If this is the case, then we can just return (0, 0)
2993 return (lmin, lmax) if lmax >= lmin else (0, 0)
2996def _linesearch_powell(func, p, xi, tol=1e-3,
2997 lower_bound=None, upper_bound=None, fval=None):
2998 """Line-search algorithm using fminbound.
3000 Find the minimium of the function ``func(x0 + alpha*direc)``.
3002 lower_bound : np.array.
3003 The lower bounds for each parameter in ``x0``. If the ``i``th
3004 parameter in ``x0`` is unbounded below, then ``lower_bound[i]``
3005 should be ``-np.inf``.
3006 Note ``np.shape(lower_bound) == (n,)``.
3007 upper_bound : np.array.
3008 The upper bounds for each parameter in ``x0``. If the ``i``th
3009 parameter in ``x0`` is unbounded above, then ``upper_bound[i]``
3010 should be ``np.inf``.
3011 Note ``np.shape(upper_bound) == (n,)``.
3012 fval : number.
3013 ``fval`` is equal to ``func(p)``, the idea is just to avoid
3014 recomputing it so we can limit the ``fevals``.
3016 """
3017 def myfunc(alpha):
3018 return func(p + alpha*xi)
3020 # if xi is zero, then don't optimize
3021 if not np.any(xi):
3022 return ((fval, p, xi) if fval is not None else (func(p), p, xi))
3023 elif lower_bound is None and upper_bound is None:
3024 # non-bounded minimization
3025 alpha_min, fret, _, _ = brent(myfunc, full_output=1, tol=tol)
3026 xi = alpha_min * xi
3027 return squeeze(fret), p + xi, xi
3028 else:
3029 bound = _line_for_search(p, xi, lower_bound, upper_bound)
3030 if np.isneginf(bound[0]) and np.isposinf(bound[1]):
3031 # equivalent to unbounded
3032 return _linesearch_powell(func, p, xi, fval=fval, tol=tol)
3033 elif not np.isneginf(bound[0]) and not np.isposinf(bound[1]):
3034 # we can use a bounded scalar minimization
3035 res = _minimize_scalar_bounded(myfunc, bound, xatol=tol / 100)
3036 xi = res.x * xi
3037 return squeeze(res.fun), p + xi, xi
3038 else:
3039 # only bounded on one side. use the tangent function to convert
3040 # the infinity bound to a finite bound. The new bounded region
3041 # is a subregion of the region bounded by -np.pi/2 and np.pi/2.
3042 bound = np.arctan(bound[0]), np.arctan(bound[1])
3043 res = _minimize_scalar_bounded(
3044 lambda x: myfunc(np.tan(x)),
3045 bound,
3046 xatol=tol / 100)
3047 xi = np.tan(res.x) * xi
3048 return squeeze(res.fun), p + xi, xi
3051def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None,
3052 maxfun=None, full_output=0, disp=1, retall=0, callback=None,
3053 direc=None):
3054 """
3055 Minimize a function using modified Powell's method.
3057 This method only uses function values, not derivatives.
3059 Parameters
3060 ----------
3061 func : callable f(x,*args)
3062 Objective function to be minimized.
3063 x0 : ndarray
3064 Initial guess.
3065 args : tuple, optional
3066 Extra arguments passed to func.
3067 xtol : float, optional
3068 Line-search error tolerance.
3069 ftol : float, optional
3070 Relative error in ``func(xopt)`` acceptable for convergence.
3071 maxiter : int, optional
3072 Maximum number of iterations to perform.
3073 maxfun : int, optional
3074 Maximum number of function evaluations to make.
3075 full_output : bool, optional
3076 If True, ``fopt``, ``xi``, ``direc``, ``iter``, ``funcalls``, and
3077 ``warnflag`` are returned.
3078 disp : bool, optional
3079 If True, print convergence messages.
3080 retall : bool, optional
3081 If True, return a list of the solution at each iteration.
3082 callback : callable, optional
3083 An optional user-supplied function, called after each
3084 iteration. Called as ``callback(xk)``, where ``xk`` is the
3085 current parameter vector.
3086 direc : ndarray, optional
3087 Initial fitting step and parameter order set as an (N, N) array, where N
3088 is the number of fitting parameters in `x0`. Defaults to step size 1.0
3089 fitting all parameters simultaneously (``np.eye((N, N))``). To
3090 prevent initial consideration of values in a step or to change initial
3091 step size, set to 0 or desired step size in the Jth position in the Mth
3092 block, where J is the position in `x0` and M is the desired evaluation
3093 step, with steps being evaluated in index order. Step size and ordering
3094 will change freely as minimization proceeds.
3096 Returns
3097 -------
3098 xopt : ndarray
3099 Parameter which minimizes `func`.
3100 fopt : number
3101 Value of function at minimum: ``fopt = func(xopt)``.
3102 direc : ndarray
3103 Current direction set.
3104 iter : int
3105 Number of iterations.
3106 funcalls : int
3107 Number of function calls made.
3108 warnflag : int
3109 Integer warning flag:
3110 1 : Maximum number of function evaluations.
3111 2 : Maximum number of iterations.
3112 3 : NaN result encountered.
3113 4 : The result is out of the provided bounds.
3114 allvecs : list
3115 List of solutions at each iteration.
3117 See also
3118 --------
3119 minimize: Interface to unconstrained minimization algorithms for
3120 multivariate functions. See the 'Powell' method in particular.
3122 Notes
3123 -----
3124 Uses a modification of Powell's method to find the minimum of
3125 a function of N variables. Powell's method is a conjugate
3126 direction method.
3128 The algorithm has two loops. The outer loop merely iterates over the inner
3129 loop. The inner loop minimizes over each current direction in the direction
3130 set. At the end of the inner loop, if certain conditions are met, the
3131 direction that gave the largest decrease is dropped and replaced with the
3132 difference between the current estimated x and the estimated x from the
3133 beginning of the inner-loop.
3135 The technical conditions for replacing the direction of greatest
3136 increase amount to checking that
3138 1. No further gain can be made along the direction of greatest increase
3139 from that iteration.
3140 2. The direction of greatest increase accounted for a large sufficient
3141 fraction of the decrease in the function value from that iteration of
3142 the inner loop.
3144 References
3145 ----------
3146 Powell M.J.D. (1964) An efficient method for finding the minimum of a
3147 function of several variables without calculating derivatives,
3148 Computer Journal, 7 (2):155-162.
3150 Press W., Teukolsky S.A., Vetterling W.T., and Flannery B.P.:
3151 Numerical Recipes (any edition), Cambridge University Press
3153 Examples
3154 --------
3155 >>> def f(x):
3156 ... return x**2
3158 >>> from scipy import optimize
3160 >>> minimum = optimize.fmin_powell(f, -1)
3161 Optimization terminated successfully.
3162 Current function value: 0.000000
3163 Iterations: 2
3164 Function evaluations: 16
3165 >>> minimum
3166 array(0.0)
3168 """
3169 opts = {'xtol': xtol,
3170 'ftol': ftol,
3171 'maxiter': maxiter,
3172 'maxfev': maxfun,
3173 'disp': disp,
3174 'direc': direc,
3175 'return_all': retall}
3177 res = _minimize_powell(func, x0, args, callback=callback, **opts)
3179 if full_output:
3180 retlist = (res['x'], res['fun'], res['direc'], res['nit'],
3181 res['nfev'], res['status'])
3182 if retall:
3183 retlist += (res['allvecs'], )
3184 return retlist
3185 else:
3186 if retall:
3187 return res['x'], res['allvecs']
3188 else:
3189 return res['x']
3192def _minimize_powell(func, x0, args=(), callback=None, bounds=None,
3193 xtol=1e-4, ftol=1e-4, maxiter=None, maxfev=None,
3194 disp=False, direc=None, return_all=False,
3195 **unknown_options):
3196 """
3197 Minimization of scalar function of one or more variables using the
3198 modified Powell algorithm.
3200 Parameters
3201 ----------
3202 fun : callable
3203 The objective function to be minimized.
3205 ``fun(x, *args) -> float``
3207 where ``x`` is a 1-D array with shape (n,) and ``args``
3208 is a tuple of the fixed parameters needed to completely
3209 specify the function.
3210 x0 : ndarray, shape (n,)
3211 Initial guess. Array of real elements of size (n,),
3212 where ``n`` is the number of independent variables.
3213 args : tuple, optional
3214 Extra arguments passed to the objective function and its
3215 derivatives (`fun`, `jac` and `hess` functions).
3216 method : str or callable, optional
3217 The present documentation is specific to ``method='powell'``, but other
3218 options are available. See documentation for `scipy.optimize.minimize`.
3219 bounds : sequence or `Bounds`, optional
3220 Bounds on decision variables. There are two ways to specify the bounds:
3222 1. Instance of `Bounds` class.
3223 2. Sequence of ``(min, max)`` pairs for each element in `x`. None
3224 is used to specify no bound.
3226 If bounds are not provided, then an unbounded line search will be used.
3227 If bounds are provided and the initial guess is within the bounds, then
3228 every function evaluation throughout the minimization procedure will be
3229 within the bounds. If bounds are provided, the initial guess is outside
3230 the bounds, and `direc` is full rank (or left to default), then some
3231 function evaluations during the first iteration may be outside the
3232 bounds, but every function evaluation after the first iteration will be
3233 within the bounds. If `direc` is not full rank, then some parameters
3234 may not be optimized and the solution is not guaranteed to be within
3235 the bounds.
3237 options : dict, optional
3238 A dictionary of solver options. All methods accept the following
3239 generic options:
3241 maxiter : int
3242 Maximum number of iterations to perform. Depending on the
3243 method each iteration may use several function evaluations.
3244 disp : bool
3245 Set to True to print convergence messages.
3247 See method-specific options for ``method='powell'`` below.
3248 callback : callable, optional
3249 Called after each iteration. The signature is:
3251 ``callback(xk)``
3253 where ``xk`` is the current parameter vector.
3255 Returns
3256 -------
3257 res : OptimizeResult
3258 The optimization result represented as a ``OptimizeResult`` object.
3259 Important attributes are: ``x`` the solution array, ``success`` a
3260 Boolean flag indicating if the optimizer exited successfully and
3261 ``message`` which describes the cause of the termination. See
3262 `OptimizeResult` for a description of other attributes.
3264 Options
3265 -------
3266 disp : bool
3267 Set to True to print convergence messages.
3268 xtol : float
3269 Relative error in solution `xopt` acceptable for convergence.
3270 ftol : float
3271 Relative error in ``fun(xopt)`` acceptable for convergence.
3272 maxiter, maxfev : int
3273 Maximum allowed number of iterations and function evaluations.
3274 Will default to ``N*1000``, where ``N`` is the number of
3275 variables, if neither `maxiter` or `maxfev` is set. If both
3276 `maxiter` and `maxfev` are set, minimization will stop at the
3277 first reached.
3278 direc : ndarray
3279 Initial set of direction vectors for the Powell method.
3280 return_all : bool, optional
3281 Set to True to return a list of the best solution at each of the
3282 iterations.
3283 """
3284 _check_unknown_options(unknown_options)
3285 maxfun = maxfev
3286 retall = return_all
3288 x = asarray(x0).flatten()
3289 if retall:
3290 allvecs = [x]
3291 N = len(x)
3292 # If neither are set, then set both to default
3293 if maxiter is None and maxfun is None:
3294 maxiter = N * 1000
3295 maxfun = N * 1000
3296 elif maxiter is None:
3297 # Convert remaining Nones, to np.inf, unless the other is np.inf, in
3298 # which case use the default to avoid unbounded iteration
3299 if maxfun == np.inf:
3300 maxiter = N * 1000
3301 else:
3302 maxiter = np.inf
3303 elif maxfun is None:
3304 if maxiter == np.inf:
3305 maxfun = N * 1000
3306 else:
3307 maxfun = np.inf
3309 # we need to use a mutable object here that we can update in the
3310 # wrapper function
3311 fcalls, func = _wrap_scalar_function_maxfun_validation(func, args, maxfun)
3313 if direc is None:
3314 direc = eye(N, dtype=float)
3315 else:
3316 direc = asarray(direc, dtype=float)
3317 if np.linalg.matrix_rank(direc) != direc.shape[0]:
3318 warnings.warn("direc input is not full rank, some parameters may "
3319 "not be optimized",
3320 OptimizeWarning, 3)
3322 if bounds is None:
3323 # don't make these arrays of all +/- inf. because
3324 # _linesearch_powell will do an unnecessary check of all the elements.
3325 # just keep them None, _linesearch_powell will not have to check
3326 # all the elements.
3327 lower_bound, upper_bound = None, None
3328 else:
3329 # bounds is standardized in _minimize.py.
3330 lower_bound, upper_bound = bounds.lb, bounds.ub
3331 if np.any(lower_bound > x0) or np.any(x0 > upper_bound):
3332 warnings.warn("Initial guess is not within the specified bounds",
3333 OptimizeWarning, 3)
3335 fval = squeeze(func(x))
3336 x1 = x.copy()
3337 iter = 0
3338 while True:
3339 try:
3340 fx = fval
3341 bigind = 0
3342 delta = 0.0
3343 for i in range(N):
3344 direc1 = direc[i]
3345 fx2 = fval
3346 fval, x, direc1 = _linesearch_powell(func, x, direc1,
3347 tol=xtol * 100,
3348 lower_bound=lower_bound,
3349 upper_bound=upper_bound,
3350 fval=fval)
3351 if (fx2 - fval) > delta:
3352 delta = fx2 - fval
3353 bigind = i
3354 iter += 1
3355 if callback is not None:
3356 callback(x)
3357 if retall:
3358 allvecs.append(x)
3359 bnd = ftol * (np.abs(fx) + np.abs(fval)) + 1e-20
3360 if 2.0 * (fx - fval) <= bnd:
3361 break
3362 if fcalls[0] >= maxfun:
3363 break
3364 if iter >= maxiter:
3365 break
3366 if np.isnan(fx) and np.isnan(fval):
3367 # Ended up in a nan-region: bail out
3368 break
3370 # Construct the extrapolated point
3371 direc1 = x - x1
3372 x1 = x.copy()
3373 # make sure that we don't go outside the bounds when extrapolating
3374 if lower_bound is None and upper_bound is None:
3375 lmax = 1
3376 else:
3377 _, lmax = _line_for_search(x, direc1, lower_bound, upper_bound)
3378 x2 = x + min(lmax, 1) * direc1
3379 fx2 = squeeze(func(x2))
3381 if (fx > fx2):
3382 t = 2.0*(fx + fx2 - 2.0*fval)
3383 temp = (fx - fval - delta)
3384 t *= temp*temp
3385 temp = fx - fx2
3386 t -= delta*temp*temp
3387 if t < 0.0:
3388 fval, x, direc1 = _linesearch_powell(
3389 func, x, direc1,
3390 tol=xtol * 100,
3391 lower_bound=lower_bound,
3392 upper_bound=upper_bound,
3393 fval=fval
3394 )
3395 if np.any(direc1):
3396 direc[bigind] = direc[-1]
3397 direc[-1] = direc1
3398 except _MaxFuncCallError:
3399 break
3401 warnflag = 0
3402 # out of bounds is more urgent than exceeding function evals or iters,
3403 # but I don't want to cause inconsistencies by changing the
3404 # established warning flags for maxfev and maxiter, so the out of bounds
3405 # warning flag becomes 3, but is checked for first.
3406 if bounds and (np.any(lower_bound > x) or np.any(x > upper_bound)):
3407 warnflag = 4
3408 msg = _status_message['out_of_bounds']
3409 elif fcalls[0] >= maxfun:
3410 warnflag = 1
3411 msg = _status_message['maxfev']
3412 if disp:
3413 warnings.warn(msg, RuntimeWarning, 3)
3414 elif iter >= maxiter:
3415 warnflag = 2
3416 msg = _status_message['maxiter']
3417 if disp:
3418 warnings.warn(msg, RuntimeWarning, 3)
3419 elif np.isnan(fval) or np.isnan(x).any():
3420 warnflag = 3
3421 msg = _status_message['nan']
3422 if disp:
3423 warnings.warn(msg, RuntimeWarning, 3)
3424 else:
3425 msg = _status_message['success']
3426 if disp:
3427 print(msg)
3428 print(" Current function value: %f" % fval)
3429 print(" Iterations: %d" % iter)
3430 print(" Function evaluations: %d" % fcalls[0])
3432 result = OptimizeResult(fun=fval, direc=direc, nit=iter, nfev=fcalls[0],
3433 status=warnflag, success=(warnflag == 0),
3434 message=msg, x=x)
3435 if retall:
3436 result['allvecs'] = allvecs
3437 return result
3440def _endprint(x, flag, fval, maxfun, xtol, disp):
3441 if flag == 0:
3442 if disp > 1:
3443 print("\nOptimization terminated successfully;\n"
3444 "The returned value satisfies the termination criteria\n"
3445 "(using xtol = ", xtol, ")")
3446 if flag == 1:
3447 if disp:
3448 print("\nMaximum number of function evaluations exceeded --- "
3449 "increase maxfun argument.\n")
3450 if flag == 2:
3451 if disp:
3452 print("\n{}".format(_status_message['nan']))
3453 return
3456def brute(func, ranges, args=(), Ns=20, full_output=0, finish=fmin,
3457 disp=False, workers=1):
3458 """Minimize a function over a given range by brute force.
3460 Uses the "brute force" method, i.e., computes the function's value
3461 at each point of a multidimensional grid of points, to find the global
3462 minimum of the function.
3464 The function is evaluated everywhere in the range with the datatype of the
3465 first call to the function, as enforced by the ``vectorize`` NumPy
3466 function. The value and type of the function evaluation returned when
3467 ``full_output=True`` are affected in addition by the ``finish`` argument
3468 (see Notes).
3470 The brute force approach is inefficient because the number of grid points
3471 increases exponentially - the number of grid points to evaluate is
3472 ``Ns ** len(x)``. Consequently, even with coarse grid spacing, even
3473 moderately sized problems can take a long time to run, and/or run into
3474 memory limitations.
3476 Parameters
3477 ----------
3478 func : callable
3479 The objective function to be minimized. Must be in the
3480 form ``f(x, *args)``, where ``x`` is the argument in
3481 the form of a 1-D array and ``args`` is a tuple of any
3482 additional fixed parameters needed to completely specify
3483 the function.
3484 ranges : tuple
3485 Each component of the `ranges` tuple must be either a
3486 "slice object" or a range tuple of the form ``(low, high)``.
3487 The program uses these to create the grid of points on which
3488 the objective function will be computed. See `Note 2` for
3489 more detail.
3490 args : tuple, optional
3491 Any additional fixed parameters needed to completely specify
3492 the function.
3493 Ns : int, optional
3494 Number of grid points along the axes, if not otherwise
3495 specified. See `Note2`.
3496 full_output : bool, optional
3497 If True, return the evaluation grid and the objective function's
3498 values on it.
3499 finish : callable, optional
3500 An optimization function that is called with the result of brute force
3501 minimization as initial guess. `finish` should take `func` and
3502 the initial guess as positional arguments, and take `args` as
3503 keyword arguments. It may additionally take `full_output`
3504 and/or `disp` as keyword arguments. Use None if no "polishing"
3505 function is to be used. See Notes for more details.
3506 disp : bool, optional
3507 Set to True to print convergence messages from the `finish` callable.
3508 workers : int or map-like callable, optional
3509 If `workers` is an int the grid is subdivided into `workers`
3510 sections and evaluated in parallel (uses
3511 `multiprocessing.Pool <multiprocessing>`).
3512 Supply `-1` to use all cores available to the Process.
3513 Alternatively supply a map-like callable, such as
3514 `multiprocessing.Pool.map` for evaluating the grid in parallel.
3515 This evaluation is carried out as ``workers(func, iterable)``.
3516 Requires that `func` be pickleable.
3518 .. versionadded:: 1.3.0
3520 Returns
3521 -------
3522 x0 : ndarray
3523 A 1-D array containing the coordinates of a point at which the
3524 objective function had its minimum value. (See `Note 1` for
3525 which point is returned.)
3526 fval : float
3527 Function value at the point `x0`. (Returned when `full_output` is
3528 True.)
3529 grid : tuple
3530 Representation of the evaluation grid. It has the same
3531 length as `x0`. (Returned when `full_output` is True.)
3532 Jout : ndarray
3533 Function values at each point of the evaluation
3534 grid, i.e., ``Jout = func(*grid)``. (Returned
3535 when `full_output` is True.)
3537 See Also
3538 --------
3539 basinhopping, differential_evolution
3541 Notes
3542 -----
3543 *Note 1*: The program finds the gridpoint at which the lowest value
3544 of the objective function occurs. If `finish` is None, that is the
3545 point returned. When the global minimum occurs within (or not very far
3546 outside) the grid's boundaries, and the grid is fine enough, that
3547 point will be in the neighborhood of the global minimum.
3549 However, users often employ some other optimization program to
3550 "polish" the gridpoint values, i.e., to seek a more precise
3551 (local) minimum near `brute's` best gridpoint.
3552 The `brute` function's `finish` option provides a convenient way to do
3553 that. Any polishing program used must take `brute's` output as its
3554 initial guess as a positional argument, and take `brute's` input values
3555 for `args` as keyword arguments, otherwise an error will be raised.
3556 It may additionally take `full_output` and/or `disp` as keyword arguments.
3558 `brute` assumes that the `finish` function returns either an
3559 `OptimizeResult` object or a tuple in the form:
3560 ``(xmin, Jmin, ... , statuscode)``, where ``xmin`` is the minimizing
3561 value of the argument, ``Jmin`` is the minimum value of the objective
3562 function, "..." may be some other returned values (which are not used
3563 by `brute`), and ``statuscode`` is the status code of the `finish` program.
3565 Note that when `finish` is not None, the values returned are those
3566 of the `finish` program, *not* the gridpoint ones. Consequently,
3567 while `brute` confines its search to the input grid points,
3568 the `finish` program's results usually will not coincide with any
3569 gridpoint, and may fall outside the grid's boundary. Thus, if a
3570 minimum only needs to be found over the provided grid points, make
3571 sure to pass in `finish=None`.
3573 *Note 2*: The grid of points is a `numpy.mgrid` object.
3574 For `brute` the `ranges` and `Ns` inputs have the following effect.
3575 Each component of the `ranges` tuple can be either a slice object or a
3576 two-tuple giving a range of values, such as (0, 5). If the component is a
3577 slice object, `brute` uses it directly. If the component is a two-tuple
3578 range, `brute` internally converts it to a slice object that interpolates
3579 `Ns` points from its low-value to its high-value, inclusive.
3581 Examples
3582 --------
3583 We illustrate the use of `brute` to seek the global minimum of a function
3584 of two variables that is given as the sum of a positive-definite
3585 quadratic and two deep "Gaussian-shaped" craters. Specifically, define
3586 the objective function `f` as the sum of three other functions,
3587 ``f = f1 + f2 + f3``. We suppose each of these has a signature
3588 ``(z, *params)``, where ``z = (x, y)``, and ``params`` and the functions
3589 are as defined below.
3591 >>> import numpy as np
3592 >>> params = (2, 3, 7, 8, 9, 10, 44, -1, 2, 26, 1, -2, 0.5)
3593 >>> def f1(z, *params):
3594 ... x, y = z
3595 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3596 ... return (a * x**2 + b * x * y + c * y**2 + d*x + e*y + f)
3598 >>> def f2(z, *params):
3599 ... x, y = z
3600 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3601 ... return (-g*np.exp(-((x-h)**2 + (y-i)**2) / scale))
3603 >>> def f3(z, *params):
3604 ... x, y = z
3605 ... a, b, c, d, e, f, g, h, i, j, k, l, scale = params
3606 ... return (-j*np.exp(-((x-k)**2 + (y-l)**2) / scale))
3608 >>> def f(z, *params):
3609 ... return f1(z, *params) + f2(z, *params) + f3(z, *params)
3611 Thus, the objective function may have local minima near the minimum
3612 of each of the three functions of which it is composed. To
3613 use `fmin` to polish its gridpoint result, we may then continue as
3614 follows:
3616 >>> rranges = (slice(-4, 4, 0.25), slice(-4, 4, 0.25))
3617 >>> from scipy import optimize
3618 >>> resbrute = optimize.brute(f, rranges, args=params, full_output=True,
3619 ... finish=optimize.fmin)
3620 >>> resbrute[0] # global minimum
3621 array([-1.05665192, 1.80834843])
3622 >>> resbrute[1] # function value at global minimum
3623 -3.4085818767
3625 Note that if `finish` had been set to None, we would have gotten the
3626 gridpoint [-1.0 1.75] where the rounded function value is -2.892.
3628 """
3629 N = len(ranges)
3630 if N > 40:
3631 raise ValueError("Brute Force not possible with more "
3632 "than 40 variables.")
3633 lrange = list(ranges)
3634 for k in range(N):
3635 if type(lrange[k]) is not type(slice(None)):
3636 if len(lrange[k]) < 3:
3637 lrange[k] = tuple(lrange[k]) + (complex(Ns),)
3638 lrange[k] = slice(*lrange[k])
3639 if (N == 1):
3640 lrange = lrange[0]
3642 grid = np.mgrid[lrange]
3644 # obtain an array of parameters that is iterable by a map-like callable
3645 inpt_shape = grid.shape
3646 if (N > 1):
3647 grid = np.reshape(grid, (inpt_shape[0], np.prod(inpt_shape[1:]))).T
3649 if not np.iterable(args):
3650 args = (args,)
3652 wrapped_func = _Brute_Wrapper(func, args)
3654 # iterate over input arrays, possibly in parallel
3655 with MapWrapper(pool=workers) as mapper:
3656 Jout = np.array(list(mapper(wrapped_func, grid)))
3657 if (N == 1):
3658 grid = (grid,)
3659 Jout = np.squeeze(Jout)
3660 elif (N > 1):
3661 Jout = np.reshape(Jout, inpt_shape[1:])
3662 grid = np.reshape(grid.T, inpt_shape)
3664 Nshape = shape(Jout)
3666 indx = argmin(Jout.ravel(), axis=-1)
3667 Nindx = np.empty(N, int)
3668 xmin = np.empty(N, float)
3669 for k in range(N - 1, -1, -1):
3670 thisN = Nshape[k]
3671 Nindx[k] = indx % Nshape[k]
3672 indx = indx // thisN
3673 for k in range(N):
3674 xmin[k] = grid[k][tuple(Nindx)]
3676 Jmin = Jout[tuple(Nindx)]
3677 if (N == 1):
3678 grid = grid[0]
3679 xmin = xmin[0]
3681 if callable(finish):
3682 # set up kwargs for `finish` function
3683 finish_args = _getfullargspec(finish).args
3684 finish_kwargs = dict()
3685 if 'full_output' in finish_args:
3686 finish_kwargs['full_output'] = 1
3687 if 'disp' in finish_args:
3688 finish_kwargs['disp'] = disp
3689 elif 'options' in finish_args:
3690 # pass 'disp' as `options`
3691 # (e.g., if `finish` is `minimize`)
3692 finish_kwargs['options'] = {'disp': disp}
3694 # run minimizer
3695 res = finish(func, xmin, args=args, **finish_kwargs)
3697 if isinstance(res, OptimizeResult):
3698 xmin = res.x
3699 Jmin = res.fun
3700 success = res.success
3701 else:
3702 xmin = res[0]
3703 Jmin = res[1]
3704 success = res[-1] == 0
3705 if not success:
3706 if disp:
3707 warnings.warn(
3708 "Either final optimization did not succeed "
3709 "or `finish` does not return `statuscode` as its last "
3710 "argument.", RuntimeWarning, 2)
3712 if full_output:
3713 return xmin, Jmin, grid, Jout
3714 else:
3715 return xmin
3718class _Brute_Wrapper:
3719 """
3720 Object to wrap user cost function for optimize.brute, allowing picklability
3721 """
3723 def __init__(self, f, args):
3724 self.f = f
3725 self.args = [] if args is None else args
3727 def __call__(self, x):
3728 # flatten needed for one dimensional case.
3729 return self.f(np.asarray(x).flatten(), *self.args)
3732def show_options(solver=None, method=None, disp=True):
3733 """
3734 Show documentation for additional options of optimization solvers.
3736 These are method-specific options that can be supplied through the
3737 ``options`` dict.
3739 Parameters
3740 ----------
3741 solver : str
3742 Type of optimization solver. One of 'minimize', 'minimize_scalar',
3743 'root', 'root_scalar', 'linprog', or 'quadratic_assignment'.
3744 method : str, optional
3745 If not given, shows all methods of the specified solver. Otherwise,
3746 show only the options for the specified method. Valid values
3747 corresponds to methods' names of respective solver (e.g., 'BFGS' for
3748 'minimize').
3749 disp : bool, optional
3750 Whether to print the result rather than returning it.
3752 Returns
3753 -------
3754 text
3755 Either None (for disp=True) or the text string (disp=False)
3757 Notes
3758 -----
3759 The solver-specific methods are:
3761 `scipy.optimize.minimize`
3763 - :ref:`Nelder-Mead <optimize.minimize-neldermead>`
3764 - :ref:`Powell <optimize.minimize-powell>`
3765 - :ref:`CG <optimize.minimize-cg>`
3766 - :ref:`BFGS <optimize.minimize-bfgs>`
3767 - :ref:`Newton-CG <optimize.minimize-newtoncg>`
3768 - :ref:`L-BFGS-B <optimize.minimize-lbfgsb>`
3769 - :ref:`TNC <optimize.minimize-tnc>`
3770 - :ref:`COBYLA <optimize.minimize-cobyla>`
3771 - :ref:`SLSQP <optimize.minimize-slsqp>`
3772 - :ref:`dogleg <optimize.minimize-dogleg>`
3773 - :ref:`trust-ncg <optimize.minimize-trustncg>`
3775 `scipy.optimize.root`
3777 - :ref:`hybr <optimize.root-hybr>`
3778 - :ref:`lm <optimize.root-lm>`
3779 - :ref:`broyden1 <optimize.root-broyden1>`
3780 - :ref:`broyden2 <optimize.root-broyden2>`
3781 - :ref:`anderson <optimize.root-anderson>`
3782 - :ref:`linearmixing <optimize.root-linearmixing>`
3783 - :ref:`diagbroyden <optimize.root-diagbroyden>`
3784 - :ref:`excitingmixing <optimize.root-excitingmixing>`
3785 - :ref:`krylov <optimize.root-krylov>`
3786 - :ref:`df-sane <optimize.root-dfsane>`
3788 `scipy.optimize.minimize_scalar`
3790 - :ref:`brent <optimize.minimize_scalar-brent>`
3791 - :ref:`golden <optimize.minimize_scalar-golden>`
3792 - :ref:`bounded <optimize.minimize_scalar-bounded>`
3794 `scipy.optimize.root_scalar`
3796 - :ref:`bisect <optimize.root_scalar-bisect>`
3797 - :ref:`brentq <optimize.root_scalar-brentq>`
3798 - :ref:`brenth <optimize.root_scalar-brenth>`
3799 - :ref:`ridder <optimize.root_scalar-ridder>`
3800 - :ref:`toms748 <optimize.root_scalar-toms748>`
3801 - :ref:`newton <optimize.root_scalar-newton>`
3802 - :ref:`secant <optimize.root_scalar-secant>`
3803 - :ref:`halley <optimize.root_scalar-halley>`
3805 `scipy.optimize.linprog`
3807 - :ref:`simplex <optimize.linprog-simplex>`
3808 - :ref:`interior-point <optimize.linprog-interior-point>`
3809 - :ref:`revised simplex <optimize.linprog-revised_simplex>`
3810 - :ref:`highs <optimize.linprog-highs>`
3811 - :ref:`highs-ds <optimize.linprog-highs-ds>`
3812 - :ref:`highs-ipm <optimize.linprog-highs-ipm>`
3814 `scipy.optimize.quadratic_assignment`
3816 - :ref:`faq <optimize.qap-faq>`
3817 - :ref:`2opt <optimize.qap-2opt>`
3819 Examples
3820 --------
3821 We can print documentations of a solver in stdout:
3823 >>> from scipy.optimize import show_options
3824 >>> show_options(solver="minimize")
3825 ...
3827 Specifying a method is possible:
3829 >>> show_options(solver="minimize", method="Nelder-Mead")
3830 ...
3832 We can also get the documentations as a string:
3834 >>> show_options(solver="minimize", method="Nelder-Mead", disp=False)
3835 Minimization of scalar function of one or more variables using the ...
3837 """
3838 import textwrap
3840 doc_routines = {
3841 'minimize': (
3842 ('bfgs', 'scipy.optimize._optimize._minimize_bfgs'),
3843 ('cg', 'scipy.optimize._optimize._minimize_cg'),
3844 ('cobyla', 'scipy.optimize._cobyla_py._minimize_cobyla'),
3845 ('dogleg', 'scipy.optimize._trustregion_dogleg._minimize_dogleg'),
3846 ('l-bfgs-b', 'scipy.optimize._lbfgsb_py._minimize_lbfgsb'),
3847 ('nelder-mead', 'scipy.optimize._optimize._minimize_neldermead'),
3848 ('newton-cg', 'scipy.optimize._optimize._minimize_newtoncg'),
3849 ('powell', 'scipy.optimize._optimize._minimize_powell'),
3850 ('slsqp', 'scipy.optimize._slsqp_py._minimize_slsqp'),
3851 ('tnc', 'scipy.optimize._tnc._minimize_tnc'),
3852 ('trust-ncg',
3853 'scipy.optimize._trustregion_ncg._minimize_trust_ncg'),
3854 ('trust-constr',
3855 'scipy.optimize._trustregion_constr.'
3856 '_minimize_trustregion_constr'),
3857 ('trust-exact',
3858 'scipy.optimize._trustregion_exact._minimize_trustregion_exact'),
3859 ('trust-krylov',
3860 'scipy.optimize._trustregion_krylov._minimize_trust_krylov'),
3861 ),
3862 'root': (
3863 ('hybr', 'scipy.optimize._minpack_py._root_hybr'),
3864 ('lm', 'scipy.optimize._root._root_leastsq'),
3865 ('broyden1', 'scipy.optimize._root._root_broyden1_doc'),
3866 ('broyden2', 'scipy.optimize._root._root_broyden2_doc'),
3867 ('anderson', 'scipy.optimize._root._root_anderson_doc'),
3868 ('diagbroyden', 'scipy.optimize._root._root_diagbroyden_doc'),
3869 ('excitingmixing', 'scipy.optimize._root._root_excitingmixing_doc'),
3870 ('linearmixing', 'scipy.optimize._root._root_linearmixing_doc'),
3871 ('krylov', 'scipy.optimize._root._root_krylov_doc'),
3872 ('df-sane', 'scipy.optimize._spectral._root_df_sane'),
3873 ),
3874 'root_scalar': (
3875 ('bisect', 'scipy.optimize._root_scalar._root_scalar_bisect_doc'),
3876 ('brentq', 'scipy.optimize._root_scalar._root_scalar_brentq_doc'),
3877 ('brenth', 'scipy.optimize._root_scalar._root_scalar_brenth_doc'),
3878 ('ridder', 'scipy.optimize._root_scalar._root_scalar_ridder_doc'),
3879 ('toms748', 'scipy.optimize._root_scalar._root_scalar_toms748_doc'),
3880 ('secant', 'scipy.optimize._root_scalar._root_scalar_secant_doc'),
3881 ('newton', 'scipy.optimize._root_scalar._root_scalar_newton_doc'),
3882 ('halley', 'scipy.optimize._root_scalar._root_scalar_halley_doc'),
3883 ),
3884 'linprog': (
3885 ('simplex', 'scipy.optimize._linprog._linprog_simplex_doc'),
3886 ('interior-point', 'scipy.optimize._linprog._linprog_ip_doc'),
3887 ('revised simplex', 'scipy.optimize._linprog._linprog_rs_doc'),
3888 ('highs-ipm', 'scipy.optimize._linprog._linprog_highs_ipm_doc'),
3889 ('highs-ds', 'scipy.optimize._linprog._linprog_highs_ds_doc'),
3890 ('highs', 'scipy.optimize._linprog._linprog_highs_doc'),
3891 ),
3892 'quadratic_assignment': (
3893 ('faq', 'scipy.optimize._qap._quadratic_assignment_faq'),
3894 ('2opt', 'scipy.optimize._qap._quadratic_assignment_2opt'),
3895 ),
3896 'minimize_scalar': (
3897 ('brent', 'scipy.optimize._optimize._minimize_scalar_brent'),
3898 ('bounded', 'scipy.optimize._optimize._minimize_scalar_bounded'),
3899 ('golden', 'scipy.optimize._optimize._minimize_scalar_golden'),
3900 ),
3901 }
3903 if solver is None:
3904 text = ["\n\n\n========\n", "minimize\n", "========\n"]
3905 text.append(show_options('minimize', disp=False))
3906 text.extend(["\n\n===============\n", "minimize_scalar\n",
3907 "===============\n"])
3908 text.append(show_options('minimize_scalar', disp=False))
3909 text.extend(["\n\n\n====\n", "root\n",
3910 "====\n"])
3911 text.append(show_options('root', disp=False))
3912 text.extend(['\n\n\n=======\n', 'linprog\n',
3913 '=======\n'])
3914 text.append(show_options('linprog', disp=False))
3915 text = "".join(text)
3916 else:
3917 solver = solver.lower()
3918 if solver not in doc_routines:
3919 raise ValueError('Unknown solver %r' % (solver,))
3921 if method is None:
3922 text = []
3923 for name, _ in doc_routines[solver]:
3924 text.extend(["\n\n" + name, "\n" + "="*len(name) + "\n\n"])
3925 text.append(show_options(solver, name, disp=False))
3926 text = "".join(text)
3927 else:
3928 method = method.lower()
3929 methods = dict(doc_routines[solver])
3930 if method not in methods:
3931 raise ValueError("Unknown method %r" % (method,))
3932 name = methods[method]
3934 # Import function object
3935 parts = name.split('.')
3936 mod_name = ".".join(parts[:-1])
3937 __import__(mod_name)
3938 obj = getattr(sys.modules[mod_name], parts[-1])
3940 # Get doc
3941 doc = obj.__doc__
3942 if doc is not None:
3943 text = textwrap.dedent(doc).strip()
3944 else:
3945 text = ""
3947 if disp:
3948 print(text)
3949 return
3950 else:
3951 return text