前言

本篇是非线性优化方法中的最后一个部分，LBFGS算法。

拟牛顿法

上篇记录了拟牛顿法的思路：通过迭代矩阵GkG_{k}Gk近似海森矩阵Hk−1H_k^{-1}Hk−1，满足拟牛顿条件：
Gk+1(J(xk+1)−J(xk))=xk+1−xkG_{k+1}(J({\bf x}_{k+1})-J({\bf x}_k))={\bf x}_{k+1} - {\bf x}_{k} \\ Gk+1(J(xk+1)−J(xk))=xk+1−xk
构造迭代表达式：
Gk+1=Gk+Pk+QkGk+1yk=sk=Gkyk+Pkyk+QkykG_{k+1}=G_k+P_k+Q_k \\ G_{k+1}y_k = s_k=G_ky_k+P_ky_k+Q_ky_k Gk+1=Gk+Pk+QkGk+1yk=sk=Gkyk+Pkyk+Qkyk

DFP算法用GkG_kGk近似Hk−1H_k^{-1}Hk−1：
Gk+1=Gk+skskTskTyk−GkykykTGkykTGkykG_{k+1} = G_k+ \frac{s_ks_k^T}{s_k^Ty_k} - \frac {G_ky_ky_k^TG_k}{y_k^TG_ky_k} \\ Gk+1=Gk+skTykskskT−ykTGkykGkykykTGk

BFGS算法用Bk−1B_k^{-1}Bk−1近似Hk−1H_k^{-1}Hk−1：
Bk+1−1=(I−skykTykTsk)Bk−1(I−ykskTykTsk)+skskTykTskB_{k+1}^{-1} = (I - \frac{s_ky_k^T}{y_k^Ts_k})B_k^{-1}(I-\frac{y_ks_k^T}{y^T_ks_k}) + \frac{s_ks_k^T}{y_k^Ts_k} Bk+1−1=(I−ykTskskykT)Bk−1(I−ykTskykskT)+ykTskskskT

迭代过程中，需要存储近似矩阵GkG_kGk或BkB_kBk，当输入维度很大时，这个近似矩阵所需的内存非常惊人，比如输入维度为10000，则存储GGG需要400MB，因此在大规模优化问题中，通常使用LBFGS算法来降低内存消耗

LBFGS算法

L-BFGS(low-memory BFGS)是BFGS的改进算法，主要思路是不存储BFGS迭代的Bk−1B_k^{-1}Bk−1，而是存储计算Bk−1B_k^{-1}Bk−1过程中的向量，通过这些向量恢复Bk−1B_k^{-1}Bk−1，从而节省内存空间。

将BFGS的迭代Bk+1−1B_{k+1}^{-1}Bk+1−1进行简写：
Bk+1−1=(I−skykTykTsk)Bk−1(I−ykskTykTsk)+skskTykTsksetρk=1ykTsk,Vk=(I−skykT),Bk−1=HkHk+1=VkHkVkT+ρkskskTB_{k+1}^{-1} = (I - \frac{s_ky_k^T}{y_k^Ts_k})B_k^{-1}(I-\frac{y_ks_k^T}{y^T_ks_k}) + \frac{s_ks_k^T}{y_k^Ts_k} \\ \quad \\ set \quad \rho_k = \frac{1}{y^T_ks_k}, \quad V_k=(I-s_ky_k^T), \quad B^{-1}_k=H_k \\ \quad \\ H_{k+1}=V_{k}H_kV^T_k+\rho_ks_ks_k^T \\ \quad \\ Bk+1−1=(I−ykTskskykT)Bk−1(I−ykTskykskT)+ykTskskskTsetρk=ykTsk1,Vk=(I−skykT),Bk−1=HkHk+1=VkHkVkT+ρkskskT

