Eigen: C++开源矩阵计算库

　　Eigen库被分为一个Core模块和几个附加的模块，每个模块有一个相关的头文件，使用该模块时需要包含该头文件，为了能便利的使用eigen的几个模块，Eigen提供了Dense和Eigen两个头文件，各个头文件和模块如下表

Module	Header file	Contents
Core	#include <Eigen/Core>	Matrix and Array classes, basic linear algebra (including triangular and selfadjoint products), array manipulation
Geometry	#include <Eigen/Geometry>	Transform, Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis)
LU	#include <Eigen/LU>	Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU)
Cholesky	#include <Eigen/Cholesky>	LLT and LDLT Cholesky factorization with solver
Householder	#include <Eigen/Householder>	Householder transformations; this module is used by several linear algebra modules
SVD	#include <Eigen/SVD>	SVD decompositions with least-squares solver (JacobiSVD, BDCSVD)
QR	#include <Eigen/QR>	QR decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR)
Eigenvalues	#include <Eigen/Eigenvalues>	Eigenvalue, eigenvector decompositions (EigenSolver, SelfAdjointEigenSolver, ComplexEigenSolver)
Sparse	#include <Eigen/Sparse>	Sparse matrix storage and related basic linear algebra (SparseMatrix, SparseVector)  (see Quick reference guide for sparse matrices for details on sparse modules)
	#include <Eigen/Dense>	Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
	#include <Eigen/Eigen>	Includes Dense and Sparse header files (the whole Eigen library)

Eigen库分为核心模块和额外模块两个部分，

Eigen和Dense头文件方便的同时包含了几个头文件可以供使用

--Core

有关矩阵和数组的类，由基本的线性代数（包含 三角形和自伴乘积相关），还有相应对数组的操作。

--Geometry

几何学的类，有关转换、平移、进位制、2d旋转、3D旋转（四元组和角轴相关）

--LU

逻辑单元的类，有关求逆，求行列式，LU分解解算器（FullPivLU,PartialPivLU）

--Cholesky

包含LLT和LDLT的Cholesky因式分解法（Cholesky分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解）

--Householder

Householder变换，这个模块供几何线性代数模块使用

--SVD

奇异值分解，最小二乘解算器解决奇异值分解

--QR

QR分解求解，三种方法：HouseholderQR、ColPivHouseholderQR,FullPivhouseholderQR

--Eigenvalues

特征值和特征向量分解的方法：EigenSolver、SelfAdjointEigenSolver、ComplexEigenSolver

--Sparse

稀疏矩阵相关类，对于系数矩阵的存储及相关线性代数

--Dense

包含：core、Geometry、LU、Cholesky、SVD、QR和Eigenvalues模块（头文件）

--Eigen

包含上述所有模块（头文件）

1、矩阵的定义

Eigen中关于矩阵类的模板函数中，共有6个参数，但是常用的只有前三个，如下所示：

template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
class Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >

......

前三个参数分别表示矩阵元素的类型，行数和列数。

矩阵定义时可以使用Dynamic 来表示矩阵的行列数是未知，例如：

typedef matrix<double, Dynamic, Dynamic> MarixXd;

在Eigen中也提供了很多常见的简化定义形式，例如

typedef Matrix< double , 3 , 1> Vector3d

注意：

（1）Eigen中无论是矩阵还是数组、向量，无论是静态矩阵还是动态矩阵都提供默认构造函数，也就是你定义这些数据结构时都可以不用提供任何参数，其大小均由运行时来确定。

（2）矩阵的构造函数中只提供行列数、元素类型的构造参数，而不提供元素值的构造，对于比较小的、固定长度的向量提供初始化元素的定义，例如：

Vector2d a(5.0, 6.0);
Vector3d b(5.0, 6.0, 7.0);
Vector4d c(5.0, 6.0, 7.0, 8.0);

Vector2f v1(x, y);

Array3i   v2(x, y, z);
Vector4d  v3(x, y, z, w);
VectorXf  v5; // empty object
ArrayXf   v6(size);

//2D objects
atrix4f  m1;
MatrixXf  m5; // empty object
MatrixXf  m6(nb_rows, nb_columns);

   (3) Eigen矩阵的赋值方法

// 单点赋值
Eigen::Matrix<double, 3, 3> srcPoints;
    srcPoints(0, 0) = 4;
    srcPoints(1, 0) = 0;
    srcPoints(2, 0) = 0;
    srcPoints(0, 1) = 2;
//...

//<<运算符赋值
Vector3f  v1;     
v1 << x, y, z;

ArrayXf   v2(4);  
v2 << 1, 2, 3, 4;

//
Matrix3f  m1;   
m1 << 1, 2, 3,
      4, 5, 6,
      7, 8, 9;

int rows=5, cols=5;
MatrixXf m(rows,cols);
m << (Matrix3f() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished(),
     MatrixXf::Zero(3,cols-3),
     MatrixXf::Zero(rows-3,3),
     MatrixXf::Identity(rows-3,cols-3);
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
（3）矩阵的预定义：

//vector x
x.setZero();
x.setOnes();
x.setConstant(value);
x.setRandom();
x.setLinSpaced(size, low, high);

VectorXf::Unit(4,1) == Vector4f(0,1,0,0)
                    == Vector4f::UnitY()

//matrix x
x.setZero(rows, cols);
x.setOnes(rows, cols);
x.setConstant(rows, cols, value);
x.setRandom(rows, cols);

x.setIdentity();//单位矩阵
x.setIdentity(rows, cols);     

（4）映射操作

　　可以将c语言类型的数组映射为矩阵或向量：

　　注意：

　　1.映射只适用于对一维数组进行操作，如果希望将二维数组映射为对应的矩阵，可以借助"mat.row(i)=Map<Vector> v(data[i],n)"进行循环来实现，其中data为二维数组名，n为数组第一维的维度。

　　2.应设置之后数组名和矩阵或向量名其实指向同一个地址，只是名字和用法不同，另外，对其中的任何一个进行修改都会修改另外一个

//连续映射
float data[] = {1,2,3,4};
Map<Vector3f> v1(data);       // uses v1 as a Vector3f object
Map<ArrayXf>  v2(data,3);     // uses v2 as a ArrayXf object
Map<Array22f> m1(data);       // uses m1 as a Array22f object
Map<MatrixXf> m2(data,2,2);   // uses m2 as a MatrixXf object

//间隔映射
float data[] = {1,2,3,4,5,6,7,8,9};
Map<VectorXf,0,InnerStride<2> >  v1(data,3);                      // = [1,3,5]
Map<VectorXf,0,InnerStride<> >   v2(data,3,InnerStride<>(3));     // = [1,4,7]
Map<MatrixXf,0,OuterStride<3> >  m2(data,2,3);                    // both lines     |1,4,7|
Map<MatrixXf,0,OuterStride<> >   m1(data,2,3,OuterStride<>(3));   // are equal to:  |2,5,8|

2、动态矩阵和静态矩阵

动态矩阵是指其大小在运行时确定，静态矩阵是指其大小在编译时确定，在Eigen中并未这样称呼矩阵。具体可见如下两段代码：

代码段1：

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
MatrixXd m = MatrixXd::Random(3,3);
m = (m + MatrixXd::Constant(3,3,1.2)) * 50;
cout << "m =" << endl << m << endl;
VectorXd v(3);
v << 1, 2, 3;
cout << "m * v =" << endl << m * v << endl;
}

