一文学会Eigen库

taotao_2016 2023-07-03 发布于辽宁

展开全文

Eigen：基于线性代数的C ++模板库，主要用于矩阵，向量，数值求解器和相关算法。SLAM中常用的Ceres、G2O等项目均是基于Eigen库。

文章图片1

Eigen库的优点：

支持整数、浮点数、复数，使用模板编程，可以为特殊的数据结构提供矩阵操作。
OpenCV自带到Eigen的接口。
支持逐元素、分块、和整体的矩阵操作。
支持使用Intel MKL加速部分功能。
支持多线程，对稀疏矩阵支持良好。
支持常用几何运算，包括旋转矩阵、四元数、矩阵变换、角轴等等。

即使不做SLAM，在3D视觉中，当处理大量数学运算时，我们也会用到Eigen库，它帮我们优化了性能。在安装完成Eigen库后，开始接下来的学习。

Eigen库的核心类是 Matrix，由6个参数构成：

Matrix< typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime, int Options = 0, // 默认（无需更改） int MaxRowsAtCompileTime = RowsAtCompileTime, // 默认（最大行数，提前知道极限） int MaxColsAtCompileTime = ColsAtCompileTime // 默认（最大列数，提前知道极限）>

其中：

前三个参数：需要我们指定
后三个参数：默认即可，无需指定

因为经常需要实例化一些方阵、向量，因此Eigen库也提供了很多直接使用的模板（利用C++的关键字：typedef），例如 Matrix4f 是的float型矩阵：

typedef Matrix<float, 4, 4> Matrix4f;

还有例如列向量：Vector3f ，其本质也是 Matrix 类：

typedef Matrix< float, 3, 1 > Vector3f;

行向量RowVector：

typedef Matrix<int, 1, 2> RowVector2i;

静态-动态-矩阵

静态矩阵：矩阵是静态的，即编译时候就知道运行结果，例如Matrix3d：表示元素类型为double大小为3*3的矩阵变量，其大小在编译时就知道。
动态矩阵：有时候运行完之后，才可以知道，这里使用MatrixXd：表示任意大小的元素类型为double的矩阵变量，其大小只有在运行被赋值之后才能知道；

数据类型

Eigen中的矩阵类型一般都是用类似MatrixNX来表示，可以根据该名字来判断其大小（2，3，4，或X，意思Dynamic）和数据类型，比如：

d：表示double类型
f：表示float类型
i：表示整数
c：表示复数；

举例：Matrix2f，表示的是一个维的，其每个元素都是float类型。

矩阵构造

默认构造，分配了大小和内存空间，但没有初始化矩阵元素（里面的数值是随机的，不能使用）：

Matrix3f a; // 3*3的元素，其中还有一个float[9]数组，其中的元素没有初始化；MatrixXf b; // 动态大小的矩阵，目前的大小是0*0，它的元素数组完全没有分配。

对于动态数组，你也可以直接分配大小（失去作用了），同样没有初始化矩阵元素：

MatrixXf a(10, 15);// 10x15动态矩阵，数组内存已经分配，但是没有初始化；VectorXf b(30);   // 大小为30的向量，数组内存已经分配，但是元素没有初始化。

或者更通用的：

Matrix< float, 3, 1 > Vector3f_def;

矩阵初始化

在构造完后，我们需要对元素进行初始化，常用的是直接赋值：

Eigen::Matrix3f m; m << 1, 2, 3,     4, 5, 6,     7, 8, 9;

它是逐行写入的，这只适用于较小的矩阵：

Eigen::MatrixXd m(3,3);m <<1,2,3, 4,5,6, 7,8,9;

对于向量，还可以在构造的时候初始化：

Vector3d v(1, 2, 3);Vector3d w(1, 0, 0);

还有一些特殊函数，函数：

MatrixXf::Zero(3,4); // 将矩阵3行4列初始化为0 MatrixXf::Ones(3,3); // 将矩阵3行3列初始化为1 Vector3f::Ones(); // 将3行的纵向量初始化为1 MatrixXi::Identity(3,3); // 单位矩阵 Matrix3d::Random(); // 随机矩阵

当前矩阵的行数、列数、大小可以通过rows()、cols()和size()来获取。遍历Eigen矩阵时最好通过rows和cols来限制访问范围，索引的方法如下：

1、矩阵访问按照先行索引、后列索引方式进行，索引下标从0开始（与Matlab不同）；

2、矩阵元素的访问可以通过**”( )”操作符完成。例如m(2, 3)**，矩阵m的第2行第3列元素；

3、针对向量还提供”**[ ]”操作符，注意矩阵则不可**如此使用。

