优化算法 - 牛顿法 and 拟牛顿法

预备知识

• 无约束优化问题

  我们以前聊过约束优化方法，将带有约束的优化问题（例如 SVM 中间隔的间隔要大于 1）通过 拉格朗日乘子法转化为凸优化问题，再对各参数求偏导从而求的最优解。

  那么对于无约束的优化问题呢？我们通常使用跟梯度有关的算法。例如 梯度下降法，在 Logistic回归 中以及神经网络中都有应用。而关于梯度的无约束优化算法还有牛顿法系列。

  在介绍牛顿法前我们先来了解两个基础知识，Hesse 矩阵以及 泰勒展开式。

• Hesse 矩阵

  简单来说，Hesse 矩阵就是某一矩阵对于其中各元素的偏导（二阶）。

  例如有矩阵：$A=\left[\begin{array}{cc}{x}_1{x}_1& x_1{x}_2\\ x_2{x}_1& x_2{x}_2\end{array}\right]$

  则，$H(x)=\left[\begin{array}{cc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}\\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}\end{array}\right]$

  Hesse矩阵的通用表达式为：$H(x)=[\frac{{\partial }^{2}f}{\partial x_i\partial x_j}{]}_{n\ast n}$

• 泰勒展开式

  高数中的经典公式，将一个函数，在某一点处，利用各阶导数计算原函数值。这里以二阶展开举例，省略高阶无穷小的余项。

  将 $f(x)$ 在 x₀ 处的二阶展开：

  $f(x)=f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{1}{2}{f}^{\prime \prime }(x_0)(x-x_0{)}^2$

  对于 x 为矩阵来讲，令 $g(x)$ 表示一阶导，$H(x)$ 表示二阶导：

  $\begin{array}{rl}f(x)=& f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{1}{2}{f}^{\prime \prime }(x_0)(x-x_0{)}^2\\ =& f(x_0)+g^T(x_0)(x-x_0)+\frac{1}{2}(x-x_0{)}^{T}H(x_0)(x-x_0)\end{array}$

牛顿法

接下来进入正题，对于无约束优化问题来讲，无法直接求解，所以需要类似梯度下降一样一步一步的迭代。牛顿法也同样基于这个理念。

• 形象化解释

  首先来看一个形象化的解释：

  $\bullet$ 目标：求 $f(x)=0$ 的解 $x^{\ast }$

  $\bullet$ 初始选择一个接近函数 $f(x)$ 零点的 $x_0$，计算相应的 $f(x_0)$ 和切线斜率 $f^{\prime }(x_0)$。

  $\bullet$ 利用“两点式”计算经过点 $(x_0,f(x_0))$ 且斜率为 $f^{\prime }(x_0)$ 的直线，与 $x$ 轴的交点 $(x_1,0)$.

  $f^{\prime }(x_0)=\frac{f(x_0)-f(x_1)}{x_0-x_1}=\frac{f(x_0)}{x_0-x_1}$

  $x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$，$x_1$ 则为下一次迭代的点，

  牛顿法又称为切线法，如图所示：
  （图片来自 https://www.cnblogs.com/shixiangwan/p/7532830.html）

  以上步骤是求解能使 $f(x)=0$ 时的解 $x^{\ast }$，但我们现在要求的是能使 $f^{\prime }(x)=0$ 的解，所以，我们的目标变为：

  $f^{\prime \prime }(x_0)=\frac{f^{\prime }(x_0)-f^{\prime }(x_1)}{x_0-x_1}$

  得到：

  $x_1=x_0-\frac{f^{\prime }(x_0)}{f^{\prime \prime }(x_0)}{}_{⟹}^k{x}_{k+1}=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$

  这便是我们牛顿法的迭代更新过程。

• 公式化解释

  对于函数 $g(x)$，其导数 $H(x)$，对 $g(x)$ 在点 $x_k$ 出使用泰勒二阶展开，得：

  $g(x)=g(x_k)+H(x_k)(x-x_k)$

  令 $x=x_{k+1}，g_k=g(x_k)，H_k=H(x_k)$，则更新迭代公式为：

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

  （$g_{k+1}=0时则与之前所推等价$）

• 算法过程

  1. $确定精度要求\epsilon ，初始化x_0，置k=0$
  2. $计算g_k=g(x_k)，若||g(x_k)||\le \epsilon 则停止计算，求得近似解x^{\ast }=x_k；否则转第3步$
  3. $计算H_k=H(x_k)，以及其逆矩阵H_k^{-1}$
  4. $计算x_{k+1}=x_k-H_k^{-1}{g}_k，令k=k+1，转第2步$

  我们可以看到，在迭代过程中，需要求 $H_k的逆矩阵H_k^{-1}$，这个计算比较复杂，效率很低，所以就产生了拟牛顿法。

拟牛顿法

拟牛顿法的思想是与牛顿法一样的，但由于求逆矩阵复杂，多以拟牛顿法寻找一个矩阵用来近似替代 $H_k^{-1}$

