牛顿法

牛顿法(Newton Method)

牛顿法其实是一个寻找函数的根的方法，而后运用到优化问题中。主要是运用泰勒级数的知识来寻找一个可导函数 $f(x)=0$ 的根。

将函数 $f(x)$ 在 $x_0$ 处展开成带皮亚诺余项的泰勒展开式：

$$f(x) = f(x_0) + f'(x_0)(x-x_0)+o(x)$$

取前两项，作为 $f(x)$ 的近似，那么就可以用

$$f(x_0) + f'(x_0)(x-x_0)=0$$

的解来近似 $f(x)$ 的解，求得的解为

$$x_1 = x_0-\frac{f(x_0)}{f'(x_0)}$$

如下图：

由于函数只有一阶展开，求解处的值离真实的根还有一定的差距，但是已经很接近结果了，我们只需要再将 $f(x)$ 在 $x_1$ 处泰勒展开，重复以上的步骤，那么就会得到比 $x_1$ 更加接近的解 $x_2$。

只要我们重复上述的步骤，循环迭代 $n$ 次，就会得到近似解 $x_n$

$$x_n =x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})}$$

利用牛顿法求解最优化问题

对于一求个无约束的最优化问题 $\min_{\mathbf x} f(\mathbf x),\mathbf x\in \mathbb R^n$ 的最优解，可以将其转换为求解其梯度等于 $0$ 时的根，即 $\nabla f(\mathbf x)=0$。你看，这就可以利用牛顿法求解最优化问题啦。

采用牛顿法求解的迭代式为：

$$\mathbf x_{k+1} = \mathbf x_k -H^{-1}(\mathbf x_k) \cdot \nabla f(\mathbf x_k)$$

这个 $H(\mathbf x)$ 是海塞矩阵，用来描述一个多元函数二阶导数。

海塞矩阵(Hessian Matrix)

对于一个多元函数 $f(x_1,x_2,\dots,x_n)$，如果函数的二阶导数都存在，则定义 $f$ 的海塞矩阵为

$$H(f)= \begin{pmatrix} \frac{\partial^2f}{\partial x_1^2} & \frac{\partial^2f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_1 \partial x_n} \\ \frac{\partial^2f}{\partial x_2 \partial x_1} & \frac{\partial^2f}{\partial x_2^2} & \cdots & \frac{\partial^2f}{\partial x_1 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2f}{\partial x_n \partial x_1} & \frac{\partial^2f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2f}{\partial x_n^2} \\ \end{pmatrix}$$

如果多元函数 $f(x_1,x_2,\dots,x_n)$ 二阶连续可导，并且在某一点 $\alpha(x_1’,x_2’,\dots,x_n’)$ 处梯度为 $0$，即$\nabla f(\alpha)=0$，那么点 $\alpha$ 为函数驻点，但不能判断该点是否为极小值。通过之前学习的知识，此时需要对函数求二阶导。

记 $f$ 在点 $\alpha$ 处的海塞矩阵为 $H(\alpha)$，因为

$$\frac{\partial^2f}{\partial x_i \partial x_j}=\frac{\partial^2f}{\partial x_j \partial x_i}$$

那么海塞矩阵是一个对称矩阵，对于 $H(\alpha)$，有以下结论

• 如果 $H(\alpha)$ 是正定矩阵，则驻点 $\alpha$ 为极小值点
• 如果 $H(\alpha)$ 是负定矩阵，则驻点 $\alpha$ 为极大值点
• 如果 $H(\alpha)$ 是不定矩阵，则驻点 $\alpha$ 不是极值点

---

因为牛顿法在进行迭代时运用了二阶导数，也就是说用二次曲线连模拟点 $\mathbf x_k$，一步便可以找到二次曲线的极值点，所以下降的速度比起梯度下降法更快，能更容易地找到答案；但由于 $H^{-1}(\mathbf x_k)$ 难于求解，难于操作，于是一种牛顿法的优化方法便出现了，就是下面提到的拟牛顿法。

