一、Eigen库介绍
简介
Eigen [1]目前最新的版本是3.4,除了C++标准库以外,不需要任何其他的依赖包。Eigen使用的CMake建立配置文件和单元测试,并自动安装。如果使用Eigen库,只需包特定模块的的头文件即可。
基本功能
Eigen适用范围广,支持包括固定大小、任意大小的所有矩阵操作,甚至是稀疏矩阵;支持所有标准的数值类型,并且可以扩展为自定义的数值类型;支持多种矩阵分解及其几何特征的求解;它不支持的模块生态系统 [2]提供了许多专门的功能,如非线性优化,矩阵功能,多项式解算器,快速傅立叶变换等。
Eigen支持多种编译环境,开发人员对库中的实例在多种编译环境下经过测试,以保证其在不同编译环境下的可靠性和实用性。
介绍来之百度百科,点这里
二、Eigen库的下载及安装
进入官方下载路径,点这里
下载eigen-3.4.0.tar.gz到本地,解压
tar -xzf eigen-3.4.0.tar.gz
编译项目
# 进入eigen源码目录
cd eigen-3.4.0
# 创建文件夹build-编译用
mkdir build
# 进入build
cd build
# 开始构建
cmake ..
安装项目
# 执行安装 - 注意这里不需要make,原因前面解释过eigen是模式式编程
# 不需要真正的编译,只需要把头文件安装到系统include目录即可
make install
上面命令执行后,默认安装 Eigen采用源码的方式提供给用户使用,在使用时只需要包含Eigen的头文件即可进行使用。之所以采用这种方式,是因为Eigen采用模板方式实现,由于模板函数不支持分离编译,所以只能提供源码而不是动态库的方式供用户使用。到/usr/local/include/
目录
三、Eigen基本使用
引入Eigen库
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
1、定义矩阵、基本访问、为矩阵赋值、重置矩阵
// 1. 定义矩阵
Eigen::MatrixXd m(2,2);
Eigen::Vector3d vec3d;
Eigen::Vector4d vec4d(1.0, 2.0, 3.0, 4.0);
// Matrix创建的矩阵默认是按列存储,Eigen在处理按列存储的矩阵时会更加高效,可通过如下方法修改优先存储方式
Matrix<int,3, 4, ColMajor> Acolmajor;
Matrix<int,3, 4, RowMajor> Arowmajor;
// 2. 定义动态矩阵、静态矩阵
Eigen::MatrixXd matrixXd; // 动态矩阵,通过mat.resise()重新修改矩阵大小
Eigen::Matrix3d matrix3d;
// 3. 矩阵元素的访问
m(0, 0) = 1;
m(0, 1) = 2;
m(1, 0) = m(0, 0) + 3;
m(1, 1) = m(0, 0) * m(0, 1);
cout << m << endl << endl;
// 4. 设置矩阵的元素
m << -5.5, 3.1,
5.9, 1.0;
cout << m << endl << endl;
int row = 4;
int col = 5;
Eigen::MatrixXf matrixXf(row, col);
matrixXf << 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20;
cout << matrixXf << endl << endl;
// 5. 重置矩阵
Eigen::MatrixXd matrixXd1(4,4);
cout << "原m的行列为:" << m.rows() << " , " << m.cols() << endl << endl;
m.resize(5,100); // 通过resize修改矩阵大小
cout << "新1的m行列为:" << m.rows() << " , " << m.cols() << endl << endl;
m = matrixXd1; // 通过直接赋值新矩阵,修改矩阵大小
cout << "新2的m行列为:" << m.rows() << " , " << m.cols() << endl << endl;
2、常用矩阵、及常用的基本矩阵运算
常用特殊矩阵:0矩阵、1矩阵、单位矩阵、随机矩阵
常用的矩阵运算:转置、共轭、共轭转置、求逆、求迹、求行列式、求模、求和、求均值、求乘积、求最小、求最大
点乘、叉乘
// 1. 其他常用特殊的矩阵
Eigen::Matrix3d m3;
m3 = Eigen::Matrix3d::Zero();
cout << "Zero() : " << endl << m3 << endl << endl;
m3 = Eigen::Matrix3d::Identity();
cout << "Identity() : " << endl << m3 << endl << endl;
m3 = Eigen::Matrix3d::Random();
cout << "Random() : " << endl << m3 << endl << endl;
// 2. 直接设置为常用特殊矩阵
m3.setRandom(); // 随机矩阵
m3.setOnes(); // 全1矩阵
m3.setZero(); // 全0矩阵
m3.setIdentity(); // 单位矩阵
// 3. 矩阵的运算1
Eigen::Matrix3d new_m3;
new_m3 = m3.transpose();
cout << "============================================" << endl << endl;
cout << "m3.transpose() - 转置: " << endl << new_m3 << endl << endl;
new_m3 = m3.conjugate();
cout << "m3.conjugate() - 共轭: " << endl << new_m3 << endl << endl;
new_m3 = m3.adjoint();
cout << "m3.adjoint() - 伴随矩阵(共轭转置): " << endl << new_m3 << endl << endl;
// 对于实数矩阵,conjugate不执行任何操作,adjoint等价于transpose
cout << "m3.inverse() - 求逆: " << endl << m3.inverse(); << endl << endl;
cout << "m3.trace() - 求迹: " << endl << m3.trace() << endl << endl;
cout << "m3.determinant() - 求行列式: " << endl << m3.determinant() << endl << endl;
cout << "m3.norm() - 求模: " << endl << m3.norm() << endl << endl;
cout << "m3.sum() - 求和: " << endl << m3.sum() << endl << endl;
cout << "m3.mean() - 求均值: " << endl << m3.mean() << endl << endl;
cout << "m3.prod() - 求乘积: " << endl << m3.prod() << endl << endl;
cout << "m3.minCoeff() - 求最小: " << endl << m3.minCoeff() << endl << endl;
cout << "m3.maxCoeff() - 求最大: " << endl << m3.maxCoeff() << endl << endl;
// 3. 矩阵的运算2
Eigen::Vector3d v1(1, 2, 3);
Eigen::Vector3d v2(2, 1, 2);
// 点乘
double a;
a = v1.dot(v2);
cout << "v1.dot(v2) - 点乘: " << endl << a << endl << endl;
// 叉乘
cout << "v1.cross(v2) - 叉乘: " << endl << v1.cross(v2) << endl << endl;
四、Eigen的应用
1、求矩阵的特征值及特征向量
// 求解特征值及特征向量
Eigen::Matrix2f matrix2f;
matrix2f << 1, 2, 3, 4;
Eigen::SelfAdjointEigenSolver<Eigen::Matrix2f> eigenSolver(matrix2f);
cout<< "SelfAdjointEigenSolver: " << eigenSolver.info() << endl;
if (eigenSolver.info() == Eigen::Success ) {
cout << eigenSolver.eigenvalues() << endl << endl;
cout << eigenSolver.eigenvectors() << endl << endl;
}
2、求解线性方程组Ax=b
方法一:直接法求线性方程组(求逆法)
const int n = 4;
Eigen::Matrix<double, n, n > A;
A = A.Random() * 10;
Eigen::Matrix<double, n, 1> b;
b = b.Random() * 10;
// 定义未知x
Eigen::Matrix<double, n, 1> x;
x = A.inverse() * b; // 直接取逆 (这种效率很低,不推荐使用)
cout << "A=" << endl << A << endl << "b=" << b << endl << endl;
cout << "x=" << x << endl << endl;
注意上面例子,直接求逆法是低效的,一般情况不推荐选择求逆法,且直接对A.inverse(),可能会取逆失败,无法完成A.inverse()*b的运算操作
方法二:QR分解法(推荐)
// const int n = 4;
// Eigen::Matrix<double, n, n > A;
// A = A.Random() * 10;
// Eigen::Matrix<double, n, 1> b;
// b = b.Random() * 10;
// case-1:使用前面的A、b矩阵,定义未知x2
Eigen::Matrix<double, n, 1> x2;
x2 = A.colPivHouseholderQr().solve(b);
cout << "QR分解法" <<endl;
cout << "x=" << x2 << endl << endl;
// case-2:当然有时为了使用QR分解对象,也可以用下面方法计算,
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qrDec(A);
x = qrDec.solve(b);
// case-3 或者采用分离式分布构造
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qrDec2; // 只声明,不待入参,并不立即构造分解
qrDec2.compute(A); // 执行计算进行分解
x = qrDec2.solve(b); // 求解
上面三种使用方法,任取一种即可
3、如何评估或计算解的精度
// 继续前面的计算,我们可以采用下面的方式计算精度
double relative_error = (A*x - b).norm() / b.norm();
std::cout << "The relative error is:\n" << relative_error << std::endl;
4、求矩阵的秩
计算秩,通常是使用各类分解算法提供的rank()方法计算,Eigen中的许多分解算法,都提供了rank()求秩
Eigen::MatrixXd A(3,3);
A << 1, 2, 5,
2, 1, 4,
3, 0, 3;
cout << "A=" << endl << A << endl << endl;
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qrDec;
// 分解
qrDec.compute(A);
// 计算秩
cout<< "A rank = " << endl << qrDec.rank() << endl << endl;
任何算法的秩的计算根据设置的阈值来进行的,也就是设置不同的阈值,对结果将造成影响;Eigen中秩的计算依赖于分解方法的阈值,并不是矩阵对角大小乘以机器的精度epsilon,可通过调用
dec.setThreshold(th_val)
设置;且注意,分解算法并不会依赖阈值的设定,也就是执行完成分解后,任可以任意设置阈值,并求解rank()得到正确的秩;
Eigen库提供了多种求解线性方程组的算法,如下所示
参考地址,点这里
在实际问题处理中,如何选择不同的方法,对求解效率将会较大的影响,可根据系数矩阵A的特征,选择不同的分解算法:
- 如果系数矩阵是对于非对称、满秩的,则最适合的分解求解方法是PartialPivLU.
- 如果系数矩阵是对称正定的,则最适合的分解方法是 LLT 或LDLT
各类算法的求解效率对比可参考如下:
参考地址,点这里