Eigen入门教程（2）矩阵的尺寸操作和子矩阵操作

myleader

矩阵的尺寸操作
想对一个矩阵或者矢量的一部分进行操作，首先必须要知道矩阵的尺寸，否则对于一个只有3*3尺寸的矩阵操作第4列是不可能的。Eigen提供了3个函数获得矩阵的尺寸

1. x.size()          // length(x)        // vector size
   
2. C.rows()          // size(C)(1)       // number of rows
   
3. C.cols()          // size(C)(2)       // number of columns

复制代码

其中size()函数能获得矢量或者矩阵的尺寸，但是笔者不推荐size(C)(dim)的用法，因为编程时容易出错，尤其是混用矢量和矩阵的时候，矢量size()返回一个值，而矩阵size()返回2个值，然后用(dim)操作符读取对应值，如果你在一个列矢量上使用size()(1)会返回1，而这很可能不是你所需要的。推荐rows()和cols()的用法，这样在阅读代码时size()就是矢量，rows()和cols()就是矩阵，容易辨别和理解。另外不管矩阵在内存中是否设定了按行存储，dim对应的1都是行数，2都是列数。还有Eigen的矩阵不支持第3维，这是不如MATLAB和Armadillo的地方。

如果要改变矩阵的尺寸也是可以的

1. Matrix<double, 3, 3> A;               // Fixed rows and cols. Same as Matrix3d.
   
2. Matrix<double, 3, Dynamic> B;         // Fixed rows, dynamic cols.
   
3. x.resize(5);       //resize a vector
   
4. A.resize(4, 4);   // Runtime error if assertions are on.
   
5. B.resize(4, 9);   // Runtime error if assertions are on.
   
6. A.resize(3, 3);   // Ok; size didn't change.
   
7. B.resize(3, 9);   // Ok; only dynamic cols changed.

复制代码

注意，如果在矩阵定义时对某个维度设定了固定尺寸，那就不能在改变那个维度的尺寸，编译通不过。

resize()函数在内存中的操作实际上很危险。

如果是扩大尺寸，那么多出来的空间会被内存中那个地方原本存储的数据填充，换句话说是不可预知的。和C中的malloc差不多。

如果是缩小尺寸，和矩阵在内存中的排列方式有关。如果是按列排列，那么当你缩小行数时，第1列中的后几个元素会被挤压到第2列，而第2列更多的元素会被挤压到第3列，直到最后的几列的一部分被删去（其实是内存被释放）以此类推；如果缩小列数，则最后面的几列会被删去。

所以在使用resize()后，推荐大家手动将其中的元素重新赋值。

相对而言conservativeResize(i,j)函数更可靠一些。如果是扩大尺寸，它的效果和resize()是一样的，你只需手工将多出来的部分重新赋值就好了。如果是缩小尺寸，那么它会把尺寸上超出的部分删去，也就是说当你缩小行数时，它会把第1列截短，然后把第2列也截短，以此类推，不会把第1列多出来的部分挤压到第2列去。当然这意味着更多的内存操作，也就意味着更大的计算开销。它的两个参数中，还支持Eigen::NoChange变量，这样它就不会更改对应维度的尺寸。

1. C.conservativeResize(Eigen::NoChange,5)

复制代码

这样C的行数不变，列数变为5，Eigen::NoChange可以减小内存分配的开销。

为了照顾MATLAB用户的习惯，还有一个resizelike()函数

1. A.resizelike(B)  // equal to A.resize(B.rows(), B.cols())
   

复制代码

可以把A的尺寸变成和B一致，不过内部使用的是resize()而不是conservativeResize()

矩阵的块操作

如果希望访问某单个元素可以用C(i,j)的方式来访问，注意在MATLAB和Fortran中数组的编号是从1开始的，而在C++中，数组的编号是从0开始的。

如果对矩阵内元素的访问都采用下标操作的话，那Eigen这个线性代数包和std::vector就没什么区别了，实际上Eigen支持从矩阵中访问一部分，也就是子矩阵或者叫块，这个块也支持Eigen数学包的各种操作，尤其是SSE指令集优化的各种计算。

如果要访问一行或者一列的话，可以使用row()和col()函数，和上面那个获得矩阵尺寸的函数很像，少了一个字母s，使用时要小心，不过这两个函数都必须要求参数，而获得尺寸的函数不需要参数，所以代码阅读时应该不会有太大的问题。

1. mat1.row(i) = mat2.col(j);          // R(i, :) = P(:, j)
   
2. VectorXd vec=mat1.row(j);
   
3. mat1.col(j1).swap(mat1.col(j2));  // R(:, [j1 j2]) = R(:, [j2, j1])

复制代码

理论上讲Vector和RowVector其实内存中都只有一列数据，只不过可以被识别为列数为1或者行数为1，而对一维数据的赋值操作中只关心内容，不关心到底哪个方向，所以上面的代码可以编译通过并正确运行。

