Levenberg-Marquardt算法(LM)的前世今生

Levenberg-Marquardt算法(LM)的前世今生

本文将从优化算法当中的梯度下降法的发展，从最基本的最速下降法逐个介绍，最后介绍LM方法，顺序如下：

1. 最速下降法
2. 牛顿法
3. 高斯牛顿法
4. LM算法（阻尼高斯牛顿法）

问题引入

优化算法的目标就是求解目标函数的最大值或最小值。我们假设要求解函数$F(x)$的最小值:
$$\underset{x}{\min }F(x)$$
假设目标函数为凸函数（随后的几种方法也是利用了最小二乘法，假设目标函数为开口朝上的二次函数，属于凸优化问题）。
直观方法：由于是凸优化问题，直接求$F^{\prime }(x)=0$的解，得到极小值。但是求解$\frac{dF}{dx}=0$过于复杂，采用逐步逼近极值点的方法（梯度下降法）求解。

最速下降法(一阶梯度下降法)

对$F(x)$作一阶泰勒展开：
$$F(x_k+\Delta x_k)\approx F(x_k)+J(x_k{)}^T\Delta x_k$$
其中$J(x_k{)}^T$为$F(x)$的雅克布矩阵（梯度矩阵，即一阶导）。要逐步逼近最小值，实际上就是确定$\Delta x_k$的取值情况（$x$应该朝哪个方向变化）。最速下降法直接令$\Delta x_k$取梯度的反向：
$$\Delta x_k=-\lambda J(x_k)$$
上式即为增量方程，$\lambda$为步长，对应模型训练当中的Learning rate。这就是反向传播当中最初的梯度下降方法，由此往后衍生出momentum，adagrad，adam等。

牛顿法(二阶梯度下降法)

对$F(x)$作二阶泰勒展开：
$$F(x_k+\Delta x_k)\approx F(x_k)+J(x_k{)}^T\Delta x_k+\frac{1}{2}\Delta x_k^{T}H(x_k)\Delta x_k$$
其中$H(x_k)$为$F(x)$的Hessian Matrix何塞矩阵。我们还是需要确定增量方程，把上式看成是以$\Delta x_k$为自变量的函数，则以$\Delta x_k$为变元对上式求导，导数为0时为极小值点：
$$\frac{\partial (F(x_k)+J(x_k{)}^T\Delta x_k+\frac{1}{2}\Delta x_k^{T}H(x_k)\Delta x_k)}{\partial \Delta x_k}=0$$
解得：
$$\Delta x_k=-\mathrm{H}({\mathrm{x}}_{\mathrm{k}}{)}^{-1}\mathrm{J}({\mathrm{x}}_{\mathrm{k}})$$
上式为牛顿法的增量方程。

缺点：Hessian计算量太大。

高斯牛顿法

针对最小二乘问题。假设$f(x)$为对数据拟合的函数（模型）。$F(x)$为目标函数：
$$\underset{\mathrm{x}}{\min }\mathrm{F}\left(\mathrm{x}\right)=||\mathrm{f}\left(\mathrm{x}\right)|{|}^2$$

不同于牛顿法直接对目标函数$\mathrm{F}\left(\mathrm{x}\right)$作泰勒展开，这里对$\mathrm{f}\left(\mathrm{x}\right)$作泰勒展开。$F(x)$在$x_k$处，具有增量$\Delta x_k$：（$x_k$为第$k$次迭代的变量）：
$$\underset{\mathrm{x}}{\min }\mathrm{F}\left({\mathrm{x}}_{\mathrm{k}}+\Delta {\mathrm{x}}_{\mathrm{k}}\right)=||\mathrm{f}({\mathrm{x}}_{\mathrm{k}}+\Delta {\mathrm{x}}_{\mathrm{k}})|{|}^2$$

对$\mathrm{f}\left(\mathrm{x}\right)$作一阶泰勒展开：
$$\mathrm{F}\left({\mathrm{x}}_{\mathrm{k}}+\Delta {\mathrm{x}}_{\mathrm{k}}\right)\approx ||\mathrm{f}\left({\mathrm{x}}_{\mathrm{k}}\right)+\mathrm{J}{\left({\mathrm{x}}_{\mathrm{k}}\right)}^{\mathrm{T}}\Delta {\mathrm{x}}_{\mathrm{k}}|{|}^2$$

此处$J(x)$为$f(x)$的雅克比矩阵。对上式以$\Delta x_k$求导，并令一阶导为0(采用分母布局)：

$$\begin{array}{r}\frac{\partial ||f(X_k)+J(X_k{)}^T\cdot \Delta X_k|{|}^2}{\partial \Delta X_k}=0\\ 2(f(X_k)+J(X_k{)}^T\cdot \Delta X_k)\cdot J(X_k)=0\\ J(X_k)\cdot (f(X_k)+J(X_k{)}^T\cdot \Delta X_k)=0\end{array}$$

$(f(X_k)+J(X_k{)}^T\cdot \Delta X_k)$为标量，交换位置：

$$\begin{array}{r}(f(X_k)+J(X_k{)}^T\cdot \Delta X_k)\cdot J(X_k)=0\\ \mathrm{J}\left({\mathrm{x}}_{\mathrm{k}}\right)\mathrm{J}{\left({\mathrm{x}}_{\mathrm{k}}\right)}^{\mathrm{T}}\mathrm{\Delta }{\mathrm{x}}_{\mathrm{k}}=-\mathrm{J}\left({\mathrm{x}}_{\mathrm{k}}\right)\mathrm{f}\left({\mathrm{x}}_{\mathrm{k}}\right)\end{array}$$

