利用ceres solver解大规模线性方程组

在工程应用中，最后要求解的线性方程组往往是原来的残差模型进行线性化后的误差方程。通常情况下，模型的线性化由人工完成，而求解误差方程则借助Eigen等矩阵运算库（参考1）来实现。现在，我们有了另一种选择，那就是ceres solver。ceres是Google公司用来解决非线性优化问题的开源库，主要是为了解决SFM问题中的光束法平差而设计的。与一般的矩阵运算库不同的是，我们只需要给ceres提供原始的残差模型就可以了，而无需自己求偏导数。ceres solver的安装和向导可以在下面给出的参考2的链接中找到，一个写的很好的中文用例可以在参考3的连接中找到。参考3给出的更多的是对参考2的翻译，最后也给出了一个空间后方交会的例子，但是未知数的数量只有6个，属于比较简单的例子。本文将结合实际应用，给出一个大规模最小二乘问题的求解实例。

现在有如上所示的一个能量函数，这个函数的第一项是观测值约束，第二项是未知数的平滑约束，这是一个工程领域十分常见的能量函数模型。我们用这个模型是要解决，利用一幅图像中稀疏的深度观测值来估计出一个密集的深度图（图1）。因此，对于一幅图像而言，每个像素对应一个未知数，比如对于1000*1000的一幅图像将会有100W个未知数。公式1是个非线性方程，但是求解这个问题却可以通过求偏导转换为一个线性方程组，具备确定解。但是由于未知数个数很多，所以利用QR分解等直接法不可能求解。好在ceres提供了丰富的选择可解决这个大规模的稀疏矩阵的求解问题（详见参考2）。

图 1. 左边为稀疏深度图，右边为密集深度图

下面，我们只阐述ceres解决该方程的具体方法，在此之前，要说的是公式1中第二项中的二阶偏导可以表示为离散形式：

这样的话，就可以直接利用ceres直接构建观测方程并求解了。

struct CostFunctor1 {CostFunctor1(double observed_depth):observed_depth(observed_depth){}template <typename T> bool operator()(const T* const depthmap, T* residual) const {residual[0] = depthmap[0] - observed_depth;return true;}static ceres::CostFunction* Create(const double observed_depth) {return (new ceres::AutoDiffCostFunction<CostFunctor1, 1, 1>(new CostFunctor1(observed_depth)));}double observed_depth;};struct CostFunctor2 {CostFunctor2(double observed_w):observed_w(observed_w){}template <typename T> bool operator()(const T* const depthleft, const T* const depth, const T* const depthright,T* residual) const {residual[0] = T(sqrt(observed_w*50.0))*(T(2.0)*depth[0] - depthleft[0] - depthright[0]);return true;}static ceres::CostFunction* Create(const double observed_w) {return (new ceres::AutoDiffCostFunction<CostFunctor2, 1, 1, 1, 1>(new CostFunctor2(observed_w)));}double observed_w;};

ceres::Problem problem;for (int i=0;i<height;i++){for (int j=0;j<width;j++){if (depthArray[i*width+j] != 0){ceres::CostFunction* cost_function = CostFunctor1::Create(depthArray[i*width+j]);problem.AddResidualBlock(cost_function,NULL,&densedepthmap[i*width+j]);}}}for (int i=0;i<height;i++){for (int j=1;j<width-1;j++){ceres::CostFunction* cost_function = CostFunctor2::Create(weightArrayX[i*width+j]);problem.AddResidualBlock(cost_function,NULL,&densedepthmap[i*width+j-1],&densedepthmap[i*width+j],&densedepthmap[i*width+j+1]);}}for (int i=1;i<height-1;i++){for (int j=0;j<width;j++){ceres::CostFunction* cost_function = CostFunctor2::Create(weightArrayY[i*width+j]);problem.AddResidualBlock(cost_function,NULL,&densedepthmap[(i-1)*width+j],&densedepthmap[i*width+j],&densedepthmap[(i+1)*width+j]);}}ceres::Solver::Options options;options.linear_solver_type = ceres::CGNR;options.num_linear_solver_threads = 6;options.function_tolerance = 0.1;options.minimizer_progress_to_stdout = true;ceres::Solver::Summary summary;ceres::Solve(options, &problem, &summary);std::cout << summary.BriefReport() << "\n";

参考1：（C++矩阵运算库推荐）http://blog.csdn.net/chenbang110/article/details/12304123

参考2：（Tutorial-Ceres Solver）http://ceres-solver.org/tutorial.html

参考3：（李民录csdn博客）http://blog.csdn.net/liminlu0314/article/details/15860677