1.模块和头文件


模块 头文件 说明
Core #include<Eigen/Core> 包含Matrix和Array类,基础的线性代数运算和数组操作
Gemetry #include<Eigen/Geometry> 包含旋转,平移,缩放,2维和3维的各种变换
LU #include<Eigen/LU> 包含求逆,行列式,LU分解
Cholesky #include<Eigen/Cholesky> 包含LLT和LDLT Cholesky分解
SVD #include<Eigen/SVD> 包含SVD分解
QR #include<Eigen/QR> 包含QR分解
Eigenvalues #include<Eigen/Eigenvalues> 包含特征值,特征向量的分解
Sparse #include<Eigen/Sparse> 包含稀疏矩阵的存储与运算
Dense #include<Eigen/Dense> 包含了Core/Geometry/LU/Cholesky/SVD/QR/Eigenvalues模块
Eigen #include<Eigen/Eigen> 包含Dense和Sparse

2.语法基础


2.1.Matrix


所有矩阵和向量都是Matrix模板类的对象,Matrix类有6个模板参数,主要使用前三个,剩下的使用默认值。


  • 默认构造时,指定大小的矩阵,只分配相应大小的空间,不进行初始化。动态大小的矩阵,则未分配空间。
  • []操作符可以用于向量元素的获取,但不能用于matrix。
  • matrix的大小可以通过rows(), cols(), size()获取,resize()可以重新调整矩阵大小。

2.2 Storage orders


  • 矩阵的条目形成一个二维网格。但是,当矩阵存储在存储器中时,必须以某种方式线性排列条目。有两种主要方法可以做到这一点,按行和按列。
  • 我们说矩阵是按行优先存储的。首先存储整个第一行,然后存储整个第二行,依此类推。

