Scipy.optimize.minimize SLSQP with linear constraints fails

Question

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
            Gradient evaluations: 6

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
            Gradient evaluations: 1

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.

Answer 1

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
            Gradient evaluations: 6

You could also use the the suggestion made in the linked article, which is to use a lambda expression with a second argument that has the desired default value:

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