牛顿法，拟牛顿法，BFGS，LBFGS

muyuu

已于 2024-07-07 17:03:59 修改

阅读量470

点赞数

分类专栏：数学工具文章标签：线性代数

于 2021-09-04 12:44:42 首次发布

本文链接：https://blog.csdn.net/muyuu/article/details/120098301

版权

数学工具专栏收录该内容

6 篇文章 0 订阅

订阅专栏

本文深入探讨了牛顿法及其在高维优化中的局限性，进而引出了拟牛顿法的概念。详细介绍了BFGS和LBFGS两种经典的拟牛顿法优化算法，阐述了它们的基本思想、计算流程及优势，并通过数学推导验证了算法的正确性。BFGS利用近似Hessian矩阵加速优化，而LBFGS则通过有限存储方式降低了内存需求，适合大规模问题求解。

摘要由CSDN通过智能技术生成

牛顿法

首先介绍一下牛顿法。给定一个优化问题
$min_x f(x)$
其中 $f (x)$ 具有连续的二阶偏导。

由Taylor展开，当 $x\rightarrow x_k$ ，
$f(x^*) = f(x_k)+\nabla f(x_k)(x^*-x_k)+\frac{1}{2}(x^*-x_k)^T\nabla^2f(x_k)(x^*-x_k)+ O(x^3)$
在极值点有 $\nabla f(x^*)=0$ ，因此 $\nabla f(x^*)=\frac{f(x^*)-f(x_k)}{x^*-x_k} = \nabla f(x_k) + \nabla^2f(x_k)(x^*-x_k)=0$
$\Longrightarrow \quad x^*=x_k-H_k^{-1}g_k$
其中 $H_k=\nabla^2f(x_k), \quad g_k=\nabla f(x_k)$ （Note：因为此处是求极小值，而当 $H_k$ 正定时， $f (x)$ 的极值为最小值，所以此处不妨假设 $H_k$ 是正定的，即 我们假设 $f (x)$ 是凸函数）

由此我们得到了牛顿法的下降方向 $H_k^{-1}g_k$ ，下面证明总 $\exists \lambda>0$ ，使得 $p_k=-\lambda H_k^{-1}g_k$ 满足： $f(x_k+p_k)<f(x_k)$

由Taylor展开， $f(x_{k+1})=f(x_k)-\lambda g_k H_k^{-1}g_k+ \frac{1}{2}\lambda^2(H_k^{-1}g_k)^T H_k (H_k^{-1}g_k)$ ，因为 $H_k$ 正定，所以 $\forall a\in \mathbb{R}^{n\times 1}, a^THa>0$ ，因此总存在一个充分小的 $\lambda_k$ ，使得 $f(x_{k+1})<f(x_k)$

Note: 牛顿法由于使用了二阶信息，可以证明它是具有二阶收敛性的，因此它的收敛速度会比作为一阶算法的梯度下降法快很多，但缺点是 $H_k^{-1}$ 的计算量很大，而且由于存在求逆的步骤，如果 $H_k$ 的条件数很大的话，容易在数值计算中出现NAN

拟牛顿法

针对上述问题，实际应用中，我们通常会找一个正定矩阵 $G_k$ 来拟合 $H_k^{-1}$ ，或者 $B_k$ 来拟合 $H_k$ ，下面我们来看牛顿法中的 $H_k$ 应该满足什么条件。

首先根据上述牛顿法，我们可以得到 $g_{k}=g_{k+1}+H_{k+1}(x_k-x_{k+1})$ ，记 $y_k=g_{k+1}-g_k ，s_k=x_{k+1}-x_k \Longrightarrow H_{k+1} s_k=y_k$
由此我们得到了拟牛顿条件： $H_{k+1} s_k=y_k$ 。并且要求 $H_{k+1}$ 是一个对称正定矩阵。

BFGS，LBFGS

BFGS