• 拟牛顿条件

  想找替代矩阵可不是随便找的，需要满足一个条件（两点式），即：

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

  此条件被称为“拟牛顿条件”。

  ---

  拟牛顿法则是使用 $G_k$ 作为 $H_k^{-1}$ 的近似，要求 $G_k$ 为正定矩阵，且满足拟牛顿条件。根据拟牛顿条件的不同，选取的正定矩阵也不一样：

  如果选择 $G_k$ 作为 $H_k^{-1}$ 的近似，就成为 DFP 算法；
  如果选择 $B_k$ 作为 $H_k$ 的近似，就成为 BFGS 算法。

  按照拟牛顿条件，在每次迭代中可以选择更新替代矩阵，用以近似 $H_k^{-1}或H_k$：

  $G_{k+1}=G_k+\Delta G_k$

  ---

  为什么要求 $G_k或B_k$ 为正定矩阵呢？那是因为只有它们为正定矩阵时，x 的搜索方向才是 $f(x)$ 下降的方向。

  具体的，根据迭代更新公式：$x_{k+1}=x_k-H_k^{-1}{g}_k$，

  引入步长因子 $\lambda$，定义 $p_k=-H_k^{-1}{g}_k$

  迭代更新公式变为：$x_{k+1}=x_k-\lambda P_k$

  ---

  $f(x_{k+1})$ 在 $x_k$ 处的泰勒一阶展开：$f(x_{k+1})=f(x_k)+g_k^T(x_{k+1}-x_k)$

  将迭代更新公式代入，得：

  $f(x_{k+1})=f(x_k)-\lambda g_k^T{H}_k^{-1}{g}_k$

  因为 $H_k$ 是正定的（$H_k^{-1}$ 也是正定的），所以有 $g_k^T{H}_k^{-1}{g}_k>0$，当 $\lambda$ 为充分小的正数时，总有 $f(x_{k+1})<f(x_k)$，就意味着 $f(x)$ 在迭代中逐步逼近最小值，也就是说 $p_k$ 的方向是函数下降方向。

  （原牛顿法中没有步长因子 $\lambda$ 的存在，当 $G_k或B_k$ 为奇异矩阵时可能会出现数值问题，即出现 $f(x_{k+1})>f(x_k)$ 的情况。此时目标函数值并没有稳定下降，而且有可能无法收敛。）

• DFP (Davidon-Fletcher-Powell)

  • 算法推导

    DFP 算法选择 $G_k$ 作为 $H_k^{-1}$ 的近似，假设每一步迭代中矩阵 $G_{k+1}$ 是由 $G_k$ 加上两个附加矩阵构成的，即：

    $G_{k+1}=G_k+P_k+Q_k$，其中 $P_k，Q_k$ 为待定矩阵，令 $y_k=g_{k+1}-g_k$，${\delta }_k=x_{k+1}-x_k$

    这时 ，$G_{k+1}{y}_k=G_k{y}_k+P_k+Q_k{y}_k$

    为使 $G_{k+1}$ 满足拟牛顿条件，可使 $P_k，Q_k$ 满足：

    $P_k{y}_k={\delta }_k$，可以找到 $P_k=\frac{{\delta }_k{\delta }_k^T}{{\delta }_k^T{y}_k}$

    $Q_k{y}_k=-G_k{y}_k$，可以找到 $Q_k=-\frac{G_k{y}_k{y}_k^T{G}_k}{y_k^T{G}_k{y}_k}$

    由以上得到矩阵 $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}$ （公式 G）

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

  • 算法过程

    1. $确定精度要求\epsilon ，初始化x_0、G_0为正定对称矩阵，置k=0$
    2. $计算g_k=g(x_k)，若||g(x_k)||\le \epsilon 则停止计算，求得近似解x^{\ast }=x_k；否则转第3步$
    3. $计算p_k=G_k{g}_k$
    4. $一维搜索：{\lambda }_k={}_{\lambda \ge 0}^{argmin}f(x_k+\lambda p_k)$
    5. $计算x_{k+1}=x_k+{\lambda }_k{p}_k$
    6. $计算g_{k+1}=g(x_{k+1})，若||g(x_k)||\le \epsilon 则停止计算，求得近似解x^{\ast }=x_k；否则，按（公式G）计算G_{k+1}，置k=k+1$

• BFGS (Broyden-Fletcher-Goldfarb-Shanno)

  BFGS 与 DFP 类似，不同的是选取的替代矩阵不同：

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

总结

对于牛顿法与拟牛顿法的计算过程的区别：

• 牛顿法：

  通过泰勒展开式近似原函数，每次迭代以计算 Hesse 的逆矩阵为核心，从而逼近 原函数的最小值。

• 拟牛顿法：

  通过初始化 $(G_0或B_0)$ ，一维搜索 ${\lambda }_k$，以及替代矩阵$(G_{k+1}或B_{k+1})$的迭代计算，代替了 Hesse 逆矩阵的计算。

  ---

  从本质上去看，牛顿法是二阶收敛，梯度下降是一阶收敛，所以牛顿法的收敛速度更快。更通俗地说，如果想找一条最短的路径走到一个山的最底部，梯度下降法每次只从当前所处位置选一个坡度最大的方向走一步，而牛顿法在选择方向时，不仅会考虑坡度是否够大，还会考虑在走了一步之后，坡度是否会变得更大。所以，可以说牛顿法比梯度下降法看得更远一点，能更快地走到最底部。（牛顿法目光更加长远，所以少走弯路；相对而言，梯度下降法只考虑了局部的最优，没有全局思想。）