ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

用Ceres求解非线性优化问题

用Ceres求解非线性优化问题 用Ceres求解非线性优化问题核心分为四步1定义参数块2构建代价函数残差3添加残差块到Problem4配置求解器并求解以ye^(mxc)为例我们有一堆观测数据点 (,)还有一个模型 e^(mxc)。目标是找到 m,c让模型预测值尽量接近真实观测。接近程度用残差衡量viyi-e^(mxic)优化目标是让所有残差平方和最小min⁡m,c∑ivi2 \min_{m,c}\sum_{i} v_i^2m,cmin​i∑​vi2​整个优化过程可以这样理解Ceres通过不断调整优化变量Parameters使得所有残差Residual的平方和最小化第一步定义参数块Parameters定义需要求解未知数并给其一个初值doublem0.0;//优化从初值点开始doublec0.0;第二步定义代价函数CostFunction定义代价函数的核心就是告诉Ceres 我给了它一个参数xi,yi后它要怎么计算残差以及残差对优化变量的导数Jacobian。Ceres实现这一目标主要有三种方式1.自动求导 (Automatic Differentiation) - 最推荐的方式这是Ceres最强大、最常用的功能。你只需要计算残差Ceres会利用C模板template和一种叫“对偶数Dual Numbers”的数学技巧自动、精确地帮你算出导数雅可比矩阵。实现步骤1定义一个仿函数Functor类在这个类中实现一个模板化的operator()。2在operator()中计算残差使用模板类型T进行所有运算。3使用AutoDiffCostFunction包装将仿函数类传递给 AutoDiffCostFunction。代码示例#includeceres/ceres.h// 1. 定义仿函数structExponentialResidual{// 构造函数传入观测到的 (x, y) 数据对ExponentialResidual(doublex,doubley):x_(x),y_(y){}// 核心计算残差templatetypenameTbooloperator()(constT*constm,// 输入参数 m数组指针constT*constc,// 输入参数 cT*residual)const// 输出残差{// residual 观测值 - 模型预测值residual[0]T(y_)-exp(m[0]*T(x_)c[0]);returntrue;// 返回 true 表示计算成功}private:constdoublex_;// 这个点的 xconstdoubley_;// 这个点的 y};intmain(){// --- 第一步生成模拟数据真实值设为 m0.3, c0.1---std::vectordoublex_data;std::vectordoubley_data;doublem_true0.3;doublec_true0.1;for(inti0;i10;i){doublexi*0.5;// x 0, 0.5, 1.0, ...doubleystd::exp(m_true*xc_true);// 根据真实模型计算 y// 可选加一点噪声使问题更真实y 0.01 * (rand() / RAND_MAX);x_data.push_back(x);y_data.push_back(y);}// --- 第二步定义待优化的变量初始值---doublem0.0;// m 的初始猜测doublec0.0;// c 的初始猜测// --- 第三步构建优化问题---ceres::Problem problem;// 遍历所有数据点为每个点添加一个残差块ResidualBlockfor(inti0;ix_data.size();i){// 使用 AutoDiffCostFunction 封装代价函数// 模板参数含义// ExponentialResidual : 我们定义的仿函数类// 1 : 每个残差块中残差的维度这里输出 1 个残差值// 1 : 第 1 个参数块m的维度标量维度为1// 1 : 第 2 个参数块c的维度标量维度为1ceres::CostFunction*cost_functionnewceres::AutoDiffCostFunctionExponentialResidual,1,1,1(newExponentialResidual(x_data[i],y_data[i]));// 将残差块添加到问题中。// 参数代价函数损失函数nullptr表示使用标准二乘以及两个优化变量的地址problem.AddResidualBlock(cost_function,nullptr,m,c);}// --- 第四步配置求解器并求解---..................................return0;}关键点解析为什么用模板 template 这是 Ceres 自动求导的关键。Ceres 会传入一种特殊类型Jet在计算残差的同时自动算出导数雅可比你完全不用手推数学公式。所以代码里所有数字都要写成 T(…)运算也用模板类型。AutoDiffCostFunction 的模板参数第一个参数是你的仿函数类名 CostFunctor。第二个参数是残差的维数这里为1。从第三个参数开始是每个参数块的维数。这里只有一个参数块 x维数为1所以是 1。如果你的问题有多个参数比如 x 是3维向量y 是4维向量则写法是 AutoDiffCostFunctionFunctor, 残差维数, 3, 4。2.数值求导 (Numeric Differentiation)当你无法或不想使用自动求导时例如调用了无法模板化的第三方库可以使用这种方法。你只需计算残差Ceres会通过有限差分法如 (f(xh)-f(x))/h来数值近似地计算导数。实现步骤1定义一个仿函数Functor类实现一个非模板化的 operator()只计算残差。2使用 NumericDiffCostFunction 包装。// 1. 定义数值求导的仿函数注意不是模板structNumericExpResidual{NumericExpResidual(doublex,doubley):x_(x),y_(y){}// 直接使用 double 计算残差不需要模板booloperator()(constdouble*constm,constdouble*constc,double*residual)const{// 计算预测值doublepredictedstd::exp(m[0]*x_c[0]);// 残差 观测值 - 预测值residual[0]y_-predicted;returntrue;}private:constdoublex_;constdoubley_;};intmain(){// 生成数据std::vectordoublex_data,y_data;GenerateData(x_data,y_data);// 待优化的初始值doublem0.0;doublec0.0;ceres::Problem problem;for(inti0;ix_data.size();i){// 2. 使用 NumericDiffCostFunction 创建代价函数// 模板参数含义// NumericExpResidual : 仿函数类// ceres::CENTRAL : 中心差分法比 FORWARD 精确// 1 : 残差维度// 1 : 第1个参数块 (m) 的维度// 1 : 第2个参数块 (c) 的维度ceres::CostFunction*cost_functionnewceres::NumericDiffCostFunctionNumericExpResidual,ceres::CENTRAL,1,1,1(newNumericExpResidual(x_data[i],y_data[i]));problem.AddResidualBlock(cost_function,nullptr,m,c);}// 求解器配置与求解ceres::Solver::Options options;options.linear_solver_typeceres::DENSE_QR;options.minimizer_progress_to_stdouttrue;ceres::Solver::Summary summary;ceres::Solve(options,problem,summary);std::coutsummary.BriefReport()\n;std::cout估计 m m, 估计 c cstd::endl;return0;}核心要点仿函数中的 operator() 不是模板函数直接使用 double 计算。使用ceres::NumericDiffCostFunction 封装。模板参数中必须指定差分方式通常选ceres::CENTRAL中心差分精度更高。精度损失因为用近似值代替精确导数计算结果会有微小的数值误差通常 10−6级别。速度较慢每次求导需要多次计算残差中心差分需要计算 2 次迭代速度比自动求导慢。无需推导这是它最大的优点当你面对极其复杂、无法模板化的第三方库函数时只能选它。3.解析求导Analytic Differentiation这种方法需要你亲自推导出残差对每个优化变量的偏导数公式然后手动填入雅可比矩阵。对于本示例vy−emxc vy-e^{mxc}vy−emxc残差v对于m的偏导数为∂v∂m−xemxc\frac{\partial v}{\partial m}-xe^{mxc}∂m∂v​−xemxc残差v对于c的偏导数为∂v∂c−emxc\frac{\partial v}{\partial c}-e^{mxc}∂c∂v​−emxc实现步骤1继承 SizedCostFunction 类在模板参数中指定残差和参数块的维数。2重写 Evaluate() 函数在里面计算残差和雅可比矩阵。代码示例// 1. 继承 SizedCostFunction// 模板参数: 残差维度, 第1参数块维度, 第2参数块维度classAnalyticExpResidual:publicceres::SizedCostFunction1,1,1{public:AnalyticExpResidual(doublex,doubley):x_(x),y_(y){}virtual~AnalyticExpResidual(){}// 2. 重写 Evaluate 函数virtualboolEvaluate(doubleconst*const*parameters,double*residuals,double**jacobians)const{// 从参数指针数组中提取 m 和 c 的值constdoublemparameters[0][0];constdoublecparameters[1][0];// 计算指数值复用避免重复计算doubleexp_valstd::exp(m*x_c);// --- 3. 计算残差 ---residuals[0]y_-exp_val;// --- 4. 计算雅可比矩阵导数---// 注意Ceres 可能不需要导数比如仅用于评估所以必须检查指针是否为空if(jacobians!nullptr){// 残差对 m 的偏导数: dv/dm -x * e^(mxc)if(jacobians[0]!nullptr){jacobians[0][0]-x_*exp_val;}// 残差对 c 的偏导数: dc/dc -e^(mxc)if(jacobians[1]!nullptr){jacobians[1][0]-exp_val;}}returntrue;}private:constdoublex_;constdoubley_;};intmain(){std::vectordoublex_data,y_data;GenerateData(x_data,y_data);doublem0.0;doublec0.0;ceres::Problem problem;for(inti0;ix_data.size();i){// 直接 new 我们自定义的代价函数类ceres::CostFunction*cost_functionnewAnalyticExpResidual(x_data[i],y_data[i]);problem.AddResidualBlock(cost_function,nullptr,m,c);}// 求解器配置与求解与之前完全一样ceres::Solver::Options options;options.linear_solver_typeceres::DENSE_QR;options.minimizer_progress_to_stdouttrue;ceres::Solver::Summary summary;ceres::Solve(options,problem,summary);std::coutsummary.BriefReport()\n;std::cout估计 m m, 估计 c cstd::endl;return0;}参数指针数组 parameters它是一个二级指针。parameters[0] 指向 mparameters[1] 指向 c。必须用 parameters[0][0] 解引用取出值。雅可比指针数组 jacobians它的结构和 parameters 对应。jacobians[0] 对应第一个参数块 m 的雅可比矩阵。由于这里 m 是标量维度1残差也是标量维度1所以矩阵大小为 1×1赋值给 jacobians[0][0]。必须检查空指针if (jacobians ! nullptr) 和 if (jacobians[0] ! nullptr) 是强制且必须的。因为 Ceres 在某些情况下比如只做损失函数评估不会分配雅可比内存如果你直接赋值程序会崩溃。第三步创建 Problem 并添加残差块Problem 是整个优化问题的容器。需要把每个数据点的残差都添加进去ceres::Problem problem;// 创建问题容器// 假设有 kNumObservations 个数据点存在 data 数组里for(inti0;ikNumObservations;i){// 3.1 创建代价函数对象ceres::CostFunction*cost_functionnewceres::AutoDiffCostFunctionExponentialResidual,1,1,1(newExponentialResidual(data_x[i],data_y[i]));// 3.2 把这个残差加入问题problem.AddResidualBlock(cost_function,nullptr,// 损失函数先不用鲁棒核m,c);// 关联的参数}第四步配置求解器并求解Ceres Solver的配置主要通过 ceres::Solver::Options 结构体来控制决定了优化过程如何执行。这些选项可以大致分为三类通用控制选项、核心算法选项和线性求解器选项。1、通用控制选项管理求解过程2.核心算法选项选择优化策略这类选项决定了Ceres使用哪种数学策略来寻找最优解。1. 选择最小化算法minimizer_typeCeres提供两种主流的非线性优化算法TRUST_REGION (信赖域方法)这是默认选项也是大多数情况下的首选。它在每一步先选择一个认为模型有效的“信任区域”然后在该区域内寻找最优步长和方向。LINE_SEARCH (线搜索方法)这种方法先确定一个下降方向然后沿着这个方向寻找最优的步长。在某些特定问题上可能更高效。2. 选择信赖域算法trust_region_strategy_type当选择 TRUST_REGION 后还可以进一步选择具体的算法LEVENBERG_MARQUARDT默认选项也是最常用、最稳定的算法适用于绝大多数问题。DOGLEG另一种信赖域算法在某些问题上可能比LM算法更快。3.线性求解器选项决定求解速度和内存这是Ceres选项中最关键也最复杂的一部分。在每次迭代中Ceres都需要求解一个大型线性方程组linear_solver_type 就用于指定求解这个方程组的方法如何选择不确定时先用默认的 DENSE_QR。参数很少如 100DENSE_QR 或 DENSE_NORMAL_CHOLESKY 都很好。参数很多但结构稀疏可以尝试 SPARSE_NORMAL_CHOLESKY。这需要你的Ceres编译时支持SuiteSparse或CXSparse等稀疏库。正在做SLAM或SfM这是BA问题应优先考虑 SPARSE_SCHUR。如果没有稀疏库DENSE_SCHUR 也可用于小规模测试。// 求解器配置与求解与之前完全一样ceres::Solver::Options options;options.linear_solver_typeceres::DENSE_QR;options.minimizer_progress_to_stdouttrue;ceres::Solver::Summary summary;ceres::Solve(options,problem,summary);//打印报告std::coutsummary.BriefReport()\n;std::cout估计 m m, 估计 c cstd::endl;综上就完成了一个最简单的非线性优化过程。在SLAM用ceres估计位姿时还会涉及到一些其他概念流形Manifold是一个数学概念直观理解就是一个看起来弯曲的空间但在每一个局部小范围内它又近似平坦像普通的欧几里得空间旋转正好处在这样的弯曲空间里不能像普通向量那样直接加减需要对它的运算进行重新定义。下面是VINS-MONO中处理位姿的流形专门处理位置3维四元数4维这种7维存储但只有6自由度的位姿定义一个流行需要实现代码中所示6个函数//Manifold是基类PoseLocalParameterization是其派生类重写了Plus、PlusJacobian、Minus、MinusJacobian等方法classPoseLocalParameterization:publicceres::Manifold{// ① 加法核心定义 x delta 怎么算virtualboolPlus(constdouble*x,constdouble*delta,double*x_plus_delta)constoverride;// ② 加法的雅可比Plus 对 delta 在 delta0 处求导virtualboolPlusJacobian(constdouble*x,double*jacobian)constoverride;// ③ 减法定义 y - x 怎么算Plus 的逆运算virtualboolMinus(constdouble*y,constdouble*x,double*y_minus_x)constoverride;// ④ 减法的雅可比virtualboolMinusJacobian(constdouble*x,double*jacobian)constoverride;// ⑤ 存储维度环境空间维度位姿7virtualintAmbientSize()constoverride{return7;}// ⑥ 自由度维度切空间维度位姿6virtualintTangentSize()constoverride{return6;}};
返回列表