牛顿法和拟牛顿法

[TOC]

考虑无约束最优化问题:

其中,$\x=(x_1,x_2,\dots,x_N)^\T \in \R^N$,假设目标函数 $f: \R^N\to \R$为凸函数,且二阶连续可微。

牛顿法

先考虑 $N=1$ 的情况(此时目标函数为$f(x)$)。牛顿法的基本思想是:在现有极小点估计值的附近对$f(x)$做二阶泰勒展开,进而找到极小点的下一个估计值。

设 $x_k$ 为当前的极小点估计值,则

表示 $f(x)$ 在 $x_k$ 附近的二阶泰勒展开式(略去高阶项)。由于求极值,可知 $\varphi(x)$ 应该满足

第一项$f(x_k)$不包含变量$x$,求导直接略去,余下后两项即

从而求得

于是,若给定初始值 $x_0$,则可以构造如下迭代式:

产生序列 $\brace{x_k}$ 来逼近 $f(x)$ 的极小点。

对 $N \gt 1$ 的情况,二阶泰勒展开式可以推广,此时

其中,$\nabla f$为$f$的梯度向量,$\nabla^2f$为$f$的海森矩阵(Hessian matrix),其定义分别为

其中,$\nabla f$和$\nabla^2f$中的元素均为关于$\x$的函数,分别简记作$\g$和$\H$。特别的,若$f$的混合偏导数可以交换次数(即对$\forall i,j$,有$\frac{\partial^2 f}{\partial x_j\partial x_j}=\frac{\partial^2 f}{\partial x_j\partial x_i}$成立),则海森矩阵$\H$为对称矩阵。 而$\nabla f(\x_k)$和$\nabla^2f(\x_k)$则表示将$\x$取为$\x_k$后得到的实值向量和矩阵,分别简记作$\g_k$和$\H_k$。

同样的,由于是求极小点,极值必要条件要求它为$\varphi(\x)$的驻点,即

亦即在二阶泰勒展开式两边作用一个梯度算子,可得

进一步,若矩阵$\H_k$非奇异,则可以解得

于是,弱给定初始值$\x_0$,则同样可以构造出迭代式

这就是原始的牛顿迭代法,其中迭代式中的搜索方向 $\d_k=-\H_k^{-1}\cdot\g_k$ 称为牛顿方向

牛顿法的完整算法描述:

  1. 给定初值$\x_0$和精度阈值$\epsilon$,令$k:=0$,
  2. 计算 $\g_k$ 和 $\H_k$,
  3. 若 $\lvert \g_k \rvert\lt\epsilon$,则停止迭代;否则确定搜索方向$\d_k$,
  4. 计算新的迭代点$\x_{k+1}:=\x_k+\d_k$,
  5. 令$k:=k+1$,转至步骤2。

当目标函数是二次函数时,由于二次泰勒展开式与原目标函数不是近似而是完全相同的二次式,$\H$退化成一个常数矩阵,从任一初始点触发,只需一步迭代即可达到$f(\x)$的极小点,因此牛顿法是一种具有二次收敛性的算法。

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

阻尼牛顿法

为了消除原始牛顿法定步长迭代的弊病,阻尼牛顿法每次迭代方向仍使用$\d_k$,但每次迭代需要沿此方向作线搜索(line search),寻求最优步长因子$\lambda_k$,即

阻尼牛顿法的算法描述:

  1. 给定初值$\x_0$和精度阈值$\epsilon$,令$k:=0$,
  2. 计算 $\g_k$ 和 $\H_k$,
  3. 若 $\lvert \g_k\rvert\lt\epsilon$,则停止迭代;否则确定搜索方向$\d_k$,
  4. 计算得到步长$\lambda_k$,并令 $\x_{k+1}:=\x_k+\lambda_k\d_k$,
  5. 令$k:=k+1$,转至步骤2。

小结

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

  • 对目标函数有较严格的要求,必须具备连续的一、二阶偏导数,海森矩阵必须正定。
  • 计算相当负责,除需计算梯度外,还需要计算二阶偏导数矩阵和它的逆矩阵,计算量、存储量均很大,且均是$O(N^2)$的复杂度。

拟牛顿法

拟牛顿法基本思想:不用二阶偏导数而构造出可以近似海森矩阵(或海森矩阵的逆)的正定对称矩阵,在拟牛顿条件下优化目标函数。

不同的构造方法产生了不同的拟牛顿法。

符号:用$B$表示对海森矩阵$H$本身的近似,用$D$表示对$H^{-1}$的近似,即$B\approx H$,$D\approx H^{-1}$。

拟牛顿条件

设经过$k+1$次迭代后得到$\x_{k+1}$,此时将目标函数$f(\x)$在$\x_{k+1}$附近做泰勒展开,取二阶近似,得到

上式两边同时作用一个梯度算子$\nabla$,可得

在上式中取$\x=\x_k$,整理可得

引入新记号

则有