MATLAB支持以逻辑矩阵或者矢量作为提取子矩阵的参数，Eigen也可以

1. mat2=mat1.col(arr1<arr2);   //the size of mat1, arr1 and arr2 must be the same
   
2. mat2=mat1.col(vec1);    //vec1 must be integer vector
   
3. mat2=mat1.col(vec1).row(vec2);
   

复制代码

如果要访问一个连续的子矩阵，那么可以使用block()函数，这个函数有2种用法

1. mat1.block(i,j,rows,cols);
   
2. mat1.block<rows,cols>(i,j);
   

复制代码

其中i,j是子矩阵左上角的位置，注意C++中从0开始。rows,cols是子矩阵的尺寸，注意不能出现内存越界，否则编译时没问题，运行时会出错，而对于大尺寸矩阵而言，很难调试。这两种方法的主要差别在于速度，第2种方法速度更快，因为它是编译时优化的，但是第2种方法要求rows,cols必须是给定数，不能是变量，哪怕是定义为常量也是不允许的，下面这段代码编译时通不过。

1. const int rows=2, cols=6;
   
2. mat1.block<rows,cols>(i,j);

复制代码

如果你实在无法确定rows,cols那就老老实实使用第1种方法吧，它有更好的通用性。

对于矢量，则有类似的函数segment()

1. vec1.segment(pos,n);
   
2. vec1.segment<n>(pos);

复制代码

这个函数也有2个版本，和矩阵的block()函数类似，第2个尖括号版本要求你预先知道子矢量的尺寸，并给定n

对于一些常用的子矩阵，Eigen还提供了一些block()和segment()的typedef，使用起来更方便

适用于矢量的head()和tail()

1. vec1.head(n);
   
2. vec1.head<n>();
   
3. vec1.tail(n)
   
4. vec1.tail<n>()

复制代码

矩阵的话类似函数稍多一些

1. mat1.topLeftCorner(rows,cols)
   
2. mat1.topLeftCorner<rows,cols>()
   
3. 
   
4. mat1.topRightCorner(rows,cols)
   
5. mat1.topRightCorner<rows,cols>()
   
6. 
   
7. mat1.bottomLeftCorner(rows,cols)
   
8. mat1.bottomLeftCorner<rows,cols>()
   
9. 
   
10. mat1.bottomRightCorner(rows,cols)
    
11. mat1.bottomRightCorner<rows,cols>()
    
12. 
    
13. mat1.topRows(rows)
    
14. mat1.topRows<rows>()
    
15. 
    
16. mat1.bottomRows(rows)
    
17. mat1.bottomRows<rows>()
    
18. 
    
19. mat1.leftCols(cols)
    
20. mat1.leftCols<cols>()
    
21. 
    
22. mat1.rightCols(cols)
    
23. mat1.rightCols<cols>()

复制代码

接下来的这部分内容是官方文档中没有的，是楼主自己在实践中遇到的，算是一点经验吧。

block()函数的返回值是一个矩阵，支持大部分矩阵操作。不过返回的不是Matrix类，不过可以自动转换为Matrix类，所以说直接处理是不行的，经过中间变量处理一下会比较好。比如说下面这段代码

1. MatrixXd function1(const &MatrixXd);
   
2. void function2(&MatrixXd);
   
3. void function3(MatrixXd);
   
4. void function4(MatrixXd*);
   
5. 
   
6. MatrixXd mat1=somevalue;
   
7. MatrixXd mat2=mat1.block(i,j,rows,cols);
   
8. 
   
9. MatrixXd mat3=function1(mat1.block(i,j,rows,cols));                             // ok to run
   
10. mat1.block(i,j,rows,cols)=function1(mat2);                                           //ok
    
11. mat1.block(i,j,rows,cols)=function1(mat1.block(i,j,rows,cols));              // ok to run
    
12. 
    
13. function2(mat1.block(i,j,rows,cols));                       // can't compile
    
14. function2(mat2);                                                    //ok
    
15. mat1.block(i,j,rows,cols)=mat3;                              //ok
    
16. 
    
17. function3(mat1.block(i,j,rows,cols));                       // ok to run, but can't change mat1 value
    
18. 
    
19. function4(&mat1.block(i,j,rows,cols));                     // can't compile

复制代码

如果不改变输入矩阵的值，那么子矩阵作为参数输入是没有问题的。也就是function1

所以如果你想写一个改变矩阵值的函数，那么其输入参数不能是子矩阵，而必须用另外一个中间量接一下。先把子矩阵赋值给中间量，再把中间量作为参数传递进去，中间量矩阵值改变以后，再把中间量通过block()函数返回给原矩阵。也就是function2或者function1的最后一种使用方式，推荐function2的使用方式，不容易出错。

这样其实就意味着更多的计算，所以推荐把原矩阵引用和块的4个参数一同作为参数传递给这个函数，这样可以不复制大矩阵，并能完成各种操作。

