Ceres求解库用法简析

ceres是C++常用的最小二乘问题求解库，最近研究colmap代码，顺便学习一下

下面举两个简单例子来说明：

问题一：一维优化，让x收敛到0

// 定义代价函数，结构体、类均可，只需重载operator()
struct CostFunctor {
      template <typename T>
      bool operator()(const T* const x, T* residual) const {
            residual[0] = T(10.0) - x[0];
            return true;
      }
};
int main(int argc, char** argv){
      // 初始值
      double x = 0.0;

      ceres::Problem problem;
      // 使用AutoDiff包装代价函数
      // 残差维度(即residual数组大小)为1
      // 参数块维度(即传给AddResidualBlock最后一项的x，也就是待优化项的大小)为1
      // nullptr处为核函数，决定如何处理离群点，传入空表示信任所有观测数据
      problem.AddResidualBlock(
            new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor),
            nullptr, &x
      );
      // 配置求解器
      // 这里使用DENSE_QR，对于小规模问题表现良好
      // colmap中使用舒尔补类型SCHUR应对视觉测量的稀疏特点，并会根据照片规模决定求解器：
      // 50张以下，使用DENSE_SCHUR
      // 1000张以下且支持稀疏计算，则使用SPARSE_SCHUR
      // 其他情况使用ITERATIVE_SCHUR，并使用预处理器SCHUR_JACOBI缩放问题
      // 以上详见bundle_adjustment.cc
      ceres::Solver::Options options;
      options.linear_solver_type = ceres::DENSE_QR;
      options.minimizer_progress_to_stdout = true; // 打印迭代情况

      ceres::Solver::Summary summary;
      ceres::Solver(options, &problem, &summary);

      std::cout << summary.BriefReport() << std::endl;
      std::cout << "Final x: " << x << std::endl;
      return 0;
}

问题二：曲线拟合y=e^{ax+b}，优化a、b

struct CostFunctor {
      CostFunctor(double x, double y) : x_(x), y_(y) {}
      template <typename T>
      bool operator()(const T* const ab, T* residual) const {
            residual[0] = y_ - exp(ab[0] * x_ + ab[1]);
            return true;
      }
      double x_, y_;
};
int main(int argc, char** argv){
      double ab[2] = {0.0, 0.0};
      ceres::Problem problem;
      // x_data和y_data是事先生成好的一组观测值
      // 对于每个观测值，增加一个约束（让ab对应的曲线去满足每一个点）
      for (int 1 = 0; i < x_data.size(); ++i){
            problem.AddResidualBlock(
                  new ceres::AutoDiffCostFunction<CostFunctor, 1, 2>(
                        new CostFunctor(x_data[i], y_data[i])
                  ),
                  nullptr, ab
            );
      }
      ceres::Solver::Options options;
      options.linear_solver_type = ceres::DENSE_QR;
      options.minimizer_progress_to_stdout = true;
      ceres::Solver::Summary summary;
      ceres::Solver(options, &problem, &summary);
      std::cout << summary.BriefReport() << std::endl;
      std::cout << "Final a, b: " << ab[0] << ", " << ab[1] << std::endl;
      return 0;
}

注意理解残差维度不等于参数维度，举例：

1. 输入x、y，输出z，残差维度=1，参数维度=2；
2. 输入x，输出坐标 (u, v)，残差维度=2，参数维度=1；

还有其他一些使用技巧，例如核函数、给problem设置固定参数和上下界，这里就不赘述了