C++数学库：Eigen & Ceres

Eigen3

线性代数库，由core模块和其他几个附加功能模块组成，只有头文件，安装非常容易。

Module	headerfile	Content
Core	Eigen/Core	数据数据类型和基本线性代数操作
Geometry	Eigen/Geometry	旋转放缩，欧拉数等几何变换相关支持
LU	Eigen/LU	LU分解以及求解器相关
Cholesky	Eigen/Cholesky	$LD{L}^T$分解和$Cholesky$分解以及求解器相关
Householder	Eigen/Householder	$Householder$变换
SVD	Eigen/SVD	SVD分解以及求解器相关
QR	Eigen/QR	QR分解以及求解器
Eigenvalues	Eigen/Eigenvalues	特征值，特征向量分解
Sparse	Eigen/Sparse	稀疏矩阵存储及相关基本线性代数
/	Eigen/Dense	Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
/	Eigen/Eigen	Includes Dense and Sparse header files (the whole Eigen library)

基础操作如下：

#include <eigen3/Eigen/Dense>
#include <iostream>

int main()
{
    Eigen::Matrix<double, 2, 2> matrix;
    matrix << 1, 2, 3, 4; // 初始化矩阵
    std::cout << "matrix = " << std::endl << matrix << std::endl;
    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_dynamic = matrix; // 矩阵赋值, Eigen::Dynamic表示矩阵的行和列都是动态的
    Eigen::MatrixXd matrix_xd = Eigen::MatrixXd::Random(2, 2); // 随机矩阵,Eigen中有很多typedef, MatrixXd是Matrix<double, Dynamic, Dynamic>的别名
    Eigen::Matrix2d matrix2d = Eigen::Matrix2d::Identity(); //单位矩阵
    Eigen::Matrix2d matrix2d2 = Eigen::Matrix2d::Zero(); //零矩阵
    Eigen::Matrix2d matrix2d3 = Eigen::Matrix2d::Ones(); //全1矩阵
    Eigen::Matrix2d matrix2d4 = Eigen::Matrix2d::Constant(2.0); //常数矩阵
    Eigen::Vector2d v(12, 12); //Vector2d就是Eigen::Matrix<double, 2, 1>
    Eigen::Matrix2d matrix_init_by_vector;
    matrix_init_by_vector << v, v/10; // 通过向量进行矩阵初始化
    matrix_init_by_vector.row(0) = v; // 通过向量赋值给矩阵的某一行
    matrix_init_by_vector.col(0) = v; // 通过向量赋值给矩阵的某一列
    matrix_init_by_vector(0, 0) = 1; // 通过下标赋值
    matrix_init_by_vector.row(0) << 1, 2; // 通过向量赋值给矩阵的某一行
    matrix_init_by_vector.row(0).tail(1) << 3; // 通过向量赋值给矩阵的某一行的某一部分
    matrix_init_by_vector.row(0).head(1) << 3; // 通过向量赋值给矩阵的某一行的某一部分
    Eigen::Matrix<double, 4, 4> matrix4x4;
    matrix4x4.setZero(); // 矩阵清零
    matrix4x4.setIdentity(); // 矩阵设为单位矩阵
    matrix4x4.setConstant(2.0); // 矩阵设为常数矩阵
    matrix4x4.setRandom(); // 矩阵设为随机矩阵
    matrix4x4.topLeftCorner(2, 2) = matrix; // 矩阵的左上角赋值
    
    Eigen::Vector2d result = matrix * v; // 矩阵和向量相乘

    Eigen::Matrix<double, 2, 2> matrix2;
    matrix2 = matrix.transpose(); // 矩阵转置

    Eigen::Matrix<double, 2, 2> matrix3;
    matrix3 = matrix.inverse(); // 矩阵求逆

    Eigen::Matrix<double, 2, 2> matrix_adjoint;
    matrix_adjoint = matrix.adjoint(); // 矩阵的伴随矩阵

    Eigen::Matrix<double, 2, 2> matrix4;
    matrix4 = matrix * matrix3 + matrix2 + matrix(1,1)*matrix - matrix3; // 矩阵加减乘除
    
    Eigen::Matrix<double, 2, 2> matrix5,matrix6;
    matrix5 = matrix.array() * matrix2.array(); // 矩阵对应元素相乘
    matrix6 = matrix.cwiseProduct(matrix2); // 矩阵对应元素相乘
    std::cout << "matrix" << std::endl << matrix << std::endl;
    std::cout << "matrix2" << std::endl << matrix2 << std::endl;
    std::cout << "matrix5:"<< std::endl << matrix5 << std::endl;
    std::cout << "matrix6" << std::endl << matrix6 << std::endl;

    Eigen::Vector2d v_cross;
    double v_dot;
    v_dot = matrix.row(0).dot(matrix.row(1)); // 矩阵行向量点乘
    v_cross = matrix.row(0).cross(matrix.row(1)); // 矩阵行向量叉乘
}

