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.5，c 取 2，d 取 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_iter=200, tol_g=1e-8, tol_step=1e-8, tol_r=1e-10): p = np.array(p0, dtype=float) 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 参数选择与调参记录

拿一组模拟数据跑一下，真值设 a=2.0、b=0.8、c=3.0、d=0.5，时间轴取 0 到 10 共 200 个点，加 2% 的高斯噪声。初值按前面说的估：a=1.5、b=0.5、c=2.0、d=0.0。初始 λ 设 1e-3，跑下来大概 18 次迭代收敛，最终参数误差都在 1% 以内，残差平方和从初值的几十降到个位数。

调参过程中有几个观察值得记下来。初始 λ 设 1e-3 是个比较通用的起点，但如果初值特别离谱，可以设大一点比如 1e-1，让前几步先走保守路线，等增益比稳定了再放开。λ 的缩放因子用 3 和 2 这组不对称的值，是因为减小 λ 时希望快一点放开、增大时希望稳一点收紧，这个不对称是有意为之的。最大迭代次数设 200 足够，正常问题几十步就收敛了，跑到 200 还没停基本说明有问题，该去查雅可比或者尺度了。

还有一个实测细节：如果数据里有个别离群点，残差平方和会被它们主导，LM 会花很多步去迁就这些点，导致整体拟合偏移。这时候可以考虑换鲁棒损失函数（比如 Huber 损失）替代纯平方和，或者先做一轮离群点剔除。LM 本身不解决鲁棒性问题，它只是把给定的目标函数优化好，目标函数选得不对，优化得再准也是白搭。

4.3 数值稳定性处理：矩阵求解那点事

LM 每步都要解一个线性方程组A·Δ = -g，A 是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 常见问题速查表

5.5 几条不写在文档里的经验

第一条，初值比算法更重要。LM 再鲁棒也架不住初值离谱到线性近似完全失效，花时间估一个好初值，比调半天 λ 划算得多。第二条，先在小规模数据上验证，拿几十个点跑通了、雅可比对上了，再上全量数据，否则出了问题根本分不清是算法错还是数据错。第三条，把每步的 λ、ρ、残差都打出来，出问题时这些日志就是破案线索，比事后猜强一百倍。第四条，如果同一个模型要反复拟合很多组数据，把雅可比的解析式固定下来、用自动微分或者符号计算生成代码，比每次数值差分稳得多也快得多。

6. 工程落地时的扩展与取舍

实际项目里很少直接拿裸的 LM 去跑，通常会做一些包装和扩展。一个常见的需求是参数边界约束，比如某个参数物理上必须为正。简单做法是在每步更新后把越界的参数夹回边界内，但这会破坏 LM 的下降性质，可能导致震荡。更稳的做法是投影到可行域后重新计算残差和增益比，或者用变量替换把有界参数映射到无界空间（比如用指数变换保证正性）。变量替换的代价是雅可比要多乘一个变换的导数，但换来的是算法在无约束空间里正常运作，我个人更倾向这种。

另一个扩展是稀疏雅可比的利用。如果残差之间相互独立（比如每个观测点只依赖少数参数），J 会非常稀疏，JᵀJ也是稀疏的。这时候用稀疏矩阵存储和稀疏 Cholesky 分解，能把大规模问题的求解速度提升一个数量级。很多成熟的优化库内部都做了这个优化，手写实现如果遇到参数上千的场景，值得花时间把稀疏性利用起来。

还有一个取舍是要不要用现成的库。如果只是做一次性的拟合，用成熟库里的 LM 实现（很多科学计算框架都有）省时省力，它们把数值稳定性、稀疏性、边界约束都处理好了。但如果你的问题有特殊结构（比如需要自定义阻尼策略、需要嵌入到实时系统里、需要对中间过程做精细监控），手写一个精简版反而更可控。我自己的判断标准是：问题规模小、需求标准，用库；问题有特殊约束、需要深度定制，手写。手写的时候把上面这些坑避开，代码量其实不大，几百行就能覆盖绝大多数场景。

最后分享一个我在实际项目里养成的习惯：每次用 LM 拟合完，不只看最终参数，还会把迭代轨迹画出来——残差随迭代的变化、λ 的变化、每个参数的收敛曲线。这些图能一眼看出收敛是健康的还是勉强压住的，如果 λ 在最后几步还在大幅波动，说明收敛得并不扎实，可能只是碰巧满足了某个停止条件，换个初值就未必了。这个习惯帮我提前发现过好几次潜在的数值问题，比只看最终结果靠谱得多。