Eigen深入剖析Matrix（进阶一）

小怪物爱吹牛

已于 2024-04-18 17:27:31 修改

阅读量4.3k

点赞数 7

分类专栏：科学计算文章标签：机器学习 python 算法

于 2023-07-05 21:41:54 首次发布

本文链接：https://blog.csdn.net/naibula/article/details/131563961

版权

科学计算专栏收录该内容

13 篇文章 44 订阅

订阅专栏

Eigen是以线性代数和矩阵理论为基础的数值计算库。在Eigen中所有的矩阵和向量都是Matrix模板类的对象，向量只是一种特殊的矩阵，默认向量是列向量。Array类的很多操作与Matrix类似，所以本章中介绍的部分操作也适用于Array。

1 模板类Matrix

矩阵类Matrix是Eigen中最基本的数据结构。Matrix类总共有六个模板参数首先介绍前三个参数，剩下的三个参数有其默认值。三个必须的参数如下：

Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

Scalar是系数的数据类型，如想要构造一个单精度的浮点类型矩阵，可以选择float。常见的Scalar类型有int、float、double、complex、complex（到底支持哪些类型还有待确认）。
RowsAtCompileTime和ColsAtCompileTime分别表示行数和列数，要求这两个参数在编译时已知。
在Eigen中提供宏定义来便捷的访问一些常用的类型，如：

typedef Matrix<float, 4, 4> Matrix4f;
typedef Matrix<float, 3, 1> Vector3f;
typedef Matrix<int, 1, 2> RowVector2i;

当然，Eigen并不局限于那些矩阵维数在编译时已知的情形。RowsAtCompileTime和ColsAtCompileTime可以取一个特殊的值Dynamic，这表示矩阵的维度在编译时是未知的，必须作为一个运行时变量来处理。Dynamic是一个const int类型的常量，其取值为-1。在Eigen术语中，Dynamic称为动态大小(dynamic size)，而运行时已知的大小称为固定大小(fixed size)。当矩阵的元素较少时，建议使用固定大小的形式，这样可以避免一些运行时错误，在运行效率上也有优势。

// 创建一个双精度的动态矩阵
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
// 创建一个整型列向量
typedef Matrix<int, Dynamic, 1> VectorXi;

Matrix的另外三个模板参数是可以选择的，完整的参数如下：

Matrix<typename Scalar,
       int RowsAtCompileTime,
       int ColsAtCompileTime,
       int Options = 0,
       int MaxRowsAtCompileTime = RowsAtCompileTime,
       int MaxColsAtCompileTime = ColsAtCompileTime>

Options是一个位域，它的值只能取RowMajor、ColMajor，分别表示行优先存储、列优先存储。默认情况下是使用列优先存储，一般地Eigen库在列优先存储的情况下效率要高。Eigen的存储形式的指定，只是影响矩阵在内存的形式，不会影响矩阵的访问习惯。
MaxRowsAtCompileTime和MaxColsAtCompileTime用于指定矩阵的维数的最大值。尽管有时候不知道矩阵的确切大小，但在编译时已经知道了矩阵维数的上界，那么可以指定这两个参数来避免动态内存分配。
这里强调一下，Eigen默认是列优先的，比如定义一个Vector3d向量，对应的维度是3行1列。PlainObjectBase既从MatrixBase派生，又从ArrayBase派生，但不是多继承。原因是PlainObjectBase是类模板，在生成不同的具体类时选择的基类不一样，实现这种选择的是使用traits技术。
类的层次结构示意图

1.1 插件技术

Eigen将很多类似的操作提取为插件，从而减少重复代码。
DenseBase插件

1.2 traits技术

