Add option to use nonlinear constraints (#346)

CHANGES.rst

@@ -33,7 +33,7 @@ Breaking changes
   ``estimagic.OptimizeLogReader``.
 - Convenience functions to create namedtuples are removed from ``estimagic.utilities``.
 
-
+- :gh:`346` Add option to use nonlinear constraints (:ghuser:`timmens`)
 - :gh:`345` Moves estimation_table to new latex functionality of pandas
   (:ghuser:`mpetrosian`)
 - :gh:`344` Adds pytree support to slice_plot (:ghuser:`janosg`)

docs/source/explanations/optimization/implementation_of_constraints.rst

@@ -24,6 +24,14 @@ constraints into constrained optimizers. Reparametrization and penalties. Below
 explain what both approaches are, why we chose the reparametrization approach over
 penalties and which reparametrizations we are using for each type of constraint.
 
+.. note::
+
+    In this text we focus on constraints that can solved by estimagic via bijective and
+    differentiable transformations. General nonlinear constraints do not fall into this
+    category. If you want to use nonlinear constraints you can still do so, but
+    estimagic will simply pass the constraints to your chosen optimizer. See
+    :ref:`constraints` for more details.
+
 
 Possible approaches
 -------------------
@@ -32,9 +40,9 @@ Possible approaches
 Reparametrizations
 ~~~~~~~~~~~~~~~~~~
 
-In the reparametrization approach need to find an invertible mapping :math:`g` such as
-well as two :math:`k'` dimensional vectors :math:`l` and :math:`u` such that:
-
+In the reparametrization approach we need to find an invertible mapping
+:math:`g : \mathbb{R}^{k'} \to \mathbb{R}^k`, and two new bounds :math:`l'` and
+:math:`u'` such that:
 
 .. math::
 
@@ -175,7 +183,7 @@ A suitable choice of :math:`\mathbf{\tilde{X}}` and :math:`\mathbf{M}` are:
 .. math::
 
     \mathbf{\tilde{X}} \equiv \{(\tilde{x}_1, \tilde{x}_2)^T \mid \mathbf{\tilde{x}}_1
-    \in \mathbb{R}^{k}$ \text{ and } \mathbf{l} \leq \mathbf{\tilde{x}}_2 \leq \mathbf{l}\}
+    \in \mathbb{R}^{k} \text{ and } \mathbf{l} \leq \mathbf{\tilde{x}}_2 \leq \mathbf{l}\}
 
     \mathbf{M} =
         \left[ {\begin{array}{cc}

docs/source/explanations/optimization/internal_optimizers.rst

@@ -43,6 +43,7 @@ this interface. The mandatory conditions for an internal optimizer function are:
      first derivative jointly
    - lower_bounds: for lower bounds in form of a 1d numpy array
    - upper_bounds: for upper bounds in form of a 1d numpy array
+   - nonlinear_constraints: for nonlinear constraints in form a list of dictionaries
 
    Of course, algorithms that do not need a certain argument (e.g. unbounded or
    derivative free ones) do not need those arguments at all.
@@ -98,6 +99,53 @@ different optimizers. However, complete transparency is possible and we try to d
 the exact meaning of all options for all optimizers.
 
 
+Nonlinear constraints
+---------------------
+
+Estimagic can pass nonlinear constraints to the internal optimizer. The internal
+interface for nonlinear constraints is as follows.
+
+A nonlinear constraint is a ``list`` of ``dict`` 's, where each ``dict`` represents a
+group of constraints. In each group the constraint function can potentially be
+multi-dimensional. We distinguish between equality and inequality constraints, which is
+signalled by a dict entry ``type`` that takes values ``"eq"`` and ``"ineq"``. The
+constraint function, which takes as input an internal parameter vector, is stored under
+the entry ``fun``, while the Jacobian of that function is stored at ``jac``. The
+tolerance for the constraints is stored under ``tol``. At last, the number of
+constraints in each group is specified under ``n_constr``. An example list with one
+constraint that would be passed to the internal optimizer is given by
+
+.. code-block::
+
+  constraints = [
+      {
+          "type": "ineq",
+          "n_constr": 1,
+          "tol": 1e-5,
+          "fun": lambda x: x**3,
+          "jac": lambda x: 3 * x**2,
+      }
+  ]
+
+
+.. note::
+
+  **Equality.** Internal equality constraints assume that the constraint is met when the function is
+  zero. That is
+
+  .. math::
+
+    0 = g(x) \in \mathbb{R}^m .
+
+  **Inequality.** Internal inequality constraints assume that the constraint is met when the function is
+  greater or equal to zero. That is
+
+  .. math::
+
+    0 \leq g(x) \in \mathbb{R}^m .
+
+
+
 Other conventions
 -----------------
 

docs/source/how_to_guides/optimization/how_to_specify_constraints.rst

@@ -10,25 +10,25 @@ Constraints vs bounds
 
 Estimagic distinguishes between bounds and constraints. Bounds are lower and upper
 bounds for parameters. In the literature they are sometimes called box constraints.
-Bounds are specified as "lower_bound" and "upper_bound" column of a params DataFrame
-or as pytrees via the ``lower_bounds`` and ``upper_bounds`` argument to ``maximize`` and
-``minimize``.
+Bounds are specified as ``lower_bounds`` and ``upper_bounds`` argument to ``maximize``
+and ``minimize``.
 
 Examples with bounds can be found in `this tutorial`_.
 
 .. _this tutorial: ../../getting_started/first_optimization_with_estimagic.ipynb
 
-Constraints are more general constraints on the parameters. This ranges from rather
-simple ones (e.g. Parameters are fixed to a value, a group of parameters is required
-to be equal) to more complex ones (like general linear constraints).
+To specify more general constraints on the parameters you use can use the argument
+``constraints``. This ranges from rather simple ones (e.g. parameters are fixed to a
+value, a group of parameters is required to be equal) to more complex ones (like general
+linear constraints, or even nonlinear constraints).
 
 Can you use constraints with all optimizers?
 ============================================
 
-We implement constraints via reparametrizations. Details are explained `here`_. This
-means that you can use all of the constraints with any optimizer that supports
-bounds. Some constraints (e.g. fixing parameters) can even be used with optimizers
-that do not support bounds.
+With the exception of general nonlinear constraints, we implement constraints via
+reparametrizations. Details are explained `here`_. This means that you can use all of
+the constraints with any optimizer that supports bounds. Some constraints (e.g. fixing
+parameters) can even be used with optimizers that do not support bounds.
 
 .. _here: ../../explanations/optimization/implementation_of_constraints.rst
 
@@ -288,8 +288,8 @@ flat numpy array are explained in the next section.
     typically it is more convenient to use the special cases instead of expressing
     them as a linear constraint. Internally, it will make no difference.
 
-    Let's impose the constraint that the sum of the average of the first four parameters
-    is at least 0.95.
+    Let's impose the constraint that the average of the first four parameters is at
+    least 0.95.
 
     .. code-block:: python
 
@@ -315,6 +315,44 @@ flat numpy array are explained in the next section.
     can also be arrays (or even pytrees) with bounds and weights for each selected
     parameter.
 
+.. tabbed:: nonlinear
+
+    .. warning::
+
+        General nonlinear constraints that are specified via a black-box constraint
+        function can only be used if you choose an optimizer that supports it.
+        The feature is currently supported by the algorithms:
+
+        * ``ipopt``
+        * ``nlopt``: ``cobyla``, ``slsqp``, ``isres``, ``mma``
+        * ``scipy``: ``cobyla``, ``slsqp``, ``trust_constr``
+
+    You can use nonlinear constraints to express restrictions of the form
+    ``lower_bound <= func(x) <= upper_bound`` or
+    ``func(x) = value`` where ``x`` are the selected parameters and ``func`` is the
+    constraint function.
+
+    Let's impose the constraint that the product of all but the last parameter is 1.
+
+    .. code-block:: python
+
+        res = minimize(
+            criterion=criterion,
+            params=np.ones(6),
+            algorithm="scipy_slsqp",
+            constraints={
+                "type": "nonlinear",
+                "selector": lambda x: x[:-1],
+                "func": lambda x: np.prod(x),
+                "value": 1.0,
+            },
+        )
+
+    This yields:
+
+    >>> array([ 1.31,  1.16,  1.01,  0.87,  0.75, -0.  ])
+
+    Where the product of the all but the last parameters is equal to 1.
 
 
 Imposing multiple constraints at once

src/estimagic/optimization/algo_options.py

@@ -179,6 +179,13 @@
 
 """
 
+CONSTRAINTS_ABSOLUTE_TOLERANCE = 1e-5
+"""float: Allowed tolerance of the equality and inequality constraints for values to be
+considered 'feasible'.
+
+"""
+
+
 """
 -------------------------
 Trust Region Parameters

src/estimagic/optimization/cyipopt_optimizers.py

@@ -24,6 +24,8 @@ def ipopt(
     lower_bounds,
     upper_bounds,
     *,
+    # nonlinear constraints
+    nonlinear_constraints=(),
     # convergence criteria
     convergence_relative_criterion_tolerance=CONVERGENCE_RELATIVE_CRITERION_TOLERANCE,
     dual_inf_tol=1.0,
@@ -494,7 +496,7 @@ def ipopt(
         x0=x,
         bounds=_get_scipy_bounds(lower_bounds, upper_bounds),
         jac=derivative,
-        constraints=(),
+        constraints=nonlinear_constraints,
         tol=convergence_relative_criterion_tolerance,
         options=options,
     )

src/estimagic/optimization/get_algorithm.py

@@ -49,6 +49,7 @@ def get_final_algorithm(
     valid_kwargs,
     lower_bounds,
     upper_bounds,
+    nonlinear_constraints,
     algo_options,
     logging,
     db_kwargs,
@@ -64,6 +65,8 @@ def get_final_algorithm(
             algorithm function.
         lower_bounds (np.ndarray): 1d numpy array with lower bounds.
         upper_bounds (np.ndarray): 1d numpy array with upper bounds.
+        nonlinear_constraints (list[dict]): List of dictionaries, each containing the
+            specification of a nonlinear constraint.
         algo_options (dict): Dictionary with additional keyword arguments for the
             algorithm. Entries that are not used by the algorithm are ignored with a
             warning.
@@ -80,6 +83,7 @@ def get_final_algorithm(
         algo_options=algo_options,
         lower_bounds=lower_bounds,
         upper_bounds=upper_bounds,
+        nonlinear_constraints=nonlinear_constraints,
         algo_name=algo_name,
         valid_kwargs=valid_kwargs,
     )
@@ -253,6 +257,7 @@ def _adjust_options_to_algorithm(
     algo_options,
     lower_bounds,
     upper_bounds,
+    nonlinear_constraints,
     algo_name,
     valid_kwargs,
 ):
@@ -291,4 +296,7 @@ def _adjust_options_to_algorithm(
     if "upper_bounds" in valid_kwargs:
         reduced["upper_bounds"] = upper_bounds
 
+    if "nonlinear_constraints" in valid_kwargs:
+        reduced["nonlinear_constraints"] = nonlinear_constraints
+
     return reduced