数值优化（六）——拟牛顿法之LBFGS理论与实战
对于无约束优化来说，线搜索算法的核心是解决这样一个迭代公式：https://www.zhihu.com/equation?tex=+x_%7Bk%2B1%7D+%3D+x_k+%2B+%5Calpha_k+p_k+%5Cquad+%280%29+
根据前几章的学习， https://www.zhihu.com/equation?tex=%5Calpha_k 可以通过Wolfe条件确定，而搜索方向则是不同算法的核心区别。
对于最速下降法来说，每次迭代都取负梯度方向，这样仅仅考虑一阶性质，不能从更宏观的范围考虑，因此收敛效率并不算高。
而牛顿法则利用黑塞矩阵计算 https://www.zhihu.com/equation?tex=p_k+%3D+-H_k%5E%7B-1%7D%5Cnabla+f_k ，能够达到非常高效的收敛，但它对问题要求非常苛刻，首先要求黑塞矩阵在迭代过程中至少是半正定的，这也就要求初始值便要非常接近最优解，其次黑塞矩阵是一个的矩阵，对于大规模问题来说如果黑塞矩阵不是稀疏矩阵，这空间代价可能无法接受，进而的计算也是问题。
拟牛顿法可以说是以上两者的折中，既能保持较为高效的收敛，同时又能保证鲁棒性。拟牛顿法不再直接计算黑塞矩阵，而是通过一定的方式构造一个近似的黑塞矩阵或者黑塞矩阵的逆，构造方式一定是简单容易的，并且一定保证近似矩阵是正定的，甚至不用显式存储这个矩阵。LBFGS是一个经典且流行的拟牛顿算法，在各个领域都表现出了强大的生命力，可以说是无约束优化领域的核心算法。
BFGS
LBFGS是由BFGS算法发展而来，BFGS是取自其四位发明者的首字母(Broyden, Fletcher, Goldfarb, Shanno)。
拟牛顿条件
根据中值定理，对于第 https://www.zhihu.com/equation?tex=k%2Ck%2B1 次迭代，有：
https://www.zhihu.com/equation?tex=+%5Cnabla+f%28x_%7Bk%2B1%7D%29+-+%5Cnabla+f%28x_k%29+%3D+%5Cnabla%5E2+f%28%5Chat%7Bx%7D%29%28x_%7Bk%2B1%7D-x_k%29+%5Cquad+%281%29+
其中 https://www.zhihu.com/equation?tex=%5Chat%7Bx%7D 是介于 https://www.zhihu.com/equation?tex=x_%7Bk%2B1%7D 和 https://www.zhihu.com/equation?tex=x_k 间的向量。
我们将 https://www.zhihu.com/equation?tex=x_%7Bk%2B1%7D-x_k 和 https://www.zhihu.com/equation?tex=%5Cnabla+f%28x_%7Bk%2B1%7D%29+-+%5Cnabla+f%28x_k%29 分别记为： https://www.zhihu.com/equation?tex=+s_k+%3D+x_%7Bk%2B1%7D-x_k+%5Cquad+%282%29+%5C%5C+y_k+%3D+%5Cnabla+f%28x_%7Bk%2B1%7D%29+-+%5Cnabla+f%28x_k%29+%5Cquad+%283%29+
那么如何让拟牛顿法构造的近似矩阵能够真的近似黑塞矩阵呢？一般会要求其满足拟牛顿条件，也就是可以带入公式(1)成立：
https://www.zhihu.com/equation?tex=+y_k+%3D+B_%7Bk%2B1%7Ds_k+%5Cquad+%284%29+
其中便是第 https://www.zhihu.com/equation?tex=k%2B1 次迭代的近似黑塞矩阵，或者直接近似黑塞矩阵的逆 https://www.zhihu.com/equation?tex=H_%7Bk%2B1%7D+%3D+B_%7Bk%2B1%7D%5E%7B-1%7D ，则有：
https://www.zhihu.com/equation?tex=+s_k+%3D+H_%7Bk%2B1%7Dy_k+%5Cquad+%285%29+
矩阵修正
对于BFGS算法，初始迭代时用单位矩阵 https://www.zhihu.com/equation?tex=I 作为黑塞矩阵的近似，然后每次迭代都是在上次迭代的近似矩阵上修正得到新的近似矩阵，修正方式也很简单：
https://www.zhihu.com/equation?tex=+B_%7Bk%2B1%7D+%3D+B_k+%2B+%5Calpha_k+u_k+u_k%5ET%2B%5Cbeta_k+v_k+v_k%5ET+%5Cquad+%286%29+
这种修正方式也叫秩2修正（秩1修正就只加一项，两种修正方式也是拟牛顿方法通用的），就是通过构造两个向量 https://www.zhihu.com/equation?tex=u_k%2C+v_k 来构造近似矩阵的变化量。而 https://www.zhihu.com/equation?tex=u_k%2Cv_k 的构造也很容易，简单来说就是将公式(6)带入拟牛顿条件然后凑出来(实际是有最优推导的，不展开介绍了)：
https://www.zhihu.com/equation?tex=+y_k+%3D+B_k+s_k+%2B+%5Calpha_k+u_k+u_k%5ET+s_k+%2B+%5Cbeta_k+v_k+v_k%5ET+s_k+%5Cquad+%287%29+
其中 https://www.zhihu.com/equation?tex=y_k%2C+B_k%2C+s_k 都是当前迭代已知量，凑的时候，不妨令 https://www.zhihu.com/equation?tex=y_k+%3D+%5Calpha_k+u_k+%5Cunderline%7Bu_k%5ET+s_k%7D ， https://www.zhihu.com/equation?tex=B_k+s_k+%2B+%5Cbeta_k+v_k+%5Cunderline%7Bv_k%5ET+s_k%7D%3D0 ，因为下划线标记部分都是标量，所以整理一下就有： https://www.zhihu.com/equation?tex=+u_k+%3D+y_k+%5C%5C+%5Calpha_k+%3D+%5Cfrac%7B1%7D%7By_k%5ET+s_k%7D+%5C%5C+v_k+%3D+B_k+s_k+%5C%5C+%5Cbeta_k+%3D+-%5Cfrac%7B1%7D%7Bs_k%5ET+B_k%5ET+s_k%7D+%5Cquad+%288%29+
带入公式(6)就有迭代公式：
https://www.zhihu.com/equation?tex=+B_%7Bk%2B1%7D%3D+B_k+%2B+%5Cfrac%7By_k+y_k%5ET%7D%7By_k%5ET+s_k%7D+-+%5Cfrac%7BB_k+s_k+s_k%5ET+B_k%5ET%7D%7Bs_k%5ET+B_k%5ET+s_k%7D+%5Cquad+%289%29+
通常用的是的逆，将公式(9)两边同时取逆，利用两次Sherman-Morrison公式( https://www.zhihu.com/equation?tex=%28A%2Buv%5ET%29%5E%7B-1%7D+%3D+A%5E%7B-1%7D-%5Cfrac%7BA%5E%7B-1%7Duv%5ETA%5E%7B-1%7D%7D%7B1%2Bv%5ETA%5E%7B-1%7Du%7D )，可以化简得到：
https://www.zhihu.com/equation?tex=+H_%7Bk%2B1%7D+%3D+%28I-%5Cfrac%7Bs_ky_k%5ET%7D%7By_k%5ETs_k%7D%29H_k%28I-%5Cfrac%7By_ks_k%5ET%7D%7By_k%5ETs_k%7D%29%2B%5Cfrac%7Bs_ks_k%5ET%7D%7By_k%5ETs_k%7D+%5Cquad+%2810%29+
这就是BFGS的近似矩阵的迭代公式，初始 https://www.zhihu.com/equation?tex=H_0+%3D+I ，就可以通过 https://www.zhihu.com/equation?tex=p_k+%3D+-H_k+%5Cnabla+f_k 求解搜索方向了。
LBFGS
根据公式(10)的近似矩阵可以看出，BFGS并没有解决牛顿法需要存储一个矩阵的缺点，因此后来提出了LBFGS算法，其中L便是Limited memory的含义，表示仅仅使用有限的内存空间。
从迭代公式可以看出，是通过前k次迭代的所有计算得到的，LBFGS的思想便是不再使用所有，而是保存最近m次的，这样仅仅需要保存2m个向量，然后每次迭代用这最近m次的计算出当前的来，并且整个计算过程是混合了的计算，避免了显式保存，因此降低了内存消耗，成就了它在数值优化上的地位。
对于公式(10)，不妨记： https://www.zhihu.com/equation?tex=+%5Crho_k+%3D+%5Cfrac%7B1%7D%7By_k%5ET+s_k%7D+%5C%5C+V_K+%3D+I-%5Crho_k+y_ks_k%5ET+
则简写为：
https://www.zhihu.com/equation?tex=+H_%7Bk%2B1%7D+%3D+%5Crho_k+s_k+s_k%5ET+%2B+V_k%5ET+H_k+V_k+%5Cquad+%2811%29+
对于右边 https://www.zhihu.com/equation?tex=H_k 可以再带入 https://www.zhihu.com/equation?tex=H_%7Bk-1%7D 一直递归到 https://www.zhihu.com/equation?tex=H_%7Bk-m%2B1%7D ，然后选 https://www.zhihu.com/equation?tex=H_%7Bk-m%2B1%7D%3D+I ，这样就通过最近m次迭代的得到了。
更巧妙的是，可以通过two-loop recursion方法，直接计算得到搜索方向，而不用显式计算矩阵。
q = -g;
for(int i = k; i > k - m; i--)
{
alpha = rho*s.dot(q); \\ 1
q = q - alpha*y;
}
r = q;
for(int i = k - m + 1; i < k + 1; i++)
{
beta = rho*y.dot(r);
r = r + s*(alpha-beta);
}
return r;
这段伪代码看起来很难理解，也很不可思议，但其实将公式(11)的递归展开，对照梳理一下就很容易理解了，公式(11)展开如下： https://www.zhihu.com/equation?tex=+-H_%7Bk%2B1%7D%5Cnabla+f_%7Bk%2B1%7D+%3D+-%5B%5Crho_k+s_k+s_k%5ET+%2B+V_k%5ET%28%5Crho_%7Bk-1%7Ds_%7Bk-1%7Ds_%7Bk-1%7D%5ET+%2B+V_%7Bk-1%7D%5ET%28%5Ccdots+%2BV_%7Bk-m%2B1%7D%5ET%28%5Crho_%7Bk-m%2B1%7Ds_%7Bk-m%2B1%7Ds_%7Bk-m%2B1%7D%5ET+%2B+V_%7Bk-m%2B1%7D%5ET+H_%7Bk-m%2B1%7DV_%7Bk-m%2B1%7D%29V_%7Bk-m%2B1%7D%29%5Ccdots+%29V_%7Bk-1%7D%29V_k%5D%5Cnabla+f_%7Bk%2B1%7D+
从右往左看，有一系列 https://www.zhihu.com/equation?tex=V_k 相乘，这其实就是代码(1)处所计算的，一直乘到 https://www.zhihu.com/equation?tex=H_%7Bk-m%2B1%7D%3DI ，然后再继续向左，就是边乘边加 https://www.zhihu.com/equation?tex=%5Crho+s+s%5ET 的部分，也就是第二个循环的内容，其实可以通过设定比如 https://www.zhihu.com/equation?tex=m%3D4 ，将整个式子完整写出来，就很清晰了。
就这样通过这个神奇的双循环，就将LBFGS的搜索方向计算出来了。可以看出，整个计算过程，只涉及到了向量的计算，避免的大规模矩阵，大大提高了算法的可用性，成为了一个经典而又流行的算法。
实战
同样的，完整的实现代码可以在github上的Chimes(https://github.com/BKHao/Chimes)中找到，其中lbfgs算法继承于线搜索算法，并且迭代框架是与最速下降法一样的，只是将迭代方向由 https://www.zhihu.com/equation?tex=p_k%3D-%5Cnabla+f_k 改为了用上述双循环求解。具体来说，有以下几点。
s = step * direction;
y = gradient - old_gradient;
old_gradient = gradient;
cursor = (cursor + 1) % Base::parameter_.lbfgs_remain_;
这里在每次迭代中及时更新最近m次的，m便是参数lbfgs_remain_。
direction = -gradient;
size_t j = cursor;
for (size_t i = 0; i < bound; ++i)
{
j = (j + Base::parameter_.lbfgs_remain_ - 1) % Base::parameter_.lbfgs_remain_;
alpha = s.dot(direction) / ys;
direction = direction - alpha * y;
}
direction = (_ys / _yy) * direction;
for (size_t i = 0; i < bound; ++i)
{
Scalar beta = y.dot(direction) / ys;
direction = direction + (alpha - beta) * s;
j = (j + 1) % Base::parameter_.lbfgs_remain_;
}
然后实现双循环求解搜索方向，这里与伪代码是完全一致的。相对最速下降法，其实就主要这两处改动，也可以看出对于线搜索算法而言，整体框架都是基本一致的，更具体的内容可以去Chimes中阅读。
比较
对于 https://www.zhihu.com/equation?tex=f%3D%28x-1%29%5E2%2B%28y-1%29%5E2 这个优化问题，最速下降法需要33次迭代才达到最优解，而LBFGS算法则仅仅需要2次迭代，可见其效率之高。对于更复杂的问题，像是优化得到下图点的分布(学术名称为CVT)，最速下降法和LBFGS的收敛过程可以从折线图中看出。
优化点的分布
最速下降法(蓝色)和LBFGS(黄色)收敛过程比较
其中蓝色线为最速下降法，黄色线为LBFGS，横轴为迭代次数，纵轴为梯度模长的对数，可以看出LBFGS的实际收敛速率明显好于最速下降法。对于理论上的收敛速率比较，可以留在之后将几类算法一起总结。 [赞][赞][赞]
页:
[1]