PlainObjectBase会根据模板参数传入类型的差别选择不同的基类。如果传入模板参数Derived是Marix类型，那么将MatrixBase作为基类；如果传入模板参数Derived是Array类型，那么将ArrayBase作为基类。这个类型推导由dense_xpr_base模板类完成。
假定PlainObjectBase<Matrix<…>>接收到一个Matrix类型的具体类型，由模板推导机制，主模板首先计算默认类型traits<Matrix<…>>::XprKind。traits实现的功能是找到Matrix类型对应的基类，那么必须要定义一个以Matrix为特化参数的traits，Eigen将Matrix对应的traits也定义在Matrix.h中。源码中traits<Matrix<…>>将XprKind设置为为MatrixXpr的别名，即主模板默认类型推导结果是XprKind=MatrixXpr。进而进行模板特化，匹配到dense_xpr_base<Derived, MatrixXpr>，最终计算出dense_xpr_base::type=MatrixBase。
这里可以看出，主模板可以引入默认类型就像函数形参可以设置默认值一样。但是模板的默认类型必须由已有的模板参数推导出来，并不能设置默认类型为一个普通类型。在例子中已有的类型是Derived，那么traits必须使用Derived进行推导，换而言之，如果dense_xpr_base主模板中使用typename Xprkind = typename float这样的具体类型是不允许的。所以，主模板默认类型推导机制和traits技术是紧密联系的，两者通常联合使用。
上面几段文字，展示了极其复杂的traits技术，有必要搞这么复杂吗？这么复杂的traits技术事实上大大简化了代码，MatrixBase和ArrayBase两个基类不一样，生成的Matrix和Array却还有一些相同的操作，如果直接派生就需要同样的代码在子类Matrix和Array各写一次。而使用traits技术，将相同的代码统一在了PlainObjectBase中，这或许就是追求极简唯美。这种理念很像中文的成语，一个成语要花很多白话才能解释，但却可以表达出内涵深远的意境。程序要为我们的生活更好的服务，它的表达能力也至关重要的一方面。C++是众所周知的令人困惑的语言，但它的生命力就在它具有极强的表达能力。
PlainObjectBase利用traits技术根据传入类型的不同选择不同基类

2 矩阵初始化与访问

2.1 多种初始化形式

Eigen::Matrix类型提供了默认初始化，它不会分配动态内存，更不会初始化任何矩阵的元素。

// 一个3*3的矩阵，矩阵的元素都没有被初始化。
Matrix3f a;
// 一个动态矩阵，它的大小是0*0，也就是说还没有为该矩阵分配内存。
MatrixXf b;

重载构造函数提供指定矩阵大小初始化。对矩阵来说，第一个参数是矩阵的行数，第二个参数是列数。对向量来说只有一个参数，默认向量是一个列向量，即指定的参数可以认为是行数。这种初始化方式会分配矩阵或向量所需的内存，但是不会初始化元素的值。

// 一个10*15的动态矩阵，内存进行了分配，但没有初始化。
MatrixXf a(10,15);
// 一个大小为30的动态数组，但没有初始化。
VectorXf b(30);
// 行向量。
RowVector3f vecRow{1, 0, 0};

为了便于输入多个数据，支持逗号表达式初始化。该方法同样可以用于修改矩阵的值，即如果在下面代码再给m输入一次值，那么m的值被第二次值所覆盖。

Matrix3i m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;

用逗号表达式初始化矩阵或者向量的方式与常规数组的初始化习惯不一样，那么矩阵或者向量可以像初始化数组一样来初始化吗？答案是有类似的形式，但略有差别，这里称为数组式初始化。如果您的Eigen在以下初始化方式时提示错误，建议更新Eigen到最新稳定版，因为这种初始化方式最近两年才支持。

MatrixXi a{
   {2, 3, 4},
   {5, 6, 7},
};
Vector2i b{{1, 2}};

矩阵初始化时允许隐式类型转换，例如一个double型的立即数可以用于初始化float型。但是Eigen中算数操作符不支持隐式类型转换，例如一个float型的矩阵不能与double型的矩阵直接相加。
对比初始化与算数操作隐式类型转换

2.2 生成特殊矩阵

Eigen提供函数生成特定类型的矩阵，如单位矩阵、常量矩阵、随机矩阵等，详见下表。生成矩阵一方面可以用于矩阵的初始化赋值，另一方面可以在程序执行过程中修改矩阵的值。适用于矩阵的生成函数往往适用于向量，而仅针对向量设计的函数却不适用于矩阵。这些函数通常也支持对象式调用和函数式调用两种方式，支持函数式调用的函数如果不指定矩阵的大小就只能用于固定大小的矩阵。零矩阵和1矩阵的实现是调用常量矩阵，也就是说零矩阵和1矩阵只是常量矩阵的特例。
在这里插入图片描述

