Eigen使用总结

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行

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（奇异值分解）、特征值分解等。其中第一部分介绍了各头文件所包含的内容。
• 对于一般形式如下的线性系统：
  $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;
}