Eigen库的介绍

Eigen库的介绍

Eigen是一个只包含头文件的库，由于在SLAM中经常需要使用Eigen，所以想要对Eigen的结构和使用方法记录总结，方便以后查询

一、Eigen的安装

Eigen的安装十分简单，只需要在终端输入
sudo apt-get install libeigen-dev
eigen的头文件默认在“/usr/include/eigen3”，如果不确定，可以输入以下命令查找

sudo uodatedb
locate eigen3

二、Eigen的头文件

模块	头文件	作用
code	#include<Eigen/Core>	核心模块，Matrix 矩阵和 Array 数组类，基本线性代数(包括三角形和自伴随乘积)，数组操作。
geometry	#include<Eigen/Geometry>	几何模块，提供 Transform 变换、Translation 平移变换、缩放、Rotation2D 2D 旋转和3D rotations3D 旋转功能使用Quaternion四元数和AngleAxis三轴旋转的方式）
LU	#include<Eigen/LU>	Inverse逆、行列式和 LU 分解的求解方法(FullPivLU, PartialPivLU)
cholesky	#inlcude<Eigen/Cholesky>	提供 LLT 和 LDLT Cholesky 因式分解
Householder	#include<Eigen/Householder>	Householder 变换，这个模块被一些线性代数模块所使用
SVD	#inlcude<Eigen/SVD>	SVD分解与最小二乘求解方法(JacobiSVD, BDCSVD)
QR	#include<Eigen/QR>	QR 分解求解方法(HouseholderQR, ColPivHouseHolderQR, FullPivHouseholderQR)
Eigenvalues	#include<Eigen/Eigenvalues>	特征值、特征向量分解(EigenSolver, SelfAdjointEigenSolver, ComplexEigenSolver)
sparse	#include<Eigen/sparse>	稀疏矩阵存储及相关的基本线性代数 (SparseMatrix, SparseVector
dense	#include<Eigen/dense>	包含Core、Geometry、LU、SVD、QR和Eigenvalues头文件
Eigen	#include<Eigen/Eigen>	包括了Dense和Spare(整合库)

三、矩阵和向量

3.1矩阵和向量的定义

Eigen支持两类对象：用模板Eigen::Matrix表示的矩阵和向量，和用Eigen::Array表示的一维和二维数组，它们的定义如下：

typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,Options>MyMatrixType;
typedef Array<Scalar,RowsAtCompileTime,ColsAtCompileTime,Options>MyArrayType;

1. Scalar 是标量的系数类型,如：float、double、int、bool等。
2. RowsAtComileTIme 和ColsAtCompileTime 表示矩阵的行列在编译时是否已知，或者是Dynamic 动态的
3. Options 可以是 ColMajor 或 RowMajor 即列主序或行主序（存储时），默认是ColMajor
   一般情况下，即可拥有固定行列的矩阵，也可拥有动态列数的矩阵,比如

Matrix<double, 6, Dynamic>                  // 固定行数，动态列数（堆分配）
Matrix<double, Dynamic, 2>                  // 动态行数，固定列数
Matrix<double, Dynamic, Dynamic, RowMajor>  // 都是动态分配，行主序
Matrix<double, 13, 3>                       // 固定大小（一般在栈上分配）

因为类模板typedef，也可以简单的使用矩阵和数组的类型定义

// Matrix
Matrix<float,Dynamic,Dynamic>   <=>   MatrixXf
Matrix<double,Dynamic,1>        <=>   VectorXd
Matrix<int,1,Dynamic>           <=>   RowVectorXi
Matrix<float,3,3>               <=>   Matrix3f
Matrix<float,4,1>               <=>   Vector4f

// Arrays
Array<float,Dynamic,Dynamic>    <=>   ArrayXXf
Array<double,Dynamic,1>         <=>   ArrayXd
Array<int,1,Dynamic>            <=>   RowArrayXi
Array<float,3,3>                <=>   Array33f
Array<float,4,1>                <=>   Array4f

具体例子如：

//声明一个2*3de1矩阵
Eigen::Matrix<float,2,3> matrix_23;
//声明一个3*3的数组
Eigen::Array<int,3,3>array_33;
//同时，Eigen通过typedef提供了许多内置类型，不过底层依然是Eigen::Matrix
//例如Vector3d实质商就是Eigen::Matrix<double,3,1>,即三维向量
Vector3d v_3d;
//这是一样的
Matrix<float,3,1>vd_3d;

3.2矩阵的基础操作

初始化

//输入矩阵
Matrix_23<<1,2,3,4,5,6;

Vextor3d v_3d;
v1<<1,2,3;

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

Arrayxf a1(4);
a1<<1,2,3,4

//输出
cout<<"matrix 2✖3 from 1 to 6"<<matrix_23<<endl;

改变矩阵，向量的大小

//vector
vector.resize(size);   //改变大小后丢失原始数据
vector.resize(other_vector);
vector.conservativeResize(size);  //改变大小后保留矩阵中的数据

//matrix
matrix.resize(ma_rows,ma_rols);
matrix.resize(Eigen::NoChange,ma_rols);
matrix.resize(ma_rows,Eigen::NoChange);
matrix.resize(other_matrix);
matrix.conservativeresize(ma_rows,ma_rols);

用（）访问矩阵的元素

