Repository: scipy/scipy
Pull request: BUG: trust-constr Bounds exclusive
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Depending on your changes, you can skip CI operations and save time and energy: 
http://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes #15093
#### What does this implement/fix?
<!--Please explain your changes.-->
The scipy docs define Bounds as inclusive but trust-constr treats bounds as exclusive. According to Nocedal, Jorge, and Stephen J. Wright. "Numerical optimization" Second Edition (2006), p. 393, interior point methods necessarily treat bounds as exclusive. To work around this limitation of interior point methods and create consistent behavior for Bounds, this PR modifies trust-constr so that Bounds passed to trust-constr are expanded by epsilon.
#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
BUG: scipy.optimize's trust-constr algorithm hangs when keep-feasible constraint is active at x0
### Describe your issue.

Scipy.optimize's trust-region constrained algorithm fails to find an optimum when the initial solution x0 lies exactly on the constraint boundary, but only when that constraint is set to "keep-feasible". Instead of converging, the solution is kept at x0.

The reason lies in the barrier subproblem (tr_interior_point.py, lines 92-95):
s[self.enforce_feasibility] = -c_ineq[self.enforce_feasibility]
log_s = [np.log(s_i) if s_i > 0 else -np.inf for s_i in s]
return f - self.barrier_parameter*np.sum(log_s)

The initial slack variable s will be zero EXACTLY, so the log_s will become infinite. The reduction ratio will then compute inf-inf and become NaN (hence the RuntimeWarning in the code example), after which all output becomes garbage and maxiter is reached.

A simple solution I see is to raise a warning if the initial solution x0 lies exactly on the constraint boundary, but I can't shake the feeling that there exists a more elegant solution.

### Reproducing Code Example

```python
from scipy.optimize import minimize, Bounds

import numpy as np

x0 = np.array([0., 0.5])


def obj(x):
    x1 = x[0]
    x2 = x[1]
    return x1**2 + x2**2

bounds = Bounds(np.array([0., 0.]), np.array([1., 1.]), keep_feasible=False)
res = minimize(
            method='trust-constr',
            fun=obj,
            x0=x0,
            options={'disp': True,
                     'verbose': 3,
                    },
            bounds=bounds)
```


### Error message

