牛顿法和拟牛顿法

Posted on |

[TOC]

\[\newcommand{\x}{\mathbf{x}} \newcommand{\R}{\mathbb{R}} \newcommand{\g}{\mathbf{g}} \newcommand{\H}{\mathit{H}} \newcommand{\d}{\mathbf{d}} \newcommand{\s}{\mathbf{s}} \newcommand{\y}{\mathbf{y}} \newcommand{\T}{\mathsf{T}}\]

考虑无约束最优化问题：

\[\min_\x f(\x)\]

其中，$\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$ 为当前的极小点估计值，则

\[\varphi(x) = f(x_k) + f'(x_k)(x-x_k) +\frac{1}{2}f''(x_k)(x-x_k)^2\]

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

\[\varphi'(x)=0\]

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

\[f'(x_k)+f''(x_k)(x-x_k)=0\]

从而求得

\[x=x_k-{f'(x_k) \over f''(x_k)}\]

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

\[x_{k+1}=x_k - {f'(x_k) \over f''(x_k)},\qquad k=0,1,\dots\]

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

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

\[\varphi(\x)=f(\x_k)+\nabla f(\x_k)\cdot(\x-\x_k)+\frac{1}{2}\cdot(\x-\x_k)^\T \cdot \nabla^2f(\x_k)\cdot(\x-\x_k)\]

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

\[\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}, \qquad \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} & \cdots & \frac{\partial^2 f}{\partial x_N^2} \end{bmatrix}_{N \times N}\]

其中，$\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)$的驻点，即

\[\nabla \varphi(\x)=0\]

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

\[\g_k+\H_k\cdot(\x-\x_k)=0\]

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

\[\x=\x_k - \H_k^{-1}\cdot \g_k\]

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

\[\x_{k+1}=\x_k-\H_k^{-1}\cdot\g_k,\qquad k=0,1,\cdots\]

这就是原始的牛顿迭代法，其中迭代式中的搜索方向 $\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$，即

\[\lambda_k= \arg \min_{\lambda \in \R} f(\x_k +\lambda \d_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}$附近做泰勒展开，取二阶近似，得到

\[f(\x)\approx f(\x_{k+1}) + \nabla f(\x_{k+1}) \cdot (\x-\x_{k+1}) + \frac{1}{2}\cdot (\x-\x_{k+1})^\T \cdot \nabla^2f(\x_{k+1})\cdot(\x-\x_{k+1})\]

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

\[\nabla f(\x) \approx \nabla f(\x_{k+1})+H_{k+1}\cdot(\x-\x_{k+1})\]

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

\[\g_{k+1} - \g_k \approx H_{k+1} \cdot (\x_{k+1} - \x_k)\]

引入新记号

\[\s_k=\x_{k+1}-\x_k, \quad \y_k=\g_{k+1}-\g_k\]

则有

\[\y_k \approx H_{k+1} \cdot \s_k, \qquad \s_k \approx H_{k+1}^{-1} \cdot \y_k\]

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

\[\y_k = B_{k+1} \cdot \s_k, \qquad \s_k = D_{k+1} \cdot \y_k\]

作为指导。

DFP

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

\[D_{k+1} = D_k+\Delta D_k, \quad k=0,1,2,\cdots\]

其中，$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完全类似。迭代式如下：

\[B_{k+1}=B_k+\Delta B_k, \quad k=0,1,2,\cdots\]

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_{k+1}^{-1}=\left(I-\frac{\s_k\y_k^\T }{\y_k^\T \s_k}\right) B_k^{-1} \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}\]

为了避免出现矩阵求逆符号，统一将$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步的迭代式

\[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}\]

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

\[D_{k+1}=V_k^\T D_kV_k + \rho_k\s_k\s_k^\T\]

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

\[\begin{eqnarray*} D_1 &=& V_0^\T D_0V_0+\rho_0\s_0\s_0^\T \\ D_2 &=& V_1^\T D_1V_1+\rho_1\s_1\s_1^\T \\ &=& V_1^\T (V_0^\T D_0V_0+\rho_0\s_0\s_0^\T )V_1+\rho_1\s_1\s_1^\T \\ &=& V_1^\T V_0^\T D_0V_0V_1+V_1^\T \rho_0\s_0\s_0^\T V_1 + \rho_1\s_1\s_1^\T \\ D_3 &=& V_2^\T D_2V_2+\rho_2\s_2\s_2^\T \\ &=& V_2^\T (V_1^\T V_0^\T D_0V_0V_1+V_1^\T \rho_0\s_0\s_0^\T V_1 + \rho_1\s_1\s_1^\T )V_2+\rho_2\s_2\s_2^\T \\ &=& V_2^\T V_1^\T V_0^\T D_0V_0V_1V_2+V_2^\T V_1^\T \rho_0\s_0\s_0^\T V_1V_2+V_2^\T \rho_1\s_1\s_1^\T V_2+\rho_2\s_2\s_2^\T \\ \cdots & & \end{eqnarray*}\]

一般的，我们有

\[\begin{eqnarray*} D_{k+1} &=& & (V_k^\T V_{k-1}^\T \cdots V_1^\T V_0^\T )D_0(V_0V_1\cdots V_{k-1}V_k) \\ &&+& (V_k^\T V_{k-1}^\T \cdots V_2^\T V_1^\T )(\rho_0\s_0\s_0^\T )(V_1V_2\cdots V_{k-1}V_k) \\ &&+& (V_k^\T V_{k-1}^\T \cdots V_3^\T V_2^\T )(\rho_1\s_1\s_1^\T )(V_2V_3\cdots V_{k-1}V_k) \\ &&+& \cdots \\ &&+& (V_k^\T V_{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}^\T )V_k \\ &&+& \rho_k\s_k\s_k^\T \end{eqnarray*}\]

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

\[\begin{eqnarray*} 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\s_0\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\s_1\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}\s_{k-2}\s_{k-2}^\T )(V_{k-1}V_k) \\ &&+& V_k^\T (\rho_{k-1}\s_{k-1}\s_{k-1}^\T )V_k \\ &&+& \rho_k\s_k\s_k^\T \end{eqnarray*}\]

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

参考

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