Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG: trust-constr evaluates function out of bounds #11712

Merged
merged 6 commits into from
Oct 27, 2022

Conversation

dpoerio
Copy link
Contributor

@dpoerio dpoerio commented Mar 22, 2020

Reference issue

Attempts to close #11649, different approach than #11650

What does this implement/fix?

Finite difference bounds are now passed to tr_interior_point, and responsibility for checking feasibility of the current solution is in the BarrierSubproblem class. If the current solution is out of bounds, the objective and constraint functions are not evaluated and are set to inf and 0, respectively. This causes the trust region to shrink until the current solution is feasible again.

@miladsade96 miladsade96 added defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.optimize labels Mar 23, 2020
@antonior92
Copy link
Member

antonior92 commented Apr 27, 2020

Ohh, Sorry for taking soo long to respond.

I not entirely sure this solution is better than #11650. At least not in principle.

#11650 -> Create a patch that in _differentiable_functions.py
this PR -> Creates a patch in minimize_trust_constr.py

I think we should try, as much as possible, to abstract away details about the differentiation from the implementation of the optimization algorithms. This make it easier to understand the code and make both the abstraction DifferentiableFunction and the optimization algorithm easier to reuse and combine in different ways.

Simply returning np.inf when outside the bounds would not fix the problem? I suggested that in #11650

@dpoerio
Copy link
Contributor Author

dpoerio commented May 18, 2020

@antonior92 Perhaps I'm overlooking something, it's been some time since I looked at this, but I don't see how to do that without modifying the _differentiable_functions, which @andyfaff said he preferred to avoid. This is basically implementing his suggestion of dealing with this in tr_interior_point.

I think, in principle, you could simply return inf when an out of bounds evaluation is attempted in the _differentiable_functions, but is that the preferred general behavior for every optimization routine that uses the _differentiable_functions? I'm not sure that would be preferred for all use cases, though maybe it is. Either way, I think you would still need to modify the _differentiable_functions, just with a slightly different approach that #11650 (returning inf instead of raising a value error).

@dpoerio
Copy link
Contributor Author

dpoerio commented Dec 13, 2020

@antonior92 I have some time to work on this bug and I'd like to get it resolved. I'm looking back over at what has been attempted thus far. As far as I can tell, what is implemented in this PR is the most straightforward way to ensure that a bound constraint violating function evaluation is not attempted. In this PR, we check for feasibility within the interior point algorithm and return inf if the bound constraints are violated. We could instead make this the responsibility of the _differentiable_functions but I think this would need to touch a lot more code. For instance, the variable in the ScalarFunction class that would allow us to check this is finite_diff_options. Consider the call stack again:

