Eigen入门教程（2）矩阵的尺寸操作和子矩阵操作
本帖最后由 myleader 于 2012-12-30 17:26 编辑矩阵的尺寸操作
想对一个矩阵或者矢量的一部分进行操作，首先必须要知道矩阵的尺寸，否则对于一个只有3*3尺寸的矩阵操作第4列是不可能的。Eigen提供了3个函数获得矩阵的尺寸
x.size() // length(x) // vector size
C.rows() // size(C)(1) // number of rows
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的地方。
如果要改变矩阵的尺寸也是可以的
Matrix<double, 3, 3> A; // Fixed rows and cols. Same as Matrix3d.
Matrix<double, 3, Dynamic> B; // Fixed rows, dynamic cols.
x.resize(5); //resize a vector
A.resize(4, 4); // Runtime error if assertions are on.
B.resize(4, 9); // Runtime error if assertions are on.
A.resize(3, 3); // Ok; size didn't change.
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变量，这样它就不会更改对应维度的尺寸。
C.conservativeResize(Eigen::NoChange,5)这样C的行数不变，列数变为5，Eigen::NoChange可以减小内存分配的开销。
为了照顾MATLAB用户的习惯，还有一个resizelike()函数
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，使用时要小心，不过这两个函数都必须要求参数，而获得尺寸的函数不需要参数，所以代码阅读时应该不会有太大的问题。
mat1.row(i) = mat2.col(j); // R(i, :) = P(:, j)
VectorXd vec=mat1.row(j);
mat1.col(j1).swap(mat1.col(j2));// R(:, ) = R(:, )理论上讲Vector和RowVector其实内存中都只有一列数据，只不过可以被识别为列数为1或者行数为1，而对一维数据的赋值操作中只关心内容，不关心到底哪个方向，所以上面的代码可以编译通过并正确运行。
MATLAB支持以逻辑矩阵或者矢量作为提取子矩阵的参数，Eigen也可以
mat2=mat1.col(arr1<arr2); //the size of mat1, arr1 and arr2 must be the same
mat2=mat1.col(vec1); //vec1 must be integer vector
mat2=mat1.col(vec1).row(vec2);
如果要访问一个连续的子矩阵，那么可以使用block()函数，这个函数有2种用法
mat1.block(i,j,rows,cols);
mat1.block<rows,cols>(i,j);
其中i,j是子矩阵左上角的位置，注意C++中从0开始。rows,cols是子矩阵的尺寸，注意不能出现内存越界，否则编译时没问题，运行时会出错，而对于大尺寸矩阵而言，很难调试。这两种方法的主要差别在于速度，第2种方法速度更快，因为它是编译时优化的，但是第2种方法要求rows,cols必须是给定数，不能是变量，哪怕是定义为常量也是不允许的，下面这段代码编译时通不过。const int rows=2, cols=6;
mat1.block<rows,cols>(i,j);
如果你实在无法确定rows,cols那就老老实实使用第1种方法吧，它有更好的通用性。
对于矢量，则有类似的函数segment()
vec1.segment(pos,n);
vec1.segment<n>(pos);
这个函数也有2个版本，和矩阵的block()函数类似，第2个尖括号版本要求你预先知道子矢量的尺寸，并给定n
对于一些常用的子矩阵，Eigen还提供了一些block()和segment()的typedef，使用起来更方便
适用于矢量的head()和tail()
vec1.head(n);
vec1.head<n>();
vec1.tail(n)
vec1.tail<n>()矩阵的话类似函数稍多一些
mat1.topLeftCorner(rows,cols)
mat1.topLeftCorner<rows,cols>()
mat1.topRightCorner(rows,cols)
mat1.topRightCorner<rows,cols>()
mat1.bottomLeftCorner(rows,cols)
mat1.bottomLeftCorner<rows,cols>()
mat1.bottomRightCorner(rows,cols)
mat1.bottomRightCorner<rows,cols>()
mat1.topRows(rows)
mat1.topRows<rows>()
mat1.bottomRows(rows)
mat1.bottomRows<rows>()
mat1.leftCols(cols)
mat1.leftCols<cols>()
mat1.rightCols(cols)
mat1.rightCols<cols>()接下来的这部分内容是官方文档中没有的，是楼主自己在实践中遇到的，算是一点经验吧。
block()函数的返回值是一个矩阵，支持大部分矩阵操作。不过返回的不是Matrix类，不过可以自动转换为Matrix类，所以说直接处理是不行的，经过中间变量处理一下会比较好。比如说下面这段代码
MatrixXd function1(const &MatrixXd);
void function2(&MatrixXd);
void function3(MatrixXd);
void function4(MatrixXd*);
MatrixXd mat1=somevalue;
MatrixXd mat2=mat1.block(i,j,rows,cols);
MatrixXd mat3=function1(mat1.block(i,j,rows,cols)); // ok to run
mat1.block(i,j,rows,cols)=function1(mat2); //ok
mat1.block(i,j,rows,cols)=function1(mat1.block(i,j,rows,cols)); // ok to run
function2(mat1.block(i,j,rows,cols)); // can't compile
function2(mat2); //ok
mat1.block(i,j,rows,cols)=mat3; //ok
function3(mat1.block(i,j,rows,cols)); // ok to run, but can't change mat1 value
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之间的数据传递。看一下下面的代码。
int main()
{
Matrix3f m;
vector<float> n(9);
m << 1, 2, 3,
4, 5, 6,
7, 8, 9;
std::cout << m <<endl;
cout<<"now map the data from matrix3f to std::vector"<<endl;
Map<Matrix3f> (n.data(),3,3)=m;
for(std::vector<float>::iterator i = n.begin(); i != n.end(); ++i )
{
cout<<*i<<endl;
}
cout<<"now map the data from vector to matrix"<<endl;
for(vector<float>::iterator i=n.begin();i!=n.end();++i)
{
*i=9-(i-n.begin());
cout<<*i<<endl;
}
Map<Matrix3f> k(n.data(),3,3);
cout<<k<<endl;
cout<<m<<endl;
for(std::vector<float>::iterator i = n.begin(); i != n.end(); ++i )
{
cout<<*i<<endl;
}
system("pause");
}不管从Matrix到vector还是vector到Matrix，都要求两者的尺寸一致，不过两层vector的方式不支持。
红色的一行把Eigen的数据传递到std::vector里，绿色的一行是把vector传递到Matrix里并创建一个新的Matrix。
注意vector的data()函数会返回实际存储的数据的第一个元素的地址，也就是返回一个指针，Map函数就是通过这个指针来完成数据传递的。Map函数完成了深拷贝，所以旧的数据删掉或者更改都没有问题的。
如果是复数矩阵，这样就可以
MatrixXcd mat(m,n);
mat.real() = Map<MatrixXd>(realData,m,n);
mat.imag() = Map<MatrixXd>(imagData,m,n);
对于复数矩阵，还有更高级的用法，因为std::complex在内存中是实部、虚部连续排列的，而std::vector在内存中也是按顺序连续排列的，所以可以将std::vector<double>指针强制转化为std::vector<std::complex>指针，然后用Map函数将其深拷贝
std::vector<double> raw = somedata;
std::complex<double> *ptr = reinterpret_cast<std::complex<double> *>(raw.data());
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不兼容的内存管理器，所以
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的算数操作，敬请期待
又是一片原创经典，谢谢分享
页:
[1]