scipy.optimize.minimize SLSQP leads to out of bounds solution

[ghost]

SLSQP algorithm goes to infinity without counting for bounds specified if local gradient in one of the directions is close to zero. This issue is found at 2D and 7D bounded constrained problems I'm running now.

objective function as:
def objective(x):
lambda _obj _x: (5._x-2)__2_sin(12.*x-4)
return _obj(numpy.linalg.norm(x)

bounds = ((-1,1),(-1,1))

[pv]

Could you give an actually runnable test script that shows this behavior?

[pv]

This works correctly:

import numpy
from numpy import *
from scipy.optimize import minimize

def objective(x):
    _obj = lambda x: (5.*x-2)*2*sin(12.*x-4)
    return _obj(numpy.linalg.norm(x))

bounds = ((-1.,1), (-1.,1))

sol = minimize(objective, [0.5, 0.75], bounds=bounds, method="SLSQP")
print sol

What sort of initial value etc are you using?

[ghost]

Good day!

Thank you for quick reply. Unfortunately I cannot recover now the state
of the code when that error happened. Now I'm running high-fidelity 7D
optimization including Ansys Fluent and I've already got this case
twice. When function is evaluated at values that are very far from
design space. I'm trying to reproduce the code when I've met this error
first time - analytic function involving Rbf interpolation. And will
email you when I'll get that issue again.

Thank you!

[pv]

Thanks! SLSQP is from ACM TOMS 733 (used in Scipy with permission from ACM), and issues like this are probably issues in the algorithm itself, and difficult to debug without a test case.

[yikelu]

Hey PV, I ran into this as well; however it is mostly due to numerical roundoff errors (for example, one of my parameters has a lower bound of 0, but SLSQP returns -1e120 or something silly). There are a few cases where it goes wildly off track, but that typically happens when the initial condition fails to respect the bounds -- I've had this happen using SLSQP as the local optimizer driving basinhopping.

[Marigold]

Happened to me too, here's minimal working example:

import scipy.optimize as optimize
import numpy as np

X = np.array([[1.020626, 1.013055], [0.989094, 1.059343]])
freq = 13.574380165289256
x_0 = [1., 1.]

def objective(b):
    def foo(r_log, freq):
        mu, sd = r_log.mean(), r_log.std()
        sd += 0.5 / freq
        return mu / sd * np.sqrt(freq)

    print(b)
    return -foo(np.log(np.maximum(np.dot(X - 1, b) + 1, 0.2)), freq=freq)

cons = ({'type': 'ineq', 'fun': lambda b: 2. - sum(b)},)
res = optimize.minimize(objective, x_0, bounds=[(0., 2.)]*len(x_0), constraints=cons, method='slsqp')
print(res)

Algorithm always gets to the correct result ([ 1.18818046 0.81181953]) and then suddenly jumps out of bounds.

[cfcohen]

It appears that I've encountered a similar issue. The following code did not produce the results that I expected. I don't know whether it's a bug or not, but I suspect that it is. Problematically (for me at least), the behavior is very different when undefined is set to None, 0.0, 1.0, etc. I expected that the bounds parameter would ensure that f(x) was not invoked with an invalid value of x. SLSQP is not the only minimizer that exhibits this problem. Sample code:

import scipy.optimize
undefined = None
def f(x):
    if x > 1.0:
        print("x > 1.0, x=%12.10f" % x)
        return undefined
    print("f(%f)=1.0" % x)
    return 1.0
# The bounds are not properly respected!
print scipy.optimize.minimize(f, 1.0, method='SLSQP', bounds=[(0.0, 1.0)])

[pv]

@cfcohen: that problem is due to the gradient numerical differentiation not respecting bounds, and the starting value lying at the edge.
@Marigold: that code works OK for me, stays within bounds.

[Marigold]

@pv yep, working like a charm now. Thanks.

[davidpasquale]

I actually still find the same issue (also with L-BFGS-B).

If the solution is the upper bound, the grandient is computed exceeding the upper bound instead of being computed inside the bound.
The workaround is to reduce the bound of a certain eps and use the same eps to compute the descrete grandient.
Pleas fix this BUG!

[mdhaber]

@davidpasquale Can you submit an example?
@Marigold's code is working without issue on Google Collab.
Tempted to close this unless we have an example.

[davidpasquale]

@mdhaber Here a basic example (the solution is also wrong!):

import numpy as np

def f(x):
  for i in range(x.shape[0]):
    if x[i] > 1.0:
        print("x", i, " > 1.0, x=%12.10f" % x[i])
        x[i] = 1.0
  print('x',  x)
  return np.sum(x**2)

lb = np.array([-1.0,-1.0])
ub = np.array([1.0, 1.0])
bnds = optimize.Bounds(lb, ub)
x0 = np.array([0.5, 1.0])
# The bounds are not properly respected!
res = optimize.minimize(f, x0, method='SLSQP', bounds=bnds, options={'eps': 1e-4})
print(res)

[davidpasquale]

@mdhaber here another example where also the result if wrong!:

from scipy import optimize
import numpy as np

def f(x):
  for i in range(x.shape[0]):
    if x[i] > 1.0:
        print("x", i, " > 1.0, x=%12.10f" % x[i])
        x[i] = 1.0
  print('x',  x)
  return x[0]*np.exp(-x[0]**2-x[1]**2)

lb = np.array([-2.0,-1.0])
ub = np.array([1.0, 1.0])
bnds = optimize.Bounds(lb, ub)
x0 = np.array([0.8, 0.0])
res = optimize.minimize(f, x0, method='SLSQP', bounds=bnds, options={'eps': 1e-6})
print(res)
print('Result:', res.x)

[mdhaber]

@davidpasquale I think that is the same sort of issue reported by @cfcohen (numerical differentiation not respecting bounds)

It can be argued that this is not a defect, as we typically consider constraints to be satisfied even when they are violated as long as the violation is within a certain tolerance. This is obvious for equality constraints, but we even allow a tolerance on inequality constraints: inequality constraints are often active at an optimal solution, "exact" equality is generally impossible to achieve, and it would require special logic (often not part of the original algorithm) to ensure that the inequality constraint is strictly respected. By the same logic, a tolerance on bounds might be allowed.

Then again, it's unlikely that this behavior is desirable, and it seems possible to change, so we might want to change it anyway.

@pv @antonior92 Is this something you think should be fixed?

[mdhaber]

@davidpasquale In the meantime you might try method='trust-constr' and as a third argument of your scipy.optimize.Bounds object specify keep_feasible=True. trust-constr does not violate the bounds while solving the second problem, at least.

But @antonior92 trust-constr fails to solve @davidpasquale's first problem with keep_feasible set to True. It hits the iteration limit with lots of warnings:

C:\Users\mhaberla\AppData\Local\Continuum\anaconda3\lib\site-packages\scipy\optimize\_trustregion_constr\equality_constrained_sqp.py:153: RuntimeWarning: invalid value encountered in double_scalars
  actual_reduction = merit_function - merit_function_next

Should I open a separate issue?

[mdhaber]

BTW @davidpasquale If you can, give analytical derivatives; all these examples solve without violating bounds when analytical derivatives are provided.

[davidpasquale]

Dear @mdhaber , thanks for your reply.
My example is just to show a simple case where the algorithm fails.
My problem is that I cannot violate the bounds ( i am dealing with parametric curves) and the optimization I am performing cannot be differenziated.
I think that the problem is just related on how you perform the gradiant approximation.
What about computing the gradient based only on left (right) side when you are at the upper (lower) bound ?

[mdhaber]

I understand, and it seems you're not alone. Yes, one-sided derivative approximation is one thing we'd need to do. There has been talk recently of overhauling the way the optimizers handle derivative approximation, so I'll add this to the conversation.

[andyfaff]

approx_derivative does take into account box bounds, and won't exceed a lower or upper limit to calculate a gradient. However, if the location at which it's asked to evaluate the gradient is outside the bounds, then it'll raise an exception.
For example, I naively used the bounds in LBFGSB in a gradient calculation, but it turns out that the current implementation asks for func and grad evaluation in locations outside the bounds, and there is therefore no point in implementing this functionality.
Edit: I was incorrect.

[mdhaber]

@andyfaff In your example, how far outside the bounds is 'L-BFGS-B' asking for the gradient to be evaluated? Considerably so, or is it within some sort of reasonable tolerance? Can you post that example?

Do you know of any such examples for TNC or SLSQP?

I've been looking for an example like that for #4916. I wrote a test to try to observe this sort of behavior, but even with an initial guess outside the bounds, every solver (except trust-constr with keep_feasible=False) brought the guess within bounds before evaluating the function. So perhaps the authors intended to stay strictly within bounds?

You wrote in #10112:

This is a simple lower/upper bound check. There should be zero tolerance for a solution not obeying box bounds.

So if possible, would you be in favor of trying to add a keep_feasible option to L-BFGS-B and others? While it might not work well for every algorithm (e.g. might affect convergence), we could try the simple thing of rounding candidate points to lie within bounds (and if necessary doing one-sided derivative approximations).

[andyfaff]

@andyfaff In your example, how far outside the bounds is 'L-BFGS-B' asking for the gradient to be evaluated? Considerably so, or is it within some sort of reasonable tolerance? Can you post that example?

In #10673 I thought I had come across such an example. It turns out that I'm wrong. The code I've introduced involves creation of a ScalarFunction object, and as part of the construction it evaluates the function and gradient with the initial guesses. If the initial guess is outside the bounds then there is an error raised. This can be avoided by clipping the initial guess.

I thought that this was the fault of the fortan optimizer loop, but it wasn't. Thankfully this issue made me look more closely and realise the mistake I was making.

The changes I'm making as part of 10673 (I've done LBFGSB, TNC, CG, BFGS) will prevent finite difference calculation going outside the bounds, so long as the guess obeys the bounds. It probably wouldn't work if the lower and upper bound are equal.
I've not looked into SLSQP yet.

