Ceres Solver 实战：从非线性最小二乘问题到高效优化解决方案-开发者社区

1. 初识Ceres Solver：非线性优化的瑞士军刀

第一次接触Ceres Solver是在处理机器人定位问题时遇到的。当时需要优化一组传感器观测数据，传统的最小二乘法在非线性场景下表现不佳，直到发现了这个由Google开源的C++库。Ceres Solver就像一把精密的瑞士军刀，专门解决各种复杂的非线性最小二乘问题。

什么是非线性最小二乘问题？简单来说，就是找到一组参数，使得模型预测值与实际观测值之间的差异平方和最小。这类问题在工程实践中无处不在：从相机标定中的参数优化，到金融领域的曲线拟合，甚至机器学习中的模型训练，本质上都是在求解最小二乘问题。

Ceres Solver的核心优势在于它的自动微分能力。传统方法需要手动推导雅可比矩阵，不仅容易出错，而且当模型变更时需要重新推导。而Ceres只需要你定义残差如何计算，它能自动计算导数，这大大提高了开发效率。我在实际项目中最喜欢的就是这个特性——修改模型后不用再头疼于求导公式的推导。

安装Ceres非常简单，在Ubuntu上一条命令即可：

sudo apt-get install libceres-dev

如果是其他系统或者需要最新版本，也可以从源码编译：

git clone https://ceres-solver.googlesource.com/ceres-solver mkdir ceres-bin && cd ceres-bin cmake ../ceres-solver make -j8 sudo make install

2. 从简单例子理解Ceres工作流程

2.1 Hello World示例：直线拟合

让我们从一个最简单的线性拟合问题开始。假设有一组二维数据点，需要拟合一条直线y = mx + c。使用Ceres解决这个问题只需要三个步骤：

定义残差计算：每个数据点的残差就是观测y值与预测y值的差
构建优化问题：将所有数据点的残差相加
配置并运行求解器

具体代码实现如下：

#include <ceres/ceres.h> #include <vector> struct LinearResidual { LinearResidual(double x, double y) : x_(x), y_(y) {} template <typename T> bool operator()(const T* const m, const T* const c, T* residual) const { residual[0] = y_ - (m[0] * x_ + c[0]); return true; } private: const double x_, y_; }; int main() { // 生成带噪声的线性数据 y = 0.5x + 2.0 std::vector<double> x_data, y_data; for (int i = 0; i < 100; ++i) { double x = i / 10.0; x_data.push_back(x); y_data.push_back(0.5 * x + 2.0 + 0.1 * rand() / RAND_MAX); } // 初始参数值 double m = 0.0, c = 0.0; // 构建问题 ceres::Problem problem; for (size_t i = 0; i < x_data.size(); ++i) { ceres::CostFunction* cost_function = new ceres::AutoDiffCostFunction<LinearResidual, 1, 1, 1>( new LinearResidual(x_data[i], y_data[i])); problem.AddResidualBlock(cost_function, nullptr, &m, &c); } // 求解配置 ceres::Solver::Options options; 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() << "\n"; std::cout << "Final m: " << m << " c: " << c << "\n"; return 0; }

这个例子虽然简单，但包含了Ceres使用的核心模式。AutoDiffCostFunction模板会自动计算导数，我们只需要关注残差的计算逻辑。在实际项目中，这种模式可以扩展到更复杂的模型而保持代码结构清晰。

2.2 非线性示例：指数曲线拟合

非线性问题的处理方式类似，只是残差计算更复杂。考虑拟合指数曲线y = exp(mx + c)：

struct ExponentialResidual { ExponentialResidual(double x, double y) : x_(x), y_(y) {} template <typename T> bool operator()(const T* const m, const T* const c, T* residual) const { residual[0] = y_ - exp(m[0] * x_ + c[0]); return true; } private: const double x_, y_; };

注意这里使用了exp函数而非std::exp，这是Ceres自动微分的关键。所有数学运算都应使用Ceres提供的版本，包括log、sin、cos等。我在项目中曾经因为误用std::sqrt导致导数计算错误，调试了很久才发现这个问题。

3. 高级特性实战：鲁棒估计与配置技巧

3.1 处理异常值：鲁棒核函数