令：
$$\mathrm{H}({\mathrm{x}}_{\mathrm{k}})=\mathrm{J}({\mathrm{x}}_{\mathrm{k}})\mathrm{J}({\mathrm{x}}_{\mathrm{k}}{)}^{\mathrm{T}},\mathrm{g}({\mathrm{x}}_{\mathrm{k}})=-\mathrm{J}({\mathrm{x}}_{\mathrm{k}})\mathrm{f}({\mathrm{x}}_{\mathrm{k}})$$

解得：
$$\Delta x_k=-\mathrm{H}({\mathrm{x}}_{\mathrm{k}}{)}^{-1}\mathrm{g}({\mathrm{x}}_{\mathrm{k}})$$

这里的$\mathrm{H}({\mathrm{x}}_{\mathrm{k}})$相当于对Hessian的近似。

病态矩阵：求解方程组时如果对数据进行较小的扰动，则得出的结果具有很大波动，这样的矩阵称为病态矩阵。

最小二乘法：又称最小平方法，通过最小化误差函数的平方和，寻找最优解。

• 泰勒展开的本质：使用以$\Delta x$为变元的多项式函数近似原函数。
• 高斯牛顿法的目标函数可以看成是一个开口朝上的抛物线。因为对$f(x)$作一阶泰勒展开，相当于把$f(x)$近似成一次函数，而${\min }_{\mathrm{x}}\mathrm{F}\left(\mathrm{x}\right)=||\mathrm{f}\left(\mathrm{x}\right)|{|}^2$，相当于$\mathrm{F}\left(\mathrm{x}\right)$为开口朝上的二次函数，导数值为0时取最小值。

LM算法（阻尼高斯牛顿法）

又称为阻尼高斯牛顿法，阻尼可以理解为系数。具体做法是：在高斯牛顿法的基础上引入置信区间，从无约束的最小二乘问题变为带约束的最小二乘问题，从而可以使用拉格朗日乘数法求解：

$$\underset{\mathrm{x}}{\min }\mathrm{F}\left({\mathrm{x}}_{\mathrm{k}}+\Delta {\mathrm{x}}_{\mathrm{k}}\right)=||\mathrm{f}({\mathrm{x}}_{\mathrm{k}}+\Delta {\mathrm{x}}_{\mathrm{k}})|{|}^2,s.t.||\mathrm{D}\Delta {\mathrm{x}}_{\mathrm{k}}|{|}^2\le \mu$$

其中$||\mathrm{D}\Delta {\mathrm{x}}_{\mathrm{k}}|{|}^2\le \mu$为约束条件，$\mathrm{D}$为系数矩阵（n维向量，系数矩阵不同对应不同的LM算法），$\mu$为置信半径。
使用拉格朗日乘数法求解：
$$\mathrm{L}(\mathrm{\Delta }{\mathrm{x}}_{\mathrm{k}},\lambda )=||f({\mathrm{x}}_{\mathrm{k}}+\mathrm{\Delta }{\mathrm{x}}_{\mathrm{k}})|{|}^2+\lambda (||D\Delta {\mathrm{x}}_{\mathrm{k}}|{|}^2-\mu )$$

对上式求导，令导数为0，可得$\Delta x_k$的解（增量方程）：
$$(\mathrm{H}+\lambda {\mathrm{D}}^{\mathrm{T}}\mathrm{D})\mathrm{\Delta }{\mathrm{x}}_{\mathrm{k}}=\mathrm{g}({\mathrm{x}}_{\mathrm{k}})$$

1. 当$\lambda$很大时，$\Delta x_k\approx -\frac{1}{\mu }g(x_k)$，相当于一阶梯度法，当$x_k$距极值点较远时效果更好。
2. 当$\lambda$很小时，$\Delta x_k\approx -{\mathrm{H}}^{-1}g(x_k)$，相当于高斯牛顿法，当$x_k$距极值点较近时效果更好。
   由此可知，对$\lambda$的取值至关重要，需要根据当前$x_k$距离极值点的远近来动态调整。

• 可以使用近似程度量化指标$\mathbf{\rho }$动态控制$\lambda$：
  $$\rho =\frac{f(x_k+\Delta x_k)-f(x_k)}{J(x_k{)}^T\cdot \Delta x_k}$$
  上式中，分子为目标函数实际减少值，分母为近似函数（泰勒一阶展开，即一次函数）的减少值。$\rho$约接近于1说明近似函数拟合得越好。

1. $\rho$很大时，说明实际减小值远大于近似函数减小值，认为此时离极值点较远，因此可以放心大胆地选用较大的$\lambda$值。
2. $\rho$很小时，说明实际减小值远小于近似函数减小值，认为此时离极值点较近，因此应当谨慎选用较小的$\lambda$值，用较小的步长逐步逼近极小值。

我们可以通过二次函数图像说明：

可以发现，越靠近极值点，二次函数的一阶导数的绝对值越小（绝对值递减）。此时发现在梯度下降的过程中，相同$\Delta x_k$的情况下，实际函数的减小值总是小于近似函数的减小值。因此这种情况认为$x_k$离极值点较近。

问题：

1. 神经网络的目标函数都是凸函数吗？为什么？
2. Hessian计算量太大。故有些算法使用Fisher Information Matrix近似Hessian。为什么能这样做？两者的联系？
3. 如何理解：Network Automatic Pruning: Start NAP and Take a Nap中的这句话：

The first approximation comes from the fact that if the training objective is the negative log-likelihood, the Hessian matrix and Fisher matrix are the expectations of second-order derivatives under the data distribution and model distribution respectively.

4. 正定矩阵？奇异矩阵为什么会导致算法不收敛，增量不稳定？

参考链接：

1. 高斯牛顿法详解 https://blog.csdn.net/qq_42138662/article/details/109289129?spm=1001.2014.3001.5506
2. Levenberg-Marquardt 最小二乘优化 https://zhuanlan.zhihu.com/p/42415718