如果你不使用引用，而是采用函数内复制对象的运算，那么它可以完成你想要的运算，但是不能改变原矩阵的值。这就是function3。在C++的大多数教材中都不推荐这种方法，所以楼主也不推荐。

而子矩阵的指针是不能传递的。所以function4编译通不过。

实际上Eigen中的矩阵类如果用指针来访问的话是很糟糕的，使用起来比较危险，因为Eigen内部实现了自己的模板和迭代器，而且指针本身想用好的话也很难，所以强烈不推荐一切指针操作。

Eigen和其他软件包之间的数据传递

很多时候我们需要和其他软件包交换数据以完成更丰富的功能，而std::vector通用性比较好，可以通过vector中转后在转化为其他软件包识别的类型。

在C++中，如果需要传递一个巨大的对象，深拷贝是非常消耗系统资源的，所以比较高效的方法是传递指针，这会带来一些危险，而Eigen解决了这个问题，不仅深拷贝速度快，而且还安全。

Eigen提供了Map函数来完成Eigen和std::vector之间的数据传递。看一下下面的代码。

1. int main()
   
2. {
   
3.     Matrix3f m;
   
4.     vector<float> n(9);
   
5. m << 1, 2, 3,
   
6.      4, 5, 6,
   
7.      7, 8, 9;
   
8. std::cout << m <<endl;
   
9. cout<<"now map the data from matrix3f to std::vector"<<endl;
   
10. Map<Matrix3f> (n.data(),3,3)=m;
    
11. for(std::vector<float>::iterator i = n.begin(); i != n.end(); ++i )
    
12. {
    
13.   cout<<*i<<endl;
    
14. }
    
15. cout<<"now map the data from vector to matrix"<<endl;
    
16. for(vector<float>::iterator i=n.begin();i!=n.end();++i)
    
17. {
    
18.     *i=9-(i-n.begin());
    
19.     cout<<*i<<endl;
    
20. }
    
21. Map<Matrix3f> k(n.data(),3,3);
    
22. cout<<k<<endl;
    
23. cout<<m<<endl;
    
24. for(std::vector<float>::iterator i = n.begin(); i != n.end(); ++i )
    
25. {
    
26.   cout<<*i<<endl;
    
27. }
    
28. system("pause");
    
29. }

复制代码

不管从Matrix到vector还是vector到Matrix，都要求两者的尺寸一致，不过两层vector的方式不支持。

红色的一行把Eigen的数据传递到std::vector里，绿色的一行是把vector传递到Matrix里并创建一个新的Matrix。

注意vector的data()函数会返回实际存储的数据的第一个元素的地址，也就是返回一个指针，Map函数就是通过这个指针来完成数据传递的。Map函数完成了深拷贝，所以旧的数据删掉或者更改都没有问题的。

如果是复数矩阵，这样就可以

1. MatrixXcd mat(m,n);
   
2. mat.real() = Map<MatrixXd>(realData,m,n);
   
3. mat.imag() = Map<MatrixXd>(imagData,m,n);

复制代码

对于复数矩阵，还有更高级的用法，因为std::complex在内存中是实部、虚部连续排列的，而std::vector在内存中也是按顺序连续排列的，所以可以将std::vector<double>指针强制转化为std::vector<std::complex>指针，然后用Map函数将其深拷贝

1. std::vector<double> raw = somedata;
   
2. std::complex<double> *ptr = reinterpret_cast<std::complex<double> *>(raw.data());
   
3. MatrixXcd mat(m,n) = Map<MatrixXcd>(ptr,m,n);

复制代码

这种方法理论上是可行的，在某些及特殊的情况下也许有用，楼主还是建议大家别这么干。理论上讲MatrixXcd复制到vector中并形成double列也是可行的，不过楼主同样不建议这么干。毕竟reinterpret_cast比较危险。

http://forum.kde.org/viewtopic.php?t=89083
说Map支持更多种类的数据传递，甚至可以访问MATLAB内存中的数据
http://eigen.tuxfamily.org/dox/TutorialMapClass.html
有更多介绍，有兴趣的可以认真阅读。

要注意的是Matrix类不能成为vector的元素，因为Matrix类内部自带了和vector不兼容的内存管理器，所以

1. std::vector<MatrixXd> vec_mat;

复制代码

是非法的，实际上编译可以通过，但是运行时会出错。
http://eigen.tuxfamily.org/dox/TopicStlContainers.html
提供了一些Eigen和stl容器联合运行的说明，就是说你在使用stl容器包含Matrix类的时候必须要指定Eigen的内存分配器，不过笔者亲自实践不太好使，运行出错，也许是楼主水平不够，如果有谁找到了好办法，请务必告知。

http://forum.kde.org/viewtopic.php?f=74&t=97365
说Eigen可以和boots::circular合用，不过楼主本人是不用boost::circular的，所以没有经验


这些内容应该足够完成矩阵的尺寸、块和批量数据操作了。下一篇准备讲Eigen的算数操作，敬请期待

caoer

又是一片原创经典，谢谢分享