Eigen使用实践

Eigen使用实践

一、安装 Eigen

1、sudo apt-get install libeigen3-dev
2、sudo updatedb
3、locate eigen3

二、Eigen简单使用

#include <iostream>
using namespace std;

#include <ctime>

//Eigen核心部分

#include <eigen3/Eigen/Dense>
//稠密矩阵的代数运算（逆，特征值等）

#include <eigen3/Eigen/Core>

using namespace Eigen;

#define MATRIX_SIZE 50

int main(int argc,char **argv){
   
    
    //Eigen中所有向量和矩阵都是Eigen::Matrix,它是一个模板类。它的前三个参数为数据类型、行、列
    
    //声明一个2＊3的float矩阵
    
    Matrix<float,2,3> matrix_23;
    
    //同时Eigen通过typedef提供了许多内置类型，不过底层仍是Eigen::Matrix
    
    //例如Vector3d 实质上是Eigen::Matrix<double,3,1>,即三维向量
    
    Vector3d v_3d;
    
    //这是一样的
    
    Matrix<float,3,1> vd_3d;
    
    //Matrix3d实质上是Eigen::Matrix<double,3,3>
    
    Matrix3d matrix_33=Matrix3d::Zero();//初始化为零
    
    //如果不确定矩阵大小，可以使用动态大小的矩阵
    
    Matrix<double,Dynamic,Dynamic> matrix_dynamic;
    
    //更简单的
    
    MatrixXd matrix_x;
    
    //这种类型还有很多，不一一列举
    //下面是对Eigen阵的操作
    //输入数据（初始化）
    
    matrix_23<<1,2,3,4,5,6;
    
    //输出
    
    cout<<"matrix 2X3 from 1 to 6: \n"<<matrix_23<<endl;
    
    //用（）访问矩阵中元素
    
    cout<<"print matrix 2X3:"<<endl;
    
    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_wrong_type=martix_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*vd_3d;
    cout<<"[1,2,3,4,5,6]*[4,5,6]:"<<result2.transpose()<<endl;
    
    //同样，你不能搞错矩阵的维度
    //试着取消下面的注释，看看Eigen会报什么错
    //Eigen::Matrix<double,2,3>reslt_wrong_dimension=matrix_23.cast<double>()*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()<<endl;     //迹
    
    cout<<"times 10: \n"<<10*matrix_33<<endl;  //数乘
    
    cout<<"inverse: \n"<<matrix_33.inverse()<<endl;  //逆
    
    cout<<"det: "<<matrix_33.determinant()<<endl;  //行列式
    
    
    //特征值
    //实对称矩阵可以保证对角化成功
    
    SelfAdjointEigenSolver<Matrix3d> eigen_solver(matrix_33.transpose()*matrix_33);
    
    cout<<"Eigen values= \n"<<eigen_solver.eigenvalues()<<endl;
    
    cout<<"Eigen vactors= \n"<<eigen_solver.eigenvectors()<<endl;
    
    //解方程
    //我们求解matrix_NN*x=v_Nd方程
    

    //N的大小在前面的宏里定义，它由随机数生成
    //直接求逆自然是最直接的，但是运算量大
    
    Matrix< double, MATRIX_SIZE, MATRIX_SIZE > matrix_NN;
    
    matrix_NN=MatrixXd::Random(MATRIX_SIZE,MATRIX_SIZE);
    
    matrix_NN=matrix_NN*matrix_NN.transpose(); //保证半正定
    
    Matrix<double,MATRIX_SIZE,1> v_Nd;
    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 Qr decomposition 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 idlt decomposition is"<<1000*(clock()-time_stt)/(double)CLOCKS_PER_SEC<<"ms"<<endl;
    cout<<"x= "<<x.transpose()<<endl;
    

    
    return 0;
    

    
}

三、运行结果

random matrix: 
 0.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
transpose: 
 0.680375 -0.211234  0.566198
  0.59688  0.823295 -0.604897
-0.329554  0.536459 -0.444451
sum:1.61307
trace:1.05922
times 10: 
 6.80375   5.9688 -3.29554
-2.11234  8.23295  5.36459
 5.66198 -6.04897 -4.44451
inverse: 
-0.198521   2.22739    2.8357
  1.00605 -0.555135  -1.41603
 -1.62213   3.59308   3.28973
det: 0.208598
Eigen values= 
0.0242899
 0.992154
  1.80558
Eigen vactors= 
-0.549013 -0.735943  0.396198
 0.253452 -0.598296 -0.760134
-0.796459  0.316906 -0.514998
time of normal inverse is：115.646ms
x=-55.7896 -298.793  130.113 -388.455 -159.312  160.654 -40.0416 -193.561  155.844  181.144  185.125 -62.7786  19.8333 -30.8772 -200.746  55.8385 -206.604  26.3559 -14.6789  122.719 -221.449   26.233  -318.95 -78.6931  50.1446  87.1986 -194.922  132.319  -171.78 -4.19736   11.876 -171.779  48.3047  84.1812 -104.958 -47.2103 -57.4502 -48.9477 -19.4237  28.9419  111.421  92.1237 -288.248 -23.3478  -275.22 -292.062  -92.698  5.96847 -93.6244  109.734
time of Qr decomposition is 5.435ms
x= -55.7896 -298.793  130.113 -388.455 -159.312  160.654 -40.0416 -193.561  155.844  181.144  185.125 -62.7786  19.8333 -30.8772 -200.746  55.8385 -206.604  26.3559 -14.6789  122.719 -221.449   26.233  -318.95 -78.6931  50.1446  87.1986 -194.922  132.319  -171.78 -4.19736   11.876 -171.779  48.3047  84.1812 -104.958 -47.2103 -57.4502 -48.9477 -19.4237  28.9419  111.421  92.1237 -288.248 -23.3478  -275.22 -292.062  -92.698  5.96847 -93.6244  109.734
time of idlt decomposition is1.209ms
x= -55.7896 -298.793  130.113 -388.455 -159.312  160.654 -40.0416 -193.561  155.844  181.144  185.125 -62.7786  19.8333 -30.8772 -200.746  55.8385 -206.604  26.3559 -14.6789  122.719 -221.449   26.233  -318.95 -78.6931  50.1446  87.1986 -194.922  132.319  -171.78 -4.19736   11.876 -171.779  48.3047  84.1812 -104.958 -47.2103 -57.4502 -48.9477 -19.4237  28.9419  111.421  92.1237 -288.248 -23.3478  -275.22 -292.062  -92.698  5.96847 -93.6244  109.734
*** Finished ***

四、初次运行出现报错信息

/home/zbz/projects/eigenMartrix/build> make -j4
[ 50%] Building CXX object CMakeFiles/eigenmartrix.dir/main.cpp.o
In file included from /usr/local/include/eigen3/Eigen/Core:348:0,
                 from /usr/local/include/eigen3/Eigen/Dense:1,
                 from /home/zbz/projects/eigenMartrix/main.cpp:9:
/usr/local/include/eigen3/Eigen/src/Core/Product.h: In instantiation of ‘struct Eigen::internal::product_result_scalar<Eigen::Matrix<float, 2, 3>, Eigen::Matrix<double, 3, 1>, Eigen::DenseShape, Eigen::DenseShape>’:
/usr/local/include/eigen3/Eigen/src/Core/Product.h:78:73:   required from ‘struct Eigen::internal::traits<Eigen::Product<Eigen::Matrix<