Eigen库的入门使用

2,842 阅读8分钟

前言

Eigen是一个基于线性代数的C++模板库,主要用于矩阵、向量、数值求解和相关算法,不论是做算法还是C++开发,使用好Eigen库会让你事半功倍。

正文

由于本人是非专业的算法工程师,所以这里也只是做一些简单入门介绍和使用。

类型定义

Eigen库中的核心类就是Matrix,它表示矩阵,是一个模板类,定义如下:

template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
class Matrix { }

其中前3个模板参数需要我们指定,后面3个参数使用默认即可。

  • 第一个参数typename _Scalar,表示该参数是一个数据类型,而不是像_Rows_Cols这样的变量。那么现在比如我想定义一个4x4float矩阵,就可以如下定义:
Matrix<float,4,4> matrix44f;

这样定义虽然简单明了,但是有点长,所以在Eigen库中提供了很多开源直接使用的模板类,比如:

  • Matrix4f,表示4x4float矩阵;
  • Vector3f,表示列向量,长度3的float向量;
  • RowVector2i,表示行向量,长度2的int向量;

所以对于一些常见矩阵可以使用内置的模板类。

除了定义已知大小的矩阵,Eigen库还可以定义动态矩阵。那什么是动态矩阵和静态矩阵呢?

静态矩阵就是编译时候就知道大小的矩阵,比如Matrix3d,在编译时我们就知道这是一个3x3double类型矩阵。

与之对应的是动态矩阵,这种矩阵需要在运行后才知道大小,比如需要从vector<vector<double>>数据类型中构建矩阵,只有当代码运行时才可以知道构建的矩阵的行列个数。这时可以使用MatrixXd来表示任意大小的元素类型为double的矩阵变量,其中X就表示未知大小,即Dynamic动态的意思。

对于矩阵类型后面的后缀表示元素类型,Eigen库定义了4种元素类型:

  • d:表示double类型;
  • f:表示float类型;
  • i:表示整数;
  • c:表示复数;

比如Matrix2c,表示一个2x2元素类型为复数的矩阵。

新建矩阵

对于默认构造函数的矩阵,分配了大小和内存空间,但是没有初始化矩阵元素。而对于基本数据类型,默认初始化的话,其值是随机的,不能使用。比如下面代码:

Eigen::Matrix2f m;
std::cout << "m:" << std::endl << m << std::endl;
Eigen::Vector4i n;
std::cout << "n:" << std::endl << n << std::endl;

//运行结果:
m:
5.88604e-039            0
9.18341e-041  2.8026e-045
n:
-517829874
        -2
1952350234
   4199631

在Eigen中重载了std::cout<<运算符,所以可以直接使用std::cout打印矩阵结果。所以定义完矩阵后,我们必须要对进行初始化。对于矩阵的初始化,有如下几种方式:

  1. 直接赋值,可以使用<<进行赋值。
Eigen::Matrix2f m;  //2x2的float矩阵
m << 1.9,2.884,     //进行初始化,期望是4个值
     4.33,5.898;
std::cout << "m:" << std::endl << m << std::endl;
​
Eigen::Vector4i n;  //长度为4的列向量
n << 3,4,2,6;
std::cout << "n:" << std::endl << n << std::endl;
​
//运行结果:
m:
1.9 2.884
4.33 5.898
n:
3
4
2
6

这种方式只适合于矩阵的大小固定的情况,这种初始化方式很像使用列表初始化的数组,只有知道矩阵的大小,其构造函数才能从左向右、从上向下来进行初始化。比如可以定义动态大小的矩阵MatrixXd,这时就无法使用<<进行直接赋值,比如:

Eigen::MatrixXi xi;     //动态大小的int矩阵
xi << 1,2,3,4,5,6;
std::cout << "xi:" << std::endl << xi << std::endl;

因为这种定义的变量xi,并不知道其行列数目,可能是3x22x3等等,所以用6个int去赋值,构造函数无法进行构造。

  1. 对于向量,可以在构造的时候进行初始化
Eigen::RowVector3d k(1,2,3);    //长度为3,类型为double的行向量
std::cout << "k:" << std::endl << k << std::endl;
  1. 一些特殊函数,比如全部初始化为0,初始化为1,随机数等。
Eigen::MatrixXi zero = MatrixXi::Zero(3,4);    //全部初始化为0的矩阵
std::cout << "------ zero ------" << std::endl << zero << std::endl;
Eigen::MatrixXi one = MatrixXi::Ones(4,4);      //全部初始化为1的矩阵
std::cout << "------ one ------" << std::endl << one << std::endl;
Eigen::MatrixXi i = MatrixXi::Identity(3,3);    //单位矩阵
std::cout << "------ i ------" << std::endl << i << std::endl;
Eigen::Matrix4f random = Matrix4f::Random();    //随机矩阵
std::cout << "------ random ------" << std::endl << random << std::endl;

运行结果:

------ zero ------
0 0 0 0
0 0 0 0
0 0 0 0
------ one ------
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
------ i ------
1 0 0
0 1 0
0 0 1
------ random ------
 -0.997497   0.170019    0.64568   0.421003
  0.127171 -0.0402539    0.49321  0.0270699
 -0.613392  -0.299417  -0.651784   -0.39201
  0.617481   0.791925   0.717887  -0.970031