然后将Hk+1H_{k+1}Hk+1展开：
Hk+1=VkVk−1Hk−1Vk−1TVkT+Vkρk−1sk−1sk−1TVkT+ρkskskT=(VkVk−1…V0)H0(VkVk−1…V0)T+(VkVk−1…V1)ρ0s0s0T(VkVk−1…V1)T+(VkVk−1…V2)ρ1s1s1T(VkVk−1…V2)T…+(Vk)ρk−1sk−1sk−1T(Vk)T+ρkskskTH_{k+1}= V_kV_{k-1}H_{k-1}V_{k-1}^TV_k^T\\+V_k\rho_{k-1}s_{k-1}s_{k-1}^TV_k^T\\+\rho_ks_ks_k^T \\ \quad \\ = (V_kV_{k-1}\dots V_{0})H_0(V_kV_{k-1}\dots V_{0})^T \\ + (V_kV_{k-1}\dots V_{1})\rho_0s_0s_0^T(V_kV_{k-1}\dots V_{1})^T \\ + (V_kV_{k-1}\dots V_{2})\rho_1s_1s_1^T(V_kV_{k-1}\dots V_{2})^T \\ \dots \\ +(V_k)\rho_{k-1}s_{k-1}s_{k-1}^T(V_k)^T \\+\rho_ks_ks_k^T \\ \quad \\ Hk+1=VkVk−1Hk−1Vk−1TVkT+Vkρk−1sk−1sk−1TVkT+ρkskskT=(VkVk−1…V0)H0(VkVk−1…V0)T+(VkVk−1…V1)ρ0s0s0T(VkVk−1…V1)T+(VkVk−1…V2)ρ1s1s1T(VkVk−1…V2)T…+(Vk)ρk−1sk−1sk−1T(Vk)T+ρkskskT
如果每次迭代的中间向量si,yis_i,y_isi,yi，随着迭代次数的增加，内存消耗逐渐增大。因此，需要丢弃一部分中间向量，只保存最近的mmm次中间向量来作近似：
Hk+1=(VkVk−1…Vk−m+1)Hk−m+1(VkVk−1…Vk−m+1)T+(VkVk−1…Vk−m+2)ρk−m+1sk−m+1sk−m+1T(VkVk−1…Vk−m+2)T+…+(Vk)ρk−1sk−1sk−1T(Vk)T+ρkskskTH_{k+1} = \\ (V_kV_{k-1}\dots V_{k-m+1})H_{k-m+1}(V_kV_{k-1}\dots V_{k-m+1})^T \\ + (V_kV_{k-1}\dots V_{k-m+2})\rho_{k-m+1}s_{k-m+1}s_{k-m+1}^T(V_kV_{k-1}\dots V_{k-m+2})^T \\ + \dots \\ + (V_k)\rho_{k-1}s_{k-1}s_{k-1}^T(V_k)^T \\+\rho_ks_ks_k^T \\ \quad \\ Hk+1=(VkVk−1…Vk−m+1)Hk−m+1(VkVk−1…Vk−m+1)T+(VkVk−1…Vk−m+2)ρk−m+1sk−m+1sk−m+1T(VkVk−1…Vk−m+2)T+…+(Vk)ρk−1sk−1sk−1T(Vk)T+ρkskskT
将Hk−m+1H_{k-m+1}Hk−m+1设置为数量矩阵：
Hk−m+1=skTykykTykIH_{k-m+1}=\frac{s_k^Ty_k}{y_k^Ty_k}I Hk−m+1=ykTykskTykI

于是，LBFGS的基本思想到这就结束了。不过以上Hk+1H_{k+1}Hk+1的计算虽然节省了矩阵存储，但过程非常复杂，计算量大。于是LBFGS还提到了一种神奇的双循环算法来简化上面这个Hk+1H_{k+1}Hk+1计算过程：

gkg_kgk表示当前的梯度；最后得到的zzz，直接就获得了迭代增量−Hk+1g-H_{k+1}g−Hk+1g。

由此可以看出，LBFGS的操作实际上是一种以时间换空间的方法，通过增加计算过程，减少需要的存储空间。

代码示例

补充代码：仍然是用python写的lbfgs：

