牛顿法与拟牛顿法

发表于 | 分类于 |

在现在的研究和学习中经常会碰到非线性最优化（nonlinear optimization）的问题，这种问题在机器学习以及系统辨识领域比较常见，之前对这一类算法的原理没有学习过，这几天就在网上查找资料，并且整理了一下，争取每一步的推导都尽可能详细。这篇博客中，从前到后分别介绍了牛顿法，拟牛顿条件，秩一校正，DFP算法，BFGS算法以及其改进的L-BFGS算法，其中特别补充了BFGS算法应用Sherman-Morrison-Woodbury公式改进迭代公式的推导。

牛顿法

牛顿法

牛顿法与一阶的简单的梯度下降法的思想一脉相承。先考虑如下的无约束最小化问题，
$$
\min_x f(x)
$$
其中$x=[x_1, x_2, \dots, x_N]^T\in\mathbb{R}^N$，继续我们假设f(x)为凸函数，且二阶连续可微。简单起见，首先考虑$N=1$的情形。牛顿法的基本思想是：在现有的极小点估计值的附近对f(x)做二阶泰勒展开，进而找到极小点的下一个估计值。设$x_k$为当前的极小点估计值，则有：
$$
\varphi(x)=f(x_{k})+f^{\prime}(x_{k})(x-x_{k})+\frac{1}{2} f^{\prime \prime}(x_{k})(x-x_{k})^{2}
$$
表示$f(x)$在极小点附近的二阶泰勒展开式，由于其求出来的是最值，由极值必要条件可知，$\varphi(x)$因满足条件
$$
\varphi^{\prime}(x)=0
$$
即
$$
f^{\prime}(x_{k})+f^{\prime \prime}(x_{k})(x-x_{k})=0
$$
于是，若给定初始值$x_0$，则可以构造如下的迭代公式
$$
x_{k+1}=x_{k}-\frac{f^{\prime}(x_{k})}{f^{\prime \prime}(x_{k})}, \quad k=0,1, \cdots
$$
使用这个迭代公式来产生序列${x_k}$来逼近$f(x)$的极小点。在一定的条件下其可以收敛到极小点。

对于N>1的情形，可以对二阶泰勒展开式做推广，此时
$$
\varphi(\mathbf{x})=f(\mathbf{x}_{k})+\nabla f(\mathbf{x}_{k}) \cdot(\mathbf{x}-\mathbf{x}_{k})+\frac{1}{2} \cdot(\mathbf{x}-\mathbf{x}_{k})^{T} \cdot \nabla^{2} f(\mathbf{x}_{k}) \cdot(\mathbf{x}-\mathbf{x}_{k})
$$
其中$\nabla f$为f的梯度向量，$\nabla^{2} f$为f的海森矩阵（Hessian），其定义分别为：
$$
\nabla f=
\begin{bmatrix}
\frac{\partial f}{\partial x_{1}} \ \frac{\partial f}{\partial x_{2}} \ {\vdots} \ \frac{\partial f}{\partial x_{N}}
\end{bmatrix},
\nabla^{2} f=
\begin{bmatrix}
\frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{N}} \\
\frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{2} \partial x_{N}}\\
& & \ddots & \\
\frac{\partial^{2} f}{\partial x_{N} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{N} \partial x_{2}} & \ldots & \frac{\partial^{2} f}{\partial x_{N}^{2}}
\end{bmatrix}_{N\times N}
$$
需要注意的是，上式中$\nabla f$和$\nabla^{2} f$中的元素均为关于$\mathbf{x}$的函数，以下分别将其记为g和H。两位，若f的混合偏导数可交换顺序，则Hessian矩阵H为对称矩阵。而$\nabla f\mathbf({x}_k)$和$\nabla^2f({x}_k)$则表示将$\mathbf{x}$取为$\mathbf{x}_k$之后得到的实值向量和矩阵，一下分别将其记为$g_k$和$H_k$（这里字母g表示gradient，H表示Hessian）

与N=1时相同，由于是求极小点，极值必要条件要求他为$\varphi(x)$的驻点，即
$$
\nabla \varphi(\mathbf{x})=0
$$
两边同时作用一个梯度算子
$$
\mathrm{g}_{k}+H_{k} \cdot(\mathbf{x}-\mathbf{x}_{k})=0
$$
若Hessian为非奇异矩阵，则
$$
\mathrm{x}=\mathbf{x}_{k}-H_{k}^{-1} \cdot \mathbf{g}_{k}
$$
于是，若给定初始值$x_0$，则同样可以构造出迭代公式
$$
\mathbf{x}_{k+1}=\mathbf{x}_{k}-H_{k}^{-1} \cdot \mathbf{g}_{k}, \quad k=0,1, \cdots
$$
这种方法就是原始牛顿法，其迭代中的搜索方向为$\mathbf{d}_{k}=-H_{k}^{-1} \cdot \mathbf{g}_{k}$称为牛顿方向。

Algorithm Newton Method
1.给定初值$x_0$和精度阈值$\epsilon$，并且令$k=0$
2.计算$\mathbf{g}_k$和$H_k$
3.若$\mathrm{norm}\mathbf{g}_k <\epsilon$，则停止迭代，否则确定搜索方向，则停止迭代，否则确定搜索方向$\mathbf{d}_k = -H_k^{-1} \cdot \mathbf{g}_k​$
4.计算新的迭代点，$x_{k+1}:=x_k+d_k$
5.令k:= k+1, go to step 2

当目标为二次函数时，则其二次泰勒展开式与原式完全相同则Hessian退化为一个常数矩阵，从任意一点出发，通过一步迭代即可以达到极小点，因此牛顿法是一种具有二次收敛性的算法。对于高次函数，牛顿法的收敛速度也比较快。

变步长牛顿法

在每一次更新的时候（牛顿法step 4），采用以下迭代公式
$$
\mathbf{x}_{k+1} :=\mathbf{x}_{k}+\lambda_{k} \mathbf{d}_{k}
$$
其中步长因子为$\lambda_k$,
$$
\lambda_{k}=\arg \min _{\lambda \in \mathbb{R}} f(\mathbf{x}_{k}+\lambda \mathbf{d}_{k})
$$

