数值计算之 Levenberg-Marquardt算法

前言

本篇是牛顿法的最后一篇，Levenberg-Marquardt算法，也就是阻尼牛顿法。

高斯牛顿法

回顾上篇中的高斯牛顿法： f ( x 0 + Δ x ) = f ( x 0 ) + J ( x 0 ) Δ x min ⁡ Δ x 1 2 ∣ ∣ f ( x ) ∣ ∣ 2 = min ⁡ Δ x 1 2 ( f ( x 0 ) + J ( x 0 ) Δ x ) T ( f ( x 0 ) + J ( x 0 ) Δ x ) J T J Δ x = − J T f f({\bf x_0+\Delta x}) = f({\bf x_0}) + J({\bf x_0}){\bf \Delta x} \\ \quad \\ \min_{\bf \Delta x} \frac{1}{2} ||f({\bf x})||^2 \\ = \min_{\bf \Delta x} \frac{1}{2} (f({\bf x_0})+J({\bf x_0}){\bf \Delta x})^T(f({\bf x_0})+J({\bf x_0}){\bf \Delta x}) \\ \quad \\ J^TJ{\bf \Delta x} = -J^Tf f(x0​+Δx)=f(x0​)+J(x0​)ΔxΔxmin​21​∣∣f(x)∣∣2=Δxmin​21​(f(x0​)+J(x0​)Δx)T(f(x0​)+J(x0​)Δx)JTJΔx=−JTf 将牛顿法中的海森矩阵 H H H，使用一阶梯度 J T J J^TJ JTJ进行了替换，因而避免了计算二阶梯度。然而， J T J J^TJ JTJ并非正定，可能是奇异阵或者病态矩阵，求逆不稳定，使得迭代的稳定性和收敛性差。

LM算法 阻尼牛顿法

为了提升迭代的稳定性，在优化函数上添加阻尼（机器学习里有时称为正则项），来降低迭代步长： min ⁡ Δ x ( 1 2 ∣ ∣ f ( x ) ∣ ∣ 2 + 1 2 λ Δ x T Δ x ) = min ⁡ Δ x 1 2 ( f ( x 0 ) + J ( x 0 ) Δ x ) T ( f ( x 0 ) + J ( x 0 ) Δ x ) + 1 2 λ Δ x T Δ x ) \min_{\bf \Delta x} (\frac{1}{2} ||f({\bf x})||^2+\frac{1}{2}\lambda{\bf \Delta x^T\Delta x}) \\ = \min_{\bf \Delta x} \frac{1}{2} (f({\bf x_0})+J({\bf x_0}){\bf \Delta x})^T(f({\bf x_0})+J({\bf x_0}){\bf \Delta x}) + \frac{1}{2} \lambda {\bf \Delta x^T\Delta x}) \\ Δxmin​(21​∣∣f(x)∣∣2+21​λΔxTΔx)=Δxmin​21​(f(x0​)+J(x0​)Δx)T(f(x0​)+J(x0​)Δx)+21​λΔxTΔx) 在 1 2 ∣ ∣ f ( x ) ∣ ∣ 2 \frac{1}{2}||f({\bf x})||^2 21​∣∣f(x)∣∣2后添加 1 2 λ Δ x T Δ x \frac {1}{2}\lambda {\bf \Delta x}^T {\bf \Delta x} 21​λΔxTΔx，当求出的 Δ x \bf \Delta x Δx比较大时，对函数值进行惩罚。

还是按照高斯牛顿法一样求一阶梯度等于零，获得迭代表达式： J T J Δ x + λ Δ x = − J T f ( J T J + λ I ) Δ x = − J T f Δ x = − ( J T J + λ I ) − 1 J T f J^TJ{\bf \Delta x} +\lambda {\bf \Delta x} = -J^Tf \\ (J^TJ+\lambda I){\bf \Delta x}=-J^Tf \\ {\bf \Delta x} = -(J^TJ+\lambda I)^{-1}J^Tf \\ JTJΔx+λΔx=−JTf(JTJ+λI)Δx=−JTfΔx=−(JTJ+λI)−1JTf 可以看出，阻尼牛顿法通过 J T J + λ I J^TJ+\lambda I JTJ+λI替换 H H H，一方面避免二阶梯度计算，另一方面保证了矩阵的可逆性。

阻尼系数 λ \lambda λ

上面的惩罚项中， λ \lambda λ用于控制 Δ x \bf \Delta x Δx对迭代步长的影响。当 f ( x 0 + Δ x ) f(\bf x_0 +\Delta x) f(x0​+Δx)与一阶泰勒展开 f ( x 0 ) + J ( x 0 ) Δ x f({\bf x_0}) +J({\bf x_0}){\bf\Delta x} f(x0​)+J(x0​)Δx足够近似，迭代步长可以大一些，加速收敛；如果函数值与泰勒展开不够近似，迭代步长要小一些，增加稳定。

因此，可以设计以下策略调整阻尼系数： ρ = f ( x 0 + Δ x ) − f ( x 0 ) J ( x 0 ) Δ x i f ρ > 3 4 : λ = 1 2 λ i f ρ < 1 2 : λ = 2 λ \rho = \frac {f({\bf x_0 + \Delta x})-f({\bf x_0})} {J({\bf x_0}){\bf \Delta x}} \\ \quad \\ if \quad \rho > \frac{3}{4}: \quad \lambda=\frac{1}{2}\lambda \\ \quad \\ if \quad \rho < \frac{1}{2}: \quad \lambda= 2\lambda ρ=J(x0​)Δxf(x0​+Δx)−f(x0​)​ifρ>43​:λ=21​λifρ 3 4 : d = 2 d i f ρ < 1 2 : d = 1 2 d \rho = \frac {f({\bf x_0 + \Delta x})-f({\bf x_0})} {J({\bf x_0}){\bf \Delta x_0}} \\ \quad \\ if \quad \rho > \frac{3}{4}: \quad d=2d \\ \quad \\ if \quad \rho < \frac{1}{2}: \quad d=\frac{1}{2}d ρ=J(x0​)Δx0​f(x0​+Δx)−f(x0​)​ifρ>43​:d=2difρ