import numpy as np
import scipy.optimize
import time
import matplotlib.pyplot as pltdef partial_derivate_xy(x, y, F):dx = (F(x + 0.001, y) - F(x, y))/0.001dy = (F(x, y + 0.001) - F(x, y))/0.001return dx, dydef partial_partial_derivate_xy(x, y, F):dx, dy = partial_derivate_xy(x, y, F)dxx = (partial_derivate_xy(x + 0.001, y, F)[0] - dx) / 0.001dyy = (partial_derivate_xy(x, y + 0.001, F)[1] - dy) / 0.001dxy = (partial_derivate_xy(x, y + 0.001, F)[0] - dx) / 0.001dyx = (partial_derivate_xy(x + 0.001, y, F)[1] - dy) / 0.001return dxx, dyy, dxy, dyxdef non_linear_func(x, y):fxy = 0.5 * (x ** 2 + y ** 2)return fxydef non_linear_func_2(x, y):fxy = x*x + 2*y*y + 2*x*y + 3*x - y - 2return fxydef non_linear_func_3(x, y):fxy = 0.5 * (x ** 2 - y ** 2)return fxydef non_linear_func_4(x, y):fxy = x**4 + 2*y**4 + 3*x**2*y**2 + 4*x*y**2 + x*y + x + 2*y + 0.5return fxydef non_linear_func_5(x, y):fxy = (np.exp(x) + np.exp(0.5 * y) + x)**2return fxydef quasi_newton_LBFGS_optim(x, y, F, g_list, lr=0.01):dx, dy = partial_derivate_xy(x, y, F)grad = np.array([dx, dy])if len(g_list) == 0:G_0 = np.eye(2)z = np.matmul(G_0, grad)else:alpha_list = []q = gradfor i in np.arange(len(g_list) - 1, -1, -1):s_i = g_list[i][0, :]y_i = g_list[i][1, :]p_i = 1. / np.dot(y_i.T, s_i)alpha_i = p_i * np.dot(s_i.T, q)q = q - alpha_i * y_ialpha_list.append(alpha_i)G_0 = np.dot(g_list[-1][0, :].T, g_list[-1][1, :]) / np.dot(g_list[-1][1, :].T, g_list[-1][1, :]) * np.eye(2)z = np.matmul(G_0, q)for i in np.arange(0, len(g_list), 1):s_i = g_list[i][0, :]y_i = g_list[i][1, :]p_i = 1. / np.dot(y_i.T, s_i)beta_i = p_i * np.dot(y_i.T, z)z = z + s_i * (alpha_list[i] - beta_i)z_direction = z / np.linalg.norm(z)vec_delta = -lr * z_directionvec_opt = np.array([x, y]) + vec_deltax_opt = vec_opt[0]y_opt = vec_opt[1]dxx, dyy = partial_derivate_xy(x_opt, y_opt, F)grad_delta = np.array([dxx - dx, dyy - dy])if len(g_list) < 5:g_list.append(np.stack([vec_delta, grad_delta]))elif len(g_list) == 5:g_list.pop(0)g_list.append(np.stack([vec_delta, grad_delta]))return x_opt, y_opt, grad, g_listdef optimizer(x0, y0, F, th=0.01):x = x0y = y0g_list = []counter = 0while True:x_opt, y_opt, grad, g_list = quasi_newton_LBFGS_optim(x, y, F, g_list)if np.linalg.norm(grad) < th:breakx = x_opty = y_optcounter = counter + 1print('iter: {}'.format(counter), 'optimized (x, y) = ({}, {})'.format(x, y))# if counter==1000: breakreturn x, ydef verify_min(x, y, F):fx = F(x, y)deltax = np.linspace(-0.1, 0.1, 100)deltay = np.linspace(-0.1, 0.1, 100)x_range = x + deltaxy_range = y + deltaycounter = 0for i in range(100):for j in range(100):f_range = F(x_range[i], y_range[j])f_delta = fx - f_rangeif f_delta < 0:counter += 1print('counter: {}'.format(counter))def scipyF(X):x, y = Xfxy = x ** 4 + 2 * y ** 4 + 3 * x ** 2 * y ** 2 + 4 * x * y ** 2 + x * y + x + 2 * y + 0.5return fxyif __name__ == '__main__':x0 = 2.y0 = 2.start = time.time()for i in range(1000):result_x, result_y = optimizer(x0, y0, non_linear_func_4)if i == 0:breakend = time.time()print(result_x, result_y, 'cost time: {}'.format(end - start))verify_min(result_x, result_y, non_linear_func_4)xaxis = np.arange(-3, 3, 0.05)yaxis = np.arange(-3, 3, 0.05)[x_, y_] = np.meshgrid(xaxis, yaxis)f = non_linear_func_4(x_, y_)plt.contour(x_, y_, f, 20)plt.show()# scipyRes = scipy.optimize.fmin_cg(scipyF, np.array([0, 0]))# print(scipyRes)

