自学笔记：Levenberg-Marquardt的改进算法

fwalker

前言

在非线性最小二乘优化的问题中，Levenberg-Marquardt算法作为高斯牛顿法与梯度下降法中的折衷方法，已经在很多问题和场合得到应用。距笔者本身的亲身体会，高斯牛顿法虽然在理论上拥有很快的收敛速度，但在实际应用过程中常因为正规矩阵(J^TJ)^{-1}的病态问题而无法投入应用；而梯度下降法虽然拥有鲁棒性能保证收敛，但收敛速度很慢，无法满足算法的实时性要求。Levenberg-Marquardt算法通过在正规矩阵上加上一个\lambda因子而变为(J^TJ+\lambda D^TD)^{-1}来获得更好的性能，但简中网站中的科普资料通常既没给出\lambda参数的选择方法，也没给出用单位矩阵I取代D矩阵的原因，因此大部分情况下沦落为初值敏感、高度依赖调参的一种算法，使用时仍然与理想状态有很大差距。笔者最近在应用该算法时也遇到了一些问题，结合查询到的文献资料而撰写此文，希望在与广大知友的讨论交流中能有所进益。
一般非线性最小二次算法的局限性





图1 简单模型$y = e^{-\theta_1 t}+e^{\theta_2 t}$的最小二乘参数拟合的代价函数等高线图

首先，有许多个参数组成的模型通常会体现出一个叫做sloppiness(稀松)的共同特征。这个特征的意思是对于这个模型来说通常只有一些参数是强相关(stiff)的，而另一些参数的相关性比较弱(sloppy)。当算法陷入参数变化对模型行为影响不敏感的区域中，即图1所示的Plateau(高原)时，算法将参数迭代到无限值却仍然找不到合适的拟合，这个现象被称为Parameter evaporation（参数蒸发）。虽然算法处于这一空间中仍然有着满足要求的特征，e.g.\nabla C = 0，但这些值要不是局部最小值，要不是无穷远处的鞍点，一般是对数据的不良拟合，没有什么实际意义。此时通常需要手动修正来引导算法找到更好的拟合，比如说加一个罚函数。但这种修正并不是完美的，而且很有可能移动最小值所在的位置。
其次，当算法沿着狭窄的峡谷寻找最佳拟合时，算法也会变得迟缓，特别是对于10个以上的参数的问题时，峡谷的纵横比会大于1000，这要求算法沿着底部槽以非常小的步长前进。而令人头痛的是，参数蒸发和缓慢收敛的解决方案相互矛盾，提升收敛速度通常提高参数蒸发的几率，约束参数蒸发的可能性往往导致收敛速度的迟缓，这个事实加剧了参数拟合的困难。
Levenberg-Marquardt算法在面对这一困难时表现出了很好的适配性。通过随着算法的演进过程动态调制\lambda因子，该方法可以在梯度下降（避免参数蒸发）和高斯牛顿法（沿峡谷快速收敛）间进行插值，故在很多问题中表现良好。但这依赖于使用者对\lambda因子的高效调教，而且它在远离最小值处仍可能无法收敛。因此通过合适的手段对Levenberg-Marquardt算法进行改进尤为重要。
Levenberg-Marquardt算法的一般步骤

Levenberg-Marquardt算法的实现包括反复重复以下五个步骤：

• 根据当前更新的参数值更新函数和相关的雅克比矩阵
  
• 更新尺度矩阵D^TD和阻尼因子\lambda
  
• 得到参数更新步长\delta \theta，并在新的参数值\theta + \delta \theta处更新函数
  
• 根据一定判据来接受或拒绝该参数更新，若接受，继续下一步骤，若拒绝，回到步骤2
  
• 若满足所有收敛标准或超过代价函数或迭代次数的限制，退出函数
  

关于阻尼因子\lambda的选择

选择阻尼因子通常有直接或间接的两种方法。直接方法即直接更改\lambda的值，间接方法则是通过首先选择一个可接受的步长\Delta后找到一个\lambda使得|\delta \theta|\leq \Delta。Marquardt最初建议只对\lambda做出简单修改就能满足要求：如果一个步骤被接受，则\lambda缩小一个固定因数（朝高斯牛顿法靠拢）；如果一个步骤被拒绝，则\lambda适当放大若干倍。一般，J^TJ的特征值在对数尺度上有很好的间隔，因此，若将\lambda缩小或放大的因数选择为与J^TJ的特征值间距相适应的幅度，Levenberg-Marquardt算法能自然获得一个不错的效果。此外，有些研究者发现，将\lambda缩小的多而放大的少也能产生更好的效果，对于中等规模的问题，一般缩小3倍而放大2倍就够了，而对于较大规模的问题，缩小5倍而放大1.5倍的效果会更好。这种方案被称为delayed gratification（延迟满足）。该策略通过修改其乘数因子同样也可以应用于间接方法。




