eigen：矩阵处理工具.docx

C+开源矩阵计算工具Eigen简单用法（一）矩阵转置函数m.transpose(); 矩阵的特征值m.eigenvalues();矩阵求逆矩阵m.inverse(); 特征向量m.eignvectors();1、 矩阵的定义Eigen中关于矩阵类的模板函数中，共有6个模板参数，但是目前常用的只有前三个，如下所示：cppview plaincopy1. template2. structtraitsMatrix3. .其前三个参数分别表示矩阵元素的类型，行数和列数。矩阵定义时可以使用Dynamic来表示矩阵的行列数为未知，例如：typedefMatrixMatrixXd;在Eigen中也提供了很多常见的简化定义形式，例如：typedef Matrix Vector3d注意：（1）Eigen中无论是矩阵还是数组、向量，无论是静态矩阵还是动态矩阵都提供默认构造函数，也就是你定义这些数据结构时都可以不用提供任何参数，其大小均由运行时来确定。（2）矩阵的构造函数中只提供行列数、元素类型的构造参数，而不提供元素值的构造，对于比较小的、固定长度的向量提供初始化元素的定义，例如：cppview plaincopy1. Vector2da(5.0,6.0);2. Vector3db(5.0,6.0,7.0);3. Vector4dc(5.0,6.0,7.0,8.0);2、动态矩阵和静态矩阵动态矩阵是指其大小在运行时确定，静态矩阵是指其大小在编译时确定，在Eigen中并未这样称呼矩阵。具体可见如下两段代码：代码段1：cppview plaincopy1. #include2. #include3. usingnamespaceEigen;4. usingnamespacestd;5. intmain()6. 7. MatrixXdm=MatrixXd:Random(3,3);8. m=(m+MatrixXd:Constant(3,3,1.2)*50;9. coutm=endlmendl;10. VectorXdv(3);11. v1,2,3;12. coutm*v=endlm*vendl;13. 代码段2：cppview plaincopy1. #include2. #include3. usingnamespaceEigen;4. usingnamespacestd;5. intmain()6. 7. Matrix3dm=Matrix3d:Random();8. m=(m+Matrix3d:Constant(1.2)*50;9. coutm=endlmendl;10. Vector3dv(1,2,3);11. coutm*v=endlm*vendl;12. 说明：1）代码段1中MatrixXd表示任意大小的元素类型为double的矩阵变量，其大小只有在运行时被赋值之后才能知道；MatrixXd:Random(3,3)表示产生一个元素类型为double的3*3的临时矩阵对象。2） 代码段2中Matrix3d表示元素类型为double大小为3*3的矩阵变量，其大小在编译时就知道；3）上例中向量的定义也是类似，不过这里的向量时列优先，在Eigen中行优先的矩阵会在其名字中包含有row，否则就是列优先。4）向量只是一个特殊的矩阵，其一个维度为1而已，如：typedef Matrix Vector3d3、矩阵元素的访问在矩阵的访问中，行索引总是作为第一个参数，需注意Eigen中遵循大家的习惯让矩阵、数组、向量的下标都是从0开始。矩阵元素的访问可以通过（）操作符完成，例如m(2,3)即是获取矩阵m的第2行第3列元素（注意行列数从0开始）。可参看如下代码：cppview plaincopy1. #include2. #include3. usingnamespaceEigen;4. intmain()5. 6. MatrixXdm(2,2);7. m(0,0)=3;8. m(1,0)=2.5;9. m(0,1)=-1;10. m(1,1)=m(1,0)+m(0,1);11. std:coutHereisthematrixm:nmstd:endl;12. VectorXdv(2);13. v(0)=4;14. v(1)=v(0)-1;15. std:coutHereisthevectorv:nvstd:endl;16. 其输出结果为：Here is the matrix m: 3 -12.5 1.5Here is the vector v:43针对向量还提供操作符，注意矩阵则不可如此使用，原因为：在C+中mi, j中逗号表达式 “i, j”的值始终都是“j”的值，即mi, j对于C+来讲就是mj;4、设置矩阵的元素在Eigen中重载了操作符，通过该操作符即可以一个一个元素的进行赋值，也可以一块一块的赋值。另外也可以使用下标进行复制，例如下面两段代码：代码段1cppview plaincopy1. Matrix3fm;2. m1,2,3,3. 4,5,6,4. 7,8,9;5. std:coutm;输出结果为：1 2 34 5 67 8 9代码段二（使用下标进行复制）cppview plaincopy1. VectorXfm_Vector_A;2. MatrixXfm_matrix_B;3. intm_iN=-1;4. 5. boolInitData(intpSrc100100,intiWidth,intiHeight)6. 7. if(NULL=pSrc|iWidth=0|iHeight=0)8. returnfalse;9. m_iN=iWidth*iHeight;10. VectorXftmp_A(m_iN);11. MatrixXftmp_B(m_iN,5);12. inti=0,j=0,iPos=0;13. while(iiWidth)14. 15. j=0;16. while(jiHeight)17. 18. tmp_A(iPos)=pSrcij*log(float)pSrcij);19. 20. tmp_B(iPos,0)=pSrcij;21. tmp_B(iPos,1)=pSrcij*i;22. tmp_B(iPos,2)=pSrcij*j;23. tmp_B(iPos,3)=pSrcij*i*i;24. tmp_B(iPos,4)=pSrcij*j*j;25. +iPos;26. +j;27. 28. +i;29. 30. m_Vector_A=tmp_A;31. m_matrix_B=tmp_B;32. 5、重置矩阵大小当前矩阵的行数、列数、大小可以通过rows()，cols()和size()来获取，对于动态矩阵可以通过resize()函数来动态修改矩阵的大小.需注意：(1)固定大小的矩阵是不能使用resize（）来修改矩阵的大小；(2)resize（）函数会析构掉原来的数据，因此调用resize（）函数之后将不能保证元素的值不改变。(3) 使用“=”操作符操作动态矩阵时，如果左右边的矩阵大小不等，则左边的动态矩阵的大小会被修改为右边的大小。例如下面的代码段：cppview plaincopy1. MatrixXfa(2,2);2. std:coutaisofsizea.rows()xa.cols()std:endl;3. MatrixXfb(3,3);4. a=b;5. std:coutaisnowofsizea.rows()xa.cols()std:endl;输出结果为：a is of size 2x2a is now of size 3x36、如何选择动态矩阵和静态矩阵？Eigen对于这问题的答案是：对于小矩阵（一般大小小于16）的使用固定大小的静态矩阵，它可以带来比较高的效率，对于大矩阵（一般大小大于32）建议使用动态矩阵。还需特别注意的是：如果特别大的矩阵使用了固定大小的静态矩阵则可能造成栈溢出的问题本文主要是Eigen中矩阵和向量的算术运算，在Eigen中的这些算术运算重载了C+的+，-，*，所以使用起来非常方便。1、矩阵的运算Eigen提供+、-、一元操作符“-”、+=、-=，例如：二元操作符+/-表示两矩阵相加（矩阵中对应元素相加/减，返回一个临时矩阵）： B+C 或B-C；一元操作符-表示对矩阵取负（矩阵中对应元素取负，返回一个临时矩阵）：-C;组合操作法+=或者-=表示（对应每隔元素都做相应操作）：A += B 或者 A-=B代码段1为矩阵的加减操作，代码如下：cppview plaincopy1. #include2. #include3. usingnamespaceEigen;4. intmain()5. 6. Matrix2da;7. a1,2,8. 3,4;9. MatrixXdb(2,2);10. b2,3,11. 1,4;12. std:couta+b=na+bstd:endl;13. std:couta-b=na-bstd:endl;14. std:coutDoinga+=b;std:endl;15. a+=b;16. std:coutNowa=nastd:endl;17. Vector3dv(1,2,3);18. Vector3dw(1,0,0);19. std:cout-v+w-v=n-v+w-vstd:endl;20. 输出结果为：a + b =3 54 8a - b =-1 -1 2 0Doing a += b;Now a =3 54 8-v + w - v =-1-4-6另外，矩阵还提供与标量（单一个数字）的乘除操作，表示每个元素都与该标量进行乘除操作。例如：二元操作符*在：A*a中表示矩阵A中的每隔元素都与数字a相乘，结果放在一个临时矩阵中，矩阵的值不会改变。对于a*A、A/a、A*=a、A /=a也是一样，例如下面的代码：cppview plaincopy1. #include2. #include3. usingnamespaceEigen;4. intmain()5. 6. Matrix2da;7. a1,2,8. 3,4;9. Vector3dv(1,2,3);10. std:couta*2.5=na*2.5std:endl;11. std:cout0.1*v=n0.1*vstd:endl;12. std:coutDoingv*=2;std:endl;13. v*=2;14. std:coutNowv=nvstd:endl;15. 输出结果为：a * 2.5 =2.5 57.5 100.1 * v =Doing v *= 2;Now v =246需要注意：在Eigen中，算术操作例如 “操作符+”并不会自己执行计算操作，他们只是返回一个“算术表达式对象”，而实际的计算则会延迟到后面的赋值时才进行。这些不影响你的使用，它只是为了方便Eigen的优化。2、求矩阵的转秩、共轭矩阵、伴随矩阵。可以通过成员函数transpose(),conjugate(),和adjoint()来完成，注意这些函数返回操作后的结果，而不会对原矩阵的元素进行直接操作，如果要让原矩阵的进行转换，则需要使用响应的InPlace函数，例如：transposeInPlace()、adjointInPlace()之类。例如下面的代码所示：cppview plaincopy1. MatrixXcfa=MatrixXcf:Random(2,2);2. coutHereisthematrixanaendl;3. coutHereisthematrixaTna.transpose()endl;4. coutHereistheconjugateofana.conjugate()endl;5. coutHereisthematrixa*na.adjoint()endl;输出结果为：Here is the matrix a (-0.211,0.68) (-0.605,0.823) (0.597,0.566) (0.536,-0.33)Here is the matrix aT(-0.211,0.68) (0.597,0.566)(-0.605,0.823) (0.536,-0.33)Here is the conjugate of a (-0.211,-0.68) (-0.605,-0.823) (0.597,-0.566) (0.536,0.33)Here is the matrix a*(-0.211,-0.68) (0.597,-0.566)(-0.605,-0.823) (0.536,0.33)3、矩阵相乘、矩阵向量相乘矩阵的相乘，矩阵与向量的相乘也是使用操作符*，共有*和*=两种操作符，其用法可以参考如下代码：cppview plaincopy1. #include2. #include3. usingnamespaceEigen;4. intmain()5. 6. Matrix2dmat;7. mat1,2,8. 3,4;9. Vector2du(-1,1),v(2,0);10. std:coutHereismat*mat:nmat*matstd:endl;11. std:coutHereismat*u:nmat*ustd:endl;12. std:coutHereisuT*mat:nu.transpose()*matstd:endl;13. std:coutHereisuT*v:nu.transpose()*vstd:endl;14. std:coutHereisu*vT:nu*v.transpose()std:endl;15. std:coutLetsmultiplymatbyitselfstd:endl;16. mat=mat*mat;17. std:coutNowmatismat:nmatstd:endl;18. 输出结果为：Here is mat*mat: 7 1015 22Here is mat*u:11Here is uT*mat:2 2Here is uT*v:-2Here is u*vT:-2 -0 2 0Lets multiply mat by itselfNow mat is mat: 7 1015 22本节主要涉及Eigen的块操作以及QR分解1、矩阵的块操作1）矩阵的块操作有两种使用方法，其定义形式为：cppview plaincopy1. matrix.block(i,j,p,q);（1）2. 3. matrix.block(i,j);（2）定义（1）表示返回从矩阵的(i, j)开始，每行取p个元素，每列取q个元素所组成的临时新矩阵对象，原矩阵的元素不变。定义（2）中block（p, q）可理解为一个p行q列的子矩阵，该定义表示从原矩阵中第(i, j)开始，获取一个p行q列的子矩阵，返回该子矩阵组成的临时 矩阵对象，原矩阵的元素不变。详细使用情况，可参考下面的代码段：cppview plaincopy1. #include2. #include3. usingnamespacestd;4. intmain()5. 6. Eigen:MatrixXfm(4,4);7. m1,2,3,4,8. 5,6,7,8,9. 9,10,11,12,10. 13,14,15,16;11. coutBlockinthemiddleendl;12. coutm.block(1,1)endlendl;13. for(inti=1;i=3;+i)14. 15. coutBlockofsizeixiendl;16. coutm.block(0,0,i,i)endlendl;17. 18. 输出的结果为：Block in the middle 6 710 11Block of size 1x11Block of size 2x21 25 6Block of size 3x3 1 2 3 5 6 7 9 10 11通过上述方式获取的子矩阵即可以作为左值也可以作为右值，也就是即可以用这个子矩阵给其他矩阵赋值，也可以给这个子矩阵对象赋值。2）矩阵也提供了获取其指定行/列的函数，其实获取某行/列也是一种特殊的获取子块。可以通过 .col()和 .row()来完成获取指定列/行的操作，参数为列/行的索引。注意：（1）需与获取矩阵的行数/列数的函数（ rows(), cols() ）的进行区别，不要弄混淆。（2）函数参数为响应行/列的索引，需注意矩阵的行列均以0开始。下面的代码段用于演示获取矩阵的指定行列：cppview plaincopy1. #include2. #include3. usingnamespacestd;4. intmain()5. 6. Eigen:MatrixXfm(3,3);7. m1,2,3,8. 4,5,6,9. 7,8,9;10. coutHereisthematrixm:endlmendl;11. cout2ndRow:m.row(1)endl;12. m.col(2)+=3*m.col(0);13. coutAfteradding3timesthefirstcolumnintothethirdcolumn,thematrixmis:n;14. coutmendl;15. 输出结果为：Here is the matrix m:1 2 34 5 67 8 92nd Row: 4 5 6After adding 3 times the first column into the third column, the matrix m is: 1 2 6 4 5 18 7 8 303）向量的块操作，其实向量只是一个特殊的矩阵，但是Eigen也为它单独提供了一些简化的块操作，如下三种形式： 获取向量的前n个元素：vector.head(n); 获取向量尾部的n个元素：vector.tail(n); 获取从向量的第i个元素开始的n个元素：vector.segment(i,n); 其用法可参考如下代码段：cppview plaincopy1. #include2. #include3. usingnamespacestd;4. intmain()5. 6. Eigen:ArrayXfv(6);7. v1,2,3,4,5,6;8. coutv.head(3)=endlv.head(3)endlendl;9. coutv.tail()=endlv.tail()endlendl;10. v.segment(1,4)*=2;11. coutafterv.segment(1,4)*=2,v=endlvendl;12. 输出结果为：v.head(3) =123v.tail() = 456after v.segment(1,4) *= 2, v =14681062、QR分解 Eigen的QR分解非常绕人，它总共提供了下面这些矩阵的分解方式：DecompositionMethodRequirements on the matrixSpeedAccuracyPartialPivLUpartialPivLu()Invertible+FullPivLUfullPivLu()None-+HouseholderQRhouseholderQr()None+ColPivHouseholderQRcolPivHouseholderQr()None+FullPivHouseholderQRfullPivHouseholderQr()None-+LLTllt()Positive definite+LDLTldlt()Positive or negative semidefinite+由于我只用到了QR分解，而且Eigen的QR分解开始使用时确实不容易入手，因此这里只提供了householderQR的分解方式的演示代码：cppview plaincopy1. voidQR2()2. 3. Matrix3dA;4. A1,1,1,5. 2,-1,-1,6. 2,-4,5;7. 8. HouseholderQRqr;9. pute(A);10. MatrixXdR=qr.matrixQR().triangularView();11. MatrixXdQ=qr.householderQ();12. std:coutQR2():HouseholderQR-std:endl;13. std:coutAstd:endlAstd:endlstd:endl;14. std:coutqr.matrixQR()std:endlqr.matrixQR()std:endlstd:endl;15. std:coutRstd:endlRstd:endlstd:endl;16. std:coutQstd:endlQstd:endlstd:endl;17. std:coutQ*Rstd:endlQ*Rstd:endlstd:endl;18. 输出结果为：3、一个矩阵使用的例子：用矩阵操作完成二维高斯拟合，并求取光斑中心下面的代码段是一个使用Eigen的矩阵操作完成二维高斯拟合求取光点的代码例子，关于二维高斯拟合求取光点的详细内容可参考：/hjx_1000/article/details/8490653cppview plaincopy1. boolGetCentrePoint(float&x0,float&y0)2. 3. if(m_iN=0)4. returnfalse;5. /QR分解6. HouseholderQRqr;7. pute(m_matrix_B);8. MatrixXfR=qr.matrixQR().triangularView();9. MatrixXfQ=qr.householderQ();10. 11. /块操作，获取向量或矩阵的局部12. VectorXfS;13. S=(Q.transpose()*m_Vector_A).head(5);14. MatrixXfR1;15. R1=R.block(0,0,5,5);16. 17. VectorXfC;18. C=R1.inverse()*S;19. 20. x0=-0.5*C1/C3;21. y0=-0.5*C2/C4;22. 23. returntrue;24. 1.首先，要把这个东东加到VS中，提供我们使用。先下载，解压缩： /index.php?title=Main_Page我的环境是 WIN7+VS2010 下载的3.2.2版本。解压缩以后有这个文件夹：eigen-eigen-1306d75b4a21恩，我只取了里面的 Eigen文件夹，放到了VS2010文件夹下的vs_2010文件夹下VC文件夹下的include文件夹中，也就是：盘符:VS2010vs_2010VCinclude然后就可以新建项目了，新建项目后，可以用：cppview plaincopyprint?1. #includeEigen/Eigen2. usingnamespaceEigen;加进来它的头文件，使用命名空间，and then 你就可以用它的函数了。2.认识它的一些头文件。Eigen这个类库，存的东西好多的，来看一下主要的几个头文件吧：最上面那段英文意思是：Eigen库分为 核心模块和额外模块两部分，每个模块都有一个用这个模块所相对应的头文件，Eigen和Dense头文件方便的同时包含了几个头文件以供使用。Core有关矩阵和数组的类，有基本的线性代数（包含 三角形 和 自伴乘积 相关），还有相应对数组的操作。Geometry几何学的类，有关转换、平移、进位制、2D旋转、3D旋转（四元组和角轴相关）LU逻辑单元的类，有关求逆，求行列式，LU分解解算器（FullPivLU,PartialPivLU)Cholesky包含LLT和LDLT的乔里斯基因式分解法。（小科普：Cholesky分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解）Householder豪斯霍尔德变换，这个模块供几个线性代数模块使用。（Householder transform: 维基百科）SVD奇异值分解，最小二乘解算器解决奇异值分解。QRQR分解求解，三种方法：HouseholderQR、ColPivHouseholderQR、FullPivHouseholderQREigenvalues特征值和特征向量分解的方法：EigenSolver、SelfAdjointEigenSolver、ComplexEigenSolverSparse稀疏矩阵相关类，对于稀疏矩阵的存储及相关基本线性代数Dense包含： Core、Gelometry、LU、Cholesky、SVD、QR和Eigenvalues模块（头文件）Eigen包含上述所有的模块（头文件）3.对矩阵的简单操作Eigen提供了两种密集的对象Matrix（矩阵）和Vector（向量）。这两者是通过矩阵模板类和一维或二维的数组模板类来实现的。这两者有几点不同：Matrix类型变量加减法，若行列数不相等，则不能做加减；Array类型的可以加减一个常数（各个元素分别加减该常数）。Matrix与Array类型变量做乘法也会有不同，Matrix是矩阵相乘，Array是对应元素相乘。但两者可以相互转换，方法为 .array() 和 .matrix()。定义（注意：定义矩阵时，默认没有初始化，必须自己初始化）Eigen的矩阵类型，一般是Matrix后面跟类型符号来表示，比如说： d 代表 double，矩阵存储的是double型的数据 f 代表float，矩阵存储的是float类型数据 c 代表complex，矩阵存数的是复数类型数据 i 代表int，矩阵存储的是整数类型相应关系为：比如：cppview plaincopyprint?1. MatrixXfm1(3,4);/建立3行4列的动态矩阵2. MatrixXfm2(3,3);3. Vector3fv1;/建立静态向量X代表动态，f代表float型初始化cppview plaincopyprint?1. m1=MatrixXf:Zero(3,4);/将矩阵3行4列初始化为02. m2=MatrixXf:Ones(3,3);/将矩阵3行3列初始化为13. v1=Vector3f:Ones();/将3行的纵向量初始化为14. 5. coutm1=nm1endl;6. coutm2=nm2endl;7. coutv1=nv1endl;运行出来效果：进一步，测试一下：cppview plaincopyprint?1. MatrixXfm3(4,5);2. m3=MatrixXf:Zero(4,5);3. coutm3_1=nm3endl;4. m3=MatrixXf:Ones(3,3);5. coutm3_2=nm3endl;6. m3=MatrixXf:Ones(6,6);7. coutm3_3=nm3endl;先定义一个，4行5列的矩阵初始化为0,4行5列，这时输出，就发现是4行5列的0初始化为3行3列的1，输出，3行3列的1初始化为6行6列的1，输出，则成为6行6列的1，具体看图：也就是说，矩阵的大小与初始化息息相关，初始化多少，它就是多少。谁让它是动态的呢？！那么，你肯定会说，定义的时候声明行列干啥？因为，下种方法初始化，就需要行列值了：cppview plaincopyprint?1. MatrixXfm3(2,3);2. m31,2,3,4,5,6;3. coutm3_1nendl;4. /为了美观点，更像个矩阵，可以换行写5. m31.3,4,-8,6. 0,0.9,2;7. coutm3_2=nm3endl;运行出来就是：访问这个就很简单了，直接就同数组的访问方式，但是不是方括号，而是圆括号：cppview plaincopyprint?1. MatrixXfm3(2,3);2. m31,2,3,4,5,6;3. coutm3_1nm3endl;4. /为了美观点，更像个矩阵，可以换行写5. m31.3,4,-8,6. 0,0.9,2;7. coutm3_2=nm3endl;8. /将第2行第3列的值改为999. m3(1,2)=99;10. coutm3_3=nm3endl;当然，同数组一样，第一行第一列的下标为（0,0）4.矩阵的基础运算代码执行了矩阵的：置0 置1 随机矩阵 单位阵 求逆 转置 数乘矩阵cppview plaincopyprint?1. MatrixXfm1(3,3);2. /矩阵全部元素置03. m1.setZero();4. coutm1_1=nm1endl;5. /矩阵全部元素置1(这里行列值不填，默认定义时候的行列，6. /若填写，则矩阵也会更改为填写的行列值)7. m1.setOnes(2,2);8. coutm1_2=nm1endl;9. /随机生成一个矩阵10. m1.setRandom();11. coutm1_3=nm1endl;12. /置单位矩阵13. m1.setIdentity(3,3);14. coutm1_4=nm1endl;15. m11,2,3,16. 5,9,10,17. 7,0,1;18. /矩阵求逆19. m1.inverse();20. coutm1_5=nm1endl;21. /矩阵转置22. m1.transpose();23. coutm1_6=nm1endl;24. /数*矩阵(数/矩阵)25. m1=2.6*m1;26. coutm1_7=nm1endl;运行结果： 啊，还有矩阵的加减乘除，这个不用写了吧？和平常的加减乘除一样的，就是矩阵乘法要注意，两个矩阵的行列值哟5.矩阵的高级运算Eigen 矩阵定义#include Matrix A; / Fixed rows and cols. Same as Matrix3d.Matrix B; / Fixed rows, dynamic cols.Matrix C; / Full dynamic. Same as MatrixXd.Matrix E; / Row major; default is column-major.Matrix3f P, Q, R; / 3x3 float matrix.Vector3f x, y, z; / 3x1 float matrix.RowVector3f a, b, c; / 1x3 float matrix.VectorXd v; / Dynamic column vector of doubles/ Eigen / Matlab / commentsx.size() / length(x) / vector sizeC.rows() / size(C,1) / number of rowsC.cols() / size(C,2) / number of columnsx(i) / x(i+1) / Mat