代码段2：

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
Matrix3d m = Matrix3d::Random();
m = (m + Matrix3d::Constant(1.2)) * 50;
cout << "m =" << endl << m << endl;
Vector3d v(1,2,3);
cout << "m * v =" << endl << m * v << endl;
}

说明：

1）代码段1中MatrixXd表示任意大小的元素类型为double的矩阵变量，其大小只有在运行时被赋值之后才能知道； MatrixXd::Random(3,3)表示产生一个元素类型为double的3*3的临时矩阵对象。

 2） 代码段2中Matrix3d表示元素类型为double大小为3*3的矩阵变量，其大小在编译时就知道；

3）上例中向量的定义也是类似，不过这里的向量时列优先，在Eigen中行优先的矩阵会在其名字中包含有row，否则就是列优先。

4）向量只是一个特殊的矩阵，其一个维度为1而已，如：typedef Matrix< double , 3 , 1> Vector3d

3、矩阵元素的访问

 在矩阵的访问中，行索引总是作为第一个参数，需注意Eigen中遵循大家的习惯让矩阵、数组、向量的下标都是从0开始。矩阵元素的访问可以通过（）操作符完成，例如m(2,3)即是获取矩阵m的第2行第3列元素（注意行列数从0开始）。可参看如下代码：

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
MatrixXd m(2,2);
m(0,0) = 3;
m(1,0) = 2.5;
m(0,1) = -1;
m(1,1) = m(1,0) + m(0,1);
std::cout << "Here is the matrix m:
" << m << std::endl;
VectorXd v(2);
v(0) = 4;
v(1) = v(0) - 1;
std::cout << "Here is the vector v:
" << v << std::endl;
}

其输出结果为：

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

针对向量还提供[]操作符，注意矩阵则不可如此使用，原因为：在C++中m[i, j]中逗号表达式 “i, j”的值始终都是“j”的值，即m[i, j]对于C++来讲就是m[j];

4、设置矩阵的元素 

在Eigen中重载了"<<"操作符，通过该操作符即可以一个一个元素的进行赋值，也可以一块一块的赋值。另外也可以使用下标进行复制，例如下面两段代码：

 代码段1:

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

输出结果为：

1 2 3
4 5 6
7 8 9

代码段二（使用下标进行复制）:

VectorXf m_Vector_A;
MatrixXf m_matrix_B;
int m_iN =-1;
 
bool InitData(int pSrc[100][100], int iWidth, int iHeight)
{
    if (NULL == pSrc || iWidth <=0 || iHeight <= 0)
        return false;
    m_iN = iWidth*iHeight;
    VectorXf tmp_A(m_iN);
    MatrixXf tmp_B(m_iN, 5);
    int i =0, j=0, iPos =0;
    while(i<iWidth)
    {
         j=0;
        while(j<iHeight)
        {
            tmp_A(iPos) = pSrc[i][j] * log((float)pSrc[i][j]);
 
            tmp_B(iPos,0) = pSrc[i][j] ;
            tmp_B(iPos,1) = pSrc[i][j] * i;
            tmp_B(iPos,2) = pSrc[i][j] * j;
            tmp_B(iPos,3) = pSrc[i][j] * i * i;
            tmp_B(iPos,4) = pSrc[i][j] * j * j;
            ++iPos;
            ++j;
        }
        ++i;
    }
    m_Vector_A = tmp_A;
    m_matrix_B = tmp_B;
}

5、重置矩阵大小

当前矩阵的行数、列数、大小可以通过rows()，cols()和size()来获取，对于动态矩阵可以通过resize()函数来动态修改矩阵的大小.

需注意：

(1) 固定大小的矩阵是不能使用resize（）来修改矩阵的大小；

(2) resize（）函数会析构掉原来的数据，因此调用resize（）函数之后将不能保证元素的值不改变。
(3) 使用“=”操作符操作动态矩阵时，如果左右边的矩阵大小不等，则左边的动态矩阵的大小会被修改为右边的大小。例如下面的代码段：

MatrixXf a(2,2);
std::cout << "a is of size " << a.rows() << "x" << a.cols() << std::endl;
MatrixXf b(3,3);
a = b;
std::cout << "a is now of size " << a.rows() << "x" << a.cols() << std::endl;

//
vector.resize(size);
vector.resizeLike(other_vector);
vector.conservativeResize(size);

//
matrix.resize(nb_rows, nb_cols);
matrix.resize(Eigen::NoChange, nb_cols);
matrix.resize(nb_rows, Eigen::NoChange);
matrix.resizeLike(other_matrix);
matrix.conservativeResize(nb_rows, nb_cols);

输出结果为：

a is of size 2x2
a is now of size 3x3

6、如何选择动态矩阵和静态矩阵

Eigen对于这问题的答案是：对于小矩阵（一般大小小于16）的使用固定大小的静态矩阵，它可以带来比较高的效率，对于大矩阵（一般大小大于32）建议使用动态矩阵。

还需特别注意的是：如果特别大的矩阵使用了固定大小的静态矩阵则可能造成栈溢出的问题.

7、矩阵的运算

Eigen重载了基础的+ - * / += -= *= /= *可以表示标量和矩阵或者矩阵和矩阵

矩阵与矩阵的运算

Eigen提供+、-、一元操作符“-”、+=、-=，例如：

二元操作符+/-表示两矩阵相加（矩阵中对应元素相加/减，返回一个临时矩阵）： B+C 或 B-C；

一元操作符-表示对矩阵取负（矩阵中对应元素取负，返回一个临时矩阵）： -C; 

组合操作法+=或者-=表示（对应每隔元素都做相应操作）：A += B 或者 A-=B

代码段1为矩阵的加减操作，代码如下：

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
a << 1, 2,
3, 4;
MatrixXd b(2,2);
b << 2, 3,
1, 4;
std::cout << "a + b =
" << a + b << std::endl;
std::cout << "a - b =
" << a - b << std::endl;
std::cout << "Doing a += b;" << std::endl;
a += b;
std::cout << "Now a =
" << a << std::endl;
Vector3d v(1,2,3);
Vector3d w(1,0,0);
std::cout << "-v + w - v =
" << -v + w - v << std::endl;
}

输出结果为：

a + b =
3 5
4 8
a - b =
-1 -1
 2  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6

矩阵与标量的运算

矩阵还提供与标量（单一个数字）的乘除操作，表示每个元素都与该标量进行乘除操作。例如：

二元操作符*在：A*a中表示矩阵A中的每个元素都与数字a相乘，结果放在一个临时矩阵中，矩阵的值不会改变。

