使用SciPy最小化满足线性等式约束条件的二次函数

Question

使用SciPy最小化满足线性等式约束条件的二次函数

pythonnumpyscipymathematical-optimizationnlopt

10

我有一个相对简单的受限制优化问题，但是根据我采取的方式会得到不同的答案。让我们先导入必要的库和漂亮打印函数:

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

以下是受限优化问题（包括雅可比矩阵）的内容。简而言之，我要最小化的目标函数只是初始值到最终值的平方百分比变化之和。线性等式约束只要求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）'进行有限制最小化（带/不带jacobian）
使用'SLSQP'进行有限制最小化（带/不带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)都是合理的解决方案，它们可以实现显著较低的目标函数值, 从直觉上讲,我们会期望具有更大起始值的变量移动得更多(无论是绝对值还是百分比方面)。
将雅可比矩阵添加到 "trust-const" 中会导致它给出错误答案（或者至少是更差的答案），并且超过最大迭代次数。也许雅可比矩阵是错误的，但是功能非常简单，我相当确定它是正确的 (?)
'SLSQP'似乎无论是否提供雅可比矩阵都不能正常工作，但速度非常快，并声称成功终止。这似乎非常令人担忧，因为得到错误答案并声称已成功终止几乎是最坏的结果。
最初我使用了非常小的起始值和目标值（只有上文提到值的千分之一），在这种情况下，上面的所有5种方法都能正常工作并给出相同的答案。我的样本数据仍然非常小，似乎很奇怪它可以处理1,2,..,5但无法处理1000,2000,..5000。
值得一提的是，所有3个错误的结果都通过将每个初始值加1,000来达到目标 - 这符合约束条件，但与最小化平方百分比差异总和相去甚远（因为具有更高起始值的变量应该比较低的变量增加更多）。

所以我的问题实际上是这里发生了什么事情，为什么只有(1)和(2b)似乎有效？

更一般地说，我想找到一种解决此类优化问题的好的基于Python的方法，并且将考虑使用除scipy之外的其他软件包的答案，尽管最好的答案理想上也应该解决这里的scipy问题（例如，这是用户错误还是我应该在github上发布的错误？）。

- JohnE

对于无约束最小化问题，如果您明确设置 fatol=1e-8，您会得到什么结果？ - user545424

1

就此而言，我使用nlopt库尝试了您的示例，并使用SLSQP获得了正确的结果，因此可以排除您的雅可比函数或局部最小值存在问题的可能性。 - user545424

@user545424 谢谢！很高兴得到确认！ - JohnE

1

由于约束是线性的，其黑塞矩阵为零。这会导致对约束施加过多的权重吗？例如，如果使用黑塞矩阵的不精确估计值乘以雅可比矩阵的逆矩阵。 - Subhaneil Lahiri

1

CVXPY下提供了更好（凸）的QP求解器。 - Erwin Kalvelagen

显示剩余6条评论

2个回答

1

这是一个部分回答，我在此放置以防止问题变得更大，但我仍然希望看到更全面和说明性的答案。这些答案基于其他两个人的评论，但他们都没有完全编写出代码，我认为明确表述会更有意义，因此在这里给出：

修复2a（具有雅可比矩阵的信任约束）

关键似乎在于指定雅可比矩阵和海森矩阵中的两者或都不指定（而不是仅指定雅可比矩阵）。@SubhaneilLahiri发表了这样的评论，并且还有一条错误消息也是这样，我最初没有注意到：

UserWarning：delta_grad == 0.0。检查近似函数是否为线性函数。如果函数是线性函数，则可以通过将海森矩阵定义为零而不是使用拟牛顿逼近来获得更好的结果。

所以我通过定义海森矩阵函数来修复它：

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 )

- JohnE

1

关于您的第二个问题，我认为如果scipy默认将ftol设置为双精度机器精度会更好。或者，当您没有设置限制时，nlopt所做的是将其默认设置为零，然后通常会出现有关舍入误差的错误警告，这迫使用户设置合理的ftol。 - user545424

1

嘿John和@user545424，你们的评论和答案刚好解决了我几天来一直苦恼的问题（一直在撞墙），我非常感激。这全都是关于ftol的问题！ - NULL

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user545424 · Accepted Answer

这个问题可以使用一个非线性优化库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。