Eigen入门教程（1）矩阵的定义和赋值

myleader

http://forum.simwe.com/thread-989245-1-1.html
曾经介绍过线性代数软件包Eigen，其速度快、无需后端运行库、C++语言写成，非常适合大型软件的开发。这里准备介绍一下Eigen的简单使用方法，希望能稍微推广一下这个优秀的软件包。原本是想直接把Eigen的手册翻译过来的，后来想想那个帮助文档其实适合有一定基础的人看，不适合从零开始学习，最后决定重新写一篇来介绍。

这是上个帖子中使用的样例代码。

1. 
   
2. #include <iostream>
   
3. #include <ctime>
   
4. #include <Eigen/Core>
   
5. 
   
6. using namespace std;
   
7. using namespace Eigen;
   
8. 
   
9. int main()
   
10. {
    
11.         int size=2048, N=10;
    
12.         clock_t start, end;
    
13.         MatrixXf A = MatrixXf::Random(size,size);
    
14.         MatrixXf B = MatrixXf::Random(size,size);
    
15.         MatrixXf Z = MatrixXf::Random(size,size);
    
16.         cout <<"start"<<endl;
    
17.         start = clock();
    
18.         for(int i=0; i<N; ++i)
    
19.         {
    
20.                 Z = A*B;  //  or Z = A+B+C ... etc
    
21.                 cout<<i+1<<endl;
    
22.         }
    
23.         end = clock();
    
24.         cout << "time taken = " << 1.0*(end-start)/CLOCKS_PER_SEC << endl;
    
25.         system("pause");
    
26.         return(0);
    
27. }

复制代码

Eigen的模块

使用Eigen首先要包含它的头文件，Eigen提供了9个模块，另外还有一系列没有正式支持的模块提供一些扩展功能。

其中最重要的是Core模块，包含了Matrix和Array两个类以及基本的线性代数运算功能。Geometry模块包含了一些基本的平移、缩放和旋转的功能，并不像VTK、CGAL等提供了更多的图形计算和显示功能。LU模块包含了求逆、行列式、LU分解求解器等功能。Cholesky模块提供了Cholesky分解的功能。Householder模块包含了辅助变换功能。SVD模块包含了SVD分解和JacobiSVD求解（最小二乘法）功能。QR模块提供了QR分解。Eigenvalues模块提供了特征值计算的功能。Sparse模块提供了专门针对稀疏矩阵的一些功能。其中Core是一切的核心，你需要那部分的功能，就把那个头文件include进去就行了。

Eigen还提供了文件名为Dense的头文件，这个头文件包含了Core、Geometry、LU、Cholesky、SVD、QR和Eigenvalues 的头文件，也就是说只要包含一个Dense就包含了这么多所有的模块。不过C++编程一般不推荐包含一个巨大的头文件，这会使编译速度降低。所以还是建议分别包含。

然后必须要声明

1. using namespace Eigen;

复制代码

这样才能正常调用Eigen的函数。

矩阵和向量的定义

线性代数中最主要的操作单元式矩阵和向量，首先需要定义矩阵。

1. Matrix<double, 3, Dynamic, RowMajor> E;

复制代码

就定义了一个矩阵E，尖括号里面的4个参数分别指：矩阵中的变量类型、行数、列数、在内存中的排列方式。上面这行代码意思就是说定义一个double型、3行、列数不固定（可按需增长）、在内存中按行存放的矩阵E。其语法和模板的调用一样。

这个定义似乎太麻烦了，Eigen针对实际使用的需求定义了一些常用的typedef

1. Matrix3f P, Q, R;

复制代码

定义了3个3*3的方阵P、Q、R，这3个方阵的类型都是float型，行数和列数都是固定不变的，在内存中的存放顺序是按列存放。Matrix3f中红色的3代表方阵的尺寸，Eigen内部代码使用了SSE2、3、4系列指令集进行优化，而指令集中支持一条指令计算4个参数，所以这个红色的3可选的固定值是2、3、4，如果你希望矩阵的尺寸是浮动可变的，需要把3换成X；蓝色的f代表数据类型，f代表float，d代表double，i代表int，cf代表complex float，cd代表complex double，ci代表complex int；不能设定内存中排列的方式。在实际使用当中如果需要如上文E那么复杂的矩阵，那还是必须得自己按模板调用标准来写。