setIdentity()并不要求矩阵是方阵，然而在数学中常称的单位矩阵是一个方阵，这也是Eigen实现中与数学定义中有差异的地方。假设一个矩阵是mn的矩阵，且m>n，那么setIdentity的作用是将前n行n列的矩阵构成一个单位矩阵，剩下的m-n行全部为0。只能说Eigen中是设计成这样，而且这样设计逻辑上也很清晰。
另外一个疑问是，Random随机矩阵生成的是一个怎样的随机矩阵？总体来看矩阵的元素服从什么样的分布，多次生成同样维度的随机矩阵又服从什么样的分布？这两个问题，这里还不能回答，后续有专门章节来讨论这个问题。当前可以说的是，setRandom()会将矩阵的元素设置为-1~1之间的数。如果只是想生成一些不一样的数字用于测试，那么可以放心的去使用这些随机数；如果要依赖随机数的统计特性，那么对统计特性应该有一个验证。
Matrix通用矩阵生成函数
Eigen中向量是特殊的矩阵，按照从特殊到一般的金子塔思维模型，在金字塔底端的向量会比矩阵有更多的细节。符合这个思路，Eigen为生成向量提供了另外几个函数，等差数列(setLinSpaced)、单位向量(setUnit)。表中还介绍了几个特殊的单位向量，UnitX()、UnitY()、UnitZ()、UnitW()分别将单位向量的第0、1、2、3号元素设置为1，而且这四个函数都是调用Unit()实现的一个封装。通常我们认识的空间是三维的，因此三维向量特殊对待，Eigen提供了Vector3f::UnitX()、Vector3f::UnitY()、Vector3f::UnitZ()，这些函数分别生成正交坐标系下坐标基。
向量专用生成函数
从类的层次上来说，生成矩阵定义在DenseBase和MatrixBase两个层次中。定义在DenseBase中的生成函数也可以在Array中调用，而定义在MatrixBase中的函数只能在Matrix中调用。就像前面所说的，生成矩阵函数支持对象式调用和函数式调用，分别对应类中的成员函数和静态成员函数。
生成矩阵类层次图

2.3 访问矩阵

下标访问与矩阵大小。Eigen库重载圆括号()访问矩阵或者向量的元素，序号从0开始。Eigen不支持使用方括号[]访问矩阵的元素(向量除外)。元素的括号访问方式，相当于读取元素的引用，既可以作为左值，又可以作为右值。

MatrixXd m(2,2);
m(0,0) = 3;
VectorXd v(2);
v(0) = 4;

可以使用matrix.rows()、matrix.cols()、matrix.size()访问矩阵当前大小，使用matrix.resize()重置矩阵的大小。如果矩阵的大小没有变化，那么resize()操作没有任何影响。如果矩阵的大小改变了，那么矩阵的值可能会改变。如果想在重设大小过程中不改变矩阵的值使用matrix.conservativeResize()。使用赋值操作符“=”，Eigen将左操作数的大小重置为右操作数的大小。
对于非方阵，在矩阵分解时往往还需要知道对角线上元素的个数，其值等于min(rows, cols)，Eigen提供了函数来获取该值matrix. diagonalSize()。
和数组一样，下标越界也是矩阵元素访问中最常见的问题。Matrix越界可以分为两类，一类是未分配内存的矩阵，任何对元素的访问都会越界；另一类是已分配内存的矩阵，超过矩阵的范围会越界。除了越界的问题，还要警惕访问未初始化的元素，在有的机器上未初始化的元素是0，但有的机器却是随机值。也就是说访问未初始化的成员，可能导致未定义行为，而且这类问题往往很难发现。
已分配内存和未分配内存元素访问对比
迭代器访问。可以通过C风格式的编码遍历矩阵的元素，也可以使用迭代器式的方式遍历元素。对于向量可以直接调用begin()和end()获取迭代器，对于Matrix和Array要先使用colwise()、rowwise()或reshaped()转换成向量后再调用C++ STL标准库。

