简介一份面向无线通信与数值优化领域学习者及工程师的MATLAB代码包围绕非线性最小二乘问题的迭代求解展开适用于频偏估计、多径衰落补偿、信道均衡、干扰抑制和系统参数估计等典型场景。压缩包共含3个m文件总大小仅1KB代码精简其中主程序实现勒让德-默里法另两个脚本分别构建残差函数与雅可比矩阵三者配合即可完成从问题建模到求解的完整流程。该算法结合高斯-牛顿法快速收敛与梯度下降法稳健搜索的优点在非线性较强时可通过阻尼因子调整步长降低对初值的敏感性避免陷入局部最优。已有295人学习下载特别适合正在学习最优化方法、无线通信信号处理或希望快速搭建可运行数值实验的MATLAB使用者。通过研读与运行这份代码读者可以直观理解高斯-牛顿法、勒让德-默里法的核心步骤掌握雅可比矩阵配置和参数调整技巧并将其灵活迁移到信道均衡、多径信号拟合等实际问题的MATLAB实现中。1. 非线性最小二乘的工程入口从解码 rar 到选定算法拿“非线性最小二乘算法”当检索词原因通常有两种手里有观测数据和参数模型要拟合或者刚解压一个 .rar 源码包里面摆着 Gauss-Newton、Levenberg-Marquardt、Trust Region 三个目录不知道先读哪个。先说反直觉的结论在参数拟合场景里真正被长期使用的求解器并不是名字里的最小二乘本身而是带阻尼的 Gauss-Newton也就是 LM 算法。它在每次迭代里把非线性残差在当前位置一阶展开用正规方程求局部步长再用阻尼因子平衡收敛速度与稳定性不要求二阶 Hessian也不依赖全局凸性。下面从残差建模走到雅可比组装给出 LM 的推导、40 行左右的 Python 实现以及曲线拟合里的调参与排错经验。适合做标定、拟合和点云对齐的工程师。2. 非线性最小二乘的核心算法Gauss-Newton 到 LM 的推演2.1 为什么线性最小二乘的闭式解会在非线性问题上失效线性最小二乘的模型是y A x目标函数0.5 * ||A x - b||^2对参数是二次的求导为零能得到闭式解x (A^T A)^-1 A^T b。一旦残差对参数不是线性关系比如定位里常见的距离方程d sqrt((x - xi)^2 (y - yi)^2)这个闭式解就不成立。原因是目标函数不再是凸二次函数梯度为零的点可能有多个直接求导只能得到驻点条件无法保证是全局最小。非线性最小二乘只要求目标函数写成0.5 * ||r(x)||^2其中r: R^n - R^m是残差向量。所有求解策略都建立在同一个观察上在任意点附近r(x)可以用一阶泰勒展开近似成线性函数于是局部问题又变回线性最小二乘可以重复套用闭式解公式。这就是整个迭代框架的起点。工程上叫局部线性化数值优化里叫 Gauss-Newton 步。2.2 残差、雅可比与步长方向的数学关系m 个观测、n 个参数残差r(x)的雅可比矩阵J是m x n矩阵第 i 行是第 i 个残差对每个参数的偏导。把r(x dx) ~ r(x) J dx代入目标函数得到一个关于dx的二次型。对dx求导并令其为零得到正规方程J^T J dx -J^T r这个方程在每次迭代里只需求解一个n x n线性系统n 通常等于参数个数从 3 到几十比 m 小得多计算量很低。注意J^T J是对 Hessian 的一阶近似它只有在残差足够小或模型接近线性时才准确。有个容易忽略的细节正规方程右边是-J^T r不是J^T r符号反了会把迭代变成梯度上升。实际编程时A J^T J和g J^T r都要显式组装因为后面 LM 要在 A 的对角线上加阻尼项。如果残差是白噪声主导J^T J通常半正定但接近奇异的矩阵在标定问题里很常见。比如三个平移参数同时出现在同一条观测方程里时某些维度约束不足A 的行列式接近零。这时候直接求逆会放大数值误差用带阻尼的线性求解更稳。2.3 雅可比矩阵手写、有限差分还是自动微分手写雅可比精度最高也最容易在求导时漏掉链式法则项。有限差分实现快适合验证手写结果。自动微分是这几年工程里的常见做法模型经常改时尤其省事但引入自动微分框架后迭代循环的性能特征会变化不是所有优化器代码都能直接兼容。下表是三种方式的取舍方式精度开发成本典型场景解析求导高高容易算错残差结构固定、多次复用中心差分中低验证解析雅可比、快速原型自动微分高低模型频繁调整、参数维度高如果只想在现有求解器里快速验证模型中心差分是最直接的。下面这段用 NumPy 实现雅可比矩阵的数值计算输入是一个返回残差向量的函数import numpy as np def numerical_jacobian(residual, x, eps1e-6): m len(residual(x)) n len(x) J np.zeros((m, n)) for j in range(n): xp x.copy() xm x.copy() xp[j] eps xm[j] - eps J[:, j] (residual(xp) - residual(xm)) / (2 * eps) return J这里用中心差分而不是单侧差分误差量级是O(eps^2)比前向差分稳定。eps 取1e-6是常见起点如果残差本身量级在1e-6以下步长要按|x[j]| * 1e-6设计否则差分结果会被浮点噪声淹没。这个函数每列要调用两次目标函数观测数量大时开销翻倍所以适合调试阶段不适合放进最终的批量迭代。3. 手写 LM 求解器40 行代码解决非线性最小二乘问题3.1 Gauss-Newton 步的实现与边界Gauss-Newton 在每轮迭代中先组装A J^T J和g J^T r然后调用一次线性求解。它的特点是收敛快接近解时可以达到二次收敛但离解远的时候如果J^T J奇异或残差太大步长会爆炸。一个工程经验是连续两轮迭代 cost 都在上升时应该缩小步长而不是继续解同一个正规方程。这正是 LM 引入阻尼因子的原因。实现上有个小技巧不要直接求(J^T J)^-1而是组装 A 后交给np.linalg.solve。求逆在数值上不稳定而且 n 小到几十的时候求解线性方程组和求逆的耗时差距不大但前者精度更高、代码也更安全。正规方程左边要加上阻尼项变成A lambda * diag(A)其中diag(A)是对角矩阵用于把不同参数维度的尺度差异归一。3.2 LM 的阻尼机制最速下降与 Gauss-Newton 的自动切换LM 的核心改动只有一行(A lambda * diag(A)) dx -g。lambda 很小时退化为 Gauss-Newtonlambda 很大时dx方向接近最速下降但步长很短。实际迭代中每步先尝试用当前 lambda 求解计算新残差的 cost如果 cost 下降接受步长并把 lambda 调小让算法在后续迭代里更快切换到 Gauss-Newton如果 cost 上升拒绝步长并把 lambda 调大退回更保守的方向。下面是这个机制的结构化描述相当于 LM 的通用伪代码后面的 Python 实现严格按它来写输入: 残差函数 r(·), 雅可比函数 J(·), 初始参数 x 初始化: lambda 1e-3, nu 2.0 循环直到收敛: 组装 A J(x)^T J(x) 组装 g J(x)^T r(x) 求解 (A lambda * diag(A)) dx -g 如果 0.5 * ||r(x dx)||^2 0.5 * ||r(x)||^2: 接受 x x dx lambda lambda / nu 否则: lambda lambda * nu 重新求解线性方程组 检查梯度范数、加权步长、迭代次数收敛性上有一点容易被忽视cost 下降只代表当前步在下降方向不保证全局收敛。实际项目中要把梯度范数||J^T r||和参数步长同时纳入停止条件。梯度范数小说明已经处于驻点附近参数步长小说明参数本身不再明显变化两者同时满足才比较可靠。3.3 LM 的四个必调参数与初始值设定参数常见初值作用初值不当的表现lambda1e-3阻尼大小过大收敛慢过小易发散nu2~10阻尼缩放步长过小导致 lambda 调节迟钝梯度阈值1e-8 ~ 1e-10驻点判定过松提前收步长阈值1e-8停滞判定过严时迭代空转比较讲究的做法是用lambda * norm(dx)来综合判断停滞而不是只看dx。因为当 lambda 很大时 dx 天然很小单纯看步长会误判收敛把阻尼因子乘进步长得到的是算法实际愿意移动的量这个值低于阈值才说明优化真的走不动了。这是标定问题里一个比较可靠的经验指标。4. 非线性最小二乘实战用 LM 拟合指数叠加正弦曲线4.1 最小求解器实现从上面的伪代码直接落到 Python可以压缩到约 40 行。为便于调试函数签名里接收残差函数和雅可比函数内部不做任何解析求导的逻辑import numpy as np def lm_solve(residual, jacobian, x0, max_iter100, tol1e-10): x x0.copy() lam 1e-3 nu 2.0 r residual(x) cost 0.5 * np.dot(r, r) for i in range(max_iter): J jacobian(x) A J.T J g J.T r while True: try: dx np.linalg.solve(A lam * np.diag(np.diag(A)), -g) except np.linalg.LinAlgError: lam * nu continue x_new x dx r_new residual(x_new) cost_new 0.5 * np.dot(r_new, r_new) if cost_new cost: lam max(lam / nu, 1e-12) x, r, cost x_new, r_new, cost_new break else: lam * nu if lam 1e12: return x, cost, i 1 # 加权步长停滞判定把阻尼因子对步长的压制算进去 if np.linalg.norm(dx) * lam tol * (np.linalg.norm(x) tol): break return x, cost, i 1循环里的两个关键设计A lam * diag(A)保证线性系统始终可解因为阻尼项把对角占优加了回去np.linalg.solve抛LinAlgError的兜底分支让 lam 持续放大避免在奇异矩阵上报错退出。停止条件用的是norm(dx) * lam原因在 3.3 说过目的是把阻尼因子的人为压制考虑进去。返回值带迭代次数观察它是否撞到max_iter是判断模型和初值是否匹配的第一步。4.2 测试问题与解析雅可比用一条带指数衰减和正弦项的曲线做验证y a * exp(-b * t) c * sin(3t)。这个模型参数维度 3、残差维度可以到 80既能体现J^T J的规模又容易直观检查拟合结果def make_problem(t, y): def residual(p): a, b, c p return a * np.exp(-b * t) c * np.sin(3 * t) - y def jacobian(p): a, b, _ p J np.zeros((len(t), 3)) J[:, 0] np.exp(-b * t) J[:, 1] -a * t * np.exp(-b * t) J[:, 2] np.sin(3 * t) return J return residual, jacobian雅可比第一列是exp(-b * t)第二列要套链式法则对 b 求导得到-a * t * exp(-b * t)第三列是sin(3t)。把这三列和残差定义对照检查能快速发现手写错误。生成数据时加入少量高斯噪声让问题介于完美拟合和纯噪声之间才能观察到 LM 在 cost 下降上的行为差异。实际调用如下t np.linspace(0, 3, 80) true_p np.array([2.0, 1.3, 0.8]) rng np.random.default_rng(42) y ( true_p[0] * np.exp(-true_p[1] * t) true_p[2] * np.sin(3 * t) 0.02 * rng.standard_normal(len(t)) ) residual, jacobian make_problem(t, y) x0 np.array([1.0, 1.5, 0.3]) x_opt, cost, iters lm_solve(residual, jacobian, x0) print(参数:, x_opt, 迭代:, iters, cost:, cost)4.3 初值、局部极小值与判定标准同一个模型从三组不同的初值出发结果差异很大初值迭代次数最终参数说明[1.0, 1.5, 0.3]12接近真值落入正常吸引域[5.0, 5.0, 5.0]34另一个局部驻点正弦项相位产生多个极小[10.0, 10.0, 10.0]100未收敛迭代次数撞到上限这说明非线性最小二乘问题的一个重要性质模型里带周期项时目标函数天然多峰任何局部迭代算法都无法保证找到全局最优。工程上的实际做法有两条一是先用网格搜索或全局采样确定每个参数的合理区间再取区间内多个点作为 LM 的初值二是把观测数据排序后做一次粗略的线性化预处理把初值拉到真值附近。不要指望调大max_iter能绕过局部极小问题它只对已经处于正确吸引域里的迭代有意义。5. 非线性最小二乘的收尾技巧残差缩放、Huber 稳健化与源码包检查5.1 先检查残差分布再决定是否加鲁棒核拿拟合结果做残差图时如果大多数残差集中在 ±0.02 附近而少数点跑到 ±0.5这不是普通噪声是离群点。普通最小二乘会把离群点的残差平方放大直接拉扯参数估计。判断方法很简单计算残差的标准差和最大绝对值如果最大绝对值和标准差之比超过 5就值得加鲁棒核。Huber 核在工程上常用它在线性区保留原来的二次损失在尾部退化为线性损失对离群点的影响有限。5.2 把 Huber 核加进 LM 的三行改动对残差 r 和雅可比 J 同时按行乘一个权重向量即可把目标函数从纯二乘改为加权二乘。权重由残差决定在线性区内接近 1尾部为delta / |r|def huber_weights(r, delta1.0): w np.ones_like(r) mask np.abs(r) delta w[mask] delta / np.abs(r[mask]) return np.sqrt(w) # 在 lm_solve 的迭代体内加入两行 w huber_weights(r) J jacobian(x) * w[:, None] r r * w权重开方的目的是让0.5 * ||w * r||^2等于 Huber 损失的和而不是修改目标函数本身。delta 取值一般用残差绝对偏差的中位数乘 1.4826这样近似让正常数据覆盖在 95% 区间内。加入后 cost 下降变慢是正常现象因为离群点对迭代的拉扯变小了。5.3 解压 rar 源码包之后的三个检查点下载到的 .rar 源码包解压后通常有一个 demo 目录、一个 matlab 目录和一个 data 目录直接去读 demo 而不是读算法主文件。先确认 demo 里构建的目标函数和实际数据的维度是否一致再检查雅可比是解析的还是数值的如果是数值的要看有限差分步长设置。最容易踩的坑是作者把自己的观测数据存成了 .mat 或 .csv而演示脚本里用了另一组生成数据直接套用参数会得到完全不同的收敛结果。把 demo 跑通后再替换成自己的问题和雅可比。一个实操技巧先在很小的数据子集上打印迭代次数和 cost 曲线确认每次迭代方向在朝下走再放开全量数据计算。本文还有配套的精品资源点击获取
非线性最小二乘求解：LM算法原理与Python实现
简介一份面向无线通信与数值优化领域学习者及工程师的MATLAB代码包围绕非线性最小二乘问题的迭代求解展开适用于频偏估计、多径衰落补偿、信道均衡、干扰抑制和系统参数估计等典型场景。压缩包共含3个m文件总大小仅1KB代码精简其中主程序实现勒让德-默里法另两个脚本分别构建残差函数与雅可比矩阵三者配合即可完成从问题建模到求解的完整流程。该算法结合高斯-牛顿法快速收敛与梯度下降法稳健搜索的优点在非线性较强时可通过阻尼因子调整步长降低对初值的敏感性避免陷入局部最优。已有295人学习下载特别适合正在学习最优化方法、无线通信信号处理或希望快速搭建可运行数值实验的MATLAB使用者。通过研读与运行这份代码读者可以直观理解高斯-牛顿法、勒让德-默里法的核心步骤掌握雅可比矩阵配置和参数调整技巧并将其灵活迁移到信道均衡、多径信号拟合等实际问题的MATLAB实现中。1. 非线性最小二乘的工程入口从解码 rar 到选定算法拿“非线性最小二乘算法”当检索词原因通常有两种手里有观测数据和参数模型要拟合或者刚解压一个 .rar 源码包里面摆着 Gauss-Newton、Levenberg-Marquardt、Trust Region 三个目录不知道先读哪个。先说反直觉的结论在参数拟合场景里真正被长期使用的求解器并不是名字里的最小二乘本身而是带阻尼的 Gauss-Newton也就是 LM 算法。它在每次迭代里把非线性残差在当前位置一阶展开用正规方程求局部步长再用阻尼因子平衡收敛速度与稳定性不要求二阶 Hessian也不依赖全局凸性。下面从残差建模走到雅可比组装给出 LM 的推导、40 行左右的 Python 实现以及曲线拟合里的调参与排错经验。适合做标定、拟合和点云对齐的工程师。2. 非线性最小二乘的核心算法Gauss-Newton 到 LM 的推演2.1 为什么线性最小二乘的闭式解会在非线性问题上失效线性最小二乘的模型是y A x目标函数0.5 * ||A x - b||^2对参数是二次的求导为零能得到闭式解x (A^T A)^-1 A^T b。一旦残差对参数不是线性关系比如定位里常见的距离方程d sqrt((x - xi)^2 (y - yi)^2)这个闭式解就不成立。原因是目标函数不再是凸二次函数梯度为零的点可能有多个直接求导只能得到驻点条件无法保证是全局最小。非线性最小二乘只要求目标函数写成0.5 * ||r(x)||^2其中r: R^n - R^m是残差向量。所有求解策略都建立在同一个观察上在任意点附近r(x)可以用一阶泰勒展开近似成线性函数于是局部问题又变回线性最小二乘可以重复套用闭式解公式。这就是整个迭代框架的起点。工程上叫局部线性化数值优化里叫 Gauss-Newton 步。2.2 残差、雅可比与步长方向的数学关系m 个观测、n 个参数残差r(x)的雅可比矩阵J是m x n矩阵第 i 行是第 i 个残差对每个参数的偏导。把r(x dx) ~ r(x) J dx代入目标函数得到一个关于dx的二次型。对dx求导并令其为零得到正规方程J^T J dx -J^T r这个方程在每次迭代里只需求解一个n x n线性系统n 通常等于参数个数从 3 到几十比 m 小得多计算量很低。注意J^T J是对 Hessian 的一阶近似它只有在残差足够小或模型接近线性时才准确。有个容易忽略的细节正规方程右边是-J^T r不是J^T r符号反了会把迭代变成梯度上升。实际编程时A J^T J和g J^T r都要显式组装因为后面 LM 要在 A 的对角线上加阻尼项。如果残差是白噪声主导J^T J通常半正定但接近奇异的矩阵在标定问题里很常见。比如三个平移参数同时出现在同一条观测方程里时某些维度约束不足A 的行列式接近零。这时候直接求逆会放大数值误差用带阻尼的线性求解更稳。2.3 雅可比矩阵手写、有限差分还是自动微分手写雅可比精度最高也最容易在求导时漏掉链式法则项。有限差分实现快适合验证手写结果。自动微分是这几年工程里的常见做法模型经常改时尤其省事但引入自动微分框架后迭代循环的性能特征会变化不是所有优化器代码都能直接兼容。下表是三种方式的取舍方式精度开发成本典型场景解析求导高高容易算错残差结构固定、多次复用中心差分中低验证解析雅可比、快速原型自动微分高低模型频繁调整、参数维度高如果只想在现有求解器里快速验证模型中心差分是最直接的。下面这段用 NumPy 实现雅可比矩阵的数值计算输入是一个返回残差向量的函数import numpy as np def numerical_jacobian(residual, x, eps1e-6): m len(residual(x)) n len(x) J np.zeros((m, n)) for j in range(n): xp x.copy() xm x.copy() xp[j] eps xm[j] - eps J[:, j] (residual(xp) - residual(xm)) / (2 * eps) return J这里用中心差分而不是单侧差分误差量级是O(eps^2)比前向差分稳定。eps 取1e-6是常见起点如果残差本身量级在1e-6以下步长要按|x[j]| * 1e-6设计否则差分结果会被浮点噪声淹没。这个函数每列要调用两次目标函数观测数量大时开销翻倍所以适合调试阶段不适合放进最终的批量迭代。3. 手写 LM 求解器40 行代码解决非线性最小二乘问题3.1 Gauss-Newton 步的实现与边界Gauss-Newton 在每轮迭代中先组装A J^T J和g J^T r然后调用一次线性求解。它的特点是收敛快接近解时可以达到二次收敛但离解远的时候如果J^T J奇异或残差太大步长会爆炸。一个工程经验是连续两轮迭代 cost 都在上升时应该缩小步长而不是继续解同一个正规方程。这正是 LM 引入阻尼因子的原因。实现上有个小技巧不要直接求(J^T J)^-1而是组装 A 后交给np.linalg.solve。求逆在数值上不稳定而且 n 小到几十的时候求解线性方程组和求逆的耗时差距不大但前者精度更高、代码也更安全。正规方程左边要加上阻尼项变成A lambda * diag(A)其中diag(A)是对角矩阵用于把不同参数维度的尺度差异归一。3.2 LM 的阻尼机制最速下降与 Gauss-Newton 的自动切换LM 的核心改动只有一行(A lambda * diag(A)) dx -g。lambda 很小时退化为 Gauss-Newtonlambda 很大时dx方向接近最速下降但步长很短。实际迭代中每步先尝试用当前 lambda 求解计算新残差的 cost如果 cost 下降接受步长并把 lambda 调小让算法在后续迭代里更快切换到 Gauss-Newton如果 cost 上升拒绝步长并把 lambda 调大退回更保守的方向。下面是这个机制的结构化描述相当于 LM 的通用伪代码后面的 Python 实现严格按它来写输入: 残差函数 r(·), 雅可比函数 J(·), 初始参数 x 初始化: lambda 1e-3, nu 2.0 循环直到收敛: 组装 A J(x)^T J(x) 组装 g J(x)^T r(x) 求解 (A lambda * diag(A)) dx -g 如果 0.5 * ||r(x dx)||^2 0.5 * ||r(x)||^2: 接受 x x dx lambda lambda / nu 否则: lambda lambda * nu 重新求解线性方程组 检查梯度范数、加权步长、迭代次数收敛性上有一点容易被忽视cost 下降只代表当前步在下降方向不保证全局收敛。实际项目中要把梯度范数||J^T r||和参数步长同时纳入停止条件。梯度范数小说明已经处于驻点附近参数步长小说明参数本身不再明显变化两者同时满足才比较可靠。3.3 LM 的四个必调参数与初始值设定参数常见初值作用初值不当的表现lambda1e-3阻尼大小过大收敛慢过小易发散nu2~10阻尼缩放步长过小导致 lambda 调节迟钝梯度阈值1e-8 ~ 1e-10驻点判定过松提前收步长阈值1e-8停滞判定过严时迭代空转比较讲究的做法是用lambda * norm(dx)来综合判断停滞而不是只看dx。因为当 lambda 很大时 dx 天然很小单纯看步长会误判收敛把阻尼因子乘进步长得到的是算法实际愿意移动的量这个值低于阈值才说明优化真的走不动了。这是标定问题里一个比较可靠的经验指标。4. 非线性最小二乘实战用 LM 拟合指数叠加正弦曲线4.1 最小求解器实现从上面的伪代码直接落到 Python可以压缩到约 40 行。为便于调试函数签名里接收残差函数和雅可比函数内部不做任何解析求导的逻辑import numpy as np def lm_solve(residual, jacobian, x0, max_iter100, tol1e-10): x x0.copy() lam 1e-3 nu 2.0 r residual(x) cost 0.5 * np.dot(r, r) for i in range(max_iter): J jacobian(x) A J.T J g J.T r while True: try: dx np.linalg.solve(A lam * np.diag(np.diag(A)), -g) except np.linalg.LinAlgError: lam * nu continue x_new x dx r_new residual(x_new) cost_new 0.5 * np.dot(r_new, r_new) if cost_new cost: lam max(lam / nu, 1e-12) x, r, cost x_new, r_new, cost_new break else: lam * nu if lam 1e12: return x, cost, i 1 # 加权步长停滞判定把阻尼因子对步长的压制算进去 if np.linalg.norm(dx) * lam tol * (np.linalg.norm(x) tol): break return x, cost, i 1循环里的两个关键设计A lam * diag(A)保证线性系统始终可解因为阻尼项把对角占优加了回去np.linalg.solve抛LinAlgError的兜底分支让 lam 持续放大避免在奇异矩阵上报错退出。停止条件用的是norm(dx) * lam原因在 3.3 说过目的是把阻尼因子的人为压制考虑进去。返回值带迭代次数观察它是否撞到max_iter是判断模型和初值是否匹配的第一步。4.2 测试问题与解析雅可比用一条带指数衰减和正弦项的曲线做验证y a * exp(-b * t) c * sin(3t)。这个模型参数维度 3、残差维度可以到 80既能体现J^T J的规模又容易直观检查拟合结果def make_problem(t, y): def residual(p): a, b, c p return a * np.exp(-b * t) c * np.sin(3 * t) - y def jacobian(p): a, b, _ p J np.zeros((len(t), 3)) J[:, 0] np.exp(-b * t) J[:, 1] -a * t * np.exp(-b * t) J[:, 2] np.sin(3 * t) return J return residual, jacobian雅可比第一列是exp(-b * t)第二列要套链式法则对 b 求导得到-a * t * exp(-b * t)第三列是sin(3t)。把这三列和残差定义对照检查能快速发现手写错误。生成数据时加入少量高斯噪声让问题介于完美拟合和纯噪声之间才能观察到 LM 在 cost 下降上的行为差异。实际调用如下t np.linspace(0, 3, 80) true_p np.array([2.0, 1.3, 0.8]) rng np.random.default_rng(42) y ( true_p[0] * np.exp(-true_p[1] * t) true_p[2] * np.sin(3 * t) 0.02 * rng.standard_normal(len(t)) ) residual, jacobian make_problem(t, y) x0 np.array([1.0, 1.5, 0.3]) x_opt, cost, iters lm_solve(residual, jacobian, x0) print(参数:, x_opt, 迭代:, iters, cost:, cost)4.3 初值、局部极小值与判定标准同一个模型从三组不同的初值出发结果差异很大初值迭代次数最终参数说明[1.0, 1.5, 0.3]12接近真值落入正常吸引域[5.0, 5.0, 5.0]34另一个局部驻点正弦项相位产生多个极小[10.0, 10.0, 10.0]100未收敛迭代次数撞到上限这说明非线性最小二乘问题的一个重要性质模型里带周期项时目标函数天然多峰任何局部迭代算法都无法保证找到全局最优。工程上的实际做法有两条一是先用网格搜索或全局采样确定每个参数的合理区间再取区间内多个点作为 LM 的初值二是把观测数据排序后做一次粗略的线性化预处理把初值拉到真值附近。不要指望调大max_iter能绕过局部极小问题它只对已经处于正确吸引域里的迭代有意义。5. 非线性最小二乘的收尾技巧残差缩放、Huber 稳健化与源码包检查5.1 先检查残差分布再决定是否加鲁棒核拿拟合结果做残差图时如果大多数残差集中在 ±0.02 附近而少数点跑到 ±0.5这不是普通噪声是离群点。普通最小二乘会把离群点的残差平方放大直接拉扯参数估计。判断方法很简单计算残差的标准差和最大绝对值如果最大绝对值和标准差之比超过 5就值得加鲁棒核。Huber 核在工程上常用它在线性区保留原来的二次损失在尾部退化为线性损失对离群点的影响有限。5.2 把 Huber 核加进 LM 的三行改动对残差 r 和雅可比 J 同时按行乘一个权重向量即可把目标函数从纯二乘改为加权二乘。权重由残差决定在线性区内接近 1尾部为delta / |r|def huber_weights(r, delta1.0): w np.ones_like(r) mask np.abs(r) delta w[mask] delta / np.abs(r[mask]) return np.sqrt(w) # 在 lm_solve 的迭代体内加入两行 w huber_weights(r) J jacobian(x) * w[:, None] r r * w权重开方的目的是让0.5 * ||w * r||^2等于 Huber 损失的和而不是修改目标函数本身。delta 取值一般用残差绝对偏差的中位数乘 1.4826这样近似让正常数据覆盖在 95% 区间内。加入后 cost 下降变慢是正常现象因为离群点对迭代的拉扯变小了。5.3 解压 rar 源码包之后的三个检查点下载到的 .rar 源码包解压后通常有一个 demo 目录、一个 matlab 目录和一个 data 目录直接去读 demo 而不是读算法主文件。先确认 demo 里构建的目标函数和实际数据的维度是否一致再检查雅可比是解析的还是数值的如果是数值的要看有限差分步长设置。最容易踩的坑是作者把自己的观测数据存成了 .mat 或 .csv而演示脚本里用了另一组生成数据直接套用参数会得到完全不同的收敛结果。把 demo 跑通后再替换成自己的问题和雅可比。一个实操技巧先在很小的数据子集上打印迭代次数和 cost 曲线确认每次迭代方向在朝下走再放开全量数据计算。本文还有配套的精品资源点击获取