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（奇异值分解）、特征值分解等。其中第一部分介绍了各头文件所包含的内容。
- 对于一般形式如下的线性系统：
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;
}
扫一扫
1825
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