矢量的定义在Eigen软件包中也是按照Matrix来的，只不过列数或者行数为1罢了。在前面如果想像E那样定义一个矢量，就把尖括号中的第2项“3”改成“1”，而如果想像定义P、Q、R那样定义一个矢量则是

1. Vector3f x, y, z;
   
2. RowVector3f a, b, c;

复制代码

Eigen中默认都是按列存储的，这与fortran保持一致，可以更方便跨语言编程或者代码的移植。所以在Vector的定义中，没有ColumnVector。
http://eigen.tuxfamily.org/dox/group__matrixtypedefs.html
有更详细的说明。

而如果你使用X来定义了一个尺寸浮动的矩阵，那还可以

1. MatrixXd A(3,16);
   
2. VectorXf B(30);
   

复制代码

这样可以直接在定义矩阵时设定矩阵的初始尺寸，当然这是不能给定矩阵中各个元素的值的。出于API一致的考虑，对于Matrix3f这种设定了固定尺寸的矩阵，下面这种用法也是合法的，不过笔者强烈建议不要这么用

1. Matrix3f a(3,3);

复制代码

实际上对于较短的矢量，也还可以直接给定值来做构造函数的参数，不过如果你要给定整形数据的矩阵，那么构造函数内的参数到底是代表矩阵尺寸呢还是代表矩阵内数据的值呢？所以笔者建议也不要使用这种方法

1. Vector2d a(5.0, 6.0);
   
2. Vector3d b(5.0, 6.0, 7.0);
   
3. Vector4d c(5.0, 6.0, 7.0, 8.0);

复制代码

矩阵的赋值和读取

首先最简单的就是针对矩阵的每个元素进行赋值，并且也可以读取，使用起来和MATLAB差不多。要注意的是，Eigen中行、列的编号是从0开始，到n-1结束，这一点和MATLAB不同。而且不管矩阵在内存中是按行排列还是按列排列，编号的第一个总是代表行，第二个代表列。

1. 
   
2. #include <iostream>
   
3. #include <Eigen/Core>
   
4. using namespace Eigen;
   
5. int main()
   
6. {
   
7. MatrixXd m(2,2);
   
8. m(0,0) = 3;
   
9. m(1,0) = 2.5;
   
10. m(0,1) = -1;
    
11. m(1,1) = m(1,0) + m(0,1);
    
12. std::cout << "Here is the matrix m:\n" << m << std::endl;
    
13. VectorXd v(2);
    
14. v(0) = 4;
    
15. v(1) = v(0) - 1;
    
16. std::cout << "Here is the vector v:\n" << v << std::endl;
    
17. }
    
18. 
    

复制代码

这段代码的输出结果是

Here is the matrix m:
  3  -1
2.5 1.5
Here is the vector v:
4
3

实际上Eigen的Matrix类还重载了流操作符，可以用于输入，但是不能用于输出。如果你希望输出，那么std::cout的流操作符是可以直接输出的。

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

复制代码

但是要注意，流操作符的读入顺序是按行读入的，这和Matrix类按列存储矛盾，所以笔者强烈建议不要使用这种方法赋值。看一下输出结果

1 2 3
4 5 6
7 8 9

其实也可以用

1. A=B；

复制代码

来进行赋值，但是前提是A和B尺寸和类型都一样或者A是浮动尺寸的。A中原本的值全部都会被覆盖，如果A是动态尺寸的，A的大小也会跟着B改变。
还有一些特殊的构造函数可以直接生成一些特殊矩阵

1. typedef {Matrix3f|Array33f} FixedXD;
   
2. FixedXD x;
   
3. x = FixedXD::Zero();//生成0矩阵
   
4. x = FixedXD::Ones();//生成1矩阵
   
5. x = FixedXD::Constant(value);//生成某固定值的矩阵
   
6. x = FixedXD::Random();//生成随机数矩阵
   
7. x = FixedXD::LinSpaced(size, low, high);//生成等差数列
   