BFGS的思路是构造 $B_k$ 来拟合 $H_k$ ，由于 $H_k$ 是正定矩阵，因此在线性空间中 $H_k$ 是对称矩阵，所以我们首先将 $B_{k+1}$ 构造为对称阵： $B_{k+1}=B_k+\alpha u u^T+\beta v v^T$ ，那么根据拟牛顿条件，有：
$\begin{aligned} B_{k+1}s_k &= B_ks_k+\alpha u u^Ts_k+\beta v v^Ts_k\\ &= B_ks_k+\alpha u^Ts_k u+\beta v^Ts_k v \end{aligned}$
不妨令： $\alpha u^Ts_k=1,\beta v^Ts_k=-1$ ，则 $\alpha=\frac{1}{u^Ts_k}, \beta= -\frac{1}{v^Ts_k}$ 。此时：
$B_ks_k+u-v=y_k$
再令 $u=y_k, v=B_k s_k$ ， $\Longrightarrow \alpha=\frac{1}{y_k^T s_k}, \beta= -\frac{1}{s_k^TB_k s_k}$
$\Longrightarrow\quad B_{k+1} = B_k + \frac{y_k y_k^T}{y_k^T s_k}-\frac{B_k s_k s_k^TB_k}{s_k^T B_k s_k}$

但注意到这个式子里依然需要求 $B_k^{-1}$ ，因此还需要利用Sherman-Morrison公式：
$(A+uv^T)^{-1} = A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}$
其中 $A$ 是n阶可逆矩阵， $A+uv^T$ 也是可逆矩阵。

引用两次Sherman-Morrison公式可得：(具体推导参考BFGS公式推导)
$G_{k+1} = (1-\frac{s_k y_k^T}{s_k^T y_k})G_k(1-\frac{s_k y_k^T}{s_k^T y_k})^T+\frac{s_k s_k^T}{s_k^T y_k}$
其中 $G_k = 1/B_k$ 。

LBFGS

下面用 $H_k$ 来代替上式的 $G_k$ 。当维度比较大的时候，矩阵 $H$ 的储存开销很大，因此实际应用中，通常使用的是LBFGS（Limited-Memory BFGS）：其基本思路是，计算过程中不储存 $H_k$ ，而是储存向量序列 ${y_k\}$ 和 ${s_k\}$ 。并且向量序列 ${y_k\}，\{s_k\}$ 也不是全部储存，而是固定存最新的 $m$ 个，需要 $H_k$ 时，再用最新的m个 $y_k$ 和 $s_k$ 计算。这样一来，空间复杂度从 $O(N^2)$ 降到了 $O (m N)$ 。

记 $\rho_k = \frac{1}{y_k^T s_k}，V_k=I-\rho_k y_k s_k^T$ ，那么上式可以写成：
$H_{k+1} = V_k^TH_kV_k + \rho_k s_k s_k^T$

$H_0$ 初始化为 $I$ ，那么：
$\begin{aligned} H_1 = & V_0^TH_0V_0 + \rho_0 s_0 s_0^T\\ \quad\\ H_2 = & V_1^TH_1V_1 + \rho_1 s_1 s_1^T\\ = & V_1^T(V_0^TH_0V_0 + \rho_0 s_0 s_0^T)V_1 + \rho_1 s_1 s_1^T\\ = & V_1^TV_0^TH_0V_0V_1 + V_1^T\rho_0 s_0 s_0^TV_1 + \rho_1 s_1 s_1^T\\ \cdots\\ H_{k+1} = &(V_k^TV_{k-1}^T\cdots V_1^T V_0^T)H_0(V_0V_1\cdots V_{k-1}V_k) \\ & + (V_k^TV_{k-1}^T\cdots V_1^T)\rho_0 s_0 s_0^T(V_1\cdots V_{k-1}V_k) \\ & \cdots \\ & + V_k^TV_{k-1}^T\rho_{k-2} s_{k-2} s_{k-2}^T V_{k-1}V_k \\ & + V_k^T \rho_{k-1} s_{k-1} s_{k-1}^TV_k \\ & + \rho_k s_k s_k^T \tag{LBFGS[1]} \end{aligned}$
只保留最近的m步：
$\begin{aligned} H_{k+1} = &(V_k^TV_{k-1}^T\cdots V_{k-m+1}^T)H_0(V_{k-m+1}\cdots V_{k-1}V_k) \\ & + (V_k^TV_{k-1}^T\cdots V_{k-m+2}^T)\rho_{k-m+1} s_{k-m+1} s_{k-m+1}^T(V_{k-m+2}\cdots V_{k-1}V_k) \\ & \cdots \\ & + V_k^T \rho_{k-1} s_{k-1} s_{k-1}^TV_k \\ & + \rho_k s_k s_k^T \tag{LBFGS[2]} \end{aligned}$