Algorithm Newton Method 2
1.给定初值$x_0$和精度阈值$\epsilon$，并且令$k=0$
2.计算$\mathbf{g}_k$和$H_k$
3.若$\mathrm{norm} \mathbf{g}_k<\epsilon$，则停止迭代，否则确定搜索方向，则停止迭代，否则确定搜索方向，则停止迭代，否则确定搜索方向，则停止迭代，否则确定搜索方向$\mathbf{d}_k = -H_k^{-1} \cdot \mathbf{g}_k$
4.计算新的迭代点，$\mathbf{x}_{k+1} :=\mathbf{x}_{k}+\lambda_{k} \mathbf{d}_{k}$
5.令k:= k+1, go to step 2

小结

牛顿法有以下缺点：

  • 对目标函数有比较严格的要求，函数必须具有连续的一，二阶偏导数，Hessian矩阵必须正定。
  • 计算复杂，除了计算梯度以外还需要计算二阶偏导数矩阵和它的你矩阵，计算量和存储量都很大。

秩一矫正

在之后要介绍的拟牛顿法中会用到Sherman-morrison定理，所以在这里先简单介绍一下秩一校正（rank-one update）。在原矩阵$A\in \mathbb{R}^{n\times n}$上增加一个秩为1的新矩阵，通常为$uv^T, u,v\in \mathbb{R}^{n\times 1}$（u与v为列向量）。我们研究增加秩一矩阵（这个过程就是秩一校正）的性质，包括逆矩阵，行列式，特征值的变化情况。

对于为什么$uv^T$为秩为1的矩阵：矩阵的每一次左乘矩阵都相当于增加了一个从原空间到某子空间的投影，因此秩越乘越小（原空间等于像空间时矩阵的秩不变）。矩阵的空间映射只可能从高维空间映射到低位空间，不可能映射成为更高维的空间。

由于秩一校正是通过加法将原矩阵与秩一矩阵连接在一起的，并且逆矩阵运算是这些主题的基础，因此，我们首先研究和式逆矩阵如何用原矩阵进行表示。考虑$n \times n$阶矩阵A与B。在一般情况下，A+B的逆矩阵并不存在有用的公式。但是如果A是可逆矩阵，同时B具有某种特殊的形态，则$(A+B)^{-1}$确实有简单的运算公式。最简单的例子是$I+\mathbf{u v}^{T}$，其中u和v是n维非零向量（即$n \times 1$矩阵），$\mathbf{u} \mathbf{v}^{T}$被称为秩一矩阵。这里导出$I+\mathbf{u} \mathbf{v}^{T}$的逆矩阵，再利用此结果推演$(A+\mathbf{u} \mathbf{v}^{T})^{-1}$的计算公式。

考虑线性方程
$$
(I+\mathbf{u v}^{T}) \mathbf{x}=\mathbf{b}
$$
解出x即可以得到$(I+\mathbf{u} \mathbf{v}^{T})^{-1}$的公式。令$k=\mathbf{v}^{T} \mathbf{x}$，乘开上式可得$\mathbf{x}+k \mathbf{u}=\mathbf{b}$。等号两边左乘$\mathbf{v}^T$，就有$k+k(\mathbf{v}^{T} \mathbf{u})=\mathbf{v}^{T} \mathbf{b}$，即得$k=(1+\mathbf{v}^{T} \mathbf{u})^{-1} \mathbf{v}^{T} \mathbf{b}$。合并以上结果，
$$
\mathbf{x}=\mathbf{b}-k \mathbf{u}=\mathbf{b}-\frac{\mathbf{v}^{T} \mathbf{b}}{1+\mathbf{v}^{T} \mathbf{u}} \mathbf{u}=\mathbf{b}-\frac{\mathbf{uv}^{T}}{1+\mathbf{v}^{T} \mathbf{u}} \mathbf{b}=(I-\frac{\mathbf{uv}^{T}}{1+\mathbf{v}^{T} \mathbf{u}}) \mathbf{b}
$$
上式中，$(\mathbf{v}^{T} \mathbf{b}) \mathbf{u}=\mathbf{u}(\mathbf{v}^{T} \mathbf{b})=(\mathbf{u} \mathbf{v}^{T}) \mathbf{b}$因为$\mathbf{v}^T\mathbf{b}$是一个纯量，并使用矩阵乘法结合律。若$1+\mathbf{v}^{T} \mathbf{u} \neq 0$，则$I+\mathbf{u v}^{T}$可逆，且
$$
(I+\mathbf{u v}^{T})^{-1}=I-\frac{\mathbf{u} \mathbf{v}^{T}}{1+\mathbf{v}^{T} \mathbf{u}}
$$
令A为一个$n\times n$阶可逆矩阵，以A取代I，利用$I+\mathbf{u} \mathbf{v}^{T}$的逆矩阵公式可导出$A+\mathbf{u} \mathbf{v}^{T}$的逆矩阵，过程如下
$$
\begin{align}
\begin{split}
(A+\mathbf{u v}^{T})^{-1} &=(A(I+A^{-1} \mathbf{u v}^{T}))^{-1}=(I+A^{-1} \mathbf{u v}^{T})^{-1} A^{-1} \\
&=(I-\frac{A^{-1} \mathbf{u v}^{T}}{1+\mathbf{v}^{T} A^{-1} \mathbf{u}}) A^{-1}=A^{-1}-\frac{A^{-1} \mathbf{u} \mathbf{v}^{T} A^{-1}}{1+\mathbf{v}^{T} A^{-1} \mathbf{u}}
\end{split}
\end{align}
$$
证明：
$$
\begin{align}
\begin{split}
&(A+u v^{T})(A^{-1}-\frac{A^{-1} u v^{T} A^{-1}}{1+v^{T} A^{-1} u})\\
=&A A^{-1}+u v^{T} A^{-1}-\frac{u v^{T} A^{-1}+u v^{T} A^{-1} u v^{T} A^{-1}}{1+v^{T} A^{-1} u} \\
=&I+u v^{T} A^{-1}-\frac{u(1+v^{T} A^{-1} u) v^{T} A^{-1}}{1+v^{T} A^{-1} u} \\
=&I+u v^{T} A^{-1}-u v^{T} A^{-1}\\
=&I
\end{split}
\end{align}
$$
上式称为Sherman-Morrison公式，$(A+\mathbf{u v}^{T})^{-1}$存在的条件是$1+\mathbf{v}^{T} A^{-1} \mathbf{u} \neq 0$。若U和V是$n \times m$阶矩阵使得$I_{m}+V^{T} A^{-1} U$可逆，则Sherman-Morrison公式可以进一步推广为
$$
(A+U V^{T})^{-1}=A^{-1}-A^{-1} U(I_{m}+V^{T} A^{-1} U)^{-1} V^{T} A^{-1}
$$
称为Sherman-Morrison-Woodbury公式。

