牛顿、高斯牛顿、LM方法

牛顿法

在数值分析中，牛顿法(Newton’s method)，也称为牛顿-拉斐森法(Newton–Raphson method)是一种求解函数根的算法。

牛顿法的基本思想：从一个接近于函数根的初始值$x_0$开始，求解函数在$x_0$点处的斜率，然后根据经过$x_0$ 函数的切线与x轴的交点，切线与x轴交点的值会比初始值接近真实的根。

优点：

• 通常情况下，收敛性是二次的：当该方法在根上收敛时，每一步的根和近似值之间的差都是平方（精确数字的数量大致翻倍）；
• 收敛速度快；

缺点：

• 牛顿法要求导数可以直接计算。导数表达式在某些情况下很难直接求得或者消耗很大；
• Overshoot：如果一阶导数在特定根的邻域中表现不佳，则该方法可能会过冲，并偏离该根；
• Stationay Point：如果遇到函数的驻点，导数为零，该方法将因除以零而终止；

Newton’s Method - Department of Mathematics at UTSA

以下以求解一维函数$f(x)$ 的根为例子：

• 目标

  求解函数 $f(x) = 0$ ；

• 推导过程

  $f:(a,b)$ 在区间$(a,b)$ 中连续可导，当前有一个近似解$x_n$ 。函数$y=f(x)$ 在 $x=x_n$ 处的切线函数为
  $$
  y =f’(x_n)(x-x_n)+f(x_n)
  $$
  $f’$ 表示函数的一阶导数，将$y$ 的解($y=0$ 时候的$x$ 值)作为下一个近似值$x_{n+1}$，$x_{n+1}$ 更接近于函数的解。可得如下推导：
  $$
  0=y =f’(x_n)(x-x_n)+f(x_n)
  $$
  则可以求解 $x_{n+1}$,
  $$
  x_{n+1} = x_n -\frac{f(x_n)}{f’(x_n)}
  $$

已经证明，如果$f’$是连续的，并且待求的$x$是孤立的，那么在零点$x$周围存在一个区域，只要初始值$x_0$位于这个邻近区域内，那么牛顿法必定收敛。并且，如果 $f’$ 不为0, 那么牛顿法将具有平方收敛的性能。粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。由于牛顿法是基于当前位置的切线来确定下一次的位置，所以牛顿法又被很形象地称为是"切线法"。

牛顿法应用于一元函数求根算法过程:

1. 初始化参数向量${x_0}$；
2. 在每次迭代中，计算函数 $ f(x)$ 和其相对于 $ x $ 的梯度（一阶导数）；
3. 求解当前迭代中参数向量的变化量 $\Delta x $ ；
4. 更新参数向量：$ x_{n+1} = x_n + \Delta x$
5. 重复步骤2-4，直到满足收敛条件。

data-science-notebook/ml/BAT_机器学习_1000_题/301-400.md at master · wizardforcel/data-science-notebook (github.com)

牛顿方法应用于最小二乘问题

考虑如下最小二乘问题：
$$
\begin{aligned}
&F(\mathbf{x})=\frac{1}{2}\sum_{i=1}^{m}\left(f_i(\mathbf{x})\right)^2
\end{aligned}
$$
其中，待拟合函数 $ f_i:\mathbf{R}^n\mapsto\mathbf{R},i=1,\ldots,m \quad m\geq n.$ 。找出使得损失函数$F(\mathbf{x})$ 最小时的$x$的值$\mathbf{x}^*$。

在最优化的问题中，线性最优化可以直接求解，但对于非线性优化问题，牛顿法提供了一种求解的办法**。任务是优化一个目标函数 $F(x)$，求函数$F$ 的极小（极大）问题，可以转化为求解函数 $F$ 的导数 $F^{\prime}(\mathbf{x}) = 0$ 的问题，这样求可以把优化问题看成方程求解问题。**

$F^{\prime}(\mathbf{x+h}) $ 的泰勒展开式为：
$$
\begin{aligned}
\mathbf{F}^{\prime}(\mathbf{x}+\mathbf{h}) & \mathbf{=}\mathbf{F}^{\prime}(\mathbf{x})+\mathbf{F}^{\prime\prime}(\mathbf{x})\mathbf{h}+O(|\mathbf{h}|^{2}) \
&\simeq\mathbf{F}^{\prime}(\mathbf{x})+\mathbf{F}^{\prime\prime}(\mathbf{x})\mathbf{h}\quad\mathrm{for}\left|\mathbf{h}\right|\text{sufficiently small}
\end{aligned}
$$
根据牛顿法可得：
$$
\mathrm{Hh}_\mathrm{n} = -\mathbf{F}^{\prime}(\mathbf{x})
$$