```shell
When keep_feasible=False: 
| niter |f evals|CG iter|  obj func   |tr radius |   opt    |  c viol  | penalty  |barrier param|CG stop|
|-------|-------|-------|-------------|----------|----------|----------|----------|-------------|-------|
|   1   |   3   |   0   | +2.5000e-01 | 1.00e+00 | 3.33e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   0   |
|   2   |   6   |   1   | +4.4923e-02 | 7.00e+00 | 2.22e-02 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   3   |   9   |   2   | +6.6160e-02 | 7.00e+00 | 5.33e-03 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   4   |
|   4   |   9   |   2   | +6.6160e-02 | 3.50e+01 | 1.36e-02 | 0.00e+00 | 1.00e+00 |  2.00e-02   |   0   |
|   5   |  12   |   3   | +2.8248e-02 | 3.50e+01 | 1.49e-03 | 0.00e+00 | 1.00e+00 |  2.00e-02   |   4   |
|   6   |  12   |   3   | +2.8248e-02 | 1.75e+02 | 3.12e-03 | 0.00e+00 | 1.00e+00 |  4.00e-03   |   0   |
|   7   |  15   |   4   | +8.9077e-03 | 1.75e+02 | 3.75e-04 | 0.00e+00 | 1.00e+00 |  4.00e-03   |   4   |
|   8   |  18   |   5   | +4.5372e-03 | 1.75e+02 | 3.88e-05 | 0.00e+00 | 1.00e+00 |  4.00e-03   |   4   |
|   9   |  18   |   5   | +4.5372e-03 | 8.75e+02 | 1.84e-04 | 0.00e+00 | 1.00e+00 |  8.00e-04   |   0   |
|  10   |  21   |   6   | +1.5528e-03 | 8.75e+02 | 2.22e-05 | 0.00e+00 | 1.00e+00 |  8.00e-04   |   4   |
|  11   |  21   |   6   | +1.5528e-03 | 4.37e+03 | 3.96e-05 | 0.00e+00 | 1.00e+00 |  1.60e-04   |   0   |
|  12   |  24   |   7   | +4.7056e-04 | 4.37e+03 | 4.89e-06 | 0.00e+00 | 1.00e+00 |  1.60e-04   |   4   |
|  13   |  27   |   8   | +2.1014e-04 | 4.37e+03 | 5.42e-07 | 0.00e+00 | 1.00e+00 |  1.60e-04   |   4   |
|  14   |  27   |   8   | +2.1014e-04 | 2.19e+04 | 1.84e-06 | 0.00e+00 | 1.00e+00 |  3.20e-05   |   0   |
|  15   |  30   |   9   | +6.9638e-05 | 2.19e+04 | 2.25e-07 | 0.00e+00 | 1.00e+00 |  3.20e-05   |   4   |
|  16   |  33   |  10   | +3.6985e-05 | 2.19e+04 | 2.22e-08 | 0.00e+00 | 1.00e+00 |  3.20e-05   |   4   |
|  17   |  33   |  10   | +3.6985e-05 | 1.09e+05 | 1.32e-07 | 0.00e+00 | 1.00e+00 |  6.40e-06   |   0   |
|  18   |  36   |  11   | +1.2713e-05 | 1.09e+05 | 1.60e-08 | 0.00e+00 | 1.00e+00 |  6.40e-06   |   4   |
|  19   |  36   |  11   | +1.2713e-05 | 5.47e+05 | 2.89e-08 | 0.00e+00 | 1.00e+00 |  1.28e-06   |   0   |
|  20   |  39   |  12   | +3.8496e-06 | 5.47e+05 | 3.57e-09 | 0.00e+00 | 1.00e+00 |  1.28e-06   |   4   |

`gtol` termination condition is satisfied.
Number of iterations: 20, function evaluations: 39, CG iterations: 12, optimality: 3.57e-09, constraint violation: 0.00e+00, execution time: 0.055 s.
When keep_feasible=True:
| niter |f evals|CG iter|  obj func   |tr radius |   opt    |  c viol  | penalty  |barrier param|CG stop|
|-------|-------|-------|-------------|----------|----------|----------|----------|-------------|-------|
|   1   |   3   |   0   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   0   |
|   2   |   6   |   1   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   3   |   9   |   3   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   4   |  12   |   5   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   5   |  15   |   7   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   6   |  18   |   9   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   7   |  21   |  11   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   8   |  24   |  12   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|   9   |  27   |  14   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  10   |  30   |  16   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  11   |  33   |  18   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
.conda\envs\python\lib\site-packages\scipy\optimize\_trustregion_constr\equality_constrained_sqp.py:152: RuntimeWarning: invalid value encountered in double_scalars
  actual_reduction = merit_function - merit_function_next
|  12   |  36   |  20   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  13   |  39   |  22   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  14   |  42   |  24   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  15   |  45   |  26   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  16   |  48   |  28   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  17   |  51   |  30   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  18   |  54   |  32   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  19   |  57   |  34   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  20   |  60   |  36   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  21   |  63   |  38   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  22   |  66   |  40   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  23   |  69   |  42   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  24   |  72   |  44   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |
|  25   |  75   |  46   | +2.5000e-01 | 1.00e+00 | 1.11e-01 | 0.00e+00 | 1.00e+00 |  1.00e-01   |   2   |

The maximum number of function evaluations is exceeded.
Number of iterations: 25, function evaluations: 75, CG iterations: 46, optimality: 1.11e-01, constraint violation: 0.00e+00, execution time: 0.03 s.
```


### SciPy/NumPy/Python version information

1.7.1 1.21.2 sys.version_info(major=3, minor=8, micro=12, releaselevel='final', serial=0)

Hints from issue discussion:

Commit message: BUG: trust-constr Bounds exclusive
Closes #15093
The scipy docs define Bounds as inclusive but trust-constr treats
bounds as exclusive. According to Nocedal, Jorge, and Stephen J.
Wright. "Numerical optimization" Second Edition (2006), p. 393,
interior point methods necessarily treat bounds as exclusive. To work
around this limitation of interior point methods and create consistent
behavior for Bounds, this PR modifies trust-constr so that Bounds
passed to trust-constr are expanded by epsilon.

Patch to review:
diff --git a/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py b/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py
index 1c86dd869b78..2835ea5445c0 100644
--- a/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py
+++ b/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py
@@ -3,7 +3,7 @@
 from scipy.sparse.linalg import LinearOperator
 from .._differentiable_functions import VectorFunction
 from .._constraints import (
-    NonlinearConstraint, LinearConstraint, PreparedConstraint, strict_bounds)
+    NonlinearConstraint, LinearConstraint, PreparedConstraint, Bounds, strict_bounds)
 from .._hessian_update_strategy import BFGS
 from .._optimize import OptimizeResult
 from .._differentiable_functions import ScalarFunction