拟牛顿法

为了克服以上提出的牛顿法的两个主要缺点，人们提出了拟牛顿法。拟牛顿法的主要思想是：不用二阶偏导数而构造出可以近似Hessian矩阵（或Hessian矩阵的逆矩阵）的正定对称阵，在“拟牛顿”的条件下优化目标函数。而不同的构造方法也就产生了不同的拟牛顿法。

简单来说，拟牛顿法就是利用不同的手段来构造Hessian矩阵的近似，而规避了计算和储存Hessian矩阵的困难。

下文中，Hessian矩阵用H来表示，Hessian矩阵的近似阵用B来表示，Hessian矩阵的逆矩阵的近似矩阵用D来表示。即$B \approx H, D \approx H^{-1}​$。

拟牛顿条件

再经过k+1次的迭代之后得到$x_{k+1}$，此时将目标函数$f(x)$在在$x_{k+1}​$附近做泰勒展开，取二阶近似值得到
$$
f(\mathbf{x}) \approx f(\mathbf{x}_{k+1})+\nabla f(\mathbf{x}_{k+1}) \cdot(\mathbf{x}-\mathbf{x}_{k+1})+\frac{1}{2} \cdot(\mathbf{x}-\mathbf{x}_{k+1})^{T} \cdot \nabla^{2} f(\mathbf{x}_{k+1}) \cdot(\mathbf{x}-\mathbf{x}_{k+1})
$$
两边同时作用一个梯度算子$\nabla$有
$$
\nabla f(\mathbf{x}) \approx \nabla f(\mathbf{x}_{k+1})+H_{k+1} \cdot(\mathbf{x}-\mathbf{x}_{k+1})
$$
在上式中取$x=x_k$，经过整理有
$$
\mathbf{g}_{k+1}-\mathbf{g}_{k} \approx H_{k+1} \cdot(\mathbf{x}_{k+1}-\mathbf{x}_{k})
$$
引入记号
$$
\mathbf{s}_{k}=\mathbf{x}_{k+1}-\mathbf{x}_{k}, \quad \mathbf{y}_{k}=\mathbf{g}_{k+1}-\mathbf{g}_{k}
$$
则
$$
\mathbf{y}_{k} \approx H_{k+1} \cdot \mathbf{s}_{k}
$$
或者
$$
\mathbf{s}_{k} \approx H_{k+1}^{-1} \cdot \mathbf{y}_{k}
$$
以上的两种形式为拟牛顿条件，而不同的拟牛顿法将迭代过程中的Hessian矩阵或者其逆矩阵近似为$B_{k+1}$或者$D_{k+1}$，有
$$
\mathbf{y}_{k}=B_{k+1} \cdot \mathbf{s}_{k}
$$
或者
$$
\mathbf{s}_{k}=D_{k+1} \cdot \mathbf{y}_{k}
$$

DFP算法

DFP法是最早的拟牛顿法，由Davidon与1959年提出。该算法的核心是通过迭代的方法，对$H_{k+1}^{-1}$做近似。其迭代格式为：
$$
D_{k+1}=D_{k}+\Delta D_{k}, \quad k=0,1,2, \cdots
$$
其中的$D_0$通常取单位矩阵I。因此，关键是每一步的校正矩阵$\Delta D_k$如何构造。

注意，迭代格式将嵌套在之前的牛顿法的算法当中。因为我们可以猜想$\Delta D_k$可能与$s_k$, $y_k$和$D_k$有关系。我们可以在这里采用“待定法”。即首先将$\Delta D_k$待定成某种特殊的形式，然后结合拟牛顿条件来进行推导。

那么$\Delta D_k$需要待定成何种格式呢？这里用到了一个比较tricky的方法，将其待定为
$$
\Delta D_{k}=\alpha \mathbf{u u}^{T}+\beta \mathbf{v} \mathbf{v}^{T}
$$
其中，$\alpha, \beta$为待定系数，$\mathbf{u}, \mathbf{v} \in \mathbb{R}^{N}$为待定向量。从形式上看，这种待定公式是少保证了矩阵$\Delta D_k$的对称性（因为向量乘以自身的转置向量均为对称矩阵）。

将上式带入之前的迭代式，结合在拟牛顿法中设定的s与y的关系式$\mathbf{s}_{k}=D_{k+1}\cdot\mathbf{y}_{k}$，有
$$
\mathbf{s}_{k}=D_{k} \mathbf{y}_{k}+\alpha \mathbf{u u}^{T} \mathbf{y}_{k}+\beta \mathbf{v} \mathbf{v}^{T} \mathbf{y}_{k}
$$
此式中将“数”单独提出来，有
$$
\begin{align}
\begin{split}
\mathbf{s}_{k} &=D_{k} \mathbf{y}_{k}+\mathbf{u}(\alpha \mathbf{u}^{T} \mathbf{y}_{k})+\mathbf{v}(\beta \mathbf{v}^{T} \mathbf{y}_{k}) \\
&=D_{k} \mathbf{y}_{k}+(\alpha \mathbf{u}^{T} \mathbf{y}_{k}) \mathbf{u}+(\beta \mathbf{v}^{T} \mathbf{y}_{k}) \mathbf{v}
\end{split}
\end{align}
$$
这里，括号中的两项均为数，我们可以做如下简单赋值
$$
\alpha \mathbf{u}^{T} \mathbf{y}_{k}=1, \quad \beta \mathbf{v}^{T} \mathbf{y}_{k}=-1
$$
即为
$$
\alpha=\frac{1}{\mathbf{u}^{T} \mathbf{y}_{k}}, \quad \beta=-\frac{1}{\mathbf{v}^{T} \mathbf{y}_{k}}
$$
其中两个向量$\mathbf{u}, \mathbf{v}$仍有待确定。