对于a*A、A/a、A*=a、A /=a也是一样，例如下面的代码：

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
a << 1, 2,
3, 4;
Vector3d v(1,2,3);
std::cout << "a * 2.5 =
" << a * 2.5 << std::endl;
std::cout << "0.1 * v =
" << 0.1 * v << std::endl;
std::cout << "Doing v *= 2;" << std::endl;
v *= 2;
std::cout << "Now v =
" << v << std::endl;
}

输出结果为：

a * 2.5 =
2.5  5
7.5 10
0.1 * v =
0.1
0.2
0.3
Doing v *= 2;
Now v =
2
4
6

需要注意：

在Eigen中，算术操作例如 “操作符+”并不会自己执行计算操作，他们只是返回一个“算术表达式对象”，而实际的计算则会延迟到后面的赋值时才进行。这些不影响你的使用，它只是为了方便Eigen的优化

算术运算方法：

add  subtract	mat3 = mat1 + mat2;           mat3 += mat1; mat3 = mat1 - mat2;           mat3 -= mat1;
scalar product	mat3 = mat1 * s1;             mat3 *= s1;           mat3 = s1 * mat1; mat3 = mat1 / s1;             mat3 /= s1;
matrix/vector  products *	col2 = mat1 * col1; row2 = row1 * mat1;           row1 *= mat1; mat3 = mat1 * mat2;           mat3 *= mat1;
transposition  adjoint *	mat1 = mat2.transpose();      mat1.transposeInPlace(); mat1 = mat2.adjoint();        mat1.adjointInPlace();
dot product  inner product *	scalar = vec1.dot(vec2); scalar = col1.adjoint() * col2; scalar = (col1.adjoint() * col2).value();
outer product *	mat = col1 * col2.transpose();
norm  normalization *	scalar = vec1.norm();         scalar = vec1.squaredNorm() vec2 = vec1.normalized();     vec1.normalize(); // inplace
cross product *	#include <Eigen/Geometry> vec3 = vec1.cross(vec2);

Coefficient-wise & Array operators

Coefficient-wise operators for matrices and vectors:

Matrix API *	Via Array conversions
mat1.cwiseMin(mat2) mat1.cwiseMax(mat2) mat1.cwiseAbs2() mat1.cwiseAbs() mat1.cwiseSqrt() mat1.cwiseProduct(mat2) mat1.cwiseQuotient(mat2)	mat1.array().min(mat2.array()) mat1.array().max(mat2.array()) mat1.array().abs2() mat1.array().abs() mat1.array().sqrt() mat1.array() * mat2.array() mat1.array() / mat2.array()

Array operators:*

Arithmetic operators	array1 * array2     array1 / array2     array1 *= array2    array1 /= array2 array1 + scalar     array1 - scalar     array1 += scalar    array1 -= scalar
Comparisons	array1 < array2     array1 > array2     array1 < scalar     array1 > scalar array1 <= array2    array1 >= array2    array1 <= scalar    array1 >= scalar array1 == array2    array1 != array2    array1 == scalar    array1 != scalar
Trigo, power, and  misc functions  and the STL variants	array1.min(array2)            array1.max(array2)            array1.abs2() array1.abs()                  abs(array1) array1.sqrt()                 sqrt(array1) array1.log()                  log(array1) array1.exp()                  exp(array1) array1.pow(exponent)          pow(array1,exponent) array1.square() array1.cube() array1.inverse() array1.sin()                  sin(array1) array1.cos()                  cos(array1) array1.tan()                  tan(array1) array1.asin()                 asin(array1) array1.acos()                 acos(array1)

Reductions

Eigen provides several reduction methods such as: minCoeff() , maxCoeff() , sum() , prod() , trace() *, norm()*, squaredNorm() *, all() , and any() . All reduction operations can be done matrix-wise, column-wise or row-wise. Usage example:

5 3 1 mat = 2 7 8       9 4 6	mat.minCoeff();	1
	mat.colwise().minCoeff();	2 3 1
	mat.rowwise().minCoeff();	1 2 4

Typical use cases of all() and any():

//Typical use cases of all() and any():

if((array1 > 0).all()) ...      // if all coefficients of array1 are greater than 0 ...
if((array1 < array2).any()) ... // if there exist a pair i,j such that array1(i,j) < array2(i,j) ...

c++矩阵运算库Eigen

　　Eigen 的简介和下载安装

最近因为要写机械臂的运动学控制程序，因此难免就要在C++中进行矩阵运算。然而C++本身不支持矩阵运算，Qt库中的矩阵运算能力也相当有限，因此考虑寻求矩阵运算库Eigen的帮助。　

Eigen是一个C++线性运算的模板库：他可以用来完成矩阵，向量，数值解等相关的运算。（Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.）

Eigen库下载： http://eigen.tuxfamily.org/index.php?title=Main_Page

Eigen库的使用相当方便，将压缩包中的Eigen文件夹拷贝到项目目录下，直接包含其中的头文件即可使用，省去了使用Cmake进行编译的烦恼。

Eigen官网有快速入门的参考文档：http://eigen.tuxfamily.org/dox/group__QuickRefPage.html

　　Eigen简单上手使用

要实现相应的功能只需要包含头相应的头文件即可：

Core	#include <Eigen/Core>	Matrix and Array classes, basic linear algebra (including triangular and selfadjoint products), array manipulation
Geometry	#include <Eigen/Geometry>	Transform, Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis)
LU	#include <Eigen/LU>	Inverse, determinant, LU decompositions with solver (FullPivLU, PartialPivLU)
Cholesky	#include <Eigen/Cholesky>	LLT and LDLT Cholesky factorization with solver
Householder	#include <Eigen/Householder>	Householder transformations; this module is used by several linear algebra modules
SVD	#include <Eigen/SVD>	SVD decomposition with least-squares solver (JacobiSVD)
QR	#include <Eigen/QR>	QR decomposition with solver (HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR)
Eigenvalues	#include <Eigen/Eigenvalues>	Eigenvalue, eigenvector decompositions (EigenSolver, SelfAdjointEigenSolver,ComplexEigenSolver)
Sparse	#include <Eigen/Sparse>	Sparse matrix storage and related basic linear algebra (SparseMatrix, DynamicSparseMatrix,SparseVector)
	#include <Eigen/Dense>	Includes Core, Geometry, LU, Cholesky, SVD, QR, and Eigenvalues header files
	#include <Eigen/Eigen>	Includes Dense and Sparse header files (the whole Eigen library)

 基本的矩阵运算只需要包含Dense即可

基本的数据类型：Array, matrix and vector

声明：

#include<Eigen/Dense>

...
//1D objects
Vector4d  v4;
Vector2f  v1(x, y);
Array3i   v2(x, y, z);
Vector4d  v3(x, y, z, w);
VectorXf  v5; // empty object
ArrayXf   v6(size);

//2D objects
atrix4f  m1;
MatrixXf  m5; // empty object
MatrixXf  m6(nb_rows, nb_columns);

赋值：

