python - L-BFGS-B 中的拉格朗日乘数

问题描述

希望做一个更难的问题，我试图从一个更简单的角度用下面的玩具最小化方程来解决它

这没有一个简单的解决方案。

为此，我需要使用具有梯度 2 和平衡点 3 的增强拉格朗日/对偶函数 1。上一个问题的增强拉格朗日版本：

拉格朗日乘数的重点是优化 [0,inf] 中的 mu，以便考虑我们问题的奇怪约束。

运行以下代码

a = 1
nbtests = 5 ; minmu = 0 ; maxmu = 5
def dual(mu) :
    x =  spo.fsolve(lambda x : 2*x  -  mu * (np.exp(x)+1) , 1)
    return   (- (x**2 - mu *(np.exp(x) + x -a )) ,- (np.exp(x) + x-a))

pl = np.empty((nbtests,2))
for i,nu in enumerate(np.linspace(minmu,maxmu,pl.shape[0]))  :
    res = spo.fmin_l_bfgs_b(dual,nu,bounds=[(0,None)],factr=1e6)
    print(nu,res[0],res[2]['task'])
    pl[i] = [nu,spo.fsolve(lambda x :  2*x  -  res[0] * (np.exp(x)+1), 1)]

plt.plot(pl[:,0],pl[:,1],label="L-BFGS-B")

plt.axhline(y=spo.fsolve(lambda x : np.exp(x)+x-a, 1)[0], label="Target",color="red")
plt.legend()
plt.show()  

plt.plot(pl[:,0],pl[:,1],label="L-BFGS-B")

plt.axhline(y=spo.fsolve(lambda x : np.exp(x)+x-a, 1)[0], label="Target",color="red")
plt.legend()
plt.show()

在这里，我们尝试从对偶问题中用 Python 中的 L-BFGS-B 优化器（这是最快的，因为我们可以访问梯度），然后使用 fsolve 恢复到原始解决方案。注意，函数和梯度内部的 - 符号是因为原始问题的最小化等于对偶问题的最大化，所以它们用于方程的实际最大化。我在这里绘制了优化器相对于 mu 的初始猜测的结果，并表明它非常依赖于初始化，根本不可靠。当将参数a更改为其他值时，该方法的收敛性更差。

相对于 mu 的优化器结果图（对于更大的 nbtests）

和程序的输出

0.0 [0.] b'ABNORMAL_TERMINATION_IN_LNSRCH'
0.5555555555555556 [0.55555556] b'ABNORMAL_TERMINATION_IN_LNSRCH'
1.1111111111111112 [3.52870269] b'ABNORMAL_TERMINATION_IN_LNSRCH'
1.6666666666666667 [3.52474085] b'ABNORMAL_TERMINATION_IN_LNSRCH'
2.2222222222222223 [3.5243099] b'ABNORMAL_TERMINATION_IN_LNSRCH'
2.7777777777777777 [3.49601967] b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
3.3333333333333335 [3.52020875] b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
3.8888888888888893 [3.88888889] b'ABNORMAL_TERMINATION_IN_LNSRCH'
4.444444444444445 [4.44444354] b'ABNORMAL_TERMINATION_IN_LNSRCH'
5.0 [5.] b'ABNORMAL_TERMINATION_IN_LNSRCH'

第一列是初始猜测，第二列是优化器之后的估计猜测，我们看到它对大多数值根本没有优化。对于 a <= 0 ，函数的域是 x < 0，对于 mu = 0，它给出了一个微不足道的最小化 x^2 = 0。所以所有解决方案在优化器返回时都应该给出 0。

错误 b'ABNORMAL_TERMINATION_IN_LNSRCH 来自错误的梯度，但在这里，它是函数的真正梯度吗...

我错过了什么？

标签： pythonpython-3.xscipymathematical-optimization

您的代码中有多个错误标志，顺便说一句，这违反了DRY 原则。对于导数，它应该是x**2 + mu * (np.exp(x) + x - a)代替和相似的。x**2 - mu * (np.exp(x) + x - a)恕我直言，类似

from scipy.optimize import fsolve, fmin_l_bfgs_b

a = 1
nbtests = 5
minmu = 0
maxmu = 5

def lagrange(x, mu):
    return x**2 + mu * (np.exp(x) + x - a)

def lagrange_grad(x, mu):
    grad_x = 2*x + mu * (np.exp(x) + 1)
    grad_mu = np.exp(x) + x - a
    return grad_x, grad_mu

def dual(mu):
    x = fsolve(lambda x: lagrange_grad(x, mu)[0], x0=1)
    obj_val = lagrange(x, mu)
    grad = lagrange_grad(x, mu)[1]
    return -1.0*obj_val, -1.0*grad

pl = np.empty((nbtests, 2))
for i, nu in enumerate(np.linspace(minmu,maxmu,nbtests)):
    res = fmin_l_bfgs_b(dual, x0=nu, bounds=[(0,None)], factr=1e6)
    mu_opt = res[0]
    x_opt = fsolve(lambda x: lagrange_grad(x, mu_opt)[0], x0=1)
    pl[i] = [nu, *x_opt]

干净多了。这给了我

array([[0.  , 0.  ],
       [1.25, 0.  ],
       [2.5 , 0.  ],
       [3.75, 0.  ],
       [5.  , 0.  ]])

如预期的。

python - L-BFGS-B 中的拉格朗日乘数

问题描述

解决方案

推荐阅读