C++中线性代数计算Eigen库的使用教程详解
作者：yuanhao
前言
Eigen是一个基于线性代数的C++模板库，主要用于矩阵、向量、数值求解和相关算法。Eigen库有如下特点：
- 支持整数、浮点数、复数，使用模板编程，可以为特殊的数据结构提供矩阵操作；
- Open-CV自带Eigen的接口；
- 支持逐元素、分块和整体的矩阵操作。
- 支持使用Intel MKL加速部分功能。
- 支持多线程，对稀疏矩阵支持良好。
- 支持常用几何运算，包括旋转矩阵、矩阵变换等。
所以不论是做算法还是C++开发，使用好Eigen库会让你事半功倍。
由于本人是非专业的算法工程师，所以这里也只是做一些简单入门介绍和使用。
类型定义
Eigen库中的核心类就是Matrix，它表示矩阵，是一个模板类，定义如下：
template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
class Matrix {
}
其中前3个模板参数需要我们指定，后面3个参数使用默认即可。
第一个参数typename _Scalar，表示该参数是一个数据类型，而不是像_Rows和_Cols这样的变量。那么现在比如我想定义一个4x4的float矩阵，就可以如下定义：
Matrix<float,4,4> matrix44f;
这样定义虽然简单明了，但是有点长，所以在Eigen库中提供了很多开源直接使用的模板类，比如：
Matrix4f，表示4x4的float矩阵；Vector3f，表示列向量，长度3的float向量；RowVector2i，表示行向量，长度2的int向量；
所以对于一些常见矩阵可以使用内置的模板类。
除了定义已知大小的矩阵，Eigen库还可以定义动态矩阵。那什么是动态矩阵和静态矩阵呢？
静态矩阵就是编译时候就知道大小的矩阵，比如Matrix3d，在编译时我们就知道这是一个3x3的double类型矩阵。
与之对应的是动态矩阵，这种矩阵需要在运行后才知道大小，比如需要从vector<vector<double>>数据类型中构建矩阵，只有当代码运行时才可以知道构建的矩阵的行列个数。这时可以使用MatrixXd来表示任意大小的元素类型为double的矩阵变量，其中X就表示未知大小，即Dynamic动态的意思。
对于矩阵类型后面的后缀表示元素类型，Eigen库定义了4种元素类型：
- d：表示
double类型； - f：表示
float类型； - i：表示整数；
- c：表示复数；
比如Matrix2c，表示一个2x2元素类型为复数的矩阵。
新建矩阵
对于默认构造函数的矩阵，分配了大小和内存空间，但是没有初始化矩阵元素。而对于基本数据类型，默认初始化的话，其值是随机的，不能使用。比如下面代码：
Eigen::Matrix2f m;
std::cout << "m:" << std::endl << m << std::endl;
Eigen::Vector4i n;
std::cout << "n:" << std::endl << n << std::endl;
//运行结果：
m:
5.88604e-039 0
9.18341e-041 2.8026e-045
n:
-517829874
-2
1952350234
4199631
在Eigen中重载了std::cout的<<运算符，所以可以直接使用std::cout打印矩阵结果。所以定义完矩阵后，我们必须要对进行初始化。对于矩阵的初始化，有如下几种方式： 0. 直接赋值，可以使用<<进行赋值。
Eigen::Matrix2f m; //2x2的float矩阵
m << 1.9,2.884, //进行初始化，期望是4个值
4.33,5.898;
std::cout << "m:" << std::endl << m << std::endl;
Eigen::Vector4i n; //长度为4的列向量
n << 3,4,2,6;
std::cout << "n:" << std::endl << n << std::endl;
//运行结果：
m:
1.9 2.884
4.33 5.898
n:
3
4
2
6
这种方式只适合于矩阵的大小固定的情况，这种初始化方式很像使用列表初始化的数组，只有知道矩阵的大小，其构造函数才能从左向右、从上向下来进行初始化。比如可以定义动态大小的矩阵MatrixXd，这时就无法使用<<进行直接赋值，比如：
Eigen::MatrixXi xi; //动态大小的int矩阵
xi << 1,2,3,4,5,6;
std::cout << "xi:" << std::endl << xi << std::endl;
因为这种定义的变量xi，并不知道其行列数目，可能是3x2、2x3等等，所以用6个int去赋值，构造函数无法进行构造。
对于向量，可以在构造的时候进行初始化。
Eigen::RowVector3d k(1,2,3); //长度为3，类型为double的行向量
std::cout << "k:" << std::endl << k << std::endl;
一些特殊函数，比如全部初始化为0，初始化为1，随机数等。
Eigen::MatrixXi zero = MatrixXi::Zero(3,4); //全部初始化为0的矩阵 std::cout << "------ zero ------" << std::endl << zero << std::endl; Eigen::MatrixXi one = MatrixXi::Ones(4,4); //全部初始化为1的矩阵 std::cout << "------ one ------" << std::endl << one << std::endl; Eigen::MatrixXi i = MatrixXi::Identity(3,3); //单位矩阵 std::cout << "------ i ------" << std::endl << i << std::endl; Eigen::Matrix4f random = Matrix4f::Random(); //随机矩阵 std::cout << "------ random ------" << std::endl << random << std::endl;
运行结果：
------ zero ------
0 0 0 0
0 0 0 0
0 0 0 0
------ one ------
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
------ i ------
1 0 0
0 1 0
0 0 1
------ random ------
-0.997497 0.170019 0.64568 0.421003
0.127171 -0.0402539 0.49321 0.0270699
-0.613392 -0.299417 -0.651784 -0.39201
0.617481 0.791925 0.717887 -0.970031
使用Eigen::Map将已有的内存块数据映射到矩阵，必须要指定需要映射的矩阵类型，是MatrixXd、VectorXd等等，在模板参数中指明。
看一下Eigen::Map函数的定义：
template< e> class Map
: public MapBase<Map<PlainObjectType, MapOptions, StrideType> >
{
}
PlanObjectType：映射数据的等效矩阵类型。MapOptions：指定指针对齐方式，以字节为单位，默认为Unaligned，即未对齐。StrideType：可选择指定步幅。默认情况下，Map采用内存布局是一个普通的、连续的数组，可以通过指定步幅来定制化。
正常我们只需指明第一个参数即可，使用Map可以在无任何其他开销的情况下，让Eigen使用非Eigen数据结构来对其进行初始化，包括多种数据结构，比如原生数组，比如double[]、float[]等，以及STL容器，比如std::vector、std::array等。
std::vector<double> vec = {1,2,3,4,5,6}; //std::vector数据类型
Eigen::Map<Eigen::VectorXd> vd(vec.data(),6); //构造成一个长度为6的列向量
std::cout << "------ vd ------" << std::endl << vd << std::endl;
Eigen::Map<Eigen::MatrixXd> xd(vec.data(),3,2); //构造成一个3x2的矩阵
std::cout << "------ xd ------" << std::endl << xd << std::endl;
//运行结果：
------ vd ------
1
2
3
4
5
6
------ xd ------
1 4
2 5
3 6
在上面代码中，我们使用vec.data()来获取数vector中的数组数据，同时必须确定矩阵或者向量的大小。
对于矩阵来说，还有一种非常常见的初始化方式，就是一行一行或者一列一列进行赋值。比如代码：
Eigen::MatrixXd mat(3,2); //创建一个3x2的double类型矩阵 mat.col(0) << 1,2,3; //对第一列进行初始化 mat.col(1) << 4,5,6; //对第二列进行初始化
这种方式在构建特殊矩阵时经常用到。
矩阵索引
当前矩阵的行数、列数和大小可以分别通过rows()、cols()和size()方法来获取，在遍历Eigen矩阵时最好获取行数、列数来限制范围。
同时索引下标从0开始，矩阵元素的访问可以通过matrix(i,j)来访问第i行第j列的元素，对于向量还可以使用类似数组取值的vector[i]来访问第i个元素。比如代码：
std::vector<double> vec = {1,2,3,4,5,6};
Eigen::Map<Eigen::VectorXd> vd(vec.data(),6); //创建一个长度为6的列向量
std::cout << "------ vd ------" << std::endl << vd << std::endl;
Eigen::Map<Eigen::MatrixXd> xd(vec.data(),3,2); //创建一个3x2的矩阵
std::cout << "------ xd ------" << std::endl << xd << std::endl;
for (int i = 0; i < xd.rows(); ++i) {
for (int j = 0; j < xd.cols(); ++j) {
std::cout << "xd.(" << i << "," << j << ")= " << xd(i,j) << std::endl;
}
}
//运行结果：
------ xd ------
1 4
2 5
3 6
xd.(0,0)= 1
xd.(0,1)= 4
xd.(1,0)= 2
xd.(1,1)= 5
xd.(2,0)= 3
xd.(2,1)= 6
还可以利用block()函数来从Matrix中取出一个小矩阵来进行处理，语法是matrix:block<p,q>(i,j)，表示从原矩阵的(i,j)位置开始，截取出pxq大小的子矩阵，比如测试代码：
Eigen::MatrixXi randomI = MatrixXi::Random(6,6); //随机值的6x6矩阵
std::cout << "------ randomI ------" << std::endl << randomI << std::endl;
Eigen::MatrixXi smallI = randomI.block<3,3>(1,1); //从(1,1)处截取一个3x3的矩阵
std::cout << "------ smallI ------" << std::endl << smallI << std::endl;
//运行结果：
------ randomI ------
-13389 -12482 3334 -14515 12319 -1243
-4442 -16231 3511 3528 7427 -8673
-11557 -16092 -10937 9283 14938 11869
-10948 -4002 5342 9915 13949 -9516
16007 1037 -1613 651 1289 9163
-1780 2332 -4846 -6490 -11720 11260
------ smallI ------
-16231 3511 3528
-16092 -10937 9283
-4002 5342 9915
数学运算
Eigen库重载了常用加减乘运算符，而对于矩阵的除法，我们一般是求逆然后再转换为乘法，使用起来比较简单，就不过多介绍了。
还有几个常用的方法，必须要介绍一下，因为在项目中非常常见。
- 矩阵转置。调用
transpose()方法即可求出一个矩阵的转置矩阵。 - 矩阵求逆。调用
inverse()方法即可求出一个矩阵的逆矩阵。 - 共轭矩阵。调用
conjugate()方法可以求出共轭矩阵。
求解线性最小二乘方程组
此外，Eigen库的最主要一个作用可以直接求解最小二乘方程组，对于方程组Ax=b，如果没有确定解，可以找到一个向量x，使得其误差平方值最小。
我们可以回顾一下在线性回归中，如何解出最佳拟合系数\theta，一般都是先构建范德蒙矩阵，然后根据公式进行计算。但是在Eigen中，可以直接使用bdcSvd类的solve()方法直接求解出最小二乘方程组的系数。
比如下面代码中使用2种方法解最小二乘方程组：
Eigen::MatrixXf X(4, 3); //创建4x3的矩阵，表示有3个系数待求解，一共4个方程
Eigen::Vector4f Y; //表示结果的长度为4的列向量
X << 0, 0, 1,
1, 1, 1,
4, 2, 1,
9, 3, 1;
Y << 2, 0, 3, 7;
std::cout << "------ X ------" << std::endl << X << std::endl;
std::cout << "------ Y ------" << std::endl << Y << std::endl;
//使用bdcSvd方法
Eigen::MatrixXf abc = X.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Y);
std::cout << "------ abc ------" << std::endl << abc << std::endl;
//使用公式法
Eigen::MatrixXf abc1 = (X.transpose() * X).inverse() * (X.transpose()) * Y;
std::cout << "------ abc1 ------" << std::endl << abc1 << std::endl;
//运行结果：
------ X ------
0 0 1
1 1 1
4 2 1
9 3 1
------ Y ------
2
0
3
7
------ abc ------
1.5
-2.7
1.8
------ abc1 ------
1.5
-2.70001
1.8
到此这篇关于C++中线性代数计算Eigen库的使用教程详解的文章就介绍到这了,更多相关C++ Eigen库内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！