简介这份资源提供Levenberg-MarquardtLM算法的Matlab实现代码面向需要解决非线性最小二乘问题的科研人员、工程师及学生可用于物理模型拟合、信号处理、机器学习参数优化等场景。LM算法融合梯度下降与牛顿法优势通过迭代调整参数最小化残差平方和在收敛速度与稳定性间取得平衡。压缩包共5个文件包含2个m脚本文件、1个txt说明、1个pdf文档和1个jpg图片整体约209KB其中m文件为核心算法与测试代码pdf与txt辅助理解算法原理及使用方式。目前已有2173人学习下载说明该实现具有一定参考价值。读者可从中获取完整的LM算法迭代流程包括残差计算、梯度求解、增广Hessian矩阵更新、线性方程组求解及参数更新等关键步骤并借助测试脚本验证算法效果适合作为非线性优化学习与项目开发的实用参考。1. 从一次拟合发散说起Levenberg–Marquardt 的 Matlab 代码到底解决什么问题做非线性最小二乘拟合的人迟早会撞上同一个场景高斯-牛顿法在初值稍偏时直接发散梯度下降法收敛又慢得让人想砸键盘。Levenberg–Marquardt简称 LM就是在这两个极端之间找平衡的算法——它用阻尼因子在“高斯-牛顿步”和“梯度下降步”之间连续切换初值差的时候偏向梯度方向保稳定接近最优解时偏向高斯-牛顿方向保速度。Matlab 里lsqcurvefit、lsqnonlin底层用的就是这类思想但很多人只会调函数不知道内部在干什么一旦拟合失败就完全无从下手。这篇笔记面向两类人一是需要自己写 LM 代码嵌入到 Simulink 或生成 C 代码的工程师二是想搞懂lsqcurvefit为什么报“局部极小值”警告的开发者。我会从雅可比矩阵怎么算、阻尼因子怎么更新、终止条件怎么设一路写到完整可复现的 Matlab 实现和踩坑记录。2. LM 算法的数学骨架与 Matlab 实现选型2.1 从高斯-牛顿到 LM阻尼因子到底在干什么非线性最小二乘的目标函数是 $S(\beta) \sum_{i1}^{m} [y_i - f(x_i, \beta)]^2$其中 $\beta$ 是待拟合参数向量。高斯-牛顿法把残差向量 $r(\beta)$ 在当前点线性化得到更新步 $\Delta\beta -(J^T J)^{-1} J^T r$其中 $J$ 是雅可比矩阵。问题在于 $J^T J$ 经常接近奇异或者初值离最优解太远导致线性化失效步长直接飞出去。LM 的做法是给 $J^T J$ 的对角线加一个阻尼项$\Delta\beta -(J^T J \lambda I)^{-1} J^T r$。$\lambda$ 大时$J^T J$ 被压制步长方向趋近于梯度下降方向 $J^T r$稳定但慢$\lambda$ 小时退化为高斯-牛顿步快但可能发散。关键在于 $\lambda$ 不是固定的而是根据每次迭代的实际效果动态调整如果新参数让残差平方和下降说明当前近似不错减小 $\lambda$更信任高斯-牛顿方向如果残差反而上升说明步子迈大了增大 $\lambda$退回梯度方向并重新计算。这个逻辑用伪代码表达就是初始化 beta, lambda 计算 r, S while 未收敛: 计算 J 求解 (JJ lambda*I) * delta -Jr 计算 beta_new beta delta 计算 S_new if S_new S: beta beta_new lambda lambda / nu % nu 通常取 10 S S_new else: lambda lambda * nu % 拒绝这一步增大阻尼Matlab 里实现这个循环核心难点不在循环本身而在雅可比矩阵的计算和 $\lambda$ 的更新策略。2.2 雅可比矩阵的三种算法与选型依据雅可比矩阵 $J_{ij} \partial r_i / \partial \beta_j$ 的计算方式直接决定代码的通用性和速度。常见做法有三种解析法手动推导偏导数公式写成 Matlab 函数。速度最快、精度最高但每换一个模型就要重新推导适合模型固定、调用频繁的场景。比如拟合 $f a \cdot e^{b x}$残差 $r_i y_i - a e^{b x_i}$偏导为 $\partial r_i/\partial a -e^{b x_i}$$\partial r_i/\partial b -a x_i e^{b x_i}$。有限差分法用 $J_{ij} \approx [r_i(\beta h e_j) - r_i(\beta)] / h$ 近似。通用性最强任何模型都能用但计算量是解析法的 $n1$ 倍$n$ 是参数个数且步长 $h$ 的选取有讲究——太大截断误差高太小舍入误差大。Matlab 的lsqcurvefit默认用有限差分步长取 $\sqrt{\text{eps}} \cdot \max(|\beta_j|, 1)$。复数步长法$J_{ij} \approx \text{Im}[r_i(\beta i h e_j)] / h$精度比有限差分高一个量级没有相减抵消误差但要求函数支持复数输入。很多工程模型含abs、max、查表不支持所以用得少。我一般会这样选模型固定且对速度有要求用解析法模型经常换、开发阶段用有限差分函数本身是解析表达式且支持复数可以试试复数步长法。下面给出有限差分法的 Matlab 实现因为它最通用读者拿去就能跑。function J jacobian_fd(func, beta, xdata, ydata, h) % jacobian_fd 用前向有限差分计算雅可比矩阵 % func: 模型函数句柄签名 func(beta, xdata) % beta: 当前参数向量 (n x 1) % xdata: 自变量 (m x 1) % ydata: 观测值 (m x 1) % h: 差分步长缺省用 sqrt(eps)*max(|beta|,1) if nargin 5 h sqrt(eps) * max(abs(beta), 1); end n numel(beta); m numel(ydata); J zeros(m, n); r0 ydata - func(beta, xdata); % 当前残差 for j 1:n beta_p beta; beta_p(j) beta_p(j) h(j); r_p ydata - func(beta_p, xdata); J(:, j) (r_p - r0) / h(j); % 前向差分 end end这段代码的逻辑很直接对每个参数分量做一次微小扰动观察残差向量的变化率。参数说明h的缺省值来自 Matlab 内部对数值微分步长的经验公式max(abs(beta),1)保证参数接近零时步长不会过小。注意h是向量每个参数分量可以有不同的步长。如果模型函数对参数尺度敏感比如一个参数是 1e-6 量级另一个是 1e6 量级建议手动传入归一化后的h。2.3 阻尼因子更新策略为什么你的 LM 收敛像蜗牛$\lambda$ 的更新策略是 LM 代码里最容易被忽视、也最影响收敛速度的部分。教科书上通常写“下降则除以 10上升则乘以 10”但实际工程中这个固定倍数在参数尺度差异大时会翻车。常见做法是引入一个缩放矩阵 $D$把阻尼项写成 $\lambda \cdot \text{diag}(J^T J)$ 而不是 $\lambda I$。这样每个参数方向的阻尼强度自动适应该方向的曲率参数尺度差异大时收敛更稳。Matlab 的lsqcurvefit内部用的就是这种缩放形式。另一个细节是 $\lambda$ 的初始值。取太小第一步就可能发散取太大前几步几乎在原地踏步。我一般取 $\lambda_0 10^{-3} \cdot \max(\text{diag}(J_0^T J_0))$让初始阻尼和问题的曲率量级匹配。如果第一次迭代就失败再乘以 10 重试。function [beta, resnorm, exitflag] lm_solve(func, beta0, xdata, ydata, opts) % lm_solve Levenberg-Marquardt 主循环 % 返回: beta 最优参数, resnorm 残差平方和, exitflag 终止原因 if nargin 5, opts struct(); end maxIter getfield_default(opts, maxIter, 200); tolX getfield_default(opts, tolX, 1e-8); tolFun getfield_default(opts, tolFun, 1e-10); nu getfield_default(opts, nu, 10); % 阻尼增减倍数 beta beta0(:); r ydata(:) - func(beta, xdata); S sum(r.^2); J jacobian_fd(func, beta, xdata, ydata); lambda 1e-3 * max(diag(J * J)); % 初始阻尼与曲率匹配 exitflag 0; for iter 1:maxIter A J * J lambda * diag(diag(J * J)); % 缩放阻尼 g J * r; if rcond(A) 1e-15 lambda lambda * nu; % 矩阵病态加大阻尼 continue; end delta -A \ g; beta_new beta delta; r_new ydata(:) - func(beta_new, xdata); S_new sum(r_new.^2); if S_new S % 接受这一步 if norm(delta) tolX * (norm(beta) tolX) || abs(S - S_new) tolFun beta beta_new; S S_new; exitflag 1; break; end beta beta_new; r r_new; S S_new; J jacobian_fd(func, beta, xdata, ydata); lambda max(lambda / nu, 1e-12); else % 拒绝这一步增大阻尼 lambda lambda * nu; if lambda 1e12 exitflag -1; break; % 阻尼过大无法继续 end end end resnorm S; end function v getfield_default(s, f, d) if isfield(s, f), v s.(f); else, v d; end end这段代码有几个关键点值得展开。第一A J*J lambda*diag(diag(J*J))用的是缩放阻尼不是单位矩阵这是和教科书版本最大的区别。第二rcond(A) 1e-15检查矩阵条件数病态时直接加大阻尼跳过这一步避免求解出天文数字的步长。第三接受步长时先检查收敛条件再更新避免多算一次雅可比。第四lambda下限设为 1e-12防止除零上限 1e12 作为发散保护。参数说明nu取 10 是经验值取 2 到 3 收敛更平滑但迭代次数多取 100 反应快但容易震荡。tolX和tolFun分别控制参数变化量和目标函数变化量的收敛阈值前者比后者更常用。maxIter设 200 对大多数曲线拟合够用如果 200 次还没收敛通常是模型本身有问题加迭代次数只是拖延时间。3. 用 LM 代码拟合实际数据的完整流程3.1 构造测试数据与模型函数光看代码不够得跑一组有噪声的数据才能看出 LM 的真实行为。我构造一个双指数衰减模型 $f a_1 e^{-b_1 x} a_2 e^{-b_2 x}$这是荧光衰减、电池弛豫等场景的常见形式也是出了名的难拟合——两个指数项容易相互补偿雅可比矩阵接近奇异。% 生成测试数据 rng(42); % 固定随机种子保证可复现 xdata linspace(0, 5, 200); true_params [3.0, 2.5, 1.2, 0.4]; % a1, b1, a2, b2 model (beta, x) beta(1)*exp(-beta(2)*x) beta(3)*exp(-beta(4)*x); ydata model(true_params, xdata) 0.02*randn(size(xdata)); % 加 2% 噪声 % 从偏离真值的初值开始拟合 beta0 [1.0, 1.0, 1.0, 1.0]; [beta_fit, resnorm, exitflag] lm_solve(model, beta0, xdata, ydata); fprintf(真实参数: %s\n, mat2str(true_params, 4)); fprintf(拟合参数: %s\n, mat2str(beta_fit, 4)); fprintf(残差平方和: %.6f, 终止标志: %d\n, resnorm, exitflag);跑这段代码你会看到拟合参数和真实参数在小数点后两位基本一致exitflag为 1。但把beta0改成[10, 10, 10, 10]再跑一次可能会发现收敛变慢甚至停在局部极小值。这就是双指数模型对初值敏感的表现也是 LM 代码需要配合多初值重启的原因。3.2 拟合结果的可视化与残差诊断拟合完不能只看参数必须画图。残差图能暴露模型系统性偏差——如果残差呈现规律性波动而不是随机分布说明模型形式选错了再怎么调 LM 参数也没用。% 可视化拟合结果与残差 figure(Position, [100 100 800 500]); subplot(2,1,1); plot(xdata, ydata, ko, MarkerSize, 3, DisplayName, 观测数据); hold on; plot(xdata, model(beta_fit, xdata), r-, LineWidth, 1.5, DisplayName, LM 拟合); xlabel(x); ylabel(y); legend(Location, northeast); title(双指数模型 LM 拟合结果); subplot(2,1,2); residuals ydata - model(beta_fit, xdata); stem(xdata, residuals, b., MarkerSize, 4); hold on; yline(0, k--); xlabel(x); ylabel(残差); title(sprintf(残差分布 (RMS %.4f), sqrt(mean(residuals.^2))));残差图如果呈现“两端大中间小”的喇叭形说明噪声是乘性的而非加性的应该拟合对数残差或者加权。如果残差在某个区间系统性偏正或偏负说明模型缺少一项。这些诊断比看resnorm数值有用得多。3.3 参数不确定度估计协方差矩阵怎么算拟合出参数只是第一步工程上还需要知道参数的置信区间。LM 收敛后协方差矩阵近似为 $\text{Cov}(\beta) \approx \sigma^2 (J^T J)^{-1}$其中 $\sigma^2 S / (m - n)$ 是残差方差$m$ 是数据点数$n$ 是参数个数。% 计算参数协方差与标准误 m numel(ydata); n numel(beta_fit); J_final jacobian_fd(model, beta_fit, xdata, ydata); sigma2 resnorm / (m - n); cov_beta sigma2 * inv(J_final * J_final); std_err sqrt(diag(cov_beta)); fprintf(参数标准误:\n); for i 1:n fprintf( beta(%d) %.4f ± %.4f\n, i, beta_fit(i), std_err(i)); end注意inv(J*J)在参数强相关时数值不稳定更稳的做法是用pinv或者对 $J$ 做 QR 分解后求逆。双指数模型里 $a_1$ 和 $a_2$ 经常强负相关标准误可能很大这本身就是模型可辨识性差的信号不是代码 bug。4. LM 代码的避坑与常见问题排查4.1 现象迭代过程中残差突然变成 NaN原因模型函数在某些参数取值下产生复数或无穷大。比如指数模型 $e^{b x}$ 在 $b$ 很大且 $x$ 很大时溢出或者对数模型参数为负。LM 的步长在初期可能很大把参数推到物理上无意义的区域。解决在模型函数内部加参数边界检查超出范围时返回一个大的惩罚残差而不是 NaN。更规范的做法是用带边界约束的 LM 变体trust-region-reflectiveMatlab 的lsqcurvefit支持lb和ub参数。自己写代码的话可以在接受步长前检查isfinite(S_new)不是有限值就拒绝这一步并增大阻尼。4.2 现象收敛到局部极小值拟合曲线只匹配了部分数据原因LM 是局部优化算法从不同初值出发会收敛到不同局部极小值。双指数、多峰高斯等模型尤其严重。解决多初值重启。用拉丁超立方采样或网格采样生成一批初值对每个初值跑一次 LM取残差最小的结果。参数少的时候3 个以内网格采样够用参数多的时候用sobolset或lhsdesign生成初值。另一个技巧是先固定部分参数拟合再放开全部参数精调。4.3 现象雅可比矩阵计算耗时占整个拟合的 80%原因有限差分法每个参数都要重新计算一次模型函数参数多的时候计算量线性增长。如果模型函数本身还包含数值积分或查表单次拟合可能几分钟。解决能推导解析雅可比的模型一定用解析法。不能推导的考虑并行计算——Matlab 的parfor可以把每个参数方向的差分计算分配到多核。另外如果模型函数支持向量化输入把xdata作为整体传入而不是循环单点计算速度能提升一个量级。4.4 现象参数标准误大得离谱但拟合曲线看起来很好原因参数之间存在强相关性模型过参数化。曲线拟合得好不代表每个参数都有物理意义可能是两个参数在相互补偿。解决看协方差矩阵的非对角元素计算相关系数矩阵。如果某两个参数的相关系数绝对值超过 0.95考虑固定其中一个参数或者重新参数化模型比如用和与差代替两个独立参数。这个问题不是 LM 代码能解决的属于模型可辨识性范畴。4.5 现象rcond(A) 1e-15频繁触发迭代几乎停滞原因$J^T J$ 本身接近奇异加阻尼后条件数仍然很差。常见于参数尺度差异巨大一个参数 1e-6另一个 1e6或者模型对某个参数不敏感。解决对参数做归一化让所有参数在同一个量级。具体做法是定义缩放因子 $s_j$拟合 $\beta_j s_j \cdot \gamma_j$对 $\gamma$ 做 LM 迭代。缩放因子取参数的典型量级或者初值的绝对值。另一个办法是增大初始 $\lambda$让前几步走梯度方向等参数进入合理区域后再减小阻尼。5. 让 LM 代码真正可用的几个进阶技巧5.1 用解析雅可比把速度提上来有限差分法通用但慢。如果模型是双指数这类形式解析雅可比不难推。以 $f a_1 e^{-b_1 x} a_2 e^{-b_2 x}$ 为例残差 $r_i y_i - f_i$偏导为function J jacobian_analytic(beta, xdata, ydata) % 双指数模型的解析雅可比 a1 beta(1); b1 beta(2); a2 beta(3); b2 beta(4); e1 exp(-b1 * xdata); e2 exp(-b2 * xdata); m numel(xdata); J zeros(m, 4); J(:,1) -e1; % dr/da1 J(:,2) a1 * xdata .* e1; % dr/db1 J(:,3) -e2; % dr/da2 J(:,4) a2 * xdata .* e2; % dr/db2 end把lm_solve里的jacobian_fd调用换成这个函数200 个数据点、4 个参数的拟合耗时能从 15 ms 降到 2 ms 左右。参数越多、模型越复杂差距越大。注意解析雅可比的符号容易推错写完一定要和有限差分的结果对比验证——用同一组参数分别算两个雅可比看最大相对误差是否在 1e-6 以内。5.2 加权最小二乘残差方差不齐时的处理如果观测数据的噪声水平随 $x$ 变化比如计数数据服从泊松分布方差正比于均值普通最小二乘的假设不成立需要加权。加权 LM 的目标函数变成 $S \sum w_i r_i^2$雅可比和梯度相应变为 $J_w W^{1/2} J$$g_w J^T W r$其中 $W \text{diag}(w_i)$。% 加权 LM 的核心修改在 lm_solve 内部 W diag(weights); % weights 是 m x 1 的权重向量 Jw sqrt(W) * J; % 加权雅可比 gw J * W * r; % 加权梯度 A Jw * Jw lambda * diag(diag(Jw * Jw)); delta -A \ gw;权重怎么取计数数据取 $w_i 1 / \max(y_i, 1)$已知测量误差取 $w_i 1 / \sigma_i^2$。权重不需要归一化因为 $\lambda$ 的初始值会根据 $J_w^T J_w$ 的尺度自动调整。5.3 从代码到部署生成不依赖 Matlab 的 C 代码工程上经常需要把拟合好的 LM 代码嵌入到嵌入式设备或 Simulink 模型里。Matlab Coder 可以把lm_solve和模型函数一起转成 C 代码但有几个限制不能用inv、rcond这类函数需要手动替换成\求解或自定义条件数估计动态内存分配要关掉所有数组大小必须编译期确定。我一般会写一个固定参数个数的专用版本把maxIter、tolX等作为宏定义雅可比用解析法矩阵求解用 LAPACK 的dgesv或者自己写一个带部分主元的高斯消元。这样生成的代码体积小、无依赖在 ARM Cortex-M 上跑 4 参数双指数拟合大约 5 ms 一次迭代。5.4 一个容易被忽略的细节终止条件的组合判断很多人只设tolX或者只设tolFun结果要么提前停止参数还在大幅变化但目标函数暂时平了要么空转参数几乎不动但目标函数还在缓慢下降。我习惯用三个条件组合参数变化量小于阈值、目标函数变化量小于阈值、梯度范数小于阈值任意两个满足才终止。梯度范数 $|J^T r|$ 是判断是否到达驻点的最直接指标但它对参数尺度敏感需要归一化后再比较。% 组合终止条件示例 grad_norm norm(J * r) / (1 norm(r)); if (norm(delta) tolX * (norm(beta) tolX)) ... (abs(S - S_new) tolFun * (1 S)) ... (grad_norm 1e-6) exitflag 1; break; end这三个条件里tolX和tolFun的相对形式乘以norm(beta)和S比绝对形式更通用因为不同问题的参数量级和残差量级差异很大。梯度范数阈值 1e-6 对大多数问题够用如果问题本身病态可以放宽到 1e-4。写 LM 代码这些年最大的教训是算法本身只有二十行但让它在你自己的模型上稳定收敛需要调的参数和排查的边界情况远超过算法核心。每次换模型先用有限差分跑通再换解析雅可比提速先用单位阻尼跑再换缩放阻尼先单初值跑再多初值重启。这个顺序能帮你把问题定位在模型本身还是优化器上。希望帮到你。本文还有配套的精品资源点击获取