​
  1. 使用Eigen::Map已有的内存块数据映射到矩阵,必须要指定需要映射的矩阵类型,是MatrixXdVectorXd等等,在模板参数中指明。

看一下Eigen::Map函数的定义:

template< e> class Map
  : public MapBase<Map<PlainObjectType, MapOptions, StrideType> >
{}
  • PlanObjectType:映射数据的等效矩阵类型。
  • MapOptions:指定指针对齐方式,以字节为单位,默认为Unaligned,即未对齐。
  • StrideType:可选择指定步幅。默认情况下,Map采用内存布局是一个普通的、连续的数组,可以通过指定步幅来定制化。

正常我们只需指明第一个参数即可,使用Map可以在无任何其他开销的情况下,让Eigen使用非Eigen数据结构来对其进行初始化,包括多种数据结构,比如原生数组,比如double[]float[]等,以及STL容器,比如std::vectorstd::array等。

std::vector<double> vec = {1,2,3,4,5,6};    //std::vector数据类型
Eigen::Map<Eigen::VectorXd> vd(vec.data(),6);   //构造成一个长度为6的列向量
std::cout << "------ vd ------" << std::endl << vd << std::endl;
Eigen::Map<Eigen::MatrixXd> xd(vec.data(),3,2); //构造成一个3x2的矩阵
std::cout << "------ xd ------" << std::endl << xd << std::endl;
    
//运行结果:
------ vd ------
1
2
3
4
5
6
------ xd ------
1 4
2 5
3 6

在上面代码中,我们使用vec.data()来获取数vector中的数组数据,同时必须确定矩阵或者向量的大小。

  1. 对于矩阵来说,还有一种非常常见的初始化方式,就是一行一行或者一列一列进行赋值。比如代码:
Eigen::MatrixXd mat(3,2); //创建一个3x2的double类型矩阵
mat.col(0) << 1,2,3;    //对第一列进行初始化
mat.col(1) << 4,5,6;    //对第二列进行初始化

这种方式在构建特殊矩阵时经常用到。

矩阵索引

当前矩阵的行数、列数和大小可以分别通过rows()cols()size()方法来获取,在遍历Eigen矩阵时最好获取行数、列数来限制范围。

同时索引下标从0开始,矩阵元素的访问可以通过matrix(i,j)来访问第i行第j列的元素,对于向量还可以使用类似数组取值的vector[i]来访问第i个元素。比如代码:

std::vector<double> vec = {1,2,3,4,5,6};
Eigen::Map<Eigen::VectorXd> vd(vec.data(),6);   //创建一个长度为6的列向量
std::cout << "------ vd ------" << std::endl << vd << std::endl;
Eigen::Map<Eigen::MatrixXd> xd(vec.data(),3,2); //创建一个3x2的矩阵
std::cout << "------ xd ------" << std::endl << xd << std::endl;
for (int i = 0; i < xd.rows(); ++i) {
        for (int j = 0; j < xd.cols(); ++j) {
            std::cout << "xd.(" << i << "," << j << ")= " << xd(i,j) << std::endl;
        }
}
​
//运行结果:
------ xd ------
1 4
2 5
3 6
xd.(0,0)= 1
xd.(0,1)= 4
xd.(1,0)= 2
xd.(1,1)= 5
xd.(2,0)= 3
xd.(2,1)= 6

还可以利用block()函数来从Matrix中取出一个小矩阵来进行处理,语法是matrix:block<p,q>(i,j),表示从原矩阵的(i,j)位置开始,截取出pxq大小的子矩阵,比如测试代码:

Eigen::MatrixXi randomI = MatrixXi::Random(6,6);    //随机值的6x6矩阵
std::cout << "------ randomI ------" << std::endl << randomI << std::endl;
Eigen::MatrixXi smallI = randomI.block<3,3>(1,1);   //从(1,1)处截取一个3x3的矩阵
std::cout << "------ smallI ------" << std::endl << smallI << std::endl;
​
//运行结果:
------ randomI ------
-13389 -12482   3334 -14515  12319  -1243
 -4442 -16231   3511   3528   7427  -8673
-11557 -16092 -10937   9283  14938  11869
-10948  -4002   5342   9915  13949  -9516
 16007   1037  -1613    651   1289   9163
 -1780   2332  -4846  -6490 -11720  11260
------ smallI ------
-16231   3511   3528
-16092 -10937   9283
 -4002   5342   9915    

数学运算

Eigen重载了常用加减乘运算符,而对于矩阵的除法,我们一般是求逆然后再转换为乘法,使用起来比较简单,就不过多介绍了。

还有几个常用的方法,必须要介绍一下,因为在项目中非常常见。

  • 矩阵转置。调用transpose()方法即可求出一个矩阵的转置矩阵。
  • 矩阵求逆。调用inverse()方法即可求出一个矩阵的逆矩阵。
  • 共轭矩阵。调用conjugate()方法可以求出共轭矩阵。
  • 求平均值。调用mean()方法求平均值。
  • 逐元素相乘。调用cwiseProduct()函数进行2个矩阵逐元素相乘。
  • 最大/小值。调用maxCoeff()minCoeff()直接求出矩阵的最大和最小值。

总结

当然在实际使用中,使用Eigen的最大场景是求解线性方程组,这个我们后面文章继续探讨。

更多算法相关总结移步专栏: juejin.cn/column/7415…