为了确定这两个待定向量，我们将刚刚的赋值式带回到原式，则有
$$
\mathbf{u}-\mathbf{v}=\mathbf{s}_{k}-D_{k} \mathbf{y}_{k}
$$
欲使上式成立，不妨直接取
$$
\mathbf{u}=\mathbf{s}_{k}, \quad \mathbf{v}=D_{k} \mathbf{y}_{k}
$$
在将其带入上面的设定式可以得到
$$
\alpha=\frac{1}{\mathrm{s}_{k}^{T} \mathrm{y}_{k}}, \quad \beta=-\frac{1}{(D_{k} \mathrm{y}_{k})^{T} \mathrm{y}_{k}}=-\frac{1}{\mathrm{y}_{k}^{T} D_{k} \mathrm{y}_{k}}
$$
其中第二个等式用到了$D_k$的对称性。其实为Hessian的对称性，如果$f$函数在区域$D$内的每个二阶导数都是连续函数，那么f的Hessian矩阵在$D$区域内为对称矩阵。

至此我们已经将校正矩阵$\Delta D_k$构造出来了，将以上设定带回到$\Delta D_k$的待定式，可以得到
$$
\Delta D_{k}=\frac{\mathbf{s}_{k} \mathbf{s}_{k}^{T}}{\mathbf{s}_{k}^{T} \mathbf{y}_{k}}-\frac{D_{k} \mathbf{y}_{k} \mathbf{y}_{k}^{T} D_{k}}{\mathbf{y}_{k}^{T} D_{k} \mathbf{y}_{k}}
$$
综上我们可以给出DFP算法的一个完整描述。

Algorithm DFP
1.给定初值$x_0$和精度阈值$\epsilon$，并且令$D_0=I$，$k:=0$
2.确定搜索方向$\mathrm{d}_{k}=-D_{k} \cdot \mathrm{g}_{k}$
3.利用$\lambda_{k}=\arg \min _{\lambda \in \mathbb{R}} f\left(\mathrm{x}_{k}+\lambda \mathrm{d}_{k}\right)$得到步长$\lambda_k$，令$\mathbf{s}_{k}=\lambda_{k} \mathbf{d}_{k}， \mathbf{x}_{k+1} :=\mathbf{x}_{k}+\mathbf{s}_{k}$
4.若$\mathrm{norm} g_{k+1} < \epsilon$，则算法结束。
5.计算$\mathbf{y}_{k}=\mathbf{g}_{k+1}-\mathbf{g}_{k}$
6.计算$D_{k+1}=D_{k}+\frac{\mathbf{s}_{k} \mathbf{s}_{k}^{T}}{\mathbf{s}_{k}^{T} \mathbf{y}_{k}}-\frac{D_{k} \mathbf{y}_{k} \mathbf{y}_{k}^{T} D_{k}}{\mathbf{y}_{k}^{T} D_{k} \mathbf{y}_{k}}$
7.令$k :=k+1$，goto Step 2

BFGS算法

BFGS算法是以其发明者Broyden, Fletcher, Goldfarb 和 Shanno 四个人的名字的首字母命名的。与DFP算法相比，BFGS算法的性能更佳。目前它已成为求解无约束非线性优化问题最常用的方法之一。BFGS算法已有较完善的局部收敛理论，对其全局收敛性的研究也取得了重要成果。

BFGS算法中核心公式的推导过程与DFP完全类似，知识互换了其中$s_k$和$y_k$的位置。需要注意的是BFGS算法是直接逼近海森矩阵，即$B_{k} \approx H_{k}$。仍然采用迭代方法，设迭代格式为
$$
B_{k+1}=B_{k}+\Delta B_{k}, \quad k=0,1,2, \cdots
$$
其中的$B_0$也经常取为单位矩阵I。因为，关键是每一步的校正矩阵$\Delta B_k$如何构造。与DFP法类似的，将其待定为
$$
\Delta B_{k}=\alpha \mathbf{u u}^{T}+\beta \mathbf{v} \mathbf{v}^{T}
$$
将其带入上式并类似于DFP算法，可将其整理为
$$
\mathbf{y}_{k}=B_{k} \mathbf{s}_{k}+(\alpha \mathbf{u}^{T} \mathbf{s}_{k}) \mathbf{u}+(\beta \mathbf{v}^{T} \mathbf{s}_{k}) \mathbf{v}
$$
通过令
$$
\alpha \mathbf{u}^{T} \mathbf{s}_{k}=1, \beta \mathbf{v}^{T} \mathbf{s}_{k}=-1
$$

$$
\mathbf{u}=\mathbf{y}_{k}, \quad \mathbf{v}=B_{k} \mathbf{s}_{k}
$$

可算得
$$
\alpha=\frac{1}{\mathrm{y}_{k}^{T} \mathrm{s}_{k}}, \quad \beta=-\frac{1}{\mathrm{s}_{k}^{T} B_{k} \mathrm{s}_{k}}
$$
综上我们便得到了校正矩阵$\Delta B_k$的构造公式
$$
\Delta B_{k}=\frac{\mathrm{y}_{k} \mathrm{y}_{k}^{T}}{\mathrm{y}_{k}^{T} \mathrm{s}_{k}}-\frac{B_{k} \mathrm{s}_{k} \mathrm{s}_{k}^{T} B_{k}}{\mathrm{s}_{k}^{T} B_{k} \mathrm{s}_{k}}
$$
现在我们可以将DFP算法与BFGS算法的两种校正矩阵$\Delta D_k$与$\Delta B_k$比较一下，可以发现，两个构造式只是将B与D相互替换，并将$s_k$与$y_k$互换位置（因为DFP算法是对Hessian矩阵的逆矩阵进行近似，而BFGS算法则是直接对Hessian矩阵进行近似计算）而已。