常见问题

• 如果使用std容器装Eigen数据类型，需要使用Eigen带的allocator，因为Eigen要求数据分配时必须按照至少16bytes对齐进行分配，所以配置分配器，Eigen自带了分配器，否则容易出bug。

对于支持C++17以后的编译器，不需要我们指定分配器，编译器会自动处理。

std::vector<Eigen::Matrix<double, 2, 2>,Eigen::aligned_allocator<Eigen::Matrix<double, 2, 2>>> matrix_vector;

• Eigen的数据类型不要通过值传递，通过引用传递，原因也是因为对齐问题。
• 编译的时候需要添加编译参数，-msse表示启动SSE指令，按照16bytes对齐，-mavx表示启动AVX指令，安装32bytes对齐，-march=native表示两种指令都开启。
  编译的时候必须所有的库统一编译参数，否则就会报错。

Eigen的对齐要求比较严格，可能引起的问题可以看官方文档这里

Ceres

非线性数学库，Google开发的，要求支持C++17，依赖Eigen和gflags库
Ceres主要是用于求解最小二乘问题，即：
$$\begin{array}{rl}\underset{x}{\min }& \frac{1}{2}\sum _i{\rho }_i\left({\Vert f_i\left(x_1,\cdots ,x_k\right)\Vert }^2\right)\\ \text{ s.t. }& l_j\le x_j\le u_j\end{array}$$
使用流程为：

• 定义问题
• 定义代价函数$f$
• 定义损失函数$\rho$

简单问题不定义损失函数，即定义损失函数为$\rho (x)=x$也没问题

具体代码如下：

#include <ceres/ceres.h>
#include <eigen3/Eigen/Dense>
#include <iostream>
#include <eigen3/Eigen/Core>

// 定义代价函数的结构体
// 代价函数主要是计算残差和求导，即残差的雅可比矩阵
// 如果要使用ceres的自动求导功能，需要将代价函数定义为结构体并重载()运算符
// 如果想要手动求导，可以继承ceres::CostFunction类，并重载Evaluate()函数
// ceres提供了自动求导，手动求导和数值求导三种方式
struct CostFunctor
{
    template <typename T>
    bool operator()(const T *const x, T *residual) const
    {
        residual[0] = T(2.0) - x[0]*x[0];
        return true;
    }
};

int main()
{
    ceres::Problem problem;
    ceres::CostFunction *cost_function = new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
    ceres::LossFunction *loss_function = new ceres::HuberLoss(0.1);
    double x = 2.0;
    problem.AddResidualBlock(cost_function, loss_function, &x); // 向问题中添加残差块和损失函数，以及待优化变量
    ceres::Solver::Options options; // 这里是设置求解器的参数
    options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY; // 选择一个线性求解器
    options.minimizer_progress_to_stdout = true; // 输出到cout

    ceres::Solver::Summary summary; // 求解过程信息
    ceres::Solve(options, &problem, &summary);
    std::cout << summary.BriefReport() << std::endl; // 输出求解结果
    std::cout << "x: " << 2 << " -> " << x << std::endl;
    return 0;
}

这篇blog讲的非常全面，可以参考。
Ceres中的自动微分功能有三种实现方式：手工求导，自动求导，数值方法求导

这里的自动求导和机器学习框架中的对tensor的自动求导一样。都是在重载了基础计算操作符，在计算的同时将导数也进行计算。

Ceres的自动求导是定义了一个Jet类型：Jet<double, number> x，Jet由数值a和扰动向量v组成，至于为什么要用一个向量，感觉可能是因为这个和机器学习中的导数存放的地方不太一样，例如对于$y=x_1+x_2$，Jet的做法是在y的v里存储$[1,1]$，两个x里分别存储$[1,0]$和$[0,1]$，而机器学习中只在计算出最后结果后进行一次backward，每个变量存储的是最终结果对自己的偏导数。

Ref

Eigen官方文档
SLAM Ceres优化库讲解
Ceres 官方文档