牛顿法和拟牛顿法

Posted on

牛顿法和拟牛顿法

牛顿法（Newton method）和拟牛顿法（quasi Newton method）是求解无约束最优化问题的常用方法，有收敛速度快的优点。

  • 牛顿法是迭代算法，每一步需要求解目标函数的hessian矩阵的逆矩阵，计算比较复杂；
  • 拟牛顿法通过正定矩阵近似hessian矩阵的逆矩阵或hessian矩阵，简化了这一计算过程。

1. 牛顿法

考虑无约束最优化问题 \[\min_{\pmb{x} \in \mathbf{R}^n} \quad f(\pmb{x})\] 其中\(\pmb{x}^*\)为目标函数的极小点。
假设\(f(\pmb{x})\)具有二阶连续偏导数，若第\(k\)次迭代值为\(\pmb{x}^{(k)}\)，则可将\(f(\pmb{x})\)在\(\pmb{x}^{(k)}\)附近进行二阶泰勒展开： \[f(\pmb{x}) = f(\pmb{x}^{(k)}) + g_k^T (\pmb{x} - \pmb{x}^{(k)}) + \frac{1}{2}(\pmb{x} - \pmb{x}^{(k)})^T H(\pmb{x}^{(k)})(\pmb{x} - \pmb{x}^{(k)}) \] 上式中\(g_k = g(\pmb{x}^{(k)}) = \nabla f(\pmb{x}^{(k)})\)是\(f(\pmb{x})\)的梯度向量在点\(\pmb{x}^{(k)}\)的值，\(H(\pmb{x}^{(k)})\)是\(f(\pmb{x})\)的hessian矩阵 \[H(\pmb{x}) = \left[ \frac{\partial^2 f}{\partial x_i \partial x_j} \right]_{n \times n}\] 在点\(\pmb{x}^{(k)}\)的值，函数\(f(\pmb{x})\)有极值的必要条件是在极值点处一阶导数为0，即梯度向量为0.特别是当\(H(\pmb{x}^{(k)})\)是正定矩阵时，函数\(f(\pmb{x})\)的极值为极小值。
牛顿法利用极小点的必要条件 \[\nabla f(\pmb{x}) = 0\] 每次迭代中从点\(\pmb{x}^{(k)}\)开始，求目标函数的极小点，作为第\(k+1\)次迭代值\(\pmb{x}^{(k+1)}\)。具体地，假如\(\pmb{x}^{(k+1)}\)满足： \[\nabla f(\pmb{x}^{(k+1)}) = 0\] 由二阶泰勒展开可得 \[\nabla f(\pmb{x}) = g_k + H_k (\pmb{x} - \pmb{x}^{(k)})\] 其中\(H_k = H(\pmb{x}^{(k)})\)，这样，\(\nabla f(\pmb{x}^{(k+1)}) = 0\)可写为 \[g_k + H_k(\pmb{x}^{(k+1)} - \pmb{x}^{(k)}) = 0\] 因此 \[\pmb{x}^{(k+1)} = \pmb{x}^{(k)} - H_k^{-1} g_k\] 或者 \[\pmb{x}^{(k+1)} = \pmb{x}^{(k)} + p_k\] 其中 \[H_k p_k = - g_k\]

算法如下：

在步骤（4）中需要求\(p_k\)，其中需要求\(H_k^{-1}\)，计算比较复杂，所以有其他改进的算法。

牛顿法和梯度下降法的比较：
梯度下降法迭代公式： \[\pmb{x}^{n+1} = \pmb{x}^n - \epsilon \nabla f(\pmb{x}^n)\] 牛顿法迭代公式： \[\pmb{x}^{n+1} = \pmb{x}^n - \frac{\nabla f(\pmb{x}^n)}{\nabla^2 f(\pmb{x}^n)}\] 与梯度下降法相比，牛顿法将学习率\(\epsilon\)替换为\((\nabla^2 f(\pmb{x}^n))^{-1}\)，可以理解为牛顿法根据代价函数的二阶导数信息，自动计算出合适的学习率，因此有更快的迭代速度。但对于hessian矩阵的计算，消耗资源大，速度慢，所以在实际应用中，牛顿法并不常用。

2. 拟牛顿法

在牛顿法的迭代中，需要计算hessian矩阵的逆矩阵\(H^{-1}\)，由于计算比较复杂，所以考虑用一个\(n\)阶矩阵\(G_k = G(\pmb{x}^{(k)})\)来近似替代\(H_k^{-1} = H^{-1}(\pmb{x}^{(k)})\)，这就是拟牛顿法的基本想法。

首先看牛顿法迭代中hessian矩阵\(H_k\)满足的条件。在式\(\nabla f(\pmb{x}) = g_k + H_k (\pmb{x} - \pmb{x}^{(k)})\)中取\(\pmb{x} = \pmb{x}^{(k+1)}\)，即得 \[g_{k+1} - g_k = H_k (\pmb{x}^{(k+1)} - \pmb{x}^{(k)})\] 即\(y_k = g_{k+1} - g_k, \quad \delta_k = \pmb{x}^{(k+1)} - \pmb{x}^{(k)}\)，则 \[y_k = H_k \delta_k\] 或 \[H_k^{-1} y_k = \delta_k\] 上述两式称为拟牛顿条件。