其中 $F^{\prime}(\mathbf{x})$ 为 $F(\mathbf{x})$ 的梯度 , ${\textrm{H}}={\textbf{F}}^{\prime\prime}({\textbf{x}})$ 为Hassian矩阵，
$$
\left.\mathbf{F}‘(\mathbf{x}) = \left[\begin{array}{c}\displaystyle\frac{\partial F}{\partial x_1}(\mathbf{x})\\
\vdots \\
\displaystyle\frac{\partial F}{\partial x_n}(\mathbf{x})\end{array}\right.\right] \ \
\mathbf{F}’'(\mathbf{x})=\left[\frac{\partial^2F}{\partial x_i\partial x_j}(\mathbf{x})\right].
$$
下一次迭代过程中：
$$
\mathbf{x}:=\mathbf{x}+\mathbf{h_{n}}.
$$
牛顿法应用于最小二乘问题算法过程：

1. 初始化参数向量$\mathbf{x_0}$；
2. 在每次迭代中，计算损失函数 $F(\mathbf{x})$ 相对于参数向量 $\mathbf{x}$ 的梯度（一阶导数）和 Hessian 矩阵（二阶导数）；
3. 求解当前迭代中参数向量的变化量 $\mathbf{h_{n}} $ ；
4. 更新参数向量：$\mathbf{x_{n+1}} =\mathbf{x_n}+\mathbf{h_{n}}.$
5. 重复步骤2-4，直到满足收敛条件。

最优化理论与算法------牛顿法（附Matlab实现）： - YuhuaStone - 博客园 (cnblogs.com)

imm3215.pdf (dtu.dk)

高斯牛顿法(Gaussian-Newton’s Method)

牛顿法（Newton’s Method）是一种用于求解数值优化问题的迭代方法，特别适用于多维场景，例如优化多元函数或求解非线性方程组。它基于泰勒级数展开，通过不断迭代来逼近函数的局部最优解。

对于多维函数，使用牛顿法进行优化时需要计算Hessian矩阵，该矩阵是一个对称矩阵，因此需要计算n*n/2次二阶偏导数，计算量相当大，所以人们为了简化计算，在牛顿法的基础上，将其发展为高斯-牛顿法。高斯牛顿法和牛顿法的区别在于，高斯牛顿法是对回归模型$f(x)$ 进行泰勒展开近似替代$f(x)$，然后通过多次迭代，多次修正回归系数，使回归系数不断逼近回归模型的最佳回归系数，最后使原模型的残差平方和达到最小。

在$x$ 点附近对 $f(x)$ 进行泰勒展开取一阶线性项近似：
$$
f\left(x+\Delta x\right)\approx f\left(x\right)+J\left(x\right)\Delta x
$$
将其带入优化目标函数：
$$
\frac{1}{2} || f(x + \Delta x) ||^2 = \frac{1}{2}(f(x)^Tf(x) + 2f(x)^TJ(x) \Delta x + \Delta x^TJ(x)^TJ(x) \Delta x)
$$
上式对$\Delta x$ 求导，并令导数为0可得：
$$
J(x)^TJ(x) \Delta x = -J(x)^Tf(x)
$$
令$H = J^TJ \quad B = -J^Tf$ 可得：
$$
H \Delta x = B
$$
$J$ 为雅各比矩阵，求解上式，便可以获得调整增量$\Delta x$。

注意,我们要求解的变量是 $\Delta x$,因此这是一个线性方程组,我们称它为增量方程,也可以称为高斯牛顿方程 (Gauss Newton equations) 或者正规方程 (Normal equations)。

这里把左侧记作 $H$ 是有意义的。对比牛顿法可见, Gauss-Newton 用 $J^TJ$ 作为牛顿法中二阶 Hessian 矩阵的近似, 从而省略了计算 $H$ 的过程。求解增量方程是整个优化问题的核心所在。

Gauss-Newton 的算法步骤可以写成:

1. 给定初始化值$x_0$ ;
2. 对第k次迭代求解雅各比矩阵$\mathbf{J}(x_k)$ 和$f(x_k)$；
3. 求解增量式方程 $H \Delta x_k= B$；
4. 更新$x_{k+1}=x_{k}+\Delta x_{k}$ ；
5. 判断$\Delta x_{k}$ 是否足够小，如果不满足则重复2-4。

