LM算法原理与MATLAB实现：阻尼最小二乘非线性拟合实战

简介资源提供基于LM算法的非线性最小二乘拟合实现面向需要处理曲线拟合与参数估计的科研、工程及数据建模人员。内容围绕Levenberg-Marquardt算法的核心思想展开包含MATLAB源函数与测试脚本并配有拟合效果图片便于直观理解算法动态调整阻尼因子、平衡梯度下降与高斯-牛顿迭代的过程。压缩包共5个文件主要为3个.m脚本和2张JPG图像总大小46KB结构简洁无需复杂依赖适合快速移植和二次开发。已有1320人学习下载。借助源码可掌握残差平方和最小化的迭代实现细节尤其是在角度测量等非线性关系数据处理中的具体用法图片结果可对照验证参数更新方向与收敛特性对初学者理解稳健拟合策略、避免局部最小值陷阱也有帮助。1. 为什么非线性拟合绕不开LM算法做过实验数据处理的人大多有这种体验模型明明是几十年前教科书里写好的但真拿自己的测量数据去拟合时梯度下降法要么慢得像蜗牛要么在某个坑里震荡高斯-牛顿法收敛快却经常因为雅可比矩阵奇异而直接发散。这时候Levenberg-Marquardt算法是绝大多数工程师的第一选择。它本质上是梯度下降与高斯-牛顿法的自适应混合通过一个阻尼因子在全局搜索和局部快速收敛之间动态切换。我最早是在一个角度传感器标定项目里接触到这套东西的当时手里只有三个MATLAB脚本lm.mLM核心迭代、angelLove.m模型函数、angleLoveForYou.m主程序。打开代码发现它把非线性最小二乘的整个流程压在了不到两百行里很值得拆开讲一讲。本文就用这套代码为线索把LM算法的数学推导、MATLAB实现、参数调节和陷阱一次说清。2. LM算法原理从梯度下降到阻尼最小二乘的数学脉络2.1 非线性最小二乘问题的目标函数与雅可比非线性最小二乘的标准形式是给定一组观测点 $(x_i, y_i)$选择一个含参数的模型 $\hat{y} f(x; \theta)$其中 $\theta \in \mathbb{R}^n$ 是待估参数向量。目标函数是残差平方和$$ S(\theta) \sum_{i1}^{m} r_i(\theta)^2 | r(\theta) |^2,\quad r_i(\theta) y_i - f(x_i; \theta) $$这里 $m$ 是样本数$n$ 是参数个数通常 $m \gg n$。注意残差向量 $r(\theta)$ 是一个从 $\mathbb{R}^n$ 映射到 $\mathbb{R}^m$ 的向量函数。求解目标就是找到 $\theta^*$ 使 $S(\theta)$ 最小。一阶最优条件是梯度为零而梯度的表达式为$$ \nabla S(\theta) 2 J(\theta)^T r(\theta) $$其中 $J(\theta)$ 是 $m \times n$ 的雅可比矩阵第 $i$ 行第 $j$ 列元素为 $\partial r_i / \partial \theta_j$。这个表达式是整个LM算法的出发点。在实际工程中如果模型形式复杂手动求雅可比容易出错所以我在看lm.m时特别注意它是不是用了数值差分。2.2 高斯牛顿法与梯度下降的缺陷高斯牛顿法的思想是对残差做一阶泰勒展开$r(\theta \delta) \approx r(\theta) J(\theta)\delta$代入目标函数后对 $\delta$ 求最小得到正规方程$$ J^T J ,\delta_{gn} -J^T r $$这个方程的解 $\delta_{gn}$ 是高斯牛顿步。优点是接近最优解时收敛极快但 $J^T J$ 可能奇异或病态导致 $\delta_{gn}$ 过大甚至发散。相反最速梯度下降法每一步沿负梯度方向走$$ \delta_{gd} -\nabla S -2J^T r $$虽然保证目标函数下降但收敛速度是线性的在接近极小值时会出现“之”字形爬行。我曾经用纯梯度下降拟合一个带衰减的正弦模型迭代两万次还差一大截换成LM后几十步就收敛了。2.3 LM算法的阻尼因子与步长控制LM算法的核心是引入一个非负阻尼因子 $\lambda$把正规方程修正为$$ (J^T J \lambda I)\delta_{lm} -J^T r $$这里 $I$ 是 $n \times n$ 单位阵。当 $\lambda$ 很大时$\delta_{lm} \approx -\frac{1}{\lambda} J^T r$相当于小步长的梯度下降保证全局搜索当 $\lambda$ 很小时$(J^T J \lambda I) \approx J^T J$退化成高斯牛顿法获得快速局部收敛。关键是 $\lambda$ 的更新策略常见的做法是如果本次迭代后 $S(\theta \delta) S(\theta)$接受新参数同时减小 $\lambda$比如 $\lambda \lambda / 3$因为这说明当前近似较好可以更激进些。如果新参数使目标函数变大拒绝更新增大 $\lambda$比如 $\lambda \lambda \times 2$让步长变小、沿梯度方向走。这个机制可以让算法在“信任域”内自适应调整。我在lm.m里看到它用的是类似Marquardt的改进版对角线上不直接加 $\lambda I$而是加 $\lambda \cdot \operatorname{diag}(J^T J)$。这样做的效果是当某个参数方向上的梯度分量很小时对应的步长不会被过度压缩数值稳定性更好。2.4 算法流程与收敛判据完整的LM迭代流程可以总结为以下步骤这也是我阅读代码时对照的标准给定初始参数 $\theta_0$初始阻尼因子 $\lambda_0$上限 $\lambda_{\max}$迭代次数上限。计算当前残差 $r(\theta_k)$ 和雅可比 $J(\theta_k)$得到近似 Hessian 矩阵 $A J^T J$ 和梯度向量 $g -J^T r$。求解线性方程组 $(A \lambda \operatorname{diag}(A))\delta g$得到参数更新量 $\delta$。计算候选新参数 $\theta_{k1} \theta_k \delta$并计算其目标函数值 $S(\theta_{k1})$。如果 $S(\theta_{k1}) S(\theta_k)$接受更新$\lambda$ 缩小否则拒绝更新$\lambda$ 扩大回到第3步。判断收敛如果 $|g|_{\infty} \epsilon_1$或者 $| \delta | \epsilon_2(|\theta_k| \epsilon_2)$或者达到最大迭代次数则停止。需要注意的是收敛判据里的 $|g|_{\infty}$ 其实反映的是梯度大小梯度过零是一阶必要条件但未必是全局最优。因为LM本质上是局部优化算法初始值选得不好可能收敛到局部极小后面第5章会专门讨论对策。3. MATLAB实现lm.m从零手写LM求解器3.1 lm.m的函数签名与输入输出这里我按照项目中lm.m的典型设计来讲解。lm.m是一个独立的MATLAB函数输入是残差函数句柄、初始参数、迭代选项输出是最优参数、残差范数、迭代信息。其基本签名如下function [theta_opt, S, info] lm(residual_fun, theta0, opt) % residual_fun : 函数句柄输入 (theta, xdata, ydata)返回残差向量 r % theta0 : 初始参数向量列向量 % opt : 结构体包含以下字段可缺省 % opt.tolX : 参数变化量的停止阈值默认 1e-6 % opt.tolFun : 目标函数变化量的停止阈值默认 1e-10 % opt.gradTol : 梯度范数阈值默认 1e-6 % opt.maxIter : 最大迭代次数默认 200 % opt.lambda0 : 初始阻尼因子默认 1e-3 % opt.lambdaMax: 阻尼因子上限默认 1e10 % opt.verbose : 是否打印迭代过程默认 false % 输出 % theta_opt : 最优参数向量 % S : 最优残差平方和 % info : 结构体包含迭代次数、停止原因、历史记录这个接口把残差函数完全交出去使得lm.m可以用于任意非线性模型。注意theta0必须是一个列向量$n \times 1$否则在后续矩阵运算中容易出维度不匹配的错误。我建议在函数内第一行就做检查加上assert(size(theta0,2) 1, theta0 必须为列向量);。3.2 雅可比矩阵的数值差分计算lm.m的核心是每次迭代计算雅可比矩阵。由于模型函数可能很复杂代码里没有让用户提供解析雅可比而是采用了中心差分公式function J numerical_jacobian(residual_fun, theta, xdata, ydata, delta) % 计算残差向量对参数的雅可比矩阵中心差分法 n length(theta); r0 residual_fun(theta, xdata, ydata); % 当前残差向量m×1 m length(r0); J zeros(m, n); for j 1:n step max(delta, delta * abs(theta(j))); % 步长与参数量级相关 theta_plus theta; theta_minus theta; theta_plus(j) theta(j) step; theta_minus(j) theta(j) - step; r_plus residual_fun(theta_plus, xdata, ydata); r_minus residual_fun(theta_minus, xdata, ydata); J(:, j) (r_plus - r_minus) / (2 * step); end end参数说明delta是差分步长的基础值我一般设为sqrt(eps)约等于 1.49e-8。这个值太小会引入浮点误差太大会使差分近似偏差变大。step max(delta, delta * abs(theta(j)))表示当参数的绝对值很大时差分步长按比例增大避免出现theta step theta的静默失效。这是数值计算里一个容易被忽略的细节尤其是在拟合幅值达数千的物理量时。中心差分比前向差分精度高但代价是每列参数要额外调用两次残差函数。如果参数很多计算量会很大。对于耗时模型可以考虑改用解析雅可比或自动微分。不过在当前这个场景下模型函数angelLove.m非常简单数值差分几乎不拖慢速度。3.3 阻尼因子的动态调整策略阻尼因子的更新是LM算法的灵魂。在lm.m中更新逻辑是这样的% 初始化 lambda opt.lambda0; theta theta0(:); S_old sum(residual_fun(theta, xdata, ydata).^2); for k 1:opt.maxIter r residual_fun(theta, xdata, ydata); J numerical_jacobian(residual_fun, theta, xdata, ydata, sqrt(eps)); % m×n A J * J; % n×n 近似Hessian g -J * r; % n×1 梯度向量 M A lambda * diag(diag(A)); % 对角缩放阻尼 % 求解正规方程用 Cholesky 分解替代直接求逆 [L, p] chol(M, lower); if p 0 % M 不是正定矩阵增大阻尼重新求解 lambda lambda * 10; continue; end delta L \ (L \ g); theta_new theta delta; S_new sum(residual_fun(theta_new, xdata, ydata).^2); if S_new S_old % 成功接受更新减小阻尼 theta theta_new; S_old S_new; lambda max(lambda / 3, 1e-12); % 下限保护 else % 失败增大阻尼本次参数不更新 lambda min(lambda * 2, opt.lambdaMax); if lambda opt.lambdaMax info.stopReason lambda上限; break; end end % 判断收敛 if norm(g, inf) opt.gradTol info.stopReason 梯度收敛; break; end if norm(delta, 2) opt.tolX * (norm(theta, 2) opt.tolX) info.stopReason 参数变化量收敛; break; end end逻辑说明用diag(diag(A))构造阻尼矩阵而不是单位阵I。这样做的理由是如果某个参数对应的 $A_{jj}$ 很小说明这个方向上的曲率很小加上一个相同的入会过度限制该方向的步长。按 $A_{jj}$ 缩放后每个参数方向上的阻尼效果相对均衡。求解线性方程组用 Cholesky 分解。M是实对称正定矩阵因为加了阻尼因此 Cholesky 分解比inv(M)*g更快速、数值更稳定。代码里利用chol的第二个返回值p判断分解是否成功。如果p 0说明当前M不正定这时直接放大阻尼重新计算避免数值异常。阻尼因子的上下限保护很关键。下限设为1e-12防止lambda变成 0 导致退化上限是用户设定的opt.lambdaMax通常给 1e10超过这个就停止说明算法已经进入死区。3.4 迭代停止条件与错误处理停止条件需要在实际操作中灵活调整。我在使用lm.m时遇到过几次不收敛的情况归纳下来有两种典型原因第一种是初始阻尼因子lambda0设置不当。如果lambda0太小算法在早期就偏向高斯牛顿法容易出现震荡如果太大收敛速度会被拖慢至纯梯度下降。经验值是取1e-3但如果是高度非线性模型比如指数乘积模型我建议试一下1或10。第二种是差分步长与模型尺度不匹配。在numerical_jacobian中当参数量级很小比如 1e-10 的系数时max(delta, delta*abs(theta))会退化成固定步长sqrt(eps)还算可用。但如果参数量级达到 1e10步长也随之变成 1e-2这时候中心差分的精度不够。所以我在实际代码里加了一个delta字段输入允许调用方根据参数范围调整。下面是一个用lm.m拟合普通指数衰减模型的最小调用示例% 构造测试数据 x linspace(0, 5, 50); y 2.5 * exp(-0.8 * x) 0.1 * randn(50, 1); % 定义残差函数 resfun (theta, x, y) y - theta(1) * exp(-theta(2) * x); % 初始猜测 theta0 [1; 0.5]; % 设置选项 opt.lambda0 0.01; opt.maxIter 100; % 调用lm求解 [theta_opt, S] lm(resfun, theta0, x, y, opt);注意resfun的签名中xdata和ydata是额外参数但在数值雅可比函数里xdata和ydata也是透传给residual_fun的因此需要保证残差函数能够正确区分数据和参数。如果你的模型包含变量权重或其他额外信息可以把这些数据打包到一个结构体中用匿名函数捕获。4. 实战angleLoveForYou角度数据的非线性拟合4.1 场景描述与模型选择项目里的angleLoveForYou.m和angelLove.m暗示了一个与角度测量相关的拟合问题。我这里构造一个典型的工程场景一个摆动物体的角度 $\theta(t)$ 随时间衰减模型为正弦衰减函数$$ \theta(t) a \cdot \sin(b \cdot t c) \cdot e^{-d \cdot t} e $$其中符号定义如表格所示参数物理意义估计值范围$a$初始振幅角度5 ~ 50$b$角频率rad/s0.1 ~ 10$c$初始相位rad-$\pi$ ~ $\pi$$d$阻尼系数0 ~ 2$e$角度偏移零漂-5 ~ 5这个模型包含了线性项$e$、非线性项$\sin$ 和 $e^{-dt}$并且参数之间存在耦合非常适合检验LM算法的稳健性。实际测量数据是带噪声的我们先用代码生成仿真数据随后再调用lm.m拟合。4.2 残差函数定义与初始参数猜测angelLove.m在这个项目里充当的就是模型残差函数。它的写法通常是这样function r angelLove(theta, t, y) % 角度衰减模型的残差 % theta [a b c d e] % t, y 为观测时间与角度 a theta(1); b theta(2); c theta(3); d theta(4); e theta(5); y_fit a * sin(b * t c) .* exp(-d * t) e; r y - y_fit; end注意. * exp(-d*t)中的点乘号t是列向量时必须逐元素运算这是MATLAB新手最容易犯的错误。如果漏掉点乘大概率会得到矩阵维度不匹配的报错或者一个错误的矩阵值。初始参数猜测上我的做法是先用快速傅里叶变换估计主频 $b$再根据前几个峰值的衰减率估算 $d$振幅 $a$ 直接用数据的标准差偏移 $e$ 用均值。对于一般数据这个方法能给出相当好的初始点。如果没有FFT条件也可以把 $b$ 粗估为 $2\pi / T$其中 $T$ 是数据的大致周期。4.3 完整调用代码与可视化主程序angleLoveForYou.m中做三件事加载数据、调用lm.m、绘制拟合曲线和残差图。以下是一个可复现的完整示例% angleLoveForYou.m 主程序 clear; close all; clc; % 仿真生成带噪声的数据 t (0:0.05:10); % 时间采样单位秒 true_theta [30, 2.0, 0.7, 0.3, 1.5]; % 真实参数 y true_theta(1) * sin(true_theta(2) * t true_theta(3)) .* exp(-true_theta(4)*t) true_theta(5); y y 0.8 * randn(size(t)); % 叠加高斯白噪声 % 定义残差函数直接使用 angelLove.m resfun (theta, x, y) angelLove(theta, x, y); % 设置LM选项 opt.lambda0 0.01; opt.maxIter 200; opt.gradTol 1e-8; opt.verbose true; % 初始参数估计按经验给一个近似 theta0 [25; 1.8; 0.5; 0.2; 1.0]; % 调用lm.m求解 [theta_opt, S, info] lm(resfun, theta0, t, y, opt); % 输出结果 fprintf(拟合参数: a%.3f b%.3f c%.3f d%.3f e%.3f\n, theta_opt); fprintf(残差平方和: %.4f, 迭代次数: %d, 停止原因: %s\n, S, info.iter, info.stopReason); % 绘制拟合曲线 y_fit theta_opt(1) * sin(theta_opt(2) * t theta_opt(3)) .* exp(-theta_opt(4)*t) theta_opt(5); figure(Position,[100,100,900,350]); subplot(1,2,1); plot(t, y, x, MarkerSize, 4); hold on; plot(t, y_fit, r-, LineWidth, 1.5); legend(观测数据,LM拟合曲线,Location,northeast); xlabel(时间 t / s); ylabel(角度 / °); title(非线性最小二乘拟合结果); subplot(1,2,2); stem(t, y - y_fit, MarkerSize, 3); xlabel(时间 t / s); ylabel(残差); title(拟合残差分布); grid on;参数说明opt.verbose true会让lm.m在命令行打印每步迭代的目标函数变化情况方便观察收敛过程。我通常建议第一次运行开启 verbose确认每步目标函数都在下降如果出现连续拒绝更新就能及早发现阻尼因子或步长问题。fprintf中%s对应的info.stopReason可能是梯度收敛或参数变化量收敛如果出现lambda上限说明算法没能收敛需要调整初始条件。这里我用stem画残差图可以直观看到残差是否随机地在0附近波动。如果残差呈现明显的正弦趋势说明模型选择有误比如漏掉了频率项需要修改angelLove.m。上述代码运行后理论上会得到与true_theta比较接近的参数值。由于加入了随机噪声每次运行的拟合结果会有小幅波动但残差平方和应该在一个较小的范围。如果你得到的结果严重偏离真实值比如 $b$ 变成了 0.5最可能的原因是初始参数落在了一个局部极小值区域第5章会给出缓解策略。5. 进阶初值选择、全局收敛与过拟合防护5.1 初始猜测的启发式策略LM算法对初始参数的敏感度比高斯牛顿低但绝不是无脑随意给初始值。我总结的三条经验是一是利用模型本身的物理约束确定参数范围。比如角度拟合中振幅 $a$ 一定不小于数据的峰值偏移 $e$ 在数据均值附近阻尼系数 $d$ 可以通过相邻峰振幅比 $\ln(A_k/A_{k1})$ 除以时间间隔来粗估。二是对周期项使用FFT或简单的过零检测估计频率对初相位利用第一个过零点处的时间反推。三是如果完全没有先验信息可以采用多起点策略在某参数空间内随机生成若干初始点每个点用LM迭代固定次数选择目标函数最小的结果作为最终答案。下面是一个多起点启动的示例代码best_theta []; best_S inf; num_starts 20; for i 1:num_starts % 在合理范围内随机生成初始参数 theta0_try [545*rand(); 0.54.5*rand(); -pi2*pi*rand(); ... 0.010.9*rand(); -24*rand()]; [theta_try, S_try] lm(resfun, theta0_try, t, y, opt); if S_try best_S best_S S_try; best_theta theta_try; end end注意随机多起点虽然能提升找到全局最优的概率但并没有数学保证。对于高维参数问题更稳妥的做法是与遗传算法等全局优化器结合先用遗传算法粗搜再用LM细化。5.2 阻尼因子的上下界限制与自适应调整lm.m中已经实现了阻尼因子上下限保护但在实际工程项目里我还会再加一组专门针对病态问题的调整规则当lambda增大到一定程度后步长会变得极小这时目标函数几乎不变但梯度仍然较大。此时应检查是否因为模型存在参数冗余比如两个参数高度相关导致 $J^TJ$ 近似奇异。解决方法是给lambda加一个上限并在超过上限时提示用户检查模型。对于序列数据可以在迭代后期减小lambda的缩放系数从lambda/3改为lambda/2以更平滑地过渡到高斯牛顿阶段避免在最优解附近来回跳动。如果发现每次迭代大量时间花在求解正规方程上可以在迭代几次后保持lambda不变只更新J和g这种“断点更新”能显著提升大规模问题的效率但需要注意收敛速度会略微下降。另外当数据存在明显的异方差性不同时间点噪声方差不同时LM算法默认的最小二乘目标不再是最优准则。此时可以对残差加权即在angelLove.m中返回残差除以标准差的权重形式$$r_i \frac{y_i - f(x_i;\theta)}{\sigma_i}$$这样LM算法直接拟合加权残差等价于最大似然估计通常能获得更好的参数精度。5.3 防止过拟合参数冗余与正则化非线性拟合同样存在过拟合问题尤其是当你的模型参数比数据中的有效结构多时。我在angleLoveForYou这个例子里增加了偏移常数 $e$但如果传感器没有系统零漂硬加这个参数就会导致过拟合。一个典型信号是拟合得到的 $e$ 值很大而残差方差并没有显著下降。这时应该用赤池信息准则AIC或贝叶斯信息准则BIC比较不同参数数量的模型$$ \mathrm{AIC} m \ln(S/m) 2n $$其中 $m$ 是样本数$n$ 是参数个数$S$ 是残差平方和。我一般会写一个循环依次去掉非核心参数重新拟合计算AIC选择AIC最小的模型。在LM框架内也可以引入正则化把目标函数改为$$ S_{\rho}(\theta) \sum_i r_i^2 \rho \sum_j (\theta_j - \theta_j^{ref})^2 $$这相当于在正规方程中额外加上一个对角项与阻尼因子的数学结构完全一致区别在于正则化系数 $\rho$ 是固定的。实现时只需要修改残差函数向量在其末尾追加正则化惩罚项。不过要谨慎使用正则化会引入偏差如果样本数足够一般用AIC选模型比正则化更直接。另一个容易混淆的概念是LM阻尼与正则化的区别阻尼因子 $\lambda$ 只在求解步长时出现不改变目标函数本身而正则化项改变了目标函数拟合出的参数实际上是带惩罚后的解。理解这一点能帮助你正确解释拟合结果。关于最后一个技巧我想推荐一个每次拟合后必做的验证步骤把拟合残差绘制成自相关图并做Ljung-Box检验。如果残差存在显著自相关说明模型没有完全提取出数据中的时间结构此时盲目相信参数估计是有风险的。具体做法是调用crosscorr(resid, resid)查看滞后期自相关是否在置信带内这个步骤虽然简单但可以避免把系统误差当成随机噪声处理。本文还有配套的精品资源点击获取