ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

多元非线性目标函数求解:从最小二乘结构到scipy工程实践

多元非线性目标函数求解:从最小二乘结构到scipy工程实践 简介这份资源针对多元非线性目标函数求解这一数学优化核心问题提供了一套MATLAB/Simulink实现样例适合正在学习最优化、需处理带约束规划问题的学生或工程师。压缩包共3个文件主要由.m函数脚本与.mdl仿真模型组成前者负责定义目标函数和约束条件后者用于搭建非线性系统的控制或仿真场景整个资源包仅13KB轻量精简便于快速下载、阅读与二次修改。已有439人学习使用可作为课程设计、科研实验或入门练习的起点。通过运行这些示例能够直观理解从约束建模、目标函数设置到调用优化算法求解的完整链路并为后续结合梯度下降法、拉格朗日乘数法、惩罚函数法或MATLAB优化工具解决更复杂问题提供参考框架。这套小资源特别适合用最小成本验证多元非线性优化方法的实际效果帮助把抽象理论落地到可运行的项目中。1. 多元非线性目标函数求解解决什么问题很多不沾边的需求——传感器标定、化学反应动力学参数拟合、机械臂参数辨识、SLAM后端图优化——落到数学上都是同一个问题多元非线性目标函数求解。共同点是目标函数对变量不是线性梯度很难看初值给得不好就掉进局部极小。能写出目标函数的工程师不少能从结构出发选对求解器的人不多。下面按一线做参数拟合的习惯把目标函数结构、求解器选型、scipy 实操、带约束与鲁棒损失、结果验证串成一条完整路径。适合有 Python 基础、正在做标定或拟合的工程师。2. 从多元非线性目标函数结构到求解器选型先做数学判断再做实现2.1 目标函数的标准形式先判断是不是最小二乘结构多元非线性目标函数求解时我第一件事不是选算法而是把问题写成标准形式对照一下。一般目标函数 minimize f(x) x ∈ R^n 非线性最小二乘 minimize F(x) 1/2 Σ r_i(x)^2如果目标函数本身就是一堆残差平方和比如拟合误差、标定重投影误差那优先交给scipy.optimize.least_squares。它会把残差向量直接展开用雅可比矩阵近似目标函数的海森矩阵而通用minimize只能把 F(x) 当作黑盒用拟牛顿法一点点试探曲率。同样的精度least_squares通常少跑一个数量级的迭代。反过来目标函数不是平方和形态——例如加了 log 鲁棒损失、L1 正则、或目标函数本身是仿真器输出——就该走minimize。这个判断决定了后面所有代码路径。认为“反正都是优化随便选”是多元非线性目标函数求解最常见的弯路。2.2 梯度、Hessian 与一阶最优性条件的作用最优解处梯度为零这是所有迭代法共同的收敛目标。差别在于下一步怎么走牛顿法要解 H(x) Δx -g(x)但显式构造 Hessian 对 n 稍大的问题代价太高。工程上最常用的是拟牛顿法BFGS 用每一步的梯度差去近似 Hessian 的逆只额外付出向量内积的开销。最小二乘结构更取巧F(x) 的 Hessian 是 J^T J Σ r_i ∇²r_i其中 J 是 m×n 雅可比矩阵。当残差接近零时J^T J 已经是对 Hessian 的可靠近似least_squares内部会大量利用这一点。这也是为什么同样的问题分别用least_squares和minimize跑前者的迭代路径会直很多。收敛时看两个判据梯度范数小于gtol或参数更新量小于xtol。调参时先放宽这两个值到 1e-8 左右别一上来就压到 1e-15否则大概率把时间耗在震荡上。2.3 求解器选型表与场景判定不同工具面向的问题结构差异很大。按我实际会用的范围选型逻辑如下工具入口适合场景主要限制scipy.optimize.least_squaresmethodtrf/lm残差形式的目标函数带边界lm不支持边界约束scipy.optimize.minimizemethodBFGS/L-BFGS-B通用目标函数、鲁棒损失、正则项无残差结构可用收敛慢scipy.optimize.differential_evolution全局优化器变量少、初值完全未知变量多时收敛极慢Ceres SolverC 接口大规模稀疏非线性最小二乘、视觉标定工程接入成本高要先组织残差块g2oC 图优化库多传感器标定、SLAM 后端基于图结构建模需要定义顶点和边判断顺序一般是能写成残差平方和先上least_squares目标函数里混入了鲁棒损失或自定义惩罚项改用minimize问题带稀疏雅可比和几十万个残差直接考虑 Ceres。注意differential_evolution不是常规武器只有你对初值完全没有概念、变量不超过十几个时才值得用。2.4 变量缩放与稀疏结构大规模场景的隐藏前提还有一个经常被忽略的预处理步骤变量缩放。如果参数里有量级 1e-6 的畸变系数也有量级 1000 的平移量trust-region 算法内部对步长的控制会非常难受轻则多跑几百次重则直接报边界不可行。我一般会在建模前把每个参数除以其量级把搜索空间归一化到相近尺度求解完成后再把结果乘回去。稀疏结构决定你是不是该离开 scipy。在机械臂运动学参数辨识这类问题里一个残差往往只依赖少数几个关节参数雅可比矩阵是高度稀疏的。least_squares默认用稠密线性代数n 到几千还能撑n 上十万就只能换 Ceres 这类稀疏求解器。判断依据很简单观察每个残差依赖的变量个数如果远小于变量总数它就是稀疏问题。3. 用 scipy.optimize 跑通多元非线性目标函数求解的最小流程3.1 构造带噪声的拟合样本用衰减正弦信号做例子。真实参数为 a、b、c、d生成合成观测整个流程替换成自己的模型和观测即可。import numpy as np from scipy.optimize import least_squares rng np.random.default_rng(42) x_data np.linspace(0, 4, 60) true_params np.array([2.5, 1.3, 3.0, 0.6]) def model(x, a, b, c, d): return a * np.exp(-b * x) * np.sin(c * x d) y_obs model(x_data, *true_params) 0.08 * rng.standard_normal(x_data.size)这段代码里model是观测模型噪声标准差 0.08 覆盖了信号峰值的约 3%。残差函数接下来直接引用model和y_obs所以替换成真实工程数据时只要保证residuals(theta)返回预测减观测的一维数组即可。3.2 用 least_squares 求最小二乘目标函数的最优参数残差函数和求解部分如下def residuals(theta): return model(x_data, *theta) - y_obs theta0 np.array([2.0, 1.0, 2.0, 0.0]) bounds ([0.1, 0.01, 0.1, -np.pi], [10.0, 5.0, 10.0, np.pi]) result least_squares( residuals, theta0, methodtrf, boundsbounds, max_nfev5000, ftol1e-12, xtol1e-12, gtol1e-12, ) print(result.x) print(result.cost)residuals返回的是每个样本的预测误差least_squares内部负责把残差平方和、雅可比矩阵和一阶最优性全部组织好。methodtrf因为设置了边界所以不能用默认的lm如果没有边界lm速度通常更快。max_nfev5000是残差函数最大调用次数作为安全阀。三个 tol 都压到 1e-12是为了让结果尽量接近真值实际工程里 1e-8 就够了压太低只会让迭代在数值噪声上来回试探。运行后result.x会回到 2.5、1.3、3.0、0.6 附近。result.cost约等于 0.5 乘以残差平方和可作为不同初值下结果的对比基准。提示如果发现result.optimality没有接近 0说明还没收敛优先加大max_nfev或放宽gtol不要急着改模型。3.3 初值、边界与目标函数形态的配合least_squares是局部优化初值的质量直接决定结果落在哪个局部极小点。我常用的策略如下场景初值策略边界设置物理参数拟合用物理量级粗略估算不用 0 做初值按量程设硬边界避免负数或奇异区外参标定用解析解或出厂标定值做初值角度给 ±π平移给量程动力学参数辨识多组随机初值并行求解无边界信息时优先加 L2 正则而不是硬边界对衰减正弦模型可以先从数据衰减趋势读出一个大概的 b用峰值位置估 c 的初始值。边界要依据参数语义给a 是幅值一定为正d 是相位角。给不出物理边界时宁可用宽边界加多起点也不要无边界裸跑否则 TRF 可能在无效区域积累数值误差。3.4 求解后的一分钟快速验证拟合参数跑出来后我不急着看 cost先看残差分布import matplotlib.pyplot as plt fit_params result.x res residuals(fit_params) print(np.mean(res)) print(np.percentile(np.abs(res), [50, 90, 99])) plt.plot(x_data, y_obs, ., labelobs) plt.plot(x_data, model(x_data, *fit_params), -, labelfit) plt.legend() plt.savefig(fit_result.png, dpi120)残差均值应该接近 0绝对值 90 分位数和 99 分位数应该大致在 2 倍和 3 倍噪声标准差附近。如果均值明显偏移说明模型有系统偏差如果高百分位过大说明存在异常样本下一步就该上鲁棒损失。这一步验证的是模型结构是否与观测匹配比任何优化器参数都重要。4. 多元非线性目标函数的通用形式、约束与梯度陷阱4.1 用 minimize 处理非最小二乘形式的目标函数真实场景里目标函数不总是干净的平方和。比如观测里有离群点时我会把残差平方和换成近似 Geman-McClure 的鲁棒损失让大残差的权重被压住from scipy.optimize import minimize def robust_misfit(theta): r residuals(theta) return np.sum(np.log1p(0.5 * r**2)) res minimize( robust_misfit, theta0, methodL-BFGS-B, boundsbounds, options{maxiter: 2000, ftol: 1e-12}, ) print(res.x)log1p(0.5 * r^2)对小的残差接近 0.5 * r^2对大的残差只按对数增长这比直接平方和鲁棒得多。这里目标函数不再是残差平方和结构least_squares的结构优势消失了只能退回到通用minimize。methodL-BFGS-B带边界且内存开销只跟迭代历史有关适合变量数几百以内的中等规模问题。4.2 解析梯度、数值梯度与自动微分怎么选minimize的jac参数不传时内部用有限差分估计梯度。对四五个变量的标定问题这完全够但变量超过 20 个后有限差分的舍入误差会明显拖慢收敛有时还会在平坦区域误判方向。手动推导解析梯度最可靠但链式法则展开很容易写错。我现在的习惯是能自动微分就用自动微分其次用手推解析梯度并配合check_grad校验最后才考虑纯有限差分。在 Python 工程里可以把目标函数写成 JAX 计算图用jax.grad得到解析梯度再传给 scipy 的求解器如果项目不允许引 JAX至少要把手推梯度写在独立函数里而不是内嵌在目标函数内部否则没法单独验证。4.3 线性约束与 QP 的边界在哪当目标函数变成二次、约束全部线性时问题已经不是通用多元非线性而是二次规划交给 QP 求解器会远比minimize快。判定标准是目标函数的 Hessian 是常数矩阵约束是 Ax ≤ b 这种形式。一个典型例子是带线性等式约束的最小二乘min 0.5 ||Ax - b||², s.t. Cx d这类问题在 QP 求解器里有专门的数据结构和预处理求解时间通常是通用内点法的零头。一般非线性等式约束才需要用trust-constr或 SLSQP非线性约束博弈的成本远高于边界约束。如果只是给参数设置上下限那不算真正意义上的约束优化least_squares的bounds和L-BFGS-B已经覆盖了。工程上尽量把约束先表达为边界表达不了再上通用约束这一条能让问题简单一个等级。4.4 多元非线性目标函数求解的典型失败排查顺序症状原因排查路径多次运行结果不同多峰初值敏感多起点扫描然后看 cost 最小值正常收敛但残差很大模型缺项或数据有离群点先画残差观察是否系统性偏移迭代震荡不收敛变量尺度差异大或目标函数不光滑做变量归一化检查模型是否出现除零已到 5000 次调用但没收敛容差太严或初值太远先放宽 ftol 到 1e-8跑通后再收紧报零主元或 singular matrix雅可比矩阵秩亏参数冗余检查哪些参数线性相关合并或固定雅可比矩阵秩亏是多元非线性目标函数求解里最难发现的坑。表面现象是求解器报 zero pivot 或者干脆 NaN本质是某些参数在数据里不可区分例如 a 与 b 同时乘以同一个常数仍能保持同样的残差。遇到这种情况我会固定其中一个参数或者往目标函数里加一个小量 L2 正则项让 Hessian 可逆。5. 多元非线性目标函数求解结果验证的三个硬技巧5.1 用 check_grad 校验解析梯度手写梯度后一定要做梯度校验。scipy.optimize.check_grad用有限差分对比解析梯度误差小于 1e-6 基本可以认为表达式正确from scipy.optimize import check_grad def f_xy(xy): x, y xy return (x - 2)**4 (x - 2*y)**2 def df_xy(xy): x, y xy return np.array([ 4*(x-2)**3 2*(x-2*y), -4*(x-2*y) ]) err check_grad(f_xy, df_xy, np.array([0.5, 0.8])) print(err)把这里的f_xy和df_xy替换成自己的目标函数与解析梯度即可。注意采样点不要选在任意分量为 0 的位置否则梯度表达式中的某些项会被掩盖。5.2 多起点扫描定位近似全局最优多元非线性目标函数求解的最优默认是局部最优通过多组随机初值扫描每组都收敛后取 cost 最小的结果可以近似逼近全局最优best None for seed in range(20): rs np.random.default_rng(seed) start np.array([ rs.uniform(0.5, 2.0), rs.uniform(1.0, 6.0), rs.uniform(2.0, 8.0), rs.uniform(-2.0, 2.0), ]) cand least_squares(residuals, start, methodtrf, boundsbounds) if best is None or cand.cost best.cost: best cand print(best.x, best.cost)如果 20 组初值里只有一两组收敛到同一个低 cost其余都落在不同高点说明目标函数多峰严重需要把扫描数量提到 50 组以上并用每个起点的方向分布判断参数空间的连通性。5.3 残差分布与一阶最优性决定何时停手最后检查残差的分布形态均值接近 0、绝对值 90 分位数在噪声的约 2 倍以内说明模型结构与数据一致如果高百分位突跳优先处理离群点而不是继续压容差。同时确认result.optimality已经低于收敛阈值这是求解器自己在告诉你一阶最优性条件已经满足。发布到自动化流水线之前把参数的边界值、残差 90 分位数和多起点扫描的初始组数一起写进验收日志。这个基线也是后续更换初值策略时判断是否变好的参照。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进