最后我们给出BFGS的完整算法描述。

Algorithm BFGS Method 1
1.给定初值$x_0$和精度阈值$\epsilon$，并且令$B_0=I，k:=0$
2.确定搜索方向$\mathbf{d}_{k} = -B_k^{-1}\cdot\mathbf{g}_k$
3.利用$\lambda_{k}=\arg \min _{\lambda \in \mathbb{R}} f\left(\mathrm{x}_{k}+\lambda \mathrm{d}_{k}\right)$得到步长得到步长$\lambda_k$，令，令$\mathbf{s}_{k}=\lambda_{k} \mathbf{d}_{k}， \mathbf{x}_{k+1} :=\mathbf{x}_{k}+\mathbf{s}_{k}$
4.若$\mathrm{norm} g_{k+1} < \epsilon$，则算法结束。
5.计算$\mathbf{y}_{k}=\mathbf{g}_{k+1}-\mathbf{g}_{k}$
6.计算$B_{k+1}=B_{k}+\frac{\mathbf{y}_{k} \mathbf{y}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}-\frac{B_{k} \mathbf{s}_{k} \mathbf{s}_{k}^{T} B_{k}}{\mathbf{s}_{k}^{T} B_{k} \mathbf{s}_{k}}$
7.令k :=k+1，goto Step 2

以上算法中的Step2通常是通过求解线性代数方程组$B_k\mathbf{d}_k = -\mathbf{g}_k$来进行的。如果真的进行这一步的话则需要求解矩阵$B_k^{-1}$，当求解参数较多时，求解逆矩阵会耗费很多算力，所以通过Sherman-Morrison公式来改善算法的性能。

更一般的做法是通过对Step6中的递推关系应用Sherman-Morrison公式，直接给出$B_{k+1}^{-1}$与$B_{k}^{-1}$之间的关系式
$$
B_{k+1}^{-1}=\left(I-\frac{\mathbf{s}_{k} \mathbf{y}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}\right) B_{k}^{-1}\left(I-\frac{\mathbf{y}_{k} \mathbf{s}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}\right)+\frac{\mathbf{s}_{k} \mathbf{s}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}
$$
或者进一步展开写成
$$
B_{k+1}^{-1}=B_{k}^{-1}+\left(\frac{1}{\mathbf{s}_{k}^{T} \mathbf{y}_{k}}+\frac{\mathbf{y}_{k}^{T} B_{k}^{-1} \mathbf{y}_{k}}{\left(\mathbf{s}_{k}^{T} \mathbf{y}_{k}\right)^{2}}\right) \mathbf{s}_{k} \mathbf{s}_{k}^{T}-\frac{1}{\mathbf{s}_{k}^{T} \mathbf{y}_{k}}\left(B_{k}^{-1} \mathbf{y}_{k} \mathbf{s}_{k}^{T}+\mathbf{s}_{k} \mathbf{y}_{k}^{T} B_{k}^{-1}\right)
$$
下面让我们来证明一下，这里我们使用Sherman-Morrison-Woodbury公式。这里先对公式进行以下定义， 有
$$
(\boldsymbol{A}+\boldsymbol{U} \boldsymbol{V}^{T})^{-1}=\boldsymbol{A}^{-1}-\boldsymbol{A}^{-1} \boldsymbol{U}(\boldsymbol{I}+\boldsymbol{V}^{T} \boldsymbol{U})^{-1} \boldsymbol{V}^{T} \boldsymbol{A}^{-1}
$$
其中
$$
U=\left[u_{1}, u_{2}, \ldots, u_{k}\right] \quad V=\left[v_{1}, v_{2}, \ldots, v_{k}\right]
$$
下面我们将这个公式展开一下
$$
\left(A+\sum_{i=1}^{k} u_{i} v_{i}^{T}\right)^{-1}=A^{-1}-A^{-1} U C^{-1} V^{T} A^{-1}
$$
其中
$$
C_{i j}=\delta_{i j}+\boldsymbol{v}_{i}^{T} \boldsymbol{u}_{j} \quad i, j=1,2, \ldots, k
$$
基于Sherman-Morrison-Woodbury的展开式，我们考虑$k=2$的情况，并且设
$$
u_{1}=v_{1}=\frac{y_{k}}{\left(s_{k}^{T} y_{k}\right)^{1 / 2}} \quad u_{2}=-v_{2}=\frac{B_{k} s_{k}}{\left(s_{k}^{T} B_{k} s_{k}\right)^{1 / 2}}
$$
为了表示上的方便我们使用$\mathbf{D}_k$来替代原本公式中的$\mathbf{B}_k^{-1}​$，则
$$
\begin{align}
\begin{split}
C_{11}&=1+\boldsymbol{v}_{1}^{T} \boldsymbol{u}_{1}=1+\frac{\boldsymbol{y}_{k}^{T} \boldsymbol{D}_{k} \boldsymbol{y}_{k}}{\boldsymbol{s}_{k}^{T} \boldsymbol{y}_{k}} \\
C_{22}&=1+\boldsymbol{v}_{2}^{T} \boldsymbol{u}_{2}=-\frac{\boldsymbol{s}_{k}^{T} \boldsymbol{B}_{k} \boldsymbol{D}_{k} \boldsymbol{B}_{k} \boldsymbol{s}_{k}}{\boldsymbol{s}_{k}^{T} \boldsymbol{B}_{k} \boldsymbol{s}_{k}}=1-1=0\\
C_{12}&=\boldsymbol{v}_{1}^{T} \boldsymbol{u}_{2}=\frac{\boldsymbol{y}_{k}^{T} \boldsymbol{B}_{k} \boldsymbol{s}_{k}}{\left(s_{k}^{T} \boldsymbol{y}_{k}\right)^{1 / 2}\left(\boldsymbol{s}_{k}^{T} \boldsymbol{B}_{k} \boldsymbol{s}_{k}\right)^{1 / 2}}=\frac{\left(\boldsymbol{s}_{k}^{T} \boldsymbol{B}_{k} \boldsymbol{s}_{k}\right)^{1 / 2}}{\left(\boldsymbol{s}_{k}^{T} \boldsymbol{y}_{k}\right)^{1 / 2}}\\
C_{21}&=\boldsymbol{v}_{2}^{T} \boldsymbol{u}_{1}=-C_{12}
\end{split}
\end{align}
$$
这时，我们可以将矩阵$C$表示称以下的形式
$$
C=\left( \begin{array}{cc}{\beta} & {\alpha} \ {-\alpha} & {0}\end{array}\right) \quad C^{-1}=\frac{1}{\alpha^{2}} \left( \begin{array}{cc}{0} & {-\alpha} \ {\alpha} & {\beta}\end{array}\right)
$$