[mdhaber]

so long as the guess obeys the bounds

Fortunately, it seems that LBFGSB, TNC, and SLSQP already clip the guess to stay within bounds, and trust-constr raises an error if the initial guess is out of bounds and keep_feasible is True, so that shouldn't be an issue for them.

You are adding bounds support for CG and BFGS? If so, maybe their guess should be clipped to match the rest.

[andyfaff]

You are adding bounds support for CG and BFGS? If so, maybe their guess should be clipped to match the rest.

I'm not adding bounds support. I'm just converting those to use approx_derivative.

[andyfaff]

25d7d56 converts SLSQP to use approx_derivative, which should respect parameter bounds. However, approx_derivative is only used for the gradient of the parameters, it's not used for calculating the gradient of the constraints.

[ghost]

I've had a lot of troubles because of this bug. I typically use SLSQP linked to external old Fortran programs like XFOIL and AVL. These codes provide output as a text file with 6-8 digits precision. My solution was to update analysis code to return 10-12 digits in exponential form (%.12e instead of %.6f). Also I always normalize bounds to [-1; 1] and set epsilon value to relatively large number from 5e-4 to 1e-3. It works in most cases, but sometimes I have doubts that solution is not a true optimum and some improvement is possible.

[andyfaff]

In #10673 the numerical differentiation code used in SLSQP, L-BFGS-B (and a few others) respects bounds. e.g. if a value is at an upper limit instead of a forward difference a reverse difference is used instead. This is implemented in optimize._numdiff.approx_derivative.
As of f865e9e all MWE presented in this issue work without problems.

