C++学习笔记——Ceres模块（最小二乘求解器）

用于求解最小二乘问题的C++模块，做一下用法的笔记。
参考书籍：《视觉SLAM十四讲-从理论到实践》——高翔

1. 安装与卸载

1.1 依赖项

安装Ceres之前，首先安装依赖项：

sudo apt-get install liblapack-dev libsuitesparse-dev libcxsparese3 libgflags-dev libgoogle-glog-dev libgtest-dev

1.2 官方地址

git官方下载地址：https://github.com/ceres-solver/ceres-solver

1.3 安装

cmake 下载文件夹路径
make
sudo make install

1.4 卸载

sudo rm -r /usr/local/lib/cmake/Ceres
sudo rm -rf /usr/local/include/ceres /usr/local/lib/libceres.a

2. 样例

样例代码改编自《视觉SLAM十四讲-从理论到实践》——高翔
做了一些更详细的傻瓜式注释以便理解。

#include <iostream>
#include <opencv2/core/core.hpp>
#include <ceres/ceres.h>
#include <chrono>

using namespace std;

// 使用Ceres求解最小二乘法优化问题

// 代价函数的计算模型
struct CURVE_FITTING_COST {
    CURVE_FITTING_COST(double x, double y): _x(x), _y(y) {}

    // 残差的计算
    // 表示T在模版实例化时可以替换任意类型，不仅包括内置类型，也包括自定义类型
    template<typename T>
    // abc是3维的模型参数
    // 用operator关键词重载了()运算符，目的是将类打造为一个可以和函数一样调用的拟函数
    // 对实例化的类对象a可以直接调用a<double>()方法
    // 接收的参数将是Ceres传递过来的雅可比矩阵，用于自动求导
    bool operator() (const T *const abc, T *residual) const {
            // y-exp(ax^2+bx+c)
            // 使用ceres::exp计算e
            residual[0] = T(_y) - ceres::exp(abc[0] * T(_x) * T(_x) + abc[1] * T(_x) + abc[2]);
            return true;
        }

    // x, y数据
    const double _x, _y;
};

int main(int argc, char **argv) {
    // 配置真实参数
    double ar  = 1.0, br = 2.0, cr = 1.0;
    // 配置估计参数的初始值
    double ae = 2.0, be = -1.0, ce = 5.0;
    // 数据点
    int N = 100;
    // 噪声Sigma值
    double w_sigma = 1.0;
    // 噪声权重，与噪声强度呈反比，用于多噪声环境，在该例中并未使用
    double inv_sigma = 1.0 / w_sigma;
    // opencv随机数产生器
    cv::RNG rng;

    // 数据生成
    vector<double> x_data, y_data;
    for (int i=0; i<N; i++) {
        // x为0到0.99，步长0.01，间隔分布的100个值
        double x = i / 100.0;
        x_data.push_back(x);
        // 在y=ax^2+bx+c的基础上加入高斯噪声作为随机偏移量
        y_data.push_back(exp(ar *x *x + br *x + cr) + rng.gaussian(w_sigma * w_sigma));
    }

    // 一个参数块/残差块
    // 放入三个用于估计的参数
    double abc[3] = {ae, be, ce};

    // 构建最小二乘问题
    ceres::Problem problem;
    for (int i=0; i<N; i++) {
        // 向问题中添加误差项
        problem.AddResidualBlock(
            // 使用自动求导
            // 模板参数：误差类型、输出维度、输入维度，维数要与前面struct中一致
            // 因为要优化的是a, b, c三个量，且都是标量，所以是1x3，配置维度为1，3
            new ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>(
                new CURVE_FITTING_COST(x_data[i], y_data[i])
            ),
            // 核函数，不使用则配置为空指针
            nullptr,
            // 待估计参数
            abc
        );
    }

    // 配置求解器
    // 这里有很多配置项可以填，包括迭代次数，步长，使用的方法等等……
    // 可以参考Ceres的Options详解：https://blog.csdn.net/DumpDoctorWang/article/details/84890792
    ceres::Solver::Options options;
    // 配置最大迭代次数（默认50），但这只是最大值，程序可能提前停止
    options.max_num_iterations = 50;
    // 求解优化问题有两类方法：信任域（TRUST_REGION）和线性搜索（LINE_SEARCH，尚不支持边界约束）
    // options.minimizer_type默认TRUST_REGION

    // 增量方程如何求解
    options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY;
    // 设置输出到cout
    options.minimizer_progress_to_stdout = true;

    // 优化信息
    ceres::Solver::Summary summary;
    // 记录开始时间
    chrono::steady_clock::time_point t1 = chrono::steady_clock::now();

    // 开始优化
    ceres::Solve(options, &problem, &summary);

    // 记录结束时间
    chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
    // 计算所用的时间并输出到屏幕
    chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
    cout << "消耗时间：" << time_used.count() << "秒。" << endl;

    // 在屏幕上输出结果
    cout << summary.BriefReport() << endl;
    cout << "a, b, c 为：";
    for (auto a:abc) cout << a << " ";
    cout << endl;

    return 0;
}

3. 其余遇到过的功能

新建边缘化次序ceres::ParameterBlockOrdering：Ceres优化

自由度约束参数LocalParameterization：Ceres LocalParameterization 理解