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

Question

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

没有平凡的解决方案。

为此，我需要使用梯度为 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 优化器（这是最快的优化器，因为我们可以访问梯度），然后恢复到原始解决方案解决。注意函数和梯度里面的 - 符号是因为原问题的最小化等于对偶问题的最大化，所以它们是为了方程的实际最大化。我在这里绘制了优化器关于 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。因此所有解都应该在优化器的 return 处给出 0。

错误b'ABNORMAL_TERMINATION_IN_LNSRCH来自一个不好的梯度，但是在这里，它是函数的真实梯度吗...

我错过了什么？

Answer 1

您的代码中有多个错误标志，顺便说一下，这违反了 DRY principle。它应该是 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.  ]])

随心所欲。

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

Lagrange multiplier in L-BFGS-B

python

python-3.x

scipy

mathematical-optimization