@@ -323,6 +323,11 @@ def _minimize_trustregion_constr(fun, x0, args, grad,
         verbose = 1
 
     if bounds is not None:
+        modified_lb = np.nextafter(bounds.lb, -np.inf, where=bounds.lb > -np.inf)
+        modified_ub = np.nextafter(bounds.ub, np.inf, where=bounds.ub < np.inf)
+        modified_lb = np.where(np.isfinite(bounds.lb), modified_lb, bounds.lb)
+        modified_ub = np.where(np.isfinite(bounds.ub), modified_ub, bounds.ub)
+        bounds = Bounds(modified_lb, modified_ub, keep_feasible=bounds.keep_feasible)
         finite_diff_bounds = strict_bounds(bounds.lb, bounds.ub,
                                            bounds.keep_feasible, n_vars)
     else:
diff --git a/scipy/optimize/tests/test_minimize_constrained.py b/scipy/optimize/tests/test_minimize_constrained.py
index 0c4a8d4676ad..587d7256c2ba 100644
--- a/scipy/optimize/tests/test_minimize_constrained.py
+++ b/scipy/optimize/tests/test_minimize_constrained.py
@@ -615,6 +615,32 @@ def callback(x, info):
         # compatibility
         assert_(result.get('niter', -1) == 1)
 
+    def test_issue_15093(self):
+        # scipy docs define bounds as inclusive, so it shouldn't be
+        # an issue to set x0 on the bounds even if keep_feasible is
+        # True. Previously, trust-constr would treat bounds as
+        # exclusive.
+
+        x0 = np.array([0., 0.5])
+
+        def obj(x):
+            x1 = x[0]
+            x2 = x[1]
+            return x1 ** 2 + x2 ** 2
+
+        bounds = Bounds(np.array([0., 0.]), np.array([1., 1.]),
+                        keep_feasible=True)
+
+        with suppress_warnings() as sup:
+            sup.filter(UserWarning, "delta_grad == 0.0")
+            result = minimize(
+                method='trust-constr',
+                fun=obj,
+                x0=x0,
+                bounds=bounds)
+
+        assert_(result.get('success'))
+
 class TestEmptyConstraint(TestCase):
     """
     Here we minimize x^2+y^2 subject to x^2-y^2>1.


Code review comments:
diff_hunk: @@ -615,6 +615,32 @@ def callback(x, info):
         # compatibility
         assert_(result.get('niter', -1) == 1)
 
+    def test_issue_15093(self):
+        # scipy docs define bounds as inclusive, so it shouldn't be
+        # an issue to set x0 on the bounds even if keep_feasible is
+        # True. Previously, trust-constr would treat bounds as
+        # exclusive.
+
+        x0 = np.array([0., 0.5])
+
+        def obj(x):
+            x1 = x[0]
+            x2 = x[1]
+            return x1 ** 2 + x2 ** 2
+
+        bounds = Bounds(np.array([0., 0.]), np.array([1., 1.]),
+                        keep_feasible=True)
+
+        with suppress_warnings() as sup:
+            sup.filter(UserWarning, "delta_grad == 0.0")
+            result = minimize(
+                method='trust-constr',
+                fun=obj,
+                x0=x0,
+                bounds=bounds)
+
+        assert_(result.get('success'))
path: scipy/optimize/tests/test_minimize_constrained.py
text: @user1:
Minor: we prefer plain `assert` in *new* test cases.

More detailed:

I checked that this fails on `main` and passes here with `main` failure summarized below:

`scipy/optimize/tests/test_minimize_constrained.py::TestTrustRegionConstr::test_issue_15093 - RuntimeWarning: invalid value encountered in scalar subtract`

The original issue reported a hang with a different reproducer, but I can no longer reproduce that hang on `main`--can you clarify what is going on there?

If I take the original reproducer, and set `keep_feasible=True`, it fails on this branch as well--is that expected? The failure is different than on `main`, but still hard for me to judge.

@author:
Hi Tyler - thanks for the review!

I'll update the `_assert` to `assert`. As for the other issues, I'm not getting the same results as you.

> The original issue reported a hang with a different reproducer, but I can no longer reproduce that hang on main--can you clarify what is going on there?

I get the hang reported in the original issue when I run the original reproducer main. If you're seeing my test fail on main because of that RuntimeWarning, it's just that the test fails on the RuntimeWarning before it gets to the `assert_(result.get('success'))`.

> If I take the original reproducer, and set keep_feasible=True, it fails on this branch as well--is that expected? The failure is different than on main, but still hard for me to judge.

The original reproducer with keep_feasible=True works on this branch for me - I get an x0 of [0, 0.0436886] with keep_feasible=True which seems reasonable to me (compared to an x0 of [0.00138686, 0.00138786] with keep_feasible=False). It does produce a warning which says `delta_grad == 0.0`, which is expected.

@user2:
Having just examined the test problem I'm surprised that `result.x` doesn't arrive very close to the true solution (it's something like 2e-3 off zero). I wonder if the default gtol could be lowered.