（LBFGS 双圈循环算法）
if $\leq m$
$\quad incr = 0, bound = k$
else
$\quad incr = k-m, bound = m$
$\quad$
$q_{bound} = \nabla f_k$
for $\cdots, 0$
$\quad j = i+incr$
$\quad \alpha_i = \rho_j s_j^Tq_{i+1}$
$\quad q_i = q_{i+1} - \alpha_i y_j$
$\quad$
$r_0 = H_0 q_0$
for $\cdots, bound-1$
$\quad j = i + incr$
$\quad \beta_i = \rho_j y_j^T r_i$
$\quad r_{i+1} = r_i + s_j(\alpha_i - \beta_i)$
$\quad$
$H_k \nabla f_k =r$ 是要求的下降方向

下面证明（LBFGS 双圈循环算法）可以正确推导出上面的 $H_{k+1}$ 表达式：
$\begin{aligned} k \leq m: \qquad \qquad \qquad \quad&\\ Loop 1: \quad q_k = & \nabla f_k\\ q_{k-i} = & q_{k-i+1} - \alpha_{k-i}y_{k-i}\\ = & q_{k-i+1} - p_{k-i}s^T_{k-i}q_{k-i+1}y_{k-i}\\ = & (I-p_{k-i}y_{k-i}s^T_{k-i})q_{k-i+1}\\ = & V_{k-i}q_{k-i+1}\\ = & V_{k-i}V_{k-i+1}\cdots V_{k-1}q_k\\ \alpha_{k-i} = & \rho_{k-i}s_{k-i}^Tq_{k-i+1} \\ = & \rho_{k-i}s_{k-i}^TV_{k-i+1}\cdots V_{k-1}q_k\\ \quad\\ Termination:\quad q_0 = & V_0V_1\cdots V_{k-1}\nabla f_k\\ \alpha_0 = & \rho_0s_0^TV_1\cdots V_{k-1}\nabla f_k\\ \quad\\ Loop 2: \quad r_0 = & H_0q_0 = H_0V_0V_1\cdots V_{k-1}\nabla f_k\\ r_{i+1} = & r_i + s_i(\alpha_i - \beta_i)& \\ = & r_i + s_i\alpha_i - \rho_is_iy_i^Tr_i\\ = & (I-\rho_is_iy_i^T) r_i + s_i\alpha_i\\ = & V_i^Tr_i + s_i\alpha_i\\ r_k = & V_{k-1}^Tr_{k-1}+ s_{k-1}\alpha_{k-1}\\ = &V_{k-1}^T(V_{k-2}^Tr_{k-2} + s_{k-1}\alpha_{k-1}) + \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k\\ = & V_{k-1}^TV_{k-2}^Tr_{k-2} + V_{k-1}^T\rho_{k-2}s_{k-2}s_{k-2}^TV_{k-1}\nabla f_k+ \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k\\ \cdots\\ = & V_{k-1}^TV_{k-2}^T\cdots V_1^Tr_0 + V_{k-1}^TV_{k-2}^T\cdots V_1^T\rho_0s_0s_0^TV_1\cdots V_{k-1}\nabla f_k\\ & + \cdots + V_{k-1}^T\rho_{k-2}s_{k-2}s_{k-2}^TV_{k-1}\nabla f_k + \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k\\ = & V_{k-1}^TV_{k-2}^T\cdots V_1^TH_0V_0V_1\cdots V_{k-1}\nabla f_k\\ & + V_{k-1}^TV_{k-2}^T\cdots V_1^T\rho_0s_0s_0^TV_1\cdots V_{k-1}\nabla f_k\\ & + \cdots \\ & + V_{k-1}^T\rho_{k-2}s_{k-2}s_{k-2}^TV_{k-1}\nabla f_k\\ & + \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k \end{aligned}$
刚好等于式(LBFGS[1])计算出的 $H_k \cdot \nabla f_k$ ，也就是LBFGS的下降方向

