数值优化（六）——拟牛顿法之LBFGS理论与实战

JoshWindsor · 发表于 2022-6-6 08:07

对于无约束优化来说，线搜索算法的核心是解决这样一个迭代公式：

$x_{k+1} = x_k + \alpha_k p_k \quad (0)$
根据前几章的学习， $\alpha_k$ 可以通过Wolfe条件确定，而搜索方向  则是不同算法的核心区别。
对于最速下降法来说，每次迭代  都取负梯度方向，这样仅仅考虑一阶性质，不能从更宏观的范围考虑，因此收敛效率并不算高。
而牛顿法则利用黑塞矩阵计算 $p_k = -H_k^{-1}\nabla f_k$ ，能够达到非常高效的收敛，但它对问题要求非常苛刻，首先要求黑塞矩阵在迭代过程中至少是半正定的，这也就要求初始值便要非常接近最优解，其次黑塞矩阵是一个  的矩阵，对于大规模问题来说如果黑塞矩阵不是稀疏矩阵，这空间代价可能无法接受，进而  的计算也是问题。
拟牛顿法可以说是以上两者的折中，既能保持较为高效的收敛，同时又能保证鲁棒性。拟牛顿法不再直接计算黑塞矩阵，而是通过一定的方式构造一个近似的黑塞矩阵或者黑塞矩阵的逆，构造方式一定是简单容易的，并且一定保证近似矩阵是正定的，甚至不用显式存储这个矩阵。LBFGS是一个经典且流行的拟牛顿算法，在各个领域都表现出了强大的生命力，可以说是无约束优化领域的核心算法。
BFGS

LBFGS是由BFGS算法发展而来，BFGS是取自其四位发明者的首字母(Broyden, Fletcher, Goldfarb, Shanno)。
拟牛顿条件

根据中值定理，对于第 $k,k+1$ 次迭代，有：

$\nabla f(x_{k+1}) - \nabla f(x_k) = \nabla^2 f(\hat{x})(x_{k+1}-x_k) \quad (1)$
其中 $\hat{x}$ 是介于 $x_{k+1}$ 和 $x_k$ 间的向量。
我们将 $x_{k+1}-x_k$ 和 $\nabla f(x_{k+1}) - \nabla f(x_k)$ 分别记为： $s_k = x_{k+1}-x_k \quad (2) \\ y_k = \nabla f(x_{k+1}) - \nabla f(x_k) \quad (3)$
那么如何让拟牛顿法构造的近似矩阵能够真的近似黑塞矩阵呢？一般会要求其满足拟牛顿条件，也就是可以带入公式(1)成立：

$y_k = B_{k+1}s_k \quad (4)$
其中  便是第 $k+1$ 次迭代的近似黑塞矩阵，或者直接近似黑塞矩阵的逆 $H_{k+1} = B_{k+1}^{-1}$ ，则有：

$s_k = H_{k+1}y_k \quad (5)$
矩阵修正

对于BFGS算法，初始迭代时用单位矩阵 $I$ 作为黑塞矩阵的近似，然后每次迭代都是在上次迭代的近似矩阵上修正得到新的近似矩阵，修正方式也很简单：

$B_{k+1} = B_k + \alpha_k u_k u_k^T+\beta_k v_k v_k^T \quad (6)$
这种修正方式也叫秩2修正（秩1修正就只加一项，两种修正方式也是拟牛顿方法通用的），就是通过构造两个向量 $u_k, v_k$ 来构造近似矩阵的变化量。而 $u_k,v_k$ 的构造也很容易，简单来说就是将公式(6)带入拟牛顿条件然后凑出来(实际是有最优推导的，不展开介绍了)：

$y_k = B_k s_k + \alpha_k u_k u_k^T s_k + \beta_k v_k v_k^T s_k \quad (7)$
其中 $y_k, B_k, s_k$ 都是当前迭代已知量，凑的时候，不妨令 $y_k = \alpha_k u_k \underline{u_k^T s_k}$ ， $B_k s_k + \beta_k v_k \underline{v_k^T s_k}=0$ ，因为下划线标记部分都是标量，所以整理一下就有： $u_k = y_k \\ \alpha_k = \frac{1}{y_k^T s_k} \\ v_k = B_k s_k \\ \beta_k = -\frac{1}{s_k^T B_k^T s_k} \quad (8)$
带入公式(6)就有迭代公式：

$B_{k+1}= B_k + \frac{y_k y_k^T}{y_k^T s_k} - \frac{B_k s_k s_k^T B_k^T}{s_k^T B_k^T s_k} \quad (9)$
通常用的是  的逆  ，将公式(9)两边同时取逆，利用两次Sherman-Morrison公式( $(A+uv^T)^{-1} = A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}$ )，可以化简得到：

