使用 SciPy 最小化受线性等式约束的二次函数
Minimize quadratic function subject to linear equality constraints with SciPy
我有一个相当简单的约束优化问题，但根据我的操作方式会得到不同的答案。让我们先把导入和一个漂亮的打印函数放在一边：
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())
样本数据很简单：
n = 5
x0 = np.arange(1,6) * 10_000
target = x0.sum() + 5_000 # increase sum from 15,000 to 20,000
这是约束优化（包括雅可比矩阵）。换句话说，我要最小化的 objective 函数只是从初始值到最终值的平方百分比变化之和。线性 equality 约束只是要求 x.sum() 等于一个常数。
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)
为了进行比较，我通过使用等式约束将 x[0] 替换为 x[1:] 的函数来重构为无约束最小化。请注意，未约束函数传递 x0[1:] 而约束函数传递 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
然后我尝试通过三种方式最小化：
- 不受'Nelder-Mead'
的约束
- 受限于 'trust-constr'（w/ & w/o jacobian）
- 受限于 'SLSQP'（w/ & w/o jacobian）
代码：
##### (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' )
这里是一些简化的输出（当然你可以 运行 代码来获得完整的输出）：
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]
备注：
(1) & (2b) 是合理的解决方案，因为它们实现了显着较低的 objective 函数值，并且直觉上我们期望具有较大起始值的变量移动更多（两者绝对和百分比）比较小的。
将 jacobian 添加到 'trust-const' 会导致它得到错误的答案（或者至少是更糟糕的答案）并且还会超过最大迭代次数。也许雅可比是错误的，但功能是如此简单，我很确定它是正确的(?)
'SLSQP' 似乎无法使用提供的 jacobian 或 w/o，但运行速度非常快，并声称可以成功终止。这似乎非常令人担忧，因为得到错误的答案并声称已成功终止几乎是最糟糕的结果。
最初我使用了非常小的起始值和目标（只有我上面的 1/1,000），在那种情况下，上面的所有 5 种方法都可以正常工作并给出相同的答案。我的样本数据仍然非常小，它处理 1,2,..,5 而不是 1000,2000,..5000.
FWIW，请注意，3 个不正确的结果都通过向每个初始值添加 1,000 来达到目标——这满足了约束条件，但离最小化 objective 函数还差得很远（b/c 具有较高初始值的变量应该比初始值较低的变量增加更多，以最小化平方百分比差异之和）。
所以我的问题实际上就是这里发生了什么，为什么只有 (1) 和 (2b) 似乎有效？
更一般地说，我想找到一个很好的基于 python 的方法来解决这个问题和类似的优化问题，并且会考虑使用除 scipy 之外的其他包的答案，尽管最好的答案也是理想的在此处解决 scipy 的问题（例如，这是用户错误还是我应该 post 到 github 的错误？）。
以下是使用 nlopt 解决此问题的方法，这是一个非线性优化库，给我留下了深刻的印象。
首先，objective函数和梯度都是用同一个函数定义的：
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 期间，最小化器最终可能 space 某处，其中下一次跳跃距离最后一个地方小于 1e-8。
这是对我提出的问题的部分回答，以防止问题变得更大，但我仍然希望看到更全面和解释性更强的答案。这些答案基于其他两个人的评论，但他们都没有完整地写出代码，我认为明确说明是有意义的，所以这里是：
修复 2a（使用 jacobian 的信任构造）
这里关于 Jacobian 和 Hessian 的关键似乎是既不指定也不指定（但不是仅指定 jacobian）。 @SubhaneilLahiri 对此发表了评论，并且还有一条我最初没有注意到的错误消息：
UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.
所以我通过定义 hessian 函数修复了它：
def constr_hess(x,v):
return np.zeros([n,n])
并将其添加到约束条件
nonlin_con = NonlinearConstraint( constr_func, 0., 0., constr_jac, constr_hess )
修复 3a 和 3b (SLSQP)
这似乎只是按照@user545424 的建议缩小公差的问题。所以我只是将 options={'ftol':1e-15} 添加到最小化：
resSQjac = minimize( obj_func, x0, method='SLSQP',
options={'ftol':1e-15},
jac = obj_jac, constraints = eq_cons )
我有一个相当简单的约束优化问题，但根据我的操作方式会得到不同的答案。让我们先把导入和一个漂亮的打印函数放在一边：
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())
样本数据很简单：
n = 5
x0 = np.arange(1,6) * 10_000
target = x0.sum() + 5_000 # increase sum from 15,000 to 20,000
这是约束优化（包括雅可比矩阵）。换句话说，我要最小化的 objective 函数只是从初始值到最终值的平方百分比变化之和。线性 equality 约束只是要求 x.sum() 等于一个常数。
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)
为了进行比较，我通过使用等式约束将 x[0] 替换为 x[1:] 的函数来重构为无约束最小化。请注意，未约束函数传递 x0[1:] 而约束函数传递 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
然后我尝试通过三种方式最小化：
- 不受'Nelder-Mead' 的约束
- 受限于 'trust-constr'（w/ & w/o jacobian）
- 受限于 'SLSQP'（w/ & w/o jacobian）
代码：
##### (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' )
这里是一些简化的输出（当然你可以 运行 代码来获得完整的输出）：
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]
备注：
(1) & (2b) 是合理的解决方案，因为它们实现了显着较低的 objective 函数值，并且直觉上我们期望具有较大起始值的变量移动更多（两者绝对和百分比）比较小的。
将 jacobian 添加到 'trust-const' 会导致它得到错误的答案（或者至少是更糟糕的答案）并且还会超过最大迭代次数。也许雅可比是错误的，但功能是如此简单，我很确定它是正确的(?)
'SLSQP' 似乎无法使用提供的 jacobian 或 w/o，但运行速度非常快，并声称可以成功终止。这似乎非常令人担忧，因为得到错误的答案并声称已成功终止几乎是最糟糕的结果。
最初我使用了非常小的起始值和目标（只有我上面的 1/1,000），在那种情况下，上面的所有 5 种方法都可以正常工作并给出相同的答案。我的样本数据仍然非常小，它处理
1,2,..,5而不是1000,2000,..5000. 似乎有点奇怪
FWIW，请注意，3 个不正确的结果都通过向每个初始值添加 1,000 来达到目标——这满足了约束条件，但离最小化 objective 函数还差得很远（b/c 具有较高初始值的变量应该比初始值较低的变量增加更多，以最小化平方百分比差异之和）。
所以我的问题实际上就是这里发生了什么，为什么只有 (1) 和 (2b) 似乎有效？
更一般地说，我想找到一个很好的基于 python 的方法来解决这个问题和类似的优化问题，并且会考虑使用除 scipy 之外的其他包的答案，尽管最好的答案也是理想的在此处解决 scipy 的问题（例如，这是用户错误还是我应该 post 到 github 的错误？）。
以下是使用 nlopt 解决此问题的方法，这是一个非线性优化库，给我留下了深刻的印象。
首先，objective函数和梯度都是用同一个函数定义的：
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 期间，最小化器最终可能 space 某处，其中下一次跳跃距离最后一个地方小于 1e-8。
这是对我提出的问题的部分回答，以防止问题变得更大，但我仍然希望看到更全面和解释性更强的答案。这些答案基于其他两个人的评论，但他们都没有完整地写出代码，我认为明确说明是有意义的，所以这里是：
修复 2a（使用 jacobian 的信任构造）
这里关于 Jacobian 和 Hessian 的关键似乎是既不指定也不指定（但不是仅指定 jacobian）。 @SubhaneilLahiri 对此发表了评论，并且还有一条我最初没有注意到的错误消息：
UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.
所以我通过定义 hessian 函数修复了它：
def constr_hess(x,v):
return np.zeros([n,n])
并将其添加到约束条件
nonlin_con = NonlinearConstraint( constr_func, 0., 0., constr_jac, constr_hess )
修复 3a 和 3b (SLSQP)
这似乎只是按照@user545424 的建议缩小公差的问题。所以我只是将 options={'ftol':1e-15} 添加到最小化：
resSQjac = minimize( obj_func, x0, method='SLSQP',
options={'ftol':1e-15},
jac = obj_jac, constraints = eq_cons )