_cobyla_py.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296
  1. """
  2. Interface to Constrained Optimization By Linear Approximation
  3. Functions
  4. ---------
  5. .. autosummary::
  6. :toctree: generated/
  7. fmin_cobyla
  8. """
  9. import numpy as np
  10. from scipy._lib._util import wrapped_inspect_signature
  11. from ._optimize import (OptimizeResult, _check_unknown_options,
  12. _prepare_scalar_function)
  13. from ._constraints import NonlinearConstraint
  14. __all__ = ['fmin_cobyla']
  15. def fmin_cobyla(func, x0, cons, args=(), consargs=None, rhobeg=1.0,
  16. rhoend=1e-4, maxfun=1000, disp=None, catol=2e-4,
  17. *, callback=None):
  18. """
  19. Minimize a function using the Constrained Optimization By Linear
  20. Approximation (COBYLA) method. This method uses the pure-python implementation
  21. of the algorithm from PRIMA.
  22. Parameters
  23. ----------
  24. func : callable
  25. Function to minimize. In the form func(x, \\*args).
  26. x0 : ndarray
  27. Initial guess.
  28. cons : sequence
  29. Constraint functions; must all be ``>=0`` (a single function
  30. if only 1 constraint). Each function takes the parameters `x`
  31. as its first argument, and it can return either a single number or
  32. an array or list of numbers.
  33. args : tuple, optional
  34. Extra arguments to pass to function.
  35. consargs : tuple, optional
  36. Extra arguments to pass to constraint functions (default of None means
  37. use same extra arguments as those passed to func).
  38. Use ``()`` for no extra arguments.
  39. rhobeg : float, optional
  40. Reasonable initial changes to the variables.
  41. rhoend : float, optional
  42. Final accuracy in the optimization (not precisely guaranteed). This
  43. is a lower bound on the size of the trust region.
  44. disp : {0, 1, 2, 3}, optional
  45. Controls the frequency of output; 0 implies no output.
  46. maxfun : int, optional
  47. Maximum number of function evaluations.
  48. catol : float, optional
  49. Absolute tolerance for constraint violations.
  50. callback : callable, optional
  51. Called after each iteration, as ``callback(x)``, where ``x`` is the
  52. current parameter vector.
  53. Returns
  54. -------
  55. x : ndarray
  56. The argument that minimises `f`.
  57. See also
  58. --------
  59. minimize: Interface to minimization algorithms for multivariate
  60. functions. See the 'COBYLA' `method` in particular.
  61. Notes
  62. -----
  63. This algorithm is based on linear approximations to the objective
  64. function and each constraint. We briefly describe the algorithm.
  65. Suppose the function is being minimized over k variables. At the
  66. jth iteration the algorithm has k+1 points v_1, ..., v_(k+1),
  67. an approximate solution x_j, and a radius RHO_j.
  68. (i.e., linear plus a constant) approximations to the objective
  69. function and constraint functions such that their function values
  70. agree with the linear approximation on the k+1 points v_1,.., v_(k+1).
  71. This gives a linear program to solve (where the linear approximations
  72. of the constraint functions are constrained to be non-negative).
  73. However, the linear approximations are likely only good
  74. approximations near the current simplex, so the linear program is
  75. given the further requirement that the solution, which
  76. will become x_(j+1), must be within RHO_j from x_j. RHO_j only
  77. decreases, never increases. The initial RHO_j is rhobeg and the
  78. final RHO_j is rhoend. In this way COBYLA's iterations behave
  79. like a trust region algorithm.
  80. Additionally, the linear program may be inconsistent, or the
  81. approximation may give poor improvement. For details about
  82. how these issues are resolved, as well as how the points v_i are
  83. updated, refer to the source code or the references below.
  84. .. versionchanged:: 1.16.0
  85. The original Powell implementation was replaced by a pure
  86. Python version from the PRIMA package, with bug fixes and
  87. improvements being made.
  88. References
  89. ----------
  90. Powell M.J.D. (1994), "A direct search optimization method that models
  91. the objective and constraint functions by linear interpolation.", in
  92. Advances in Optimization and Numerical Analysis, eds. S. Gomez and
  93. J-P Hennart, Kluwer Academic (Dordrecht), pp. 51-67
  94. Powell M.J.D. (1998), "Direct search algorithms for optimization
  95. calculations", Acta Numerica 7, 287-336
  96. Powell M.J.D. (2007), "A view of algorithms for optimization without
  97. derivatives", Cambridge University Technical Report DAMTP 2007/NA03
  98. Zhang Z. (2023), "PRIMA: Reference Implementation for Powell's Methods with
  99. Modernization and Amelioration", https://www.libprima.net,
  100. :doi:`10.5281/zenodo.8052654`
  101. Examples
  102. --------
  103. Minimize the objective function f(x,y) = x*y subject
  104. to the constraints x**2 + y**2 < 1 and y > 0::
  105. >>> def objective(x):
  106. ... return x[0]*x[1]
  107. ...
  108. >>> def constr1(x):
  109. ... return 1 - (x[0]**2 + x[1]**2)
  110. ...
  111. >>> def constr2(x):
  112. ... return x[1]
  113. ...
  114. >>> from scipy.optimize import fmin_cobyla
  115. >>> fmin_cobyla(objective, [0.0, 0.1], [constr1, constr2], rhoend=1e-7)
  116. array([-0.70710685, 0.70710671])
  117. The exact solution is (-sqrt(2)/2, sqrt(2)/2).
  118. """
  119. err = "cons must be a sequence of callable functions or a single"\
  120. " callable function."
  121. try:
  122. len(cons)
  123. except TypeError as e:
  124. if callable(cons):
  125. cons = [cons]
  126. else:
  127. raise TypeError(err) from e
  128. else:
  129. for thisfunc in cons:
  130. if not callable(thisfunc):
  131. raise TypeError(err)
  132. if consargs is None:
  133. consargs = args
  134. # build constraints
  135. nlcs = []
  136. for con in cons:
  137. # Use default argument, otherwise the last `con` is captured by all wrapped_con
  138. def wrapped_con(x, confunc=con):
  139. return confunc(x, *consargs)
  140. nlcs.append(NonlinearConstraint(wrapped_con, 0, np.inf))
  141. # options
  142. opts = {'rhobeg': rhobeg,
  143. 'tol': rhoend,
  144. 'disp': disp,
  145. 'maxiter': maxfun,
  146. 'catol': catol,
  147. 'callback': callback}
  148. sol = _minimize_cobyla(func, x0, args, constraints=nlcs,
  149. **opts)
  150. if disp and not sol['success']:
  151. print(f"COBYLA failed to find a solution: {sol.message}")
  152. return sol['x']
  153. def _minimize_cobyla(fun, x0, args=(), constraints=(),
  154. rhobeg=1.0, tol=1e-4, maxiter=1000,
  155. disp=0, catol=None, f_target=-np.inf,
  156. callback=None, bounds=None, **unknown_options):
  157. """
  158. Minimize a scalar function of one or more variables using the
  159. Constrained Optimization BY Linear Approximation (COBYLA) algorithm.
  160. This method uses the pure-python implementation of the algorithm from PRIMA.
  161. Options
  162. -------
  163. rhobeg : float
  164. Reasonable initial changes to the variables.
  165. tol : float
  166. Final accuracy in the optimization (not precisely guaranteed).
  167. This is a lower bound on the size of the trust region.
  168. disp : int
  169. Controls the frequency of output:
  170. 0. (default) There will be no printing
  171. 1. A message will be printed to the screen at the end of iteration, showing
  172. the best vector of variables found and its objective function value
  173. 2. in addition to 1, each new value of RHO is printed to the screen,
  174. with the best vector of variables so far and its objective function
  175. value.
  176. 3. in addition to 2, each function evaluation with its variables will
  177. be printed to the screen.
  178. maxiter : int
  179. Maximum number of function evaluations.
  180. catol : float
  181. Tolerance (absolute) for constraint violations
  182. f_target : float
  183. Stop if the objective function is less than `f_target`.
  184. .. versionchanged:: 1.16.0
  185. The original Powell implementation was replaced by a pure
  186. Python version from the PRIMA package, with bug fixes and
  187. improvements being made.
  188. References
  189. ----------
  190. Zhang Z. (2023), "PRIMA: Reference Implementation for Powell's Methods with
  191. Modernization and Amelioration", https://www.libprima.net,
  192. :doi:`10.5281/zenodo.8052654`
  193. """
  194. from .._lib.pyprima import minimize
  195. from .._lib.pyprima.common.infos import SMALL_TR_RADIUS, FTARGET_ACHIEVED
  196. from .._lib.pyprima.common.message import get_info_string
  197. _check_unknown_options(unknown_options)
  198. rhoend = tol
  199. iprint = disp if disp is not None else 0
  200. if iprint != 0 and iprint != 1 and iprint != 2 and iprint != 3:
  201. raise ValueError(f'disp argument to minimize must be 0, 1, 2, or 3,\
  202. received {iprint}')
  203. # create the ScalarFunction, cobyla doesn't require derivative function
  204. def _jac(x, *args):
  205. return None
  206. sf = _prepare_scalar_function(fun, x0, args=args, jac=_jac)
  207. if callback is not None:
  208. sig = wrapped_inspect_signature(callback)
  209. if set(sig.parameters) == {"intermediate_result"}:
  210. def wrapped_callback_intermediate(x, f, nf, tr, cstrv, nlconstrlist):
  211. intermediate_result = OptimizeResult(x=np.copy(x), fun=f, nfev=nf,
  212. nit=tr, maxcv=cstrv)
  213. callback(intermediate_result=intermediate_result)
  214. else:
  215. def wrapped_callback_intermediate(x, f, nf, tr, cstrv, nlconstrlist):
  216. callback(np.copy(x))
  217. def wrapped_callback(x, f, nf, tr, cstrv, nlconstrlist):
  218. try:
  219. wrapped_callback_intermediate(x, f, nf, tr, cstrv, nlconstrlist)
  220. return False
  221. except StopIteration:
  222. return True
  223. else:
  224. wrapped_callback = None
  225. ctol = catol if catol is not None else np.sqrt(np.finfo(float).eps)
  226. options = {
  227. 'rhobeg': rhobeg,
  228. 'rhoend': rhoend,
  229. 'maxfev': maxiter,
  230. 'iprint': iprint,
  231. 'ctol': ctol,
  232. 'ftarget': f_target,
  233. }
  234. result = minimize(sf.fun, x0, method='cobyla', bounds=bounds,
  235. constraints=constraints, callback=wrapped_callback,
  236. options=options)
  237. if result.cstrv > ctol:
  238. success = False
  239. message = ('Did not converge to a solution satisfying the constraints. See '
  240. '`maxcv` for the magnitude of the violation.')
  241. else:
  242. success = result.info == SMALL_TR_RADIUS or result.info == FTARGET_ACHIEVED
  243. message = get_info_string('COBYLA', result.info)
  244. return OptimizeResult(x=result.x,
  245. status=result.info,
  246. success=success,
  247. message=message,
  248. nfev=result.nf,
  249. fun=result.f,
  250. maxcv=result.cstrv)