图2 Minpack-2项目中一个涂层厚度标准化的问题的Hessian矩阵的特征值分布，该问题由134个参数和252个残差组成。尽管该特征值横跨8个数量级，但他们在这些范围内的分布并不均匀，大部分特征值集中在$10^3$和1~10这个数量级中；故对于该问题，使用间接方法更新$\lambda$因子更为合适。

关于尺度矩阵D^TD的选择

关于尺度矩阵D，绝大多数人在使用该算法时会依着提出者Levenberg的建议，直接令D^TD=I。之后有研究者发现，令D^TD为对角矩阵，且其元素与J^TJ的对角矩阵相等能使该尺度矩阵更精确地捕获参数的缩放规律。然而，这种保持尺度不变性的方法的问题在于它大大增加了参数蒸发的敏感性。特别是，当一个参数开始蒸发时，模型对该参数的敏感度降低，因此它在J^TJ对角线的对应元素变小，进而降低了该参数的阻尼。然而为了抑制参数蒸发，阻尼应该增加，因此尺度不变的选择在此时适得其反。然而在另一方面，若算法处于峡谷区域前进，采用该尺度矩阵可以大大加快算法的前进速度。
因此有些人对对角D^TD的选择作出一定改进。它仍然令D^TD为对角矩阵，但其对角上的元素的值为历次更新的J^TJ对角元素对应值上的最大值。该方法同样保留了缩放的不变性，但仍然比单位矩阵缩放更容易出现参数蒸发现象，这是因为初值可能位于不产生足够阻尼的区域。
一个很好的折衷办法是在D^TD中指定阻尼项的最小值，这样既提高了鲁棒性，又能使算法在峡谷区域微调缩放的方向。
收敛判据

首先，必须区分收敛判据和停止判据。收敛判据是表面算法确实找到了局部极小值的判据，而停止判据则是令算法放弃继续执行的判据。关于前者，Bate & Watts给出了一种优雅的收敛准则，它缘起于最小二乘问题的几何解释。这种方法重点检测残差向量和切平面之间的角度:
\cos \phi = \frac{|P^Tr|}{|r|}
其中P^T是投影矩阵，将残差向量投影到模型对应的流形的切平面。一般投影矩阵由雅克比矩阵的SVD分解获得:J = U \Sigma V^T，则P^T = U U^T。当\cos \phi小于10^{-2}或10^{-3}数量级的某个值时，可认为算法收敛。
（这里需要对该SVD分解做出一定额外解释。在常规的SVD分解，U和V分别为m阶和n阶的酉矩阵， \Sigma 为 m\times n (这里假设m<n)的矩阵，但在此处，V和$\Sigma$均为n阶方阵，而U为 m\times n 矩阵，此时的U虽然有着近似酉矩阵的性质，但还不是酉矩阵）
该方法为一种无量纲的收敛标准，指示了算法结果与最小值的接近程度。它同时是对解的准确性的一种统计解释。然而，当魔性流行的边界很窄时，该方法存在严重缺陷。同时最佳拟合若正好具有参数蒸发的特性（通常在含噪音数据的大型模型场景中出现），那么即使算法已实质上收敛，但$\cos \phi$仍然可能很大。
若参数蒸发，模型对特定的参数组合将不会特别敏感，雅克比矩阵的奇异值会变得非常小。当它足够小时，这些参数的更新方向应位于J的零空间中。即使奇异值在形式上可能是非零的，但出于计算目的，我们知道算法不会通过在这些方向上的参数移动来获得任何进益，它也不会对残差的切向分量作出贡献。为了修正这个这个问题，需要将投影矩阵修正为P = \tilde{U} \tilde{U}^T，其中\tilde{U}为剔除了小奇异值对应的左奇异向量的奇异向量组成的矩阵。若设定算法的精度为\epsilon，则这个小奇异值的判定标准为是否小于\sqrt{\epsilon}\max{\Sigma}，其中\max{\Sigma}指代为最大的奇异值。
另一种修正的方法是当代价函数的梯度低于某个阈值时，判定为算法收敛。
除了上述收敛标准，还需要提供一个停止标准。停止标准除了与算法的最大迭代次数相关外，也应该提供何时达到最大残差和雅克比值评估次数的停止标准，或者是代价函数的梯度下降到某个阈值、阻尼项变得太大或代价函数值本身达到某个可接受的值。
测地线加速度