[fuglede]

Just ran into this one for L-BFGS-B (SciPy 1.6.2):

  File "lib\site-packages\scipy\optimize\_minimize.py", line 619, in minimize
    return _minimize_lbfgsb(fun, x0, args, jac, bounds,
  File "lib\site-packages\scipy\optimize\lbfgsb.py", line 360, in _minimize_lbfgsb
    f, g = func_and_grad(x)
  File "lib\site-packages\scipy\optimize\_differentiable_functions.py", line 261, in fun_and_grad
    self._update_grad()
  File "lib\site-packages\scipy\optimize\_differentiable_functions.py", line 231, in _update_grad
    self._update_grad_impl()
  File "lib\site-packages\scipy\optimize\_differentiable_functions.py", line 151, in update_grad
    self.g = approx_derivative(fun_wrapped, self.x, f0=self.f,
  File "lib\site-packages\scipy\optimize\_numdiff.py", line 451, in approx_derivative
    raise ValueError("`x0` violates bound constraints.")

In my case, what I'm doing is using res = minimize(f1, x0, bounds=[(-20, 20)]*n), followed by minimize(f2, res.x, bounds=[(-20, 20)]*n); that is, using the result of one optimization as the initial guess for another one with the same bounds. If res.x is within bounds (as it should be), this should never fail.

I don't have a short repro, and even if I did, it wouldn't be helpful, since I have been unable to reproduce this even after quite a few attempts (RNG seems friendly to me).

[andyfaff]

@fuglede all the minimize methods are deterministic, so if you do find a MWE that demonstrates the issue please post it in a new issue (closed issues are not likely to be curated).

When an LBFGSB minimization is started the initial vector is clipped to bounds, so the initial guess should be within the bounds.

The numerical differentiation strictly obeys bounds now. This causes issues where an upper bound is equal to a lower bound, but that's not the case for the [(-20, 20) * n]. My guess is that the Fortran code for LBFGSB sometimes (albeit very rarely) proposes a new x where one of the values is ever so slightly outside one of the bounds, I'm talking a ULP or something along those lines. We've seen this happen in SLSQP (#11403), but not in LBFGSB. The fix is quite simple (https://github.com/scipy/scipy/blob/master/scipy/optimize/slsqp.py#L380), but it might not be worth doing until we can reproducibly create the problem.

[fuglede]

Ah, both of my f1 and f2 were non-deterministic. Perhaps then what's happening is simply that for such a function minimize might produce results outside of the bounds (and maybe I'm being overly optimistic in using such a non-deterministic thing with minimize in the first place (although it generally handles that situation really well!)? But still, you're saying that that result would be clipped, so even then it shouldn't really be a problem. I'll play around a bit more but doubt I'll find a useful reproducing example.

[andyfaff]

minimize shouldn't produce a solution that is out of bounds. Can you check that you're not altering x in your user function?

[fuglede]

Can you check that you're not altering x in your user function?

I'm not. (But would I expect that to make a difference here -- are there function calls taking place between the result being returned from one minimize, and the input validated in the next one?)

But again, I'm still unable to reproduce the issue, so I probably wouldn't pay more attention to me myself, unless something like this happens to show up elsewhere.