李航统计学习方法-附录B 牛顿法与拟牛顿法

牛顿法(Newton method)和拟牛顿法(quasi-Newton method)也是求解无约束最优化的常用方法，有收敛速度快的优点。牛顿法是迭代算法，每一步需要求解目标函数的黑塞矩阵的逆矩阵，计算比较复杂。拟牛顿法通过正定矩阵近似黑塞矩阵的逆矩阵或黑塞矩阵，简化了这一计算过程。

1. 牛顿法

   考虑无约束最优化问题
   $$\underset{x\in R^n}{\min }f(x)\text{\hspace{1em}(B.1)}$$

其中$x^{\ast }$为目标函数的极小值。

​ 假设f(x)具有二阶连续偏导数，若第k次迭代值为$x^{(k)}$,则可将f(x)在$x^{(k)}$附近进行二阶泰勒展开：
$$f(x)=f(x^{(k)})+g_k^T(x-x^{(k)})+\frac{1}{2}(x-x^{(k)}{)}^{T}H(x^{(k)})(x-x^{(k)})\text{\hspace{1em}(B.2)}$$
这里，$g_k=g(x^{(k)})=\nabla f(x^{(k)})$是f(x)的梯度向量在点$x^{(k)}$的值，$H(x^{(k)})$是f(x)的黑塞矩阵（Hessian matrix）
$$H(x)=[\frac{{\partial }^{2}f}{\partial x_i\partial x_j}{]}_{n\times n}\text{\hspace{1em}(B.3)}$$
在点$(x^{(k)}$的值，函数f(x)有极值的必要条件是在极值点初一阶导数为0，即梯度向量为0。特别是当$H((x^{(k)})$是正定矩阵时，函数f(x)的极值为极小值。

​ 牛顿法利用极小值点的必要条件
$$\nabla f(x)=0\text{\hspace{1em}(B.4)}$$
每次迭代中从点$(x^{(k)})$开始，求目标函数的极小点，作为第k+1次迭代值$x^{(k+1)}$。具体地，假设$x^{(k+1)}$满足
$$\nabla f(x^{(k+1)})=0\text{\hspace{1em}(B.5)}$$
由式(B.2)，对$\nabla f(x)$在$x^{(k)}$处进行一阶泰勒展开。
$$\nabla f(x)=g_k+H_k(x-x^{(k)})\text{\hspace{1em}(B.6)}$$
其中$H_k=H(x^{(k)})$。这样，式(B.5)成为
$$g_k+H_k(x^{(k+1)}-x^{(k)})=0\text{\hspace{1em}(B.7)}$$
因此，
$$x^{(k+1)}=x^{(k)}-H_k^{-1}{g}_k\text{\hspace{1em}(B.8)}$$
或者
$$x^{(k+1)}=x^{(k)}+p_k\text{\hspace{1em}(B.9)}$$
其中，
$$H_k{p}_k=-g_k$$
用式(B.8)作为迭代公式的算法就是牛顿法。

算法B.1(牛顿法)

输入：目标函数f(x)，梯度g(x)=$\nabla f(x)$,黑塞矩阵H(x),精度要求$ϵ$;

输出：f(x)的极小值$x^{\ast }$。

(1)取初始点$x^{(}0)$，置k=0。

(2)计算$g_k=g(x^{(k)})$。

(3)若$||g_k||<ϵ$,则停止计算，得近似解$x^{\ast }=x^{(k)}$。

(4)计算$H_k=H(x^{(k)})$,并求$p_k$
$$H_k{p}_k=-g_k$$
(5)置$x^{(k+1)}=x^{(k)}+p_k$。

(6)置k=k+1,转(2)。

步骤(4)求$p_k,p_k=-H_k^{(-1)}{g}_k$,要求$H_k^{(-1)}$,计算比较复杂，所以有其他改进得方法。。

2. 拟牛顿法得思路

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

   先看牛顿法迭代中黑塞矩阵$H_k$满足的条件。首先，$H_k$满足以下关系。在式(B.6)中取$x=x^{(k+1)}$,即得
   $$g_{(k+1)}-g_k=H_k(x^{(k+1)}-x^{(k)})\text{\hspace{1em}(B.11)}$$
   记$y_k=g_{k+1}-g_k,{\delta }_k=x^{(k+1)}-x^{(k)},$则
   $$y_k=H_k{\delta }_k\text{\hspace{1em}(B.12)}$$
   或
   $$H_k^{-1}{y}_k={\delta }_k\text{\hspace{1em}(B.13)}$$
   式(B.12)或式(B.13)称为拟牛顿条件。

   ​ 如果$H_k$是正定的($H_k^{-1}$也是正定的)，那么可以保证牛顿法搜索方向$p_k$是下降方向。这是因为搜索方向是$p_k=-H_k^{-1}{g}_k$,由式(B.8)有
   $$x=x^{(k)}+\lambda p_k=x^{(k)}-\lambda H_k^{-1}{g}_k\text{\hspace{1em}(B.14)}$$
   所以f(x)在$x^{(k)}$的泰勒展开式(B.2)可以去掉2阶项可以近似写成：
   $$f(x)=f(x^{(k)}-\lambda g_k^T{H}_k^{-1}{g}_k\text{\hspace{1em}(B.15)}$$
   因为$H_k^{-1}$正定，故有$g_k^T{H}_k^{-1}{g}_k>0$。当$\lambda$为一个充分小的正数时，总有$f(x)<f(x^{(k)}),$也就是说$p_k$是下降方向。
   $$G_{k+1}{y}_k={\delta }_k\text{\hspace{1em}(B.16)}$$
   ​ 拟牛顿法将$G_k$作为$H_k^{-1}$的近似或选择$B_k$作为$H_k$的近似的算法称为拟牛顿法。

   ​ 按照拟牛顿法条件，在每次迭代中可以选择更新矩阵$G_{k+1}$;
   $$G_{k+1}=G_k+△G_k\text{\hspace{1em}(B.17)}$$
   这种选择有一定的灵活性，因此有多种具体实现方法。下面介绍Broyden类拟牛顿法。

   3. DFP(Davidon-Flecher-Powell)算法

      DFP算法选择$G_{k+1}$的方法是，假设每一步迭代中矩阵$G_{k+1}$是由$G_k$加上两个附加项构成的，即
      $$G_{k+1}=G_k+P_k+Q_k\text{\hspace{1em}(B.18)}$$
      其中$P_k,Q_k$是待定矩阵，这时。
      $$G_{k+1}{y}_k=G_k{y}_k+P_k{y}_k+Q_k{y}_k\text{\hspace{1em}(B.19)}$$
      为使$G_{k+1}$满足拟牛顿条件，可以$P_k$和$Q_k$满足：
      $$\begin{gathered} P_k{y}_k={\delta }_k \\ \text{\hspace{1em}(B.20)} \end{gathered}$$

      $$Q_k{y}_k=-G_k{y}_k\text{\hspace{1em}(B.21)}$$

      事实上，不难找出这样的$P_k$和$Q_k$，例如取
      $$P_k=\frac{{\delta }_k{\delta }_k^T}{{\delta }_k^T{y}_k}\text{\hspace{1em}(B.22)}$$

      $$Q_k=-\frac{G_k{y}_k{y}_k^T{G}_k}{y_k^T{G}_k{y}_k}\text{\hspace{1em}(B.23)}$$

      这样就可得到矩阵$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}\text{\hspace{1em}(B.24)}$$
      称为DFP算法：

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

      DFP算法如下：

      算法B.2(DFP算法)

      输入：目标函数f(x),梯度g(x)=$\nabla f(x)$,精度要求$ϵ$;

      输出：f(x)的极小点$x^{\ast }$。

      (1)选定初始点$x^{(0)}$,取$G_0$为正定对称矩阵，置k=0。

      (2)计算$g_k=g(x^{(k)})$。若$||g_k||<ϵ$，则停止计算，得近似解$x^{\ast }=x^k$;否则转(3)。

      (3)置$p_k=-G_k{g}_k$。

      (4)一维搜索：求${\lambda }_k$使得
      $$f(x^{(k)}+{\lambda }_k{p}_k)=\underset{\lambda \ge 0}{\min }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_{k+1}||<ϵ$,则停止计算，得近似解$x^{\ast }=x^{(k+1)};$否则，按式(B.24)算出$G_{k+1}$。

      (7)置k=k+1,转(3)。

3. BFGS(Broyden-Fletcher-Goldfarb-Shanno)算法(BFGS algorithm)

   BFGS算法是最流行得拟牛顿算法。

   可以考虑用$G_k$逼近黑塞矩阵得逆矩阵$H^{-1}$，也可以考虑用$B_k$逼近黑塞矩阵H。这时，相应得拟牛顿条件是
   $$B_{k+1}{\delta }_k=y_k\text{\hspace{1em}(B.25)}$$
   可以用同样的方法得到另一迭代公式。首先令
   $$B_{k+1}=B_k+P_k+Q_K\text{\hspace{1em}(B.26)}$$

   $$B_{k+1}{\delta }_k=B_k{\delta }_k+P_k{\delta }_k+Q_k{\delta }_k\text{\hspace{1em}(B.27)}$$

   考虑使$P_k$和$Q_k$满足：
   $$P_k{\delta }_k=y_k\text{\hspace{1em}(B.28)}$$

   $$Q_k{\delta }_k=-B_k{\delta }_k\text{\hspace{1em}(B.29)}$$

   找出适合条件的$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}\text{\hspace{1em}(B.30)}$$
   可以证明，如果初始矩阵$B_0$是正定的，则迭代过程中的每个矩阵$B_k$都是正定的。下面写出BFGS拟牛顿法。

   算法B.3(BFGS算法)

   输入：目标函数f(x)，g(x)=$\nabla f(x)$,精度要求$ϵ$;

   输出：f(x)的极小点$x^{\ast }$。

   (1)选定初始点$x^{(0)}$,取$G_0$为正定对称矩阵，置k=0。

   (2)计算$g_k=g(x^{(k)})$。若$||g_k||<ϵ$，则停止计算，得近似解$x^{\ast }=x^k$;否则转(3)。

   (3)置$B_k{p}_k=-g_k$,求出$p_k$。

   (4)一维搜索：求${\lambda }_k$使得
   $$f(x^{(k)}+{\lambda }_k{p}_k)=\underset{\lambda \ge 0}{\min }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_{k+1}||<ϵ$,则停止计算，得近似解$x^{\ast }=x^{(k+1)};$否则，按式(B.30)算出$B_{k+1}$。

   (7)置k=k+1,转(3)。

4. Broyden类算法(Broyden’s algorithm)

   我们可以从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}=(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}\text{\hspace{1em}(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}\text{\hspace{1em}(B.32)}$$
   也满足拟牛顿条件式，而且是正定的。其中$0\le \alpha \le 1$。这样就得到了一类拟牛顿法，称为Broyden类算法。
   [1] 李航. 统计学习方法-2版.北京：清华大学出版社，2019