$$
\beta=1+\frac{\boldsymbol{y}_{k}^{T} \boldsymbol{D}_{k} \boldsymbol{y}_{k}}{\boldsymbol{s}_{k}^{T} \boldsymbol{y}_{k}} \quad \alpha=\frac{\left(\boldsymbol{s}_{k}^{T} \boldsymbol{B}_{k} \boldsymbol{s}_{k}\right)^{1 / 2}}{\left(\boldsymbol{s}_{k}^{T} \boldsymbol{y}_{k}\right)^{1 / 2}}
$$

我们设定$\tilde{\mathbf{U}}=\mathbf{D}_{k} \mathbf{U}$并且$\tilde{\mathbf{V}}=\mathbf{D}_{k} \mathbf{V}$，这里
$$
\tilde{\mathbf{u}}_{i}=\mathbf{D}_{k} \mathbf{u}_{i}, \quad \tilde{\mathbf{v}}_{i}=\mathbf{D}_{k} \mathbf{v}_{i} \quad i=1,2
$$
这时我们就可以很直观的得到
$$
\begin{align}
\begin{split}
\boldsymbol{D}_{k+1} &\leftarrow \boldsymbol{D}_{k}-\boldsymbol{D}_{k} \boldsymbol{U} \boldsymbol{C}^{-1} \boldsymbol{V}^{T} \boldsymbol{D}_{k}=\boldsymbol{D}_{k}-\tilde{\boldsymbol{U}} \boldsymbol{C}^{-1} \tilde{\boldsymbol{V}}^{T}\\
&=\boldsymbol{D}_{k}+\frac{1}{\alpha}\left(-\widetilde{\boldsymbol{u}}_{1} \widetilde{\boldsymbol{v}}_{2}^{T}+\widetilde{\boldsymbol{u}}_{2} \widetilde{\boldsymbol{v}}_{1}^{T}\right)-\frac{\beta}{\alpha^{2}} \widetilde{\boldsymbol{u}}_{2} \widetilde{\boldsymbol{v}}_{2}^{T}
\end{split}
\end{align}
$$
再将之前设定的$\alpha, \beta, \mathbf{U}, \mathbf{V}​$带回到式中可以得到
$$
\boldsymbol{D}_{k+1} \leftarrow \boldsymbol{D}_{k}-\frac{\boldsymbol{D}_{k} \boldsymbol{y}_{k} \boldsymbol{s}_{k}^{T}+\boldsymbol{s}_{k} \boldsymbol{y}_{k}^{T} \boldsymbol{D}_{k}}{\boldsymbol{s}_{k}^{T} \boldsymbol{y}_{k}}+\frac{\boldsymbol{s}_{k} \boldsymbol{s}_{k}^{T}}{\boldsymbol{s}_{k}^{T} \boldsymbol{y}_{k}}\left(1+\frac{\boldsymbol{y}_{k}^{T} \boldsymbol{D}_{k} \boldsymbol{y}_{k}}{\boldsymbol{s}_{k}^{T} \boldsymbol{y}_{k}}\right)
$$
Q.E.D

利用上式我们可以简单的得到$B_{k+1}^{-1}$与$B_{k}^{-1}$之间的关系式，由此我们可以得到改进后的BFGS算法。

Algorithm BFGS Method 2
1.给定初值$x_0$和精度阈值$\epsilon$，并且令$D_0=I，k:=0$
2.确定搜索方向$\mathbf{d}_{k}=-D_{k} \cdot \mathbf{g}_{k}$
3.利用$\lambda_{k}=\arg \min _{\lambda \in \mathbb{R}} f\left(\mathrm{x}_{k}+\lambda \mathrm{d}_{k}\right)$得到步长得到步长$\lambda_k$，令，令$\mathbf{s}_{k}=\lambda_{k} \mathbf{d}_{k}， \mathbf{x}_{k+1} :=\mathbf{x}_{k}+\mathbf{s}_{k}$
4.若$\mathrm{norm} g_{k+1} < \epsilon$，则算法结束。
5.计算$\mathbf{y}_{k}=\mathbf{g}_{k+1}-\mathbf{g}_{k}$
6.计算$D_{k+1}=\left(I-\frac{\mathbf{s}_{k} \mathbf{y}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}\right) D_{k}\left(I-\frac{\mathbf{y}_{k} \mathbf{s}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}\right)+\frac{\mathbf{s}_{k} \mathbf{s}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}$
7.令k :=k+1，goto Step 2

其实从一开始到现在，我们在这些算法中使用的计算步长$\lambda_k$的方法都是一种精确搜索。在实际应用中还有很多像Wolfe型搜索，Armijo型搜索以及满足Goldstrin条件的非精确搜索。带非精确搜索的拟牛顿法是从1976年Powell的工作开始的，它证明了带Wolfe搜索的BFGS算法的全局收敛性和超线性收敛性。

L-BFGS算法

