# Scipy.optimize.minimize SLSQP with linear constraints fails

Posted on

### Question :

Scipy.optimize.minimize SLSQP with linear constraints fails

Consider the following (convex) optimization problem:

``````minimize 0.5 * y.T * y
s.t.     A*x - b == y
``````

where the optimization (vector) variables are `x` and `y` and `A`, `b` are a matrix and vector, respectively, of appropriate dimensions.

The code below finds a solution easily using the `SLSQP` method from Scipy:

``````import numpy as np
from scipy.optimize import minimize

# problem dimensions:
n = 10 # arbitrary integer set by user
m = 2 * n

# generate parameters A, b:
np.random.seed(123) # for reproducibility of results
A = np.random.randn(m,n)
b = np.random.randn(m)

# objective function:
def obj(z):
vy = z[n:]
return 0.5 * vy.dot(vy)

# constraint function:
def cons(z):
vx = z[:n]
vy = z[n:]
return A.dot(vx) - b - vy

# constraints input for SLSQP:
cons = ({'type': 'eq','fun': cons})

# generate a random initial estimate:
z0 = np.random.randn(n+m)

sol = minimize(obj, x0 = z0, constraints = cons, method = 'SLSQP',  options={'disp': True})
``````
``````Optimization terminated successfully.    (Exit mode 0)
Current function value: 2.12236220865
Iterations: 6
Function evaluations: 192
``````

Note that the constraint function is a convenient ‘array-output’ function.

Now, instead of an array-output function for the constraint, one could in principle use an equivalent set of ‘scalar-output’ constraint functions (actually, the scipy.optimize documentation discusses only this type of constraint functions as input to `minimize`).

Here is the equivalent constraint set followed by the output of `minimize` (same `A`, `b`, and initial value as the above listing):

``````# this is the i-th element of cons(z):
def cons_i(z, i):
vx = z[:n]
vy = z[n:]
return A[i].dot(vx) - b[i] - vy[i]

# listable of scalar-output constraints input for SLSQP:
cons_per_i = [{'type':'eq', 'fun': lambda z: cons_i(z, i)} for i in np.arange(m)]

sol2 = minimize(obj, x0 = z0, constraints = cons_per_i, method = 'SLSQP', options={'disp': True})
``````
``````Singular matrix C in LSQ subproblem    (Exit mode 6)
Current function value: 6.87999270692
Iterations: 1
Function evaluations: 32
``````

Evidently, the algorithm fails (the returning objective value is actually the objective value for the given initialization), which I find a bit weird. Note that running `[cons_per_i[i]['fun'](sol.x) for i in np.arange(m)]` shows that `sol.x`, obtained using the array-output constraint formulation, satisfies all scalar-output constraints of `cons_per_i` as expected (within numerical tolerance).

I would appreciate if anyone has some explanation for this issue.

You’ve run into the “late binding closures” gotcha. All the calls to `cons_i` are being made with the second argument equal to 19.

A fix is to use the `args` dictionary element in the dictionary that defines the constraints instead of the lambda function closures:

``````cons_per_i = [{'type':'eq', 'fun': cons_i, 'args': (i,)} for i in np.arange(m)]
``````

With this, the minimization works:

``````In [417]: sol2 = minimize(obj, x0 = z0, constraints = cons_per_i, method = 'SLSQP', options={'disp': True})
Optimization terminated successfully.    (Exit mode 0)
Current function value: 2.1223622086
Iterations: 6
Function evaluations: 192
``````cons_per_i = [{'type':'eq', 'fun': lambda z, i=i: cons_i(z, i)} for i in np.arange(m)]