LM算法深度解析：非线性最小二乘拟合的Python实现与工程实践

简介：面向数值计算与数据拟合学习者，提供基于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 对角元也大，λ 乘以这个对角元后不会把这个维度的步长完全压死；而量级很小的参数维度，阻尼对其步长的影响也同步缩小，不会因为统一缩放被忽略。

λ 的取值直接决定算法落在行为谱系的哪一端，我把这个对应关系整理成一张便于调参时对照的表：

工程实现时，我一般把 λ 的初始值取 0.001，同时引入一个缩放因子 ν，通常取 10。每次迭代先按当前 λ 求解正规方程得到候选步长，再用候选参数计算目标函数值：如果目标函数下降了，说明当前方向可信，λ 除以 ν，让算法更激进地转向高斯牛顿；如果目标函数没有下降甚至上升，说明步长跨过了信任区域，λ 乘以 ν，把方向拉回更保守的梯度下降行为。这种以目标函数升降为信号的调节方式，本质上是一个简单的信任域控制，可靠性足够覆盖绝大多数拟合问题。

2.3 LM 的完整迭代流程与收敛判据

把上面的公式串起来，LM 算法的一次完整迭代可以拆成下面几个步骤：

1. 计算残差 r(θ_k)、雅可比 J(θ_k)，进而得到 A = JᵀJ 和梯度向量 g = -Jᵀr。
2. 对当前 λ 求解 (A + λ·diag(A))·Δ = g，得到候选步长 Δ。
3. 用 θ_k + Δ 计算候选残差和目标函数值 S_candidate。
4. 若 S_candidate < S_k，接受步长，λ 缩小，进入下一次外层迭代。
5. 若 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, step=1e-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, lam0=1e-3, nu=10.0, tol_grad=1e-8, tol_step=1e-8, max_iter=200): """ Levenberg-Marquardt 拟合入口。 返回: p : 最优参数向量 history : 每次成功迭代的目标函数值列表 converged : 是否正常收敛 """ p = np.asarray(p0, dtype=float) 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, ord=np.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 关键参数的含义与合理取值范围

参数调整往往是实际拟合中最耗时间的部分，把各个参数单独拆开说明：

有一个常见的误导性认知：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(size=xdata.size) # 自定义实现 p_custom, _, conv = lm_fit(model, [1.0, 1.0, 0.0], xdata, ydata) # SciPy 实现，method='lm' 对应 LM 算法 res = least_squares( lambda p: ydata - model(p, xdata), [1.0, 1.0, 0.0], method='lm' ) 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_fit：

rng = 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(size=xdata.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.02，b 的偏差在 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% 计算预测误差。如果训练集误差远小于验证集误差，说明模型复杂度超过了数据能支撑的信息量。

当噪声方差为 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, grid=3): """ 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)

这里隐含的假设是残差独立同分布且服从正态分布。当数据存在自相关或异方差时，比如噪声幅度随信号大小变化，这个估计会有偏差，正确处理是改为加权最小二乘，把每个样本的权重设为方差的倒数，再用同样的协方差公式。报告拟合结果时，把参数值和标准差一起给出，比只给残差平方和更能体现拟合质量的边界。

本文还有配套的精品资源，点击获取