在这之前介绍的两种拟牛顿法（DFP算法以及BFGS算法）都解决了原本牛顿法中需要计算Hessian矩阵的问题，使其更容易在计算机上实现，但随着仿真（Simulation）技术以及最优化（Optimization）技术以及机器学习（Machine Learning）技术的发展我们所计算的参数越来越多，当参数多到一定程度的时候，那么在这两种方法中用以存储近似Hessian矩阵的内存就不够用了。假设我们有一个10000参数的评价函数需要进行最优化，那么就需要存储一个$10000 \times 10000$的矩阵$D_k$，假设我们使用Double型来存储数据，那么我们需要
$$
\frac{10000 \times 10000 \times 8}{2^{10} \times 2^{10} \times 2^{10}} = 0.745 (GB)
$$
对于一般的个人电脑来说这已经是一个很大的内存空间了，我们可以想象，例如说在我的研究中，经常会有100多万个数据点需要进行类似的运算，你可以想象那将会需要多大的内存。暂且估算一下，长宽都拉伸100倍，容量需要这个例子的10000倍，也就是说只是存储这个矩阵就需要大约7450GB的内存，就算把握裤子当了也买不起这么多内存。基于以上事实我们就需要对BFGS算法进行改造，从而减少其迭代过程中所需的内存开销。

L-BFGS（Limited-memory BFGS or Limited-storage BFGS）算法就是这样一种改进算法。它对BFGS算法及行了近似，其基本思想是：不再存储完整的矩阵$D_k$，而是存储计算过程中的向量序列$\{\mathbf{s}_i\}, \{\mathbf{y}_i\}$，需要矩阵$D_k$时，利用向量序列$\{\mathbf{s}_i\}, \{\mathbf{y}_i\}$的计算来代替。而且也并不是所有的向量序列都需要储存，而是固定存储最新的m个（参数m可以由用户根据自己的计算机的内存进行自行指定）。每次计算$D_k$时，只利用最新的m个$\{\mathbf{s}_i\}, \{\mathbf{y}_i\}$。显然这样一来，我们的存储空间的空间复杂度由原来的$O(N^2)$降到了$O(mN)$。

首先我们先考虑在上一节改进之后的BFGS算法中$D_k$的迭代式
$$
D_{k+1}=\left(I-\frac{\mathbf{s}_{k} \mathbf{y}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}\right) D_{k}\left(I-\frac{\mathbf{y}_{k} \mathbf{s}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}\right)+\frac{\mathbf{s}_{k} \mathbf{s}_{k}^{T}}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}
$$
若我们记$\rho_{k}=\frac{1}{\mathbf{y}_{k}^{T} \mathbf{s}_{k}}, V_{k}=I-\rho_{k} \mathbf{y}_{k} \mathbf{s}_{k}^{T}$，则上面的迭代式可以写成
$$
D_{k+1}=V_{k}^{T} D_{k} V_{k}+\rho_{k} \mathbf{s}_{k} \mathbf{s}_{k}^{T}
$$
如果设给定的初始矩阵为$D_0$（通常为正定的对角矩阵，如$D_0=I$），则利用上式，依次可得
$$
\begin{align}
\begin{split}
D_{1} &=V_{0}^{T} D_{0} V_{0}+\rho_{0} \mathbf{s}_{0} \mathbf{s}_{0}^{T} \\
D_{2} &=V_{1}^{T} D_{1} V_{1}+\rho_{1} \mathrm{s}_{1} \mathrm{s}_{1}^{T} \\
&=V_{1}^{T}\left(V_{0}^{T} D_{0} V_{0}+\rho_{0} \mathrm{s}_{0} \mathrm{s}_{0}^{T}\right) V_{1}+\rho_{1} \mathrm{s}_{1} \mathrm{s}_{1}^{T}\\
&=V_{1}^{T} V_{0}^{T} D_{0} V_{0} V_{1}+V_{1}^{T} \rho_{0} \mathrm{s}_{0} \mathrm{s}_{0}^{T} V_{1}+\rho_{1} \mathrm{s}_{1} \mathrm{s}_{1}^{T} \\
D_{3} &=V_{2}^{T} D_{2} V_{2}+\rho_{2} \mathrm{s}_{2} \mathrm{s}_{2}^{T}\\
&=V_{2}^{T}\left(V_{1}^{T} V_{0}^{T} D_{0} V_{0} V_{1}+V_{1}^{T} \rho_{0} \mathrm{s}_{0} \mathrm{s}_{0}^{T} V_{1}+\rho_{1} \mathrm{s}_{1} \mathrm{s}_{1}^{T}\right) V_{2}+\rho_{2} \mathrm{s}_{2} \mathrm{s}_{2}^{T}\\
&=V_{2}^{T} V_{1}^{T} V_{0}^{T} D_{0} V_{0} V_{1} V_{2}+V_{2}^{T} V_{1}^{T} \rho_{0} \mathrm{s}_{0} \mathrm{s}_{0}^{T} V_{1} V_{2}+V_{2}^{T} \rho_{1} \mathrm{s}_{1} \mathrm{s}_{1}^{T} V_{2}++\rho_{2} \mathrm{s}_{2} \mathrm{s}_{2}^{T}\\
\cdots
\end{split}
\end{align}
$$
更一般的
$$
\begin{align}
\begin{split}
D_{k+1}=&(V_{k}^{T} V_{k-1}^{T} \cdots V_{1}^{T} V_{0}^{T}) D_{0}(V_{0} V_{1} \cdots V_{k-1} V_{k})\\
+&(V_{k}^{T} V_{k-1}^{T} \cdots V_{2}^{T} V_{1}^{T})(\rho_{0} \mathbf{s}_{0} \mathbf{s}_{0}^{T})(V_{1} V_{2} \cdots V_{k-1} V_{k})\\
+&(V_{k}^{T} V_{k-1}^{T} \cdots V_{3}^{T} V_{2}^{T})(\rho_{1} \mathbf{s}_{1} \mathbf{s}_{1}^{T})(V_{2} V_{3} \cdots V_{k-1} V_{k})\\
+&\cdots \\
+&(V_{k}^{T} V_{k-1}^{T})(\rho_{k-2} \mathbf{S}_{k-2} \mathbf{s}_{k-2}^{T})(V_{k-1} V_{k})\\
+&V_{k}^{T}(\rho_{k-1} \mathbf{S}_{k-1} \mathbf{s}_{k-1}^{T}) V_{k}\\
+&\rho_{k} \mathbf{S}_{k} \mathbf{s}_{k}^{T}
\end{split}
\end{align}
$$
由上式可知，计算$D_{k+1}$需要用到$\left\{s_{i}, y_{i}\right\}_{i=0}^{k}$，因此，若从$\mathbf{s}_{0}, \mathbf{y}_{0}$开始连续的存储m组的话，只能存储到$\mathbf{s}_{m-1}, \mathbf{y}_{m-1}$，亦即，只能依次计算$D_{1}, D_{2}, \cdots$，直到$D_m$。那么我们如何计算$D_{m+1}, D_{m+2}$？

