牛顿法和拟牛顿法

牛顿法（Newton method）和拟牛顿法（quasi Newton method）也是求解无约束最优化问题的常用方法。方法使用函数 $f(x)$的一二阶导数。具有收敛速度快、迭代次数少、计算量大的特点。

问题描述

考虑如下无约束极小化问题：

$$\min_{x\in R^n} \; f(x)$$

最优点记为$x^\ast$

牛顿法

• 输入：目标函数$f(x)$，梯度$g(x)=\bigtriangledown f(x)$，海森矩阵$H(x)$，精度要求$\varepsilon $
• 输出：$f(x)$的极小点$x^\ast$
  1. 确定初始点$x_0$，置$k=0$。
  2. 计算 $g_k=g(x^{(k)})$
  3. 若$\left | g_k \right |<\varepsilon $，则停止计算，得近似解$x^\ast=x^{(k)}$
  4. 计算$H_k=H(x^{(k)})$
  5. 置$x^{(k+1)}=x^{(k)}-H_k^{-1}g_k$
  6. 置$k=k+1$，转第二步骤

其中$g_k=\bigtriangledown f(x^{(k)})$，$H_k=\bigtriangledown^2 f(x^{(k)})=\left [\frac{\partial^2 f}{\partial x_i\partial x_j} \right ]_{n\times n}$。$H_{k}^{-1}$为$x^{(k)}$的Hessian的逆矩阵,

$$g(f)=\bigtriangledown f=\begin{bmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \\ \vdots \\ \frac{\partial f}{\partial x_n} \\ \end{bmatrix} \; \; \; \; \; H(f)= \bigtriangledown^2 f= \begin{bmatrix} \frac{\partial^2 f}{\partial x^2_1} & \frac{\partial^2 f}{\partial x_1 \partial x_2} \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x^2_2} \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\\ \vdots & \ddots & \vdots \\\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} \cdots & \frac{\partial^2 f}{\partial x^2_n} \end{bmatrix}$$

代数推导

假设函数 $f(x)$ 二次可微，则在当前第$k$次迭代点 $x^{(k)}$ 对 $f(x)$ 进行二阶泰勒展开

$$f \left ( x \right ) \approx f \left ( x^{(k)} \right )+{f}'\left ( x^{(k)} \right )\left ( x-x^{(k)} \right )+\frac{1}{2}{f}''\left ( x^{(k)} \right )\left ( x-x^{(k)} \right )^2$$

求函数 $f(x)$ 极值则可以转化为对 $f(x)$ 求导并令其为0，

$$\frac{\partial f}{\partial x} \approx f'\left ( x^{(k)} \right )+{f}''\left ( x^{(k)} \right )\left ( x-x^{(k)} \right )=0$$

得到