$H_{k+1} = (I-\frac{s_ky_k^T}{y_k^Ts_k})H_k(I-\frac{y_ks_k^T}{y_k^Ts_k})+\frac{s_ks_k^T}{y_k^Ts_k} \quad (10)$
这就是BFGS的近似矩阵的迭代公式，初始 $H_0 = I$ ，就可以通过 $p_k = -H_k \nabla f_k$ 求解搜索方向了。
LBFGS

根据公式(10)的近似矩阵可以看出，BFGS并没有解决牛顿法需要存储一个  矩阵的缺点，因此后来提出了LBFGS算法，其中L便是Limited memory的含义，表示仅仅使用有限的内存空间。
从迭代公式可以看出，  是通过前k次迭代的所有  计算得到的，LBFGS的思想便是不再使用所有  ，而是保存最近m次的  ，这样仅仅需要保存2m个向量，然后每次迭代用这最近m次的  计算出当前的  来，并且整个计算过程是混合了  的计算，避免了显式保存  ，因此降低了内存消耗，成就了它在数值优化上的地位。
对于公式(10)，不妨记： $\rho_k = \frac{1}{y_k^T s_k} \\ V_K = I-\rho_k y_ks_k^T$
则简写为：

$H_{k+1} = \rho_k s_k s_k^T + V_k^T H_k V_k \quad (11)$
对于右边 $H_k$ 可以再带入 $H_{k-1}$ 一直递归到 $H_{k-m+1}$ ，然后选 $H_{k-m+1}= I$ ，这样就通过最近m次迭代的  得到了  。
更巧妙的是，可以通过two-loop recursion方法，直接计算得到搜索方向，而不用显式计算矩阵。
q = -g[k+1];
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)展开如下： $-H_{k+1}\nabla f_{k+1} = -[\rho_k s_k s_k^T + V_k^T(\rho_{k-1}s_{k-1}s_{k-1}^T + V_{k-1}^T(\cdots +V_{k-m+1}^T(\rho_{k-m+1}s_{k-m+1}s_{k-m+1}^T + V_{k-m+1}^T H_{k-m+1}V_{k-m+1})V_{k-m+1})\cdots )V_{k-1})V_k]\nabla f_{k+1}$
从右往左看，有一系列 $V_k$ 相乘，这其实就是代码(1)处所计算的，一直乘到 $H_{k-m+1}=I$ ，然后再继续向左，就是边乘边加 $\rho s s^T$ 的部分，也就是第二个循环的内容，其实可以通过设定比如 $m=4$ ，将整个式子完整写出来，就很清晰了。
就这样通过这个神奇的双循环，就将LBFGS的搜索方向计算出来了。可以看出，整个计算过程，只涉及到了向量的计算，避免的大规模矩阵，大大提高了算法的可用性，成为了一个经典而又流行的算法。
实战

同样的，完整的实现代码可以在github上的Chimes(https://github.com/BKHao/Chimes)中找到，其中lbfgs算法继承于线搜索算法，并且迭代框架是与最速下降法一样的，只是将迭代方向由 $p_k=-\nabla f_k$ 改为了用上述双循环求解。具体来说，有以下几点。
s[cursor] = step * direction;
y[cursor] = 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[j] = s[j].dot(direction) / ys[j];
direction = direction - alpha[j] * y[j];
}
direction = (_ys / _yy) * direction;
for (size_t i = 0; i < bound; ++i)
{
Scalar beta = y[j].dot(direction) / ys[j];
direction = direction + (alpha[j] - beta) * s[j];
j = (j + 1) % Base::parameter_.lbfgs_remain_;
}
然后实现双循环求解搜索方向，这里与伪代码是完全一致的。相对最速下降法，其实就主要这两处改动，也可以看出对于线搜索算法而言，整体框架都是基本一致的，更具体的内容可以去Chimes中阅读。
比较

对于 $f=(x-1)^2+(y-1)^2$ 这个优化问题，最速下降法需要33次迭代才达到最优解，而LBFGS算法则仅仅需要2次迭代，可见其效率之高。对于更复杂的问题，像是优化得到下图点的分布(学术名称为CVT)，最速下降法和LBFGS的收敛过程可以从折线图中看出。

优化点的分布

最速下降法(蓝色)和LBFGS(黄色)收敛过程比较

其中蓝色线为最速下降法，黄色线为LBFGS，横轴为迭代次数，纵轴为梯度模长的对数，可以看出LBFGS的实际收敛速率明显好于最速下降法。对于理论上的收敛速率比较，可以留在之后将几类算法一起总结。

acecase · 发表于 2022-6-6 08:12

[赞][赞][赞]

		自动登录	找回密码
密码			立即注册

数值优化（六）——拟牛顿法之LBFGS理论与实战

本帖子中包含更多资源