Traceback (most recent call last):
  File "/home/dvp/.conda/envs/scipydev/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3418, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-5-d298ce39ec79>", line 9, in <module>
    tst = opt.minimize(fun =  of, x0 = x0, method = 'trust-constr', bounds = bnds, constraints = nlcs)
  File "/home/dvp/Projects/scipy/scipy/optimize/_minimize.py", line 630, in minimize
    return _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py", line 509, in _minimize_trustregion_constr
    _, result = tr_interior_point(
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/tr_interior_point.py", line 322, in tr_interior_point
    z, state = equality_constrained_sqp(
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py", line 147, in equality_constrained_sqp
    f_next, b_next = fun_and_constr(x_next)
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/tr_interior_point.py", line 83, in function_and_constraints
    f = self.fun(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 242, in fun
    self._update_x_impl(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 215, in update_x
    self._update_hess()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 236, in _update_hess
    self._update_hess_impl()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 200, in update_hess
    self._update_grad()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 231, in _update_grad
    self._update_grad_impl()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 151, in update_grad
    self.g = approx_derivative(fun_wrapped, self.x, f0=self.f,
  File "/home/dvp/Projects/scipy/scipy/optimize/_numdiff.py", line 451, in approx_derivative
    raise ValueError("`x0` violates bound constraints.")
ValueError: `x0` violates bound constraints.

This PR implements the check in tr_interior_point.py, line 82. Once we go into _differentiable_functions.py, finite_diff_options is not in scope (as currently implemented) until approx_derivative is called. It doesn't seem to me there is a graceful way to exit at that point, so the class would have to be reconfigured a bit. Also, would we want the same behavior for all other _differentiable_functions? What do you think?

@andyfaff
Copy link
Contributor

One needs to be careful when approx_derivative or differentiable functions are modified. Modifications have to make sense from a context in which all of the minimizers in optimize are considered, because the aforementioned code underlie all of them.

@andyfaff
Copy link
Contributor

Approx_derivative will obey any bounds for numerical differentiation, but requesting gradients for an x outside the bounds is a definite no-no

@dpoerio
Copy link
Contributor Author

dpoerio commented Dec 13, 2020

One needs to be careful when approx_derivative or differentiable functions are modified. Modifications have to make sense from a context in which all of the minimizers in optimize are considered, because the aforementioned code underlie all of them.

Agreed, I think this is why it makes sense to keep this fix local to the interior point algorithm, at least for now. I've not observed the bug with other optimizers.

@antonior92
Copy link
Member

Hi all, it has been some time since I last checked this thread. So correct me if I miss something.

But I think I have never suggested changing approx_derivative. As I said, in this comment here I believe the problem is due to the quasi-newton approximation, which does ask for approx_derivative to compute gradients out of bounds...

In this sense, this is very much a problem of _differentiable_functions.py in general. I think the reason for not observing it in other optimizers is because different optimizers rely on _differentiable_functions.py to a different extent. Some of them have implementations of quasi-newton approximations embedded in them...

Anyway, I still think trying to fix it here in minimize_trust_constr.py doesn't attack the real problem (which is on the quasi-newton implementation) . It also breaks some useful abstractions. As I said before:

I think we should try, as much as possible, to abstract away details about the differentiation from the implementation of the optimization algorithms. This make it easier to understand the code and make both the abstraction DifferentiableFunction and the optimization algorithm easier to reuse and combine in different ways.

Have you tried, just changing:

        elif isinstance(hess, HessianUpdateStrategy):
            self.H = hess
            self.H.initialize(self.n, 'hess')
            self.H_updated = True
            self.x_prev = None
            self.g_prev = None

            def update_hess():
                self._update_grad()
                self.H.update(self.x - self.x_prev, self.g - self.g_prev)

        self._update_hess_impl = update_hess

to:

        elif isinstance(hess, HessianUpdateStrategy):
            self.H = hess
            self.H.initialize(self.n, 'hess')
            self.H_updated = True
            self.x_prev = None
            self.g_prev = None

            def update_hess():
                if ``x inside bounds``:
                    self._update_grad()
                    self.H.update(self.x - self.x_prev, self.g - self.g_prev)

        self._update_hess_impl = update_hess

I totally agree with @andyfaff that we need to be very careful when changing the ``_differentiable_functions.py`. But I think a punctual change is probably the best option here...

@dpoerio
Copy link
Contributor Author

dpoerio commented Dec 23, 2020

@antonior92 Thanks for your comments. I briefly tried your suggested change, and it is not working as-is. First, as noted above, the variable that would allow us to check the bounds is not an instance variable of the ScalarFunction class and is not in scope through all the calls in differentiable_functions.py. So first I assigned it as an instance variable:

        finite_diff_options = {}
        self.finite_diff_options = finite_diff_options

Then I implement your suggested change:

        elif isinstance(hess, HessianUpdateStrategy):
            self.H = hess
            self.H.initialize(self.n, 'hess')
            self.H_updated = True
            self.x_prev = None
            self.g_prev = None

            def update_hess():
                if not np.any((self.x < self.finite_diff_options['bounds'][0])
                           | (self.x > self.finite_diff_options['bounds'][1])):
                    self._update_grad()
                    self.H.update(self.x - self.x_prev, self.g - self.g_prev)

We still error out in the next line in tr_interior_point.py, the call to self.constr(). Here's the traceback:

Traceback (most recent call last):
  File "/home/dvp/.conda/envs/scipydev/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3418, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-d298ce39ec79>", line 9, in <module>
    tst = opt.minimize(fun =  of, x0 = x0, method = 'trust-constr', bounds = bnds, constraints = nlcs)
  File "/home/dvp/Projects/scipy/scipy/optimize/_minimize.py", line 630, in minimize
    return _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py", line 509, in _minimize_trustregion_constr
    _, result = tr_interior_point(
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/tr_interior_point.py", line 321, in tr_interior_point
    z, state = equality_constrained_sqp(
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py", line 147, in equality_constrained_sqp
    f_next, b_next = fun_and_constr(x_next)
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/tr_interior_point.py", line 83, in function_and_constraints
    c_eq, c_ineq = self.constr(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/canonical_constraint.py", line 104, in fun
    *[c.fun(x) for c in canonical_constraints])
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/canonical_constraint.py", line 104, in <listcomp>
    *[c.fun(x) for c in canonical_constraints])
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/canonical_constraint.py", line 238, in fun
    return empty_fun, lb - cfun.fun(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 512, in fun
    self._update_x(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 494, in _update_x
    self._update_x_impl(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 477, in update_x
    self._update_hess()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 508, in _update_hess
    self._update_hess_impl()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 458, in update_hess
    self._update_jac()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 503, in _update_jac
    self._update_jac_impl()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 406, in update_jac
    approx_derivative(fun_wrapped, self.x, f0=self.f,
  File "/home/dvp/Projects/scipy/scipy/optimize/_numdiff.py", line 451, in approx_derivative
    raise ValueError("`x0` violates bound constraints.")
ValueError: `x0` violates bound constraints.

In the call to the constraint evaluation, it looks like we are using the VectorFunction class. I can implement the same change as above for this class:

        finite_diff_options = {}
        self.finite_diff_options = finite_diff_options

and

        elif isinstance(hess, HessianUpdateStrategy):
            self.H = hess
            self.H.initialize(self.n, 'hess')
            self.H_updated = True
            self.x_prev = None
            self.J_prev = None

            def update_hess():
                if not np.any((self.x < self.finite_diff_options['bounds'][0])
                           | (self.x > self.finite_diff_options['bounds'][1])):
                    self._update_jac()
                    # When v is updated before x was updated, then x_prev and
                    # J_prev are None and we need this check.
                    if self.x_prev is not None and self.J_prev is not None:
                        delta_x = self.x - self.x_prev
                        delta_g = self.J.T.dot(self.v) - self.J_prev.T.dot(self.v)
                        self.H.update(delta_x, delta_g)

        self._update_hess_impl = update_hess

This allows us to get past both initial failing calls to self.fun() and self.constr(). However, we error out at a different point in the next sequence of calls:

Traceback (most recent call last):
  File "/home/dvp/.conda/envs/scipydev/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3418, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-d298ce39ec79>", line 9, in <module>
    tst = opt.minimize(fun =  of, x0 = x0, method = 'trust-constr', bounds = bnds, constraints = nlcs)
  File "/home/dvp/Projects/scipy/scipy/optimize/_minimize.py", line 630, in minimize
    return _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py", line 509, in _minimize_trustregion_constr
    _, result = tr_interior_point(
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/tr_interior_point.py", line 322, in tr_interior_point
    z, state = equality_constrained_sqp(
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py", line 147, in equality_constrained_sqp
    f_next, b_next = fun_and_constr(x_next)
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/tr_interior_point.py", line 83, in function_and_constraints
    f = self.fun(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/canonical_constraint.py", line 104, in fun
    *[c.fun(x) for c in canonical_constraints])
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/canonical_constraint.py", line 104, in <listcomp>
    *[c.fun(x) for c in canonical_constraints])
  File "/home/dvp/Projects/scipy/scipy/optimize/_trustregion_constr/canonical_constraint.py", line 238, in fun
    return empty_fun, lb - cfun.fun(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 519, in fun
    self._update_x(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 501, in _update_x
    self._update_x_impl(x)
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 477, in update_x
    self._update_jac()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 510, in _update_jac
    self._update_jac_impl()
  File "/home/dvp/Projects/scipy/scipy/optimize/_differentiable_functions.py", line 411, in update_jac
    approx_derivative(fun_wrapped, self.x, f0=self.f,
  File "/home/dvp/Projects/scipy/scipy/optimize/_numdiff.py", line 451, in approx_derivative
    raise ValueError("`x0` violates bound constraints.")
ValueError: `x0` violates bound constraints.

Do we want to just add checks for constraint violations at every possible point in _differentiable_functions.py? Adding this check in many places does not strike me as the most satisfying solution.

@antonior92
Copy link
Member

antonior92 commented Dec 23, 2020

Thank you for trying it out so quickly @dpoerio.

Do we want to just add checks for constraint violations at every possible point in _differentiable_functions.py? Adding this check in many places does not strike me as the most satisfying solution.

I agree that adding checks everywhere is definitely not a good solution. I still think the problem is general from _differentiable_functions.py, but I can't think of some easy fix here.

I agree that your solution is a clean and easy one. Since it solves the problem where it is most apparent I think we can merge it. @andyfaff what do you think?

@mdhaber
Copy link
Contributor

mdhaber commented Oct 8, 2022

I went ahead and resolved the merge conflicts, but I noticed that PEP8 issues will need to be resolved. I can do that.

Otherwise, @andyfaff did this look good to you? It seemed to have Antonio's approval. If you're happy with the approach - at least if it's better than not having it at all - it might be worth merging.

@mdhaber mdhaber requested a review from andyfaff October 15, 2022 20:10
@antonior92
Copy link
Member

Hi, I agree with merging it

@mdhaber
Copy link
Contributor

mdhaber commented Oct 27, 2022

Thanks @antonior92. Normally I'd hesitate to make the objective/constraints non-smooth like this, but it sounds like it is compatible with the algorithm and fixes the bug - plus it has the support of the trust-constr author! If there is a better way of handling this somewhere out there, we'd welcome the additional enhancement. In the meantime, it is important to respect keep_feasible=True if that's something we offer, so I'll merge when tests pass.

The CI failure looked unrelated, so I went ahead and merged.

@mdhaber mdhaber merged commit 48c3dc8 into scipy:main Oct 27, 2022
@dpoerio
Copy link
Contributor Author

dpoerio commented Oct 28, 2022

@mdhaber It looks like this is causing a failure at the CI for ARM Macs. I guess this is not tested for PRs? I don't have the hardware handy to debug this (without provisioning a cloud machine at least). On first pass, it looks like just a warning, maybe caused by numerous function evaluations returning inf at multiple points beyond the bounds, making the function appear constant/linear.

FAILED scipy/optimize/tests/test_minimize_constrained.py::test_gh11649 - UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.

Perhaps we should revert until this can be debugged/investigated?

mdhaber added a commit to mdhaber/scipy that referenced this pull request Oct 30, 2022
WarrenWeckesser added a commit that referenced this pull request Nov 6, 2022
MAINT: optimize.minimize: revert gh-11712 to fix macos_arm64_test
@mdhaber mdhaber added this to the 1.10.0 milestone Nov 29, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.optimize
Projects
None yet
Development

Successfully merging this pull request may close these issues.

trust-constr error when attempting to keep bound constrained solutions feasible
5 participants