$$\begin{align*} x&=x^{(k)}-\frac{f'\left ( x^{(k)} \right )}{f''\left ( x^{(k)} \right )} \\ &= x^{(k)}-H_{k}^{-1}g_k \\ x^{(k+1)} &= x^{(k)}-H_{k}^{-1}g_k \end{align*}$$

阻尼牛顿法

原始牛顿法由于迭代公式中没有步长因子，而是固定步长，有时会使迭代超过极值点，即$f(x^{(k+1)})>f(x^k)$，在牛顿方向上附加了步长因子，每次调整时会在搜索空间，在该方向找到最优步长，然后调整

$$\lambda _k=arg\; \underset{\lambda \in \mathbb{R}}{min}\; f(x_k+\lambda H_k^{-1}g_k)\;\;\;\;\;\;\; \tag 2$$

• 输入：目标函数$f(x)$，梯度$g(x)=\bigtriangledown f(x)$，海森矩阵$H(x)$，精度要求$\varepsilon $
• 输出：$f(x)$的极小点$x^\ast$
  1. 确定初始点$x_0$，置$k=0$。
  2. 计算 $g_k=g(x^{(k)})$
  3. 若$\left | g_k \right |<\varepsilon $，则停止计算，得近似解$x^\ast=x^{(k)}$
  4. 计算$H_k=H(x^{(k)})$
  5. 利用$(2)$式子得到步长$\lambda_k$置$x^{(k+1)}=x^{(k)}-\lambda_k H_k^{-1}g_k$
  6. 置$k=k+1$，转第二步骤

拟牛顿算法

拟牛顿思路

因为Hessian矩阵的计算量大而且无法保持正定所以采用拟牛顿算法。拟牛顿算法的基本思想是不用二阶偏导数而构造出可以近似Hessian矩阵(或者Hessian矩阵的逆矩阵)的正定对称矩阵，在拟牛顿的条件下优化目标函数

拟牛顿条件

$$\begin{align*} \frac{\partial f}{\partial x} &\approx f'\left ( x^{(k)} \right )+{f}''\left ( x^{(k)} \right )\left ( x-x^{(k)} \right )\\ f'\left ( x^{(k+1)} \right ) &\approx f'\left ( x^{(k)} \right )+{f}''\left ( x^{(k)} \right )\left ( x^{(k+1)}-x^{(k)} \right ) \\ g_{k+1}-g_k &\approx H_k\left ( x^{(k+1)}-x^{(k)} \right ) \end{align*}$$

拟牛顿法对$H_k$或$H_k^{-1}$取近似值，可减少计算量。记$B\approx H$，$D\approx H^{-1}$，$y_k=g_{k+1}-g_k$，$s_k=x_{k+1}-x_k$。根据拟牛顿条件，可得近似公式：

$$\color{red}{B_{k+1}=\frac{y_k}{s_k} \;\;\;\;\;\;\;\;D_{k+1}=\frac{s_k}{y_k}}$$

DFP算法

DFP算法采用的是$D$，但并不直接计算$D$，而是计算每一步DD的增量$ΔD$来间接的求出$D$。这也是很多优化算法的做法，因为一般上一步的中间结果对下一步的计算仍有价值，若直接抛弃重新计算耗时耗力耗内存，重新发明了轮子。

$$D_{k+1}=D_k+\Delta D_k$$

$D_0$通常取单位矩阵I，关键导出每一步的$ΔD_k$。
通过一系列艰苦而又卓绝的推导计算假设取便，最终的导出结果为：

$$\color{red}{\begin{align*} D_{k+1}&=D_k+\Delta D_k\; \; \; \; k=0,1,2,\cdots \\ \Delta D_k &=\frac{s_ks_k^T}{s_k^Ty_k} -\frac{D_ky_ky_k^TD_k}{y_k^TD_ky_k} \end{align*}}$$

一般来说，在进行中间增量计算时，都要经过这一步艰苦而又卓绝的推导计算。

BFGS算法

BFGS算法与DFP算法类似，只是采用的BB来近似HH。最终的公式为：

$$\color{red}{\begin{align*} \Delta B_k &=\frac{y_ky_k^T}{s_k^Ty_k} -\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k} \end{align*}}$$

L-BFGS

L-BFGS算法对BFGS算法进行改进，不再存储矩阵$D_k$，因为$D_k$有时候比较大，计算机的肚子盛不下。但是我们用到$D_k$的时候怎么办呢？答案是根据公式求出来。

从上面的算法推导可知，$D_k$只跟$D_0$和序列$\{s_k\}$和$\{y_k\}$有关。即我们知道了后者，即可以求得前者。进一步近似，我们只需要序列$\{s_k\}$和$\{y_k\}$的最近$m$个值即可。这样说来，我们的计算机内存中只需要存储这两个序列即可，瞬间卸掉了很多东西，正是春风得意马蹄轻。当然，这样cpu的计算量也会相应的增加，这是可以接受的，马，也是可以接受的。

最终的递推关系为

$$D_{k+1}=V^T_kD_kV_k+\rho_k s_ks^T_k$$

其中

$$\rho_k=\frac{1}{y^T_ks_k},V_k=I-\rho_ky_ks^T_k$$