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

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

学习李航的《统计学习方法》时，关于牛顿法和拟牛顿法的笔记。

牛顿法（Newton method）和拟牛顿法（quasi-Newton method）时求解无约束优化问题的常用方法，有收敛速度快的优点。牛顿法时迭代算法，每一步需要求解目标函数的 Hession 矩阵的逆矩阵，计算较为复杂；拟牛顿法通过正定矩阵近似 Hession 矩阵或其逆矩阵，简化了计算过程。

牛顿法

牛顿法的推导：考虑无约束最优化问题：
$$\underset{x\in {\mathbb{R}}^n}{\min }f(x)$$
设 $x^{\ast }$ 是目标函数的极小点。

假设 $f(x)$ 具有二姐连续偏导数，若第 $k$ 次迭代的值为 $x^{(k)}$ ，则对 $f(x)$ 在 $x^{(k)}$ 附近二阶泰勒展开：
$$f(x)=f(x^{(k)})+g_k^{\mathrm{T}}(x-x^{(k)})+\frac{1}{2}(x-x^{(k)})H_k(x-x^{(k)})$$
其中 $g_k$ 为 $f(x)$ 在 $x^{(k)}$ 处的梯度向量，$H_k$ 为 $f(x)$ 的 Hession 矩阵在 $x^{(k)}$ 处的值：
$$g(x)=\nabla f(x)={\left[\frac{\partial f}{\partial x_i}\right]}_{n\times 1},H(x)={\left[\frac{{\partial }^{2}f}{\partial x_i\partial x_j}\right]}_{n\times n}$$
函数 $f(x)$ 有极值的必要条件是在极值点处一阶导数（即梯度向量）为 $0$ 。特别地，当 $H(x)$ 在极值点处是正定矩阵时，函数 $f(x)$ 的极值为极小值。

牛顿法利用极小点的必要条件，每次从上一次迭代得到的极小点 $x^{(k)}$ 开始，得到当前目标函数的极小点 $x^{(k+1)}$ ，即：
$$\nabla f(x^{(k+1)})=g_k+(x^{(k+1)}-x^{(k)})=0$$
解得：
$$x^{(k+1)}=x^{(k)}-H_k^{-1}{g}_k$$
假设 $p_k$ 为一 $n\times 1$ 的列向量，且满足：
$$H_k{p}_k=-g_k$$
则迭代公式可写为：
$$x^{(k+1)}=x^{(k)}+p_k$$
算法：牛顿法

• 输入：目标函数 $f(x)$ ，梯度 $g(x)=\nabla f(x)$ ，Hession 矩阵 $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$ ，$k=k+1$ ，跳转至 2；

拟牛顿法

牛顿法中需要计算 $H_k^{-1}$ ，比较耗时，我们考虑使用一个与 Hession 矩阵具有类似性质的 $n$ 阶矩阵 $G_k=G(x^{(k)})$ 来近似代替 $H_k^{-1}$ 。

条件一：首先，对 $x^{(k)}$ 处的二阶泰勒展开进行求导，取 $x{=}^{(k+1)}$ ，得：
$$g_{k+1}-g_k=H_k(x^{(k+1)}-x^{(k)})$$
记 $y_k=g_{k+1}-g_k$ ，${\delta }_k=x^{(k+1)}-x^{(k)}$ ，则满足：
$$y_k=H_k{\delta }_k$$
或：
$$H_k^{-1}{y}_k={\delta }_k$$
该条件称为拟牛顿条件。

条件二：如果 $H_k$ 是正定的（$H_k^{-1}$ 也是正定的），那么可以保证牛顿法搜索方向 $p_k$ 是下降方向，因为：
$$x^{(k+1)}=x^{(k)}+\lambda p_k=x^{(k)}-\lambda H_k^{-1}{g}_k$$
带入前面的泰勒展开式，为：
$$f(x^{(k+1)})=f(x^k)-\lambda g_k^{\mathrm{T}}{H}_k^{-1}{g}_k+\frac{1}{2}{\lambda }^2{g}_k^{\mathrm{T}}{H}_k^{-\mathrm{T}}{H}_k{H}_k^{-1}{g}_k\approx f(x^k)-\lambda g_k^{\mathrm{T}}{H}_k^{-1}{g}_k$$
当 $\lambda$ 为一个足够小的正数时，由于 $H_k$ 正定，即 $\lambda g_k^{\mathrm{T}}{H}_k^{-1}{g}_k>0$ ，因此有 $f(x^{(k+1)})<f(x^{(k)})$ ，也就是说 $p_k$ 时下降方向。