8. x.setZero();
   
9. x.setOnes();
   
10. x.setConstant(value);
    
11. x.setRandom();
    
12. x.setLinSpaced(size, low, high);

复制代码

还有另外一些不太常用的赋值函数

1. A.fill(10);       // Fill A with all 10's.
   
2. A.setRandom();    // Fill A with uniform random numbers in (-1, 1).
   
3.                   // Requires #include <Eigen/Array>.
   
4. A.setIdentity();  // Fill A with the identity.

复制代码

注意这些函数不能在定义矩阵时使用。

也可以对矩阵的子矩阵进行操作来赋值，不过分块的操作留待下篇文章再讲，所以这里只是演示一下

1. int rows=5, cols=5;
   
2. MatrixXf m(rows,cols);
   
3. m << (Matrix3f() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished(),
   
4. MatrixXf::Zero(3,cols-3),
   
5. MatrixXf::Zero(rows-3,3),
   
6. MatrixXf::Identity(rows-3,cols-3);
   
7. cout << m;

复制代码

这段代码的输出结果是

1 2 3 0 0
4 5 6 0 0
7 8 9 0 0
0 0 0 1 0
0 0 0 0 1

这个赋值函数太复杂了，楼主一般是不用的。

对于复数矩阵，要求输入参数为std::complex，如果复数矩阵是double型的，那么也要求std::complex是double型的，其他float和int以此类推。

1. Matrix2cd matcx1;
   
2. matcx1<<std::complex<double>(0,1), std::complex<double>(2,4.5),
   
3.                  std::complex<double>(100.3, 250), std::complex<double>(-3.1415926, -2.71828);

复制代码

复数矩阵也可以而通过下标赋值，不过也仍然要求输入为std::complex类型，使用下标访问复数矩阵，得到的也是一个std::complex型数据，在内存中是实部在前、虚步在后。实际上复数矩阵类包含两个函数real()和imag()分别访问实部和虚步，不过返回值是不可赋值的。换句话说，下面这段代码无法编译通过

1. matcx1(1,1).real()=-3.1415926;
   
2. matcx1(1,1).imag()=-2.71828;
   

复制代码

但是复数矩阵可以通过两个实数矩阵来构造，就是说real()和imag()允许等号右边的矩阵对象却不允许单纯的数字

1. matxc1.real()=Matrix2d::Random();
   
2. matxc1.imag()=Matrix2d::Random();

复制代码

就没问题

到这里请回过头看一下帖子最上面的代码，应该能看懂每一行了。

矩阵的定义和赋值就到这里，下一篇准备讲尺寸和子矩阵（块）的操作，敬请期待。

caoer

leader大作又出现了。谢谢分享

myleader

caoer 发表于 2012-11-26 01:40
leader大作又出现了。谢谢分享

公司效益不太好，再加上年底，不忙了。抽空把早就想写完的东西弄完，算是实现自己的许诺。

caoer

the first script i got this,
start ...
1
2
3
4
5
6
7
8
9
10
time taken = 26.815
请按任意键继续. . .

from windows 7

myleader

caoer 发表于 2012-12-9 08:42
the first script i got this,
start ...
1

我第一次使用也遇到了这个问题，不要忘了关闭debug模式，改用release，这样才会真的快起来。

还需要调整编译的优化参数，必须要启用SSE2/3/4指令集，而且还要开启OpenMP多线程，这样就可以大大加快。如果使用mingw编译器

-mmmx -msse -msse2 -msse3 -msse4a -fopenmp -lgomp

当然CPU的温度也会更高

我在奔腾G620、win7x64上release单线程只用了9.7s，

如果你使用MSVC，我不知道该怎么开启多线程支持，也不知道该怎么添加指令集支持，爱莫能助。

caoer

谢谢提示。
我用的是 codeblocks上的g++,只开了 o3, 也就是release mode.
windows上的g++ 是比 msvs慢。我的是 i3770k, 呵呵，看来编译器和flag的确很重要

caoer

add flags of MMX and SSE, it is interesting, proformance goes up..

start ...
1
2
3
4
5
6
7
8
9
10
time taken = 9.461
请按任意键继续. . .