牛顿法（Newton Methods）、阻尼牛顿法和拟牛顿法

令 $X=(x_1,x_2,\cdots ,x_N{)}^T\in {\mathbf{R}}^N$，目标函数 $f:{\mathbf{R}}^N\to \mathbf{R}$， $f$ 为凸函数，且二阶连续可微，我们希望求解如下的无约束极小化问题：

$$\underset{X}{\min }f(X)$$

1 牛顿法

1.1 $N=1$ 时的迭代公式

为了简单起见，这里先考虑 $N=1$ 的情形，此时目标函数 $f(X)$ 变为 $f(x)$。

牛顿法的基本思想是： 在现有极小点估计值得附近对 $f(x)$ 做二阶泰勒展开，进而找到极小点的下一个估计值。假设 $x_k$ 是当前的极小点估计值，则：

$$\phi (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)$ 在 $x_k$ 附近的二阶泰勒展开式（其中略去了关于 $x-x_k$ 的高阶项）。因为我们的目标是求最值，由极值的必要条件可知，$\phi (x)$ 应该满足：

$${\phi }^{\prime }(x)=f^{\prime }(x_k)+f^{\prime \prime }(x_k)(x-x_k)=0$$

从而有：

$$x=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$$

于是，若给定初始值 $x_0$，则可以按照下面的迭代公式

$$x_{k+1}=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)},k=0,1,\cdots$$

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

在一定的条件下，{xk}\{x_k\}{xk​} 可以收敛到 $f(x)$ 的极小值点。

1.2 $N>1$ 时的迭代公式

当 $N>1$ 时，二阶泰勒展示式写作：

$$\phi (X)=f(X_k)+\nabla f(X_k)\cdot (X-X_k)+\frac{1}{2}\cdot (X-X_k{)}^T\cdot {\nabla }^{2}f(X_k)\cdot (X-X_k)$$

其中，$\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\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{array}\right]}_{N\times N}$$

$\nabla f$ 与 ${\nabla }^{2}f$ 中的函数均为关于 $X$ 的函数，分别记为 $g$ 和 $H$（ $g$ 表示 gradient，$H$ 表示 Hessian）。特别地，若 $f$ 的混合偏导数可交换次序（即对 $\forall 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(X_k)$ 与 ${\nabla }^{2}f(X_k)$ 表示将 $X$ 取为 $X_k$ 后得到的实值向量和矩阵，以下分别记为 $g_k$ 和 $H_k$。

同样，因为是求极小点，极值的必要条件要求它是 $\phi (X)$ 的 驻点，即：

$$\nabla \phi (X)=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,k=0,1,\cdots \text{\hspace{1em}(1-1)}$$

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

1.3 完整的算法描述

算法的伪代码描述为：

1. 给定初始值 $X_0$ 和精度阈值 $\epsilon$，并令 $k:=0$
2. 计算 $g_k$ 和 $H_k$
3. 若 $\Vert g_k\Vert <\epsilon$，则停止迭代；否则计算搜索方向 $d_k=-H_k^{-1}\cdot g_k$
4. 计算新的迭代点 $X_{k+1}:=X_k+d_k$
5. 令 $k:=k+1$，跳到第2步

优点

当目标函数是二次函数时，由于二次泰勒展开函数是与原目标函数完全相同的二次式，海森矩阵退化成一个常数矩阵，从任一初始点出发，利用公式（1-1）只需一步迭代就可以达到 $f(x)$ 的极小点 $x^{\ast }$，因此牛顿法是一种具有 二次收敛性 的算法。

对于非二次函数，若函数的二次性态较强，或迭代点已进入极小点的邻域，则其收敛速度也是很快的，这是牛顿法的主要优点。

缺点

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

2 阻尼牛顿法

为了消除原始牛顿法的缺点，人们提出了阻尼牛顿法。

对于牛顿法，确定了迭代方向之后，迭代步长默认为1，但是这个迭代方向并不一定是朝着函数值下降的方向。阻尼牛顿法每次的迭代方向仍采用 $d_k$，但每次迭代需沿此方向做一维搜索（line search），寻找最优的步长因子 ${\lambda }_k$，即：

$${\lambda }_k=arg\underset{\lambda \in R}{\min }f(X_k+\lambda d_k)\text{\hspace{1em}(2-1)}$$

算法的伪代码描述为：

1. 给定初始值 $X_0$ 和精度阈值 $\epsilon$，并令 $k:=0$
2. 计算 $g_k$ 和 $H_k$
3. 若 $\Vert g_k\Vert <\epsilon$，则停止迭代；否则计算搜索方向 $d_k=-H_k^{-1}\cdot g_k$
4. 利用公式（2-1）得到步长${\lambda }_k$，并计算新的迭代点 $X_{k+1}:=X_k+{\lambda }_k{d}_k$
5. 令 $k:=k+1$，跳到第2步

