Joh*_*hnE 9 python numpy mathematical-optimization scipy quadratic-programming
I have a reasonably simple constrained optimization problem but get different answers depending on how I do it. Let's get the import and a pretty print function out of the way first:
import numpy as np
from scipy.optimize import minimize, LinearConstraint, NonlinearConstraint, SR1
def print_res( res, label ):
print("\n\n ***** ", label, " ***** \n")
print(res.message)
print("obj func value at solution", obj_func(res.x))
print("starting values: ", x0)
print("ending values: ", res.x.astype(int) )
print("% diff", (100.*(res.x-x0)/x0).astype(int) )
print("target achieved?",target,res.x.sum())
The sample data is very simple:
n = 5
x0 = np.arange(1,6) * 10_000
target = x0.sum() + 5_000 # increase sum from 15,000 to 20,000
Here's the constrained optimization (including jacobians). In words, the objective function I want to minimize is just the sum of squared percentage changes from the initial values to final values. The linear equality constraint is simply requiring x.sum() to equal a constant.
def obj_func(x):
return ( ( ( x - x0 ) / x0 ) ** 2 ).sum()
def obj_jac(x):
return 2. * ( x - x0 ) / x0 ** 2
def constr_func(x):
return x.sum() - target
def constr_jac(x):
return np.ones(n)
And for comparison, I've re-factored as an unconstrained minimization by using the equality constraint to replace x[0] with a function of x[1:]. Note that the unconstrained function is passed x0[1:] whereas the constrained function is passed x0.
def unconstr_func(x):
x_one = target - x.sum()
first_term = ( ( x_one - x0[0] ) / x0[0] ) ** 2
second_term = ( ( ( x - x0[1:] ) / x0[1:] ) ** 2 ).sum()
return first_term + second_term
I then try to minimize in three ways:
Code:
##### (1) unconstrained
res0 = minimize( unconstr_func, x0[1:], method='Nelder-Mead') # OK, but weird note
res0.x = np.hstack( [target - res0.x.sum(), res0.x] )
print_res( res0, 'unconstrained' )
##### (2a) constrained -- trust-constr w/ jacobian
nonlin_con = NonlinearConstraint( constr_func, 0., 0., constr_jac )
resTCjac = minimize( obj_func, x0, method='trust-constr',
jac='2-point', hess=SR1(), constraints = nonlin_con )
print_res( resTCjac, 'trust-const w/ jacobian' )
##### (2b) constrained -- trust-constr w/o jacobian
nonlin_con = NonlinearConstraint( constr_func, 0., 0. )
resTC = minimize( obj_func, x0, method='trust-constr',
jac='2-point', hess=SR1(), constraints = nonlin_con )
print_res( resTC, 'trust-const w/o jacobian' )
##### (3a) constrained -- SLSQP w/ jacobian
eq_cons = { 'type': 'eq', 'fun' : constr_func, 'jac' : constr_jac }
resSQjac = minimize( obj_func, x0, method='SLSQP',
jac = obj_jac, constraints = eq_cons )
print_res( resSQjac, 'SLSQP w/ jacobian' )
##### (3b) constrained -- SLSQP w/o jacobian
eq_cons = { 'type': 'eq', 'fun' : constr_func }
resSQ = minimize( obj_func, x0, method='SLSQP',
jac = obj_jac, constraints = eq_cons )
print_res( resSQ, 'SLSQP w/o jacobian' )
Here is some simplified output (and of course you can run the code to get the full output):
starting values: [10000 20000 30000 40000 50000]
***** (1) unconstrained *****
Optimization terminated successfully.
obj func value at solution 0.0045454545454545305
ending values: [10090 20363 30818 41454 52272]
***** (2a) trust-const w/ jacobian *****
The maximum number of function evaluations is exceeded.
obj func value at solution 0.014635854609684874
ending values: [10999 21000 31000 41000 51000]
***** (2b) trust-const w/o jacobian *****
`gtol` termination condition is satisfied.
obj func value at solution 0.0045454545462939935
ending values: [10090 20363 30818 41454 52272]
***** (3a) SLSQP w/ jacobian *****
Optimization terminated successfully.
obj func value at solution 0.014636111111111114
ending values: [11000 21000 31000 41000 51000]
***** (3b) SLSQP w/o jacobian *****
Optimization terminated successfully.
obj func value at solution 0.014636111111111114
ending values: [11000 21000 31000 41000 51000]
Notes:
(1) & (2b) are plausible solutions in that they achieve significantly lower objective function values and intuitively we'd expect the variables with larger starting values to move more (both absolutely and in percentage terms) than the smaller ones.
Adding the jacobian to 'trust-const' causes it to get the wrong answer (or at least a worse answer) and also to exceed max iterations. Maybe the jacobian is wrong, but the function is so simple that I'm pretty sure it's correct (?)
'SLSQP' doesn't seem to work w/ or w/o the jacobian supplied, but works very fast and claims to terminate successfully. This seems very worrisome in that getting the wrong answer and claiming to have terminated successfully is pretty much the worst possible outcome.
Initially I used very small starting values and targets (just 1/1,000 of what I have above) and in that case all 5 approaches above work fine and give the same answers. My sample data is still extremely small, and it seems kinda bizarre for it to handle 1,2,..,5 but not 1000,2000,..5000.
FWIW, note that the 3 incorrect results all hit the target by adding 1,000 to each initial value -- this satisfies the constraint but comes nowhere near minimizing the objective function (b/c variables with higher initial values should be increased more than lower ones to minimize the sum of squared percentage differences).
So my question is really just what is happening here and why do only (1) and (2b) seem to work?
More generally, I'd like to find a good python-based approach to this and similar optimization problems and will consider answers using other packages besides scipy although the best answer would ideally also address what is going on with scipy here (e.g. is this user error or a bug I should post to github?).
这是使用nlopt非线性解决方案优化库解决的方法，对此我印象深刻。
首先，使用相同的函数定义目标函数和梯度：
def obj_func(x, grad):
if grad.size > 0:
grad[:] = obj_jac(x)
return ( ( ( x/x0 - 1 )) ** 2 ).sum()
def obj_jac(x):
return 2. * ( x - x0 ) / x0 ** 2
def constr_func(x, grad):
if grad.size > 0:
grad[:] = constr_jac(x)
return x.sum() - target
def constr_jac(x):
return np.ones(n)
然后，使用Nelder-Mead和SLSQP运行最小化：
opt = nlopt.opt(nlopt.LN_NELDERMEAD,len(x0)-1)
opt.set_min_objective(unconstr_func)
opt.set_ftol_abs(1e-15)
xopt = opt.optimize(x0[1:].copy())
xopt = np.hstack([target - xopt.sum(), xopt])
fval = opt.last_optimum_value()
print_res(xopt,fval,"Nelder-Mead");
opt = nlopt.opt(nlopt.LD_SLSQP,len(x0))
opt.set_min_objective(obj_func)
opt.add_equality_constraint(constr_func)
opt.set_ftol_abs(1e-15)
xopt = opt.optimize(x0.copy())
fval = opt.last_optimum_value()
print_res(xopt,fval,"SLSQP w/ jacobian");
结果如下：
***** Nelder-Mead *****
obj func value at solution 0.00454545454546
result: 3
starting values: [ 10000. 20000. 30000. 40000. 50000.]
ending values: [10090 20363 30818 41454 52272]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0
***** SLSQP w/ jacobian *****
obj func value at solution 0.00454545454545
result: 3
starting values: [ 10000. 20000. 30000. 40000. 50000.]
ending values: [10090 20363 30818 41454 52272]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0
当进行测试时，我想我发现了最初尝试的问题所在。如果我1e-8对scipy函数默认设置的功能设置绝对公差，则会得到：
***** Nelder-Mead *****
obj func value at solution 0.0045454580693
result: 3
starting values: [ 10000. 20000. 30000. 40000. 50000.]
ending values: [10090 20363 30816 41454 52274]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0
***** SLSQP w/ jacobian *****
obj func value at solution 0.0146361108503
result: 3
starting values: [ 10000. 20000. 30000. 40000. 50000.]
ending values: [10999 21000 31000 41000 51000]
% diff [9 5 3 2 2]
target achieved? 155000.0 155000.0
这正是您所看到的。因此，我的猜测是，在SLSQP期间，最小化器最终会出现在似然空间中的某个位置，其中下一个跳转小于1e-8最后一个跳转。