梯度类算法原理：最速下降法、牛顿法和拟牛顿法

算法结构

梯度类算法有很多，本文主要学习最常见的3个算法：最速下降法、牛顿法和拟牛顿法。算法名称虽多，但是他们的算法结构都是一样的，可以描述为

（1）选定初始点$x_0$。

（2）迭代公式为$x_{k+1}=x_k+{\alpha }_k{d}_k$，其中${\alpha }_k$和$d_k$分别为第$k$次的迭代步长和迭代方向。

人们常说，选择大于努力。对优化算法来说，也是类似的：如果$d_k$设计的不好，可能无论${\alpha }_k$怎么努力，都得不到最优解。所以，我们应该重点关注$d_k$的构造方法。

最速下降法

最朴素的想法，就是如果针对任意初始点，都能直接算出让函数下降最快的方向，就好了。下面我们看看怎么能够做到。

假设$x_k=[x_{k1},x_{k2},..,x_{kn}]$，给定方向$s$上的一个增量为
$$\Delta s=[\Delta x_1,\Delta x_2,...,\Delta x_n]$$
定义在点$x_k$处沿方向$s$的变化率为
$$\frac{\partial f(x)}{\partial s}{|}_{x_k}=\underset{\Delta s\to 0}{\lim }\frac{f(x_k+\Delta s)-f(x_k)}{\Vert \Delta s\Vert }$$
当$\Delta s$足够小时，上式右侧的分子变为全微分形式
$$f(x_k+\Delta s)-f(x_k)=\frac{\partial f}{\partial x_{k1}}{|}_{x_k}\Delta x_1+\frac{\partial f}{\partial x_{k2}}{|}_{x_k}\Delta x_2+\cdot \cdot \cdot +\frac{\partial f}{\partial x_{kn}}{|}_{x_k}\Delta x_n$$
将上式带入上上式
$$\begin{gathered} \frac{\partial f(x)}{\partial s}{|}_{x_k}=\underset{\Delta s\to 0}{\lim }(\frac{\partial f}{\partial x_{k1}}{|}_{x_k}\frac{\Delta x_1}{\Delta s}+\frac{\partial f}{\partial x_{k2}}{|}_{x_k}\frac{\Delta x_2}{\Delta s}+\cdot \cdot \cdot +\frac{\partial f}{\partial x_{kn}}{|}_{x_k}\frac{\Delta x_n}{\Delta s}) \\ =\frac{\partial f}{\partial x_1}{|}_{x_k}\cos {\alpha }_1+\frac{\partial f}{\partial x_2}{|}_{x_k}\cos {\alpha }_2+\cdot \cdot \cdot +\frac{\partial f}{\partial x_n}{|}_{x_k}\cos {\alpha }_n \\ =[\nabla f(x_k){]}^{T}e \end{gathered}$$
此处，$\cos {\alpha }_i=\underset{\Delta s\to 0}{\lim }\frac{\Delta x_i}{\Delta s}(i=1,2,\cdot \cdot \cdot ,n)$，表示方向$s$与坐标轴$x_i夹$角的余弦值。$[\nabla f(x_k)]$是函数$f$在$x_k$处的梯度向量，表达式为
$$[\nabla f(x_k){]}^T=[\frac{\partial f}{\partial x_{k1}}{|}_{x_k},\frac{\partial f}{\partial x_{k2}}{|}_{x_k},\cdot \cdot \cdot ,\frac{\partial f}{\partial x_{kn}}{|}_{x_k}]$$
$e=[\cos {\alpha }_1,\cos {\alpha }_2,\cdot \cdot \cdot ,\cos {\alpha }_n]$，且$\Vert e\Vert =1$，所以$e$是$s$方向上的单位向量。

在$x_k$处，由于梯度向量是固定的，而$e$的模值为1，所以变化率的值仅随梯度向量和$e$之间的相对空间关系变化。当两者为同一方向时，变化率为最大正值；当两者为相反方向时，变化率为最小负值。

在最速下降法中，就设定
$$d_k=-[\nabla f(x_k)]$$

牛顿法

看起来，最速下降法已经非常优秀了，毕竟都已经沿着让函数下降最快的方向了，还能更好的设计方案嘛？

答案是有的。单从下一步来看，负梯度方向确实是最好的；但是如果看两步呢？两次负梯度方向的矢量和能否通过一次计算得到？甚至是，出现比两次负梯度方向矢量和更优的方向？

下面看一下牛顿法是如何解决以上问题的。