$\begin{aligned} k > m: \qquad \qquad \qquad \quad&\\ Loop 1: \quad q_m = & \nabla f_k\\ q_{m-i} = & q_{m-i+1} - \alpha_{m-i}y_{k-i}\\ = & q_{m-i+1} - p_{k-i}s^T_{k-i}q_{m-i+1}y_{k-i}\\ = & (I-p_{k-i}y_{k-i}s^T_{k-i})q_{m-i+1}\\ = & V_{k-i}q_{m-i+1}\\ = & V_{k-i}V_{k-i+1}\cdots V_{k-1}q_m\\ \alpha_{m-i} = & \rho_{k-i}s_{k-i}^Tq_{m-i+1} \\ = & \rho_{k-i}s_{k-i}^TV_{k-i+1}\cdots V_{k-1}q_m\\ \quad\\ Termination:\quad q_0 = & V_{k-m}V_{k-m+1}\cdots V_{k-1}\nabla f_k\\ \alpha_0 = & \rho_{k-m}s_{k-m}^TV_{k-m+1}\cdots V_{k-1}\nabla f_k\\ \quad\\ Loop 2: \quad r_0 = & H_0q_0 = H_0V_{k-m}V_{k-m+1}\cdots V_{k-1}\nabla f_k\\ r_{i+1} = & r_i + s_{k-m+i}(\alpha_i - \beta_i)& \\ = & r_i + s_{k-m+i}\alpha_i - \rho_{k-m+i}s_{k-m+i}y_{k-m+i}^Tr_i\\ = & (I-\rho_{k-m+i}s_{k-m+i}y_{k-m+i}^T) r_i + s_{k-m+i}\alpha_i\\ = & V_{k-m+i}^Tr_i + s_{k-m+i}\alpha_i\\ r_m = & V_{k-1}^Tr_{m-1}+ s_{k-1}\alpha_{m-1}\\ = &V_{k-1}^T(V_{k-2}^Tr_{m-2} + s_{k-1}\alpha_{m-1}) + \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k\\ = & V_{k-1}^TV_{k-2}^Tr_{m-2} + V_{k-1}^T\rho_{k-2}s_{k-2}s_{k-2}^TV_{k-1}\nabla f_k+ \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k\\ \cdots\\ = & V_{k-1}^TV_{k-2}^T\cdots V_{k-m}^Tr_0 + V_{k-1}^TV_{k-2}^T\cdots V_{k-m+1}^T\rho_{k-m}s_{k-m}s_{k-m}^TV_{k-m+1}\cdots V_{k-1}\nabla f_k\\ & + \cdots + V_{k-1}^T\rho_{k-2}s_{k-2}s_{k-2}^TV_{k-1}\nabla f_k + \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k\\ = & V_{k-1}^TV_{k-2}^T\cdots V_{k-m}^TH_0V_{k-m}V_{k-m+1}\cdots V_{k-1}\nabla f_k\\ & + V_{k-1}^TV_{k-2}^T\cdots V_{k-m+1}^T\rho_{k-m}s_{k-m}s_{k-m}^TV_{k-m+1}\cdots V_{k-1}\nabla f_k\\ & + \cdots \\ & + V_{k-1}^T\rho_{k-2}s_{k-2}s_{k-2}^TV_{k-1}\nabla f_k\\ & + \rho_{k-1}s_{k-1}s_{k-1}^T\nabla f_k \end{aligned}$
刚好等于式(LBFGS[2])计算出的 $H_k \cdot \nabla f_k$ ，也就只保存m个 $s_i,y_i$ 时计算的LBFGS下降方向，得证。