for(int i=0;i<2;i++)
	for(int j=0;j<3;j++)
	cout<<matrix_23(i,j)<<"\t";
cout<<endl;

矩阵和向量相乘

v_3d<<3,2,1;
vd_3d<<4,5,6;
//注意在Eigen里不能混合两种不同类型的矩阵，像这样就是错的
//Matrix<double,2,1>result=matrix_23*v_3d
Matrix<double,2,1>result=matrix_23.cast<double>()*v_3d
cout<<"[1,2,3;4,5,6]*[3,2,1]="result.transpose()<<endl;//输出结果的转置

Matrix<float,2,1>result2=matrix_23*v_3d;

四则运算

matrix_33=Matrix3d::random();//随机数矩阵
cout<<"random matrix:\n"<<matrix_33<<endl;
cout<<"transpose:\n"<<matrix_33.transpose()<<endl;
cout<<"sum:"<<matrix_33.sum()<<endl;   //各元素的和
cout<<"trace:"<<matrix_33.trace(0<<endl;   //迹
cout<<"inverse:\n"<<matrix_33.inverse()<<endl;  //逆
cout<<"times 10:\n"<<10*matrix_33<<endl;        //数乘
cout<<"det:\n"<<matrix_33.determinant()<<endl;    //行列式

四、geometry模块

Eigen::gemotry模块提供了各种旋转和平移的表示，和它们之间的变换方式。

4.1初始化

Eigen中四种旋转的双精度表示

Eigen::Matrix3d//旋转矩阵（3*3）
Eigen::AngleAxisd//旋转向量（3*1）
Eigen::Vector3d//欧拉角（3*1）
Eigen::quaterniond//四元数（4*1）
Eigen::Isometry  //欧式变换矩阵（4*4）

//3D旋转矩阵可以直接使用Matrix3d或matrix3f
matrix3d rotation_matrix=matrix3d::Identitf();//初始化为单位矩阵
//旋转向量使用AngleAxis
AngleAxisd ratation_vector(M_PI/4,vector3d(0,0,1));//沿z轴旋转45度
//将旋转向量转化为旋转矩阵
cout.precision(3)//控制输出的精度
cout<<"rotation matrix=\n"<<rotation_vector.matrix()<<endl;
//同理，矩阵可以直接赋值
rotation_matrix=rotation_vector.toRotationMatrix();

//使用AngleAxis实现矩阵旋转
Vector3d v(1,0,0);
Vector3d v_rotated=rotation_vector*v;

//使用旋转矩阵实现旋转
v_rotated=rotation_matrix*v;

4.2矩阵变换

欧拉角：将旋转矩阵变为欧拉角

vector3d euler_angles=rotation_matrix.eulerAngles(2,1,0);   //ZYX顺序
cout<<"yaw pitch roll="<<euler_angle.transpose()<<endl;

欧式变换矩阵的初始化

//欧式变换矩阵可以使用Eigen::Isometry
Isometry3d T=Isometry::Identity();//虽然是3D，但是实际上是4*4的矩阵
T.rotate(rotation_vector);           //旋转矩阵
T.pretranslate(Vector3d(1,3,4));     //平移向量
cout<<"Transform matrix=\n"<<T.matrix()<<endl;

用变换矩阵进行坐标变换

Vector3d v_transformed=T*v;     //相当于R*V+T
cout<<"v transform="<<v_transform.transpose()<<endl;

四元数的初始化

//法1：直接把AngleAxis复赋值给四元数。反之亦然
Quaterniond q=Quaterniond(rotation_vector);
//请注意coffes的顺序为（x,y,z,w),w为实部，前三者为虚部
cout<<"quaternion form rotation vector="<<q.coeffs().transpose()<<endl;
//法2：也可以把旋转矩阵赋给它
q=Quaterniond(rotation_matrix);

用四元数表示旋转

v_rotated=q*v  //注意数学上是qvq^（-1）

五、Ax=b的求解方法

	//解方程
    //求解matrix_NN*X=V_nd方程
    Matrix<double,Matrix_size,Matrix_size>matrix_NN=MatrixXd::Random(Matrix_size,Matrix_size);
    matrix_NN=matrix_NN*matrix_NN.transpose();
    Matrix<double,Matrix_size,1>v_Nd=MatrixXd::Random(Matrix_size,1);

    clock_t time_stt=clock();
    //直接求逆
    Matrix<double,Matrix_size,1> x=matrix_NN.inverse()*v_Nd;
    cout<<"time of normal inverse is"<<1000*(clock()-time_stt)/(double)CLOCKS_PER_SEC<<"ms"<<endl;
    cout<<"x="<<x.transpose()<<endl;

    //使用QR分解
    time_stt=clock();
    x=matrix_NN.colPivHouseholderQr().solve(v_Nd);
    cout<<"time of normal inverse is"<<1000*(clock()-time_stt)/(double)CLOCKS_PER_SEC<<"ms"<<endl;
    cout<<"x="<<x.transpose()<<endl;

    //使用cholesky分解
    time_stt=clock();
    x=matrix_NN.ldlt().solve(v_Nd);
    cout<<"time of normal inverse is"<<1000*(clock()-time_stt)/(double)CLOCKS_PER_SEC<<"ms"<<endl;
    cout<<"x="<<x.transpose()<<endl;
    return 0;