SciPy's trust-constr is ignoring my constraints

Question

I am doing some statistics-related optimization. I am trying to use scipy's minimize. However, the solution that it gives me is not feasible. My code is

import numpy as np
from scipy.optimize import minimize, NonlinearConstraint
from math import e

def get_expectation(f,mu):
    return f.dot(mu)

def get_neg_subgaus_coeff(f,mu):
    return -(e**((f/10)**2)).dot(mu)

def grad_neg_subgaus_coeff(f,mu):
    return -(2*e**((f/10)**2)*f)*mu

 def my_mu(n):
    return 1/(np.abs((np.arange(n)-n//2))+1)**2

def get_constraints(n):
    constr_list = []
    for i in range(n-1):

        constr = lambda f: f[i]- f[i+1]
        constr_list.append(NonlinearConstraint(constr,.1, .1))

    return constr_list

def callback(x,_): # not important, only debug info
    print(get_neg_subgaus_coeff(x, mu))
    return False

def get_solution(mu,n):
    constraints = get_constraints(n) + [NonlinearConstraint(lambda f: get_expectation(f,mu), 0, 0)]
    return minimize(get_neg_subgaus_coeff, np.zeros(n)-1, args=(mu),jac=grad_neg_subgaus_coeff, method="trust-constr", bounds=[(-100, 100) for _ in range(n)], constraints=constraints, tol=1e-8, options={"maxiter":1000, "disp":True}, callback=callback)

if __name__ == "__main__":
    n=10
    mu = my_mu(n)
    sol = get_solution(mu,n).x
    print(sol)

Now, the problem is that when the optimization finishes (with "xtol termination condition is satisfied." and "constraint violation: 8.33e-17"), it gives me the following solution

[-0.93668183 -0.63135599 -0.35261352 -0.54613932 -0.52013279  0.21760397
  1.25553898 -0.48837217 -2.04095876 -2.14095876]

which clearly does not satisfy the conditions produced by get_constraints which say that two adjacent elements of the array may differ by at most 0.1.

It behaves the same (it only runs slower) when I remove the Jacobian. There are several similar-sounding questions on stack overflow but those are about SLSQP and do not seem to be very relevant.

Any ideas why this is the case? Or how to solve it?

Answer 1

I think this related to a well-known "feature" of lambdas. If you look at:

 flist = [lambda x: x+i for i in range(3)]

then when executing:

  [flist[i](1) for i in range(3)]

you may expect [1+0,1+1,1+2] = [1,2,3]. However, you will see:

  [3, 3, 3]

Basically, there is just one lambda in this setup (or rather: late binding at end of loop).

One possible fix is:

 from functools import partial
 flist = [partial(lambda x,j: x+j, j=i) for i in range(3)]
 [flist[i](1) for i in range(3)]

This will show: [1,2,3].

When applying this technique to your problem, we can write:

def get_constraints(n):
    flist = [partial(lambda f,j: f[j]-f[j+1], j=i) for i in range(n-1)]
    return [NonlinearConstraint(flist[i],.1,.1) for i in range(n-1)]

Now the solution looks like:

x: array([ 0.49289571,  0.39289571,  0.29289571,  0.19289571,  0.09289571,
   -0.00710429, -0.10710429, -0.20710429, -0.30710429, -0.40710429])

Conclusion: nothing to do with the solver, but rather with a misunderstanding of how lambdas work.

PS. To be precise this is due to "late binding closures" (the same thing can happen without the use of a lambda function, but there is where I typically see this behavior). See https://docs.python-guide.org/writing/gotchas/.