1. 从一次拟合翻车说起LM算法到底在解决什么问题几年前接手一个传感器标定的小项目手里有一组实验数据理论上完全符合一个非线性模型参数就那么五六个我心想随便找个优化器跑一跑就完事了。结果梯度下降法在接近最优解的地方开始疯狂震荡步长稍微调大一点就发散调小一点又慢得像蜗牛换成高斯-牛顿法迭代几步直接报矩阵奇异残差平方和反而越跑越大。折腾了大半天最后把Levenberg-Marquardt 算法下文简称 LM 算法搬出来设了个初始阻尼因子十几步迭代就稳稳收敛了。那次经历让我彻底记住了这个算法——它不是什么新潮的东西诞生至今已经几十年但在非线性最小二乘这个领域里它依然是工业界和学术界最常被拿出来用的那一档工具。LM 算法本质上是一种求解非线性最小二乘问题的迭代优化算法。说人话就是你有一堆观测数据也有一个带未知参数的数学模型模型算出来的值和实际观测值之间有偏差你想找到一组参数让这些偏差的平方和最小。这个问题在工程里太常见了——曲线拟合、相机标定、机器人位姿估计、神经网络训练里的某些环节、甚至经济学里的参数估计背后都可能是同一个数学结构。LM 算法能做的事情就是在高斯-牛顿法和梯度下降法之间做一个自适应切换离最优解远的时候它表现得像梯度下降步子稳、不容易跑飞离最优解近的时候它又切换成高斯-牛顿的形态收敛快、精度高。这个“切换”靠的是一个叫阻尼因子的东西来调节后面会详细拆。这篇文章适合谁看如果你写过最小二乘拟合但被发散问题折磨过如果你在做视觉、机器人、信号处理相关的工作需要手写优化器或者你只是好奇为什么那么多库的底层都藏着一个叫 LM 的东西那这篇内容应该能给你一些直接能抄的作业。我会从设计思路讲起把核心公式拆开给出可以直接跑的代码再把我踩过的坑和排查经验整理成表。全程尽量说人话复杂的数学用生活化的类比带过去但关键的参数和步骤不会含糊。2. 算法整体设计与思路拆解2.1 为什么偏偏要在两种方法之间“骑墙”要理解 LM 的设计得先搞清楚它调和的两个极端。梯度下降法的思路很朴素沿着当前残差下降最快的方向走一小步再重新算方向。它的优点是几乎不会发散只要步长够小残差平方和一定在降缺点是收敛慢尤其在接近最优解的区域梯度趋近于零方向信息变得不可靠走起来像在泥地里挪步。高斯-牛顿法则激进得多它用一阶泰勒展开去近似非线性模型把问题转化成一系列线性最小二乘问题来解每一步都试图直接跳到线性化后的最优点。它在最优解附近收敛极快通常是二次收敛但代价是如果当前点离最优解很远线性近似根本不成立算出来的步长可能巨大无比直接把参数甩到天边去甚至导致残差增大、矩阵奇异。LM 算法的聪明之处在于它不选边站而是引入一个阻尼因子 λ有些文献写作 μ把高斯-牛顿的法方程改造成一个带阻尼的形式。当 λ 很大时阻尼项主导更新方向退化成梯度下降方向的一个缩放版步子小且稳当 λ 很小时阻尼项可以忽略算法就退化成纯高斯-牛顿。迭代过程中λ 会根据每一步的实际效果动态调整如果这一步让残差下降了说明当前线性近似还算靠谱就减小 λ让算法更激进一点如果这一步让残差上升了说明冲太猛了就增大 λ退回来稳一稳。这个机制用一句话概括就是用残差的实际变化来投票决定下一步该保守还是该激进。这种设计的优势非常实际。它避免了高斯-牛顿在初值不好时直接崩盘的问题同时又不像梯度下降那样在最优解附近磨蹭。代价是每一步要多维护一个阻尼因子、多判断一次残差变化计算量略增但换来的是鲁棒性的大幅提升。在工程实践中初值往往靠经验或者粗估计给出不可能保证落在最优解附近这时候 LM 的“骑墙”特性就是救命稻草。2.2 阻尼项背后的几何直觉很多人第一次看 LM 的公式会觉得那个 λ 加得莫名其妙为什么是在对角线上加这里给一个几何上的解释。高斯-牛顿的法方程是(JᵀJ) Δ -Jᵀr其中 J 是残差对参数的雅可比矩阵r 是残差向量。矩阵JᵀJ近似了目标函数的曲率信息海森矩阵的一部分。当某些参数方向上的曲率很小甚至接近零时JᵀJ接近奇异解出来的 Δ 会非常大这就是发散和震荡的根源。在对角线上加上λ严格说是λ·diag(JᵀJ)或λ·I两种变体后面会讲相当于给每个参数方向的曲率“垫”了一个下限。你可以把它想象成给一个凹凸不平的碗底铺了一层软垫原来某些方向几乎是平的解会在上面滑得很远垫了之后每个方向都有了最低限度的“坡度”解就不会跑得没边。λ 越大垫得越厚解越靠近原点方向即步子越小λ 越小垫得越薄越接近原始的高斯-牛顿解。这个几何图像比纯公式好记多了调参的时候心里也有底残差震荡就加厚垫子收敛太慢就减薄垫子。2.3 两种阻尼变体的取舍实际实现里阻尼项的写法有两种主流变体选哪种会直接影响数值表现。第一种是(JᵀJ λI) Δ -Jᵀr即在所有对角线上加同一个 λ这叫加性阻尼或直接叫 Levenberg 形式。第二种是(JᵀJ λ·diag(JᵀJ)) Δ -Jᵀr即按对角线元素的比例来加这叫乘性阻尼是 Marquardt 后来提出的改进。两者的区别在于对参数尺度的敏感性。如果不同参数的量级差异很大比如一个参数是几千另一个是零点几加性阻尼对大量级参数几乎没影响对小量级参数却可能过度压制导致收敛不均衡。乘性阻尼按各自的对角元素缩放相当于做了某种尺度归一化对量级差异大的问题更友好。我个人的习惯是默认用乘性阻尼只有在参数量级比较统一、且需要严格的理论分析时才用加性形式。这个选择没有绝对的对错但知道背后的理由遇到收敛异常时就知道该往哪个方向排查。3. 核心细节解析与实操要点3.1 雅可比矩阵整个算法的发动机LM 算法每一步的核心计算量都花在雅可比矩阵 J上。J 的第 i 行第 j 列元素是第 i 个残差对第 j 个参数的偏导数。这个矩阵刻画了“参数动一点点残差会怎么变”是整个算法获取方向信息的唯一来源。J 算得准不准、快不快直接决定了 LM 的表现。获取 J 有两条路。第一条是解析求导手动推导每个残差对每个参数的偏导表达式写成代码。优点是精度高、速度快缺点是推导容易出错模型一改就得重新推。第二条是数值差分用(r(x h·eⱼ) - r(x)) / h来近似第 j 列的偏导其中 eⱼ 是第 j 个单位向量h 是一个小步长。优点是通用、不用推导缺点是每个参数都要多算一次残差参数多的时候开销成倍增长而且 h 的选取很讲究——太大精度不够太小会被浮点误差淹没。我的经验是参数少于十个、模型不太复杂时优先解析求导一次推导省下后面无数次差分的开销参数多或者模型频繁变动时用数值差分图省事但 h 要选在sqrt(机器精度) * max(|xⱼ|, 1)这个量级附近大约是 1e-8 左右实测比较稳。如果模型本身支持自动微分比如用某些科学计算框架写的那是最省心的精度和速度都兼顾直接拿来用就行。注意数值差分时如果某个参数当前值恰好是零max(|xⱼ|, 1)里的 1 就是防止 h 变成零的保护这个细节很多手写实现会漏掉导致该列雅可比全为零算法在那个方向上彻底失去信息。3.2 阻尼因子的更新策略增益比是裁判λ 怎么调是 LM 实现里最见功力的地方。主流做法是引入一个增益比 ρgain ratio定义为实际残差下降量与线性模型预测的下降量之比。具体来说当前残差平方和是F(x)按 LM 步长算出候选点x_new后实际下降是F(x) - F(x_new)而线性化模型预测的下降量可以用Δᵀ(λ·D·Δ - Jᵀr)这类表达式算出来不同文献写法略有差异核心是用当前线性近似估计“理论上该降多少”。两者一比得到 ρ。ρ 的含义很直观如果 ρ 接近 1说明线性近似很准实际降的和预测的差不多可以放心减小 λ、加大步子如果 ρ 是负数或者很小说明线性近似完全失灵实际残差反而涨了必须增大 λ、退回保守模式。常见的更新规则是ρ 大于某个阈值比如 0.75时 λ 除以一个因子比如 3ρ 小于某个阈值比如 0.25时 λ 乘以一个因子比如 2中间区域保持不变。这些数字不是金科玉律但经过大量实践检验作为起点很靠谱。还有一个细节容易被忽略λ 要有上下界。下界防止它小到零导致矩阵奇异上界防止它大到更新步长趋近于零、算法原地踏步。我一般设下界 1e-12、上界 1e12配合增益比策略绝大多数问题都能覆盖。3.3 收敛判据什么时候该停迭代不能无限跑得有停下来的条件。LM 常用的收敛判据有三类通常组合使用。第一类是梯度判据如果Jᵀr的无穷范数小于某个容差说明已经接近驻点可以停。第二类是步长判据如果更新量 Δ 的范数相对于参数范数小到可以忽略说明再走也走不动了。第三类是残差判据如果连续若干步残差平方和的相对下降量小于阈值说明已经榨不出更多改进了。三类判据各有盲区。梯度判据在参数尺度差异大时可能误判步长判据在 λ 被顶到上界时可能假性触发残差判据在目标函数有平台区时可能过早停止。我的做法是三个都设任意一个满足就停同时加一个最大迭代次数兜底一般设 100 到 200 次。容差方面梯度容差 1e-8、步长容差 1e-8、残差相对容差 1e-10这套组合在多数中小规模问题上表现稳定。如果你的问题残差量级特别大或特别小这些绝对值要相应缩放别死记。3.4 参数尺度归一化被低估的前置工作这一条是我踩坑最多的地方。LM 虽然比纯高斯-牛顿鲁棒但它对参数尺度依然敏感尤其是用加性阻尼的时候。假设你在拟合一个模型其中一个参数是频率量级几千另一个是相位量级零点几如果不做归一化直接扔进去阻尼因子对两个参数的实际影响天差地别收敛会非常别扭——频率方向早就收敛了相位方向还在慢慢爬。解决办法是在优化前对参数做尺度变换让所有参数大致落在同一个量级。最简单的做法是给每个参数估计一个典型尺度 sⱼ令xⱼ sⱼ · yⱼ对 y 做优化雅可比相应除以 sⱼ。这样阻尼因子对所有方向的作用就均衡了。尺度怎么估用参数的物理量级、初值的绝对值、或者先跑几步粗优化看参数变化范围都行。这一步花不了多少时间但能省下后面大量调 λ 的功夫。我现在的习惯是只要参数超过三个且量级不统一先归一化再说几乎成了肌肉记忆。4. 实操过程与核心环节实现4.1 从一个曲线拟合例子走完整流程光说理论没意思拿一个具体的非线性拟合问题把整个流程走一遍。假设有一组观测数据理论上服从y a · exp(-b · t) · cos(c · t d)这样的衰减振荡模型四个参数 a、b、c、d 待定。数据是带噪声的初值靠肉眼估a 取数据最大值b 取 0.5c 取 2d 取 0。这个初值离真值不算近正好用来检验 LM 的鲁棒性。残差向量定义为每个数据点的rᵢ y_model(tᵢ) - y_obs(tᵢ)目标是最小化Σ rᵢ²。雅可比矩阵是 4 列第 i 行分别是∂rᵢ/∂a exp(-b·tᵢ)·cos(c·tᵢd)、∂rᵢ/∂b -a·tᵢ·exp(-b·tᵢ)·cos(c·tᵢd)、∂rᵢ/∂c -a·tᵢ·exp(-b·tᵢ)·sin(c·tᵢd)、∂rᵢ/∂d -a·exp(-b·tᵢ)·sin(c·tᵢd)。这四个偏导手推一遍不算难但 c 和 d 在三角函数里符号容易搞错推完最好用数值差分交叉验证一下对不上就说明推错了。下面是核心迭代的 Python 实现用乘性阻尼增益比控制 λ参数做了简单归一化。代码不依赖任何优化库纯手写方便你直接改吧改吧用到自己项目里。import numpy as np def lm_fit(t, y_obs, p0, max_iter200, tol_g1e-8, tol_step1e-8, tol_r1e-10): p np.array(p0, dtypefloat) n len(p) lam 1e-3 lam_min, lam_max 1e-12, 1e12 # 参数尺度归一化用初值绝对值的量级作为尺度 scale np.maximum(np.abs(p), 1.0) p p / scale def residuals(p_scaled): a, b, c, d p_scaled * scale return a * np.exp(-b * t) * np.cos(c * t d) - y_obs def jacobian(p_scaled): a, b, c, d p_scaled * scale e np.exp(-b * t) cosv np.cos(c * t d) sinv np.sin(c * t d) J np.zeros((len(t), n)) J[:, 0] e * cosv * scale[0] J[:, 1] -a * t * e * cosv * scale[1] J[:, 2] -a * t * e * sinv * scale[2] J[:, 3] -a * e * sinv * scale[3] return J r residuals(p) F 0.5 * np.dot(r, r) for it in range(max_iter): J jacobian(p) g J.T r if np.linalg.norm(g, np.inf) tol_g: break # 乘性阻尼 JtJ J.T J D np.diag(np.diag(JtJ)) A JtJ lam * D try: delta np.linalg.solve(A, -g) except np.linalg.LinAlgError: lam min(lam * 10, lam_max) continue p_new p delta r_new residuals(p_new) F_new 0.5 * np.dot(r_new, r_new) # 增益比 predicted 0.5 * delta (lam * D delta - g) rho (F - F_new) / (predicted 1e-30) if rho 0.75: lam max(lam / 3.0, lam_min) elif rho 0.25: lam min(lam * 2.0, lam_max) if rho 0: p p_new r r_new F F_new if np.linalg.norm(delta) tol_step * (np.linalg.norm(p) tol_step): break if it 0 and abs(F - F_new) tol_r * F: break return p * scale, F这段代码有几个地方值得单独说。第一scale用初值绝对值的量级简单粗暴但有效避免了参数尺度悬殊的问题。第二雅可比里每个偏导都乘了对应的scale这是因为我们对归一化后的参数求导链式法则要求补上这个因子漏了就会导致方向错误。第三增益比计算里predicted的表达式用了0.5 * delta (lam * D delta - g)这是从线性模型预测下降量推出来的分母加了个 1e-30 防止除零。第四rho 0才接受新点否则保持原参数、只增大 λ 重试这个判断保证了残差单调不增至少不会因为一次坏步长而变差。4.2 参数选择与调参记录拿一组模拟数据跑一下真值设 a2.0、b0.8、c3.0、d0.5时间轴取 0 到 10 共 200 个点加 2% 的高斯噪声。初值按前面说的估a1.5、b0.5、c2.0、d0.0。初始 λ 设 1e-3跑下来大概 18 次迭代收敛最终参数误差都在 1% 以内残差平方和从初值的几十降到个位数。调参过程中有几个观察值得记下来。初始 λ 设 1e-3 是个比较通用的起点但如果初值特别离谱可以设大一点比如 1e-1让前几步先走保守路线等增益比稳定了再放开。λ 的缩放因子用 3 和 2 这组不对称的值是因为减小 λ 时希望快一点放开、增大时希望稳一点收紧这个不对称是有意为之的。最大迭代次数设 200 足够正常问题几十步就收敛了跑到 200 还没停基本说明有问题该去查雅可比或者尺度了。还有一个实测细节如果数据里有个别离群点残差平方和会被它们主导LM 会花很多步去迁就这些点导致整体拟合偏移。这时候可以考虑换鲁棒损失函数比如 Huber 损失替代纯平方和或者先做一轮离群点剔除。LM 本身不解决鲁棒性问题它只是把给定的目标函数优化好目标函数选得不对优化得再准也是白搭。4.3 数值稳定性处理矩阵求解那点事LM 每步都要解一个线性方程组A·Δ -gA 是JᵀJ λD。这个矩阵理论上在 λ 大于零时是正定的但数值上未必。如果 J 的某些列几乎线性相关参数之间强耦合JᵀJ的条件数会非常大加上阻尼后可能依然病态直接求逆或者用普通的高斯消元会放大误差。稳妥的做法是用Cholesky 分解来解因为 A 是对称正定的Cholesky 比一般 LU 快一倍且数值更稳。如果 Cholesky 失败说明 A 不够正定就增大 λ 重试这本身就是一种保护机制。我上面的代码用了np.linalg.solve底层是 LU对中小规模问题够用如果参数上百、对稳定性要求高建议换成scipy.linalg.cho_factor和cho_solve。另外JᵀJ的累加最好用高精度或者分块累加避免大矩阵相乘时的舍入误差累积这个在参数多的时候影响明显。提示如果迭代中频繁出现 Cholesky 失败或者 λ 被顶到上界八成是雅可比有列全零或者参数之间存在完全共线。先去检查雅可比别急着调 λ那是治标不治本。5. 常见问题与排查技巧实录5.1 残差震荡不收敛先查尺度和雅可比残差在几步之内反复上下跳、λ 忽大忽小、迭代次数用完了还没停这是最常见的症状。排查顺序我总结成三步。第一步查参数尺度如果不同参数量级差了两个数量级以上先做归一化十有八九能缓解。第二步查雅可比正确性用数值差分和解析求导对一遍看每一列的相对误差超过 1e-4 就说明有问题常见的是符号错了或者漏了链式法则的因子。第三步查阻尼变体如果用的是加性阻尼且参数尺度不统一换成乘性阻尼试试。还有一种隐蔽的情况残差函数本身在某个参数区域不连续或者有尖点比如用了绝对值、分段函数导致雅可比在那里不存在或者跳变LM 的线性近似彻底失效。这种问题不是调参能解决的要么把模型改光滑要么换用不依赖梯度的优化方法。5.2 收敛太慢λ 下界和增益比阈值在作怪迭代能收敛但步数特别多每次残差只降一点点。先看 λ 是不是被卡在下界附近下不来如果是检查增益比的计算是不是有问题——预测下降量算错会导致 ρ 一直偏小λ 就只增不减。再一个可能是增益比的阈值设得太苛刻比如把减小 λ 的阈值设成 0.9那只有线性近似极准的时候才放开自然慢。把阈值放宽到 0.75 左右通常能提速。另外如果问题本身条件数很差收敛慢是固有的这时候可以考虑换更好的初值或者对参数做更精细的归一化从源头改善条件数。5.3 矩阵奇异报错共线性和零列的识别LinAlgError或者解出来的 Δ 大得离谱基本是矩阵奇异或接近奇异。先算JᵀJ的条件数超过 1e12 就要警惕。然后看 J 的列有没有某一列全零某个参数对残差完全没影响说明模型里这个参数是冗余的有没有两列高度相关两个参数在起同样的作用模型过参数化。这两种情况都要回到模型层面解决——删冗余参数、合并共线参数、或者加正则化项。单纯靠增大 λ 只能压住症状压不住病根而且 λ 太大收敛会慢得没法用。5.4 常见问题速查表症状可能原因排查动作解决方向残差震荡、λ 反复跳参数尺度悬殊检查各参数量级做尺度归一化换乘性阻尼迭代步数多、收敛慢λ 下界过高或增益比阈值苛刻打印 λ 和 ρ 序列放宽阈值检查预测下降量计算矩阵奇异报错雅可比列共线或全零算条件数看列范数删冗余参数加正则化残差降不下去、停在平台模型与数据不匹配或离群点主导看残差分布换模型或换鲁棒损失参数跑到物理上不合理的值缺少约束检查迭代轨迹加参数边界或惩罚项数值差分雅可比不准步长 h 选取不当对比解析求导h 取 1e-8 量级或改自动微分5.5 几条不写在文档里的经验第一条初值比算法更重要。LM 再鲁棒也架不住初值离谱到线性近似完全失效花时间估一个好初值比调半天 λ 划算得多。第二条先在小规模数据上验证拿几十个点跑通了、雅可比对上了再上全量数据否则出了问题根本分不清是算法错还是数据错。第三条把每步的 λ、ρ、残差都打出来出问题时这些日志就是破案线索比事后猜强一百倍。第四条如果同一个模型要反复拟合很多组数据把雅可比的解析式固定下来、用自动微分或者符号计算生成代码比每次数值差分稳得多也快得多。6. 工程落地时的扩展与取舍实际项目里很少直接拿裸的 LM 去跑通常会做一些包装和扩展。一个常见的需求是参数边界约束比如某个参数物理上必须为正。简单做法是在每步更新后把越界的参数夹回边界内但这会破坏 LM 的下降性质可能导致震荡。更稳的做法是投影到可行域后重新计算残差和增益比或者用变量替换把有界参数映射到无界空间比如用指数变换保证正性。变量替换的代价是雅可比要多乘一个变换的导数但换来的是算法在无约束空间里正常运作我个人更倾向这种。另一个扩展是稀疏雅可比的利用。如果残差之间相互独立比如每个观测点只依赖少数参数J 会非常稀疏JᵀJ也是稀疏的。这时候用稀疏矩阵存储和稀疏 Cholesky 分解能把大规模问题的求解速度提升一个数量级。很多成熟的优化库内部都做了这个优化手写实现如果遇到参数上千的场景值得花时间把稀疏性利用起来。还有一个取舍是要不要用现成的库。如果只是做一次性的拟合用成熟库里的 LM 实现很多科学计算框架都有省时省力它们把数值稳定性、稀疏性、边界约束都处理好了。但如果你的问题有特殊结构比如需要自定义阻尼策略、需要嵌入到实时系统里、需要对中间过程做精细监控手写一个精简版反而更可控。我自己的判断标准是问题规模小、需求标准用库问题有特殊约束、需要深度定制手写。手写的时候把上面这些坑避开代码量其实不大几百行就能覆盖绝大多数场景。最后分享一个我在实际项目里养成的习惯每次用 LM 拟合完不只看最终参数还会把迭代轨迹画出来——残差随迭代的变化、λ 的变化、每个参数的收敛曲线。这些图能一眼看出收敛是健康的还是勉强压住的如果 λ 在最后几步还在大幅波动说明收敛得并不扎实可能只是碰巧满足了某个停止条件换个初值就未必了。这个习惯帮我提前发现过好几次潜在的数值问题比只看最终结果靠谱得多。