综合前面两个条件，我们在选择近似矩阵 $G_k$ 时，首先要求每次迭代时 $G_k$ 是正定的，其次 $G_k$ 需要满足以下的拟牛顿条件：
$$G_{k+1}{y}_k={\delta }_k$$
同时，按照拟牛顿条件，每次迭代可以选择更新矩阵 $G_{k+1}$ ：
$$G_{k+1}=G_k+\Delta G_k$$
有多种选择近似矩阵的方法：

DFP 算法

DFP （Davidon-Fletcher-Powell）算法对 $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={\delta }_k$$
为了满足拟牛顿条件，可以令：
$$P_k{y}_k={\delta }_k,Q_k{y}_k=-G_k{y}_k$$
书上给了一个例子：
$$P_k=\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k},Q_k=-\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{G}_k{y}_k}$$
可以证明（咱也不知道怎么证明），如果 $G_0$ 是正定的，则迭代过程中的每个矩阵 $G_k$ 都是正定的。

算法：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(x^{(k)})\Vert <\epsilon$ ，则停止计算，返回近似解 $x^{\ast }=x^{(k)}$ ，否则继续；
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}$ ，置 $k=k+1$ ，转 2；

$$G_{k+1}=G_k+\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}-\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{G}_k{y}_k}$$

BFGS 算法

BFGS（Broyden-Fletcher-Goldfarb-Shanno）算法是最流行的拟牛顿算法。此时考虑使用 $B_k$ 近似 $H_k$ ，更新方法类似：
$$B_{k+1}=B_k+P_k+Q_k$$
为了满足拟牛顿条件：
$$B_{k+1}{\delta }_k=B_k{\delta }_k+P_k{\delta }_k+Q_k{\delta }_k=y_k$$
取：
$$P_k{\delta }_k=y_k,Q_k{\delta }_k=-B_k{\delta }_k$$
同样可令：
$$P_k=\frac{y_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k},Q_k=-\frac{B_k{\delta }_k{\delta }_k^{\mathrm{T}}{B}_k}{{\delta }_k^{\mathrm{T}}{B}_k{\delta }_k}$$
算法：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(x^{(k)})\Vert <\epsilon$ ，则停止计算，返回近似解 $x^{\ast }=x^{(k)}$ ，否则继续；
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. 计算 $B_{k+1}$ ，置 $k=k+1$ ，转 2；

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

Broyden 类算法