这就是所谓的拟牛顿条件,也叫割线条件(secant condition),它对迭代过程中的海森矩阵$H_{k+1}$作约束。因此,对$H_{k+1}$做近似的$B_{k+1}$,以及对$H_{k+1}^{-1}$ 做近似的$D_{k+1}$可以将

作为指导。

DFP

算法核心是通过迭代的方法,对$H_{k+1}^{-1}$做近似,迭代格式为

其中,$D_0$通常取单位矩阵$I$,关键是校正矩阵$\Delta D_k$的构造方法。

DFP算法:

  1. 给定初值$\x_0$和精度阈值$\epsilon$,并令$D_0=I$,$k:=0$,
  2. 确定搜索方向 $\d_k = -D_k\cdot \g_k$,
  3. 线搜索步长 $\lambda_k$,令$\s_k=\lambda_k \d_k$,$\x_{k+1}:=\x_k+\s_k$,
  4. 若 $\lvert\g_{k+1}\rvert\lt \epsilon$,则算法结束,
  5. 计算 $\y_k=\g_{k+1}-\g_k$,
  6. 计算 $D_{k+1}=D_k+{\s_k\s_k^\T \over \s_k^\T \y_k} - {D_k\y_k\y_k^\T D_k \over \y_k^\T D_k\y_k}$,
  7. 令$k:=k+1$,转至步骤2。

BFGS

BFGS和DFP区别是直接逼近$H$,而不是$H^{-1}$,其核心公式推导过程和DFP完全类似。迭代式如下:

BGGS算法:

  1. 给定初值$\x_0$和精度阈值$\epsilon$,并令$B_0=I$,$k:=0$,
  2. 确定搜索方向 $\d_k = -B_k^{-1}\cdot \g_k$,
  3. 线搜索步长 $\lambda_k$,令$\s_k=\lambda_k \d_k$,$\x_{k+1}:=\x_k+\s_k$,
  4. 若 $\lvert\g_{k+1}\rvert\lt \epsilon$,则算法结束,
  5. 计算 $\y_k=\g_{k+1}-\g_k$,
  6. 计算 $B_{k+1}=B_k+{\y_k\y_k^\T \over \y_k^\T \s_k} - {B_k\s_k\s_k^\T B_k \over \s_k^\T B_k\s_k}$,
  7. 令$k:=k+1$,转至步骤2。

以上算法第2步通常通过求解线性代数方程组 $B_k\d_k=-\g_k$来进行。然而,更一般的做法是通过对第6步的递推关系应用Sherman-Morrison公式,直接给出$B_{k+1}^{-1}$与$B_k^{-1}$之间的关系式:

为了避免出现矩阵求逆符号,统一将$B_i^{-1}$换为$D_i$,这样整个算法中不再需要求解线性代数方程组,由矩阵-向量运算就可以完成。

BFGS算法II:

  1. 给定初值$\x_0$和精度阈值$\epsilon$,并令$D_0=I$,$k:=0$,
  2. 确定搜索方向 $\d_k = -D_k\cdot \g_k$,
  3. 线搜索步长 $\lambda_k$,令$\s_k=\lambda_k \d_k$,$\x_{k+1}:=\x_k+\s_k$,
  4. 若 $\lvert\g_{k+1}\rvert\lt \epsilon$,则算法结束,
  5. 计算 $\y_k=\g_{k+1}-\g_k$,
  6. 计算 $D_{k+1}=\left(I-\frac{\s_k\y_k^\T }{\y_k^\T \s_k}\right) D_k \left(I-\frac{\y_k\s_k^\T }{\y_k^\T \s_k}\right) + \frac{\s_k\s_k^\T }{\y_k^\T \s_k}$,
  7. 令$k:=k+1$,转至步骤2。

L-BFGS

在BFGS算法中,需要用到$N\times N$矩阵$D_k$,当$N$很大时,该矩阵很耗资源。L-BFGS对BFGS算法进行了近似,其基本思想是:不再存储完整的矩阵$D_k$,而是存储计算过程中最近的$m$个向量序列$\brace{\s_i}$、$\brace{\y_i}$,需要矩阵$D_k$时,利用该向量序列的计算来替代。这样存储由$O(N^2)$降为$O(mN)$。

出发点是BFGS算法II第6步的迭代式

若记$\rho_k=\frac{1}{\y_k^\T \s_k}$,$V_k=I-\rho_k\y_k\s_k^\T $,则上式可写成

如果给定初始矩阵$D_0$(通常为正定对角矩阵,如$D_0=I$),利用上式依次可得

一般的,我们有

由上式可见,计算$D_{k+1}$需用到$\brace{s_i,y_i}{i=0}^k$。但我们可以丢掉一些最早的向量做近似计算,只由最近的$m$个向量计算$D{k+1}$,近似公式如下:

在文献Nocedal J. Updating quasi-Newton matrices with limited storage中,上式被称为special BFGS matrices。

参考

  • http://blog.csdn.net/itplus/article/details/21896453