Eigen 简介

2019-08-11
5 min read
  1. Short tutorial / The Matrix class / Quick reference guide /
  2. 矩阵计算 / Eigen 的数组操作 / Modules and Header files /
  3. CSCI2240 /

因为 Eigen 是 header-only 的,所以直接下载 Eigen 头文件即可使用,CMake 中引入 Eigen 的方法官方也提供了说明

Pitfalls

细节请参考官网 Common pitfalls

Matrix Class

Eigen 中,所有的矩阵和向量都是 Matrix 模板类。向量是特殊的矩阵,有着一行或者一列数据

矩阵元素引用

Eigen 中主要的矩阵元素引用方式是使用重载后的小括号。对于矩阵,行索引总是第一个参数。对于向量而言,只需要一个索引值就可以了。所有索引值都从 0 开始。只有一个索引值的元素引用语法 m(index)并不仅限于向量,也适用于一般的矩阵,意味着读取一行(列)元素。然而这个特性十分依赖Eigen的存储顺序。所有的 Eigen 矩阵默认按列存储在内存中,不过你可以手动修改存储顺序,可参考 Storage orders

运算符 [] 也被重载用于引用向量的元素,但是因为 C++ 语言的特性,[]不能用于矩阵元素的引用,因为在C++中[i,j][j]的计算值是相同的,都是[j]

逗号初始化

矩阵和向量都可以使用逗号初始化语法(其他初始化方式):

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

维度调整

Matrix 的尺寸信息可以通过 rows()、cols() 和 size() 获取,分别返回矩阵的行数、列数和所有元素的个数。矩阵尺寸调整可以通过方法 resize()

如果矩阵的实际尺寸没有发生变化,那么resize() 方法将不进行任何实际操作。如果使用 resize() 调整了矩阵的尺寸,那么矩阵元素的值可能发生变化。如果你需要一个保守的不改变矩阵元素值的resize方法,可以参考 conservativeResize()

为了实现API的统一,所有调整矩阵尺寸的方法都适用于固定尺寸的矩阵,不过大部分时候这些函数都会报错或者不做任何实际操作。对于固定尺寸的矩阵 Matrix4d m,做这样的调整m.resize(4,4)是不会报错的,因为m并没有改变尺寸

定维 vs 动态

什么时候使用定维矩阵?什么时候使用动态矩阵?比较好的答案是:当矩阵尺寸比较小的时候使用定维矩阵,当矩阵尺寸比较大或者你不能使用定维矩阵时使用动态矩阵。对于定维矩阵而言,当尺寸较小时其执行效率要高于动态矩阵。定维矩阵,Eigen 将不会为其分配动态内存。在 Eigen 中定维矩阵在内存中的实现方式是普通的数组,Matrix4f m 的实现方式和float m[16]是非常类似的,而 MatrixXf mx 的实现方式和 float mx*=new float[rows*cols] 是类似的。

当矩阵元素个数大于32或者更多时,固定尺寸的性能优势将不再明显。创建一个非常大的固定尺寸矩阵很可能造成栈溢出,因为Eigen在栈中创建定维矩阵。根据不同的编译与运行环境,Eigen可能使用更激进的向量化手段(例如使用SIMD指令),这些特性可以参考 Vectorization

内置类型

Eigen 预定义了以下 Matrix 类型

  1. MatrixNt ,等价于 Matrix<type, N, N>
  2. VectorNt,等价于 Matrix<type, N, 1>
  3. RowVectorNt,等价 Matrix<type, 1, N>

N 的取值可以是 2、3、4、或者 X(X表示动态);t 可以是 i、f、d、cf(complex<float>)、cd 或者用户自定义类型

矩阵计算

  1. 矩阵计算,主要对矩阵、向量和标量之间的计算做一些简要介绍
  2. Eigen 的数组操作 /

Array

Array 的加减运算是 element-wise 的,也就是相同尺寸 Array 对应元素之间的加减。Array 也提供了 Array 和标量之间的加减运算(这是与 Matrix 的不同之处):array+scalar,Array 每个元素都与 scalar 相加

Matrix&Array

Eigen 中 Matrix 和 Array 的职能不同,前者完成线性代数相关算法,后者完成元素相关的计算。有些时候我们需要同时对数据进行线性代数相关的操作和元素相关的操作,这个时候就需要在 Matrix 和 Array 之间进行转换

Matrix 有一个 array() 方法可将 Matrix 转化为 Array 对象。同样,Array 有一个 matrix() 方法可将 Array 对象转化为 Matrix 对象。和 Eigen 表达式优化一样,这种转换没有运行时消耗(假设你允许编译器优化)。.array() 方法和 .matrix() 方法的返回值均可以当做左值或者右值进行使用

Eigen 不允许在同一个表达式中同时使用 Array 和 Matrix 对象,也就是说你不能直接将 Array 和 Matrix 对象相加。重载的运算符 +,其两侧要么全是 Array,要么全是 Matrix 对象。赋值运算符是个例外,你可以将 Matrix d对象赋予 Array 对象,同样你也可以将 Array 对象赋予 Matrix对象

线性方程组

Eigen 提供了稠密矩阵和稀疏矩阵,对应由线性方程组生成的系数矩阵也有稀疏和稠密之别。不同矩阵有不同的特性,所以 Eigen 同时提供了稠密和稀疏线性方程组的算法。细节请参考官网

高性能计算

延迟计算

延迟计算(Lazy Evaluation是 Eigen 运行时特性,以减少计算量。在 Eigen 中算术运算符并不做实际的计算,这些运算符只是返回一个被标识要做相关计算的对象,真实的计算发生在对整个表达式求值的时候,一般是遇见赋值等号的时候。这样做有利于编译器做优化。举个例子:

VectorXf a(50), b(50), c(50), d(50);
...
a = 3*b + 4*c + 5*d;

Eigen 将上面的表达式编译为一个 for 循环,在不考虑其他优化(如 SIMD)的情况下可以等价为如下代码:

for(int i = 0; i < 50; ++i)
  a[i] = 3*b[i] + 4*c[i] + 5*d[i];

译者注:如果所有的实际计算发生在运算符出现的地方,那么上面的代码最少需要三个 for循环分别迭代b、c、d这三个矩阵,和一个for循环比起来在循环次数上后者是前者的1/3

所以不要怕使用非常长的算术表达式,长的算术表达式便于编译器做优化

并行计算

  1. Eigen 支持 CPU 的 SIMD 指令,例如 SSE2
  2. Eigen 支持多线程加速,相关细节可以参考官网
  3. Eigen 支持更换底层线性代数库,例如 BLAS、MKL 等,Eigen 3.3 后 Eigen 也支持 CUDA,实现 GPU 计算