根据 BFGS 算法中 $B_k$ 矩阵的迭代公式，我们可以得到对应的 $G_k$ 矩阵。有 $G_k=B_k^{-1}$ ，$G_{k+1}^{-1}=B_{k+1}^{-1}$ ；有：
$$\begin{array}{rl}{G}_{k+1}^{-1}=& G_k^{-1}+\frac{y_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}-\frac{G_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k}\\ \Rightarrow G_{k+1}=& {\left(G_k^{-1}+\frac{y_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}+\frac{-G_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k}\right)}^{-1}\end{array}$$
需要用到 Shermax-Morrison 公式：
$$(A+u{v}^{\mathrm{T}}{)}^{-1}=A^{-1}-\frac{A^{-1}u{v}^{\mathrm{T}}{A}^{-1}}{1+v^{\mathrm{T}}{A}^{-1}u}$$
其中 $u$ 和 $v$ 为 $n$ 维向量，$A$ 为可逆矩阵；记：
$$A=G_k^{-1}+\frac{y_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k},u=-G_k^{-1}{\delta }_k,v^{\mathrm{T}}=\frac{{\delta }_k^{\mathrm{T}}{G}_k^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k}$$
则：
$$\begin{array}{rl}{G}_{k+1}=& {\left(A+\frac{-G_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k}\right)}^{-1}\\ =& A^{-1}+\frac{A^{-1}\frac{G_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k}{A}^{-1}}{1-\frac{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{A}^{-1}{G}_k^{-1}{\delta }_k}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k}}\\ =& A^{-1}+\frac{A^{-1}{G}_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}{A}^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k-{\delta }_k^{\mathrm{T}}{G}_k^{-1}{A}^{-1}{G}_k^{-1}{\delta }_k}\end{array}$$
对 $A$ 再次使用公式，此时 $u=y_k$ ，$v^{\mathrm{T}}=\frac{y_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}$ ，得到：
$$A^{-1}={\left(G_k^{-1}+\frac{y_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}\right)}^{-1}=G_k-\frac{G_k\frac{y_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}{G}_k}{1+\frac{y_k^{\mathrm{T}}{G}_k{y}_k}{y_k^{\mathrm{T}}{\delta }_k}}=G_k-\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}$$
将 $A^{-1}$ 逐步代入 $G_{k+1}$ ，得到：
$$\begin{array}{rl}{G}_{k+1}=& A^{-1}+\frac{A^{-1}{G}_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}{A}^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k-{\delta }_k^{\mathrm{T}}{G}_k^{-1}{A}^{-1}{G}_k^{-1}{\delta }_k}\\ =& A^{-1}+\frac{A^{-1}{G}_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}{A}^{-1}}{{\delta }_k^{\mathrm{T}}{G}_k^{-1}{\delta }_k-{\delta }_k^{\mathrm{T}}{G}_k^{-1}\left(G_k-\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}\right)G_k^{-1}{\delta }_k}\\ =& A^{-1}+\frac{A^{-1}{G}_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}{A}^{-1}}{\frac{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}}\\ =& A^{-1}+\left(G_k-\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}\right)\left(\frac{G_k^{-1}{\delta }_k{\delta }_k^{\mathrm{T}}{G}_k^{-1}}{\frac{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}}\right)\left(G_k-\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}\right)\\ =& A^{-1}+\left(I-\frac{G_k{y}_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}\right)\left(\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{\frac{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}}\right)\left(I-\frac{y_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}\right)\\ =& A^{-1}+\left(\frac{({\delta }_k{\delta }_k^{\mathrm{T}})(y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k)}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)-\left(\frac{{\delta }_k{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{G}_k}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)-\left(\frac{G_k{y}_k{y}_k^{\mathrm{T}}{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)+\left(\frac{G_k{y}_k{y}_k^{\mathrm{T}}{\delta }_k{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{G}_k}{({\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k)(y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k)}\right)\\ =& G_k-\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}\\ +& \left(\frac{({\delta }_k{\delta }_k^{\mathrm{T}})(y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k)}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)-\left(\frac{{\delta }_k{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{G}_k}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)-\left(\frac{G_k{y}_k{y}_k^{\mathrm{T}}{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)+\left(\frac{G_k{y}_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k+y_k^{\mathrm{T}}{G}_k{y}_k}\right)\\ =& G_k+\left(\frac{({\delta }_k{\delta }_k^{\mathrm{T}})(y_k^{\mathrm{T}}{\delta }_k)}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)+\left(\frac{({\delta }_k{\delta }_k^{\mathrm{T}})(y_k^{\mathrm{T}}{G}_k{y}_k)}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)-\left(\frac{{\delta }_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k}\right)-\left(\frac{G_k{y}_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)\\ =& G_k-\left(\frac{G_k{y}_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)-\left(\frac{{\delta }_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k}\right)+\left(\frac{{\delta }_k(y_k^{\mathrm{T}}{G}_k{y}_k){\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k{y}_k^{\mathrm{T}}{\delta }_k}\right)+\left(\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)\\ =& G_k\left(I-\frac{y_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)-\left(\frac{{\delta }_k{y}_k^{\mathrm{T}}{G}_k}{y_k^{\mathrm{T}}{\delta }_k}\right)\left(I-\frac{y_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)+\left(\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)\\ =& \left(I-\frac{{\delta }_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}\right)G_k\left(I-\frac{y_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)+\left(\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)\\ =& \left(I-\frac{{\delta }_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}\right)G_k{\left(I-\frac{{\delta }_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}\right)}^{\mathrm{T}}+\left(\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}\right)\end{array}$$
因此得到 BFGS 算法的 $G_k$ 的迭代公式：
$$G_{k+1}^{\text{BFGS}}=\left(I-\frac{{\delta }_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}\right)G_k^{\text{BFGS}}{\left(I-\frac{{\delta }_k{y}_k^{\mathrm{T}}}{y_k^{\mathrm{T}}{\delta }_k}\right)}^{\mathrm{T}}+\frac{{\delta }_k{\delta }_k^{\mathrm{T}}}{{\delta }_k^{\mathrm{T}}{y}_k}$$
由 DFP 算法得到的 $G_k$ 记作 $G_k^{\text{DFP}}$ ，可知二者的线性组合也满足拟牛顿条件式，而且是正定的：
$$G_{k+1}=\alpha G_k^{\text{DFP}}+(1-\alpha )G_k^{\text{BFGS}},0\le \alpha \le 1$$
这样就得到了一类拟牛顿算法，称为 Broyden 类算法。