可以看到，阻尼牛顿法相比于牛顿法，在每次参数更新之前，利用一维搜索法计算更新步长，确保优化方向为下降方向。

3 拟牛顿法

原始牛顿法虽然收敛速度快，但是需要计算海森矩阵的逆矩阵 $H^{-1}$，而且有时目标函数的海森矩阵无法保持正定，从而使得原始牛顿法失效。为了克服这两个问题，人们提出了拟牛顿法。这个方法的 基本思想 是：不用二阶偏导数，而是构造出一个可以近似海森矩阵（或海森矩阵的逆）的正定对称阵。不同的构造方法就产生了不同的拟牛顿法。

下面我们先推导一下拟牛顿条件，它给“对海森矩阵（或海森矩阵的逆）做近似”提供了理论指导，指出了用来近似的矩阵应该满足的条件。

3.1 拟牛顿条件

对 $\nabla f(x)$ 在 $x_k$ 做泰勒展开我们可以得到以下近似：

$$\nabla f(x)=g_k+H_k(x-x_k)$$

取 $x=x_{k+1}$，则有：

$$\nabla f(x_{k+1})=g_{k+1}=g_k+H_k(x_{k+1}-x_k)$$

即：

$$g_{k+1}-g_k=H_k(x_{k+1}-x_k)$$

记 $y_k=g_{k+1}-g_k$，${\delta }_k=x_{k+1}-x_k$，则有：

$$y_k=H_k{\delta }_k\text{\hspace{1em}(3-1)}$$

或

$$H_k^{-1}{y}_k={\delta }_k\text{\hspace{1em}(3-2)}$$

以上即为拟牛顿条件。

常用的拟牛顿法有DFP、BFGS、L-BFGS，区别在于如何选取替代矩阵。

3.2 DFP算法

GFP算法用于近似拟牛顿条件（3-2），这里用 $G_k$ 代表对 $H_k^{-1}$ 的近似，下面直接给出计算公式：

$$G_{k+1}=G_k+\frac{{\delta }_k{\delta }_k^T}{{\delta }_k^T{y}_k}-\frac{G_k{y}_k{y}_k^T{G}_k}{y_k^T{G}_k{y}_k}$$

可以证明，如果初始矩阵 $G_0$ 是正定对称的，则迭代过程中的每个矩阵 $G_k$ 都是正定对称的，一般取 $G_0=I$。

3.3 BFGS算法

BFGS 算法用于近似拟牛顿条件（3-1），与 DFP 相比，BFGS 的性能更佳。这里用 $B_k$ 代表对 $H_k$ 的近似，下面直接给出计算公式：

$$B_{k+1}=B_k+\frac{y_k{y}_k^T}{y_k^T{\delta }_k}-\frac{B_k{\delta }_k{\delta }_k^T{B}_k}{{\delta }_k^T{B}_k{\delta }_k}$$

可以证明，如果初始矩阵 $B_0$ 是正定对称的，则迭代过程中的每个矩阵 $B_k$ 都是正定对称的，一般取 $B_0=I$。

若记 $G_k=B_k^{-1}$， $G_{k+1}=B_{k+1}^{-1}$ ，那么应用Sherman-Morrison公式可以将上述迭代公式改写为：

$$G_{k+1}=(I-\frac{{\delta }_k{y}_k^T}{{\delta }_k^T{y}_k})G_k(I-\frac{{\delta }_k{y}_k^T}{{\delta }_k^T{y}_k}{)}^T+\frac{{\delta }_k{\delta }_k^T}{{\delta }_k^T{y}_k}$$

上式就是 BFGS 算法关于 $G_k$ 的迭代公式。

3.4 L-BFGS算法

在BFGS中，需要用到一个 $N$ 阶矩阵 $G_k$ ，当 $N$ 很大时，存储这个矩阵将消耗大量的计算机资源。为了解决这个问题，减少 BFGS 迭代过程中所需的内存开销，就有了 L-BFGS。

L-BFGS（Limited-memory BFGS 或 Limited-storage BFGS）对 BFGS 进行了近似，其 基本思想 是：不再存储完整的矩阵 $G_k$ ，而是存储计算过程中的向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk​}{yk​}，需要矩阵 $G_k$ 时，利用向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk​}{yk​} 的计算来代替。而且，向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk​}{yk​} 也不是所有的都存储，而是保留最新的 $m$ 个，每次计算 $G_k$ 时，只利用最新的 $m$ 个向量序列 {δk}{yk}\{\delta_k\}\{y_k\}{δk​}{yk​}。这样一来，存储空间由 $O(N^2)$ 降至 $O(mN)$。