//    
Vector3f  v1;     v1 << x, y, z;
ArrayXf   v2(4);  v2 << 1, 2, 3, 4;

//
Matrix3f  m1;   m1 << 1, 2, 3,
                      4, 5, 6,
                      7, 8, 9;

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

Resize矩阵：

//
vector.resize(size);
vector.resizeLike(other_vector);
vector.conservativeResize(size);

//
matrix.resize(nb_rows, nb_cols);
matrix.resize(Eigen::NoChange, nb_cols);
matrix.resize(nb_rows, Eigen::NoChange);
matrix.resizeLike(other_matrix);
matrix.conservativeResize(nb_rows, nb_cols);

元素读取：

vector(i)；    
vector[i]；

matrix(i,j)；

矩阵的预定义：

//vector x
x.setZero();
x.setOnes();
x.setConstant(value);
x.setRandom();
x.setLinSpaced(size, low, high);

VectorXf::Unit(4,1) == Vector4f(0,1,0,0)
                    == Vector4f::UnitY()

//matrix x
x.setZero(rows, cols);
x.setOnes(rows, cols);
x.setConstant(rows, cols, value);
x.setRandom(rows, cols);

x.setIdentity();//单位矩阵
x.setIdentity(rows, cols);     

映射操作，可以将c语言类型的数组映射为矩阵或向量：

（注意：

1.映射只适用于对一维数组进行操作，如果希望将二维数组映射为对应的矩阵，可以借助"mat.row(i)=Map<Vector> v(data[i],n)"进行循环来实现，其中data为二维数组名，n为数组第一维的维度。

2.应设置之后数组名和矩阵或向量名其实指向同一个地址，只是名字和用法不同，另外，对其中的任何一个进行修改都会修改另外一个）

//连续映射
float data[] = {1,2,3,4};
Map<Vector3f> v1(data);       // uses v1 as a Vector3f object
Map<ArrayXf>  v2(data,3);     // uses v2 as a ArrayXf object
Map<Array22f> m1(data);       // uses m1 as a Array22f object
Map<MatrixXf> m2(data,2,2);   // uses m2 as a MatrixXf object

//间隔映射
float data[] = {1,2,3,4,5,6,7,8,9};
Map<VectorXf,0,InnerStride<2> >  v1(data,3);                      // = [1,3,5]
Map<VectorXf,0,InnerStride<> >   v2(data,3,InnerStride<>(3));     // = [1,4,7]
Map<MatrixXf,0,OuterStride<3> >  m2(data,2,3);                    // both lines     |1,4,7|
Map<MatrixXf,0,OuterStride<> >   m1(data,2,3,OuterStride<>(3));   // are equal to:  |2,5,8|

算术运算：

add  subtract	mat3 = mat1 + mat2;           mat3 += mat1; mat3 = mat1 - mat2;           mat3 -= mat1;
scalar product	mat3 = mat1 * s1;             mat3 *= s1;           mat3 = s1 * mat1; mat3 = mat1 / s1;             mat3 /= s1;
matrix/vector  products *	col2 = mat1 * col1; row2 = row1 * mat1;           row1 *= mat1; mat3 = mat1 * mat2;           mat3 *= mat1;
transposition  adjoint *	mat1 = mat2.transpose();      mat1.transposeInPlace(); mat1 = mat2.adjoint();        mat1.adjointInPlace();
dot product  inner product *	scalar = vec1.dot(vec2); scalar = col1.adjoint() * col2; scalar = (col1.adjoint() * col2).value();
outer product *	mat = col1 * col2.transpose();
norm  normalization *	scalar = vec1.norm();         scalar = vec1.squaredNorm() vec2 = vec1.normalized();     vec1.normalize(); // inplace
cross product *	#include <Eigen/Geometry> vec3 = vec1.cross(vec2);

Coefficient-wise & Array operators

Coefficient-wise operators for matrices and vectors:

Matrix API *	Via Array conversions
mat1.cwiseMin(mat2) mat1.cwiseMax(mat2) mat1.cwiseAbs2() mat1.cwiseAbs() mat1.cwiseSqrt() mat1.cwiseProduct(mat2) mat1.cwiseQuotient(mat2)	mat1.array().min(mat2.array()) mat1.array().max(mat2.array()) mat1.array().abs2() mat1.array().abs() mat1.array().sqrt() mat1.array() * mat2.array() mat1.array() / mat2.array()

Array operators:*

Arithmetic operators	array1 * array2     array1 / array2     array1 *= array2    array1 /= array2 array1 + scalar     array1 - scalar     array1 += scalar    array1 -= scalar
Comparisons	array1 < array2     array1 > array2     array1 < scalar     array1 > scalar array1 <= array2    array1 >= array2    array1 <= scalar    array1 >= scalar array1 == array2    array1 != array2    array1 == scalar    array1 != scalar
Trigo, power, and  misc functions  and the STL variants	array1.min(array2)            array1.max(array2)            array1.abs2() array1.abs()                  abs(array1) array1.sqrt()                 sqrt(array1) array1.log()                  log(array1) array1.exp()                  exp(array1) array1.pow(exponent)          pow(array1,exponent) array1.square() array1.cube() array1.inverse() array1.sin()                  sin(array1) array1.cos()                  cos(array1) array1.tan()                  tan(array1) array1.asin()                 asin(array1) array1.acos()                 acos(array1)

Reductions

Eigen provides several reduction methods such as: minCoeff() , maxCoeff() , sum() , prod() , trace() *, norm()*, squaredNorm() *, all() , and any() . All reduction operations can be done matrix-wise, column-wise or row-wise. Usage example:

5 3 1 mat = 2 7 8       9 4 6	mat.minCoeff();	1
	mat.colwise().minCoeff();	2 3 1
	mat.rowwise().minCoeff();	1 2 4

Typical use cases of all() and any():

//Typical use cases of all() and any():

if((array1 > 0).all()) ...      // if all coefficients of array1 are greater than 0 ...
if((array1 < array2).any()) ... // if there exist a pair i,j such that array1(i,j) < array2(i,j) ...

 

Sub-matrices

Read-write access to a column or a row of a matrix (or array):

mat1.row(i) = mat2.col(j);

mat1.col(j1).swap(mat1.col(j2));

Read-write access to sub-vectors:

Default versions	Optimized versions when the size  is known at compile time
vec1.head(n)	vec1.head<n>()	the first n coeffs
vec1.tail(n)	vec1.tail<n>()	the last n coeffs
vec1.segment(pos,n)	vec1.segment<n>(pos)	the n coeffs in the  range [pos : pos + n - 1]
Read-write access to sub-matrices:
mat1.block(i,j,rows,cols) (more)	mat1.block<rows,cols>(i,j) (more)	the rows x cols sub-matrix  starting from position (i,j)
mat1.topLeftCorner(rows,cols) mat1.topRightCorner(rows,cols) mat1.bottomLeftCorner(rows,cols) mat1.bottomRightCorner(rows,cols)	mat1.topLeftCorner<rows,cols>() mat1.topRightCorner<rows,cols>() mat1.bottomLeftCorner<rows,cols>() mat1.bottomRightCorner<rows,cols>()	the rows x cols sub-matrix  taken in one of the four corners
mat1.topRows(rows) mat1.bottomRows(rows) mat1.leftCols(cols) mat1.rightCols(cols)	mat1.topRows<rows>() mat1.bottomRows<rows>() mat1.leftCols<cols>() mat1.rightCols<cols>()	specialized versions of block()  when the block fit two corners