真实数据中常有异常值，普通最小二乘对异常很敏感。Ceres提供了多种鲁棒核函数（Loss Function）来降低异常值影响。例如使用Huber损失：

problem.AddResidualBlock(cost_function, new ceres::HuberLoss(1.0), // 鲁棒核函数 &m, &c);

Huber损失在残差较小时保持平方特性，较大时变为线性，避免单个异常点主导优化。其他可选核函数包括：

CauchyLoss：对重尾分布更鲁棒
SoftLOneLoss：平滑过渡到线性
TukeyLoss：完全忽略大残差

选择核函数需要根据数据特性。在我的一个视觉SLAM项目中，使用CauchyLoss有效降低了动态物体对位姿估计的影响。

3.2 求解器配置的艺术

Ceres的求解器配置对性能影响很大，主要选项包括：

ceres::Solver::Options options; options.linear_solver_type = ceres::DENSE_QR; // 小规模稠密问题 options.max_num_iterations = 100; // 最大迭代次数 options.function_tolerance = 1e-6; // 成本函数变化阈值 options.gradient_tolerance = 1e-10; // 梯度阈值 options.parameter_tolerance = 1e-8; // 参数变化阈值 options.minimizer_progress_to_stdout = true; // 打印进度 options.num_threads = 4; // 使用多线程

对于不同规模的问题，线性求解器的选择很关键：

DENSE_QR：参数少（<100）的稠密问题
DENSE_NORMAL_CHOLESKY：中等规模，比QR更快
SPARSE_NORMAL_CHOLESKY：大规模稀疏问题（如BA）
ITERATIVE_SCHUR：超大规模问题

我曾经优化一个包含3D重建的问题，将线性求解器从DENSE_QR改为SPARSE_NORMAL_CHOLESKY后，求解时间从15分钟降到了30秒。

4. 实战案例：Bundle Adjustment实现

Bundle Adjustment（光束法平差）是Ceres的典型应用场景，它通过优化相机位姿和3D点位置，最小化重投影误差。以下是简化的实现：

4.1 定义重投影误差

struct SnavelyReprojectionError { SnavelyReprojectionError(double observed_x, double observed_y) : observed_x(observed_x), observed_y(observed_y) {} template <typename T> bool operator()(const T* const camera, const T* const point, T* residuals) const { // camera[0,1,2]为旋转向量，[3,4,5]为平移 T p[3]; ceres::AngleAxisRotatePoint(camera, point, p); p[0] += camera[3]; p[1] += camera[4]; p[2] += camera[5]; // 投影到归一化平面 T xp = p[0] / p[2]; T yp = p[1] / p[2]; // 计算预测位置 T predicted_x = focal * xp; T predicted_y = focal * yp; // 残差 residuals[0] = predicted_x - observed_x; residuals[1] = predicted_y - observed_y; return true; } static ceres::CostFunction* Create(double observed_x, double observed_y) { return new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>( new SnavelyReprojectionError(observed_x, observed_y)); } double observed_x, observed_y; };

4.2 构建BA问题

void BundleAdjustment(std::vector<Feature>& features, std::vector<Camera>& cameras, std::vector<Point3D>& points) { ceres::Problem problem; for (auto& f : features) { ceres::CostFunction* cost_function = SnavelyReprojectionError::Create(f.x, f.y); problem.AddResidualBlock(cost_function, new ceres::HuberLoss(1.0), cameras[f.camera_idx].params, points[f.point_idx].coords); } ceres::Solver::Options options; options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY; options.max_num_iterations = 100; ceres::Solver::Summary summary; ceres::Solve(options, &problem, &summary); std::cout << summary.FullReport() << "\n"; }

在实际的SLAM系统中，还需要考虑：

参数块参数化（如旋转使用四元数）
固定某些参数（如第一帧相机位姿）
使用舒尔补加速求解

5. 性能优化与调试技巧

5.1 分析求解过程

Ceres的输出日志包含丰富信息：

iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter 0 4.512500e+01 0.00e+00 9.50e+00 0.00e+00 0.00e+00 1.00e+04 0 1 4.511598e-07 4.51e+01 9.50e-04 9.50e+00 1.00e+00 3.00e+04 1

关键指标解读：