优化算法（一）：牛顿法与拟牛顿法

机器学习算法中经常碰到非线性优化问题，如 Sparse Filtering 算法，其主要工作在于求解一个非线性极小化问题。在具体实现中，大多调用的是成熟的软件包做支撑，其中最常用的一个算法是 L-BFGS。为了解这个算法的数学机理，这几天做了一些调研，现把学习过程中理解的一些东西整理出来。

拟牛顿法（Quasi-Newton Methods）是求解非线性优化问题最有效的方法之一，在20世纪50年代由美国Argonne国家实验室的物理学家W. C. Davidon提出。Davidon设计的这种算法在当时看来是非线性优化领域最具有创造性的发明之一。不就R. Fletcher和M. J. D. Powell证实了这种新的算法远比其他方法快速和可靠，使得非线性优化这门学科在一夜之间突飞猛进。在之后的20年里，拟牛顿方法得到了蓬勃发展，出现了大量的变形公式以及数以百计的相关论文。

DFP方法、BFGS方法以及L-BFGS算法都是重要的拟牛顿法。本文将对这些方法进行简要介绍，当然，在介绍拟牛顿法之前，我们先看看什么是牛顿法，以及拟牛顿法和牛顿法之间有什么关系，为此，考虑如下无约束的极小化问题：
$$mi{n}_{x}f(x)\text{\hspace{1em}(0.1)}$$
其中$\mathbf{x}=(x_1,x_2,...,x_N{)}^T\in {\mathbb{R}}^N$.由于本文不准备对收敛性进行讨论，因此不妨对目标函数$f:{\mathbb{R}}^N\to \mathbb{R}$作一个比较苛刻的假设，这里我们假定$f$是凸函数，且两阶连续可微。此外，记绩效问题(0.1)的解为$x^{\ast }$

牛顿法

原始牛顿法

为简单起见，先考虑$N=1$的简单情形，此时目标函数$f(\mathbf{x})变为$$f(x)$。牛顿法的基本思想是：在现有极小点估计值的附近对$f(x)$做二阶泰勒展开，进而找到极小点的下一个估计值。设$x_k$为当前的极小点估计值，则
$$\phi (x)\approx f(x_k)+f^{\prime }(x_k)(x-x_k)+\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2\text{\hspace{1em}(1.2)}$$
表示$f(x)$在$x_k$附近的二阶泰勒展开式(略去了关于$x-x_k$的高阶项). 由于求得是最值，由极值必要条件可知，$\phi (x)$应该满足：
$${\phi }^{\prime }(x)=0\text{\hspace{1em}(1.3)}$$
即
$$f^{\prime }(x_k)+f^{\prime \prime }(x_k)(x-x_k)=0\text{\hspace{1em}(1.4)}$$
从而求得
$$x=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}\text{\hspace{1em}(1.5)}$$

