ENH: add cobyqa to scipy.optimize. (scipy#19994)

Co-authored-by: Tom M. Ragonneau <t.ragonneau@gmail.com>
Co-authored-by: Ralf Gommers <ralf.gommers@gmail.com>
Co-authored-by: Lucas Colley <lucas.colley8@gmail.com>

Author: 3 people

.gitmodules

@@ -23,3 +23,6 @@
 [submodule "scipy/_lib/pocketfft"]
 	path = scipy/_lib/pocketfft
 	url = https://github.com/scipy/pocketfft
+[submodule "scipy/_lib/cobyqa"]
+	path = scipy/_lib/cobyqa
+	url = https://github.com/cobyqa/cobyqa.git

benchmarks/benchmarks/optimize.py

@@ -273,8 +273,8 @@ def bench_run(self, x0, methods=None, **minimizer_kwargs):
 
         # L-BFGS-B, BFGS, trust-constr, SLSQP can use gradients, but examine
         # performance when numerical differentiation is used.
-        fonly_methods = ["COBYLA", 'Powell', 'nelder-mead', 'L-BFGS-B', 'BFGS',
-                         'trust-constr', 'SLSQP']
+        fonly_methods = ["COBYLA", 'COBYQA', 'Powell', 'nelder-mead',
+                         'L-BFGS-B', 'BFGS', 'trust-constr', 'SLSQP']
         for method in fonly_methods:
             if method not in methods:
                 continue
@@ -316,7 +316,7 @@ class BenchSmoothUnbounded(Benchmark):
         ['rosenbrock_slow', 'rosenbrock_nograd', 'rosenbrock', 'rosenbrock_tight',
          'simple_quadratic', 'asymmetric_quadratic',
          'sin_1d', 'booth', 'beale', 'LJ'],
-        ["COBYLA", 'Powell', 'nelder-mead',
+        ["COBYLA", 'COBYQA', 'Powell', 'nelder-mead',
          'L-BFGS-B', 'BFGS', 'CG', 'TNC', 'SLSQP',
          "Newton-CG", 'dogleg', 'trust-ncg', 'trust-exact',
          'trust-krylov', 'trust-constr'],
@@ -599,7 +599,8 @@ class BenchDFO(Benchmark):
 
     params = [
         list(range(53)),  # adjust which problems to solve
-        ["COBYLA", "SLSQP", "Powell", "nelder-mead", "L-BFGS-B", "BFGS",
+        ["COBYLA", "COBYQA", "SLSQP", "Powell", "nelder-mead", "L-BFGS-B",
+         "BFGS",
          "trust-constr"],  # note: methods must also be listed in bench_run
         ["mean_nfev", "min_obj"],  # defined in average_results
     ]

doc/source/reference/optimize.minimize-cobyqa.rst

@@ -0,0 +1,8 @@
+.. _optimize.minimize-cobyqa:
+
+minimize(method='COBYQA')
+-------------------------
+
+.. scipy-optimize:function:: scipy.optimize.minimize
+   :impl: scipy.optimize._cobyqa_py._minimize_cobyqa
+   :method: COBYQA

doc/source/tutorial/optimize.rst

@@ -541,8 +541,8 @@ Constrained minimization of multivariate scalar functions (:func:`minimize`)
 ----------------------------------------------------------------------------
 
 The :func:`minimize` function provides algorithms for constrained minimization,
-namely ``'trust-constr'`` ,  ``'SLSQP'`` and ``'COBYLA'``. They require the constraints
-to be defined using slightly different structures. The method ``'trust-constr'`` requires
+namely ``'trust-constr'`` ,  ``'SLSQP'``, ``'COBYLA'``, and ``'COBYQA'``. They require the constraints
+to be defined using slightly different structures. The methods ``'trust-constr'`` and ``'COBYQA'`` require
 the  constraints to be defined as a sequence of objects :func:`LinearConstraint` and
 :func:`NonlinearConstraint`. Methods ``'SLSQP'`` and ``'COBYLA'``, on the other hand,
 require constraints to be defined  as a sequence of dictionaries, with keys

pytest.ini

@@ -1,6 +1,6 @@
 [pytest]
 addopts = -l
-norecursedirs = doc tools scipy/_lib/array_api_compat scipy/_lib/highs
+norecursedirs = doc tools scipy/_lib/array_api_compat scipy/_lib/cobyqa scipy/_lib/highs
 junit_family=xunit2
 
 filterwarnings =

scipy/_lib/cobyqa

@@ -14,6 +14,9 @@ endif
 if not fs.exists('pocketfft/README.md')
   error('Missing the `pocketfft` submodule! Run `git submodule update --init` to fix this.')
 endif
+if not fs.exists('cobyqa/README.rst')
+  error('Missing the `cobyqa` submodule! Run `git submodule update --init` to fix this.')
+endif
 
 _lib_pxd = [
   fs.copyfile('__init__.py'),
@@ -199,5 +202,36 @@ py3.install_sources(
   subdir: 'scipy/_lib/array_api_compat/torch',
 )
 
+py3.install_sources(
+  [
+    'cobyqa/cobyqa/__init__.py',
+    'cobyqa/cobyqa/framework.py',
+    'cobyqa/cobyqa/main.py',
+    'cobyqa/cobyqa/models.py',
+    'cobyqa/cobyqa/problem.py',
+    'cobyqa/cobyqa/settings.py',
+  ],
+  subdir: 'scipy/_lib/cobyqa',
+)
+
+py3.install_sources(
+  [
+    'cobyqa/cobyqa/subsolvers/__init__.py',
+    'cobyqa/cobyqa/subsolvers/geometry.py',
+    'cobyqa/cobyqa/subsolvers/optim.py',
+  ],
+  subdir: 'scipy/_lib/cobyqa/subsolvers',
+)
+
+py3.install_sources(
+  [
+    'cobyqa/cobyqa/utils/__init__.py',
+    'cobyqa/cobyqa/utils/exceptions.py',
+    'cobyqa/cobyqa/utils/math.py',
+    'cobyqa/cobyqa/utils/versions.py',
+  ],
+  subdir: 'scipy/_lib/cobyqa/utils',
+)
+
 subdir('_uarray')
 subdir('tests')

scipy/_lib/meson.build

@@ -65,6 +65,7 @@
    optimize.minimize-lbfgsb
    optimize.minimize-tnc
    optimize.minimize-cobyla
+   optimize.minimize-cobyqa
    optimize.minimize-slsqp
    optimize.minimize-trustconstr
    optimize.minimize-dogleg

scipy/optimize/__init__.py

@@ -0,0 +1,57 @@
+import numpy as np
+
+from ._optimize import _check_unknown_options
+
+
+def _minimize_cobyqa(fun, x0, args=(), bounds=None, constraints=(),
+                     callback=None, disp=False, maxfev=None, maxiter=None,
+                     f_target=-np.inf, feasibility_tol=1e-8,
+                     initial_tr_radius=1.0, final_tr_radius=1e-6, scale=False,
+                     **unknown_options):
+    """
+    Minimize a scalar function of one or more variables using the
+    Constrained Optimization BY Quadratic Approximations (COBYQA) algorithm.
+
+    .. versionadded:: 1.14.0
+
+    Options
+    -------
+    disp : bool
+        Set to True to print information about the optimization procedure.
+    maxfev : int
+        Maximum number of function evaluations.
+    maxiter : int
+        Maximum number of iterations.
+    f_target : float
+        Target value for the objective function. The optimization procedure is
+        terminated when the objective function value of a feasible point (see
+        `feasibility_tol` below) is less than or equal to this target.
+    feasibility_tol : float
+        Absolute tolerance for the constraint violation.
+    initial_tr_radius : float
+        Initial trust-region radius. Typically, this value should be in the
+        order of one tenth of the greatest expected change to the variables.
+    final_tr_radius : float
+        Final trust-region radius. It should indicate the accuracy required in
+        the final values of the variables. If provided, this option overrides
+        the value of `tol` in the `minimize` function.
+    scale : bool
+        Set to True to scale the variables according to the bounds. If True and
+        if all the lower and upper bounds are finite, the variables are scaled
+        to be within the range :math:`[-1, 1]`. If any of the lower or upper
+        bounds is infinite, the variables are not scaled.
+    """
+    from .._lib.cobyqa import minimize  # import here to avoid circular imports
+
+    _check_unknown_options(unknown_options)
+    options = {
+        'disp': bool(disp),
+        'maxfev': int(maxfev) if maxfev is not None else 500 * len(x0),
+        'maxiter': int(maxiter) if maxiter is not None else 1000 * len(x0),
+        'target': float(f_target),
+        'feasibility_tol': float(feasibility_tol),
+        'radius_init': float(initial_tr_radius),
+        'radius_final': float(final_tr_radius),
+        'scale': bool(scale),
+    }
+    return minimize(fun, x0, args, bounds, constraints, callback, options)

scipy/optimize/_cobyqa_py.py

@@ -30,6 +30,7 @@
 from ._lbfgsb_py import _minimize_lbfgsb
 from ._tnc import _minimize_tnc
 from ._cobyla_py import _minimize_cobyla
+from ._cobyqa_py import _minimize_cobyqa
 from ._slsqp_py import _minimize_slsqp
 from ._constraints import (old_bound_to_new, new_bounds_to_old,
                            old_constraint_to_new, new_constraint_to_old,
@@ -38,13 +39,14 @@
 from ._differentiable_functions import FD_METHODS
 
 MINIMIZE_METHODS = ['nelder-mead', 'powell', 'cg', 'bfgs', 'newton-cg',
-                    'l-bfgs-b', 'tnc', 'cobyla', 'slsqp', 'trust-constr',
-                    'dogleg', 'trust-ncg', 'trust-exact', 'trust-krylov']
+                    'l-bfgs-b', 'tnc', 'cobyla', 'cobyqa', 'slsqp',
+                    'trust-constr', 'dogleg', 'trust-ncg', 'trust-exact',
+                    'trust-krylov']
 
 # These methods support the new callback interface (passed an OptimizeResult)
 MINIMIZE_METHODS_NEW_CB = ['nelder-mead', 'powell', 'cg', 'bfgs', 'newton-cg',
                            'l-bfgs-b', 'trust-constr', 'dogleg', 'trust-ncg',
-                           'trust-exact', 'trust-krylov']
+                           'trust-exact', 'trust-krylov', 'cobyqa']
 
 MINIMIZE_SCALAR_METHODS = ['brent', 'bounded', 'golden']
 
@@ -80,6 +82,7 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
             - 'L-BFGS-B'    :ref:`(see here) <optimize.minimize-lbfgsb>`
             - 'TNC'         :ref:`(see here) <optimize.minimize-tnc>`
             - 'COBYLA'      :ref:`(see here) <optimize.minimize-cobyla>`
+            - 'COBYQA'      :ref:`(see here) <optimize.minimize-cobyqa>`
             - 'SLSQP'       :ref:`(see here) <optimize.minimize-slsqp>`
             - 'trust-constr':ref:`(see here) <optimize.minimize-trustconstr>`
             - 'dogleg'      :ref:`(see here) <optimize.minimize-dogleg>`
@@ -146,15 +149,15 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
         parameters.
     bounds : sequence or `Bounds`, optional
         Bounds on variables for Nelder-Mead, L-BFGS-B, TNC, SLSQP, Powell,
-        trust-constr, and COBYLA methods. There are two ways to specify the
-        bounds:
+        trust-constr, COBYLA, and COBYQA methods. There are two ways to specify
+        the bounds:
 
             1. Instance of `Bounds` class.
             2. Sequence of ``(min, max)`` pairs for each element in `x`. None
                is used to specify no bound.
 
     constraints : {Constraint, dict} or List of {Constraint, dict}, optional
-        Constraints definition. Only for COBYLA, SLSQP and trust-constr.
+        Constraints definition. Only for COBYLA, COBYQA, SLSQP and trust-constr.
 
         Constraints for 'trust-constr' are defined as a single object or a
         list of objects specifying constraints to the optimization problem.
@@ -178,6 +181,8 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
         Equality constraint means that the constraint function result is to
         be zero whereas inequality means that it is to be non-negative.
         Note that COBYLA only supports inequality constraints.
+
+        Constraints for COBYQA are defined as any of the above.
     tol : float, optional
         Tolerance for termination. When `tol` is specified, the selected
         minimization algorithm sets some relevant solver-specific tolerance(s)
@@ -335,6 +340,13 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
     constraints functions 'fun' may return either a single number
     or an array or list of numbers.
 
+    Method :ref:`COBYQA <optimize.minimize-cobyqa>` uses the Constrained
+    Optimization BY Quadratic Approximations (COBYQA) method [18]_. The
+    algorithm is a derivative-free trust-region SQP method based on quadratic
+    approximations to the objective function and each nonlinear constraint. The
+    bounds are treated as unrelaxable constraints, in the sense that the
+    algorithm always respects them throughout the optimization process.
+
     Method :ref:`SLSQP <optimize.minimize-slsqp>` uses Sequential
     Least SQuares Programming to minimize a function of several
     variables with any combination of bounds, equality and inequality
@@ -466,6 +478,10 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
     .. [17] Lalee, Marucha, Jorge Nocedal, and Todd Plantega. 1998. On the
         implementation of an algorithm for large-scale equality constrained
         optimization. SIAM Journal on Optimization 8.3: 682-706.
+    .. [18] Ragonneau, T. M. *Model-Based Derivative-Free Optimization Methods
+        and Software*. PhD thesis, Department of Applied Mathematics, The Hong
+        Kong Polytechnic University, Hong Kong, China, 2022. URL:
+        https://theses.lib.polyu.edu.hk/handle/200/12294.
 
     Examples
     --------
@@ -558,7 +574,7 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
         options = {}
     # check if optional parameters are supported by the selected method
     # - jac
-    if meth in ('nelder-mead', 'powell', 'cobyla') and bool(jac):
+    if meth in ('nelder-mead', 'powell', 'cobyla', 'cobyqa') and bool(jac):
         warn('Method %s does not use gradient information (jac).' % method,
              RuntimeWarning, stacklevel=2)
     # - hess
@@ -574,16 +590,17 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
              'information (hessp).' % method,
              RuntimeWarning, stacklevel=2)
     # - constraints or bounds
-    if (meth not in ('cobyla', 'slsqp', 'trust-constr', '_custom') and
+    if (meth not in ('cobyla', 'cobyqa', 'slsqp', 'trust-constr', '_custom') and
             np.any(constraints)):
         warn('Method %s cannot handle constraints.' % method,
              RuntimeWarning, stacklevel=2)
-    if meth not in ('nelder-mead', 'powell', 'l-bfgs-b', 'cobyla', 'slsqp',
-                    'tnc', 'trust-constr', '_custom') and bounds is not None:
+    if meth not in (
+            'nelder-mead', 'powell', 'l-bfgs-b', 'cobyla', 'cobyqa', 'slsqp',
+            'tnc', 'trust-constr', '_custom') and bounds is not None:
         warn('Method %s cannot handle bounds.' % method,
              RuntimeWarning, stacklevel=2)
     # - return_all
-    if (meth in ('l-bfgs-b', 'tnc', 'cobyla', 'slsqp') and
+    if (meth in ('l-bfgs-b', 'tnc', 'cobyla', 'cobyqa', 'slsqp') and
             options.get('return_all', False)):
         warn('Method %s does not support the return_all option.' % method,
              RuntimeWarning, stacklevel=2)
@@ -624,6 +641,8 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
             options.setdefault('gtol', tol)
         if meth in ('cobyla', '_custom'):
             options.setdefault('tol', tol)
+        if meth == 'cobyqa':
+            options.setdefault('final_tr_radius', tol)
         if meth == 'trust-constr':
             options.setdefault('xtol', tol)
             options.setdefault('gtol', tol)
@@ -718,6 +737,9 @@ def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
     elif meth == 'cobyla':
         res = _minimize_cobyla(fun, x0, args, constraints, callback=callback,
                                bounds=bounds, **options)
+    elif meth == 'cobyqa':
+        res = _minimize_cobyqa(fun, x0, args, bounds, constraints, callback,
+                               **options)
     elif meth == 'slsqp':
         res = _minimize_slsqp(fun, x0, args, jac, bounds,
                               constraints, callback=callback, **options)
@@ -1016,7 +1038,8 @@ def _validate_bounds(bounds, x0, meth):
 
 def standardize_bounds(bounds, x0, meth):
     """Converts bounds to the form required by the solver."""
-    if meth in {'trust-constr', 'powell', 'nelder-mead', 'cobyla', 'new'}:
+    if meth in {'trust-constr', 'powell', 'nelder-mead', 'cobyla', 'cobyqa',
+                'new'}:
         if not isinstance(bounds, Bounds):
             lb, ub = old_bound_to_new(bounds)
             bounds = Bounds(lb, ub)