Python实战：利用scipy.optimize.least_squares实现鲁棒最小二乘拟合

1. 为什么需要鲁棒最小二乘拟合在实际工程和科研中我们经常会遇到这样的场景采集到的实验数据里混入了异常值或者某些测量点明显偏离正常范围。这时候如果直接用传统的最小二乘法进行拟合结果往往会受到这些坏点的严重影响。我去年在做传感器标定时就遇到过这种情况——明明标定板摆放得很正但总有那么几个点的测量值明显异常。传统最小二乘法的目标是最小化残差的平方和这意味着它对大残差异常值非常敏感。举个例子假设我们有10个数据点其中9个都很准但有1个偏离了10倍。这个异常点对最终拟合结果的影响会是正常点的100倍因为平方关系。这显然不是我们想要的结果。鲁棒最小二乘法Robust Least Squares就是为了解决这个问题而生的。它的核心思想是通过修改损失函数降低大残差对整体优化的影响。在Python的scipy.optimize模块中least_squares函数提供了多种鲁棒损失函数选项让我们能够轻松实现这种抗干扰拟合。2. least_squares函数核心参数解析2.1 基本调用方式scipy.optimize.least_squares的基本调用形式如下from scipy.optimize import least_squares result least_squares( fun, # 残差计算函数 x0, # 初始参数猜测 methodtrf, # 优化算法 losslinear, # 损失函数类型 args(), # 传递给fun的额外参数 kwargs{} # 传递给fun的额外关键字参数 )其中最关键的是loss参数它决定了如何处理残差。默认的linear就是普通最小二乘而要实现鲁棒拟合我们需要改用其他选项。2.2 鲁棒损失函数选项least_squares提供了几种常用的鲁棒损失函数huber对中等大小的残差使用平方损失对大残差使用线性损失soft_l1平滑版的L1损失对大残差的抑制比huber更强cauchy使用对数损失对异常值非常不敏感arctan使用反正切函数对大残差的抑制最温和我在处理工业传感器数据时发现对于中等程度的异常值huber通常就能取得不错的效果而当数据中有明显离群点时cauchy表现更好。下面这个表格总结了各损失函数的特性损失函数对大残差的敏感度适用场景linear非常高干净数据huber中等少量异常值soft_l1较低较多异常值cauchy非常低严重离群点arctan最低极端异常值3. 实战带噪声数据的曲线拟合3.1 生成模拟数据让我们通过一个具体例子来看看鲁棒拟合的效果。假设我们要拟合一个简单的线性模型y a*x b但数据中包含一些异常值import numpy as np import matplotlib.pyplot as plt from scipy.optimize import least_squares # 生成干净数据 np.random.seed(42) x np.linspace(0, 10, 50) y_true 2.5 * x 1.0 # 添加高斯噪声 y_noise y_true np.random.normal(scale1.0, sizelen(x)) # 人为添加异常值 outlier_indices [5, 15, 25, 35, 45] y_noise[outlier_indices] 15 * np.random.rand(len(outlier_indices))3.2 定义残差函数我们需要定义一个计算残差的函数这个函数将作为least_squares的第一个参数def residual(params, x, y): 计算模型预测值与实际值的残差 a, b params return y - (a * x b)3.3 比较不同损失函数的效果现在我们可以比较不同loss参数下的拟合效果# 初始参数猜测 x0 [1.0, 0.0] # 普通最小二乘 res_linear least_squares(residual, x0, losslinear, args(x, y_noise)) # Huber损失 res_huber least_squares(residual, x0, losshuber, args(x, y_noise)) # Cauchy损失 res_cauchy least_squares(residual, x0, losscauchy, args(x, y_noise)) # 绘制结果 plt.figure(figsize(10, 6)) plt.scatter(x, y_noise, label带噪声数据, alpha0.6) plt.plot(x, y_true, k--, label真实模型, linewidth2) plt.plot(x, res_linear.x[0]*x res_linear.x[1], labellinear) plt.plot(x, res_huber.x[0]*x res_huber.x[1], labelhuber) plt.plot(x, res_cauchy.x[0]*x res_cauchy.x[1], labelcauchy) plt.legend() plt.show()从图中可以明显看出使用默认linear损失的拟合线被异常值严重拉偏而huber和cauchy损失则能较好地抵抗异常值的影响其中cauchy的表现尤为出色。4. 进阶技巧与实战建议4.1 损失函数参数调优除了选择损失函数类型我们还可以通过f_scale参数进一步调整鲁棒性。这个参数控制着损失函数从二次行为转变为线性行为的阈值。例如对于huber损失# 调整huber损失的阈值 res_huber_tuned least_squares( residual, x0, losshuber, f_scale0.1, # 默认是1.0 args(x, y_noise) )较小的f_scale会使算法对异常值更敏感较大的值则会更鲁棒。在实际项目中我通常会尝试几个不同的f_scale值然后选择使最终残差分布最合理的那个。4.2 结合其他优化算法least_squares提供了多种优化算法通过method参数指定不同的算法在处理不同问题时效果可能不同trfTrust Region Reflective默认算法适用于边界约束问题lmLevenberg-Marquardt适用于无约束问题对小残差问题效率高dogbox适用于变量较多的问题我的经验是对于大多数鲁棒拟合问题默认的trf就已经足够好。但在处理非常大的数据集时dogbox可能会更快一些。4.3 处理实际问题时的注意事项在实际项目中应用鲁棒最小二乘时有几个容易踩的坑值得注意初始值选择虽然鲁棒方法对初始值不如传统方法敏感但好的初始值仍然能加快收敛。我通常会先用普通最小二乘快速拟合一次用其结果作为鲁棒拟合的初始值。参数缩放当不同参数的数值尺度差异很大时比如一个参数在0.1左右另一个在1000左右最好先对数据进行归一化处理否则可能会影响优化效果。结果验证拟合完成后一定要检查残差分布。好的拟合应该使残差随机分布没有明显的模式。我习惯绘制残差直方图和残差-预测值散点图来验证。异常值检测鲁棒拟合完成后那些残差特别大的点很可能就是异常值。可以设置一个阈值比如3倍标准差之外的点来自动标记这些异常值供后续分析使用。5. 真实案例传感器标定应用去年我在做一个工业机器人的力传感器标定项目时就深刻体会到了鲁棒最小二乘的价值。我们需要通过一组已知力和对应的传感器读数来建立力-电压关系的数学模型。由于现场环境复杂采集的数据中总会有一些异常点——可能是由于瞬时电磁干扰、机械振动或其他未知因素造成的。最初使用普通最小二乘时标定结果总是不稳定机器人表现时好时坏。后来改用huber损失的鲁棒拟合后问题迎刃而解。具体实现代码如下def sensor_residual(params, known_forces, sensor_readings): 传感器标定残差计算 scale, offset params return sensor_readings - (scale * known_forces offset) # 初始猜测根据传感器规格书 x0 [0.1, 0.0] # 鲁棒拟合 result least_squares( sensor_residual, x0, losshuber, args(calibration_forces, calibration_readings), methodtrf ) # 提取标定参数 scale, offset result.x这个案例中鲁棒拟合不仅提高了标定精度还使得我们可以自动识别出那些异常的标定点残差特别大的点帮助我们发现了一些传感器安装问题。