use of uninitialized variable in SLSQP when n equality constraints (Trac #1231)

[scipy-gitbot]

Original ticket http://projects.scipy.org/scipy/ticket/1231 on 2010-07-12 by trac user stevenj, assigned to unknown.

I was adapting the SLSQP code for use in another library (NLopt, ab-initio.mit.edu/nlopt), and I noticed a bug that you might want to fix. In the case where the number of equality constraints equals the dimension of the problem, valgrind reported the use of an uninitialized variable in a conditional statement. I tracked the source of the uninitialized value down to line 813 of slsqp_optmz.f in SUBROUTINE lsei:

  CALL dcopy_ (mg-mc,w(mc1),0,w(mc1),1)

I believe that this should be:

  CALL dcopy_ (mg,w(mc1),0,w(mc1),1)

This initializes part of the work array w to 0. If the number of inequality constraints differs from the number of dimensions (mc.NE.n in lsei), then this statement doesn't matter because w(mc1...mc1+mg-1) get overwritten anyway by the Lagrange multipliers when lsi is called a little later. But if mc.EQ.n, then the GOTO statment on line 815 jumps to the end, where w(mc1...mc1+mg-1) is used in a dot product.

I've verified that, with the above fix, it correctly solves a test problem with various numbers of equality and/or inequality constraints, and valgrind no longer complains.

PS. I should also mention that, because of rounding errors, SLSQP will occasionally evaluate the objective/constraints slightly outside of the bounding box. I'm not sure if you care about this; my NLopt library guarantees that the bound constraints are strictly honored (unlike nonlinear constraints), so I had to tweak SLSQP in a couple of places to check that rounding errors don't push the solution out of bounds.

[scipy-gitbot]

@charris wrote on 2010-07-12

Ah, test problems, my precious, we needs them.

[scipy-gitbot]

trac user stevenj wrote on 2010-07-15

I was running the code directly, not via SciPy, so my test case is not written in the SciPy syntax. The test problem I was using was the simple 2d problem described here: http://jdj.mit.edu/wiki/index.php/NLopt_Tutorial ...with the two nonlinear inequality constraints changed into equality constraints.

Of course, the thing about uninitialized-memory bugs is that it may still succeed as long as the uninitialized memory takes on innocuous enough values, so I wouldn't have caught this without valgrind. Of course, "there's a good chance this will succeed" is not particularly reassuring, which is why I tracked down and fixed the bug.

[scipy-gitbot]

trac user stevenj wrote on 2010-07-15

Correction: because SLSQP does not support unbounded variables, I changed that test problem to bound x1 in [-10,10] and x2 in [0,10].

(I've since modified my SLSQP version to support infinite bounds, but that's a separate patch.)

[rgommers]

OK, so http://jdj.mit.edu/wiki/index.php/NLopt_Tutorial#Example_in_Python needs to be translated in a unit test.

[andyfaff]

This looks as if it can be closed now. @rgommers, the suggested test (expressed as two equality constraints) runs on my machine:

# https://nlopt.readthedocs.io/en/latest/NLopt_Tutorial/
# implement two equality constraints, in R^2.

import numpy as np
from scipy.optimize import minimize

def fun(x):
    return np.sqrt(x[1])

def f_eqcon(x):
    """ Equality constraint """
    return x[1] - (2 * x[0])**3

def f_eqcon2(x):
    """ Equality constraint """
    return x[1] - (-x[0] + 1)**3

c1 = {'type': 'eq', 'fun': f_eqcon}
c2 = {'type': 'eq', 'fun': f_eqcon2}

res = minimize(fun, [8, 0.25], method='SLSQP',
                       constraints=[c1, c2], bounds=[(-0.5, 1), (0, 8)])

np.testing.assert_allclose(res.fun, 0.5443310539518)
np.testing.assert_allclose(res.x, [0.33333333, 0.2962963])
assert res.success

Can this be closed as is, or does the code above need to be included in a unit test?

[andyfaff]

Added the suggested test in PR #10708 . Seems to work.

[rgommers]

PR merged, closing. Thanks @andyfaff

[andyfaff]

From what I can tell of the NLopt slsqp.c code the dcopy_ call is made with the first argument simply being mg, rather than (mg-mc) (as we have in scipy). Thus, while the suggested example seems to work, the underlying Fortran code remains unchanged.

[ilayn]

After spending some time with this bug for a week for the C-rewrite, I can confirm that the correction mg-mc to mg is required. Otherwise, part of the multiplier vector remains uninitialized and causes wrong output. What is saving us on Python side is the initial buffer is sent with zeros in it

scipy/scipy/optimize/_slsqp_py.py

Line 351 in ab84560

w = zeros(len_w)

But in my C tests, I used uninitialized buffer and this triggers the bug. I changed this anyways instead of relying on the caller zero input.