resize：不同于matlab、Python，对于动态矩阵虽然可以通过resize()函数来动态修改矩阵的大小，但是需要说明的是，在Eigen中：
不能用：固定大小的矩阵是不能使用resize()来修改矩阵的大小；
数据会变：resize()函数会析构掉原来的数据，变为0.，因此最好使用：conservativeResize()函数
大小修改：使用”=”操作符操作动态矩阵时，如果左右两边的矩阵大小不等，则左边的动态矩阵的大小会被修改为右边的大小。

利用block()函数，可以从Matrix中取出一个小矩阵来进行处理，使用的语法为：

matrix.block<p,q>(i,j);

例如：

Eigen::MatrixXf m(4, 4);m << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16;cout << 'Block in the middle' << endl;cout << m.block<2, 2>(1, 1) << endl << endl;for (int i = 1; i <= 3; ++i){ cout << 'Block of size ' << i << 'x' << i << endl; cout << m.block(0, 0, i, i) << endl << endl;}// Output is:// Block in the middle// 6 7// 10 11// Block of size 1x1// 1// Block of size 2x2// 1 2// 5 6// Block of size 3x3// 1 2 3// 5 6 7// 9 10 11

单独的列和行是块的特殊情况。Eigen提供了可以轻松解决它们的方法：.col（）和.row（）：

Eigen::MatrixXi m(2, 2);m << 1, 2, 3, 4;cout << m.col(0) << endl;//  1 3

Eigen帮我们重载了，直接运算：

Vector3d v(1, 2, 3);Vector3d w(1, 0, 0);cout << v + w << endl;

除法：通常我们是除以标量。对于矩阵除法，我们是求它的逆，再转换为矩阵乘法。因此较为简单：

Vector3d v(1, 2, 3);Vector3d r = v / 3;cout << r << endl;

矩阵乘法：*

乘法，标量非常简单：

cout << v * 2 << endl;v *= 2; // 原地操作

Matrix2d mat;mat << 1, 2,    3, 4;Vector2d u(-1, 1), v(2, 0);// 矩阵乘法 乘以矩阵std::cout << 'Here is mat*mat:'          << mat * mat << std::endl;// 矩阵乘法 乘以向量std::cout << 'Here is mat*u:'          << mat * u << std::endl;// 转置之后，再矩阵乘法std::cout << 'Here is u^T*mat:'          << u.transpose() * mat << std::endl;// 转置之后，向量的矩阵乘法std::cout << 'Here is u^T*v:'          << u.transpose() * v << std::endl;std::cout << 'Here is u*v^T:'          << u * v.transpose() << std::endl;// 矩阵乘法std::cout << 'Let's multiply mat by itself' << std::endl;mat = mat * mat;std::cout << 'Now mat is mat:'          << mat << std::endl;//Output is:// Here is mat*mat://  7 10// 15 22// Here is mat*u:// 1// 1// Here is u^T*mat:// 2 2// Here is u^T*v:// -2// Here is u*v^T:// -2 -0//  2  0// Let's multiply mat by itself// Now mat is mat://  7 10// 15 22

补充：转置

向量、矩阵的乘法，因为需要size一致，因此需要用到转置：

MatrixXcf a = MatrixXcf::Random(2, 2); //MatrixXcf 为复数矩阵cout << 'Here is the matrix a' << a << endl;// 矩阵转置cout << 'Here is the matrix a^T' << a.transpose() << endl;// 共轭矩阵cout << 'Here is the conjugate of a' << a.conjugate() << endl;// 共轭转置矩阵cout << 'Here is the matrix a^*' << a.adjoint() << endl;

需要说明的是，在Eigen中，对于自身的操作，都有专门的函数，例如对自身的转置：

 a.transposeInPlace(); // 直接在a上操作

点乘和叉乘

Vector3d v(1, 2, 3);Vector3d w(0, 1, 2);// 点乘cout << 'Dot product: ' << v.dot(w) << endl;// 叉乘cout << 'Cross product:' << v.cross(w) << endl;// 点成结果Dot product: 8 // 1 * 0 + 2 * 1 + 3 * 2=8 Cross product: 1 // 2 * 2 - 1 * 3 = 1-2 // 3 * 0 - 1 * 2 = -2 1 // 1 * 1 - 0 * 2 = 1

在Eigen中，向量的叉乘只支持三维的向量，这是因为叉乘通常用于计算方向、夹角等，它的计算规则如下：