LM 修正牛顿法

Levenberg-Marquardt 修正牛顿法主要着手于处理 $H^{-1}(\mathbf x_k)$ 为不可逆的情况，令

$$G(\mathbb x_k) = H(\mathbf x_k) + \mu I$$

使 $\mu$ 充分大，则 $G(\mathbf x_k)$ 的特征值 $\lambda_i + \mu$ 均大于 $0$。

当 $\mu \to \infty$ 时， $H^{-1}(\mathbf x_k) \cdot \nabla f(\mathbf x_k)$ 趋近于 $\frac 1 \mu \nabla f(\mathbf x_k)$，那么该法就趋近于梯度下降法。

高斯-牛顿法(Gauss-Newton Method)

高斯-牛顿法是用来解非线性的最小二乘问题。设有 $m$ 个样本点 ${ (x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\cdots,(x^{(m)},y^{(m)}) }$，$x^{(i)}\in \mathbb R^n$，用于训练一个可以用来预测的函数 $f_\theta(x)$。

有别于其他模型，这里我们特别地记 $f_i(x)=f(x^{(i)};\theta)$ 。有损失函数：

$$L(\theta) = \sum_{i=1}^m(f_i(\theta) - y^{(i)})^2$$

上式对 $\theta_j$ 进行求导，有

$$\frac{\partial L}{\partial \theta_j} = \sum_{i=1}^m\left[2(f_i(\theta) - y^{(i)})\frac{\partial f_i}{\partial\theta_j}\right]$$

写成向量和矩阵形式：

$$\begin{align} \nabla L(\theta) &= \begin{pmatrix} \frac{\partial L}{\partial \theta_1}\\ \frac{\partial L}{\partial \theta_2}\\ \vdots\\ \frac{\partial L}{\partial \theta_n} \end{pmatrix} =2 \begin{pmatrix} \frac{\partial f_1}{\partial \theta_1} & \frac{\partial f_2}{\partial \theta_1} & \cdots&\frac{\partial f_m}{\partial \theta_1}\\ \frac{\partial f_1}{\partial \theta_2} & \frac{\partial f_2}{\partial \theta_2} & \cdots&\frac{\partial f_m}{\partial \theta_2}\\ \vdots & \vdots & \ddots & \vdots\\ \frac{\partial f_1}{\partial \theta_n} & \frac{\partial f_2}{\partial \theta_n} & \cdots&\frac{\partial f_m}{\partial \theta_n} \end{pmatrix} \begin{pmatrix} (f_1(\theta) - y^{(1)})\\ (f_2(\theta) - y^{(2)})\\ \vdots\\ (f_m(\theta) - y^{(m)}) \end{pmatrix} \\ &=2J^Tr \end{align}$$

其中 $J$ 为 $f(x;\theta)$ 关于 $\theta$ 的雅可比矩阵，$\mathbf r$ 为样本与模型的残差。

接下来我们根据已经求得的梯度，更进一步求海森矩阵的元素 $h_{jk}$

$$\begin{align} h_{jk} &= \frac{\partial^2 L}{\partial\theta_j\theta_k}\\ &=\frac{\partial}{\partial\theta_j}\left(\sum_{i=1}^m\left[2(f_i(\theta) - y^{(i)})\frac{\partial f_i}{\partial\theta_k}\right]\right)\\ &= 2\sum_{i=1}^m \left( \frac{\partial f_i}{\partial\theta_j}\frac{\partial f_i}{\partial\theta_k} + (f_i(\theta) - y^{(i)})\frac{\partial^2 f_i}{\partial\theta_j\theta_k} \right) \end{align}$$

同样的，写成矩阵形式

$$H(\theta) = 2(J^TJ+S)$$

其中

$$s_{ij} = \sum_{i=1}^m(f_i(\theta) - y^{(i)})\frac{\partial^2 f_i}{\partial\theta_j\theta_k}$$

在迭代过程中，由于残差项很小，往往（?）在迭代过程中可以忽略，则我们找到了一个新的矩阵 $J^TJ$ 来代替海塞矩阵 $H$，回归到求解 $\mathbf x$ 的问题，得到迭代式

$$\mathbf x_{k+1} = \mathbf x_k -(J^TJ)^{-1} J^T\mathbf r$$

拟牛顿法(Quasi-Newton Methods)

由于海塞矩阵矩阵不一定可逆，而且求矩阵逆的复杂度较高，所以引入了拟牛顿法。拟牛顿法其实就是准找可以替代矩阵 $H^{-1}(\mathbf x_k)$ 的方法。

我们对 $\nabla f(\mathbf x)$ 在 $x_k$ 点处做一阶泰勒展开：

$$\nabla f(\mathbf x) = \nabla f(\mathbf x_k) + H(\mathbf x_{k})\cdot(\mathbf x-\mathbf x_k)$$

或者说对 $f(\mathbf x)$ 在 $x_k$ 点处做二阶泰勒展开：

$$f(\mathbf x) = f(\mathbf x_k) + (\mathbf x-\mathbf x_k)^T\nabla f(\mathbf x_k) + \frac 1 2 (\mathbf x-\mathbf x_k)^TH(\mathbf x_{k})(\mathbf x-\mathbf x_k)$$

并对其求导。

取 $\mathbf x = \mathbf x_{k+1}$ 有：

$$\nabla f(\mathbf x_{k+1})-\nabla f(\mathbf x_k) = H(\mathbf x_k)\cdot(\mathbf x_{k+1}-\mathbf x_k)$$

记 $\mathbf yk = \nabla f(\mathbf x{k+1})-\nabla f(\mathbf xk)$，$\delta_k = \mathbf x{k+1}-\mathbf x_k$，有：

$$\mathbf y_k = H(\mathbf x_k) \cdot \delta_k$$

上式被称为拟牛顿条件。根据对矩阵 $H^{-1}(\mathbf x_k)$ 或者 $H(\mathbf x_k)$ 估计方法的不同，拟牛顿法可以分为以下的几种：

DFP 变尺度法

DFP 方法采用矩阵 $G(\mathbf x_k)$ 来对 $H^{-1}(\mathbf x_k)$ 进行近似，最早是由 Davidon 提出，后经 Fletcher 和 Powell 解释和改进，在命名时以三个人名字的首字母命名。

在 DFP 方法中，假设迭代式：

$$G(\mathbf x_{k+1}) = G(\mathbf x_{k}) + \Delta G(\mathbf x_{k})$$

其中 $G(\mathbf x_{k})$ 是正定的，令

$$\Delta G_k = \alpha u_ku_k^T + \beta v_kv_k^T, \quad u_k,v_k \in \mathbb R^n$$

代入上上式可得

$$G(\mathbf x_{k+1}) = G(\mathbf x_{k}) + \alpha u_ku_k^T + \beta v_kv_k^T$$

已知有 $ H^{-1}(\mathbf x_k) \cdot \mathbf y_k =\delta_k$，则

$$\begin{align} \delta_k =& G(\mathbf x_{k})\mathbf y_k + \alpha u_ku_k^T\mathbf y_k + \beta v_kv_k^T\mathbf y_k\\ =& G(\mathbf x_{k})\mathbf y_k + \alpha (u_k^T\mathbf y_k)u_k + \beta (v_k^T\mathbf y_k)v_k \end{align}$$

可知 $u_k^T\mathbf y_k,v_k^T\mathbf y_k$ 为数，我们设 $u_k = r\delta_k$， $v_k = \theta G(\mathbf x_k)\mathbf y_k$，有

$$\delta_k = G(\mathbf x_{k})\mathbf y_k+\alpha( r\delta_k^T\mathbf y_k) r\delta_k + \beta (\theta \mathbf y_k^TG(\mathbf x_k)\mathbf y_k)\theta G(\mathbf x_k)\mathbf y_k\$$

整理得：

$$(\alpha r^2(\delta_k^T\mathbf y_k)-1)\delta_k+(\beta\theta^2(\mathbf y_k^TG(\mathbf x_k)\mathbf y_k)+1) G(\mathbf x_k)\mathbf y_k=0$$

令

$$\begin{cases} \alpha r^2(\delta_k^T\mathbf y_k) = 1\\ \beta\theta^2(\mathbf y_k^TG(\mathbf x_k)\mathbf y_k)=-1 \end{cases}$$

最终代入可得

$$G(\mathbf x_{k+1}) = G(\mathbf x_k) + \frac{\delta_k\delta_k^T}{\delta_k^T\mathbf y_k} - \frac{G(\mathbf x_k)\mathbf y_k\mathbf y_k^TG(\mathbf x_k)}{\mathbf y_k^TG(\mathbf x_k)\mathbf y_k}$$

BFGS 方法

BFGS 算法采用矩阵 $B(\mathbf x_k)$ 来对 $H(\mathbf x_k)$ 进行近似，采取上节 DFP 方法的推导方法，可得

$$B(\mathbf x_{k+1}) = B(\mathbf x_k) + \frac{\mathbf y_k\mathbf y_k^T}{\mathbf y_k^T\delta_k} - \frac{B(\mathbf x_k)\delta_k \delta_k^TB(\mathbf x_k)}{\delta_k^TB(\mathbf x_k)\delta_k}$$

可记 $G(\mathbf x{k})=B(\mathbf x{k})^{-1}$，那么有

$$G(\mathbf x_{k+1})=\left(B(\mathbf x_k) + \frac{\mathbf y_k\mathbf y_k^T}{\mathbf y_k^T\delta_k} - \frac{B(\mathbf x_k)\delta_k \delta_k^TB(\mathbf x_k)}{\delta_k^TB(\mathbf x_k)\delta_k}\right)^{-1}$$

本式可以运用 Sherman-Morrison-Woodbury 公式求解

Sherman-Morrison-Woodbury 公式是一种矩阵求逆的方法，公式如下：

$$(A+UCV)^{-1} = A^{-1}-\frac{A^{-1}UVA^{-1}}{C^{-1}+VA^{-1}U}$$

还有一个向量版的，叫 Sherman-Morrison 公式：

$$(A+uv^{T})^{-1} = A^{-1}-\frac{A^{-1}uv^{T}A^{-1}}{1+v^{T}A^{-1}u}$$

经过上面两个计算的不停的折腾（完整推倒过程可以 点这里 ），求解得：

$$G(\mathbf x_{k+1}) = (I-\frac{\delta_k\mathbf y_k^T}{\mathbf y^T_k\delta_k}) G(\mathbf x_{k}) (I-\frac{\mathbf y_k\delta_k^T}{\mathbf y^T_k\delta_k})+\frac{\delta_k\delta_k^T}{\mathbf y^T_k\delta_k}$$

令 $\rho = \frac 1 {\mathbf y^T_k\delta_k}$，可以得出最终的 BFGS 方法的迭代方程：

$$G(\mathbf x_{k+1}) = (I-\rho{\delta_k\mathbf y_k^T}) G(\mathbf x_{k}) (I-\rho{\mathbf y_k\delta_k^T})+\rho{\delta_k\delta_k^T}$$

L-BFGS 算法

这个算法的英文全称为Limited-memory BFGS（或 Limited-storage BFGS），意思就是上面的 BFGS 方法同样需要消耗计算机的很多计算能力，于是需要一个简洁版的算法。这个 L-BFGS 算法就是这个意思。

但是我现在还不想写，想看的看下面这个博客

https://liuxiaofei.com.cn/blog/lbfgs%e6%96%b9%e6%b3%95%e6%8e%a8%e5%af%bc/