统计学习方法--牛顿法和拟牛顿法

与公众号同步更新，详细内容及相关ipynb文件在公众号中，公众号：AI入门小白

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

黑塞矩阵是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。

牛顿法

考虑无约束最优化问题
$$\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)})=▽f(x^{(k)})$是$f(x)$的梯度向量在点$x^{(k)}$的值，$H(x^{(k)})$是$f(x)$的黑塞矩阵
$$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)$的极值为极小值。

牛顿法利用极小点的必要条件
$$▽f(x)=0\text{\hspace{1em}(B.4)}$$
每次选代中从点$x^{(k)}$开始，求目标函数的极小点，作为第$k+1$次迭代值$x^{(k+1)}$。具体地，假设$x^{(k+1)}$满足:
$$▽f(x^{(k+1)})=0\text{\hspace{1em}(B.5)}$$
由式(B.2) 有
$$▽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\text{\hspace{1em}(B.10)}$$
用式(B.8)作为迭代公式的算法就是牛顿法。

算法B.1(牛顿法)
输入：目标函数$f(x)$，梯度$g(x)=▽f(x)$，黑塞矩阵$H(x)$，精度要求$\epsilon$；
输出：$f(x)$的极小点$x^{\ast }$。
(1)取初始点$x^{(0)}$，置$k=0$。
(2)计算$g_k=g(x^{(k)})$。
(3)若$\Vert g_k\Vert <\epsilon$，则停止计算，得近似解$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}$，计算比较复杂，所以有其它改进的方法。

拟牛顿法的思路

在牛顿法的迭代中，需要计算黑塞矩阵的逆矩阵$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) 可以近似写成:
$$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$作为$H_k^{-1}$的近似，要求矩阵$G_k$满足同样的条件。首先，每次选代矩阵$G_k$是正定的。同时，$G_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+\Delta G_k\text{\hspace{1em}(B.17)}$$
这种选择有一定的灵活性，因此有多种具体实现方法。下面介绍Broyden 类拟牛顿法。

DFP (Davidon-Fletcher- Powell) 算法(DFP algorithm)

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$满足:
$$P_k{y}_k={\delta }_k\text{\hspace{1em}(B.20)}$$
$$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$都是正定的。

算法B.2(DFP算法)
输入：目标函数$f(x)$，梯度$g(x)=\nabla f(x)$，精度要求$\epsilon$；
输出：$f(x)$的极小点$x^{\ast }$。
(1)选定初始点$x^{(0)}$，取$G_0$为正定对称矩阵，置$k=0$。
(2)计算$g_k=g(x^{(k)})$。若$\Vert g_k\Vert <\epsilon$，则停止计算，得近似解$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)})$，若$\Vert g_{k+1}\Vert <\epsilon$，则停止计算，得近似解$x^{\ast }=x^{(k+1)}$；否则，按式(B.24)算出$G_{k+1}$。
(7)置$k=k+1$，转(3)。

BFGS (Broyden-Fletcher-Goldfarl-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$都是正定的。

算法B.3(BFGS算法)
输入：目标函数$f(x),g(x)=\nabla f(x)$，精度要求$\epsilon$；
输出：$f(x)$的极小点$x^{\ast }$。
(1)选定初始点$x^{(0)}$，取$B_0$为正定对称矩阵，置$k=0$。
(2)计算$g_k=g(x^{(k)})$。若$\Vert g_k\Vert <\epsilon$，则停止计算，得近似解$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)})$，若$\Vert g_{k+1}\Vert <\epsilon$，则停止计算，得近似解$x^{\ast }=x^{(k+1)}$；否则，按式(B.30)算出$B_{k+1}$。
(7)置$k=k+1$，转(3)。

Broyden 类算法(Broyden’s algorithm)

Sherman-Morrison 公式:假设$A$是$n$阶可逆矩阵，$u,v$是$n$维向量，且$A+u{v}^T$也是可逆矩阵，则
$$(A+u{v}^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}=(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 类算法。

数据来源：统计学习方法(第二版)