TEST(IterSTL, Vector)
{
    VectorXf vec(10);
    vec.setRandom();
    cout << "Here is a random vector vec" << endl;
    for (auto x : vec)
        cout << x << "\t";
    cout << endl;
    cout << "Here is sorted vec using std::sort" << endl;
    std::sort(vec.begin(), vec.end());
    cout << vec << endl;
}

2.4 块操作

块操作是从矩阵中生成子矩阵，既可以作为左值，又可以作为右值。块操作不仅适用于Matrix也适用于Array。需要传入参数的块操作需要用户保证，传入的参数不会使矩阵越界。
块是矩阵的一个矩形区域，块表达式可以是左值也可以是右值。最通用的是block()，它有两个版本的调用方式：

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

i、j表示块的起始位置，即块的左上角元素在矩阵中的角标，p、q表示块的大小。矩阵是固定大小时，块操作运行时效率要高。从图中可以看出，下标从0开始标号，且子矩阵包含指定位置(i, j)。
block操作示意图
矩阵计算中常常需要访问某一行或者某一列，Eigen提供了row()和col()来实现这一功能。图中展示了这两个函数的使用方法和产生的作用。
取1行或者1列元素
从边角框选一些元素
Eigen也提供了从某边某角开始框选元素的函数，如图所示。这些操作和block()一样，都支持两种形式。topLeftCorner(rows, cols)、bottomLeftCorner(rows, cols)、topRightCorner(rows, cols)、bottomRightCorner(rows, cols)分别是从左上角、左下角、右上角、右下角框选rows*cols维子矩阵。topRows(rows)、bottomRows(rows)表示分别从顶部、底部取rows行。leftCols(cols)、rightCols(cols)表示分别从左侧、右侧取cols列向量。middleCols(i, cols)表示从第i列起取cols列，middleRows(i, rows )表示从i列取rows行。
与特殊矩阵生成一样，对向量也提供了特殊的块操作函数。如图，vec.head(n)取向量的头n个元素，vec.tail(n)取向量的后n个元素，vec.segment(i, n)从第i个元素起，取n个元素。
向量支持的块操作

2.5 切片

如果Eigen现有的块操作函数不能满足需求，可以考虑使用切片自定义一些块操作函数。切片是从矩阵或者数组中截取特定的元素构成新的矩阵或数组。切片调用的核心函数是operator()(const RowIndices& rowIndices, const ColIndices& colIndices)，DenseBase通过包含IndexedViewMethods.h而获得该方法。
序列生成是切片实现的基础。ArithmeticSequence.h中实现了seqN函数生成序列，进一步封装了seq函数给开发者使用。seq本质上是返回一个ArithmeticSequence的模板类对象，里面只包含了描述序列的参数，并没有真正生成序列的所有元素。此外，Eigen还提供了空对象all用于占位，其占位含义是取所有的行或者列。还有last表示矩阵最后一行或者最后一列的下标。
在这里插入图片描述
在表中，介绍了seq实现与block API等价的例子。这几个例子只是为了说明block API可以由seq实现，而且当前block API没有实现的功能也可以使用seq实现。

2.6 Reshape

从Eigen 3.4开始，提供便捷的方法a.reshaped()，在保持元素值不变的情况下，改变矩阵的形状。reshaped()默认按列遍历元素，支持指定模板参数RowMajor按行遍历元素。从图中可以看出，reshaped(2,8)从原始矩阵的第一列开始取元素，然后放在新矩阵的第一列，不论是原始矩阵还是新矩阵，一列到达末尾后就进入新的一列。而reshaped(2,8)则是按行取元素，一行到达末尾跳到新的一行的起始位置。在不指定行数和列数时，reshaped()按列展开为1列，reshaped()按行展开为1列。
Reshaped示意图
reshaped()并不会改变调用对象，而只是返回一个副本。如果想原位改变形状可以调用resize()，但一个固定大小的向量是不允许resize()成一个矩阵的。与转置和求逆类似，a = a.reshaped(2,8)这种赋值给“自己”的形式是不允许的。

3 Matrix逻辑运算

3.1逻辑运算符

Matrix类型支持的逻辑运算操作符仅有“==”和“！=”，返回值为bool类型。如图，“>、>=、<、<=”等操作符都不支持，而且不支持矩阵与数字比较大小。
Matrix支持的逻辑运算符

