简介：基于列文伯格-马夸尔特（Levenberg-Marquardt，简称LM）算法的Matlab实现资源包，专为需要利用非线性最小二乘方法完成复杂模型参数估计的科研人员、算法工程师及高年级学生设计。LM算法兼具梯度下降法的全局探索能力与牛顿法的局部快速收敛特性，通过自适应调整阻尼因子，在迭代过程中动态平衡收敛速度与稳定性，尤其适合处理Hessian矩阵近似病态或初值偏离最优解较远的拟合问题。压缩包仅209KB却包含5个文件，2个.m源文件分别提供LMFnlsq2核心函数和配套测试脚本，可直接运行验证；txt文件补充算法基础说明，PDF文档对误差定义、Hessian矩阵近似、步长更新策略和病态情况处理等关键实现进行了详细解释，jpg示意图则帮助直观理解算法效果。读者可以跟随PDF文档中的测试案例完整走通LM算法流程，通过自定义目标函数与观测数据将代码迁移至物理模型拟合、信号处理、滤波器系数估计、神经网络权重优化等实际任务，有效缩短算法开发周期。该资源自发布以来已有2171人学习/下载，以极其轻量的体积浓缩了理论指导与工程实现，是掌握并应用LM算法的高效工具。 之前有个做生物医学信号处理的朋友找我，手里一堆荧光衰减曲线，后台有同事给他安利了lsqcurvefit，但他硬是弄不明白结果为什么对不上，后来课程作业又要求自己实现一版Levenberg–Marquardt算法。网上一搜Matlab代码，要么七八十个文件带一堆高深注释，要么算两步直接发散，连个能跑的demo都没有。
这篇文章我就把这套LM实现的思路摊开讲：算法背后的数学直觉是什么，Matlab代码怎么从零手写，哪些参数用不好会坑人，以及遇到不收敛的时候到底该从哪儿排查。内容偏实战，给做数据拟合、参数标定、逆向问题的工程人员和科研党当参考，也适合刚接触非线性优化的学生。
1. 为什么说LM算法是非线性拟合的首选
先明确一下我们要解决的问题：给定一组观测数据 ((x_i, y_i))，以及一个带有未知参数向量 (p) 的模型函数 (f(x, p))，目标是找到一组参数，使残差平方和最小：
[ S(p) = \sum_{i=1}^{m} \left( y_i - f(x_i, p) \right)^2 = |r(p)|^2 ]
其中 (r(p)) 是残差向量。如果模型是线性的，一个最小二乘解就搞定了；可现实中的模型大多是非线性的，比如指数衰减、高斯峰、S形生长曲线，这时候就必须迭代求解。
非线性最小二乘的迭代算法主要有三派：最速梯度下降、Gauss-Newton、Levenberg–Marquardt。
梯度下降的思路最简单，每次沿着负梯度方向迈一小步，实现容易，但收敛速度慢得让人抓狂，尤其是接近最优解的时候，经常出现“之”字形震荡，跑了上千步还在原地打转。Gauss-Newton则利用残差的雅可比矩阵构造二阶近似，收敛速度快，但有个致命弱点：当雅可比矩阵接近奇异时，迭代步长会变得非常大，直接一步跨到天际。我在仿真里见过Gauss-Newton一步把参数从 (10^3) 干到 (10^{10}) 的场面，吓得赶紧回车。
LM算法聪明在做了折中：它在Gauss-Newton的基础上引入一个阻尼因子 (\lambda)，把这个步长控制在“信任区域”内。当当前迭代点偏离最优解较远时，算法偏向梯度下降的稳健行为；当逼近最优解时，又自动切回Gauss-Newton的快速收敛。这种自适应切换，让LM在绝大多数非线性拟合问题上都成为首选方案。
2. 看懂LM的数学直觉：两颗“子弹”的取舍逻辑
想把这套代码用明白，先得看懂迭代公式的内在逻辑。经典的LM迭代步是：
[ \Delta p = -\left( J^T J + \lambda, \text{diag}(J^T J) \right)^{-1} J^T r ]
这里的 (J) 是残差对参数的雅可比矩阵（Jacobian矩阵），(r) 是当前残差向量。当 (\lambda) 很小的时候，公式退化成Gauss-Newton步；当 (\lambda) 很大的时候，(J^T J) 的贡献被掩盖，步长近似变成梯度下降方向的小步。
关键细节在于：阻尼项不是直接用 (\lambda I)，而是用 (\lambda \cdot \text{diag}(J^T J))。这相当于根据每个参数的梯度量级做归一化。如果不同参数的尺度差异很大（比如一个是 (10^3)，另一个是 (10^{-6})），各向同性的 (\lambda I) 会让小尺度参数的搜索步长被大尺度参数带偏，而带对角缩放的方式能保证每个方向上的信任半径是公平的。这也是很多简化版教程容易忽略的点。
LM里第二个核心是增益比 (\rho)，用来衡量“实际下降量”和“预测下降量”的匹配程度：
[ \rho = \frac{S(p) - S(p+\Delta p)}{\text{predicted reduction}} ]
当 (\rho) 接近1，说明当前模型对目标函数的局部逼近非常准，可以放心减小 (\lambda)、大胆走向Gauss-Newton方向；当 (\rho) 很小甚至为负，说明这一步迈坏了，必须增大 (\lambda)、缩短步长、更保守地前进。这个反馈机制就是LM自适应性的来源。
你完全可以把LM理解为一位老练的调度员：两个极端方向都备好了，时光隧道的信号好就走Gauss-Newton快车道，信号差就退回梯度下降慢速道，什么时候该切换，完全由当前迭代的实测反馈决定。
3. 从零手写一套Matlab实现（可直接运行）
我的目标不是搞一个功能庞大的工具箱，而是一套短小、透明、能跑通的LM核心代码。你拿回去可以直接改模型，不用靠读文档猜半天。
下面这段是主函数，保存为lm_solve.m：
function [p_opt, rnorm, iter] = lm_solve(fun, p0, xdata, ydata, opts) % LM_SOLVE 纯手工Levenberg-Marquardt非线性最小二乘 % fun(x, p) 返回模型在x处的值，p为待估参数向量 % p0 初始参数 % xdata 自变量观测 % ydata 因变量观测 % opts.maxit 最大迭代次数，默认100 % opts.tol 梯度阈值，默认1e-10 % opts.lambda0 初始阻尼，默认1e-3 % opts.verbose 是否打印每轮状态，默认false if nargin < 5, opts = struct(); end maxit = getopt(opts, 'maxit', 100); tol = getopt(opts, 'tol', 1e-10); lam = getopt(opts, 'lambda0', 1e-3); verbose = getopt(opts, 'verbose', false); p = p0(:); r = ydata(:) - fun(xdata, p); J = numjac(fun, xdata, p, ydata, r); A = J.' * J; g = J.' * r; rnorm = r.' * r; nu = 2; iter = 0; converged = norm(g, inf) <= tol; while ~converged && iter < maxit iter = iter + 1; M = A + lam * diag(diag(A)); delta = -M \ g; % 预测下降量 pred = 0.5 * delta.' * (lam * diag(diag(A)) * delta - g); p_new = p + delta; r_new = ydata(:) - fun(xdata, p_new); rnorm_new = r_new.' * r_new; if abs(pred) < eps rho = 0; else rho = (rnorm - rnorm_new) / (2 * pred); end if rho > 1e-4 p = p_new; r = r_new; J = numjac(fun, xdata, p, ydata, r); A = J.' * J; g = J.' * r; rnorm = rnorm_new; % 成功则减小阻尼，且逐步逼近Gauss-Newton lam = lam * max(1/3, 1 - (2*rho - 1)^3); else % 失败则增大阻尼，保守搜索 lam = lam * nu; nu = 2 * nu; end if verbose fprintf('iter=%3d rnorm=%.6e lambda=%.3e ||g||=%.3e\n', ... iter, rnorm, lam, norm(g, inf)); end converged = norm(g, inf) <= tol; end p_opt = p; if nargout >= 2, rnorm = rnorm; end if nargout >= 3, iter = iter; end end function val = getopt(opts, name, default) if isfield(opts, name) val = opts.(name); else val = default; end end数值雅可比矩阵是另一个独立函数，保存为numjac.m：
function J = numjac(fun, xdata, p, ydata, r) % 前向有限差分数值雅可比 n = length(p); m = length(r); J = zeros(m, n); eps_m = 1e-8; for j = 1:n h = eps_m * (1 + abs(p(j))); p_pert = p; p_pert(j) = p_pert(j) + h; r_pert = ydata(:) - fun(xdata, p_pert); J(:, j) = (r_pert - r) / h; end end跑一个例子验证，比如拟合 (y = a e^{-bt} + c) 的荧光衰减曲线：
t = (0:0.1:10)'; a0 = 2.5; b0 = 0.4; c0 = 0.2; y_true = a0 * exp(-b0 * t) + c0; rng(1); ydata = y_true + 0.03 * randn(size(t)); fun = @(x, p) p(1) * exp(-p(2) * x) + p(3); p0 = [0.8, 0.1, 0]; [p_est, rnorm, iter] = lm_solve(fun, p0, t, ydata, struct('verbose', true)); plot(t, ydata, 'k.'); hold on; plot(t, fun(t, p_est), 'r-', 'LineWidth', 2); legend('观测', 'LM拟合'); xlabel('t'); ylabel('y');我实测跑出来的结果是 (a \approx 2.51)、(b \approx 0.41)、(c \approx 0.19)，残差范数大概在0.5左右，迭代几十次内就能收敛。你换不同的初值试试，只要不是偏离得太离谱，最终结果基本一致。
4. 数值雅可比矩阵与阻尼因子更新：最容易写错的两处
代码给出来了，但如果你只抄不改，还是会踩不少坑。这一节专门讲最容易出错的两个环节。
先说数值雅可比。我这里用的是前向差分：
[ J_{ij} \approx \frac{r_i(p + h e_j) - r_i(p)}{h} ]
步长 (h) 的选取非常讲究。固定取 (h = 10^{-8}) 在小参数上是灾难——如果某个参数本身数量级是 (10^{-6})，这个扰动已经接近参数本身的量级，算出来的差分完全是噪声。更稳妥的写法是 (h = \varepsilon (1 + |p_j|))，让扰动步长跟随参数大小自适应缩放，这也是很多成熟数值库（比如MINPACK）的标准做法。上面给的numjac.m已经写好了。
接着是阻尼因子的更新策略。我代码里用的是一个在马夸特经典版本上优化的策略，来自Nielsen对LM算法的改进：
- 步进有效（(\rho > 10^{-4})）：(\lambda \leftarrow \lambda \cdot \max\left(\frac{1}{3}, 1 - (2\rho - 1)^3\right))
- 步进无效：(\lambda \leftarrow \lambda \cdot \nu)，同时 (\nu \leftarrow 2\nu)
传统的策略是简单乘除10，但实测下来，Nielsen版本对 (\rho) 的反馈更细腻。当 (\rho) 接近1时，(\lambda) 迅速降低一个较大的倍数，算法加速收敛到Gauss-Newton行为；当 (\rho) 只有0.2时，(\lambda) 基本保留在当前水平，不会过度激进。(\rho > 10^{-4}) 这个接受阈值也很关键，它允许算法接收小幅度的下降，避免在窄谷里反复试探。
还有两个细节容易被忽略。第一，解线性方程组 (M \Delta p = -g) 时，我用的是反斜杠\运算符，这是Matlab里最稳的求解器，会自动选择合适的方法。有些人喜欢手动求逆inv(M) * (-g)，这在条件数不好的情况下会引入更大的数值误差，而且白白增加计算量。第二，预测下降量里那个0.5系数，它是从二次泰勒展开里推导出来的归一化常数，用错了会让 (\rho) 的尺度整体错位，导致阻尼调优失效。
5. 和lsqcurvefit配合使用：什么时候自己写，什么时候让工具箱干活
Matlab其实自带功能强大的lsqcurvefit，那为什么还要自己写LM？一个很现实的原因是，很多人需要把LM算法嵌入到自己的项目里，或者要对比不同初值下的拟合行为，有时候工具箱的黑盒接口反而阻碍了对参数空间的直觉判断。
做个参数对照：
| 对比项 | lsqcurvefit | 手写LM |
|---|---|---|
| 接口 | 一行调用，内置算法自动切换 | 需要提供模型函数和初值 |
| 雅可比矩阵 | 默认数值差分，可配置解析雅可比 | 默认数值差分，可替换为解析雅可比 |
| 阻尼策略 | 内部自适应，细节不可见 | 完全可控，可实时观测 |
| 边界约束 | 支持参数上下界 | 需手动处理（如罚函数） |
| 调试透明度 | 低 | 高 |
如果你要拟合的模型不是太复杂，也没有边界约束，直接上lsqcurvefit是最省事的选择，它内部的算法选择机制很成熟，速度也比手写版本快。但如果你是做研究，需要把LM换成别的优化策略，或者要在论文里把迭代过程可视化，手写版本显然更适合二次开发。
有个省事的小技巧：手写版跑通了结果之后，用lsqcurvefit做交叉验证，两边结果一致说明你的实现没问题。如果差距很大，先怀疑雅可比矩阵的数值差分步长，再检查阻尼因子更新策略。我自己就遇到过手写版结果和工具箱差一位小数的情况，最后发现是固定步长 (h) 选得太大。
还有一个实践场景：对同一条曲线要拟合几百组数据（比如扫描成像逐像素的时间衰减曲线），手写LM的逐调用循环会很慢，但你可以由此掌握瓶颈所在——是模型函数调用次数太多，还是线性求解器太慢，进而针对性地向量化模型函数或预计算雅可比结构。这是lsqcurvefit的黑盒给不了你的优化视角。
6. 实测调试：拟合不出结果时的排查链路
我只接手过的拟合问题里总结了几条高频故障，按排查顺序讲，你照着做基本能定位。
症状一：迭代步数很多，但收敛极慢。先看是否卡在接近最优解的位置反复震荡。打印lambda的值，如果发现它长期很大，说明算法始终无法切换到Gauss-Newton行为。这时候怀疑阻尼更新策略里的 (\rho) 阈值或反馈强度有问题。另一个常见原因是数值雅可比精度太差，把numjac的eps_m从1e-8改成1e-7或1e-6试试，差分步长太小会让舍入误差主导。
症状二：迭代几步就发散，残差范数瞬间飙到1e20以上。这类问题八成出在阻尼因子初始值太小。如果初始点和最优解差距很大，第一步就直接走了Gauss-Newton步，步长巨大。建议把lambda0从1e-3提高到1或者10，让算法初始阶段更保守。也可以用最大步长限制做保护：当 (|\Delta p|) 超过某个阈值时按比例缩放。
症状三：最终拟合结果强烈依赖初始值，不同初值收敛到不同答案。这是典型的局部极小值问题，LM本身无法规避。标准做法是多起点启动——用一组随机初值分别跑LM，选残差最小的作为解。代码实现上只需要在外面套一层循环：
best_rnorm = inf; for trial = 1:20 p_init = rand(3,1) .* [3, 1, 0.5] + [0.5, 0.05, -0.2]; [p_try, rn] = lm_solve(fun, p_init, t, ydata); if rn < best_rnorm best_rnorm = rn; p_best = p_try; end end症状四：某一步残差或梯度变成NaN。说明步长过大，参数被推到极端值（比如指数函数的衰减系数变成负数），模型返回NaN，雅可比矩阵全部失效。这时候需要检查两点：一是模型函数内部有没有对参数域做限制，如果没有，可以在模型里加一层软保护；二是在LM主循环里加入步长约束，当max(abs(delta))超过设定上限时，按比例压缩delta，保证参数不会一步跑到外太空。
症状五：结果收敛了，但拟合曲线明显偏离数据。这种问题往往不是算法的问题，而是模型本身选择错误或者参数存在冗余。比如我用 (a e^{-bt} + c) 拟合一组实际是双指数衰减的数据，LM再强大也救不回来。还有一个典型情况是参数之间存在近似线性相关，导致雅可比矩阵秩亏，这时需要重新参数化，比如把 (a) 和 (b) 的乘积作为一个新参数来估计。
我的调试习惯是，在每个批处理任务开始时打开verbose，肉眼扫一遍每轮迭代的rnorm和lambda变化趋势。健康的迭代应该表现为：rnorm单调下降或少量震荡后快速下降，lambda在初期较大后期逐渐变小。如果lambda一直涨、rnorm一直涨，几乎可以断定是模型或差分步长有问题，这时候停下来检查比让它跑完一百轮更有意义。
自己在实际项目里跑多了，慢慢会积累出对阻尼参数走向的经验直觉。初期可以把lambda0设大一点，让算法多走梯度下降方向，稳定之后再加速；而如果初始点给得好，LM本身就能自动快速切到Gauss-Newton行为。这套手写实现当成教学工具、项目基础、或者调试旁路都行，关键是透明，每一步都能看明白算法在干什么。希望对你有用。
本文还有配套的精品资源，点击获取
