★ FEATURED ARTICLE 数值优化这门手艺做久了会发现一个规律真正天天被调用的算法往往不是教科书里最优雅的那个而是最“抗造”的那个。Levenberg-Marquardt 算法下称 LM 算法就是典型代表。它既不是梯度下降的纯粹版也不是高斯-牛顿的纯粹版而是把两者揉在一起用一个阻尼因子在“稳”和“快”之间来回切换。搞计算机视觉的用它做光束法平差搞机器人的用它做位姿图优化搞曲线拟合的用它把一堆散点压成一条光滑曲线搞工程反演的用它从观测数据里反推材料参数。可以说只要你的问题最后归结为“找一组参数让模型输出和实测数据的残差平方和最小”LM 算法基本都能上场。这篇文章面向的是真正要动手写代码、调参数、排查不收敛问题的人。我会从设计思路讲起把阻尼因子怎么起作用、雅可比矩阵怎么组织、迭代什么时候该停这些关键点掰开揉碎然后给出一套可以直接照着写的实现流程和参数取值参考最后把我这些年踩过的坑整理成一张速查表。全文不依赖任何特定框架伪代码和 Python 示例都有数学基础只需要你懂矩阵乘法和最小二乘的基本概念。读完之后你应该能独立实现一个能跑、能收敛、能解释每一步在干什么的 LM 求解器。1. 算法整体设计与思路拆解1.1 为什么要在梯度下降和高斯-牛顿之间“骑墙”要理解 LM 算法得先明白它想解决谁的痛点。假设我们要最小化一个残差平方和函数 $F(\mathbf{x}) \frac{1}{2}|\mathbf{r}(\mathbf{x})|^2$其中 $\mathbf{r}$ 是残差向量$\mathbf{x}$ 是待求参数。梯度下降的做法是沿着负梯度方向走一小步优点是无论从哪儿出发都大概率能往下走缺点是接近最优点时步子越来越小收敛慢得让人抓狂而且对参数尺度很敏感——一个参数是毫米级、另一个是万级同一个学习率根本没法用。高斯-牛顿则是另一套思路把残差在当前点做一阶线性化直接解一个线性最小二乘问题求出步长。它在最优点附近收敛极快二阶收敛不是梦但前提是初始点得离最优解足够近而且雅可比矩阵得满秩、条件数不能太差。一旦初始点选偏了或者模型本身有强非线性高斯-牛顿给出的步长可能大得离谱直接把迭代带到更差的点上甚至发散。LM 算法的聪明之处在于它不选边站而是引入一个阻尼因子 $\lambda$把高斯-牛顿的法方程改造成$$(\mathbf{J}^T\mathbf{J} \lambda \mathbf{I})\Delta\mathbf{x} -\mathbf{J}^T\mathbf{r}$$这里的 $\mathbf{J}$ 是残差对参数的雅可比矩阵。当 $\lambda \to 0$ 时式子退化成纯高斯-牛顿当 $\lambda$ 很大时$\lambda \mathbf{I}$ 主导步长方向趋近于负梯度方向严格说是梯度除以 $\lambda$步长被压得很小。所以 $\lambda$ 本质上是一个“信任调节旋钮”模型可信、当前点不错就把 $\lambda$ 调小放开步子走模型不可信、刚走了一步反而变差了就把 $\lambda$ 调大缩回梯度下降的保守姿态。这个设计背后的逻辑很朴素我们没法事先知道当前点离最优解多远那就让算法自己根据每步的实际效果来动态判断。走得好就大胆走得差就谨慎。这种“用结果反馈调节策略”的思路比固定学习率或者固定阻尼要稳健得多也是 LM 算法几十年来一直是非线性最小二乘默认选择的核心原因。1.2 阻尼项为什么加在 $\mathbf{J}^T\mathbf{J}$ 上而不是别处有人会问阻尼为什么不直接加在步长上或者加在残差上加在 $\mathbf{J}^T\mathbf{J}$ 的对角线上有几个很实在的好处。第一$\mathbf{J}^T\mathbf{J}$ 本身是半正定的加上 $\lambda \mathbf{I}$$\lambda 0$之后一定是正定的这保证了线性方程组有唯一解数值上也更稳定。第二这个形式等价于在每一步求解一个带正则的最小二乘问题正则项惩罚步长过大天然抑制了高斯-牛顿那种“一步冲过头”的毛病。第三当某个参数在当前点几乎不影响残差时$\mathbf{J}^T\mathbf{J}$ 对应行列接近零矩阵接近奇异此时阻尼项正好把这个方向“托住”避免求解器报奇异矩阵错误。还有一个细节值得说严格意义上的 LM 用的是 $\lambda \cdot \text{diag}(\mathbf{J}^T\mathbf{J})$ 而不是 $\lambda \mathbf{I}$也就是按每个参数自身的敏感度来缩放阻尼。这个变体叫“尺度自适应阻尼”在参数量纲差异大的问题里效果明显更好。举个例子你同时拟合一个长度参数量级 $10^{-3}$和一个温度参数量级 $10^2$如果统一加 $\lambda \mathbf{I}$阻尼对长度参数的相对影响会被放大几万倍导致那个方向几乎不动。用 $\text{diag}(\mathbf{J}^T\mathbf{J})$ 做缩放每个参数受到的阻尼强度就和它自身的敏感度匹配了。我在实际项目里默认用这个变体除非问题规模特别小、参数尺度统一才会退回简单的 $\lambda \mathbf{I}$。1.3 收敛判据的组合设计迭代什么时候停这件事比想象中重要。只用一个判据往往会出问题只看步长 $|\Delta\mathbf{x}|$ 小于阈值可能在参数几乎不动但残差还很大的平坦区域误判收敛只看残差变化量可能在残差下降极慢的长尾阶段空转几千次只看梯度范数又可能因为尺度问题设不好阈值。稳妥的做法是组合多个判据任意一个满足就停同时设一个最大迭代次数兜底。常用的组合是梯度范数 $|\mathbf{J}^T\mathbf{r}|_\infty$ 小于设定容差说明已经到驻点附近、步长相对变化 $|\Delta\mathbf{x}|/(|\mathbf{x}|\epsilon)$ 小于容差说明参数不再显著变化、残差相对下降量小于容差说明目标函数已经榨不出更多改进。三个判据的阈值一般取 $10^{-8}$ 到 $10^{-12}$ 量级具体看你的残差量级和精度需求。我习惯把梯度判据设得最严因为它最能反映“是否真的到了极值点”另外两个作为辅助防止在病态问题上无限循环。2. 核心细节解析与实操要点2.1 雅可比矩阵的组织方式与稀疏性利用LM 算法每一步的计算大头都在构造雅可比矩阵 $\mathbf{J}$ 和解线性方程组。$\mathbf{J}$ 是 $m \times n$ 的矩阵$m$ 是残差个数$n$ 是参数个数。对于小规模问题$n$ 几十到几百直接稠密存储、用 Cholesky 分解解 $\mathbf{J}^T\mathbf{J} \lambda\mathbf{I}$ 就够了。但很多实际问题的 $n$ 能到几千甚至上万比如大规模光束法平差里每张影像的位姿加每个三维点参数量轻松破万这时候稠密矩阵根本存不下。关键洞察是$\mathbf{J}$ 通常是稀疏的因为每个残差只依赖于少数几个参数。比如一个重投影误差只涉及一个相机位姿和一个三维点那这一行雅可比里非零元素就那么几个。利用这个结构$\mathbf{J}^T\mathbf{J}$ 也是稀疏的可以用稀疏 Cholesky 或者共轭梯度类迭代解法来解。工程上常见的做法是先分析问题的关联结构把参数分块用 Schur 补消元把其中一块消掉再解剩下的块。这套路子在视觉 SLAM 里是标配但原理对所有稀疏非线性最小二乘都通用。实操中还有一个容易忽略的点雅可比的计算方式。解析求导当然最准最快但模型复杂时手推导数容易出错而且改模型就得改导数代码维护成本高。自动微分现在很成熟对于中等规模问题完全可用代价是每步计算量增加但换来的是正确性和开发效率。我的经验是模型定型前用自动微分快速迭代模型稳定后如果性能吃紧再把热点部分的雅可比换成解析形式。数值差分有限差分尽量别用步长选不好误差大而且每列要多算一次残差$n$ 大时开销惊人。2.2 阻尼因子的更新策略Marquardt 规则与 Nielsen 改进$\lambda$ 怎么变是 LM 算法里最需要拿捏的地方。最经典的 Marquardt 规则很直接算出一个候选步长后比较新旧残差。如果新残差更小接受这一步并把 $\lambda$ 除以一个因子比如 10表示“这步走得不错下次更大胆”如果新残差更大拒绝这一步把 $\lambda$ 乘以一个因子比如 10重新解方程再试直到走出一个能降低残差的步或者 $\lambda$ 大到步长小到可以忽略。这个规则简单有效但有两个毛病。一是乘除因子固定问题尺度不同时表现差异大二是 $\lambda$ 的调整是“跳跃式”的可能在一个本来合适的值附近来回震荡。Nielsen 的改进版用了一个更平滑的更新公式引入增益比 $\rho$$$\rho \frac{F(\mathbf{x}) - F(\mathbf{x} \Delta\mathbf{x})}{ \text{预测下降量} }$$预测下降量用线性化模型算出来是 $\Delta\mathbf{x}^T(\lambda \Delta\mathbf{x} - \mathbf{J}^T\mathbf{r})$ 这类形式。$\rho$ 接近 1 说明线性化模型预测得很准可以大幅减小 $\lambda$$\rho$ 很小甚至为负说明模型完全不可信要大幅增大 $\lambda$。Nielsen 还给出了 $\lambda$ 更新的具体系数涉及 $\nu$ 参数初始取 2根据 $\rho$ 调整实测下来比固定因子版本收敛更稳尤其在残差曲面比较扭曲的问题上。我自己的习惯是先用经典 Marquardt 规则把流程跑通确认雅可比和残差计算没问题再换成 Nielsen 版本做精细调优。因为 Nielsen 版本对 $\rho$ 的计算依赖线性化预测如果雅可比本身有 bug$\rho$ 会给出误导性的反馈反而掩盖问题。2.3 参数尺度归一化被低估的收敛关键前面提过用 $\text{diag}(\mathbf{J}^T\mathbf{J})$ 做阻尼缩放但这只解决了一部分问题。如果参数本身的物理量纲差异巨大比如一个参数在 $10^{-6}$ 量级、另一个在 $10^{6}$ 量级即使有阻尼缩放数值精度和步长控制还是会出问题。更彻底的做法是在优化前对参数做归一化把每个参数除以其典型尺度让所有参数在优化过程中都处于 $O(1)$ 量级优化完再乘回去。典型尺度怎么定如果参数有物理意义用它的预期量级如果没有用初始值的绝对值加个下限防止除零或者用雅可比列范数估计。这个预处理步骤看起来不起眼但我见过太多“算法不收敛”的案例最后发现就是参数尺度没处理好。一个直观的判断方法如果你发现 $\mathbf{J}^T\mathbf{J}$ 的条件数超过 $10^{10}$基本可以确定需要做尺度归一化或者检查模型参数化方式。注意归一化会改变雅可比的值如果你在优化过程中还要用雅可比做其他分析比如不确定性估计记得把归一化的影响换算回去否则协方差矩阵的量纲会错。3. 实操过程与核心环节实现3.1 从零实现一个可用的 LM 求解器下面给一个结构清晰的 Python 实现不依赖 scipy 的 least_squares目的是让你看清每一步。实际项目里当然可以直接调库但自己写一遍对排查问题极有帮助。import numpy as np def lm_optimize(residual_fn, jacobian_fn, x0, max_iter200, lambda_init1e-3, tol_grad1e-10, tol_step1e-10, tol_cost1e-12): x np.array(x0, dtypefloat) n x.size lam lambda_init r residual_fn(x) cost 0.5 * np.dot(r, r) J jacobian_fn(x) for it in range(max_iter): # 构造法方程 JtJ J.T J Jtr J.T r # 尺度自适应阻尼 diag np.diag(JtJ).copy() diag[diag 1e-12] 1e-12 A JtJ lam * np.diag(diag) # 解线性方程组 try: delta np.linalg.solve(A, -Jtr) except np.linalg.LinAlgError: lam * 10.0 continue # 试探步 x_new x delta r_new residual_fn(x_new) cost_new 0.5 * np.dot(r_new, r_new) # 增益比 predicted -0.5 * delta (lam * diag * delta Jtr) rho (cost - cost_new) / (predicted 1e-30) if cost_new cost: # 接受步 x x_new r r_new cost cost_new J jacobian_fn(x) lam * max(1.0/3.0, 1.0 - (2.0*rho - 1.0)**3) lam max(lam, 1e-12) else: # 拒绝步 lam * 2.0 if rho 0.25 else 10.0 # 收敛判据 grad_inf np.max(np.abs(Jtr)) step_rel np.linalg.norm(delta) / (np.linalg.norm(x) 1e-12) cost_rel abs(cost - cost_new) / (cost 1e-12) if grad_inf tol_grad or step_rel tol_step or cost_rel tol_cost: break return x, cost, it这段代码有几个地方值得展开。第一predicted的计算用的是线性化模型的预测下降量公式里那个lam * diag * delta项对应阻尼带来的额外惩罚少了它 $\rho$ 会算偏。第二接受步之后 $\lambda$ 的更新用了 Nielsen 的立方公式max(1/3, ...)保证 $\lambda$ 不会一次降得太狠。第三拒绝步时根据 $\rho$ 分两档增大 $\lambda$$\rho$ 为负说明方向都错了要更激进地增大。第四收敛判据里Jtr用的是接受步之前的梯度严格说应该在更新后重算但为了省一次雅可比计算用旧值判断在大多数情况下够用精度要求高的话可以调整。3.2 一个曲线拟合的完整案例拿一个指数衰减加正弦的模型练手$y a e^{-bt} c \sin(dt)$参数 $\mathbf{x} (a, b, c, d)$。生成带噪声的合成数据从偏离真值的初始点出发看 LM 能不能收敛回去。def model(x, t): a, b, c, d x return a * np.exp(-b * t) c * np.sin(d * t) def residual_fn(x): return model(x, t_data) - y_data def jacobian_fn(x): a, b, c, d x e np.exp(-b * t_data) s np.sin(d * t_data) co np.cos(d * t_data) J np.zeros((t_data.size, 4)) J[:, 0] e J[:, 1] -a * t_data * e J[:, 2] s J[:, 3] c * t_data * co return J真值取 $a2.0, b0.5, c1.0, d3.0$初始点故意设成 $a1.0, b1.0, c0.5, d2.0$噪声标准差 0.05。跑下来一般 20 到 40 次迭代收敛最终参数误差在噪声水平内。这个案例的价值在于$b$ 和 $d$ 分别控制衰减和振荡频率两者在数据覆盖不充分时容易耦合是检验阻尼策略是否合理的好场景。如果你发现迭代在 $b$ 和 $d$ 之间来回拉扯、残差下降很慢多半是 $\lambda$ 的初始值或更新因子需要调。参数取值上lambda_init我一般取 $10^{-3}$问题非线性强的话取 $10^{-2}$ 甚至 $10^{-1}$让算法一开始就保守些。max_iter设 200 是保险值正常问题几十次就收敛了如果跑满 200 次还没停基本可以断定有问题要么雅可比错了要么模型和数据的尺度不匹配。3.3 大规模稀疏问题的处理路径当参数上千稠密解法就顶不住了。这时候要做的第一件事是分析雅可比的稀疏结构。以位姿图优化为例残差是每条边上的相对位姿误差只涉及两个节点的参数所以雅可比每行只有两个非零块。$\mathbf{J}^T\mathbf{J}$ 的结构对应图的邻接结构可以用稀疏矩阵存储解方程用稀疏 Cholesky比如 CHOLMOD 这类库或者预条件共轭梯度。如果问题有“块状”结构比如同时优化一批局部参数和一批全局参数Schur 补消元是标准操作把 $\mathbf{J}^T\mathbf{J} \lambda \text{diag}$ 按参数分块写成 $\begin{pmatrix} A B \ B^T C \end{pmatrix}$先消掉 $C$ 块通常是对角块容易求逆得到关于 $A$ 的 Schur 补 $A - B C^{-1} B^T$解出第一块参数再回代求第二块。这个技巧能把一个巨大的稠密问题拆成若干小问题计算量降一个量级。实现上要注意 Schur 补矩阵可能变稠密需要根据具体结构决定消元顺序。提示稀疏实现里$\lambda$ 的尺度自适应阻尼要作用在 $\mathbf{J}^T\mathbf{J}$ 的对角线上稀疏矩阵取对角线元素有专门接口别不小心转成稠密再取内存会爆。4. 常见问题与排查技巧实录4.1 迭代不收敛的典型症状与对策LM 不收敛的表现有好几种对症下药才有效。第一种是残差一开始下降然后停滞梯度范数降不下去。这通常是模型本身在最优解附近不可辨识比如两个参数效果几乎相同共线性或者某个参数对残差几乎没影响。排查方法是看 $\mathbf{J}^T\mathbf{J}$ 的特征值谱如果有接近零的特征值对应的特征向量就指出了不可辨识的方向。对策要么改模型参数化去掉冗余参数要么加更强的正则增大 $\lambda$ 的下限要么固定那个不敏感的参数。第二种是残差震荡接受步和拒绝步交替出现$\lambda$ 忽大忽小。这多半是 $\lambda$ 的更新因子太激进或者增益比计算有误。检查predicted的公式有没有漏项把 $\lambda$ 的乘除因子从 10 改成 3 到 5 试试通常能缓解。第三种是直接发散残差越来越大。先查雅可比是不是符号错了或者漏了项——用有限差分对一遍解析雅可比误差应该在 $10^{-6}$ 相对量级以内。再查残差函数的输入输出维度对不对参数顺序有没有搞混。这两类低级错误占了发散案例的一大半。症状可能原因排查动作对策残差停滞不降参数不可辨识/共线性看 $\mathbf{J}^T\mathbf{J}$ 最小特征值改参数化或加正则接受拒绝交替震荡$\lambda$ 更新过激/增益比算错打印 $\rho$ 序列减小更新因子核对预测下降公式残差持续增大雅可比错误/维度错位有限差分对比雅可比修正导数或维度步长极小但梯度大参数尺度失衡看 $\mathbf{J}^T\mathbf{J}$ 条件数参数归一化迭代次数打满收敛判据过严打印各判据实际值放宽阈值或检查残差量级4.2 雅可比验证的实操手法解析雅可比写错是最高频的 bug而且往往错得很隐蔽——程序能跑就是收敛慢或者收敛到错的地方。最可靠的验证手段是中心差分对比对每个参数 $x_i$取步长 $h 10^{-6} \cdot \max(1, |x_i|)$算 $(r(x h e_i) - r(x - h e_i)) / (2h)$和解析雅可比的第 $i$ 列逐元素比。相对误差超过 $10^{-4}$ 就要警惕超过 $10^{-2}$ 基本可以确定有错。这个检查要放在单元测试里每次改模型都跑一遍。我吃过一次亏模型里加了一个平滑项忘了同步更新雅可比结果优化结果总是差那么一点点查了两天才发现。从那以后雅可比对比成了我所有非线性最小二乘代码的标配测试。4.3 收敛精度的现实预期很多人对 LM 的收敛精度有不切实际的期待。理论上高斯-牛顿在最优解附近二阶收敛但那是理想情况。实际数据有噪声、模型有偏差残差不可能降到零参数估计的精度受限于数据的信噪比不是算法能突破的。判断是否“收敛到位”的实用标准是看参数估计的协方差 $\sigma^2 (\mathbf{J}^T\mathbf{J})^{-1}$$\sigma^2$ 是残差方差估计如果连续几步的参数变化已经远小于这个标准差再迭代也没意义了。另外浮点精度也是天花板。残差平方和降到 $10^{-16}$ 相对量级附近时数值误差开始主导继续迭代可能反而让结果变差。所以收敛判据里的tol_cost别设得比机器精度还小$10^{-12}$ 到 $10^{-14}$ 是合理下限。4.4 初值选择的经验法则LM 对初值比高斯-牛顿宽容但不是无限宽容。强非线性问题比如多峰、周期性强对初值依然敏感。实用的做法是先用粗略的网格搜索或者随机多起点跑几轮挑残差最小的几个作为 LM 的初值如果问题有物理约束把初值限制在合理范围内对于序列数据比如轨迹优化用上一帧的结果做初值通常效果很好这也是增量式优化的常见套路。还有一个技巧如果发现 LM 从某个初值出发总是收敛到局部极小可以先把 $\lambda$ 设大跑几十步相当于先做梯度下降“热身”再让 $\lambda$ 自然降下来进入高斯-牛顿模式。这种“先稳后快”的策略在残差曲面崎岖的问题上很管用。我在实际项目里养成的习惯是每次接入一个新的非线性最小二乘问题先不急着调算法参数而是花时间把残差函数和雅可比在几个随机点上验证一遍再画一画残差随参数变化的切片图对曲面的“地形”有个直观感受。这个前置工作看起来费时但能省下后面大量盲调参数的时间。LM 算法本身足够稳健大多数不收敛的问题根子都在模型表达、雅可比正确性和参数尺度这三件事上算法参数反而是最后才需要动的。
阅读完成 · 觉得有帮助？
