Ceres Solver实战:从最小二乘到流形优化的完整上手笔记
做三维视觉这行手头没一个好用的非线性优化库就像修车没有扳手。Ceres Solver不是最炫的库但胜在稳、快、资料齐全尤其适合处理目标函数能写成误差平方和的问题——相机标定、点云配准、激光建图、视觉惯性里程计里的姿态图优化甚至纯函数曲线拟合都能往里塞。这篇笔记是我自己从零啃Ceres Solver的记录尽量把公式翻译成人话。目标是让被高深术语劝退的读者有一条更平缓的上手路径先搞懂它到底在解决什么问题再拆核心概念然后跑通一个真实例子最后聊调参和排错。不管你是刚接触SLAM的学生还是要做标定或拟合工程的开发者按这个顺序读下来至少能把代码跑通并且知道下一步该去查哪里。1. Ceres Solver到底在解决什么问题一个公式看懂非线性最小二乘1.1 最小二乘问题从“算得差不多”到“算得最准”很多工程问题的本质可以概括成一句话有一堆带噪声的观测数据想求一组参数让模型尽量符合这些数据。用数学语言写出来Ceres要解决的优化问题长这样minimize 0.5 * Σ ρᵢ( || fᵢ(x) ||² )先别被符号吓到。这里的 x 是我们要估计的参数比如相机的内外参、机器人的位姿、拟合曲线的系数。fᵢ(x) 是“误差项”表示模型预测值和观测值之间的差距。整个式子的意思就是找一组 x让所有误差的平方和最小。为什么是平方而不是绝对值因为平方函数在零点附近光滑可导数学上处理方便而且对小幅噪声的统计假设很自然。这就是“最小二乘”二字的由来。有了这个框架你会发现好多任务都能套进来。比如给定一堆三维点和它们对应的二维像素坐标估计相机位姿PnP问题给定一组里程计数据和回环检测结果优化机器人的轨迹姿态图优化给出一组温度和电阻的测量记录拟合成一个多项式。1.2 怎么去求最小从暴力搜索到梯度下山有了目标函数关键问题是怎么找最小值最容易想到的办法是网格搜索把参数空间每个组合都试一遍。但参数一多组合数指数爆炸根本不现实。比如位姿有6个自由度哪怕每个自由度只试100个值就是10的12次方种组合。工程上更聪明的方法是迭代求解。想象你站在山坡上想走到山谷最低点。最朴素的做法是每走一步看看周围哪个方向是下坡然后朝那个方向走一步不断重复。这就是梯度下降法的直观思想。但纯梯度下降有两个问题一是接近最优值时收敛越来越慢二是它只用了目标函数的一阶信息没有考虑曲率步子迈大了容易震荡迈小了又太慢。后来人们发现如果目标函数本身是平方和形式可以利用误差项的结构——误差函数的一阶导Jacobian和二阶近似Hessian——构造出更聪明的步长。这就是高斯牛顿法和列文伯格-马夸尔特法Levenberg-Marquardt。Ceres内部默认用的就是LM策略配合各种线性求解器能把几十万维的优化问题在几秒内解决。1.3 Ceres在这个生态里的位置写优化算法不难难的是把数值稳定性、稀疏矩阵、自动求导、并行计算这些细节都处理好。Ceres把这些脏活累活全包了你只需要定义误差函数告诉它误差怎么算。告诉它哪些是待优化的参数。调几个求解选项然后喊“开始”。它内部会自动求导或使用你给的导数选合适的线性求解器迭代到满足精度条件为止。它的典型应用场景非常集中场景待优化参数误差来源多项式拟合多项式系数拟合值与观测值之差相机标定内参、畸变系数、外参重投影误差点云配准旋转矩阵平移向量对应点距离误差SLAM后端优化位姿序列、路标点观测误差和回环约束误差传感器外参标定两传感器间的位姿变换同步观测数据之间的偏差一句话总结凡是能写成“误差平方和最小”的问题都可以先想到Ceres。2. 核心概念拆解Problem、CostFunction、LossFunction 与参数块2.1 用组装积木的方式理解这四个组件Ceres最关键的抽象可以分四块参数块ParameterBlock一个待优化的变量数组。比如一个位姿通常是6维或7维一个相机内参可能是4维或5维。代价函数CostFunction给定参数值计算出误差向量。这是整个优化问题的核心你可以用自动求导、数值求导或解析求导来定义它。残差块ResidualBlock把代价函数、参数块、可选的核函数绑在一起形成目标函数中的一个项。一个优化问题通常有几百甚至几十万个残差块。损失函数LossFunction也叫核函数用来调节误差的影响。普通最小二乘对所有误差一视同仁但数据里往往有离群点一个错误匹配就能把结果带偏核函数能压低大误差的权重。2.2 数据流先建问题再求解Ceres的使用流程几乎固定不变构造一个 ceres::Problem 对象。循环添加残差块AddResidualBlock(cost_function, loss_function, 参数块指针...)。设置参数块的可优化状态SetParameterBlockConstant固定某些参数或 SetParameterization/SetManifold处理特殊参数空间。构造 Solver::Options配置求解器偏好。调用 Solve(options, problem, summary) 求解。从summary里读迭代信息从参数块指针里取优化结果。有一件事很多新手容易搞混Ceres不会返回一个新的参数数组它直接原地更新你传给它的指针指向的内存。所以求解完参数就自动变成优化后的值了。2.3 为什么必须要有LossFunction聊聊野值我最初用Ceres做拟合时觉得LossFunction是多余的。普通最小二乘不是挺好的吗后来被真实数据狠狠教育了一回。假如你在拟合一条曲线数据里混入了一个偏移极大的坏点。普通最小二乘会对这个点的误差平方项施加巨大惩罚整个拟合结果被它拽走看起来就像一根好好直尺被一个钉子顶弯了。LossFunction的作用就是给某些误差加了上限把大误差的影响从“平方级”降成“线性级”甚至更弱。常用的核包括Huber核和Cauchy核。Huber核在小误差区段保持平方惩罚大误差区段变成线性惩罚相当于给野值设置了一个“忍耐阈值”。在标定和SLAM里不加核函数的优化几乎不可用——因为特征匹配永远有错误匹配。2.4 参数空间不是普通的加减法空间对普通数组比如多项式的三个系数加减法天然成立。但对旋转矩阵、四元数这类参数情况就麻烦了。四元数可以表示旋转但它有一个约束模长必须是1。如果你把四元数当普通4维向量去加减更新后的模长就不再是1它就不再表示一个合法的旋转。这就需要用LocalParameterization新版叫Manifold来告诉Ceres你的参数虽然存在4维数组里但实际运动只有3个自由度更新的时候要用特殊的运算规则。这个点特别重要我在第5节会详细展开先在这里建立一个概念不是所有参数都适合直接加一个增量。3. 写第一个能跑的Ceres程序曲线拟合实战3.1 环境准备与编译工程的搭建Ceres本身是C库依赖线性代数库Eigen。编译安装的工作建议放在一个独立环境里做避免污染系统。最省心的方式是用系统包管理器直接安装预编译版本。如果你用的是常见Linux发行版或macOS一条命令就能装好包名通常类似libceres-dev。装好之后头文件在include/ceres库在lib目录。Windows下推荐用vcpkg同样能直接拉取。我自己的习惯是源码编译下载源码安装依赖Eigen、一些命令行解析、日志相关库然后依次跑cmake、make、install。源码编译的好处是可以开启更多测试用例也能跟进新特性。代价是第一次配置时会碰到各种依赖报错耐心一点即可。新建工程时CMakeLists.txt几乎是模板写好一次到处复用cmake_minimum_required(VERSION 3.10) project(curve_fitting) set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) find_package(Ceres REQUIRED) add_executable(curve_fitting main.cpp) target_link_libraries(curve_fitting Ceres::ceres)3.2 用自动求导定义误差函数我们要拟合的目标是二次曲线 y ax² bx c。先构造带噪声的观测数据然后用Ceres找回系数。用AutoDiff自动求导定义代价函数核心是写一个继承自 ceres::AutoDiffCostFunction 的类型。重点在于 operator() 必须是模板函数因为Ceres要拿不同的标量类型double和Jet分别计算值和导数#include ceres/ceres.h #include cmath #include random #include vector struct QuadraticCostFunctor { QuadraticCostFunctor(double x, double y) : x_(x), y_(y) {} template typename T bool operator()(const T* const params, T* residual) const { residual[0] y_ - (params[0] * x_ * x_ params[1] * x_ params[2]); return true; } private: double x_; double y_; };error 真实观测值 - 预测值。取谁的符号其实不影响最小化因为后面会对误差平方正负无所谓。3.3 构造Problem并求解主函数逻辑如下生成一批模拟数据。设真值为 a1.0, b2.0, c0.5在 x 从-1到1均匀采样y按真实公式计算再加一点高斯噪声。初始化参数为接近零的随机值比如 a0, b0, c0。构造problem。对每个数据点新建一个代价函数对象加入残差块。loss这里先用NULL等价于不加核等加了野值再换Huber。设置求解选项。第一次跑用默认配置就行但把输出打开看迭代过程。完整代码如下int main(int argc, char** argv) { double true_a 1.0, true_b 2.0, true_c 0.5; double init_a 0.0, init_b 0.0, init_c 0.0; std::vectordouble x_data, y_data; std::mt19937 rng(42); std::normal_distributiondouble noise(0.0, 0.2); for (int i 0; i 100; i) { double x -1.0 2.0 * i / 99.0; double y true_a * x * x true_b * x true_c noise(rng); x_data.push_back(x); y_data.push_back(y); } ceres::Problem problem; double params[3] {init_a, init_b, init_c}; for (int i 0; i x_data.size(); i) { ceres::CostFunction* cost new ceres::AutoDiffCostFunctionQuadraticCostFunctor, 1, 3( new QuadraticCostFunctor(x_data[i], y_data[i])); problem.AddResidualBlock(cost, nullptr, params); } ceres::Solver::Options options; options.max_num_iterations 100; options.linear_solver_type ceres::DENSE_QR; options.minimizer_progress_to_stdout true; ceres::Solver::Summary summary; ceres::Solve(options, problem, summary); std::cout summary.BriefReport() std::endl; std::cout estimated a,b,c: params[0] , params[1] , params[2] std::endl; return 0; }这里选 DENSE_QR 是因为只有3个参数稠密求解器足够。如果参数数量成千上万就得换稀疏求解器。3.4 看求解日志别忽略Summary这个宝库跑完程序控制台会打印迭代日志和BriefReport。日志里每一行对应一次迭代显示iter第几次迭代cost当前总代价梯度模长、步长、信赖域半径等内部指标很多新人只看“estimated”结果不看迭代过程这是不对的。如果最终迭代停在“max iterations reached”而不是“converged”说明你设的迭代上限太小或者参数设置有问题。如果出现cost剧烈震荡甚至NaN就要回头检查误差函数和初始值了。4. 三种求导方式怎么选AutoDiff、数值求导和解析求导4.1 为什么要求导无论哪种优化算法都需要知道目标函数在当前位置的梯度信息。Ceres把优化框架搭好了但梯度信息怎么来这就是求导方式要解决的问题。Ceres支持三种方式自动求导AutoDiff你只写原函数它用模板魔法算出精确导数。数值求导NumericDiff用有限差分近似导数。解析求导AnalyticDiff你手写出导数的解析表达式。隐藏的一点是同样一个问题求导方式的差异会直接体现在收敛速度和稳定上。4.2 自动求导的工作原理Jet是个好东西AutoDiff不是差分它是把普通double替换成一种叫Jet的类型。Jet同时存着一个值和一组导数分量。运算符重载让所有数学操作都自动携带导数信息。误差函数里的乘法和加法执行完后残差值有了关于所有参数的偏导数也顺带算出来了。好处是导数精确到机器精度没有截断误差也不会出现数值求导那种步长怎么选的纠结。代价是模板代码写起来有点绕特别是嵌套复杂函数时模板层级变深编译时间会变长。我个人的经验是默认优先用AutoDiff。代码直观不容易出错适合绝大多数场合。4.3 数值求导调试友好但精度有限数值求导有现成封装你只需要提供普通的非模板operator()。它的原理是给参数加微小扰动计算差分。好处是写起来最简单坏处是步长选择很敏感步长太大截断误差大步长太小浮点数舍入误差大。还有一个麻烦对于几十万维参数加残差的组合数值求导意味着每个残差块都要额外计算多次原函数计算量成倍增加。所以在生产代码里我很少用数值求导只会在调试某个拟合问题时临时用一下。4.4 解析求导性能王者代价是易错手写导数就像手写汇编。能拿到最快的执行速度和最小的内存占用在实时性要求极高的场合比如实时SLAM前端值得考虑。但手推Jacobian不仅费时间还容易推错。一个符号错了结果就在误差平面上乱飘。实践中建议这样定策略求导方式编码难度精度性能推荐程度自动求导中精确良好首选数值求导低有近似误差差仅调试解析求导高精确最优性能瓶颈时再用如果刚开始学先把AutoDiff用熟后面觉得单次求解耗时压不下来再考虑针对最热点的残差块手写解析导数。5. 姿态优化与四元数中文笔记里最容易翻车的部分5.1 一个数值上合法、几何上荒谬的四元数在SLAM里待优化变量经常是世界坐标系下的机器人位姿。表示旋转用的通常是四元数。四元数要满足模长为1的约束。问题是Ceres的迭代总是给参数加增量比如 x_new x delta。如果直接对四元数的4个分量做这种加法更新后的模长大概率不再是1。一个模长不是1的四元数对应着一个不存在的旋转。这种错误很阴险计算过程中所有数值都正常没有NaN但旋转矩阵变得畸形重投影误差爆炸或者优化结果完全离谱。5.2 LocalParameterization的核心接口老版本Ceres用 ceres::LocalParameterization 来定义“在切空间上更新”。假设参数块是4维的四元数但真正的自由度只有3维。Ceres在迭代时不会直接加一个4维增量而是把当前4维向量和增量3维交给Plus函数。Plus函数内部做3维到4维的映射执行四元数乘法再做归一化。你得自己实现Plus和Jacobian两个函数。简单说Plus定义了“如何在流形上迈出一步”Jacobian定义了“这步的导数怎么传到残差”。一段极简的四元数流形实现通常长这样class QuaternionParameterization : public ceres::LocalParameterization { public: bool Plus(const double* x, const double* delta, double* x_plus_delta) const override { // 将delta的3维小角增量转为四元数 Eigen::Quaterniond dq; double theta delta[0] * delta[0] delta[1] * delta[1] delta[2] * delta[2]; if (theta 0.0) { theta std::sqrt(theta); dq Eigen::Quaterniond( std::cos(theta / 2.0), delta[0] / theta * std::sin(theta / 2.0), delta[1] / theta * std::sin(theta / 2.0), delta[2] / theta * std::sin(theta / 2.0)); } else { dq Eigen::Quaterniond(1.0, 0.0, 0.0, 0.0); } Eigen::Quaterniond orig(x[3], x[0], x[1], x[2]); Eigen::Quaterniond result (dq * orig).normalized(); x_plus_delta[0] result.x(); x_plus_delta[1] result.y(); x_plus_delta[2] result.z(); x_plus_delta[3] result.w(); return true; } bool ComputeJacobian(const double* x, double* jacobian) const override { // 通常可以让Ceres用自动求导来数值雅可比或者手动填4x3矩阵 // 这里简化处理大多数情况下填近似恒等映射对应的矩阵也能工作 Eigen::MapEigen::Matrixdouble, 4, 3, Eigen::RowMajor J(jacobian); Eigen::Quaterniond q(x[3], x[0], x[1], x[2]); // 实际工程中这里需要推导完整Jacobian这里只做示意 J.setZero(); return true; } int GlobalSize() const override { return 4; } int LocalSize() const override { return 3; } };需要说明的是手动实现ComputeJacobian很容易写错。在Ceres较新的版本里推荐改用Manifold接口并且可以配合AutoDiff来自动生成流形上的Jacobian大幅降低手推公式的工作量。5.3 新版Manifold接口接口变化怎么办Ceres新版本把 LocalParameterization 标记为废弃改用 ceres::Manifold。Manifold的核心接口更清晰Plus提供从切空间到流形空间的更新Minus提供从流形空间到切空间的差此外还有两个Jacobian接口。对于四元数新版本里已经内置了 QuaternionManifold 这个东西。也就是说你大概率不用自己写四元数流形了直接用内置的就行ceres::Manifold* quat_manifold new ceres::QuaternionManifold(); problem.SetManifold(q_ptr, quat_manifold);这是我在实际项目里最常做的事情。内置版本经过大量测试比自己实现的版本稳得多。只有当你要处理更复杂的空间比如SE(3)、SO(3)尺度的联合空间时才有必要写自定义Manifold。5.4 姿态图优化的简化思路姿态图优化是SLAM后端常见任务。问题是有一串位姿节点相邻节点之间有相对位姿约束比如来自里程计回环处有绝对约束来自回环检测想让全体位姿尽可能满足这些约束。用Ceres描述这个问题每个位姿可以拆成两部分平移向量3维普通加法更新和旋转四元数4维套Manifold更新。误差函数根据两个节点的相对位姿预测值与观测值之差来构造。核心就在于四元数残差不能简单做四则运算而要用旋转的“差”语义。写这个例子的完整代码篇幅太长不在此展开。但关键路径很清晰先定义位姿对应的残差块再给旋转部分配Manifold最后让求解器迭代。如果日志显示cost持续下降但最终姿态图有零空间漂移往往是因为没有给某个节点固定约束SetParameterBlockConstant。6. 求解器配置与调参Options里这些参数值得细调6.1 先认识这些常用选项Ceres Solver::Options里的参数非常多绝大多数默认值已经不错。新手不需要全部理解记牢这几个就够用linear_solver_type线性求解器直接影响求解速度和内存占用。trust_region_strategy信赖域策略一般选LEVENBERG_MARQUARDT。max_num_iterations最大迭代次数收敛失败经常是这里太小。function_tolerance / gradient_tolerance收敛判断阈值默认值较小调大可以提前结束。num_threads并行线程数残差块之间相互独立适合开多线程。minimizer_progress_to_stdout调试时设为true看实时迭代信息。6.2 线性求解器怎么选DENSE还是SPARSE每次迭代里核心工作是解一个线性方程组。参数少时可以用稠密求解器简单粗暴。参数多时矩阵是稀疏的每个残差只连到少数几个参数就必须用稀疏求解器。经验法则问题规模推荐求解器参数少于几十DENSE_QR 或 DENSE_SCHUR参数中等且结构明显SPARSE_SCHUR大规模无特殊结构SPARSE_NORMAL_CHOLESKY几十万参数的BA问题SPARSE_SCHUR 选项调整我曾在一个规模很大的BA问题上用了DENSE_QR结果内存直接爆炸——稠密矩阵要存几十万乘几十万的元素。换成SPARSE_SCHUR以后内存占用下降了一个数量级速度也快得多。选求解器之前先估算问题矩阵的稀疏度比盲目调参有效得多。6.3 LM还是Dogleg什么时候换一种信赖域策略LM列文伯格-马夸尔特是默认且最通用的选择。它在梯度下降和高斯牛顿之间自适应插值数值稳定性好。Dogleg是另一种信赖域方法在某些问题中收敛更快但对初始值更敏感。绝大多数情况下直接用LM只有当你发现收敛太慢且问题比较平滑时可以试试Dogleg。注意一点LM里的初始信赖域半径不是越大越好。半径太大第一步可能直接飞出去半径太小收敛极慢。Ceres会根据代价函数的前几步自动调整所以不用太担心初始值。6.4 从日志判断收敛状态cost下降曲线会说话我用Ceres调试问题时一定会打开minimizer_progress_to_stdout。迭代日志里有几个信号要关注每步cost都在稳定下降最后梯度接近0是健康状态。cost前几轮大幅下降后面小幅波动说明初始值还行但存在轻微病态。cost突然增大说明信任域半径给大了或代价函数里有不连续处。出现nan或inf优先级最高立刻停止查找误差函数、参数更新逻辑或求导定义。很多时候“不收敛”不是Ceres的问题而是误差函数本身定义得不合理或者参数空间更新方式不对。先看日志再检查模型最后才调求解器参数。7. 踩坑与排错速查我在实战中遇到的那些问题7.1 编译期错误模板地狱与Eigen对齐用AutoDiff写复杂代价函数时编译报错往往是一长串模板实例化日志看着像天书。我遇到最多的几类残差维度或参数块维度写错Ceres在编译期用静态断言检查维数错误信息会提示expected和actual。参数块数量不一致代价函数模板声明了3个参数块AddResidualBlock却传了两个。忘写模板关键字operator()里有时需要显式写typename或template。Eigen固定大小类型的内存对齐在结构体里包含Eigen::Quaterniond或Matrix时可能要加EIGEN_MAKE_ALIGNED_OPERATOR_NEW。7.2 运行期NaN先查三个地方求解过程出现NaN按下面顺序排查命中率极高除零错误。检查残差里有没有除以接近0的变量。根号或对数语义错误。比如对负数开根号。参数流形实现错误。四元数Plus之后没有归一化几轮迭代后数值迅速膨胀。我调试NaN时常用的手段是在代价函数里加printf把输入参数和中间量都打出来定位到具体哪个残差块先出问题。7.3 收敛到错误极值野值、初始值、不可观拟合结果明显不对但日志显示已收敛。这时多半是下面几种情况数据里存在强野值且你没加LossFunction。初始值离真值太远陷入了局部极小值。问题本身有不可观方向。比如单目SLAM存在尺度漂移你不固定某些参数优化结果自然在不可观方向上随意漂移。对应策略加Huber核用更可靠的初值比如先用PnP粗算一遍固定零空间SetParameterBlockConstant把关键节点钉住。7.4 性能上不去残差块过大或并行没开优化一个大规模问题时即使参数不多总耗时也可能惊人。常见原因是单个残差块计算太重或者线程没开满。Ceres的并行是基于残差块级别的残差块数量越大于线程数并行收益越明显。如果你构造了一个巨大的残差块包含几万项并行基本上失效应该拆分成成千上万个小的残差块。另外如果代价函数里有大量重复计算比如多个残差块共享同一个中间量可以用代价函数的构造函数把这些中间量先算好存下来避免每次求值都重算。7.5 踩坑快查表现象可能原因建议动作编译报模板错误维度不匹配检查ResidualBlock参数个数运行出NaN除零/根号负数/流形未归一化在代价函数里打印中间值不收敛但cost降幅小初始值差用其他方法提供更接近的初值结果被野值带偏没有LossFunction换HuberLoss出现不可观漂移参数自由度语义有冗余固定部分参数单次迭代很慢残差块过大/求解器不匹配拆分残差块换稀疏求解器四元数结果异常更新时直接加了4维增量使用内置QuaternionManifold最后补一句处理问题的通用思路不要让Ceres替你回答“模型对不对”它只能告诉你“在你给的模型下什么参数让误差最小”。误差函数写错再好的求解器也救不回来。把这句刻在脑子里遇到玄学问题先回去审模型本身。我个人最常用的一条调参技巧是第一次跑任何新问题先开max_num_iterations到比较大的值再把minimizer_progress_to_stdout打开全程盯日志。如果cost降得顺利再逐步收紧收敛阈值。别一上来就追求一次调到完美参数先让问题跑通再谈优化。Ceres这套框架一旦用熟你会发现大部分体力活都被它兜住了剩下的精力可以花在真正有挑战的建模和特征处理上。