要看两步，就需要用到$x$的二次项。针对任意函数$f(x)$，使用泰勒公式展开，并保留二次项

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

其中，$H(x_k)$为海森矩阵。

上式的一阶导数为
$$\nabla f(x)=[\nabla f(x_k)]+H(x_k)[x-x_k]$$
令导数值为0，并写成迭代形式，得到
$$x_{k+1}=x_k-[H(x_k){]}^{-1}[\nabla f(x_k)]$$
以上即为牛顿法。相比最速下降法，牛顿法在$[\nabla f(x_k)]$的前面，多了一项$[H(x_k){]}^{-1}$。

牛顿法因为利用了泰勒展开式，所以当$x_k$接近（局部）最优解时，收敛速度非常快。

拟牛顿法

既然牛顿法相比最速下降法，可以理解为多看了一步，那是不是照葫芦画瓢，再多看几步，然后不断更新迭代方向呢？原则上是可以的，但是牛顿法已经有自己的问题了：（1）函数$f(x)$必须二阶可导；（2）$H$和$H^{-1}$的计算过程较为复杂。再多看几步，对函数$f(x)$的要求只会更高，所以更靠谱的优化方案是去降低$H^{-1}$的计算复杂度。

接下来要介绍的拟牛顿法，其主要思路是通过构造复杂度低的函数来替代$H^{-1}$，从而降低复杂度。

把牛顿法推导过程中用到的一阶导数公式再抄一遍，并写成迭代形式

$$\nabla f(x_{k+1})=[\nabla f(x_k)]+H(x_k)[x_{k+1}-x_k]$$
令$y_k=[\nabla f(x_{k+1})]-[\nabla f(x_k)]$，$s_k=x_{k+1}-x_k$，上式可以化简为
$$s_k=[H(x_k){]}^{-1}{y}_k$$

令$[H(x_k){]}^{-1}=G(x_k)$，并构造迭代公式
$$G_{k+1}=G_k+E_k$$
如果设定$G_0$为单位阵，只需要确定$E_k$，便可以替代海森矩阵了。

令$E_k=\alpha {\mu }_k{\mu }_k^T+\beta {\nu }_k{\nu }_k^T$，其中${\mu }_k$和${\nu }_k$均为$n\times 1$向量(构造成两项的原因会在后面描述)。
$$s_k=(G_k+\alpha {\mu }_k{\mu }_k^T+\beta {\nu }_k{\nu }_k^T)y_k$$
上式做一下变换
$$\alpha ({\mu }_k^T{y}_k){\mu }_k+\beta ({\nu }_k^T{y}_k){\nu }_k=s_k-G_k{y}_k$$
其中，${\mu }_k^T{y}_k$和${\nu }_k^T{y}_k$为实数，$s_k-G_k$为$n\times 1$向量，$\alpha$和$\beta$可以任意选取，我们取特殊的一组：${\mu }_k=r{G}_k{y}_k$，${\nu }_k=\theta s_k$，此时$E_k$变为
$$E_k=\alpha r^2{G}_k{y}_k{y}_k^T{G}_k^T+\beta {\theta }^2{s}_k{s}_k^T$$

将${\mu }_k$和${\nu }_k$带入上上式
$$\alpha r^2(y_k^T{G}_k^T{y}_k)G_k{y}_k+\beta {\theta }^2(s_k^T{y}_k)s_k=s_k-G_k{y}_k$$
化简一下
$$[\alpha r^2(y_k^T{G}_k^T{y}_k)G_k+1]G_k{y}_k+[\beta {\theta }^2{s}_k^T{y}_k-1]s_k=0$$
由于$G_k{y}_k$和$s_k$是任意的，所以需要$\alpha r^2(y_k^T{G}_k^T{y}_k)G_k+1=0$和$\beta {\theta }^2{s}_k^T{y}_k-1=0$，得到
$$\alpha r^2=-\frac{1}{y_k^T{G}_k^T{y}_k},\beta {\theta }^2=\frac{1}{s_k^T{y}_k}$$
从上上式可以看出，如果$E_k$在构造时只有一项，是无法保证恒等式的。

所以，最终的拟牛顿法迭代公式为
$$G_{k+1}=G_k-\frac{G_k{y}_k{y}_k^T{G}_k^T}{y_k^T{G}_k^T{y}_k}+\frac{s_k{s}_k^T}{s_k^T{y}_k}$$
相比$H^{-1}$，上式的计算不仅复杂度低，而且不需要$f(x)$为二阶可导。