top

Miscellaneous operations

Reverse

Vectors, rows, and/or columns of a matrix can be reversed (see DenseBase::reverse(), DenseBase::reverseInPlace(),VectorwiseOp::reverse()).

vec.reverse()           mat.colwise().reverse()   mat.rowwise().reverse()

vec.reverseInPlace()

Replicate

Vectors, matrices, rows, and/or columns can be replicated in any direction (see DenseBase::replicate(),VectorwiseOp::replicate())

vec.replicate(times)                                          vec.replicate<Times>

mat.replicate(vertical_times, horizontal_times)               mat.replicate<VerticalTimes, HorizontalTimes>()

mat.colwise().replicate(vertical_times, horizontal_times)     mat.colwise().replicate<VerticalTimes, HorizontalTimes>()

mat.rowwise().replicate(vertical_times, horizontal_times)     mat.rowwise().replicate<VerticalTimes, HorizontalTimes>()

top

Diagonal, Triangular, and Self-adjoint matrices

(matrix world *)

Diagonal matrices

Operation	Code
view a vector as a diagonal matrix	mat1 = vec1.asDiagonal();
Declare a diagonal matrix	DiagonalMatrix<Scalar,SizeAtCompileTime> diag1(size); diag1.diagonal() = vector;
Access the diagonal and super/sub diagonals of a matrix as a vector (read/write)	vec1 = mat1.diagonal();        mat1.diagonal() = vec1;      // main diagonal vec1 = mat1.diagonal(+n);      mat1.diagonal(+n) = vec1;    // n-th super diagonal vec1 = mat1.diagonal(-n);      mat1.diagonal(-n) = vec1;    // n-th sub diagonal vec1 = mat1.diagonal<1>();     mat1.diagonal<1>() = vec1;   // first super diagonal vec1 = mat1.diagonal<-2>();    mat1.diagonal<-2>() = vec1;  // second sub diagonal
Optimized products and inverse	mat3  = scalar * diag1 * mat1; mat3 += scalar * mat1 * vec1.asDiagonal(); mat3 = vec1.asDiagonal().inverse() * mat1 mat3 = mat1 * diag1.inverse()

Triangular views

TriangularView gives a view on a triangular part of a dense matrix and allows to perform optimized operations on it. The opposite triangular part is never referenced and can be used to store other information.

Note

The .triangularView() template member function requires the template keyword if it is used on an object of a type that depends on a template parameter; see The template and typename keywords in C++ for details.

Operation	Code
Reference to a triangular with optional  unit or null diagonal (read/write):	m.triangularView<Xxx>() Xxx = Upper, Lower, StrictlyUpper, StrictlyLower, UnitUpper,UnitLower
Writing to a specific triangular part: (only the referenced triangular part is evaluated)	m1.triangularView<Eigen::Lower>() = m2 + m3
Conversion to a dense matrix setting the opposite triangular part to zero:	m2 = m1.triangularView<Eigen::UnitUpper>()
Products:	m3 += s1 * m1.adjoint().triangularView<Eigen::UnitUpper>() * m2 m3 -= s1 * m2.conjugate() * m1.adjoint().triangularView<Eigen::Lower>()
Solving linear equations:	L1.triangularView<Eigen::UnitLower>().solveInPlace(M2) L1.triangularView<Eigen::Lower>().adjoint().solveInPlace(M3) U1.triangularView<Eigen::Upper>().solveInPlace<OnTheRight>(M4)

8、求矩阵的转置、共轭矩阵、伴随矩阵、行列式、逆矩阵。

可以通过 成员函数transpose(), conjugate(),determinant(), adjoint()和inverse()来完成，注意这些函数返回操作后的结果，而不会对原矩阵的元素进行直接操作，如果要让原矩阵的进行转换，则需要使用响应的InPlace函数，例如：transposeInPlace() 、 adjointInPlace() 之类。

例如下面的代码所示：

MatrixXcf a = MatrixXcf::Random(2,2);
cout << "Here is the matrix a
" << a << endl;
cout << "Here is the matrix a^T
" << a.transpose() << endl;
cout << "Here is the conjugate of a
" << a.conjugate() << endl;
cout << "Here is the matrix a^*
" << a.adjoint() << endl;
cout << "determinant： " << c.determinant() << endl;
cout << "inverse:" << c.inverse() << 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 a^T
(-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)
......

矩阵的特征值: m.eigenvalues();

特征向量: m.eigenvector();

9、矩阵相乘、矩阵向量相乘

矩阵的相乘，矩阵与向量的相乘也是使用操作符*，共有*和*=两种操作符，其用法可以参考如下代码：

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d mat;
mat << 1, 2,
3, 4;
Vector2d u(-1,1), v(2,0);
std::cout << "Here is mat*mat:
" << mat*mat << std::endl;
std::cout << "Here is mat*u:
" << mat*u << std::endl;
std::cout << "Here is u^T*mat:
" << u.transpose()*mat << std::endl;
std::cout << "Here is u^T*v:
" << u.transpose()*v << std::endl;
std::cout << "Here is u*v^T:
" << u*v.transpose() << std::endl;
std::cout << "Let's multiply mat by itself" << std::endl;
mat = mat*mat;
std::cout << "Now mat is mat:
" << mat << std::endl;
}

输出结果为：

Here is mat*mat:
 7 10
15 22
Here is mat*u:
1
1
Here is u^T*mat:
2 2
Here is u^T*v:
-2
Here is u*v^T:
-2 -0
 2  0
Let's multiply mat by itself
Now mat is mat:
 7 10
15 22

10、点积和叉积

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
    //点积、叉积(针对向量的)
    Vector3d v(1,2,3);
    Vector3d w(0,1,2);
    std::cout<<v.dot(w)<<std::endl<<std::endl;
    std::cout<<w.cross(v)<<std::endl<<std::endl;
}
*/

11、矩阵的块操作

   1）矩阵的块操作有两种使用方法，其定义形式为：

matrix.block(i,j,p,q);      （1）
 
matrix.block<p,q>(i,j);    （2）

定义（1）表示返回从矩阵的(i, j)开始，每行取p个元素，每列取q个元素所组成的临时新矩阵对象，原矩阵的元素不变。

定义（2）中block（p, q）可理解为一个p行q列的子矩阵，该定义表示从原矩阵中第(i, j)开始，获取一个p行q列的子矩阵，返回该子矩阵组成的临时 矩阵对象，原矩阵的元素不变。

