(0)模块和头文件
模块 | 头文件 | 内容 |
---|---|---|
Core | #include | 矩阵和数组类,基本线性代数,数组操作 |
Geometry | #include | 变换、平移、缩放、2D旋转 和 3D 旋转(四元数、AngleAxis) |
LU | #include | 使用求解器求方阵的逆与行列式,以及对矩阵进行 LU 分解 (FullPivLU, PartialPivLU) |
Cholesky | #include | 使用求解器进行 LLT 和 LDLT Cholesky 分解 |
Householder | #include | Householder transformations; this module is used by several linear algebra modules |
SVD | #include | 使用最小二乘求解器进行 SVD 分解(JacobiSVD、BDCSVD) |
QR | #include | 使用求解器进行 QR 分解(HouseholderQR、ColPivHouseholderQR、FullPivHouseholderQR) |
Eigenvalues | #include | 特征值、特征向量分解(EigenSolver、SelfAdjointEigenSolver、ComplexEigenSolver) |
Sparse | #include | 稀疏矩阵存储及相关基础线性代数(SparseMatrix、SparseVector) |
#include | 包括 Core、Geometry、LU、Cholesky、SVD、QR 和 Eigenvalues 头文件 | |
#include | 包括 Dense 和 Sparse 头文件(整个 Eigen 库) |
(1)矩阵类
1)Matrix模板类定义:
Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>typedef Matrix<double, Dynamic, Dynamic> MatrixXd;typedef Matrix<int, Dynamic, 1> VectorXi;typedef Matrix<float, 3, 1> Vector3f;typedef Matrix<float, 4, 4> Matrix4f;
2)构造
Matrix3f a; //3*3MatrixXf b;MatrixXf a(10,15); //10*15VectorXf b(30); //30*1
固定大小的向量提供了初始化操作:
Vector2d a(5.0, 6.0);Vector3d b(5.0, 6.0, 7.0);Vector4d c(5.0, 6.0, 7.0, 8.0);
定义矩阵大小,但未初始化,该方式在堆上分配内存
MatrixXd m(2, 2); //2*2的矩阵// 矩阵元素赋值,index 从0开始m(0, 0) = 3;m(1, 0) = 2.5;m(0, 1) = -1;m(1, 1) = m(1, 0) + m(0, 1);
矩阵和矢量系数可以使用所谓的逗号初始化语法方便地设置
Matrix3f m;m << 1, 2, 3,4, 5, 6,7, 8, 9;
MatrixXi a { // construct a 2x2 matrix{1, 2}, // first row{3, 4} // second row};Matrix<double, 2, 3> b {{2, 3, 4},{5, 6, 7},};
3)动态调整
不能调整固定大小的矩阵的大小
MatrixXd m(2, 5);m.resize(4, 3);std::cout << "The matrix m is of size "<< m.rows() << "x" << m.cols() << std::endl;std::cout << "It has " << m.size() << " coefficients" << std::endl;VectorXd v(2);v.resize(5);std::cout << "The vector v is of size " << v.size() << std::endl;std::cout << "As a matrix, v is of size "<< v.rows() << "x" << v.cols() << std::endl;// Output is:// The matrix m is of size 4x3// It has 12 coefficients// The vector v is of size 5// As a matrix, v is of size 5x1
//operator= 将矩阵复制到另一个矩阵中的操作。Eigen自动调整左侧矩阵的大小,使其与右侧大小的矩阵大小匹配。例如:MatrixXf a(2, 2);std::cout << "a is of size " << a.rows() << "x" << a.cols() << std::endl;MatrixXf b(3, 3);a = b;std::cout << "a is now of size " << a.rows() << "x" << a.cols() << std::endl;// Output://a is of size 2x2//a is now of size 3x3
4)固定尺寸与动态尺寸
// 固定大小的矩阵or向量在栈上分配内存,因为它在编译时期就可以确定大小// Matrix4f mymatrix ≈ float mymatrix[16];// MatrixX表示运行时才确定矩阵的大小,因为它在堆上分配内存// MatrixXf mymatrix(rows,columns) ≈ float *mymatrix = new float[rows*columns];MatrixXd m = MatrixXd::Random(3, 3); // random values between -1 and 1m = (m + MatrixXd::Constant(3, 3, 1.2)) * 50; // MatrixXd::Constant(3, 3, 1.2) represents a 3-by-3 matrix expression having all coefficients equal to 1.2cout << "m =" << endl<< m << endl;VectorXd v(3); // 尚未初始化v << 1, 2, 3; // uses the so-called comma-initializercout << "m * v =" << endl<< m * v << endl;
什么时候应该使用固定尺寸(例如Matrix4f),什么时候应该使用动态尺寸(例如MatrixXf)?
简单的答案是:在可能的地方使用固定尺寸来显示非常小的尺寸,在需要的地方使用动态尺寸来显示较大的尺寸。
对于小尺寸,尤其是对于小于(大约)16的尺寸,使用固定尺寸对性能有极大的好处,因为它使Eigen避免了动态内存分配并展开了循环。
在内部,固定大小的本征矩阵只是一个简单的数组,即 Matrix4f mymatrix;
真的等于只是在做 float[16];
因此这确实具有零运行时间成本。相比之下,动态大小矩阵的数组始终分配在堆上,因此 MatrixXf mymatrix(行,列);
等于做 float * mymatrix = new [行*列];
除此之外,MatrixXf对象还将其行数和列数存储为成员变量。当然,使用固定大小的限制是,只有当您在编译时知道大小时,才有可能这样做。同样,对于足够大的尺寸(例如,对于大于(大约)32的尺寸),使用固定尺寸的性能优势变得可以忽略不计。更糟糕的是,尝试使用函数内部的固定大小创建非常大的矩阵可能会导致堆栈溢出,因为Eigen会尝试自动将数组分配为局部变量,而这通常是在堆栈上完成的。最后,视情况而定,当使用动态尺寸时,Eigen还可尝试进行矢量化(使用SIMD指令),请参见参考矢量化。
5)可选模板参数
在页面开始时提到Matrix类采用六个模板参数,但到目前为止,我们仅讨论了前三个。其余三个参数是可选的。这是模板参数的完整列表:
Matrix<typename Scalar,int RowsAtCompileTime,int ColsAtCompileTime,int Options = 0,int MaxRowsAtCompileTime = RowsAtCompileTime,int MaxColsAtCompileTime = ColsAtCompileTime>
Options是位字段。在这里,我们只讨论一点:RowMajor。它指定这种类型的矩阵使用行优先存储顺序。默认情况下,存储顺序为“按列的顺序存储”。例如,此类型表示行优先存储的3×3矩阵:Matrix
MaxRowsAtCompileTime并且MaxColsAtCompileTime在您希望指定时很有用,即使在编译时不知道矩阵的确切大小,在编译时也知道固定的上限。可能要这样做的最大原因是避免动态内存分配。例如,以下矩阵类型使用12个浮点数的普通数组,而不分配动态内存: Matrix
6)Matrix typedef
MatrixNt for Matrix
For example, MatrixXi for Matrix
VectorNt for Matrix
For example, Vector2f for Matrix
RowVectorNt for Matrix
For example, RowVector3d for Matrix
其中:
N可以是任何一个2,3,4,或X(意思Dynamic)。
t可以是i(表示int),f(表示float),d(表示double),cf(表示complex )或cd(表示complex )的任何一种。
仅针对这五个类型定义typedef的事实并不意味着它们是唯一受支持的标量类型。例如,支持所有标准整数类型,请参阅标量类型。
(2)矩阵与向量运算
1)加减运算
#include #include int main(){Eigen::Matrix2d a;a << 1, 2,3, 4;Eigen::MatrixXd b(2,2);b << 2, 3,1, 4;std::cout << "a + b =\n" << a + b << std::endl;std::cout << "a - b =\n" << a - b << std::endl;std::cout << "Doing a += b;" << std::endl;a += b;std::cout << "Now a =\n" << a << std::endl;Eigen::Vector3d v(1,2,3);Eigen::Vector3d w(1,0,0);std::cout << "-v + w - v =\n" << -v + w - v << std::endl;}
2)标量乘法除法
#include #include int main(){Eigen::Matrix2d a;a << 1, 2,3, 4;Eigen::Vector3d v(1,2,3);std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;std::cout << "0.1 * v =\n" << 0.1 * v << std::endl;std::cout << "Doing v *= 2;" << std::endl;v *= 2;std::cout << "Now v =\n" << v << std::endl;}
3)转置和共轭
MatrixXcf a = MatrixXcf::Random(2,2);cout << "Here is the matrix a\n" << a << endl;cout << "Here is the matrix a^T\n" << a.transpose() << endl; //转置矩阵cout << "Here is the conjugate of a\n" << a.conjugate() << endl; //共轭矩阵cout << "Here is the matrix a^*\n" << a.adjoint() << endl; //伴随矩阵
注意:对于一个矩阵自身的转置,应该使用.transposeInPlace()
Matrix2i a; a << 1, 2, 3, 4;cout << "Here is the matrix a:\n" << a << endl;//错误做法a = a.transpose(); // !!! do NOT do this !!!cout << "and the result of the aliasing effect:\n" << a << endl;//正确做法MatrixXf a(2,3);a << 1, 2, 3, 4, 5, 6;cout << "Here is the initial matrix a:\n" << a << endl;a.transposeInPlace();cout << "and after being transposed:\n" << a << endl;
4)矩阵乘法
#include #include int main(){Eigen::Matrix2d mat;mat << 1, 2,3, 4;Eigen::Vector2d u(-1,1), v(2,0);std::cout << "Here is mat*mat:\n" << mat*mat << std::endl;std::cout << "Here is mat*u:\n" << mat*u << std::endl;std::cout << "Here is u^T*mat:\n" << u.transpose()*mat << std::endl;std::cout << "Here is u^T*v:\n" << u.transpose()*v << std::endl;std::cout << "Here is u*v^T:\n" << u*v.transpose() << std::endl;std::cout << "Let's multiply mat by itself" << std::endl;mat = mat*mat;std::cout << "Now mat is mat:\n" << mat << std::endl;}
5)点乘和叉乘
叉积只适用于大小为3的向量,点积适用于任意向量
#include #include int main(){Eigen::Vector3d v(1,2,3);Eigen::Vector3d w(0,1,2);std::cout << "Dot product: " << v.dot(w) << std::endl;double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalarstd::cout << "Dot product via a matrix product: " << dp << std::endl;std::cout << "Cross product:\n" << v.cross(w) << std::endl;}
6)基本算术操作
#include #include using namespace std;int main(){Eigen::Matrix2d mat;mat << 1, 2,3, 4;cout << "Here is mat.sum(): " << mat.sum() << endl; //元素和cout << "Here is mat.prod(): " << mat.prod() << endl; //元素乘积cout << "Here is mat.mean(): " << mat.mean() << endl; //元素均值cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; //最小系数cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; //最大系数cout << "Here is mat.trace(): " << mat.trace() << endl; //迹}
可以输出元素位置
Matrix3f m = Matrix3f::Random();std::ptrdiff_t i, j;float minOfM = m.minCoeff(&i,&j);cout << "Here is the matrix m:\n" << m << endl;cout << "Its minimum coefficient (" << minOfM<< ") is at position (" << i << "," << j << ")\n\n";RowVector4i v = RowVector4i::Random();int maxOfV = v.maxCoeff(&i);cout << "Here is the vector v: " << v << endl;cout << "Its maximum coefficient (" << maxOfV<< ") is at position " << i << endl;
(3)Array类与系数操作
1)Array
与用于线性代数的Matrix类相反,Array类提供了通用数组。此外,Array类提供了一种执行逐系数运算的简便方法,该运算可能没有线性代数含义,例如将常数添加到数组中的每个系数或按系数乘两个数组。
Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
Eigen还为某些常见情况提供了typedef,其方式类似于Matrix typedef,但有一些细微的差异,因为单词 “array”既用于一维数组,也用于二维数组。
一维数组采用ArrayNt形式的typedef代表,其中N和t是大小和标量类型
二维数组采用ArrayNNt形式的typedef
Array<float,Dynamic,1> //ArrayXf Array<float,3,1> //Array3f Array<double,Dynamic,Dynamic> //ArrayXXd Array<double,3,3> //Array33d
2)访问数组中的值
#include #include int main(){Eigen::ArrayXXf m(2,2);// assign some values coefficient by coefficientm(0,0) = 1.0; m(0,1) = 2.0;m(1,0) = 3.0; m(1,1) = m(0,1) + m(1,0);// print values to standard outputstd::cout << m << std::endl << std::endl;// using the comma-initializer is also allowedm << 1.0,2.0,3.0,4.0;// print values to standard outputstd::cout << m << std::endl;}
3)加减法
两个数组的加减法与矩阵相同。如果两个数组的大小相同,并且该加法或减法是按系数进行的,则此操作有效。数组还支持以下形式的表达式,该表达式array + scalar将标量添加到数组中的每个系数。这提供了不能直接用于Matrix对象的功能。
#include #include int main(){Eigen::ArrayXXf a(3,3);Eigen::ArrayXXf b(3,3);a << 1,2,3,4,5,6,7,8,9;b << 1,2,3,1,2,3,1,2,3;// Adding two arraysstd::cout << "a + b = " << std::endl << a + b << std::endl << std::endl;// Subtracting a scalar from an arraystd::cout << "a - 2 = " << std::endl << a - 2 << std::endl;}
4)数组乘法
矩阵将乘法解释为矩阵乘积,而数组将乘法解释为按系数乘积
#include #include int main(){Eigen::ArrayXXf a(2,2);Eigen::ArrayXXf b(2,2);a << 1,2,3,4;b << 5,6,7,8;std::cout << "a * b = " << std::endl << a * b << std::endl;}// Output is:// a * b =// 5 12// 21 32
5) 其他运算
数组类定义除上述加法,减法和乘法运算符其他系数为单位的运算。例如,.abs()方法获取每个系数的绝对值,而.sqrt()计算系数的平方根。如果有两个大小相同的数组,则可以调用.min()来构造其系数是两个给定数组中对应系数的最小值的数组。注意:.array()和.matrix() 转换为相同的维数,但是不同的对象具有不同的方法
#include #include int main(){Eigen::ArrayXf a = Eigen::ArrayXf::Random(5);a *= 2;std::cout << "a =" << std::endl<< a << std::endl;std::cout << "a.abs() =" << std::endl<< a.abs() << std::endl;std::cout << "a.abs().sqrt() =" << std::endl<< a.abs().sqrt() << std::endl;std::cout << "a.min(a.abs().sqrt()) =" << std::endl<< a.min(a.abs().sqrt()) << std::endl;}// Output is:// a =// 1.36// -0.422// 1.13// 1.19// 1.65// a.abs() =// 1.36// 0.422// 1.13// 1.19// 1.65// a.abs().sqrt() =// 1.17// 0.65// 1.06// 1.09// 1.28// a.min(a.abs().sqrt()) =// 1.17// -0.422// 1.06// 1.09// 1.28
6)数组与矩阵之间的转换
什么时候应该使用Matrix类的对象,什么时候应该使用Array类的对象?
简单来说,不能对数组应用矩阵运算,也不能对矩阵应用数组运算。因此,如果需要进行线性代数运算(例如矩阵乘法),则应使用矩阵。如果需要进行系数运算,则应使用数组。但是,有时并不是那么简单,但是您需要同时使用Matrix和Array操作。在这种情况下,需要将矩阵转换为数组或反向转换。无论选择将对象声明为数组还是矩阵,都可以访问所有操作。
矩阵表达式具有.array()方法,可以将它们“转换”为数组表达式,因此可以轻松地应用按系数进行运算。相反,数组表达式具有.matrix()方法。与所有Eigen表达式抽象一样,这没有任何运行时开销(只要您让编译器进行优化).array()和.matrix() 可被用作右值和作为左值。
Eigen禁止在表达式中混合矩阵和数组。 例如,不能直接矩阵和数组相加。运算符+的操作数要么都是矩阵,要么都是数组。但是,使用.array()和.matrix()可以轻松地将其转换为另一种。注意,此规则的例外是赋值运算符=:允许将矩阵表达式分配给数组变量,或将数组表达式分配给矩阵变量。下面的示例演示如何通过使用.array()方法对Matrix对象使用数组操作。
例如,语句需要两个矩阵和,他们两个转换为阵列,用来将它们相乘系数明智并将结果指定给矩阵变量(这是合法的,因为本征允许分配数组表达式到矩阵的变量)。result = m.array() * n.array()mnresult。实际上,这种使用情况非常普遍,以至于Eigen为矩阵提供了const .cwiseProduct()方法来计算系数乘积。
#include #include using Eigen::MatrixXf;int main(){MatrixXf m(2,2);MatrixXf n(2,2);MatrixXf result(2,2);m << 1,2,3,4;n << 5,6,7,8;result = m * n;std::cout << "-- Matrix m*n: --\n" << result << "\n\n";result = m.array() * n.array();std::cout << "-- Array m*n: --\n" << result << "\n\n";result = m.cwiseProduct(n);std::cout << "-- With cwiseProduct: --\n" << result << "\n\n";result = m.array() + 4;std::cout << "-- Array m + 4: --\n" << result << "\n\n";}// Output is:// -- Matrix m*n: --// 19 22// 43 50// -- Array m*n: --// 5 12// 21 32// -- With cwiseProduct: --// 5 12// 21 32// -- Array m + 4: --// 5 6// 7 8
(4)块操作
1)使用块
块是矩阵或数组的某一部分。块表达式既可以用作右值,也可以用作左值。
Block of size (p,q), starting at (i,j)
1 起始位置+行列数 matrix.block(i,j,p,q);
2 模板参数为行列数,函数参数为起始位置 matrix.block
(i,j);
两种版本均可用于固定大小和动态大小的矩阵和数组。这两个表达式在语义上是等效的。唯一的区别是,如果块大小较小,则固定大小版本通常会为您提供更快的代码,但是需要在编译时知道此大小。
#include #include using namespace std;int main(){Eigen::MatrixXf m(4,4);m << 1, 2, 3, 4,5, 6, 7, 8,9,10,11,12,13,14,15,16;cout << "Block in the middle" << endl;cout << m.block<2,2>(1,1) << endl << endl;for (int i = 1; i <= 3; ++i){cout << "Block of size " << i << "x" << i << endl;cout << m.block(0,0,i,i) << endl << endl;}}
#include #include int main(){Eigen::Array22f m;m << 1,2,3,4;Eigen::Array44f a = Eigen::Array44f::Constant(0.6);std::cout << "Here is the array a:\n" << a << "\n\n";a.block<2,2>(1,1) = m; //可作为左值std::cout << "Here is now a with m copied into its central 2x2 block:\n" << a << "\n\n";a.block(0,0,2,3) = a.block(2,1,2,3);std::cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x3 block:\n" << a << "\n\n";}
2)列(columns)和行(rows)
#include #include using namespace std;int main(){Eigen::MatrixXf m(3,3);m << 1,2,3,4,5,6,7,8,9;cout << "Here is the matrix m:" << endl << m << endl;cout << "2nd Row: " << m.row(1) << endl;m.col(2) += 3 * m.col(0);cout << "After adding 3 times the first column into the third column, the matrix m is:\n";cout << m << endl;}
3)一系列矩阵的边角块操作
#include #include using namespace std;int main(){Eigen::Matrix4f m;m << 1, 2, 3, 4,5, 6, 7, 8,9, 10,11,12,13,14,15,16;cout << "m.leftCols(2) =" << endl << m.leftCols(2) << endl << endl; //前两列cout << "m.bottomRows() =" << endl << m.bottomRows<2>() << endl << endl; //下两行m.topLeftCorner(1,3) = m.bottomRightCorner(3,1).transpose();cout << "After assignment, m = " << endl << m << endl;}
4)向量的边角块操作
#include #include using namespace std;int main(){Eigen::ArrayXf v(6);v << 1, 2, 3, 4, 5, 6;cout << "v.head(3) =" << endl << v.head(3) << endl << endl;cout << "v.tail() = " << endl << v.tail<3>() << endl << endl;v.segment(1,4) *= 2;cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;}
(5)高级初始化方法
1)逗号初始化器
需要预先指定对象的大小。如果列出的系数太少或太多,就会报错。
Matrix3f m;m << 1, 2, 3,4, 5, 6,7, 8, 9;std::cout << m;
此外,初始化列表的元素本身可以是向量或矩阵。通常的用途是将向量或矩阵连接在一起。例如,这是将两个行向量连接在一起的方法。 注意:必须先设置大小,然后才能使用逗号初始化程序。
RowVectorXd vec1(3);vec1 << 1, 2, 3;std::cout << "vec1 = " << vec1 << std::endl;RowVectorXd vec2(4);vec2 << 1, 4, 9, 16;std::cout << "vec2 = " << vec2 << std::endl;RowVectorXd joined(7);joined << vec1, vec2;std::cout << "joined = " << joined << std::endl;
类似分块矩阵的初始化方式
MatrixXf matA(2, 2);matA << 1, 2, 3, 4;MatrixXf matB(4, 4);matB << matA, matA/10, matA/10, matA;std::cout << matB << std::endl;
对矩阵的某一块赋值
Matrix3f m;m.row(0) << 1, 2, 3;m.block(1,0,2,2) << 4, 5, 7, 8;m.col(2).tail(2) << 6, 9;std::cout << m;
2)特殊矩阵和数组
模板类Matrix和Array有静态方法,可以帮助初始化;
有三种变体:
第一个变体不带参数,只能用于固定大小的对象。如果要将动态尺寸对象初始化为零,则需要指定尺寸。
第二个变体需要一个参数,并且可以用于一维动态尺寸对象,
第三个变体需要两个参数,并且可以用于二维对象。
std::cout << "A fixed-size array:\n";Array33f a1 = Array33f::Zero();std::cout << a1 << "\n\n";std::cout << "A one-dimensional dynamic-size array:\n";ArrayXf a2 = ArrayXf::Zero(3);std::cout << a2 << "\n\n";std::cout << "A two-dimensional dynamic-size array:\n";ArrayXXf a3 = ArrayXXf::Zero(3, 4);std::cout << a3 << "\n";//Output//A fixed-size array://0 0 0//0 0 0//0 0 0//A one-dimensional dynamic-size array://0//0//0//A two-dimensional dynamic-size array://0 0 0 0//0 0 0 0//0 0 0 0
同样,静态方法Constant(value)会将所有系数设置为value。如果需要指定对象的大小,则附加参数放在value参数之前,如
MatrixXd::Constant(rows, cols, value)
。
Random()用随机系数填充矩阵或数组。Identity()获得单位矩阵, 此方法仅适用于Matrix,不适用于Array,因为“单位矩阵”是线性代数概念。该方法LinSpaced(尺寸,低,高)是仅可用于载体和一维数组; 它产生一个指定大小的向量,其系数在low和之间平均间隔high。方法LinSpaced()以下示例说明了该示例,该示例打印一张表格,其中包含以度为单位的角度,以弧度为单位的相应角度以及它们的正弦和余弦值。
ArrayXXf table(10, 4);table.col(0) = ArrayXf::LinSpaced(10, 0, 90);table.col(1) = M_PI / 180 * table.col(0);table.col(2) = table.col(1).sin();table.col(3) = table.col(1).cos();std::cout << " Degrees Radians Sine Cosine\n";std::cout << table << std::endl;
Eigen定义了诸如setZero()
,MatrixBase :: setIdentity()
和DenseBase :: setLinSpaced()
之类的实用程序函数来方便地执行此操作。即,可以采用对象的成员函数设置初始值。
const int size = 6;MatrixXd mat1(size, size);mat1.topLeftCorner(size/2, size/2) = MatrixXd::Zero(size/2, size/2);mat1.topRightCorner(size/2, size/2) = MatrixXd::Identity(size/2, size/2);mat1.bottomLeftCorner(size/2, size/2) = MatrixXd::Identity(size/2, size/2);mat1.bottomRightCorner(size/2, size/2) = MatrixXd::Zero(size/2, size/2);std::cout << mat1 << std::endl << std::endl;MatrixXd mat2(size, size);mat2.topLeftCorner(size/2, size/2).setZero();mat2.topRightCorner(size/2, size/2).setIdentity();mat2.bottomLeftCorner(size/2, size/2).setIdentity();mat2.bottomRightCorner(size/2, size/2).setZero();std::cout << mat2 << std::endl << std::endl;MatrixXd mat3(size, size);mat3 << MatrixXd::Zero(size/2, size/2), MatrixXd::Identity(size/2, size/2),MatrixXd::Identity(size/2, size/2), MatrixXd::Zero(size/2, size/2);std::cout << mat3 << std::endl;
(6)归约,访问者和广播
1)归约
在Eigen中,约简是一个采用矩阵或数组并返回单个标量值的函数。最常用的归约方法之一是.sum()
,它返回给定矩阵或数组内所有系数的总和。
#include #include using namespace std;int main(){Eigen::Matrix2d mat;mat << 1, 2,3, 4;cout << "Here is mat.sum(): " << mat.sum() << endl; //所有元素之和cout << "Here is mat.prod(): " << mat.prod() << endl; //所有元素乘积cout << "Here is mat.mean(): " << mat.mean() << endl; //和/个数cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; //最小值cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; //最大值cout << "Here is mat.trace(): " << mat.trace() << endl; //迹}
2)范数计算
2范数的平方可以通过squaredNorm()
获得。它本身等于矢量的点积,并且等效于其系数的平方绝对值的总和。norm()
方法返回squaredNorm()
的平方根。
这些运算也可以在矩阵上运算。在这种情况下,n×p矩阵被视为大小(n * p)的向量,因此,例如norm()方法返回“ Frobenius”或“ Hilbert-Schmidt”范数。
我们避免谈论l2矩阵的范数,因为那可能意味着不同的事情。如果需要其他按系数分配的lp范数,请使用lpNorm
()。如果需要无穷范数,则模板参数p可以采用特殊值Infinity,这是系数绝对值的最大值。
#include #include int main(){Eigen::VectorXf v(2);Eigen::MatrixXf m(2,2), n(2,2);v << -1,2;m << 1,-2,-3,4;//向量范数std::cout << "v.squaredNorm() = " << v.squaredNorm() << std::endl; //2-范数的平方std::cout << "v.norm() = " << v.norm() << std::endl; //2-范数std::cout << "v.lpNorm() = " << v.lpNorm<1>() << std::endl; //1-范数std::cout << "v.lpNorm() = " << v.lpNorm<Eigen::Infinity>() << std::endl; //无穷范数std::cout << std::endl;//矩阵范数std::cout << "m.squaredNorm() = " << m.squaredNorm() << std::endl; //m2-范数的平方std::cout << "m.norm() = " << m.norm() << std::endl; //m2-范数 std::cout << "m.lpNorm() = " << m.lpNorm<1>() << std::endl; //m1-范数std::cout << "m.lpNorm() = " << m.lpNorm<Eigen::Infinity>() << std::endl; //?既不是m无穷范数也不是无穷范数}
#include #include int main(){Eigen::MatrixXf m(2,2);m << 1,-2,-3,4;std::cout << "1-norm(m) = " << m.cwiseAbs().colwise().sum().maxCoeff()<< " == " << m.colwise().lpNorm<1>().maxCoeff() << std::endl; //1-范数(列和范数)std::cout << "infty-norm(m) = " << m.cwiseAbs().rowwise().sum().maxCoeff()<< " == " << m.rowwise().lpNorm<1>().maxCoeff() << std::endl; //无穷范数(行和范数)}
3)布尔
如果给定 Matrix 或 Array 中的所有系数都为 true ,则 all() 返回 true 。
如果给定 Matrix 或 Array 中的至少一个系数计算结果为 true ,则 any() 返回 true 。
count() 返回给定矩阵或数组中计算结果为真的系数的数量。
#include #include int main(){Eigen::ArrayXXf a(2,2);a << 1,2,3,4;std::cout << "(a > 0).all() = " << (a > 0).all() << std::endl;std::cout << "(a > 0).any() = " << (a > 0).any() << std::endl;std::cout << "(a > 0).count() = " << (a > 0).count() << std::endl;std::cout << std::endl;std::cout << "(a > 2).all() = " << (a > 2).all() << std::endl;std::cout << "(a > 2).any() = " << (a > 2).any() << std::endl;std::cout << "(a > 2).count() = " << (a > 2).count() << std::endl;}
4)访问者
在矩阵和数组的所有元素中,想要获得一个系数在Matrix或Array中的位置时,访问者很有用。最简单的示例是maxCoeff(&x,&y)和minCoeff(&x,&y),可用于查找Matrix或Array中最大或最小系数的位置。传递给访问者的参数是指向要存储行和列位置的变量的指针。这些变量应为Index类型。
#include #include <Eigen/Denseint main(){Eigen::MatrixXf m(2,2);m << 1, 2,3, 4;//get location of maximumEigen::Index maxRow, maxCol;float max = m.maxCoeff(&maxRow, &maxCol);//get location of minimumEigen::Index minRow, minCol;float min = m.minCoeff(&minRow, &minCol);std::cout << "Max: " << max << ", at: " <<maxRow << "," << maxCol << std::endl;std:: cout << "Min: " << min << ", at: " <<minRow << "," << minCol << std::endl;}
5)部分归约
在矩阵或数组的列向量和行向量中,element-wise是按元素的,那么colwise()或rowwise()表示按列或行的。部分归约是可以在Matrix或Array上按列或按行操作的归约,对每个列或行应用归约运算并返回具有相应值的列或行向量。部分缩减适用于colwis()或rowwise()
一个简单的示例是获取给定矩阵中每一列中元素的最大值,并将结果存储在行向量中:
#include #include using namespace std;int main(){Eigen::MatrixXf mat(2,4);mat << 1, 2, 6, 9,3, 1, 7, 2;std::cout << "Column's maximum: " << std::endl<< mat.colwise().maxCoeff() << std::endl; // 对于矩阵mat的每一列,取最大系数值}
#include #include using namespace std;int main(){Eigen::MatrixXf mat(2,4);mat << 1, 2, 6, 9,3, 1, 7, 2;std::cout << "Row's maximum: " << std::endl<< mat.rowwise().maxCoeff() << std::endl; // 对于矩阵mat的每一行,取最大系数值}
6)部分归约与其他操作结合
#include #include int main(){Eigen::MatrixXf mat(2,4);mat << 1, 2, 6, 9,3, 1, 7, 2;Eigen::Index maxIndex;float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex); //对于矩阵的每一列中的元素求和,结果的最大系数在第2列std::cout << "Maximum sum at position " << maxIndex << std::endl;std::cout << "The corresponding vector is: " << std::endl;std::cout << mat.col( maxIndex ) << std::endl;std::cout << "And its sum is is: " << maxNorm << std::endl;}
7)广播
#include #include using namespace std;int main(){//广播背后的概念类似于部分归约,区别在于广播构造了一个表达式,其中向量(列或行)通过在一个方向上复制而被解释为矩阵。//一个简单的示例是将某个列向量添加到矩阵中的每一列。这可以通过以下方式完成:Eigen::MatrixXf mat(2, 4);Eigen::VectorXf v(2);mat << 1, 2, 6, 9,3, 1, 7, 2;v << 0,1;//add v to each column of m//mat.colwise() += v用两种等效的方式解释指令。//它将向量添加v到矩阵的每一列。或者,可以将其解释为重复向量v四次以形成四乘二矩阵,然后将其加到matmat.colwise() += v;std::cout << "Broadcasting result: " << std::endl;std::cout << mat << std::endl;}
(7)Map类
1)Map类型与声明Map变量
Map 类 实现C++中的数组内存和Eigen对象的交互
Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >
在这种默认情况下,Map仅需要一个模板参数。
要构造Map变量,还需要其他两条信息:指向定义系数数组的内存区域的指针,以及所需的矩阵或矢量形状。(注意区分模板参数和函数形参)
例如,要定义一个float在编译时确定大小的矩阵,您可以执行以下操作:Map mf(pf,rows,columns);
其中pf是一个float类型的指针,指向内存中的数组。固定大小的整数只读向量可能会声明为Map mi(pi);
其中pi是int *。在这种情况下,不必将大小传递给构造函数,因为它已经由Matrix / Array类型指定。
Map足够灵活,可以容纳各种不同的数据表示形式。还有其他两个(可选)模板参数:Map
,MapOptions指定指针是Aligned还是Unaligned。默认值为Unaligned。StrideType允许使用Stride类为内存阵列指定自定义布局。一个示例是指定数据数组以行优先格式进行组织MapConstruct()。
int array[8];for(int i = 0; i < 8; ++i) array[i] = i;cout << "Column-major:\n" << Map<Matrix<int,2,4>>(array) << endl;cout << "Row-major:\n" << Map<Matrix<int,2,4,RowMajor>>(array) << endl;cout << "Row-major using stride:\n"<< Map<Matrix<int,2,4>, Unaligned, Stride<1,4>>(array) << endl;
2)使用Map参数
有点复杂了。。。
typedef Matrix<float,1,Dynamic> MatrixType;typedef Map<MatrixType> MapType;typedef Map<const MatrixType> MapTypeConst; // a read-only mapconst int n_dims = 5;MatrixType m1(n_dims), m2(n_dims);m1.setRandom();m2.setRandom();float *p = &m2(0); // get the address storing the data for m2MapType m2map(p,m2.size()); // m2map shares data with m2MapTypeConst m2mapconst(p,m2.size()); // a read-only accessor for m2cout << "m1: " << m1 << endl;cout << "m2: " << m2 << endl;cout << "Squared euclidean distance: " << (m1-m2).squaredNorm() << endl;cout << "Squared euclidean distance, using map: " <<(m1-m2map).squaredNorm() << endl;m2map(3) = 7; // this will change m2, since they share the same arraycout << "Updated m2: " << m2 << endl;cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << endl;/* m2mapconst(2) = 5; */ // this yields a compile-time error
(8)对齐问题
1)对齐错误
Eigen::internal::matrix_array::internal::matrix_array()[with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:Assertion `(reinterpret_cast(array) & (sizemask)) == 0 && "this assertionis explained here: http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.htmlREAD THIS WEB PAGE !!! ****"' failed.
2)四种原因
原因1:结构体中具有Eigen对象成员
注意,此处Eigen :: Vector2d仅用作示例,更一般而言,所有固定大小的可矢量化Eigen类型都会出现此问题,固定大小的可矢量化Eigen类型是如果它具有固定的大小并且大小是16字节的倍数。
Eigen::Vector2dEigen::Vector4dEigen::Vector4fEigen::Matrix2dEigen::Matrix2fEigen::Matrix4dEigen::Matrix4fEigen::Affine3dEigen::Affine3fEigen::QuaterniondEigen::Quaternionf
首先, “固定大小”应该清楚:如果在编译时,Eigen对象的行数和列数是固定的,则其固定大小。因此,例如Matrix3f具有固定大小,但MatrixXf没有(固定大小的对立是动态大小)。固定大小的Eigen对象的系数数组是普通的“静态数组”,不会动态分配。例如,Matrix4f后面的数据只是一个“float array[16]”。固定大小的对象通常很小,这意味着我们要以零的运行时开销(在内存使用和速度方面)来处理它们。现在,矢量化(SSE和AltiVec)都可以处理128位数据包。此外,出于性能原因,这些数据包必须具有128位对齐。因此,事实证明,固定大小的Eigen对象可以向量化的唯一方法是,如果它们的大小是128位或16个字节的倍数。然后,Eigen将为这些对象请求16字节对齐,并且此后将依赖于这些对象进行对齐,因此不会执行运行时检查以进行对齐。
class Foo{...Eigen::Vector2d v;...};...Foo *foo = new Foo;
Eigen需要Eigen :: Vector2d的数组(2个双精度)的128位对齐。对于GCC,这是通过属性((aligned(16)))完成的。Eigen重载了Eigen :: Vector2d的“ operator new”,因此它将始终返回128位对齐的指针。因此,通常情况下,您不必担心任何事情,Eigen会为您处理对齐。
除了一种情况。当您具有上述的Foo类,并且如上所述动态分配新的Foo时,则由于Foo没有对齐“ operator new”,因此返回的指针foo不一定是128位对齐的。然后,成员v的alignment属性相对于类的开头foo。如果foo指针未对齐,则foo-> v也不会对齐!解决方案是让Foo类具有一致的“Operator new”
解决方法:
如果定义的结构具有固定大小的可矢量化Eigen类型的成员,则必须重载其“ operator new”,以便它生成16字节对齐的指针。幸运的是,Eigen提供了一个宏EIGEN_MAKE_ALIGNED_OPERATOR_NEW
来执行此操作。换句话说:有一个类,该类具有固定大小的可矢量化Eigen对象作为成员,然后动态创建该类的对象。
class Foo{//...Eigen::Vector4d v;//...public:EIGEN_MAKE_ALIGNED_OPERATOR_NEW};//... Foo *foo = new Foo;
该宏使“ new Foo”始终返回对齐的指针。一个 Eigen::Vector2d有两个double类型,一个double为8字节=64位,则一个Eigen::Vector2d为128位,这恰好是SSE数据包的大小,这使得可以使用SSE对该向量执行各种操作。但是SSE指令(至少Eigen使用的是快速指令)需要128位对齐。否则会出现段错误。
原因2:STL容器或手动内存分配
原因3:通过值传递Eigen对象
原因4:编译器对堆栈对齐做出错误假设(例如Windows上的GCC)
此章节中还有少部分未写:Reshape,Slicing and Indexing,Aliasing,Storage orders,一是用不到,二是感觉今后的矩阵运算基本上都是线性代数的内容,出错率较低,错了之后再查也是可以的。