基于线性代数的C ++模板库

3D视觉工坊 2023-06-25 354

描述

00 Eigen简介

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

Eigen库的优点：

支持整数、浮点数、复数，使用模板编程，可以为特殊的数据结构提供矩阵操作。

OpenCV自带到Eigen的接口。

支持逐元素、分块、和整体的矩阵操作。

支持使用Intel MKL加速部分功能。

支持多线程，对稀疏矩阵支持良好。

支持常用几何运算，包括旋转矩阵、四元数、矩阵变换、角轴等等。

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

01 数据类型

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 Matrix4f;

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

typedef Matrix< float, 3, 1 > Vector3f;

行向量RowVector：

typedef Matrix RowVector2i;

静态-动态-矩阵

静态矩阵：矩阵是静态的，即编译时候就知道运行结果，例如Matrix3d：表示元素类型为double大小为3*3的矩阵变量，其大小在编译时就知道。

动态矩阵：有时候运行完之后，才可以知道，这里使用MatrixXd：表示任意大小的元素类型为double的矩阵变量，其大小只有在运行被赋值之后才能知道；

数据类型

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

d：表示double类型

f：表示float类型

i：表示整数

c：表示复数；

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

02 新建矩阵

矩阵构造

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

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();      // 随机矩阵

03 矩阵索引

当前矩阵的行数、列数、大小可以通过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(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

04 数学运算

4.1 加减法

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

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

4.2 乘除法

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

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中，向量的叉乘只支持三维的向量，这是因为叉乘通常用于计算方向、夹角等，它的计算规则如下：

4.3 特征运算

// 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)

05 通用数组

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

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

m1.array().atan()；

常见数据类型

Array                   ArrayXf
Array                         Array3f
Array            ArrayXXd
Array                        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函数

06 更多操作

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

打开APP阅读更多精彩内容