如果\(H_k\)是正定的（\(H_k^{-1}\)也是正定的），那么可以保证牛顿法搜索方向\(p_k\)是下降方向。这是因为搜索方向是\(p_k = -H_k^{-1} g_k\)，由式\(\pmb{x}^{(k+1)} = \pmb{x}^{(k)} - H_k^{-1} g_k\)，有 \[\pmb{x} = \pmb{x}^{(k)} + \lambda p_k = \pmb{x}^{(k)} - \lambda H_k^{-1} g_k\] 所以\(f(\pmb{x})\)在\(\pmb{x}^{(k)}\)的泰勒展开式可以近似写成： \[f(\pmb{x}) = f(\pmb{x}^{(k)}) - \lambda g_k^T H_k^{-1} g_k\] 因为\(H_k^{-1}\)正定，所以有\(g_k^T H_k^{-1} g_k > 0\)，当\(\lambda\)是一个充分小的正数时，总有\(f(\pmb{x}) < f(\pmb{x}^{(k)})\)，也就是说\(p_k\)是下降方向。
拟牛顿法将\(G_k\)作为\(H_k^{-1}\)的近似，要求矩阵\(G_k\)满足同样的条件。首先，每次迭代矩阵\(G_k\)是正定的。同时，\(G_k\)满足下面的拟牛顿条件： \[G_k y_k = \delta_k\] 按照拟牛顿条件选择\(G_k\)作为\(H_k^{-1}\)的近似或选择\(B_k\)作为\(H_k\)的近似的算法称为拟牛顿法。
按照拟牛顿条件，在每次迭代中可以选择更新矩阵\(G_{k+1}\): \[G_{k+1} = G_k + \Delta G_k\] 这种选择有一定的灵活性，因此有多种具体实现方法。

2.1. DFP（Davidon-Fletcher-Powell）算法

DFP算法选择\(G_{k+1}\)的方法是，假设每一步迭代中矩阵\(G_{k+1}\)是由\(G_k\)加上两个附加项构成的，即 \[G_{k+1} = G_k + P_k + Q_k\] 其中\(P_k\)和\(Q_k\)是待定矩阵，这时 \[G_{k+1} y_k = G_k y_k + P_k y_k + Q_k y_k\] 为使\(G_{k+1}\)满足拟牛顿条件，可使\(P_k\)和\(Q_k\)满足： \[\begin{equation} P_k y_k = \delta_k \nonumber \\ Q_k y_k = - G_k y_k \nonumber \end{equation}\] 事实上找出\(P_k\)和\(Q_k\)并不困难，例如取 \[\begin{equation} P_k = \frac{\delta_k \delta_k^T}{\delta_k^T y_k} \nonumber \\ Q_k = - \frac{G_k y_k y_k^T G_k}{y_k^T G_k y_k} \nonumber \end{equation}\] 这样就可以得到矩阵\(G_{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} \tag{B.23}\] 称为DFP算法。
可以证明，如果初始矩阵\(G_0\)是正定的，则迭代过程中的每个矩阵\(G_k\)都是正定的。

算法如下：

2.2. BFGS（Broyden-Fletcher-Goldfarb-Shanno）算法

除了用\(G_k\)逼近hessian矩阵的逆矩阵\(H_k^{-1}\)，也可以用\(B_k\)逼近hessian矩阵\(H_k\)，相应的拟牛顿条件为 \[B_{k+1} \delta_k = y_k\] 可以用同样的方法得到另一条迭代公式，令 \[\begin{equation} B_{k+1} = B_k + P_k + Q_k \nonumber \\ B_{k+1} \delta_k = B_k \delta_k + P_k \delta_k + Q_k \delta_k \nonumber \end{equation}\] 使\(P_k\)和\(Q_k\)满足： \[\begin{equation} P_k \delta_k = y_k \nonumber \\ Q_k \delta_k = - B_k \delta_k \nonumber \end{equation}\] 找到合适条件的\(P_k\)和\(Q_k\)，得到BFGS算法矩阵\(B_{k+1}\)的迭代公式： \[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} \tag{B.30}\] 可以证明，如果初始矩阵\(B_0\)是正定的，则迭代过程中的每个矩阵\(B_k\)都是正定的。

算法如下：

2.3. Broyden类算法

Sherman-Morrison公式：假设\(A\)是\(n\)阶可逆矩阵，\(u,v\)是\(n\)维向量，且\(A+uv^T\)也是可逆矩阵，则 \[(A+uv^T)^{-1} = A^{-1} - \frac{A^{-1} u v^T A^{-1}}{1+v^T A^{-1} u}\]

可以从BFGS算法矩阵\(B_k\)的迭代公式（B.30）得到BFGS算法关于\(G_k\)的迭代公式。事实上，若记\(G_k = B_k^{-1}\)，\(G_{k+1} = B_{k+1}^{-1}\)，那么对式（B.30）两次应用Sherman-Morrison公式即得 \[G_{k+1} = \left( I - \frac{\delta_k y_k^T}{\delta_k^T y_k} \right) G_k \left( I - \frac{\delta_k y_k^T}{\delta_k^T y_k} \right)^T + \frac{\delta_k \delta_k^T}{\delta_k^T y_k} \tag{B.31}\] 称为BFGS算法关于\(G_k\)的迭代公式。
有DFP算法\(G_k\)的迭代公式（B.23）得到的\(G_{k+1}\)记作\(G^{DFP}\)，有BFGS算法\(G_k\)的迭代公式（B.31）得到的\(G_{k+1}\)记作\(G^{BFGS}\)，它们均满足方程拟牛顿条件式，所以它们的线性组合 \[G_{k+1} = \alpha G^{DFP} + (1 - \alpha) G^{BFGS}\] 也满足拟牛顿条件式，而且是正定的。其中\(0 \leqslant \alpha \leqslant 1\)。这样就得到一类拟牛顿法，称为Broyden算法。