于是，若给定初始值$x_0$,则可以构造如下的迭代格式：
$$x_{k+1}=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)},k=0,1,...\text{\hspace{1em}(1.6)}$$
产生序列$x_k$来逼近$f(x)$的极小点。在一定条件下，$x_k$可以收敛到$f(x)$的极小点。
对于$N>1$的情形，二阶泰勒展开式(1.2)可以做推广，此时
$$\phi (\mathbf{x})\approx f({\mathbf{x}}_k)+\nabla f({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^T{\nabla }^{2}f({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)\text{\hspace{1em}(1.7)}$$
其中$\nabla f$为$f$的梯度向量，${\nabla }^{2}f$为$f$的海森矩阵(Hessian Matrix)，其定义分别为：
$$\nabla f=\left[\begin{array}{c}\frac{\partial f}{\partial x_1}\\ \frac{\partial f}{\partial x_2}\\ ...\\ \frac{\partial f}{\partial x_N}\end{array}\right],{\nabla }^{2}f={\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1{x}_2}& ...& \frac{{\partial }^{2}f}{\partial x_1{x}_N}\\ \frac{{\partial }^{2}f}{\partial x_2{x}_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& ...& \frac{{\partial }^{2}f}{\partial x_2{x}_N}\\ ...& ...& ...& ...\\ \frac{{\partial }^{2}f}{\partial x_N{x}_1}& \frac{{\partial }^{2}f}{\partial x_N{x}_2}& ...& \frac{{\partial }^{2}f}{\partial x_N^2}\end{array}\right]}_{N\times N}\text{\hspace{1em}(1.8)}$$
注意，$\nabla f$和${\nabla }^{2}f$中的元素均为关于$\mathbf{x}$的函数，以下分别将其记为$\mathbf{g}$和$H$。特别的，若$f$的混合偏导数可以交换次序（即对任意$i,j$，成立$\frac{{\partial }^{2}f}{\partial x_i\partial x_j}=\frac{{\partial }^{2}f}{\partial x_j\partial x_i}$）,则海森矩阵$H$为对称矩阵，而$\nabla f({\mathbf{x}}_k)$和${\nabla }^{2}f({\mathbf{x}}_k)$则表示将$\mathbf{x}$取为${\mathbf{x}}_k$后得到的实值向量和矩阵，以下分别将其记为${\mathbf{g}}_k$和$H_k$(这里字母$\mathbf{g}$表示gradient，$H$表示Hessian)

同样的，由于是求极小点，极值必要条件要求它为$\phi (x)$的驻点，即
$$\nabla \phi (x)=0\text{\hspace{1em}(1.9)}$$
亦即（通过在(1.7)两边作用一个梯度算子）
$${\mathbf{g}}_k+H_k\cdot (\mathbf{x}-{\mathbf{x}}_k)=0\text{\hspace{1em}(1.10)}$$
进一步，若矩阵$H_k$非奇异，则可解得
$$\mathbf{x}={\mathbf{x}}_k-H_k^{-1}\cdot {\mathbf{g}}_k\text{\hspace{1em}(1.11)}$$
于是，若给定初始值${\mathbf{x}}_0$，则同样可以构造出迭代格式
$${\mathbf{x}}_{k+1}={\mathbf{x}}_k-H_k^{-1}\cdot {\mathbf{g}}_k,k=0,1,...\text{\hspace{1em}(1.12)}$$
这就是原始的牛顿迭代法，其迭代格式中的搜索方向${\mathbf{d}}_k=H_k^{-1}\cdot {\mathbf{g}}_k$称为牛顿方向。下面给出一个牛顿法的完整算法描述：

---

算法 1.1 （牛顿法）

1、给定初始值${\mathbf{x}}_0$和精度阈值$ϵ$，并令$k:=0$

2、计算${\mathbf{g}}_k$和$H_k$

3、若$||{\mathbf{g}}_k||<ϵ$，则停止迭代；否则确定搜索方向${\mathbf{d}}_k=-H_k^{-1}\cdot {\mathbf{g}}_k$

4、计算新的迭代点${\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\mathbf{d}}_k$

5、令$k:=k+1$，转至第2步

---

当目标函数是二次函数时，由于二次泰勒展开函数与原目标函数不是近似而是完全相同的二次式，海森矩阵退化成一个常数矩阵，从任一初始点出发，利用(1.12)只需要一步迭代即可到达$f(\mathbf{x})$的极小点$x^{\ast }$，因此牛顿法是一种具有二次收敛性的算法，对于非二次函数，若函数的二次性态较强，或迭代点已进入极小点的领域，其收敛速度也是很快的，这是牛顿法的主要优点。

但是牛顿法由于迭代公式中没有步长因子，而是定步长迭代，对于非二次性目标函数有时会使得函数值上升，即出现$f({\mathbf{x}}_{k+1})>f({\mathbf{x}}_k)$的情况，这表明原始牛顿法不能保证函数值稳定的下降，在严重的情况下甚至可能造成迭代点列 { x k } \{\mathbf{x}_k\} {xk​}的发散而导致计算失败。

阻尼牛顿法

为了消除牛顿法中的弊病，人们提出了“阻尼牛顿法”。阻尼牛顿法每次迭代的方向任然是${\mathbf{d}}_k$，但每次迭代需沿着此方向作一维搜索(line search),寻求最优的步长因子${\lambda }_k$,即
$${\lambda }_k=argmi{n}_{\lambda \in \mathbb{R}}f({\mathbf{x}}_k+\lambda {\mathbf{d}}_k)\text{\hspace{1em}(1.13)}$$

下面给出一个阻尼牛顿法的完整算法描述。

---

算法 1.2 （阻尼牛顿法）

1、给定初始值${\mathbf{x}}_0$和精度阈值$ϵ$，并令$k:=0$

2、计算${\mathbf{g}}_k$和$H_k$

3、若$||{\mathbf{g}}_k||<ϵ$，则停止迭代；否则确定搜索方向${\mathbf{d}}_k=-H_k^{-1}\cdot {\mathbf{g}}_k$

4、利用（1.13）得到步长${\lambda }_k$，计算新的迭代点${\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k$

5、令$k:=k+1$，转至第2步

---

注1.1 算法1.3的步3中设计到$H_k^{-1}$的计算，实际应用中，通常并不直接对$H_k$进行求逆，而是将其转化为求解线性代数方程组$H_k{\mathbf{d}}_k=-{\mathbf{g}}_k$，此时可根据系数矩阵$H_k^{-1}$的性态来选择适合的迭代法，如预条件共轭梯度法（PCG）、代数多重网格法（AMG）等。

注1.2 有些文献资料里将算法1.2步3中的搜索方向写成${\mathbf{d}}_k=H_k^{-1}\cdot {\mathbf{g}}_k$，这里是没问题的，原因是在4中有一个求步长的过程，将搜索方向写成了${\mathbf{d}}_k=H_k^{-1}\cdot {\mathbf{g}}_k$，无非是求得的最佳步长和原来求得的最佳步长相差一个符号罢了。

至此已完成了牛顿法的算法介绍，接下来对其做个小结：

牛顿法是梯度（下降）法的进一步发展，梯度法利用目标函数的一阶偏导数信息，以负梯度方向作为搜索方向，只考虑目标函数在迭代点的局部性质；而牛顿法不仅使用目标函数的一阶偏导数，还进一步利用目标函数的二阶偏导数，这样就考虑了梯度变化的趋势，因而能更全面的确定合适的搜索方向以加快收敛，它具有二阶收敛性。但牛顿法主要存在以下两个缺点：

• 对目标函数有严格的要求。函数必须具有连续的一、二阶偏导数，海森矩阵必须正定。
• 计算相当复杂，除需计算梯度而外，还需计算二阶偏导矩阵和它的逆矩阵。计算量、存储量均很大，且均已维数$N$的平方比增加，当N很大时这个问题更加突出。

拟牛顿法

如上节所说，牛顿法虽然收敛速度很快，但是计算过程中需要计算目标函数的二阶偏导数，计算复杂度较大。而且有时目标函数的海森矩阵无法保持正定，从而使得牛顿法失效。为了克服这两个问题，人们提出了拟牛顿法。这个方法的基本思想是：不用二阶偏导数而构造出可以近似还森矩阵（或海森矩阵的逆）的正定对称阵，在“拟牛顿”的条件下优化目标函数。不同的构造方法就产生了不同的拟牛顿法。

也有人把“拟牛顿法”翻译成“准牛顿法”。

在介绍具体的拟牛顿法之前，我们先推导一个拟牛顿条件，或者叫做拟牛顿方程，还有的叫做割线条件，因为对海森矩阵（或海森矩阵的逆）做近似总不能随便近似吧，我们也需要理论指导，而拟牛顿条件则是用来提供理论指导的，它指出了用来近似的矩阵应该满足的条件。

为明确起见，下文中用$B$表示对海森矩阵$H$本身的近似，而用$D$表示对还森矩阵的逆$H^{-1}$的近似，即$B\approx H,D\approx H^{-1}$

拟牛顿条件

设经过$k+1$次迭代后得到${\mathbf{x}}_{k+1}$，此时将目标函数$f(\mathbf{x})$在${\mathbf{x}}_{k+1}$附近作泰勒展开，取二阶近似，得到
$$f(\mathbf{x})\approx f({\mathbf{x}}_{k+1})+\nabla f({\mathbf{x}}_{k+1})(\mathbf{x}-{\mathbf{x}}_{k+1})$$

$$+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_{k+1}{)}^T{\nabla }^{2}f({\mathbf{x}}_{k+1})(\mathbf{x}-{\mathbf{x}}_{k+1})\text{\hspace{1em}(2.14)}$$

在（2.14）两边同时作用一个梯度算子$\nabla$，可得
$$\nabla f(\mathbf{x})\approx \nabla f({\mathbf{x}}_{k+1})+H_{k+1}\cdot (\mathbf{x}-{\mathbf{x}}_{k+1})\text{\hspace{1em}(2.15)}$$

在（2.15）中取$\mathbf{x}={\mathbf{x}}_k$，并整理，可得
$${\mathbf{g}}_{k+1}-{\mathbf{g}}_k\approx H_{k+1}\cdot ({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)\text{\hspace{1em}(2.16)}$$
若引入记号
$${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k,{\mathbf{y}}_k={\mathbf{g}}_{k+1}-{\mathbf{g}}_k\text{\hspace{1em}(2.17)}$$
则（2.16）可紧凑地写成
$${\mathbf{y}}_k\approx H_{k+1}\cdot {\mathbf{s}}_k\text{\hspace{1em}(2.18)}$$
或者
$${\mathbf{s}}_k\approx H_{k+1}^{-1}\cdot {\mathbf{y}}_k\text{\hspace{1em}(2.19)}$$

这就是所谓的拟牛顿条件，它对迭代过程中的海森矩阵$H_{k+1}$作约束，因此。对$H_{k+1}$做近似的$B_{k+1}$，以及对$H_{k+1}^{-1}$做近似的$D_{k+1}$可将
$${\mathbf{y}}_k\approx B_{k+1}\cdot {\mathbf{s}}_k\text{\hspace{1em}(2.20)}$$
或者
$${\mathbf{s}}_k\approx D_{k+1}\cdot {\mathbf{y}}_k\text{\hspace{1em}(2.21)}$$
作为指导。

接下来，我们依次介绍几种常见的拟牛顿法。

DFP算法

DFP算法是以William C. Davidon、Roger Fletcher、Micha J. D. Powell三个人的名字的首字母命名的，它由Davidon于1959年首先提出，后经Fletcher和Powell加以发展和完善，是最早的拟牛顿法。该算法的核心是：通过迭代的方法，对$H_{k+1}^{-1}$做近似，迭代格式为
$$D_{k+1}=D_k+\Delta D_k,k=0,1,2,...\text{\hspace{1em}(2.22)}$$
其中的$D_0$通常取单位矩阵$I$，因此，关键是每一步的校正矩阵$\Delta D_k$如何构造。

注意，迭代格式（2.22）将嵌套在算法1.2中，因此，我们猜想$\Delta D_k$可能与${\mathbf{s}}_k,{\mathbf{y}}_k$和$D_k$发生关联。这里，我们采用“待定法”，即首先将$\Delta D_k$待定成某种形式，然后结合拟牛顿条件（2.21）来进行推导。

那将$\Delta D_k$待定成什么形式呢？这个说起来比较tricky，我们将其待定为
$$\Delta D_k=\alpha \mathbf{u}{\mathbf{u}}^T+\beta \mathbf{v}{\mathbf{v}}^T\text{\hspace{1em}(2.23)}$$
其中$\alpha$和$\beta$为待定系数，$\mathbf{u},\mathbf{v}\in {\mathbb{R}}^N$为待定向量。从形式上看，这种待定公式至少保证了矩阵$\Delta D_k$的对称性（因为$\mathbf{u}{\mathbf{u}}^T$和$\mathbf{v}{\mathbf{v}}^T$均为对称矩阵）

将（2.23）带入（2.22），并结合指导条件（2.21），可得
$${\mathbf{s}}_k=D_k{\mathbf{y}}_k+\alpha \mathbf{u}{\mathbf{u}}^T{\mathbf{y}}_k+\beta \mathbf{v}{\mathbf{v}}^T{\mathbf{y}}_k\text{\hspace{1em}(2.24)}$$

从（2.24）似乎也看不出什么啊！别急，我们将其改写一下：
$${\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}\text{\hspace{1em}(2.25)}$$

看到了吧？括号中的$\alpha {\mathbf{u}}^T{\mathbf{y}}_k$和$\beta {\mathbf{v}}^T{\mathbf{y}}_k$是两个数，既然是数，我们不妨做如下简单赋值
$$\alpha {\mathbf{u}}^T{\mathbf{y}}_k=1,\beta {\mathbf{v}}^T{\mathbf{y}}_k=-1\text{\hspace{1em}(2.26)}$$
即
$$\alpha =\frac{1}{{\mathbf{u}}^T{\mathbf{y}}_k},\beta =-\frac{1}{{\mathbf{v}}^T{\mathbf{y}}_k}\text{\hspace{1em}(2.27)}$$

其中向量$\mathbf{u},\mathbf{v}$仍有待确定。

那么$\mathbf{u},\mathbf{v}$如何确定呢？将（2.26）带入（2.25），得到
$$\mathbf{u}-\mathbf{v}={\mathbf{s}}_k-D_k{\mathbf{y}}_k\text{\hspace{1em}(2.28)}$$

要上式成立，不放直接取
$$\mathbf{u}={\mathbf{s}}_k,\mathbf{v}=D_k{\mathbf{y}}_k\text{\hspace{1em}(2.29)}$$

再将（2.29）带入（2.27），便得到
$$\alpha =\frac{1}{{\mathbf{s}}_k^T{\mathbf{y}}_k},\beta =-\frac{1}{(D_k{\mathbf{y}}_k{)}^T{\mathbf{y}}_k}=-\frac{1}{{\mathbf{y}}_k^T{D}_k{\mathbf{y}}_k}\text{\hspace{1em}(2.30)}$$

其中第二个等式用到了$D_k$的对称性。
至此，我们已经将校正矩阵$\Delta D_k$构造出来了，将（2.29）和（2.30）带入（2.23），便得
$$\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}\text{\hspace{1em}(2.31)}$$
综上，我们给出了DFP算法的一个完整算法描述。

---

算法 2.1 （DFP算法）

1、给定初始值${\mathbf{x}}_0$和精度阈值$ϵ$，并令$D_0=I,k:=0$

2、确定搜索方向${\mathbf{d}}_k=-D_k\cdot {\mathbf{g}}_k$

3、利用（1.13）得到步长${\lambda }_k$,令${\mathbf{s}}_k={\lambda }_k{\mathbf{d}}_k,{\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\mathbf{s}}_k$

4、若$||{\mathbf{g}}_{k+1}||<ϵ$，则算法结束；

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$，转至第2步

---

BFGS算法

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

BFGS算法中核心公式的推导过程和DFP完全类似，只是互换了其中${\mathbf{s}}_k$和${\mathbf{y}}_k$的位置。为了方便自己以后查阅，我打算将上一节的推导过程再重写一遍，已经领会该过程的读者不妨直接跳过以下推导，直接看结论。

需要注意的是，BFGS算法是直接逼近海森矩阵，即$B_k\approx H_k$.仍采用迭代方法，设迭代格式为
$$B_{k+1}=B_k+\Delta B_k,k=0,1,2,...\text{\hspace{1em}(2.32)}$$

其中的$B_0$也取单位矩阵$I$，因此，关键是每一步的校正矩阵$\Delta B_k$如何构造。同样的，将其待定为：

$$\Delta B_k=\alpha \mathbf{u}{\mathbf{u}}^T+\beta \mathbf{v}{\mathbf{v}}^T\text{\hspace{1em}(2.33)}$$

将（2.33）带入（2.32），并结合指导条件（2.20），可得
$${\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}\text{\hspace{1em}(2.34)}$$

通过令$\alpha {\mathbf{u}}^T{\mathbf{s}}_k=1,\beta {\mathbf{v}}^T{\mathbf{s}}_k=-1$,以及
$$\mathbf{u}={\mathbf{y}}_k,\mathbf{v}=B_k{\mathbf{s}}_k\text{\hspace{1em}(2.35)}$$
可算得
$$\alpha =\frac{1}{{\mathbf{y}}_k^T{\mathbf{s}}_k},\beta =-\frac{1}{{\mathbf{s}}_k^T{B}_k{\mathbf{s}}_k}\text{\hspace{1em}(2.36)}$$
综上，便得到了如下校正矩阵$B_k\approx H_k$的公式
$$\Delta 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}\text{\hspace{1em}(2.31)}$$

好了，现在把矩阵$B_k\approx H_k$和$D_k\approx H_k$拿出来对比一下，是不是除了$D$换成$B$外，其它只是将${\mathbf{s}}_k$和${\mathbf{y}}_k$互换了一下位置呢？

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

---

算法 2.2 （BFGS算法（1））

1、给定初始值${\mathbf{x}}_0$和精度阈值$ϵ$，并令$B_0=I,k:=0$

2、确定搜索方向${\mathbf{d}}_k=-B_k^{-1}\cdot {\mathbf{g}}_k$

3、利用（1.13）得到步长${\lambda }_k$,令${\mathbf{s}}_k={\lambda }_k{\mathbf{d}}_k,{\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\mathbf{s}}_k$

4、若$||{\mathbf{g}}_{k+1}||<ϵ$，则算法结束；

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$，转至第2步

---

算法2.2中的步2通常是哦通过求解线性代数方程组$B_k{\mathbf{d}}_k=-{\mathbf{g}}_k$来进行。然而，更一般的做法是，通过对步6中的递推关系应用Sherman-Morrion公式，直接给出$B_{k+1}^{-1}$与$B_k^{-1}$之间的关系式
$$B_{k+1}^{-1}=(I-\frac{{\mathbf{s}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k})B_k^{-1}(I-\frac{{\mathbf{y}}_k{\mathbf{s}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k})+\frac{{\mathbf{s}}_k{\mathbf{s}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}\text{\hspace{1em}(2.38)}$$

或者进一步展开写成
$$B_{k+1}^{-1}=B_k^{-1}$$

$$+(\frac{1}{{\mathbf{s}}_k^T{\mathbf{y}}_k}+\frac{{\mathbf{y}}_k^T{B}_k^{-1}{\mathbf{y}}_k}{({\mathbf{s}}_k^T{\mathbf{y}}_k{)}^2}){\mathbf{s}}_k{\mathbf{s}}_k^T-\frac{1}{{\mathbf{s}}_k^T{\mathbf{y}}_k}(B_k^{-1}{\mathbf{y}}_k{\mathbf{s}}_k^T+{\mathbf{s}}_k{\mathbf{y}}_k^T{B}_k^{-1})\text{\hspace{1em}(2.39)}$$

其中$(\frac{1}{{\mathbf{s}}_k^T{\mathbf{y}}_k}+\frac{{\mathbf{y}}_k^T{B}_k^{-1}{\mathbf{y}}_k}{({\mathbf{s}}_k^T{\mathbf{y}}_k{)}^2})$和$\frac{1}{{\mathbf{s}}_k^T{\mathbf{y}}_k}$是实数。

注2.1 关于Sherman-Morrison公式
设$A\in {\mathbb{R}}^n$为非奇异方阵，$\mathbf{u},\mathbf{v}\in {\mathbb{R}}^n$，若$1+{\mathbf{v}}^T{A}^{-1}\mathbf{u}\ne 0$，则有
$$(A+\mathbf{u}{\mathbf{v}}^T{)}^{-1}=A^{-1}-\frac{A^{-1}\mathbf{u}{\mathbf{v}}^T{A}^{-1}}{1+{\mathbf{v}}^T{A}^{-1}\mathbf{u}}\text{\hspace{1em}(2.40)}$$*

利用（2.38），我们很容易将算法2.2改写成2.3. 注意，为了避免出现矩阵求你符号，我们统一将$B_i^{-1}$换成$D_i$(这样做仅仅只是符号上看起来舒服起见)。这样，整个算法中不再需要求解线性代数方程组，由矩阵-向量运算就可以完成了。

---

算法 2.3 （BFGS算法（2））

1、给定初始值${\mathbf{x}}_0$和精度阈值$ϵ$，并令$D_0=I,k:=0$

2、确定搜索方向${\mathbf{d}}_k=-D_k\cdot {\mathbf{g}}_k$

3、利用（1.13）得到步长${\lambda }_k$,令${\mathbf{s}}_k={\lambda }_k{\mathbf{d}}_k,{\mathbf{x}}_{k+1}:={\mathbf{x}}_k+{\mathbf{s}}_k$

4、若$||{\mathbf{g}}_{k+1}||<ϵ$，则算法结束；

5、计算${\mathbf{y}}_k={\mathbf{g}}_{k+1}-{\mathbf{g}}_k$

6、计算$D_{k+1}=(I-\frac{{\mathbf{s}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k})D_k(I-\frac{{\mathbf{y}}_k{\mathbf{s}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k})+\frac{{\mathbf{s}}_k{\mathbf{s}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}$

7、令$k:=k+1$，转至第2步

---

至此，关于DFP算法和BFGS算法的介绍就完成了，回过头来，我们对比一下算法2.1和算法2.3，容易发现，这两个算法的唯一不同仅在于$D_{k+1}$的迭代公式不同罢了。

最后，再补充谈谈一维搜索（line search）的问题，在之前几个算法描述中，为简单起见，均采用（1.13）来计算步长${\lambda }_k$，其实这是一种精确搜索。实际应用中，还有像Wolfe型搜索，Armijo搜索以及满足Goldstein条件的非精准搜索。这里我们以Wolfe搜索为例，简单做个介绍。

设$\tilde{\beta }\in (0,\frac{1}{2}),\beta \in (\tilde{\beta },1)$，所谓的Wolfe搜索是指${\lambda }_k$满足如下Wolfe条件
$$\{\begin{array}{rl}f({\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k)& \le f({\mathbf{x}}_k)+\tilde{\beta }{\lambda }_k{\mathbf{d}}_k^T{\mathbf{g}}_k\\ {\mathbf{d}}_k^T\mathbf{g}({\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k)& \ge \beta {\mathbf{d}}_k^T{\mathbf{g}}_k\end{array}\text{\hspace{1em}(2.41)}$$

带非精确搜索的拟牛顿法的研究是从1976年Powell的工作开始的，他证明了带Wolfe搜索的BFGS算法的全局收敛性和超线性收敛性。

L-BFGS算法

在BFGS算法中，需要用到一个$N\times N$的矩阵$D_k$,当$N$很大时，存储这个矩阵将变得很耗计算机资源。例如，考虑$N$个10万的情形，且用double型（8字节）在存储$D_k$，需要多大的内存呢？我们来计算一下
$$\frac{N阶矩阵的字节数}{1GB的字节数}=\frac{10^5\times 10^5\times 8}{2^{10}\times 2^{10}\times 2^{10}}=74.5(GB)\text{\hspace{1em}(2.42)}$$

74.5GB，很惊人是吧，这对于一般的服务器是很难承受的。当然，考虑到矩阵$D_k$的对称性，内存还可以降一半，但是，在机器学习问题中，像10万这样的规模还只能算是中小规模。那么，是否可以通过对BFGS算法进行改造，从而减少其迭代过程中所需的内存开销呢？

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

接下来，讨论L-BFGS算法的具体实现过程。我们的出发点是算法2.3步6的迭代式
$$D_{k+1}=(I-\frac{{\mathbf{s}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k})D_k(I-\frac{{\mathbf{y}}_k{\mathbf{s}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k})+\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{s}}_k{\mathbf{s}}_k^T$，则上式可写成
$$D_k+1=V_k^T{D}_k{V}_k+{\rho }_k{\mathbf{s}}_k{\mathbf{s}}_k^T\text{\hspace{1em}(2.43)}$$

若给定初始矩阵$D_0$(通常为正定的对角矩阵，如$D_0=I$)，则利用（2.43），依次可得

$$\begin{array}{rl}{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{\mathbf{s}}_1{\mathbf{s}}_1^T\\ & =V_1^T(V_0^T{D}_0{V}_0+{\rho }_0{\mathbf{s}}_0{\mathbf{s}}_0^T)V_1+{\rho }_1{\mathbf{s}}_1{\mathbf{s}}_1^T\\ & =V_1^T{V}_0^T{D}_0{V}_0{V}_1+V_1^T{\rho }_0{\mathbf{s}}_0{\mathbf{s}}_0^T{V}_1+{\rho }_1{\mathbf{s}}_1{\mathbf{s}}_1^T\\ D_3& =V_2^T{D}_2{V}_2+{\rho }_2{\mathbf{s}}_2{\mathbf{s}}_2^T\\ & =V_2^T(V_1^T{V}_0^T{D}_0{V}_0{V}_1+V_1^T{\rho }_0{\mathbf{s}}_0{\mathbf{s}}_0^T{V}_1+{\rho }_1{\mathbf{s}}_1{\mathbf{s}}_1^T)V_2+{\rho }_2{\mathbf{s}}_2{\mathbf{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{\mathbf{s}}_0{\mathbf{s}}_0^T{V}_1{V}_2+V_2^T{\rho }_1{\mathbf{s}}_1{\mathbf{s}}_1^T{V}_2+{\rho }_2{\mathbf{s}}_2{\mathbf{s}}_2^T\end{array}$$

一般的，我们有
$$\begin{array}{rl}{D}_{k+1}& =(V_k^T{V}_{k-1}^T...V_1^T{V}_0^T)D_0(V_0{V}_1...V_{k-1}{V}_k)\\ & +(V_k^T{V}_{k-1}^T...V_2^T{V}_1^T)({\rho }_0{\mathbf{s}}_0{\mathbf{s}}_0^T)(V_1{V}_2...V_{k-1}{V}_k)\\ & +(V_k^T{V}_{k-1}^T...V_3^T{V}_2^T)({\rho }_1{\mathbf{s}}_1{\mathbf{s}}_1^T)(V_2{V}_3...V_{k-1}{V}_k)\\ & +...\\ & +(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{array}\text{\hspace{1em}(2.44)}$$

由上式可见，计算$D_{k+1}$需要用到 { s i , y i } i = 0 k \{s_i,y_i\}^k_{i=0} {si​,yi​}i=0k​，因此，若从$s_0,y_0$开始连续的存储$m$组的话，智能存储到$s_{m-1},y_{m-1}$，亦即，只能依次计算$D_1,D_2,...,D_m$。那么$D_{m+1},D_{m+2}$该如何计算呢？

自然的，如果一定要丢掉一些向量，那么肯定优先考虑那些最早生成的向量。具体来说，计算$D_{m+1}$时，我们保存 { s i , y i } i = 1 m \{s_i,y_i\}^m_{i=1} {si​,yi​}i=1m​，丢掉了 { s 0 , y 0 } \{s_0,y_0\} {s0​,y0​};计算$D_{m+2}$时，我们保存 { s i , y i } i = 2 m + 1 \{s_i,y_i\}^{m+1}_{i=2} {si​,yi​}i=2m+1​，丢掉了 { s 0 , y 0 } i = 0 1 \{s_0,y_0\}^1_{i=0} {s0​,y0​}i=01​；…

但是舍弃掉一些向量后，就只能近似计算了。当$k+1>m$时，仿照（2.44），可以构造近似计算公式

$$\begin{array}{rl}{D}_{k+1}& =(V_k^T{V}_{k-1}^T...V_{k-m+2}^T{V}_{k-m+1}^T)D_0(V_{k-m+1}{V}_{k-m+2}...V_{k-1}{V}_k)\\ & +(V_k^T{V}_{k-1}^T...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}...V_{k-1}{V}_k)\\ & +(V_k^T{V}_{k-1}^T...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}...V_{k-1}{V}_k)\\ & +...\\ & +(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{array}\text{\hspace{1em}(2.45)}$$

（2.44）和（2.45）被称为Special BFGS Matrics。若引入 m ^ = m i n { k , m − 1 } \hat{m}=min\{k,m-1\} m^=min{k,m−1}，则还可以将（2.44）和（2.45）合并的写成

$$\begin{array}{rl}{D}_{k+1}& =(V_k^T{V}_{k-1}^T...V_{k-\hat{m}+1}^T{V}_{k-\hat{m}}^T)D_0(V_{k-\hat{m}}{V}_{k-\hat{m}+1}...V_{k-1}{V}_k)\\ & +(V_k^T{V}_{k-1}^T...V_{k-\hat{m}+2}^T{V}_{k-\hat{m}+1}^T)({\rho }_0{\mathbf{s}}_0{\mathbf{s}}_0^T)(V_{k-\hat{m}+1}{V}_{k-\hat{m}+2}...V_{k-1}{V}_k)\\ & +(V_k^T{V}_{k-1}^T...V_{k-\hat{m}+3}^T{V}_{k-\hat{m}+2}^T)({\rho }_1{\mathbf{s}}_1{\mathbf{s}}_1^T)(V_{k-\hat{m}+2}{V}_{k-\hat{m}+3}...V_{k-1}{V}_k)\\ & +...\\ & +(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{array}\text{\hspace{1em}(2.46)}$$

看到这里，千万不要被（2.46）冗长复杂的形式吓到，事实上，由BFGS算法流程易知，$D_k$的作用仅用来计算$D_k{\mathbf{g}}_k$获取搜索方向，因此，若能利用表达式（2.46）设计出一种计算$D_k{\mathbf{g}}_k$的快速算法，则大功告成。具体算法如下：

---

算法 2.4 （$D_k\cdot {\mathbf{g}}_k$）的快速算法

Step 1 初始化

$$\delta =\{\begin{array}{rlr}0,& & 若k\le m\\ k-m,& & 若k>m\end{array};L=\{\begin{array}{rlr}k,& & 若k\le m\\ m,& & 若k>m\end{array};{\mathbf{q}}_L={\mathbf{g}}_k$$

Step 2 后向循环$

  For $i=L-1,L-2,...,1,0$ DO
  {
   $j=i+\delta$;
   ${\alpha }_i={\rho }_j{\mathbf{s}}_j^T{\mathbf{q}}_{i+1}$; //${\alpha }_i$需要寻下来，前向循环要用！
   ${\mathbf{q}}_i={\mathbf{q}}_{i+1}-{\alpha }_i{\mathbf{y}}_j$.
  }

Step 3 前项循环

  ${\mathbf{r}}_0=D_0\cdot {\mathbf{q}}_0$;
  For $i=0,1,...,L-2,L-1$ DO
  {
   $j=i+\delta$;
   ${\beta }_i={\rho }_j{\mathbf{y}}_j^T{\mathbf{r}}_i$;
   ${\mathbf{r}}_{i+1}={\mathbf{r}}_i+({\alpha }_i-{\beta }_i){\mathbf{s}}_j$;

  }

---

最后计算出的${\mathbf{r}}_L$即为$H_k\cdot {\mathbf{g}}_k$的值

参考文献

牛顿法与拟牛顿法学习笔记（一）牛顿法
牛顿法与拟牛顿法学习笔记（二）拟牛顿条件
牛顿法与拟牛顿法学习笔记（三）DFP 算法
牛顿法与拟牛顿法学习笔记（四）BFGS 算法
牛顿法与拟牛顿法学习笔记（五）L-BFGS 算法