_minimize.py 45 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038
  1. """
  2. Unified interfaces to minimization algorithms.
  3. Functions
  4. ---------
  5. - minimize : minimization of a function of several variables.
  6. - minimize_scalar : minimization of a function of one variable.
  7. """
  8. __all__ = ['minimize', 'minimize_scalar']
  9. from warnings import warn
  10. import numpy as np
  11. # unconstrained minimization
  12. from ._optimize import (_minimize_neldermead, _minimize_powell, _minimize_cg,
  13. _minimize_bfgs, _minimize_newtoncg,
  14. _minimize_scalar_brent, _minimize_scalar_bounded,
  15. _minimize_scalar_golden, MemoizeJac, OptimizeResult)
  16. from ._trustregion_dogleg import _minimize_dogleg
  17. from ._trustregion_ncg import _minimize_trust_ncg
  18. from ._trustregion_krylov import _minimize_trust_krylov
  19. from ._trustregion_exact import _minimize_trustregion_exact
  20. from ._trustregion_constr import _minimize_trustregion_constr
  21. # constrained minimization
  22. from ._lbfgsb_py import _minimize_lbfgsb
  23. from ._tnc import _minimize_tnc
  24. from ._cobyla_py import _minimize_cobyla
  25. from ._slsqp_py import _minimize_slsqp
  26. from ._constraints import (old_bound_to_new, new_bounds_to_old,
  27. old_constraint_to_new, new_constraint_to_old,
  28. NonlinearConstraint, LinearConstraint, Bounds,
  29. PreparedConstraint)
  30. from ._differentiable_functions import FD_METHODS
  31. MINIMIZE_METHODS = ['nelder-mead', 'powell', 'cg', 'bfgs', 'newton-cg',
  32. 'l-bfgs-b', 'tnc', 'cobyla', 'slsqp', 'trust-constr',
  33. 'dogleg', 'trust-ncg', 'trust-exact', 'trust-krylov']
  34. MINIMIZE_SCALAR_METHODS = ['brent', 'bounded', 'golden']
  35. def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
  36. hessp=None, bounds=None, constraints=(), tol=None,
  37. callback=None, options=None):
  38. """Minimization of scalar function of one or more variables.
  39. Parameters
  40. ----------
  41. fun : callable
  42. The objective function to be minimized.
  43. ``fun(x, *args) -> float``
  44. where ``x`` is a 1-D array with shape (n,) and ``args``
  45. is a tuple of the fixed parameters needed to completely
  46. specify the function.
  47. x0 : ndarray, shape (n,)
  48. Initial guess. Array of real elements of size (n,),
  49. where ``n`` is the number of independent variables.
  50. args : tuple, optional
  51. Extra arguments passed to the objective function and its
  52. derivatives (`fun`, `jac` and `hess` functions).
  53. method : str or callable, optional
  54. Type of solver. Should be one of
  55. - 'Nelder-Mead' :ref:`(see here) <optimize.minimize-neldermead>`
  56. - 'Powell' :ref:`(see here) <optimize.minimize-powell>`
  57. - 'CG' :ref:`(see here) <optimize.minimize-cg>`
  58. - 'BFGS' :ref:`(see here) <optimize.minimize-bfgs>`
  59. - 'Newton-CG' :ref:`(see here) <optimize.minimize-newtoncg>`
  60. - 'L-BFGS-B' :ref:`(see here) <optimize.minimize-lbfgsb>`
  61. - 'TNC' :ref:`(see here) <optimize.minimize-tnc>`
  62. - 'COBYLA' :ref:`(see here) <optimize.minimize-cobyla>`
  63. - 'SLSQP' :ref:`(see here) <optimize.minimize-slsqp>`
  64. - 'trust-constr':ref:`(see here) <optimize.minimize-trustconstr>`
  65. - 'dogleg' :ref:`(see here) <optimize.minimize-dogleg>`
  66. - 'trust-ncg' :ref:`(see here) <optimize.minimize-trustncg>`
  67. - 'trust-exact' :ref:`(see here) <optimize.minimize-trustexact>`
  68. - 'trust-krylov' :ref:`(see here) <optimize.minimize-trustkrylov>`
  69. - custom - a callable object, see below for description.
  70. If not given, chosen to be one of ``BFGS``, ``L-BFGS-B``, ``SLSQP``,
  71. depending on whether or not the problem has constraints or bounds.
  72. jac : {callable, '2-point', '3-point', 'cs', bool}, optional
  73. Method for computing the gradient vector. Only for CG, BFGS,
  74. Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov,
  75. trust-exact and trust-constr.
  76. If it is a callable, it should be a function that returns the gradient
  77. vector:
  78. ``jac(x, *args) -> array_like, shape (n,)``
  79. where ``x`` is an array with shape (n,) and ``args`` is a tuple with
  80. the fixed parameters. If `jac` is a Boolean and is True, `fun` is
  81. assumed to return a tuple ``(f, g)`` containing the objective
  82. function and the gradient.
  83. Methods 'Newton-CG', 'trust-ncg', 'dogleg', 'trust-exact', and
  84. 'trust-krylov' require that either a callable be supplied, or that
  85. `fun` return the objective and gradient.
  86. If None or False, the gradient will be estimated using 2-point finite
  87. difference estimation with an absolute step size.
  88. Alternatively, the keywords {'2-point', '3-point', 'cs'} can be used
  89. to select a finite difference scheme for numerical estimation of the
  90. gradient with a relative step size. These finite difference schemes
  91. obey any specified `bounds`.
  92. hess : {callable, '2-point', '3-point', 'cs', HessianUpdateStrategy}, optional
  93. Method for computing the Hessian matrix. Only for Newton-CG, dogleg,
  94. trust-ncg, trust-krylov, trust-exact and trust-constr.
  95. If it is callable, it should return the Hessian matrix:
  96. ``hess(x, *args) -> {LinearOperator, spmatrix, array}, (n, n)``
  97. where ``x`` is a (n,) ndarray and ``args`` is a tuple with the fixed
  98. parameters.
  99. The keywords {'2-point', '3-point', 'cs'} can also be used to select
  100. a finite difference scheme for numerical estimation of the hessian.
  101. Alternatively, objects implementing the `HessianUpdateStrategy`
  102. interface can be used to approximate the Hessian. Available
  103. quasi-Newton methods implementing this interface are:
  104. - `BFGS`;
  105. - `SR1`.
  106. Not all of the options are available for each of the methods; for
  107. availability refer to the notes.
  108. hessp : callable, optional
  109. Hessian of objective function times an arbitrary vector p. Only for
  110. Newton-CG, trust-ncg, trust-krylov, trust-constr.
  111. Only one of `hessp` or `hess` needs to be given. If `hess` is
  112. provided, then `hessp` will be ignored. `hessp` must compute the
  113. Hessian times an arbitrary vector:
  114. ``hessp(x, p, *args) -> ndarray shape (n,)``
  115. where ``x`` is a (n,) ndarray, ``p`` is an arbitrary vector with
  116. dimension (n,) and ``args`` is a tuple with the fixed
  117. parameters.
  118. bounds : sequence or `Bounds`, optional
  119. Bounds on variables for Nelder-Mead, L-BFGS-B, TNC, SLSQP, Powell, and
  120. trust-constr methods. There are two ways to specify the bounds:
  121. 1. Instance of `Bounds` class.
  122. 2. Sequence of ``(min, max)`` pairs for each element in `x`. None
  123. is used to specify no bound.
  124. constraints : {Constraint, dict} or List of {Constraint, dict}, optional
  125. Constraints definition. Only for COBYLA, SLSQP and trust-constr.
  126. Constraints for 'trust-constr' are defined as a single object or a
  127. list of objects specifying constraints to the optimization problem.
  128. Available constraints are:
  129. - `LinearConstraint`
  130. - `NonlinearConstraint`
  131. Constraints for COBYLA, SLSQP are defined as a list of dictionaries.
  132. Each dictionary with fields:
  133. type : str
  134. Constraint type: 'eq' for equality, 'ineq' for inequality.
  135. fun : callable
  136. The function defining the constraint.
  137. jac : callable, optional
  138. The Jacobian of `fun` (only for SLSQP).
  139. args : sequence, optional
  140. Extra arguments to be passed to the function and Jacobian.
  141. Equality constraint means that the constraint function result is to
  142. be zero whereas inequality means that it is to be non-negative.
  143. Note that COBYLA only supports inequality constraints.
  144. tol : float, optional
  145. Tolerance for termination. When `tol` is specified, the selected
  146. minimization algorithm sets some relevant solver-specific tolerance(s)
  147. equal to `tol`. For detailed control, use solver-specific
  148. options.
  149. options : dict, optional
  150. A dictionary of solver options. All methods except `TNC` accept the
  151. following generic options:
  152. maxiter : int
  153. Maximum number of iterations to perform. Depending on the
  154. method each iteration may use several function evaluations.
  155. For `TNC` use `maxfun` instead of `maxiter`.
  156. disp : bool
  157. Set to True to print convergence messages.
  158. For method-specific options, see :func:`show_options()`.
  159. callback : callable, optional
  160. Called after each iteration. For 'trust-constr' it is a callable with
  161. the signature:
  162. ``callback(xk, OptimizeResult state) -> bool``
  163. where ``xk`` is the current parameter vector. and ``state``
  164. is an `OptimizeResult` object, with the same fields
  165. as the ones from the return. If callback returns True
  166. the algorithm execution is terminated.
  167. For all the other methods, the signature is:
  168. ``callback(xk)``
  169. where ``xk`` is the current parameter vector.
  170. Returns
  171. -------
  172. res : OptimizeResult
  173. The optimization result represented as a ``OptimizeResult`` object.
  174. Important attributes are: ``x`` the solution array, ``success`` a
  175. Boolean flag indicating if the optimizer exited successfully and
  176. ``message`` which describes the cause of the termination. See
  177. `OptimizeResult` for a description of other attributes.
  178. See also
  179. --------
  180. minimize_scalar : Interface to minimization algorithms for scalar
  181. univariate functions
  182. show_options : Additional options accepted by the solvers
  183. Notes
  184. -----
  185. This section describes the available solvers that can be selected by the
  186. 'method' parameter. The default method is *BFGS*.
  187. **Unconstrained minimization**
  188. Method :ref:`CG <optimize.minimize-cg>` uses a nonlinear conjugate
  189. gradient algorithm by Polak and Ribiere, a variant of the
  190. Fletcher-Reeves method described in [5]_ pp.120-122. Only the
  191. first derivatives are used.
  192. Method :ref:`BFGS <optimize.minimize-bfgs>` uses the quasi-Newton
  193. method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) [5]_
  194. pp. 136. It uses the first derivatives only. BFGS has proven good
  195. performance even for non-smooth optimizations. This method also
  196. returns an approximation of the Hessian inverse, stored as
  197. `hess_inv` in the OptimizeResult object.
  198. Method :ref:`Newton-CG <optimize.minimize-newtoncg>` uses a
  199. Newton-CG algorithm [5]_ pp. 168 (also known as the truncated
  200. Newton method). It uses a CG method to the compute the search
  201. direction. See also *TNC* method for a box-constrained
  202. minimization with a similar algorithm. Suitable for large-scale
  203. problems.
  204. Method :ref:`dogleg <optimize.minimize-dogleg>` uses the dog-leg
  205. trust-region algorithm [5]_ for unconstrained minimization. This
  206. algorithm requires the gradient and Hessian; furthermore the
  207. Hessian is required to be positive definite.
  208. Method :ref:`trust-ncg <optimize.minimize-trustncg>` uses the
  209. Newton conjugate gradient trust-region algorithm [5]_ for
  210. unconstrained minimization. This algorithm requires the gradient
  211. and either the Hessian or a function that computes the product of
  212. the Hessian with a given vector. Suitable for large-scale problems.
  213. Method :ref:`trust-krylov <optimize.minimize-trustkrylov>` uses
  214. the Newton GLTR trust-region algorithm [14]_, [15]_ for unconstrained
  215. minimization. This algorithm requires the gradient
  216. and either the Hessian or a function that computes the product of
  217. the Hessian with a given vector. Suitable for large-scale problems.
  218. On indefinite problems it requires usually less iterations than the
  219. `trust-ncg` method and is recommended for medium and large-scale problems.
  220. Method :ref:`trust-exact <optimize.minimize-trustexact>`
  221. is a trust-region method for unconstrained minimization in which
  222. quadratic subproblems are solved almost exactly [13]_. This
  223. algorithm requires the gradient and the Hessian (which is
  224. *not* required to be positive definite). It is, in many
  225. situations, the Newton method to converge in fewer iterations
  226. and the most recommended for small and medium-size problems.
  227. **Bound-Constrained minimization**
  228. Method :ref:`Nelder-Mead <optimize.minimize-neldermead>` uses the
  229. Simplex algorithm [1]_, [2]_. This algorithm is robust in many
  230. applications. However, if numerical computation of derivative can be
  231. trusted, other algorithms using the first and/or second derivatives
  232. information might be preferred for their better performance in
  233. general.
  234. Method :ref:`L-BFGS-B <optimize.minimize-lbfgsb>` uses the L-BFGS-B
  235. algorithm [6]_, [7]_ for bound constrained minimization.
  236. Method :ref:`Powell <optimize.minimize-powell>` is a modification
  237. of Powell's method [3]_, [4]_ which is a conjugate direction
  238. method. It performs sequential one-dimensional minimizations along
  239. each vector of the directions set (`direc` field in `options` and
  240. `info`), which is updated at each iteration of the main
  241. minimization loop. The function need not be differentiable, and no
  242. derivatives are taken. If bounds are not provided, then an
  243. unbounded line search will be used. If bounds are provided and
  244. the initial guess is within the bounds, then every function
  245. evaluation throughout the minimization procedure will be within
  246. the bounds. If bounds are provided, the initial guess is outside
  247. the bounds, and `direc` is full rank (default has full rank), then
  248. some function evaluations during the first iteration may be
  249. outside the bounds, but every function evaluation after the first
  250. iteration will be within the bounds. If `direc` is not full rank,
  251. then some parameters may not be optimized and the solution is not
  252. guaranteed to be within the bounds.
  253. Method :ref:`TNC <optimize.minimize-tnc>` uses a truncated Newton
  254. algorithm [5]_, [8]_ to minimize a function with variables subject
  255. to bounds. This algorithm uses gradient information; it is also
  256. called Newton Conjugate-Gradient. It differs from the *Newton-CG*
  257. method described above as it wraps a C implementation and allows
  258. each variable to be given upper and lower bounds.
  259. **Constrained Minimization**
  260. Method :ref:`COBYLA <optimize.minimize-cobyla>` uses the
  261. Constrained Optimization BY Linear Approximation (COBYLA) method
  262. [9]_, [10]_, [11]_. The algorithm is based on linear
  263. approximations to the objective function and each constraint. The
  264. method wraps a FORTRAN implementation of the algorithm. The
  265. constraints functions 'fun' may return either a single number
  266. or an array or list of numbers.
  267. Method :ref:`SLSQP <optimize.minimize-slsqp>` uses Sequential
  268. Least SQuares Programming to minimize a function of several
  269. variables with any combination of bounds, equality and inequality
  270. constraints. The method wraps the SLSQP Optimization subroutine
  271. originally implemented by Dieter Kraft [12]_. Note that the
  272. wrapper handles infinite values in bounds by converting them into
  273. large floating values.
  274. Method :ref:`trust-constr <optimize.minimize-trustconstr>` is a
  275. trust-region algorithm for constrained optimization. It swiches
  276. between two implementations depending on the problem definition.
  277. It is the most versatile constrained minimization algorithm
  278. implemented in SciPy and the most appropriate for large-scale problems.
  279. For equality constrained problems it is an implementation of Byrd-Omojokun
  280. Trust-Region SQP method described in [17]_ and in [5]_, p. 549. When
  281. inequality constraints are imposed as well, it swiches to the trust-region
  282. interior point method described in [16]_. This interior point algorithm,
  283. in turn, solves inequality constraints by introducing slack variables
  284. and solving a sequence of equality-constrained barrier problems
  285. for progressively smaller values of the barrier parameter.
  286. The previously described equality constrained SQP method is
  287. used to solve the subproblems with increasing levels of accuracy
  288. as the iterate gets closer to a solution.
  289. **Finite-Difference Options**
  290. For Method :ref:`trust-constr <optimize.minimize-trustconstr>`
  291. the gradient and the Hessian may be approximated using
  292. three finite-difference schemes: {'2-point', '3-point', 'cs'}.
  293. The scheme 'cs' is, potentially, the most accurate but it
  294. requires the function to correctly handle complex inputs and to
  295. be differentiable in the complex plane. The scheme '3-point' is more
  296. accurate than '2-point' but requires twice as many operations. If the
  297. gradient is estimated via finite-differences the Hessian must be
  298. estimated using one of the quasi-Newton strategies.
  299. **Method specific options for the** `hess` **keyword**
  300. +--------------+------+----------+-------------------------+-----+
  301. | method/Hess | None | callable | '2-point/'3-point'/'cs' | HUS |
  302. +==============+======+==========+=========================+=====+
  303. | Newton-CG | x | (n, n) | x | x |
  304. | | | LO | | |
  305. +--------------+------+----------+-------------------------+-----+
  306. | dogleg | | (n, n) | | |
  307. +--------------+------+----------+-------------------------+-----+
  308. | trust-ncg | | (n, n) | x | x |
  309. +--------------+------+----------+-------------------------+-----+
  310. | trust-krylov | | (n, n) | x | x |
  311. +--------------+------+----------+-------------------------+-----+
  312. | trust-exact | | (n, n) | | |
  313. +--------------+------+----------+-------------------------+-----+
  314. | trust-constr | x | (n, n) | x | x |
  315. | | | LO | | |
  316. | | | sp | | |
  317. +--------------+------+----------+-------------------------+-----+
  318. where LO=LinearOperator, sp=Sparse matrix, HUS=HessianUpdateStrategy
  319. **Custom minimizers**
  320. It may be useful to pass a custom minimization method, for example
  321. when using a frontend to this method such as `scipy.optimize.basinhopping`
  322. or a different library. You can simply pass a callable as the ``method``
  323. parameter.
  324. The callable is called as ``method(fun, x0, args, **kwargs, **options)``
  325. where ``kwargs`` corresponds to any other parameters passed to `minimize`
  326. (such as `callback`, `hess`, etc.), except the `options` dict, which has
  327. its contents also passed as `method` parameters pair by pair. Also, if
  328. `jac` has been passed as a bool type, `jac` and `fun` are mangled so that
  329. `fun` returns just the function values and `jac` is converted to a function
  330. returning the Jacobian. The method shall return an `OptimizeResult`
  331. object.
  332. The provided `method` callable must be able to accept (and possibly ignore)
  333. arbitrary parameters; the set of parameters accepted by `minimize` may
  334. expand in future versions and then these parameters will be passed to
  335. the method. You can find an example in the scipy.optimize tutorial.
  336. References
  337. ----------
  338. .. [1] Nelder, J A, and R Mead. 1965. A Simplex Method for Function
  339. Minimization. The Computer Journal 7: 308-13.
  340. .. [2] Wright M H. 1996. Direct search methods: Once scorned, now
  341. respectable, in Numerical Analysis 1995: Proceedings of the 1995
  342. Dundee Biennial Conference in Numerical Analysis (Eds. D F
  343. Griffiths and G A Watson). Addison Wesley Longman, Harlow, UK.
  344. 191-208.
  345. .. [3] Powell, M J D. 1964. An efficient method for finding the minimum of
  346. a function of several variables without calculating derivatives. The
  347. Computer Journal 7: 155-162.
  348. .. [4] Press W, S A Teukolsky, W T Vetterling and B P Flannery.
  349. Numerical Recipes (any edition), Cambridge University Press.
  350. .. [5] Nocedal, J, and S J Wright. 2006. Numerical Optimization.
  351. Springer New York.
  352. .. [6] Byrd, R H and P Lu and J. Nocedal. 1995. A Limited Memory
  353. Algorithm for Bound Constrained Optimization. SIAM Journal on
  354. Scientific and Statistical Computing 16 (5): 1190-1208.
  355. .. [7] Zhu, C and R H Byrd and J Nocedal. 1997. L-BFGS-B: Algorithm
  356. 778: L-BFGS-B, FORTRAN routines for large scale bound constrained
  357. optimization. ACM Transactions on Mathematical Software 23 (4):
  358. 550-560.
  359. .. [8] Nash, S G. Newton-Type Minimization Via the Lanczos Method.
  360. 1984. SIAM Journal of Numerical Analysis 21: 770-778.
  361. .. [9] Powell, M J D. A direct search optimization method that models
  362. the objective and constraint functions by linear interpolation.
  363. 1994. Advances in Optimization and Numerical Analysis, eds. S. Gomez
  364. and J-P Hennart, Kluwer Academic (Dordrecht), 51-67.
  365. .. [10] Powell M J D. Direct search algorithms for optimization
  366. calculations. 1998. Acta Numerica 7: 287-336.
  367. .. [11] Powell M J D. A view of algorithms for optimization without
  368. derivatives. 2007.Cambridge University Technical Report DAMTP
  369. 2007/NA03
  370. .. [12] Kraft, D. A software package for sequential quadratic
  371. programming. 1988. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace
  372. Center -- Institute for Flight Mechanics, Koln, Germany.
  373. .. [13] Conn, A. R., Gould, N. I., and Toint, P. L.
  374. Trust region methods. 2000. Siam. pp. 169-200.
  375. .. [14] F. Lenders, C. Kirches, A. Potschka: "trlib: A vector-free
  376. implementation of the GLTR method for iterative solution of
  377. the trust region problem", :arxiv:`1611.04718`
  378. .. [15] N. Gould, S. Lucidi, M. Roma, P. Toint: "Solving the
  379. Trust-Region Subproblem using the Lanczos Method",
  380. SIAM J. Optim., 9(2), 504--525, (1999).
  381. .. [16] Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal. 1999.
  382. An interior point algorithm for large-scale nonlinear programming.
  383. SIAM Journal on Optimization 9.4: 877-900.
  384. .. [17] Lalee, Marucha, Jorge Nocedal, and Todd Plantega. 1998. On the
  385. implementation of an algorithm for large-scale equality constrained
  386. optimization. SIAM Journal on Optimization 8.3: 682-706.
  387. Examples
  388. --------
  389. Let us consider the problem of minimizing the Rosenbrock function. This
  390. function (and its respective derivatives) is implemented in `rosen`
  391. (resp. `rosen_der`, `rosen_hess`) in the `scipy.optimize`.
  392. >>> from scipy.optimize import minimize, rosen, rosen_der
  393. A simple application of the *Nelder-Mead* method is:
  394. >>> x0 = [1.3, 0.7, 0.8, 1.9, 1.2]
  395. >>> res = minimize(rosen, x0, method='Nelder-Mead', tol=1e-6)
  396. >>> res.x
  397. array([ 1., 1., 1., 1., 1.])
  398. Now using the *BFGS* algorithm, using the first derivative and a few
  399. options:
  400. >>> res = minimize(rosen, x0, method='BFGS', jac=rosen_der,
  401. ... options={'gtol': 1e-6, 'disp': True})
  402. Optimization terminated successfully.
  403. Current function value: 0.000000
  404. Iterations: 26
  405. Function evaluations: 31
  406. Gradient evaluations: 31
  407. >>> res.x
  408. array([ 1., 1., 1., 1., 1.])
  409. >>> print(res.message)
  410. Optimization terminated successfully.
  411. >>> res.hess_inv
  412. array([[ 0.00749589, 0.01255155, 0.02396251, 0.04750988, 0.09495377], # may vary
  413. [ 0.01255155, 0.02510441, 0.04794055, 0.09502834, 0.18996269],
  414. [ 0.02396251, 0.04794055, 0.09631614, 0.19092151, 0.38165151],
  415. [ 0.04750988, 0.09502834, 0.19092151, 0.38341252, 0.7664427 ],
  416. [ 0.09495377, 0.18996269, 0.38165151, 0.7664427, 1.53713523]])
  417. Next, consider a minimization problem with several constraints (namely
  418. Example 16.4 from [5]_). The objective function is:
  419. >>> fun = lambda x: (x[0] - 1)**2 + (x[1] - 2.5)**2
  420. There are three constraints defined as:
  421. >>> cons = ({'type': 'ineq', 'fun': lambda x: x[0] - 2 * x[1] + 2},
  422. ... {'type': 'ineq', 'fun': lambda x: -x[0] - 2 * x[1] + 6},
  423. ... {'type': 'ineq', 'fun': lambda x: -x[0] + 2 * x[1] + 2})
  424. And variables must be positive, hence the following bounds:
  425. >>> bnds = ((0, None), (0, None))
  426. The optimization problem is solved using the SLSQP method as:
  427. >>> res = minimize(fun, (2, 0), method='SLSQP', bounds=bnds,
  428. ... constraints=cons)
  429. It should converge to the theoretical solution (1.4 ,1.7).
  430. """
  431. x0 = np.atleast_1d(np.asarray(x0))
  432. if x0.ndim != 1:
  433. message = ('Use of `minimize` with `x0.ndim != 1` is deprecated. '
  434. 'Currently, singleton dimensions will be removed from '
  435. '`x0`, but an error will be raised in SciPy 1.11.0.')
  436. warn(message, DeprecationWarning, stacklevel=2)
  437. x0 = np.atleast_1d(np.squeeze(x0))
  438. if x0.dtype.kind in np.typecodes["AllInteger"]:
  439. x0 = np.asarray(x0, dtype=float)
  440. if not isinstance(args, tuple):
  441. args = (args,)
  442. if method is None:
  443. # Select automatically
  444. if constraints:
  445. method = 'SLSQP'
  446. elif bounds is not None:
  447. method = 'L-BFGS-B'
  448. else:
  449. method = 'BFGS'
  450. if callable(method):
  451. meth = "_custom"
  452. else:
  453. meth = method.lower()
  454. if options is None:
  455. options = {}
  456. # check if optional parameters are supported by the selected method
  457. # - jac
  458. if meth in ('nelder-mead', 'powell', 'cobyla') and bool(jac):
  459. warn('Method %s does not use gradient information (jac).' % method,
  460. RuntimeWarning)
  461. # - hess
  462. if meth not in ('newton-cg', 'dogleg', 'trust-ncg', 'trust-constr',
  463. 'trust-krylov', 'trust-exact', '_custom') and hess is not None:
  464. warn('Method %s does not use Hessian information (hess).' % method,
  465. RuntimeWarning)
  466. # - hessp
  467. if meth not in ('newton-cg', 'trust-ncg', 'trust-constr',
  468. 'trust-krylov', '_custom') \
  469. and hessp is not None:
  470. warn('Method %s does not use Hessian-vector product '
  471. 'information (hessp).' % method, RuntimeWarning)
  472. # - constraints or bounds
  473. if (meth not in ('cobyla', 'slsqp', 'trust-constr', '_custom') and
  474. np.any(constraints)):
  475. warn('Method %s cannot handle constraints.' % method,
  476. RuntimeWarning)
  477. if meth not in ('nelder-mead', 'powell', 'l-bfgs-b', 'tnc', 'slsqp',
  478. 'trust-constr', '_custom') and bounds is not None:
  479. warn('Method %s cannot handle bounds.' % method,
  480. RuntimeWarning)
  481. # - return_all
  482. if (meth in ('l-bfgs-b', 'tnc', 'cobyla', 'slsqp') and
  483. options.get('return_all', False)):
  484. warn('Method %s does not support the return_all option.' % method,
  485. RuntimeWarning)
  486. # check gradient vector
  487. if callable(jac):
  488. pass
  489. elif jac is True:
  490. # fun returns func and grad
  491. fun = MemoizeJac(fun)
  492. jac = fun.derivative
  493. elif (jac in FD_METHODS and
  494. meth in ['trust-constr', 'bfgs', 'cg', 'l-bfgs-b', 'tnc', 'slsqp']):
  495. # finite differences with relative step
  496. pass
  497. elif meth in ['trust-constr']:
  498. # default jac calculation for this method
  499. jac = '2-point'
  500. elif jac is None or bool(jac) is False:
  501. # this will cause e.g. LBFGS to use forward difference, absolute step
  502. jac = None
  503. else:
  504. # default if jac option is not understood
  505. jac = None
  506. # set default tolerances
  507. if tol is not None:
  508. options = dict(options)
  509. if meth == 'nelder-mead':
  510. options.setdefault('xatol', tol)
  511. options.setdefault('fatol', tol)
  512. if meth in ('newton-cg', 'powell', 'tnc'):
  513. options.setdefault('xtol', tol)
  514. if meth in ('powell', 'l-bfgs-b', 'tnc', 'slsqp'):
  515. options.setdefault('ftol', tol)
  516. if meth in ('bfgs', 'cg', 'l-bfgs-b', 'tnc', 'dogleg',
  517. 'trust-ncg', 'trust-exact', 'trust-krylov'):
  518. options.setdefault('gtol', tol)
  519. if meth in ('cobyla', '_custom'):
  520. options.setdefault('tol', tol)
  521. if meth == 'trust-constr':
  522. options.setdefault('xtol', tol)
  523. options.setdefault('gtol', tol)
  524. options.setdefault('barrier_tol', tol)
  525. if meth == '_custom':
  526. # custom method called before bounds and constraints are 'standardised'
  527. # custom method should be able to accept whatever bounds/constraints
  528. # are provided to it.
  529. return method(fun, x0, args=args, jac=jac, hess=hess, hessp=hessp,
  530. bounds=bounds, constraints=constraints,
  531. callback=callback, **options)
  532. constraints = standardize_constraints(constraints, x0, meth)
  533. remove_vars = False
  534. if bounds is not None:
  535. if meth in {"tnc", "slsqp", "l-bfgs-b"}:
  536. # These methods can't take the finite-difference derivatives they
  537. # need when a variable is fixed by the bounds. To avoid this issue,
  538. # remove fixed variables from the problem.
  539. # NOTE: if this list is expanded, then be sure to update the
  540. # accompanying tests and test_optimize.eb_data. Consider also if
  541. # default OptimizeResult will need updating.
  542. # convert to new-style bounds so we only have to consider one case
  543. bounds = standardize_bounds(bounds, x0, 'new')
  544. # determine whether any variables are fixed
  545. i_fixed = (bounds.lb == bounds.ub)
  546. if np.all(i_fixed):
  547. # all the parameters are fixed, a minimizer is not able to do
  548. # anything
  549. return _optimize_result_for_equal_bounds(
  550. fun, bounds, meth, args=args, constraints=constraints
  551. )
  552. # determine whether finite differences are needed for any grad/jac
  553. fd_needed = (not callable(jac))
  554. for con in constraints:
  555. if not callable(con.get('jac', None)):
  556. fd_needed = True
  557. # If finite differences are ever used, remove all fixed variables
  558. # Always remove fixed variables for TNC; see gh-14565
  559. remove_vars = i_fixed.any() and (fd_needed or meth == "tnc")
  560. if remove_vars:
  561. x_fixed = (bounds.lb)[i_fixed]
  562. x0 = x0[~i_fixed]
  563. bounds = _remove_from_bounds(bounds, i_fixed)
  564. fun = _remove_from_func(fun, i_fixed, x_fixed)
  565. if callable(callback):
  566. callback = _remove_from_func(callback, i_fixed, x_fixed)
  567. if callable(jac):
  568. jac = _remove_from_func(jac, i_fixed, x_fixed, remove=1)
  569. # make a copy of the constraints so the user's version doesn't
  570. # get changed. (Shallow copy is ok)
  571. constraints = [con.copy() for con in constraints]
  572. for con in constraints: # yes, guaranteed to be a list
  573. con['fun'] = _remove_from_func(con['fun'], i_fixed,
  574. x_fixed, min_dim=1,
  575. remove=0)
  576. if callable(con.get('jac', None)):
  577. con['jac'] = _remove_from_func(con['jac'], i_fixed,
  578. x_fixed, min_dim=2,
  579. remove=1)
  580. bounds = standardize_bounds(bounds, x0, meth)
  581. if meth == 'nelder-mead':
  582. res = _minimize_neldermead(fun, x0, args, callback, bounds=bounds,
  583. **options)
  584. elif meth == 'powell':
  585. res = _minimize_powell(fun, x0, args, callback, bounds, **options)
  586. elif meth == 'cg':
  587. res = _minimize_cg(fun, x0, args, jac, callback, **options)
  588. elif meth == 'bfgs':
  589. res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
  590. elif meth == 'newton-cg':
  591. res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
  592. **options)
  593. elif meth == 'l-bfgs-b':
  594. res = _minimize_lbfgsb(fun, x0, args, jac, bounds,
  595. callback=callback, **options)
  596. elif meth == 'tnc':
  597. res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,
  598. **options)
  599. elif meth == 'cobyla':
  600. res = _minimize_cobyla(fun, x0, args, constraints, callback=callback,
  601. **options)
  602. elif meth == 'slsqp':
  603. res = _minimize_slsqp(fun, x0, args, jac, bounds,
  604. constraints, callback=callback, **options)
  605. elif meth == 'trust-constr':
  606. res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
  607. bounds, constraints,
  608. callback=callback, **options)
  609. elif meth == 'dogleg':
  610. res = _minimize_dogleg(fun, x0, args, jac, hess,
  611. callback=callback, **options)
  612. elif meth == 'trust-ncg':
  613. res = _minimize_trust_ncg(fun, x0, args, jac, hess, hessp,
  614. callback=callback, **options)
  615. elif meth == 'trust-krylov':
  616. res = _minimize_trust_krylov(fun, x0, args, jac, hess, hessp,
  617. callback=callback, **options)
  618. elif meth == 'trust-exact':
  619. res = _minimize_trustregion_exact(fun, x0, args, jac, hess,
  620. callback=callback, **options)
  621. else:
  622. raise ValueError('Unknown solver %s' % method)
  623. if remove_vars:
  624. res.x = _add_to_array(res.x, i_fixed, x_fixed)
  625. res.jac = _add_to_array(res.jac, i_fixed, np.nan)
  626. if "hess_inv" in res:
  627. res.hess_inv = None # unknown
  628. return res
  629. def minimize_scalar(fun, bracket=None, bounds=None, args=(),
  630. method=None, tol=None, options=None):
  631. """Minimization of scalar function of one variable.
  632. Parameters
  633. ----------
  634. fun : callable
  635. Objective function.
  636. Scalar function, must return a scalar.
  637. bracket : sequence, optional
  638. For methods 'brent' and 'golden', `bracket` defines the bracketing
  639. interval and can either have three items ``(a, b, c)`` so that
  640. ``a < b < c`` and ``fun(b) < fun(a), fun(c)`` or two items ``a`` and
  641. ``c`` which are assumed to be a starting interval for a downhill
  642. bracket search (see `bracket`); it doesn't always mean that the
  643. obtained solution will satisfy ``a <= x <= c``.
  644. bounds : sequence, optional
  645. For method 'bounded', `bounds` is mandatory and must have two finite
  646. items corresponding to the optimization bounds.
  647. args : tuple, optional
  648. Extra arguments passed to the objective function.
  649. method : str or callable, optional
  650. Type of solver. Should be one of:
  651. - :ref:`Brent <optimize.minimize_scalar-brent>`
  652. - :ref:`Bounded <optimize.minimize_scalar-bounded>`
  653. - :ref:`Golden <optimize.minimize_scalar-golden>`
  654. - custom - a callable object (added in version 0.14.0), see below
  655. Default is "Bounded" if bounds are provided and "Brent" otherwise.
  656. See the 'Notes' section for details of each solver.
  657. tol : float, optional
  658. Tolerance for termination. For detailed control, use solver-specific
  659. options.
  660. options : dict, optional
  661. A dictionary of solver options.
  662. maxiter : int
  663. Maximum number of iterations to perform.
  664. disp : bool
  665. Set to True to print convergence messages.
  666. See :func:`show_options()` for solver-specific options.
  667. Returns
  668. -------
  669. res : OptimizeResult
  670. The optimization result represented as a ``OptimizeResult`` object.
  671. Important attributes are: ``x`` the solution array, ``success`` a
  672. Boolean flag indicating if the optimizer exited successfully and
  673. ``message`` which describes the cause of the termination. See
  674. `OptimizeResult` for a description of other attributes.
  675. See also
  676. --------
  677. minimize : Interface to minimization algorithms for scalar multivariate
  678. functions
  679. show_options : Additional options accepted by the solvers
  680. Notes
  681. -----
  682. This section describes the available solvers that can be selected by the
  683. 'method' parameter. The default method is the ``"Bounded"`` Brent method if
  684. `bounds` are passed and unbounded ``"Brent"`` otherwise.
  685. Method :ref:`Brent <optimize.minimize_scalar-brent>` uses Brent's
  686. algorithm [1]_ to find a local minimum. The algorithm uses inverse
  687. parabolic interpolation when possible to speed up convergence of
  688. the golden section method.
  689. Method :ref:`Golden <optimize.minimize_scalar-golden>` uses the
  690. golden section search technique [1]_. It uses analog of the bisection
  691. method to decrease the bracketed interval. It is usually
  692. preferable to use the *Brent* method.
  693. Method :ref:`Bounded <optimize.minimize_scalar-bounded>` can
  694. perform bounded minimization [2]_ [3]_. It uses the Brent method to find a
  695. local minimum in the interval x1 < xopt < x2.
  696. **Custom minimizers**
  697. It may be useful to pass a custom minimization method, for example
  698. when using some library frontend to minimize_scalar. You can simply
  699. pass a callable as the ``method`` parameter.
  700. The callable is called as ``method(fun, args, **kwargs, **options)``
  701. where ``kwargs`` corresponds to any other parameters passed to `minimize`
  702. (such as `bracket`, `tol`, etc.), except the `options` dict, which has
  703. its contents also passed as `method` parameters pair by pair. The method
  704. shall return an `OptimizeResult` object.
  705. The provided `method` callable must be able to accept (and possibly ignore)
  706. arbitrary parameters; the set of parameters accepted by `minimize` may
  707. expand in future versions and then these parameters will be passed to
  708. the method. You can find an example in the scipy.optimize tutorial.
  709. .. versionadded:: 0.11.0
  710. References
  711. ----------
  712. .. [1] Press, W., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery.
  713. Numerical Recipes in C. Cambridge University Press.
  714. .. [2] Forsythe, G.E., M. A. Malcolm, and C. B. Moler. "Computer Methods
  715. for Mathematical Computations." Prentice-Hall Series in Automatic
  716. Computation 259 (1977).
  717. .. [3] Brent, Richard P. Algorithms for Minimization Without Derivatives.
  718. Courier Corporation, 2013.
  719. Examples
  720. --------
  721. Consider the problem of minimizing the following function.
  722. >>> def f(x):
  723. ... return (x - 2) * x * (x + 2)**2
  724. Using the *Brent* method, we find the local minimum as:
  725. >>> from scipy.optimize import minimize_scalar
  726. >>> res = minimize_scalar(f)
  727. >>> res.fun
  728. -9.9149495908
  729. The minimizer is:
  730. >>> res.x
  731. 1.28077640403
  732. Using the *Bounded* method, we find a local minimum with specified
  733. bounds as:
  734. >>> res = minimize_scalar(f, bounds=(-3, -1), method='bounded')
  735. >>> res.fun # minimum
  736. 3.28365179850e-13
  737. >>> res.x # minimizer
  738. -2.0000002026
  739. """
  740. if not isinstance(args, tuple):
  741. args = (args,)
  742. if callable(method):
  743. meth = "_custom"
  744. elif method is None:
  745. meth = 'brent' if bounds is None else 'bounded'
  746. else:
  747. meth = method.lower()
  748. if options is None:
  749. options = {}
  750. if bounds is not None and meth in {'brent', 'golden'}:
  751. message = f"Use of `bounds` is incompatible with 'method={method}'."
  752. raise ValueError(message)
  753. if tol is not None:
  754. options = dict(options)
  755. if meth == 'bounded' and 'xatol' not in options:
  756. warn("Method 'bounded' does not support relative tolerance in x; "
  757. "defaulting to absolute tolerance.", RuntimeWarning)
  758. options['xatol'] = tol
  759. elif meth == '_custom':
  760. options.setdefault('tol', tol)
  761. else:
  762. options.setdefault('xtol', tol)
  763. # replace boolean "disp" option, if specified, by an integer value.
  764. disp = options.get('disp')
  765. if isinstance(disp, bool):
  766. options['disp'] = 2 * int(disp)
  767. if meth == '_custom':
  768. return method(fun, args=args, bracket=bracket, bounds=bounds, **options)
  769. elif meth == 'brent':
  770. return _minimize_scalar_brent(fun, bracket, args, **options)
  771. elif meth == 'bounded':
  772. if bounds is None:
  773. raise ValueError('The `bounds` parameter is mandatory for '
  774. 'method `bounded`.')
  775. return _minimize_scalar_bounded(fun, bounds, args, **options)
  776. elif meth == 'golden':
  777. return _minimize_scalar_golden(fun, bracket, args, **options)
  778. else:
  779. raise ValueError('Unknown solver %s' % method)
  780. def _remove_from_bounds(bounds, i_fixed):
  781. """Removes fixed variables from a `Bounds` instance"""
  782. lb = bounds.lb[~i_fixed]
  783. ub = bounds.ub[~i_fixed]
  784. return Bounds(lb, ub) # don't mutate original Bounds object
  785. def _remove_from_func(fun_in, i_fixed, x_fixed, min_dim=None, remove=0):
  786. """Wraps a function such that fixed variables need not be passed in"""
  787. def fun_out(x_in, *args, **kwargs):
  788. x_out = np.zeros_like(i_fixed, dtype=x_in.dtype)
  789. x_out[i_fixed] = x_fixed
  790. x_out[~i_fixed] = x_in
  791. y_out = fun_in(x_out, *args, **kwargs)
  792. y_out = np.array(y_out)
  793. if min_dim == 1:
  794. y_out = np.atleast_1d(y_out)
  795. elif min_dim == 2:
  796. y_out = np.atleast_2d(y_out)
  797. if remove == 1:
  798. y_out = y_out[..., ~i_fixed]
  799. elif remove == 2:
  800. y_out = y_out[~i_fixed, ~i_fixed]
  801. return y_out
  802. return fun_out
  803. def _add_to_array(x_in, i_fixed, x_fixed):
  804. """Adds fixed variables back to an array"""
  805. i_free = ~i_fixed
  806. if x_in.ndim == 2:
  807. i_free = i_free[:, None] @ i_free[None, :]
  808. x_out = np.zeros_like(i_free, dtype=x_in.dtype)
  809. x_out[~i_free] = x_fixed
  810. x_out[i_free] = x_in.ravel()
  811. return x_out
  812. def standardize_bounds(bounds, x0, meth):
  813. """Converts bounds to the form required by the solver."""
  814. if meth in {'trust-constr', 'powell', 'nelder-mead', 'new'}:
  815. if not isinstance(bounds, Bounds):
  816. lb, ub = old_bound_to_new(bounds)
  817. bounds = Bounds(lb, ub)
  818. elif meth in ('l-bfgs-b', 'tnc', 'slsqp', 'old'):
  819. if isinstance(bounds, Bounds):
  820. bounds = new_bounds_to_old(bounds.lb, bounds.ub, x0.shape[0])
  821. return bounds
  822. def standardize_constraints(constraints, x0, meth):
  823. """Converts constraints to the form required by the solver."""
  824. all_constraint_types = (NonlinearConstraint, LinearConstraint, dict)
  825. new_constraint_types = all_constraint_types[:-1]
  826. if constraints is None:
  827. constraints = []
  828. elif isinstance(constraints, all_constraint_types):
  829. constraints = [constraints]
  830. else:
  831. constraints = list(constraints) # ensure it's a mutable sequence
  832. if meth in ['trust-constr', 'new']:
  833. for i, con in enumerate(constraints):
  834. if not isinstance(con, new_constraint_types):
  835. constraints[i] = old_constraint_to_new(i, con)
  836. else:
  837. # iterate over copy, changing original
  838. for i, con in enumerate(list(constraints)):
  839. if isinstance(con, new_constraint_types):
  840. old_constraints = new_constraint_to_old(con, x0)
  841. constraints[i] = old_constraints[0]
  842. constraints.extend(old_constraints[1:]) # appends 1 if present
  843. return constraints
  844. def _optimize_result_for_equal_bounds(
  845. fun, bounds, method, args=(), constraints=()
  846. ):
  847. """
  848. Provides a default OptimizeResult for when a bounded minimization method
  849. has (lb == ub).all().
  850. Parameters
  851. ----------
  852. fun: callable
  853. bounds: Bounds
  854. method: str
  855. constraints: Constraint
  856. """
  857. success = True
  858. message = 'All independent variables were fixed by bounds.'
  859. # bounds is new-style
  860. x0 = bounds.lb
  861. if constraints:
  862. message = ("All independent variables were fixed by bounds at values"
  863. " that satisfy the constraints.")
  864. constraints = standardize_constraints(constraints, x0, 'new')
  865. maxcv = 0
  866. for c in constraints:
  867. pc = PreparedConstraint(c, x0)
  868. violation = pc.violation(x0)
  869. if np.sum(violation):
  870. maxcv = max(maxcv, np.max(violation))
  871. success = False
  872. message = (f"All independent variables were fixed by bounds, but "
  873. f"the independent variables do not satisfy the "
  874. f"constraints exactly. (Maximum violation: {maxcv}).")
  875. return OptimizeResult(
  876. x=x0, fun=fun(x0, *args), success=success, message=message, nfev=1,
  877. njev=0, nhev=0,
  878. )