在测试时我发现，如果像之前我写的牛顿法、拟牛顿法代码，不严格约束增量步长，LBFGS算法几乎都会发散。因此，下篇还将介绍一个梯度法使用中的步长约束原则。

后记

在无约束非线性优化中，我仅学习记录了梯度下降、牛顿法、拟牛顿法等各种方法的基本思想和具体操作，而没有记录它们的收敛条件推导以及步长条件。

等我后续在《最优化理论》中学习到相应部分后，再回来做更加深入的学习探索。

数值计算之 LBFGS相关推荐

数值计算之拟牛顿法
数值计算之拟牛顿法前言拟牛顿条件几何解释 DFP算法 BFGS算法 Broyden类算法代码示例后记前言无论是牛顿法,高斯牛顿法,还是LM算法,都需要计算海森矩阵或其近似,然后求海森矩 ...
大规模优化算法 - LBFGS算法
L-BFGS算法比较适合在大规模的数值计算中,具备牛顿法收敛速度快的特点,但不需要牛顿法那样存储Hesse矩阵,因此节省了大量的空间以及计算资源.本文主要通过对于无约束最优化问题的一些常用算法总结,一 ...
L-BFGS算法简介
参考大规模优化算法 - LBFGS算法上面文献的CSDN转载版因为在CRF的学习中需要用到LBFGS算法,所以就学习下. L-BFGS算法比较适合在大规模的数值计算中,具备牛顿法收敛速度快的特点 ...
优化算法SGD/ASGD/AdaGrad/Adadelta/RMSprop/Adam/Adamax/SparseAdam/L-BFGS/Rprop
机器学习界有一群炼丹师,他们每天的日常是: 拿来药材(数据),架起八卦炉(模型),点着六味真火(优化算法),就摇着蒲扇等着丹药出炉了. 不过,当过厨子的都知道,同样的食材,同样的菜谱,但火候不一样了, ...
【opencv】(2) 图像处理：边界填充、图像融合、图像阈值、数值计算
主要内容有:边界填充 cv2.copyMakeBorder(),数值计算 cv2.add(),改变尺寸 cv2.resize(),图像融合 cv2.addWeighted(),图像阈值 cv2.thr ...
Linux下数值计算
数值计算: Linux,在Linux执行最简单的加减乘除运算还是可以的,涉及到更高级的运算,我建议你还是算了吧在Linux中进行数值计算是非常不理想的,在默认情况下不支持数值计算,比如a=1;b=2 ...
数值计算领域的“圣经”，图灵出了新版本 | 11月书讯
简报本月图灵虽然只出品了 5 本新书,但其中有一本数学领域的重量级巨著,就是 Gene H. Golub 的<矩阵计算>(第4版). 此外,其他4本好书分别是:<C#并发编程经典 ...
送你38个常用的Python库，数值计算、可视化、机器学习等8大领域都有了
来源:大数据DT(ID:bigdatadt) 作者:李明江张良均周东平张尚佳内容摘编自<Python3智能数据分析快速入门> 本文约5200字,建议阅读10分钟. 本文为你总结了常 ...
sklearn警告：ConvergenceWarning: lbfgs failed to converge (status=1):
问题这个警告是训练逻辑回归模型的时候出来的. model=LogisticRegression() train_model("logistic regression",model ...

数值计算之 LBFGS

数值计算之 LBFGS

前言

拟牛顿法

LBFGS算法

代码示例

后记

数值计算之 LBFGS相关推荐

最新文章

热门文章