Eigen 使用总结
- 1.模块和头文件
- 2.语法基础
- 2.1.Matrix类
- 2.2 Storage orders
- 2.2.矩阵与向量的运算
- 2.3.Array类
- 2.4.矩阵初始化
- 2.5. 归约，迭代器，广播
- 范数计算
- 布尔归约
- 迭代器，获取某元素位置
- 2.6 Eigen 矩阵单个元素操作
- 2.6.1 vector operator
- 2.6.2 matrix operator
- 2.7 Eigen 矩阵化简
- 2.8 Eigen 求解线性方程组 Ax = b
- 2.9 矩阵的特殊生成
- 3.Eigen线性方程与矩阵分解
1.模块和头文件
| 模块 | 头文件 | 说明 |
|---|---|---|
| Core | #include |
包含Matrix和Array类，基础的线性代数运算和数组操作 |
| Gemetry | #include |
包含旋转，平移，缩放，2维和3维的各种变换 |
| LU | #include |
包含求逆，行列式，LU分解 |
| Cholesky | #include |
包含LLT和LDLT Cholesky分解 |
| SVD | #include |
包含SVD分解 |
| QR | #include |
包含QR分解 |
| Eigenvalues | #include |
包含特征值，特征向量的分解 |
| Sparse | #include |
包含稀疏矩阵的存储与运算 |
| Dense | #include |
包含了Core/Geometry/LU/Cholesky/SVD/QR/Eigenvalues模块 |
| Eigen | #include |
包含Dense和Sparse |
- 默认构造时，指定大小的矩阵，只分配相应大小的空间，不进行初始化。动态大小的矩阵，则未分配空间。
- []操作符可以用于向量元素的获取，但不能用于matrix。
- matrix的大小可以通过rows(), cols(), size()获取，resize()可以重新调整矩阵大小。
- 矩阵的条目形成一个二维网格。但是，当矩阵存储在存储器中时，必须以某种方式线性排列条目。有两种主要方法可以做到这一点，按行和按列。
- 我们说矩阵是按行优先存储的。首先存储整个第一行，然后存储整个第二行，依此类推。
Matrix 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 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 5In memory (column-major):
893215244945In memory (row-major):
822991443545
那么，您应该在程序中使用哪个存储顺序？这个问题没有简单的答案。这取决于您的应用程序。请记住以下几点：
- 您的用户可能希望您使用特定的存储顺序。或者，您可以使用Eigen以外的其他库，并且这些其他库可能需要一定的存储顺序。在这些情况下，在整个程序中使用此存储顺序可能是最简单，最快的。
- 当矩阵以行优先顺序存储时，逐行遍历矩阵的算法会更快，因为数据位置更好。同样，对于主要列矩阵，逐列遍历更快。可能需要尝试一下以找出对您的特定应用程序更快的方法。
- Eigen中的默认值是column-major。自然，对Eigen库的大多数开发和测试都是使用列主矩阵完成的。这意味着，即使我们旨在透明地支持列主存储和行主存储顺序，Eigen库也可能很好地与列主存储矩阵配合使用。
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 &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()函数
- 逗号初始化：为矩阵元素赋值，顺序是从左到右，从上到下，数目必须匹配。
- 特殊矩阵
-
- 零阵：类静态成员函数Zero()
-
- 常量矩阵：Constant(rows, cols, value)
-
- 随机矩阵：Random()
-
- 单位矩阵：Identity()
- LinSpaced(size, low, high)：构建从low到high等间距的size长度的序列，适用于vector和一维数组。
- 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
using namespace std;
#include
// Eigen 部分
#include
// 稠密矩阵的代数运算（逆，特征值等）
#include #define MATRIX_SIZE 50
****************************
* 本程序演示了 Eigen 基本类型的使用
****************************/int main( int argc, char** argv )
{
// Eigen 中所有向量和矩阵都是Eigen::Matrix，它是一个模板类。它的前三个参数为：数据类型，行，列
// 声明一个2*3的float矩阵
Eigen::Matrix matrix_23;
// 同时，Eigen 通过 typedef 提供了许多内置类型，不过底层仍是Eigen::Matrix
// 例如 Vector3d 实质上是 Eigen::Matrix，即三维向量
Eigen::Vector3d v_3d;
// 这是一样的
Eigen::Matrix vd_3d;
// Matrix3d 实质上是 Eigen::Matrix
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< result_wrong_type = matrix_23 * v_3d;
// 应该显式转换
Eigen::Matrix result = matrix_23.cast() * v_3d;
cout << result << endl;
Eigen::Matrix result2 = matrix_23 * vd_3d;
cout << result2 << endl;
// 同样你不能搞错矩阵的维度
// 试着取消下面的注释，看看Eigen会报什么错
// Eigen::Matrix result_wrong_dimension = matrix_23.cast() * 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_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 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 LDLT分解法实际上是Cholesky分解法的改进
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()
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进行分解，当然最基本的方法是高斯消元法。
#include
#include 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 < 【Eigen使用总结】使用这些接口也可以解决矩阵相乘的问题：
#include
#include 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"< Eigen也提供了计算特征值和特征向量的算法：
#include
#include using namespace std;
using namespace Eigen;
int main()
{
Matrix2f A;
A << 1,2,2,3;
cout<<"Here is the matrix A:\n"< eigensolver(A);
if( eigensolver.info() != Success ) abort();
cout<<" The eigenvalues of A are:\n"< Eigen 也提供了求逆矩阵和求矩阵行列式的算法，但是这两种算法对于大型矩阵来说都是非常不经济的算法，当需要对大型矩阵做这种的操作时，需要自己判断到底需不需这样做。但是对于小型矩阵 则可以没有顾虑地使用。
#include
#include 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"< Eigen也提供了解最小二乘问题的解法，并给出两种实现，分别是BDCSVD和JacobiSVD,其中推荐使用的一种是BDCSVD。
#include
#include using namespace std;
using namespace Eigen;
int main()
{
MatrixXf A = MatrixXf::Random(3,2);
cout<<"Here is the matrix A:\n"< 