// Eigen also provides some reduction operations to reduce a given matrix or vector to a single value// such as the sum (computed by sum()), product (prod()), or the maximum (maxCoeff()) and minimum (minCoeff()) of all its coefficients.Eigen::Matrix2d mat;mat << 1, 2,       3, 4;//元素和，元素乘积，元素均值，最小系数，最大系数，踪cout << 'Here is mat.sum():       ' << mat.sum() << endl;cout << 'Here is mat.prod():      ' << mat.prod() << endl;cout << 'Here is mat.mean():      ' << mat.mean() << endl;cout << 'Here is mat.minCoeff():  ' << mat.minCoeff() << endl;cout << 'Here is mat.maxCoeff():  ' << mat.maxCoeff() << endl;cout << 'Here is mat.trace():     ' << mat.trace() << endl;// 可以返回元素位置Matrix3f m = Matrix3f::Random();std::ptrdiff_t i, j;  // std::ptrdiff_t 是二个指针相减结果的有符号整数类型float minOfM = m.minCoeff(&i, &j);cout << 'Here is the matrix m:'     << m << endl;cout << 'Its minimum coefficient (' << minOfM     << ') is at position (' << i << ',' << j << ')';RowVector4i v = RowVector4i::Random();int maxOfV = v.maxCoeff(&i);cout << 'Here is the vector v: ' << v << endl;cout << 'Its maximum coefficient (' << maxOfV     << ') is at position ' << i << endl;// Output is:// Here is mat.sum():       10// Here is mat.prod():      24// Here is mat.mean():      2.5// Here is mat.minCoeff():  1// Here is mat.maxCoeff():  4// Here is mat.trace():     5// Here is the matrix m://  -0.444451   0.257742   0.904459//    0.10794  -0.270431    0.83239// -0.0452059  0.0268018   0.271423// Its minimum coefficient (-0.444451) is at position (0,0)

Array类提供了通用数组。此外，Array类提供了一种执行逐系数运算的简便方法，该运算可能没有线性代数含义，例如将常数添加到数组中的每个系数或按系数乘两个数组。

注：Eigen计算三角函数等，Matrix并不支持，需要通过.array() 转换到Array类，再计算！
m1.array().atan()；

常见数据类型

Array<float,Dynamic,1>                   ArrayXfArray<float,3,1>                         Array3fArray<double,Dynamic,Dynamic>            ArrayXXdArray<double,3,3>                        Array

常见操作：

// 逐元素操作Vectorized operations on each element independently // Eigen // Matlab //注释 R = P.cwiseProduct(Q); // R = P .* Q //逐元素乘法 R = P.array() * s.array(); // R = P .* s //逐元素乘法（s为标量） R = P.cwiseQuotient(Q); // R = P ./ Q //逐元素除法 R = P.array() / Q.array(); // R = P ./ Q //逐元素除法 R = P.array() + s.array(); // R = P + s //逐元素加法（s为标量） R = P.array() - s.array(); // R = P - s //逐元素减法（s为标量） R.array() += s; // R = R + s //逐元素加法（s为标量） R.array() -= s; // R = R - s //逐元素减法（s为标量） R.array() < Q.array(); // R < Q //逐元素比较运算 R.array() <= Q.array(); // R <= Q //逐元素比较运算 R.cwiseInverse(); // 1 ./ P //逐元素取倒数 R.array().inverse(); // 1 ./ P //逐元素取倒数 R.array().sin() // sin(P) //逐元素计算正弦函数 R.array().cos() // cos(P) //逐元素计算余弦函数 R.array().pow(s) // P .^ s //逐元素计算幂函数 R.array().square() // P .^ 2 //逐元素计算平方 R.array().cube() // P .^ 3 //逐元素计算立方 R.cwiseSqrt() // sqrt(P) //逐元素计算平方根 R.array().sqrt() // sqrt(P) //逐元素计算平方根 R.array().exp() // exp(P) //逐元素计算指数函数 R.array().log() // log(P) //逐元素计算对数函数 R.cwiseMax(P) // max(R, P) //逐元素计算R和P的最大值 R.array().max(P.array()) // max(R, P) //逐元素计算R和P的最大值 R.cwiseMin(P) // min(R, P) //逐元素计算R和P的最小值 R.array().min(P.array()) // min(R, P) //逐元素计算R和P的最小值 R.cwiseAbs(P) // abs(P) //逐元素计算R和P的绝对值 R.array().abs() // abs(P) //逐元素计算绝对值 R.cwiseAbs2() // abs(P.^2) //逐元素计算平方 R.array().abs2() // abs(P.^2) //逐元素计算平方 (R.array() < s).select(P,Q); // (R < s ? P : Q) //根据R的元素值是否小于s，选择P和Q的对应元素 R = (Q.array()==0).select(P,A) // R(Q==0) = P(Q==0) R(Q!=0) = P(Q!=0) //根据Q中元素等于零的位置选择P中元素 R = P.unaryExpr(ptr_fun(func)) // R = arrayfun(func, P) // 对P中的每个元素应用func函数

对于Eigen，它适合一个简单的数值计算库，并没有什么实用技巧。其实大多数时候，你只需要利用Google和百度去查询你需要的操作即可！对于更多的操作，可以参考：Eigen 常用函数查询，对比MatLab操作。