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?
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/.