详细使用情况，可参考下面的代码段：

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::MatrixXf m(4,4);
m << 1, 2, 3, 4,
5, 6, 7, 8,
9,10,11,12,
13,14,15,16;
cout << "Block in the middle" << endl;
cout << m.block<2,2>(1,1) << endl << endl;
for (int i = 1; i <= 3; ++i)
{
cout << "Block of size " << i << "x" << i << endl;
cout << m.block(0,0,i,i) << endl << endl;
}
}

输出的结果为：

Block in the middle
6  7
10 11
Block of size 1x1
1
Block of size 2x2
1 2
5 6
Block of size 3x3
1  2  3
5  6  7
9 10 11

通过上述方式获取的子矩阵即可以作为左值也可以作为右值，也就是即可以用这个子矩阵给其他矩阵赋值，也可以给这个子矩阵对象赋值。

2）矩阵也提供了获取其指定行/列的函数，其实获取某行/列也是一种特殊的获取子块。可以通过 .col()和 .row()来完成获取指定列/行的操作，参数为列/行的索引。
注意：
（1）需与获取矩阵的行数/列数的函数（ rows(), cols() ）的进行区别，不要弄混淆。
（2）函数参数为响应行/列的索引，需注意矩阵的行列均以0开始。
下面的代码段用于演示获取矩阵的指定行列：

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::MatrixXf m(3,3);
m << 1,2,3,
4,5,6,
7,8,9;
cout << "Here is the matrix m:" << endl << m << endl;
cout << "2nd Row: " << m.row(1) << endl;
m.col(2) += 3 * m.col(0);
cout << "After adding 3 times the first column into the third column, the matrix m is:
";
cout << m << endl;
}

输出结果为：

Here is the matrix m:
1 2 3
4 5 6
7 8 9
2nd Row: 4 5 6
After adding 3 times the first column into the third column, the matrix m is:
 1  2  6
 4  5 18
 7  8 30

3）向量的块操作，其实向量只是一个特殊的矩阵，但是Eigen也为它单独提供了一些简化的块操作，如下三种形式：
     获取向量的前n个元素：vector.head(n); 
     获取向量尾部的n个元素：vector.tail(n);
     获取从向量的第i个元素开始的n个元素：vector.segment(i,n);
     其用法可参考如下代码段：

#include <Eigen/Dense>
#include <iostream>
using namespace std;
int main()
{
Eigen::ArrayXf v(6);
v << 1, 2, 3, 4, 5, 6;
cout << "v.head(3) =" << endl << v.head(3) << endl << endl;
cout << "v.tail<3>() = " << endl << v.tail<3>() << endl << endl;
v.segment(1,4) *= 2;
cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
}

输出结果为：

v.head(3) =
1
2
3
v.tail<3>() = 
4
5
6
after 'v.segment(1,4) *= 2', v =
1
4
6
8
10
6

12、计算特征值和特征向量

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    //特征向量、特征值
    std::cout << "Here is the matrix A:
" << a << std::endl;
    SelfAdjointEigenSolver<Matrix2d> eigensolver(a);
    if (eigensolver.info() != Success) abort();
     std::cout << "特征值:
" << eigensolver.eigenvalues() << std::endl;
     std::cout << "Here's a matrix whose columns are eigenvectors of A 
"
     << "corresponding to these eigenvalues:
"
     << eigensolver.eigenvectors() << std::endl;
}

13、矩阵分解

矩阵分解 (decomposition, factorization)是将矩阵拆解为数个矩阵的乘积，可分为三角分解、满秩分解、QR分解、Jordan分解和SVD（奇异值）分解等，常见的有三种：1)三角分解法 (Triangular Factorization)，2)QR 分解法 (QR Factorization)，3)奇异值分解法 (Singular Value Decompostion)。

Eigen总共提供了下面这些矩阵的分解方式：

Decomposition	Method	Requirements on the matrix	Speed	Accuracy
PartialPivLU	partialPivLu()	Invertible(可逆)	++	+
FullPivLU	fullPivLu()	None	-	+++
HouseholderQR	householderQr()	None	++	+
ColPivHouseholderQR	colPivHouseholderQr()	None	+	++
FullPivHouseholderQR	fullPivHouseholderQr()	None	-	+++
LLT	llt()	Positive definite(正定)	+++	+
LDLT	ldlt()	Positive or negative semidefinite(正或负半定)	+++	++

QR分解

   QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵,所以称为QR分解法,与此正规正交矩阵的通用符号Q有关。

householderQR的分解方式的演示代码：

void QR2()
{
    Matrix3d A;
    A<<1,1,1,
       2,-1,-1,
       2,-4,5;

    HouseholderQR<Matrix3d> qr;
    qr.compute(A);
    MatrixXd R = qr.matrixQR().triangularView<Upper>();
    MatrixXd Q =  qr.householderQ();
    std::cout << "QR2(): HouseholderQR---------------------------------------------"<< std::endl;
    std::cout << "A "<< std::endl <<A << std::endl << std::endl;
    std::cout <<"qr.matrixQR()"<< std::endl << qr.matrixQR() << std::endl << std::endl;
    std::cout << "R"<< std::endl <<R << std::endl << std::endl;
    std::cout << "Q "<< std::endl <<Q << std::endl << std::endl;
    std::cout <<"Q*R" << std::endl <<Q*R << std::endl << std::endl;
}

输出为：

QR2(): HouseholderQR---------------------------------------------
A 
 1  1  1
 2 -1 -1
 2 -4  5

qr.matrixQR()
 -3   3  -3
0.5  -3   3
0.5  -1  -3

R
-3  3 -3
 0 -3  3
 0  0 -3

Q 
-0.333333 -0.666667 -0.666667
-0.666667 -0.333333  0.666667
-0.666667  0.666667 -0.333333

Q*R
 1  1  1
 2 -1 -1
 2 -4  5

矩阵的LU三角分解 

矩阵的SVD（奇异值分解）分解

奇异值分解 (singular value decomposition,SVD) 是另一种正交矩阵分解法；SVD是最可靠的分解法，但是它比QR 分解法要花上近十倍的计算时间。[U,S,V]=svd(A)，其中U和V分别代表两个正交矩阵，而S代表一对角矩阵。 和QR分解法相同， 原矩阵A不必为正方矩阵。使用SVD分解法的用途是解最小平方误差法和数据压缩。
 

矩阵的LL^T分解 

Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的（LU三角分解法的变形）。

LDLT分解

对于对称正定矩阵A，可以将A分解为A=LDL^T，并且该分解是唯一的，其中L为一下三角形单位矩阵（即主对角线元素皆为1），D为一对角矩阵（只在主对角线上有元素且不为0，其余皆为零），L^T为L的转置矩阵。LDL^T 是Cholesky分解的变形：

LDL^T分解法实际上是Cholesky分解法的改进，因为Cholesky分解法虽然不需要选主元，但其运算过程中涉及到开方问题，而LDL^T分解法则避免了这一问题，可用于求解线性方程组。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    //线性方程求解 Ax =B;
    Matrix4d A;
    A << 2,-1,-1,1,
        1,1,-2,1,
        4,-6,2,-2,
        3,6,-9,7;

    Vector4d b(2,4,4,9);

    Vector4d x1 = A.colPivHouseholderQr().solve(b);
   

