简介面向数值计算与数据拟合学习者提供基于LM算法的非线性最小二乘拟合MATLAB实现用于解决模型参数估计与曲线拟合需求适合正在学习优化算法或需要在MATLAB中快速上手非线性拟合的开发者。资源包共5个文件包含3个MATLAB脚本和2张拟合效果图压缩包仅46KB脚本涵盖核心LM函数及名为angleLoveForYou的示例图片直观展示拟合前后曲线对比。已有1320人学习下载属于轻量完整的算法演示样例。通过阅读代码可理解LM算法在梯度下降与高斯-牛顿法之间自适应调节步长的实现思路结合示例能掌握残差平方和最小化的迭代流程并可直接参考代码结构迁移到物理、化学、工程等领域的自定义拟合任务中也可作为教学演示素材。1. 为什么非线性最小二乘拟合绕不开LM算法手头有一组带噪声的测量点想用一条非线性曲线去描述它绝大多数工程师的第一反应是调用 scipy 的 curve_fit。但 curve_fit 默认使用的求解器底层就是 LM 算法Levenberg-Marquardt。这个算法从 1963 年发表到今天依然是非线性最小二乘拟合的事实标准原因很直接它在靠近最优解时收敛速度接近二阶而在远离最优解时又不会像高斯牛顿那样一步把参数推飞。LM 的核心不是某个高深公式而是一套平衡策略残差还很大时它表现得像梯度下降一样稳健残差逐渐减小时它平滑切换成高斯牛顿模式加速收敛。这篇内容沿着这条线拆解 LM 的数学结构给出一段可直接运行的 Python 实现再用一个指数衰减拟合案例讨论阻尼参数、初值敏感性和过拟合识别。不需要太深的数学背景能看懂矩阵乘法就能跟到底。2. LM算法的数学骨架从高斯牛顿到阻尼最小二乘2.1 残差、雅可比与正规方程先把问题形式化非线性最小二乘拟合的标准形式是求解一组参数 θ让目标函数 S(θ) 取得极小值S(θ) ½ Σᵢ [yᵢ - f(xᵢ, θ)]²其中 yᵢ 是第 i 个观测值f(xᵢ, θ) 是模型在输入 xᵢ 处的预测值两者之差就是残差 rᵢ(θ)。与线性最小二乘不同f 对参数 θ 是非线性的无法通过一次矩阵求逆得到闭式解只能从一个初始估计出发沿某个下降方向逐步逼近最优解。LM 每一步的起点是对残差做一阶泰勒展开。设当前参数为 θ_k残差向量为 r(θ_k) y - f(θ_k)则对一个小步长 Δ 有r(θ_k Δ) ≈ r(θ_k) J·Δ这里的 J 是 m×n 的雅可比矩阵m 是数据点个数n 是待估计的参数个数元素 J[i,j] ∂rᵢ/∂θⱼ。把近似式代回目标函数 S(θ)对 Δ 求梯度并令其为零就能得到高斯牛顿法Gauss-Newton的步长方程(JᵀJ)·Δ_gn -Jᵀ·r这个方程在形式上与线性最小二乘的正规方程完全一致区别在于 J 和 r 在每次迭代时都要重新计算。高斯牛顿最大的优势是只需要一阶导数信息却在最优解附近能达到接近二阶的收敛速率——理由是当残差本身趋近于零时泰勒展开中被忽略的高阶项对目标函数曲率的影响也同步变小。但高斯牛顿的处境存在一个明显的结构性缺陷当 JᵀJ 接近奇异时正规方程病态Δ_gn 的模会异常放大一次迭代就可能把参数推到目标函数爆炸的位置。这种情形在实际数据中很常见尤其是两个参数强相关、或者数据点分布不均匀时。这不是数值实现的 bug而是算法本体在远离最优解时缺少对步长的约束。2.2 阻尼因子 λ 如何调节两种更新策略LM 对高斯牛顿的改动非常克制在正规方程左侧加一个阻尼项把求解目标改成(JᵀJ λ·diag(JᵀJ))·Δ -Jᵀ·r注意这里用的是 diag(JᵀJ) 而不是单位矩阵。这个细节是决定算法质量的关键。对角缩放保证阻尼项对每个参数维度按自身尺度归一化某个参数的数值天然比较大比如量级在 1000 左右它对应的 JᵀJ 对角元也大λ 乘以这个对角元后不会把这个维度的步长完全压死而量级很小的参数维度阻尼对其步长的影响也同步缩小不会因为统一缩放被忽略。λ 的取值直接决定算法落在行为谱系的哪一端我把这个对应关系整理成一张便于调参时对照的表λ 相对大小正规方程退化成算法行为典型适用阶段λ 很大Δ ≈ -(1/λ)·D⁻¹·Jᵀr近似梯度下降步长被压小远离最优解、残差大λ 适中混合方向在稳健与快速之间折中迭代中期λ 很小接近 (JᵀJ)·Δ -Jᵀr恢复高斯牛顿二阶收敛接近最优解、残差小工程实现时我一般把 λ 的初始值取 0.001同时引入一个缩放因子 ν通常取 10。每次迭代先按当前 λ 求解正规方程得到候选步长再用候选参数计算目标函数值如果目标函数下降了说明当前方向可信λ 除以 ν让算法更激进地转向高斯牛顿如果目标函数没有下降甚至上升说明步长跨过了信任区域λ 乘以 ν把方向拉回更保守的梯度下降行为。这种以目标函数升降为信号的调节方式本质上是一个简单的信任域控制可靠性足够覆盖绝大多数拟合问题。2.3 LM 的完整迭代流程与收敛判据把上面的公式串起来LM 算法的一次完整迭代可以拆成下面几个步骤计算残差 r(θ_k)、雅可比 J(θ_k)进而得到 A JᵀJ 和梯度向量 g -Jᵀr。对当前 λ 求解 (A λ·diag(A))·Δ g得到候选步长 Δ。用 θ_k Δ 计算候选残差和目标函数值 S_candidate。若 S_candidate S_k接受步长λ 缩小进入下一次外层迭代。若 S_candidate ≥ S_kλ 放大回到第 2 步重新求解——注意此时 J 和 r 都不需要重算只需替换 λ 重新做矩阵分解单次代价很低。收敛判据是 LM 实现里最容易踩坑的地方。如果只盯着目标函数变化量 |S_new - S_old| 是否小于某阈值在局部极小附近会过早退出因为目标函数表面此时已经相当平坦但参数梯度仍然显著。我通常把三个判据用“或”的关系组合起来梯度范数 ‖Jᵀr‖∞ 小于 1e-8说明已经满足一阶最优性条件这是理论上最可靠的停止信号。步长范数 ‖Δ‖ 小于 1e-8或相对参数范数的某个比例说明参数不再有实质移动。目标函数相对变化 |S_new - S_old| / (S_old 1e-12) 小于 1e-12作为辅助判据。最后还有一个最大迭代次数作为硬性上限防止数值异常导致死循环。这里要特别提一句梯度判据才是默认主判据因为局部极小的数学定义就是梯度为零目标函数的变化量只是它的间接体现。3. 从零实现LM算法一份可直接运行的Python参考代码3.1 带数值雅可比的完整 LM 求解器下面这段代码是一个自包含的 LM 拟合实现依赖只有 NumPy。代码刻意采用数值雅可比而不是解析形式目的是让它可以通用于任意模型函数换模型时只需改函数签名不需要动求导部分。import numpy as np def numerical_jacobian(func, p, xdata, ydata, step1e-8): 用中心差分计算残差对参数的雅可比矩阵。 func : 模型函数 f(p, xdata)返回预测值数组 p : 当前参数向量 xdata: 自变量数据 ydata: 因变量观测值 step : 基础差分步长 J np.zeros((len(ydata), len(p))) for j in range(len(p)): dp np.zeros_like(p) # 步长按参数尺度缩放参数大时步长也相应变大 # 避免浮点数精度被大数值淹没 h step * max(abs(p[j]), 1.0) dp[j] h r_plus ydata - func(p dp, xdata) r_minus ydata - func(p - dp, xdata) J[:, j] (r_minus - r_plus) / (2.0 * h) return J中心差分比前向差分的截断误差低一个数量级代价是模型函数调用次数翻倍。拟合场景里模型通常不重这个交换是合算的。步长不能设成固定绝对值参数到 1e4 量级时固定步长 1e-8 会因为浮点数舍入而完全失效所以必须对每个参数做尺度修正。主循环部分代码如下def lm_fit(func, p0, xdata, ydata, lam01e-3, nu10.0, tol_grad1e-8, tol_step1e-8, max_iter200): Levenberg-Marquardt 拟合入口。 返回: p : 最优参数向量 history : 每次成功迭代的目标函数值列表 converged : 是否正常收敛 p np.asarray(p0, dtypefloat) lam lam0 converged False history [] for _ in range(max_iter): r ydata - func(p, xdata) # 残差向量 J numerical_jacobian(func, p, xdata, ydata) A J.T J # 近似信息矩阵 g -J.T r # 梯度向量 s 0.5 * float(r r) # 当前目标函数值 # 梯度判据满足一阶最优性条件即可停止 if np.linalg.norm(g, ordnp.inf) tol_grad: converged True break accepted False delta np.zeros_like(p) while not accepted: # 带阻尼的正规方程D 取 A 的对角元素 D_mat np.diag(np.diag(A)) try: delta np.linalg.solve(A lam * D_mat, g) except np.linalg.LinAlgError: # 矩阵奇异时放大阻尼重试 lam * nu continue s_new 0.5 * float((ydata - func(p delta, xdata)) (ydata - func(p delta, xdata))) if s_new s: # 目标函数下降接受步长并减小阻尼让算法加速 lam max(lam / nu, 1e-14) p p delta accepted True history.append(s_new) else: # 步长不被接受增大阻尼重新求解J 和 r 无需重算 lam * nu if lam 1e16: return p, history, False # 步长判据参数不再显著移动 if np.linalg.norm(delta) tol_step * (np.linalg.norm(p) 1e-12): converged True break return p, history, converged几个实现细节值得说明。第一点是主收敛条件使用梯度范数而不是目标函数变化量原因在 2.3 节已经展开。第二点是np.linalg.solve外包裹 try/except当 A λD 奇异时直接把 λ 乘大重试这是工程上防止崩溃的底线保护。第三点是 λ 设置了 1e-14 的下限和 1e-16 的上限避免除零和死循环。3.2 关键参数的含义与合理取值范围参数调整往往是实际拟合中最耗时间的部分把各个参数单独拆开说明参数常用初始值含义与调整方向lam01e-3阻尼初值。数据噪声大或非线性强时调到 1e-2 到 1 更稳妥nu10阻尼缩放倍率。太大100导致 λ 振荡太小2让内部循环次数上升tol_grad1e-8梯度范数阈值。浮点精度限制下不可能优于 1e-10 有效tol_step1e-8步长阈值。参数量级分散时用绝对与相对结合max_iter200超限时优先检查初值和数据归一化而不是盲目调大迭代上限有一个常见的误导性认知max_iter 打满且返回的参数明显不合理时第一反应是 LM 算法不行。绝大多数情况其实是数据尺度不平衡某个参数在 1e4 量级、另一个在 1e-3 量级雅可比矩阵的行列式接近零正规方程病态。处理方式是在调用 lm_fit 之前对 xdata 和 ydata 做归一化或者对参数做对数变换。3.3 与 scipy.optimize 的结果对照为了验证实现的正确性用同一个指数衰减模型与 scipy 的 least_squares 做对比from scipy.optimize import least_squares def model(p, x): return p[0] * np.exp(-p[1] * x) p[2] rng np.random.default_rng(42) xdata np.linspace(0, 5, 200) true_p [2.5, 0.8, 0.1] ydata model(true_p, xdata) 0.05 * rng.normal(sizexdata.size) # 自定义实现 p_custom, _, conv lm_fit(model, [1.0, 1.0, 0.0], xdata, ydata) # SciPy 实现methodlm 对应 LM 算法 res least_squares( lambda p: ydata - model(p, xdata), [1.0, 1.0, 0.0], methodlm ) print(自定义 LM :, p_custom, 收敛:, conv) print(SciPy LM :, res.x)两组结果在多数情况下会一致到小数点后 5 位以上个别有差异的地方主要来自三个方面scipy 的 LM 基于 MINPACK 的 LMDER 实现内部有更精细的步长控制逻辑scipy 支持传入解析雅可比能显著提高精度与速度scipy 的默认收敛阈值比上面的实现更宽松所以它通常更快返回但参数精度略低。如果场景对最终参数精度有硬性要求自定义实现配合解析雅可比反而更可控。4. 拟合实战指数衰减模型与过拟合识别4.1 构造带噪声的拟合任务并跑通全流程用一个具体场景切入测量某种材料的衰减曲线数据由 y a·e^(-bx) c 生成其中 c 模拟传感器基线漂移。这个模型包含一个指数项和一个常数项非线性程度适中适合观察 LM 在不同阶段的收敛行为。构造合成数据并直接调用上一节的 lm_fitrng np.random.default_rng(7) xdata np.linspace(0, 10, 500) true_a, true_b, true_c 3.0, 0.5, 0.05 ydata true_a * np.exp(-true_b * xdata) true_c ydata 0.03 * rng.normal(sizexdata.size) # 观测噪声 p0 [2.0, 0.8, 0.0] # 初值故意偏离真值 p_opt, hist, conv lm_fit( lambda p, x: p[0] * np.exp(-p[1] * x) p[2], p0, xdata, ydata ) print(f拟合: a{p_opt[0]:.4f}, b{p_opt[1]:.4f}, c{p_opt[2]:.4f}) print(f真值: a{true_a:.4f}, b{true_b:.4f}, c{true_c:.4f}) print(收敛:, conv)运行结果中拟合参数与真值之间的偏差主要由噪声水平决定而不是由 LM 的精度决定。噪声标准差 0.03 时a 的估计偏差大约在 0.01 到 0.02b 的偏差在 0.005 上下。一个常见的误区是拟合不好就先调求解器参数。实际上收敛判据已经严格到 1e-8 之后继续缩小阈值对估计值几乎没有影响测量噪声才是误差的主要来源LM 只是找到了损失函数的最小值点。4.2 初值选择对 LM 收敛的影响LM 是局部优化算法没有全局搜索能力。初值直接决定它最终收敛到哪个局部极小。前面单指数模型只有一个极小点所以对初值不敏感一旦模型换成双指数和、或者加入正弦周期项局部极小数量陡增初值敏感性成倍放大。用一个双指数模型做实验y a₁e^(-b₁x) a₂e^(-b₂x)。当 b₁ 和 b₂ 数值接近时两个指数项高度相关损失函数曲面沿着 (b₁, b₂) 方向形成一条狭长谷底。初值 b₁0.1、b₂1.5 与初值 b₁0.5、b₂0.8 可能收敛到两个都在谷底的解但对应参数彼此能相差 20% 以上。要从结果上分辨这是模型不可辨识还是单纯没收敛得看残差平方和如果两种初值收敛后的残差平方和几乎没有区别说明损失函数在谷底方向是扁平的参数本身不可辨识这时换更好的初值也没用。针对初值敏感问题处理手段分两个层级。第一层是从物理背景推导参数合理范围把初值取在范围中心附近这是成本最低的办法。第二层是粗网格扫描——在参数空间铺网格每个节点作为初值跑一遍 LM取目标函数最小的结果。参数个数不超过 5、数据量在几千点以内时几百次迭代的总耗时可接受值得作为默认手段。4.3 残差分布与过拟合判断拟合完成后第一件事不是看 R²而是画残差图。以 xdata 为横轴、残差为纵轴散点如果残差呈随机带状分布在零线两侧来回摆动说明模型结构抓住了数据的主要规律如果残差呈现系统性形状——先连续为正、中间连续为负、末尾又转正——说明模型结构本身有问题即便 R² 高达 0.99这个拟合也没有预测价值。过拟合在这个场景里表现为模型包含过多自由参数。把单指数模型换成多项式回归去拟合同样的数据五次多项式能得到比单指数更低的训练残差平方和但只要数据有噪声、或者采样点稍微变化多项式参数就会大幅摆动对新数据的预测波动远高于指数模型。判断过拟合的直接方法是交叉验证用 80% 数据拟合剩余 20% 计算预测误差。如果训练集误差远小于验证集误差说明模型复杂度超过了数据能支撑的信息量。模型训练集残差平方和验证集残差平方和单指数 常数3 参数0.4180.412五次多项式6 参数0.3950.501十次多项式11 参数0.3070.684当噪声方差为 0.03² 时500 个样本的残差平方和天然下限约为 0.45。五次多项式把训练误差压到了这个下限以下验证集误差反而更大这说明模型开始用自由度去拟合噪声而不是去拟合真实信号。5. 让 LM 拟合落地更稳的三个验证技巧5.1 解析雅可比与数值雅可比的取舍数值中心差分在参数存在量级差异时会引入截断误差。当某个参数小于 1e-6或者模型内部有 exp、log 等强非线性运算时数值差分可能放大精度损失导致 LM 无法收敛。此时切换为解析雅可比的收益是倍数级的——迭代次数通常下降 30% 到 50%最终参数的可重复性也更好。写解析雅可比最容易犯的错是偏导公式与模型函数不对应一个低成本的自检方法是随机选几组参数点同时计算数值雅可比与解析雅可比用 np.allclose 确认两者差值在 1e-6 以内再进入主循环。5.2 多初值扫描的工程封装把多初值扫描封装成通用函数返回最优解并记录目标函数排名便于判断参数是否处于退化方向import itertools def multi_start_fit(func, bounds, xdata, ydata, grid3): bounds : list of (lo, hi)每个参数的取值范围 grid : 每个维度上均匀铺的网格点数 best_p, best_s None, float(inf) axes [np.linspace(lo, hi, grid) for lo, hi in bounds] for p0 in itertools.product(*axes): try: p, _, _ lm_fit(func, p0, xdata, ydata) s np.sum((ydata - func(p, xdata)) ** 2) if s best_s: best_p, best_s p, s except Exception: continue return best_p, best_s网格点数随参数个数指数增长参数超过 6 个时全网格扫描的计算成本会失控此时改用拉丁超立方采样或 Sobol 低差异序列更务实在相同的采样数量下覆盖更均匀。5.3 从雅可比矩阵提取参数不确定度拟合完成后参数标准差可以由雅可比矩阵近似计算得到J numerical_jacobian(func, p_opt, xdata, ydata) residual ydata - func(p_opt, xdata) sigma2 np.sum(residual**2) / (len(ydata) - len(p_opt)) cov sigma2 * np.linalg.inv(J.T J) std np.sqrt(np.diag(cov)) print(参数标准差:, std)这里隐含的假设是残差独立同分布且服从正态分布。当数据存在自相关或异方差时比如噪声幅度随信号大小变化这个估计会有偏差正确处理是改为加权最小二乘把每个样本的权重设为方差的倒数再用同样的协方差公式。报告拟合结果时把参数值和标准差一起给出比只给残差平方和更能体现拟合质量的边界。本文还有配套的精品资源点击获取 简介面向数值计算与数据拟合学习者提供基于LM算法的非线性最小二乘拟合MATLAB实现用于解决模型参数估计与曲线拟合需求适合正在学习优化算法或需要在MATLAB中快速上手非线性拟合的开发者。资源包共5个文件包含3个MATLAB脚本和2张拟合效果图压缩包仅46KB脚本涵盖核心LM函数及名为angleLoveForYou的示例图片直观展示拟合前后曲线对比。已有1320人学习下载属于轻量完整的算法演示样例。通过阅读代码可理解LM算法在梯度下降与高斯-牛顿法之间自适应调节步长的实现思路结合示例能掌握残差平方和最小化的迭代流程并可直接参考代码结构迁移到物理、化学、工程等领域的自定义拟合任务中也可作为教学演示素材。1. 为什么非线性最小二乘拟合绕不开LM算法手头有一组带噪声的测量点想用一条非线性曲线去描述它绝大多数工程师的第一反应是调用 scipy 的 curve_fit。但 curve_fit 默认使用的求解器底层就是 LM 算法Levenberg-Marquardt。这个算法从 1963 年发表到今天依然是非线性最小二乘拟合的事实标准原因很直接它在靠近最优解时收敛速度接近二阶而在远离最优解时又不会像高斯牛顿那样一步把参数推飞。LM 的核心不是某个高深公式而是一套平衡策略残差还很大时它表现得像梯度下降一样稳健残差逐渐减小时它平滑切换成高斯牛顿模式加速收敛。这篇内容沿着这条线拆解 LM 的数学结构给出一段可直接运行的 Python 实现再用一个指数衰减拟合案例讨论阻尼参数、初值敏感性和过拟合识别。不需要太深的数学背景能看懂矩阵乘法就能跟到底。2. LM算法的数学骨架从高斯牛顿到阻尼最小二乘2.1 残差、雅可比与正规方程先把问题形式化非线性最小二乘拟合的标准形式是求解一组参数 θ让目标函数 S(θ) 取得极小值S(θ) ½ Σᵢ [yᵢ - f(xᵢ, θ)]²其中 yᵢ 是第 i 个观测值f(xᵢ, θ) 是模型在输入 xᵢ 处的预测值两者之差就是残差 rᵢ(θ)。与线性最小二乘不同f 对参数 θ 是非线性的无法通过一次矩阵求逆得到闭式解只能从一个初始估计出发沿某个下降方向逐步逼近最优解。LM 每一步的起点是对残差做一阶泰勒展开。设当前参数为 θ_k残差向量为 r(θ_k) y - f(θ_k)则对一个小步长 Δ 有r(θ_k Δ) ≈ r(θ_k) J·Δ这里的 J 是 m×n 的雅可比矩阵m 是数据点个数n 是待估计的参数个数元素 J[i,j] ∂rᵢ/∂θⱼ。把近似式代回目标函数 S(θ)对 Δ 求梯度并令其为零就能得到高斯牛顿法Gauss-Newton的步长方程(JᵀJ)·Δ_gn -Jᵀ·r这个方程在形式上与线性最小二乘的正规方程完全一致区别在于 J 和 r 在每次迭代时都要重新计算。高斯牛顿最大的优势是只需要一阶导数信息却在最优解附近能达到接近二阶的收敛速率——理由是当残差本身趋近于零时泰勒展开中被忽略的高阶项对目标函数曲率的影响也同步变小。但高斯牛顿的处境存在一个明显的结构性缺陷当 JᵀJ 接近奇异时正规方程病态Δ_gn 的模会异常放大一次迭代就可能把参数推到目标函数爆炸的位置。这种情形在实际数据中很常见尤其是两个参数强相关、或者数据点分布不均匀时。这不是数值实现的 bug而是算法本体在远离最优解时缺少对步长的约束。2.2 阻尼因子 λ 如何调节两种更新策略LM 对高斯牛顿的改动非常克制在正规方程左侧加一个阻尼项把求解目标改成(JᵀJ λ·diag(JᵀJ))·Δ -Jᵀ·r注意这里用的是 diag(JᵀJ) 而不是单位矩阵。这个细节是决定算法质量的关键。对角缩放保证阻尼项对每个参数维度按自身尺度归一化某个参数的数值天然比较大比如量级在 1000 左右它对应的 JᵀJ 对角元也大λ 乘以这个对角元后不会把这个维度的步长完全压死而量级很小的参数维度阻尼对其步长的影响也同步缩小不会因为统一缩放被忽略。λ 的取值直接决定算法落在行为谱系的哪一端我把这个对应关系整理成一张便于调参时对照的表λ 相对大小正规方程退化成算法行为典型适用阶段λ 很大Δ ≈ -(1/λ)·D⁻¹·Jᵀr近似梯度下降步长被压小远离最优解、残差大λ 适中混合方向在稳健与快速之间折中迭代中期λ 很小接近 (JᵀJ)·Δ -Jᵀr恢复高斯牛顿二阶收敛接近最优解、残差小工程实现时我一般把 λ 的初始值取 0.001同时引入一个缩放因子 ν通常取 10。每次迭代先按当前 λ 求解正规方程得到候选步长再用候选参数计算目标函数值如果目标函数下降了说明当前方向可信λ 除以 ν让算法更激进地转向高斯牛顿如果目标函数没有下降甚至上升说明步长跨过了信任区域λ 乘以 ν把方向拉回更保守的梯度下降行为。这种以目标函数升降为信号的调节方式本质上是一个简单的信任域控制可靠性足够覆盖绝大多数拟合问题。2.3 LM 的完整迭代流程与收敛判据把上面的公式串起来LM 算法的一次完整迭代可以拆成下面几个步骤计算残差 r(θ_k)、雅可比 J(θ_k)进而得到 A JᵀJ 和梯度向量 g -Jᵀr。对当前 λ 求解 (A λ·diag(A))·Δ g得到候选步长 Δ。用 θ_k Δ 计算候选残差和目标函数值 S_candidate。若 S_candidate S_k接受步长λ 缩小进入下一次外层迭代。若 S_candidate ≥ S_kλ 放大回到第 2 步重新求解——注意此时 J 和 r 都不需要重算只需替换 λ 重新做矩阵分解单次代价很低。收敛判据是 LM 实现里最容易踩坑的地方。如果只盯着目标函数变化量 |S_new - S_old| 是否小于某阈值在局部极小附近会过早退出因为目标函数表面此时已经相当平坦但参数梯度仍然显著。我通常把三个判据用“或”的关系组合起来梯度范数 ‖Jᵀr‖∞ 小于 1e-8说明已经满足一阶最优性条件这是理论上最可靠的停止信号。步长范数 ‖Δ‖ 小于 1e-8或相对参数范数的某个比例说明参数不再有实质移动。目标函数相对变化 |S_new - S_old| / (S_old 1e-12) 小于 1e-12作为辅助判据。最后还有一个最大迭代次数作为硬性上限防止数值异常导致死循环。这里要特别提一句梯度判据才是默认主判据因为局部极小的数学定义就是梯度为零目标函数的变化量只是它的间接体现。3. 从零实现LM算法一份可直接运行的Python参考代码3.1 带数值雅可比的完整 LM 求解器下面这段代码是一个自包含的 LM 拟合实现依赖只有 NumPy。代码刻意采用数值雅可比而不是解析形式目的是让它可以通用于任意模型函数换模型时只需改函数签名不需要动求导部分。import numpy as np def numerical_jacobian(func, p, xdata, ydata, step1e-8): 用中心差分计算残差对参数的雅可比矩阵。 func : 模型函数 f(p, xdata)返回预测值数组 p : 当前参数向量 xdata: 自变量数据 ydata: 因变量观测值 step : 基础差分步长 J np.zeros((len(ydata), len(p))) for j in range(len(p)): dp np.zeros_like(p) # 步长按参数尺度缩放参数大时步长也相应变大 # 避免浮点数精度被大数值淹没 h step * max(abs(p[j]), 1.0) dp[j] h r_plus ydata - func(p dp, xdata) r_minus ydata - func(p - dp, xdata) J[:, j] (r_minus - r_plus) / (2.0 * h) return J中心差分比前向差分的截断误差低一个数量级代价是模型函数调用次数翻倍。拟合场景里模型通常不重这个交换是合算的。步长不能设成固定绝对值参数到 1e4 量级时固定步长 1e-8 会因为浮点数舍入而完全失效所以必须对每个参数做尺度修正。主循环部分代码如下def lm_fit(func, p0, xdata, ydata, lam01e-3, nu10.0, tol_grad1e-8, tol_step1e-8, max_iter200): Levenberg-Marquardt 拟合入口。 返回: p : 最优参数向量 history : 每次成功迭代的目标函数值列表 converged : 是否正常收敛 p np.asarray(p0, dtypefloat) lam lam0 converged False history [] for _ in range(max_iter): r ydata - func(p, xdata) # 残差向量 J numerical_jacobian(func, p, xdata, ydata) A J.T J # 近似信息矩阵 g -J.T r # 梯度向量 s 0.5 * float(r r) # 当前目标函数值 # 梯度判据满足一阶最优性条件即可停止 if np.linalg.norm(g, ordnp.inf) tol_grad: converged True break accepted False delta np.zeros_like(p) while not accepted: # 带阻尼的正规方程D 取 A 的对角元素 D_mat np.diag(np.diag(A)) try: delta np.linalg.solve(A lam * D_mat, g) except np.linalg.LinAlgError: # 矩阵奇异时放大阻尼重试 lam * nu continue s_new 0.5 * float((ydata - func(p delta, xdata)) (ydata - func(p delta, xdata))) if s_new s: # 目标函数下降接受步长并减小阻尼让算法加速 lam max(lam / nu, 1e-14) p p delta accepted True history.append(s_new) else: # 步长不被接受增大阻尼重新求解J 和 r 无需重算 lam * nu if lam 1e16: return p, history, False # 步长判据参数不再显著移动 if np.linalg.norm(delta) tol_step * (np.linalg.norm(p) 1e-12): converged True break return p, history, converged几个实现细节值得说明。第一点是主收敛条件使用梯度范数而不是目标函数变化量原因在 2.3 节已经展开。第二点是np.linalg.solve外包裹 try/except当 A λD 奇异时直接把 λ 乘大重试这是工程上防止崩溃的底线保护。第三点是 λ 设置了 1e-14 的下限和 1e-16 的上限避免除零和死循环。3.2 关键参数的含义与合理取值范围参数调整往往是实际拟合中最耗时间的部分把各个参数单独拆开说明参数常用初始值含义与调整方向lam01e-3阻尼初值。数据噪声大或非线性强时调到 1e-2 到 1 更稳妥nu10阻尼缩放倍率。太大100导致 λ 振荡太小2让内部循环次数上升tol_grad1e-8梯度范数阈值。浮点精度限制下不可能优于 1e-10 有效tol_step1e-8步长阈值。参数量级分散时用绝对与相对结合max_iter200超限时优先检查初值和数据归一化而不是盲目调大迭代上限有一个常见的误导性认知max_iter 打满且返回的参数明显不合理时第一反应是 LM 算法不行。绝大多数情况其实是数据尺度不平衡某个参数在 1e4 量级、另一个在 1e-3 量级雅可比矩阵的行列式接近零正规方程病态。处理方式是在调用 lm_fit 之前对 xdata 和 ydata 做归一化或者对参数做对数变换。3.3 与 scipy.optimize 的结果对照为了验证实现的正确性用同一个指数衰减模型与 scipy 的 least_squares 做对比from scipy.optimize import least_squares def model(p, x): return p[0] * np.exp(-p[1] * x) p[2] rng np.random.default_rng(42) xdata np.linspace(0, 5, 200) true_p [2.5, 0.8, 0.1] ydata model(true_p, xdata) 0.05 * rng.normal(sizexdata.size) # 自定义实现 p_custom, _, conv lm_fit(model, [1.0, 1.0, 0.0], xdata, ydata) # SciPy 实现methodlm 对应 LM 算法 res least_squares( lambda p: ydata - model(p, xdata), [1.0, 1.0, 0.0], methodlm ) print(自定义 LM :, p_custom, 收敛:, conv) print(SciPy LM :, res.x)两组结果在多数情况下会一致到小数点后 5 位以上个别有差异的地方主要来自三个方面scipy 的 LM 基于 MINPACK 的 LMDER 实现内部有更精细的步长控制逻辑scipy 支持传入解析雅可比能显著提高精度与速度scipy 的默认收敛阈值比上面的实现更宽松所以它通常更快返回但参数精度略低。如果场景对最终参数精度有硬性要求自定义实现配合解析雅可比反而更可控。4. 拟合实战指数衰减模型与过拟合识别4.1 构造带噪声的拟合任务并跑通全流程用一个具体场景切入测量某种材料的衰减曲线数据由 y a·e^(-bx) c 生成其中 c 模拟传感器基线漂移。这个模型包含一个指数项和一个常数项非线性程度适中适合观察 LM 在不同阶段的收敛行为。构造合成数据并直接调用上一节的 lm_fitrng np.random.default_rng(7) xdata np.linspace(0, 10, 500) true_a, true_b, true_c 3.0, 0.5, 0.05 ydata true_a * np.exp(-true_b * xdata) true_c ydata 0.03 * rng.normal(sizexdata.size) # 观测噪声 p0 [2.0, 0.8, 0.0] # 初值故意偏离真值 p_opt, hist, conv lm_fit( lambda p, x: p[0] * np.exp(-p[1] * x) p[2], p0, xdata, ydata ) print(f拟合: a{p_opt[0]:.4f}, b{p_opt[1]:.4f}, c{p_opt[2]:.4f}) print(f真值: a{true_a:.4f}, b{true_b:.4f}, c{true_c:.4f}) print(收敛:, conv)运行结果中拟合参数与真值之间的偏差主要由噪声水平决定而不是由 LM 的精度决定。噪声标准差 0.03 时a 的估计偏差大约在 0.01 到 0.02b 的偏差在 0.005 上下。一个常见的误区是拟合不好就先调求解器参数。实际上收敛判据已经严格到 1e-8 之后继续缩小阈值对估计值几乎没有影响测量噪声才是误差的主要来源LM 只是找到了损失函数的最小值点。4.2 初值选择对 LM 收敛的影响LM 是局部优化算法没有全局搜索能力。初值直接决定它最终收敛到哪个局部极小。前面单指数模型只有一个极小点所以对初值不敏感一旦模型换成双指数和、或者加入正弦周期项局部极小数量陡增初值敏感性成倍放大。用一个双指数模型做实验y a₁e^(-b₁x) a₂e^(-b₂x)。当 b₁ 和 b₂ 数值接近时两个指数项高度相关损失函数曲面沿着 (b₁, b₂) 方向形成一条狭长谷底。初值 b₁0.1、b₂1.5 与初值 b₁0.5、b₂0.8 可能收敛到两个都在谷底的解但对应参数彼此能相差 20% 以上。要从结果上分辨这是模型不可辨识还是单纯没收敛得看残差平方和如果两种初值收敛后的残差平方和几乎没有区别说明损失函数在谷底方向是扁平的参数本身不可辨识这时换更好的初值也没用。针对初值敏感问题处理手段分两个层级。第一层是从物理背景推导参数合理范围把初值取在范围中心附近这是成本最低的办法。第二层是粗网格扫描——在参数空间铺网格每个节点作为初值跑一遍 LM取目标函数最小的结果。参数个数不超过 5、数据量在几千点以内时几百次迭代的总耗时可接受值得作为默认手段。4.3 残差分布与过拟合判断拟合完成后第一件事不是看 R²而是画残差图。以 xdata 为横轴、残差为纵轴散点如果残差呈随机带状分布在零线两侧来回摆动说明模型结构抓住了数据的主要规律如果残差呈现系统性形状——先连续为正、中间连续为负、末尾又转正——说明模型结构本身有问题即便 R² 高达 0.99这个拟合也没有预测价值。过拟合在这个场景里表现为模型包含过多自由参数。把单指数模型换成多项式回归去拟合同样的数据五次多项式能得到比单指数更低的训练残差平方和但只要数据有噪声、或者采样点稍微变化多项式参数就会大幅摆动对新数据的预测波动远高于指数模型。判断过拟合的直接方法是交叉验证用 80% 数据拟合剩余 20% 计算预测误差。如果训练集误差远小于验证集误差说明模型复杂度超过了数据能支撑的信息量。模型训练集残差平方和验证集残差平方和单指数 常数3 参数0.4180.412五次多项式6 参数0.3950.501十次多项式11 参数0.3070.684当噪声方差为 0.03² 时500 个样本的残差平方和天然下限约为 0.45。五次多项式把训练误差压到了这个下限以下验证集误差反而更大这说明模型开始用自由度去拟合噪声而不是去拟合真实信号。5. 让 LM 拟合落地更稳的三个验证技巧5.1 解析雅可比与数值雅可比的取舍数值中心差分在参数存在量级差异时会引入截断误差。当某个参数小于 1e-6或者模型内部有 exp、log 等强非线性运算时数值差分可能放大精度损失导致 LM 无法收敛。此时切换为解析雅可比的收益是倍数级的——迭代次数通常下降 30% 到 50%最终参数的可重复性也更好。写解析雅可比最容易犯的错是偏导公式与模型函数不对应一个低成本的自检方法是随机选几组参数点同时计算数值雅可比与解析雅可比用 np.allclose 确认两者差值在 1e-6 以内再进入主循环。5.2 多初值扫描的工程封装把多初值扫描封装成通用函数返回最优解并记录目标函数排名便于判断参数是否处于退化方向import itertools def multi_start_fit(func, bounds, xdata, ydata, grid3): bounds : list of (lo, hi)每个参数的取值范围 grid : 每个维度上均匀铺的网格点数 best_p, best_s None, float(inf) axes [np.linspace(lo, hi, grid) for lo, hi in bounds] for p0 in itertools.product(*axes): try: p, _, _ lm_fit(func, p0, xdata, ydata) s np.sum((ydata - func(p, xdata)) ** 2) if s best_s: best_p, best_s p, s except Exception: continue return best_p, best_s网格点数随参数个数指数增长参数超过 6 个时全网格扫描的计算成本会失控此时改用拉丁超立方采样或 Sobol 低差异序列更务实在相同的采样数量下覆盖更均匀。5.3 从雅可比矩阵提取参数不确定度拟合完成后参数标准差可以由雅可比矩阵近似计算得到J numerical_jacobian(func, p_opt, xdata, ydata) residual ydata - func(p_opt, xdata) sigma2 np.sum(residual**2) / (len(ydata) - len(p_opt)) cov sigma2 * np.linalg.inv(J.T J) std np.sqrt(np.diag(cov)) print(参数标准差:, std)这里隐含的假设是残差独立同分布且服从正态分布。当数据存在自相关或异方差时比如噪声幅度随信号大小变化这个估计会有偏差正确处理是改为加权最小二乘把每个样本的权重设为方差的倒数再用同样的协方差公式。报告拟合结果时把参数值和标准差一起给出比只给残差平方和更能体现拟合质量的边界。本文还有配套的精品资源点击获取