Matrix<int, 3, 4, ColMajor> Acolmajor;
Acolmajor << 8, 2, 2, 9,
9, 1, 4, 4,
3, 5, 4, 5;
cout << “The matrix A:” << endl;
cout << Acolmajor << endl << endl;
cout << “In memory (column-major):” << endl;
for (int i = 0; i < Acolmajor.size(); i++)
cout << _(Acolmajor.data() + i) << “ “;
cout << endl << endl;
Matrix<int, 3, 4, RowMajor> Arowmajor = Acolmajor;
cout << “In memory (row-major):” << endl;
for (int i = 0; i < Arowmajor.size(); i++)
cout << _(Arowmajor.data() + i) << “ “;
cout << endl;

    output:


    The matrix A:
    8 2 2 9
    9 1 4 4
    3 5 4 5

    In memory (column-major):
    8 9 3 2 1 5 2 4 4 9 4 5

    In memory (row-major):
    8 2 2 9 9 1 4 4 3 5 4 5

      那么,您应该在程序中使用哪个存储顺序?这个问题没有简单的答案。这取决于您的应用程序。请记住以下几点:


      • 您的用户可能希望您使用特定的存储顺序。或者,您可以使用Eigen以外的其他库,并且这些其他库可能需要一定的存储顺序。在这些情况下,在整个程序中使用此存储顺序可能是最简单,最快的。
      • 当矩阵以行优先顺序存储时,逐行遍历矩阵的算法会更快,因为数据位置更好。同样,对于主要列矩阵,逐列遍历更快。可能需要尝试一下以找出对您的特定应用程序更快的方法。
      • Eigen中的默认值是column-major。自然,对Eigen库的大多数开发和测试都是使用列主矩阵完成的。这意味着,即使我们旨在透明地支持列主存储和行主存储顺序,Eigen库也可能很好地与列主存储矩阵配合使用。

      2.2.矩阵与向量的运算


      MatrixXcf a = MatrixXcf::Random(3,3);
      a.transpose(); # 转置
      a.conjugate(); # 共轭
      a.adjoint(); # 共轭转置(伴随矩阵)
      # 对于实数矩阵,conjugate不执行任何操作,adjoint等价于transpose
      a.transposeInPlace() #原地转置

      Vector3d v(1,2,3);
      Vector3d w(4,5,6);
      v.dot(w); # 点积
      v.cross(w); # 叉积

      Matrix2d a;
      a << 1, 2, 3, 4;
      a.sum(); # 所有元素求和
      a.prod(); # 所有元素乘积
      a.mean(); # 所有元素求平均
      a.minCoeff(); # 所有元素中最小元素
      a.maxCoeff(); # 所有元素中最大元素
      a.trace(); # 迹,对角元素的和
      #minCoeff和maxCoeff还可以返回结果元素的位置信息
      int i, j;
      a.minCoeff(&i, &j);

      行列相取:row 取整行


      Eigen::Matrix<double, 3, 4> &Pose0
      Pose0.row(2) //取第三行
      Pose0.row(1) //取第2行
      Pose0.row(0) //取第1行
      • 1
      • 2
      • 3
      • 4

      2.3.Array类


      Array是个类模板,前三个参数必须指定,后三个参数可选。


      • 当执行array_array时,执行的是相应元素的乘积,所以两个array必须具有相同的尺寸。
      • Matrix对象——>Array对象:.array()函数
      • Array对象——>Matrix对象:.matrix()函数

      2.4.矩阵初始化


      • 逗号初始化:为矩阵元素赋值,顺序是从左到右,从上到下,数目必须匹配。
      • 特殊矩阵

        • 零阵:类静态成员函数Zero()

        • 常量矩阵:Constant(rows, cols, value)

        • 随机矩阵:Random()

        • 单位矩阵:Identity()
      • LinSpaced(size, low, high):构建从low到high等间距的size长度的序列,适用于vector和一维数组。

      2.5. 归约,迭代器,广播


      范数计算


      • squareNorm():L2范数,等价于计算vector自身点积
      • norm():返回`squareNorm的开方根
      • .lpNorm

        ():p范数,p可以取Infinity,表无穷范数


      布尔归约


      • all()=true: matrix或array中所有元素为true
      • any()=true: 到少有一个为true
      • count(): 返回true元素个数

      迭代器,获取某元素位置


      // sample
      Eigen::MatrixXf m(2,2);
      m << 1,2,3,4;
      MatrixXf::Index maxRow, maxCol;
      float max = m.maxCoeff(&minRow, &minCol);

        2.6 Eigen 矩阵单个元素操作


        2.6.1 vector operator


        // 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全加s
        R.array() -= s; // R = R - s全减s
        R.array() < Q.array(); // R < Q 以下的都是针对矩阵的单个元素的操作
        R.array() <= Q.array(); // R <= Q矩阵元素比较,会在相应位置置0或1
        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)这个也是单个元素的操作

        2.6.2 matrix operator


        #include <iostream>
        using namespace std;
        #include <ctime>
        // Eigen 部分
        #include <Eigen/Core>
        // 稠密矩阵的代数运算(逆,特征值等)
        #include <Eigen/Dense>

        #define MATRIX_SIZE 50
        ____________________________
        _ 本程序演示了 Eigen 基本类型的使用
        ____________________________/

        int main( int argc, char__ argv )
        {
        // Eigen 中所有向量和矩阵都是Eigen::Matrix,它是一个模板类。它的前三个参数为:数据类型,行,列
        // 声明一个2_3的float矩阵
        Eigen::Matrix<float, 2, 3> matrix_23;

        // 同时,Eigen 通过 typedef 提供了许多内置类型,不过底层仍是Eigen::Matrix
        // 例如 Vector3d 实质上是 Eigen::Matrix<double, 3, 1>,即三维向量
        Eigen::Vector3d v_3d;
        // 这是一样的
        Eigen::Matrix<float,3,1> vd_3d;

        // Matrix3d 实质上是 Eigen::Matrix<double, 3, 3>
        Eigen::Matrix3d matrix_33 = Eigen::Matrix3d::Zero(); //初始化为零
        // 如果不确定矩阵大小,可以使用动态大小的矩阵
        Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_dynamic;
        // 更简单的
        Eigen::MatrixXd matrix_x;
        // 这种类型还有很多,我们不一一列举

        // 下面是对Eigen阵的操作
        // 输入数据(初始化)
        matrix_23 << 1, 2, 3, 4, 5, 6;
        // 输出
        cout << matrix_23 << 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里你不能混合两种不同类型的矩阵,像这样是错的
        // Eigen::Matrix<double, 2, 1> result_wrong_type = matrix_23 _ v_3d;
        // 应该显式转换
        Eigen::Matrix<double, 2, 1> result = matrix_23.cast<double>() _ v_3d;
        cout << result << endl;

        Eigen::Matrix<float, 2, 1> result2 = matrix_23 _ vd_3d;
        cout << result2 << endl;

        // 同样你不能搞错矩阵的维度
        // 试着取消下面的注释,看看Eigen会报什么错
        // Eigen::Matrix<double, 2, 3> result_wrong_dimension = matrix_23.cast<double>() _ v_3d;

        // 一些矩阵运算
        // 四则运算就不演示了,直接用+-_/即可。
        matrix_33 = Eigen::Matrix3d::Random(); // 随机数矩阵
        cout << matrix_33 << endl << endl;

        cout << matrix_33.transpose() << endl; // 转置
        cout << matrix_33.sum() << endl; // 各元素和
        cout << matrix_33.trace() << endl; // 迹
        cout << 10_matrix_33 << endl; // 数乘
        cout << matrix_33.inverse() << endl; // 逆
        cout << matrix_33.determinant() << endl; // 行列式

        // 特征值
        // 实对称矩阵可以保证对角化成功
        Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigen_solver ( matrix_33.transpose()_matrix_33 );
        cout << “Eigen values = \n” << eigen_solver.eigenvalues() << endl;
        cout << “Eigen vectors = \n” << eigen_solver.eigenvectors() << endl;

        // 解方程
        // 我们求解 matrix_NN _ x = v_Nd 这个方程
        // N的大小在前边的宏里定义,它由随机数生成
        // 直接求逆自然是最直接的,但是求逆运算量大

        Eigen::Matrix< double, MATRIX_SIZE, MATRIX_SIZE > matrix_NN;
        matrix_NN = Eigen::MatrixXd::Random( MATRIX_SIZE, MATRIX_SIZE );
        Eigen::Matrix< double, MATRIX_SIZE, 1> v_Nd;
        v_Nd = Eigen::MatrixXd::Random( MATRIX_SIZE,1 );

        clock_t time_stt = clock(); // 计时
        // 直接求逆
        Eigen::Matrix<double,MATRIX_SIZE,1> x = matrix_NN.inverse()_v_Nd;
        cout <<“time use in normal inverse is “ << 1000_ (clock() - time_stt)/(double)CLOCKS_PER_SEC << “ms”<< endl;

        // 通常用矩阵分解来求,例如QR分解,速度会快很多
        time_stt = clock();
        x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
        cout <<“time use in Qr decomposition is “ <<1000_ (clock() - time_stt)/(double)CLOCKS_PER_SEC <<“ms” << endl;

        return 0;
        }




        2.7 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)列求和1×N
        R.rowwise().sum() // sum(R, 2) or sum(R’)’行求和N×1
        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) 或运算

        2.8 Eigen 求解线性方程组 Ax = b


        // Solve Ax = b. Result stored in x. Matlab: x = A \ b.
        x = A.ldlt().solve(b)); // #include <Eigen/Cholesky>LDLT分解法实际上是Cholesky分解法的改进
        x = A.llt() .solve(b)); // A sym. p.d. #include <Eigen/Cholesky>
        x = A.lu() .solve(b)); // Stable and fast. #include <Eigen/LU>
        x = A.qr() .solve(b)); // No pivoting. #include <Eigen/QR>
        x = 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()

          2.9 矩阵的特殊生成


          // 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 // 元素随机在-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)’线性分布的数组

            3.Eigen线性方程与矩阵分解


            • Eigen提供了解线性方程的计算方法,包括LU分解法,QR分解法,SVD(奇异值分解)、特征值分解等。其中第一部分介绍了各头文件所包含的内容。
            • 对于一般形式如下的线性系统:





              A


              x


              =


              b



              {Ax=b}


              Ax=b
            • 解决上述方程的方式一般是将矩阵A进行分解,当然最基本的方法是高斯消元法。

            Eigen 官方给的一个案例:


            #include <iostream>
            #include <Eigen/Dense>

            using namespace std;
            using namespace Eigen;

            int main()
            {
            Matrix3f A;
            Vector3f b;
            A << 1,2,3, 4,5,6, 7,8,10;
            b << 3,3,4;
            cout<<“Here is the Matrix A:\n”<< A <<endl;
            cout<<“ Here is the vector b:\n”<< b <<endl;
            Vector3f x = A.colPivHouseholderQr().solve(b);
            cout<<“The solution is:\n”<<x<<endl;
            return 0;
            }

            使用这些接口也可以解决矩阵相乘的问题:


            #include <iostream>
            #include <Eigen/Dense>

            using namespace std;
            using namespace Eigen;

            int main()
            {
            Matrix2f A,b;
            A << 2,-1,-1,3;
            b << 1,2,3,1;
            cout<<“Here is the matrix A:\n”<<A<<endl;
            cout<<“Here is the right hand side b:\n”<<b<<endl;
            Matrix2f x = A.ldlt().solve(b);
            cout<<“The solution is:\n”<<x<<endl;
            return 0;
            }

            Eigen也提供了计算特征值和特征向量的算法:


            #include <iostream>
            #include <Eigen/Dense>

            using namespace std;
            using namespace Eigen;

            int main()
            {
            Matrix2f A;
            A << 1,2,2,3;
            cout<<“Here is the matrix A:\n”<<A<<endl;
            SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
            if( eigensolver.info() != Success ) abort();
            cout<<“ The eigenvalues of A are:\n”<<eigensolver.eigenvalues()<<endl;
            cout<<“ Here is a matrix whose columns are eigenvectors of A\n”
            <<“ corresponding to these eigenvalues:\n”
            <<eigensolver.eigenvectors()<<endl;
            return 0;
            }

            Eigen 也提供了求逆矩阵和求矩阵行列式的算法,但是这两种算法对于大型矩阵来说都是非常不经济的算法,当需要对大型矩阵做这种的操作时,需要自己判断到底需不需这样做。但是对于小型矩阵 则可以没有顾虑地使用。


            #include <iostream>
            #include <Eigen/Dense>

            using namespace std;
            using namespace Eigen;

            int main()
            {
            Matrix3f A;
            A << 1,2,1,
            2,1,0,
            -1,1,2;

            cout<<“Here is the matrix A:\n”<<A<<endl;
            cout<<“The determinant of A is “<<A.determinant()<<endl;
            cout<<“The inverse of A is:\n”<<A.inverse()<<endl;
            return 0;
            }

            Eigen也提供了解最小二乘问题的解法,并给出两种实现,分别是BDCSVD和JacobiSVD,其中推荐使用的一种是BDCSVD。


            #include <iostream>
            #include <Eigen/Dense>

            using namespace std;
            using namespace Eigen;

            int main()
            {
            MatrixXf A = MatrixXf::Random(3,2);
            cout<<“Here is the matrix A:\n”<<A<<endl;
            VectorXf b = VectorXf::Random(3);
            cout<<“Here is the right hand side b:\n”<<b<<endl;
            cout<<“The least-squares solution is:\n”
            <<A.bdcSvd(ComputeThinU|ComputeThinV).solve(b)<<endl;
            return 0;
            }