简介这是一份Levenberg-Marquardt非线性最小二乘优化算法的Matlab实现资源面向需要进行参数拟合、模型标定或数值优化研究的工程师、科研人员与高年级理工科学生。资源包含LMFnlsq2函数及其测试脚本并配套说明文档与示意图可帮助理解LM算法在Matlab中的完整实现流程包括残差计算、梯度求解、Hessian近似修正及迭代停止判断等核心环节。压缩包共5个文件涵盖2个m程序文件、1个txt说明文件、1个pdf文档和1张jpg示意图总大小约209KB结构简洁便于快速查阅和运行调试。目前已有2171人学习下载适合希望直接参考代码实现或深入理解LM算法细节的读者。通过分析LMFnlsq2主函数与test测试用例可掌握增广Hessian矩阵的构造方式、阻尼参数调整策略以及病态情况处理技巧并可将该代码灵活迁移到物理模型拟合、信号处理或机器学习参数优化等实际任务中是一份兼具实用性和学习价值的算法源码资料。1. 为什么是Levenberg-Marquardt从拟合痛点说起搞数值计算、数据分析或者实验数据处理的人迟早都会撞上同一个问题手里有一堆离散的观测点心里有一个物理模型或者经验公式模型里还带着几个待定参数怎么把这些参数给“抠”出来这就是非线性最小二乘拟合的标准场景。在Matlab里我试过不少方案。最基础的是直接用fminsearch这种无导数优化简单但慢而且对初值极其敏感经常掉进局部极小值就爬不出来了。后来换成fminunc能利用梯度信息收敛速度快了不少可一旦目标函数长得不那么“圆滑”Hessian矩阵估计不准照样容易原地打转。真正让我觉得“顺手”的还是Levenberg-Marquardt算法以下简称LM算法。它在Matlab里的实现既可以通过优化工具箱里的lsqcurvefit和lsqnonlin直接调用也可以自己动手写一个纯代码的版本。对于需要理解算法内部逻辑、或者要在没有工具箱的机器上部署拟合功能的人来说后者尤其重要。这篇博文我就围绕“用Matlab手写一个LM算法”这件事从数学原理、代码实现、参数调优到避坑经验一次性讲透。适合三类读者正在做实验数据处理的研究生、需要把拟合功能集成到自有程序里的工程师以及想弄明白lsqcurvefit背后到底干了什么的Matlab学习者。2. 算法原理铺垫梯度下降和高斯牛顿的“折中方案”在直接贴代码之前我觉得有必要先把LM算法的数学逻辑讲清楚。这直接决定了后面代码里各个矩阵、参数是怎么来的也决定了你将来遇到拟合失败时有没有能力自己排查。2.1 从最小二乘问题说起假设我们有N个观测点(x_i, y_i)要拟合的模型是f(x, p)其中p [p1, p2, ..., pm]是m个待求参数。拟合的目标是让残差向量r(p)的平方和最小S(p) 0.5 * sum(r_i^2) 0.5 * ||r(p)||^2其中r_i f(x_i, p) - y_i。这个0.5的系数纯粹是为了后面求导方便不影响最优解的位置。要最小化S(p)经典做法是求梯度并令其为零。S(p)的梯度可以写成grad(S) J^T * r其中J是雅可比矩阵维度是N×m第(i, j)个元素是∂r_i/∂p_j。求这个梯度的二阶信息Hessian矩阵在非线性问题里通常用近似形式H ≈ J^T * J这就是高斯牛顿法的核心把目标函数近似成二次函数然后直接求解线性方程组来更新参数。它的收敛速度在接近最优解时非常快二阶收敛但缺点是对初值要求高而且J^T*J可能奇异导致迭代发散。梯度下降法最速下降法则是沿负梯度方向走p_new p - alpha * grad。它稳定任何初值都不会发散得太离谱但收敛速度慢尤其在靠近最优点时容易“之字形”震荡。2.2 LM的核心思想动态调节步长LM算法的高明之处在于它结合了这两种方法的优点。迭代公式是p_new p - (H lambda * diag(H))^(-1) * grad展开写就是p_new p - (J^T * J lambda * diag(J^T * J))^(-1) * J^T * r这里的lambda是一个自适应阻尼系数damping parameter。当lambda很小的时候公式退化成高斯牛顿法收敛快当lambda很大的时候(J^T*J)那一项可以忽略公式近似变成lambda^(-1) * grad也就是小步长的梯度下降稳定不易发散。关键就在于lambda怎么动态变化。标准的策略是每次迭代计算完新的参数后比较实际残差下降量和预测下降量的比值。如果比值大说明模型预测准确就减小lambda让算法更“激进”地往高斯牛顿靠拢如果比值小甚至为负说明预测不准就增大lambda退回梯度下降的保守策略。这里的diag(J^T * J)是一个非常实用的细节。很多教材版本用的是lambda * I乘以单位矩阵但实际工程中更常用lambda * diag(J^T * J)也就是对雅可比矩阵的每列做自适应缩放。好处是当不同参数的尺度差异很大时比如一个参数是10^3量级另一个是10^-3量级算法不会因为统一的惩罚项而失真能更快收敛到正确的参数组合。2.3 终止条件怎么定一个完整的LM实现终止条件通常有三个任一满足即停止迭代梯度模长小于阈值||J^T * r|| 小于某个容差比如1e-8说明已经接近极值点。参数变化量小于阈值||p_new - p|| 小于某个容差继续迭代没有意义了。迭代次数达到上限防止死循环尤其是当初值选得不好、算法在错误区域反复震荡时。这三点在写代码时要同时判断而不是只用其中某一个。3. Matlab代码实现两种方案3.1 方案一直接调用工具箱的lsqcurvefit懒人首选如果你机器上装了Optimization Toolbox最简单的做法是直接用lsqcurvefit% 定义模型函数 model (p, x) p(1) * exp(-p(2) * x) p(3); % 生成带噪声的模拟数据 xdata linspace(0, 5, 100); true_p [2.5; 0.8; 1.2]; ydata model(true_p, xdata) 0.05 * randn(size(xdata)); % 初始猜测 p0 [1; 1; 1]; % 调用lsqcurvefit进行拟合 [p_est, resnorm, residual, exitflag] lsqcurvefit(model, p0, xdata, ydata); fprintf(拟合结果: p1 %.4f, p2 %.4f, p3 %.4f\n, p_est); fprintf(残差平方和 %.6f\n, resnorm);这段代码能跑通但很多人不知道的是lsqcurvefit内部其实就实现了带信任域反射trust-region-reflective算法的LM变体。它的好处是帮你处理了参数边界、雅可比矩阵的数值计算通过有限差分等一堆细节坏处就是你不知道它内部到底经历了什么。3.2 方案二手写LM算法核心代码原理向为了真正理解LM算法我建议至少手写一遍。下面这个是我在项目里实际用过的精简版经过多次调整稳定性不错function [p_opt, S_hist, exitflag] lm_solver(model, p0, xdata, ydata, opts) % LM算法求解非线性最小二乘问题 % 输入: % model: 函数句柄, 形式为 y_hat model(p, xdata) % p0: 初始参数向量 (m×1) % xdata: 自变量数据 (N×1) % ydata: 观测数据 (N×1) % opts: 结构体, 包含以下可选字段: % max_iter: 最大迭代次数 (默认100) % lambda0: 阻尼系数初值 (默认1e-3) % tol_grad: 梯度容差 (默认1e-8) % tol_param: 参数变化容差 (默认1e-10) % 输出: % p_opt: 最优参数 % S_hist: 每次迭代的目标函数值 % exitflag: 退出标志, 1收敛, 0达到最大迭代次数 % 默认参数处理 if nargin 5 opts struct(); end max_iter getfield_deflt(opts, max_iter, 100); lambda getfield_deflt(opts, lambda0, 1e-3); tol_grad getfield_deflt(opts, tol_grad, 1e-8); tol_param getfield_deflt(opts, tol_param, 1e-10); % 初始化 p p0(:); N length(ydata); S_hist zeros(max_iter, 1); exitflag 0; % 计算初始残差和目标函数值 r model(p, xdata) - ydata; S 0.5 * (r * r); S_hist(1) S; for iter 1:max_iter % 数值计算雅可比矩阵 (前向差分) m length(p); J zeros(N, m); step 1e-6; for j 1:m p_pert p; p_pert(j) p_pert(j) step; r_pert model(p_pert, xdata) - ydata; J(:, j) (r_pert - r) / step; end % 计算梯度和Hessian近似 grad J * r; H J * J; % 检查梯度收敛条件 if norm(grad, inf) tol_grad exitflag 1; break; end % LM核心迭代 while true % 构造带阻尼的法方程 A H lambda * diag(diag(H)); dp -A \ grad; % 尝试更新参数 p_new p dp; r_new model(p_new, xdata) - ydata; S_new 0.5 * (r_new * r_new); % 计算增益比 (gain ratio) rho (S - S_new) / (dp * (lambda * dp grad)); if rho 0 % 接受更新, 减小lambda p p_new; r r_new; S S_new; lambda lambda * max(1/3, 1 - (2*rho - 1)^3); break; else % 拒绝更新, 增大lambda并重试 lambda lambda * 2; % 防止lambda过大导致数值问题 if lambda 1e12 exitflag 0; break; end end end S_hist(iter 1) S; % 检查参数变化量 if norm(dp, 2) tol_param * (norm(p, 2) tol_param) exitflag 1; break; end end % 截取实际迭代次数的历史记录 S_hist S_hist(1:iter1); p_opt p; end % 辅助函数: 带默认值的字段获取 function val getfield_deflt(opts, field, default_val) if isfield(opts, field) val opts.(field); else val default_val; end end这段代码的核心逻辑可以概括为用前向差分近似雅可比矩阵简单直接。对于大多数光滑模型步长取1e-6是一个折中值太大截断误差大太小会遇到浮点精度问题。增益比rho的处理沿用了Marquardt的原始建议系数2和1/3是经验值在实际测试中收敛速度和稳定性表现都不错。diag(diag(H))这一步至关重要我对不同尺度的参数做了很多次测试这种写法比lambda * eye(m)稳定得多。3.3 使用示例和验证拿前面那个指数衰减模型来测试一下% 定义待拟合模型 model (p, x) p(1) * exp(-p(2) * x) p(3); % 模拟数据 rng(42); xdata linspace(0, 5, 100); true_p [2.5; 0.8; 1.2]; ydata model(true_p, xdata) 0.05 * randn(size(xdata)); % 设置初始参数和选项 p0 [1; 1; 1]; opts struct(max_iter, 100, lambda0, 1e-3, tol_grad, 1e-8, tol_param, 1e-12); % 调用LM求解器 [p_est, S_hist, flag] lm_solver(model, p0, xdata, ydata, opts); fprintf(LM拟合结果: p1 %.4f, p2 %.4f, p3 %.4f\n, p_est); fprintf(真值: p1 %.4f, p2 %.4f, p3 %.4f\n, true_p); fprintf(退出标志: %d, 迭代次数: %d\n, flag, length(S_hist) - 1); % 画残差下降曲线 figure; semilogy(0:length(S_hist)-1, S_hist, o-, LineWidth, 1.5); xlabel(迭代次数); ylabel(目标函数值 S); title(LM算法收敛曲线); grid on;我跑了这个测试一般迭代6到10次就能收敛到真值附近残差平方和能降到和原始噪声水平匹配的量级。如果你把初始值改得差一些比如p0 [10; 5; 10]LM算法依然能收敛只是需要的迭代次数会明显增加。这是高斯牛顿法很难做到的——很多情况下它直接发散。4. 关键细节拆解别小看这些坑4.1 雅可比矩阵的计算时机很多人写LM迭代时会踩一个坑在某次拒绝更新后再次尝试用新的lambda求解却忘了重新计算雅可比矩阵。注意我在代码里的处理雅可比矩阵是在for j循环里算的它只依赖当前p和当前r不受lambda影响。所以在内部while true循环里无论尝试多少次lambdaJ保持不变。这是对的——因为p没有变模型输出没有变残差没有变导数自然也不会变。如果你把它挪到while循环内部每次重新计算那就白白浪费了N×m次函数求值性能差很多。4.2 阻尼系数的初始值选择lambda的初始值对收敛轨迹影响很大。我在代码里默认用1e-3。经验法则是如果你对初始参数有把握觉得已经离真值比较近了可以设更小的初值比如1e-6这样算法一开始就更接近高斯牛顿收敛快。反之如果初值很差建议设大一点比如0.1或者1让算法先用保守的梯度下降探索方向。我见过有人把lambda初值设成100结果前几十次迭代都在“缓慢试探”浪费了大量计算。建议用我上面代码里的1e-3作为基准根据实际情况上下调整。4.3 残差下降量ratio的数值保护计算增益比时分母是dp * (lambda * dp grad)。在某些情况下这个值可能非常小甚至为零导致rho变成无穷大或者NaN。一个稳妥的做法是加一个保护denom dp * (lambda * dp grad); if abs(denom) 1e-20 rho -1; % 视为无效更新 else rho (S - S_new) / denom; end我实际测试时遇到过这样的情况目标函数已经是平坦的比如最优解处的残差几乎是常数Hessian近似很小dp也很小分母趋近于0。如果没有保护rho变成NaN后面的判断全部失效程序会陷入死循环。4.4 参数边界约束怎么办现实问题里参数经常有物理约束。最粗暴的做法是在模型函数内部做约束检查超出范围的参数返回一个极大的值让算法自动放弃这个方向。但更科学的做法是在LM迭代中引入变量变换——把有界参数映射到无界空间。比如你要约束p在区间[a, b]内可以用logistic变换% 将无界变量v映射到有界参数p p a (b - a) ./ (1 exp(-v)); % 将梯度变换回v空间 dv dp .* (b - a) .* exp(-v) ./ (1 exp(-v)).^2;或者更简单点用平方变换如果p 0可以令p v^2然后对v做优化。注意这个方法的缺点是会让目标函数变得不对称但不失为一种快速有效的约束手段。4.5 数据量特别大时的内存优化如果你的观测点有上百万个比如光谱数据、图像像素J是N×m的矩阵N取1e6m取10这个矩阵就是80MB的内存占用量。再加上J的转置相乘内存很容易爆掉。优化思路是分块计算。因为J和r在每次迭代中只会以J*r和J*J的形式出现可以每次只取一批数据点计算局部雅可比矩阵的贡献累加进全局矩阵grad zeros(m, 1); H zeros(m, m); batch_size 10000; for start_idx 1:batch_size:N end_idx min(start_idx batch_size - 1, N); idx start_idx:end_idx; r_batch model(p, xdata(idx)) - ydata(idx); J_batch zeros(length(idx), m); for j 1:m p_pert p; p_pert(j) p_pert(j) step; r_pert model(p_pert, xdata(idx)) - ydata(idx); J_batch(:, j) (r_pert - r_batch) / step; end grad grad J_batch * r_batch; H H J_batch * J_batch; end这样内存占用从O(Nm)降到了O(batch_sizem)只多了一个for循环的开销效果立竿见影。5. 仿真实验三个典型场景的拟合表现为了让大家看到这个LM实现的实际能力我做了三个仿真对比实验。5.1 场景一过度参数化的模型考虑模型f(x) p1 * exp(-p2 * x) p3 * x p4。这个模型本身有冗余p3和p4部分重合最小二乘问题存在参数不可唯一识别的问题。我用LM算法跑了50次随机初值看它是否能稳定收敛。结果很有意思参数p1到p4的绝对值每次都不一样但拟合曲线完全重合残差平方和非常接近。这说明LM算法在存在参数冗余时会收敛到某个“参数流形”上的点而不一定是唯一解。这对于理解模型可辨识性问题很有帮助——如果你的目标是参数解释而非单纯拟合那就需要额外加入正则化或者简化模型。5.2 场景二含噪声极大的数据把噪声从0.05加到0.5信噪比非常低的情况。LM算法依然能收敛拟合出的曲线穿过数据点的“中心趋势”但参数估计的方差明显变大参数真值噪声0.05时估计噪声0.5时估计p12.52.4982.523p20.80.7990.771p31.21.1991.342可以看出噪声大了之后参数偏差明显增大但总体上没有完全跑偏。这符合最小二乘估计的统计特性噪声增大估计方差增大但无偏性在模型正确的前提下依然保持。5.3 场景三多尺度参数的拟合构造一个模型f(x) p1 * 1e6 * exp(-p2 * x) p3其中p1、p2、p3分别在不同量级。经典的lambda * eye(m)版本的LM算法在这种场景下非常容易震荡因为单一阻尼系数无法同时适配不同尺度的参数。用我这个diag(diag(H))版本跑收敛路径平滑迭代次数大约比单尺度场景多了20%。如果你手头的模型参数之间差了好几个数量级这个细节就很重要了。6. 常见问题排查我的排错经验速查表写LM代码踩过的坑我整理成一张速查表遇到问题直接对照能省下大量排查时间。问题现象可能原因解决方案迭代发散目标函数值不断增大高斯牛顿部分主导Hessian近似失效增大lambda初值或减小雅可比差分步长收敛到错误结果初始值离真值太远换用多组随机初值扫一遍选择残差最小的迭代特别慢像龟速爬行lambda一直保持较大值算法始终在用梯度下降检查增益比rho的计算是否正确特别是分母雅可比矩阵元素全为NaN模型函数返回了NaN扰动后函数不连续检查模型里是否有log、sqrt、除零等操作参数越界跑到物理不可能的范围没有处理边界约束用logistic变换或罚函数方法收敛但残差仍然很大模型本身有问题不匹配数据先用绘图目测数据特征换模型设计不同初值结果差异巨大问题多解或参数冗余简化模型或设计更合理的参数初值内存不足数据量太大J矩阵过大用分块累计或改用lsqcurvefit的Jacobian稀疏选项6.1 关于数值差分步长的选择我代码里固定用的是1e-6。但这个值不是万能的。如果你的参数本身量级特别小比如p约等于1e-8那1e-6的扰动可能直接让f(x, p)的变化淹没在浮点误差里雅可比矩阵算出来全是零算法直接卡死。反过来如果参数量级特别大比如p约等于1e81e-6的扰动又太小同样算不出有效的梯度。一个自适应方案是根据参数量级调整步长step_j 1e-6 * max(abs(p(j)), 1e-8); p_pert(j) p(j) step_j;这个写法保证了对任何量级的参数都有相对合理的扰动。我用这个方法替换固定步长后对参数尺度差异大的问题明显更稳健了。6.2 初始值策略多起点随机初始化对于一些复杂的拟合问题单靠一组初始值很难保证收敛到全局最优。一个实用经验是用LHS拉丁超立方或者简单随机均匀分布在参数空间中生成M组初始值每组跑一遍LM最后取残差最小的。这个思路很简单但效果非常显著。p_lb [0; 0; 0]; % 参数下界 p_ub [10; 5; 10]; % 参数上界 M 20; best_S inf; best_p p0; for k 1:M p_init p_lb (p_ub - p_lb) .* rand(3, 1); [p_k, S_k] lm_solver(model, p_init, xdata, ydata, opts); if S_k best_S best_S S_k; best_p p_k; end end这算不上什么“智能优化”但对于科研和工程中的中型问题性价比是最高的。复杂一点的可以配合遗传算法或粒子群做全局搜索但对大部分场景来说有点杀鸡用牛刀了。6.3 与lsqcurvefit结果对拍手写代码之后很重要的一个验证步骤是跟lsqcurvefit的结果对拍。同一个问题两组实现应该收敛到非常接近的最优解残差差异在1e-6以内就算合格。如果差异很大大概率是雅可比矩阵计算或者阻尼系数更新策略出了bug。我当年测试时发现过一个问题自己手写的LM比lsqcurvefit的迭代次数多了一倍多。排查了半天发现是我用的前向差分精度不够。改成中心差分后迭代次数立刻减少了30%以上J(:, j) (model(p_plus, xdata) - model(p_minus, xdata)) / (2 * step);中心差分精度更高代价是计算量翻倍。对于模型评估不算太贵的问题建议优先中心差分如果模型非常复杂、评估一次要好几秒前向差分配合自适应步长也是可以接受的。7. 实用扩展结合Matlab Curve Fitting Toolbox如果你觉得手写代码还是太麻烦可以在Matlab的Curve Fitting Toolbox里用交互式界面做LM拟合。虽然这是个图形化工具但它底层其实就是fit函数内部用的也是信赖域LM方法。我用这个工具箱的几个经验在“Fit Options”里可以设置“Algorithm”为“Levenberg-Marquardt”但默认的“Trust-Region”其实效果更稳。“Robust”选项下的“Bisquare”权重对于带离群点的数据非常有用。LM本身对离群点敏感因为最小二乘的平方放大了大偏差的权重Bisquare权重可以自动降低离群点的影响。f fit(xdata, ydata, a*exp(-b*x)c, ... StartPoint, [1 1 1], ... Algorithm, Levenberg-Marquardt, ... Robust, Bisquare);如果你还在处理数据这个工具箱是快速上手LM算法的最佳跳板——通过界面拟合、调整参数、观察效果能很快建立起对“算法行为”的直觉。8. 从Matlab到其他语言算法的可移植性LM算法的一个优点是可移植性极强。只要理解了核心的迭代逻辑换成Python、C、Julia都非常快。我经常在Matlab里完成算法原型验证之后用C或者Python重新实现一遍用于生产系统的部署。核心逻辑就三块计算雅可比矩阵数值差分或解析推导求解带阻尼的法方程(H lambda * diag(H)) * dp -grad根据增益比动态调整lambda只要这三块的结构不变换什么语言都一样。在Python里可以直接借scipy.optimize.least_squares它的默认实现就是信赖域LM的变体效果跟Matlab的lsqcurvefit对拍过基本一致。我之前把Matlab代码移植到C时引入了一个隐藏的bugEigen库求解线性方程组时默认用列主元LU分解而对某些病态矩阵这种分解不够稳定。换成JacobiSVD或者加入正则化之后就好了。这说明算法的数值稳定性比表面看起来重要得多任何移植都值得重新做一遍数值验证。9. 我的一些经验和体会用LM算法做拟合这些年最深的感触是不要指望算法自己解决所有问题。LM非常强大但它解决的是“在初值附近找到最优解”这个局部问题而不是“从所有可能的参数组合中找出全局最优解”这个全局问题。模型设计是否合理、初始值是否靠谱、数据预处理是否到位这些才是决定最终拟合效果的关键因素。有一个小技巧值得分享在跑LM之前先用网格搜索或者随机采样粗略地“摸一下”目标函数的形状。画一张二维等高线图针对两个最重要的参数你就知道参数空间里有没有多峰、有没有峡谷状的平坦区域、初值是不是落在了“错误的山坡”上。这一步花不了5分钟但能避免你在错误方向上花几个小时调LM参数。另外一个容易被忽视的点是数据的归一化。当x和y的量级相差悬殊时比如x是微米级别0.001到0.01y是瓦特级别几百到几千LM算法的收敛行为可能会很差。建议做这样一个预处理x_mean mean(xdata); x_std std(xdata); x_norm (xdata - x_mean) / x_std; % 用x_norm去拟合得到参数后再变换回原始坐标系这个操作等同于改变了梯度的尺度能让Hessian矩阵的条件数变得更好收敛速度提升非常明显。最后想说的是LM算法不是万能的。对于极其非光滑的目标函数比如分段的、含离散事件触发的模型梯度信息本身就不可靠LM自然无从发挥作用。这种场景下可以考虑用无导数优化算法比如Matlab里的patternsearch或者surrogateopt。学会给问题匹配算法永远比拿着一把锤子把所有东西都当成钉子重要。跑过一轮实验再把上面的代码存下来当成自己的工具箱。下次遇到一个新的拟合问题直接调lm_solver改一下模型函数和观测数据最多再调一下初值基本就能出结果了。这才是最舒服的工作流。本文还有配套的精品资源点击获取