Vector4d x2

 = A.ldlt().solve(b));  // A sym. p.s.d.    #include <Eigen/Cholesky>

Vector4d x3

 = A.llt().solve(b));  // A sym. p.d.      #include <Eigen/Cholesky>

Vector4d x4

 = A.lu().solve(b));  // Stable and fast. #include <Eigen/LU>

Vector4d x5

 = A.qr().solve(b));  // No pivoting.     #include <Eigen/QR>

Vector4d x6

 = A.svd().solve(b));  // Stable, slowest. #include <Eigen/SVD>

// .ldlt() -> .matrixL() and .matrixD()
// .llt()  -> .matrixL()
// .lu()   -> .matrixL() and .matrixU()
// .qr()   -> .matrixQ() and .matrixR()
// .svd()  -> .matrixU(), .singularValues(), and .matrixV()

    std::cout << "The solution is:
" << x <<"

"<<x2<<"

"<<x3 <<std::endl;
}

输出为： 

 0
-1
-4
-3

-289.143
 448.714
 29.9082
 3.97959

  1.52903
   0.1758
-0.340206
0.0423223

结果不一致，这是为什么？

14、特征值分解

最小二乘法求解

最小二乘求解有两种方式，jacobiSvd或者colPivHouseholderQr，4*4以下的小矩阵速度没有区别，jacobiSvd可能更快，大矩阵最好用colPivHouseholderQr

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    MatrixXf A1 = MatrixXf::Random(3, 2);
    std::cout << "Here is the matrix A:
" << A1 << std::endl;
    VectorXf b1 = VectorXf::Random(3);
    std::cout << "Here is the right hand side b:
" << b1 << std::endl;
    //jacobiSvd 方式:Slow (but fast for small matrices)
    std::cout << "The least-squares solution is:
"
    << A1.jacobiSvd(ComputeThinU | ComputeThinV).solve(b1) << std::endl;
    //colPivHouseholderQr方法:fast
    std::cout << "The least-squares solution is:
"
    << A1.colPivHouseholderQr().solve(b1) << std::endl;
}

输出为：

Here is the matrix A:
 0.680375   0.59688
-0.211234  0.823295
 0.566198 -0.604897
Here is the right hand side b:
-0.329554
 0.536459
-0.444451
The least-squares solution is:
-0.669626
 0.314253
The least-squares solution is:
-0.669626
 0.314253

15、稀疏矩阵

稀疏矩阵的头文件包括：

#include <Eigen/SparseCore>
#include <Eigen/SparseCholesky>
#include <Eigen/IterativeLinearSolvers>
#include <Eigen/Sparse>

初始化有两种方式：

1.使用三元组插入

typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;
triplets.reserve(estimation_of_entries); //estimation_of_entries是预估的条目
for(...)
{
    tripletList.push_back(T(i,j,v_ij));//第 i,j个有值的位置的值
}
SparseMatrixType mat(rows,cols);
mat.setFromTriplets(tripletList.begin(), tripletList.end());
// mat is ready to go!

2.直接将已知的非0值插入

SparseMatrix<double> mat(rows,cols);
mat.reserve(VectorXi::Constant(cols,6));
for(...)
{
    // i,j 个非零值 v_ij != 0
    mat.insert(i,j) = v_ij;
}
mat.makeCompressed(); // optional

稀疏矩阵支持大部分一元和二元运算:

sm1.real() sm1.imag() -sm1 0.5*sm1
sm1+sm2 sm1-sm2 sm1.cwiseProduct(sm2)

二元运算中，稀疏矩阵和普通矩阵可以混合使用

//dm表示普通矩阵
dm2 = sm1 + dm1;

16、C++数组和矩阵转换

 使用Map函数，可以实现Eigen的矩阵和c++中的数组直接转换,语法如下：

//@param MatrixType 矩阵类型
//@param MapOptions 可选参数，指的是指针是否对齐，Aligned, or Unaligned. The default is Unaligned.
//@param StrideType 可选参数，步长
/*
    Map<typename MatrixType,
        int MapOptions,
        typename StrideType>
*/
    int i;
    //数组转矩阵
    double *aMat = new double[20];
    for(i =0;i<20;i++)
    {
        aMat[i] = rand()%11;
    }
    //静态矩阵，编译时确定维数 Matrix<double,4,5> 
    Eigen:Map<Matrix<double,4,5> > staMat(aMat);
 
 
    //输出
    for (int i = 0; i < staMat.size(); i++)
        std::cout << *(staMat.data() + i) << " ";
    std::cout << std::endl << std::endl;
 
 
    //动态矩阵，运行时确定 MatrixXd
    Map<MatrixXd> dymMat(aMat,4,5);
 
 
    //输出，应该和上面一致
    for (int i = 0; i < dymMat.size(); i++)
        std::cout << *(dymMat.data() + i) << " ";
    std::cout << std::endl << std::endl;
 
    //Matrix为列优先，如下返回指针
    dymMat.data();

矩阵的应用实例：

矩阵的定义 

#include<Eigen/Dense>

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的基础使用

// Basic usage

// Eigen // Matlab // comments

x.size() // length(x) // vector size

C.rows() // size(C,1) // number of rows

C.cols() // size(C,2) // number of columns

x(i) // x(i+1) // Matlab is 1-based

C(i, j) // C(i+1,j+1) //

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.

A <<1, 2, 3, // Initialize A. The elements can also be

4, 5, 6, // matrices, which are stacked along cols

7, 8, 9; // and then the rows are stacked.

B << A, A, A; // B is three horizontally stacked A's.

A.fill(10); // Fill A with all 10's.

 Eigen特殊矩阵生成

// Eigen // Matlab

MatrixXd::Identity(rows,cols) // eye(rows,cols)

C.setIdentity(rows,cols) // C = eye(rows,cols)

MatrixXd::Zero(rows,cols) // zeros(rows,cols)

C.setZero(rows,cols) // C = ones(rows,cols)

MatrixXd::Ones(rows,cols) // ones(rows,cols)

C.setOnes(rows,cols) // C = ones(rows,cols)

MatrixXd::Random(rows,cols) // rand(rows,cols)*2-1 // MatrixXd::Random returns uniform random numbers in (-1, 1).

C.setRandom(rows,cols) // C = rand(rows,cols)*2-1

VectorXd::LinSpaced(size,low,high) // linspace(low,high,size)'

v.setLinSpaced(size,low,high) // v = linspace(low,high,size)'

 Eigen矩阵分块

// Matrix slicing and blocks. All expressions listed here are read/write.

// Templated size versions are faster. Note that Matlab is 1-based (a size N // vector is x(1)...x(N)).

// Eigen // Matlab

x.head(n) // x(1:n)

x.head() // x(1:n)