除了修改阻尼因子和尺度矩阵，也可以通过将每次的迭代增量修正为包含更高阶相关的项来提高算法效率。在此略去该校正的推导过程，仅呈现具体校正步骤和校正后的结果供大家参考。
令每次的迭代增量\delta \theta为一阶量和二阶量的和：
\delta \theta = \delta \theta_1 + \delta \theta_2 \equiv v\delta t+ 1/2a \delta t^2
一般来说，Levenberg-Marquardt 迭代增量表现为该非线性模型的一阶展开（该展开被称为速度项），具体为：
\delta \theta_1 = -(J^TJ+\lambda D^TD)^{-1}J^Tr
而二阶量（加速度项）为：
\delta \theta_2 \approx -\frac{1}{2}(J^TJ+\lambda D^TD)^{-1}J^Tr''
其中我们使用了方向二阶导数：r_m'' = K_{m\mu v}\delta \theta_1^\mu \delta\theta_1^v(为了避免歧义，这里使用了哑下标表示法,不了解可以翻我专栏的前一篇文章)，其中K_{m\mu v} = \frac{\partial f_m}{\partial \theta_{\mu}\partial\theta_v}。但这并不意味这需要对模型本身求完整的二阶导数，而只要计算沿着一阶校正\delta \theta_1方向的方向二阶导数，故上式可近似为：
K_{m\mu v}\delta \theta_1^\mu \delta\theta_1^v \approx \frac{2}{h}(\frac{r_m(\theta+h\delta \theta_1)-r_m(\theta)}{h}-J_{m\mu}\delta \theta_1^{\mu})
在实际使用中，一般选择h = 0.1。
此外，利用加速度校正迭代增量时，需要一个额外的条件来判断该增量是否被接受：
\frac{2|\delta \theta_2|}{|\delta \theta_1|}\leq \alpha
其中\alpha为人为设定的常量，其最优值不定。这个条件的动机是理想中的迭代增量代表了一个被截断的扰动序列，故这个值应该减小来保证收敛。对于大多数问题\alpha = 0.75已经是一个很好的设定值了，但对于那些收敛困难的问题，可以将其设定为0.1。而如果不设这个步骤，算法的结果是不可预测的，可能会出现很多大的、失控的迭代增量，导致参数蒸发。
使用测地线加速度的一个好处是，当算法迭代过程将处于狭谷区域时，可以通过抛物线获取近似路径，这样相比起原先的折线可以更准确地跟随蜿蜒路径的朝向。




图3 当处于峡谷位置朝着最佳拟合点前进时，通过用抛物线近似路径，可以用更少的迭代找到最佳拟合



图4 图1的log化后的图，如图所示，从高原区域开始迭代算法，纯LM法虽然能减少cost，但会使算法原理最佳拟合并卡在一个平坦区域，但这个点与数据的拟合度很差。引入加速度项后，算法会识别出$|\delta \theta_2|>|\delta \theta_1|$,从而增大$\lambda$来减小增量

是否“上坡”

通常来说 Levenberg-Marquardt方法的标准是接受能降低cost的迭代增量而拒绝可能导致cost增加的增量。这是一个自然并且安全的选择，但通常不是最有效的。当算法沿着一个纵横比非常大的峡谷进行拟合时，通常只有一个很小的迭代增量能降低成本，但如果允许一定程度的“上坡”，可以使算法朝着最合适的方向前进。这里我们对第i次迭代计算：
\beta_i = cos(\delta \theta_1^{new},\delta \theta_1^{old})
它表示建议的速度增量与上一次迭代被接受的速度增量的夹角，若该夹角为锐角，则接受上坡运动，随着角度的减小乃至消失，接受度增加。更具体的，如果有：
(1-\beta_i)^bC_{i+1}\leq C_i
或更保守的：
(1-\beta_i)^bC_{i+1}\leq \min(C1,...,Ci)
则接受上坡运动。其中C_i第i次迭代的cost.b的合理取值为1或2，其中2允许更大程度的上坡。
采用该方法除了能以更少的迭代次数完成优化意外，该方法对于有多个极小值的系统，能平均找到更好的极小值，即可以通过不绝对地下降到固定来避免一些沿途的坑洼处。不过该方法的高效是以牺牲稳定性为前提的，不能保证收敛。
参考文献

Transtrum M K, Sethna J P. Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization[J]. arXiv preprint arXiv:1201.5885, 2012.

本文使用 Zhihu On VSCode 创作并发布