自然的，我们需要丢到一些向量，那么肯定优先考虑丢掉那些最早生成的向量。具体地说，计算$D_{m+1}$时，我们保存$\left\{s_{i}, y_{i}\right\}_{i=1}^{m}$，丢掉了$\left\{s_{0}, y_{0}\right\}$；计算$D_{m+2}$时，我们保存$\left\{s_{i}, y_{i}\right\}_{i=2}^{m+1}$丢掉$\left\{s_{1}, y_{1}\right\}$。

但是舍弃掉一些向量后，我们就无法进行精确计算，只能够对其进行近似，当$k+1>m$时，我们可以构造近似计算公式
$$
\begin{align}
\begin{split}
D_{k+1}=&(V_{k}^{T} V_{k-1}^{T} \cdots V_{k-m+2}^{T} V_{k-m+1}^{T}) D_{0}(V_{k-m+1} V_{k-m+2} \cdots V_{k-1} V_{k})\\
+&(V_{k}^{T} V_{k-1}^{T} \cdots V_{k-m+3}^{T} V_{k-m+2}^{T})(\rho_{0} \mathbf{S}_{0} \mathbf{S}_{0}^{T})(V_{k-m+2} V_{k-m+3} \cdots V_{k-1} V_{k})\\
+&(V_{k}^{T} V_{k-1}^{T} \cdots V_{k-m+4}^{T} V_{k-m+3}^{T})(\rho_{1} \mathbf{s}_{1} \mathbf{s}_{1}^{T})(V_{k-m+3} V_{k-m+4} \cdots V_{k-1} V_{k})\\
+&\cdots\\
+&(V_{k}^{T} V_{k-1}^{T})(\rho_{k-2} \mathbf{s}_{k-2} \mathbf{s}_{k-2}^{T})(V_{k-1} V_{k})\\
+&V_{k}^{T}(\rho_{k-1} \mathbf{S}_{k-1} \mathbf{s}_{k-1}^{T}) V_{k}\\
+&\rho_{k} \mathbf{s}_{k} \mathbf{s}_{k}^{T}
\end{split}
\end{align}
$$

Algorithm L-BFGS Method
input starting point $x_0$, integer history size $m>0$(generally m=10), k=1
output the position $x$ with a minimal objective function
1: while no converge do
2: Calculate gradient $\nabla f(x_k)$ at point $x_k$
3: Compute direction $d^k$ using $\mathbf{d}_{k}=-D_{k} \cdot \mathbf{g}_{k}$
4: Compute $x_{k+1} = x_k +\lambda_kd_k$ where $\lambda_k$ is chose to satisfy the condition $f(x_k +\lambda_kd_k) = \min_{\lambda \geq 0}f(x_k +\lambda d_k)$
5: if $k>m$ then
6: Discard vector pair $p_{k-m}, q_{k-m}$ from memory storage
7: end if
8: update $p^{k}=x^{k+1}-x^{k}, q^{k}=\nabla f\left(x^{k+1}\right)-\nabla f\left(x^{k}\right), k=k+1$
9: end while

现在还有一个问题，那就是怎么计算下降方向$d_k$，当然我们可以直接根据算法求，矩阵$D_k$只用于计算$D_k\mathbf{g}_k$从而获取搜索方向。因此只要能够利用上式设计出一种计算$D_k\mathbf{g}_k$的算法就好。一种更为常用，更为方便的方式是利用two-loop recursion方法来计算。详细算法如下

Algorithm L-BFGS two-loop recursion
input $\nabla f(x_k), p_{i}, q_{i} \text { where } i=k-m, \ldots, k-1$
output new direction $d$
1: $d=-\nabla f\left(x_{k}\right)$
2: for $i=k-1 \text { to } k-m$ do
3: $\alpha^{i}=\frac{\left(p_{i}\right)^{T} d}{\left(q_{i}\right)^{T} p_{i}}$
4: $d=d-\alpha_{i} q_{i}$
5: end for
6: $d=\frac{\left(p_{k-1}\right)^{T} q_{k-1}}{\left(q_{k-1}\right)^{T} q_{k-1}} d$
7: for $i=k-m \text { to } k-1$ do
8: $\beta=\frac{\left(q_{i}\right)^{T} d}{\left(q_{i}\right)^{T} p_{i}}$
9: $d=d+\left(\alpha_{i}-\beta\right)$
10: end for

这种算法中在step 6中使用的公式就是利用了如下经验公式
$$
r_{k}=\frac{p_{k-1}^{T} q_{k-1}}{q_{k-1}^{T} q_{k-1}}
$$
最后计算出的$\mathbf{r}_L$即为$D_{k} \cdot \mathbf{g}_{k}$的值。其实在这个算法中的$\mathbf{r}_k$表示比例系数，它利用最近一次的曲率信息来估计真是海森矩阵的大小，这就使得当前步的搜索方向较为理想。这样就省去了步长搜索的步骤。

至此，简略的牛顿法与主要的拟牛顿法已经全部介绍完毕。

参考文献