x.tail(n) // x(end - n + 1: end)

x.tail() // x(end - n + 1: end)

x.segment(i, n) // x(i+1 : i+n)

x.segment(i) // x(i+1 : i+n)

P.block(i, j, rows, cols) // P(i+1 : i+rows, j+1 : j+cols)

P.block(i, j) // P(i+1 : i+rows, j+1 : j+cols)

P.row(i) // P(i+1, :)

P.col(j) // P(:, j+1)

P.leftCols() // P(:, 1:cols)

P.leftCols(cols) // P(:, 1:cols)

P.middleCols(j) // P(:, j+1:j+cols)

P.middleCols(j, cols) // P(:, j+1:j+cols)

P.rightCols() // P(:, end-cols+1:end)

P.rightCols(cols) // P(:, end-cols+1:end)

P.topRows() // P(1:rows, :)

P.topRows(rows) // P(1:rows, :)

P.middleRows(i) // P(i+1:i+rows, :)

P.middleRows(i, rows) // P(i+1:i+rows, :)

P.bottomRows() // P(end-rows+1:end, :)

P.bottomRows(rows) // P(end-rows+1:end, :)

P.topLeftCorner(rows, cols) // P(1:rows, 1:cols)

P.topRightCorner(rows, cols) // P(1:rows, end-cols+1:end)

P.bottomLeftCorner(rows, cols) // P(end-rows+1:end, 1:cols)

P.bottomRightCorner(rows, cols) // P(end-rows+1:end, end-cols+1:end) P.topLeftCorner<rows,cols>() // P(1:rows, 1:cols)

P.topRightCorner<rows,cols>() // P(1:rows, end-cols+1:end)

P.bottomLeftCorner<rows,cols>() // P(end-rows+1:end, 1:cols)

P.bottomRightCorner<rows,cols>() // P(end-rows+1:end, end-cols+1:end)

 Eigen矩阵元素交换

// Of particular note is Eigen's swap function which is highly optimized. 

// Eigen 　　　　　　　　　　// Matlab

R.row(i) = P.col(j); 　　// R(i, :) = P(:, i)

R.col(j1).swap(mat1.col(j2)); // R(:, [j1 j2]) = R(:, [j2, j1])

Eigen矩阵转置

// Views, transpose, etc; all read-write except for .adjoint().

// Eigen // Matlab

R.adjoint() // R'

R.transpose() // R.' or conj(R')

R.diagonal() // diag(R)

x.asDiagonal() // diag(x)

R.transpose().colwise().reverse(); // rot90(R)

R.conjugate() // conj(R)

 Eigen矩阵乘积

// All the same as Matlab, but matlab doesn't have *= style operators.

// Matrix-vector. Matrix-matrix.Matrix-scalar.

y = M*x; R = P*Q; R = P*s;

a = b*M; R = P - Q; R = s*P;

a *= M; R = P + Q; R = P/s;

R *= Q; R = s*P;

R += Q; R *= s;

R -= Q; R /= s;

Eigen矩阵单个元素操作

// Vectorized operations on each element independently // Eigen // Matlab

R = P.cwiseProduct(Q); // R = P .* Q

R = P.array() * s.array();// R = P .* s

R = P.cwiseQuotient(Q); // R = P ./ Q

R = P.array() / Q.array();// R = P ./ Q

R = P.array() + s.array();// R = P + s

R = P.array() - s.array();// R = P - s

R.array() += s; // R = R + s

R.array() -= s; // R = R - s

R.array()

R.array() <= Q.array(); // R <= Q

R.cwiseInverse(); // 1 ./ P

R.array().inverse(); // 1 ./ P

R.array().sin() // sin(P)

R.array().cos() // cos(P)

R.array().pow(s) // P .^ s

R.array().square() // P .^ 2

R.array().cube() // P .^ 3

R.cwiseSqrt() // sqrt(P)

R.array().sqrt() // sqrt(P)

R.array().exp() // exp(P)

R.array().log() // log(P)

R.cwiseMax(P) // max(R, P)

R.array().max(P.array()) // max(R, P)

R.cwiseMin(P) // min(R, P)

R.array().min(P.array()) // min(R, P)

R.cwiseAbs() // abs(P)

R.array().abs() // abs(P)

R.cwiseAbs2() // abs(P.^2)

R.array().abs2() // abs(P.^2)

(R.array() < s).select(P,Q); // (R < s ? P : Q)

 

Eigen矩阵简化 

// Reductions.

int r, c;

// Eigen // Matlab

R.minCoeff() // min(R(:))

R.maxCoeff() // max(R(:))

s = R.minCoeff(&r, &c) // [s, i] = min(R(:)); [r, c] = ind2sub(size(R), i);

s = R.maxCoeff(&r, &c) // [s, i] = max(R(:)); [r, c] = ind2sub(size(R), i);

R.sum() // sum(R(:))

R.colwise().sum() // sum(R)

R.rowwise().sum() // sum(R, 2) or sum(R')'

R.prod() // prod(R(:))

R.colwise().prod() // prod(R)

R.rowwise().prod() // prod(R, 2) or prod(R')'

R.trace() // trace(R)

R.all() // all(R(:))

R.colwise().all() // all(R)

R.rowwise().all() // all(R, 2)

R.any() // any(R(:))

R.colwise().any() // any(R)

R.rowwise().any() // any(R, 2)

Eigen矩阵点乘

// Dot products, norms, etc.

// Eigen // Matlab
x.norm() // norm(x). Note that norm(R) doesn't work in Eigen. 
x.squaredNorm() // dot(x, x) Note the equivalence is not true for complex
x.dot(y) // dot(x, y)
x.cross(y) // cross(x, y) Requires #include

矩阵类型转换

//// Type conversion

// Eigen // Matlab

A.cast(); // double(A)

A.cast(); // single(A)

A.cast(); // int32(A)

A.real(); // real(A)

A.imag(); // imag(A)

// if the original type equals destination type, no work is done

  

Eigen求解线性方程组Ax = b

// Solve Ax = b. Result stored in x. Matlab: x = A  b.
x = A.ldlt().solve(b)); // A sym. p.s.d. #include 

x = A.llt() .solve(b)); // A sym. p.d. #include 

x = A.lu() .solve(b)); // Stable and fast. #include

x = A.qr() .solve(b)); // No pivoting. #include

x = A.svd() .solve(b)); // Stable, slowest. #include

// .ldlt() -> .matrixL() and .matrixD()

// .llt() -> .matrixL()

// .lu() -> .matrixL() and .matrixU()

// .qr() -> .matrixQ() and .matrixR()

// .svd() -> .matrixU(), .singularValues(), and .matrixV()

 Eigen矩阵特征值 

// Eigenvalue problems

// Eigen // Matlab

A.eigenvalues(); // eig(A);

EigenSolvereig(A); // [vecval] = eig(A)

eig.eigenvalues(); // diag(val)

eig.eigenvectors(); // vec

// For self-adjoint matrices use SelfAdjointEigenSolver<>