4 矩阵向量代数

矩阵向量乘法需满足数学中对维度的约束，Eigen使用大量assert断言来表示这些约束，当约束不满足时调用abort()退出程序。

4.1 加减与数乘除

重载C++中“+”、“-”、“+=”、“-=”操作符，要求左右操作数的维度相同。不允许一个向量加上或者减去一个数。只要维度一样，向量和矩阵可以相加减，进一步说明向量是特殊的矩阵。
在这里插入图片描述

重载C++中“”、“/”、“=”、“/=”操作符，支持矩阵和向量乘以或者除以一个数。再次强调，比较浮点数大小，需使用isApprox()。

4.2 转置共轭倒序

转置是将位置在(i, j)处的元素与(j, i)的元素交换，共轭是将复数的虚部变为原来的反数，共轭转置是共轭与转置的结合。转置aT、共轭、共轭转置aH分别通过transpose()、conjugate()、adjoint()实现。调用格式a.transpose()，a.conjugate()，a.adjoint()。对于实数而言，共轭没有任何影响，共轭转置等价于转置。使用a = a.transpose()会出现错误，这是因为Eigen在进行转置或者共轭操作时，会同时写左操作数，从而得到意想不到的结果。要实现原位转置可以使用a.transposeInPlace()。类似的，也支持adjointInPlace()。
在这里插入图片描述
a.reverse()的行为与转置共轭类似，所以就放在本节中介绍。该函数是将一个矩阵的元素顺序反过来排列。

4.3 矩阵乘法

由于在Eigen中向量只是特殊的矩阵，因此只需重载“”、“=”即可实现矩阵和向量的乘法。如果你担心m=mm会导致错误，现在可以消除这个疑虑，因为Eigen以一种特殊的方式处理矩阵乘法，编译m=mm时，实际会构造一个中间变量tmp。

TEST(Matrix, Mul)
{
    Matrix3f a{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
    Matrix3f b{{0, 1, 0}, {1, 0, 0}, {0, 0, 1}};
    Matrix3f colSwap{{2, 1, 3}, {5, 4, 6}, {8, 7, 9}};
    Matrix3f rowSwap{{4, 5, 6}, {1, 2, 3}, {7, 8, 9}};
    Vector3f vec{1, 0, 0};
    Vector3f colExtract{1, 4, 7};
    RowVector3f vecRow{1, 0, 0};
    RowVector3f rowExtract{1, 2, 3};

    cout << "Here is the matrix a\n"
         << a << endl;
    cout << "Here is the matrix b\n"
         << b << endl;
    cout << "Here is a * b\n"
         << a * b << endl;
    EXPECT_TRUE(a * b == colSwap);
    cout << "Here is b * a\n"
         << b * a << endl;
    EXPECT_TRUE(b * a == rowSwap);

    ASSERT_EQ(vec.rows(), 3);
    ASSERT_EQ(vecRow.rows(), 1);

    cout << "Here is column major vector vec\n"
         << vec << endl;
    cout << "Here is row major vector vecRow\n"
         << vecRow << endl;
    cout << "Here is a * vec\n"
         << a * vec << endl;
    EXPECT_TRUE(a * vec == colExtract);
    cout << "Here is vecRow * a\n"
         << vecRow * a << endl;
    EXPECT_TRUE(vecRow * a == rowExtract);

    a *= b;
    cout << "Here is a *= b\n"
         << a << endl;
    EXPECT_TRUE(a == colSwap);
}

4.4 内积向量积

内积又可以称为点积，Eigen分别使用dot()和cross()来实现内积和向量积，向量积只适用于三维向量。

TEST(Matrix, Dot)
{
    Vector3d a{1, 2, 3};
    Vector3d b{2, 3, 4};

    cout << "Here is vector a\n"
         << a << endl;
    cout << "Here is vector b\n"
         << b << endl;
    cout << "Here is return type of Dot\n"
         << typeid(a.dot(b)).name() << endl;
    cout << "Here is a.dot(b)\n"
         << a.dot(b) << endl;
    EXPECT_EQ(a.dot(b), 20);
}