线性代数在3D游戏建模、深度学习、机器人工程、视频图像处理、航空航天等领域都有着举足轻重的位置。在使用C++编程时,如果从零编写矩阵的乘法、求逆、求特征值等运算,将会是件很麻烦的事情。而且,不管是从计算正确性,还是从效率方面都很难得到保障。因此,成熟稳定的线性代数库显得尤为重要。
Eigen是线性代数的C++模板库,它包含矩阵、向量、四元数、数值求解器和相关算法,其拥有良好的稳定性和计算效率,是我们搭建C++工程研究算法的必备利器。
通过官网下载,然后解压到自己指定的目录即可。在进行VS工程配置时,只需要如下图那样包含Eigen库所在目录即可。
右击工程名——》属性——》VC++目录——》常规——》包含目录
一、向量
常用的向量类型有 Vector3d Vector4d VectorXd VectorXcd等。其中VectorXd 表示实数向量,VectorXcd 表示复数向量,前者应用在许多科学计算和工程中,如物理模拟、信号处理等,后者在需要处理复数时应用,如电磁学、量子物理等。
1.1初始化
Vector3d v3d_a = { 1,2,3 }; cout << "Vector3d构造一:\n" << v3d_a << endl; Vector3d v3d_b; v3d_b << 4, 5, 6; cout << "Vector3d构造二:\n" << v3d_b << endl;
Vector4d v4d = { 1.1,2.2,3.3,4.4 }; cout << "Vector4d:\n" << v4d << endl; VectorXd vxd(8); //这里一定要指定向量长度,否则报错! vxd << 3, 1, 4, 1, 5, 9, 2, 6; cout << "VectorXd:\n" << vxd << endl;
初始化,给矩阵元素赋值,这里有两种方式
一种是"<<"操作符,用该操作符一个一个的赋值(注意顺序是从左到右,从上到下),或者一块一块的赋值。另外还有一些特殊矩阵(零矩阵,单位矩阵,常数矩阵)。
// 挨个赋值
Matrix3f m1;
m1 << 1, 2, 3,
4, 5, 6,
7, 8, 9;
cout << m1 << endl;
// 分块儿赋值,这对两个矩阵合并有用
int rows = 5, cols = 5;
MatrixXf m(rows, cols);
m << (Matrix3f() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished(),
MatrixXf::Zero(3, cols - 3),
MatrixXf::Zero(rows - 3, 3),
MatrixXf::Identity(rows - 3, cols - 3); // Zero()、Identity() 分别是0矩阵和单位矩阵。
cout << m << endl;
输出 :
一种是按下标赋值
int rows = 5, cols = 5;
MatrixXf m(rows,cols);
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
m(i, j) = 2;
}
}
cout << m << endl;
输出:
特殊矩阵:
Constant(rows,cols,value); //常量矩阵
Random();//随机
Zeor();//0
Identity();//单位
===================================
VectorXcd xcd(3, 1); xcd.setZero(); xcd[0]._Val[0] = 1.2; xcd[2]._Val[1] = 1.3; cout << "VectorXcd:\n" << xcd << endl; cout << "VectorXcd real:\n" << xcd.real() << endl;
1.2二范数
Vector3d v3d_a = { 1,2,3 }; cout << "二范数:\n" << v3d_a.norm() << endl;
1.3点积
Vector3d v3d_a = { 1,2,3 }; Vector3d v3d_b; v3d_b << 4, 5, 6; cout << "点积:\n" << v3d_a.dot(v3d_b) << endl;
1.4叉积
Vector3d v3d_a = { 1,2,3 }; Vector3d v3d_b; v3d_b << 4, 5, 6; cout << "叉积:\n" << v3d_a.cross(v3d_b) << endl;
二、矩阵
2.1单位矩阵
MatrixXd I = MatrixXd::Identity(4, 4); cout << "单位矩阵:\n" << I << endl;
2.2矩阵的迹
MatrixXd MatA(3, 3); MatA << 1, 2, 3, 4, 5, 6, 7, 8, 9; cout << "矩阵的迹:\n" << MatA.trace() << endl;
2.3转置
cout << "矩阵的转置:\n" << MatA.transpose() << endl;
2.4块操作
大小为 (p,q)
, 起始于 (i,j)
的块
动态大小的块表达式 matrix.block(i,j,p,q);
固定大小的块表达式matrix.block<p,q>(i,j);
二者区别是,固定大小的版本会在块比较小的时候快一点,但要求块大小在编译的时候就知道。
MatrixXd MatB(3, 3); Vector3d va = { 0.1,0.4,0.7 }; Vector3d vb = { 0.2,0.5,0.8 }; Vector3d vc = { 0.3,0.6,0.9 }; MatB.block<3, 1>(0, 0) = va; //固定 从0编码,以0,0为起点,存入3行1列的向量 MatB.block(0, 1, 3, 1) = vb; //动态 从0编码,以0,1为起点,存入3行1列的向量 MatB.block<3, 1>(0, 2) = vc; cout << "块操作:\n" << MatB << endl;
2.5乘法
cout << "矩阵乘法:\n" << MatA * MatB << endl;
2.6求逆
MatrixXd MatC(3, 3); MatC << 1, 2, 3, 2, 5, 6, 3, 4, 8; cout << "矩阵求逆:\n" << MatC.inverse() << endl;
2.7特征值
MatrixXd K(4, 4); K << -0.504129, 0.0190935, -0.836998, -0.137961, 0.0190935, -0.534845, 0.114843, -0.813407, -0.836998, 0.114843, 0.526649, -0.0167767, -0.137961, -0.813407, -0.0167767, 0.512325; VectorXd eigVals = K.eigenvalues().real(); cout << "特征值:\n" << eigVals << endl;
2.8特征向量
// 求解矩阵的最大特征值和对应特征向量 EigenSolver<MatrixXd> es(K); MatrixXcd eigenvectors = es.eigenvectors(); int index_of_max = MaxElementIndex(eigVals); // 找到最大特征值对应的下标 VectorXcd max_eigenvector = eigenvectors.col(index_of_max); Vector4d vec4 = max_eigenvector.real(); cout << "最大特征值对应的特征向量:\n" << vec4 << endl;
三、四元数
3.1初始化
Quaterniond q(4, 1, 2, 3);//要注意Eigen中四元数赋值的顺序,实数w在首; cout << "四元数:\n" << q << endl; //但是实际上它的内部存储顺序是[x y z w] cout << q.w() << " " << q.x() << " " << q.y() << " " << q.z() << endl; //单个取值
3.2旋转矩阵
q = vec4; cout << "旋转矩阵:\n" << q.toRotationMatrix() << endl; //将四元数转换为旋转矩阵
四、调试技巧
在用VS断点调试一般变量时,只需要鼠标悬停在变量上即可看见变量值,但由于矩阵和向量都是一个抽象对象,若要观察其具体变量就需要使用“快速监视“窗口,并且编辑表达式,使数据内容完全展开。
4.1向量的调试
双击选中变量
Shift+F9
一直展开到array,就可以看见向量全部内容
4.2矩阵的调试
同理,选中变量,Shift+F9,打开快速监视
此时展开到m_data变量,发现只能看到一个数值,这是因为该变量是个double指针类型
为了观察整个矩阵,在表达式栏,输入”,16”,表示展开指针后面16个数值,如下
敲击两次回车
认真比对存储顺序发现——矩阵元素存储按列依次排列,而非按行依次排列!
五、整体代码
#include <iostream> #include <Eigen/Core> #include <Eigen/Geometry> #include <Eigen/Eigenvalues> #include <Eigen/Dense> using namespace std; using namespace Eigen; int MaxElementIndex(VectorXd v) { int index = 0; for (int i = 1; i < v.size(); i++) { if (v(index) < v(i)) { index = i; } } return index; } int main() { Vector3d v3d_a = { 1,2,3 }; cout << "Vector3d构造一:\n" << v3d_a << endl; Vector3d v3d_b; v3d_b << 4, 5, 6; cout << "Vector3d构造二:\n" << v3d_b << endl; cout << "二范数:\n" << v3d_a.norm() << endl; cout << "点积:\n" << v3d_a.dot(v3d_b) << endl; cout << "叉积:\n" << v3d_a.cross(v3d_b) << endl; Vector4d v4d = { 1.1,2.2,3.3,4.4 }; cout << "Vector4d:\n" << v4d << endl; VectorXd vxd(8); //这里一定要指定向量长度,否则报错! vxd << 3, 1, 4, 1, 5, 9, 2, 6; cout << "VectorXd:\n" << vxd << endl; VectorXcd xcd(3, 1); xcd.setZero(); xcd[0]._Val[0] = 1.2; xcd[2]._Val[1] = 1.3; cout << "VectorXcd:\n" << xcd << endl; cout << "VectorXcd real:\n" << xcd.real() << endl; MatrixXd I = MatrixXd::Identity(4, 4); cout << "单位矩阵:\n" << I << endl; MatrixXd MatA(3, 3); MatA << 1, 2, 3, 4, 5, 6, 7, 8, 9; cout << "矩阵的迹:\n" << MatA.trace() << endl; cout << "矩阵的转置:\n" << MatA.transpose() << endl; MatrixXd MatB(3, 3); Vector3d va = { 0.1,0.4,0.7 }; Vector3d vb = { 0.2,0.5,0.8 }; Vector3d vc = { 0.3,0.6,0.9 }; MatB.block<3, 1>(0, 0) = va; //固定 从0编码,以0,0为起点,存入3行1列的向量 MatB.block(0, 1, 3, 1) = vb; //动态 从0编码,以0,1为起点,存入3行1列的向量 MatB.block<3, 1>(0, 2) = vc; cout << "块操作:\n" << MatB << endl; cout << "矩阵乘法:\n" << MatA * MatB << endl; MatrixXd MatC(3, 3); MatC << 1, 2, 3, 2, 5, 6, 3, 4, 8; cout << "矩阵求逆:\n" << MatC.inverse() << endl; MatrixXd K(4, 4); K << -0.504129, 0.0190935, -0.836998, -0.137961, 0.0190935, -0.534845, 0.114843, -0.813407, -0.836998, 0.114843, 0.526649, -0.0167767, -0.137961, -0.813407, -0.0167767, 0.512325; VectorXd eigVals = K.eigenvalues().real(); cout << "特征值:\n" << eigVals << endl; // 求解矩阵的最大特征值和对应特征向量 EigenSolver<MatrixXd> es(K); MatrixXcd eigenvectors = es.eigenvectors(); int index_of_max = MaxElementIndex(eigVals); // 找到最大特征值对应的下标 VectorXcd max_eigenvector = eigenvectors.col(index_of_max); Vector4d vec4 = max_eigenvector.real(); cout << "最大特征值对应的特征向量:\n" << vec4 << endl; Quaterniond q(4, 1, 2, 3);//要注意Eigen中四元数赋值的顺序,实数w在首; cout << "四元数:\n" << q << endl; //但是实际上它的内部存储顺序是[x y z w] cout << q.w() << " " << q.x() << " " << q.y() << " " << q.z() << endl; //单个取值 q = vec4; cout << "旋转矩阵:\n" << q.toRotationMatrix() << endl; //将四元数转换为旋转矩阵 }