求解增量方程时，原则上，要求我们所用的近似 H 矩阵是可逆的(而且是正定的)，但实际数据中计算得到的 $H = J^TJ $ 却只有半正定性。也就是说，在使用 Gauss Newton 方法时，可能出现 JTJ为奇异矩阵或者病态 (illcondition) 的情况，此时增量的稳定性较差，导致算法不收敛。更严重的是，就算我们假设 H 非奇异也非病态，如果我们求出来的步长 $\Delta x_{k}$ 太大，也会导致我们采用的局部近似不够准确，这样一来我们甚至都无法保证它的迭代收敛，哪怕是让目标函数变得更大都是有可能的。

牛顿法和高斯牛顿法对比

视觉slam14讲

LM(Levenberg–Marquardt)法

由于 Gauss-Newton 方法中采用的近似泰勒展开只能在展开点附近有较好的近似效果，通过信赖区域 (Trust Region)的方法，在信赖区域里边，认为近似是有效的; 出了这个区域，近似可能会出问题。

那么如何确定这个信赖区域的范围呢？一个比较好的方法是根据我们的近似模型跟实际函数之间的差异来确定这个范围：如果差异小，就让范围尽可能大；如果差异大，就缩小这个近似范围。具体数学表达式为：
$$
\rho=\frac{f\left(\boldsymbol{x}+\Delta\boldsymbol{x}\right)-f\left(\boldsymbol{x}\right)}{\boldsymbol{J}\left(\boldsymbol{x}\right)\Delta\boldsymbol{x}}.
$$

• $\rho$ 接近1，近似是好的；
• $\rho$ 太小，实际减小的值远少于近似减小的值，认为近似比较差，需要缩小近似范围；
• $\rho$ 太大，实际减小的值远大于近似减小的值，需要增大近似范围；

则L-M算法框架

1. 给定初始值 $x_0$, 以及初始优化半径 $\mu$。
2. 对于第 $k$ 次迭代, 求解:

$$
\min_{\Delta\boldsymbol{x}_k}\frac{1}{2}|f\left(\boldsymbol{x}_k\right)+\boldsymbol{J}\left(\boldsymbol{x}_k\right)\Delta\boldsymbol{x}_k|^2,\quad s.t.|\boldsymbol{D}\Delta\boldsymbol{x}_k|^2\leq\mu,
$$

​ 这里 $\mu$ 是信赖区域的半径。

3. 计算 $\rho$。
4. 若$\rho>\frac{3}{4},$ 则$\mu=2\mu;$
5. 若$\rho<\frac{1}{4},:\text{则}\mu=0.5\mu;$
6. 如果 $\rho$ 大于某阈值, 认为近似可行。令 $x_{k+1}=x_k+\Delta x_k$。
7. 判断算法是否收敛。如不收敛则返回 2, 否则结束。

近似范围扩大的倍数和阈值都是经验值。在 Levenberg 提出的优化方法中，把 $D $取成单位阵 $I$，相当于直接把 $\Delta x$ 约束在一个球中。随后，Marqaurdt 提出将 $D $ 取成非负数对角阵，实际中通常用$ J ^T J $的对角元素平方根，使得在梯度小的维度上约束范围更大一些。不论如何，在 L-M 优化中，我们都需要解一个子问题来获得梯度。这个子问题是带不等式约束的优化问题，我们用 Lagrange 乘子将它转化为一个无约束优化问题:
$$
\min\limits_{\Delta\boldsymbol{x}_k}\frac{1}{2}\left|f\left(\boldsymbol{x}_k\right)+\boldsymbol{J}\left(\boldsymbol{x}_k\right)\Delta\boldsymbol{x}_k\right|^2+\frac{\lambda}{2}\left|D\Delta\boldsymbol{x}\right|^2.
$$
$\lambda$ 为Lagrange 乘子，类似的可得计算增量的线性方程：
$$
\left(H+\lambda D^TD\right)\Delta x=g
$$
考虑其简化形式，将$D=I$，
$$
(\boldsymbol{H}+\lambda\boldsymbol{I})\Delta x=\boldsymbol{g}
$$
从上式子看到,当参数 $\lambda$ 比较小时, $H$ 占主要地位,这说明二次近似模型在该范围内是比较好的, L-M 方法更接近于高斯牛顿法。另一方面,当 $\lambda$ 比较大时, $\lambda\boldsymbol{I} $ 占据主要地位, L-M更接近于一阶梯度下降法(即最速下降),这说明附近的二次近似不够好。L-M 的求解方式,可在一定程度上避免线性方程组的系数矩阵的非奇异和病态问题, 提供更稳定更准确的增量 $\Delta x$。

LMA.pdf (ucsb.edu)

imm3215.pdf (dtu.dk)