BUG: Corruption/segfault in scipy.optimize.minimize with constraints

[i404788]

Describe your issue.

When minimizing with scipy.optimize.minimize it crashes python (kernel), with a corruption error.
This only happens when constraints are active with SLSQP.

I also tried scipy 1.7-1.6, but it has the same issue. Also found this issue which might be related: #14159

Reproducing Code Example

import numpy as np
import scipy.optimize as sco

x = np.random.rand(22,365)
target = np.linspace(0.9, 4.0, 50)

def metric(v, weights):
    return [[0, 0],[1, 1]]

def efficient_metric(v, target):
    def metric_a(weights):
        return metric(v, weights)[1][0]
    
    def metric_b(weights, v):
        return metric(v, weights)[0][0]

    constraints = ({'type': 'eq', 'fun': lambda x: metric_a(x) - target},
                   {'type': 'eq', 'fun': lambda x: np.sum(x) - 1})
    weights = np.array([len(v)*[1./len(v)]])
    result = sco.minimize(metric_b, weights, args=(v,),
                          method='SLSQP', constraints=constraints)
    return result

efficient_metric(x, target)

Error message

corrupted size vs. prev_size
OR
corrupted double-linked list
OR
double free or corruption (!prev)

SciPy/NumPy/Python version information

1.7.1 1.21.3 sys.version_info(major=3, minor=9, micro=2, releaselevel='final', serial=0)

[i404788]

Update I've located the issue, I think.
Changing target to 1 will result in a graceful failure, so a constraint with np.ndarray return type doesn't work.

[tupui]

Hi @i404788, I believe the constraints need to return a scalar value. I don't think there is an issue here.

[andyfaff]

It shouldn't crash python though.

[tupui]

This yes 👍 I am not sure if we can do proper input validation of the constraints before it goes to compiled code.

[andyfaff]

Probably

[i404788]

Yeah I kept it open so we can prevent it from crashing python (even if it's user error). I looked into it out of interest, there is already some input validation in place:

scipy/scipy/optimize/slsqp.py

Line 468 in 62eb424

c_eq = concatenate([atleast_1d(con['fun'](x, *con['args']))

, but it only checks if it's at least 1d, while we want it to be exactly one element.

It seems like the function is also wrapped depending on if it has a jacobian defined:

scipy/scipy/optimize/slsqp.py

Line 310 in 62eb424

cons[ctype] += ({'fun': con['fun'],

And it's used a few more times with different contexts, so maybe it's easiest to add another wrapper for input validation. For example:

def add_input_validation(f):
    def _f(*args, **kwargs):
        v = f(*args,**kwargs)
        v = atleast1d(v)
        assert np.prod(v.shape) == 1, "Constraints need to return a single scalar"
        return v
   return _f

[andyfaff]

For future reference, the constraints functions can return multiple values. For example:

import numpy as np
from scipy.optimize import minimize, rosen, NonlinearConstraint
from scipy.optimize._constraints import new_constraint_to_old

N = 4
rng = np.random.default_rng(1)
x0 = rng.uniform(size=N) * 10
bounds = [(0, 10)] * N

def con(x):
    # there are 3 inequality constraints here, and one equality
    return x[0] + x[1], x[2] + x[3], np.sum(x), x[0] - x[1]

nlc = NonlinearConstraint(con, [0.4, 2, 4, 0], [0.6, 7, 8, 0])

# warning: this is a private function and is subject to change
old_constraint = new_constraint_to_old(nlc, x0)

for c in old_constraint:
    cev = c['fun'](x0)
    print(len(cev), cev)

res0 = minimize(rosen, x0, bounds=bounds, constraints=nlc)
res1 = minimize(rosen, x0, bounds=bounds, constraints=old_constraint)
assert res0.success
assert res1.success
np.testing.assert_allclose(res0.x, res1.x)

gives

1 [-4.38642072]
6 [ 14.22285321   8.9280906   21.55094381 -14.02285321  -3.9280906
 -17.55094381]

[andyfaff]

Changing target to 1 will result in a graceful failure, so a constraint with np.ndarray return type doesn't work.

This is not True, constraint functions can return multiple values.

The original example is poor:

• a seed should be provided if there is random number generation
• metric_b will only ever return 0. Thus SLSQP is working in a perfectly flat landscape, there is no variation of metric_b with any variation of weights.
• the only function that would change value with different inputs is the second equality constraint.

Having said all that, it still causes a crash, and that shouldn't happen.

[i404788]

Hey @andyfaff,
The example was simplified for the bug report as the original was quite a bit larger and messier, but I verified the same issue occurred. All of the issues mentioned don't occur in the original code, so I don't think those are a cause.